granicus.if.org Git - imagemagick/blob - MagickCore/morphology.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %    M   M    OOO    RRRR   PPPP   H   H   OOO   L       OOO    GGGG  Y   Y   %
   7 %    MM MM   O   O   R   R  P   P  H   H  O   O  L      O   O  G       Y Y    %
   8 %    M M M   O   O   RRRR   PPPP   HHHHH  O   O  L      O   O  G GGG    Y     %
   9 %    M   M   O   O   R R    P      H   H  O   O  L      O   O  G   G    Y     %
  10 %    M   M    OOO    R  R   P      H   H   OOO   LLLLL   OOO    GGG     Y     %
  11 %                                                                             %
  12 %                                                                             %
  13 %                        MagickCore Morphology Methods                        %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                              Anthony Thyssen                                %
  17 %                               January 2010                                  %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2015 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 % Morpology is the the application of various kernels, of any size and even
  37 % shape, to a image in various ways (typically binary, but not always).
  38 %
  39 % Convolution (weighted sum or average) is just one specific type of
  40 % morphology. Just one that is very common for image bluring and sharpening
  41 % effects.  Not only 2D Gaussian blurring, but also 2-pass 1D Blurring.
  42 %
  43 % This module provides not only a general morphology function, and the ability
  44 % to apply more advanced or iterative morphologies, but also functions for the
  45 % generation of many different types of kernel arrays from user supplied
  46 % arguments. Prehaps even the generation of a kernel from a small image.
  47 */
  48 \f
  49 /*
  50   Include declarations.
  51 */
  52 #include "MagickCore/studio.h"
  53 #include "MagickCore/artifact.h"
  54 #include "MagickCore/cache-view.h"
  55 #include "MagickCore/channel.h"
  56 #include "MagickCore/color-private.h"
  57 #include "MagickCore/enhance.h"
  58 #include "MagickCore/exception.h"
  59 #include "MagickCore/exception-private.h"
  60 #include "MagickCore/gem.h"
  61 #include "MagickCore/gem-private.h"
  62 #include "MagickCore/hashmap.h"
  63 #include "MagickCore/image.h"
  64 #include "MagickCore/image-private.h"
  65 #include "MagickCore/list.h"
  66 #include "MagickCore/magick.h"
  67 #include "MagickCore/memory_.h"
  68 #include "MagickCore/memory-private.h"
  69 #include "MagickCore/monitor-private.h"
  70 #include "MagickCore/morphology.h"
  71 #include "MagickCore/morphology-private.h"
  72 #include "MagickCore/option.h"
  73 #include "MagickCore/pixel-accessor.h"
  74 #include "MagickCore/pixel-private.h"
  75 #include "MagickCore/prepress.h"
  76 #include "MagickCore/quantize.h"
  77 #include "MagickCore/resource_.h"
  78 #include "MagickCore/registry.h"
  79 #include "MagickCore/semaphore.h"
  80 #include "MagickCore/splay-tree.h"
  81 #include "MagickCore/statistic.h"
  82 #include "MagickCore/string_.h"
  83 #include "MagickCore/string-private.h"
  84 #include "MagickCore/thread-private.h"
  85 #include "MagickCore/token.h"
  86 #include "MagickCore/utility.h"
  87 #include "MagickCore/utility-private.h"
  88 \f
  89 /*
  90   Other global definitions used by module.
  91 */
  92 static inline double MagickMin(const double x,const double y)
  93 {
  94   return( x < y ? x : y);
  95 }
  96 static inline double MagickMax(const double x,const double y)
  97 {
  98   return( x > y ? x : y);
  99 }
 100 #define Minimize(assign,value) assign=MagickMin(assign,value)
 101 #define Maximize(assign,value) assign=MagickMax(assign,value)
 102
 103 /* Integer Factorial Function - for a Binomial kernel */
 104 #if 1
 105 static inline size_t fact(size_t n)
 106 {
 107   size_t f,l;
 108   for(f=1, l=2; l <= n; f=f*l, l++);
 109   return(f);
 110 }
 111 #elif 1 /* glibc floating point alternatives */
 112 #define fact(n) ((size_t)tgamma((double)n+1))
 113 #else
 114 #define fact(n) ((size_t)lgamma((double)n+1))
 115 #endif
 116
 117
 118 /* Currently these are only internal to this module */
 119 static void
 120   CalcKernelMetaData(KernelInfo *),
 121   ExpandMirrorKernelInfo(KernelInfo *),
 122   ExpandRotateKernelInfo(KernelInfo *, const double),
 123   RotateKernelInfo(KernelInfo *, double);
 124 \f
 125
 126 /* Quick function to find last kernel in a kernel list */
 127 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
 128 {
 129   while (kernel->next != (KernelInfo *) NULL)
 130     kernel=kernel->next;
 131   return(kernel);
 132 }
 133
 134 /*
 135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 136 %                                                                             %
 137 %                                                                             %
 138 %                                                                             %
 139 %     A c q u i r e K e r n e l I n f o                                       %
 140 %                                                                             %
 141 %                                                                             %
 142 %                                                                             %
 143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 144 %
 145 %  AcquireKernelInfo() takes the given string (generally supplied by the
 146 %  user) and converts it into a Morphology/Convolution Kernel.  This allows
 147 %  users to specify a kernel from a number of pre-defined kernels, or to fully
 148 %  specify their own kernel for a specific Convolution or Morphology
 149 %  Operation.
 150 %
 151 %  The kernel so generated can be any rectangular array of floating point
 152 %  values (doubles) with the 'control point' or 'pixel being affected'
 153 %  anywhere within that array of values.
 154 %
 155 %  Previously IM was restricted to a square of odd size using the exact
 156 %  center as origin, this is no longer the case, and any rectangular kernel
 157 %  with any value being declared the origin. This in turn allows the use of
 158 %  highly asymmetrical kernels.
 159 %
 160 %  The floating point values in the kernel can also include a special value
 161 %  known as 'nan' or 'not a number' to indicate that this value is not part
 162 %  of the kernel array. This allows you to shaped the kernel within its
 163 %  rectangular area. That is 'nan' values provide a 'mask' for the kernel
 164 %  shape.  However at least one non-nan value must be provided for correct
 165 %  working of a kernel.
 166 %
 167 %  The returned kernel should be freed using the DestroyKernelInfo() when you
 168 %  are finished with it.  Do not free this memory yourself.
 169 %
 170 %  Input kernel defintion strings can consist of any of three types.
 171 %
 172 %    "name:args[[@><]"
 173 %         Select from one of the built in kernels, using the name and
 174 %         geometry arguments supplied.  See AcquireKernelBuiltIn()
 175 %
 176 %    "WxH[+X+Y][@><]:num, num, num ..."
 177 %         a kernel of size W by H, with W*H floating point numbers following.
 178 %         the 'center' can be optionally be defined at +X+Y (such that +0+0
 179 %         is top left corner). If not defined the pixel in the center, for
 180 %         odd sizes, or to the immediate top or left of center for even sizes
 181 %         is automatically selected.
 182 %
 183 %    "num, num, num, num, ..."
 184 %         list of floating point numbers defining an 'old style' odd sized
 185 %         square kernel.  At least 9 values should be provided for a 3x3
 186 %         square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
 187 %         Values can be space or comma separated.  This is not recommended.
 188 %
 189 %  You can define a 'list of kernels' which can be used by some morphology
 190 %  operators A list is defined as a semi-colon separated list kernels.
 191 %
 192 %     " kernel ; kernel ; kernel ; "
 193 %
 194 %  Any extra ';' characters, at start, end or between kernel defintions are
 195 %  simply ignored.
 196 %
 197 %  The special flags will expand a single kernel, into a list of rotated
 198 %  kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
 199 %  cyclic rotations, while a '>' will generate a list of 90-degree rotations.
 200 %  The '<' also exands using 90-degree rotates, but giving a 180-degree
 201 %  reflected kernel before the +/- 90-degree rotations, which can be important
 202 %  for Thinning operations.
 203 %
 204 %  Note that 'name' kernels will start with an alphabetic character while the
 205 %  new kernel specification has a ':' character in its specification string.
 206 %  If neither is the case, it is assumed an old style of a simple list of
 207 %  numbers generating a odd-sized square kernel has been given.
 208 %
 209 %  The format of the AcquireKernal method is:
 210 %
 211 %      KernelInfo *AcquireKernelInfo(const char *kernel_string)
 212 %
 213 %  A description of each parameter follows:
 214 %
 215 %    o kernel_string: the Morphology/Convolution kernel wanted.
 216 %
 217 */
 218
 219 /* This was separated so that it could be used as a separate
 220 ** array input handling function, such as for -color-matrix
 221 */
 222 static KernelInfo *ParseKernelArray(const char *kernel_string)
 223 {
 224   KernelInfo
 225     *kernel;
 226
 227   char
 228     token[MaxTextExtent];
 229
 230   const char
 231     *p,
 232     *end;
 233
 234   register ssize_t
 235     i;
 236
 237   double
 238     nan = sqrt((double)-1.0);  /* Special Value : Not A Number */
 239
 240   MagickStatusType
 241     flags;
 242
 243   GeometryInfo
 244     args;
 245
 246   kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
 247   if (kernel == (KernelInfo *) NULL)
 248     return(kernel);
 249   (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
 250   kernel->minimum = kernel->maximum = kernel->angle = 0.0;
 251   kernel->negative_range = kernel->positive_range = 0.0;
 252   kernel->type = UserDefinedKernel;
 253   kernel->next = (KernelInfo *) NULL;
 254   kernel->signature = MagickSignature;
 255   if (kernel_string == (const char *) NULL)
 256     return(kernel);
 257
 258   /* find end of this specific kernel definition string */
 259   end = strchr(kernel_string, ';');
 260   if ( end == (char *) NULL )
 261     end = strchr(kernel_string, '\0');
 262
 263   /* clear flags - for Expanding kernel lists thorugh rotations */
 264    flags = NoValue;
 265
 266   /* Has a ':' in argument - New user kernel specification
 267      FUTURE: this split on ':' could be done by StringToken()
 268    */
 269   p = strchr(kernel_string, ':');
 270   if ( p != (char *) NULL && p < end)
 271     {
 272       /* ParseGeometry() needs the geometry separated! -- Arrgghh */
 273       memcpy(token, kernel_string, (size_t) (p-kernel_string));
 274       token[p-kernel_string] = '\0';
 275       SetGeometryInfo(&args);
 276       flags = ParseGeometry(token, &args);
 277
 278       /* Size handling and checks of geometry settings */
 279       if ( (flags & WidthValue) == 0 ) /* if no width then */
 280         args.rho = args.sigma;         /* then  width = height */
 281       if ( args.rho < 1.0 )            /* if width too small */
 282          args.rho = 1.0;               /* then  width = 1 */
 283       if ( args.sigma < 1.0 )          /* if height too small */
 284         args.sigma = args.rho;         /* then  height = width */
 285       kernel->width = (size_t)args.rho;
 286       kernel->height = (size_t)args.sigma;
 287
 288       /* Offset Handling and Checks */
 289       if ( args.xi  < 0.0 || args.psi < 0.0 )
 290         return(DestroyKernelInfo(kernel));
 291       kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
 292                                         : (ssize_t) (kernel->width-1)/2;
 293       kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
 294                                         : (ssize_t) (kernel->height-1)/2;
 295       if ( kernel->x >= (ssize_t) kernel->width ||
 296            kernel->y >= (ssize_t) kernel->height )
 297         return(DestroyKernelInfo(kernel));
 298
 299       p++; /* advance beyond the ':' */
 300     }
 301   else
 302     { /* ELSE - Old old specification, forming odd-square kernel */
 303       /* count up number of values given */
 304       p=(const char *) kernel_string;
 305       while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
 306         p++;  /* ignore "'" chars for convolve filter usage - Cristy */
 307       for (i=0; p < end; i++)
 308       {
 309         GetMagickToken(p,&p,token);
 310         if (*token == ',')
 311           GetMagickToken(p,&p,token);
 312       }
 313       /* set the size of the kernel - old sized square */
 314       kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
 315       kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
 316       p=(const char *) kernel_string;
 317       while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
 318         p++;  /* ignore "'" chars for convolve filter usage - Cristy */
 319     }
 320
 321   /* Read in the kernel values from rest of input string argument */
 322   kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
 323     kernel->width,kernel->height*sizeof(*kernel->values)));
 324   if (kernel->values == (MagickRealType *) NULL)
 325     return(DestroyKernelInfo(kernel));
 326   kernel->minimum=MagickMaximumValue;
 327   kernel->maximum=(-MagickMaximumValue);
 328   kernel->negative_range = kernel->positive_range = 0.0;
 329   for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
 330   {
 331     GetMagickToken(p,&p,token);
 332     if (*token == ',')
 333       GetMagickToken(p,&p,token);
 334     if (    LocaleCompare("nan",token) == 0
 335         || LocaleCompare("-",token) == 0 ) {
 336       kernel->values[i] = nan; /* this value is not part of neighbourhood */
 337     }
 338     else {
 339       kernel->values[i] = StringToDouble(token,(char **) NULL);
 340       ( kernel->values[i] < 0)
 341           ?  ( kernel->negative_range += kernel->values[i] )
 342           :  ( kernel->positive_range += kernel->values[i] );
 343       Minimize(kernel->minimum, kernel->values[i]);
 344       Maximize(kernel->maximum, kernel->values[i]);
 345     }
 346   }
 347
 348   /* sanity check -- no more values in kernel definition */
 349   GetMagickToken(p,&p,token);
 350   if ( *token != '\0' && *token != ';' && *token != '\'' )
 351     return(DestroyKernelInfo(kernel));
 352
 353 #if 0
 354   /* this was the old method of handling a incomplete kernel */
 355   if ( i < (ssize_t) (kernel->width*kernel->height) ) {
 356     Minimize(kernel->minimum, kernel->values[i]);
 357     Maximize(kernel->maximum, kernel->values[i]);
 358     for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
 359       kernel->values[i]=0.0;
 360   }
 361 #else
 362   /* Number of values for kernel was not enough - Report Error */
 363   if ( i < (ssize_t) (kernel->width*kernel->height) )
 364     return(DestroyKernelInfo(kernel));
 365 #endif
 366
 367   /* check that we recieved at least one real (non-nan) value! */
 368   if (kernel->minimum == MagickMaximumValue)
 369     return(DestroyKernelInfo(kernel));
 370
 371   if ( (flags & AreaValue) != 0 )         /* '@' symbol in kernel size */
 372     ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
 373   else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
 374     ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
 375   else if ( (flags & LessValue) != 0 )    /* '<' symbol in kernel args */
 376     ExpandMirrorKernelInfo(kernel);       /* 90 degree mirror rotate */
 377
 378   return(kernel);
 379 }
 380
 381 static KernelInfo *ParseKernelName(const char *kernel_string,
 382   ExceptionInfo *exception)
 383 {
 384   char
 385     token[MaxTextExtent];
 386
 387   const char
 388     *p,
 389     *end;
 390
 391   GeometryInfo
 392     args;
 393
 394   KernelInfo
 395     *kernel;
 396
 397   MagickStatusType
 398     flags;
 399
 400   ssize_t
 401     type;
 402
 403   /* Parse special 'named' kernel */
 404   GetMagickToken(kernel_string,&p,token);
 405   type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
 406   if ( type < 0 || type == UserDefinedKernel )
 407     return((KernelInfo *) NULL);  /* not a valid named kernel */
 408
 409   while (((isspace((int) ((unsigned char) *p)) != 0) ||
 410           (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
 411     p++;
 412
 413   end = strchr(p, ';'); /* end of this kernel defintion */
 414   if ( end == (char *) NULL )
 415     end = strchr(p, '\0');
 416
 417   /* ParseGeometry() needs the geometry separated! -- Arrgghh */
 418   memcpy(token, p, (size_t) (end-p));
 419   token[end-p] = '\0';
 420   SetGeometryInfo(&args);
 421   flags = ParseGeometry(token, &args);
 422
 423 #if 0
 424   /* For Debugging Geometry Input */
 425   (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
 426     flags, args.rho, args.sigma, args.xi, args.psi );
 427 #endif
 428
 429   /* special handling of missing values in input string */
 430   switch( type ) {
 431     /* Shape Kernel Defaults */
 432     case UnityKernel:
 433       if ( (flags & WidthValue) == 0 )
 434         args.rho = 1.0;    /* Default scale = 1.0, zero is valid */
 435       break;
 436     case SquareKernel:
 437     case DiamondKernel:
 438     case OctagonKernel:
 439     case DiskKernel:
 440     case PlusKernel:
 441     case CrossKernel:
 442       if ( (flags & HeightValue) == 0 )
 443         args.sigma = 1.0;    /* Default scale = 1.0, zero is valid */
 444       break;
 445     case RingKernel:
 446       if ( (flags & XValue) == 0 )
 447         args.xi = 1.0;       /* Default scale = 1.0, zero is valid */
 448       break;
 449     case RectangleKernel:    /* Rectangle - set size defaults */
 450       if ( (flags & WidthValue) == 0 ) /* if no width then */
 451         args.rho = args.sigma;         /* then  width = height */
 452       if ( args.rho < 1.0 )            /* if width too small */
 453           args.rho = 3;                /* then  width = 3 */
 454       if ( args.sigma < 1.0 )          /* if height too small */
 455         args.sigma = args.rho;         /* then  height = width */
 456       if ( (flags & XValue) == 0 )     /* center offset if not defined */
 457         args.xi = (double)(((ssize_t)args.rho-1)/2);
 458       if ( (flags & YValue) == 0 )
 459         args.psi = (double)(((ssize_t)args.sigma-1)/2);
 460       break;
 461     /* Distance Kernel Defaults */
 462     case ChebyshevKernel:
 463     case ManhattanKernel:
 464     case OctagonalKernel:
 465     case EuclideanKernel:
 466       if ( (flags & HeightValue) == 0 )           /* no distance scale */
 467         args.sigma = 100.0;                       /* default distance scaling */
 468       else if ( (flags & AspectValue ) != 0 )     /* '!' flag */
 469         args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
 470       else if ( (flags & PercentValue ) != 0 )    /* '%' flag */
 471         args.sigma *= QuantumRange/100.0;         /* percentage of color range */
 472       break;
 473     default:
 474       break;
 475   }
 476
 477   kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args, exception);
 478   if ( kernel == (KernelInfo *) NULL )
 479     return(kernel);
 480
 481   /* global expand to rotated kernel list - only for single kernels */
 482   if ( kernel->next == (KernelInfo *) NULL ) {
 483     if ( (flags & AreaValue) != 0 )         /* '@' symbol in kernel args */
 484       ExpandRotateKernelInfo(kernel, 45.0);
 485     else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
 486       ExpandRotateKernelInfo(kernel, 90.0);
 487     else if ( (flags & LessValue) != 0 )    /* '<' symbol in kernel args */
 488       ExpandMirrorKernelInfo(kernel);
 489   }
 490
 491   return(kernel);
 492 }
 493
 494 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string,
 495   ExceptionInfo *exception)
 496 {
 497   KernelInfo
 498     *kernel,
 499     *new_kernel;
 500
 501   char
 502     *kernel_cache,
 503     token[MaxTextExtent];
 504
 505   const char
 506     *p;
 507
 508   if (kernel_string == (const char *) NULL)
 509     return(ParseKernelArray(kernel_string));
 510   p=kernel_string;
 511   kernel_cache=(char *) NULL;
 512   if (*kernel_string == '@')
 513     {
 514       kernel_cache=FileToString(kernel_string+1,~0UL,exception);
 515       if (kernel_cache == (char *) NULL)
 516         return((KernelInfo *) NULL);
 517       p=(const char *) kernel_cache;
 518     }
 519   kernel=NULL;
 520   while (GetMagickToken(p,NULL,token), *token != '\0')
 521   {
 522     /* ignore extra or multiple ';' kernel separators */
 523     if (*token != ';')
 524       {
 525         /* tokens starting with alpha is a Named kernel */
 526         if (isalpha((int) ((unsigned char) *token)) != 0)
 527           new_kernel=ParseKernelName(p,exception);
 528         else /* otherwise a user defined kernel array */
 529           new_kernel=ParseKernelArray(p);
 530
 531         /* Error handling -- this is not proper error handling! */
 532         if (new_kernel == (KernelInfo *) NULL)
 533           {
 534             if (kernel != (KernelInfo *) NULL)
 535               kernel=DestroyKernelInfo(kernel);
 536             return((KernelInfo *) NULL);
 537           }
 538
 539         /* initialise or append the kernel list */
 540         if (kernel == (KernelInfo *) NULL)
 541           kernel=new_kernel;
 542         else
 543           LastKernelInfo(kernel)->next=new_kernel;
 544       }
 545
 546     /* look for the next kernel in list */
 547     p=strchr(p,';');
 548     if (p == (char *) NULL)
 549       break;
 550     p++;
 551   }
 552   if (kernel_cache != (char *) NULL)
 553     kernel_cache=DestroyString(kernel_cache);
 554   return(kernel);
 555 }
 556 \f
 557 /*
 558 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 559 %                                                                             %
 560 %                                                                             %
 561 %                                                                             %
 562 %     A c q u i r e K e r n e l B u i l t I n                                 %
 563 %                                                                             %
 564 %                                                                             %
 565 %                                                                             %
 566 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 567 %
 568 %  AcquireKernelBuiltIn() returned one of the 'named' built-in types of
 569 %  kernels used for special purposes such as gaussian blurring, skeleton
 570 %  pruning, and edge distance determination.
 571 %
 572 %  They take a KernelType, and a set of geometry style arguments, which were
 573 %  typically decoded from a user supplied string, or from a more complex
 574 %  Morphology Method that was requested.
 575 %
 576 %  The format of the AcquireKernalBuiltIn method is:
 577 %
 578 %      KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
 579 %           const GeometryInfo args)
 580 %
 581 %  A description of each parameter follows:
 582 %
 583 %    o type: the pre-defined type of kernel wanted
 584 %
 585 %    o args: arguments defining or modifying the kernel
 586 %
 587 %  Convolution Kernels
 588 %
 589 %    Unity
 590 %       The a No-Op or Scaling single element kernel.
 591 %
 592 %    Gaussian:{radius},{sigma}
 593 %       Generate a two-dimensional gaussian kernel, as used by -gaussian.
 594 %       The sigma for the curve is required.  The resulting kernel is
 595 %       normalized,
 596 %
 597 %       If 'sigma' is zero, you get a single pixel on a field of zeros.
 598 %
 599 %       NOTE: that the 'radius' is optional, but if provided can limit (clip)
 600 %       the final size of the resulting kernel to a square 2*radius+1 in size.
 601 %       The radius should be at least 2 times that of the sigma value, or
 602 %       sever clipping and aliasing may result.  If not given or set to 0 the
 603 %       radius will be determined so as to produce the best minimal error
 604 %       result, which is usally much larger than is normally needed.
 605 %
 606 %    LoG:{radius},{sigma}
 607 %        "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
 608 %        The supposed ideal edge detection, zero-summing kernel.
 609 %
 610 %        An alturnative to this kernel is to use a "DoG" with a sigma ratio of
 611 %        approx 1.6 (according to wikipedia).
 612 %
 613 %    DoG:{radius},{sigma1},{sigma2}
 614 %        "Difference of Gaussians" Kernel.
 615 %        As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
 616 %        from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
 617 %        The result is a zero-summing kernel.
 618 %
 619 %    Blur:{radius},{sigma}[,{angle}]
 620 %       Generates a 1 dimensional or linear gaussian blur, at the angle given
 621 %       (current restricted to orthogonal angles).  If a 'radius' is given the
 622 %       kernel is clipped to a width of 2*radius+1.  Kernel can be rotated
 623 %       by a 90 degree angle.
 624 %
 625 %       If 'sigma' is zero, you get a single pixel on a field of zeros.
 626 %
 627 %       Note that two convolutions with two "Blur" kernels perpendicular to
 628 %       each other, is equivalent to a far larger "Gaussian" kernel with the
 629 %       same sigma value, However it is much faster to apply. This is how the
 630 %       "-blur" operator actually works.
 631 %
 632 %    Comet:{width},{sigma},{angle}
 633 %       Blur in one direction only, much like how a bright object leaves
 634 %       a comet like trail.  The Kernel is actually half a gaussian curve,
 635 %       Adding two such blurs in opposite directions produces a Blur Kernel.
 636 %       Angle can be rotated in multiples of 90 degrees.
 637 %
 638 %       Note that the first argument is the width of the kernel and not the
 639 %       radius of the kernel.
 640 %
 641 %    Binomial:[{radius}]
 642 %       Generate a discrete kernel using a 2 dimentional Pascel's Triangle
 643 %       of values. Used for special forma of image filters.
 644 %
 645 %    # Still to be implemented...
 646 %    #
 647 %    # Filter2D
 648 %    # Filter1D
 649 %    #    Set kernel values using a resize filter, and given scale (sigma)
 650 %    #    Cylindrical or Linear.   Is this possible with an image?
 651 %    #
 652 %
 653 %  Named Constant Convolution Kernels
 654 %
 655 %  All these are unscaled, zero-summing kernels by default. As such for
 656 %  non-HDRI version of ImageMagick some form of normalization, user scaling,
 657 %  and biasing the results is recommended, to prevent the resulting image
 658 %  being 'clipped'.
 659 %
 660 %  The 3x3 kernels (most of these) can be circularly rotated in multiples of
 661 %  45 degrees to generate the 8 angled varients of each of the kernels.
 662 %
 663 %    Laplacian:{type}
 664 %      Discrete Lapacian Kernels, (without normalization)
 665 %        Type 0 :  3x3 with center:8 surounded by -1  (8 neighbourhood)
 666 %        Type 1 :  3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
 667 %        Type 2 :  3x3 with center:4 edge:1 corner:-2
 668 %        Type 3 :  3x3 with center:4 edge:-2 corner:1
 669 %        Type 5 :  5x5 laplacian
 670 %        Type 7 :  7x7 laplacian
 671 %        Type 15 : 5x5 LoG (sigma approx 1.4)
 672 %        Type 19 : 9x9 LoG (sigma approx 1.4)
 673 %
 674 %    Sobel:{angle}
 675 %      Sobel 'Edge' convolution kernel (3x3)
 676 %          | -1, 0, 1 |
 677 %          | -2, 0,-2 |
 678 %          | -1, 0, 1 |
 679 %
 680 %    Roberts:{angle}
 681 %      Roberts convolution kernel (3x3)
 682 %          |  0, 0, 0 |
 683 %          | -1, 1, 0 |
 684 %          |  0, 0, 0 |
 685 %
 686 %    Prewitt:{angle}
 687 %      Prewitt Edge convolution kernel (3x3)
 688 %          | -1, 0, 1 |
 689 %          | -1, 0, 1 |
 690 %          | -1, 0, 1 |
 691 %
 692 %    Compass:{angle}
 693 %      Prewitt's "Compass" convolution kernel (3x3)
 694 %          | -1, 1, 1 |
 695 %          | -1,-2, 1 |
 696 %          | -1, 1, 1 |
 697 %
 698 %    Kirsch:{angle}
 699 %      Kirsch's "Compass" convolution kernel (3x3)
 700 %          | -3,-3, 5 |
 701 %          | -3, 0, 5 |
 702 %          | -3,-3, 5 |
 703 %
 704 %    FreiChen:{angle}
 705 %      Frei-Chen Edge Detector is based on a kernel that is similar to
 706 %      the Sobel Kernel, but is designed to be isotropic. That is it takes
 707 %      into account the distance of the diagonal in the kernel.
 708 %
 709 %          |   1,     0,   -1     |
 710 %          | sqrt(2), 0, -sqrt(2) |
 711 %          |   1,     0,   -1     |
 712 %
 713 %    FreiChen:{type},{angle}
 714 %
 715 %      Frei-Chen Pre-weighted kernels...
 716 %
 717 %        Type 0:  default un-nomalized version shown above.
 718 %
 719 %        Type 1: Orthogonal Kernel (same as type 11 below)
 720 %          |   1,     0,   -1     |
 721 %          | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
 722 %          |   1,     0,   -1     |
 723 %
 724 %        Type 2: Diagonal form of Kernel...
 725 %          |   1,     sqrt(2),    0     |
 726 %          | sqrt(2),   0,     -sqrt(2) | / 2*sqrt(2)
 727 %          |   0,    -sqrt(2)    -1     |
 728 %
 729 %      However this kernel is als at the heart of the FreiChen Edge Detection
 730 %      Process which uses a set of 9 specially weighted kernel.  These 9
 731 %      kernels not be normalized, but directly applied to the image. The
 732 %      results is then added together, to produce the intensity of an edge in
 733 %      a specific direction.  The square root of the pixel value can then be
 734 %      taken as the cosine of the edge, and at least 2 such runs at 90 degrees
 735 %      from each other, both the direction and the strength of the edge can be
 736 %      determined.
 737 %
 738 %        Type 10: All 9 of the following pre-weighted kernels...
 739 %
 740 %        Type 11: |   1,     0,   -1     |
 741 %                 | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
 742 %                 |   1,     0,   -1     |
 743 %
 744 %        Type 12: | 1, sqrt(2), 1 |
 745 %                 | 0,   0,     0 | / 2*sqrt(2)
 746 %                 | 1, sqrt(2), 1 |
 747 %
 748 %        Type 13: | sqrt(2), -1,    0     |
 749 %                 |  -1,      0,    1     | / 2*sqrt(2)
 750 %                 |   0,      1, -sqrt(2) |
 751 %
 752 %        Type 14: |    0,     1, -sqrt(2) |
 753 %                 |   -1,     0,     1    | / 2*sqrt(2)
 754 %                 | sqrt(2), -1,     0    |
 755 %
 756 %        Type 15: | 0, -1, 0 |
 757 %                 | 1,  0, 1 | / 2
 758 %                 | 0, -1, 0 |
 759 %
 760 %        Type 16: |  1, 0, -1 |
 761 %                 |  0, 0,  0 | / 2
 762 %                 | -1, 0,  1 |
 763 %
 764 %        Type 17: |  1, -2,  1 |
 765 %                 | -2,  4, -2 | / 6
 766 %                 | -1, -2,  1 |
 767 %
 768 %        Type 18: | -2, 1, -2 |
 769 %                 |  1, 4,  1 | / 6
 770 %                 | -2, 1, -2 |
 771 %
 772 %        Type 19: | 1, 1, 1 |
 773 %                 | 1, 1, 1 | / 3
 774 %                 | 1, 1, 1 |
 775 %
 776 %      The first 4 are for edge detection, the next 4 are for line detection
 777 %      and the last is to add a average component to the results.
 778 %
 779 %      Using a special type of '-1' will return all 9 pre-weighted kernels
 780 %      as a multi-kernel list, so that you can use them directly (without
 781 %      normalization) with the special "-set option:morphology:compose Plus"
 782 %      setting to apply the full FreiChen Edge Detection Technique.
 783 %
 784 %      If 'type' is large it will be taken to be an actual rotation angle for
 785 %      the default FreiChen (type 0) kernel.  As such  FreiChen:45  will look
 786 %      like a  Sobel:45  but with 'sqrt(2)' instead of '2' values.
 787 %
 788 %      WARNING: The above was layed out as per
 789 %          http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
 790 %      But rotated 90 degrees so direction is from left rather than the top.
 791 %      I have yet to find any secondary confirmation of the above. The only
 792 %      other source found was actual source code at
 793 %          http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
 794 %      Neigher paper defineds the kernels in a way that looks locical or
 795 %      correct when taken as a whole.
 796 %
 797 %  Boolean Kernels
 798 %
 799 %    Diamond:[{radius}[,{scale}]]
 800 %       Generate a diamond shaped kernel with given radius to the points.
 801 %       Kernel size will again be radius*2+1 square and defaults to radius 1,
 802 %       generating a 3x3 kernel that is slightly larger than a square.
 803 %
 804 %    Square:[{radius}[,{scale}]]
 805 %       Generate a square shaped kernel of size radius*2+1, and defaulting
 806 %       to a 3x3 (radius 1).
 807 %
 808 %    Octagon:[{radius}[,{scale}]]
 809 %       Generate octagonal shaped kernel of given radius and constant scale.
 810 %       Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
 811 %       in "Diamond" kernel.
 812 %
 813 %    Disk:[{radius}[,{scale}]]
 814 %       Generate a binary disk, thresholded at the radius given, the radius
 815 %       may be a float-point value. Final Kernel size is floor(radius)*2+1
 816 %       square. A radius of 5.3 is the default.
 817 %
 818 %       NOTE: That a low radii Disk kernels produce the same results as
 819 %       many of the previously defined kernels, but differ greatly at larger
 820 %       radii.  Here is a table of equivalences...
 821 %          "Disk:1"    => "Diamond", "Octagon:1", or "Cross:1"
 822 %          "Disk:1.5"  => "Square"
 823 %          "Disk:2"    => "Diamond:2"
 824 %          "Disk:2.5"  => "Octagon"
 825 %          "Disk:2.9"  => "Square:2"
 826 %          "Disk:3.5"  => "Octagon:3"
 827 %          "Disk:4.5"  => "Octagon:4"
 828 %          "Disk:5.4"  => "Octagon:5"
 829 %          "Disk:6.4"  => "Octagon:6"
 830 %       All other Disk shapes are unique to this kernel, but because a "Disk"
 831 %       is more circular when using a larger radius, using a larger radius is
 832 %       preferred over iterating the morphological operation.
 833 %
 834 %    Rectangle:{geometry}
 835 %       Simply generate a rectangle of 1's with the size given. You can also
 836 %       specify the location of the 'control point', otherwise the closest
 837 %       pixel to the center of the rectangle is selected.
 838 %
 839 %       Properly centered and odd sized rectangles work the best.
 840 %
 841 %  Symbol Dilation Kernels
 842 %
 843 %    These kernel is not a good general morphological kernel, but is used
 844 %    more for highlighting and marking any single pixels in an image using,
 845 %    a "Dilate" method as appropriate.
 846 %
 847 %    For the same reasons iterating these kernels does not produce the
 848 %    same result as using a larger radius for the symbol.
 849 %
 850 %    Plus:[{radius}[,{scale}]]
 851 %    Cross:[{radius}[,{scale}]]
 852 %       Generate a kernel in the shape of a 'plus' or a 'cross' with
 853 %       a each arm the length of the given radius (default 2).
 854 %
 855 %       NOTE: "plus:1" is equivalent to a "Diamond" kernel.
 856 %
 857 %    Ring:{radius1},{radius2}[,{scale}]
 858 %       A ring of the values given that falls between the two radii.
 859 %       Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
 860 %       This is the 'edge' pixels of the default "Disk" kernel,
 861 %       More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
 862 %
 863 %  Hit and Miss Kernels
 864 %
 865 %    Peak:radius1,radius2
 866 %       Find any peak larger than the pixels the fall between the two radii.
 867 %       The default ring of pixels is as per "Ring".
 868 %    Edges
 869 %       Find flat orthogonal edges of a binary shape
 870 %    Corners
 871 %       Find 90 degree corners of a binary shape
 872 %    Diagonals:type
 873 %       A special kernel to thin the 'outside' of diagonals
 874 %    LineEnds:type
 875 %       Find end points of lines (for pruning a skeletion)
 876 %       Two types of lines ends (default to both) can be searched for
 877 %         Type 0: All line ends
 878 %         Type 1: single kernel for 4-conneected line ends
 879 %         Type 2: single kernel for simple line ends
 880 %    LineJunctions
 881 %       Find three line junctions (within a skeletion)
 882 %         Type 0: all line junctions
 883 %         Type 1: Y Junction kernel
 884 %         Type 2: Diagonal T Junction kernel
 885 %         Type 3: Orthogonal T Junction kernel
 886 %         Type 4: Diagonal X Junction kernel
 887 %         Type 5: Orthogonal + Junction kernel
 888 %    Ridges:type
 889 %       Find single pixel ridges or thin lines
 890 %         Type 1: Fine single pixel thick lines and ridges
 891 %         Type 2: Find two pixel thick lines and ridges
 892 %    ConvexHull
 893 %       Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
 894 %    Skeleton:type
 895 %       Traditional skeleton generating kernels.
 896 %         Type 1: Tradional Skeleton kernel (4 connected skeleton)
 897 %         Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
 898 %         Type 3: Thinning skeleton based on a ressearch paper by
 899 %                 Dan S. Bloomberg (Default Type)
 900 %    ThinSE:type
 901 %       A huge variety of Thinning Kernels designed to preserve conectivity.
 902 %       many other kernel sets use these kernels as source definitions.
 903 %       Type numbers are 41-49, 81-89, 481, and 482 which are based on
 904 %       the super and sub notations used in the source research paper.
 905 %
 906 %  Distance Measuring Kernels
 907 %
 908 %    Different types of distance measuring methods, which are used with the
 909 %    a 'Distance' morphology method for generating a gradient based on
 910 %    distance from an edge of a binary shape, though there is a technique
 911 %    for handling a anti-aliased shape.
 912 %
 913 %    See the 'Distance' Morphological Method, for information of how it is
 914 %    applied.
 915 %
 916 %    Chebyshev:[{radius}][x{scale}[%!]]
 917 %       Chebyshev Distance (also known as Tchebychev or Chessboard distance)
 918 %       is a value of one to any neighbour, orthogonal or diagonal. One why
 919 %       of thinking of it is the number of squares a 'King' or 'Queen' in
 920 %       chess needs to traverse reach any other position on a chess board.
 921 %       It results in a 'square' like distance function, but one where
 922 %       diagonals are given a value that is closer than expected.
 923 %
 924 %    Manhattan:[{radius}][x{scale}[%!]]
 925 %       Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
 926 %       Cab distance metric), it is the distance needed when you can only
 927 %       travel in horizontal or vertical directions only.  It is the
 928 %       distance a 'Rook' in chess would have to travel, and results in a
 929 %       diamond like distances, where diagonals are further than expected.
 930 %
 931 %    Octagonal:[{radius}][x{scale}[%!]]
 932 %       An interleving of Manhatten and Chebyshev metrics producing an
 933 %       increasing octagonally shaped distance.  Distances matches those of
 934 %       the "Octagon" shaped kernel of the same radius.  The minimum radius
 935 %       and default is 2, producing a 5x5 kernel.
 936 %
 937 %    Euclidean:[{radius}][x{scale}[%!]]
 938 %       Euclidean distance is the 'direct' or 'as the crow flys' distance.
 939 %       However by default the kernel size only has a radius of 1, which
 940 %       limits the distance to 'Knight' like moves, with only orthogonal and
 941 %       diagonal measurements being correct.  As such for the default kernel
 942 %       you will get octagonal like distance function.
 943 %
 944 %       However using a larger radius such as "Euclidean:4" you will get a
 945 %       much smoother distance gradient from the edge of the shape. Especially
 946 %       if the image is pre-processed to include any anti-aliasing pixels.
 947 %       Of course a larger kernel is slower to use, and not always needed.
 948 %
 949 %    The first three Distance Measuring Kernels will only generate distances
 950 %    of exact multiples of {scale} in binary images. As such you can use a
 951 %    scale of 1 without loosing any information.  However you also need some
 952 %    scaling when handling non-binary anti-aliased shapes.
 953 %
 954 %    The "Euclidean" Distance Kernel however does generate a non-integer
 955 %    fractional results, and as such scaling is vital even for binary shapes.
 956 %
 957 */
 958
 959 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
 960   const GeometryInfo *args,ExceptionInfo *exception)
 961 {
 962   KernelInfo
 963     *kernel;
 964
 965   register ssize_t
 966     i;
 967
 968   register ssize_t
 969     u,
 970     v;
 971
 972   double
 973     nan = sqrt((double)-1.0);  /* Special Value : Not A Number */
 974
 975   /* Generate a new empty kernel if needed */
 976   kernel=(KernelInfo *) NULL;
 977   switch(type) {
 978     case UndefinedKernel:    /* These should not call this function */
 979     case UserDefinedKernel:
 980       assert("Should not call this function" != (char *) NULL);
 981       break;
 982     case LaplacianKernel:   /* Named Descrete Convolution Kernels */
 983     case SobelKernel:       /* these are defined using other kernels */
 984     case RobertsKernel:
 985     case PrewittKernel:
 986     case CompassKernel:
 987     case KirschKernel:
 988     case FreiChenKernel:
 989     case EdgesKernel:       /* Hit and Miss kernels */
 990     case CornersKernel:
 991     case DiagonalsKernel:
 992     case LineEndsKernel:
 993     case LineJunctionsKernel:
 994     case RidgesKernel:
 995     case ConvexHullKernel:
 996     case SkeletonKernel:
 997     case ThinSEKernel:
 998       break;               /* A pre-generated kernel is not needed */
 999 #if 0
1000     /* set to 1 to do a compile-time check that we haven't missed anything */
1001     case UnityKernel:
1002     case GaussianKernel:
1003     case DoGKernel:
1004     case LoGKernel:
1005     case BlurKernel:
1006     case CometKernel:
1007     case BinomialKernel:
1008     case DiamondKernel:
1009     case SquareKernel:
1010     case RectangleKernel:
1011     case OctagonKernel:
1012     case DiskKernel:
1013     case PlusKernel:
1014     case CrossKernel:
1015     case RingKernel:
1016     case PeaksKernel:
1017     case ChebyshevKernel:
1018     case ManhattanKernel:
1019     case OctangonalKernel:
1020     case EuclideanKernel:
1021 #else
1022     default:
1023 #endif
1024       /* Generate the base Kernel Structure */
1025       kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1026       if (kernel == (KernelInfo *) NULL)
1027         return(kernel);
1028       (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1029       kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1030       kernel->negative_range = kernel->positive_range = 0.0;
1031       kernel->type = type;
1032       kernel->next = (KernelInfo *) NULL;
1033       kernel->signature = MagickSignature;
1034       break;
1035   }
1036
1037   switch(type) {
1038     /*
1039       Convolution Kernels
1040     */
1041     case UnityKernel:
1042       {
1043         kernel->height = kernel->width = (size_t) 1;
1044         kernel->x = kernel->y = (ssize_t) 0;
1045         kernel->values=(MagickRealType *) MagickAssumeAligned(
1046           AcquireAlignedMemory(1,sizeof(*kernel->values)));
1047         if (kernel->values == (MagickRealType *) NULL)
1048           return(DestroyKernelInfo(kernel));
1049         kernel->maximum = kernel->values[0] = args->rho;
1050         break;
1051       }
1052       break;
1053     case GaussianKernel:
1054     case DoGKernel:
1055     case LoGKernel:
1056       { double
1057           sigma = fabs(args->sigma),
1058           sigma2 = fabs(args->xi),
1059           A, B, R;
1060
1061         if ( args->rho >= 1.0 )
1062           kernel->width = (size_t)args->rho*2+1;
1063         else if ( (type != DoGKernel) || (sigma >= sigma2) )
1064           kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1065         else
1066           kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1067         kernel->height = kernel->width;
1068         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1069         kernel->values=(MagickRealType *) MagickAssumeAligned(
1070           AcquireAlignedMemory(kernel->width,kernel->height*
1071           sizeof(*kernel->values)));
1072         if (kernel->values == (MagickRealType *) NULL)
1073           return(DestroyKernelInfo(kernel));
1074
1075         /* WARNING: The following generates a 'sampled gaussian' kernel.
1076          * What we really want is a 'discrete gaussian' kernel.
1077          *
1078          * How to do this is I don't know, but appears to be basied on the
1079          * Error Function 'erf()' (intergral of a gaussian)
1080          */
1081
1082         if ( type == GaussianKernel || type == DoGKernel )
1083           { /* Calculate a Gaussian,  OR positive half of a DoG */
1084             if ( sigma > MagickEpsilon )
1085               { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
1086                 B = (double) (1.0/(Magick2PI*sigma*sigma));
1087                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1088                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1089                       kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1090               }
1091             else /* limiting case - a unity (normalized Dirac) kernel */
1092               { (void) ResetMagickMemory(kernel->values,0, (size_t)
1093                   kernel->width*kernel->height*sizeof(*kernel->values));
1094                 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1095               }
1096           }
1097
1098         if ( type == DoGKernel )
1099           { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1100             if ( sigma2 > MagickEpsilon )
1101               { sigma = sigma2;                /* simplify loop expressions */
1102                 A = 1.0/(2.0*sigma*sigma);
1103                 B = (double) (1.0/(Magick2PI*sigma*sigma));
1104                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1105                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1106                     kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1107               }
1108             else /* limiting case - a unity (normalized Dirac) kernel */
1109               kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1110           }
1111
1112         if ( type == LoGKernel )
1113           { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1114             if ( sigma > MagickEpsilon )
1115               { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
1116                 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1117                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1118                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1119                     { R = ((double)(u*u+v*v))*A;
1120                       kernel->values[i] = (1-R)*exp(-R)*B;
1121                     }
1122               }
1123             else /* special case - generate a unity kernel */
1124               { (void) ResetMagickMemory(kernel->values,0, (size_t)
1125                   kernel->width*kernel->height*sizeof(*kernel->values));
1126                 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1127               }
1128           }
1129
1130         /* Note the above kernels may have been 'clipped' by a user defined
1131         ** radius, producing a smaller (darker) kernel.  Also for very small
1132         ** sigma's (> 0.1) the central value becomes larger than one, and thus
1133         ** producing a very bright kernel.
1134         **
1135         ** Normalization will still be needed.
1136         */
1137
1138         /* Normalize the 2D Gaussian Kernel
1139         **
1140         ** NB: a CorrelateNormalize performs a normal Normalize if
1141         ** there are no negative values.
1142         */
1143         CalcKernelMetaData(kernel);  /* the other kernel meta-data */
1144         ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1145
1146         break;
1147       }
1148     case BlurKernel:
1149       { double
1150           sigma = fabs(args->sigma),
1151           alpha, beta;
1152
1153         if ( args->rho >= 1.0 )
1154           kernel->width = (size_t)args->rho*2+1;
1155         else
1156           kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1157         kernel->height = 1;
1158         kernel->x = (ssize_t) (kernel->width-1)/2;
1159         kernel->y = 0;
1160         kernel->negative_range = kernel->positive_range = 0.0;
1161         kernel->values=(MagickRealType *) MagickAssumeAligned(
1162           AcquireAlignedMemory(kernel->width,kernel->height*
1163           sizeof(*kernel->values)));
1164         if (kernel->values == (MagickRealType *) NULL)
1165           return(DestroyKernelInfo(kernel));
1166
1167 #if 1
1168 #define KernelRank 3
1169         /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1170         ** It generates a gaussian 3 times the width, and compresses it into
1171         ** the expected range.  This produces a closer normalization of the
1172         ** resulting kernel, especially for very low sigma values.
1173         ** As such while wierd it is prefered.
1174         **
1175         ** I am told this method originally came from Photoshop.
1176         **
1177         ** A properly normalized curve is generated (apart from edge clipping)
1178         ** even though we later normalize the result (for edge clipping)
1179         ** to allow the correct generation of a "Difference of Blurs".
1180         */
1181
1182         /* initialize */
1183         v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1184         (void) ResetMagickMemory(kernel->values,0, (size_t)
1185           kernel->width*kernel->height*sizeof(*kernel->values));
1186         /* Calculate a Positive 1D Gaussian */
1187         if ( sigma > MagickEpsilon )
1188           { sigma *= KernelRank;               /* simplify loop expressions */
1189             alpha = 1.0/(2.0*sigma*sigma);
1190             beta= (double) (1.0/(MagickSQ2PI*sigma ));
1191             for ( u=-v; u <= v; u++) {
1192               kernel->values[(u+v)/KernelRank] +=
1193                               exp(-((double)(u*u))*alpha)*beta;
1194             }
1195           }
1196         else /* special case - generate a unity kernel */
1197           kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1198 #else
1199         /* Direct calculation without curve averaging
1200            This is equivelent to a KernelRank of 1 */
1201
1202         /* Calculate a Positive Gaussian */
1203         if ( sigma > MagickEpsilon )
1204           { alpha = 1.0/(2.0*sigma*sigma);    /* simplify loop expressions */
1205             beta = 1.0/(MagickSQ2PI*sigma);
1206             for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1207               kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1208           }
1209         else /* special case - generate a unity kernel */
1210           { (void) ResetMagickMemory(kernel->values,0, (size_t)
1211               kernel->width*kernel->height*sizeof(*kernel->values));
1212             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1213           }
1214 #endif
1215         /* Note the above kernel may have been 'clipped' by a user defined
1216         ** radius, producing a smaller (darker) kernel.  Also for very small
1217         ** sigma's (> 0.1) the central value becomes larger than one, as a
1218         ** result of not generating a actual 'discrete' kernel, and thus
1219         ** producing a very bright 'impulse'.
1220         **
1221         ** Becuase of these two factors Normalization is required!
1222         */
1223
1224         /* Normalize the 1D Gaussian Kernel
1225         **
1226         ** NB: a CorrelateNormalize performs a normal Normalize if
1227         ** there are no negative values.
1228         */
1229         CalcKernelMetaData(kernel);  /* the other kernel meta-data */
1230         ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1231
1232         /* rotate the 1D kernel by given angle */
1233         RotateKernelInfo(kernel, args->xi );
1234         break;
1235       }
1236     case CometKernel:
1237       { double
1238           sigma = fabs(args->sigma),
1239           A;
1240
1241         if ( args->rho < 1.0 )
1242           kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1243         else
1244           kernel->width = (size_t)args->rho;
1245         kernel->x = kernel->y = 0;
1246         kernel->height = 1;
1247         kernel->negative_range = kernel->positive_range = 0.0;
1248         kernel->values=(MagickRealType *) MagickAssumeAligned(
1249           AcquireAlignedMemory(kernel->width,kernel->height*
1250           sizeof(*kernel->values)));
1251         if (kernel->values == (MagickRealType *) NULL)
1252           return(DestroyKernelInfo(kernel));
1253
1254         /* A comet blur is half a 1D gaussian curve, so that the object is
1255         ** blurred in one direction only.  This may not be quite the right
1256         ** curve to use so may change in the future. The function must be
1257         ** normalised after generation, which also resolves any clipping.
1258         **
1259         ** As we are normalizing and not subtracting gaussians,
1260         ** there is no need for a divisor in the gaussian formula
1261         **
1262         ** It is less comples
1263         */
1264         if ( sigma > MagickEpsilon )
1265           {
1266 #if 1
1267 #define KernelRank 3
1268             v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1269             (void) ResetMagickMemory(kernel->values,0, (size_t)
1270               kernel->width*sizeof(*kernel->values));
1271             sigma *= KernelRank;            /* simplify the loop expression */
1272             A = 1.0/(2.0*sigma*sigma);
1273             /* B = 1.0/(MagickSQ2PI*sigma); */
1274             for ( u=0; u < v; u++) {
1275               kernel->values[u/KernelRank] +=
1276                   exp(-((double)(u*u))*A);
1277               /*  exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1278             }
1279             for (i=0; i < (ssize_t) kernel->width; i++)
1280               kernel->positive_range += kernel->values[i];
1281 #else
1282             A = 1.0/(2.0*sigma*sigma);     /* simplify the loop expression */
1283             /* B = 1.0/(MagickSQ2PI*sigma); */
1284             for ( i=0; i < (ssize_t) kernel->width; i++)
1285               kernel->positive_range +=
1286                 kernel->values[i] = exp(-((double)(i*i))*A);
1287                 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1288 #endif
1289           }
1290         else /* special case - generate a unity kernel */
1291           { (void) ResetMagickMemory(kernel->values,0, (size_t)
1292               kernel->width*kernel->height*sizeof(*kernel->values));
1293             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1294             kernel->positive_range = 1.0;
1295           }
1296
1297         kernel->minimum = 0.0;
1298         kernel->maximum = kernel->values[0];
1299         kernel->negative_range = 0.0;
1300
1301         ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1302         RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1303         break;
1304       }
1305     case BinomialKernel:
1306       {
1307         size_t
1308           order_f;
1309
1310         if (args->rho < 1.0)
1311           kernel->width = kernel->height = 3;  /* default radius = 1 */
1312         else
1313           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1314         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1315
1316         order_f = fact(kernel->width-1);
1317
1318         kernel->values=(MagickRealType *) MagickAssumeAligned(
1319           AcquireAlignedMemory(kernel->width,kernel->height*
1320           sizeof(*kernel->values)));
1321         if (kernel->values == (MagickRealType *) NULL)
1322           return(DestroyKernelInfo(kernel));
1323
1324         /* set all kernel values within diamond area to scale given */
1325         for ( i=0, v=0; v < (ssize_t)kernel->height; v++)
1326           { size_t
1327               alpha = order_f / ( fact((size_t) v) * fact(kernel->height-v-1) );
1328             for ( u=0; u < (ssize_t)kernel->width; u++, i++)
1329               kernel->positive_range += kernel->values[i] = (double)
1330                 (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-u-1) ));
1331           }
1332         kernel->minimum = 1.0;
1333         kernel->maximum = kernel->values[kernel->x+kernel->y*kernel->width];
1334         kernel->negative_range = 0.0;
1335         break;
1336       }
1337
1338     /*
1339       Convolution Kernels - Well Known Named Constant Kernels
1340     */
1341     case LaplacianKernel:
1342       { switch ( (int) args->rho ) {
1343           case 0:
1344           default: /* laplacian square filter -- default */
1345             kernel=ParseKernelArray("3: -1,-1,-1  -1,8,-1  -1,-1,-1");
1346             break;
1347           case 1:  /* laplacian diamond filter */
1348             kernel=ParseKernelArray("3: 0,-1,0  -1,4,-1  0,-1,0");
1349             break;
1350           case 2:
1351             kernel=ParseKernelArray("3: -2,1,-2  1,4,1  -2,1,-2");
1352             break;
1353           case 3:
1354             kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  1,-2,1");
1355             break;
1356           case 5:   /* a 5x5 laplacian */
1357             kernel=ParseKernelArray(
1358               "5: -4,-1,0,-1,-4  -1,2,3,2,-1  0,3,4,3,0  -1,2,3,2,-1  -4,-1,0,-1,-4");
1359             break;
1360           case 7:   /* a 7x7 laplacian */
1361             kernel=ParseKernelArray(
1362               "7:-10,-5,-2,-1,-2,-5,-10 -5,0,3,4,3,0,-5 -2,3,6,7,6,3,-2 -1,4,7,8,7,4,-1 -2,3,6,7,6,3,-2 -5,0,3,4,3,0,-5 -10,-5,-2,-1,-2,-5,-10" );
1363             break;
1364           case 15:  /* a 5x5 LoG (sigma approx 1.4) */
1365             kernel=ParseKernelArray(
1366               "5: 0,0,-1,0,0  0,-1,-2,-1,0  -1,-2,16,-2,-1  0,-1,-2,-1,0  0,0,-1,0,0");
1367             break;
1368           case 19:  /* a 9x9 LoG (sigma approx 1.4) */
1369             /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1370             kernel=ParseKernelArray(
1371               "9: 0,-1,-1,-2,-2,-2,-1,-1,0  -1,-2,-4,-5,-5,-5,-4,-2,-1  -1,-4,-5,-3,-0,-3,-5,-4,-1  -2,-5,-3,12,24,12,-3,-5,-2  -2,-5,-0,24,40,24,-0,-5,-2  -2,-5,-3,12,24,12,-3,-5,-2  -1,-4,-5,-3,-0,-3,-5,-4,-1  -1,-2,-4,-5,-5,-5,-4,-2,-1  0,-1,-1,-2,-2,-2,-1,-1,0");
1372             break;
1373         }
1374         if (kernel == (KernelInfo *) NULL)
1375           return(kernel);
1376         kernel->type = type;
1377         break;
1378       }
1379     case SobelKernel:
1380       { /* Simple Sobel Kernel */
1381         kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1382         if (kernel == (KernelInfo *) NULL)
1383           return(kernel);
1384         kernel->type = type;
1385         RotateKernelInfo(kernel, args->rho);
1386         break;
1387       }
1388     case RobertsKernel:
1389       {
1390         kernel=ParseKernelArray("3: 0,0,0  1,-1,0  0,0,0");
1391         if (kernel == (KernelInfo *) NULL)
1392           return(kernel);
1393         kernel->type = type;
1394         RotateKernelInfo(kernel, args->rho);
1395         break;
1396       }
1397     case PrewittKernel:
1398       {
1399         kernel=ParseKernelArray("3: 1,0,-1  1,0,-1  1,0,-1");
1400         if (kernel == (KernelInfo *) NULL)
1401           return(kernel);
1402         kernel->type = type;
1403         RotateKernelInfo(kernel, args->rho);
1404         break;
1405       }
1406     case CompassKernel:
1407       {
1408         kernel=ParseKernelArray("3: 1,1,-1  1,-2,-1  1,1,-1");
1409         if (kernel == (KernelInfo *) NULL)
1410           return(kernel);
1411         kernel->type = type;
1412         RotateKernelInfo(kernel, args->rho);
1413         break;
1414       }
1415     case KirschKernel:
1416       {
1417         kernel=ParseKernelArray("3: 5,-3,-3  5,0,-3  5,-3,-3");
1418         if (kernel == (KernelInfo *) NULL)
1419           return(kernel);
1420         kernel->type = type;
1421         RotateKernelInfo(kernel, args->rho);
1422         break;
1423       }
1424     case FreiChenKernel:
1425       /* Direction is set to be left to right positive */
1426       /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1427       /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1428       { switch ( (int) args->rho ) {
1429           default:
1430           case 0:
1431             kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1432             if (kernel == (KernelInfo *) NULL)
1433               return(kernel);
1434             kernel->type = type;
1435             kernel->values[3] = +(MagickRealType) MagickSQ2;
1436             kernel->values[5] = -(MagickRealType) MagickSQ2;
1437             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1438             break;
1439           case 2:
1440             kernel=ParseKernelArray("3: 1,2,0  2,0,-2  0,-2,-1");
1441             if (kernel == (KernelInfo *) NULL)
1442               return(kernel);
1443             kernel->type = type;
1444             kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2;
1445             kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2;
1446             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1447             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1448             break;
1449           case 10:
1450           {
1451             kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19",exception);
1452             if (kernel == (KernelInfo *) NULL)
1453               return(kernel);
1454             break;
1455           }
1456           case 1:
1457           case 11:
1458             kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1459             if (kernel == (KernelInfo *) NULL)
1460               return(kernel);
1461             kernel->type = type;
1462             kernel->values[3] = +(MagickRealType) MagickSQ2;
1463             kernel->values[5] = -(MagickRealType) MagickSQ2;
1464             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1465             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1466             break;
1467           case 12:
1468             kernel=ParseKernelArray("3: 1,2,1  0,0,0  1,2,1");
1469             if (kernel == (KernelInfo *) NULL)
1470               return(kernel);
1471             kernel->type = type;
1472             kernel->values[1] = +(MagickRealType) MagickSQ2;
1473             kernel->values[7] = +(MagickRealType) MagickSQ2;
1474             CalcKernelMetaData(kernel);
1475             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1476             break;
1477           case 13:
1478             kernel=ParseKernelArray("3: 2,-1,0  -1,0,1  0,1,-2");
1479             if (kernel == (KernelInfo *) NULL)
1480               return(kernel);
1481             kernel->type = type;
1482             kernel->values[0] = +(MagickRealType) MagickSQ2;
1483             kernel->values[8] = -(MagickRealType) MagickSQ2;
1484             CalcKernelMetaData(kernel);
1485             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1486             break;
1487           case 14:
1488             kernel=ParseKernelArray("3: 0,1,-2  -1,0,1  2,-1,0");
1489             if (kernel == (KernelInfo *) NULL)
1490               return(kernel);
1491             kernel->type = type;
1492             kernel->values[2] = -(MagickRealType) MagickSQ2;
1493             kernel->values[6] = +(MagickRealType) MagickSQ2;
1494             CalcKernelMetaData(kernel);
1495             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1496             break;
1497           case 15:
1498             kernel=ParseKernelArray("3: 0,-1,0  1,0,1  0,-1,0");
1499             if (kernel == (KernelInfo *) NULL)
1500               return(kernel);
1501             kernel->type = type;
1502             ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1503             break;
1504           case 16:
1505             kernel=ParseKernelArray("3: 1,0,-1  0,0,0  -1,0,1");
1506             if (kernel == (KernelInfo *) NULL)
1507               return(kernel);
1508             kernel->type = type;
1509             ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1510             break;
1511           case 17:
1512             kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  -1,-2,1");
1513             if (kernel == (KernelInfo *) NULL)
1514               return(kernel);
1515             kernel->type = type;
1516             ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1517             break;
1518           case 18:
1519             kernel=ParseKernelArray("3: -2,1,-2  1,4,1  -2,1,-2");
1520             if (kernel == (KernelInfo *) NULL)
1521               return(kernel);
1522             kernel->type = type;
1523             ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1524             break;
1525           case 19:
1526             kernel=ParseKernelArray("3: 1,1,1  1,1,1  1,1,1");
1527             if (kernel == (KernelInfo *) NULL)
1528               return(kernel);
1529             kernel->type = type;
1530             ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1531             break;
1532         }
1533         if ( fabs(args->sigma) >= MagickEpsilon )
1534           /* Rotate by correctly supplied 'angle' */
1535           RotateKernelInfo(kernel, args->sigma);
1536         else if ( args->rho > 30.0 || args->rho < -30.0 )
1537           /* Rotate by out of bounds 'type' */
1538           RotateKernelInfo(kernel, args->rho);
1539         break;
1540       }
1541
1542     /*
1543       Boolean or Shaped Kernels
1544     */
1545     case DiamondKernel:
1546       {
1547         if (args->rho < 1.0)
1548           kernel->width = kernel->height = 3;  /* default radius = 1 */
1549         else
1550           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1551         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1552
1553         kernel->values=(MagickRealType *) MagickAssumeAligned(
1554           AcquireAlignedMemory(kernel->width,kernel->height*
1555           sizeof(*kernel->values)));
1556         if (kernel->values == (MagickRealType *) NULL)
1557           return(DestroyKernelInfo(kernel));
1558
1559         /* set all kernel values within diamond area to scale given */
1560         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1561           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1562             if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1563               kernel->positive_range += kernel->values[i] = args->sigma;
1564             else
1565               kernel->values[i] = nan;
1566         kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1567         break;
1568       }
1569     case SquareKernel:
1570     case RectangleKernel:
1571       { double
1572           scale;
1573         if ( type == SquareKernel )
1574           {
1575             if (args->rho < 1.0)
1576               kernel->width = kernel->height = 3;  /* default radius = 1 */
1577             else
1578               kernel->width = kernel->height = (size_t) (2*args->rho+1);
1579             kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1580             scale = args->sigma;
1581           }
1582         else {
1583             /* NOTE: user defaults set in "AcquireKernelInfo()" */
1584             if ( args->rho < 1.0 || args->sigma < 1.0 )
1585               return(DestroyKernelInfo(kernel));    /* invalid args given */
1586             kernel->width = (size_t)args->rho;
1587             kernel->height = (size_t)args->sigma;
1588             if ( args->xi  < 0.0 || args->xi  > (double)kernel->width ||
1589                  args->psi < 0.0 || args->psi > (double)kernel->height )
1590               return(DestroyKernelInfo(kernel));    /* invalid args given */
1591             kernel->x = (ssize_t) args->xi;
1592             kernel->y = (ssize_t) args->psi;
1593             scale = 1.0;
1594           }
1595         kernel->values=(MagickRealType *) MagickAssumeAligned(
1596           AcquireAlignedMemory(kernel->width,kernel->height*
1597           sizeof(*kernel->values)));
1598         if (kernel->values == (MagickRealType *) NULL)
1599           return(DestroyKernelInfo(kernel));
1600
1601         /* set all kernel values to scale given */
1602         u=(ssize_t) (kernel->width*kernel->height);
1603         for ( i=0; i < u; i++)
1604             kernel->values[i] = scale;
1605         kernel->minimum = kernel->maximum = scale;   /* a flat shape */
1606         kernel->positive_range = scale*u;
1607         break;
1608       }
1609       case OctagonKernel:
1610         {
1611           if (args->rho < 1.0)
1612             kernel->width = kernel->height = 5;  /* default radius = 2 */
1613           else
1614             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1615           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1616
1617           kernel->values=(MagickRealType *) MagickAssumeAligned(
1618             AcquireAlignedMemory(kernel->width,kernel->height*
1619             sizeof(*kernel->values)));
1620           if (kernel->values == (MagickRealType *) NULL)
1621             return(DestroyKernelInfo(kernel));
1622
1623           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1624             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1625               if ( (labs((long) u)+labs((long) v)) <=
1626                         ((long)kernel->x + (long)(kernel->x/2)) )
1627                 kernel->positive_range += kernel->values[i] = args->sigma;
1628               else
1629                 kernel->values[i] = nan;
1630           kernel->minimum = kernel->maximum = args->sigma;  /* a flat shape */
1631           break;
1632         }
1633       case DiskKernel:
1634         {
1635           ssize_t
1636             limit = (ssize_t)(args->rho*args->rho);
1637
1638           if (args->rho < 0.4)           /* default radius approx 4.3 */
1639             kernel->width = kernel->height = 9L, limit = 18L;
1640           else
1641             kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1642           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1643
1644           kernel->values=(MagickRealType *) MagickAssumeAligned(
1645             AcquireAlignedMemory(kernel->width,kernel->height*
1646             sizeof(*kernel->values)));
1647           if (kernel->values == (MagickRealType *) NULL)
1648             return(DestroyKernelInfo(kernel));
1649
1650           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1651             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1652               if ((u*u+v*v) <= limit)
1653                 kernel->positive_range += kernel->values[i] = args->sigma;
1654               else
1655                 kernel->values[i] = nan;
1656           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1657           break;
1658         }
1659       case PlusKernel:
1660         {
1661           if (args->rho < 1.0)
1662             kernel->width = kernel->height = 5;  /* default radius 2 */
1663           else
1664             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1665           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1666
1667           kernel->values=(MagickRealType *) MagickAssumeAligned(
1668             AcquireAlignedMemory(kernel->width,kernel->height*
1669             sizeof(*kernel->values)));
1670           if (kernel->values == (MagickRealType *) NULL)
1671             return(DestroyKernelInfo(kernel));
1672
1673           /* set all kernel values along axises to given scale */
1674           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1675             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1676               kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1677           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1678           kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1679           break;
1680         }
1681       case CrossKernel:
1682         {
1683           if (args->rho < 1.0)
1684             kernel->width = kernel->height = 5;  /* default radius 2 */
1685           else
1686             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1687           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1688
1689           kernel->values=(MagickRealType *) MagickAssumeAligned(
1690             AcquireAlignedMemory(kernel->width,kernel->height*
1691             sizeof(*kernel->values)));
1692           if (kernel->values == (MagickRealType *) NULL)
1693             return(DestroyKernelInfo(kernel));
1694
1695           /* set all kernel values along axises to given scale */
1696           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1697             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1698               kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1699           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1700           kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1701           break;
1702         }
1703       /*
1704         HitAndMiss Kernels
1705       */
1706       case RingKernel:
1707       case PeaksKernel:
1708         {
1709           ssize_t
1710             limit1,
1711             limit2,
1712             scale;
1713
1714           if (args->rho < args->sigma)
1715             {
1716               kernel->width = ((size_t)args->sigma)*2+1;
1717               limit1 = (ssize_t)(args->rho*args->rho);
1718               limit2 = (ssize_t)(args->sigma*args->sigma);
1719             }
1720           else
1721             {
1722               kernel->width = ((size_t)args->rho)*2+1;
1723               limit1 = (ssize_t)(args->sigma*args->sigma);
1724               limit2 = (ssize_t)(args->rho*args->rho);
1725             }
1726           if ( limit2 <= 0 )
1727             kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1728
1729           kernel->height = kernel->width;
1730           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1731           kernel->values=(MagickRealType *) MagickAssumeAligned(
1732             AcquireAlignedMemory(kernel->width,kernel->height*
1733             sizeof(*kernel->values)));
1734           if (kernel->values == (MagickRealType *) NULL)
1735             return(DestroyKernelInfo(kernel));
1736
1737           /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1738           scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1739           for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1740             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1741               { ssize_t radius=u*u+v*v;
1742                 if (limit1 < radius && radius <= limit2)
1743                   kernel->positive_range += kernel->values[i] = (double) scale;
1744                 else
1745                   kernel->values[i] = nan;
1746               }
1747           kernel->minimum = kernel->maximum = (double) scale;
1748           if ( type == PeaksKernel ) {
1749             /* set the central point in the middle */
1750             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1751             kernel->positive_range = 1.0;
1752             kernel->maximum = 1.0;
1753           }
1754           break;
1755         }
1756       case EdgesKernel:
1757         {
1758           kernel=AcquireKernelInfo("ThinSE:482",exception);
1759           if (kernel == (KernelInfo *) NULL)
1760             return(kernel);
1761           kernel->type = type;
1762           ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1763           break;
1764         }
1765       case CornersKernel:
1766         {
1767           kernel=AcquireKernelInfo("ThinSE:87",exception);
1768           if (kernel == (KernelInfo *) NULL)
1769             return(kernel);
1770           kernel->type = type;
1771           ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1772           break;
1773         }
1774       case DiagonalsKernel:
1775         {
1776           switch ( (int) args->rho ) {
1777             case 0:
1778             default:
1779               { KernelInfo
1780                   *new_kernel;
1781                 kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
1782                 if (kernel == (KernelInfo *) NULL)
1783                   return(kernel);
1784                 kernel->type = type;
1785                 new_kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
1786                 if (new_kernel == (KernelInfo *) NULL)
1787                   return(DestroyKernelInfo(kernel));
1788                 new_kernel->type = type;
1789                 LastKernelInfo(kernel)->next = new_kernel;
1790                 ExpandMirrorKernelInfo(kernel);
1791                 return(kernel);
1792               }
1793             case 1:
1794               kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
1795               break;
1796             case 2:
1797               kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
1798               break;
1799           }
1800           if (kernel == (KernelInfo *) NULL)
1801             return(kernel);
1802           kernel->type = type;
1803           RotateKernelInfo(kernel, args->sigma);
1804           break;
1805         }
1806       case LineEndsKernel:
1807         { /* Kernels for finding the end of thin lines */
1808           switch ( (int) args->rho ) {
1809             case 0:
1810             default:
1811               /* set of kernels to find all end of lines */
1812               return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>",exception));
1813             case 1:
1814               /* kernel for 4-connected line ends - no rotation */
1815               kernel=ParseKernelArray("3: 0,0,-  0,1,1  0,0,-");
1816               break;
1817           case 2:
1818               /* kernel to add for 8-connected lines - no rotation */
1819               kernel=ParseKernelArray("3: 0,0,0  0,1,0  0,0,1");
1820               break;
1821           case 3:
1822               /* kernel to add for orthogonal line ends - does not find corners */
1823               kernel=ParseKernelArray("3: 0,0,0  0,1,1  0,0,0");
1824               break;
1825           case 4:
1826               /* traditional line end - fails on last T end */
1827               kernel=ParseKernelArray("3: 0,0,0  0,1,-  0,0,-");
1828               break;
1829           }
1830           if (kernel == (KernelInfo *) NULL)
1831             return(kernel);
1832           kernel->type = type;
1833           RotateKernelInfo(kernel, args->sigma);
1834           break;
1835         }
1836       case LineJunctionsKernel:
1837         { /* kernels for finding the junctions of multiple lines */
1838           switch ( (int) args->rho ) {
1839             case 0:
1840             default:
1841               /* set of kernels to find all line junctions */
1842               return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>",exception));
1843             case 1:
1844               /* Y Junction */
1845               kernel=ParseKernelArray("3: 1,-,1  -,1,-  -,1,-");
1846               break;
1847             case 2:
1848               /* Diagonal T Junctions */
1849               kernel=ParseKernelArray("3: 1,-,-  -,1,-  1,-,1");
1850               break;
1851             case 3:
1852               /* Orthogonal T Junctions */
1853               kernel=ParseKernelArray("3: -,-,-  1,1,1  -,1,-");
1854               break;
1855             case 4:
1856               /* Diagonal X Junctions */
1857               kernel=ParseKernelArray("3: 1,-,1  -,1,-  1,-,1");
1858               break;
1859             case 5:
1860               /* Orthogonal X Junctions - minimal diamond kernel */
1861               kernel=ParseKernelArray("3: -,1,-  1,1,1  -,1,-");
1862               break;
1863           }
1864           if (kernel == (KernelInfo *) NULL)
1865             return(kernel);
1866           kernel->type = type;
1867           RotateKernelInfo(kernel, args->sigma);
1868           break;
1869         }
1870       case RidgesKernel:
1871         { /* Ridges - Ridge finding kernels */
1872           KernelInfo
1873             *new_kernel;
1874           switch ( (int) args->rho ) {
1875             case 1:
1876             default:
1877               kernel=ParseKernelArray("3x1:0,1,0");
1878               if (kernel == (KernelInfo *) NULL)
1879                 return(kernel);
1880               kernel->type = type;
1881               ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1882               break;
1883             case 2:
1884               kernel=ParseKernelArray("4x1:0,1,1,0");
1885               if (kernel == (KernelInfo *) NULL)
1886                 return(kernel);
1887               kernel->type = type;
1888               ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1889
1890               /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1891               /* Unfortunatally we can not yet rotate a non-square kernel */
1892               /* But then we can't flip a non-symetrical kernel either */
1893               new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1894               if (new_kernel == (KernelInfo *) NULL)
1895                 return(DestroyKernelInfo(kernel));
1896               new_kernel->type = type;
1897               LastKernelInfo(kernel)->next = new_kernel;
1898               new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1899               if (new_kernel == (KernelInfo *) NULL)
1900                 return(DestroyKernelInfo(kernel));
1901               new_kernel->type = type;
1902               LastKernelInfo(kernel)->next = new_kernel;
1903               new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1904               if (new_kernel == (KernelInfo *) NULL)
1905                 return(DestroyKernelInfo(kernel));
1906               new_kernel->type = type;
1907               LastKernelInfo(kernel)->next = new_kernel;
1908               new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1909               if (new_kernel == (KernelInfo *) NULL)
1910                 return(DestroyKernelInfo(kernel));
1911               new_kernel->type = type;
1912               LastKernelInfo(kernel)->next = new_kernel;
1913               new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1914               if (new_kernel == (KernelInfo *) NULL)
1915                 return(DestroyKernelInfo(kernel));
1916               new_kernel->type = type;
1917               LastKernelInfo(kernel)->next = new_kernel;
1918               new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1919               if (new_kernel == (KernelInfo *) NULL)
1920                 return(DestroyKernelInfo(kernel));
1921               new_kernel->type = type;
1922               LastKernelInfo(kernel)->next = new_kernel;
1923               new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1924               if (new_kernel == (KernelInfo *) NULL)
1925                 return(DestroyKernelInfo(kernel));
1926               new_kernel->type = type;
1927               LastKernelInfo(kernel)->next = new_kernel;
1928               new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1929               if (new_kernel == (KernelInfo *) NULL)
1930                 return(DestroyKernelInfo(kernel));
1931               new_kernel->type = type;
1932               LastKernelInfo(kernel)->next = new_kernel;
1933               break;
1934           }
1935           break;
1936         }
1937       case ConvexHullKernel:
1938         {
1939           KernelInfo
1940             *new_kernel;
1941           /* first set of 8 kernels */
1942           kernel=ParseKernelArray("3: 1,1,-  1,0,-  1,-,0");
1943           if (kernel == (KernelInfo *) NULL)
1944             return(kernel);
1945           kernel->type = type;
1946           ExpandRotateKernelInfo(kernel, 90.0);
1947           /* append the mirror versions too - no flip function yet */
1948           new_kernel=ParseKernelArray("3: 1,1,1  1,0,-  -,-,0");
1949           if (new_kernel == (KernelInfo *) NULL)
1950             return(DestroyKernelInfo(kernel));
1951           new_kernel->type = type;
1952           ExpandRotateKernelInfo(new_kernel, 90.0);
1953           LastKernelInfo(kernel)->next = new_kernel;
1954           break;
1955         }
1956       case SkeletonKernel:
1957         {
1958           switch ( (int) args->rho ) {
1959             case 1:
1960             default:
1961               /* Traditional Skeleton...
1962               ** A cyclically rotated single kernel
1963               */
1964               kernel=AcquireKernelInfo("ThinSE:482",exception);
1965               if (kernel == (KernelInfo *) NULL)
1966                 return(kernel);
1967               kernel->type = type;
1968               ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1969               break;
1970             case 2:
1971               /* HIPR Variation of the cyclic skeleton
1972               ** Corners of the traditional method made more forgiving,
1973               ** but the retain the same cyclic order.
1974               */
1975               kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;",exception);
1976               if (kernel == (KernelInfo *) NULL)
1977                 return(kernel);
1978               if (kernel->next == (KernelInfo *) NULL)
1979                 return(DestroyKernelInfo(kernel));
1980               kernel->type = type;
1981               kernel->next->type = type;
1982               ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1983               break;
1984             case 3:
1985               /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1986               ** "Connectivity-Preserving Morphological Image Thransformations"
1987               ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1988               **   http://www.leptonica.com/papers/conn.pdf
1989               */
1990               kernel=AcquireKernelInfo("ThinSE:41; ThinSE:42; ThinSE:43",
1991                 exception);
1992               if (kernel == (KernelInfo *) NULL)
1993                 return(kernel);
1994               kernel->type = type;
1995               kernel->next->type = type;
1996               kernel->next->next->type = type;
1997               ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1998               break;
1999            }
2000           break;
2001         }
2002       case ThinSEKernel:
2003         { /* Special kernels for general thinning, while preserving connections
2004           ** "Connectivity-Preserving Morphological Image Thransformations"
2005           ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
2006           **   http://www.leptonica.com/papers/conn.pdf
2007           ** And
2008           **   http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2009           **
2010           ** Note kernels do not specify the origin pixel, allowing them
2011           ** to be used for both thickening and thinning operations.
2012           */
2013           switch ( (int) args->rho ) {
2014             /* SE for 4-connected thinning */
2015             case 41: /* SE_4_1 */
2016               kernel=ParseKernelArray("3: -,-,1  0,-,1  -,-,1");
2017               break;
2018             case 42: /* SE_4_2 */
2019               kernel=ParseKernelArray("3: -,-,1  0,-,1  -,0,-");
2020               break;
2021             case 43: /* SE_4_3 */
2022               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,-,1");
2023               break;
2024             case 44: /* SE_4_4 */
2025               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,-");
2026               break;
2027             case 45: /* SE_4_5 */
2028               kernel=ParseKernelArray("3: -,0,1  0,-,1  -,0,-");
2029               break;
2030             case 46: /* SE_4_6 */
2031               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,1");
2032               break;
2033             case 47: /* SE_4_7 */
2034               kernel=ParseKernelArray("3: -,1,1  0,-,1  -,0,-");
2035               break;
2036             case 48: /* SE_4_8 */
2037               kernel=ParseKernelArray("3: -,-,1  0,-,1  0,-,1");
2038               break;
2039             case 49: /* SE_4_9 */
2040               kernel=ParseKernelArray("3: 0,-,1  0,-,1  -,-,1");
2041               break;
2042             /* SE for 8-connected thinning - negatives of the above */
2043             case 81: /* SE_8_0 */
2044               kernel=ParseKernelArray("3: -,1,-  0,-,1  -,1,-");
2045               break;
2046             case 82: /* SE_8_2 */
2047               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,-,-");
2048               break;
2049             case 83: /* SE_8_3 */
2050               kernel=ParseKernelArray("3: 0,-,-  0,-,1  -,1,-");
2051               break;
2052             case 84: /* SE_8_4 */
2053               kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,-");
2054               break;
2055             case 85: /* SE_8_5 */
2056               kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,-");
2057               break;
2058             case 86: /* SE_8_6 */
2059               kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,1");
2060               break;
2061             case 87: /* SE_8_7 */
2062               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,0,-");
2063               break;
2064             case 88: /* SE_8_8 */
2065               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,1,-");
2066               break;
2067             case 89: /* SE_8_9 */
2068               kernel=ParseKernelArray("3: 0,1,-  0,-,1  -,1,-");
2069               break;
2070             /* Special combined SE kernels */
2071             case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2072               kernel=ParseKernelArray("3: -,-,1  0,-,-  -,0,-");
2073               break;
2074             case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2075               kernel=ParseKernelArray("3: -,1,-  -,-,1  0,-,-");
2076               break;
2077             case 481: /* SE_48_1 - General Connected Corner Kernel */
2078               kernel=ParseKernelArray("3: -,1,1  0,-,1  0,0,-");
2079               break;
2080             default:
2081             case 482: /* SE_48_2 - General Edge Kernel */
2082               kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,1");
2083               break;
2084           }
2085           if (kernel == (KernelInfo *) NULL)
2086             return(kernel);
2087           kernel->type = type;
2088           RotateKernelInfo(kernel, args->sigma);
2089           break;
2090         }
2091       /*
2092         Distance Measuring Kernels
2093       */
2094       case ChebyshevKernel:
2095         {
2096           if (args->rho < 1.0)
2097             kernel->width = kernel->height = 3;  /* default radius = 1 */
2098           else
2099             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2100           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2101
2102           kernel->values=(MagickRealType *) MagickAssumeAligned(
2103             AcquireAlignedMemory(kernel->width,kernel->height*
2104             sizeof(*kernel->values)));
2105           if (kernel->values == (MagickRealType *) NULL)
2106             return(DestroyKernelInfo(kernel));
2107
2108           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2109             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2110               kernel->positive_range += ( kernel->values[i] =
2111                   args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2112           kernel->maximum = kernel->values[0];
2113           break;
2114         }
2115       case ManhattanKernel:
2116         {
2117           if (args->rho < 1.0)
2118             kernel->width = kernel->height = 3;  /* default radius = 1 */
2119           else
2120             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2121           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2122
2123           kernel->values=(MagickRealType *) MagickAssumeAligned(
2124             AcquireAlignedMemory(kernel->width,kernel->height*
2125             sizeof(*kernel->values)));
2126           if (kernel->values == (MagickRealType *) NULL)
2127             return(DestroyKernelInfo(kernel));
2128
2129           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2130             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2131               kernel->positive_range += ( kernel->values[i] =
2132                   args->sigma*(labs((long) u)+labs((long) v)) );
2133           kernel->maximum = kernel->values[0];
2134           break;
2135         }
2136       case OctagonalKernel:
2137       {
2138         if (args->rho < 2.0)
2139           kernel->width = kernel->height = 5;  /* default/minimum radius = 2 */
2140         else
2141           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2142         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2143
2144         kernel->values=(MagickRealType *) MagickAssumeAligned(
2145           AcquireAlignedMemory(kernel->width,kernel->height*
2146           sizeof(*kernel->values)));
2147         if (kernel->values == (MagickRealType *) NULL)
2148           return(DestroyKernelInfo(kernel));
2149
2150         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2151           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2152             {
2153               double
2154                 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2155                 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2156               kernel->positive_range += kernel->values[i] =
2157                         args->sigma*MagickMax(r1,r2);
2158             }
2159         kernel->maximum = kernel->values[0];
2160         break;
2161       }
2162     case EuclideanKernel:
2163       {
2164         if (args->rho < 1.0)
2165           kernel->width = kernel->height = 3;  /* default radius = 1 */
2166         else
2167           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2168         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2169
2170         kernel->values=(MagickRealType *) MagickAssumeAligned(
2171           AcquireAlignedMemory(kernel->width,kernel->height*
2172           sizeof(*kernel->values)));
2173         if (kernel->values == (MagickRealType *) NULL)
2174           return(DestroyKernelInfo(kernel));
2175
2176         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2177           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2178             kernel->positive_range += ( kernel->values[i] =
2179               args->sigma*sqrt((double)(u*u+v*v)) );
2180         kernel->maximum = kernel->values[0];
2181         break;
2182       }
2183     default:
2184       {
2185         /* No-Op Kernel - Basically just a single pixel on its own */
2186         kernel=ParseKernelArray("1:1");
2187         if (kernel == (KernelInfo *) NULL)
2188           return(kernel);
2189         kernel->type = UndefinedKernel;
2190         break;
2191       }
2192       break;
2193   }
2194   return(kernel);
2195 }
2196 \f
2197 /*
2198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2199 %                                                                             %
2200 %                                                                             %
2201 %                                                                             %
2202 %     C l o n e K e r n e l I n f o                                           %
2203 %                                                                             %
2204 %                                                                             %
2205 %                                                                             %
2206 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2207 %
2208 %  CloneKernelInfo() creates a new clone of the given Kernel List so that its
2209 %  can be modified without effecting the original.  The cloned kernel should
2210 %  be destroyed using DestoryKernelInfo() when no longer needed.
2211 %
2212 %  The format of the CloneKernelInfo method is:
2213 %
2214 %      KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2215 %
2216 %  A description of each parameter follows:
2217 %
2218 %    o kernel: the Morphology/Convolution kernel to be cloned
2219 %
2220 */
2221 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2222 {
2223   register ssize_t
2224     i;
2225
2226   KernelInfo
2227     *new_kernel;
2228
2229   assert(kernel != (KernelInfo *) NULL);
2230   new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2231   if (new_kernel == (KernelInfo *) NULL)
2232     return(new_kernel);
2233   *new_kernel=(*kernel); /* copy values in structure */
2234
2235   /* replace the values with a copy of the values */
2236   new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2237     AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2238   if (new_kernel->values == (MagickRealType *) NULL)
2239     return(DestroyKernelInfo(new_kernel));
2240   for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2241     new_kernel->values[i]=kernel->values[i];
2242
2243   /* Also clone the next kernel in the kernel list */
2244   if ( kernel->next != (KernelInfo *) NULL ) {
2245     new_kernel->next = CloneKernelInfo(kernel->next);
2246     if ( new_kernel->next == (KernelInfo *) NULL )
2247       return(DestroyKernelInfo(new_kernel));
2248   }
2249
2250   return(new_kernel);
2251 }
2252 \f
2253 /*
2254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2255 %                                                                             %
2256 %                                                                             %
2257 %                                                                             %
2258 %     D e s t r o y K e r n e l I n f o                                       %
2259 %                                                                             %
2260 %                                                                             %
2261 %                                                                             %
2262 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2263 %
2264 %  DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2265 %  kernel.
2266 %
2267 %  The format of the DestroyKernelInfo method is:
2268 %
2269 %      KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2270 %
2271 %  A description of each parameter follows:
2272 %
2273 %    o kernel: the Morphology/Convolution kernel to be destroyed
2274 %
2275 */
2276 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2277 {
2278   assert(kernel != (KernelInfo *) NULL);
2279   if (kernel->next != (KernelInfo *) NULL)
2280     kernel->next=DestroyKernelInfo(kernel->next);
2281   kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2282   kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2283   return(kernel);
2284 }
2285 \f
2286 /*
2287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2288 %                                                                             %
2289 %                                                                             %
2290 %                                                                             %
2291 +     E x p a n d M i r r o r K e r n e l I n f o                             %
2292 %                                                                             %
2293 %                                                                             %
2294 %                                                                             %
2295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2296 %
2297 %  ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2298 %  sequence of 90-degree rotated kernels but providing a reflected 180
2299 %  rotatation, before the -/+ 90-degree rotations.
2300 %
2301 %  This special rotation order produces a better, more symetrical thinning of
2302 %  objects.
2303 %
2304 %  The format of the ExpandMirrorKernelInfo method is:
2305 %
2306 %      void ExpandMirrorKernelInfo(KernelInfo *kernel)
2307 %
2308 %  A description of each parameter follows:
2309 %
2310 %    o kernel: the Morphology/Convolution kernel
2311 %
2312 % This function is only internel to this module, as it is not finalized,
2313 % especially with regard to non-orthogonal angles, and rotation of larger
2314 % 2D kernels.
2315 */
2316
2317 #if 0
2318 static void FlopKernelInfo(KernelInfo *kernel)
2319     { /* Do a Flop by reversing each row. */
2320       size_t
2321         y;
2322       register ssize_t
2323         x,r;
2324       register double
2325         *k,t;
2326
2327       for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2328         for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2329           t=k[x],  k[x]=k[r],  k[r]=t;
2330
2331       kernel->x = kernel->width - kernel->x - 1;
2332       angle = fmod(angle+180.0, 360.0);
2333     }
2334 #endif
2335
2336 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2337 {
2338   KernelInfo
2339     *clone,
2340     *last;
2341
2342   last = kernel;
2343
2344   clone = CloneKernelInfo(last);
2345   RotateKernelInfo(clone, 180);   /* flip */
2346   LastKernelInfo(last)->next = clone;
2347   last = clone;
2348
2349   clone = CloneKernelInfo(last);
2350   RotateKernelInfo(clone, 90);   /* transpose */
2351   LastKernelInfo(last)->next = clone;
2352   last = clone;
2353
2354   clone = CloneKernelInfo(last);
2355   RotateKernelInfo(clone, 180);  /* flop */
2356   LastKernelInfo(last)->next = clone;
2357
2358   return;
2359 }
2360 \f
2361 /*
2362 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2363 %                                                                             %
2364 %                                                                             %
2365 %                                                                             %
2366 +     E x p a n d R o t a t e K e r n e l I n f o                             %
2367 %                                                                             %
2368 %                                                                             %
2369 %                                                                             %
2370 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2371 %
2372 %  ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2373 %  incrementally by the angle given, until the kernel repeats.
2374 %
2375 %  WARNING: 45 degree rotations only works for 3x3 kernels.
2376 %  While 90 degree roatations only works for linear and square kernels
2377 %
2378 %  The format of the ExpandRotateKernelInfo method is:
2379 %
2380 %      void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2381 %
2382 %  A description of each parameter follows:
2383 %
2384 %    o kernel: the Morphology/Convolution kernel
2385 %
2386 %    o angle: angle to rotate in degrees
2387 %
2388 % This function is only internel to this module, as it is not finalized,
2389 % especially with regard to non-orthogonal angles, and rotation of larger
2390 % 2D kernels.
2391 */
2392
2393 /* Internal Routine - Return true if two kernels are the same */
2394 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2395      const KernelInfo *kernel2)
2396 {
2397   register size_t
2398     i;
2399
2400   /* check size and origin location */
2401   if (    kernel1->width != kernel2->width
2402        || kernel1->height != kernel2->height
2403        || kernel1->x != kernel2->x
2404        || kernel1->y != kernel2->y )
2405     return MagickFalse;
2406
2407   /* check actual kernel values */
2408   for (i=0; i < (kernel1->width*kernel1->height); i++) {
2409     /* Test for Nan equivalence */
2410     if ( IfNaN(kernel1->values[i]) && !IfNaN(kernel2->values[i]) )
2411       return MagickFalse;
2412     if ( IfNaN(kernel2->values[i]) && !IfNaN(kernel1->values[i]) )
2413       return MagickFalse;
2414     /* Test actual values are equivalent */
2415     if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2416       return MagickFalse;
2417   }
2418
2419   return MagickTrue;
2420 }
2421
2422 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2423 {
2424   KernelInfo
2425     *clone,
2426     *last;
2427
2428   last = kernel;
2429 DisableMSCWarning(4127)
2430   while(1) {
2431 RestoreMSCWarning
2432     clone = CloneKernelInfo(last);
2433     RotateKernelInfo(clone, angle);
2434     if ( SameKernelInfo(kernel, clone) != MagickFalse )
2435       break;
2436     LastKernelInfo(last)->next = clone;
2437     last = clone;
2438   }
2439   clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2440   return;
2441 }
2442 \f
2443 /*
2444 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2445 %                                                                             %
2446 %                                                                             %
2447 %                                                                             %
2448 +     C a l c M e t a K e r n a l I n f o                                     %
2449 %                                                                             %
2450 %                                                                             %
2451 %                                                                             %
2452 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2453 %
2454 %  CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2455 %  using the kernel values.  This should only ne used if it is not possible to
2456 %  calculate that meta-data in some easier way.
2457 %
2458 %  It is important that the meta-data is correct before ScaleKernelInfo() is
2459 %  used to perform kernel normalization.
2460 %
2461 %  The format of the CalcKernelMetaData method is:
2462 %
2463 %      void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2464 %
2465 %  A description of each parameter follows:
2466 %
2467 %    o kernel: the Morphology/Convolution kernel to modify
2468 %
2469 %  WARNING: Minimum and Maximum values are assumed to include zero, even if
2470 %  zero is not part of the kernel (as in Gaussian Derived kernels). This
2471 %  however is not true for flat-shaped morphological kernels.
2472 %
2473 %  WARNING: Only the specific kernel pointed to is modified, not a list of
2474 %  multiple kernels.
2475 %
2476 % This is an internal function and not expected to be useful outside this
2477 % module.  This could change however.
2478 */
2479 static void CalcKernelMetaData(KernelInfo *kernel)
2480 {
2481   register size_t
2482     i;
2483
2484   kernel->minimum = kernel->maximum = 0.0;
2485   kernel->negative_range = kernel->positive_range = 0.0;
2486   for (i=0; i < (kernel->width*kernel->height); i++)
2487     {
2488       if ( fabs(kernel->values[i]) < MagickEpsilon )
2489         kernel->values[i] = 0.0;
2490       ( kernel->values[i] < 0)
2491           ?  ( kernel->negative_range += kernel->values[i] )
2492           :  ( kernel->positive_range += kernel->values[i] );
2493       Minimize(kernel->minimum, kernel->values[i]);
2494       Maximize(kernel->maximum, kernel->values[i]);
2495     }
2496
2497   return;
2498 }
2499 \f
2500 /*
2501 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2502 %                                                                             %
2503 %                                                                             %
2504 %                                                                             %
2505 %     M o r p h o l o g y A p p l y                                           %
2506 %                                                                             %
2507 %                                                                             %
2508 %                                                                             %
2509 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2510 %
2511 %  MorphologyApply() applies a morphological method, multiple times using
2512 %  a list of multiple kernels.  This is the method that should be called by
2513 %  other 'operators' that internally use morphology operations as part of
2514 %  their processing.
2515 %
2516 %  It is basically equivalent to as MorphologyImage() (see below) but without
2517 %  any user controls.  This allows internel programs to use this method to
2518 %  perform a specific task without possible interference by any API user
2519 %  supplied settings.
2520 %
2521 %  It is MorphologyImage() task to extract any such user controls, and
2522 %  pass them to this function for processing.
2523 %
2524 %  More specifically all given kernels should already be scaled, normalised,
2525 %  and blended appropriatally before being parred to this routine. The
2526 %  appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2527 %
2528 %  The format of the MorphologyApply method is:
2529 %
2530 %      Image *MorphologyApply(const Image *image,MorphologyMethod method,
2531 %        const ssize_t iterations,const KernelInfo *kernel,
2532 %        const CompositeMethod compose,const double bias,
2533 %        ExceptionInfo *exception)
2534 %
2535 %  A description of each parameter follows:
2536 %
2537 %    o image: the source image
2538 %
2539 %    o method: the morphology method to be applied.
2540 %
2541 %    o iterations: apply the operation this many times (or no change).
2542 %                  A value of -1 means loop until no change found.
2543 %                  How this is applied may depend on the morphology method.
2544 %                  Typically this is a value of 1.
2545 %
2546 %    o channel: the channel type.
2547 %
2548 %    o kernel: An array of double representing the morphology kernel.
2549 %
2550 %    o compose: How to handle or merge multi-kernel results.
2551 %          If 'UndefinedCompositeOp' use default for the Morphology method.
2552 %          If 'NoCompositeOp' force image to be re-iterated by each kernel.
2553 %          Otherwise merge the results using the compose method given.
2554 %
2555 %    o bias: Convolution Output Bias.
2556 %
2557 %    o exception: return any errors or warnings in this structure.
2558 %
2559 */
2560 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2561   const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2562   ExceptionInfo *exception)
2563 {
2564 #define MorphologyTag  "Morphology/Image"
2565
2566   CacheView
2567     *image_view,
2568     *morphology_view;
2569
2570   OffsetInfo
2571     offset;
2572
2573   register ssize_t
2574     i;
2575
2576   ssize_t
2577     y;
2578
2579   size_t
2580     *changes,
2581     changed,
2582     width;
2583
2584   MagickBooleanType
2585     status;
2586
2587   MagickOffsetType
2588     progress;
2589
2590   assert(image != (Image *) NULL);
2591   assert(image->signature == MagickSignature);
2592   assert(morphology_image != (Image *) NULL);
2593   assert(morphology_image->signature == MagickSignature);
2594   assert(kernel != (KernelInfo *) NULL);
2595   assert(kernel->signature == MagickSignature);
2596   assert(exception != (ExceptionInfo *) NULL);
2597   assert(exception->signature == MagickSignature);
2598   status=MagickTrue;
2599   progress=0;
2600   image_view=AcquireVirtualCacheView(image,exception);
2601   morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2602   width=image->columns+kernel->width-1;
2603   offset.x=0.0;
2604   offset.y=0.0;
2605   switch (method)
2606   {
2607     case ConvolveMorphology:
2608     case DilateMorphology:
2609     case DilateIntensityMorphology:
2610     case IterativeDistanceMorphology:
2611     {
2612       /*
2613         Kernel needs to used with reflection about origin.
2614       */
2615       offset.x=(ssize_t) kernel->width-kernel->x-1;
2616       offset.y=(ssize_t) kernel->height-kernel->y-1;
2617       break;
2618     }
2619     case ErodeMorphology:
2620     case ErodeIntensityMorphology:
2621     case HitAndMissMorphology:
2622     case ThinningMorphology:
2623     case ThickenMorphology:
2624     {
2625       offset.x=kernel->x;
2626       offset.y=kernel->y;
2627       break;
2628     }
2629     default:
2630     {
2631       assert("Not a Primitive Morphology Method" != (char *) NULL);
2632       break;
2633     }
2634   }
2635   changed=0;
2636   changes=(size_t *) AcquireQuantumMemory(GetOpenMPMaximumThreads(),
2637     sizeof(*changes));
2638   if (changes == (size_t *) NULL)
2639     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2640   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2641     changes[i]=0;
2642   if ((method == ConvolveMorphology) && (kernel->width == 1))
2643     {
2644       register ssize_t
2645         x;
2646
2647       /*
2648         Special handling (for speed) of vertical (blur) kernels.  This performs
2649         its handling in columns rather than in rows.  This is only done
2650         for convolve as it is the only method that generates very large 1-D
2651         vertical kernels (such as a 'BlurKernel')
2652      */
2653 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2654      #pragma omp parallel for schedule(static,4) shared(progress,status) \
2655        magick_threads(image,morphology_image,image->columns,1)
2656 #endif
2657       for (x=0; x < (ssize_t) image->columns; x++)
2658       {
2659         const int
2660           id = GetOpenMPThreadId();
2661
2662         register const Quantum
2663           *restrict p;
2664
2665         register Quantum
2666           *restrict q;
2667
2668         register ssize_t
2669           y;
2670
2671         ssize_t
2672           center;
2673
2674         if (status == MagickFalse)
2675           continue;
2676         p=GetCacheViewVirtualPixels(image_view,x,-offset.y,1,image->rows+
2677           kernel->height-1,exception);
2678         q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2679           morphology_image->rows,exception);
2680         if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2681           {
2682             status=MagickFalse;
2683             continue;
2684           }
2685         center=(ssize_t) GetPixelChannels(image)*offset.y;
2686         for (y=0; y < (ssize_t) image->rows; y++)
2687         {
2688           register ssize_t
2689             i;
2690
2691           for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2692           {
2693             double
2694               alpha,
2695               gamma,
2696               pixel;
2697
2698             PixelChannel
2699               channel;
2700
2701             PixelTrait
2702               morphology_traits,
2703               traits;
2704
2705             register const MagickRealType
2706               *restrict k;
2707
2708             register const Quantum
2709               *restrict pixels;
2710
2711             register ssize_t
2712               u;
2713
2714             size_t
2715               count;
2716
2717             ssize_t
2718               v;
2719
2720             channel=GetPixelChannelChannel(image,i);
2721             traits=GetPixelChannelTraits(image,channel);
2722             morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2723             if ((traits == UndefinedPixelTrait) ||
2724                 (morphology_traits == UndefinedPixelTrait))
2725               continue;
2726             if (((traits & CopyPixelTrait) != 0) ||
2727                 (GetPixelReadMask(image,p+center) == 0))
2728               {
2729                 SetPixelChannel(morphology_image,channel,p[center+i],q);
2730                 continue;
2731               }
2732             k=(&kernel->values[kernel->width*kernel->height-1]);
2733             pixels=p;
2734             pixel=bias;
2735             gamma=0.0;
2736             count=0;
2737             if ((morphology_traits & BlendPixelTrait) == 0)
2738               for (v=0; v < (ssize_t) kernel->height; v++)
2739               {
2740                 for (u=0; u < (ssize_t) kernel->width; u++)
2741                 {
2742                   if (IfNaN(*k) == MagickFalse)
2743                     {
2744                       pixel+=(*k)*pixels[i];
2745                       gamma+=(*k);
2746                       count++;
2747                     }
2748                   k--;
2749                   pixels+=GetPixelChannels(image);
2750                 }
2751               }
2752             else
2753               for (v=0; v < (ssize_t) kernel->height; v++)
2754               {
2755                 for (u=0; u < (ssize_t) kernel->width; u++)
2756                 {
2757                   if (IfNaN(*k) == MagickFalse)
2758                     {
2759                       alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2760                       pixel+=alpha*(*k)*pixels[i];
2761                       gamma+=alpha*(*k);
2762                       count++;
2763                     }
2764                   k--;
2765                   pixels+=GetPixelChannels(image);
2766                 }
2767               }
2768             if (fabs(pixel-p[center+i]) > MagickEpsilon)
2769               changes[id]++;
2770             gamma=PerceptibleReciprocal(gamma);
2771             if (count != 0)
2772               gamma*=(double) kernel->height*kernel->width/count;
2773             SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2774               pixel),q);
2775           }
2776           p+=GetPixelChannels(image);
2777           q+=GetPixelChannels(morphology_image);
2778         }
2779         if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2780           status=MagickFalse;
2781         if (image->progress_monitor != (MagickProgressMonitor) NULL)
2782           {
2783             MagickBooleanType
2784               proceed;
2785
2786 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2787             #pragma omp critical (MagickCore_MorphologyPrimitive)
2788 #endif
2789             proceed=SetImageProgress(image,MorphologyTag,progress++,
2790               image->rows);
2791             if (proceed == MagickFalse)
2792               status=MagickFalse;
2793           }
2794       }
2795       morphology_image->type=image->type;
2796       morphology_view=DestroyCacheView(morphology_view);
2797       image_view=DestroyCacheView(image_view);
2798       for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2799         changed+=changes[i];
2800       changes=(size_t *) RelinquishMagickMemory(changes);
2801       return(status ? (ssize_t) changed : 0);
2802     }
2803   /*
2804     Normal handling of horizontal or rectangular kernels (row by row).
2805   */
2806 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2807   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2808     magick_threads(image,morphology_image,image->rows,1)
2809 #endif
2810   for (y=0; y < (ssize_t) image->rows; y++)
2811   {
2812     const int
2813       id = GetOpenMPThreadId();
2814
2815     register const Quantum
2816       *restrict p;
2817
2818     register Quantum
2819       *restrict q;
2820
2821     register ssize_t
2822       x;
2823
2824     ssize_t
2825       center;
2826
2827     if (status == MagickFalse)
2828       continue;
2829     p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,
2830       kernel->height,exception);
2831     q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2832       1,exception);
2833     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2834       {
2835         status=MagickFalse;
2836         continue;
2837       }
2838     center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
2839       GetPixelChannels(image)*offset.x);
2840     for (x=0; x < (ssize_t) image->columns; x++)
2841     {
2842       register ssize_t
2843         i;
2844
2845       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2846       {
2847         double
2848           alpha,
2849           gamma,
2850           maximum,
2851           minimum,
2852           pixel;
2853
2854         PixelChannel
2855           channel;
2856
2857         PixelTrait
2858           morphology_traits,
2859           traits;
2860
2861         register const MagickRealType
2862           *restrict k;
2863
2864         register const Quantum
2865           *restrict pixels;
2866
2867         register ssize_t
2868           u;
2869
2870         size_t
2871           count;
2872
2873         ssize_t
2874           v;
2875
2876         channel=GetPixelChannelChannel(image,i);
2877         traits=GetPixelChannelTraits(image,channel);
2878         morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2879         if ((traits == UndefinedPixelTrait) ||
2880             (morphology_traits == UndefinedPixelTrait))
2881           continue;
2882         if (((traits & CopyPixelTrait) != 0) ||
2883             (GetPixelReadMask(image,p+center) == 0))
2884           {
2885             SetPixelChannel(morphology_image,channel,p[center+i],q);
2886             continue;
2887           }
2888         pixels=p;
2889         maximum=0.0;
2890         minimum=(double) QuantumRange;
2891         count=kernel->width*kernel->height;
2892         switch (method)
2893         {
2894           case ConvolveMorphology: pixel=bias; break;
2895           case HitAndMissMorphology: pixel=(double) QuantumRange; break;
2896           case ThinningMorphology: pixel=(double) QuantumRange; break;
2897           case ThickenMorphology: pixel=(double) QuantumRange; break;
2898           case ErodeMorphology: pixel=(double) QuantumRange; break;
2899           case DilateMorphology: pixel=0.0; break;
2900           case ErodeIntensityMorphology:
2901           case DilateIntensityMorphology:
2902           case IterativeDistanceMorphology:
2903           {
2904             pixel=(double) p[center+i];
2905             break;
2906           }
2907           default: pixel=0; break;
2908         }
2909         gamma=1.0;
2910         switch (method)
2911         {
2912           case ConvolveMorphology:
2913           {
2914             /*
2915                Weighted Average of pixels using reflected kernel
2916
2917                For correct working of this operation for asymetrical kernels,
2918                the kernel needs to be applied in its reflected form.  That is
2919                its values needs to be reversed.
2920
2921                Correlation is actually the same as this but without reflecting
2922                the kernel, and thus 'lower-level' that Convolution.  However as
2923                Convolution is the more common method used, and it does not
2924                really cost us much in terms of processing to use a reflected
2925                kernel, so it is Convolution that is implemented.
2926
2927                Correlation will have its kernel reflected before calling this
2928                function to do a Convolve.
2929
2930                For more details of Correlation vs Convolution see
2931                  http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2932             */
2933             k=(&kernel->values[kernel->width*kernel->height-1]);
2934             count=0;
2935             if ((morphology_traits & BlendPixelTrait) == 0)
2936               {
2937                 /*
2938                   No alpha blending.
2939                 */
2940                 for (v=0; v < (ssize_t) kernel->height; v++)
2941                 {
2942                   for (u=0; u < (ssize_t) kernel->width; u++)
2943                   {
2944                     if (IfNaN(*k) == MagickFalse)
2945                       {
2946                         pixel+=(*k)*pixels[i];
2947                         count++;
2948                       }
2949                     k--;
2950                     pixels+=GetPixelChannels(image);
2951                   }
2952                   pixels+=(image->columns-1)*GetPixelChannels(image);
2953                 }
2954                 break;
2955               }
2956             /*
2957               Alpha blending.
2958             */
2959             gamma=0.0;
2960             for (v=0; v < (ssize_t) kernel->height; v++)
2961             {
2962               for (u=0; u < (ssize_t) kernel->width; u++)
2963               {
2964                 if (IfNaN(*k) == MagickFalse)
2965                   {
2966                     alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2967                     pixel+=alpha*(*k)*pixels[i];
2968                     gamma+=alpha*(*k);
2969                     count++;
2970                   }
2971                 k--;
2972                 pixels+=GetPixelChannels(image);
2973               }
2974               pixels+=(image->columns-1)*GetPixelChannels(image);
2975             }
2976             break;
2977           }
2978           case ErodeMorphology:
2979           {
2980             /*
2981               Minimum value within kernel neighbourhood.
2982
2983               The kernel is not reflected for this operation.  In normal
2984               Greyscale Morphology, the kernel value should be added
2985               to the real value, this is currently not done, due to the
2986               nature of the boolean kernels being used.
2987             */
2988             k=kernel->values;
2989             for (v=0; v < (ssize_t) kernel->height; v++)
2990             {
2991               for (u=0; u < (ssize_t) kernel->width; u++)
2992               {
2993                 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
2994                   {
2995                     if ((double) pixels[i] < pixel)
2996                       pixel=(double) pixels[i];
2997                   }
2998                 k++;
2999                 pixels+=GetPixelChannels(image);
3000               }
3001               pixels+=(image->columns-1)*GetPixelChannels(image);
3002             }
3003             break;
3004           }
3005           case DilateMorphology:
3006           {
3007             /*
3008                Maximum value within kernel neighbourhood.
3009
3010                For correct working of this operation for asymetrical kernels,
3011                the kernel needs to be applied in its reflected form.  That is
3012                its values needs to be reversed.
3013
3014                In normal Greyscale Morphology, the kernel value should be
3015                added to the real value, this is currently not done, due to the
3016                nature of the boolean kernels being used.
3017             */
3018             count=0;
3019             k=(&kernel->values[kernel->width*kernel->height-1]);
3020             for (v=0; v < (ssize_t) kernel->height; v++)
3021             {
3022               for (u=0; u < (ssize_t) kernel->width; u++)
3023               {
3024                 if ((IfNaN(*k) == MagickFalse) && (*k > 0.5))
3025                   {
3026                     if ((double) pixels[i] > pixel)
3027                       pixel=(double) pixels[i];
3028                   }
3029                 k--;
3030                 pixels+=GetPixelChannels(image);
3031               }
3032               pixels+=(image->columns-1)*GetPixelChannels(image);
3033             }
3034             break;
3035           }
3036           case HitAndMissMorphology:
3037           case ThinningMorphology:
3038           case ThickenMorphology:
3039           {
3040             /*
3041                Minimum of foreground pixel minus maxumum of background pixels.
3042
3043                The kernel is not reflected for this operation, and consists
3044                of both foreground and background pixel neighbourhoods, 0.0 for
3045                background, and 1.0 for foreground with either Nan or 0.5 values
3046                for don't care.
3047
3048                This never produces a meaningless negative result.  Such results
3049                cause Thinning/Thicken to not work correctly when used against a
3050                greyscale image.
3051             */
3052             count=0;
3053             k=kernel->values;
3054             for (v=0; v < (ssize_t) kernel->height; v++)
3055             {
3056               for (u=0; u < (ssize_t) kernel->width; u++)
3057               {
3058                 if (IfNaN(*k) == MagickFalse)
3059                   {
3060                     if (*k > 0.7)
3061                       {
3062                         if ((double) pixels[i] < pixel)
3063                           pixel=(double) pixels[i];
3064                       }
3065                     else
3066                       if (*k < 0.3)
3067                         {
3068                           if ((double) pixels[i] > maximum)
3069                             maximum=(double) pixels[i];
3070                         }
3071                     count++;
3072                   }
3073                 k++;
3074                 pixels+=GetPixelChannels(image);
3075               }
3076               pixels+=(image->columns-1)*GetPixelChannels(image);
3077             }
3078             pixel-=maximum;
3079             if (pixel < 0.0)
3080               pixel=0.0;
3081             if (method ==  ThinningMorphology)
3082               pixel=(double) p[center+i]-pixel;
3083             else
3084               if (method ==  ThickenMorphology)
3085                 pixel+=(double) p[center+i]+pixel;
3086             break;
3087           }
3088           case ErodeIntensityMorphology:
3089           {
3090             /*
3091               Select pixel with minimum intensity within kernel neighbourhood.
3092
3093               The kernel is not reflected for this operation.
3094             */
3095             count=0;
3096             k=kernel->values;
3097             for (v=0; v < (ssize_t) kernel->height; v++)
3098             {
3099               for (u=0; u < (ssize_t) kernel->width; u++)
3100               {
3101                 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3102                   {
3103                     if (GetPixelIntensity(image,pixels) < minimum)
3104                       {
3105                         pixel=(double) pixels[i];
3106                         minimum=GetPixelIntensity(image,pixels);
3107                       }
3108                     count++;
3109                   }
3110                 k++;
3111                 pixels+=GetPixelChannels(image);
3112               }
3113               pixels+=(image->columns-1)*GetPixelChannels(image);
3114             }
3115             break;
3116           }
3117           case DilateIntensityMorphology:
3118           {
3119             /*
3120               Select pixel with maximum intensity within kernel neighbourhood.
3121
3122               The kernel is not reflected for this operation.
3123             */
3124             count=0;
3125             k=(&kernel->values[kernel->width*kernel->height-1]);
3126             for (v=0; v < (ssize_t) kernel->height; v++)
3127             {
3128               for (u=0; u < (ssize_t) kernel->width; u++)
3129               {
3130                 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3131                   {
3132                     if (GetPixelIntensity(image,pixels) > maximum)
3133                       {
3134                         pixel=(double) pixels[i];
3135                         maximum=GetPixelIntensity(image,pixels);
3136                       }
3137                     count++;
3138                   }
3139                 k--;
3140                 pixels+=GetPixelChannels(image);
3141               }
3142               pixels+=(image->columns-1)*GetPixelChannels(image);
3143             }
3144             break;
3145           }
3146           case IterativeDistanceMorphology:
3147           {
3148             /*
3149                Compute th iterative distance from black edge of a white image
3150                shape.  Essentually white values are decreased to the smallest
3151                'distance from edge' it can find.
3152
3153                It works by adding kernel values to the neighbourhood, and and
3154                select the minimum value found. The kernel is rotated before
3155                use, so kernel distances match resulting distances, when a user
3156                provided asymmetric kernel is applied.
3157
3158                This code is nearly identical to True GrayScale Morphology but
3159                not quite.
3160
3161                GreyDilate Kernel values added, maximum value found Kernel is
3162                rotated before use.
3163
3164                GrayErode:  Kernel values subtracted and minimum value found No
3165                kernel rotation used.
3166
3167                Note the the Iterative Distance method is essentially a
3168                GrayErode, but with negative kernel values, and kernel rotation
3169                applied.
3170             */
3171             count=0;
3172             k=(&kernel->values[kernel->width*kernel->height-1]);
3173             for (v=0; v < (ssize_t) kernel->height; v++)
3174             {
3175               for (u=0; u < (ssize_t) kernel->width; u++)
3176               {
3177                 if (IfNaN(*k) == MagickFalse)
3178                   {
3179                     if ((pixels[i]+(*k)) < pixel)
3180                       pixel=(double) pixels[i]+(*k);
3181                     count++;
3182                   }
3183                 k--;
3184                 pixels+=GetPixelChannels(image);
3185               }
3186               pixels+=(image->columns-1)*GetPixelChannels(image);
3187             }
3188             break;
3189           }
3190           case UndefinedMorphology:
3191           default:
3192             break;
3193         }
3194         if (fabs(pixel-p[center+i]) > MagickEpsilon)
3195           changes[id]++;
3196         gamma=PerceptibleReciprocal(gamma);
3197         if (count != 0)
3198           gamma*=(double) kernel->height*kernel->width/count;
3199         SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*pixel),q);
3200       }
3201       p+=GetPixelChannels(image);
3202       q+=GetPixelChannels(morphology_image);
3203     }
3204     if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3205       status=MagickFalse;
3206     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3207       {
3208         MagickBooleanType
3209           proceed;
3210
3211 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3212         #pragma omp critical (MagickCore_MorphologyPrimitive)
3213 #endif
3214         proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3215         if (proceed == MagickFalse)
3216           status=MagickFalse;
3217       }
3218   }
3219   morphology_view=DestroyCacheView(morphology_view);
3220   image_view=DestroyCacheView(image_view);
3221   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
3222     changed+=changes[i];
3223   changes=(size_t *) RelinquishMagickMemory(changes);
3224   return(status ? (ssize_t) changed : -1);
3225 }
3226
3227 /*
3228   This is almost identical to the MorphologyPrimative() function above, but
3229   applies the primitive directly to the actual image using two passes, once in
3230   each direction, with the results of the previous (and current) row being
3231   re-used.
3232
3233   That is after each row is 'Sync'ed' into the image, the next row makes use of
3234   those values as part of the calculation of the next row.  It repeats, but
3235   going in the oppisite (bottom-up) direction.
3236
3237   Because of this 're-use of results' this function can not make use of multi-
3238   threaded, parellel processing.
3239 */
3240 static ssize_t MorphologyPrimitiveDirect(Image *image,
3241   const MorphologyMethod method,const KernelInfo *kernel,
3242   ExceptionInfo *exception)
3243 {
3244   CacheView
3245     *morphology_view,
3246     *image_view;
3247
3248   MagickBooleanType
3249     status;
3250
3251   MagickOffsetType
3252     progress;
3253
3254   OffsetInfo
3255     offset;
3256
3257   size_t
3258     width,
3259     changed;
3260
3261   ssize_t
3262     y;
3263
3264   assert(image != (Image *) NULL);
3265   assert(image->signature == MagickSignature);
3266   assert(kernel != (KernelInfo *) NULL);
3267   assert(kernel->signature == MagickSignature);
3268   assert(exception != (ExceptionInfo *) NULL);
3269   assert(exception->signature == MagickSignature);
3270   status=MagickTrue;
3271   changed=0;
3272   progress=0;
3273   switch(method)
3274   {
3275     case DistanceMorphology:
3276     case VoronoiMorphology:
3277     {
3278       /*
3279         Kernel reflected about origin.
3280       */
3281       offset.x=(ssize_t) kernel->width-kernel->x-1;
3282       offset.y=(ssize_t) kernel->height-kernel->y-1;
3283       break;
3284     }
3285     default:
3286     {
3287       offset.x=kernel->x;
3288       offset.y=kernel->y;
3289       break;
3290     }
3291   }
3292   /*
3293     Two views into same image, do not thread.
3294   */
3295   image_view=AcquireVirtualCacheView(image,exception);
3296   morphology_view=AcquireAuthenticCacheView(image,exception);
3297   width=image->columns+kernel->width-1;
3298   for (y=0; y < (ssize_t) image->rows; y++)
3299   {
3300     register const Quantum
3301       *restrict p;
3302
3303     register Quantum
3304       *restrict q;
3305
3306     register ssize_t
3307       x;
3308
3309     ssize_t
3310       center;
3311
3312     /*
3313       Read virtual pixels, and authentic pixels, from the same image!  We read
3314       using virtual to get virtual pixel handling, but write back into the same
3315       image.
3316
3317       Only top half of kernel is processed as we do a single pass downward
3318       through the image iterating the distance function as we go.
3319     */
3320     if (status == MagickFalse)
3321       continue;
3322     p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t)
3323       offset.y+1,exception);
3324     q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3325       exception);
3326     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3327       {
3328         status=MagickFalse;
3329         continue;
3330       }
3331     center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
3332       GetPixelChannels(image)*offset.x);
3333     for (x=0; x < (ssize_t) image->columns; x++)
3334     {
3335       register ssize_t
3336         i;
3337
3338       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3339       {
3340         double
3341           pixel;
3342
3343         PixelTrait
3344           traits;
3345
3346         register const MagickRealType
3347           *restrict k;
3348
3349         register const Quantum
3350           *restrict pixels;
3351
3352         register ssize_t
3353           u;
3354
3355         ssize_t
3356           v;
3357
3358         traits=GetPixelChannelTraits(image,(PixelChannel) i);
3359         if (traits == UndefinedPixelTrait)
3360           continue;
3361         if (((traits & CopyPixelTrait) != 0) ||
3362             (GetPixelReadMask(image,p+center) == 0))
3363           continue;
3364         pixels=p;
3365         pixel=(double) QuantumRange;
3366         switch (method)
3367         {
3368           case DistanceMorphology:
3369           {
3370             k=(&kernel->values[kernel->width*kernel->height-1]);
3371             for (v=0; v <= offset.y; v++)
3372             {
3373               for (u=0; u < (ssize_t) kernel->width; u++)
3374               {
3375                 if (IfNaN(*k) == MagickFalse)
3376                   {
3377                     if ((pixels[i]+(*k)) < pixel)
3378                       pixel=(double) pixels[i]+(*k);
3379                   }
3380                 k--;
3381                 pixels+=GetPixelChannels(image);
3382               }
3383               pixels+=(image->columns-1)*GetPixelChannels(image);
3384             }
3385             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3386             pixels=q-offset.x*GetPixelChannels(image);
3387             for (u=0; u < offset.x; u++)
3388             {
3389               if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3390                 {
3391                   if ((pixels[i]+(*k)) < pixel)
3392                     pixel=(double) pixels[i]+(*k);
3393                 }
3394               k--;
3395               pixels+=GetPixelChannels(image);
3396             }
3397             break;
3398           }
3399           case VoronoiMorphology:
3400           {
3401             k=(&kernel->values[kernel->width*kernel->height-1]);
3402             for (v=0; v < offset.y; v++)
3403             {
3404               for (u=0; u < (ssize_t) kernel->width; u++)
3405               {
3406                 if (IfNaN(*k) == MagickFalse)
3407                   {
3408                     if ((pixels[i]+(*k)) < pixel)
3409                       pixel=(double) pixels[i]+(*k);
3410                   }
3411                 k--;
3412                 pixels+=GetPixelChannels(image);
3413               }
3414               pixels+=(image->columns-1)*GetPixelChannels(image);
3415             }
3416             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3417             pixels=q-offset.x*GetPixelChannels(image);
3418             for (u=0; u < offset.x; u++)
3419             {
3420               if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3421                 {
3422                   if ((pixels[i]+(*k)) < pixel)
3423                     pixel=(double) pixels[i]+(*k);
3424                 }
3425               k--;
3426               pixels+=GetPixelChannels(image);
3427             }
3428             break;
3429           }
3430           default:
3431             break;
3432         }
3433         if (fabs(pixel-q[i]) > MagickEpsilon)
3434           changed++;
3435         q[i]=ClampToQuantum(pixel);
3436       }
3437       p+=GetPixelChannels(image);
3438       q+=GetPixelChannels(image);
3439     }
3440     if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3441       status=MagickFalse;
3442     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3443       {
3444         MagickBooleanType
3445           proceed;
3446
3447         proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3448         if (proceed == MagickFalse)
3449           status=MagickFalse;
3450       }
3451   }
3452   morphology_view=DestroyCacheView(morphology_view);
3453   image_view=DestroyCacheView(image_view);
3454   /*
3455     Do the reverse pass through the image.
3456   */
3457   image_view=AcquireVirtualCacheView(image,exception);
3458   morphology_view=AcquireAuthenticCacheView(image,exception);
3459   for (y=(ssize_t) image->rows-1; y >= 0; y--)
3460   {
3461     register const Quantum
3462       *restrict p;
3463
3464     register Quantum
3465       *restrict q;
3466
3467     register ssize_t
3468       x;
3469
3470     ssize_t
3471       center;
3472
3473     /*
3474        Read virtual pixels, and authentic pixels, from the same image.  We
3475        read using virtual to get virtual pixel handling, but write back
3476        into the same image.
3477
3478        Only the bottom half of the kernel is processed as we up the image.
3479     */
3480     if (status == MagickFalse)
3481       continue;
3482     p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t)
3483       kernel->y+1,exception);
3484     q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3485       exception);
3486     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3487       {
3488         status=MagickFalse;
3489         continue;
3490       }
3491     p+=(image->columns-1)*GetPixelChannels(image);
3492     q+=(image->columns-1)*GetPixelChannels(image);
3493     center=(ssize_t) (offset.x*GetPixelChannels(image));
3494     for (x=(ssize_t) image->columns-1; x >= 0; x--)
3495     {
3496       register ssize_t
3497         i;
3498
3499       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3500       {
3501         double
3502           pixel;
3503
3504         PixelTrait
3505           traits;
3506
3507         register const MagickRealType
3508           *restrict k;
3509
3510         register const Quantum
3511           *restrict pixels;
3512
3513         register ssize_t
3514           u;
3515
3516         ssize_t
3517           v;
3518
3519         traits=GetPixelChannelTraits(image,(PixelChannel) i);
3520         if (traits == UndefinedPixelTrait)
3521           continue;
3522         if (((traits & CopyPixelTrait) != 0) ||
3523             (GetPixelReadMask(image,p+center) == 0))
3524           continue;
3525         pixels=p;
3526         pixel=(double) QuantumRange;
3527         switch (method)
3528         {
3529           case DistanceMorphology:
3530           {
3531             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3532             for (v=offset.y; v < (ssize_t) kernel->height; v++)
3533             {
3534               for (u=0; u < (ssize_t) kernel->width; u++)
3535               {
3536                 if (IfNaN(*k) == MagickFalse)
3537                   {
3538                     if ((pixels[i]+(*k)) < pixel)
3539                       pixel=(double) pixels[i]+(*k);
3540                   }
3541                 k--;
3542                 pixels+=GetPixelChannels(image);
3543               }
3544               pixels+=(image->columns-1)*GetPixelChannels(image);
3545             }
3546             k=(&kernel->values[kernel->width*kernel->y+kernel->x-1]);
3547             pixels=q;
3548             for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3549             {
3550               pixels+=GetPixelChannels(image);
3551               if ((IfNaN(*k) == MagickFalse) &&
3552                   ((x+u-offset.x) < (ssize_t) image->columns))
3553                 {
3554                   if ((pixels[i]+(*k)) < pixel)
3555                     pixel=(double) pixels[i]+(*k);
3556                 }
3557               k--;
3558             }
3559             break;
3560           }
3561           case VoronoiMorphology:
3562           {
3563             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3564             for (v=offset.y; v < (ssize_t) kernel->height; v++)
3565             {
3566               for (u=0; u < (ssize_t) kernel->width; u++)
3567               {
3568                 if (IfNaN(*k) == MagickFalse)
3569                   {
3570                     if ((pixels[i]+(*k)) < pixel)
3571                       pixel=(double) pixels[i]+(*k);
3572                   }
3573                 k--;
3574                 pixels+=GetPixelChannels(image);
3575               }
3576               pixels+=(image->columns-1)*GetPixelChannels(image);
3577             }
3578             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3579             pixels=q;
3580             for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3581             {
3582               pixels+=GetPixelChannels(image);
3583               if ((IfNaN(*k) == MagickFalse) &&
3584                   ((x+u-offset.x) < (ssize_t) image->columns))
3585                 {
3586                   if ((pixels[i]+(*k)) < pixel)
3587                     pixel=(double) pixels[i]+(*k);
3588                 }
3589               k--;
3590             }
3591             break;
3592           }
3593           default:
3594             break;
3595         }
3596         if (fabs(pixel-q[i]) > MagickEpsilon)
3597           changed++;
3598         q[i]=ClampToQuantum(pixel);
3599       }
3600       p-=GetPixelChannels(image);
3601       q-=GetPixelChannels(image);
3602     }
3603     if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3604       status=MagickFalse;
3605     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3606       {
3607         MagickBooleanType
3608           proceed;
3609
3610         proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3611         if (proceed == MagickFalse)
3612           status=MagickFalse;
3613       }
3614   }
3615   morphology_view=DestroyCacheView(morphology_view);
3616   image_view=DestroyCacheView(image_view);
3617   return(status ? (ssize_t) changed : -1);
3618 }
3619
3620 /*
3621   Apply a Morphology by calling one of the above low level primitive
3622   application functions.  This function handles any iteration loops,
3623   composition or re-iteration of results, and compound morphology methods that
3624   is based on multiple low-level (staged) morphology methods.
3625
3626   Basically this provides the complex glue between the requested morphology
3627   method and raw low-level implementation (above).
3628 */
3629 MagickPrivate Image *MorphologyApply(const Image *image,
3630   const MorphologyMethod method, const ssize_t iterations,
3631   const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3632   ExceptionInfo *exception)
3633 {
3634   CompositeOperator
3635     curr_compose;
3636
3637   Image
3638     *curr_image,    /* Image we are working with or iterating */
3639     *work_image,    /* secondary image for primitive iteration */
3640     *save_image,    /* saved image - for 'edge' method only */
3641     *rslt_image;    /* resultant image - after multi-kernel handling */
3642
3643   KernelInfo
3644     *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3645     *norm_kernel,      /* the current normal un-reflected kernel */
3646     *rflt_kernel,      /* the current reflected kernel (if needed) */
3647     *this_kernel;      /* the kernel being applied */
3648
3649   MorphologyMethod
3650     primitive;      /* the current morphology primitive being applied */
3651
3652   CompositeOperator
3653     rslt_compose;   /* multi-kernel compose method for results to use */
3654
3655   MagickBooleanType
3656     special,        /* do we use a direct modify function? */
3657     verbose;        /* verbose output of results */
3658
3659   size_t
3660     method_loop,    /* Loop 1: number of compound method iterations (norm 1) */
3661     method_limit,   /*         maximum number of compound method iterations */
3662     kernel_number,  /* Loop 2: the kernel number being applied */
3663     stage_loop,     /* Loop 3: primitive loop for compound morphology */
3664     stage_limit,    /*         how many primitives are in this compound */
3665     kernel_loop,    /* Loop 4: iterate the kernel over image */
3666     kernel_limit,   /*         number of times to iterate kernel */
3667     count,          /* total count of primitive steps applied */
3668     kernel_changed, /* total count of changed using iterated kernel */
3669     method_changed; /* total count of changed over method iteration */
3670
3671   ssize_t
3672     changed;        /* number pixels changed by last primitive operation */
3673
3674   char
3675     v_info[MaxTextExtent];
3676
3677   assert(image != (Image *) NULL);
3678   assert(image->signature == MagickSignature);
3679   assert(kernel != (KernelInfo *) NULL);
3680   assert(kernel->signature == MagickSignature);
3681   assert(exception != (ExceptionInfo *) NULL);
3682   assert(exception->signature == MagickSignature);
3683
3684   count = 0;      /* number of low-level morphology primitives performed */
3685   if ( iterations == 0 )
3686     return((Image *) NULL);   /* null operation - nothing to do! */
3687
3688   kernel_limit = (size_t) iterations;
3689   if ( iterations < 0 )  /* negative interations = infinite (well alomst) */
3690      kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3691
3692   verbose = IsStringTrue(GetImageArtifact(image,"debug"));
3693
3694   /* initialise for cleanup */
3695   curr_image = (Image *) image;
3696   curr_compose = image->compose;
3697   (void) curr_compose;
3698   work_image = save_image = rslt_image = (Image *) NULL;
3699   reflected_kernel = (KernelInfo *) NULL;
3700
3701   /* Initialize specific methods
3702    * + which loop should use the given iteratations
3703    * + how many primitives make up the compound morphology
3704    * + multi-kernel compose method to use (by default)
3705    */
3706   method_limit = 1;       /* just do method once, unless otherwise set */
3707   stage_limit = 1;        /* assume method is not a compound */
3708   special = MagickFalse;   /* assume it is NOT a direct modify primitive */
3709   rslt_compose = compose; /* and we are composing multi-kernels as given */
3710   switch( method ) {
3711     case SmoothMorphology:  /* 4 primitive compound morphology */
3712       stage_limit = 4;
3713       break;
3714     case OpenMorphology:    /* 2 primitive compound morphology */
3715     case OpenIntensityMorphology:
3716     case TopHatMorphology:
3717     case CloseMorphology:
3718     case CloseIntensityMorphology:
3719     case BottomHatMorphology:
3720     case EdgeMorphology:
3721       stage_limit = 2;
3722       break;
3723     case HitAndMissMorphology:
3724       rslt_compose = LightenCompositeOp;  /* Union of multi-kernel results */
3725       /* FALL THUR */
3726     case ThinningMorphology:
3727     case ThickenMorphology:
3728       method_limit = kernel_limit;  /* iterate the whole method */
3729       kernel_limit = 1;             /* do not do kernel iteration  */
3730       break;
3731     case DistanceMorphology:
3732     case VoronoiMorphology:
3733       special = MagickTrue;         /* use special direct primative */
3734       break;
3735     default:
3736       break;
3737   }
3738
3739   /* Apply special methods with special requirments
3740   ** For example, single run only, or post-processing requirements
3741   */
3742   if ( special != MagickFalse )
3743     {
3744       rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3745       if (rslt_image == (Image *) NULL)
3746         goto error_cleanup;
3747       if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3748         goto error_cleanup;
3749
3750       changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception);
3751
3752       if ( IfMagickTrue(verbose) )
3753         (void) (void) FormatLocaleFile(stderr,
3754           "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3755           CommandOptionToMnemonic(MagickMorphologyOptions, method),
3756           1.0,0.0,1.0, (double) changed);
3757
3758       if ( changed < 0 )
3759         goto error_cleanup;
3760
3761       if ( method == VoronoiMorphology ) {
3762         /* Preserve the alpha channel of input image - but turned it off */
3763         (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3764           exception);
3765         (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3766           MagickTrue,0,0,exception);
3767         (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3768           exception);
3769       }
3770       goto exit_cleanup;
3771     }
3772
3773   /* Handle user (caller) specified multi-kernel composition method */
3774   if ( compose != UndefinedCompositeOp )
3775     rslt_compose = compose;  /* override default composition for method */
3776   if ( rslt_compose == UndefinedCompositeOp )
3777     rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3778
3779   /* Some methods require a reflected kernel to use with primitives.
3780    * Create the reflected kernel for those methods. */
3781   switch ( method ) {
3782     case CorrelateMorphology:
3783     case CloseMorphology:
3784     case CloseIntensityMorphology:
3785     case BottomHatMorphology:
3786     case SmoothMorphology:
3787       reflected_kernel = CloneKernelInfo(kernel);
3788       if (reflected_kernel == (KernelInfo *) NULL)
3789         goto error_cleanup;
3790       RotateKernelInfo(reflected_kernel,180);
3791       break;
3792     default:
3793       break;
3794   }
3795
3796   /* Loops around more primitive morpholgy methods
3797   **  erose, dilate, open, close, smooth, edge, etc...
3798   */
3799   /* Loop 1:  iterate the compound method */
3800   method_loop = 0;
3801   method_changed = 1;
3802   while ( method_loop < method_limit && method_changed > 0 ) {
3803     method_loop++;
3804     method_changed = 0;
3805
3806     /* Loop 2:  iterate over each kernel in a multi-kernel list */
3807     norm_kernel = (KernelInfo *) kernel;
3808     this_kernel = (KernelInfo *) kernel;
3809     rflt_kernel = reflected_kernel;
3810
3811     kernel_number = 0;
3812     while ( norm_kernel != NULL ) {
3813
3814       /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3815       stage_loop = 0;          /* the compound morphology stage number */
3816       while ( stage_loop < stage_limit ) {
3817         stage_loop++;   /* The stage of the compound morphology */
3818
3819         /* Select primitive morphology for this stage of compound method */
3820         this_kernel = norm_kernel; /* default use unreflected kernel */
3821         primitive = method;        /* Assume method is a primitive */
3822         switch( method ) {
3823           case ErodeMorphology:      /* just erode */
3824           case EdgeInMorphology:     /* erode and image difference */
3825             primitive = ErodeMorphology;
3826             break;
3827           case DilateMorphology:     /* just dilate */
3828           case EdgeOutMorphology:    /* dilate and image difference */
3829             primitive = DilateMorphology;
3830             break;
3831           case OpenMorphology:       /* erode then dialate */
3832           case TopHatMorphology:     /* open and image difference */
3833             primitive = ErodeMorphology;
3834             if ( stage_loop == 2 )
3835               primitive = DilateMorphology;
3836             break;
3837           case OpenIntensityMorphology:
3838             primitive = ErodeIntensityMorphology;
3839             if ( stage_loop == 2 )
3840               primitive = DilateIntensityMorphology;
3841             break;
3842           case CloseMorphology:      /* dilate, then erode */
3843           case BottomHatMorphology:  /* close and image difference */
3844             this_kernel = rflt_kernel; /* use the reflected kernel */
3845             primitive = DilateMorphology;
3846             if ( stage_loop == 2 )
3847               primitive = ErodeMorphology;
3848             break;
3849           case CloseIntensityMorphology:
3850             this_kernel = rflt_kernel; /* use the reflected kernel */
3851             primitive = DilateIntensityMorphology;
3852             if ( stage_loop == 2 )
3853               primitive = ErodeIntensityMorphology;
3854             break;
3855           case SmoothMorphology:         /* open, close */
3856             switch ( stage_loop ) {
3857               case 1: /* start an open method, which starts with Erode */
3858                 primitive = ErodeMorphology;
3859                 break;
3860               case 2:  /* now Dilate the Erode */
3861                 primitive = DilateMorphology;
3862                 break;
3863               case 3:  /* Reflect kernel a close */
3864                 this_kernel = rflt_kernel; /* use the reflected kernel */
3865                 primitive = DilateMorphology;
3866                 break;
3867               case 4:  /* Finish the Close */
3868                 this_kernel = rflt_kernel; /* use the reflected kernel */
3869                 primitive = ErodeMorphology;
3870                 break;
3871             }
3872             break;
3873           case EdgeMorphology:        /* dilate and erode difference */
3874             primitive = DilateMorphology;
3875             if ( stage_loop == 2 ) {
3876               save_image = curr_image;      /* save the image difference */
3877               curr_image = (Image *) image;
3878               primitive = ErodeMorphology;
3879             }
3880             break;
3881           case CorrelateMorphology:
3882             /* A Correlation is a Convolution with a reflected kernel.
3883             ** However a Convolution is a weighted sum using a reflected
3884             ** kernel.  It may seem stange to convert a Correlation into a
3885             ** Convolution as the Correlation is the simplier method, but
3886             ** Convolution is much more commonly used, and it makes sense to
3887             ** implement it directly so as to avoid the need to duplicate the
3888             ** kernel when it is not required (which is typically the
3889             ** default).
3890             */
3891             this_kernel = rflt_kernel; /* use the reflected kernel */
3892             primitive = ConvolveMorphology;
3893             break;
3894           default:
3895             break;
3896         }
3897         assert( this_kernel != (KernelInfo *) NULL );
3898
3899         /* Extra information for debugging compound operations */
3900         if ( IfMagickTrue(verbose) ) {
3901           if ( stage_limit > 1 )
3902             (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
3903              CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
3904              method_loop,(double) stage_loop);
3905           else if ( primitive != method )
3906             (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
3907               CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
3908               method_loop);
3909           else
3910             v_info[0] = '\0';
3911         }
3912
3913         /* Loop 4: Iterate the kernel with primitive */
3914         kernel_loop = 0;
3915         kernel_changed = 0;
3916         changed = 1;
3917         while ( kernel_loop < kernel_limit && changed > 0 ) {
3918           kernel_loop++;     /* the iteration of this kernel */
3919
3920           /* Create a clone as the destination image, if not yet defined */
3921           if ( work_image == (Image *) NULL )
3922             {
3923               work_image=CloneImage(image,0,0,MagickTrue,exception);
3924               if (work_image == (Image *) NULL)
3925                 goto error_cleanup;
3926               if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
3927                 goto error_cleanup;
3928             }
3929
3930           /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3931           count++;
3932           changed = MorphologyPrimitive(curr_image, work_image, primitive,
3933                        this_kernel, bias, exception);
3934           if ( IfMagickTrue(verbose) ) {
3935             if ( kernel_loop > 1 )
3936               (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
3937             (void) (void) FormatLocaleFile(stderr,
3938               "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3939               v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
3940               primitive),(this_kernel == rflt_kernel ) ? "*" : "",
3941               (double) (method_loop+kernel_loop-1),(double) kernel_number,
3942               (double) count,(double) changed);
3943           }
3944           if ( changed < 0 )
3945             goto error_cleanup;
3946           kernel_changed += changed;
3947           method_changed += changed;
3948
3949           /* prepare next loop */
3950           { Image *tmp = work_image;   /* swap images for iteration */
3951             work_image = curr_image;
3952             curr_image = tmp;
3953           }
3954           if ( work_image == image )
3955             work_image = (Image *) NULL; /* replace input 'image' */
3956
3957         } /* End Loop 4: Iterate the kernel with primitive */
3958
3959         if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed )
3960           (void) FormatLocaleFile(stderr, "   Total %.20g",(double) kernel_changed);
3961         if ( IfMagickTrue(verbose) && stage_loop < stage_limit )
3962           (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
3963
3964 #if 0
3965     (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3966     (void) FormatLocaleFile(stderr, "      curr =0x%lx\n", (unsigned long)curr_image);
3967     (void) FormatLocaleFile(stderr, "      work =0x%lx\n", (unsigned long)work_image);
3968     (void) FormatLocaleFile(stderr, "      save =0x%lx\n", (unsigned long)save_image);
3969     (void) FormatLocaleFile(stderr, "      union=0x%lx\n", (unsigned long)rslt_image);
3970 #endif
3971
3972       } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3973
3974       /*  Final Post-processing for some Compound Methods
3975       **
3976       ** The removal of any 'Sync' channel flag in the Image Compositon
3977       ** below ensures the methematical compose method is applied in a
3978       ** purely mathematical way, and only to the selected channels.
3979       ** Turn off SVG composition 'alpha blending'.
3980       */
3981       switch( method ) {
3982         case EdgeOutMorphology:
3983         case EdgeInMorphology:
3984         case TopHatMorphology:
3985         case BottomHatMorphology:
3986           if ( IfMagickTrue(verbose) )
3987             (void) FormatLocaleFile(stderr,
3988               "\n%s: Difference with original image",CommandOptionToMnemonic(
3989               MagickMorphologyOptions, method) );
3990           (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
3991             MagickTrue,0,0,exception);
3992           break;
3993         case EdgeMorphology:
3994           if ( IfMagickTrue(verbose) )
3995             (void) FormatLocaleFile(stderr,
3996               "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
3997               MagickMorphologyOptions, method) );
3998           (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
3999             MagickTrue,0,0,exception);
4000           save_image = DestroyImage(save_image); /* finished with save image */
4001           break;
4002         default:
4003           break;
4004       }
4005
4006       /* multi-kernel handling:  re-iterate, or compose results */
4007       if ( kernel->next == (KernelInfo *) NULL )
4008         rslt_image = curr_image;   /* just return the resulting image */
4009       else if ( rslt_compose == NoCompositeOp )
4010         { if ( IfMagickTrue(verbose) ) {
4011             if ( this_kernel->next != (KernelInfo *) NULL )
4012               (void) FormatLocaleFile(stderr, " (re-iterate)");
4013             else
4014               (void) FormatLocaleFile(stderr, " (done)");
4015           }
4016           rslt_image = curr_image; /* return result, and re-iterate */
4017         }
4018       else if ( rslt_image == (Image *) NULL)
4019         { if ( IfMagickTrue(verbose) )
4020             (void) FormatLocaleFile(stderr, " (save for compose)");
4021           rslt_image = curr_image;
4022           curr_image = (Image *) image;  /* continue with original image */
4023         }
4024       else
4025         { /* Add the new 'current' result to the composition
4026           **
4027           ** The removal of any 'Sync' channel flag in the Image Compositon
4028           ** below ensures the methematical compose method is applied in a
4029           ** purely mathematical way, and only to the selected channels.
4030           ** IE: Turn off SVG composition 'alpha blending'.
4031           */
4032           if ( IfMagickTrue(verbose) )
4033             (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4034               CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4035           (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4036             0,0,exception);
4037           curr_image = DestroyImage(curr_image);
4038           curr_image = (Image *) image;  /* continue with original image */
4039         }
4040       if ( IfMagickTrue(verbose) )
4041         (void) FormatLocaleFile(stderr, "\n");
4042
4043       /* loop to the next kernel in a multi-kernel list */
4044       norm_kernel = norm_kernel->next;
4045       if ( rflt_kernel != (KernelInfo *) NULL )
4046         rflt_kernel = rflt_kernel->next;
4047       kernel_number++;
4048     } /* End Loop 2: Loop over each kernel */
4049
4050   } /* End Loop 1: compound method interation */
4051
4052   goto exit_cleanup;
4053
4054   /* Yes goto's are bad, but it makes cleanup lot more efficient */
4055 error_cleanup:
4056   if ( curr_image == rslt_image )
4057     curr_image = (Image *) NULL;
4058   if ( rslt_image != (Image *) NULL )
4059     rslt_image = DestroyImage(rslt_image);
4060 exit_cleanup:
4061   if ( curr_image == rslt_image || curr_image == image )
4062     curr_image = (Image *) NULL;
4063   if ( curr_image != (Image *) NULL )
4064     curr_image = DestroyImage(curr_image);
4065   if ( work_image != (Image *) NULL )
4066     work_image = DestroyImage(work_image);
4067   if ( save_image != (Image *) NULL )
4068     save_image = DestroyImage(save_image);
4069   if ( reflected_kernel != (KernelInfo *) NULL )
4070     reflected_kernel = DestroyKernelInfo(reflected_kernel);
4071   return(rslt_image);
4072 }
4073
4074 \f
4075 /*
4076 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4077 %                                                                             %
4078 %                                                                             %
4079 %                                                                             %
4080 %     M o r p h o l o g y I m a g e                                           %
4081 %                                                                             %
4082 %                                                                             %
4083 %                                                                             %
4084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4085 %
4086 %  MorphologyImage() applies a user supplied kernel to the image according to
4087 %  the given mophology method.
4088 %
4089 %  This function applies any and all user defined settings before calling
4090 %  the above internal function MorphologyApply().
4091 %
4092 %  User defined settings include...
4093 %    * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4094 %    * Kernel Scale/normalize settings             ("-define convolve:scale=??")
4095 %      This can also includes the addition of a scaled unity kernel.
4096 %    * Show Kernel being applied                   ("-define showkernel=1")
4097 %
4098 %  Other operators that do not want user supplied options interfering,
4099 %  especially "convolve:bias" and "showkernel" should use MorphologyApply()
4100 %  directly.
4101 %
4102 %  The format of the MorphologyImage method is:
4103 %
4104 %      Image *MorphologyImage(const Image *image,MorphologyMethod method,
4105 %        const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4106 %
4107 %  A description of each parameter follows:
4108 %
4109 %    o image: the image.
4110 %
4111 %    o method: the morphology method to be applied.
4112 %
4113 %    o iterations: apply the operation this many times (or no change).
4114 %                  A value of -1 means loop until no change found.
4115 %                  How this is applied may depend on the morphology method.
4116 %                  Typically this is a value of 1.
4117 %
4118 %    o kernel: An array of double representing the morphology kernel.
4119 %              Warning: kernel may be normalized for the Convolve method.
4120 %
4121 %    o exception: return any errors or warnings in this structure.
4122 %
4123 */
4124 MagickExport Image *MorphologyImage(const Image *image,
4125   const MorphologyMethod method,const ssize_t iterations,
4126   const KernelInfo *kernel,ExceptionInfo *exception)
4127 {
4128   KernelInfo
4129     *curr_kernel;
4130
4131   CompositeOperator
4132     compose;
4133
4134   Image
4135     *morphology_image;
4136
4137   double
4138     bias;
4139
4140   curr_kernel = (KernelInfo *) kernel;
4141   bias=0.0;
4142   compose = UndefinedCompositeOp;  /* use default for method */
4143
4144   /* Apply Convolve/Correlate Normalization and Scaling Factors.
4145    * This is done BEFORE the ShowKernelInfo() function is called so that
4146    * users can see the results of the 'option:convolve:scale' option.
4147    */
4148   if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4149       const char
4150         *artifact;
4151
4152       /* Get the bias value as it will be needed */
4153       artifact = GetImageArtifact(image,"convolve:bias");
4154       if ( artifact != (const char *) NULL) {
4155         if (IfMagickFalse(IsGeometry(artifact)))
4156           (void) ThrowMagickException(exception,GetMagickModule(),
4157                OptionWarning,"InvalidSetting","'%s' '%s'",
4158                "convolve:bias",artifact);
4159         else
4160           bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4161       }
4162
4163       /* Scale kernel according to user wishes */
4164       artifact = GetImageArtifact(image,"convolve:scale");
4165       if ( artifact != (const char *) NULL ) {
4166         if (IfMagickFalse(IsGeometry(artifact)))
4167           (void) ThrowMagickException(exception,GetMagickModule(),
4168                OptionWarning,"InvalidSetting","'%s' '%s'",
4169                "convolve:scale",artifact);
4170         else {
4171           if ( curr_kernel == kernel )
4172             curr_kernel = CloneKernelInfo(kernel);
4173           if (curr_kernel == (KernelInfo *) NULL)
4174             return((Image *) NULL);
4175           ScaleGeometryKernelInfo(curr_kernel, artifact);
4176         }
4177       }
4178     }
4179
4180   /* display the (normalized) kernel via stderr */
4181   if ( IfStringTrue(GetImageArtifact(image,"showkernel"))
4182     || IfStringTrue(GetImageArtifact(image,"convolve:showkernel"))
4183     || IfStringTrue(GetImageArtifact(image,"morphology:showkernel")) )
4184     ShowKernelInfo(curr_kernel);
4185
4186   /* Override the default handling of multi-kernel morphology results
4187    * If 'Undefined' use the default method
4188    * If 'None' (default for 'Convolve') re-iterate previous result
4189    * Otherwise merge resulting images using compose method given.
4190    * Default for 'HitAndMiss' is 'Lighten'.
4191    */
4192   { const char
4193       *artifact;
4194     ssize_t
4195       parse;
4196
4197     artifact = GetImageArtifact(image,"morphology:compose");
4198     if ( artifact != (const char *) NULL) {
4199       parse=ParseCommandOption(MagickComposeOptions,
4200         MagickFalse,artifact);
4201       if ( parse < 0 )
4202         (void) ThrowMagickException(exception,GetMagickModule(),
4203              OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4204              "morphology:compose",artifact);
4205       else
4206         compose=(CompositeOperator)parse;
4207     }
4208   }
4209   /* Apply the Morphology */
4210   morphology_image = MorphologyApply(image,method,iterations,
4211     curr_kernel,compose,bias,exception);
4212
4213   /* Cleanup and Exit */
4214   if ( curr_kernel != kernel )
4215     curr_kernel=DestroyKernelInfo(curr_kernel);
4216   return(morphology_image);
4217 }
4218 \f
4219 /*
4220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4221 %                                                                             %
4222 %                                                                             %
4223 %                                                                             %
4224 +     R o t a t e K e r n e l I n f o                                         %
4225 %                                                                             %
4226 %                                                                             %
4227 %                                                                             %
4228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4229 %
4230 %  RotateKernelInfo() rotates the kernel by the angle given.
4231 %
4232 %  Currently it is restricted to 90 degree angles, of either 1D kernels
4233 %  or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4234 %  It will ignore usless rotations for specific 'named' built-in kernels.
4235 %
4236 %  The format of the RotateKernelInfo method is:
4237 %
4238 %      void RotateKernelInfo(KernelInfo *kernel, double angle)
4239 %
4240 %  A description of each parameter follows:
4241 %
4242 %    o kernel: the Morphology/Convolution kernel
4243 %
4244 %    o angle: angle to rotate in degrees
4245 %
4246 % This function is currently internal to this module only, but can be exported
4247 % to other modules if needed.
4248 */
4249 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4250 {
4251   /* angle the lower kernels first */
4252   if ( kernel->next != (KernelInfo *) NULL)
4253     RotateKernelInfo(kernel->next, angle);
4254
4255   /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4256   **
4257   ** TODO: expand beyond simple 90 degree rotates, flips and flops
4258   */
4259
4260   /* Modulus the angle */
4261   angle = fmod(angle, 360.0);
4262   if ( angle < 0 )
4263     angle += 360.0;
4264
4265   if ( 337.5 < angle || angle <= 22.5 )
4266     return;   /* Near zero angle - no change! - At least not at this time */
4267
4268   /* Handle special cases */
4269   switch (kernel->type) {
4270     /* These built-in kernels are cylindrical kernels, rotating is useless */
4271     case GaussianKernel:
4272     case DoGKernel:
4273     case LoGKernel:
4274     case DiskKernel:
4275     case PeaksKernel:
4276     case LaplacianKernel:
4277     case ChebyshevKernel:
4278     case ManhattanKernel:
4279     case EuclideanKernel:
4280       return;
4281
4282     /* These may be rotatable at non-90 angles in the future */
4283     /* but simply rotating them in multiples of 90 degrees is useless */
4284     case SquareKernel:
4285     case DiamondKernel:
4286     case PlusKernel:
4287     case CrossKernel:
4288       return;
4289
4290     /* These only allows a +/-90 degree rotation (by transpose) */
4291     /* A 180 degree rotation is useless */
4292     case BlurKernel:
4293       if ( 135.0 < angle && angle <= 225.0 )
4294         return;
4295       if ( 225.0 < angle && angle <= 315.0 )
4296         angle -= 180;
4297       break;
4298
4299     default:
4300       break;
4301   }
4302   /* Attempt rotations by 45 degrees  -- 3x3 kernels only */
4303   if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4304     {
4305       if ( kernel->width == 3 && kernel->height == 3 )
4306         { /* Rotate a 3x3 square by 45 degree angle */
4307           double t  = kernel->values[0];
4308           kernel->values[0] = kernel->values[3];
4309           kernel->values[3] = kernel->values[6];
4310           kernel->values[6] = kernel->values[7];
4311           kernel->values[7] = kernel->values[8];
4312           kernel->values[8] = kernel->values[5];
4313           kernel->values[5] = kernel->values[2];
4314           kernel->values[2] = kernel->values[1];
4315           kernel->values[1] = t;
4316           /* rotate non-centered origin */
4317           if ( kernel->x != 1 || kernel->y != 1 ) {
4318             ssize_t x,y;
4319             x = (ssize_t) kernel->x-1;
4320             y = (ssize_t) kernel->y-1;
4321                  if ( x == y  ) x = 0;
4322             else if ( x == 0  ) x = -y;
4323             else if ( x == -y ) y = 0;
4324             else if ( y == 0  ) y = x;
4325             kernel->x = (ssize_t) x+1;
4326             kernel->y = (ssize_t) y+1;
4327           }
4328           angle = fmod(angle+315.0, 360.0);  /* angle reduced 45 degrees */
4329           kernel->angle = fmod(kernel->angle+45.0, 360.0);
4330         }
4331       else
4332         perror("Unable to rotate non-3x3 kernel by 45 degrees");
4333     }
4334   if ( 45.0 < fmod(angle, 180.0)  && fmod(angle,180.0) <= 135.0 )
4335     {
4336       if ( kernel->width == 1 || kernel->height == 1 )
4337         { /* Do a transpose of a 1 dimensional kernel,
4338           ** which results in a fast 90 degree rotation of some type.
4339           */
4340           ssize_t
4341             t;
4342           t = (ssize_t) kernel->width;
4343           kernel->width = kernel->height;
4344           kernel->height = (size_t) t;
4345           t = kernel->x;
4346           kernel->x = kernel->y;
4347           kernel->y = t;
4348           if ( kernel->width == 1 ) {
4349             angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
4350             kernel->angle = fmod(kernel->angle+90.0, 360.0);
4351           } else {
4352             angle = fmod(angle+90.0, 360.0);   /* angle increased 90 degrees */
4353             kernel->angle = fmod(kernel->angle+270.0, 360.0);
4354           }
4355         }
4356       else if ( kernel->width == kernel->height )
4357         { /* Rotate a square array of values by 90 degrees */
4358           { register ssize_t
4359               i,j,x,y;
4360
4361             register MagickRealType
4362               *k,t;
4363
4364             k=kernel->values;
4365             for( i=0, x=(ssize_t) kernel->width-1;  i<=x;   i++, x--)
4366               for( j=0, y=(ssize_t) kernel->height-1;  j<y;   j++, y--)
4367                 { t                    = k[i+j*kernel->width];
4368                   k[i+j*kernel->width] = k[j+x*kernel->width];
4369                   k[j+x*kernel->width] = k[x+y*kernel->width];
4370                   k[x+y*kernel->width] = k[y+i*kernel->width];
4371                   k[y+i*kernel->width] = t;
4372                 }
4373           }
4374           /* rotate the origin - relative to center of array */
4375           { register ssize_t x,y;
4376             x = (ssize_t) (kernel->x*2-kernel->width+1);
4377             y = (ssize_t) (kernel->y*2-kernel->height+1);
4378             kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4379             kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4380           }
4381           angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
4382           kernel->angle = fmod(kernel->angle+90.0, 360.0);
4383         }
4384       else
4385         perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4386     }
4387   if ( 135.0 < angle && angle <= 225.0 )
4388     {
4389       /* For a 180 degree rotation - also know as a reflection
4390        * This is actually a very very common operation!
4391        * Basically all that is needed is a reversal of the kernel data!
4392        * And a reflection of the origon
4393        */
4394       MagickRealType
4395         t;
4396
4397       register MagickRealType
4398         *k;
4399
4400       ssize_t
4401         i,
4402         j;
4403
4404       k=kernel->values;
4405       j=(ssize_t) (kernel->width*kernel->height-1);
4406       for (i=0;  i < j;  i++, j--)
4407         t=k[i],  k[i]=k[j],  k[j]=t;
4408
4409       kernel->x = (ssize_t) kernel->width  - kernel->x - 1;
4410       kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4411       angle = fmod(angle-180.0, 360.0);   /* angle+180 degrees */
4412       kernel->angle = fmod(kernel->angle+180.0, 360.0);
4413     }
4414   /* At this point angle should at least between -45 (315) and +45 degrees
4415    * In the future some form of non-orthogonal angled rotates could be
4416    * performed here, posibily with a linear kernel restriction.
4417    */
4418
4419   return;
4420 }
4421 \f
4422 /*
4423 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4424 %                                                                             %
4425 %                                                                             %
4426 %                                                                             %
4427 %     S c a l e G e o m e t r y K e r n e l I n f o                           %
4428 %                                                                             %
4429 %                                                                             %
4430 %                                                                             %
4431 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4432 %
4433 %  ScaleGeometryKernelInfo() takes a geometry argument string, typically
4434 %  provided as a  "-set option:convolve:scale {geometry}" user setting,
4435 %  and modifies the kernel according to the parsed arguments of that setting.
4436 %
4437 %  The first argument (and any normalization flags) are passed to
4438 %  ScaleKernelInfo() to scale/normalize the kernel.  The second argument
4439 %  is then passed to UnityAddKernelInfo() to add a scled unity kernel
4440 %  into the scaled/normalized kernel.
4441 %
4442 %  The format of the ScaleGeometryKernelInfo method is:
4443 %
4444 %      void ScaleGeometryKernelInfo(KernelInfo *kernel,
4445 %        const double scaling_factor,const MagickStatusType normalize_flags)
4446 %
4447 %  A description of each parameter follows:
4448 %
4449 %    o kernel: the Morphology/Convolution kernel to modify
4450 %
4451 %    o geometry:
4452 %             The geometry string to parse, typically from the user provided
4453 %             "-set option:convolve:scale {geometry}" setting.
4454 %
4455 */
4456 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4457   const char *geometry)
4458 {
4459   MagickStatusType
4460     flags;
4461
4462   GeometryInfo
4463     args;
4464
4465   SetGeometryInfo(&args);
4466   flags = ParseGeometry(geometry, &args);
4467
4468 #if 0
4469   /* For Debugging Geometry Input */
4470   (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4471        flags, args.rho, args.sigma, args.xi, args.psi );
4472 #endif
4473
4474   if ( (flags & PercentValue) != 0 )      /* Handle Percentage flag*/
4475     args.rho *= 0.01,  args.sigma *= 0.01;
4476
4477   if ( (flags & RhoValue) == 0 )          /* Set Defaults for missing args */
4478     args.rho = 1.0;
4479   if ( (flags & SigmaValue) == 0 )
4480     args.sigma = 0.0;
4481
4482   /* Scale/Normalize the input kernel */
4483   ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4484
4485   /* Add Unity Kernel, for blending with original */
4486   if ( (flags & SigmaValue) != 0 )
4487     UnityAddKernelInfo(kernel, args.sigma);
4488
4489   return;
4490 }
4491 /*
4492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4493 %                                                                             %
4494 %                                                                             %
4495 %                                                                             %
4496 %     S c a l e K e r n e l I n f o                                           %
4497 %                                                                             %
4498 %                                                                             %
4499 %                                                                             %
4500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4501 %
4502 %  ScaleKernelInfo() scales the given kernel list by the given amount, with or
4503 %  without normalization of the sum of the kernel values (as per given flags).
4504 %
4505 %  By default (no flags given) the values within the kernel is scaled
4506 %  directly using given scaling factor without change.
4507 %
4508 %  If either of the two 'normalize_flags' are given the kernel will first be
4509 %  normalized and then further scaled by the scaling factor value given.
4510 %
4511 %  Kernel normalization ('normalize_flags' given) is designed to ensure that
4512 %  any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4513 %  morphology methods will fall into -1.0 to +1.0 range.  Note that for
4514 %  non-HDRI versions of IM this may cause images to have any negative results
4515 %  clipped, unless some 'bias' is used.
4516 %
4517 %  More specifically.  Kernels which only contain positive values (such as a
4518 %  'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4519 %  ensuring a 0.0 to +1.0 output range for non-HDRI images.
4520 %
4521 %  For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4522 %  the kernel will be scaled by the absolute of the sum of kernel values, so
4523 %  that it will generally fall within the +/- 1.0 range.
4524 %
4525 %  For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4526 %  will be scaled by just the sum of the postive values, so that its output
4527 %  range will again fall into the  +/- 1.0 range.
4528 %
4529 %  For special kernels designed for locating shapes using 'Correlate', (often
4530 %  only containing +1 and -1 values, representing foreground/brackground
4531 %  matching) a special normalization method is provided to scale the positive
4532 %  values separately to those of the negative values, so the kernel will be
4533 %  forced to become a zero-sum kernel better suited to such searches.
4534 %
4535 %  WARNING: Correct normalization of the kernel assumes that the '*_range'
4536 %  attributes within the kernel structure have been correctly set during the
4537 %  kernels creation.
4538 %
4539 %  NOTE: The values used for 'normalize_flags' have been selected specifically
4540 %  to match the use of geometry options, so that '!' means NormalizeValue, '^'
4541 %  means CorrelateNormalizeValue.  All other GeometryFlags values are ignored.
4542 %
4543 %  The format of the ScaleKernelInfo method is:
4544 %
4545 %      void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4546 %               const MagickStatusType normalize_flags )
4547 %
4548 %  A description of each parameter follows:
4549 %
4550 %    o kernel: the Morphology/Convolution kernel
4551 %
4552 %    o scaling_factor:
4553 %             multiply all values (after normalization) by this factor if not
4554 %             zero.  If the kernel is normalized regardless of any flags.
4555 %
4556 %    o normalize_flags:
4557 %             GeometryFlags defining normalization method to use.
4558 %             specifically: NormalizeValue, CorrelateNormalizeValue,
4559 %                           and/or PercentValue
4560 %
4561 */
4562 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4563   const double scaling_factor,const GeometryFlags normalize_flags)
4564 {
4565   register double
4566     pos_scale,
4567     neg_scale;
4568
4569   register ssize_t
4570     i;
4571
4572   /* do the other kernels in a multi-kernel list first */
4573   if ( kernel->next != (KernelInfo *) NULL)
4574     ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4575
4576   /* Normalization of Kernel */
4577   pos_scale = 1.0;
4578   if ( (normalize_flags&NormalizeValue) != 0 ) {
4579     if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4580       /* non-zero-summing kernel (generally positive) */
4581       pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4582     else
4583       /* zero-summing kernel */
4584       pos_scale = kernel->positive_range;
4585   }
4586   /* Force kernel into a normalized zero-summing kernel */
4587   if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4588     pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4589                  ? kernel->positive_range : 1.0;
4590     neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4591                  ? -kernel->negative_range : 1.0;
4592   }
4593   else
4594     neg_scale = pos_scale;
4595
4596   /* finialize scaling_factor for positive and negative components */
4597   pos_scale = scaling_factor/pos_scale;
4598   neg_scale = scaling_factor/neg_scale;
4599
4600   for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4601     if ( ! IfNaN(kernel->values[i]) )
4602       kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4603
4604   /* convolution output range */
4605   kernel->positive_range *= pos_scale;
4606   kernel->negative_range *= neg_scale;
4607   /* maximum and minimum values in kernel */
4608   kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4609   kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4610
4611   /* swap kernel settings if user's scaling factor is negative */
4612   if ( scaling_factor < MagickEpsilon ) {
4613     double t;
4614     t = kernel->positive_range;
4615     kernel->positive_range = kernel->negative_range;
4616     kernel->negative_range = t;
4617     t = kernel->maximum;
4618     kernel->maximum = kernel->minimum;
4619     kernel->minimum = 1;
4620   }
4621
4622   return;
4623 }
4624 \f
4625 /*
4626 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4627 %                                                                             %
4628 %                                                                             %
4629 %                                                                             %
4630 %     S h o w K e r n e l I n f o                                             %
4631 %                                                                             %
4632 %                                                                             %
4633 %                                                                             %
4634 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4635 %
4636 %  ShowKernelInfo() outputs the details of the given kernel defination to
4637 %  standard error, generally due to a users 'showkernel' option request.
4638 %
4639 %  The format of the ShowKernel method is:
4640 %
4641 %      void ShowKernelInfo(const KernelInfo *kernel)
4642 %
4643 %  A description of each parameter follows:
4644 %
4645 %    o kernel: the Morphology/Convolution kernel
4646 %
4647 */
4648 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4649 {
4650   const KernelInfo
4651     *k;
4652
4653   size_t
4654     c, i, u, v;
4655
4656   for (c=0, k=kernel;  k != (KernelInfo *) NULL;  c++, k=k->next ) {
4657
4658     (void) FormatLocaleFile(stderr, "Kernel");
4659     if ( kernel->next != (KernelInfo *) NULL )
4660       (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4661     (void) FormatLocaleFile(stderr, " \"%s",
4662           CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4663     if ( fabs(k->angle) >= MagickEpsilon )
4664       (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4665     (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4666       k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4667     (void) FormatLocaleFile(stderr,
4668           " with values from %.*lg to %.*lg\n",
4669           GetMagickPrecision(), k->minimum,
4670           GetMagickPrecision(), k->maximum);
4671     (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4672           GetMagickPrecision(), k->negative_range,
4673           GetMagickPrecision(), k->positive_range);
4674     if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4675       (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4676     else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4677       (void) FormatLocaleFile(stderr, " (Normalized)\n");
4678     else
4679       (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4680           GetMagickPrecision(), k->positive_range+k->negative_range);
4681     for (i=v=0; v < k->height; v++) {
4682       (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4683       for (u=0; u < k->width; u++, i++)
4684         if ( IfNaN(k->values[i]) )
4685           (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4686         else
4687           (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4688               GetMagickPrecision(), (double) k->values[i]);
4689       (void) FormatLocaleFile(stderr,"\n");
4690     }
4691   }
4692 }
4693 \f
4694 /*
4695 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4696 %                                                                             %
4697 %                                                                             %
4698 %                                                                             %
4699 %     U n i t y A d d K e r n a l I n f o                                     %
4700 %                                                                             %
4701 %                                                                             %
4702 %                                                                             %
4703 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4704 %
4705 %  UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4706 %  to the given pre-scaled and normalized Kernel.  This in effect adds that
4707 %  amount of the original image into the resulting convolution kernel.  This
4708 %  value is usually provided by the user as a percentage value in the
4709 %  'convolve:scale' setting.
4710 %
4711 %  The resulting effect is to convert the defined kernels into blended
4712 %  soft-blurs, unsharp kernels or into sharpening kernels.
4713 %
4714 %  The format of the UnityAdditionKernelInfo method is:
4715 %
4716 %      void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4717 %
4718 %  A description of each parameter follows:
4719 %
4720 %    o kernel: the Morphology/Convolution kernel
4721 %
4722 %    o scale:
4723 %             scaling factor for the unity kernel to be added to
4724 %             the given kernel.
4725 %
4726 */
4727 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4728   const double scale)
4729 {
4730   /* do the other kernels in a multi-kernel list first */
4731   if ( kernel->next != (KernelInfo *) NULL)
4732     UnityAddKernelInfo(kernel->next, scale);
4733
4734   /* Add the scaled unity kernel to the existing kernel */
4735   kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4736   CalcKernelMetaData(kernel);  /* recalculate the meta-data */
4737
4738   return;
4739 }
4740 \f
4741 /*
4742 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4743 %                                                                             %
4744 %                                                                             %
4745 %                                                                             %
4746 %     Z e r o K e r n e l N a n s                                             %
4747 %                                                                             %
4748 %                                                                             %
4749 %                                                                             %
4750 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4751 %
4752 %  ZeroKernelNans() replaces any special 'nan' value that may be present in
4753 %  the kernel with a zero value.  This is typically done when the kernel will
4754 %  be used in special hardware (GPU) convolution processors, to simply
4755 %  matters.
4756 %
4757 %  The format of the ZeroKernelNans method is:
4758 %
4759 %      void ZeroKernelNans (KernelInfo *kernel)
4760 %
4761 %  A description of each parameter follows:
4762 %
4763 %    o kernel: the Morphology/Convolution kernel
4764 %
4765 */
4766 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4767 {
4768   register size_t
4769     i;
4770
4771   /* do the other kernels in a multi-kernel list first */
4772   if ( kernel->next != (KernelInfo *) NULL)
4773     ZeroKernelNans(kernel->next);
4774
4775   for (i=0; i < (kernel->width*kernel->height); i++)
4776     if ( IfNaN(kernel->values[i]) )
4777       kernel->values[i] = 0.0;
4778
4779   return;
4780 }