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