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

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