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