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

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