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