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

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