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

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