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