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