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-2011 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/hashmap.h"
  61 #include "MagickCore/image.h"
  62 #include "MagickCore/image-private.h"
  63 #include "MagickCore/list.h"
  64 #include "MagickCore/magick.h"
  65 #include "MagickCore/memory_.h"
  66 #include "MagickCore/monitor-private.h"
  67 #include "MagickCore/morphology.h"
  68 #include "MagickCore/morphology-private.h"
  69 #include "MagickCore/option.h"
  70 #include "MagickCore/pixel-accessor.h"
  71 #include "MagickCore/prepress.h"
  72 #include "MagickCore/quantize.h"
  73 #include "MagickCore/registry.h"
  74 #include "MagickCore/semaphore.h"
  75 #include "MagickCore/splay-tree.h"
  76 #include "MagickCore/statistic.h"
  77 #include "MagickCore/string_.h"
  78 #include "MagickCore/string-private.h"
  79 #include "MagickCore/token.h"
  80 #include "MagickCore/utility.h"
  81 \f
  82
  83 /*
  84 ** The following test is for special floating point numbers of value NaN (not
  85 ** a number), that may be used within a Kernel Definition.  NaN's are defined
  86 ** as part of the IEEE standard for floating point number representation.
  87 **
  88 ** These are used as a Kernel value to mean that this kernel position is not
  89 ** part of the kernel neighbourhood for convolution or morphology processing,
  90 ** and thus should be ignored.  This allows the use of 'shaped' kernels.
  91 **
  92 ** The special properity that two NaN's are never equal, even if they are from
  93 ** the same variable allow you to test if a value is special NaN value.
  94 **
  95 ** This macro  IsNaN() is thus is only true if the value given is NaN.
  96 */
  97 #define IsNan(a)   ((a)!=(a))
  98
  99 /*
 100   Other global definitions used by module.
 101 */
 102 static inline double MagickMin(const double x,const double y)
 103 {
 104   return( x < y ? x : y);
 105 }
 106 static inline double MagickMax(const double x,const double y)
 107 {
 108   return( x > y ? x : y);
 109 }
 110 #define Minimize(assign,value) assign=MagickMin(assign,value)
 111 #define Maximize(assign,value) assign=MagickMax(assign,value)
 112
 113 /* Currently these are only internal to this module */
 114 static void
 115   CalcKernelMetaData(KernelInfo *),
 116   ExpandMirrorKernelInfo(KernelInfo *),
 117   ExpandRotateKernelInfo(KernelInfo *, const double),
 118   RotateKernelInfo(KernelInfo *, double);
 119 \f
 120
 121 /* Quick function to find last kernel in a kernel list */
 122 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
 123 {
 124   while (kernel->next != (KernelInfo *) NULL)
 125     kernel = kernel->next;
 126   return(kernel);
 127 }
 128
 129 /*
 130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 131 %                                                                             %
 132 %                                                                             %
 133 %                                                                             %
 134 %     A c q u i r e K e r n e l I n f o                                       %
 135 %                                                                             %
 136 %                                                                             %
 137 %                                                                             %
 138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 139 %
 140 %  AcquireKernelInfo() takes the given string (generally supplied by the
 141 %  user) and converts it into a Morphology/Convolution Kernel.  This allows
 142 %  users to specify a kernel from a number of pre-defined kernels, or to fully
 143 %  specify their own kernel for a specific Convolution or Morphology
 144 %  Operation.
 145 %
 146 %  The kernel so generated can be any rectangular array of floating point
 147 %  values (doubles) with the 'control point' or 'pixel being affected'
 148 %  anywhere within that array of values.
 149 %
 150 %  Previously IM was restricted to a square of odd size using the exact
 151 %  center as origin, this is no longer the case, and any rectangular kernel
 152 %  with any value being declared the origin. This in turn allows the use of
 153 %  highly asymmetrical kernels.
 154 %
 155 %  The floating point values in the kernel can also include a special value
 156 %  known as 'nan' or 'not a number' to indicate that this value is not part
 157 %  of the kernel array. This allows you to shaped the kernel within its
 158 %  rectangular area. That is 'nan' values provide a 'mask' for the kernel
 159 %  shape.  However at least one non-nan value must be provided for correct
 160 %  working of a kernel.
 161 %
 162 %  The returned kernel should be freed using the DestroyKernelInfo() when you
 163 %  are finished with it.  Do not free this memory yourself.
 164 %
 165 %  Input kernel defintion strings can consist of any of three types.
 166 %
 167 %    "name:args[[@><]"
 168 %         Select from one of the built in kernels, using the name and
 169 %         geometry arguments supplied.  See AcquireKernelBuiltIn()
 170 %
 171 %    "WxH[+X+Y][@><]:num, num, num ..."
 172 %         a kernel of size W by H, with W*H floating point numbers following.
 173 %         the 'center' can be optionally be defined at +X+Y (such that +0+0
 174 %         is top left corner). If not defined the pixel in the center, for
 175 %         odd sizes, or to the immediate top or left of center for even sizes
 176 %         is automatically selected.
 177 %
 178 %    "num, num, num, num, ..."
 179 %         list of floating point numbers defining an 'old style' odd sized
 180 %         square kernel.  At least 9 values should be provided for a 3x3
 181 %         square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
 182 %         Values can be space or comma separated.  This is not recommended.
 183 %
 184 %  You can define a 'list of kernels' which can be used by some morphology
 185 %  operators A list is defined as a semi-colon separated list kernels.
 186 %
 187 %     " kernel ; kernel ; kernel ; "
 188 %
 189 %  Any extra ';' characters, at start, end or between kernel defintions are
 190 %  simply ignored.
 191 %
 192 %  The special flags will expand a single kernel, into a list of rotated
 193 %  kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
 194 %  cyclic rotations, while a '>' will generate a list of 90-degree rotations.
 195 %  The '<' also exands using 90-degree rotates, but giving a 180-degree
 196 %  reflected kernel before the +/- 90-degree rotations, which can be important
 197 %  for Thinning operations.
 198 %
 199 %  Note that 'name' kernels will start with an alphabetic character while the
 200 %  new kernel specification has a ':' character in its specification string.
 201 %  If neither is the case, it is assumed an old style of a simple list of
 202 %  numbers generating a odd-sized square kernel has been given.
 203 %
 204 %  The format of the AcquireKernal method is:
 205 %
 206 %      KernelInfo *AcquireKernelInfo(const char *kernel_string)
 207 %
 208 %  A description of each parameter follows:
 209 %
 210 %    o kernel_string: the Morphology/Convolution kernel wanted.
 211 %
 212 */
 213
 214 /* This was separated so that it could be used as a separate
 215 ** array input handling function, such as for -color-matrix
 216 */
 217 static KernelInfo *ParseKernelArray(const char *kernel_string)
 218 {
 219   KernelInfo
 220     *kernel;
 221
 222   char
 223     token[MaxTextExtent];
 224
 225   const char
 226     *p,
 227     *end;
 228
 229   register ssize_t
 230     i;
 231
 232   double
 233     nan = sqrt((double)-1.0);  /* Special Value : Not A Number */
 234
 235   MagickStatusType
 236     flags;
 237
 238   GeometryInfo
 239     args;
 240
 241   kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
 242   if (kernel == (KernelInfo *)NULL)
 243     return(kernel);
 244   (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
 245   kernel->minimum = kernel->maximum = kernel->angle = 0.0;
 246   kernel->negative_range = kernel->positive_range = 0.0;
 247   kernel->type = UserDefinedKernel;
 248   kernel->next = (KernelInfo *) NULL;
 249   kernel->signature = MagickSignature;
 250   if (kernel_string == (const char *) NULL)
 251     return(kernel);
 252
 253   /* find end of this specific kernel definition string */
 254   end = strchr(kernel_string, ';');
 255   if ( end == (char *) NULL )
 256     end = strchr(kernel_string, '\0');
 257
 258   /* clear flags - for Expanding kernel lists thorugh rotations */
 259    flags = NoValue;
 260
 261   /* Has a ':' in argument - New user kernel specification */
 262   p = strchr(kernel_string, ':');
 263   if ( p != (char *) NULL && p < end)
 264     {
 265       /* ParseGeometry() needs the geometry separated! -- Arrgghh */
 266       memcpy(token, kernel_string, (size_t) (p-kernel_string));
 267       token[p-kernel_string] = '\0';
 268       SetGeometryInfo(&args);
 269       flags = ParseGeometry(token, &args);
 270
 271       /* Size handling and checks of geometry settings */
 272       if ( (flags & WidthValue) == 0 ) /* if no width then */
 273         args.rho = args.sigma;         /* then  width = height */
 274       if ( args.rho < 1.0 )            /* if width too small */
 275          args.rho = 1.0;               /* then  width = 1 */
 276       if ( args.sigma < 1.0 )          /* if height too small */
 277         args.sigma = args.rho;         /* then  height = width */
 278       kernel->width = (size_t)args.rho;
 279       kernel->height = (size_t)args.sigma;
 280
 281       /* Offset Handling and Checks */
 282       if ( args.xi  < 0.0 || args.psi < 0.0 )
 283         return(DestroyKernelInfo(kernel));
 284       kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
 285                                                : (ssize_t) (kernel->width-1)/2;
 286       kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
 287                                                : (ssize_t) (kernel->height-1)/2;
 288       if ( kernel->x >= (ssize_t) kernel->width ||
 289            kernel->y >= (ssize_t) kernel->height )
 290         return(DestroyKernelInfo(kernel));
 291
 292       p++; /* advance beyond the ':' */
 293     }
 294   else
 295     { /* ELSE - Old old specification, forming odd-square kernel */
 296       /* count up number of values given */
 297       p=(const char *) kernel_string;
 298       while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
 299         p++;  /* ignore "'" chars for convolve filter usage - Cristy */
 300       for (i=0; p < end; i++)
 301       {
 302         GetMagickToken(p,&p,token);
 303         if (*token == ',')
 304           GetMagickToken(p,&p,token);
 305       }
 306       /* set the size of the kernel - old sized square */
 307       kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
 308       kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
 309       p=(const char *) kernel_string;
 310       while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
 311         p++;  /* ignore "'" chars for convolve filter usage - Cristy */
 312     }
 313
 314   /* Read in the kernel values from rest of input string argument */
 315   kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 316                         kernel->height*sizeof(double));
 317   if (kernel->values == (double *) NULL)
 318     return(DestroyKernelInfo(kernel));
 319
 320   kernel->minimum = +MagickHuge;
 321   kernel->maximum = -MagickHuge;
 322   kernel->negative_range = kernel->positive_range = 0.0;
 323
 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] = InterpretLocaleValue(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=(double *) AcquireQuantumMemory(1,sizeof(double));
1032         if (kernel->values == (double *) NULL)
1033           return(DestroyKernelInfo(kernel));
1034         kernel->maximum = kernel->values[0] = args->rho;
1035         break;
1036       }
1037       break;
1038     case GaussianKernel:
1039     case DoGKernel:
1040     case LoGKernel:
1041       { double
1042           sigma = fabs(args->sigma),
1043           sigma2 = fabs(args->xi),
1044           A, B, R;
1045
1046         if ( args->rho >= 1.0 )
1047           kernel->width = (size_t)args->rho*2+1;
1048         else if ( (type != DoGKernel) || (sigma >= sigma2) )
1049           kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1050         else
1051           kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1052         kernel->height = kernel->width;
1053         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1054         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1055                               kernel->height*sizeof(double));
1056         if (kernel->values == (double *) NULL)
1057           return(DestroyKernelInfo(kernel));
1058
1059         /* WARNING: The following generates a 'sampled gaussian' kernel.
1060          * What we really want is a 'discrete gaussian' kernel.
1061          *
1062          * How to do this is I don't know, but appears to be basied on the
1063          * Error Function 'erf()' (intergral of a gaussian)
1064          */
1065
1066         if ( type == GaussianKernel || type == DoGKernel )
1067           { /* Calculate a Gaussian,  OR positive half of a DoG */
1068             if ( sigma > MagickEpsilon )
1069               { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
1070                 B = (double) (1.0/(Magick2PI*sigma*sigma));
1071                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1072                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1073                       kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1074               }
1075             else /* limiting case - a unity (normalized Dirac) kernel */
1076               { (void) ResetMagickMemory(kernel->values,0, (size_t)
1077                             kernel->width*kernel->height*sizeof(double));
1078                 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1079               }
1080           }
1081
1082         if ( type == DoGKernel )
1083           { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1084             if ( sigma2 > MagickEpsilon )
1085               { sigma = sigma2;                /* simplify loop expressions */
1086                 A = 1.0/(2.0*sigma*sigma);
1087                 B = (double) (1.0/(Magick2PI*sigma*sigma));
1088                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1089                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1090                     kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1091               }
1092             else /* limiting case - a unity (normalized Dirac) kernel */
1093               kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1094           }
1095
1096         if ( type == LoGKernel )
1097           { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1098             if ( sigma > MagickEpsilon )
1099               { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
1100                 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1101                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1102                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1103                     { R = ((double)(u*u+v*v))*A;
1104                       kernel->values[i] = (1-R)*exp(-R)*B;
1105                     }
1106               }
1107             else /* special case - generate a unity kernel */
1108               { (void) ResetMagickMemory(kernel->values,0, (size_t)
1109                             kernel->width*kernel->height*sizeof(double));
1110                 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1111               }
1112           }
1113
1114         /* Note the above kernels may have been 'clipped' by a user defined
1115         ** radius, producing a smaller (darker) kernel.  Also for very small
1116         ** sigma's (> 0.1) the central value becomes larger than one, and thus
1117         ** producing a very bright kernel.
1118         **
1119         ** Normalization will still be needed.
1120         */
1121
1122         /* Normalize the 2D Gaussian Kernel
1123         **
1124         ** NB: a CorrelateNormalize performs a normal Normalize if
1125         ** there are no negative values.
1126         */
1127         CalcKernelMetaData(kernel);  /* the other kernel meta-data */
1128         ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1129
1130         break;
1131       }
1132     case BlurKernel:
1133       { double
1134           sigma = fabs(args->sigma),
1135           alpha, beta;
1136
1137         if ( args->rho >= 1.0 )
1138           kernel->width = (size_t)args->rho*2+1;
1139         else
1140           kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1141         kernel->height = 1;
1142         kernel->x = (ssize_t) (kernel->width-1)/2;
1143         kernel->y = 0;
1144         kernel->negative_range = kernel->positive_range = 0.0;
1145         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1146                               kernel->height*sizeof(double));
1147         if (kernel->values == (double *) NULL)
1148           return(DestroyKernelInfo(kernel));
1149
1150 #if 1
1151 #define KernelRank 3
1152         /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1153         ** It generates a gaussian 3 times the width, and compresses it into
1154         ** the expected range.  This produces a closer normalization of the
1155         ** resulting kernel, especially for very low sigma values.
1156         ** As such while wierd it is prefered.
1157         **
1158         ** I am told this method originally came from Photoshop.
1159         **
1160         ** A properly normalized curve is generated (apart from edge clipping)
1161         ** even though we later normalize the result (for edge clipping)
1162         ** to allow the correct generation of a "Difference of Blurs".
1163         */
1164
1165         /* initialize */
1166         v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1167         (void) ResetMagickMemory(kernel->values,0, (size_t)
1168                      kernel->width*kernel->height*sizeof(double));
1169         /* Calculate a Positive 1D Gaussian */
1170         if ( sigma > MagickEpsilon )
1171           { sigma *= KernelRank;               /* simplify loop expressions */
1172             alpha = 1.0/(2.0*sigma*sigma);
1173             beta= (double) (1.0/(MagickSQ2PI*sigma ));
1174             for ( u=-v; u <= v; u++) {
1175               kernel->values[(u+v)/KernelRank] +=
1176                               exp(-((double)(u*u))*alpha)*beta;
1177             }
1178           }
1179         else /* special case - generate a unity kernel */
1180           kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1181 #else
1182         /* Direct calculation without curve averaging */
1183
1184         /* Calculate a Positive Gaussian */
1185         if ( sigma > MagickEpsilon )
1186           { alpha = 1.0/(2.0*sigma*sigma);    /* simplify loop expressions */
1187             beta = 1.0/(MagickSQ2PI*sigma);
1188             for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1189               kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1190           }
1191         else /* special case - generate a unity kernel */
1192           { (void) ResetMagickMemory(kernel->values,0, (size_t)
1193                          kernel->width*kernel->height*sizeof(double));
1194             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1195           }
1196 #endif
1197         /* Note the above kernel may have been 'clipped' by a user defined
1198         ** radius, producing a smaller (darker) kernel.  Also for very small
1199         ** sigma's (> 0.1) the central value becomes larger than one, and thus
1200         ** producing a very bright kernel.
1201         **
1202         ** Normalization will still be needed.
1203         */
1204
1205         /* Normalize the 1D Gaussian Kernel
1206         **
1207         ** NB: a CorrelateNormalize performs a normal Normalize if
1208         ** there are no negative values.
1209         */
1210         CalcKernelMetaData(kernel);  /* the other kernel meta-data */
1211         ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1212
1213         /* rotate the 1D kernel by given angle */
1214         RotateKernelInfo(kernel, args->xi );
1215         break;
1216       }
1217     case CometKernel:
1218       { double
1219           sigma = fabs(args->sigma),
1220           A;
1221
1222         if ( args->rho < 1.0 )
1223           kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1224         else
1225           kernel->width = (size_t)args->rho;
1226         kernel->x = kernel->y = 0;
1227         kernel->height = 1;
1228         kernel->negative_range = kernel->positive_range = 0.0;
1229         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1230                               kernel->height*sizeof(double));
1231         if (kernel->values == (double *) NULL)
1232           return(DestroyKernelInfo(kernel));
1233
1234         /* A comet blur is half a 1D gaussian curve, so that the object is
1235         ** blurred in one direction only.  This may not be quite the right
1236         ** curve to use so may change in the future. The function must be
1237         ** normalised after generation, which also resolves any clipping.
1238         **
1239         ** As we are normalizing and not subtracting gaussians,
1240         ** there is no need for a divisor in the gaussian formula
1241         **
1242         ** It is less comples
1243         */
1244         if ( sigma > MagickEpsilon )
1245           {
1246 #if 1
1247 #define KernelRank 3
1248             v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1249             (void) ResetMagickMemory(kernel->values,0, (size_t)
1250                           kernel->width*sizeof(double));
1251             sigma *= KernelRank;            /* simplify the loop expression */
1252             A = 1.0/(2.0*sigma*sigma);
1253             /* B = 1.0/(MagickSQ2PI*sigma); */
1254             for ( u=0; u < v; u++) {
1255               kernel->values[u/KernelRank] +=
1256                   exp(-((double)(u*u))*A);
1257               /*  exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1258             }
1259             for (i=0; i < (ssize_t) kernel->width; i++)
1260               kernel->positive_range += kernel->values[i];
1261 #else
1262             A = 1.0/(2.0*sigma*sigma);     /* simplify the loop expression */
1263             /* B = 1.0/(MagickSQ2PI*sigma); */
1264             for ( i=0; i < (ssize_t) kernel->width; i++)
1265               kernel->positive_range +=
1266                 kernel->values[i] =
1267                   exp(-((double)(i*i))*A);
1268                 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1269 #endif
1270           }
1271         else /* special case - generate a unity kernel */
1272           { (void) ResetMagickMemory(kernel->values,0, (size_t)
1273                          kernel->width*kernel->height*sizeof(double));
1274             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1275             kernel->positive_range = 1.0;
1276           }
1277
1278         kernel->minimum = 0.0;
1279         kernel->maximum = kernel->values[0];
1280         kernel->negative_range = 0.0;
1281
1282         ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1283         RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1284         break;
1285       }
1286
1287     /*
1288       Convolution Kernels - Well Known Named Constant Kernels
1289     */
1290     case LaplacianKernel:
1291       { switch ( (int) args->rho ) {
1292           case 0:
1293           default: /* laplacian square filter -- default */
1294             kernel=ParseKernelArray("3: -1,-1,-1  -1,8,-1  -1,-1,-1");
1295             break;
1296           case 1:  /* laplacian diamond filter */
1297             kernel=ParseKernelArray("3: 0,-1,0  -1,4,-1  0,-1,0");
1298             break;
1299           case 2:
1300             kernel=ParseKernelArray("3: -2,1,-2  1,4,1  -2,1,-2");
1301             break;
1302           case 3:
1303             kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  1,-2,1");
1304             break;
1305           case 5:   /* a 5x5 laplacian */
1306             kernel=ParseKernelArray(
1307               "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");
1308             break;
1309           case 7:   /* a 7x7 laplacian */
1310             kernel=ParseKernelArray(
1311               "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" );
1312             break;
1313           case 15:  /* a 5x5 LoG (sigma approx 1.4) */
1314             kernel=ParseKernelArray(
1315               "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");
1316             break;
1317           case 19:  /* a 9x9 LoG (sigma approx 1.4) */
1318             /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1319             kernel=ParseKernelArray(
1320               "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");
1321             break;
1322         }
1323         if (kernel == (KernelInfo *) NULL)
1324           return(kernel);
1325         kernel->type = type;
1326         break;
1327       }
1328     case SobelKernel:
1329       { /* Simple Sobel Kernel */
1330         kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1331         if (kernel == (KernelInfo *) NULL)
1332           return(kernel);
1333         kernel->type = type;
1334         RotateKernelInfo(kernel, args->rho);
1335         break;
1336       }
1337     case RobertsKernel:
1338       {
1339         kernel=ParseKernelArray("3: 0,0,0  1,-1,0  0,0,0");
1340         if (kernel == (KernelInfo *) NULL)
1341           return(kernel);
1342         kernel->type = type;
1343         RotateKernelInfo(kernel, args->rho);
1344         break;
1345       }
1346     case PrewittKernel:
1347       {
1348         kernel=ParseKernelArray("3: 1,0,-1  1,0,-1  1,0,-1");
1349         if (kernel == (KernelInfo *) NULL)
1350           return(kernel);
1351         kernel->type = type;
1352         RotateKernelInfo(kernel, args->rho);
1353         break;
1354       }
1355     case CompassKernel:
1356       {
1357         kernel=ParseKernelArray("3: 1,1,-1  1,-2,-1  1,1,-1");
1358         if (kernel == (KernelInfo *) NULL)
1359           return(kernel);
1360         kernel->type = type;
1361         RotateKernelInfo(kernel, args->rho);
1362         break;
1363       }
1364     case KirschKernel:
1365       {
1366         kernel=ParseKernelArray("3: 5,-3,-3  5,0,-3  5,-3,-3");
1367         if (kernel == (KernelInfo *) NULL)
1368           return(kernel);
1369         kernel->type = type;
1370         RotateKernelInfo(kernel, args->rho);
1371         break;
1372       }
1373     case FreiChenKernel:
1374       /* Direction is set to be left to right positive */
1375       /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1376       /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1377       { switch ( (int) args->rho ) {
1378           default:
1379           case 0:
1380             kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1381             if (kernel == (KernelInfo *) NULL)
1382               return(kernel);
1383             kernel->type = type;
1384             kernel->values[3] = +MagickSQ2;
1385             kernel->values[5] = -MagickSQ2;
1386             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1387             break;
1388           case 2:
1389             kernel=ParseKernelArray("3: 1,2,0  2,0,-2  0,-2,-1");
1390             if (kernel == (KernelInfo *) NULL)
1391               return(kernel);
1392             kernel->type = type;
1393             kernel->values[1] = kernel->values[3] = +MagickSQ2;
1394             kernel->values[5] = kernel->values[7] = -MagickSQ2;
1395             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1396             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1397             break;
1398           case 10:
1399             kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1400             if (kernel == (KernelInfo *) NULL)
1401               return(kernel);
1402             break;
1403           case 1:
1404           case 11:
1405             kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1406             if (kernel == (KernelInfo *) NULL)
1407               return(kernel);
1408             kernel->type = type;
1409             kernel->values[3] = +MagickSQ2;
1410             kernel->values[5] = -MagickSQ2;
1411             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1412             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1413             break;
1414           case 12:
1415             kernel=ParseKernelArray("3: 1,2,1  0,0,0  1,2,1");
1416             if (kernel == (KernelInfo *) NULL)
1417               return(kernel);
1418             kernel->type = type;
1419             kernel->values[1] = +MagickSQ2;
1420             kernel->values[7] = +MagickSQ2;
1421             CalcKernelMetaData(kernel);
1422             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1423             break;
1424           case 13:
1425             kernel=ParseKernelArray("3: 2,-1,0  -1,0,1  0,1,-2");
1426             if (kernel == (KernelInfo *) NULL)
1427               return(kernel);
1428             kernel->type = type;
1429             kernel->values[0] = +MagickSQ2;
1430             kernel->values[8] = -MagickSQ2;
1431             CalcKernelMetaData(kernel);
1432             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1433             break;
1434           case 14:
1435             kernel=ParseKernelArray("3: 0,1,-2  -1,0,1  2,-1,0");
1436             if (kernel == (KernelInfo *) NULL)
1437               return(kernel);
1438             kernel->type = type;
1439             kernel->values[2] = -MagickSQ2;
1440             kernel->values[6] = +MagickSQ2;
1441             CalcKernelMetaData(kernel);
1442             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1443             break;
1444           case 15:
1445             kernel=ParseKernelArray("3: 0,-1,0  1,0,1  0,-1,0");
1446             if (kernel == (KernelInfo *) NULL)
1447               return(kernel);
1448             kernel->type = type;
1449             ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1450             break;
1451           case 16:
1452             kernel=ParseKernelArray("3: 1,0,-1  0,0,0  -1,0,1");
1453             if (kernel == (KernelInfo *) NULL)
1454               return(kernel);
1455             kernel->type = type;
1456             ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1457             break;
1458           case 17:
1459             kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  -1,-2,1");
1460             if (kernel == (KernelInfo *) NULL)
1461               return(kernel);
1462             kernel->type = type;
1463             ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1464             break;
1465           case 18:
1466             kernel=ParseKernelArray("3: -2,1,-2  1,4,1  -2,1,-2");
1467             if (kernel == (KernelInfo *) NULL)
1468               return(kernel);
1469             kernel->type = type;
1470             ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1471             break;
1472           case 19:
1473             kernel=ParseKernelArray("3: 1,1,1  1,1,1  1,1,1");
1474             if (kernel == (KernelInfo *) NULL)
1475               return(kernel);
1476             kernel->type = type;
1477             ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1478             break;
1479         }
1480         if ( fabs(args->sigma) > MagickEpsilon )
1481           /* Rotate by correctly supplied 'angle' */
1482           RotateKernelInfo(kernel, args->sigma);
1483         else if ( args->rho > 30.0 || args->rho < -30.0 )
1484           /* Rotate by out of bounds 'type' */
1485           RotateKernelInfo(kernel, args->rho);
1486         break;
1487       }
1488
1489     /*
1490       Boolean or Shaped Kernels
1491     */
1492     case DiamondKernel:
1493       {
1494         if (args->rho < 1.0)
1495           kernel->width = kernel->height = 3;  /* default radius = 1 */
1496         else
1497           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1498         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1499
1500         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1501                               kernel->height*sizeof(double));
1502         if (kernel->values == (double *) NULL)
1503           return(DestroyKernelInfo(kernel));
1504
1505         /* set all kernel values within diamond area to scale given */
1506         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1507           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1508             if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1509               kernel->positive_range += kernel->values[i] = args->sigma;
1510             else
1511               kernel->values[i] = nan;
1512         kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1513         break;
1514       }
1515     case SquareKernel:
1516     case RectangleKernel:
1517       { double
1518           scale;
1519         if ( type == SquareKernel )
1520           {
1521             if (args->rho < 1.0)
1522               kernel->width = kernel->height = 3;  /* default radius = 1 */
1523             else
1524               kernel->width = kernel->height = (size_t) (2*args->rho+1);
1525             kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1526             scale = args->sigma;
1527           }
1528         else {
1529             /* NOTE: user defaults set in "AcquireKernelInfo()" */
1530             if ( args->rho < 1.0 || args->sigma < 1.0 )
1531               return(DestroyKernelInfo(kernel));    /* invalid args given */
1532             kernel->width = (size_t)args->rho;
1533             kernel->height = (size_t)args->sigma;
1534             if ( args->xi  < 0.0 || args->xi  > (double)kernel->width ||
1535                  args->psi < 0.0 || args->psi > (double)kernel->height )
1536               return(DestroyKernelInfo(kernel));    /* invalid args given */
1537             kernel->x = (ssize_t) args->xi;
1538             kernel->y = (ssize_t) args->psi;
1539             scale = 1.0;
1540           }
1541         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1542                               kernel->height*sizeof(double));
1543         if (kernel->values == (double *) NULL)
1544           return(DestroyKernelInfo(kernel));
1545
1546         /* set all kernel values to scale given */
1547         u=(ssize_t) (kernel->width*kernel->height);
1548         for ( i=0; i < u; i++)
1549             kernel->values[i] = scale;
1550         kernel->minimum = kernel->maximum = scale;   /* a flat shape */
1551         kernel->positive_range = scale*u;
1552         break;
1553       }
1554       case OctagonKernel:
1555         {
1556           if (args->rho < 1.0)
1557             kernel->width = kernel->height = 5;  /* default radius = 2 */
1558           else
1559             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1560           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1561
1562           kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1563                                 kernel->height*sizeof(double));
1564           if (kernel->values == (double *) NULL)
1565             return(DestroyKernelInfo(kernel));
1566
1567           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1568             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1569               if ( (labs((long) u)+labs((long) v)) <=
1570                         ((long)kernel->x + (long)(kernel->x/2)) )
1571                 kernel->positive_range += kernel->values[i] = args->sigma;
1572               else
1573                 kernel->values[i] = nan;
1574           kernel->minimum = kernel->maximum = args->sigma;  /* a flat shape */
1575           break;
1576         }
1577       case DiskKernel:
1578         {
1579           ssize_t
1580             limit = (ssize_t)(args->rho*args->rho);
1581
1582           if (args->rho < 0.4)           /* default radius approx 4.3 */
1583             kernel->width = kernel->height = 9L, limit = 18L;
1584           else
1585             kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1586           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1587
1588           kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1589                                 kernel->height*sizeof(double));
1590           if (kernel->values == (double *) NULL)
1591             return(DestroyKernelInfo(kernel));
1592
1593           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1594             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1595               if ((u*u+v*v) <= limit)
1596                 kernel->positive_range += kernel->values[i] = args->sigma;
1597               else
1598                 kernel->values[i] = nan;
1599           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1600           break;
1601         }
1602       case PlusKernel:
1603         {
1604           if (args->rho < 1.0)
1605             kernel->width = kernel->height = 5;  /* default radius 2 */
1606           else
1607             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1608           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1609
1610           kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1611                                 kernel->height*sizeof(double));
1612           if (kernel->values == (double *) NULL)
1613             return(DestroyKernelInfo(kernel));
1614
1615           /* set all kernel values along axises to given scale */
1616           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1617             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1618               kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1619           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1620           kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1621           break;
1622         }
1623       case CrossKernel:
1624         {
1625           if (args->rho < 1.0)
1626             kernel->width = kernel->height = 5;  /* default radius 2 */
1627           else
1628             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1629           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1630
1631           kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1632                                 kernel->height*sizeof(double));
1633           if (kernel->values == (double *) NULL)
1634             return(DestroyKernelInfo(kernel));
1635
1636           /* set all kernel values along axises to given scale */
1637           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1638             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1639               kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1640           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1641           kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1642           break;
1643         }
1644       /*
1645         HitAndMiss Kernels
1646       */
1647       case RingKernel:
1648       case PeaksKernel:
1649         {
1650           ssize_t
1651             limit1,
1652             limit2,
1653             scale;
1654
1655           if (args->rho < args->sigma)
1656             {
1657               kernel->width = ((size_t)args->sigma)*2+1;
1658               limit1 = (ssize_t)(args->rho*args->rho);
1659               limit2 = (ssize_t)(args->sigma*args->sigma);
1660             }
1661           else
1662             {
1663               kernel->width = ((size_t)args->rho)*2+1;
1664               limit1 = (ssize_t)(args->sigma*args->sigma);
1665               limit2 = (ssize_t)(args->rho*args->rho);
1666             }
1667           if ( limit2 <= 0 )
1668             kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1669
1670           kernel->height = kernel->width;
1671           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1672           kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1673                                 kernel->height*sizeof(double));
1674           if (kernel->values == (double *) NULL)
1675             return(DestroyKernelInfo(kernel));
1676
1677           /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1678           scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1679           for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1680             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1681               { ssize_t radius=u*u+v*v;
1682                 if (limit1 < radius && radius <= limit2)
1683                   kernel->positive_range += kernel->values[i] = (double) scale;
1684                 else
1685                   kernel->values[i] = nan;
1686               }
1687           kernel->minimum = kernel->maximum = (double) scale;
1688           if ( type == PeaksKernel ) {
1689             /* set the central point in the middle */
1690             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1691             kernel->positive_range = 1.0;
1692             kernel->maximum = 1.0;
1693           }
1694           break;
1695         }
1696       case EdgesKernel:
1697         {
1698           kernel=AcquireKernelInfo("ThinSE:482");
1699           if (kernel == (KernelInfo *) NULL)
1700             return(kernel);
1701           kernel->type = type;
1702           ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1703           break;
1704         }
1705       case CornersKernel:
1706         {
1707           kernel=AcquireKernelInfo("ThinSE:87");
1708           if (kernel == (KernelInfo *) NULL)
1709             return(kernel);
1710           kernel->type = type;
1711           ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1712           break;
1713         }
1714       case DiagonalsKernel:
1715         {
1716           switch ( (int) args->rho ) {
1717             case 0:
1718             default:
1719               { KernelInfo
1720                   *new_kernel;
1721                 kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
1722                 if (kernel == (KernelInfo *) NULL)
1723                   return(kernel);
1724                 kernel->type = type;
1725                 new_kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
1726                 if (new_kernel == (KernelInfo *) NULL)
1727                   return(DestroyKernelInfo(kernel));
1728                 new_kernel->type = type;
1729                 LastKernelInfo(kernel)->next = new_kernel;
1730                 ExpandMirrorKernelInfo(kernel);
1731                 return(kernel);
1732               }
1733             case 1:
1734               kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
1735               break;
1736             case 2:
1737               kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
1738               break;
1739           }
1740           if (kernel == (KernelInfo *) NULL)
1741             return(kernel);
1742           kernel->type = type;
1743           RotateKernelInfo(kernel, args->sigma);
1744           break;
1745         }
1746       case LineEndsKernel:
1747         { /* Kernels for finding the end of thin lines */
1748           switch ( (int) args->rho ) {
1749             case 0:
1750             default:
1751               /* set of kernels to find all end of lines */
1752               return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1753             case 1:
1754               /* kernel for 4-connected line ends - no rotation */
1755               kernel=ParseKernelArray("3: 0,0,-  0,1,1  0,0,-");
1756               break;
1757           case 2:
1758               /* kernel to add for 8-connected lines - no rotation */
1759               kernel=ParseKernelArray("3: 0,0,0  0,1,0  0,0,1");
1760               break;
1761           case 3:
1762               /* kernel to add for orthogonal line ends - does not find corners */
1763               kernel=ParseKernelArray("3: 0,0,0  0,1,1  0,0,0");
1764               break;
1765           case 4:
1766               /* traditional line end - fails on last T end */
1767               kernel=ParseKernelArray("3: 0,0,0  0,1,-  0,0,-");
1768               break;
1769           }
1770           if (kernel == (KernelInfo *) NULL)
1771             return(kernel);
1772           kernel->type = type;
1773           RotateKernelInfo(kernel, args->sigma);
1774           break;
1775         }
1776       case LineJunctionsKernel:
1777         { /* kernels for finding the junctions of multiple lines */
1778           switch ( (int) args->rho ) {
1779             case 0:
1780             default:
1781               /* set of kernels to find all line junctions */
1782               return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1783             case 1:
1784               /* Y Junction */
1785               kernel=ParseKernelArray("3: 1,-,1  -,1,-  -,1,-");
1786               break;
1787             case 2:
1788               /* Diagonal T Junctions */
1789               kernel=ParseKernelArray("3: 1,-,-  -,1,-  1,-,1");
1790               break;
1791             case 3:
1792               /* Orthogonal T Junctions */
1793               kernel=ParseKernelArray("3: -,-,-  1,1,1  -,1,-");
1794               break;
1795             case 4:
1796               /* Diagonal X Junctions */
1797               kernel=ParseKernelArray("3: 1,-,1  -,1,-  1,-,1");
1798               break;
1799             case 5:
1800               /* Orthogonal X Junctions - minimal diamond kernel */
1801               kernel=ParseKernelArray("3: -,1,-  1,1,1  -,1,-");
1802               break;
1803           }
1804           if (kernel == (KernelInfo *) NULL)
1805             return(kernel);
1806           kernel->type = type;
1807           RotateKernelInfo(kernel, args->sigma);
1808           break;
1809         }
1810       case RidgesKernel:
1811         { /* Ridges - Ridge finding kernels */
1812           KernelInfo
1813             *new_kernel;
1814           switch ( (int) args->rho ) {
1815             case 1:
1816             default:
1817               kernel=ParseKernelArray("3x1:0,1,0");
1818               if (kernel == (KernelInfo *) NULL)
1819                 return(kernel);
1820               kernel->type = type;
1821               ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1822               break;
1823             case 2:
1824               kernel=ParseKernelArray("4x1:0,1,1,0");
1825               if (kernel == (KernelInfo *) NULL)
1826                 return(kernel);
1827               kernel->type = type;
1828               ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1829
1830               /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1831               /* Unfortunatally we can not yet rotate a non-square kernel */
1832               /* But then we can't flip a non-symetrical kernel either */
1833               new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1834               if (new_kernel == (KernelInfo *) NULL)
1835                 return(DestroyKernelInfo(kernel));
1836               new_kernel->type = type;
1837               LastKernelInfo(kernel)->next = new_kernel;
1838               new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1839               if (new_kernel == (KernelInfo *) NULL)
1840                 return(DestroyKernelInfo(kernel));
1841               new_kernel->type = type;
1842               LastKernelInfo(kernel)->next = new_kernel;
1843               new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1844               if (new_kernel == (KernelInfo *) NULL)
1845                 return(DestroyKernelInfo(kernel));
1846               new_kernel->type = type;
1847               LastKernelInfo(kernel)->next = new_kernel;
1848               new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1849               if (new_kernel == (KernelInfo *) NULL)
1850                 return(DestroyKernelInfo(kernel));
1851               new_kernel->type = type;
1852               LastKernelInfo(kernel)->next = new_kernel;
1853               new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1854               if (new_kernel == (KernelInfo *) NULL)
1855                 return(DestroyKernelInfo(kernel));
1856               new_kernel->type = type;
1857               LastKernelInfo(kernel)->next = new_kernel;
1858               new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1859               if (new_kernel == (KernelInfo *) NULL)
1860                 return(DestroyKernelInfo(kernel));
1861               new_kernel->type = type;
1862               LastKernelInfo(kernel)->next = new_kernel;
1863               new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1864               if (new_kernel == (KernelInfo *) NULL)
1865                 return(DestroyKernelInfo(kernel));
1866               new_kernel->type = type;
1867               LastKernelInfo(kernel)->next = new_kernel;
1868               new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1869               if (new_kernel == (KernelInfo *) NULL)
1870                 return(DestroyKernelInfo(kernel));
1871               new_kernel->type = type;
1872               LastKernelInfo(kernel)->next = new_kernel;
1873               break;
1874           }
1875           break;
1876         }
1877       case ConvexHullKernel:
1878         {
1879           KernelInfo
1880             *new_kernel;
1881           /* first set of 8 kernels */
1882           kernel=ParseKernelArray("3: 1,1,-  1,0,-  1,-,0");
1883           if (kernel == (KernelInfo *) NULL)
1884             return(kernel);
1885           kernel->type = type;
1886           ExpandRotateKernelInfo(kernel, 90.0);
1887           /* append the mirror versions too - no flip function yet */
1888           new_kernel=ParseKernelArray("3: 1,1,1  1,0,-  -,-,0");
1889           if (new_kernel == (KernelInfo *) NULL)
1890             return(DestroyKernelInfo(kernel));
1891           new_kernel->type = type;
1892           ExpandRotateKernelInfo(new_kernel, 90.0);
1893           LastKernelInfo(kernel)->next = new_kernel;
1894           break;
1895         }
1896       case SkeletonKernel:
1897         {
1898           switch ( (int) args->rho ) {
1899             case 1:
1900             default:
1901               /* Traditional Skeleton...
1902               ** A cyclically rotated single kernel
1903               */
1904               kernel=AcquireKernelInfo("ThinSE:482");
1905               if (kernel == (KernelInfo *) NULL)
1906                 return(kernel);
1907               kernel->type = type;
1908               ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1909               break;
1910             case 2:
1911               /* HIPR Variation of the cyclic skeleton
1912               ** Corners of the traditional method made more forgiving,
1913               ** but the retain the same cyclic order.
1914               */
1915               kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1916               if (kernel == (KernelInfo *) NULL)
1917                 return(kernel);
1918               if (kernel->next == (KernelInfo *) NULL)
1919                 return(DestroyKernelInfo(kernel));
1920               kernel->type = type;
1921               kernel->next->type = type;
1922               ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1923               break;
1924             case 3:
1925               /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1926               ** "Connectivity-Preserving Morphological Image Thransformations"
1927               ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1928               **   http://www.leptonica.com/papers/conn.pdf
1929               */
1930               kernel=AcquireKernelInfo(
1931                             "ThinSE:41; ThinSE:42; ThinSE:43");
1932               if (kernel == (KernelInfo *) NULL)
1933                 return(kernel);
1934               kernel->type = type;
1935               kernel->next->type = type;
1936               kernel->next->next->type = type;
1937               ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1938               break;
1939            }
1940           break;
1941         }
1942       case ThinSEKernel:
1943         { /* Special kernels for general thinning, while preserving connections
1944           ** "Connectivity-Preserving Morphological Image Thransformations"
1945           ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1946           **   http://www.leptonica.com/papers/conn.pdf
1947           ** And
1948           **   http://tpgit.github.com/Leptonica/ccthin_8c_source.html
1949           **
1950           ** Note kernels do not specify the origin pixel, allowing them
1951           ** to be used for both thickening and thinning operations.
1952           */
1953           switch ( (int) args->rho ) {
1954             /* SE for 4-connected thinning */
1955             case 41: /* SE_4_1 */
1956               kernel=ParseKernelArray("3: -,-,1  0,-,1  -,-,1");
1957               break;
1958             case 42: /* SE_4_2 */
1959               kernel=ParseKernelArray("3: -,-,1  0,-,1  -,0,-");
1960               break;
1961             case 43: /* SE_4_3 */
1962               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,-,1");
1963               break;
1964             case 44: /* SE_4_4 */
1965               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,-");
1966               break;
1967             case 45: /* SE_4_5 */
1968               kernel=ParseKernelArray("3: -,0,1  0,-,1  -,0,-");
1969               break;
1970             case 46: /* SE_4_6 */
1971               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,1");
1972               break;
1973             case 47: /* SE_4_7 */
1974               kernel=ParseKernelArray("3: -,1,1  0,-,1  -,0,-");
1975               break;
1976             case 48: /* SE_4_8 */
1977               kernel=ParseKernelArray("3: -,-,1  0,-,1  0,-,1");
1978               break;
1979             case 49: /* SE_4_9 */
1980               kernel=ParseKernelArray("3: 0,-,1  0,-,1  -,-,1");
1981               break;
1982             /* SE for 8-connected thinning - negatives of the above */
1983             case 81: /* SE_8_0 */
1984               kernel=ParseKernelArray("3: -,1,-  0,-,1  -,1,-");
1985               break;
1986             case 82: /* SE_8_2 */
1987               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,-,-");
1988               break;
1989             case 83: /* SE_8_3 */
1990               kernel=ParseKernelArray("3: 0,-,-  0,-,1  -,1,-");
1991               break;
1992             case 84: /* SE_8_4 */
1993               kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,-");
1994               break;
1995             case 85: /* SE_8_5 */
1996               kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,-");
1997               break;
1998             case 86: /* SE_8_6 */
1999               kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,1");
2000               break;
2001             case 87: /* SE_8_7 */
2002               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,0,-");
2003               break;
2004             case 88: /* SE_8_8 */
2005               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,1,-");
2006               break;
2007             case 89: /* SE_8_9 */
2008               kernel=ParseKernelArray("3: 0,1,-  0,-,1  -,1,-");
2009               break;
2010             /* Special combined SE kernels */
2011             case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2012               kernel=ParseKernelArray("3: -,-,1  0,-,-  -,0,-");
2013               break;
2014             case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2015               kernel=ParseKernelArray("3: -,1,-  -,-,1  0,-,-");
2016               break;
2017             case 481: /* SE_48_1 - General Connected Corner Kernel */
2018               kernel=ParseKernelArray("3: -,1,1  0,-,1  0,0,-");
2019               break;
2020             default:
2021             case 482: /* SE_48_2 - General Edge Kernel */
2022               kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,1");
2023               break;
2024           }
2025           if (kernel == (KernelInfo *) NULL)
2026             return(kernel);
2027           kernel->type = type;
2028           RotateKernelInfo(kernel, args->sigma);
2029           break;
2030         }
2031       /*
2032         Distance Measuring Kernels
2033       */
2034       case ChebyshevKernel:
2035         {
2036           if (args->rho < 1.0)
2037             kernel->width = kernel->height = 3;  /* default radius = 1 */
2038           else
2039             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2040           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2041
2042           kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2043                                 kernel->height*sizeof(double));
2044           if (kernel->values == (double *) NULL)
2045             return(DestroyKernelInfo(kernel));
2046
2047           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2048             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2049               kernel->positive_range += ( kernel->values[i] =
2050                   args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2051           kernel->maximum = kernel->values[0];
2052           break;
2053         }
2054       case ManhattanKernel:
2055         {
2056           if (args->rho < 1.0)
2057             kernel->width = kernel->height = 3;  /* default radius = 1 */
2058           else
2059             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2060           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2061
2062           kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2063                                 kernel->height*sizeof(double));
2064           if (kernel->values == (double *) NULL)
2065             return(DestroyKernelInfo(kernel));
2066
2067           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2068             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2069               kernel->positive_range += ( kernel->values[i] =
2070                   args->sigma*(labs((long) u)+labs((long) v)) );
2071           kernel->maximum = kernel->values[0];
2072           break;
2073         }
2074       case OctagonalKernel:
2075       {
2076         if (args->rho < 2.0)
2077           kernel->width = kernel->height = 5;  /* default/minimum radius = 2 */
2078         else
2079           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2080         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2081
2082         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2083                               kernel->height*sizeof(double));
2084         if (kernel->values == (double *) NULL)
2085           return(DestroyKernelInfo(kernel));
2086
2087         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2088           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2089             {
2090               double
2091                 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2092                 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2093               kernel->positive_range += kernel->values[i] =
2094                         args->sigma*MagickMax(r1,r2);
2095             }
2096         kernel->maximum = kernel->values[0];
2097         break;
2098       }
2099     case EuclideanKernel:
2100       {
2101         if (args->rho < 1.0)
2102           kernel->width = kernel->height = 3;  /* default radius = 1 */
2103         else
2104           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2105         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2106
2107         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2108                               kernel->height*sizeof(double));
2109         if (kernel->values == (double *) NULL)
2110           return(DestroyKernelInfo(kernel));
2111
2112         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2113           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2114             kernel->positive_range += ( kernel->values[i] =
2115                  args->sigma*sqrt((double)(u*u+v*v)) );
2116         kernel->maximum = kernel->values[0];
2117         break;
2118       }
2119     default:
2120       {
2121         /* No-Op Kernel - Basically just a single pixel on its own */
2122         kernel=ParseKernelArray("1:1");
2123         if (kernel == (KernelInfo *) NULL)
2124           return(kernel);
2125         kernel->type = UndefinedKernel;
2126         break;
2127       }
2128       break;
2129   }
2130   return(kernel);
2131 }
2132 \f
2133 /*
2134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2135 %                                                                             %
2136 %                                                                             %
2137 %                                                                             %
2138 %     C l o n e K e r n e l I n f o                                           %
2139 %                                                                             %
2140 %                                                                             %
2141 %                                                                             %
2142 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2143 %
2144 %  CloneKernelInfo() creates a new clone of the given Kernel List so that its
2145 %  can be modified without effecting the original.  The cloned kernel should
2146 %  be destroyed using DestoryKernelInfo() when no longer needed.
2147 %
2148 %  The format of the CloneKernelInfo method is:
2149 %
2150 %      KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2151 %
2152 %  A description of each parameter follows:
2153 %
2154 %    o kernel: the Morphology/Convolution kernel to be cloned
2155 %
2156 */
2157 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2158 {
2159   register ssize_t
2160     i;
2161
2162   KernelInfo
2163     *new_kernel;
2164
2165   assert(kernel != (KernelInfo *) NULL);
2166   new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2167   if (new_kernel == (KernelInfo *) NULL)
2168     return(new_kernel);
2169   *new_kernel=(*kernel); /* copy values in structure */
2170
2171   /* replace the values with a copy of the values */
2172   new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2173     kernel->height*sizeof(double));
2174   if (new_kernel->values == (double *) NULL)
2175     return(DestroyKernelInfo(new_kernel));
2176   for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2177     new_kernel->values[i]=kernel->values[i];
2178
2179   /* Also clone the next kernel in the kernel list */
2180   if ( kernel->next != (KernelInfo *) NULL ) {
2181     new_kernel->next = CloneKernelInfo(kernel->next);
2182     if ( new_kernel->next == (KernelInfo *) NULL )
2183       return(DestroyKernelInfo(new_kernel));
2184   }
2185
2186   return(new_kernel);
2187 }
2188 \f
2189 /*
2190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2191 %                                                                             %
2192 %                                                                             %
2193 %                                                                             %
2194 %     D e s t r o y K e r n e l I n f o                                       %
2195 %                                                                             %
2196 %                                                                             %
2197 %                                                                             %
2198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2199 %
2200 %  DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2201 %  kernel.
2202 %
2203 %  The format of the DestroyKernelInfo method is:
2204 %
2205 %      KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2206 %
2207 %  A description of each parameter follows:
2208 %
2209 %    o kernel: the Morphology/Convolution kernel to be destroyed
2210 %
2211 */
2212 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2213 {
2214   assert(kernel != (KernelInfo *) NULL);
2215
2216   if ( kernel->next != (KernelInfo *) NULL )
2217     kernel->next = DestroyKernelInfo(kernel->next);
2218
2219   kernel->values = (double *)RelinquishMagickMemory(kernel->values);
2220   kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
2221   return(kernel);
2222 }
2223 \f
2224 /*
2225 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2226 %                                                                             %
2227 %                                                                             %
2228 %                                                                             %
2229 +     E x p a n d M i r r o r K e r n e l I n f o                             %
2230 %                                                                             %
2231 %                                                                             %
2232 %                                                                             %
2233 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2234 %
2235 %  ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2236 %  sequence of 90-degree rotated kernels but providing a reflected 180
2237 %  rotatation, before the -/+ 90-degree rotations.
2238 %
2239 %  This special rotation order produces a better, more symetrical thinning of
2240 %  objects.
2241 %
2242 %  The format of the ExpandMirrorKernelInfo method is:
2243 %
2244 %      void ExpandMirrorKernelInfo(KernelInfo *kernel)
2245 %
2246 %  A description of each parameter follows:
2247 %
2248 %    o kernel: the Morphology/Convolution kernel
2249 %
2250 % This function is only internel to this module, as it is not finalized,
2251 % especially with regard to non-orthogonal angles, and rotation of larger
2252 % 2D kernels.
2253 */
2254
2255 #if 0
2256 static void FlopKernelInfo(KernelInfo *kernel)
2257     { /* Do a Flop by reversing each row. */
2258       size_t
2259         y;
2260       register ssize_t
2261         x,r;
2262       register double
2263         *k,t;
2264
2265       for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2266         for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2267           t=k[x],  k[x]=k[r],  k[r]=t;
2268
2269       kernel->x = kernel->width - kernel->x - 1;
2270       angle = fmod(angle+180.0, 360.0);
2271     }
2272 #endif
2273
2274 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2275 {
2276   KernelInfo
2277     *clone,
2278     *last;
2279
2280   last = kernel;
2281
2282   clone = CloneKernelInfo(last);
2283   RotateKernelInfo(clone, 180);   /* flip */
2284   LastKernelInfo(last)->next = clone;
2285   last = clone;
2286
2287   clone = CloneKernelInfo(last);
2288   RotateKernelInfo(clone, 90);   /* transpose */
2289   LastKernelInfo(last)->next = clone;
2290   last = clone;
2291
2292   clone = CloneKernelInfo(last);
2293   RotateKernelInfo(clone, 180);  /* flop */
2294   LastKernelInfo(last)->next = clone;
2295
2296   return;
2297 }
2298 \f
2299 /*
2300 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2301 %                                                                             %
2302 %                                                                             %
2303 %                                                                             %
2304 +     E x p a n d R o t a t e K e r n e l I n f o                             %
2305 %                                                                             %
2306 %                                                                             %
2307 %                                                                             %
2308 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2309 %
2310 %  ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2311 %  incrementally by the angle given, until the kernel repeats.
2312 %
2313 %  WARNING: 45 degree rotations only works for 3x3 kernels.
2314 %  While 90 degree roatations only works for linear and square kernels
2315 %
2316 %  The format of the ExpandRotateKernelInfo method is:
2317 %
2318 %      void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2319 %
2320 %  A description of each parameter follows:
2321 %
2322 %    o kernel: the Morphology/Convolution kernel
2323 %
2324 %    o angle: angle to rotate in degrees
2325 %
2326 % This function is only internel to this module, as it is not finalized,
2327 % especially with regard to non-orthogonal angles, and rotation of larger
2328 % 2D kernels.
2329 */
2330
2331 /* Internal Routine - Return true if two kernels are the same */
2332 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2333      const KernelInfo *kernel2)
2334 {
2335   register size_t
2336     i;
2337
2338   /* check size and origin location */
2339   if (    kernel1->width != kernel2->width
2340        || kernel1->height != kernel2->height
2341        || kernel1->x != kernel2->x
2342        || kernel1->y != kernel2->y )
2343     return MagickFalse;
2344
2345   /* check actual kernel values */
2346   for (i=0; i < (kernel1->width*kernel1->height); i++) {
2347     /* Test for Nan equivalence */
2348     if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2349       return MagickFalse;
2350     if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2351       return MagickFalse;
2352     /* Test actual values are equivalent */
2353     if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
2354       return MagickFalse;
2355   }
2356
2357   return MagickTrue;
2358 }
2359
2360 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2361 {
2362   KernelInfo
2363     *clone,
2364     *last;
2365
2366   last = kernel;
2367   while(1) {
2368     clone = CloneKernelInfo(last);
2369     RotateKernelInfo(clone, angle);
2370     if ( SameKernelInfo(kernel, clone) == MagickTrue )
2371       break;
2372     LastKernelInfo(last)->next = clone;
2373     last = clone;
2374   }
2375   clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2376   return;
2377 }
2378 \f
2379 /*
2380 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2381 %                                                                             %
2382 %                                                                             %
2383 %                                                                             %
2384 +     C a l c M e t a K e r n a l I n f o                                     %
2385 %                                                                             %
2386 %                                                                             %
2387 %                                                                             %
2388 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2389 %
2390 %  CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2391 %  using the kernel values.  This should only ne used if it is not possible to
2392 %  calculate that meta-data in some easier way.
2393 %
2394 %  It is important that the meta-data is correct before ScaleKernelInfo() is
2395 %  used to perform kernel normalization.
2396 %
2397 %  The format of the CalcKernelMetaData method is:
2398 %
2399 %      void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2400 %
2401 %  A description of each parameter follows:
2402 %
2403 %    o kernel: the Morphology/Convolution kernel to modify
2404 %
2405 %  WARNING: Minimum and Maximum values are assumed to include zero, even if
2406 %  zero is not part of the kernel (as in Gaussian Derived kernels). This
2407 %  however is not true for flat-shaped morphological kernels.
2408 %
2409 %  WARNING: Only the specific kernel pointed to is modified, not a list of
2410 %  multiple kernels.
2411 %
2412 % This is an internal function and not expected to be useful outside this
2413 % module.  This could change however.
2414 */
2415 static void CalcKernelMetaData(KernelInfo *kernel)
2416 {
2417   register size_t
2418     i;
2419
2420   kernel->minimum = kernel->maximum = 0.0;
2421   kernel->negative_range = kernel->positive_range = 0.0;
2422   for (i=0; i < (kernel->width*kernel->height); i++)
2423     {
2424       if ( fabs(kernel->values[i]) < MagickEpsilon )
2425         kernel->values[i] = 0.0;
2426       ( kernel->values[i] < 0)
2427           ?  ( kernel->negative_range += kernel->values[i] )
2428           :  ( kernel->positive_range += kernel->values[i] );
2429       Minimize(kernel->minimum, kernel->values[i]);
2430       Maximize(kernel->maximum, kernel->values[i]);
2431     }
2432
2433   return;
2434 }
2435 \f
2436 /*
2437 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2438 %                                                                             %
2439 %                                                                             %
2440 %                                                                             %
2441 %     M o r p h o l o g y A p p l y                                           %
2442 %                                                                             %
2443 %                                                                             %
2444 %                                                                             %
2445 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2446 %
2447 %  MorphologyApply() applies a morphological method, multiple times using
2448 %  a list of multiple kernels.
2449 %
2450 %  It is basically equivalent to as MorphologyImage() (see below) but
2451 %  without any user controls.  This allows internel programs to use this
2452 %  function, to actually perform a specific task without possible interference
2453 %  by any API user supplied settings.
2454 %
2455 %  It is MorphologyImage() task to extract any such user controls, and
2456 %  pass them to this function for processing.
2457 %
2458 %  More specifically kernels are not normalized/scaled/blended by the
2459 %  'convolve:scale' Image Artifact (setting), nor is the convolve bias
2460 %  (-bias setting or image->bias) loooked at, but must be supplied from the
2461 %  function arguments.
2462 %
2463 %  The format of the MorphologyApply method is:
2464 %
2465 %      Image *MorphologyApply(const Image *image,MorphologyMethod method,
2466 %        const ssize_t iterations,const KernelInfo *kernel,
2467 %        const CompositeMethod compose,const double bias,
2468 %        ExceptionInfo *exception)
2469 %
2470 %  A description of each parameter follows:
2471 %
2472 %    o image: the source image
2473 %
2474 %    o method: the morphology method to be applied.
2475 %
2476 %    o iterations: apply the operation this many times (or no change).
2477 %                  A value of -1 means loop until no change found.
2478 %                  How this is applied may depend on the morphology method.
2479 %                  Typically this is a value of 1.
2480 %
2481 %    o channel: the channel type.
2482 %
2483 %    o kernel: An array of double representing the morphology kernel.
2484 %
2485 %    o compose: How to handle or merge multi-kernel results.
2486 %          If 'UndefinedCompositeOp' use default for the Morphology method.
2487 %          If 'NoCompositeOp' force image to be re-iterated by each kernel.
2488 %          Otherwise merge the results using the compose method given.
2489 %
2490 %    o bias: Convolution Output Bias.
2491 %
2492 %    o exception: return any errors or warnings in this structure.
2493 %
2494 */
2495
2496 /* Apply a Morphology Primative to an image using the given kernel.
2497 ** Two pre-created images must be provided, and no image is created.
2498 ** It returns the number of pixels that changed between the images
2499 ** for result convergence determination.
2500 */
2501 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2502   const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2503   ExceptionInfo *exception)
2504 {
2505 #define MorphologyTag  "Morphology/Image"
2506
2507   CacheView
2508     *image_view,
2509     *morphology_view;
2510
2511   ssize_t
2512     y, offx, offy;
2513
2514   size_t
2515     virt_width,
2516     changed;
2517
2518   MagickBooleanType
2519     status;
2520
2521   MagickOffsetType
2522     progress;
2523
2524   assert(image != (Image *) NULL);
2525   assert(image->signature == MagickSignature);
2526   assert(morphology_image != (Image *) NULL);
2527   assert(morphology_image->signature == MagickSignature);
2528   assert(kernel != (KernelInfo *) NULL);
2529   assert(kernel->signature == MagickSignature);
2530   assert(exception != (ExceptionInfo *) NULL);
2531   assert(exception->signature == MagickSignature);
2532
2533   status=MagickTrue;
2534   changed=0;
2535   progress=0;
2536
2537   image_view=AcquireCacheView(image);
2538   morphology_view=AcquireCacheView(morphology_image);
2539   virt_width=image->columns+kernel->width-1;
2540
2541   /* Some methods (including convolve) needs use a reflected kernel.
2542    * Adjust 'origin' offsets to loop though kernel as a reflection.
2543    */
2544   offx = kernel->x;
2545   offy = kernel->y;
2546   switch(method) {
2547     case ConvolveMorphology:
2548     case DilateMorphology:
2549     case DilateIntensityMorphology:
2550     /*case DistanceMorphology:*/
2551       /* kernel needs to used with reflection about origin */
2552       offx = (ssize_t) kernel->width-offx-1;
2553       offy = (ssize_t) kernel->height-offy-1;
2554       break;
2555     case ErodeMorphology:
2556     case ErodeIntensityMorphology:
2557     case HitAndMissMorphology:
2558     case ThinningMorphology:
2559     case ThickenMorphology:
2560       /* kernel is used as is, without reflection */
2561       break;
2562     default:
2563       assert("Not a Primitive Morphology Method" != (char *) NULL);
2564       break;
2565   }
2566
2567
2568   if ( method == ConvolveMorphology && kernel->width == 1 )
2569   { /* Special handling (for speed) of vertical (blur) kernels.
2570     ** This performs its handling in columns rather than in rows.
2571     ** This is only done for convolve as it is the only method that
2572     ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
2573     **
2574     ** Timing tests (on single CPU laptop)
2575     ** Using a vertical 1-d Blue with normal row-by-row (below)
2576     **   time convert logo: -morphology Convolve Blur:0x10+90 null:
2577     **      0.807u
2578     ** Using this column method
2579     **   time convert logo: -morphology Convolve Blur:0x10+90 null:
2580     **      0.620u
2581     **
2582     ** Anthony Thyssen, 14 June 2010
2583     */
2584     register ssize_t
2585       x;
2586
2587 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2588 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2589 #endif
2590     for (x=0; x < (ssize_t) image->columns; x++)
2591     {
2592       register const Quantum
2593         *restrict p;
2594
2595       register Quantum
2596         *restrict q;
2597
2598       register ssize_t
2599         y;
2600
2601       ssize_t
2602         r;
2603
2604       if (status == MagickFalse)
2605         continue;
2606       p=GetCacheViewVirtualPixels(image_view,x,-offy,1,image->rows+
2607         kernel->height-1,exception);
2608       q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2609         morphology_image->rows,exception);
2610       if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2611         {
2612           status=MagickFalse;
2613           continue;
2614         }
2615       /* offset to origin in 'p'. while 'q' points to it directly */
2616       r = offy;
2617
2618       for (y=0; y < (ssize_t) image->rows; y++)
2619       {
2620         PixelInfo
2621           result;
2622
2623         register ssize_t
2624           v;
2625
2626         register const double
2627           *restrict k;
2628
2629         register const Quantum
2630           *restrict k_pixels;
2631
2632         /* Copy input image to the output image for unused channels
2633         * This removes need for 'cloning' a new image every iteration
2634         */
2635         SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2636           GetPixelChannels(image)),q);
2637         SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2638           GetPixelChannels(image)),q);
2639         SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2640           GetPixelChannels(image)),q);
2641         if (image->colorspace == CMYKColorspace)
2642           SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2643             GetPixelChannels(image)),q);
2644
2645         /* Set the bias of the weighted average output */
2646         result.red     =
2647         result.green   =
2648         result.blue    =
2649         result.alpha =
2650         result.black   = bias;
2651
2652
2653         /* Weighted Average of pixels using reflected kernel
2654         **
2655         ** NOTE for correct working of this operation for asymetrical
2656         ** kernels, the kernel needs to be applied in its reflected form.
2657         ** That is its values needs to be reversed.
2658         */
2659         k = &kernel->values[ kernel->height-1 ];
2660         k_pixels = p;
2661         if ( (image->sync == MagickFalse) || (image->matte == MagickFalse) )
2662           { /* No 'Sync' involved.
2663             ** Convolution is simple greyscale channel operation
2664             */
2665             for (v=0; v < (ssize_t) kernel->height; v++) {
2666               if ( IsNan(*k) ) continue;
2667               result.red     += (*k)*GetPixelRed(image,k_pixels);
2668               result.green   += (*k)*GetPixelGreen(image,k_pixels);
2669               result.blue    += (*k)*GetPixelBlue(image,k_pixels);
2670               if (image->colorspace == CMYKColorspace)
2671                 result.black+=(*k)*GetPixelBlack(image,k_pixels);
2672               result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2673               k--;
2674               k_pixels+=GetPixelChannels(image);
2675             }
2676             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2677               SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
2678             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2679               SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
2680             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2681               SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
2682             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2683                 (image->colorspace == CMYKColorspace))
2684               SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
2685             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2686                 (image->matte == MagickTrue))
2687               SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2688           }
2689         else
2690           { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2691             ** Weight the color channels with Alpha Channel so that
2692             ** transparent pixels are not part of the results.
2693             */
2694             MagickRealType
2695               alpha,  /* alpha weighting of colors : kernel*alpha  */
2696               gamma;  /* divisor, sum of color weighting values */
2697
2698             gamma=0.0;
2699             for (v=0; v < (ssize_t) kernel->height; v++) {
2700               if ( IsNan(*k) ) continue;
2701               alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels));
2702               gamma += alpha;
2703               result.red     += alpha*GetPixelRed(image,k_pixels);
2704               result.green   += alpha*GetPixelGreen(image,k_pixels);
2705               result.blue    += alpha*GetPixelBlue(image,k_pixels);
2706               if (image->colorspace == CMYKColorspace)
2707                 result.black += alpha*GetPixelBlack(image,k_pixels);
2708               result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2709               k--;
2710               k_pixels+=GetPixelChannels(image);
2711             }
2712             /* Sync'ed channels, all channels are modified */
2713             gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2714             SetPixelRed(morphology_image,
2715               ClampToQuantum(gamma*result.red),q);
2716             SetPixelGreen(morphology_image,
2717               ClampToQuantum(gamma*result.green),q);
2718             SetPixelBlue(morphology_image,
2719               ClampToQuantum(gamma*result.blue),q);
2720             if (image->colorspace == CMYKColorspace)
2721               SetPixelBlack(morphology_image,
2722                 ClampToQuantum(gamma*result.black),q);
2723             SetPixelAlpha(morphology_image,
2724               ClampToQuantum(result.alpha),q);
2725           }
2726
2727         /* Count up changed pixels */
2728         if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q))
2729             || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q))
2730             || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q))
2731             || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q))
2732             || ((image->colorspace == CMYKColorspace) &&
2733                 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
2734           changed++;  /* The pixel was changed in some way! */
2735         p+=GetPixelChannels(image);
2736         q+=GetPixelChannels(morphology_image);
2737       } /* y */
2738       if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2739         status=MagickFalse;
2740       if (image->progress_monitor != (MagickProgressMonitor) NULL)
2741         {
2742           MagickBooleanType
2743             proceed;
2744
2745 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2746   #pragma omp critical (MagickCore_MorphologyImage)
2747 #endif
2748           proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2749           if (proceed == MagickFalse)
2750             status=MagickFalse;
2751         }
2752     } /* x */
2753     morphology_image->type=image->type;
2754     morphology_view=DestroyCacheView(morphology_view);
2755     image_view=DestroyCacheView(image_view);
2756     return(status ? (ssize_t) changed : 0);
2757   }
2758
2759   /*
2760   ** Normal handling of horizontal or rectangular kernels (row by row)
2761   */
2762 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2763   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2764 #endif
2765   for (y=0; y < (ssize_t) image->rows; y++)
2766   {
2767     register const Quantum
2768       *restrict p;
2769
2770     register Quantum
2771       *restrict q;
2772
2773     register ssize_t
2774       x;
2775
2776     size_t
2777       r;
2778
2779     if (status == MagickFalse)
2780       continue;
2781     p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
2782       kernel->height,  exception);
2783     q=GetCacheViewAuthenticPixels(morphology_view,0,y,
2784       morphology_image->columns,1,exception);
2785     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2786       {
2787         status=MagickFalse;
2788         continue;
2789       }
2790     /* offset to origin in 'p'. while 'q' points to it directly */
2791     r = virt_width*offy + offx;
2792
2793     for (x=0; x < (ssize_t) image->columns; x++)
2794     {
2795        ssize_t
2796         v;
2797
2798       register ssize_t
2799         u;
2800
2801       register const double
2802         *restrict k;
2803
2804       register const Quantum
2805         *restrict k_pixels;
2806
2807       PixelInfo
2808         result,
2809         min,
2810         max;
2811
2812       /* Copy input image to the output image for unused channels
2813        * This removes need for 'cloning' a new image every iteration
2814        */
2815       SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2816         GetPixelChannels(image)),q);
2817       SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2818         GetPixelChannels(image)),q);
2819       SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2820         GetPixelChannels(image)),q);
2821       if (image->colorspace == CMYKColorspace)
2822         SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2823           GetPixelChannels(image)),q);
2824
2825       /* Defaults */
2826       min.red     =
2827       min.green   =
2828       min.blue    =
2829       min.alpha =
2830       min.black   = (MagickRealType) QuantumRange;
2831       max.red     =
2832       max.green   =
2833       max.blue    =
2834       max.alpha =
2835       max.black   = (MagickRealType) 0;
2836       /* default result is the original pixel value */
2837       result.red     = (MagickRealType) GetPixelRed(image,p+r*GetPixelChannels(image));
2838       result.green   = (MagickRealType) GetPixelGreen(image,p+r*GetPixelChannels(image));
2839       result.blue    = (MagickRealType) GetPixelBlue(image,p+r*GetPixelChannels(image));
2840       result.black   = 0.0;
2841       if (image->colorspace == CMYKColorspace)
2842         result.black = (MagickRealType) GetPixelBlack(image,p+r*GetPixelChannels(image));
2843       result.alpha=(MagickRealType) GetPixelAlpha(image,p+r*GetPixelChannels(image));
2844
2845       switch (method) {
2846         case ConvolveMorphology:
2847           /* Set the bias of the weighted average output */
2848           result.red     =
2849           result.green   =
2850           result.blue    =
2851           result.alpha =
2852           result.black   = bias;
2853           break;
2854         case DilateIntensityMorphology:
2855         case ErodeIntensityMorphology:
2856           /* use a boolean flag indicating when first match found */
2857           result.red = 0.0;  /* result is not used otherwise */
2858           break;
2859         default:
2860           break;
2861       }
2862
2863       switch ( method ) {
2864         case ConvolveMorphology:
2865             /* Weighted Average of pixels using reflected kernel
2866             **
2867             ** NOTE for correct working of this operation for asymetrical
2868             ** kernels, the kernel needs to be applied in its reflected form.
2869             ** That is its values needs to be reversed.
2870             **
2871             ** Correlation is actually the same as this but without reflecting
2872             ** the kernel, and thus 'lower-level' that Convolution.  However
2873             ** as Convolution is the more common method used, and it does not
2874             ** really cost us much in terms of processing to use a reflected
2875             ** kernel, so it is Convolution that is implemented.
2876             **
2877             ** Correlation will have its kernel reflected before calling
2878             ** this function to do a Convolve.
2879             **
2880             ** For more details of Correlation vs Convolution see
2881             **   http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2882             */
2883             k = &kernel->values[ kernel->width*kernel->height-1 ];
2884             k_pixels = p;
2885             if ( (image->sync == MagickFalse) ||
2886                  (image->matte == MagickFalse) )
2887               { /* No 'Sync' involved.
2888                 ** Convolution is simple greyscale channel operation
2889                 */
2890                 for (v=0; v < (ssize_t) kernel->height; v++) {
2891                   for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2892                     if ( IsNan(*k) ) continue;
2893                     result.red     += (*k)*
2894                       GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2895                     result.green   += (*k)*
2896                       GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2897                     result.blue    += (*k)*
2898                       GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2899                     if (image->colorspace == CMYKColorspace)
2900                       result.black += (*k)*
2901                         GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2902                     result.alpha += (*k)*
2903                       GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2904                   }
2905                   k_pixels += virt_width*GetPixelChannels(image);
2906                 }
2907                 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2908                   SetPixelRed(morphology_image,ClampToQuantum(result.red),
2909                     q);
2910                 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2911                   SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2912                     q);
2913                 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2914                   SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2915                     q);
2916                 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2917                     (image->colorspace == CMYKColorspace))
2918                   SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2919                     q);
2920                 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2921                     (image->matte == MagickTrue))
2922                   SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
2923                     q);
2924               }
2925             else
2926               { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2927                 ** Weight the color channels with Alpha Channel so that
2928                 ** transparent pixels are not part of the results.
2929                 */
2930                 MagickRealType
2931                   alpha,  /* alpha weighting of colors : kernel*alpha  */
2932                   gamma;  /* divisor, sum of color weighting values */
2933
2934                 gamma=0.0;
2935                 for (v=0; v < (ssize_t) kernel->height; v++) {
2936                   for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2937                     if ( IsNan(*k) ) continue;
2938                     alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels+u*
2939                       GetPixelChannels(image)));
2940                     gamma += alpha;
2941                     result.red     += alpha*
2942                       GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2943                     result.green   += alpha*
2944                       GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2945                     result.blue    += alpha*
2946                       GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2947                     if (image->colorspace == CMYKColorspace)
2948                       result.black+=alpha*
2949                         GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2950                     result.alpha += (*k)*
2951                       GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2952                   }
2953                   k_pixels += virt_width*GetPixelChannels(image);
2954                 }
2955                 /* Sync'ed channels, all channels are modified */
2956                 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2957                 SetPixelRed(morphology_image,
2958                   ClampToQuantum(gamma*result.red),q);
2959                 SetPixelGreen(morphology_image,
2960                   ClampToQuantum(gamma*result.green),q);
2961                 SetPixelBlue(morphology_image,
2962                   ClampToQuantum(gamma*result.blue),q);
2963                 if (image->colorspace == CMYKColorspace)
2964                   SetPixelBlack(morphology_image,
2965                     ClampToQuantum(gamma*result.black),q);
2966                 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2967               }
2968             break;
2969
2970         case ErodeMorphology:
2971             /* Minimum Value within kernel neighbourhood
2972             **
2973             ** NOTE that the kernel is not reflected for this operation!
2974             **
2975             ** NOTE: in normal Greyscale Morphology, the kernel value should
2976             ** be added to the real value, this is currently not done, due to
2977             ** the nature of the boolean kernels being used.
2978             */
2979             k = kernel->values;
2980             k_pixels = p;
2981             for (v=0; v < (ssize_t) kernel->height; v++) {
2982               for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2983                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2984                 Minimize(min.red,     (double)
2985                   GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
2986                 Minimize(min.green,   (double)
2987                   GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
2988                 Minimize(min.blue,    (double)
2989                   GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
2990                 if (image->colorspace == CMYKColorspace)
2991                   Minimize(min.black,(double)
2992                     GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
2993                 Minimize(min.alpha,(double)
2994                   GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
2995               }
2996               k_pixels += virt_width*GetPixelChannels(image);
2997             }
2998             break;
2999
3000         case DilateMorphology:
3001             /* Maximum Value within kernel neighbourhood
3002             **
3003             ** NOTE for correct working of this operation for asymetrical
3004             ** kernels, the kernel needs to be applied in its reflected form.
3005             ** That is its values needs to be reversed.
3006             **
3007             ** NOTE: in normal Greyscale Morphology, the kernel value should
3008             ** be added to the real value, this is currently not done, due to
3009             ** the nature of the boolean kernels being used.
3010             **
3011             */
3012             k = &kernel->values[ kernel->width*kernel->height-1 ];
3013             k_pixels = p;
3014             for (v=0; v < (ssize_t) kernel->height; v++) {
3015               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3016                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3017                 Maximize(max.red,     (double)
3018                   GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3019                 Maximize(max.green,   (double)
3020                   GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3021                 Maximize(max.blue,    (double)
3022                   GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3023                 if (image->colorspace == CMYKColorspace)
3024                   Maximize(max.black,   (double)
3025                     GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3026                 Maximize(max.alpha,(double)
3027                   GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3028               }
3029               k_pixels += virt_width*GetPixelChannels(image);
3030             }
3031             break;
3032
3033         case HitAndMissMorphology:
3034         case ThinningMorphology:
3035         case ThickenMorphology:
3036             /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3037             **
3038             ** NOTE that the kernel is not reflected for this operation,
3039             ** and consists of both foreground and background pixel
3040             ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3041             ** with either Nan or 0.5 values for don't care.
3042             **
3043             ** Note that this will never produce a meaningless negative
3044             ** result.  Such results can cause Thinning/Thicken to not work
3045             ** correctly when used against a greyscale image.
3046             */
3047             k = kernel->values;
3048             k_pixels = p;
3049             for (v=0; v < (ssize_t) kernel->height; v++) {
3050               for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3051                 if ( IsNan(*k) ) continue;
3052                 if ( (*k) > 0.7 )
3053                 { /* minimim of foreground pixels */
3054                   Minimize(min.red,     (double)
3055                     GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3056                   Minimize(min.green,   (double)
3057                     GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3058                   Minimize(min.blue,    (double)
3059                     GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3060                   if ( image->colorspace == CMYKColorspace)
3061                     Minimize(min.black,(double)
3062                       GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3063                   Minimize(min.alpha,(double)
3064                     GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3065                 }
3066                 else if ( (*k) < 0.3 )
3067                 { /* maximum of background pixels */
3068                   Maximize(max.red,     (double)
3069                     GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3070                   Maximize(max.green,   (double)
3071                     GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3072                   Maximize(max.blue,    (double)
3073                     GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3074                   if (image->colorspace == CMYKColorspace)
3075                     Maximize(max.black,   (double)
3076                       GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3077                   Maximize(max.alpha,(double)
3078                     GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3079                 }
3080               }
3081               k_pixels += virt_width*GetPixelChannels(image);
3082             }
3083             /* Pattern Match if difference is positive */
3084             min.red     -= max.red;     Maximize( min.red,     0.0 );
3085             min.green   -= max.green;   Maximize( min.green,   0.0 );
3086             min.blue    -= max.blue;    Maximize( min.blue,    0.0 );
3087             min.black   -= max.black;   Maximize( min.black,   0.0 );
3088             min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3089             break;
3090
3091         case ErodeIntensityMorphology:
3092             /* Select Pixel with Minimum Intensity within kernel neighbourhood
3093             **
3094             ** WARNING: the intensity test fails for CMYK and does not
3095             ** take into account the moderating effect of the alpha channel
3096             ** on the intensity.
3097             **
3098             ** NOTE that the kernel is not reflected for this operation!
3099             */
3100             k = kernel->values;
3101             k_pixels = p;
3102             for (v=0; v < (ssize_t) kernel->height; v++) {
3103               for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3104                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3105                 if ( result.red == 0.0 ||
3106                      GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) < GetPixelIntensity(morphology_image,q) ) {
3107                   /* copy the whole pixel - no channel selection */
3108                   SetPixelRed(morphology_image,GetPixelRed(image,
3109                     k_pixels+u*GetPixelChannels(image)),q);
3110                   SetPixelGreen(morphology_image,GetPixelGreen(image,
3111                     k_pixels+u*GetPixelChannels(image)),q);
3112                   SetPixelBlue(morphology_image,GetPixelBlue(image,
3113                     k_pixels+u*GetPixelChannels(image)),q);
3114                   SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3115                     k_pixels+u*GetPixelChannels(image)),q);
3116                   if ( result.red > 0.0 ) changed++;
3117                   result.red = 1.0;
3118                 }
3119               }
3120               k_pixels += virt_width*GetPixelChannels(image);
3121             }
3122             break;
3123
3124         case DilateIntensityMorphology:
3125             /* Select Pixel with Maximum Intensity within kernel neighbourhood
3126             **
3127             ** WARNING: the intensity test fails for CMYK and does not
3128             ** take into account the moderating effect of the alpha channel
3129             ** on the intensity (yet).
3130             **
3131             ** NOTE for correct working of this operation for asymetrical
3132             ** kernels, the kernel needs to be applied in its reflected form.
3133             ** That is its values needs to be reversed.
3134             */
3135             k = &kernel->values[ kernel->width*kernel->height-1 ];
3136             k_pixels = p;
3137             for (v=0; v < (ssize_t) kernel->height; v++) {
3138               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3139                 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3140                 if ( result.red == 0.0 ||
3141                      GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) > GetPixelIntensity(morphology_image,q) ) {
3142                   /* copy the whole pixel - no channel selection */
3143                   SetPixelRed(morphology_image,GetPixelRed(image,
3144                     k_pixels+u*GetPixelChannels(image)),q);
3145                   SetPixelGreen(morphology_image,GetPixelGreen(image,
3146                     k_pixels+u*GetPixelChannels(image)),q);
3147                   SetPixelBlue(morphology_image,GetPixelBlue(image,
3148                     k_pixels+u*GetPixelChannels(image)),q);
3149                   SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3150                     k_pixels+u*GetPixelChannels(image)),q);
3151                   if ( result.red > 0.0 ) changed++;
3152                   result.red = 1.0;
3153                 }
3154               }
3155               k_pixels += virt_width*GetPixelChannels(image);
3156             }
3157             break;
3158 #if 0
3159   This code has been obsoleted by the MorphologyPrimitiveDirect() function.
3160   However it is still (almost) correct coding for Grayscale Morphology.
3161   That is...
3162
3163   GrayErode    is equivalent but with kernel values subtracted from pixels
3164                without the kernel rotation
3165   GreyDilate   is equivalent but using Maximum() instead of Minimum()
3166                useing kernel rotation
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 double
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       SetPixelInfo(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                     SetPixelInfo(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                     SetPixelInfo(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           SetPixelPixelInfo(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 double
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       SetPixelInfo(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                     SetPixelInfo(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                     SetPixelInfo(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           SetPixelPixelInfo(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 MagickExport 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) == MagickFalse)
3849         {
3850           InheritException(exception,&rslt_image->exception);
3851           goto error_cleanup;
3852         }
3853
3854       changed = MorphologyPrimitiveDirect(rslt_image, method,
3855          kernel, exception);
3856
3857       if ( verbose == MagickTrue )
3858         (void) (void) FormatLocaleFile(stderr,
3859           "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3860           CommandOptionToMnemonic(MagickMorphologyOptions, method),
3861           1.0,0.0,1.0, (double) changed);
3862
3863       if ( changed < 0 )
3864         goto error_cleanup;
3865
3866       if ( method == VoronoiMorphology ) {
3867         /* Preserve the alpha channel of input image - but turned off */
3868         (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel);
3869         (void) CompositeImage(rslt_image, CopyOpacityCompositeOp, image, 0, 0);
3870         (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel);
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) == MagickFalse)
4029                 {
4030                   InheritException(exception,&work_image->exception);
4031                   goto error_cleanup;
4032                 }
4033               /* work_image->type=image->type; ??? */
4034             }
4035
4036           /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4037           count++;
4038           changed = MorphologyPrimitive(curr_image, work_image, primitive,
4039                        this_kernel, bias, exception);
4040
4041           if ( verbose == MagickTrue ) {
4042             if ( kernel_loop > 1 )
4043               (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4044             (void) (void) FormatLocaleFile(stderr,
4045               "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4046               v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4047               primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4048               (double) (method_loop+kernel_loop-1),(double) kernel_number,
4049               (double) count,(double) changed);
4050           }
4051           if ( changed < 0 )
4052             goto error_cleanup;
4053           kernel_changed += changed;
4054           method_changed += changed;
4055
4056           /* prepare next loop */
4057           { Image *tmp = work_image;   /* swap images for iteration */
4058             work_image = curr_image;
4059             curr_image = tmp;
4060           }
4061           if ( work_image == image )
4062             work_image = (Image *) NULL; /* replace input 'image' */
4063
4064         } /* End Loop 4: Iterate the kernel with primitive */
4065
4066         if ( verbose == MagickTrue && kernel_changed != (size_t)changed )
4067           (void) FormatLocaleFile(stderr, "   Total %.20g",(double) kernel_changed);
4068         if ( verbose == MagickTrue && stage_loop < stage_limit )
4069           (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4070
4071 #if 0
4072     (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4073     (void) FormatLocaleFile(stderr, "      curr =0x%lx\n", (unsigned long)curr_image);
4074     (void) FormatLocaleFile(stderr, "      work =0x%lx\n", (unsigned long)work_image);
4075     (void) FormatLocaleFile(stderr, "      save =0x%lx\n", (unsigned long)save_image);
4076     (void) FormatLocaleFile(stderr, "      union=0x%lx\n", (unsigned long)rslt_image);
4077 #endif
4078
4079       } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4080
4081       /*  Final Post-processing for some Compound Methods
4082       **
4083       ** The removal of any 'Sync' channel flag in the Image Compositon
4084       ** below ensures the methematical compose method is applied in a
4085       ** purely mathematical way, and only to the selected channels.
4086       ** Turn off SVG composition 'alpha blending'.
4087       */
4088       switch( method ) {
4089         case EdgeOutMorphology:
4090         case EdgeInMorphology:
4091         case TopHatMorphology:
4092         case BottomHatMorphology:
4093           if ( verbose == MagickTrue )
4094             (void) FormatLocaleFile(stderr, "\n%s: Difference with original image",
4095                  CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4096           curr_image->sync=MagickFalse;
4097           (void) CompositeImage(curr_image,DifferenceCompositeOp,image,0,0);
4098           break;
4099         case EdgeMorphology:
4100           if ( verbose == MagickTrue )
4101             (void) FormatLocaleFile(stderr, "\n%s: Difference of Dilate and Erode",
4102                  CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4103           curr_image->sync=MagickFalse;
4104           (void) CompositeImage(curr_image,DifferenceCompositeOp,save_image,0,
4105             0);
4106           save_image = DestroyImage(save_image); /* finished with save image */
4107           break;
4108         default:
4109           break;
4110       }
4111
4112       /* multi-kernel handling:  re-iterate, or compose results */
4113       if ( kernel->next == (KernelInfo *) NULL )
4114         rslt_image = curr_image;   /* just return the resulting image */
4115       else if ( rslt_compose == NoCompositeOp )
4116         { if ( verbose == MagickTrue ) {
4117             if ( this_kernel->next != (KernelInfo *) NULL )
4118               (void) FormatLocaleFile(stderr, " (re-iterate)");
4119             else
4120               (void) FormatLocaleFile(stderr, " (done)");
4121           }
4122           rslt_image = curr_image; /* return result, and re-iterate */
4123         }
4124       else if ( rslt_image == (Image *) NULL)
4125         { if ( verbose == MagickTrue )
4126             (void) FormatLocaleFile(stderr, " (save for compose)");
4127           rslt_image = curr_image;
4128           curr_image = (Image *) image;  /* continue with original image */
4129         }
4130       else
4131         { /* Add the new 'current' result to the composition
4132           **
4133           ** The removal of any 'Sync' channel flag in the Image Compositon
4134           ** below ensures the methematical compose method is applied in a
4135           ** purely mathematical way, and only to the selected channels.
4136           ** IE: Turn off SVG composition 'alpha blending'.
4137           */
4138           if ( verbose == MagickTrue )
4139             (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4140                  CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4141           rslt_image->sync=MagickFalse;
4142           (void) CompositeImage(rslt_image, rslt_compose, curr_image, 0, 0);
4143           curr_image = DestroyImage(curr_image);
4144           curr_image = (Image *) image;  /* continue with original image */
4145         }
4146       if ( verbose == MagickTrue )
4147         (void) FormatLocaleFile(stderr, "\n");
4148
4149       /* loop to the next kernel in a multi-kernel list */
4150       norm_kernel = norm_kernel->next;
4151       if ( rflt_kernel != (KernelInfo *) NULL )
4152         rflt_kernel = rflt_kernel->next;
4153       kernel_number++;
4154     } /* End Loop 2: Loop over each kernel */
4155
4156   } /* End Loop 1: compound method interation */
4157
4158   goto exit_cleanup;
4159
4160   /* Yes goto's are bad, but it makes cleanup lot more efficient */
4161 error_cleanup:
4162   if ( curr_image == rslt_image )
4163     curr_image = (Image *) NULL;
4164   if ( rslt_image != (Image *) NULL )
4165     rslt_image = DestroyImage(rslt_image);
4166 exit_cleanup:
4167   if ( curr_image == rslt_image || curr_image == image )
4168     curr_image = (Image *) NULL;
4169   if ( curr_image != (Image *) NULL )
4170     curr_image = DestroyImage(curr_image);
4171   if ( work_image != (Image *) NULL )
4172     work_image = DestroyImage(work_image);
4173   if ( save_image != (Image *) NULL )
4174     save_image = DestroyImage(save_image);
4175   if ( reflected_kernel != (KernelInfo *) NULL )
4176     reflected_kernel = DestroyKernelInfo(reflected_kernel);
4177   return(rslt_image);
4178 }
4179
4180 \f
4181 /*
4182 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4183 %                                                                             %
4184 %                                                                             %
4185 %                                                                             %
4186 %     M o r p h o l o g y I m a g e                                           %
4187 %                                                                             %
4188 %                                                                             %
4189 %                                                                             %
4190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4191 %
4192 %  MorphologyImage() applies a user supplied kernel to the image
4193 %  according to the given mophology method.
4194 %
4195 %  This function applies any and all user defined settings before calling
4196 %  the above internal function MorphologyApply().
4197 %
4198 %  User defined settings include...
4199 %    * Output Bias for Convolution and correlation   ("-bias")
4200 %    * Kernel Scale/normalize settings     ("-set 'option:convolve:scale'")
4201 %      This can also includes the addition of a scaled unity kernel.
4202 %    * Show Kernel being applied           ("-set option:showkernel 1")
4203 %
4204 %  The format of the MorphologyImage method is:
4205 %
4206 %      Image *MorphologyImage(const Image *image,MorphologyMethod method,
4207 %        const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4208 %
4209 %      Image *MorphologyImage(const Image *image, const ChannelType
4210 %        channel,MorphologyMethod method,const ssize_t iterations,
4211 %        KernelInfo *kernel,ExceptionInfo *exception)
4212 %
4213 %  A description of each parameter follows:
4214 %
4215 %    o image: the image.
4216 %
4217 %    o method: the morphology method to be applied.
4218 %
4219 %    o iterations: apply the operation this many times (or no change).
4220 %                  A value of -1 means loop until no change found.
4221 %                  How this is applied may depend on the morphology method.
4222 %                  Typically this is a value of 1.
4223 %
4224 %    o kernel: An array of double representing the morphology kernel.
4225 %              Warning: kernel may be normalized for the Convolve method.
4226 %
4227 %    o exception: return any errors or warnings in this structure.
4228 %
4229 */
4230 MagickExport Image *MorphologyImage(const Image *image,
4231   const MorphologyMethod method,const ssize_t iterations,
4232   const KernelInfo *kernel,ExceptionInfo *exception)
4233 {
4234   KernelInfo
4235     *curr_kernel;
4236
4237   CompositeOperator
4238     compose;
4239
4240   Image
4241     *morphology_image;
4242
4243
4244   /* Apply Convolve/Correlate Normalization and Scaling Factors.
4245    * This is done BEFORE the ShowKernelInfo() function is called so that
4246    * users can see the results of the 'option:convolve:scale' option.
4247    */
4248   curr_kernel = (KernelInfo *) kernel;
4249   if ( method == ConvolveMorphology ||  method == CorrelateMorphology )
4250     {
4251       const char
4252         *artifact;
4253       artifact = GetImageArtifact(image,"convolve:scale");
4254       if ( artifact != (const char *)NULL ) {
4255         if ( curr_kernel == kernel )
4256           curr_kernel = CloneKernelInfo(kernel);
4257         if (curr_kernel == (KernelInfo *) NULL) {
4258           curr_kernel=DestroyKernelInfo(curr_kernel);
4259           return((Image *) NULL);
4260         }
4261         ScaleGeometryKernelInfo(curr_kernel, artifact);
4262       }
4263     }
4264
4265   /* display the (normalized) kernel via stderr */
4266   if ( IsMagickTrue(GetImageArtifact(image,"showkernel"))
4267     || IsMagickTrue(GetImageArtifact(image,"convolve:showkernel"))
4268     || IsMagickTrue(GetImageArtifact(image,"morphology:showkernel")) )
4269     ShowKernelInfo(curr_kernel);
4270
4271   /* Override the default handling of multi-kernel morphology results
4272    * If 'Undefined' use the default method
4273    * If 'None' (default for 'Convolve') re-iterate previous result
4274    * Otherwise merge resulting images using compose method given.
4275    * Default for 'HitAndMiss' is 'Lighten'.
4276    */
4277   { const char
4278       *artifact;
4279     artifact = GetImageArtifact(image,"morphology:compose");
4280     compose = UndefinedCompositeOp;  /* use default for method */
4281     if ( artifact != (const char *) NULL)
4282       compose=(CompositeOperator) ParseCommandOption(MagickComposeOptions,
4283         MagickFalse,artifact);
4284   }
4285   /* Apply the Morphology */
4286   morphology_image = MorphologyApply(image, method, iterations,
4287     curr_kernel, compose, image->bias, exception);
4288
4289   /* Cleanup and Exit */
4290   if ( curr_kernel != kernel )
4291     curr_kernel=DestroyKernelInfo(curr_kernel);
4292   return(morphology_image);
4293 }
4294 \f
4295 /*
4296 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4297 %                                                                             %
4298 %                                                                             %
4299 %                                                                             %
4300 +     R o t a t e K e r n e l I n f o                                         %
4301 %                                                                             %
4302 %                                                                             %
4303 %                                                                             %
4304 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4305 %
4306 %  RotateKernelInfo() rotates the kernel by the angle given.
4307 %
4308 %  Currently it is restricted to 90 degree angles, of either 1D kernels
4309 %  or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4310 %  It will ignore usless rotations for specific 'named' built-in kernels.
4311 %
4312 %  The format of the RotateKernelInfo method is:
4313 %
4314 %      void RotateKernelInfo(KernelInfo *kernel, double angle)
4315 %
4316 %  A description of each parameter follows:
4317 %
4318 %    o kernel: the Morphology/Convolution kernel
4319 %
4320 %    o angle: angle to rotate in degrees
4321 %
4322 % This function is currently internal to this module only, but can be exported
4323 % to other modules if needed.
4324 */
4325 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4326 {
4327   /* angle the lower kernels first */
4328   if ( kernel->next != (KernelInfo *) NULL)
4329     RotateKernelInfo(kernel->next, angle);
4330
4331   /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4332   **
4333   ** TODO: expand beyond simple 90 degree rotates, flips and flops
4334   */
4335
4336   /* Modulus the angle */
4337   angle = fmod(angle, 360.0);
4338   if ( angle < 0 )
4339     angle += 360.0;
4340
4341   if ( 337.5 < angle || angle <= 22.5 )
4342     return;   /* Near zero angle - no change! - At least not at this time */
4343
4344   /* Handle special cases */
4345   switch (kernel->type) {
4346     /* These built-in kernels are cylindrical kernels, rotating is useless */
4347     case GaussianKernel:
4348     case DoGKernel:
4349     case LoGKernel:
4350     case DiskKernel:
4351     case PeaksKernel:
4352     case LaplacianKernel:
4353     case ChebyshevKernel:
4354     case ManhattanKernel:
4355     case EuclideanKernel:
4356       return;
4357
4358     /* These may be rotatable at non-90 angles in the future */
4359     /* but simply rotating them in multiples of 90 degrees is useless */
4360     case SquareKernel:
4361     case DiamondKernel:
4362     case PlusKernel:
4363     case CrossKernel:
4364       return;
4365
4366     /* These only allows a +/-90 degree rotation (by transpose) */
4367     /* A 180 degree rotation is useless */
4368     case BlurKernel:
4369     case RectangleKernel:
4370       if ( 135.0 < angle && angle <= 225.0 )
4371         return;
4372       if ( 225.0 < angle && angle <= 315.0 )
4373         angle -= 180;
4374       break;
4375
4376     default:
4377       break;
4378   }
4379   /* Attempt rotations by 45 degrees */
4380   if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4381     {
4382       if ( kernel->width == 3 && kernel->height == 3 )
4383         { /* Rotate a 3x3 square by 45 degree angle */
4384           MagickRealType t  = kernel->values[0];
4385           kernel->values[0] = kernel->values[3];
4386           kernel->values[3] = kernel->values[6];
4387           kernel->values[6] = kernel->values[7];
4388           kernel->values[7] = kernel->values[8];
4389           kernel->values[8] = kernel->values[5];
4390           kernel->values[5] = kernel->values[2];
4391           kernel->values[2] = kernel->values[1];
4392           kernel->values[1] = t;
4393           /* rotate non-centered origin */
4394           if ( kernel->x != 1 || kernel->y != 1 ) {
4395             ssize_t x,y;
4396             x = (ssize_t) kernel->x-1;
4397             y = (ssize_t) kernel->y-1;
4398                  if ( x == y  ) x = 0;
4399             else if ( x == 0  ) x = -y;
4400             else if ( x == -y ) y = 0;
4401             else if ( y == 0  ) y = x;
4402             kernel->x = (ssize_t) x+1;
4403             kernel->y = (ssize_t) y+1;
4404           }
4405           angle = fmod(angle+315.0, 360.0);  /* angle reduced 45 degrees */
4406           kernel->angle = fmod(kernel->angle+45.0, 360.0);
4407         }
4408       else
4409         perror("Unable to rotate non-3x3 kernel by 45 degrees");
4410     }
4411   if ( 45.0 < fmod(angle, 180.0)  && fmod(angle,180.0) <= 135.0 )
4412     {
4413       if ( kernel->width == 1 || kernel->height == 1 )
4414         { /* Do a transpose of a 1 dimensional kernel,
4415           ** which results in a fast 90 degree rotation of some type.
4416           */
4417           ssize_t
4418             t;
4419           t = (ssize_t) kernel->width;
4420           kernel->width = kernel->height;
4421           kernel->height = (size_t) t;
4422           t = kernel->x;
4423           kernel->x = kernel->y;
4424           kernel->y = t;
4425           if ( kernel->width == 1 ) {
4426             angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
4427             kernel->angle = fmod(kernel->angle+90.0, 360.0);
4428           } else {
4429             angle = fmod(angle+90.0, 360.0);   /* angle increased 90 degrees */
4430             kernel->angle = fmod(kernel->angle+270.0, 360.0);
4431           }
4432         }
4433       else if ( kernel->width == kernel->height )
4434         { /* Rotate a square array of values by 90 degrees */
4435           { register size_t
4436               i,j,x,y;
4437             register MagickRealType
4438               *k,t;
4439             k=kernel->values;
4440             for( i=0, x=kernel->width-1;  i<=x;   i++, x--)
4441               for( j=0, y=kernel->height-1;  j<y;   j++, y--)
4442                 { t                    = k[i+j*kernel->width];
4443                   k[i+j*kernel->width] = k[j+x*kernel->width];
4444                   k[j+x*kernel->width] = k[x+y*kernel->width];
4445                   k[x+y*kernel->width] = k[y+i*kernel->width];
4446                   k[y+i*kernel->width] = t;
4447                 }
4448           }
4449           /* rotate the origin - relative to center of array */
4450           { register ssize_t x,y;
4451             x = (ssize_t) (kernel->x*2-kernel->width+1);
4452             y = (ssize_t) (kernel->y*2-kernel->height+1);
4453             kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4454             kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4455           }
4456           angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
4457           kernel->angle = fmod(kernel->angle+90.0, 360.0);
4458         }
4459       else
4460         perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4461     }
4462   if ( 135.0 < angle && angle <= 225.0 )
4463     {
4464       /* For a 180 degree rotation - also know as a reflection
4465        * This is actually a very very common operation!
4466        * Basically all that is needed is a reversal of the kernel data!
4467        * And a reflection of the origon
4468        */
4469       size_t
4470         i,j;
4471       register double
4472         *k,t;
4473
4474       k=kernel->values;
4475       for ( i=0, j=kernel->width*kernel->height-1;  i<j;  i++, j--)
4476         t=k[i],  k[i]=k[j],  k[j]=t;
4477
4478       kernel->x = (ssize_t) kernel->width  - kernel->x - 1;
4479       kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4480       angle = fmod(angle-180.0, 360.0);   /* angle+180 degrees */
4481       kernel->angle = fmod(kernel->angle+180.0, 360.0);
4482     }
4483   /* At this point angle should at least between -45 (315) and +45 degrees
4484    * In the future some form of non-orthogonal angled rotates could be
4485    * performed here, posibily with a linear kernel restriction.
4486    */
4487
4488   return;
4489 }
4490 \f
4491 /*
4492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4493 %                                                                             %
4494 %                                                                             %
4495 %                                                                             %
4496 %     S c a l e G e o m e t r y K e r n e l I n f o                           %
4497 %                                                                             %
4498 %                                                                             %
4499 %                                                                             %
4500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4501 %
4502 %  ScaleGeometryKernelInfo() takes a geometry argument string, typically
4503 %  provided as a  "-set option:convolve:scale {geometry}" user setting,
4504 %  and modifies the kernel according to the parsed arguments of that setting.
4505 %
4506 %  The first argument (and any normalization flags) are passed to
4507 %  ScaleKernelInfo() to scale/normalize the kernel.  The second argument
4508 %  is then passed to UnityAddKernelInfo() to add a scled unity kernel
4509 %  into the scaled/normalized kernel.
4510 %
4511 %  The format of the ScaleGeometryKernelInfo method is:
4512 %
4513 %      void ScaleGeometryKernelInfo(KernelInfo *kernel,
4514 %        const double scaling_factor,const MagickStatusType normalize_flags)
4515 %
4516 %  A description of each parameter follows:
4517 %
4518 %    o kernel: the Morphology/Convolution kernel to modify
4519 %
4520 %    o geometry:
4521 %             The geometry string to parse, typically from the user provided
4522 %             "-set option:convolve:scale {geometry}" setting.
4523 %
4524 */
4525 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4526      const char *geometry)
4527 {
4528   GeometryFlags
4529     flags;
4530   GeometryInfo
4531     args;
4532
4533   SetGeometryInfo(&args);
4534   flags = (GeometryFlags) ParseGeometry(geometry, &args);
4535
4536 #if 0
4537   /* For Debugging Geometry Input */
4538   (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4539        flags, args.rho, args.sigma, args.xi, args.psi );
4540 #endif
4541
4542   if ( (flags & PercentValue) != 0 )      /* Handle Percentage flag*/
4543     args.rho *= 0.01,  args.sigma *= 0.01;
4544
4545   if ( (flags & RhoValue) == 0 )          /* Set Defaults for missing args */
4546     args.rho = 1.0;
4547   if ( (flags & SigmaValue) == 0 )
4548     args.sigma = 0.0;
4549
4550   /* Scale/Normalize the input kernel */
4551   ScaleKernelInfo(kernel, args.rho, flags);
4552
4553   /* Add Unity Kernel, for blending with original */
4554   if ( (flags & SigmaValue) != 0 )
4555     UnityAddKernelInfo(kernel, args.sigma);
4556
4557   return;
4558 }
4559 /*
4560 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4561 %                                                                             %
4562 %                                                                             %
4563 %                                                                             %
4564 %     S c a l e K e r n e l I n f o                                           %
4565 %                                                                             %
4566 %                                                                             %
4567 %                                                                             %
4568 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4569 %
4570 %  ScaleKernelInfo() scales the given kernel list by the given amount, with or
4571 %  without normalization of the sum of the kernel values (as per given flags).
4572 %
4573 %  By default (no flags given) the values within the kernel is scaled
4574 %  directly using given scaling factor without change.
4575 %
4576 %  If either of the two 'normalize_flags' are given the kernel will first be
4577 %  normalized and then further scaled by the scaling factor value given.
4578 %
4579 %  Kernel normalization ('normalize_flags' given) is designed to ensure that
4580 %  any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4581 %  morphology methods will fall into -1.0 to +1.0 range.  Note that for
4582 %  non-HDRI versions of IM this may cause images to have any negative results
4583 %  clipped, unless some 'bias' is used.
4584 %
4585 %  More specifically.  Kernels which only contain positive values (such as a
4586 %  'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4587 %  ensuring a 0.0 to +1.0 output range for non-HDRI images.
4588 %
4589 %  For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4590 %  the kernel will be scaled by the absolute of the sum of kernel values, so
4591 %  that it will generally fall within the +/- 1.0 range.
4592 %
4593 %  For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4594 %  will be scaled by just the sum of the postive values, so that its output
4595 %  range will again fall into the  +/- 1.0 range.
4596 %
4597 %  For special kernels designed for locating shapes using 'Correlate', (often
4598 %  only containing +1 and -1 values, representing foreground/brackground
4599 %  matching) a special normalization method is provided to scale the positive
4600 %  values separately to those of the negative values, so the kernel will be
4601 %  forced to become a zero-sum kernel better suited to such searches.
4602 %
4603 %  WARNING: Correct normalization of the kernel assumes that the '*_range'
4604 %  attributes within the kernel structure have been correctly set during the
4605 %  kernels creation.
4606 %
4607 %  NOTE: The values used for 'normalize_flags' have been selected specifically
4608 %  to match the use of geometry options, so that '!' means NormalizeValue, '^'
4609 %  means CorrelateNormalizeValue.  All other GeometryFlags values are ignored.
4610 %
4611 %  The format of the ScaleKernelInfo method is:
4612 %
4613 %      void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4614 %               const MagickStatusType normalize_flags )
4615 %
4616 %  A description of each parameter follows:
4617 %
4618 %    o kernel: the Morphology/Convolution kernel
4619 %
4620 %    o scaling_factor:
4621 %             multiply all values (after normalization) by this factor if not
4622 %             zero.  If the kernel is normalized regardless of any flags.
4623 %
4624 %    o normalize_flags:
4625 %             GeometryFlags defining normalization method to use.
4626 %             specifically: NormalizeValue, CorrelateNormalizeValue,
4627 %                           and/or PercentValue
4628 %
4629 */
4630 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4631   const double scaling_factor,const GeometryFlags normalize_flags)
4632 {
4633   register ssize_t
4634     i;
4635
4636   register double
4637     pos_scale,
4638     neg_scale;
4639
4640   /* do the other kernels in a multi-kernel list first */
4641   if ( kernel->next != (KernelInfo *) NULL)
4642     ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4643
4644   /* Normalization of Kernel */
4645   pos_scale = 1.0;
4646   if ( (normalize_flags&NormalizeValue) != 0 ) {
4647     if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4648       /* non-zero-summing kernel (generally positive) */
4649       pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4650     else
4651       /* zero-summing kernel */
4652       pos_scale = kernel->positive_range;
4653   }
4654   /* Force kernel into a normalized zero-summing kernel */
4655   if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4656     pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4657                  ? kernel->positive_range : 1.0;
4658     neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
4659                  ? -kernel->negative_range : 1.0;
4660   }
4661   else
4662     neg_scale = pos_scale;
4663
4664   /* finialize scaling_factor for positive and negative components */
4665   pos_scale = scaling_factor/pos_scale;
4666   neg_scale = scaling_factor/neg_scale;
4667
4668   for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4669     if ( ! IsNan(kernel->values[i]) )
4670       kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4671
4672   /* convolution output range */
4673   kernel->positive_range *= pos_scale;
4674   kernel->negative_range *= neg_scale;
4675   /* maximum and minimum values in kernel */
4676   kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4677   kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4678
4679   /* swap kernel settings if user's scaling factor is negative */
4680   if ( scaling_factor < MagickEpsilon ) {
4681     double t;
4682     t = kernel->positive_range;
4683     kernel->positive_range = kernel->negative_range;
4684     kernel->negative_range = t;
4685     t = kernel->maximum;
4686     kernel->maximum = kernel->minimum;
4687     kernel->minimum = 1;
4688   }
4689
4690   return;
4691 }
4692 \f
4693 /*
4694 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4695 %                                                                             %
4696 %                                                                             %
4697 %                                                                             %
4698 %     S h o w K e r n e l I n f o                                             %
4699 %                                                                             %
4700 %                                                                             %
4701 %                                                                             %
4702 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4703 %
4704 %  ShowKernelInfo() outputs the details of the given kernel defination to
4705 %  standard error, generally due to a users 'showkernel' option request.
4706 %
4707 %  The format of the ShowKernel method is:
4708 %
4709 %      void ShowKernelInfo(KernelInfo *kernel)
4710 %
4711 %  A description of each parameter follows:
4712 %
4713 %    o kernel: the Morphology/Convolution kernel
4714 %
4715 */
4716 MagickExport void ShowKernelInfo(KernelInfo *kernel)
4717 {
4718   KernelInfo
4719     *k;
4720
4721   size_t
4722     c, i, u, v;
4723
4724   for (c=0, k=kernel;  k != (KernelInfo *) NULL;  c++, k=k->next ) {
4725
4726     (void) FormatLocaleFile(stderr, "Kernel");
4727     if ( kernel->next != (KernelInfo *) NULL )
4728       (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4729     (void) FormatLocaleFile(stderr, " \"%s",
4730           CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4731     if ( fabs(k->angle) > MagickEpsilon )
4732       (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4733     (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4734       k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4735     (void) FormatLocaleFile(stderr,
4736           " with values from %.*lg to %.*lg\n",
4737           GetMagickPrecision(), k->minimum,
4738           GetMagickPrecision(), k->maximum);
4739     (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4740           GetMagickPrecision(), k->negative_range,
4741           GetMagickPrecision(), k->positive_range);
4742     if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4743       (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4744     else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4745       (void) FormatLocaleFile(stderr, " (Normalized)\n");
4746     else
4747       (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4748           GetMagickPrecision(), k->positive_range+k->negative_range);
4749     for (i=v=0; v < k->height; v++) {
4750       (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4751       for (u=0; u < k->width; u++, i++)
4752         if ( IsNan(k->values[i]) )
4753           (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4754         else
4755           (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4756               GetMagickPrecision(), k->values[i]);
4757       (void) FormatLocaleFile(stderr,"\n");
4758     }
4759   }
4760 }
4761 \f
4762 /*
4763 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4764 %                                                                             %
4765 %                                                                             %
4766 %                                                                             %
4767 %     U n i t y A d d K e r n a l I n f o                                     %
4768 %                                                                             %
4769 %                                                                             %
4770 %                                                                             %
4771 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4772 %
4773 %  UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4774 %  to the given pre-scaled and normalized Kernel.  This in effect adds that
4775 %  amount of the original image into the resulting convolution kernel.  This
4776 %  value is usually provided by the user as a percentage value in the
4777 %  'convolve:scale' setting.
4778 %
4779 %  The resulting effect is to convert the defined kernels into blended
4780 %  soft-blurs, unsharp kernels or into sharpening kernels.
4781 %
4782 %  The format of the UnityAdditionKernelInfo method is:
4783 %
4784 %      void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4785 %
4786 %  A description of each parameter follows:
4787 %
4788 %    o kernel: the Morphology/Convolution kernel
4789 %
4790 %    o scale:
4791 %             scaling factor for the unity kernel to be added to
4792 %             the given kernel.
4793 %
4794 */
4795 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4796   const double scale)
4797 {
4798   /* do the other kernels in a multi-kernel list first */
4799   if ( kernel->next != (KernelInfo *) NULL)
4800     UnityAddKernelInfo(kernel->next, scale);
4801
4802   /* Add the scaled unity kernel to the existing kernel */
4803   kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4804   CalcKernelMetaData(kernel);  /* recalculate the meta-data */
4805
4806   return;
4807 }
4808 \f
4809 /*
4810 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4811 %                                                                             %
4812 %                                                                             %
4813 %                                                                             %
4814 %     Z e r o K e r n e l N a n s                                             %
4815 %                                                                             %
4816 %                                                                             %
4817 %                                                                             %
4818 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4819 %
4820 %  ZeroKernelNans() replaces any special 'nan' value that may be present in
4821 %  the kernel with a zero value.  This is typically done when the kernel will
4822 %  be used in special hardware (GPU) convolution processors, to simply
4823 %  matters.
4824 %
4825 %  The format of the ZeroKernelNans method is:
4826 %
4827 %      void ZeroKernelNans (KernelInfo *kernel)
4828 %
4829 %  A description of each parameter follows:
4830 %
4831 %    o kernel: the Morphology/Convolution kernel
4832 %
4833 */
4834 MagickExport void ZeroKernelNans(KernelInfo *kernel)
4835 {
4836   register size_t
4837     i;
4838
4839   /* do the other kernels in a multi-kernel list first */
4840   if ( kernel->next != (KernelInfo *) NULL)
4841     ZeroKernelNans(kernel->next);
4842
4843   for (i=0; i < (kernel->width*kernel->height); i++)
4844     if ( IsNan(kernel->values[i]) )
4845       kernel->values[i] = 0.0;
4846
4847   return;
4848 }