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