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