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

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