granicus.if.org Git - imagemagick/blob - MagickCore/morphology.c

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