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

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