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