/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % M M OOO RRRR PPPP H H OOO L OOO GGGG Y Y % % MM MM O O R R P P H H O O L O O G Y Y % % M M M O O RRRR PPPP HHHHH O O L O O G GGG Y % % M M O O R R P H H O O L O O G G Y % % M M OOO R R P H H OOO LLLLL OOO GGG Y % % % % % % MagickCore Morphology Methods % % % % Software Design % % Anthony Thyssen % % September 2009 % % % % % % Copyright 1999-2010 ImageMagick Studio LLC, a non-profit organization % % dedicated to making software imaging solutions freely available. % % % % You may not use this file except in compliance with the License. You may % % obtain a copy of the License at % % % % http://www.imagemagick.org/script/license.php % % % % Unless required by applicable law or agreed to in writing, software % % distributed under the License is distributed on an "AS IS" BASIS, % % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. % % See the License for the specific language governing permissions and % % limitations under the License. % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Morpology is the the application of various kernals, of any size and even % shape, to a image in various ways (typically binary, but not always). % % Convolution (weighted sum or average) is just one specific type of % morphology. Just one that is very common for image bluring and sharpening % effects. Not only 2D Gaussian blurring, but also 2-pass 1D Blurring. % % This module provides not only a general morphology function, and the ability % to apply more advanced or iterative morphologies, but also functions for the % generation of many different types of kernel arrays from user supplied % arguments. Prehaps even the generation of a kernel from a small image. */ /* Include declarations. */ #include "magick/studio.h" #include "magick/artifact.h" #include "magick/cache-view.h" #include "magick/color-private.h" #include "magick/enhance.h" #include "magick/exception.h" #include "magick/exception-private.h" #include "magick/gem.h" #include "magick/hashmap.h" #include "magick/image.h" #include "magick/image-private.h" #include "magick/list.h" #include "magick/memory_.h" #include "magick/monitor-private.h" #include "magick/morphology.h" #include "magick/option.h" #include "magick/pixel-private.h" #include "magick/prepress.h" #include "magick/quantize.h" #include "magick/registry.h" #include "magick/semaphore.h" #include "magick/splay-tree.h" #include "magick/statistic.h" #include "magick/string_.h" #include "magick/string-private.h" #include "magick/token.h" /* * The following are assignments and tests for special floating point numbers * of value NaN (not a number), that may be used within a Kernel Definition. * NaN's are defined as part of the IEEE standard for floating point number * representation. * * These are used a Kernel value of NaN means that that kernal position * is not part of the normal convolution or morphology process, and thus * allowing the use of 'shaped' kernels. * * Special Properities * Two NaN's are never equal, even if they are from the same variable * That is the IsNaN() macro is only true if the value is NaN. */ #define IsNan(a) ((a)!=(a)) /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % A c q u i r e K e r n e l F r o m S t r i n g % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % AcquireKernelFromString() takes the given string (generally supplied by the % user) and converts it into a Morphology/Convolution Kernel. This allows % users to specify a kernel from a number of pre-defined kernels, or to fully % specify their own kernel for a specific Convolution or Morphology % Operation. % % The kernel so generated can be any rectangular array of floating point % values (doubles) with the 'control point' or 'pixel being affected' % anywhere within that array of values. % % ASIDE: Previously IM was restricted to a square of odd size using the exact % center. % % The floating point values in the kernel can also include a special value % known as 'NaN' or 'not a number' to indicate that this value is not part % of the kernel array. This allows you to specify a non-rectangular shaped % kernel, for use in Morphological operators, without the need for some type % of kernal mask. % % The returned kernel should be freed using the DestroyKernel() when you are % finished with it. % % Input kernel defintion strings can consist of any of three types. % % "num, num, num, num, ..." % list of floating point numbers defining an 'old style' odd sized % square kernel. At least 9 values should be provided for a 3x3 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc. % Values can be space or comma separated. % % "WxH[+X+Y]:num, num, num ..." % a kernal of size W by H, with W*H floating point numbers following. % the 'center' can be optionally be defined at +X+Y (such that +0+0 % is top left corner). If not defined a pixel closest to the center % of the array is automatically defined. % % "name:args" % Select from one of the built in kernels. See AcquireKernelBuiltIn() % % Note that 'name' kernels will start with an alphabetic character % while the new kernel specification has a ':' character in its % specification. % % TODO: bias and auto-scale handling of the kernel % The given kernel is assumed to have been pre-scaled appropriatally, usally % by the kernel generator. % % The format of the AcquireKernal method is: % % MagickKernel *AcquireKernelFromString(const char *kernel_string) % % A description of each parameter follows: % % o kernel_string: the Morphology/Convolution kernel wanted. % */ MagickExport MagickKernel *AcquireKernelFromString(const char *kernel_string) { MagickKernel *kernel; char token[MaxTextExtent]; register unsigned long i; const char *p; MagickStatusType flags; GeometryInfo args; assert(kernel_string != (const char *) NULL); SetGeometryInfo(&args); /* does it start with an alpha - Return a builtin kernel */ GetMagickToken(kernel_string,&p,token); if ( isalpha((int)token[0]) ) { long type; type=ParseMagickOption(MagickKernelOptions,MagickFalse,token); if ( type < 0 || type == UserDefinedKernel ) return((MagickKernel *)NULL); while (((isspace((int) ((unsigned char) *p)) != 0) || (*p == ',') || (*p == ':' )) && (*p != '\0')) p++; flags = ParseGeometry(p, &args); /* special handling of missing values in input string */ if ( type == RectangleKernel ) { if ( (flags & WidthValue) == 0 ) /* if no width then */ args.rho = args.sigma; /* then width = height */ if ( args.rho < 1.0 ) /* if width too small */ args.rho = 3; /* then width = 3 */ if ( args.sigma < 1.0 ) /* if height too small */ args.sigma = args.rho; /* then height = width */ if ( (flags & XValue) == 0 ) /* center offset if not defined */ args.xi = (double)(((long)args.rho-1)/2); if ( (flags & YValue) == 0 ) args.psi = (double)(((long)args.sigma-1)/2); } return(AcquireKernelBuiltIn((MagickKernelType)type, &args)); } kernel=(MagickKernel *) AcquireMagickMemory(sizeof(*kernel)); if (kernel == (MagickKernel *)NULL) return(kernel); (void) ResetMagickMemory(kernel,0,sizeof(*kernel)); kernel->type = UserDefinedKernel; /* Has a ':' in argument - New user kernel specification */ p = strchr(kernel_string, ':'); if ( p != (char *) NULL) { #if 1 /* ParseGeometry() needs the geometry separated! -- Arrgghh */ memcpy(token, kernel_string, p-kernel_string); token[p-kernel_string] = '\0'; flags = ParseGeometry(token, &args); #else flags = ParseGeometry(kernel_string, &args); #endif /* Size Handling and Checks */ if ( (flags & WidthValue) == 0 ) /* if no width then */ args.rho = args.sigma; /* then width = height */ if ( args.rho < 1.0 ) /* if width too small */ args.rho = 1.0; /* then width = 1 */ if ( args.sigma < 1.0 ) /* if height too small */ args.sigma = args.rho; /* then height = width */ kernel->width = (unsigned long)args.rho; kernel->height = (unsigned long)args.sigma; /* Offset Handling and Checks */ if ( args.xi < 0.0 || args.psi < 0.0 ) return(DestroyKernel(kernel)); kernel->offset_x = ((flags & XValue)!=0) ? (unsigned long)args.xi : (kernel->width-1)/2; kernel->offset_y = ((flags & YValue)!=0) ? (unsigned long)args.psi : (kernel->height-1)/2; if ( kernel->offset_x >= kernel->width || kernel->offset_y >= kernel->height ) return(DestroyKernel(kernel)); p++; /* advance beyond the ':' */ } else { /* ELSE - Old old kernel specification, forming odd-square kernel */ /* count up number of values given */ p=(const char *) kernel_string; for (i=0; *p != '\0'; i++) { GetMagickToken(p,&p,token); if (*token == ',') GetMagickToken(p,&p,token); } /* set the size of the kernel - old sized square */ kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0); kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; p=(const char *) kernel_string; } /* Read in the kernel values from rest of input string argument */ kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); kernel->range_neg = kernel->range_pos = 0.0; for (i=0; (i < kernel->width*kernel->height) && (*p != '\0'); i++) { GetMagickToken(p,&p,token); if (*token == ',') GetMagickToken(p,&p,token); (( kernel->values[i] = StringToDouble(token) ) < 0) ? ( kernel->range_neg += kernel->values[i] ) : ( kernel->range_pos += kernel->values[i] ); } for ( ; i < kernel->width*kernel->height; i++) kernel->values[i]=0.0; return(kernel); } /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % A c q u i r e K e r n e l B u i l t I n % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % AcquireKernelBuiltIn() returned one of the 'named' built-in types of % kernels used for special purposes such as gaussian blurring, skeleton % pruning, and edge distance determination. % % They take a KernelType, and a set of geometry style arguments, which were % typically decoded from a user supplied string, or from a more complex % Morphology Method that was requested. % % The format of the AcquireKernalBuiltIn method is: % % MagickKernel *AcquireKernelBuiltIn(const MagickKernelType type, % const GeometryInfo args) % % A description of each parameter follows: % % o type: the pre-defined type of kernel wanted % % o args: arguments defining or modifying the kernel % % Convolution Kernels % % Gaussian "[{radius}]x{sigma}" % Generate a two-dimentional gaussian kernel, as used by -gaussian % A sigma is required, (with the 'x'), due to historical reasons. % % NOTE: that the 'radius' is optional, but if provided can limit (clip) % the final size of the resulting kernel to a square 2*radius+1 in size. % The radius should be at least 2 times that of the sigma value, or % sever clipping and aliasing may result. If not given or set to 0 the % radius will be determined so as to produce the best minimal error % result, which is usally much larger than is normally needed. % % Blur "[{radius}]x{sigma}[+angle]" % As per Gaussian, but generates a 1 dimensional or linear gaussian % blur, at the angle given (current restricted to orthogonal angles). % If a 'radius' is given the kernel is clipped to a width of 2*radius+1. % % NOTE that two such blurs perpendicular to each other is equivelent to % -blur and the previous gaussian, but is often 10 or more times faster. % % Comet "[{width}]x{sigma}[+angle]" % Blur in one direction only, mush like how a bright object leaves % a comet like trail. The Kernel is actually half a gaussian curve, % Adding two such blurs in oppiste directions produces a Linear Blur. % % NOTE: that the first argument is the width of the kernel and not the % radius of the kernel. % % # Still to be implemented... % # % # Laplacian "{radius}x{sigma}" % # Laplacian (a mexican hat like) Function % # % # LOG "{radius},{sigma1},{sigma2} % # Laplacian of Gaussian % # % # DOG "{radius},{sigma1},{sigma2} % # Difference of Gaussians % % Boolean Kernels % % Rectangle "{geometry}" % Simply generate a rectangle of 1's with the size given. You can also % specify the location of the 'control point', otherwise the closest % pixel to the center of the rectangle is selected. % % Properly centered and odd sized rectangles work the best. % % Diamond "[{radius}]" % Generate a diamond shaped kernal with given radius to the points. % Kernel size will again be radius*2+1 square and defaults to radius 1, % generating a 3x3 kernel that is slightly larger than a square. % % Square "[{radius}]" % Generate a square shaped kernel of size radius*2+1, and defaulting % to a 3x3 (radius 1). % % Note that using a larger radius for the "Square" or the "Diamond" % is also equivelent to iterating the basic morphological method % that many times. However However iterating with the smaller radius 1 % default is actually faster than using a larger kernel radius. % % Disk "[{radius}] % Generate a binary disk of the radius given, radius may be a float. % Kernel size will be ceil(radius)*2+1 square. % NOTE: Here are some disk shapes of specific interest % "disk:1" => "diamond" or "cross:1" % "disk:1.5" => "square" % "disk:2" => "diamond:2" % "disk:2.5" => default - radius 2 disk shape % "disk:2.9" => "square:2" % "disk:3.5" => octagonal/disk shape of radius 3 % "disk:4.2" => roughly octagonal shape of radius 4 % "disk:4.3" => disk shape of radius 4 % After this all the kernel shape becomes more and more circular. % % Because a "disk" is more circular when using a larger radius, using a % larger radius is preferred over iterating the morphological operation. % % Plus "[{radius}]" % Generate a kernel in the shape of a 'plus' sign. The length of each % arm is also the radius, which defaults to 2. % % This kernel is not a good general morphological kernel, but is used % more for highlighting and marking any single pixels in an image using, % a "Dilate" or "Erode" method as appropriate. # % NOTE: "plus:1" is equivelent to a "Diamond" kernel. % % Note that unlike other kernels iterating a plus does not produce the % same result as using a larger radius for the cross. % % Distance Measuring Kernels % % Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1) % Manhatten "[{radius}][x{scale}]" square grid distance (default r=1) % Knight "[{radius}][x{scale}]" octagonal distance (default r=1) % Euclidean "[{radius}][x{scale}]" direct distance (default r=4) % % Different types of distance measuring methods, which are used with the % a 'Distance' morphology method for generating a gradient based on % distance from an edge of a binary shape, though there is a technique % for handling a anti-aliased shape. % % The first 3 are simplifications that alow the use of a small kernel % which is iterated. The lest is more accurate but requires a larger % kernel to produce a accurate distance measure. The larger the better. % % The actual distance is scaled the size give, which while unnecessary % for a "Chebyshev" or "Manhatten" distance, is needed to allow for % correct handling of fractional distances in "Knight" and "Euclidean" % distance formulas. If no scale is provided it is set to a value of % 100, allowing for a maximum distance measurement of 655 pixels from % any edge, using a Q16 version of IM. % % See the 'Distance' Morphological Method, for information of how it % is applied. % */ static void KernelRotate(MagickKernel *kernel, double angle) { /* Rotate a kernel appropriately for the angle given ** ** Currently assumes the kernel (rightly) horizontally is symetrical ** ** TODO: expand beyond simple 90 degree rotates, flips and flops */ /* Modulus the angle */ angle = fmod(angle, 360.0); if ( angle < 0 ) angle += 360.0; if ( 315.0 < angle || angle <= 45.0 ) return; /* no change! - At least at this time */ switch (kernel->type) { /* These kernels are cylindrical kernel, rotating is useless */ case GaussianKernel: case LaplacianKernel: case LOGKernel: case DOGKernel: case DiskKernel: case ChebyshevKernel: case ManhattenKernel: case KnightKernel: case EuclideanKernel: return; /* These may be rotatable at non-90 angles in the future */ /* but simply rotating them 90 degrees is useless */ case SquareKernel: case DiamondKernel: case PlusKernel: return; /* These only allows a +/-90 degree rotation (transpose) */ case BlurKernel: case RectangleKernel: if ( 135.0 < angle && angle <= 225.0 ) return; if ( 225.0 < angle && angle <= 315.0 ) angle -= 180; break; /* these are freely rotatable in 90 degree units */ case CometKernel: case UndefinedKernel: case UserDefinedKernel: break; } fprintf(stderr, "angle2 = %lf\n", angle); if ( 135.0 < angle && angle <= 315.0 ) { /* Do a flop, this assumes kernel is horizontally symetrical. */ /* Each kernel data row need to be reversed! */ unsigned long y; register unsigned long x,r; register double *k,t; for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width) { for ( x=0, r=kernel->width-1; xwidth/2; x++, r--) t=k[x], k[x]=k[r], k[r]=t; } kernel->offset_x = kernel->width - kernel->offset_x - 1; angle = fmod(angle+180.0, 360.0); } if ( 45.0 < angle && angle <= 135.0 ) { /* Do a transpose, this assumes the kernel is orthoginally symetrical */ /* The data is the same, just the size and offsets needs to be swapped. */ unsigned long t; t = kernel->width; kernel->width = kernel->height; kernel->height = t; t = kernel->offset_x; kernel->offset_x = kernel->offset_y; kernel->offset_y = t; angle = fmod(450.0 - angle, 360.0); } /* at this point angle should be between +45 and -45 (315) degrees */ return; } MagickExport MagickKernel *AcquireKernelBuiltIn(const MagickKernelType type, const GeometryInfo *args) { MagickKernel *kernel; register unsigned long i; register long u, v; double nan = sqrt((double)-1.0); /* Special Value : Not A Number */ kernel=(MagickKernel *) AcquireMagickMemory(sizeof(*kernel)); if (kernel == (MagickKernel *) NULL) return(kernel); (void) ResetMagickMemory(kernel,0,sizeof(*kernel)); kernel->range_neg = kernel->range_pos = 0.0; kernel->type = type; switch(type) { /* Convolution Kernels */ case GaussianKernel: { double sigma = fabs(args->sigma); sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma; kernel->width = kernel->height = GetOptimalKernelWidth2D(args->rho,sigma); kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; kernel->range_neg = kernel->range_pos = 0.0; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); sigma = 2.0*sigma*sigma; /* simplify the expression */ for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++) for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++) kernel->range_pos += ( kernel->values[i] = exp(-((double)(u*u+v*v))/sigma) /* / (MagickPI*sigma) */ ); /* Normalize the Kernel - see notes in BlurKernel, below */ u=kernel->width*kernel->height; for (i=0; i < (unsigned long)u; i++) kernel->values[i] /= kernel->range_pos; kernel->range_pos=1.0; break; } case BlurKernel: { double sigma = fabs(args->sigma); sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma; kernel->width = GetOptimalKernelWidth1D(args->rho,sigma); kernel->offset_x = (kernel->width-1)/2; kernel->height = 1; kernel->offset_y = 0; kernel->range_neg = kernel->range_pos = 0.0; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); #if 1 #define KernelRank 3 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix). ** It generates a gaussian 3 times the width, and compresses it into ** the expected range. This produces a closer normalization of the ** resulting kernel, especially for very low sigma values. ** As such while wierd it is prefered. ** ** I am told this method originally came from Photoshop. */ sigma *= KernelRank; /* simplify expanded curve */ v = (kernel->width*KernelRank-1)/2; /* start/end points to fit range */ (void) ResetMagickMemory(kernel->values,0, (size_t) kernel->width*sizeof(double)); for ( u=-v; u <= v; u++) { kernel->values[(u+v)/KernelRank] += exp(-((double)(u*u))/(2.0*sigma*sigma)) /* / (MagickSQ2PI*sigma/KernelRank) */ ; } for (i=0; i < kernel->width; i++) kernel->range_pos += kernel->values[i]; #else for ( i=0, u=-kernel->offset_x; i < kernel->width; i++, u++) kernel->range_pos += ( kernel->values[i] = exp(-((double)(u*u))/(2.0*sigma*sigma)) /* / (MagickSQ2PI*sigma) */ ); #endif /* Note that both the above methods do not generate a normalized ** kernel, though it gets close. The kernel may be 'clipped' by a user ** defined radius, producing a smaller (darker) kernel. Also for very ** small sigma's (> 0.1) the central value becomes larger than one, ** and thus producing a bright kernel. */ #if 1 /* Normalize the 1D Gaussian Kernel ** ** Because of this the divisor in the above kernel generator is ** not needed, and is taken care of here. */ for (i=0; i < kernel->width; i++) kernel->values[i] /= kernel->range_pos; kernel->range_pos=1.0; #endif /* rotate the kernel by given angle */ KernelRotate(kernel, args->xi); break; } case CometKernel: { double sigma = fabs(args->sigma); sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma; if ( args->rho < 1.0 ) kernel->width = GetOptimalKernelWidth1D(args->rho,sigma); else kernel->width = (unsigned long)args->rho; kernel->offset_x = kernel->offset_y = 0; kernel->height = 1; kernel->range_neg = kernel->range_pos = 0.0; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); /* A comet blur is half a gaussian curve, so that the object is ** blurred in one direction only. This may not be quite the right ** curve so may change in the future. The function must be normalised. */ #if 1 #define KernelRank 3 sigma *= KernelRank; /* simplify expanded curve */ v = kernel->width*KernelRank; /* start/end points to fit range */ (void) ResetMagickMemory(kernel->values,0, (size_t) kernel->width*sizeof(double)); for ( u=0; u < v; u++) { kernel->values[u/KernelRank] += exp(-((double)(u*u))/(2.0*sigma*sigma)) /* / (MagickSQ2PI*sigma/KernelRank) */ ; } for (i=0; i < kernel->width; i++) kernel->range_pos += kernel->values[i]; #else for ( i=0; i < kernel->width; i++) kernel->range_pos += ( kernel->values[i] = exp(-((double)(i*i))/(2.0*sigma*sigma)) /* / (MagickSQ2PI*sigma) */ ); #endif /* Normalize the Kernel - see notes in BlurKernel */ for (i=0; i < kernel->width; i++) kernel->values[i] /= kernel->range_pos; kernel->range_pos=1.0; /* rotate the kernel by given angle */ KernelRotate(kernel, args->xi); break; } /* Boolean Kernels */ case RectangleKernel: case SquareKernel: { if ( type == SquareKernel ) { if (args->rho < 1.0) kernel->width = kernel->height = 3; /* radius 1 */ else kernel->width = kernel->height = 2*(long)args->rho+1; kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; } else { if ( args->rho < 1.0 || args->sigma < 1.0 ) return(DestroyKernel(kernel)); kernel->width = (unsigned long)args->rho; kernel->height = (unsigned long)args->sigma; if ( args->xi < 0.0 || args->xi > (double)kernel->width || args->psi < 0.0 || args->psi > (double)kernel->height ) return(DestroyKernel(kernel)); kernel->offset_x = (unsigned long)args->xi; kernel->offset_y = (unsigned long)args->psi; } kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); u=kernel->width*kernel->height; for ( i=0; i < (unsigned long)u; i++) kernel->values[i] = 1.0; break; } case DiamondKernel: { if (args->rho < 1.0) kernel->width = kernel->height = 3; /* radius 1 */ else kernel->width = kernel->height = ((unsigned long)args->rho)*2+1; kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++) for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++) if ((labs(u)+labs(v)) <= (long)kernel->offset_x) kernel->range_pos += kernel->values[i] = 1.0; else kernel->values[i] = nan; break; } case DiskKernel: { long limit; limit = (long)(args->rho*args->rho); if (args->rho < 1.0) /* default: ~2.5 radius disk */ kernel->width = kernel->height = 5L, limit = 5L; else kernel->width = kernel->height = ((unsigned long)args->rho)*2+1; kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++) for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++) if ((u*u+v*v) <= limit) kernel->range_pos += kernel->values[i] = 1.0; else kernel->values[i] = nan; break; } case PlusKernel: { if (args->rho < 1.0) kernel->width = kernel->height = 5; /* radius 2 */ else kernel->width = kernel->height = ((unsigned long)args->rho)*2+1; kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++) for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++) kernel->values[i] = (u == 0 || v == 0) ? 1.0 : nan; kernel->range_pos = kernel->width*2.0 - 1.0; break; } /* Distance Measuring Kernels */ case ChebyshevKernel: { double scale; if (args->rho < 1.0) kernel->width = kernel->height = 3; else kernel->width = kernel->height = ((unsigned long)args->rho)*2+1; kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); scale = (args->sigma < 1.0) ? 100.0 : args->sigma; for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++) for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++) kernel->range_pos += ( kernel->values[i] = scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) ); break; } case ManhattenKernel: { double scale; if (args->rho < 1.0) kernel->width = kernel->height = 3; else kernel->width = kernel->height = ((unsigned long)args->rho)*2+1; kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); scale = (args->sigma < 1.0) ? 100.0 : args->sigma; for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++) for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++) kernel->range_pos += ( kernel->values[i] = scale*(labs(u)+labs(v)) ); break; } case KnightKernel: { double scale; if (args->rho < 1.0) kernel->width = kernel->height = 3; else kernel->width = kernel->height = ((unsigned long)args->rho)*2+1; kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); scale = (args->sigma < 1.0) ? 100.0 : args->sigma; for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++) for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++) kernel->range_pos += ( kernel->values[i] = scale*((labs(u)rho < 1.0) kernel->width = kernel->height = 9; else kernel->width = kernel->height = ((unsigned long)args->rho)*2+1; kernel->offset_x = kernel->offset_y = (kernel->width-1)/2; kernel->values=(double *) AcquireQuantumMemory(kernel->width, kernel->height*sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); scale = (args->sigma < 1.0) ? 100.0 : args->sigma; for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++) for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++) kernel->range_pos += ( kernel->values[i] = scale*sqrt((double)(u*u+v*v)) ); break; } /* Undefined Kernels */ case LaplacianKernel: case LOGKernel: case DOGKernel: assert("Kernel Type has not been defined yet"); /* FALL THRU */ default: /* Generate a No-Op minimal kernel - 1x1 pixel */ kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double)); if (kernel->values == (double *) NULL) return(DestroyKernel(kernel)); kernel->range_pos = kernel->values[0] = 1.0; kernel->width = kernel->height = 1; kernel->offset_x = kernel->offset_x = 0; kernel->type = UndefinedKernel; break; } return(kernel); } /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % D e s t r o y K e r n e l % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % DestroyKernel() frees the memory used by a Convolution/Morphology kernel. % % The format of the DestroyKernel method is: % % MagickKernel *DestroyKernel(MagickKernel *kernel) % % A description of each parameter follows: % % o kernel: the Morphology/Convolution kernel to be destroyed % */ MagickExport MagickKernel *DestroyKernel(MagickKernel *kernel) { assert(kernel != (MagickKernel *) NULL); kernel->values=(double *)RelinquishMagickMemory(kernel->values); kernel=(MagickKernel *) RelinquishMagickMemory(kernel); return(kernel); } /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % % % % M o r p h o l o g y I m a g e % % % % % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % MorphologyImage() applies a user supplied kernel to the image according to % the given mophology method. % % The given kernel is assumed to have been pre-scaled appropriatally, usally % by the kernel generator. % % The format of the MorphologyImage method is: % % Image *MorphologyImage(const Image *image, const MorphologyMethod % method, const long iterations, const ChannelType channel, % const MagickKernel *kernel, ExceptionInfo *exception) % % A description of each parameter follows: % % o image: the image. % % o method: the morphology method to be applied. % % o iterations: apply the operation this many times (or no change). % A value of -1 means loop until no change found. % How this is applied may depend on the morphology method. % Typically this is a value of 1. % % o channel: the channel type. % % o kernel: An array of double representing the morphology kernel. % This is assumed to have been pre-scaled (normalized). % % o exception: return any errors or warnings in this structure. % % % TODO: bias and auto-scale handling of the kernel for convolution % The given kernel is assumed to have been pre-scaled appropriatally, usally % by the kernel generator. % */ static inline double MagickMin(const MagickRealType x,const MagickRealType y) { return( x < y ? x : y); } static inline double MagickMax(const MagickRealType x,const MagickRealType y) { return( x > y ? x : y); } #define Minimize(assign,value) assign=MagickMin(assign,value) #define Maximize(assign,value) assign=MagickMax(assign,value) /* incr change if the value being assigned changed */ #define Assign(channel,value) \ { q->channel = RoundToQuantum(value); \ if ( p[r].channel != q->channel ) changed++; \ } #define AssignIndex(value) \ { q_indexes[x] = RoundToQuantum(value); \ if ( p_indexes[r] != q_indexes[x] ) changed++; \ } /* Internal function * Apply the Morphology method with the given Kernel * And return the number of values changed. */ static unsigned long MorphologyApply(const Image *image, Image *result_image, const MorphologyMethod method, const ChannelType channel, const MagickKernel *kernel, ExceptionInfo *exception) { #define MorphologyTag "Morphology/Image" long progress, y; unsigned long changed; MagickBooleanType status; MagickPixelPacket bias; CacheView *p_view, *q_view; /* Apply Morphology to Image. */ status=MagickTrue; changed=0; progress=0; GetMagickPixelPacket(image,&bias); SetMagickPixelPacketBias(image,&bias); p_view=AcquireCacheView(image); q_view=AcquireCacheView(result_image); #if defined(MAGICKCORE_OPENMP_SUPPORT) #pragma omp parallel for schedule(dynamic,4) shared(progress,status) #endif for (y=0; y < (long) image->rows; y++) { MagickBooleanType sync; register const PixelPacket *restrict p; register const IndexPacket *restrict p_indexes; register PixelPacket *restrict q; register IndexPacket *restrict q_indexes; register long x; long r; if (status == MagickFalse) continue; p=GetCacheViewVirtualPixels(p_view, -kernel->offset_x, y-kernel->offset_y, image->columns+kernel->width, kernel->height, exception); q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1, exception); if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL)) { status=MagickFalse; continue; } p_indexes=GetCacheViewVirtualIndexQueue(p_view); q_indexes=GetCacheViewAuthenticIndexQueue(q_view); r = (image->columns+kernel->width)*kernel->offset_y+kernel->offset_x; for (x=0; x < (long) image->columns; x++) { long v; register long u; register const double *restrict k; register const PixelPacket *restrict k_pixels; register const IndexPacket *restrict k_indexes; MagickPixelPacket result; /* Copy input to ouput image - removes need for 'cloning' new images */ *q = p[r]; if (image->colorspace == CMYKColorspace) q_indexes[x] = p_indexes[r]; result.index=0; switch (method) { case ConvolveMorphology: result=bias; break; /* default result is the convolution bias */ case DialateIntensityMorphology: case ErodeIntensityMorphology: /* result is the pixel as is */ result.red = p[r].red; result.green = p[r].green; result.blue = p[r].blue; result.opacity = p[r].opacity; if ( image->colorspace == CMYKColorspace) result.index = p_indexes[r]; break; default: /* most need to handle transparency as alpha */ result.red = p[r].red; result.green = p[r].green; result.blue = p[r].blue; result.opacity = QuantumRange - p[r].opacity; if ( image->colorspace == CMYKColorspace) result.index = p_indexes[r]; break; } switch ( method ) { case ConvolveMorphology: /* Weighted Average of pixels */ if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse)) { /* Kernel Weighted Convolution (no transparency) */ k = kernel->values; k_pixels = p; k_indexes = p_indexes; for (v=0; v < (long) kernel->height; v++) { for (u=0; u < (long) kernel->width; u++, k++) { if ( IsNan(*k) ) continue; result.red += (*k)*k_pixels[u].red; result.green += (*k)*k_pixels[u].green; result.blue += (*k)*k_pixels[u].blue; /* result.opacity += no involvment */ if ( image->colorspace == CMYKColorspace) result.index += (*k)*k_indexes[u]; } k_pixels += image->columns+kernel->width; k_indexes += image->columns+kernel->width; } if ((channel & RedChannel) != 0) Assign(red,result.red); if ((channel & GreenChannel) != 0) Assign(green,result.green); if ((channel & BlueChannel) != 0) Assign(blue,result.blue); /* no transparency involved */ if ((channel & IndexChannel) != 0 && image->colorspace == CMYKColorspace) AssignIndex(result.index); } else { /* Kernel & Alpha weighted Convolution */ MagickRealType alpha, /* alpha value * kernel weighting */ gamma; /* weighting divisor */ gamma=0.0; k = kernel->values; k_pixels = p; k_indexes = p_indexes; for (v=0; v < (long) kernel->height; v++) { for (u=0; u < (long) kernel->width; u++, k++) { if ( IsNan(*k) ) continue; alpha=(*k)*(QuantumScale*(QuantumRange- k_pixels[u].opacity)); gamma += alpha; result.red += alpha*k_pixels[u].red; result.green += alpha*k_pixels[u].green; result.blue += alpha*k_pixels[u].blue; result.opacity += (*k)*k_pixels[u].opacity; if ( image->colorspace == CMYKColorspace) result.index += alpha*k_indexes[u]; } k_pixels += image->columns+kernel->width; k_indexes += image->columns+kernel->width; } gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma); if ((channel & RedChannel) != 0) Assign(red,gamma*result.red); if ((channel & GreenChannel) != 0) Assign(green,gamma*result.green); if ((channel & BlueChannel) != 0) Assign(blue,gamma*result.blue); if ((channel & OpacityChannel) != 0 && image->matte == MagickTrue ) Assign(opacity,result.opacity); if ((channel & IndexChannel) != 0 && image->colorspace == CMYKColorspace) AssignIndex(gamma*result.index); } break; case DialateMorphology: /* Maximize Value - Kernel should be boolean */ k = kernel->values; k_pixels = p; k_indexes = p_indexes; for (v=0; v < (long) kernel->height; v++) { for (u=0; u < (long) kernel->width; u++, k++) { if ( IsNan(*k) || (*k) < 0.5 ) continue; Maximize(result.red, k_pixels[u].red); Maximize(result.green, k_pixels[u].green); Maximize(result.blue, k_pixels[u].blue); Maximize(result.opacity, QuantumRange-k_pixels[u].opacity); if ( image->colorspace == CMYKColorspace) Maximize(result.index, k_indexes[u]); } k_pixels += image->columns+kernel->width; k_indexes += image->columns+kernel->width; } if ((channel & RedChannel) != 0) Assign(red,result.red); if ((channel & GreenChannel) != 0) Assign(green,result.green); if ((channel & BlueChannel) != 0) Assign(blue,result.blue); if ((channel & OpacityChannel) != 0 && image->matte == MagickTrue ) Assign(opacity,QuantumRange-result.opacity); if ((channel & IndexChannel) != 0 && image->colorspace == CMYKColorspace) AssignIndex(result.index); break; case ErodeMorphology: /* Minimize Value - Kernel should be boolean */ k = kernel->values; k_pixels = p; k_indexes = p_indexes; for (v=0; v < (long) kernel->height; v++) { for (u=0; u < (long) kernel->width; u++, k++) { if ( IsNan(*k) || (*k) < 0.5 ) continue; Minimize(result.red, k_pixels[u].red); Minimize(result.green, k_pixels[u].green); Minimize(result.blue, k_pixels[u].blue); Minimize(result.opacity, QuantumRange-k_pixels[u].opacity); if ( image->colorspace == CMYKColorspace) Minimize(result.index, k_indexes[u]); } k_pixels += image->columns+kernel->width; k_indexes += image->columns+kernel->width; } if ((channel & RedChannel) != 0) Assign(red,result.red); if ((channel & GreenChannel) != 0) Assign(green,result.green); if ((channel & BlueChannel) != 0) Assign(blue,result.blue); if ((channel & OpacityChannel) != 0 && image->matte == MagickTrue ) Assign(opacity,QuantumRange-result.opacity); if ((channel & IndexChannel) != 0 && image->colorspace == CMYKColorspace) AssignIndex(result.index); break; case DialateIntensityMorphology: /* Maximum Intensity Pixel - Kernel should be boolean */ k = kernel->values; k_pixels = p; k_indexes = p_indexes; for (v=0; v < (long) kernel->height; v++) { for (u=0; u < (long) kernel->width; u++, k++) { if ( IsNan(*k) || (*k) < 0.5 ) continue; if ( PixelIntensity(&p[r]) > PixelIntensity(&(k_pixels[u])) ) continue; result.red = k_pixels[u].red; result.green = k_pixels[u].green; result.blue = k_pixels[u].blue; result.opacity = k_pixels[u].opacity; if ( image->colorspace == CMYKColorspace) result.index = k_indexes[u]; } k_pixels += image->columns+kernel->width; k_indexes += image->columns+kernel->width; } if ((channel & RedChannel) != 0) Assign(red,result.red); if ((channel & GreenChannel) != 0) Assign(green,result.green); if ((channel & BlueChannel) != 0) Assign(blue,result.blue); if ((channel & OpacityChannel) != 0 && image->matte == MagickTrue ) Assign(opacity,result.opacity); if ((channel & IndexChannel) != 0 && image->colorspace == CMYKColorspace) AssignIndex(result.index); break; case ErodeIntensityMorphology: /* Minimum Intensity Pixel - Kernel should be boolean */ k = kernel->values; k_pixels = p; k_indexes = p_indexes; for (v=0; v < (long) kernel->height; v++) { for (u=0; u < (long) kernel->width; u++, k++) { if ( IsNan(*k) || (*k) < 0.5 ) continue; if ( PixelIntensity(&p[r]) < PixelIntensity(&(k_pixels[u])) ) continue; result.red = k_pixels[u].red; result.green = k_pixels[u].green; result.blue = k_pixels[u].blue; result.opacity = k_pixels[u].opacity; if ( image->colorspace == CMYKColorspace) result.index = k_indexes[u]; } k_pixels += image->columns+kernel->width; k_indexes += image->columns+kernel->width; } if ((channel & RedChannel) != 0) Assign(red,result.red); if ((channel & GreenChannel) != 0) Assign(green,result.green); if ((channel & BlueChannel) != 0) Assign(blue,result.blue); if ((channel & OpacityChannel) != 0 && image->matte == MagickTrue ) Assign(opacity,result.opacity); if ((channel & IndexChannel) != 0 && image->colorspace == CMYKColorspace) AssignIndex(result.index); break; case DistanceMorphology: #if 0 /* No need to do distance morphology if all values are zero */ /* Unfortunatally I have not been able to get this right! */ if ( ((channel & RedChannel) == 0 && p[r].red == 0) || ((channel & GreenChannel) == 0 && p[r].green == 0) || ((channel & BlueChannel) == 0 && p[r].blue == 0) || ((channel & OpacityChannel) == 0 && p[r].opacity == 0) || (( (channel & IndexChannel) == 0 || image->colorspace != CMYKColorspace ) && p_indexes[x] ==0 ) ) break; #endif k = kernel->values; k_pixels = p; k_indexes = p_indexes; for (v=0; v < (long) kernel->height; v++) { for (u=0; u < (long) kernel->width; u++, k++) { if ( IsNan(*k) ) continue; Minimize(result.red, (*k)+k_pixels[u].red); Minimize(result.green, (*k)+k_pixels[u].green); Minimize(result.blue, (*k)+k_pixels[u].blue); Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity); if ( image->colorspace == CMYKColorspace) Minimize(result.index, (*k)+k_indexes[u]); } k_pixels += image->columns+kernel->width; k_indexes += image->columns+kernel->width; } #if 1 if ((channel & RedChannel) != 0) Assign(red,result.red); if ((channel & GreenChannel) != 0) Assign(green,result.green); if ((channel & BlueChannel) != 0) Assign(blue,result.blue); if ((channel & OpacityChannel) != 0 && image->matte == MagickTrue ) Assign(opacity,QuantumRange-result.opacity); if ((channel & IndexChannel) != 0 && image->colorspace == CMYKColorspace) AssignIndex(result.index); #else /* By returning the number of 'maximum' values still to process ** we can get the Distance iteration to finish faster. ** BUT this may cause an infinite loop on very large shapes, ** which may have a distance that reachs a maximum gradient. */ if ((channel & RedChannel) != 0) { q->red = RoundToQuantum(result.red); if ( q->red == QuantumRange ) changed++; /* more to do */ } if ((channel & GreenChannel) != 0) { q->green = RoundToQuantum(result.green); if ( q->green == QuantumRange ) changed++; /* more to do */ } if ((channel & BlueChannel) != 0) { q->blue = RoundToQuantum(result.blue); if ( q->blue == QuantumRange ) changed++; /* more to do */ } if ((channel & OpacityChannel) != 0) { q->opacity = RoundToQuantum(QuantumRange-result.opacity); if ( q->opacity == 0 ) changed++; /* more to do */ } if (((channel & IndexChannel) != 0) && (image->colorspace == CMYKColorspace)) { q_indexes[x] = RoundToQuantum(result.index); if ( q_indexes[x] == QuantumRange ) changed++; } #endif break; case UndefinedMorphology: default: break; /* Do nothing */ } p++; q++; } sync=SyncCacheViewAuthenticPixels(q_view,exception); if (sync == MagickFalse) status=MagickFalse; if (image->progress_monitor != (MagickProgressMonitor) NULL) { MagickBooleanType proceed; #if defined(MAGICKCORE_OPENMP_SUPPORT) #pragma omp critical (MagickCore_MorphologyImage) #endif proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows); if (proceed == MagickFalse) status=MagickFalse; } } result_image->type=image->type; q_view=DestroyCacheView(q_view); p_view=DestroyCacheView(p_view); return(status ? changed : 0); } MagickExport Image *MorphologyImage(const Image *image, MorphologyMethod method, const long iterations, const ChannelType channel, const MagickKernel *kernel, ExceptionInfo *exception) { unsigned long count, limit, changed; Image *new_image, *old_image; assert(image != (Image *) NULL); assert(image->signature == MagickSignature); assert(exception != (ExceptionInfo *) NULL); assert(exception->signature == MagickSignature); if ( GetImageArtifact(image,"showkernel") != (const char *) NULL) { /* Show the Kernel that was input by the user */ unsigned long i, u, v; fprintf(stderr, "Kernel \"%s\" size %lux%lu%+ld%+ld scaling %+lg to %+lg\n", MagickOptionToMnemonic(MagickKernelOptions, kernel->type), kernel->width, kernel->height, kernel->offset_x, kernel->offset_y, kernel->range_neg, kernel->range_pos); for (i=v=0; v < kernel->height; v++) { fprintf(stderr,"%2ld: ",v); for (u=0; u < kernel->width; u++, i++) fprintf(stderr,"%5.3lf ",kernel->values[i]); fprintf(stderr,"\n"); } } if ( iterations == 0 ) return((Image *)NULL); /* null operation - nothing to do! */ /* kernel must be valid at this point * (except maybe for posible future morphology methods like "Prune" */ assert(kernel != (MagickKernel *)NULL); count = 0; limit = iterations; if ( iterations < 0 ) limit = image->columns > image->rows ? image->columns : image->rows; /* Special morphology cases */ changed=MagickFalse; switch( method ) { case CloseMorphology: new_image = MorphologyImage(image, DialateMorphology, iterations, channel, kernel, exception); if (new_image == (Image *) NULL) return((Image *) NULL); method = ErodeMorphology; break; case OpenMorphology: new_image = MorphologyImage(image, ErodeMorphology, iterations, channel, kernel, exception); if (new_image == (Image *) NULL) return((Image *) NULL); method = DialateMorphology; break; case CloseIntensityMorphology: new_image = MorphologyImage(image, DialateIntensityMorphology, iterations, channel, kernel, exception); if (new_image == (Image *) NULL) return((Image *) NULL); method = ErodeIntensityMorphology; break; case OpenIntensityMorphology: new_image = MorphologyImage(image, ErodeIntensityMorphology, iterations, channel, kernel, exception); if (new_image == (Image *) NULL) return((Image *) NULL); method = DialateIntensityMorphology; break; default: /* Do a morphology once! This ensures a new_image has been generated, but allows us to skip the creation of 'old_image' if it wasn't needed. */ new_image=CloneImage(image,0,0,MagickTrue,exception); if (new_image == (Image *) NULL) return((Image *) NULL); if (SetImageStorageClass(new_image,DirectClass) == MagickFalse) { InheritException(exception,&new_image->exception); new_image=DestroyImage(new_image); return((Image *) NULL); } changed = MorphologyApply(image,new_image,method,channel,kernel, exception); count++; if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) fprintf(stderr, "Morphology %s:%lu => Changed %lu\n", MagickOptionToMnemonic(MagickMorphologyOptions, method), count, changed); } /* Repeat the interative morphology until count or no change */ if ( count < limit && changed > 0 ) { old_image = CloneImage(new_image,0,0,MagickTrue,exception); if (old_image == (Image *) NULL) return(DestroyImage(new_image)); if (SetImageStorageClass(old_image,DirectClass) == MagickFalse) { InheritException(exception,&old_image->exception); old_image=DestroyImage(old_image); return(DestroyImage(new_image)); } while( count < limit && changed != 0 ) { Image *tmp = old_image; old_image = new_image; new_image = tmp; changed = MorphologyApply(old_image,new_image,method,channel,kernel, exception); count++; if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) fprintf(stderr, "Morphology %s:%lu => Changed %lu\n", MagickOptionToMnemonic(MagickMorphologyOptions, method), count, changed); } DestroyImage(old_image); } return(new_image); }