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

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