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

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