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

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