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