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