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

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %    M   M    OOO    RRRR   PPPP   H   H   OOO   L       OOO    GGGG  Y   Y   %
   7 %    MM MM   O   O   R   R  P   P  H   H  O   O  L      O   O  G       Y Y    %
   8 %    M M M   O   O   RRRR   PPPP   HHHHH  O   O  L      O   O  G GGG    Y     %
   9 %    M   M   O   O   R R    P      H   H  O   O  L      O   O  G   G    Y     %
  10 %    M   M    OOO    R  R   P      H   H   OOO   LLLLL   OOO    GGG     Y     %
  11 %                                                                             %
  12 %                                                                             %
  13 %                        MagickCore Morphology Methods                        %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                              Anthony Thyssen                                %
  17 %                               January 2010                                  %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 % Morpology is the the application of various kernels, of any size and even
  37 % shape, to a image in various ways (typically binary, but not always).
  38 %
  39 % Convolution (weighted sum or average) is just one specific type of
  40 % morphology. Just one that is very common for image bluring and sharpening
  41 % effects.  Not only 2D Gaussian blurring, but also 2-pass 1D Blurring.
  42 %
  43 % This module provides not only a general morphology function, and the ability
  44 % to apply more advanced or iterative morphologies, but also functions for the
  45 % generation of many different types of kernel arrays from user supplied
  46 % arguments. Prehaps even the generation of a kernel from a small image.
  47 */
  48 \f
  49 /*
  50   Include declarations.
  51 */
  52 #include "MagickCore/studio.h"
  53 #include "MagickCore/artifact.h"
  54 #include "MagickCore/cache-view.h"
  55 #include "MagickCore/color-private.h"
  56 #include "MagickCore/enhance.h"
  57 #include "MagickCore/exception.h"
  58 #include "MagickCore/exception-private.h"
  59 #include "MagickCore/gem.h"
  60 #include "MagickCore/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 *) AcquireMagickMemory(sizeof(*kernel));
 244   if (kernel == (KernelInfo *)NULL)
 245     return(kernel);
 246   (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
 247   kernel->minimum = kernel->maximum = kernel->angle = 0.0;
 248   kernel->negative_range = kernel->positive_range = 0.0;
 249   kernel->type = UserDefinedKernel;
 250   kernel->next = (KernelInfo *) NULL;
 251   kernel->signature = MagickSignature;
 252   if (kernel_string == (const char *) NULL)
 253     return(kernel);
 254
 255   /* find end of this specific kernel definition string */
 256   end = strchr(kernel_string, ';');
 257   if ( end == (char *) NULL )
 258     end = strchr(kernel_string, '\0');
 259
 260   /* clear flags - for Expanding kernel lists thorugh rotations */
 261    flags = NoValue;
 262
 263   /* Has a ':' in argument - New user kernel specification */
 264   p = strchr(kernel_string, ':');
 265   if ( p != (char *) NULL && p < end)
 266     {
 267       /* ParseGeometry() needs the geometry separated! -- Arrgghh */
 268       memcpy(token, kernel_string, (size_t) (p-kernel_string));
 269       token[p-kernel_string] = '\0';
 270       SetGeometryInfo(&args);
 271       flags = ParseGeometry(token, &args);
 272
 273       /* Size handling and checks of geometry settings */
 274       if ( (flags & WidthValue) == 0 ) /* if no width then */
 275         args.rho = args.sigma;         /* then  width = height */
 276       if ( args.rho < 1.0 )            /* if width too small */
 277          args.rho = 1.0;               /* then  width = 1 */
 278       if ( args.sigma < 1.0 )          /* if height too small */
 279         args.sigma = args.rho;         /* then  height = width */
 280       kernel->width = (size_t)args.rho;
 281       kernel->height = (size_t)args.sigma;
 282
 283       /* Offset Handling and Checks */
 284       if ( args.xi  < 0.0 || args.psi < 0.0 )
 285         return(DestroyKernelInfo(kernel));
 286       kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
 287                                                : (ssize_t) (kernel->width-1)/2;
 288       kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
 289                                                : (ssize_t) (kernel->height-1)/2;
 290       if ( kernel->x >= (ssize_t) kernel->width ||
 291            kernel->y >= (ssize_t) kernel->height )
 292         return(DestroyKernelInfo(kernel));
 293
 294       p++; /* advance beyond the ':' */
 295     }
 296   else
 297     { /* ELSE - Old old specification, forming odd-square kernel */
 298       /* count up number of values given */
 299       p=(const char *) kernel_string;
 300       while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
 301         p++;  /* ignore "'" chars for convolve filter usage - Cristy */
 302       for (i=0; p < end; i++)
 303       {
 304         GetMagickToken(p,&p,token);
 305         if (*token == ',')
 306           GetMagickToken(p,&p,token);
 307       }
 308       /* set the size of the kernel - old sized square */
 309       kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
 310       kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
 311       p=(const char *) kernel_string;
 312       while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
 313         p++;  /* ignore "'" chars for convolve filter usage - Cristy */
 314     }
 315
 316   /* Read in the kernel values from rest of input string argument */
 317   kernel->values=(double *) AcquireAlignedMemory(kernel->width,
 318                         kernel->height*sizeof(double));
 319   if (kernel->values == (double *) NULL)
 320     return(DestroyKernelInfo(kernel));
 321
 322   kernel->minimum = +MagickHuge;
 323   kernel->maximum = -MagickHuge;
 324   kernel->negative_range = kernel->positive_range = 0.0;
 325
 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; /* do not include this value in kernel */
 334     }
 335     else {
 336       kernel->values[i] = InterpretLocaleValue(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=(double *) AcquireAlignedMemory(1,sizeof(double));
1034         if (kernel->values == (double *) NULL)
1035           return(DestroyKernelInfo(kernel));
1036         kernel->maximum = kernel->values[0] = args->rho;
1037         break;
1038       }
1039       break;
1040     case GaussianKernel:
1041     case DoGKernel:
1042     case LoGKernel:
1043       { double
1044           sigma = fabs(args->sigma),
1045           sigma2 = fabs(args->xi),
1046           A, B, R;
1047
1048         if ( args->rho >= 1.0 )
1049           kernel->width = (size_t)args->rho*2+1;
1050         else if ( (type != DoGKernel) || (sigma >= sigma2) )
1051           kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1052         else
1053           kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1054         kernel->height = kernel->width;
1055         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1056         kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1057                               kernel->height*sizeof(double));
1058         if (kernel->values == (double *) NULL)
1059           return(DestroyKernelInfo(kernel));
1060
1061         /* WARNING: The following generates a 'sampled gaussian' kernel.
1062          * What we really want is a 'discrete gaussian' kernel.
1063          *
1064          * How to do this is I don't know, but appears to be basied on the
1065          * Error Function 'erf()' (intergral of a gaussian)
1066          */
1067
1068         if ( type == GaussianKernel || type == DoGKernel )
1069           { /* Calculate a Gaussian,  OR positive half of a DoG */
1070             if ( sigma > MagickEpsilon )
1071               { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
1072                 B = (double) (1.0/(Magick2PI*sigma*sigma));
1073                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1074                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1075                       kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1076               }
1077             else /* limiting case - a unity (normalized Dirac) kernel */
1078               { (void) ResetMagickMemory(kernel->values,0, (size_t)
1079                             kernel->width*kernel->height*sizeof(double));
1080                 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1081               }
1082           }
1083
1084         if ( type == DoGKernel )
1085           { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1086             if ( sigma2 > MagickEpsilon )
1087               { sigma = sigma2;                /* simplify loop expressions */
1088                 A = 1.0/(2.0*sigma*sigma);
1089                 B = (double) (1.0/(Magick2PI*sigma*sigma));
1090                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1091                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1092                     kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1093               }
1094             else /* limiting case - a unity (normalized Dirac) kernel */
1095               kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1096           }
1097
1098         if ( type == LoGKernel )
1099           { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1100             if ( sigma > MagickEpsilon )
1101               { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
1102                 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1103                 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1104                   for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1105                     { R = ((double)(u*u+v*v))*A;
1106                       kernel->values[i] = (1-R)*exp(-R)*B;
1107                     }
1108               }
1109             else /* special case - generate a unity kernel */
1110               { (void) ResetMagickMemory(kernel->values,0, (size_t)
1111                             kernel->width*kernel->height*sizeof(double));
1112                 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1113               }
1114           }
1115
1116         /* Note the above kernels may have been 'clipped' by a user defined
1117         ** radius, producing a smaller (darker) kernel.  Also for very small
1118         ** sigma's (> 0.1) the central value becomes larger than one, and thus
1119         ** producing a very bright kernel.
1120         **
1121         ** Normalization will still be needed.
1122         */
1123
1124         /* Normalize the 2D Gaussian Kernel
1125         **
1126         ** NB: a CorrelateNormalize performs a normal Normalize if
1127         ** there are no negative values.
1128         */
1129         CalcKernelMetaData(kernel);  /* the other kernel meta-data */
1130         ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1131
1132         break;
1133       }
1134     case BlurKernel:
1135       { double
1136           sigma = fabs(args->sigma),
1137           alpha, beta;
1138
1139         if ( args->rho >= 1.0 )
1140           kernel->width = (size_t)args->rho*2+1;
1141         else
1142           kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1143         kernel->height = 1;
1144         kernel->x = (ssize_t) (kernel->width-1)/2;
1145         kernel->y = 0;
1146         kernel->negative_range = kernel->positive_range = 0.0;
1147         kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1148                               kernel->height*sizeof(double));
1149         if (kernel->values == (double *) NULL)
1150           return(DestroyKernelInfo(kernel));
1151
1152 #if 1
1153 #define KernelRank 3
1154         /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1155         ** It generates a gaussian 3 times the width, and compresses it into
1156         ** the expected range.  This produces a closer normalization of the
1157         ** resulting kernel, especially for very low sigma values.
1158         ** As such while wierd it is prefered.
1159         **
1160         ** I am told this method originally came from Photoshop.
1161         **
1162         ** A properly normalized curve is generated (apart from edge clipping)
1163         ** even though we later normalize the result (for edge clipping)
1164         ** to allow the correct generation of a "Difference of Blurs".
1165         */
1166
1167         /* initialize */
1168         v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1169         (void) ResetMagickMemory(kernel->values,0, (size_t)
1170                      kernel->width*kernel->height*sizeof(double));
1171         /* Calculate a Positive 1D Gaussian */
1172         if ( sigma > MagickEpsilon )
1173           { sigma *= KernelRank;               /* simplify loop expressions */
1174             alpha = 1.0/(2.0*sigma*sigma);
1175             beta= (double) (1.0/(MagickSQ2PI*sigma ));
1176             for ( u=-v; u <= v; u++) {
1177               kernel->values[(u+v)/KernelRank] +=
1178                               exp(-((double)(u*u))*alpha)*beta;
1179             }
1180           }
1181         else /* special case - generate a unity kernel */
1182           kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1183 #else
1184         /* Direct calculation without curve averaging */
1185
1186         /* Calculate a Positive Gaussian */
1187         if ( sigma > MagickEpsilon )
1188           { alpha = 1.0/(2.0*sigma*sigma);    /* simplify loop expressions */
1189             beta = 1.0/(MagickSQ2PI*sigma);
1190             for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1191               kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1192           }
1193         else /* special case - generate a unity kernel */
1194           { (void) ResetMagickMemory(kernel->values,0, (size_t)
1195                          kernel->width*kernel->height*sizeof(double));
1196             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1197           }
1198 #endif
1199         /* Note the above kernel may have been 'clipped' by a user defined
1200         ** radius, producing a smaller (darker) kernel.  Also for very small
1201         ** sigma's (> 0.1) the central value becomes larger than one, and thus
1202         ** producing a very bright kernel.
1203         **
1204         ** Normalization will still be needed.
1205         */
1206
1207         /* Normalize the 1D Gaussian Kernel
1208         **
1209         ** NB: a CorrelateNormalize performs a normal Normalize if
1210         ** there are no negative values.
1211         */
1212         CalcKernelMetaData(kernel);  /* the other kernel meta-data */
1213         ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1214
1215         /* rotate the 1D kernel by given angle */
1216         RotateKernelInfo(kernel, args->xi );
1217         break;
1218       }
1219     case CometKernel:
1220       { double
1221           sigma = fabs(args->sigma),
1222           A;
1223
1224         if ( args->rho < 1.0 )
1225           kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1226         else
1227           kernel->width = (size_t)args->rho;
1228         kernel->x = kernel->y = 0;
1229         kernel->height = 1;
1230         kernel->negative_range = kernel->positive_range = 0.0;
1231         kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1232                               kernel->height*sizeof(double));
1233         if (kernel->values == (double *) NULL)
1234           return(DestroyKernelInfo(kernel));
1235
1236         /* A comet blur is half a 1D gaussian curve, so that the object is
1237         ** blurred in one direction only.  This may not be quite the right
1238         ** curve to use so may change in the future. The function must be
1239         ** normalised after generation, which also resolves any clipping.
1240         **
1241         ** As we are normalizing and not subtracting gaussians,
1242         ** there is no need for a divisor in the gaussian formula
1243         **
1244         ** It is less comples
1245         */
1246         if ( sigma > MagickEpsilon )
1247           {
1248 #if 1
1249 #define KernelRank 3
1250             v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1251             (void) ResetMagickMemory(kernel->values,0, (size_t)
1252                           kernel->width*sizeof(double));
1253             sigma *= KernelRank;            /* simplify the loop expression */
1254             A = 1.0/(2.0*sigma*sigma);
1255             /* B = 1.0/(MagickSQ2PI*sigma); */
1256             for ( u=0; u < v; u++) {
1257               kernel->values[u/KernelRank] +=
1258                   exp(-((double)(u*u))*A);
1259               /*  exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1260             }
1261             for (i=0; i < (ssize_t) kernel->width; i++)
1262               kernel->positive_range += kernel->values[i];
1263 #else
1264             A = 1.0/(2.0*sigma*sigma);     /* simplify the loop expression */
1265             /* B = 1.0/(MagickSQ2PI*sigma); */
1266             for ( i=0; i < (ssize_t) kernel->width; i++)
1267               kernel->positive_range +=
1268                 kernel->values[i] =
1269                   exp(-((double)(i*i))*A);
1270                 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1271 #endif
1272           }
1273         else /* special case - generate a unity kernel */
1274           { (void) ResetMagickMemory(kernel->values,0, (size_t)
1275                          kernel->width*kernel->height*sizeof(double));
1276             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1277             kernel->positive_range = 1.0;
1278           }
1279
1280         kernel->minimum = 0.0;
1281         kernel->maximum = kernel->values[0];
1282         kernel->negative_range = 0.0;
1283
1284         ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1285         RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1286         break;
1287       }
1288
1289     /*
1290       Convolution Kernels - Well Known Named Constant Kernels
1291     */
1292     case LaplacianKernel:
1293       { switch ( (int) args->rho ) {
1294           case 0:
1295           default: /* laplacian square filter -- default */
1296             kernel=ParseKernelArray("3: -1,-1,-1  -1,8,-1  -1,-1,-1");
1297             break;
1298           case 1:  /* laplacian diamond filter */
1299             kernel=ParseKernelArray("3: 0,-1,0  -1,4,-1  0,-1,0");
1300             break;
1301           case 2:
1302             kernel=ParseKernelArray("3: -2,1,-2  1,4,1  -2,1,-2");
1303             break;
1304           case 3:
1305             kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  1,-2,1");
1306             break;
1307           case 5:   /* a 5x5 laplacian */
1308             kernel=ParseKernelArray(
1309               "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");
1310             break;
1311           case 7:   /* a 7x7 laplacian */
1312             kernel=ParseKernelArray(
1313               "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" );
1314             break;
1315           case 15:  /* a 5x5 LoG (sigma approx 1.4) */
1316             kernel=ParseKernelArray(
1317               "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");
1318             break;
1319           case 19:  /* a 9x9 LoG (sigma approx 1.4) */
1320             /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1321             kernel=ParseKernelArray(
1322               "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");
1323             break;
1324         }
1325         if (kernel == (KernelInfo *) NULL)
1326           return(kernel);
1327         kernel->type = type;
1328         break;
1329       }
1330     case SobelKernel:
1331       { /* Simple Sobel Kernel */
1332         kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1333         if (kernel == (KernelInfo *) NULL)
1334           return(kernel);
1335         kernel->type = type;
1336         RotateKernelInfo(kernel, args->rho);
1337         break;
1338       }
1339     case RobertsKernel:
1340       {
1341         kernel=ParseKernelArray("3: 0,0,0  1,-1,0  0,0,0");
1342         if (kernel == (KernelInfo *) NULL)
1343           return(kernel);
1344         kernel->type = type;
1345         RotateKernelInfo(kernel, args->rho);
1346         break;
1347       }
1348     case PrewittKernel:
1349       {
1350         kernel=ParseKernelArray("3: 1,0,-1  1,0,-1  1,0,-1");
1351         if (kernel == (KernelInfo *) NULL)
1352           return(kernel);
1353         kernel->type = type;
1354         RotateKernelInfo(kernel, args->rho);
1355         break;
1356       }
1357     case CompassKernel:
1358       {
1359         kernel=ParseKernelArray("3: 1,1,-1  1,-2,-1  1,1,-1");
1360         if (kernel == (KernelInfo *) NULL)
1361           return(kernel);
1362         kernel->type = type;
1363         RotateKernelInfo(kernel, args->rho);
1364         break;
1365       }
1366     case KirschKernel:
1367       {
1368         kernel=ParseKernelArray("3: 5,-3,-3  5,0,-3  5,-3,-3");
1369         if (kernel == (KernelInfo *) NULL)
1370           return(kernel);
1371         kernel->type = type;
1372         RotateKernelInfo(kernel, args->rho);
1373         break;
1374       }
1375     case FreiChenKernel:
1376       /* Direction is set to be left to right positive */
1377       /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1378       /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1379       { switch ( (int) args->rho ) {
1380           default:
1381           case 0:
1382             kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1383             if (kernel == (KernelInfo *) NULL)
1384               return(kernel);
1385             kernel->type = type;
1386             kernel->values[3] = +MagickSQ2;
1387             kernel->values[5] = -MagickSQ2;
1388             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1389             break;
1390           case 2:
1391             kernel=ParseKernelArray("3: 1,2,0  2,0,-2  0,-2,-1");
1392             if (kernel == (KernelInfo *) NULL)
1393               return(kernel);
1394             kernel->type = type;
1395             kernel->values[1] = kernel->values[3] = +MagickSQ2;
1396             kernel->values[5] = kernel->values[7] = -MagickSQ2;
1397             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1398             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1399             break;
1400           case 10:
1401             kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1402             if (kernel == (KernelInfo *) NULL)
1403               return(kernel);
1404             break;
1405           case 1:
1406           case 11:
1407             kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
1408             if (kernel == (KernelInfo *) NULL)
1409               return(kernel);
1410             kernel->type = type;
1411             kernel->values[3] = +MagickSQ2;
1412             kernel->values[5] = -MagickSQ2;
1413             CalcKernelMetaData(kernel);     /* recalculate meta-data */
1414             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1415             break;
1416           case 12:
1417             kernel=ParseKernelArray("3: 1,2,1  0,0,0  1,2,1");
1418             if (kernel == (KernelInfo *) NULL)
1419               return(kernel);
1420             kernel->type = type;
1421             kernel->values[1] = +MagickSQ2;
1422             kernel->values[7] = +MagickSQ2;
1423             CalcKernelMetaData(kernel);
1424             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1425             break;
1426           case 13:
1427             kernel=ParseKernelArray("3: 2,-1,0  -1,0,1  0,1,-2");
1428             if (kernel == (KernelInfo *) NULL)
1429               return(kernel);
1430             kernel->type = type;
1431             kernel->values[0] = +MagickSQ2;
1432             kernel->values[8] = -MagickSQ2;
1433             CalcKernelMetaData(kernel);
1434             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1435             break;
1436           case 14:
1437             kernel=ParseKernelArray("3: 0,1,-2  -1,0,1  2,-1,0");
1438             if (kernel == (KernelInfo *) NULL)
1439               return(kernel);
1440             kernel->type = type;
1441             kernel->values[2] = -MagickSQ2;
1442             kernel->values[6] = +MagickSQ2;
1443             CalcKernelMetaData(kernel);
1444             ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1445             break;
1446           case 15:
1447             kernel=ParseKernelArray("3: 0,-1,0  1,0,1  0,-1,0");
1448             if (kernel == (KernelInfo *) NULL)
1449               return(kernel);
1450             kernel->type = type;
1451             ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1452             break;
1453           case 16:
1454             kernel=ParseKernelArray("3: 1,0,-1  0,0,0  -1,0,1");
1455             if (kernel == (KernelInfo *) NULL)
1456               return(kernel);
1457             kernel->type = type;
1458             ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1459             break;
1460           case 17:
1461             kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  -1,-2,1");
1462             if (kernel == (KernelInfo *) NULL)
1463               return(kernel);
1464             kernel->type = type;
1465             ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1466             break;
1467           case 18:
1468             kernel=ParseKernelArray("3: -2,1,-2  1,4,1  -2,1,-2");
1469             if (kernel == (KernelInfo *) NULL)
1470               return(kernel);
1471             kernel->type = type;
1472             ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1473             break;
1474           case 19:
1475             kernel=ParseKernelArray("3: 1,1,1  1,1,1  1,1,1");
1476             if (kernel == (KernelInfo *) NULL)
1477               return(kernel);
1478             kernel->type = type;
1479             ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1480             break;
1481         }
1482         if ( fabs(args->sigma) > MagickEpsilon )
1483           /* Rotate by correctly supplied 'angle' */
1484           RotateKernelInfo(kernel, args->sigma);
1485         else if ( args->rho > 30.0 || args->rho < -30.0 )
1486           /* Rotate by out of bounds 'type' */
1487           RotateKernelInfo(kernel, args->rho);
1488         break;
1489       }
1490
1491     /*
1492       Boolean or Shaped Kernels
1493     */
1494     case DiamondKernel:
1495       {
1496         if (args->rho < 1.0)
1497           kernel->width = kernel->height = 3;  /* default radius = 1 */
1498         else
1499           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1500         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1501
1502         kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1503                               kernel->height*sizeof(double));
1504         if (kernel->values == (double *) NULL)
1505           return(DestroyKernelInfo(kernel));
1506
1507         /* set all kernel values within diamond area to scale given */
1508         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1509           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1510             if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1511               kernel->positive_range += kernel->values[i] = args->sigma;
1512             else
1513               kernel->values[i] = nan;
1514         kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1515         break;
1516       }
1517     case SquareKernel:
1518     case RectangleKernel:
1519       { double
1520           scale;
1521         if ( type == SquareKernel )
1522           {
1523             if (args->rho < 1.0)
1524               kernel->width = kernel->height = 3;  /* default radius = 1 */
1525             else
1526               kernel->width = kernel->height = (size_t) (2*args->rho+1);
1527             kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1528             scale = args->sigma;
1529           }
1530         else {
1531             /* NOTE: user defaults set in "AcquireKernelInfo()" */
1532             if ( args->rho < 1.0 || args->sigma < 1.0 )
1533               return(DestroyKernelInfo(kernel));    /* invalid args given */
1534             kernel->width = (size_t)args->rho;
1535             kernel->height = (size_t)args->sigma;
1536             if ( args->xi  < 0.0 || args->xi  > (double)kernel->width ||
1537                  args->psi < 0.0 || args->psi > (double)kernel->height )
1538               return(DestroyKernelInfo(kernel));    /* invalid args given */
1539             kernel->x = (ssize_t) args->xi;
1540             kernel->y = (ssize_t) args->psi;
1541             scale = 1.0;
1542           }
1543         kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1544                               kernel->height*sizeof(double));
1545         if (kernel->values == (double *) NULL)
1546           return(DestroyKernelInfo(kernel));
1547
1548         /* set all kernel values to scale given */
1549         u=(ssize_t) (kernel->width*kernel->height);
1550         for ( i=0; i < u; i++)
1551             kernel->values[i] = scale;
1552         kernel->minimum = kernel->maximum = scale;   /* a flat shape */
1553         kernel->positive_range = scale*u;
1554         break;
1555       }
1556       case OctagonKernel:
1557         {
1558           if (args->rho < 1.0)
1559             kernel->width = kernel->height = 5;  /* default radius = 2 */
1560           else
1561             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1562           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1563
1564           kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1565                                 kernel->height*sizeof(double));
1566           if (kernel->values == (double *) NULL)
1567             return(DestroyKernelInfo(kernel));
1568
1569           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1570             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1571               if ( (labs((long) u)+labs((long) v)) <=
1572                         ((long)kernel->x + (long)(kernel->x/2)) )
1573                 kernel->positive_range += kernel->values[i] = args->sigma;
1574               else
1575                 kernel->values[i] = nan;
1576           kernel->minimum = kernel->maximum = args->sigma;  /* a flat shape */
1577           break;
1578         }
1579       case DiskKernel:
1580         {
1581           ssize_t
1582             limit = (ssize_t)(args->rho*args->rho);
1583
1584           if (args->rho < 0.4)           /* default radius approx 4.3 */
1585             kernel->width = kernel->height = 9L, limit = 18L;
1586           else
1587             kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1588           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1589
1590           kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1591                                 kernel->height*sizeof(double));
1592           if (kernel->values == (double *) NULL)
1593             return(DestroyKernelInfo(kernel));
1594
1595           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1596             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1597               if ((u*u+v*v) <= limit)
1598                 kernel->positive_range += kernel->values[i] = args->sigma;
1599               else
1600                 kernel->values[i] = nan;
1601           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1602           break;
1603         }
1604       case PlusKernel:
1605         {
1606           if (args->rho < 1.0)
1607             kernel->width = kernel->height = 5;  /* default radius 2 */
1608           else
1609             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1610           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1611
1612           kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1613                                 kernel->height*sizeof(double));
1614           if (kernel->values == (double *) NULL)
1615             return(DestroyKernelInfo(kernel));
1616
1617           /* set all kernel values along axises to given scale */
1618           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1619             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1620               kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1621           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1622           kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1623           break;
1624         }
1625       case CrossKernel:
1626         {
1627           if (args->rho < 1.0)
1628             kernel->width = kernel->height = 5;  /* default radius 2 */
1629           else
1630             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1631           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1632
1633           kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1634                                 kernel->height*sizeof(double));
1635           if (kernel->values == (double *) NULL)
1636             return(DestroyKernelInfo(kernel));
1637
1638           /* set all kernel values along axises to given scale */
1639           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1640             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1641               kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1642           kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
1643           kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1644           break;
1645         }
1646       /*
1647         HitAndMiss Kernels
1648       */
1649       case RingKernel:
1650       case PeaksKernel:
1651         {
1652           ssize_t
1653             limit1,
1654             limit2,
1655             scale;
1656
1657           if (args->rho < args->sigma)
1658             {
1659               kernel->width = ((size_t)args->sigma)*2+1;
1660               limit1 = (ssize_t)(args->rho*args->rho);
1661               limit2 = (ssize_t)(args->sigma*args->sigma);
1662             }
1663           else
1664             {
1665               kernel->width = ((size_t)args->rho)*2+1;
1666               limit1 = (ssize_t)(args->sigma*args->sigma);
1667               limit2 = (ssize_t)(args->rho*args->rho);
1668             }
1669           if ( limit2 <= 0 )
1670             kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1671
1672           kernel->height = kernel->width;
1673           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1674           kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1675                                 kernel->height*sizeof(double));
1676           if (kernel->values == (double *) NULL)
1677             return(DestroyKernelInfo(kernel));
1678
1679           /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1680           scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1681           for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1682             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1683               { ssize_t radius=u*u+v*v;
1684                 if (limit1 < radius && radius <= limit2)
1685                   kernel->positive_range += kernel->values[i] = (double) scale;
1686                 else
1687                   kernel->values[i] = nan;
1688               }
1689           kernel->minimum = kernel->maximum = (double) scale;
1690           if ( type == PeaksKernel ) {
1691             /* set the central point in the middle */
1692             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1693             kernel->positive_range = 1.0;
1694             kernel->maximum = 1.0;
1695           }
1696           break;
1697         }
1698       case EdgesKernel:
1699         {
1700           kernel=AcquireKernelInfo("ThinSE:482");
1701           if (kernel == (KernelInfo *) NULL)
1702             return(kernel);
1703           kernel->type = type;
1704           ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1705           break;
1706         }
1707       case CornersKernel:
1708         {
1709           kernel=AcquireKernelInfo("ThinSE:87");
1710           if (kernel == (KernelInfo *) NULL)
1711             return(kernel);
1712           kernel->type = type;
1713           ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1714           break;
1715         }
1716       case DiagonalsKernel:
1717         {
1718           switch ( (int) args->rho ) {
1719             case 0:
1720             default:
1721               { KernelInfo
1722                   *new_kernel;
1723                 kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
1724                 if (kernel == (KernelInfo *) NULL)
1725                   return(kernel);
1726                 kernel->type = type;
1727                 new_kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
1728                 if (new_kernel == (KernelInfo *) NULL)
1729                   return(DestroyKernelInfo(kernel));
1730                 new_kernel->type = type;
1731                 LastKernelInfo(kernel)->next = new_kernel;
1732                 ExpandMirrorKernelInfo(kernel);
1733                 return(kernel);
1734               }
1735             case 1:
1736               kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
1737               break;
1738             case 2:
1739               kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
1740               break;
1741           }
1742           if (kernel == (KernelInfo *) NULL)
1743             return(kernel);
1744           kernel->type = type;
1745           RotateKernelInfo(kernel, args->sigma);
1746           break;
1747         }
1748       case LineEndsKernel:
1749         { /* Kernels for finding the end of thin lines */
1750           switch ( (int) args->rho ) {
1751             case 0:
1752             default:
1753               /* set of kernels to find all end of lines */
1754               return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1755             case 1:
1756               /* kernel for 4-connected line ends - no rotation */
1757               kernel=ParseKernelArray("3: 0,0,-  0,1,1  0,0,-");
1758               break;
1759           case 2:
1760               /* kernel to add for 8-connected lines - no rotation */
1761               kernel=ParseKernelArray("3: 0,0,0  0,1,0  0,0,1");
1762               break;
1763           case 3:
1764               /* kernel to add for orthogonal line ends - does not find corners */
1765               kernel=ParseKernelArray("3: 0,0,0  0,1,1  0,0,0");
1766               break;
1767           case 4:
1768               /* traditional line end - fails on last T end */
1769               kernel=ParseKernelArray("3: 0,0,0  0,1,-  0,0,-");
1770               break;
1771           }
1772           if (kernel == (KernelInfo *) NULL)
1773             return(kernel);
1774           kernel->type = type;
1775           RotateKernelInfo(kernel, args->sigma);
1776           break;
1777         }
1778       case LineJunctionsKernel:
1779         { /* kernels for finding the junctions of multiple lines */
1780           switch ( (int) args->rho ) {
1781             case 0:
1782             default:
1783               /* set of kernels to find all line junctions */
1784               return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1785             case 1:
1786               /* Y Junction */
1787               kernel=ParseKernelArray("3: 1,-,1  -,1,-  -,1,-");
1788               break;
1789             case 2:
1790               /* Diagonal T Junctions */
1791               kernel=ParseKernelArray("3: 1,-,-  -,1,-  1,-,1");
1792               break;
1793             case 3:
1794               /* Orthogonal T Junctions */
1795               kernel=ParseKernelArray("3: -,-,-  1,1,1  -,1,-");
1796               break;
1797             case 4:
1798               /* Diagonal X Junctions */
1799               kernel=ParseKernelArray("3: 1,-,1  -,1,-  1,-,1");
1800               break;
1801             case 5:
1802               /* Orthogonal X Junctions - minimal diamond kernel */
1803               kernel=ParseKernelArray("3: -,1,-  1,1,1  -,1,-");
1804               break;
1805           }
1806           if (kernel == (KernelInfo *) NULL)
1807             return(kernel);
1808           kernel->type = type;
1809           RotateKernelInfo(kernel, args->sigma);
1810           break;
1811         }
1812       case RidgesKernel:
1813         { /* Ridges - Ridge finding kernels */
1814           KernelInfo
1815             *new_kernel;
1816           switch ( (int) args->rho ) {
1817             case 1:
1818             default:
1819               kernel=ParseKernelArray("3x1:0,1,0");
1820               if (kernel == (KernelInfo *) NULL)
1821                 return(kernel);
1822               kernel->type = type;
1823               ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1824               break;
1825             case 2:
1826               kernel=ParseKernelArray("4x1:0,1,1,0");
1827               if (kernel == (KernelInfo *) NULL)
1828                 return(kernel);
1829               kernel->type = type;
1830               ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1831
1832               /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1833               /* Unfortunatally we can not yet rotate a non-square kernel */
1834               /* But then we can't flip a non-symetrical kernel either */
1835               new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1836               if (new_kernel == (KernelInfo *) NULL)
1837                 return(DestroyKernelInfo(kernel));
1838               new_kernel->type = type;
1839               LastKernelInfo(kernel)->next = new_kernel;
1840               new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1841               if (new_kernel == (KernelInfo *) NULL)
1842                 return(DestroyKernelInfo(kernel));
1843               new_kernel->type = type;
1844               LastKernelInfo(kernel)->next = new_kernel;
1845               new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1846               if (new_kernel == (KernelInfo *) NULL)
1847                 return(DestroyKernelInfo(kernel));
1848               new_kernel->type = type;
1849               LastKernelInfo(kernel)->next = new_kernel;
1850               new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1851               if (new_kernel == (KernelInfo *) NULL)
1852                 return(DestroyKernelInfo(kernel));
1853               new_kernel->type = type;
1854               LastKernelInfo(kernel)->next = new_kernel;
1855               new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1856               if (new_kernel == (KernelInfo *) NULL)
1857                 return(DestroyKernelInfo(kernel));
1858               new_kernel->type = type;
1859               LastKernelInfo(kernel)->next = new_kernel;
1860               new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1861               if (new_kernel == (KernelInfo *) NULL)
1862                 return(DestroyKernelInfo(kernel));
1863               new_kernel->type = type;
1864               LastKernelInfo(kernel)->next = new_kernel;
1865               new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1866               if (new_kernel == (KernelInfo *) NULL)
1867                 return(DestroyKernelInfo(kernel));
1868               new_kernel->type = type;
1869               LastKernelInfo(kernel)->next = new_kernel;
1870               new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1871               if (new_kernel == (KernelInfo *) NULL)
1872                 return(DestroyKernelInfo(kernel));
1873               new_kernel->type = type;
1874               LastKernelInfo(kernel)->next = new_kernel;
1875               break;
1876           }
1877           break;
1878         }
1879       case ConvexHullKernel:
1880         {
1881           KernelInfo
1882             *new_kernel;
1883           /* first set of 8 kernels */
1884           kernel=ParseKernelArray("3: 1,1,-  1,0,-  1,-,0");
1885           if (kernel == (KernelInfo *) NULL)
1886             return(kernel);
1887           kernel->type = type;
1888           ExpandRotateKernelInfo(kernel, 90.0);
1889           /* append the mirror versions too - no flip function yet */
1890           new_kernel=ParseKernelArray("3: 1,1,1  1,0,-  -,-,0");
1891           if (new_kernel == (KernelInfo *) NULL)
1892             return(DestroyKernelInfo(kernel));
1893           new_kernel->type = type;
1894           ExpandRotateKernelInfo(new_kernel, 90.0);
1895           LastKernelInfo(kernel)->next = new_kernel;
1896           break;
1897         }
1898       case SkeletonKernel:
1899         {
1900           switch ( (int) args->rho ) {
1901             case 1:
1902             default:
1903               /* Traditional Skeleton...
1904               ** A cyclically rotated single kernel
1905               */
1906               kernel=AcquireKernelInfo("ThinSE:482");
1907               if (kernel == (KernelInfo *) NULL)
1908                 return(kernel);
1909               kernel->type = type;
1910               ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1911               break;
1912             case 2:
1913               /* HIPR Variation of the cyclic skeleton
1914               ** Corners of the traditional method made more forgiving,
1915               ** but the retain the same cyclic order.
1916               */
1917               kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1918               if (kernel == (KernelInfo *) NULL)
1919                 return(kernel);
1920               if (kernel->next == (KernelInfo *) NULL)
1921                 return(DestroyKernelInfo(kernel));
1922               kernel->type = type;
1923               kernel->next->type = type;
1924               ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1925               break;
1926             case 3:
1927               /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1928               ** "Connectivity-Preserving Morphological Image Thransformations"
1929               ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1930               **   http://www.leptonica.com/papers/conn.pdf
1931               */
1932               kernel=AcquireKernelInfo(
1933                             "ThinSE:41; ThinSE:42; ThinSE:43");
1934               if (kernel == (KernelInfo *) NULL)
1935                 return(kernel);
1936               kernel->type = type;
1937               kernel->next->type = type;
1938               kernel->next->next->type = type;
1939               ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1940               break;
1941            }
1942           break;
1943         }
1944       case ThinSEKernel:
1945         { /* Special kernels for general thinning, while preserving connections
1946           ** "Connectivity-Preserving Morphological Image Thransformations"
1947           ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1948           **   http://www.leptonica.com/papers/conn.pdf
1949           ** And
1950           **   http://tpgit.github.com/Leptonica/ccthin_8c_source.html
1951           **
1952           ** Note kernels do not specify the origin pixel, allowing them
1953           ** to be used for both thickening and thinning operations.
1954           */
1955           switch ( (int) args->rho ) {
1956             /* SE for 4-connected thinning */
1957             case 41: /* SE_4_1 */
1958               kernel=ParseKernelArray("3: -,-,1  0,-,1  -,-,1");
1959               break;
1960             case 42: /* SE_4_2 */
1961               kernel=ParseKernelArray("3: -,-,1  0,-,1  -,0,-");
1962               break;
1963             case 43: /* SE_4_3 */
1964               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,-,1");
1965               break;
1966             case 44: /* SE_4_4 */
1967               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,-");
1968               break;
1969             case 45: /* SE_4_5 */
1970               kernel=ParseKernelArray("3: -,0,1  0,-,1  -,0,-");
1971               break;
1972             case 46: /* SE_4_6 */
1973               kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,1");
1974               break;
1975             case 47: /* SE_4_7 */
1976               kernel=ParseKernelArray("3: -,1,1  0,-,1  -,0,-");
1977               break;
1978             case 48: /* SE_4_8 */
1979               kernel=ParseKernelArray("3: -,-,1  0,-,1  0,-,1");
1980               break;
1981             case 49: /* SE_4_9 */
1982               kernel=ParseKernelArray("3: 0,-,1  0,-,1  -,-,1");
1983               break;
1984             /* SE for 8-connected thinning - negatives of the above */
1985             case 81: /* SE_8_0 */
1986               kernel=ParseKernelArray("3: -,1,-  0,-,1  -,1,-");
1987               break;
1988             case 82: /* SE_8_2 */
1989               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,-,-");
1990               break;
1991             case 83: /* SE_8_3 */
1992               kernel=ParseKernelArray("3: 0,-,-  0,-,1  -,1,-");
1993               break;
1994             case 84: /* SE_8_4 */
1995               kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,-");
1996               break;
1997             case 85: /* SE_8_5 */
1998               kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,-");
1999               break;
2000             case 86: /* SE_8_6 */
2001               kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,1");
2002               break;
2003             case 87: /* SE_8_7 */
2004               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,0,-");
2005               break;
2006             case 88: /* SE_8_8 */
2007               kernel=ParseKernelArray("3: -,1,-  0,-,1  0,1,-");
2008               break;
2009             case 89: /* SE_8_9 */
2010               kernel=ParseKernelArray("3: 0,1,-  0,-,1  -,1,-");
2011               break;
2012             /* Special combined SE kernels */
2013             case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2014               kernel=ParseKernelArray("3: -,-,1  0,-,-  -,0,-");
2015               break;
2016             case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2017               kernel=ParseKernelArray("3: -,1,-  -,-,1  0,-,-");
2018               break;
2019             case 481: /* SE_48_1 - General Connected Corner Kernel */
2020               kernel=ParseKernelArray("3: -,1,1  0,-,1  0,0,-");
2021               break;
2022             default:
2023             case 482: /* SE_48_2 - General Edge Kernel */
2024               kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,1");
2025               break;
2026           }
2027           if (kernel == (KernelInfo *) NULL)
2028             return(kernel);
2029           kernel->type = type;
2030           RotateKernelInfo(kernel, args->sigma);
2031           break;
2032         }
2033       /*
2034         Distance Measuring Kernels
2035       */
2036       case ChebyshevKernel:
2037         {
2038           if (args->rho < 1.0)
2039             kernel->width = kernel->height = 3;  /* default radius = 1 */
2040           else
2041             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2042           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2043
2044           kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2045                                 kernel->height*sizeof(double));
2046           if (kernel->values == (double *) NULL)
2047             return(DestroyKernelInfo(kernel));
2048
2049           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2050             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2051               kernel->positive_range += ( kernel->values[i] =
2052                   args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2053           kernel->maximum = kernel->values[0];
2054           break;
2055         }
2056       case ManhattanKernel:
2057         {
2058           if (args->rho < 1.0)
2059             kernel->width = kernel->height = 3;  /* default radius = 1 */
2060           else
2061             kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2062           kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2063
2064           kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2065                                 kernel->height*sizeof(double));
2066           if (kernel->values == (double *) NULL)
2067             return(DestroyKernelInfo(kernel));
2068
2069           for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2070             for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2071               kernel->positive_range += ( kernel->values[i] =
2072                   args->sigma*(labs((long) u)+labs((long) v)) );
2073           kernel->maximum = kernel->values[0];
2074           break;
2075         }
2076       case OctagonalKernel:
2077       {
2078         if (args->rho < 2.0)
2079           kernel->width = kernel->height = 5;  /* default/minimum radius = 2 */
2080         else
2081           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2082         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2083
2084         kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2085                               kernel->height*sizeof(double));
2086         if (kernel->values == (double *) NULL)
2087           return(DestroyKernelInfo(kernel));
2088
2089         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2090           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2091             {
2092               double
2093                 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2094                 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2095               kernel->positive_range += kernel->values[i] =
2096                         args->sigma*MagickMax(r1,r2);
2097             }
2098         kernel->maximum = kernel->values[0];
2099         break;
2100       }
2101     case EuclideanKernel:
2102       {
2103         if (args->rho < 1.0)
2104           kernel->width = kernel->height = 3;  /* default radius = 1 */
2105         else
2106           kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2107         kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2108
2109         kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2110                               kernel->height*sizeof(double));
2111         if (kernel->values == (double *) NULL)
2112           return(DestroyKernelInfo(kernel));
2113
2114         for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2115           for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2116             kernel->positive_range += ( kernel->values[i] =
2117                  args->sigma*sqrt((double)(u*u+v*v)) );
2118         kernel->maximum = kernel->values[0];
2119         break;
2120       }
2121     default:
2122       {
2123         /* No-Op Kernel - Basically just a single pixel on its own */
2124         kernel=ParseKernelArray("1:1");
2125         if (kernel == (KernelInfo *) NULL)
2126           return(kernel);
2127         kernel->type = UndefinedKernel;
2128         break;
2129       }
2130       break;
2131   }
2132   return(kernel);
2133 }
2134 \f
2135 /*
2136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2137 %                                                                             %
2138 %                                                                             %
2139 %                                                                             %
2140 %     C l o n e K e r n e l I n f o                                           %
2141 %                                                                             %
2142 %                                                                             %
2143 %                                                                             %
2144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2145 %
2146 %  CloneKernelInfo() creates a new clone of the given Kernel List so that its
2147 %  can be modified without effecting the original.  The cloned kernel should
2148 %  be destroyed using DestoryKernelInfo() when no longer needed.
2149 %
2150 %  The format of the CloneKernelInfo method is:
2151 %
2152 %      KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2153 %
2154 %  A description of each parameter follows:
2155 %
2156 %    o kernel: the Morphology/Convolution kernel to be cloned
2157 %
2158 */
2159 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2160 {
2161   register ssize_t
2162     i;
2163
2164   KernelInfo
2165     *new_kernel;
2166
2167   assert(kernel != (KernelInfo *) NULL);
2168   new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2169   if (new_kernel == (KernelInfo *) NULL)
2170     return(new_kernel);
2171   *new_kernel=(*kernel); /* copy values in structure */
2172
2173   /* replace the values with a copy of the values */
2174   new_kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2175     kernel->height*sizeof(double));
2176   if (new_kernel->values == (double *) NULL)
2177     return(DestroyKernelInfo(new_kernel));
2178   for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2179     new_kernel->values[i]=kernel->values[i];
2180
2181   /* Also clone the next kernel in the kernel list */
2182   if ( kernel->next != (KernelInfo *) NULL ) {
2183     new_kernel->next = CloneKernelInfo(kernel->next);
2184     if ( new_kernel->next == (KernelInfo *) NULL )
2185       return(DestroyKernelInfo(new_kernel));
2186   }
2187
2188   return(new_kernel);
2189 }
2190 \f
2191 /*
2192 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2193 %                                                                             %
2194 %                                                                             %
2195 %                                                                             %
2196 %     D e s t r o y K e r n e l I n f o                                       %
2197 %                                                                             %
2198 %                                                                             %
2199 %                                                                             %
2200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2201 %
2202 %  DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2203 %  kernel.
2204 %
2205 %  The format of the DestroyKernelInfo method is:
2206 %
2207 %      KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2208 %
2209 %  A description of each parameter follows:
2210 %
2211 %    o kernel: the Morphology/Convolution kernel to be destroyed
2212 %
2213 */
2214 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2215 {
2216   assert(kernel != (KernelInfo *) NULL);
2217
2218   if ( kernel->next != (KernelInfo *) NULL )
2219     kernel->next = DestroyKernelInfo(kernel->next);
2220
2221   kernel->values = (double *)RelinquishMagickMemory(kernel->values);
2222   kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
2223   return(kernel);
2224 }
2225 \f
2226 /*
2227 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2228 %                                                                             %
2229 %                                                                             %
2230 %                                                                             %
2231 +     E x p a n d M i r r o r K e r n e l I n f o                             %
2232 %                                                                             %
2233 %                                                                             %
2234 %                                                                             %
2235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2236 %
2237 %  ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2238 %  sequence of 90-degree rotated kernels but providing a reflected 180
2239 %  rotatation, before the -/+ 90-degree rotations.
2240 %
2241 %  This special rotation order produces a better, more symetrical thinning of
2242 %  objects.
2243 %
2244 %  The format of the ExpandMirrorKernelInfo method is:
2245 %
2246 %      void ExpandMirrorKernelInfo(KernelInfo *kernel)
2247 %
2248 %  A description of each parameter follows:
2249 %
2250 %    o kernel: the Morphology/Convolution kernel
2251 %
2252 % This function is only internel to this module, as it is not finalized,
2253 % especially with regard to non-orthogonal angles, and rotation of larger
2254 % 2D kernels.
2255 */
2256
2257 #if 0
2258 static void FlopKernelInfo(KernelInfo *kernel)
2259     { /* Do a Flop by reversing each row. */
2260       size_t
2261         y;
2262       register ssize_t
2263         x,r;
2264       register double
2265         *k,t;
2266
2267       for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2268         for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2269           t=k[x],  k[x]=k[r],  k[r]=t;
2270
2271       kernel->x = kernel->width - kernel->x - 1;
2272       angle = fmod(angle+180.0, 360.0);
2273     }
2274 #endif
2275
2276 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2277 {
2278   KernelInfo
2279     *clone,
2280     *last;
2281
2282   last = kernel;
2283
2284   clone = CloneKernelInfo(last);
2285   RotateKernelInfo(clone, 180);   /* flip */
2286   LastKernelInfo(last)->next = clone;
2287   last = clone;
2288
2289   clone = CloneKernelInfo(last);
2290   RotateKernelInfo(clone, 90);   /* transpose */
2291   LastKernelInfo(last)->next = clone;
2292   last = clone;
2293
2294   clone = CloneKernelInfo(last);
2295   RotateKernelInfo(clone, 180);  /* flop */
2296   LastKernelInfo(last)->next = clone;
2297
2298   return;
2299 }
2300 \f
2301 /*
2302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2303 %                                                                             %
2304 %                                                                             %
2305 %                                                                             %
2306 +     E x p a n d R o t a t e K e r n e l I n f o                             %
2307 %                                                                             %
2308 %                                                                             %
2309 %                                                                             %
2310 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2311 %
2312 %  ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2313 %  incrementally by the angle given, until the kernel repeats.
2314 %
2315 %  WARNING: 45 degree rotations only works for 3x3 kernels.
2316 %  While 90 degree roatations only works for linear and square kernels
2317 %
2318 %  The format of the ExpandRotateKernelInfo method is:
2319 %
2320 %      void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2321 %
2322 %  A description of each parameter follows:
2323 %
2324 %    o kernel: the Morphology/Convolution kernel
2325 %
2326 %    o angle: angle to rotate in degrees
2327 %
2328 % This function is only internel to this module, as it is not finalized,
2329 % especially with regard to non-orthogonal angles, and rotation of larger
2330 % 2D kernels.
2331 */
2332
2333 /* Internal Routine - Return true if two kernels are the same */
2334 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2335      const KernelInfo *kernel2)
2336 {
2337   register size_t
2338     i;
2339
2340   /* check size and origin location */
2341   if (    kernel1->width != kernel2->width
2342        || kernel1->height != kernel2->height
2343        || kernel1->x != kernel2->x
2344        || kernel1->y != kernel2->y )
2345     return MagickFalse;
2346
2347   /* check actual kernel values */
2348   for (i=0; i < (kernel1->width*kernel1->height); i++) {
2349     /* Test for Nan equivalence */
2350     if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2351       return MagickFalse;
2352     if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2353       return MagickFalse;
2354     /* Test actual values are equivalent */
2355     if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
2356       return MagickFalse;
2357   }
2358
2359   return MagickTrue;
2360 }
2361
2362 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2363 {
2364   KernelInfo
2365     *clone,
2366     *last;
2367
2368   last = kernel;
2369   while(1) {
2370     clone = CloneKernelInfo(last);
2371     RotateKernelInfo(clone, angle);
2372     if ( SameKernelInfo(kernel, clone) == MagickTrue )
2373       break;
2374     LastKernelInfo(last)->next = clone;
2375     last = clone;
2376   }
2377   clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2378   return;
2379 }
2380 \f
2381 /*
2382 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2383 %                                                                             %
2384 %                                                                             %
2385 %                                                                             %
2386 +     C a l c M e t a K e r n a l I n f o                                     %
2387 %                                                                             %
2388 %                                                                             %
2389 %                                                                             %
2390 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2391 %
2392 %  CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2393 %  using the kernel values.  This should only ne used if it is not possible to
2394 %  calculate that meta-data in some easier way.
2395 %
2396 %  It is important that the meta-data is correct before ScaleKernelInfo() is
2397 %  used to perform kernel normalization.
2398 %
2399 %  The format of the CalcKernelMetaData method is:
2400 %
2401 %      void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2402 %
2403 %  A description of each parameter follows:
2404 %
2405 %    o kernel: the Morphology/Convolution kernel to modify
2406 %
2407 %  WARNING: Minimum and Maximum values are assumed to include zero, even if
2408 %  zero is not part of the kernel (as in Gaussian Derived kernels). This
2409 %  however is not true for flat-shaped morphological kernels.
2410 %
2411 %  WARNING: Only the specific kernel pointed to is modified, not a list of
2412 %  multiple kernels.
2413 %
2414 % This is an internal function and not expected to be useful outside this
2415 % module.  This could change however.
2416 */
2417 static void CalcKernelMetaData(KernelInfo *kernel)
2418 {
2419   register size_t
2420     i;
2421
2422   kernel->minimum = kernel->maximum = 0.0;
2423   kernel->negative_range = kernel->positive_range = 0.0;
2424   for (i=0; i < (kernel->width*kernel->height); i++)
2425     {
2426       if ( fabs(kernel->values[i]) < MagickEpsilon )
2427         kernel->values[i] = 0.0;
2428       ( kernel->values[i] < 0)
2429           ?  ( kernel->negative_range += kernel->values[i] )
2430           :  ( kernel->positive_range += kernel->values[i] );
2431       Minimize(kernel->minimum, kernel->values[i]);
2432       Maximize(kernel->maximum, kernel->values[i]);
2433     }
2434
2435   return;
2436 }
2437 \f
2438 /*
2439 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2440 %                                                                             %
2441 %                                                                             %
2442 %                                                                             %
2443 %     M o r p h o l o g y A p p l y                                           %
2444 %                                                                             %
2445 %                                                                             %
2446 %                                                                             %
2447 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2448 %
2449 %  MorphologyApply() applies a morphological method, multiple times using
2450 %  a list of multiple kernels.
2451 %
2452 %  It is basically equivalent to as MorphologyImage() (see below) but
2453 %  without any user controls.  This allows internel programs to use this
2454 %  function, to actually perform a specific task without possible interference
2455 %  by any API user supplied settings.
2456 %
2457 %  It is MorphologyImage() task to extract any such user controls, and
2458 %  pass them to this function for processing.
2459 %
2460 %  More specifically kernels are not normalized/scaled/blended by the
2461 %  'convolve:scale' Image Artifact (setting), nor is the convolve bias
2462 %  (-bias setting or image->bias) loooked at, but must be supplied from the
2463 %  function arguments.
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 DistanceMorphology:*/
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(dynamic,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->sync == MagickFalse) || (image->matte == MagickFalse) )
2663           { /* No 'Sync' involved.
2664             ** Convolution is simple greyscale channel operation
2665             */
2666             for (v=0; v < (ssize_t) kernel->height; v++) {
2667               if ( IsNan(*k) ) continue;
2668               result.red     += (*k)*GetPixelRed(image,k_pixels);
2669               result.green   += (*k)*GetPixelGreen(image,k_pixels);
2670               result.blue    += (*k)*GetPixelBlue(image,k_pixels);
2671               if (image->colorspace == CMYKColorspace)
2672                 result.black+=(*k)*GetPixelBlack(image,k_pixels);
2673               result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2674               k--;
2675               k_pixels+=GetPixelChannels(image);
2676             }
2677             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2678               SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
2679             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2680               SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
2681             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2682               SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
2683             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2684                 (image->colorspace == CMYKColorspace))
2685               SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
2686             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2687                 (image->matte == MagickTrue))
2688               SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2689           }
2690         else
2691           { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2692             ** Weight the color channels with Alpha Channel so that
2693             ** transparent pixels are not part of the results.
2694             */
2695             MagickRealType
2696               alpha,  /* alpha weighting of colors : kernel*alpha  */
2697               gamma;  /* divisor, sum of color weighting values */
2698
2699             gamma=0.0;
2700             for (v=0; v < (ssize_t) kernel->height; v++) {
2701               if ( IsNan(*k) ) continue;
2702               alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels));
2703               gamma += alpha;
2704               result.red     += alpha*GetPixelRed(image,k_pixels);
2705               result.green   += alpha*GetPixelGreen(image,k_pixels);
2706               result.blue    += alpha*GetPixelBlue(image,k_pixels);
2707               if (image->colorspace == CMYKColorspace)
2708                 result.black += alpha*GetPixelBlack(image,k_pixels);
2709               result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2710               k--;
2711               k_pixels+=GetPixelChannels(image);
2712             }
2713             /* Sync'ed channels, all channels are modified */
2714             gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2715             SetPixelRed(morphology_image,
2716               ClampToQuantum(gamma*result.red),q);
2717             SetPixelGreen(morphology_image,
2718               ClampToQuantum(gamma*result.green),q);
2719             SetPixelBlue(morphology_image,
2720               ClampToQuantum(gamma*result.blue),q);
2721             if (image->colorspace == CMYKColorspace)
2722               SetPixelBlack(morphology_image,
2723                 ClampToQuantum(gamma*result.black),q);
2724             SetPixelAlpha(morphology_image,
2725               ClampToQuantum(result.alpha),q);
2726           }
2727
2728         /* Count up changed pixels */
2729         if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q))
2730             || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q))
2731             || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q))
2732             || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q))
2733             || ((image->colorspace == CMYKColorspace) &&
2734                 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
2735           changed++;  /* The pixel was changed in some way! */
2736         p+=GetPixelChannels(image);
2737         q+=GetPixelChannels(morphology_image);
2738       } /* y */
2739       if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2740         status=MagickFalse;
2741       if (image->progress_monitor != (MagickProgressMonitor) NULL)
2742         {
2743           MagickBooleanType
2744             proceed;
2745
2746 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2747   #pragma omp critical (MagickCore_MorphologyImage)
2748 #endif
2749           proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2750           if (proceed == MagickFalse)
2751             status=MagickFalse;
2752         }
2753     } /* x */
2754     morphology_image->type=image->type;
2755     morphology_view=DestroyCacheView(morphology_view);
2756     image_view=DestroyCacheView(image_view);
2757     return(status ? (ssize_t) changed : 0);
2758   }
2759
2760   /*
2761   ** Normal handling of horizontal or rectangular kernels (row by row)
2762   */
2763 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2764   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2765 #endif
2766   for (y=0; y < (ssize_t) image->rows; y++)
2767   {
2768     register const Quantum
2769       *restrict p;
2770
2771     register Quantum
2772       *restrict q;
2773
2774     register ssize_t
2775       x;
2776
2777     size_t
2778       r;
2779
2780     if (status == MagickFalse)
2781       continue;
2782     p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
2783       kernel->height,  exception);
2784     q=GetCacheViewAuthenticPixels(morphology_view,0,y,
2785       morphology_image->columns,1,exception);
2786     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2787       {
2788         status=MagickFalse;
2789         continue;
2790       }
2791     /* offset to origin in 'p'. while 'q' points to it directly */
2792     r = virt_width*offy + offx;
2793
2794     for (x=0; x < (ssize_t) image->columns; x++)
2795     {
2796        ssize_t
2797         v;
2798
2799       register ssize_t
2800         u;
2801
2802       register const double
2803         *restrict k;
2804
2805       register const Quantum
2806         *restrict k_pixels;
2807
2808       PixelInfo
2809         result,
2810         min,
2811         max;
2812
2813       /* Copy input image to the output image for unused channels
2814        * This removes need for 'cloning' a new image every iteration
2815        */
2816       SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2817         GetPixelChannels(image)),q);
2818       SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2819         GetPixelChannels(image)),q);
2820       SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2821         GetPixelChannels(image)),q);
2822       if (image->colorspace == CMYKColorspace)
2823         SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2824           GetPixelChannels(image)),q);
2825
2826       /* Defaults */
2827       min.red     =
2828       min.green   =
2829       min.blue    =
2830       min.alpha =
2831       min.black   = (MagickRealType) QuantumRange;
2832       max.red     =
2833       max.green   =
2834       max.blue    =
2835       max.alpha =
2836       max.black   = (MagickRealType) 0;
2837       /* default result is the original pixel value */
2838       result.red     = (MagickRealType) GetPixelRed(image,p+r*GetPixelChannels(image));
2839       result.green   = (MagickRealType) GetPixelGreen(image,p+r*GetPixelChannels(image));
2840       result.blue    = (MagickRealType) GetPixelBlue(image,p+r*GetPixelChannels(image));
2841       result.black   = 0.0;
2842       if (image->colorspace == CMYKColorspace)
2843         result.black = (MagickRealType) GetPixelBlack(image,p+r*GetPixelChannels(image));
2844       result.alpha=(MagickRealType) GetPixelAlpha(image,p+r*GetPixelChannels(image));
2845
2846       switch (method) {
2847         case ConvolveMorphology:
2848           /* Set the bias of the weighted average output */
2849           result.red     =
2850           result.green   =
2851           result.blue    =
2852           result.alpha =
2853           result.black   = bias;
2854           break;
2855         case DilateIntensityMorphology:
2856         case ErodeIntensityMorphology:
2857           /* use a boolean flag indicating when first match found */
2858           result.red = 0.0;  /* result is not used otherwise */
2859           break;
2860         default:
2861           break;
2862       }
2863
2864       switch ( method ) {
2865         case ConvolveMorphology:
2866             /* Weighted Average of pixels using reflected kernel
2867             **
2868             ** NOTE for correct working of this operation for asymetrical
2869             ** kernels, the kernel needs to be applied in its reflected form.
2870             ** That is its values needs to be reversed.
2871             **
2872             ** Correlation is actually the same as this but without reflecting
2873             ** the kernel, and thus 'lower-level' that Convolution.  However
2874             ** as Convolution is the more common method used, and it does not
2875             ** really cost us much in terms of processing to use a reflected
2876             ** kernel, so it is Convolution that is implemented.
2877             **
2878             ** Correlation will have its kernel reflected before calling
2879             ** this function to do a Convolve.
2880             **
2881             ** For more details of Correlation vs Convolution see
2882             **   http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2883             */
2884             k = &kernel->values[ kernel->width*kernel->height-1 ];
2885             k_pixels = p;
2886             if ( (image->sync == MagickFalse) ||
2887                  (image->matte == MagickFalse) )
2888               { /* No 'Sync' involved.
2889                 ** Convolution is simple greyscale channel operation
2890                 */
2891                 for (v=0; v < (ssize_t) kernel->height; v++) {
2892                   for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2893                     if ( IsNan(*k) ) continue;
2894                     result.red     += (*k)*
2895                       GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2896                     result.green   += (*k)*
2897                       GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2898                     result.blue    += (*k)*
2899                       GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2900                     if (image->colorspace == CMYKColorspace)
2901                       result.black += (*k)*
2902                         GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2903                     result.alpha += (*k)*
2904                       GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2905                   }
2906                   k_pixels += virt_width*GetPixelChannels(image);
2907                 }
2908                 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2909                   SetPixelRed(morphology_image,ClampToQuantum(result.red),
2910                     q);
2911                 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2912                   SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2913                     q);
2914                 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2915                   SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2916                     q);
2917                 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2918                     (image->colorspace == CMYKColorspace))
2919                   SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2920                     q);
2921                 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2922                     (image->matte == MagickTrue))
2923                   SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
2924                     q);
2925               }
2926             else
2927               { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2928                 ** Weight the color channels with Alpha Channel so that
2929                 ** transparent pixels are not part of the results.
2930                 */
2931                 MagickRealType
2932                   alpha,  /* alpha weighting of colors : kernel*alpha  */
2933                   gamma;  /* divisor, sum of color weighting values */
2934
2935                 gamma=0.0;
2936                 for (v=0; v < (ssize_t) kernel->height; v++) {
2937                   for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2938                     if ( IsNan(*k) ) continue;
2939                     alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels+u*
2940                       GetPixelChannels(image)));
2941                     gamma += alpha;
2942                     result.red     += alpha*
2943                       GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2944                     result.green   += alpha*
2945                       GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2946                     result.blue    += alpha*
2947                       GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2948                     if (image->colorspace == CMYKColorspace)
2949                       result.black+=alpha*
2950                         GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2951                     result.alpha += (*k)*
2952                       GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2953                   }
2954                   k_pixels += virt_width*GetPixelChannels(image);
2955                 }
2956                 /* Sync'ed channels, all channels are modified */
2957                 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2958                 SetPixelRed(morphology_image,
2959                   ClampToQuantum(gamma*result.red),q);
2960                 SetPixelGreen(morphology_image,
2961                   ClampToQuantum(gamma*result.green),q);
2962                 SetPixelBlue(morphology_image,
2963                   ClampToQuantum(gamma*result.blue),q);
2964                 if (image->colorspace == CMYKColorspace)
2965                   SetPixelBlack(morphology_image,
2966                     ClampToQuantum(gamma*result.black),q);
2967                 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2968               }
2969             break;
2970
2971         case ErodeMorphology:
2972             /* Minimum Value within kernel neighbourhood
2973             **
2974             ** NOTE that the kernel is not reflected for this operation!
2975             **
2976             ** NOTE: in normal Greyscale Morphology, the kernel value should
2977             ** be added to the real value, this is currently not done, due to
2978             ** the nature of the boolean kernels being used.
2979             */
2980             k = kernel->values;
2981             k_pixels = p;
2982             for (v=0; v < (ssize_t) kernel->height; v++) {
2983               for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2984                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2985                 Minimize(min.red,     (double)
2986                   GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
2987                 Minimize(min.green,   (double)
2988                   GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
2989                 Minimize(min.blue,    (double)
2990                   GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
2991                 if (image->colorspace == CMYKColorspace)
2992                   Minimize(min.black,(double)
2993                     GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
2994                 Minimize(min.alpha,(double)
2995                   GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
2996               }
2997               k_pixels += virt_width*GetPixelChannels(image);
2998             }
2999             break;
3000
3001         case DilateMorphology:
3002             /* Maximum Value within kernel neighbourhood
3003             **
3004             ** NOTE for correct working of this operation for asymetrical
3005             ** kernels, the kernel needs to be applied in its reflected form.
3006             ** That is its values needs to be reversed.
3007             **
3008             ** NOTE: in normal Greyscale Morphology, the kernel value should
3009             ** be added to the real value, this is currently not done, due to
3010             ** the nature of the boolean kernels being used.
3011             **
3012             */
3013             k = &kernel->values[ kernel->width*kernel->height-1 ];
3014             k_pixels = p;
3015             for (v=0; v < (ssize_t) kernel->height; v++) {
3016               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3017                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3018                 Maximize(max.red,     (double)
3019                   GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3020                 Maximize(max.green,   (double)
3021                   GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3022                 Maximize(max.blue,    (double)
3023                   GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3024                 if (image->colorspace == CMYKColorspace)
3025                   Maximize(max.black,   (double)
3026                     GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3027                 Maximize(max.alpha,(double)
3028                   GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3029               }
3030               k_pixels += virt_width*GetPixelChannels(image);
3031             }
3032             break;
3033
3034         case HitAndMissMorphology:
3035         case ThinningMorphology:
3036         case ThickenMorphology:
3037             /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3038             **
3039             ** NOTE that the kernel is not reflected for this operation,
3040             ** and consists of both foreground and background pixel
3041             ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3042             ** with either Nan or 0.5 values for don't care.
3043             **
3044             ** Note that this will never produce a meaningless negative
3045             ** result.  Such results can cause Thinning/Thicken to not work
3046             ** correctly when used against a greyscale image.
3047             */
3048             k = kernel->values;
3049             k_pixels = p;
3050             for (v=0; v < (ssize_t) kernel->height; v++) {
3051               for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3052                 if ( IsNan(*k) ) continue;
3053                 if ( (*k) > 0.7 )
3054                 { /* minimim of foreground pixels */
3055                   Minimize(min.red,     (double)
3056                     GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3057                   Minimize(min.green,   (double)
3058                     GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3059                   Minimize(min.blue,    (double)
3060                     GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3061                   if ( image->colorspace == CMYKColorspace)
3062                     Minimize(min.black,(double)
3063                       GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3064                   Minimize(min.alpha,(double)
3065                     GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3066                 }
3067                 else if ( (*k) < 0.3 )
3068                 { /* maximum of background pixels */
3069                   Maximize(max.red,     (double)
3070                     GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3071                   Maximize(max.green,   (double)
3072                     GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3073                   Maximize(max.blue,    (double)
3074                     GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3075                   if (image->colorspace == CMYKColorspace)
3076                     Maximize(max.black,   (double)
3077                       GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3078                   Maximize(max.alpha,(double)
3079                     GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3080                 }
3081               }
3082               k_pixels += virt_width*GetPixelChannels(image);
3083             }
3084             /* Pattern Match if difference is positive */
3085             min.red     -= max.red;     Maximize( min.red,     0.0 );
3086             min.green   -= max.green;   Maximize( min.green,   0.0 );
3087             min.blue    -= max.blue;    Maximize( min.blue,    0.0 );
3088             min.black   -= max.black;   Maximize( min.black,   0.0 );
3089             min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3090             break;
3091
3092         case ErodeIntensityMorphology:
3093             /* Select Pixel with Minimum Intensity within kernel neighbourhood
3094             **
3095             ** WARNING: the intensity test fails for CMYK and does not
3096             ** take into account the moderating effect of the alpha channel
3097             ** on the intensity.
3098             **
3099             ** NOTE that the kernel is not reflected for this operation!
3100             */
3101             k = kernel->values;
3102             k_pixels = p;
3103             for (v=0; v < (ssize_t) kernel->height; v++) {
3104               for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3105                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3106                 if ( result.red == 0.0 ||
3107                      GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) < GetPixelIntensity(morphology_image,q) ) {
3108                   /* copy the whole pixel - no channel selection */
3109                   SetPixelRed(morphology_image,GetPixelRed(image,
3110                     k_pixels+u*GetPixelChannels(image)),q);
3111                   SetPixelGreen(morphology_image,GetPixelGreen(image,
3112                     k_pixels+u*GetPixelChannels(image)),q);
3113                   SetPixelBlue(morphology_image,GetPixelBlue(image,
3114                     k_pixels+u*GetPixelChannels(image)),q);
3115                   SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3116                     k_pixels+u*GetPixelChannels(image)),q);
3117                   if ( result.red > 0.0 ) changed++;
3118                   result.red = 1.0;
3119                 }
3120               }
3121               k_pixels += virt_width*GetPixelChannels(image);
3122             }
3123             break;
3124
3125         case DilateIntensityMorphology:
3126             /* Select Pixel with Maximum Intensity within kernel neighbourhood
3127             **
3128             ** WARNING: the intensity test fails for CMYK and does not
3129             ** take into account the moderating effect of the alpha channel
3130             ** on the intensity (yet).
3131             **
3132             ** NOTE for correct working of this operation for asymetrical
3133             ** kernels, the kernel needs to be applied in its reflected form.
3134             ** That is its values needs to be reversed.
3135             */
3136             k = &kernel->values[ kernel->width*kernel->height-1 ];
3137             k_pixels = p;
3138             for (v=0; v < (ssize_t) kernel->height; v++) {
3139               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3140                 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3141                 if ( result.red == 0.0 ||
3142                      GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) > GetPixelIntensity(morphology_image,q) ) {
3143                   /* copy the whole pixel - no channel selection */
3144                   SetPixelRed(morphology_image,GetPixelRed(image,
3145                     k_pixels+u*GetPixelChannels(image)),q);
3146                   SetPixelGreen(morphology_image,GetPixelGreen(image,
3147                     k_pixels+u*GetPixelChannels(image)),q);
3148                   SetPixelBlue(morphology_image,GetPixelBlue(image,
3149                     k_pixels+u*GetPixelChannels(image)),q);
3150                   SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3151                     k_pixels+u*GetPixelChannels(image)),q);
3152                   if ( result.red > 0.0 ) changed++;
3153                   result.red = 1.0;
3154                 }
3155               }
3156               k_pixels += virt_width*GetPixelChannels(image);
3157             }
3158             break;
3159 #if 0
3160   This code has been obsoleted by the MorphologyPrimitiveDirect() function.
3161   However it is still (almost) correct coding for Grayscale Morphology.
3162   That is...
3163
3164   GrayErode    is equivalent but with kernel values subtracted from pixels
3165                without the kernel rotation
3166   GreyDilate   is equivalent but using Maximum() instead of Minimum()
3167                using kernel rotation
3168
3169   It has thus been preserved for future implementation of those methods.
3170
3171         case DistanceMorphology:
3172             /* Add kernel Value and select the minimum value found.
3173             ** The result is a iterative distance from edge of image shape.
3174             **
3175             ** All Distance Kernels are symetrical, but that may not always
3176             ** be the case. For example how about a distance from left edges?
3177             ** To work correctly with asymetrical kernels the reflected kernel
3178             ** needs to be applied.
3179             */
3180             k = &kernel->values[ kernel->width*kernel->height-1 ];
3181             k_pixels = p;
3182             for (v=0; v < (ssize_t) kernel->height; v++) {
3183               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3184                 if ( IsNan(*k) ) continue;
3185                 Minimize(result.red,     (*k)+k_pixels[u].red);
3186                 Minimize(result.green,   (*k)+k_pixels[u].green);
3187                 Minimize(result.blue,    (*k)+k_pixels[u].blue);
3188                 Minimize(result.alpha, (*k)+k_pixels[u].alpha);
3189                 if ( image->colorspace == CMYKColorspace)
3190                   Minimize(result.black,(*k)+GetPixelBlack(p_image,k_indexes+u));
3191               }
3192               k_pixels += virt_width*GetPixelChannels(image);
3193             }
3194             break;
3195 #endif
3196         case UndefinedMorphology:
3197         default:
3198             break; /* Do nothing */
3199       }
3200       /* Final mathematics of results (combine with original image?)
3201       **
3202       ** NOTE: Difference Morphology operators Edge* and *Hat could also
3203       ** be done here but works better with iteration as a image difference
3204       ** in the controling function (below).  Thicken and Thinning however
3205       ** should be done here so thay can be iterated correctly.
3206       */
3207       switch ( method ) {
3208         case HitAndMissMorphology:
3209         case ErodeMorphology:
3210           result = min;    /* minimum of neighbourhood */
3211           break;
3212         case DilateMorphology:
3213           result = max;    /* maximum of neighbourhood */
3214           break;
3215         case ThinningMorphology:
3216           /* subtract pattern match from original */
3217           result.red     -= min.red;
3218           result.green   -= min.green;
3219           result.blue    -= min.blue;
3220           result.black   -= min.black;
3221           result.alpha -= min.alpha;
3222           break;
3223         case ThickenMorphology:
3224           /* Add the pattern matchs to the original */
3225           result.red     += min.red;
3226           result.green   += min.green;
3227           result.blue    += min.blue;
3228           result.black   += min.black;
3229           result.alpha += min.alpha;
3230           break;
3231         default:
3232           /* result directly calculated or assigned */
3233           break;
3234       }
3235       /* Assign the resulting pixel values - Clamping Result */
3236       switch ( method ) {
3237         case UndefinedMorphology:
3238         case ConvolveMorphology:
3239         case DilateIntensityMorphology:
3240         case ErodeIntensityMorphology:
3241           break;  /* full pixel was directly assigned - not a channel method */
3242         default:
3243           if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3244             SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
3245           if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3246             SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
3247           if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3248             SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
3249           if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3250               (image->colorspace == CMYKColorspace))
3251             SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
3252           if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3253               (image->matte == MagickTrue))
3254             SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3255           break;
3256       }
3257       /* Count up changed pixels */
3258       if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q)) ||
3259           (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q)) ||
3260           (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q)) ||
3261           (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q)) ||
3262           ((image->colorspace == CMYKColorspace) &&
3263            (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
3264         changed++;  /* The pixel was changed in some way! */
3265       p+=GetPixelChannels(image);
3266       q+=GetPixelChannels(morphology_image);
3267     } /* x */
3268     if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3269       status=MagickFalse;
3270     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3271       {
3272         MagickBooleanType
3273           proceed;
3274
3275 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3276   #pragma omp critical (MagickCore_MorphologyImage)
3277 #endif
3278         proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3279         if (proceed == MagickFalse)
3280           status=MagickFalse;
3281       }
3282   } /* y */
3283   morphology_view=DestroyCacheView(morphology_view);
3284   image_view=DestroyCacheView(image_view);
3285   return(status ? (ssize_t)changed : -1);
3286 }
3287
3288 /* This is almost identical to the MorphologyPrimative() function above,
3289 ** but will apply the primitive directly to the image in two passes.
3290 **
3291 ** That is after each row is 'Sync'ed' into the image, the next row will
3292 ** make use of those values as part of the calculation of the next row.
3293 ** It then repeats, but going in the oppisite (bottom-up) direction.
3294 **
3295 ** Because of this 'iterative' handling this function can not make use
3296 ** of multi-threaded, parellel processing.
3297 */
3298 static ssize_t MorphologyPrimitiveDirect(Image *image,
3299   const MorphologyMethod method,const KernelInfo *kernel,
3300   ExceptionInfo *exception)
3301 {
3302   CacheView
3303     *auth_view,
3304     *virt_view;
3305
3306   MagickBooleanType
3307     status;
3308
3309   MagickOffsetType
3310     progress;
3311
3312   ssize_t
3313     y, offx, offy;
3314
3315   size_t
3316     virt_width,
3317     changed;
3318
3319   status=MagickTrue;
3320   changed=0;
3321   progress=0;
3322
3323   assert(image != (Image *) NULL);
3324   assert(image->signature == MagickSignature);
3325   assert(kernel != (KernelInfo *) NULL);
3326   assert(kernel->signature == MagickSignature);
3327   assert(exception != (ExceptionInfo *) NULL);
3328   assert(exception->signature == MagickSignature);
3329
3330   /* Some methods (including convolve) needs use a reflected kernel.
3331    * Adjust 'origin' offsets to loop though kernel as a reflection.
3332    */
3333   offx = kernel->x;
3334   offy = kernel->y;
3335   switch(method) {
3336     case DistanceMorphology:
3337     case VoronoiMorphology:
3338       /* kernel needs to used with reflection about origin */
3339       offx = (ssize_t) kernel->width-offx-1;
3340       offy = (ssize_t) kernel->height-offy-1;
3341       break;
3342 #if 0
3343     case ?????Morphology:
3344       /* kernel is used as is, without reflection */
3345       break;
3346 #endif
3347     default:
3348       assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
3349       break;
3350   }
3351
3352   /* DO NOT THREAD THIS CODE! */
3353   /* two views into same image (virtual, and actual) */
3354   virt_view=AcquireCacheView(image);
3355   auth_view=AcquireCacheView(image);
3356   virt_width=image->columns+kernel->width-1;
3357
3358   for (y=0; y < (ssize_t) image->rows; y++)
3359   {
3360     register const Quantum
3361       *restrict p;
3362
3363     register Quantum
3364       *restrict q;
3365
3366     register ssize_t
3367       x;
3368
3369     ssize_t
3370       r;
3371
3372     /* NOTE read virtual pixels, and authentic pixels, from the same image!
3373     ** we read using virtual to get virtual pixel handling, but write back
3374     ** into the same image.
3375     **
3376     ** Only top half of kernel is processed as we do a single pass downward
3377     ** through the image iterating the distance function as we go.
3378     */
3379     if (status == MagickFalse)
3380       break;
3381     p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
3382       offy+1,exception);
3383     q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3384       exception);
3385     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3386       status=MagickFalse;
3387     if (status == MagickFalse)
3388       break;
3389
3390     /* offset to origin in 'p'. while 'q' points to it directly */
3391     r = (ssize_t) virt_width*offy + offx;
3392
3393     for (x=0; x < (ssize_t) image->columns; x++)
3394     {
3395       ssize_t
3396         v;
3397
3398       register ssize_t
3399         u;
3400
3401       register const double
3402         *restrict k;
3403
3404       register const Quantum
3405         *restrict k_pixels;
3406
3407       PixelInfo
3408         result;
3409
3410       /* Starting Defaults */
3411       GetPixelInfo(image,&result);
3412       SetPixelInfo(image,q,&result);
3413       if ( method != VoronoiMorphology )
3414         result.alpha = QuantumRange - result.alpha;
3415
3416       switch ( method ) {
3417         case DistanceMorphology:
3418             /* Add kernel Value and select the minimum value found. */
3419             k = &kernel->values[ kernel->width*kernel->height-1 ];
3420             k_pixels = p;
3421             for (v=0; v <= (ssize_t) offy; v++) {
3422               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3423                 if ( IsNan(*k) ) continue;
3424                 Minimize(result.red,     (*k)+
3425                   GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3426                 Minimize(result.green,   (*k)+
3427                   GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3428                 Minimize(result.blue,    (*k)+
3429                   GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3430                 if (image->colorspace == CMYKColorspace)
3431                   Minimize(result.black,(*k)+
3432                     GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3433                 Minimize(result.alpha, (*k)+
3434                   GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3435               }
3436               k_pixels += virt_width*GetPixelChannels(image);
3437             }
3438             /* repeat with the just processed pixels of this row */
3439             k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3440             k_pixels = q-offx*GetPixelChannels(image);
3441               for (u=0; u < (ssize_t) offx; u++, k--) {
3442                 if ( x+u-offx < 0 ) continue;  /* off the edge! */
3443                 if ( IsNan(*k) ) continue;
3444                 Minimize(result.red,     (*k)+
3445                   GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3446                 Minimize(result.green,   (*k)+
3447                   GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3448                 Minimize(result.blue,    (*k)+
3449                   GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3450                 if (image->colorspace == CMYKColorspace)
3451                   Minimize(result.black,(*k)+
3452                     GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3453                 Minimize(result.alpha,(*k)+
3454                   GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3455               }
3456             break;
3457         case VoronoiMorphology:
3458             /* Apply Distance to 'Matte' channel, coping the closest color.
3459             **
3460             ** This is experimental, and realy the 'alpha' component should
3461             ** be completely separate 'masking' channel.
3462             */
3463             k = &kernel->values[ kernel->width*kernel->height-1 ];
3464             k_pixels = p;
3465             for (v=0; v <= (ssize_t) offy; v++) {
3466               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3467                 if ( IsNan(*k) ) continue;
3468                 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3469                   {
3470                     SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3471                       &result);
3472                     result.alpha += *k;
3473                   }
3474               }
3475               k_pixels += virt_width*GetPixelChannels(image);
3476             }
3477             /* repeat with the just processed pixels of this row */
3478             k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3479             k_pixels = q-offx*GetPixelChannels(image);
3480               for (u=0; u < (ssize_t) offx; u++, k--) {
3481                 if ( x+u-offx < 0 ) continue;  /* off the edge! */
3482                 if ( IsNan(*k) ) continue;
3483                 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3484                   {
3485                     SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3486                       &result);
3487                     result.alpha += *k;
3488                   }
3489               }
3490             break;
3491         default:
3492           /* result directly calculated or assigned */
3493           break;
3494       }
3495       /* Assign the resulting pixel values - Clamping Result */
3496       switch ( method ) {
3497         case VoronoiMorphology:
3498           SetPixelPixelInfo(image,&result,q);
3499           break;
3500         default:
3501           if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3502             SetPixelRed(image,ClampToQuantum(result.red),q);
3503           if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3504             SetPixelGreen(image,ClampToQuantum(result.green),q);
3505           if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3506             SetPixelBlue(image,ClampToQuantum(result.blue),q);
3507           if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3508               (image->colorspace == CMYKColorspace))
3509             SetPixelBlack(image,ClampToQuantum(result.black),q);
3510           if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3511               (image->matte == MagickTrue))
3512             SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3513           break;
3514       }
3515       /* Count up changed pixels */
3516       if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q)) ||
3517           (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q)) ||
3518           (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q)) ||
3519           (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q)) ||
3520           ((image->colorspace == CMYKColorspace) &&
3521            (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3522         changed++;  /* The pixel was changed in some way! */
3523
3524       p+=GetPixelChannels(image); /* increment pixel buffers */
3525       q+=GetPixelChannels(image);
3526     } /* x */
3527
3528     if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3529       status=MagickFalse;
3530     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3531       if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3532                 == MagickFalse )
3533         status=MagickFalse;
3534
3535   } /* y */
3536
3537   /* Do the reversed pass through the image */
3538   for (y=(ssize_t)image->rows-1; y >= 0; y--)
3539   {
3540     register const Quantum
3541       *restrict p;
3542
3543     register Quantum
3544       *restrict q;
3545
3546     register ssize_t
3547       x;
3548
3549     ssize_t
3550       r;
3551
3552     if (status == MagickFalse)
3553       break;
3554     /* NOTE read virtual pixels, and authentic pixels, from the same image!
3555     ** we read using virtual to get virtual pixel handling, but write back
3556     ** into the same image.
3557     **
3558     ** Only the bottom half of the kernel will be processes as we
3559     ** up the image.
3560     */
3561     p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
3562       kernel->y+1,exception);
3563     q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3564       exception);
3565     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3566       status=MagickFalse;
3567     if (status == MagickFalse)
3568       break;
3569
3570     /* adjust positions to end of row */
3571     p += (image->columns-1)*GetPixelChannels(image);
3572     q += (image->columns-1)*GetPixelChannels(image);
3573
3574     /* offset to origin in 'p'. while 'q' points to it directly */
3575     r = offx;
3576
3577     for (x=(ssize_t)image->columns-1; x >= 0; x--)
3578     {
3579       ssize_t
3580         v;
3581
3582       register ssize_t
3583         u;
3584
3585       register const double
3586         *restrict k;
3587
3588       register const Quantum
3589         *restrict k_pixels;
3590
3591       PixelInfo
3592         result;
3593
3594       /* Default - previously modified pixel */
3595       GetPixelInfo(image,&result);
3596       SetPixelInfo(image,q,&result);
3597       if ( method != VoronoiMorphology )
3598         result.alpha = QuantumRange - result.alpha;
3599
3600       switch ( method ) {
3601         case DistanceMorphology:
3602             /* Add kernel Value and select the minimum value found. */
3603             k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3604             k_pixels = p;
3605             for (v=offy; v < (ssize_t) kernel->height; v++) {
3606               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3607                 if ( IsNan(*k) ) continue;
3608                 Minimize(result.red,     (*k)+
3609                   GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3610                 Minimize(result.green,   (*k)+
3611                   GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3612                 Minimize(result.blue,    (*k)+
3613                   GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3614                 if ( image->colorspace == CMYKColorspace)
3615                   Minimize(result.black,(*k)+
3616                     GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3617                 Minimize(result.alpha, (*k)+
3618                   GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3619               }
3620               k_pixels += virt_width*GetPixelChannels(image);
3621             }
3622             /* repeat with the just processed pixels of this row */
3623             k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3624             k_pixels = q-offx*GetPixelChannels(image);
3625               for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3626                 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3627                 if ( IsNan(*k) ) continue;
3628                 Minimize(result.red,     (*k)+
3629                   GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3630                 Minimize(result.green,   (*k)+
3631                   GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3632                 Minimize(result.blue,    (*k)+
3633                   GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3634                 if ( image->colorspace == CMYKColorspace)
3635                   Minimize(result.black,   (*k)+
3636                     GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3637                 Minimize(result.alpha, (*k)+
3638                   GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3639               }
3640             break;
3641         case VoronoiMorphology:
3642             /* Apply Distance to 'Matte' channel, coping the closest color.
3643             **
3644             ** This is experimental, and realy the 'alpha' component should
3645             ** be completely separate 'masking' channel.
3646             */
3647             k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3648             k_pixels = p;
3649             for (v=offy; v < (ssize_t) kernel->height; v++) {
3650               for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3651                 if ( IsNan(*k) ) continue;
3652                 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3653                   {
3654                     SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3655                       &result);
3656                     result.alpha += *k;
3657                   }
3658               }
3659               k_pixels += virt_width*GetPixelChannels(image);
3660             }
3661             /* repeat with the just processed pixels of this row */
3662             k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3663             k_pixels = q-offx*GetPixelChannels(image);
3664               for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3665                 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3666                 if ( IsNan(*k) ) continue;
3667                 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3668                   {
3669                     SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3670                       &result);
3671                     result.alpha += *k;
3672                   }
3673               }
3674             break;
3675         default:
3676           /* result directly calculated or assigned */
3677           break;
3678       }
3679       /* Assign the resulting pixel values - Clamping Result */
3680       switch ( method ) {
3681         case VoronoiMorphology:
3682           SetPixelPixelInfo(image,&result,q);
3683           break;
3684         default:
3685           if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3686             SetPixelRed(image,ClampToQuantum(result.red),q);
3687           if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3688             SetPixelGreen(image,ClampToQuantum(result.green),q);
3689           if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3690             SetPixelBlue(image,ClampToQuantum(result.blue),q);
3691           if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3692               (image->colorspace == CMYKColorspace))
3693             SetPixelBlack(image,ClampToQuantum(result.black),q);
3694           if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3695               (image->matte == MagickTrue))
3696             SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3697           break;
3698       }
3699       /* Count up changed pixels */
3700       if (   (GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q))
3701           || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q))
3702           || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q))
3703           || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q))
3704           || ((image->colorspace == CMYKColorspace) &&
3705               (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3706         changed++;  /* The pixel was changed in some way! */
3707
3708       p-=GetPixelChannels(image); /* go backward through pixel buffers */
3709       q-=GetPixelChannels(image);
3710     } /* x */
3711     if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3712       status=MagickFalse;
3713     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3714       if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3715                 == MagickFalse )
3716         status=MagickFalse;
3717
3718   } /* y */
3719
3720   auth_view=DestroyCacheView(auth_view);
3721   virt_view=DestroyCacheView(virt_view);
3722   return(status ? (ssize_t) changed : -1);
3723 }
3724
3725 /* Apply a Morphology by calling theabove low level primitive application
3726 ** functions.  This function handles any iteration loops, composition or
3727 ** re-iteration of results, and compound morphology methods that is based
3728 ** on multiple low-level (staged) morphology methods.
3729 **
3730 ** Basically this provides the complex grue between the requested morphology
3731 ** method and raw low-level implementation (above).
3732 */
3733 MagickPrivate Image *MorphologyApply(const Image *image,
3734   const MorphologyMethod method, const ssize_t iterations,
3735   const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3736   ExceptionInfo *exception)
3737 {
3738   CompositeOperator
3739     curr_compose;
3740
3741   Image
3742     *curr_image,    /* Image we are working with or iterating */
3743     *work_image,    /* secondary image for primitive iteration */
3744     *save_image,    /* saved image - for 'edge' method only */
3745     *rslt_image;    /* resultant image - after multi-kernel handling */
3746
3747   KernelInfo
3748     *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3749     *norm_kernel,      /* the current normal un-reflected kernel */
3750     *rflt_kernel,      /* the current reflected kernel (if needed) */
3751     *this_kernel;      /* the kernel being applied */
3752
3753   MorphologyMethod
3754     primitive;      /* the current morphology primitive being applied */
3755
3756   CompositeOperator
3757     rslt_compose;   /* multi-kernel compose method for results to use */
3758
3759   MagickBooleanType
3760     special,        /* do we use a direct modify function? */
3761     verbose;        /* verbose output of results */
3762
3763   size_t
3764     method_loop,    /* Loop 1: number of compound method iterations (norm 1) */
3765     method_limit,   /*         maximum number of compound method iterations */
3766     kernel_number,  /* Loop 2: the kernel number being applied */
3767     stage_loop,     /* Loop 3: primitive loop for compound morphology */
3768     stage_limit,    /*         how many primitives are in this compound */
3769     kernel_loop,    /* Loop 4: iterate the kernel over image */
3770     kernel_limit,   /*         number of times to iterate kernel */
3771     count,          /* total count of primitive steps applied */
3772     kernel_changed, /* total count of changed using iterated kernel */
3773     method_changed; /* total count of changed over method iteration */
3774
3775   ssize_t
3776     changed;        /* number pixels changed by last primitive operation */
3777
3778   char
3779     v_info[80];
3780
3781   assert(image != (Image *) NULL);
3782   assert(image->signature == MagickSignature);
3783   assert(kernel != (KernelInfo *) NULL);
3784   assert(kernel->signature == MagickSignature);
3785   assert(exception != (ExceptionInfo *) NULL);
3786   assert(exception->signature == MagickSignature);
3787
3788   count = 0;      /* number of low-level morphology primitives performed */
3789   if ( iterations == 0 )
3790     return((Image *)NULL);   /* null operation - nothing to do! */
3791
3792   kernel_limit = (size_t) iterations;
3793   if ( iterations < 0 )  /* negative interations = infinite (well alomst) */
3794      kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3795
3796   verbose = IsMagickTrue(GetImageArtifact(image,"verbose"));
3797
3798   /* initialise for cleanup */
3799   curr_image = (Image *) image;
3800   curr_compose = image->compose;
3801   (void) curr_compose;
3802   work_image = save_image = rslt_image = (Image *) NULL;
3803   reflected_kernel = (KernelInfo *) NULL;
3804
3805   /* Initialize specific methods
3806    * + which loop should use the given iteratations
3807    * + how many primitives make up the compound morphology
3808    * + multi-kernel compose method to use (by default)
3809    */
3810   method_limit = 1;       /* just do method once, unless otherwise set */
3811   stage_limit = 1;        /* assume method is not a compound */
3812   special = MagickFalse;   /* assume it is NOT a direct modify primitive */
3813   rslt_compose = compose; /* and we are composing multi-kernels as given */
3814   switch( method ) {
3815     case SmoothMorphology:  /* 4 primitive compound morphology */
3816       stage_limit = 4;
3817       break;
3818     case OpenMorphology:    /* 2 primitive compound morphology */
3819     case OpenIntensityMorphology:
3820     case TopHatMorphology:
3821     case CloseMorphology:
3822     case CloseIntensityMorphology:
3823     case BottomHatMorphology:
3824     case EdgeMorphology:
3825       stage_limit = 2;
3826       break;
3827     case HitAndMissMorphology:
3828       rslt_compose = LightenCompositeOp;  /* Union of multi-kernel results */
3829       /* FALL THUR */
3830     case ThinningMorphology:
3831     case ThickenMorphology:
3832       method_limit = kernel_limit;  /* iterate the whole method */
3833       kernel_limit = 1;             /* do not do kernel iteration  */
3834       break;
3835     case DistanceMorphology:
3836     case VoronoiMorphology:
3837       special = MagickTrue;
3838       break;
3839     default:
3840       break;
3841   }
3842
3843   /* Apply special methods with special requirments
3844   ** For example, single run only, or post-processing requirements
3845   */
3846   if ( special == MagickTrue )
3847     {
3848       rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3849       if (rslt_image == (Image *) NULL)
3850         goto error_cleanup;
3851       if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3852         goto error_cleanup;
3853
3854       changed = MorphologyPrimitiveDirect(rslt_image, method,
3855          kernel, exception);
3856
3857       if ( verbose == MagickTrue )
3858         (void) (void) FormatLocaleFile(stderr,
3859           "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3860           CommandOptionToMnemonic(MagickMorphologyOptions, method),
3861           1.0,0.0,1.0, (double) changed);
3862
3863       if ( changed < 0 )
3864         goto error_cleanup;
3865
3866       if ( method == VoronoiMorphology ) {
3867         /* Preserve the alpha channel of input image - but turned off */
3868         (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3869           exception);
3870         (void) CompositeImage(rslt_image, CopyOpacityCompositeOp, image, 0, 0);
3871         (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3872           exception);
3873       }
3874       goto exit_cleanup;
3875     }
3876
3877   /* Handle user (caller) specified multi-kernel composition method */
3878   if ( compose != UndefinedCompositeOp )
3879     rslt_compose = compose;  /* override default composition for method */
3880   if ( rslt_compose == UndefinedCompositeOp )
3881     rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3882
3883   /* Some methods require a reflected kernel to use with primitives.
3884    * Create the reflected kernel for those methods. */
3885   switch ( method ) {
3886     case CorrelateMorphology:
3887     case CloseMorphology:
3888     case CloseIntensityMorphology:
3889     case BottomHatMorphology:
3890     case SmoothMorphology:
3891       reflected_kernel = CloneKernelInfo(kernel);
3892       if (reflected_kernel == (KernelInfo *) NULL)
3893         goto error_cleanup;
3894       RotateKernelInfo(reflected_kernel,180);
3895       break;
3896     default:
3897       break;
3898   }
3899
3900   /* Loops around more primitive morpholgy methods
3901   **  erose, dilate, open, close, smooth, edge, etc...
3902   */
3903   /* Loop 1:  iterate the compound method */
3904   method_loop = 0;
3905   method_changed = 1;
3906   while ( method_loop < method_limit && method_changed > 0 ) {
3907     method_loop++;
3908     method_changed = 0;
3909
3910     /* Loop 2:  iterate over each kernel in a multi-kernel list */
3911     norm_kernel = (KernelInfo *) kernel;
3912     this_kernel = (KernelInfo *) kernel;
3913     rflt_kernel = reflected_kernel;
3914
3915     kernel_number = 0;
3916     while ( norm_kernel != NULL ) {
3917
3918       /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3919       stage_loop = 0;          /* the compound morphology stage number */
3920       while ( stage_loop < stage_limit ) {
3921         stage_loop++;   /* The stage of the compound morphology */
3922
3923         /* Select primitive morphology for this stage of compound method */
3924         this_kernel = norm_kernel; /* default use unreflected kernel */
3925         primitive = method;        /* Assume method is a primitive */
3926         switch( method ) {
3927           case ErodeMorphology:      /* just erode */
3928           case EdgeInMorphology:     /* erode and image difference */
3929             primitive = ErodeMorphology;
3930             break;
3931           case DilateMorphology:     /* just dilate */
3932           case EdgeOutMorphology:    /* dilate and image difference */
3933             primitive = DilateMorphology;
3934             break;
3935           case OpenMorphology:       /* erode then dialate */
3936           case TopHatMorphology:     /* open and image difference */
3937             primitive = ErodeMorphology;
3938             if ( stage_loop == 2 )
3939               primitive = DilateMorphology;
3940             break;
3941           case OpenIntensityMorphology:
3942             primitive = ErodeIntensityMorphology;
3943             if ( stage_loop == 2 )
3944               primitive = DilateIntensityMorphology;
3945             break;
3946           case CloseMorphology:      /* dilate, then erode */
3947           case BottomHatMorphology:  /* close and image difference */
3948             this_kernel = rflt_kernel; /* use the reflected kernel */
3949             primitive = DilateMorphology;
3950             if ( stage_loop == 2 )
3951               primitive = ErodeMorphology;
3952             break;
3953           case CloseIntensityMorphology:
3954             this_kernel = rflt_kernel; /* use the reflected kernel */
3955             primitive = DilateIntensityMorphology;
3956             if ( stage_loop == 2 )
3957               primitive = ErodeIntensityMorphology;
3958             break;
3959           case SmoothMorphology:         /* open, close */
3960             switch ( stage_loop ) {
3961               case 1: /* start an open method, which starts with Erode */
3962                 primitive = ErodeMorphology;
3963                 break;
3964               case 2:  /* now Dilate the Erode */
3965                 primitive = DilateMorphology;
3966                 break;
3967               case 3:  /* Reflect kernel a close */
3968                 this_kernel = rflt_kernel; /* use the reflected kernel */
3969                 primitive = DilateMorphology;
3970                 break;
3971               case 4:  /* Finish the Close */
3972                 this_kernel = rflt_kernel; /* use the reflected kernel */
3973                 primitive = ErodeMorphology;
3974                 break;
3975             }
3976             break;
3977           case EdgeMorphology:        /* dilate and erode difference */
3978             primitive = DilateMorphology;
3979             if ( stage_loop == 2 ) {
3980               save_image = curr_image;      /* save the image difference */
3981               curr_image = (Image *) image;
3982               primitive = ErodeMorphology;
3983             }
3984             break;
3985           case CorrelateMorphology:
3986             /* A Correlation is a Convolution with a reflected kernel.
3987             ** However a Convolution is a weighted sum using a reflected
3988             ** kernel.  It may seem stange to convert a Correlation into a
3989             ** Convolution as the Correlation is the simplier method, but
3990             ** Convolution is much more commonly used, and it makes sense to
3991             ** implement it directly so as to avoid the need to duplicate the
3992             ** kernel when it is not required (which is typically the
3993             ** default).
3994             */
3995             this_kernel = rflt_kernel; /* use the reflected kernel */
3996             primitive = ConvolveMorphology;
3997             break;
3998           default:
3999             break;
4000         }
4001         assert( this_kernel != (KernelInfo *) NULL );
4002
4003         /* Extra information for debugging compound operations */
4004         if ( verbose == MagickTrue ) {
4005           if ( stage_limit > 1 )
4006             (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
4007              CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
4008              method_loop,(double) stage_loop);
4009           else if ( primitive != method )
4010             (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
4011               CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
4012               method_loop);
4013           else
4014             v_info[0] = '\0';
4015         }
4016
4017         /* Loop 4: Iterate the kernel with primitive */
4018         kernel_loop = 0;
4019         kernel_changed = 0;
4020         changed = 1;
4021         while ( kernel_loop < kernel_limit && changed > 0 ) {
4022           kernel_loop++;     /* the iteration of this kernel */
4023
4024           /* Create a clone as the destination image, if not yet defined */
4025           if ( work_image == (Image *) NULL )
4026             {
4027               work_image=CloneImage(image,0,0,MagickTrue,exception);
4028               if (work_image == (Image *) NULL)
4029                 goto error_cleanup;
4030               if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
4031                 goto error_cleanup;
4032               /* work_image->type=image->type; ??? */
4033             }
4034
4035           /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4036           count++;
4037           changed = MorphologyPrimitive(curr_image, work_image, primitive,
4038                        this_kernel, bias, exception);
4039
4040           if ( verbose == MagickTrue ) {
4041             if ( kernel_loop > 1 )
4042               (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4043             (void) (void) FormatLocaleFile(stderr,
4044               "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4045               v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4046               primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4047               (double) (method_loop+kernel_loop-1),(double) kernel_number,
4048               (double) count,(double) changed);
4049           }
4050           if ( changed < 0 )
4051             goto error_cleanup;
4052           kernel_changed += changed;
4053           method_changed += changed;
4054
4055           /* prepare next loop */
4056           { Image *tmp = work_image;   /* swap images for iteration */
4057             work_image = curr_image;
4058             curr_image = tmp;
4059           }
4060           if ( work_image == image )
4061             work_image = (Image *) NULL; /* replace input 'image' */
4062
4063         } /* End Loop 4: Iterate the kernel with primitive */
4064
4065         if ( verbose == MagickTrue && kernel_changed != (size_t)changed )
4066           (void) FormatLocaleFile(stderr, "   Total %.20g",(double) kernel_changed);
4067         if ( verbose == MagickTrue && stage_loop < stage_limit )
4068           (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4069
4070 #if 0
4071     (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4072     (void) FormatLocaleFile(stderr, "      curr =0x%lx\n", (unsigned long)curr_image);
4073     (void) FormatLocaleFile(stderr, "      work =0x%lx\n", (unsigned long)work_image);
4074     (void) FormatLocaleFile(stderr, "      save =0x%lx\n", (unsigned long)save_image);
4075     (void) FormatLocaleFile(stderr, "      union=0x%lx\n", (unsigned long)rslt_image);
4076 #endif
4077
4078       } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4079
4080       /*  Final Post-processing for some Compound Methods
4081       **
4082       ** The removal of any 'Sync' channel flag in the Image Compositon
4083       ** below ensures the methematical compose method is applied in a
4084       ** purely mathematical way, and only to the selected channels.
4085       ** Turn off SVG composition 'alpha blending'.
4086       */
4087       switch( method ) {
4088         case EdgeOutMorphology:
4089         case EdgeInMorphology:
4090         case TopHatMorphology:
4091         case BottomHatMorphology:
4092           if ( verbose == MagickTrue )
4093             (void) FormatLocaleFile(stderr, "\n%s: Difference with original image",
4094                  CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4095           curr_image->sync=MagickFalse;
4096           (void) CompositeImage(curr_image,DifferenceCompositeOp,image,0,0);
4097           break;
4098         case EdgeMorphology:
4099           if ( verbose == MagickTrue )
4100             (void) FormatLocaleFile(stderr, "\n%s: Difference of Dilate and Erode",
4101                  CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4102           curr_image->sync=MagickFalse;
4103           (void) CompositeImage(curr_image,DifferenceCompositeOp,save_image,0,
4104             0);
4105           save_image = DestroyImage(save_image); /* finished with save image */
4106           break;
4107         default:
4108           break;
4109       }
4110
4111       /* multi-kernel handling:  re-iterate, or compose results */
4112       if ( kernel->next == (KernelInfo *) NULL )
4113         rslt_image = curr_image;   /* just return the resulting image */
4114       else if ( rslt_compose == NoCompositeOp )
4115         { if ( verbose == MagickTrue ) {
4116             if ( this_kernel->next != (KernelInfo *) NULL )
4117               (void) FormatLocaleFile(stderr, " (re-iterate)");
4118             else
4119               (void) FormatLocaleFile(stderr, " (done)");
4120           }
4121           rslt_image = curr_image; /* return result, and re-iterate */
4122         }
4123       else if ( rslt_image == (Image *) NULL)
4124         { if ( verbose == MagickTrue )
4125             (void) FormatLocaleFile(stderr, " (save for compose)");
4126           rslt_image = curr_image;
4127           curr_image = (Image *) image;  /* continue with original image */
4128         }
4129       else
4130         { /* Add the new 'current' result to the composition
4131           **
4132           ** The removal of any 'Sync' channel flag in the Image Compositon
4133           ** below ensures the methematical compose method is applied in a
4134           ** purely mathematical way, and only to the selected channels.
4135           ** IE: Turn off SVG composition 'alpha blending'.
4136           */
4137           if ( verbose == MagickTrue )
4138             (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4139                  CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4140           rslt_image->sync=MagickFalse;
4141           (void) CompositeImage(rslt_image, rslt_compose, curr_image, 0, 0);
4142           curr_image = DestroyImage(curr_image);
4143           curr_image = (Image *) image;  /* continue with original image */
4144         }
4145       if ( verbose == MagickTrue )
4146         (void) FormatLocaleFile(stderr, "\n");
4147
4148       /* loop to the next kernel in a multi-kernel list */
4149       norm_kernel = norm_kernel->next;
4150       if ( rflt_kernel != (KernelInfo *) NULL )
4151         rflt_kernel = rflt_kernel->next;
4152       kernel_number++;
4153     } /* End Loop 2: Loop over each kernel */
4154
4155   } /* End Loop 1: compound method interation */
4156
4157   goto exit_cleanup;
4158
4159   /* Yes goto's are bad, but it makes cleanup lot more efficient */
4160 error_cleanup:
4161   if ( curr_image == rslt_image )
4162     curr_image = (Image *) NULL;
4163   if ( rslt_image != (Image *) NULL )
4164     rslt_image = DestroyImage(rslt_image);
4165 exit_cleanup:
4166   if ( curr_image == rslt_image || curr_image == image )
4167     curr_image = (Image *) NULL;
4168   if ( curr_image != (Image *) NULL )
4169     curr_image = DestroyImage(curr_image);
4170   if ( work_image != (Image *) NULL )
4171     work_image = DestroyImage(work_image);
4172   if ( save_image != (Image *) NULL )
4173     save_image = DestroyImage(save_image);
4174   if ( reflected_kernel != (KernelInfo *) NULL )
4175     reflected_kernel = DestroyKernelInfo(reflected_kernel);
4176   return(rslt_image);
4177 }
4178
4179 \f
4180 /*
4181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4182 %                                                                             %
4183 %                                                                             %
4184 %                                                                             %
4185 %     M o r p h o l o g y I m a g e                                           %
4186 %                                                                             %
4187 %                                                                             %
4188 %                                                                             %
4189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4190 %
4191 %  MorphologyImage() applies a user supplied kernel to the image
4192 %  according to the given mophology method.
4193 %
4194 %  This function applies any and all user defined settings before calling
4195 %  the above internal function MorphologyApply().
4196 %
4197 %  User defined settings include...
4198 %    * Output Bias for Convolution and correlation   ("-bias")
4199 %    * Kernel Scale/normalize settings     ("-set 'option:convolve:scale'")
4200 %      This can also includes the addition of a scaled unity kernel.
4201 %    * Show Kernel being applied           ("-set option:showkernel 1")
4202 %
4203 %  The format of the MorphologyImage method is:
4204 %
4205 %      Image *MorphologyImage(const Image *image,MorphologyMethod method,
4206 %        const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4207 %
4208 %      Image *MorphologyImage(const Image *image, const ChannelType
4209 %        channel,MorphologyMethod method,const ssize_t iterations,
4210 %        KernelInfo *kernel,ExceptionInfo *exception)
4211 %
4212 %  A description of each parameter follows:
4213 %
4214 %    o image: the image.
4215 %
4216 %    o method: the morphology method to be applied.
4217 %
4218 %    o iterations: apply the operation this many times (or no change).
4219 %                  A value of -1 means loop until no change found.
4220 %                  How this is applied may depend on the morphology method.
4221 %                  Typically this is a value of 1.
4222 %
4223 %    o kernel: An array of double representing the morphology kernel.
4224 %              Warning: kernel may be normalized for the Convolve method.
4225 %
4226 %    o exception: return any errors or warnings in this structure.
4227 %
4228 */
4229 MagickExport Image *MorphologyImage(const Image *image,
4230   const MorphologyMethod method,const ssize_t iterations,
4231   const KernelInfo *kernel,ExceptionInfo *exception)
4232 {
4233   KernelInfo
4234     *curr_kernel;
4235
4236   CompositeOperator
4237     compose;
4238
4239   Image
4240     *morphology_image;
4241
4242
4243   /* Apply Convolve/Correlate Normalization and Scaling Factors.
4244    * This is done BEFORE the ShowKernelInfo() function is called so that
4245    * users can see the results of the 'option:convolve:scale' option.
4246    */
4247   curr_kernel = (KernelInfo *) kernel;
4248   if ( method == ConvolveMorphology ||  method == CorrelateMorphology )
4249     {
4250       const char
4251         *artifact;
4252       artifact = GetImageArtifact(image,"convolve:scale");
4253       if ( artifact != (const char *)NULL ) {
4254         if ( curr_kernel == kernel )
4255           curr_kernel = CloneKernelInfo(kernel);
4256         if (curr_kernel == (KernelInfo *) NULL) {
4257           curr_kernel=DestroyKernelInfo(curr_kernel);
4258           return((Image *) NULL);
4259         }
4260         ScaleGeometryKernelInfo(curr_kernel, artifact);
4261       }
4262     }
4263
4264   /* display the (normalized) kernel via stderr */
4265   if ( IsMagickTrue(GetImageArtifact(image,"showkernel"))
4266     || IsMagickTrue(GetImageArtifact(image,"convolve:showkernel"))
4267     || IsMagickTrue(GetImageArtifact(image,"morphology:showkernel")) )
4268     ShowKernelInfo(curr_kernel);
4269
4270   /* Override the default handling of multi-kernel morphology results
4271    * If 'Undefined' use the default method
4272    * If 'None' (default for 'Convolve') re-iterate previous result
4273    * Otherwise merge resulting images using compose method given.
4274    * Default for 'HitAndMiss' is 'Lighten'.
4275    */
4276   { const char
4277       *artifact;
4278     artifact = GetImageArtifact(image,"morphology:compose");
4279     compose = UndefinedCompositeOp;  /* use default for method */
4280     if ( artifact != (const char *) NULL)
4281       compose=(CompositeOperator) ParseCommandOption(MagickComposeOptions,
4282         MagickFalse,artifact);
4283   }
4284   /* Apply the Morphology */
4285   morphology_image = MorphologyApply(image, method, iterations,
4286     curr_kernel, compose, image->bias, exception);
4287
4288   /* Cleanup and Exit */
4289   if ( curr_kernel != kernel )
4290     curr_kernel=DestroyKernelInfo(curr_kernel);
4291   return(morphology_image);
4292 }
4293 \f
4294 /*
4295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4296 %                                                                             %
4297 %                                                                             %
4298 %                                                                             %
4299 +     R o t a t e K e r n e l I n f o                                         %
4300 %                                                                             %
4301 %                                                                             %
4302 %                                                                             %
4303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4304 %
4305 %  RotateKernelInfo() rotates the kernel by the angle given.
4306 %
4307 %  Currently it is restricted to 90 degree angles, of either 1D kernels
4308 %  or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4309 %  It will ignore usless rotations for specific 'named' built-in kernels.
4310 %
4311 %  The format of the RotateKernelInfo method is:
4312 %
4313 %      void RotateKernelInfo(KernelInfo *kernel, double angle)
4314 %
4315 %  A description of each parameter follows:
4316 %
4317 %    o kernel: the Morphology/Convolution kernel
4318 %
4319 %    o angle: angle to rotate in degrees
4320 %
4321 % This function is currently internal to this module only, but can be exported
4322 % to other modules if needed.
4323 */
4324 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4325 {
4326   /* angle the lower kernels first */
4327   if ( kernel->next != (KernelInfo *) NULL)
4328     RotateKernelInfo(kernel->next, angle);
4329
4330   /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4331   **
4332   ** TODO: expand beyond simple 90 degree rotates, flips and flops
4333   */
4334
4335   /* Modulus the angle */
4336   angle = fmod(angle, 360.0);
4337   if ( angle < 0 )
4338     angle += 360.0;
4339
4340   if ( 337.5 < angle || angle <= 22.5 )
4341     return;   /* Near zero angle - no change! - At least not at this time */
4342
4343   /* Handle special cases */
4344   switch (kernel->type) {
4345     /* These built-in kernels are cylindrical kernels, rotating is useless */
4346     case GaussianKernel:
4347     case DoGKernel:
4348     case LoGKernel:
4349     case DiskKernel:
4350     case PeaksKernel:
4351     case LaplacianKernel:
4352     case ChebyshevKernel:
4353     case ManhattanKernel:
4354     case EuclideanKernel:
4355       return;
4356
4357     /* These may be rotatable at non-90 angles in the future */
4358     /* but simply rotating them in multiples of 90 degrees is useless */
4359     case SquareKernel:
4360     case DiamondKernel:
4361     case PlusKernel:
4362     case CrossKernel:
4363       return;
4364
4365     /* These only allows a +/-90 degree rotation (by transpose) */
4366     /* A 180 degree rotation is useless */
4367     case BlurKernel:
4368       if ( 135.0 < angle && angle <= 225.0 )
4369         return;
4370       if ( 225.0 < angle && angle <= 315.0 )
4371         angle -= 180;
4372       break;
4373
4374     default:
4375       break;
4376   }
4377   /* Attempt rotations by 45 degrees  -- 3x3 kernels only */
4378   if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4379     {
4380       if ( kernel->width == 3 && kernel->height == 3 )
4381         { /* Rotate a 3x3 square by 45 degree angle */
4382           MagickRealType t  = kernel->values[0];
4383           kernel->values[0] = kernel->values[3];
4384           kernel->values[3] = kernel->values[6];
4385           kernel->values[6] = kernel->values[7];
4386           kernel->values[7] = kernel->values[8];
4387           kernel->values[8] = kernel->values[5];
4388           kernel->values[5] = kernel->values[2];
4389           kernel->values[2] = kernel->values[1];
4390           kernel->values[1] = t;
4391           /* rotate non-centered origin */
4392           if ( kernel->x != 1 || kernel->y != 1 ) {
4393             ssize_t x,y;
4394             x = (ssize_t) kernel->x-1;
4395             y = (ssize_t) kernel->y-1;
4396                  if ( x == y  ) x = 0;
4397             else if ( x == 0  ) x = -y;
4398             else if ( x == -y ) y = 0;
4399             else if ( y == 0  ) y = x;
4400             kernel->x = (ssize_t) x+1;
4401             kernel->y = (ssize_t) y+1;
4402           }
4403           angle = fmod(angle+315.0, 360.0);  /* angle reduced 45 degrees */
4404           kernel->angle = fmod(kernel->angle+45.0, 360.0);
4405         }
4406       else
4407         perror("Unable to rotate non-3x3 kernel by 45 degrees");
4408     }
4409   if ( 45.0 < fmod(angle, 180.0)  && fmod(angle,180.0) <= 135.0 )
4410     {
4411       if ( kernel->width == 1 || kernel->height == 1 )
4412         { /* Do a transpose of a 1 dimensional kernel,
4413           ** which results in a fast 90 degree rotation of some type.
4414           */
4415           ssize_t
4416             t;
4417           t = (ssize_t) kernel->width;
4418           kernel->width = kernel->height;
4419           kernel->height = (size_t) t;
4420           t = kernel->x;
4421           kernel->x = kernel->y;
4422           kernel->y = t;
4423           if ( kernel->width == 1 ) {
4424             angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
4425             kernel->angle = fmod(kernel->angle+90.0, 360.0);
4426           } else {
4427             angle = fmod(angle+90.0, 360.0);   /* angle increased 90 degrees */
4428             kernel->angle = fmod(kernel->angle+270.0, 360.0);
4429           }
4430         }
4431       else if ( kernel->width == kernel->height )
4432         { /* Rotate a square array of values by 90 degrees */
4433           { register size_t
4434               i,j,x,y;
4435             register MagickRealType
4436               *k,t;
4437             k=kernel->values;
4438             for( i=0, x=kernel->width-1;  i<=x;   i++, x--)
4439               for( j=0, y=kernel->height-1;  j<y;   j++, y--)
4440                 { t                    = k[i+j*kernel->width];
4441                   k[i+j*kernel->width] = k[j+x*kernel->width];
4442                   k[j+x*kernel->width] = k[x+y*kernel->width];
4443                   k[x+y*kernel->width] = k[y+i*kernel->width];
4444                   k[y+i*kernel->width] = t;
4445                 }
4446           }
4447           /* rotate the origin - relative to center of array */
4448           { register ssize_t x,y;
4449             x = (ssize_t) (kernel->x*2-kernel->width+1);
4450             y = (ssize_t) (kernel->y*2-kernel->height+1);
4451             kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4452             kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4453           }
4454           angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
4455           kernel->angle = fmod(kernel->angle+90.0, 360.0);
4456         }
4457       else
4458         perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4459     }
4460   if ( 135.0 < angle && angle <= 225.0 )
4461     {
4462       /* For a 180 degree rotation - also know as a reflection
4463        * This is actually a very very common operation!
4464        * Basically all that is needed is a reversal of the kernel data!
4465        * And a reflection of the origon
4466        */
4467       size_t
4468         i,j;
4469       register double
4470         *k,t;
4471
4472       k=kernel->values;
4473       for ( i=0, j=kernel->width*kernel->height-1;  i<j;  i++, j--)
4474         t=k[i],  k[i]=k[j],  k[j]=t;
4475
4476       kernel->x = (ssize_t) kernel->width  - kernel->x - 1;
4477       kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4478       angle = fmod(angle-180.0, 360.0);   /* angle+180 degrees */
4479       kernel->angle = fmod(kernel->angle+180.0, 360.0);
4480     }
4481   /* At this point angle should at least between -45 (315) and +45 degrees
4482    * In the future some form of non-orthogonal angled rotates could be
4483    * performed here, posibily with a linear kernel restriction.
4484    */
4485
4486   return;
4487 }
4488 \f
4489 /*
4490 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4491 %                                                                             %
4492 %                                                                             %
4493 %                                                                             %
4494 %     S c a l e G e o m e t r y K e r n e l I n f o                           %
4495 %                                                                             %
4496 %                                                                             %
4497 %                                                                             %
4498 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4499 %
4500 %  ScaleGeometryKernelInfo() takes a geometry argument string, typically
4501 %  provided as a  "-set option:convolve:scale {geometry}" user setting,
4502 %  and modifies the kernel according to the parsed arguments of that setting.
4503 %
4504 %  The first argument (and any normalization flags) are passed to
4505 %  ScaleKernelInfo() to scale/normalize the kernel.  The second argument
4506 %  is then passed to UnityAddKernelInfo() to add a scled unity kernel
4507 %  into the scaled/normalized kernel.
4508 %
4509 %  The format of the ScaleGeometryKernelInfo method is:
4510 %
4511 %      void ScaleGeometryKernelInfo(KernelInfo *kernel,
4512 %        const double scaling_factor,const MagickStatusType normalize_flags)
4513 %
4514 %  A description of each parameter follows:
4515 %
4516 %    o kernel: the Morphology/Convolution kernel to modify
4517 %
4518 %    o geometry:
4519 %             The geometry string to parse, typically from the user provided
4520 %             "-set option:convolve:scale {geometry}" setting.
4521 %
4522 */
4523 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4524      const char *geometry)
4525 {
4526   GeometryFlags
4527     flags;
4528   GeometryInfo
4529     args;
4530
4531   SetGeometryInfo(&args);
4532   flags = (GeometryFlags) ParseGeometry(geometry, &args);
4533
4534 #if 0
4535   /* For Debugging Geometry Input */
4536   (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4537        flags, args.rho, args.sigma, args.xi, args.psi );
4538 #endif
4539
4540   if ( (flags & PercentValue) != 0 )      /* Handle Percentage flag*/
4541     args.rho *= 0.01,  args.sigma *= 0.01;
4542
4543   if ( (flags & RhoValue) == 0 )          /* Set Defaults for missing args */
4544     args.rho = 1.0;
4545   if ( (flags & SigmaValue) == 0 )
4546     args.sigma = 0.0;
4547
4548   /* Scale/Normalize the input kernel */
4549   ScaleKernelInfo(kernel, args.rho, flags);
4550
4551   /* Add Unity Kernel, for blending with original */
4552   if ( (flags & SigmaValue) != 0 )
4553     UnityAddKernelInfo(kernel, args.sigma);
4554
4555   return;
4556 }
4557 /*
4558 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4559 %                                                                             %
4560 %                                                                             %
4561 %                                                                             %
4562 %     S c a l e K e r n e l I n f o                                           %
4563 %                                                                             %
4564 %                                                                             %
4565 %                                                                             %
4566 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4567 %
4568 %  ScaleKernelInfo() scales the given kernel list by the given amount, with or
4569 %  without normalization of the sum of the kernel values (as per given flags).
4570 %
4571 %  By default (no flags given) the values within the kernel is scaled
4572 %  directly using given scaling factor without change.
4573 %
4574 %  If either of the two 'normalize_flags' are given the kernel will first be
4575 %  normalized and then further scaled by the scaling factor value given.
4576 %
4577 %  Kernel normalization ('normalize_flags' given) is designed to ensure that
4578 %  any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4579 %  morphology methods will fall into -1.0 to +1.0 range.  Note that for
4580 %  non-HDRI versions of IM this may cause images to have any negative results
4581 %  clipped, unless some 'bias' is used.
4582 %
4583 %  More specifically.  Kernels which only contain positive values (such as a
4584 %  'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4585 %  ensuring a 0.0 to +1.0 output range for non-HDRI images.
4586 %
4587 %  For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4588 %  the kernel will be scaled by the absolute of the sum of kernel values, so
4589 %  that it will generally fall within the +/- 1.0 range.
4590 %
4591 %  For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4592 %  will be scaled by just the sum of the postive values, so that its output
4593 %  range will again fall into the  +/- 1.0 range.
4594 %
4595 %  For special kernels designed for locating shapes using 'Correlate', (often
4596 %  only containing +1 and -1 values, representing foreground/brackground
4597 %  matching) a special normalization method is provided to scale the positive
4598 %  values separately to those of the negative values, so the kernel will be
4599 %  forced to become a zero-sum kernel better suited to such searches.
4600 %
4601 %  WARNING: Correct normalization of the kernel assumes that the '*_range'
4602 %  attributes within the kernel structure have been correctly set during the
4603 %  kernels creation.
4604 %
4605 %  NOTE: The values used for 'normalize_flags' have been selected specifically
4606 %  to match the use of geometry options, so that '!' means NormalizeValue, '^'
4607 %  means CorrelateNormalizeValue.  All other GeometryFlags values are ignored.
4608 %
4609 %  The format of the ScaleKernelInfo method is:
4610 %
4611 %      void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4612 %               const MagickStatusType normalize_flags )
4613 %
4614 %  A description of each parameter follows:
4615 %
4616 %    o kernel: the Morphology/Convolution kernel
4617 %
4618 %    o scaling_factor:
4619 %             multiply all values (after normalization) by this factor if not
4620 %             zero.  If the kernel is normalized regardless of any flags.
4621 %
4622 %    o normalize_flags:
4623 %             GeometryFlags defining normalization method to use.
4624 %             specifically: NormalizeValue, CorrelateNormalizeValue,
4625 %                           and/or PercentValue
4626 %
4627 */
4628 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4629   const double scaling_factor,const GeometryFlags normalize_flags)
4630 {
4631   register ssize_t
4632     i;
4633
4634   register double
4635     pos_scale,
4636     neg_scale;
4637
4638   /* do the other kernels in a multi-kernel list first */
4639   if ( kernel->next != (KernelInfo *) NULL)
4640     ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4641
4642   /* Normalization of Kernel */
4643   pos_scale = 1.0;
4644   if ( (normalize_flags&NormalizeValue) != 0 ) {
4645     if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4646       /* non-zero-summing kernel (generally positive) */
4647       pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4648     else
4649       /* zero-summing kernel */
4650       pos_scale = kernel->positive_range;
4651   }
4652   /* Force kernel into a normalized zero-summing kernel */
4653   if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4654     pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4655                  ? kernel->positive_range : 1.0;
4656     neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
4657                  ? -kernel->negative_range : 1.0;
4658   }
4659   else
4660     neg_scale = pos_scale;
4661
4662   /* finialize scaling_factor for positive and negative components */
4663   pos_scale = scaling_factor/pos_scale;
4664   neg_scale = scaling_factor/neg_scale;
4665
4666   for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4667     if ( ! IsNan(kernel->values[i]) )
4668       kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4669
4670   /* convolution output range */
4671   kernel->positive_range *= pos_scale;
4672   kernel->negative_range *= neg_scale;
4673   /* maximum and minimum values in kernel */
4674   kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4675   kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4676
4677   /* swap kernel settings if user's scaling factor is negative */
4678   if ( scaling_factor < MagickEpsilon ) {
4679     double t;
4680     t = kernel->positive_range;
4681     kernel->positive_range = kernel->negative_range;
4682     kernel->negative_range = t;
4683     t = kernel->maximum;
4684     kernel->maximum = kernel->minimum;
4685     kernel->minimum = 1;
4686   }
4687
4688   return;
4689 }
4690 \f
4691 /*
4692 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4693 %                                                                             %
4694 %                                                                             %
4695 %                                                                             %
4696 %     S h o w K e r n e l I n f o                                             %
4697 %                                                                             %
4698 %                                                                             %
4699 %                                                                             %
4700 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4701 %
4702 %  ShowKernelInfo() outputs the details of the given kernel defination to
4703 %  standard error, generally due to a users 'showkernel' option request.
4704 %
4705 %  The format of the ShowKernel method is:
4706 %
4707 %      void ShowKernelInfo(const KernelInfo *kernel)
4708 %
4709 %  A description of each parameter follows:
4710 %
4711 %    o kernel: the Morphology/Convolution kernel
4712 %
4713 */
4714 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4715 {
4716   const KernelInfo
4717     *k;
4718
4719   size_t
4720     c, i, u, v;
4721
4722   for (c=0, k=kernel;  k != (KernelInfo *) NULL;  c++, k=k->next ) {
4723
4724     (void) FormatLocaleFile(stderr, "Kernel");
4725     if ( kernel->next != (KernelInfo *) NULL )
4726       (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4727     (void) FormatLocaleFile(stderr, " \"%s",
4728           CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4729     if ( fabs(k->angle) > MagickEpsilon )
4730       (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4731     (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4732       k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4733     (void) FormatLocaleFile(stderr,
4734           " with values from %.*lg to %.*lg\n",
4735           GetMagickPrecision(), k->minimum,
4736           GetMagickPrecision(), k->maximum);
4737     (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4738           GetMagickPrecision(), k->negative_range,
4739           GetMagickPrecision(), k->positive_range);
4740     if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4741       (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4742     else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4743       (void) FormatLocaleFile(stderr, " (Normalized)\n");
4744     else
4745       (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4746           GetMagickPrecision(), k->positive_range+k->negative_range);
4747     for (i=v=0; v < k->height; v++) {
4748       (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4749       for (u=0; u < k->width; u++, i++)
4750         if ( IsNan(k->values[i]) )
4751           (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4752         else
4753           (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4754               GetMagickPrecision(), k->values[i]);
4755       (void) FormatLocaleFile(stderr,"\n");
4756     }
4757   }
4758 }
4759 \f
4760 /*
4761 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4762 %                                                                             %
4763 %                                                                             %
4764 %                                                                             %
4765 %     U n i t y A d d K e r n a l I n f o                                     %
4766 %                                                                             %
4767 %                                                                             %
4768 %                                                                             %
4769 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4770 %
4771 %  UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4772 %  to the given pre-scaled and normalized Kernel.  This in effect adds that
4773 %  amount of the original image into the resulting convolution kernel.  This
4774 %  value is usually provided by the user as a percentage value in the
4775 %  'convolve:scale' setting.
4776 %
4777 %  The resulting effect is to convert the defined kernels into blended
4778 %  soft-blurs, unsharp kernels or into sharpening kernels.
4779 %
4780 %  The format of the UnityAdditionKernelInfo method is:
4781 %
4782 %      void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4783 %
4784 %  A description of each parameter follows:
4785 %
4786 %    o kernel: the Morphology/Convolution kernel
4787 %
4788 %    o scale:
4789 %             scaling factor for the unity kernel to be added to
4790 %             the given kernel.
4791 %
4792 */
4793 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4794   const double scale)
4795 {
4796   /* do the other kernels in a multi-kernel list first */
4797   if ( kernel->next != (KernelInfo *) NULL)
4798     UnityAddKernelInfo(kernel->next, scale);
4799
4800   /* Add the scaled unity kernel to the existing kernel */
4801   kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4802   CalcKernelMetaData(kernel);  /* recalculate the meta-data */
4803
4804   return;
4805 }
4806 \f
4807 /*
4808 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4809 %                                                                             %
4810 %                                                                             %
4811 %                                                                             %
4812 %     Z e r o K e r n e l N a n s                                             %
4813 %                                                                             %
4814 %                                                                             %
4815 %                                                                             %
4816 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4817 %
4818 %  ZeroKernelNans() replaces any special 'nan' value that may be present in
4819 %  the kernel with a zero value.  This is typically done when the kernel will
4820 %  be used in special hardware (GPU) convolution processors, to simply
4821 %  matters.
4822 %
4823 %  The format of the ZeroKernelNans method is:
4824 %
4825 %      void ZeroKernelNans (KernelInfo *kernel)
4826 %
4827 %  A description of each parameter follows:
4828 %
4829 %    o kernel: the Morphology/Convolution kernel
4830 %
4831 */
4832 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4833 {
4834   register size_t
4835     i;
4836
4837   /* do the other kernels in a multi-kernel list first */
4838   if ( kernel->next != (KernelInfo *) NULL)
4839     ZeroKernelNans(kernel->next);
4840
4841   for (i=0; i < (kernel->width*kernel->height); i++)
4842     if ( IsNan(kernel->values[i]) )
4843       kernel->values[i] = 0.0;
4844
4845   return;
4846 }