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 %                              September 2009                                 %
  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 kernals, 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/memory_.h"
  65 #include "magick/monitor-private.h"
  66 #include "magick/morphology.h"
  67 #include "magick/option.h"
  68 #include "magick/pixel-private.h"
  69 #include "magick/prepress.h"
  70 #include "magick/quantize.h"
  71 #include "magick/registry.h"
  72 #include "magick/semaphore.h"
  73 #include "magick/splay-tree.h"
  74 #include "magick/statistic.h"
  75 #include "magick/string_.h"
  76 #include "magick/string-private.h"
  77 #include "magick/token.h"
  78
  79 /*
  80  * The following are assignments and tests for special floating point numbers
  81  * of value NaN (not a number), that may be used within a Kernel Definition.
  82  * NaN's are defined as part of the IEEE standard for floating point number
  83  * representation.
  84  *
  85  * These are used a Kernel value of NaN means that that kernal position
  86  * is not part of the normal convolution or morphology process, and thus
  87  * allowing the use of 'shaped' kernels.
  88  *
  89  * Special Properities
  90  *   Two NaN's are never equal, even if they are from the same variable
  91  * That is the IsNaN() macro is only true if the value is NaN.
  92  */
  93 #define IsNan(a)   ((a)!=(a))
  94
  95 \f
  96 /*
  97 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  98 %                                                                             %
  99 %                                                                             %
 100 %                                                                             %
 101 %     A c q u i r e K e r n e l F r o m S t r i n g                           %
 102 %                                                                             %
 103 %                                                                             %
 104 %                                                                             %
 105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 106 %
 107 %  AcquireKernelFromString() takes the given string (generally supplied by the
 108 %  user) and converts it into a Morphology/Convolution Kernel.  This allows
 109 %  users to specify a kernel from a number of pre-defined kernels, or to fully
 110 %  specify their own kernel for a specific Convolution or Morphology
 111 %  Operation.
 112 %
 113 %  The kernel so generated can be any rectangular array of floating point
 114 %  values (doubles) with the 'control point' or 'pixel being affected'
 115 %  anywhere within that array of values.
 116 %
 117 %  ASIDE: Previously IM was restricted to a square of odd size using the exact
 118 %  center.
 119 %
 120 %  The floating point values in the kernel can also include a special value
 121 %  known as 'NaN' or 'not a number' to indicate that this value is not part
 122 %  of the kernel array. This allows you to specify a non-rectangular shaped
 123 %  kernel, for use in Morphological operators, without the need for some type
 124 %  of kernal mask.
 125 %
 126 %  The returned kernel should be freed using the DestroyKernel() when you are
 127 %  finished with it.
 128 %
 129 %  Input kernel defintion strings can consist of any of three types.
 130 %
 131 %    "num, num, num, num, ..."
 132 %         list of floating point numbers defining an 'old style' odd sized
 133 %         square kernel.  At least 9 values should be provided for a 3x3
 134 %         square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
 135 %         Values can be space or comma separated.
 136 %
 137 %    "WxH[+X+Y]:num, num, num ..."
 138 %         a kernal of size W by H, with W*H floating point numbers following.
 139 %         the 'center' can be optionally be defined at +X+Y (such that +0+0
 140 %         is top left corner). If not defined a pixel closest to the center
 141 %         of the array is automatically defined.
 142 %
 143 %    "name:args"
 144 %         Select from one of the built in kernels. See AcquireKernelBuiltIn()
 145 %
 146 %  Note that 'name' kernels will start with an alphabetic character
 147 %  while the new kernel specification has a ':' character in its
 148 %  specification.
 149 %
 150 % TODO: bias and auto-scale handling of the kernel
 151 %     The given kernel is assumed to have been pre-scaled appropriatally, usally
 152 %     by the kernel generator.
 153 %
 154 %  The format of the AcquireKernal method is:
 155 %
 156 %      MagickKernel *AcquireKernelFromString(const char *kernel_string)
 157 %
 158 %  A description of each parameter follows:
 159 %
 160 %    o kernel_string: the Morphology/Convolution kernel wanted.
 161 %
 162 */
 163
 164 MagickExport MagickKernel *AcquireKernelFromString(const char *kernel_string)
 165 {
 166   MagickKernel
 167     *kernel;
 168
 169   char
 170     token[MaxTextExtent];
 171
 172   register unsigned long
 173     i;
 174
 175   const char
 176     *p;
 177
 178   MagickStatusType
 179     flags;
 180
 181   GeometryInfo
 182     args;
 183
 184   assert(kernel_string != (const char *) NULL);
 185   SetGeometryInfo(&args);
 186
 187   /* does it start with an alpha - Return a builtin kernel */
 188   GetMagickToken(kernel_string,&p,token);
 189   if ( isalpha((int)token[0]) )
 190   {
 191     long
 192       type;
 193
 194     type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
 195     if ( type < 0 || type == UserDefinedKernel )
 196       return((MagickKernel *)NULL);
 197
 198     while (((isspace((int) ((unsigned char) *p)) != 0) ||
 199            (*p == ',') || (*p == ':' )) && (*p != '\0'))
 200       p++;
 201     flags = ParseGeometry(p, &args);
 202
 203     /* special handling of missing values in input string */
 204     if ( type == RectangleKernel ) {
 205       if ( (flags & WidthValue) == 0 ) /* if no width then */
 206         args.rho = args.sigma;         /* then  width = height */
 207       if ( args.rho < 1.0 )            /* if width too small */
 208          args.rho = 3;                 /* then  width = 3 */
 209       if ( args.sigma < 1.0 )          /* if height too small */
 210         args.sigma = args.rho;         /* then  height = width */
 211       if ( (flags & XValue) == 0 )     /* center offset if not defined */
 212         args.xi = (double)(((long)args.rho-1)/2);
 213       if ( (flags & YValue) == 0 )
 214         args.psi = (double)(((long)args.sigma-1)/2);
 215     }
 216
 217     return(AcquireKernelBuiltIn((MagickKernelType)type, &args));
 218   }
 219
 220   kernel=(MagickKernel *) AcquireMagickMemory(sizeof(*kernel));
 221   if (kernel == (MagickKernel *)NULL)
 222     return(kernel);
 223   (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
 224   kernel->type = UserDefinedKernel;
 225
 226   /* Has a ':' in argument - New user kernel specification */
 227   p = strchr(kernel_string, ':');
 228   if ( p != (char *) NULL)
 229     {
 230 #if 1
 231       /* ParseGeometry() needs the geometry separated! -- Arrgghh */
 232       memcpy(token, kernel_string, p-kernel_string);
 233       token[p-kernel_string] = '\0';
 234       flags = ParseGeometry(token, &args);
 235 #else
 236       flags = ParseGeometry(kernel_string, &args);
 237 #endif
 238
 239       /* Size Handling and Checks */
 240       if ( (flags & WidthValue) == 0 ) /* if no width then */
 241         args.rho = args.sigma;         /* then  width = height */
 242       if ( args.rho < 1.0 )            /* if width too small */
 243          args.rho = 1.0;               /* then  width = 1 */
 244       if ( args.sigma < 1.0 )          /* if height too small */
 245         args.sigma = args.rho;         /* then  height = width */
 246       kernel->width = (unsigned long)args.rho;
 247       kernel->height = (unsigned long)args.sigma;
 248
 249       /* Offset Handling and Checks */
 250       if ( args.xi  < 0.0 || args.psi < 0.0 )
 251         return(DestroyKernel(kernel));
 252       kernel->offset_x = ((flags & XValue)!=0) ? (unsigned long)args.xi
 253                                                : (kernel->width-1)/2;
 254       kernel->offset_y = ((flags & YValue)!=0) ? (unsigned long)args.psi
 255                                                : (kernel->height-1)/2;
 256       if ( kernel->offset_x >= kernel->width ||
 257            kernel->offset_y >= kernel->height )
 258         return(DestroyKernel(kernel));
 259
 260       p++; /* advance beyond the ':' */
 261     }
 262   else
 263     { /* ELSE - Old old kernel specification, forming odd-square kernel */
 264       /* count up number of values given */
 265       p=(const char *) kernel_string;
 266       for (i=0; *p != '\0'; i++)
 267       {
 268         GetMagickToken(p,&p,token);
 269         if (*token == ',')
 270           GetMagickToken(p,&p,token);
 271       }
 272       /* set the size of the kernel - old sized square */
 273       kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
 274       kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 275       p=(const char *) kernel_string;
 276     }
 277
 278   /* Read in the kernel values from rest of input string argument */
 279   kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 280                         kernel->height*sizeof(double));
 281   if (kernel->values == (double *) NULL)
 282     return(DestroyKernel(kernel));
 283
 284   kernel->range_neg = kernel->range_pos = 0.0;
 285   for (i=0; (i < kernel->width*kernel->height) && (*p != '\0'); i++)
 286   {
 287     GetMagickToken(p,&p,token);
 288     if (*token == ',')
 289       GetMagickToken(p,&p,token);
 290     (( kernel->values[i] = StringToDouble(token) ) < 0)
 291         ?  ( kernel->range_neg += kernel->values[i] )
 292         :  ( kernel->range_pos += kernel->values[i] );
 293   }
 294   for ( ; i < kernel->width*kernel->height; i++)
 295     kernel->values[i]=0.0;
 296
 297   return(kernel);
 298 }
 299 \f
 300 /*
 301 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 302 %                                                                             %
 303 %                                                                             %
 304 %                                                                             %
 305 %     A c q u i r e K e r n e l B u i l t I n                                 %
 306 %                                                                             %
 307 %                                                                             %
 308 %                                                                             %
 309 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 310 %
 311 %  AcquireKernelBuiltIn() returned one of the 'named' built-in types of
 312 %  kernels used for special purposes such as gaussian blurring, skeleton
 313 %  pruning, and edge distance determination.
 314 %
 315 %  They take a KernelType, and a set of geometry style arguments, which were
 316 %  typically decoded from a user supplied string, or from a more complex
 317 %  Morphology Method that was requested.
 318 %
 319 %  The format of the AcquireKernalBuiltIn method is:
 320 %
 321 %      MagickKernel *AcquireKernelBuiltIn(const MagickKernelType type,
 322 %           const GeometryInfo args)
 323 %
 324 %  A description of each parameter follows:
 325 %
 326 %    o type: the pre-defined type of kernel wanted
 327 %
 328 %    o args: arguments defining or modifying the kernel
 329 %
 330 %  Convolution Kernels
 331 %
 332 %    Gaussian  "[{radius}]x{sigma}"
 333 %       Generate a two-dimentional gaussian kernel, as used by -gaussian
 334 %       A sigma is required, (with the 'x'), due to historical reasons.
 335 %
 336 %       NOTE: that the 'radius' is optional, but if provided can limit (clip)
 337 %       the final size of the resulting kernel to a square 2*radius+1 in size.
 338 %       The radius should be at least 2 times that of the sigma value, or
 339 %       sever clipping and aliasing may result.  If not given or set to 0 the
 340 %       radius will be determined so as to produce the best minimal error
 341 %       result, which is usally much larger than is normally needed.
 342 %
 343 %    Blur  "[{radius}]x{sigma}[+angle]"
 344 %       As per Gaussian, but generates a 1 dimensional or linear gaussian
 345 %       blur, at the angle given (current restricted to orthogonal angles).
 346 %       If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
 347 %
 348 %       NOTE that two such blurs perpendicular to each other is equivelent to
 349 %       -blur and the previous gaussian, but is often 10 or more times faster.
 350 %
 351 %    Comet  "[{width}]x{sigma}[+angle]"
 352 %       Blur in one direction only, mush like how a bright object leaves
 353 %       a comet like trail.  The Kernel is actually half a gaussian curve,
 354 %       Adding two such blurs in oppiste directions produces a Linear Blur.
 355 %
 356 %       NOTE: that the first argument is the width of the kernel and not the
 357 %       radius of the kernel.
 358 %
 359 %    # Still to be implemented...
 360 %    #
 361 %    # Laplacian  "{radius}x{sigma}"
 362 %    #    Laplacian (a mexican hat like) Function
 363 %    #
 364 %    # LOG  "{radius},{sigma1},{sigma2}
 365 %    #    Laplacian of Gaussian
 366 %    #
 367 %    # DOG  "{radius},{sigma1},{sigma2}
 368 %    #    Difference of Gaussians
 369 %
 370 %  Boolean Kernels
 371 %
 372 %    Rectangle "{geometry}"
 373 %       Simply generate a rectangle of 1's with the size given. You can also
 374 %       specify the location of the 'control point', otherwise the closest
 375 %       pixel to the center of the rectangle is selected.
 376 %
 377 %       Properly centered and odd sized rectangles work the best.
 378 %
 379 %    Diamond  "[{radius}]"
 380 %       Generate a diamond shaped kernal with given radius to the points.
 381 %       Kernel size will again be radius*2+1 square and defaults to radius 1,
 382 %       generating a 3x3 kernel that is slightly larger than a square.
 383 %
 384 %    Square  "[{radius}]"
 385 %       Generate a square shaped kernel of size radius*2+1, and defaulting
 386 %       to a 3x3 (radius 1).
 387 %
 388 %       Note that using a larger radius for the "Square" or the "Diamond"
 389 %       is also equivelent to iterating the basic morphological method
 390 %       that many times. However However iterating with the smaller radius 1
 391 %       default is actually faster than using a larger kernel radius.
 392 %
 393 %    Disk   "[{radius}]
 394 %       Generate a binary disk of the radius given, radius may be a float.
 395 %       Kernel size will be ceil(radius)*2+1 square.
 396 %       NOTE: Here are some disk shapes of specific interest
 397 %          "disk:1"    => "diamond" or "cross:1"
 398 %          "disk:1.5"  => "square"
 399 %          "disk:2"    => "diamond:2"
 400 %          "disk:2.5"  => default - radius 2 disk shape
 401 %          "disk:2.9"  => "square:2"
 402 %          "disk:3.5"  => octagonal/disk shape of radius 3
 403 %          "disk:4.2"  => roughly octagonal shape of radius 4
 404 %          "disk:4.3"  => disk shape of radius 4
 405 %       After this all the kernel shape becomes more and more circular.
 406 %
 407 %       Because a "disk" is more circular when using a larger radius, using a
 408 %       larger radius is preferred over iterating the morphological operation.
 409 %
 410 %    Plus  "[{radius}]"
 411 %       Generate a kernel in the shape of a 'plus' sign. The length of each
 412 %       arm is also the radius, which defaults to 2.
 413 %
 414 %       This kernel is not a good general morphological kernel, but is used
 415 %       more for highlighting and marking any single pixels in an image using,
 416 %       a "Dilate" or "Erode" method as appropriate.
 417 #
 418 %       NOTE: "plus:1" is equivelent to a "Diamond" kernel.
 419 %
 420 %       Note that unlike other kernels iterating a plus does not produce the
 421 %       same result as using a larger radius for the cross.
 422 %
 423 %  Distance Measuring Kernels
 424 %
 425 %    Chebyshev "[{radius}][x{scale}]"   largest x or y distance (default r=1)
 426 %    Manhatten "[{radius}][x{scale}]"   square grid distance    (default r=1)
 427 %    Knight    "[{radius}][x{scale}]"   octagonal distance      (default r=1)
 428 %    Euclidean "[{radius}][x{scale}]"   direct distance         (default r=4)
 429 %
 430 %       Different types of distance measuring methods, which are used with the
 431 %       a 'Distance' morphology method for generating a gradient based on
 432 %       distance from an edge of a binary shape, though there is a technique
 433 %       for handling a anti-aliased shape.
 434 %
 435 %       The first 3 are simplifications that alow the use of a small kernel
 436 %       which is iterated. The lest is more accurate but requires a larger
 437 %       kernel to produce a accurate distance measure. The larger the better.
 438 %
 439 %       The actual distance is scaled the size give, which while unnecessary
 440 %       for a "Chebyshev" or "Manhatten" distance, is needed to allow for
 441 %       correct handling of fractional distances in "Knight" and "Euclidean"
 442 %       distance formulas.  If no scale is provided it is set to a value of
 443 %       100, allowing for a maximum distance measurement of 655 pixels from
 444 %       any edge, using a Q16 version of IM.
 445 %
 446 %       See the 'Distance' Morphological Method, for information of how it
 447 %       is applied.
 448 %
 449 */
 450
 451 static void KernelRotate(MagickKernel *kernel, double angle)
 452 {
 453   /* Rotate a kernel appropriately for the angle given
 454   **
 455   ** Currently assumes the kernel (rightly) horizontally is symetrical
 456   **
 457   ** TODO: expand beyond simple 90 degree rotates, flips and flops
 458   */
 459
 460   /* Modulus the angle */
 461   angle = fmod(angle, 360.0);
 462   if ( angle < 0 )
 463     angle += 360.0;
 464
 465   if ( 315.0 < angle || angle <= 45.0 )
 466     return;   /* no change! - At least at this time */
 467
 468   switch (kernel->type) {
 469     /* These kernels are cylindrical kernel, rotating is useless */
 470     case GaussianKernel:
 471     case LaplacianKernel:
 472     case LOGKernel:
 473     case DOGKernel:
 474     case DiskKernel:
 475     case ChebyshevKernel:
 476     case ManhattenKernel:
 477     case KnightKernel:
 478     case EuclideanKernel:
 479       return;
 480
 481     /* These may be rotatable at non-90 angles in the future */
 482     /* but simply rotating them 90 degrees is useless */
 483     case SquareKernel:
 484     case DiamondKernel:
 485     case PlusKernel:
 486       return;
 487
 488     /* These only allows a +/-90 degree rotation (transpose) */
 489     case BlurKernel:
 490     case RectangleKernel:
 491       if ( 135.0 < angle && angle <= 225.0 )
 492         return;
 493       if ( 225.0 < angle && angle <= 315.0 )
 494         angle -= 180;
 495       break;
 496
 497     /* these are freely rotatable in 90 degree units */
 498     case CometKernel:
 499     case UndefinedKernel:
 500     case UserDefinedKernel:
 501       break;
 502   }
 503
 504 fprintf(stderr, "angle2 = %lf\n", angle);
 505
 506   if ( 135.0 < angle && angle <= 315.0 )
 507     {
 508       /* Do a flop, this assumes kernel is horizontally symetrical. */
 509       /* Each kernel data row need to be reversed! */
 510       unsigned long
 511         y;
 512       register unsigned long
 513         x,r;
 514       register double
 515         *k,t;
 516       for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width) {
 517         for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
 518           t=k[x],  k[x]=k[r],  k[r]=t;
 519       }
 520       kernel->offset_x = kernel->width - kernel->offset_x - 1;
 521       angle = fmod(angle+180.0, 360.0);
 522     }
 523   if ( 45.0 < angle && angle <= 135.0 )
 524     {
 525       /* Do a transpose, this assumes the kernel is orthoginally symetrical */
 526       /* The data is the same, just the size and offsets needs to be swapped. */
 527       unsigned long
 528         t;
 529       t = kernel->width;
 530       kernel->width = kernel->height;
 531       kernel->height = t;
 532       t = kernel->offset_x;
 533       kernel->offset_x = kernel->offset_y;
 534       kernel->offset_y = t;
 535       angle = fmod(450.0 - angle, 360.0);
 536     }
 537   /* at this point angle should be between +45 and -45 (315) degrees */
 538   return;
 539 }
 540
 541 MagickExport MagickKernel *AcquireKernelBuiltIn(const MagickKernelType type,
 542    const GeometryInfo *args)
 543 {
 544   MagickKernel
 545     *kernel;
 546
 547   register unsigned long
 548     i;
 549
 550   register long
 551     u,
 552     v;
 553
 554   double
 555     nan = sqrt((double)-1.0);  /* Special Value : Not A Number */
 556
 557   kernel=(MagickKernel *) AcquireMagickMemory(sizeof(*kernel));
 558   if (kernel == (MagickKernel *) NULL)
 559     return(kernel);
 560   (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
 561   kernel->range_neg = kernel->range_pos = 0.0;
 562   kernel->type = type;
 563
 564   switch(type) {
 565     /* Convolution Kernels */
 566     case GaussianKernel:
 567       { double
 568           sigma = fabs(args->sigma);
 569
 570         sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
 571
 572         kernel->width = kernel->height =
 573                             GetOptimalKernelWidth2D(args->rho,sigma);
 574         kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 575         kernel->range_neg = kernel->range_pos = 0.0;
 576         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 577                               kernel->height*sizeof(double));
 578         if (kernel->values == (double *) NULL)
 579           return(DestroyKernel(kernel));
 580
 581         sigma = 2.0*sigma*sigma; /* simplify the expression */
 582         for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
 583           for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
 584             kernel->range_pos += (
 585               kernel->values[i] =
 586                  exp(-((double)(u*u+v*v))/sigma)
 587                        /*  / (MagickPI*sigma)  */ );
 588
 589         /* Normalize the Kernel - see notes in BlurKernel, below */
 590         u=kernel->width*kernel->height;
 591         for (i=0; i < (unsigned long)u; i++)
 592           kernel->values[i] /= kernel->range_pos;
 593         kernel->range_pos=1.0;
 594
 595         break;
 596       }
 597     case BlurKernel:
 598       { double
 599           sigma = fabs(args->sigma);
 600
 601         sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
 602
 603         kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
 604         kernel->offset_x = (kernel->width-1)/2;
 605         kernel->height = 1;
 606         kernel->offset_y = 0;
 607         kernel->range_neg = kernel->range_pos = 0.0;
 608         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 609                               kernel->height*sizeof(double));
 610         if (kernel->values == (double *) NULL)
 611           return(DestroyKernel(kernel));
 612
 613 #if 1
 614 #define KernelRank 3
 615         /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
 616         ** It generates a gaussian 3 times the width, and compresses it into
 617         ** the expected range.  This produces a closer normalization of the
 618         ** resulting kernel, especially for very low sigma values.
 619         ** As such while wierd it is prefered.
 620         **
 621         ** I am told this method originally came from Photoshop.
 622         */
 623         sigma *= KernelRank;                /* simplify expanded curve */
 624         v = (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
 625         (void) ResetMagickMemory(kernel->values,0, (size_t)
 626                        kernel->width*sizeof(double));
 627         for ( u=-v; u <= v; u++) {
 628           kernel->values[(u+v)/KernelRank] +=
 629                 exp(-((double)(u*u))/(2.0*sigma*sigma))
 630                        /*   / (MagickSQ2PI*sigma/KernelRank)  */ ;
 631         }
 632         for (i=0; i < kernel->width; i++)
 633           kernel->range_pos += kernel->values[i];
 634 #else
 635         for ( i=0, u=-kernel->offset_x; i < kernel->width; i++, u++)
 636           kernel->range_pos += (
 637               kernel->values[i] =
 638                 exp(-((double)(u*u))/(2.0*sigma*sigma))
 639                        /*  / (MagickSQ2PI*sigma)  */ );
 640 #endif
 641         /* Note that both the above methods do not generate a normalized
 642         ** kernel, though it gets close. The kernel may be 'clipped' by a user
 643         ** defined radius, producing a smaller (darker) kernel.  Also for very
 644         ** small sigma's (> 0.1) the central value becomes larger than one,
 645         ** and thus producing a bright kernel.
 646         */
 647 #if 1
 648         /* Normalize the 1D Gaussian Kernel
 649         **
 650         ** Because of this the divisor in the above kernel generator is
 651         ** not needed, and is taken care of here.
 652         */
 653         for (i=0; i < kernel->width; i++)
 654           kernel->values[i] /= kernel->range_pos;
 655         kernel->range_pos=1.0;
 656 #endif
 657         /* rotate the kernel by given angle */
 658         KernelRotate(kernel, args->xi);
 659         break;
 660       }
 661     case CometKernel:
 662       { double
 663           sigma = fabs(args->sigma);
 664
 665         sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
 666
 667         if ( args->rho < 1.0 )
 668           kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
 669         else
 670           kernel->width = (unsigned long)args->rho;
 671         kernel->offset_x = kernel->offset_y = 0;
 672         kernel->height = 1;
 673         kernel->range_neg = kernel->range_pos = 0.0;
 674         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 675                               kernel->height*sizeof(double));
 676         if (kernel->values == (double *) NULL)
 677           return(DestroyKernel(kernel));
 678
 679         /* A comet blur is half a gaussian curve, so that the object is
 680         ** blurred in one direction only.  This may not be quite the right
 681         ** curve so may change in the future. The function must be normalised.
 682         */
 683 #if 1
 684 #define KernelRank 3
 685         sigma *= KernelRank;                /* simplify expanded curve */
 686         v = kernel->width*KernelRank; /* start/end points to fit range */
 687         (void) ResetMagickMemory(kernel->values,0, (size_t)
 688                        kernel->width*sizeof(double));
 689         for ( u=0; u < v; u++) {
 690           kernel->values[u/KernelRank] +=
 691                exp(-((double)(u*u))/(2.0*sigma*sigma))
 692                        /*   / (MagickSQ2PI*sigma/KernelRank)  */ ;
 693         }
 694         for (i=0; i < kernel->width; i++)
 695           kernel->range_pos += kernel->values[i];
 696 #else
 697         for ( i=0; i < kernel->width; i++)
 698           kernel->range_pos += (
 699             kernel->values[i] =
 700                exp(-((double)(i*i))/(2.0*sigma*sigma))
 701                        /*  / (MagickSQ2PI*sigma)  */ );
 702 #endif
 703         /* Normalize the Kernel - see notes in BlurKernel */
 704         for (i=0; i < kernel->width; i++)
 705           kernel->values[i] /= kernel->range_pos;
 706         kernel->range_pos=1.0;
 707
 708         /* rotate the kernel by given angle */
 709         KernelRotate(kernel, args->xi);
 710         break;
 711       }
 712     /* Boolean Kernels */
 713     case RectangleKernel:
 714     case SquareKernel:
 715       {
 716         if ( type == SquareKernel )
 717           {
 718             if (args->rho < 1.0)
 719               kernel->width = kernel->height = 3; /* radius 1 */
 720             else
 721               kernel->width = kernel->height = 2*(long)args->rho+1;
 722             kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 723           }
 724         else {
 725             if ( args->rho < 1.0 || args->sigma < 1.0 )
 726               return(DestroyKernel(kernel));
 727             kernel->width = (unsigned long)args->rho;
 728             kernel->height = (unsigned long)args->sigma;
 729             if ( args->xi  < 0.0 || args->xi  > (double)kernel->width ||
 730                  args->psi < 0.0 || args->psi > (double)kernel->height )
 731               return(DestroyKernel(kernel));
 732             kernel->offset_x = (unsigned long)args->xi;
 733             kernel->offset_y = (unsigned long)args->psi;
 734           }
 735         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 736                               kernel->height*sizeof(double));
 737         if (kernel->values == (double *) NULL)
 738           return(DestroyKernel(kernel));
 739
 740         u=kernel->width*kernel->height;
 741         for ( i=0; i < (unsigned long)u; i++)
 742             kernel->values[i] = 1.0;
 743         break;
 744       }
 745     case DiamondKernel:
 746       {
 747         if (args->rho < 1.0)
 748           kernel->width = kernel->height = 3; /* radius 1 */
 749         else
 750           kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
 751         kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 752
 753         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 754                               kernel->height*sizeof(double));
 755         if (kernel->values == (double *) NULL)
 756           return(DestroyKernel(kernel));
 757
 758         for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
 759           for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
 760             if ((labs(u)+labs(v)) <= (long)kernel->offset_x)
 761               kernel->range_pos += kernel->values[i] = 1.0;
 762             else
 763               kernel->values[i] = nan;
 764         break;
 765       }
 766     case DiskKernel:
 767       {
 768         long
 769           limit;
 770
 771         limit = (long)(args->rho*args->rho);
 772         if (args->rho < 1.0)        /* default: ~2.5 radius disk */
 773           kernel->width = kernel->height = 5L, limit = 5L;
 774         else
 775            kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
 776         kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 777
 778         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 779                               kernel->height*sizeof(double));
 780         if (kernel->values == (double *) NULL)
 781           return(DestroyKernel(kernel));
 782
 783         for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
 784           for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
 785             if ((u*u+v*v) <= limit)
 786               kernel->range_pos += kernel->values[i] = 1.0;
 787             else
 788               kernel->values[i] = nan;
 789         break;
 790       }
 791     case PlusKernel:
 792       {
 793         if (args->rho < 1.0)
 794           kernel->width = kernel->height = 5;  /* radius 2 */
 795         else
 796            kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
 797         kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 798
 799         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 800                               kernel->height*sizeof(double));
 801         if (kernel->values == (double *) NULL)
 802           return(DestroyKernel(kernel));
 803
 804         for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
 805           for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
 806             kernel->values[i] = (u == 0 || v == 0) ? 1.0 : nan;
 807         kernel->range_pos = kernel->width*2.0 - 1.0;
 808         break;
 809       }
 810     /* Distance Measuring Kernels */
 811     case ChebyshevKernel:
 812       {
 813         double
 814           scale;
 815
 816         if (args->rho < 1.0)
 817           kernel->width = kernel->height = 3;
 818         else
 819           kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
 820         kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 821
 822         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 823                               kernel->height*sizeof(double));
 824         if (kernel->values == (double *) NULL)
 825           return(DestroyKernel(kernel));
 826
 827         scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
 828         for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
 829           for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
 830             kernel->range_pos += ( kernel->values[i] =
 831                  scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
 832         break;
 833       }
 834     case ManhattenKernel:
 835       {
 836         double
 837           scale;
 838
 839         if (args->rho < 1.0)
 840           kernel->width = kernel->height = 3;
 841         else
 842            kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
 843         kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 844
 845         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 846                               kernel->height*sizeof(double));
 847         if (kernel->values == (double *) NULL)
 848           return(DestroyKernel(kernel));
 849
 850         scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
 851         for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
 852           for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
 853             kernel->range_pos += ( kernel->values[i] =
 854                  scale*(labs(u)+labs(v)) );
 855         break;
 856       }
 857     case KnightKernel:
 858       {
 859         double
 860           scale;
 861
 862         if (args->rho < 1.0)
 863           kernel->width = kernel->height = 3;
 864         else
 865            kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
 866         kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 867
 868         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 869                               kernel->height*sizeof(double));
 870         if (kernel->values == (double *) NULL)
 871           return(DestroyKernel(kernel));
 872
 873         scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
 874         for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
 875           for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
 876             kernel->range_pos += ( kernel->values[i] =
 877                scale*((labs(u)<labs(v)) ? (MagickSQ2-1.0)*labs(u)+labs(v)
 878                                         : (MagickSQ2-1.0)*labs(v)+labs(u) ) );
 879         break;
 880       }
 881     case EuclideanKernel:
 882       {
 883         double
 884           scale;
 885
 886         if (args->rho < 1.0)
 887           kernel->width = kernel->height = 9;
 888         else
 889            kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
 890         kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
 891
 892         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
 893                               kernel->height*sizeof(double));
 894         if (kernel->values == (double *) NULL)
 895           return(DestroyKernel(kernel));
 896
 897         scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
 898         for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
 899           for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
 900             kernel->range_pos += ( kernel->values[i] =
 901                  scale*sqrt((double)(u*u+v*v)) );
 902         break;
 903       }
 904     /* Undefined Kernels */
 905     case LaplacianKernel:
 906     case LOGKernel:
 907     case DOGKernel:
 908       assert("Kernel Type has not been defined yet");
 909       /* FALL THRU */
 910     default:
 911       /* Generate a No-Op minimal kernel - 1x1 pixel */
 912       kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
 913       if (kernel->values == (double *) NULL)
 914         return(DestroyKernel(kernel));
 915       kernel->range_pos = kernel->values[0] = 1.0;
 916       kernel->width = kernel->height = 1;
 917       kernel->offset_x = kernel->offset_x = 0;
 918       kernel->type = UndefinedKernel;
 919       break;
 920   }
 921
 922   return(kernel);
 923 }
 924
 925 /*
 926 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 927 %                                                                             %
 928 %                                                                             %
 929 %                                                                             %
 930 %     D e s t r o y K e r n e l                                               %
 931 %                                                                             %
 932 %                                                                             %
 933 %                                                                             %
 934 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 935 %
 936 %  DestroyKernel() frees the memory used by a Convolution/Morphology kernel.
 937 %
 938 %  The format of the DestroyKernel method is:
 939 %
 940 %      MagickKernel *DestroyKernel(MagickKernel *kernel)
 941 %
 942 %  A description of each parameter follows:
 943 %
 944 %    o kernel: the Morphology/Convolution kernel to be destroyed
 945 %
 946 */
 947
 948 MagickExport MagickKernel *DestroyKernel(MagickKernel *kernel)
 949 {
 950   assert(kernel != (MagickKernel *) NULL);
 951   kernel->values=(double *)RelinquishMagickMemory(kernel->values);
 952   kernel=(MagickKernel *) RelinquishMagickMemory(kernel);
 953   return(kernel);
 954 }
 955
 956 \f
 957 /*
 958 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 959 %                                                                             %
 960 %                                                                             %
 961 %                                                                             %
 962 %     M o r p h o l o g y I m a g e                                           %
 963 %                                                                             %
 964 %                                                                             %
 965 %                                                                             %
 966 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 967 %
 968 %  MorphologyImage() applies a user supplied kernel to the image according to
 969 %  the given mophology method.
 970 %
 971 %  The given kernel is assumed to have been pre-scaled appropriatally, usally
 972 %  by the kernel generator.
 973 %
 974 %  The format of the MorphologyImage method is:
 975 %
 976 %      Image *MorphologyImage(const Image *image, const MorphologyMethod
 977 %        method, const long iterations, const ChannelType channel,
 978 %        const MagickKernel *kernel, ExceptionInfo *exception)
 979 %
 980 %  A description of each parameter follows:
 981 %
 982 %    o image: the image.
 983 %
 984 %    o method: the morphology method to be applied.
 985 %
 986 %    o iterations: apply the operation this many times (or no change).
 987 %                  A value of -1 means loop until no change found.
 988 %                  How this is applied may depend on the morphology method.
 989 %                  Typically this is a value of 1.
 990 %
 991 %    o channel: the channel type.
 992 %
 993 %    o kernel: An array of double representing the morphology kernel.
 994 %              This is assumed to have been pre-scaled (normalized).
 995 %
 996 %    o exception: return any errors or warnings in this structure.
 997 %
 998 %
 999 % TODO: bias and auto-scale handling of the kernel for convolution
1000 %     The given kernel is assumed to have been pre-scaled appropriatally, usally
1001 %     by the kernel generator.
1002 %
1003 */
1004
1005 static inline double MagickMin(const MagickRealType x,const MagickRealType y)
1006 {
1007   return( x < y ? x : y);
1008 }
1009 static inline double MagickMax(const MagickRealType x,const MagickRealType y)
1010 {
1011   return( x > y ? x : y);
1012 }
1013 #define Minimize(assign,value) assign=MagickMin(assign,value)
1014 #define Maximize(assign,value) assign=MagickMax(assign,value)
1015
1016 /* incr change if the value being assigned changed */
1017 #define Assign(channel,value) \
1018   { q->channel = RoundToQuantum(value); \
1019     if ( p[r].channel != q->channel ) changed++; \
1020   }
1021 #define AssignIndex(value) \
1022   { q_indexes[x] = RoundToQuantum(value); \
1023     if ( p_indexes[r] != q_indexes[x] ) changed++; \
1024   }
1025
1026 /* Internal function
1027  * Apply the Morphology method with the given Kernel
1028  * And return the number of values changed.
1029  */
1030 static unsigned long MorphologyApply(const Image *image, Image
1031      *result_image, const MorphologyMethod method, const ChannelType channel,
1032      const MagickKernel *kernel, ExceptionInfo *exception)
1033 {
1034   #define MorphologyTag  "Morphology/Image"
1035
1036   long
1037     progress,
1038     y;
1039
1040   unsigned long
1041     changed;
1042
1043   MagickBooleanType
1044     status;
1045
1046   MagickPixelPacket
1047     bias;
1048
1049   CacheView
1050     *p_view,
1051     *q_view;
1052
1053   /*
1054     Apply Morphology to Image.
1055   */
1056   status=MagickTrue;
1057   changed=0;
1058   progress=0;
1059
1060   GetMagickPixelPacket(image,&bias);
1061   SetMagickPixelPacketBias(image,&bias);
1062
1063   p_view=AcquireCacheView(image);
1064   q_view=AcquireCacheView(result_image);
1065 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1066   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1067 #endif
1068   for (y=0; y < (long) image->rows; y++)
1069   {
1070     MagickBooleanType
1071       sync;
1072
1073     register const PixelPacket
1074       *restrict p;
1075
1076     register const IndexPacket
1077       *restrict p_indexes;
1078
1079     register PixelPacket
1080       *restrict q;
1081
1082     register IndexPacket
1083       *restrict q_indexes;
1084
1085     register long
1086       x;
1087
1088     long
1089       r;
1090
1091     if (status == MagickFalse)
1092       continue;
1093     p=GetCacheViewVirtualPixels(p_view, -kernel->offset_x,  y-kernel->offset_y,
1094          image->columns+kernel->width, kernel->height, exception);
1095     q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1096          exception);
1097     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1098       {
1099         status=MagickFalse;
1100         continue;
1101       }
1102     p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1103     q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1104     r = (image->columns+kernel->width)*kernel->offset_y+kernel->offset_x;
1105     for (x=0; x < (long) image->columns; x++)
1106     {
1107       long
1108         v;
1109
1110       register long
1111         u;
1112
1113       register const double
1114         *restrict k;
1115
1116       register const PixelPacket
1117         *restrict k_pixels;
1118
1119       register const IndexPacket
1120         *restrict k_indexes;
1121
1122       MagickPixelPacket
1123         result;
1124
1125       /* Copy input to ouput image - removes need for 'cloning' new images */
1126       *q = p[r];
1127       if (image->colorspace == CMYKColorspace)
1128         q_indexes[x] = p_indexes[r];
1129
1130       result.index=0;
1131       switch (method) {
1132         case ConvolveMorphology:
1133           result=bias;
1134           break;  /* default result is the convolution bias */
1135         case DialateIntensityMorphology:
1136         case ErodeIntensityMorphology:
1137           /* result is the pixel as is */
1138           result.red     = p[r].red;
1139           result.green   = p[r].green;
1140           result.blue    = p[r].blue;
1141           result.opacity = p[r].opacity;
1142           if ( image->colorspace == CMYKColorspace)
1143              result.index   = p_indexes[r];
1144           break;
1145         default:
1146           /* most need to handle transparency as alpha */
1147           result.red     = p[r].red;
1148           result.green   = p[r].green;
1149           result.blue    = p[r].blue;
1150           result.opacity = QuantumRange - p[r].opacity;
1151           if ( image->colorspace == CMYKColorspace)
1152              result.index   = p_indexes[r];
1153           break;
1154       }
1155
1156       switch ( method ) {
1157         case ConvolveMorphology:
1158             /* Weighted Average of pixels */
1159             if (((channel & OpacityChannel) == 0) ||
1160                       (image->matte == MagickFalse))
1161               {
1162                 /* Kernel Weighted Convolution (no transparency) */
1163                 k = kernel->values;
1164                 k_pixels = p;
1165                 k_indexes = p_indexes;
1166                 for (v=0; v < (long) kernel->height; v++) {
1167                   for (u=0; u < (long) kernel->width; u++, k++) {
1168                     if ( IsNan(*k) ) continue;
1169                     result.red     += (*k)*k_pixels[u].red;
1170                     result.green   += (*k)*k_pixels[u].green;
1171                     result.blue    += (*k)*k_pixels[u].blue;
1172                     /* result.opacity += no involvment */
1173                     if ( image->colorspace == CMYKColorspace)
1174                       result.index   += (*k)*k_indexes[u];
1175                   }
1176                   k_pixels += image->columns+kernel->width;
1177                   k_indexes += image->columns+kernel->width;
1178                 }
1179                 if ((channel & RedChannel) != 0)
1180                   Assign(red,result.red);
1181                 if ((channel & GreenChannel) != 0)
1182                   Assign(green,result.green);
1183                 if ((channel & BlueChannel) != 0)
1184                   Assign(blue,result.blue);
1185                 /* no transparency involved */
1186                 if ((channel & IndexChannel) != 0
1187                     && image->colorspace == CMYKColorspace)
1188                   AssignIndex(result.index);
1189               }
1190             else
1191               { /* Kernel & Alpha weighted Convolution */
1192                 MagickRealType
1193                   alpha,  /* alpha value * kernel weighting */
1194                   gamma;  /* weighting divisor */
1195
1196                 gamma=0.0;
1197                 k = kernel->values;
1198                 k_pixels = p;
1199                 k_indexes = p_indexes;
1200                 for (v=0; v < (long) kernel->height; v++) {
1201                   for (u=0; u < (long) kernel->width; u++, k++) {
1202                     if ( IsNan(*k) ) continue;
1203                     alpha=(*k)*(QuantumScale*(QuantumRange-
1204                                           k_pixels[u].opacity));
1205                     gamma += alpha;
1206                     result.red     += alpha*k_pixels[u].red;
1207                     result.green   += alpha*k_pixels[u].green;
1208                     result.blue    += alpha*k_pixels[u].blue;
1209                     result.opacity += (*k)*k_pixels[u].opacity;
1210                     if ( image->colorspace == CMYKColorspace)
1211                       result.index   += alpha*k_indexes[u];
1212                   }
1213                   k_pixels += image->columns+kernel->width;
1214                   k_indexes += image->columns+kernel->width;
1215                 }
1216                 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1217                 if ((channel & RedChannel) != 0)
1218                   Assign(red,gamma*result.red);
1219                 if ((channel & GreenChannel) != 0)
1220                   Assign(green,gamma*result.green);
1221                 if ((channel & BlueChannel) != 0)
1222                   Assign(blue,gamma*result.blue);
1223                 if ((channel & OpacityChannel) != 0
1224                     && image->matte == MagickTrue )
1225                   Assign(opacity,result.opacity);
1226                 if ((channel & IndexChannel) != 0
1227                     && image->colorspace == CMYKColorspace)
1228                   AssignIndex(gamma*result.index);
1229               }
1230             break;
1231
1232         case DialateMorphology:
1233             /* Maximize Value - Kernel should be boolean */
1234             k = kernel->values;
1235             k_pixels = p;
1236             k_indexes = p_indexes;
1237             for (v=0; v < (long) kernel->height; v++) {
1238               for (u=0; u < (long) kernel->width; u++, k++) {
1239                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1240                 Maximize(result.red,     k_pixels[u].red);
1241                 Maximize(result.green,   k_pixels[u].green);
1242                 Maximize(result.blue,    k_pixels[u].blue);
1243                 Maximize(result.opacity, QuantumRange-k_pixels[u].opacity);
1244                 if ( image->colorspace == CMYKColorspace)
1245                   Maximize(result.index,   k_indexes[u]);
1246               }
1247               k_pixels += image->columns+kernel->width;
1248               k_indexes += image->columns+kernel->width;
1249             }
1250             if ((channel & RedChannel) != 0)
1251               Assign(red,result.red);
1252             if ((channel & GreenChannel) != 0)
1253               Assign(green,result.green);
1254             if ((channel & BlueChannel) != 0)
1255               Assign(blue,result.blue);
1256             if ((channel & OpacityChannel) != 0
1257                  && image->matte == MagickTrue )
1258               Assign(opacity,QuantumRange-result.opacity);
1259             if ((channel & IndexChannel) != 0
1260                 && image->colorspace == CMYKColorspace)
1261               AssignIndex(result.index);
1262             break;
1263
1264         case ErodeMorphology:
1265             /* Minimize Value - Kernel should be boolean */
1266             k = kernel->values;
1267             k_pixels = p;
1268             k_indexes = p_indexes;
1269             for (v=0; v < (long) kernel->height; v++) {
1270               for (u=0; u < (long) kernel->width; u++, k++) {
1271                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1272                 Minimize(result.red,     k_pixels[u].red);
1273                 Minimize(result.green,   k_pixels[u].green);
1274                 Minimize(result.blue,    k_pixels[u].blue);
1275                 Minimize(result.opacity, QuantumRange-k_pixels[u].opacity);
1276                 if ( image->colorspace == CMYKColorspace)
1277                   Minimize(result.index,   k_indexes[u]);
1278               }
1279               k_pixels += image->columns+kernel->width;
1280               k_indexes += image->columns+kernel->width;
1281             }
1282             if ((channel & RedChannel) != 0)
1283               Assign(red,result.red);
1284             if ((channel & GreenChannel) != 0)
1285               Assign(green,result.green);
1286             if ((channel & BlueChannel) != 0)
1287               Assign(blue,result.blue);
1288             if ((channel & OpacityChannel) != 0
1289                  && image->matte == MagickTrue )
1290               Assign(opacity,QuantumRange-result.opacity);
1291             if ((channel & IndexChannel) != 0
1292                 && image->colorspace == CMYKColorspace)
1293               AssignIndex(result.index);
1294             break;
1295
1296         case DialateIntensityMorphology:
1297             /* Maximum Intensity Pixel - Kernel should be boolean */
1298             k = kernel->values;
1299             k_pixels = p;
1300             k_indexes = p_indexes;
1301             for (v=0; v < (long) kernel->height; v++) {
1302               for (u=0; u < (long) kernel->width; u++, k++) {
1303                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1304                 if ( PixelIntensity(&p[r]) >
1305                      PixelIntensity(&(k_pixels[u])) ) continue;
1306                 result.red     = k_pixels[u].red;
1307                 result.green   = k_pixels[u].green;
1308                 result.blue    = k_pixels[u].blue;
1309                 result.opacity = k_pixels[u].opacity;
1310                 if ( image->colorspace == CMYKColorspace)
1311                   result.index   = k_indexes[u];
1312               }
1313               k_pixels += image->columns+kernel->width;
1314               k_indexes += image->columns+kernel->width;
1315             }
1316             if ((channel & RedChannel) != 0)
1317               Assign(red,result.red);
1318             if ((channel & GreenChannel) != 0)
1319               Assign(green,result.green);
1320             if ((channel & BlueChannel) != 0)
1321               Assign(blue,result.blue);
1322             if ((channel & OpacityChannel) != 0
1323                  && image->matte == MagickTrue )
1324               Assign(opacity,result.opacity);
1325             if ((channel & IndexChannel) != 0
1326                 && image->colorspace == CMYKColorspace)
1327               AssignIndex(result.index);
1328             break;
1329
1330         case ErodeIntensityMorphology:
1331             /* Minimum Intensity Pixel - Kernel should be boolean */
1332             k = kernel->values;
1333             k_pixels = p;
1334             k_indexes = p_indexes;
1335             for (v=0; v < (long) kernel->height; v++) {
1336               for (u=0; u < (long) kernel->width; u++, k++) {
1337                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1338                 if ( PixelIntensity(&p[r]) <
1339                      PixelIntensity(&(k_pixels[u])) ) continue;
1340                 result.red     = k_pixels[u].red;
1341                 result.green   = k_pixels[u].green;
1342                 result.blue    = k_pixels[u].blue;
1343                 result.opacity = k_pixels[u].opacity;
1344                 if ( image->colorspace == CMYKColorspace)
1345                   result.index   = k_indexes[u];
1346               }
1347               k_pixels += image->columns+kernel->width;
1348               k_indexes += image->columns+kernel->width;
1349             }
1350             if ((channel & RedChannel) != 0)
1351               Assign(red,result.red);
1352             if ((channel & GreenChannel) != 0)
1353               Assign(green,result.green);
1354             if ((channel & BlueChannel) != 0)
1355               Assign(blue,result.blue);
1356             if ((channel & OpacityChannel) != 0
1357                  && image->matte == MagickTrue )
1358               Assign(opacity,result.opacity);
1359             if ((channel & IndexChannel) != 0
1360                 && image->colorspace == CMYKColorspace)
1361               AssignIndex(result.index);
1362             break;
1363
1364         case DistanceMorphology:
1365 #if 0
1366           /* No need to do distance morphology if all values are zero */
1367           /* Unfortunatally I have not been able to get this right! */
1368           if (   ((channel & RedChannel) == 0 && p[r].red == 0)
1369               || ((channel & GreenChannel) == 0 && p[r].green == 0)
1370               || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1371               || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1372               || (( (channel & IndexChannel) == 0
1373                   || image->colorspace != CMYKColorspace
1374                                                ) && p_indexes[x] ==0 )
1375              )
1376             break;
1377 #endif
1378             k = kernel->values;
1379             k_pixels = p;
1380             k_indexes = p_indexes;
1381             for (v=0; v < (long) kernel->height; v++) {
1382               for (u=0; u < (long) kernel->width; u++, k++) {
1383                 if ( IsNan(*k) ) continue;
1384                 Minimize(result.red,     (*k)+k_pixels[u].red);
1385                 Minimize(result.green,   (*k)+k_pixels[u].green);
1386                 Minimize(result.blue,    (*k)+k_pixels[u].blue);
1387                 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1388                 if ( image->colorspace == CMYKColorspace)
1389                   Minimize(result.index,   (*k)+k_indexes[u]);
1390               }
1391               k_pixels += image->columns+kernel->width;
1392               k_indexes += image->columns+kernel->width;
1393             }
1394 #if 1
1395             if ((channel & RedChannel) != 0)
1396               Assign(red,result.red);
1397             if ((channel & GreenChannel) != 0)
1398               Assign(green,result.green);
1399             if ((channel & BlueChannel) != 0)
1400               Assign(blue,result.blue);
1401             if ((channel & OpacityChannel) != 0
1402                  && image->matte == MagickTrue )
1403               Assign(opacity,QuantumRange-result.opacity);
1404             if ((channel & IndexChannel) != 0
1405                 && image->colorspace == CMYKColorspace)
1406               AssignIndex(result.index);
1407 #else
1408             /* By returning the number of 'maximum' values still to process
1409             ** we can get the Distance iteration to finish faster.
1410             ** BUT this may cause an infinite loop on very large shapes,
1411             ** which may have a distance that reachs a maximum gradient.
1412             */
1413             if ((channel & RedChannel) != 0)
1414               { q->red = RoundToQuantum(result.red);
1415                 if ( q->red == QuantumRange ) changed++; /* more to do */
1416               }
1417             if ((channel & GreenChannel) != 0)
1418               { q->green = RoundToQuantum(result.green);
1419                 if ( q->green == QuantumRange ) changed++; /* more to do */
1420               }
1421             if ((channel & BlueChannel) != 0)
1422               { q->blue = RoundToQuantum(result.blue);
1423                 if ( q->blue == QuantumRange ) changed++; /* more to do */
1424               }
1425             if ((channel & OpacityChannel) != 0)
1426               { q->opacity = RoundToQuantum(QuantumRange-result.opacity);
1427                 if ( q->opacity == 0 ) changed++; /* more to do */
1428               }
1429             if (((channel & IndexChannel) != 0) &&
1430                 (image->colorspace == CMYKColorspace))
1431               { q_indexes[x] = RoundToQuantum(result.index);
1432                 if ( q_indexes[x] == QuantumRange ) changed++;
1433               }
1434 #endif
1435             break;
1436
1437         case UndefinedMorphology:
1438         default:
1439             break; /* Do nothing */
1440         }
1441       p++;
1442       q++;
1443     }
1444     sync=SyncCacheViewAuthenticPixels(q_view,exception);
1445     if (sync == MagickFalse)
1446       status=MagickFalse;
1447     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1448       {
1449         MagickBooleanType
1450           proceed;
1451
1452 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1453   #pragma omp critical (MagickCore_MorphologyImage)
1454 #endif
1455         proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1456         if (proceed == MagickFalse)
1457           status=MagickFalse;
1458       }
1459   }
1460   result_image->type=image->type;
1461   q_view=DestroyCacheView(q_view);
1462   p_view=DestroyCacheView(p_view);
1463   return(status ? changed : 0);
1464 }
1465
1466
1467 MagickExport Image *MorphologyImage(const Image *image, MorphologyMethod
1468   method, const long iterations, const ChannelType channel,
1469   const MagickKernel *kernel, ExceptionInfo *exception)
1470 {
1471   unsigned long
1472     count,
1473     limit,
1474     changed;
1475
1476   Image
1477     *new_image,
1478     *old_image;
1479
1480   assert(image != (Image *) NULL);
1481   assert(image->signature == MagickSignature);
1482   assert(exception != (ExceptionInfo *) NULL);
1483   assert(exception->signature == MagickSignature);
1484
1485   if ( GetImageArtifact(image,"showkernel") != (const char *) NULL)
1486     {
1487       /* Show the Kernel that was input by the user */
1488       unsigned long
1489         i, u, v;
1490
1491       fprintf(stderr, "Kernel \"%s\" size %lux%lu%+ld%+ld scaling %+lg to %+lg\n",
1492             MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
1493             kernel->width, kernel->height,
1494             kernel->offset_x, kernel->offset_y,
1495             kernel->range_neg, kernel->range_pos);
1496       for (i=v=0; v < kernel->height; v++) {
1497         fprintf(stderr,"%2ld: ",v);
1498         for (u=0; u < kernel->width; u++, i++)
1499           fprintf(stderr,"%5.3lf ",kernel->values[i]);
1500         fprintf(stderr,"\n");
1501       }
1502     }
1503
1504   if ( iterations == 0 )
1505     return((Image *)NULL); /* null operation - nothing to do! */
1506
1507   /* kernel must be valid at this point
1508    * (except maybe for posible future morphology methods like "Prune"
1509    */
1510   assert(kernel != (MagickKernel *)NULL);
1511
1512   count = 0;
1513   limit = iterations;
1514   if ( iterations < 0 )
1515     limit = image->columns > image->rows ? image->columns : image->rows;
1516
1517   /* Special morphology cases */
1518   changed=MagickFalse;
1519   switch( method ) {
1520     case CloseMorphology:
1521       new_image = MorphologyImage(image, DialateMorphology, iterations, channel,
1522             kernel, exception);
1523       if (new_image == (Image *) NULL)
1524         return((Image *) NULL);
1525       method = ErodeMorphology;
1526       break;
1527     case OpenMorphology:
1528       new_image = MorphologyImage(image, ErodeMorphology, iterations, channel,
1529             kernel, exception);
1530       if (new_image == (Image *) NULL)
1531         return((Image *) NULL);
1532       method = DialateMorphology;
1533       break;
1534     case CloseIntensityMorphology:
1535       new_image = MorphologyImage(image, DialateIntensityMorphology,
1536             iterations, channel, kernel, exception);
1537       if (new_image == (Image *) NULL)
1538         return((Image *) NULL);
1539       method = ErodeIntensityMorphology;
1540       break;
1541     case OpenIntensityMorphology:
1542       new_image = MorphologyImage(image, ErodeIntensityMorphology,
1543             iterations, channel, kernel, exception);
1544       if (new_image == (Image *) NULL)
1545         return((Image *) NULL);
1546       method = DialateIntensityMorphology;
1547       break;
1548
1549     default:
1550       /* Do a morphology once!
1551         This ensures a new_image has been generated, but allows us
1552         to skip the creation of 'old_image' if it wasn't needed.
1553       */
1554       new_image=CloneImage(image,0,0,MagickTrue,exception);
1555       if (new_image == (Image *) NULL)
1556         return((Image *) NULL);
1557       if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1558         {
1559           InheritException(exception,&new_image->exception);
1560           new_image=DestroyImage(new_image);
1561           return((Image *) NULL);
1562         }
1563       changed = MorphologyApply(image,new_image,method,channel,kernel,
1564             exception);
1565       count++;
1566       if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1567         fprintf(stderr, "Morphology %s:%lu => Changed %lu\n",
1568               MagickOptionToMnemonic(MagickMorphologyOptions, method),
1569               count, changed);
1570   }
1571
1572   /* Repeat the interative morphology until count or no change */
1573   if ( count < limit && changed > 0 ) {
1574     old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1575     if (old_image == (Image *) NULL)
1576         return(DestroyImage(new_image));
1577     if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1578       {
1579         InheritException(exception,&old_image->exception);
1580         old_image=DestroyImage(old_image);
1581         return(DestroyImage(new_image));
1582       }
1583     while( count < limit && changed != 0 )
1584       {
1585         Image *tmp = old_image;
1586         old_image = new_image;
1587         new_image = tmp;
1588         changed = MorphologyApply(old_image,new_image,method,channel,kernel,
1589               exception);
1590         count++;
1591         if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1592           fprintf(stderr, "Morphology %s:%lu => Changed %lu\n",
1593                 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1594                 count, changed);
1595       }
1596     DestroyImage(old_image);
1597   }
1598
1599   return(new_image);
1600 }
1601