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