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