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

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