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