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