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 (((morphology_traits & CopyPixelTrait) != 0) ||
2714                 (GetPixelReadMask(image,p+center) == 0))
2715               {
2716                 SetPixelChannel(morphology_image,channel,p[center+i],q);
2717                 continue;
2718               }
2719             k=(&kernel->values[kernel->width*kernel->height-1]);
2720             pixels=p;
2721             pixel=bias;
2722             count=0;
2723             if ((morphology_traits & BlendPixelTrait) == 0)
2724               {
2725                 /*
2726                   No alpha blending.
2727                 */
2728                 for (v=0; v < (ssize_t) kernel->height; v++)
2729                 {
2730                   for (u=0; u < (ssize_t) kernel->width; u++)
2731                   {
2732                     if (IfNaN(*k) == MagickFalse)
2733                       {
2734                         pixel+=(*k)*pixels[i];
2735                         count++;
2736                       }
2737                     k--;
2738                     pixels+=GetPixelChannels(image);
2739                   }
2740                 }
2741                 if (fabs(pixel-p[center+i]) > MagickEpsilon)
2742                   changes[id]++;
2743                 gamma=(double) kernel->height*kernel->width/count;
2744                 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2745                   pixel),q);
2746                 continue;
2747               }
2748             /*
2749               Alpha blending.
2750             */
2751             gamma=0.0;
2752             for (v=0; v < (ssize_t) kernel->height; v++)
2753             {
2754               for (u=0; u < (ssize_t) kernel->width; u++)
2755               {
2756                 if (IfNaN(*k) == MagickFalse)
2757                   {
2758                     alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2759                     pixel+=(*k)*alpha*pixels[i];
2760                     gamma+=(*k)*alpha;
2761                     count++;
2762                   }
2763                 k--;
2764                 pixels+=GetPixelChannels(image);
2765               }
2766             }
2767             if (fabs(pixel-p[center+i]) > MagickEpsilon)
2768               changes[id]++;
2769             gamma=PerceptibleReciprocal(gamma);
2770             gamma*=(double) kernel->height*kernel->width/count;
2771             SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2772               pixel),q);
2773           }
2774           p+=GetPixelChannels(image);
2775           q+=GetPixelChannels(morphology_image);
2776         }
2777         if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2778           status=MagickFalse;
2779         if (image->progress_monitor != (MagickProgressMonitor) NULL)
2780           {
2781             MagickBooleanType
2782               proceed;
2783
2784 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2785             #pragma omp critical (MagickCore_MorphologyPrimitive)
2786 #endif
2787             proceed=SetImageProgress(image,MorphologyTag,progress++,
2788               image->rows);
2789             if (proceed == MagickFalse)
2790               status=MagickFalse;
2791           }
2792       }
2793       morphology_image->type=image->type;
2794       morphology_view=DestroyCacheView(morphology_view);
2795       image_view=DestroyCacheView(image_view);
2796       for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2797         changed+=changes[i];
2798       changes=(size_t *) RelinquishMagickMemory(changes);
2799       return(status ? (ssize_t) changed : 0);
2800     }
2801   /*
2802     Normal handling of horizontal or rectangular kernels (row by row).
2803   */
2804 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2805   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2806     magick_threads(image,morphology_image,image->rows,1)
2807 #endif
2808   for (y=0; y < (ssize_t) image->rows; y++)
2809   {
2810     const int
2811       id = GetOpenMPThreadId();
2812
2813     register const Quantum
2814       *restrict p;
2815
2816     register Quantum
2817       *restrict q;
2818
2819     register ssize_t
2820       x;
2821
2822     ssize_t
2823       center;
2824
2825     if (status == MagickFalse)
2826       continue;
2827     p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,
2828       kernel->height,exception);
2829     q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2830       1,exception);
2831     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2832       {
2833         status=MagickFalse;
2834         continue;
2835       }
2836     center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
2837       GetPixelChannels(image)*offset.x);
2838     for (x=0; x < (ssize_t) image->columns; x++)
2839     {
2840       register ssize_t
2841         i;
2842
2843       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2844       {
2845         double
2846           alpha,
2847           gamma,
2848           maximum,
2849           minimum,
2850           pixel;
2851
2852         PixelChannel
2853           channel;
2854
2855         PixelTrait
2856           morphology_traits,
2857           traits;
2858
2859         register const MagickRealType
2860           *restrict k;
2861
2862         register const Quantum
2863           *restrict pixels;
2864
2865         register ssize_t
2866           u;
2867
2868         size_t
2869           count;
2870
2871         ssize_t
2872           v;
2873
2874         channel=GetPixelChannelChannel(image,i);
2875         traits=GetPixelChannelTraits(image,channel);
2876         morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2877         if ((traits == UndefinedPixelTrait) ||
2878             (morphology_traits == UndefinedPixelTrait))
2879           continue;
2880         if (((morphology_traits & CopyPixelTrait) != 0) ||
2881             (GetPixelReadMask(image,p+center) == 0))
2882           {
2883             SetPixelChannel(morphology_image,channel,p[center+i],q);
2884             continue;
2885           }
2886         pixels=p;
2887         maximum=0.0;
2888         minimum=(double) QuantumRange;
2889         count=kernel->width*kernel->height;
2890         switch (method)
2891         {
2892           case ConvolveMorphology: pixel=bias; break;
2893           case HitAndMissMorphology: pixel=(double) QuantumRange; break;
2894           case ThinningMorphology: pixel=(double) QuantumRange; break;
2895           case ThickenMorphology: pixel=(double) QuantumRange; break;
2896           case ErodeMorphology: pixel=(double) QuantumRange; break;
2897           case DilateMorphology: pixel=0.0; break;
2898           case ErodeIntensityMorphology:
2899           case DilateIntensityMorphology:
2900           case IterativeDistanceMorphology:
2901           {
2902             pixel=(double) p[center+i];
2903             break;
2904           }
2905           default: pixel=0; break;
2906         }
2907         gamma=1.0;
2908         switch (method)
2909         {
2910           case ConvolveMorphology:
2911           {
2912             /*
2913                Weighted Average of pixels using reflected kernel
2914
2915                For correct working of this operation for asymetrical
2916                kernels, the kernel needs to be applied in its reflected form.
2917                That is its values needs to be reversed.
2918
2919                Correlation is actually the same as this but without reflecting
2920                the kernel, and thus 'lower-level' that Convolution.  However
2921                as Convolution is the more common method used, and it does not
2922                really cost us much in terms of processing to use a reflected
2923                kernel, so it is Convolution that is implemented.
2924
2925                Correlation will have its kernel reflected before calling
2926                this function to do a Convolve.
2927
2928                For more details of Correlation vs Convolution see
2929                  http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2930             */
2931             k=(&kernel->values[kernel->width*kernel->height-1]);
2932             count=0;
2933             if ((morphology_traits & BlendPixelTrait) == 0)
2934               {
2935                 /*
2936                   No alpha blending.
2937                 */
2938                 for (v=0; v < (ssize_t) kernel->height; v++)
2939                 {
2940                   for (u=0; u < (ssize_t) kernel->width; u++)
2941                   {
2942                     if (IfNaN(*k) == MagickFalse)
2943                       {
2944                         pixel+=(*k)*pixels[i];
2945                         count++;
2946                       }
2947                     k--;
2948                     pixels+=GetPixelChannels(image);
2949                   }
2950                   pixels+=(image->columns-1)*GetPixelChannels(image);
2951                 }
2952                 break;
2953               }
2954             /*
2955               Alpha blending.
2956             */
2957             for (v=0; v < (ssize_t) kernel->height; v++)
2958             {
2959               for (u=0; u < (ssize_t) kernel->width; u++)
2960               {
2961                 if (IfNaN(*k) == MagickFalse)
2962                   {
2963                     alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2964                     pixel+=(*k)*alpha*pixels[i];
2965                     gamma+=(*k)*alpha;
2966                     count++;
2967                   }
2968                 k--;
2969                 pixels+=GetPixelChannels(image);
2970               }
2971               pixels+=(image->columns-1)*GetPixelChannels(image);
2972             }
2973             break;
2974           }
2975           case ErodeMorphology:
2976           {
2977             /*
2978               Minimum value within kernel neighbourhood.
2979
2980               The kernel is not reflected for this operation.  In normal
2981               Greyscale Morphology, the kernel value should be added
2982               to the real value, this is currently not done, due to the
2983               nature of the boolean kernels being used.
2984             */
2985             k=kernel->values;
2986             for (v=0; v < (ssize_t) kernel->height; v++)
2987             {
2988               for (u=0; u < (ssize_t) kernel->width; u++)
2989               {
2990                 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
2991                   {
2992                     if ((double) pixels[i] < pixel)
2993                       pixel=(double) pixels[i];
2994                   }
2995                 k++;
2996                 pixels+=GetPixelChannels(image);
2997               }
2998               pixels+=(image->columns-1)*GetPixelChannels(image);
2999             }
3000             break;
3001           }
3002           case DilateMorphology:
3003           {
3004             /*
3005                Maximum value within kernel neighbourhood.
3006
3007                For correct working of this operation for asymetrical kernels,
3008                the kernel needs to be applied in its reflected form.  That is
3009                its values needs to be reversed.
3010
3011                In normal Greyscale Morphology, the kernel value should be
3012                added to the real value, this is currently not done, due to the
3013                nature of the boolean kernels being used.
3014             */
3015             count=0;
3016             k=(&kernel->values[kernel->width*kernel->height-1]);
3017             for (v=0; v < (ssize_t) kernel->height; v++)
3018             {
3019               for (u=0; u < (ssize_t) kernel->width; u++)
3020               {
3021                 if ((IfNaN(*k) == MagickFalse) && (*k > 0.5))
3022                   {
3023                     if ((double) pixels[i] > pixel)
3024                       pixel=(double) pixels[i];
3025                     count++;
3026                   }
3027                 k--;
3028                 pixels+=GetPixelChannels(image);
3029               }
3030               pixels+=(image->columns-1)*GetPixelChannels(image);
3031             }
3032             break;
3033           }
3034           case HitAndMissMorphology:
3035           case ThinningMorphology:
3036           case ThickenMorphology:
3037           {
3038             /*
3039                Minimum of foreground pixel minus maxumum of background pixels.
3040
3041                The kernel is not reflected for this operation, and consists
3042                of both foreground and background pixel neighbourhoods, 0.0 for
3043                background, and 1.0 for foreground with either Nan or 0.5 values
3044                for don't care.
3045
3046                This never produces a meaningless negative result.  Such results
3047                cause Thinning/Thicken to not work correctly when used against a
3048                greyscale image.
3049             */
3050             count=0;
3051             k=kernel->values;
3052             for (v=0; v < (ssize_t) kernel->height; v++)
3053             {
3054               for (u=0; u < (ssize_t) kernel->width; u++)
3055               {
3056                 if (IfNaN(*k) == MagickFalse)
3057                   {
3058                     if (*k > 0.7)
3059                       {
3060                         if ((double) pixels[i] < pixel)
3061                           pixel=(double) pixels[i];
3062                       }
3063                     else
3064                       if (*k < 0.3)
3065                         {
3066                           if ((double) pixels[i] > maximum)
3067                             maximum=(double) pixels[i];
3068                         }
3069                     count++;
3070                   }
3071                 k++;
3072                 pixels+=GetPixelChannels(image);
3073               }
3074               pixels+=(image->columns-1)*GetPixelChannels(image);
3075             }
3076             pixel-=maximum;
3077             if (pixel < 0.0)
3078               pixel=0.0;
3079             if (method ==  ThinningMorphology)
3080               pixel=(double) p[center+i]-pixel;
3081             else
3082               if (method ==  ThickenMorphology)
3083                 pixel+=(double) p[center+i]+pixel;
3084             break;
3085           }
3086           case ErodeIntensityMorphology:
3087           {
3088             /*
3089               Select pixel with minimum intensity within kernel neighbourhood.
3090
3091               The kernel is not reflected for this operation.
3092             */
3093             count=0;
3094             k=kernel->values;
3095             for (v=0; v < (ssize_t) kernel->height; v++)
3096             {
3097               for (u=0; u < (ssize_t) kernel->width; u++)
3098               {
3099                 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3100                   {
3101                     if (GetPixelIntensity(image,pixels) < minimum)
3102                       {
3103                         pixel=(double) pixels[i];
3104                         minimum=GetPixelIntensity(image,pixels);
3105                       }
3106                     count++;
3107                   }
3108                 k++;
3109                 pixels+=GetPixelChannels(image);
3110               }
3111               pixels+=(image->columns-1)*GetPixelChannels(image);
3112             }
3113             break;
3114           }
3115           case DilateIntensityMorphology:
3116           {
3117             /*
3118               Select pixel with maximum intensity within kernel neighbourhood.
3119
3120               The kernel is not reflected for this operation.
3121             */
3122             count=0;
3123             k=(&kernel->values[kernel->width*kernel->height-1]);
3124             for (v=0; v < (ssize_t) kernel->height; v++)
3125             {
3126               for (u=0; u < (ssize_t) kernel->width; u++)
3127               {
3128                 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3129                   {
3130                     if (GetPixelIntensity(image,pixels) > maximum)
3131                       {
3132                         pixel=(double) pixels[i];
3133                         maximum=GetPixelIntensity(image,pixels);
3134                       }
3135                     count++;
3136                   }
3137                 k--;
3138                 pixels+=GetPixelChannels(image);
3139               }
3140               pixels+=(image->columns-1)*GetPixelChannels(image);
3141             }
3142             break;
3143           }
3144           case IterativeDistanceMorphology:
3145           {
3146             /*
3147                Compute th iterative distance from black edge of a white image
3148                shape.  Essentually white values are decreased to the smallest
3149                'distance from edge' it can find.
3150
3151                It works by adding kernel values to the neighbourhood, and and
3152                select the minimum value found. The kernel is rotated before
3153                use, so kernel distances match resulting distances, when a user
3154                provided asymmetric kernel is applied.
3155
3156                This code is nearly identical to True GrayScale Morphology but
3157                not quite.
3158
3159                GreyDilate Kernel values added, maximum value found Kernel is
3160                rotated before use.
3161
3162                GrayErode:  Kernel values subtracted and minimum value found No
3163                kernel rotation used.
3164
3165                Note the the Iterative Distance method is essentially a
3166                GrayErode, but with negative kernel values, and kernel rotation
3167                applied.
3168             */
3169             count=0;
3170             k=(&kernel->values[kernel->width*kernel->height-1]);
3171             for (v=0; v < (ssize_t) kernel->height; v++)
3172             {
3173               for (u=0; u < (ssize_t) kernel->width; u++)
3174               {
3175                 if (IfNaN(*k) == MagickFalse)
3176                   {
3177                     if ((pixels[i]+(*k)) < pixel)
3178                       pixel=(double) pixels[i]+(*k);
3179                     count++;
3180                   }
3181                 k--;
3182                 pixels+=GetPixelChannels(image);
3183               }
3184               pixels+=(image->columns-1)*GetPixelChannels(image);
3185             }
3186             break;
3187           }
3188           case UndefinedMorphology:
3189           default:
3190             break;
3191         }
3192         if (fabs(pixel-p[center+i]) > MagickEpsilon)
3193           changes[id]++;
3194         gamma=PerceptibleReciprocal(gamma);
3195         gamma*=(double) kernel->height*kernel->width/count;
3196         SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*pixel),q);
3197       }
3198       p+=GetPixelChannels(image);
3199       q+=GetPixelChannels(morphology_image);
3200     }
3201     if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3202       status=MagickFalse;
3203     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3204       {
3205         MagickBooleanType
3206           proceed;
3207
3208 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3209         #pragma omp critical (MagickCore_MorphologyPrimitive)
3210 #endif
3211         proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3212         if (proceed == MagickFalse)
3213           status=MagickFalse;
3214       }
3215   }
3216   morphology_view=DestroyCacheView(morphology_view);
3217   image_view=DestroyCacheView(image_view);
3218   for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
3219     changed+=changes[i];
3220   changes=(size_t *) RelinquishMagickMemory(changes);
3221   return(status ? (ssize_t) changed : -1);
3222 }
3223
3224 /*
3225   This is almost identical to the MorphologyPrimative() function above, but
3226   applies the primitive directly to the actual image using two passes, once in
3227   each direction, with the results of the previous (and current) row being
3228   re-used.
3229
3230   That is after each row is 'Sync'ed' into the image, the next row makes use of
3231   those values as part of the calculation of the next row.  It repeats, but
3232   going in the oppisite (bottom-up) direction.
3233
3234   Because of this 're-use of results' this function can not make use of multi-
3235   threaded, parellel processing.
3236 */
3237 static ssize_t MorphologyPrimitiveDirect(Image *image,
3238   const MorphologyMethod method,const KernelInfo *kernel,
3239   ExceptionInfo *exception)
3240 {
3241   CacheView
3242     *morphology_view,
3243     *image_view;
3244
3245   MagickBooleanType
3246     status;
3247
3248   MagickOffsetType
3249     progress;
3250
3251   OffsetInfo
3252     offset;
3253
3254   size_t
3255     width,
3256     changed;
3257
3258   ssize_t
3259     y;
3260
3261   assert(image != (Image *) NULL);
3262   assert(image->signature == MagickSignature);
3263   assert(kernel != (KernelInfo *) NULL);
3264   assert(kernel->signature == MagickSignature);
3265   assert(exception != (ExceptionInfo *) NULL);
3266   assert(exception->signature == MagickSignature);
3267   status=MagickTrue;
3268   changed=0;
3269   progress=0;
3270   switch(method)
3271   {
3272     case DistanceMorphology:
3273     case VoronoiMorphology:
3274     {
3275       /*
3276         Kernel reflected about origin.
3277       */
3278       offset.x=(ssize_t) kernel->width-kernel->x-1;
3279       offset.y=(ssize_t) kernel->height-kernel->y-1;
3280       break;
3281     }
3282     default:
3283     {
3284       offset.x=kernel->x;
3285       offset.y=kernel->y;
3286       break;
3287     }
3288   }
3289   /*
3290     Two views into same image, do not thread.
3291   */
3292   image_view=AcquireVirtualCacheView(image,exception);
3293   morphology_view=AcquireAuthenticCacheView(image,exception);
3294   width=image->columns+kernel->width-1;
3295   for (y=0; y < (ssize_t) image->rows; y++)
3296   {
3297     register const Quantum
3298       *restrict p;
3299
3300     register Quantum
3301       *restrict q;
3302
3303     register ssize_t
3304       x;
3305
3306     ssize_t
3307       center;
3308
3309     /*
3310       Read virtual pixels, and authentic pixels, from the same image!  We read
3311       using virtual to get virtual pixel handling, but write back into the same
3312       image.
3313
3314       Only top half of kernel is processed as we do a single pass downward
3315       through the image iterating the distance function as we go.
3316     */
3317     if (status == MagickFalse)
3318       continue;
3319     p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t)
3320       offset.y+1,exception);
3321     q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3322       exception);
3323     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3324       {
3325         status=MagickFalse;
3326         continue;
3327       }
3328     center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
3329       GetPixelChannels(image)*offset.x);
3330     for (x=0; x < (ssize_t) image->columns; x++)
3331     {
3332       register ssize_t
3333         i;
3334
3335       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3336       {
3337         double
3338           pixel;
3339
3340         PixelTrait
3341           traits;
3342
3343         register const MagickRealType
3344           *restrict k;
3345
3346         register const Quantum
3347           *restrict pixels;
3348
3349         register ssize_t
3350           u;
3351
3352         ssize_t
3353           v;
3354
3355         traits=GetPixelChannelTraits(image,(PixelChannel) i);
3356         if (traits == UndefinedPixelTrait)
3357           continue;
3358         if (((traits & CopyPixelTrait) != 0) ||
3359             (GetPixelReadMask(image,p+center) == 0))
3360           continue;
3361         pixels=p;
3362         pixel=(double) QuantumRange;
3363         switch (method)
3364         {
3365           case DistanceMorphology:
3366           {
3367             k=(&kernel->values[kernel->width*kernel->height-1]);
3368             for (v=0; v <= offset.y; v++)
3369             {
3370               for (u=0; u < (ssize_t) kernel->width; u++)
3371               {
3372                 if (IfNaN(*k) == MagickFalse)
3373                   {
3374                     if ((pixels[i]+(*k)) < pixel)
3375                       pixel=(double) pixels[i]+(*k);
3376                   }
3377                 k--;
3378                 pixels+=GetPixelChannels(image);
3379               }
3380               pixels+=(image->columns-1)*GetPixelChannels(image);
3381             }
3382             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3383             pixels=q-offset.x*GetPixelChannels(image);
3384             for (u=0; u < offset.x; u++)
3385             {
3386               if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3387                 {
3388                   if ((pixels[i]+(*k)) < pixel)
3389                     pixel=(double) pixels[i]+(*k);
3390                 }
3391               k--;
3392               pixels+=GetPixelChannels(image);
3393             }
3394             break;
3395           }
3396           case VoronoiMorphology:
3397           {
3398             k=(&kernel->values[kernel->width*kernel->height-1]);
3399             for (v=0; v < offset.y; v++)
3400             {
3401               for (u=0; u < (ssize_t) kernel->width; u++)
3402               {
3403                 if (IfNaN(*k) == MagickFalse)
3404                   {
3405                     if ((pixels[i]+(*k)) < pixel)
3406                       pixel=(double) pixels[i]+(*k);
3407                   }
3408                 k--;
3409                 pixels+=GetPixelChannels(image);
3410               }
3411               pixels+=(image->columns-1)*GetPixelChannels(image);
3412             }
3413             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3414             pixels=q-offset.x*GetPixelChannels(image);
3415             for (u=0; u < offset.x; u++)
3416             {
3417               if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3418                 {
3419                   if ((pixels[i]+(*k)) < pixel)
3420                     pixel=(double) pixels[i]+(*k);
3421                 }
3422               k--;
3423               pixels+=GetPixelChannels(image);
3424             }
3425             break;
3426           }
3427           default:
3428             break;
3429         }
3430         if (fabs(pixel-q[i]) > MagickEpsilon)
3431           changed++;
3432         q[i]=ClampToQuantum(pixel);
3433       }
3434       p+=GetPixelChannels(image);
3435       q+=GetPixelChannels(image);
3436     }
3437     if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3438       status=MagickFalse;
3439     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3440       {
3441         MagickBooleanType
3442           proceed;
3443
3444         proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3445         if (proceed == MagickFalse)
3446           status=MagickFalse;
3447       }
3448   }
3449   morphology_view=DestroyCacheView(morphology_view);
3450   image_view=DestroyCacheView(image_view);
3451   /*
3452     Do the reverse pass through the image.
3453   */
3454   image_view=AcquireVirtualCacheView(image,exception);
3455   morphology_view=AcquireAuthenticCacheView(image,exception);
3456   for (y=(ssize_t) image->rows-1; y >= 0; y--)
3457   {
3458     register const Quantum
3459       *restrict p;
3460
3461     register Quantum
3462       *restrict q;
3463
3464     register ssize_t
3465       x;
3466
3467     ssize_t
3468       center;
3469
3470     /*
3471        Read virtual pixels, and authentic pixels, from the same image.  We
3472        read using virtual to get virtual pixel handling, but write back
3473        into the same image.
3474
3475        Only the bottom half of the kernel is processed as we up the image.
3476     */
3477     if (status == MagickFalse)
3478       continue;
3479     p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t)
3480       kernel->y+1,exception);
3481     q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3482       exception);
3483     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3484       {
3485         status=MagickFalse;
3486         continue;
3487       }
3488     p+=(image->columns-1)*GetPixelChannels(image);
3489     q+=(image->columns-1)*GetPixelChannels(image);
3490     center=(ssize_t) (offset.x*GetPixelChannels(image));
3491     for (x=(ssize_t) image->columns-1; x >= 0; x--)
3492     {
3493       register ssize_t
3494         i;
3495
3496       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3497       {
3498         double
3499           pixel;
3500
3501         PixelTrait
3502           traits;
3503
3504         register const MagickRealType
3505           *restrict k;
3506
3507         register const Quantum
3508           *restrict pixels;
3509
3510         register ssize_t
3511           u;
3512
3513         ssize_t
3514           v;
3515
3516         traits=GetPixelChannelTraits(image,(PixelChannel) i);
3517         if (traits == UndefinedPixelTrait)
3518           continue;
3519         if (((traits & CopyPixelTrait) != 0) ||
3520             (GetPixelReadMask(image,p+center) == 0))
3521           continue;
3522         pixels=p;
3523         pixel=(double) QuantumRange;
3524         switch (method)
3525         {
3526           case DistanceMorphology:
3527           {
3528             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3529             for (v=offset.y; v < (ssize_t) kernel->height; v++)
3530             {
3531               for (u=0; u < (ssize_t) kernel->width; u++)
3532               {
3533                 if (IfNaN(*k) == MagickFalse)
3534                   {
3535                     if ((pixels[i]+(*k)) < pixel)
3536                       pixel=(double) pixels[i]+(*k);
3537                   }
3538                 k--;
3539                 pixels+=GetPixelChannels(image);
3540               }
3541               pixels+=(image->columns-1)*GetPixelChannels(image);
3542             }
3543             k=(&kernel->values[kernel->width*kernel->y+kernel->x-1]);
3544             pixels=q-offset.x*GetPixelChannels(image);
3545             for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3546             {
3547               if ((IfNaN(*k) == MagickFalse) &&
3548                   ((x+u-offset.x) < (ssize_t) image->columns))
3549                 {
3550                   if ((pixels[i]+(*k)) < pixel)
3551                     pixel=(double) pixels[i]+(*k);
3552                 }
3553               k--;
3554               pixels+=GetPixelChannels(image);
3555             }
3556             break;
3557           }
3558           case VoronoiMorphology:
3559           {
3560             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3561             for (v=offset.y; v < (ssize_t) kernel->height; v++)
3562             {
3563               for (u=0; u < (ssize_t) kernel->width; u++)
3564               {
3565                 if (IfNaN(*k) == MagickFalse)
3566                   {
3567                     if ((pixels[i]+(*k)) < pixel)
3568                       pixel=(double) pixels[i]+(*k);
3569                   }
3570                 k--;
3571                 pixels+=GetPixelChannels(image);
3572               }
3573               pixels+=(image->columns-1)*GetPixelChannels(image);
3574             }
3575             k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3576             pixels=q-offset.x*GetPixelChannels(image);
3577             for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3578             {
3579               if ((IfNaN(*k) == MagickFalse) &&
3580                   ((x+u-offset.x) < (ssize_t) image->columns))
3581                 {
3582                   if ((pixels[i]+(*k)) < pixel)
3583                     pixel=(double) pixels[i]+(*k);
3584                 }
3585               k--;
3586               pixels+=GetPixelChannels(image);
3587             }
3588             break;
3589           }
3590           default:
3591             break;
3592         }
3593         if (fabs(pixel-q[i]) > MagickEpsilon)
3594           changed++;
3595         q[i]=ClampToQuantum(pixel);
3596       }
3597       p-=GetPixelChannels(image);
3598       q-=GetPixelChannels(image);
3599     }
3600     if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3601       status=MagickFalse;
3602     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3603       {
3604         MagickBooleanType
3605           proceed;
3606
3607         proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3608         if (proceed == MagickFalse)
3609           status=MagickFalse;
3610       }
3611   }
3612   morphology_view=DestroyCacheView(morphology_view);
3613   image_view=DestroyCacheView(image_view);
3614   return(status ? (ssize_t) changed : -1);
3615 }
3616
3617 /*
3618   Apply a Morphology by calling one of the above low level primitive
3619   application functions.  This function handles any iteration loops,
3620   composition or re-iteration of results, and compound morphology methods that
3621   is based on multiple low-level (staged) morphology methods.
3622
3623   Basically this provides the complex glue between the requested morphology
3624   method and raw low-level implementation (above).
3625 */
3626 MagickPrivate Image *MorphologyApply(const Image *image,
3627   const MorphologyMethod method, const ssize_t iterations,
3628   const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3629   ExceptionInfo *exception)
3630 {
3631   CompositeOperator
3632     curr_compose;
3633
3634   Image
3635     *curr_image,    /* Image we are working with or iterating */
3636     *work_image,    /* secondary image for primitive iteration */
3637     *save_image,    /* saved image - for 'edge' method only */
3638     *rslt_image;    /* resultant image - after multi-kernel handling */
3639
3640   KernelInfo
3641     *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3642     *norm_kernel,      /* the current normal un-reflected kernel */
3643     *rflt_kernel,      /* the current reflected kernel (if needed) */
3644     *this_kernel;      /* the kernel being applied */
3645
3646   MorphologyMethod
3647     primitive;      /* the current morphology primitive being applied */
3648
3649   CompositeOperator
3650     rslt_compose;   /* multi-kernel compose method for results to use */
3651
3652   MagickBooleanType
3653     special,        /* do we use a direct modify function? */
3654     verbose;        /* verbose output of results */
3655
3656   size_t
3657     method_loop,    /* Loop 1: number of compound method iterations (norm 1) */
3658     method_limit,   /*         maximum number of compound method iterations */
3659     kernel_number,  /* Loop 2: the kernel number being applied */
3660     stage_loop,     /* Loop 3: primitive loop for compound morphology */
3661     stage_limit,    /*         how many primitives are in this compound */
3662     kernel_loop,    /* Loop 4: iterate the kernel over image */
3663     kernel_limit,   /*         number of times to iterate kernel */
3664     count,          /* total count of primitive steps applied */
3665     kernel_changed, /* total count of changed using iterated kernel */
3666     method_changed; /* total count of changed over method iteration */
3667
3668   ssize_t
3669     changed;        /* number pixels changed by last primitive operation */
3670
3671   char
3672     v_info[80];
3673
3674   assert(image != (Image *) NULL);
3675   assert(image->signature == MagickSignature);
3676   assert(kernel != (KernelInfo *) NULL);
3677   assert(kernel->signature == MagickSignature);
3678   assert(exception != (ExceptionInfo *) NULL);
3679   assert(exception->signature == MagickSignature);
3680
3681   count = 0;      /* number of low-level morphology primitives performed */
3682   if ( iterations == 0 )
3683     return((Image *)NULL);   /* null operation - nothing to do! */
3684
3685   kernel_limit = (size_t) iterations;
3686   if ( iterations < 0 )  /* negative interations = infinite (well alomst) */
3687      kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3688
3689   verbose = IsStringTrue(GetImageArtifact(image,"debug"));
3690
3691   /* initialise for cleanup */
3692   curr_image = (Image *) image;
3693   curr_compose = image->compose;
3694   (void) curr_compose;
3695   work_image = save_image = rslt_image = (Image *) NULL;
3696   reflected_kernel = (KernelInfo *) NULL;
3697
3698   /* Initialize specific methods
3699    * + which loop should use the given iteratations
3700    * + how many primitives make up the compound morphology
3701    * + multi-kernel compose method to use (by default)
3702    */
3703   method_limit = 1;       /* just do method once, unless otherwise set */
3704   stage_limit = 1;        /* assume method is not a compound */
3705   special = MagickFalse;   /* assume it is NOT a direct modify primitive */
3706   rslt_compose = compose; /* and we are composing multi-kernels as given */
3707   switch( method ) {
3708     case SmoothMorphology:  /* 4 primitive compound morphology */
3709       stage_limit = 4;
3710       break;
3711     case OpenMorphology:    /* 2 primitive compound morphology */
3712     case OpenIntensityMorphology:
3713     case TopHatMorphology:
3714     case CloseMorphology:
3715     case CloseIntensityMorphology:
3716     case BottomHatMorphology:
3717     case EdgeMorphology:
3718       stage_limit = 2;
3719       break;
3720     case HitAndMissMorphology:
3721       rslt_compose = LightenCompositeOp;  /* Union of multi-kernel results */
3722       /* FALL THUR */
3723     case ThinningMorphology:
3724     case ThickenMorphology:
3725       method_limit = kernel_limit;  /* iterate the whole method */
3726       kernel_limit = 1;             /* do not do kernel iteration  */
3727       break;
3728     case DistanceMorphology:
3729     case VoronoiMorphology:
3730       special = MagickTrue;         /* use special direct primative */
3731       break;
3732     default:
3733       break;
3734   }
3735
3736   /* Apply special methods with special requirments
3737   ** For example, single run only, or post-processing requirements
3738   */
3739   if ( special != MagickFalse )
3740     {
3741       rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3742       if (rslt_image == (Image *) NULL)
3743         goto error_cleanup;
3744       if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3745         goto error_cleanup;
3746
3747       changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception);
3748
3749       if ( IfMagickTrue(verbose) )
3750         (void) (void) FormatLocaleFile(stderr,
3751           "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3752           CommandOptionToMnemonic(MagickMorphologyOptions, method),
3753           1.0,0.0,1.0, (double) changed);
3754
3755       if ( changed < 0 )
3756         goto error_cleanup;
3757
3758       if ( method == VoronoiMorphology ) {
3759         /* Preserve the alpha channel of input image - but turned it off */
3760         (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3761           exception);
3762         (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3763           MagickTrue,0,0,exception);
3764         (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3765           exception);
3766       }
3767       goto exit_cleanup;
3768     }
3769
3770   /* Handle user (caller) specified multi-kernel composition method */
3771   if ( compose != UndefinedCompositeOp )
3772     rslt_compose = compose;  /* override default composition for method */
3773   if ( rslt_compose == UndefinedCompositeOp )
3774     rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3775
3776   /* Some methods require a reflected kernel to use with primitives.
3777    * Create the reflected kernel for those methods. */
3778   switch ( method ) {
3779     case CorrelateMorphology:
3780     case CloseMorphology:
3781     case CloseIntensityMorphology:
3782     case BottomHatMorphology:
3783     case SmoothMorphology:
3784       reflected_kernel = CloneKernelInfo(kernel);
3785       if (reflected_kernel == (KernelInfo *) NULL)
3786         goto error_cleanup;
3787       RotateKernelInfo(reflected_kernel,180);
3788       break;
3789     default:
3790       break;
3791   }
3792
3793   /* Loops around more primitive morpholgy methods
3794   **  erose, dilate, open, close, smooth, edge, etc...
3795   */
3796   /* Loop 1:  iterate the compound method */
3797   method_loop = 0;
3798   method_changed = 1;
3799   while ( method_loop < method_limit && method_changed > 0 ) {
3800     method_loop++;
3801     method_changed = 0;
3802
3803     /* Loop 2:  iterate over each kernel in a multi-kernel list */
3804     norm_kernel = (KernelInfo *) kernel;
3805     this_kernel = (KernelInfo *) kernel;
3806     rflt_kernel = reflected_kernel;
3807
3808     kernel_number = 0;
3809     while ( norm_kernel != NULL ) {
3810
3811       /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3812       stage_loop = 0;          /* the compound morphology stage number */
3813       while ( stage_loop < stage_limit ) {
3814         stage_loop++;   /* The stage of the compound morphology */
3815
3816         /* Select primitive morphology for this stage of compound method */
3817         this_kernel = norm_kernel; /* default use unreflected kernel */
3818         primitive = method;        /* Assume method is a primitive */
3819         switch( method ) {
3820           case ErodeMorphology:      /* just erode */
3821           case EdgeInMorphology:     /* erode and image difference */
3822             primitive = ErodeMorphology;
3823             break;
3824           case DilateMorphology:     /* just dilate */
3825           case EdgeOutMorphology:    /* dilate and image difference */
3826             primitive = DilateMorphology;
3827             break;
3828           case OpenMorphology:       /* erode then dialate */
3829           case TopHatMorphology:     /* open and image difference */
3830             primitive = ErodeMorphology;
3831             if ( stage_loop == 2 )
3832               primitive = DilateMorphology;
3833             break;
3834           case OpenIntensityMorphology:
3835             primitive = ErodeIntensityMorphology;
3836             if ( stage_loop == 2 )
3837               primitive = DilateIntensityMorphology;
3838             break;
3839           case CloseMorphology:      /* dilate, then erode */
3840           case BottomHatMorphology:  /* close and image difference */
3841             this_kernel = rflt_kernel; /* use the reflected kernel */
3842             primitive = DilateMorphology;
3843             if ( stage_loop == 2 )
3844               primitive = ErodeMorphology;
3845             break;
3846           case CloseIntensityMorphology:
3847             this_kernel = rflt_kernel; /* use the reflected kernel */
3848             primitive = DilateIntensityMorphology;
3849             if ( stage_loop == 2 )
3850               primitive = ErodeIntensityMorphology;
3851             break;
3852           case SmoothMorphology:         /* open, close */
3853             switch ( stage_loop ) {
3854               case 1: /* start an open method, which starts with Erode */
3855                 primitive = ErodeMorphology;
3856                 break;
3857               case 2:  /* now Dilate the Erode */
3858                 primitive = DilateMorphology;
3859                 break;
3860               case 3:  /* Reflect kernel a close */
3861                 this_kernel = rflt_kernel; /* use the reflected kernel */
3862                 primitive = DilateMorphology;
3863                 break;
3864               case 4:  /* Finish the Close */
3865                 this_kernel = rflt_kernel; /* use the reflected kernel */
3866                 primitive = ErodeMorphology;
3867                 break;
3868             }
3869             break;
3870           case EdgeMorphology:        /* dilate and erode difference */
3871             primitive = DilateMorphology;
3872             if ( stage_loop == 2 ) {
3873               save_image = curr_image;      /* save the image difference */
3874               curr_image = (Image *) image;
3875               primitive = ErodeMorphology;
3876             }
3877             break;
3878           case CorrelateMorphology:
3879             /* A Correlation is a Convolution with a reflected kernel.
3880             ** However a Convolution is a weighted sum using a reflected
3881             ** kernel.  It may seem stange to convert a Correlation into a
3882             ** Convolution as the Correlation is the simplier method, but
3883             ** Convolution is much more commonly used, and it makes sense to
3884             ** implement it directly so as to avoid the need to duplicate the
3885             ** kernel when it is not required (which is typically the
3886             ** default).
3887             */
3888             this_kernel = rflt_kernel; /* use the reflected kernel */
3889             primitive = ConvolveMorphology;
3890             break;
3891           default:
3892             break;
3893         }
3894         assert( this_kernel != (KernelInfo *) NULL );
3895
3896         /* Extra information for debugging compound operations */
3897         if ( IfMagickTrue(verbose) ) {
3898           if ( stage_limit > 1 )
3899             (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
3900              CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
3901              method_loop,(double) stage_loop);
3902           else if ( primitive != method )
3903             (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
3904               CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
3905               method_loop);
3906           else
3907             v_info[0] = '\0';
3908         }
3909
3910         /* Loop 4: Iterate the kernel with primitive */
3911         kernel_loop = 0;
3912         kernel_changed = 0;
3913         changed = 1;
3914         while ( kernel_loop < kernel_limit && changed > 0 ) {
3915           kernel_loop++;     /* the iteration of this kernel */
3916
3917           /* Create a clone as the destination image, if not yet defined */
3918           if ( work_image == (Image *) NULL )
3919             {
3920               work_image=CloneImage(image,0,0,MagickTrue,exception);
3921               if (work_image == (Image *) NULL)
3922                 goto error_cleanup;
3923               if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
3924                 goto error_cleanup;
3925             }
3926
3927           /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3928           count++;
3929           changed = MorphologyPrimitive(curr_image, work_image, primitive,
3930                        this_kernel, bias, exception);
3931
3932           if ( IfMagickTrue(verbose) ) {
3933             if ( kernel_loop > 1 )
3934               (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
3935             (void) (void) FormatLocaleFile(stderr,
3936               "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3937               v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
3938               primitive),(this_kernel == rflt_kernel ) ? "*" : "",
3939               (double) (method_loop+kernel_loop-1),(double) kernel_number,
3940               (double) count,(double) changed);
3941           }
3942           if ( changed < 0 )
3943             goto error_cleanup;
3944           kernel_changed += changed;
3945           method_changed += changed;
3946
3947           /* prepare next loop */
3948           { Image *tmp = work_image;   /* swap images for iteration */
3949             work_image = curr_image;
3950             curr_image = tmp;
3951           }
3952           if ( work_image == image )
3953             work_image = (Image *) NULL; /* replace input 'image' */
3954
3955         } /* End Loop 4: Iterate the kernel with primitive */
3956
3957         if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed )
3958           (void) FormatLocaleFile(stderr, "   Total %.20g",(double) kernel_changed);
3959         if ( IfMagickTrue(verbose) && stage_loop < stage_limit )
3960           (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
3961
3962 #if 0
3963     (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3964     (void) FormatLocaleFile(stderr, "      curr =0x%lx\n", (unsigned long)curr_image);
3965     (void) FormatLocaleFile(stderr, "      work =0x%lx\n", (unsigned long)work_image);
3966     (void) FormatLocaleFile(stderr, "      save =0x%lx\n", (unsigned long)save_image);
3967     (void) FormatLocaleFile(stderr, "      union=0x%lx\n", (unsigned long)rslt_image);
3968 #endif
3969
3970       } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3971
3972       /*  Final Post-processing for some Compound Methods
3973       **
3974       ** The removal of any 'Sync' channel flag in the Image Compositon
3975       ** below ensures the methematical compose method is applied in a
3976       ** purely mathematical way, and only to the selected channels.
3977       ** Turn off SVG composition 'alpha blending'.
3978       */
3979       switch( method ) {
3980         case EdgeOutMorphology:
3981         case EdgeInMorphology:
3982         case TopHatMorphology:
3983         case BottomHatMorphology:
3984           if ( IfMagickTrue(verbose) )
3985             (void) FormatLocaleFile(stderr,
3986               "\n%s: Difference with original image",CommandOptionToMnemonic(
3987               MagickMorphologyOptions, method) );
3988           (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
3989             MagickTrue,0,0,exception);
3990           break;
3991         case EdgeMorphology:
3992           if ( IfMagickTrue(verbose) )
3993             (void) FormatLocaleFile(stderr,
3994               "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
3995               MagickMorphologyOptions, method) );
3996           (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
3997             MagickTrue,0,0,exception);
3998           save_image = DestroyImage(save_image); /* finished with save image */
3999           break;
4000         default:
4001           break;
4002       }
4003
4004       /* multi-kernel handling:  re-iterate, or compose results */
4005       if ( kernel->next == (KernelInfo *) NULL )
4006         rslt_image = curr_image;   /* just return the resulting image */
4007       else if ( rslt_compose == NoCompositeOp )
4008         { if ( IfMagickTrue(verbose) ) {
4009             if ( this_kernel->next != (KernelInfo *) NULL )
4010               (void) FormatLocaleFile(stderr, " (re-iterate)");
4011             else
4012               (void) FormatLocaleFile(stderr, " (done)");
4013           }
4014           rslt_image = curr_image; /* return result, and re-iterate */
4015         }
4016       else if ( rslt_image == (Image *) NULL)
4017         { if ( IfMagickTrue(verbose) )
4018             (void) FormatLocaleFile(stderr, " (save for compose)");
4019           rslt_image = curr_image;
4020           curr_image = (Image *) image;  /* continue with original image */
4021         }
4022       else
4023         { /* Add the new 'current' result to the composition
4024           **
4025           ** The removal of any 'Sync' channel flag in the Image Compositon
4026           ** below ensures the methematical compose method is applied in a
4027           ** purely mathematical way, and only to the selected channels.
4028           ** IE: Turn off SVG composition 'alpha blending'.
4029           */
4030           if ( IfMagickTrue(verbose) )
4031             (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4032               CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4033           (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4034             0,0,exception);
4035           curr_image = DestroyImage(curr_image);
4036           curr_image = (Image *) image;  /* continue with original image */
4037         }
4038       if ( IfMagickTrue(verbose) )
4039         (void) FormatLocaleFile(stderr, "\n");
4040
4041       /* loop to the next kernel in a multi-kernel list */
4042       norm_kernel = norm_kernel->next;
4043       if ( rflt_kernel != (KernelInfo *) NULL )
4044         rflt_kernel = rflt_kernel->next;
4045       kernel_number++;
4046     } /* End Loop 2: Loop over each kernel */
4047
4048   } /* End Loop 1: compound method interation */
4049
4050   goto exit_cleanup;
4051
4052   /* Yes goto's are bad, but it makes cleanup lot more efficient */
4053 error_cleanup:
4054   if ( curr_image == rslt_image )
4055     curr_image = (Image *) NULL;
4056   if ( rslt_image != (Image *) NULL )
4057     rslt_image = DestroyImage(rslt_image);
4058 exit_cleanup:
4059   if ( curr_image == rslt_image || curr_image == image )
4060     curr_image = (Image *) NULL;
4061   if ( curr_image != (Image *) NULL )
4062     curr_image = DestroyImage(curr_image);
4063   if ( work_image != (Image *) NULL )
4064     work_image = DestroyImage(work_image);
4065   if ( save_image != (Image *) NULL )
4066     save_image = DestroyImage(save_image);
4067   if ( reflected_kernel != (KernelInfo *) NULL )
4068     reflected_kernel = DestroyKernelInfo(reflected_kernel);
4069   return(rslt_image);
4070 }
4071
4072 \f
4073 /*
4074 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4075 %                                                                             %
4076 %                                                                             %
4077 %                                                                             %
4078 %     M o r p h o l o g y I m a g e                                           %
4079 %                                                                             %
4080 %                                                                             %
4081 %                                                                             %
4082 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4083 %
4084 %  MorphologyImage() applies a user supplied kernel to the image according to
4085 %  the given mophology method.
4086 %
4087 %  This function applies any and all user defined settings before calling
4088 %  the above internal function MorphologyApply().
4089 %
4090 %  User defined settings include...
4091 %    * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4092 %    * Kernel Scale/normalize settings             ("-define convolve:scale=??")
4093 %      This can also includes the addition of a scaled unity kernel.
4094 %    * Show Kernel being applied                   ("-define showkernel=1")
4095 %
4096 %  Other operators that do not want user supplied options interfering,
4097 %  especially "convolve:bias" and "showkernel" should use MorphologyApply()
4098 %  directly.
4099 %
4100 %  The format of the MorphologyImage method is:
4101 %
4102 %      Image *MorphologyImage(const Image *image,MorphologyMethod method,
4103 %        const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4104 %
4105 %  A description of each parameter follows:
4106 %
4107 %    o image: the image.
4108 %
4109 %    o method: the morphology method to be applied.
4110 %
4111 %    o iterations: apply the operation this many times (or no change).
4112 %                  A value of -1 means loop until no change found.
4113 %                  How this is applied may depend on the morphology method.
4114 %                  Typically this is a value of 1.
4115 %
4116 %    o kernel: An array of double representing the morphology kernel.
4117 %              Warning: kernel may be normalized for the Convolve method.
4118 %
4119 %    o exception: return any errors or warnings in this structure.
4120 %
4121 */
4122 MagickExport Image *MorphologyImage(const Image *image,
4123   const MorphologyMethod method,const ssize_t iterations,
4124   const KernelInfo *kernel,ExceptionInfo *exception)
4125 {
4126   KernelInfo
4127     *curr_kernel;
4128
4129   CompositeOperator
4130     compose;
4131
4132   Image
4133     *morphology_image;
4134
4135   double
4136     bias;
4137
4138   curr_kernel = (KernelInfo *) kernel;
4139   bias=0.0;
4140   compose = UndefinedCompositeOp;  /* use default for method */
4141
4142   /* Apply Convolve/Correlate Normalization and Scaling Factors.
4143    * This is done BEFORE the ShowKernelInfo() function is called so that
4144    * users can see the results of the 'option:convolve:scale' option.
4145    */
4146   if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4147       const char
4148         *artifact;
4149
4150       /* Get the bias value as it will be needed */
4151       artifact = GetImageArtifact(image,"convolve:bias");
4152       if ( artifact != (const char *) NULL) {
4153         if (IfMagickFalse(IsGeometry(artifact)))
4154           (void) ThrowMagickException(exception,GetMagickModule(),
4155                OptionWarning,"InvalidSetting","'%s' '%s'",
4156                "convolve:bias",artifact);
4157         else
4158           bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4159       }
4160
4161       /* Scale kernel according to user wishes */
4162       artifact = GetImageArtifact(image,"convolve:scale");
4163       if ( artifact != (const char *)NULL ) {
4164         if (IfMagickFalse(IsGeometry(artifact)))
4165           (void) ThrowMagickException(exception,GetMagickModule(),
4166                OptionWarning,"InvalidSetting","'%s' '%s'",
4167                "convolve:scale",artifact);
4168         else {
4169           if ( curr_kernel == kernel )
4170             curr_kernel = CloneKernelInfo(kernel);
4171           if (curr_kernel == (KernelInfo *) NULL)
4172             return((Image *) NULL);
4173           ScaleGeometryKernelInfo(curr_kernel, artifact);
4174         }
4175       }
4176     }
4177
4178   /* display the (normalized) kernel via stderr */
4179   if ( IfStringTrue(GetImageArtifact(image,"showkernel"))
4180     || IfStringTrue(GetImageArtifact(image,"convolve:showkernel"))
4181     || IfStringTrue(GetImageArtifact(image,"morphology:showkernel")) )
4182     ShowKernelInfo(curr_kernel);
4183
4184   /* Override the default handling of multi-kernel morphology results
4185    * If 'Undefined' use the default method
4186    * If 'None' (default for 'Convolve') re-iterate previous result
4187    * Otherwise merge resulting images using compose method given.
4188    * Default for 'HitAndMiss' is 'Lighten'.
4189    */
4190   { const char
4191       *artifact;
4192     ssize_t
4193       parse;
4194
4195     artifact = GetImageArtifact(image,"morphology:compose");
4196     if ( artifact != (const char *) NULL) {
4197       parse=ParseCommandOption(MagickComposeOptions,
4198         MagickFalse,artifact);
4199       if ( parse < 0 )
4200         (void) ThrowMagickException(exception,GetMagickModule(),
4201              OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4202              "morphology:compose",artifact);
4203       else
4204         compose=(CompositeOperator)parse;
4205     }
4206   }
4207   /* Apply the Morphology */
4208   morphology_image = MorphologyApply(image,method,iterations,
4209     curr_kernel,compose,bias,exception);
4210
4211   /* Cleanup and Exit */
4212   if ( curr_kernel != kernel )
4213     curr_kernel=DestroyKernelInfo(curr_kernel);
4214   return(morphology_image);
4215 }
4216 \f
4217 /*
4218 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4219 %                                                                             %
4220 %                                                                             %
4221 %                                                                             %
4222 +     R o t a t e K e r n e l I n f o                                         %
4223 %                                                                             %
4224 %                                                                             %
4225 %                                                                             %
4226 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4227 %
4228 %  RotateKernelInfo() rotates the kernel by the angle given.
4229 %
4230 %  Currently it is restricted to 90 degree angles, of either 1D kernels
4231 %  or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4232 %  It will ignore usless rotations for specific 'named' built-in kernels.
4233 %
4234 %  The format of the RotateKernelInfo method is:
4235 %
4236 %      void RotateKernelInfo(KernelInfo *kernel, double angle)
4237 %
4238 %  A description of each parameter follows:
4239 %
4240 %    o kernel: the Morphology/Convolution kernel
4241 %
4242 %    o angle: angle to rotate in degrees
4243 %
4244 % This function is currently internal to this module only, but can be exported
4245 % to other modules if needed.
4246 */
4247 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4248 {
4249   /* angle the lower kernels first */
4250   if ( kernel->next != (KernelInfo *) NULL)
4251     RotateKernelInfo(kernel->next, angle);
4252
4253   /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4254   **
4255   ** TODO: expand beyond simple 90 degree rotates, flips and flops
4256   */
4257
4258   /* Modulus the angle */
4259   angle = fmod(angle, 360.0);
4260   if ( angle < 0 )
4261     angle += 360.0;
4262
4263   if ( 337.5 < angle || angle <= 22.5 )
4264     return;   /* Near zero angle - no change! - At least not at this time */
4265
4266   /* Handle special cases */
4267   switch (kernel->type) {
4268     /* These built-in kernels are cylindrical kernels, rotating is useless */
4269     case GaussianKernel:
4270     case DoGKernel:
4271     case LoGKernel:
4272     case DiskKernel:
4273     case PeaksKernel:
4274     case LaplacianKernel:
4275     case ChebyshevKernel:
4276     case ManhattanKernel:
4277     case EuclideanKernel:
4278       return;
4279
4280     /* These may be rotatable at non-90 angles in the future */
4281     /* but simply rotating them in multiples of 90 degrees is useless */
4282     case SquareKernel:
4283     case DiamondKernel:
4284     case PlusKernel:
4285     case CrossKernel:
4286       return;
4287
4288     /* These only allows a +/-90 degree rotation (by transpose) */
4289     /* A 180 degree rotation is useless */
4290     case BlurKernel:
4291       if ( 135.0 < angle && angle <= 225.0 )
4292         return;
4293       if ( 225.0 < angle && angle <= 315.0 )
4294         angle -= 180;
4295       break;
4296
4297     default:
4298       break;
4299   }
4300   /* Attempt rotations by 45 degrees  -- 3x3 kernels only */
4301   if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4302     {
4303       if ( kernel->width == 3 && kernel->height == 3 )
4304         { /* Rotate a 3x3 square by 45 degree angle */
4305           double t  = kernel->values[0];
4306           kernel->values[0] = kernel->values[3];
4307           kernel->values[3] = kernel->values[6];
4308           kernel->values[6] = kernel->values[7];
4309           kernel->values[7] = kernel->values[8];
4310           kernel->values[8] = kernel->values[5];
4311           kernel->values[5] = kernel->values[2];
4312           kernel->values[2] = kernel->values[1];
4313           kernel->values[1] = t;
4314           /* rotate non-centered origin */
4315           if ( kernel->x != 1 || kernel->y != 1 ) {
4316             ssize_t x,y;
4317             x = (ssize_t) kernel->x-1;
4318             y = (ssize_t) kernel->y-1;
4319                  if ( x == y  ) x = 0;
4320             else if ( x == 0  ) x = -y;
4321             else if ( x == -y ) y = 0;
4322             else if ( y == 0  ) y = x;
4323             kernel->x = (ssize_t) x+1;
4324             kernel->y = (ssize_t) y+1;
4325           }
4326           angle = fmod(angle+315.0, 360.0);  /* angle reduced 45 degrees */
4327           kernel->angle = fmod(kernel->angle+45.0, 360.0);
4328         }
4329       else
4330         perror("Unable to rotate non-3x3 kernel by 45 degrees");
4331     }
4332   if ( 45.0 < fmod(angle, 180.0)  && fmod(angle,180.0) <= 135.0 )
4333     {
4334       if ( kernel->width == 1 || kernel->height == 1 )
4335         { /* Do a transpose of a 1 dimensional kernel,
4336           ** which results in a fast 90 degree rotation of some type.
4337           */
4338           ssize_t
4339             t;
4340           t = (ssize_t) kernel->width;
4341           kernel->width = kernel->height;
4342           kernel->height = (size_t) t;
4343           t = kernel->x;
4344           kernel->x = kernel->y;
4345           kernel->y = t;
4346           if ( kernel->width == 1 ) {
4347             angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
4348             kernel->angle = fmod(kernel->angle+90.0, 360.0);
4349           } else {
4350             angle = fmod(angle+90.0, 360.0);   /* angle increased 90 degrees */
4351             kernel->angle = fmod(kernel->angle+270.0, 360.0);
4352           }
4353         }
4354       else if ( kernel->width == kernel->height )
4355         { /* Rotate a square array of values by 90 degrees */
4356           { register ssize_t
4357               i,j,x,y;
4358
4359             register MagickRealType
4360               *k,t;
4361
4362             k=kernel->values;
4363             for( i=0, x=(ssize_t) kernel->width-1;  i<=x;   i++, x--)
4364               for( j=0, y=(ssize_t) kernel->height-1;  j<y;   j++, y--)
4365                 { t                    = k[i+j*kernel->width];
4366                   k[i+j*kernel->width] = k[j+x*kernel->width];
4367                   k[j+x*kernel->width] = k[x+y*kernel->width];
4368                   k[x+y*kernel->width] = k[y+i*kernel->width];
4369                   k[y+i*kernel->width] = t;
4370                 }
4371           }
4372           /* rotate the origin - relative to center of array */
4373           { register ssize_t x,y;
4374             x = (ssize_t) (kernel->x*2-kernel->width+1);
4375             y = (ssize_t) (kernel->y*2-kernel->height+1);
4376             kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4377             kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4378           }
4379           angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
4380           kernel->angle = fmod(kernel->angle+90.0, 360.0);
4381         }
4382       else
4383         perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4384     }
4385   if ( 135.0 < angle && angle <= 225.0 )
4386     {
4387       /* For a 180 degree rotation - also know as a reflection
4388        * This is actually a very very common operation!
4389        * Basically all that is needed is a reversal of the kernel data!
4390        * And a reflection of the origon
4391        */
4392       MagickRealType
4393         t;
4394
4395       register MagickRealType
4396         *k;
4397
4398       ssize_t
4399         i,
4400         j;
4401
4402       k=kernel->values;
4403       j=(ssize_t) (kernel->width*kernel->height-1);
4404       for (i=0;  i < j;  i++, j--)
4405         t=k[i],  k[i]=k[j],  k[j]=t;
4406
4407       kernel->x = (ssize_t) kernel->width  - kernel->x - 1;
4408       kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4409       angle = fmod(angle-180.0, 360.0);   /* angle+180 degrees */
4410       kernel->angle = fmod(kernel->angle+180.0, 360.0);
4411     }
4412   /* At this point angle should at least between -45 (315) and +45 degrees
4413    * In the future some form of non-orthogonal angled rotates could be
4414    * performed here, posibily with a linear kernel restriction.
4415    */
4416
4417   return;
4418 }
4419 \f
4420 /*
4421 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4422 %                                                                             %
4423 %                                                                             %
4424 %                                                                             %
4425 %     S c a l e G e o m e t r y K e r n e l I n f o                           %
4426 %                                                                             %
4427 %                                                                             %
4428 %                                                                             %
4429 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4430 %
4431 %  ScaleGeometryKernelInfo() takes a geometry argument string, typically
4432 %  provided as a  "-set option:convolve:scale {geometry}" user setting,
4433 %  and modifies the kernel according to the parsed arguments of that setting.
4434 %
4435 %  The first argument (and any normalization flags) are passed to
4436 %  ScaleKernelInfo() to scale/normalize the kernel.  The second argument
4437 %  is then passed to UnityAddKernelInfo() to add a scled unity kernel
4438 %  into the scaled/normalized kernel.
4439 %
4440 %  The format of the ScaleGeometryKernelInfo method is:
4441 %
4442 %      void ScaleGeometryKernelInfo(KernelInfo *kernel,
4443 %        const double scaling_factor,const MagickStatusType normalize_flags)
4444 %
4445 %  A description of each parameter follows:
4446 %
4447 %    o kernel: the Morphology/Convolution kernel to modify
4448 %
4449 %    o geometry:
4450 %             The geometry string to parse, typically from the user provided
4451 %             "-set option:convolve:scale {geometry}" setting.
4452 %
4453 */
4454 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4455   const char *geometry)
4456 {
4457   MagickStatusType
4458     flags;
4459
4460   GeometryInfo
4461     args;
4462
4463   SetGeometryInfo(&args);
4464   flags = ParseGeometry(geometry, &args);
4465
4466 #if 0
4467   /* For Debugging Geometry Input */
4468   (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4469        flags, args.rho, args.sigma, args.xi, args.psi );
4470 #endif
4471
4472   if ( (flags & PercentValue) != 0 )      /* Handle Percentage flag*/
4473     args.rho *= 0.01,  args.sigma *= 0.01;
4474
4475   if ( (flags & RhoValue) == 0 )          /* Set Defaults for missing args */
4476     args.rho = 1.0;
4477   if ( (flags & SigmaValue) == 0 )
4478     args.sigma = 0.0;
4479
4480   /* Scale/Normalize the input kernel */
4481   ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4482
4483   /* Add Unity Kernel, for blending with original */
4484   if ( (flags & SigmaValue) != 0 )
4485     UnityAddKernelInfo(kernel, args.sigma);
4486
4487   return;
4488 }
4489 /*
4490 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4491 %                                                                             %
4492 %                                                                             %
4493 %                                                                             %
4494 %     S c a l e K e r n e l I n f o                                           %
4495 %                                                                             %
4496 %                                                                             %
4497 %                                                                             %
4498 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4499 %
4500 %  ScaleKernelInfo() scales the given kernel list by the given amount, with or
4501 %  without normalization of the sum of the kernel values (as per given flags).
4502 %
4503 %  By default (no flags given) the values within the kernel is scaled
4504 %  directly using given scaling factor without change.
4505 %
4506 %  If either of the two 'normalize_flags' are given the kernel will first be
4507 %  normalized and then further scaled by the scaling factor value given.
4508 %
4509 %  Kernel normalization ('normalize_flags' given) is designed to ensure that
4510 %  any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4511 %  morphology methods will fall into -1.0 to +1.0 range.  Note that for
4512 %  non-HDRI versions of IM this may cause images to have any negative results
4513 %  clipped, unless some 'bias' is used.
4514 %
4515 %  More specifically.  Kernels which only contain positive values (such as a
4516 %  'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4517 %  ensuring a 0.0 to +1.0 output range for non-HDRI images.
4518 %
4519 %  For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4520 %  the kernel will be scaled by the absolute of the sum of kernel values, so
4521 %  that it will generally fall within the +/- 1.0 range.
4522 %
4523 %  For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4524 %  will be scaled by just the sum of the postive values, so that its output
4525 %  range will again fall into the  +/- 1.0 range.
4526 %
4527 %  For special kernels designed for locating shapes using 'Correlate', (often
4528 %  only containing +1 and -1 values, representing foreground/brackground
4529 %  matching) a special normalization method is provided to scale the positive
4530 %  values separately to those of the negative values, so the kernel will be
4531 %  forced to become a zero-sum kernel better suited to such searches.
4532 %
4533 %  WARNING: Correct normalization of the kernel assumes that the '*_range'
4534 %  attributes within the kernel structure have been correctly set during the
4535 %  kernels creation.
4536 %
4537 %  NOTE: The values used for 'normalize_flags' have been selected specifically
4538 %  to match the use of geometry options, so that '!' means NormalizeValue, '^'
4539 %  means CorrelateNormalizeValue.  All other GeometryFlags values are ignored.
4540 %
4541 %  The format of the ScaleKernelInfo method is:
4542 %
4543 %      void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4544 %               const MagickStatusType normalize_flags )
4545 %
4546 %  A description of each parameter follows:
4547 %
4548 %    o kernel: the Morphology/Convolution kernel
4549 %
4550 %    o scaling_factor:
4551 %             multiply all values (after normalization) by this factor if not
4552 %             zero.  If the kernel is normalized regardless of any flags.
4553 %
4554 %    o normalize_flags:
4555 %             GeometryFlags defining normalization method to use.
4556 %             specifically: NormalizeValue, CorrelateNormalizeValue,
4557 %                           and/or PercentValue
4558 %
4559 */
4560 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4561   const double scaling_factor,const GeometryFlags normalize_flags)
4562 {
4563   register double
4564     pos_scale,
4565     neg_scale;
4566
4567   register ssize_t
4568     i;
4569
4570   /* do the other kernels in a multi-kernel list first */
4571   if ( kernel->next != (KernelInfo *) NULL)
4572     ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4573
4574   /* Normalization of Kernel */
4575   pos_scale = 1.0;
4576   if ( (normalize_flags&NormalizeValue) != 0 ) {
4577     if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4578       /* non-zero-summing kernel (generally positive) */
4579       pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4580     else
4581       /* zero-summing kernel */
4582       pos_scale = kernel->positive_range;
4583   }
4584   /* Force kernel into a normalized zero-summing kernel */
4585   if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4586     pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4587                  ? kernel->positive_range : 1.0;
4588     neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4589                  ? -kernel->negative_range : 1.0;
4590   }
4591   else
4592     neg_scale = pos_scale;
4593
4594   /* finialize scaling_factor for positive and negative components */
4595   pos_scale = scaling_factor/pos_scale;
4596   neg_scale = scaling_factor/neg_scale;
4597
4598   for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4599     if ( ! IfNaN(kernel->values[i]) )
4600       kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4601
4602   /* convolution output range */
4603   kernel->positive_range *= pos_scale;
4604   kernel->negative_range *= neg_scale;
4605   /* maximum and minimum values in kernel */
4606   kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4607   kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4608
4609   /* swap kernel settings if user's scaling factor is negative */
4610   if ( scaling_factor < MagickEpsilon ) {
4611     double t;
4612     t = kernel->positive_range;
4613     kernel->positive_range = kernel->negative_range;
4614     kernel->negative_range = t;
4615     t = kernel->maximum;
4616     kernel->maximum = kernel->minimum;
4617     kernel->minimum = 1;
4618   }
4619
4620   return;
4621 }
4622 \f
4623 /*
4624 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4625 %                                                                             %
4626 %                                                                             %
4627 %                                                                             %
4628 %     S h o w K e r n e l I n f o                                             %
4629 %                                                                             %
4630 %                                                                             %
4631 %                                                                             %
4632 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4633 %
4634 %  ShowKernelInfo() outputs the details of the given kernel defination to
4635 %  standard error, generally due to a users 'showkernel' option request.
4636 %
4637 %  The format of the ShowKernel method is:
4638 %
4639 %      void ShowKernelInfo(const KernelInfo *kernel)
4640 %
4641 %  A description of each parameter follows:
4642 %
4643 %    o kernel: the Morphology/Convolution kernel
4644 %
4645 */
4646 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4647 {
4648   const KernelInfo
4649     *k;
4650
4651   size_t
4652     c, i, u, v;
4653
4654   for (c=0, k=kernel;  k != (KernelInfo *) NULL;  c++, k=k->next ) {
4655
4656     (void) FormatLocaleFile(stderr, "Kernel");
4657     if ( kernel->next != (KernelInfo *) NULL )
4658       (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4659     (void) FormatLocaleFile(stderr, " \"%s",
4660           CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4661     if ( fabs(k->angle) >= MagickEpsilon )
4662       (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4663     (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4664       k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4665     (void) FormatLocaleFile(stderr,
4666           " with values from %.*lg to %.*lg\n",
4667           GetMagickPrecision(), k->minimum,
4668           GetMagickPrecision(), k->maximum);
4669     (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4670           GetMagickPrecision(), k->negative_range,
4671           GetMagickPrecision(), k->positive_range);
4672     if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4673       (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4674     else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4675       (void) FormatLocaleFile(stderr, " (Normalized)\n");
4676     else
4677       (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4678           GetMagickPrecision(), k->positive_range+k->negative_range);
4679     for (i=v=0; v < k->height; v++) {
4680       (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4681       for (u=0; u < k->width; u++, i++)
4682         if ( IfNaN(k->values[i]) )
4683           (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4684         else
4685           (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4686               GetMagickPrecision(), (double) k->values[i]);
4687       (void) FormatLocaleFile(stderr,"\n");
4688     }
4689   }
4690 }
4691 \f
4692 /*
4693 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4694 %                                                                             %
4695 %                                                                             %
4696 %                                                                             %
4697 %     U n i t y A d d K e r n a l I n f o                                     %
4698 %                                                                             %
4699 %                                                                             %
4700 %                                                                             %
4701 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4702 %
4703 %  UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4704 %  to the given pre-scaled and normalized Kernel.  This in effect adds that
4705 %  amount of the original image into the resulting convolution kernel.  This
4706 %  value is usually provided by the user as a percentage value in the
4707 %  'convolve:scale' setting.
4708 %
4709 %  The resulting effect is to convert the defined kernels into blended
4710 %  soft-blurs, unsharp kernels or into sharpening kernels.
4711 %
4712 %  The format of the UnityAdditionKernelInfo method is:
4713 %
4714 %      void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4715 %
4716 %  A description of each parameter follows:
4717 %
4718 %    o kernel: the Morphology/Convolution kernel
4719 %
4720 %    o scale:
4721 %             scaling factor for the unity kernel to be added to
4722 %             the given kernel.
4723 %
4724 */
4725 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4726   const double scale)
4727 {
4728   /* do the other kernels in a multi-kernel list first */
4729   if ( kernel->next != (KernelInfo *) NULL)
4730     UnityAddKernelInfo(kernel->next, scale);
4731
4732   /* Add the scaled unity kernel to the existing kernel */
4733   kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4734   CalcKernelMetaData(kernel);  /* recalculate the meta-data */
4735
4736   return;
4737 }
4738 \f
4739 /*
4740 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4741 %                                                                             %
4742 %                                                                             %
4743 %                                                                             %
4744 %     Z e r o K e r n e l N a n s                                             %
4745 %                                                                             %
4746 %                                                                             %
4747 %                                                                             %
4748 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4749 %
4750 %  ZeroKernelNans() replaces any special 'nan' value that may be present in
4751 %  the kernel with a zero value.  This is typically done when the kernel will
4752 %  be used in special hardware (GPU) convolution processors, to simply
4753 %  matters.
4754 %
4755 %  The format of the ZeroKernelNans method is:
4756 %
4757 %      void ZeroKernelNans (KernelInfo *kernel)
4758 %
4759 %  A description of each parameter follows:
4760 %
4761 %    o kernel: the Morphology/Convolution kernel
4762 %
4763 */
4764 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4765 {
4766   register size_t
4767     i;
4768
4769   /* do the other kernels in a multi-kernel list first */
4770   if ( kernel->next != (KernelInfo *) NULL)
4771     ZeroKernelNans(kernel->next);
4772
4773   for (i=0; i < (kernel->width*kernel->height); i++)
4774     if ( IfNaN(kernel->values[i]) )
4775       kernel->values[i] = 0.0;
4776
4777   return;
4778 }