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