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