2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %
13 % MagickCore Morphology Methods %
20 % Copyright 1999-2010 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % http://www.imagemagick.org/script/license.php %
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. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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).
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.
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.
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"
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.
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.
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.
92 #define IsNan(a) ((a)!=(a))
96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100 % A c q u i r e K e r n e l F r o m S t r i n g %
104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 % AcquireKernelFromString() takes the given string (generally supplied by the
107 % user) and converts it into a Morphology/Convolution Kernel. This allows
108 % users to specify a kernel from a number of pre-defined kernels, or to fully
109 % specify their own kernel for a specific Convolution or Morphology
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.
116 % ASIDE: Previously IM was restricted to a square of odd size using the exact
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
125 % The returned kernel should be freed using the DestroyKernel() when you are
128 % Input kernel defintion strings can consist of any of three types.
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.
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.
143 % Select from one of the built in kernels. See AcquireKernelBuiltIn()
145 % Note that 'name' kernels will start with an alphabetic character
146 % while the new kernel specification has a ':' character in its
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.
153 % The format of the AcquireKernal method is:
155 % MagickKernel *AcquireKernelFromString(const char *kernel_string)
157 % A description of each parameter follows:
159 % o kernel_string: the Morphology/Convolution kernel wanted.
163 MagickExport MagickKernel *AcquireKernelFromString(const char *kernel_string)
169 token[MaxTextExtent];
171 register unsigned long
183 assert(kernel_string != (const char *) NULL);
184 SetGeometryInfo(&args);
186 /* does it start with an alpha - Return a builtin kernel */
187 GetMagickToken(kernel_string,&p,token);
188 if ( isalpha((int)token[0]) )
193 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
194 if ( type < 0 || type == UserDefinedKernel )
195 return((MagickKernel *)NULL);
197 while (((isspace((int) ((unsigned char) *p)) != 0) ||
198 (*p == ',') || (*p == ':' )) && (*p != '\0'))
200 flags = ParseGeometry(p, &args);
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);
216 return(AcquireKernelBuiltIn((MagickKernelType)type, &args));
219 kernel=(MagickKernel *) AcquireMagickMemory(sizeof(*kernel));
220 if (kernel == (MagickKernel *)NULL)
222 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
223 kernel->type = UserDefinedKernel;
224 kernel->signature = MagickSignature;
226 /* Has a ':' in argument - New user kernel specification */
227 p = strchr(kernel_string, ':');
228 if ( p != (char *) NULL)
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);
236 flags = ParseGeometry(kernel_string, &args);
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;
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));
260 p++; /* advance beyond the ':' */
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 == '\''))
268 for (i=0; *p != '\0'; i++)
270 GetMagickToken(p,&p,token);
272 GetMagickToken(p,&p,token);
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;
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));
286 kernel->range_neg = kernel->range_pos = 0.0;
287 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
289 for (i=0; (i < kernel->width*kernel->height) && (*p != '\0'); i++)
291 GetMagickToken(p,&p,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] );
298 for ( ; i < kernel->width*kernel->height; i++)
299 kernel->values[i]=0.0;
305 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
309 % A c q u i r e K e r n e l B u i l t I n %
313 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
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.
323 % The format of the AcquireKernalBuiltIn method is:
325 % MagickKernel *AcquireKernelBuiltIn(const MagickKernelType type,
326 % const GeometryInfo args)
328 % A description of each parameter follows:
330 % o type: the pre-defined type of kernel wanted
332 % o args: arguments defining or modifying the kernel
334 % Convolution Kernels
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.
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.
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.
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.
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.
360 % NOTE: that the first argument is the width of the kernel and not the
361 % radius of the kernel.
363 % # Still to be implemented...
365 % # Laplacian "{radius}x{sigma}"
366 % # Laplacian (a mexican hat like) Function
368 % # LOG "{radius},{sigma1},{sigma2}
369 % # Laplacian of Gaussian
371 % # DOG "{radius},{sigma1},{sigma2}
372 % # Difference of Gaussians
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.
381 % Properly centered and odd sized rectangles work the best.
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.
388 % Square "[{radius}]"
389 % Generate a square shaped kernel of size radius*2+1, and defaulting
390 % to a 3x3 (radius 1).
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.
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.
411 % Because a "disk" is more circular when using a larger radius, using a
412 % larger radius is preferred over iterating the morphological operation.
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.
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.
422 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
424 % Note that unlike other kernels iterating a plus does not produce the
425 % same result as using a larger radius for the cross.
427 % Distance Measuring Kernels
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)
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.
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
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.
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
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.
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.
471 % See the 'Distance' Morphological Method, for information of how it is
476 MagickExport MagickKernel *AcquireKernelBuiltIn(const MagickKernelType type,
477 const GeometryInfo *args)
482 register unsigned long
490 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
492 kernel=(MagickKernel *) AcquireMagickMemory(sizeof(*kernel));
493 if (kernel == (MagickKernel *) NULL)
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;
499 kernel->signature = MagickSignature;
502 /* Convolution Kernels */
505 sigma = fabs(args->sigma);
507 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
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));
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 += (
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 ];
529 KernelNormalize(kernel);
535 sigma = fabs(args->sigma);
537 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
539 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
540 kernel->offset_x = (kernel->width-1)/2;
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));
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.
557 ** I am told this method originally came from Photoshop.
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) */ ;
568 for (i=0; i < kernel->width; i++)
569 kernel->range_pos += kernel->values[i];
571 for ( i=0, u=-kernel->offset_x; i < kernel->width; i++, u++)
572 kernel->range_pos += (
574 exp(-((double)(u*u))/(2.0*sigma*sigma))
575 /* / (MagickSQ2PI*sigma) */ );
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.
586 /* Normalize the 1D Gaussian Kernel
588 ** Because of this the divisor in the above kernel generator is
589 ** not needed, so is not done above.
591 KernelNormalize(kernel);
593 /* rotate the kernel by given angle */
594 KernelRotate(kernel, args->xi);
599 sigma = fabs(args->sigma);
601 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
603 if ( args->rho < 1.0 )
604 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
606 kernel->width = (unsigned long)args->rho;
607 kernel->offset_x = kernel->offset_y = 0;
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));
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.
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) */ ;
630 for (i=0; i < kernel->width; i++)
631 kernel->range_pos += kernel->values[i];
633 for ( i=0; i < kernel->width; i++)
634 kernel->range_pos += (
636 exp(-((double)(i*i))/(2.0*sigma*sigma))
637 /* / (MagickSQ2PI*sigma) */ );
639 kernel->value_min = 0;
640 kernel->value_max = kernel->values[0];
642 KernelNormalize(kernel);
643 KernelRotate(kernel, args->xi);
646 /* Boolean Kernels */
647 case RectangleKernel:
650 if ( type == SquareKernel )
653 kernel->width = kernel->height = 3; /* default radius = 1 */
655 kernel->width = kernel->height = 2*(long)args->rho+1;
656 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
659 /* NOTE: user defaults set in "AcquireKernelFromString()" */
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;
670 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
671 kernel->height*sizeof(double));
672 if (kernel->values == (double *) NULL)
673 return(DestroyKernel(kernel));
675 u=kernel->width*kernel->height;
676 for ( i=0; i < (unsigned long)u; i++)
677 kernel->values[i] = 1.0;
679 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
680 kernel->range_pos = (double) u;
685 kernel->width = kernel->height = 3; /* default radius = 1 */
687 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
688 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
690 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
691 kernel->height*sizeof(double));
692 if (kernel->values == (double *) NULL)
693 return(DestroyKernel(kernel));
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;
700 kernel->values[i] = nan;
701 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
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;
713 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
714 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
716 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
717 kernel->height*sizeof(double));
718 if (kernel->values == (double *) NULL)
719 return(DestroyKernel(kernel));
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;
726 kernel->values[i] = nan;
727 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
733 kernel->width = kernel->height = 5; /* default radius 2 */
735 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
736 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
738 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
739 kernel->height*sizeof(double));
740 if (kernel->values == (double *) NULL)
741 return(DestroyKernel(kernel));
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;
750 /* Distance Measuring Kernels */
751 case ChebyshevKernel:
757 kernel->width = kernel->height = 3; /* default radius = 1 */
759 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
760 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
762 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
763 kernel->height*sizeof(double));
764 if (kernel->values == (double *) NULL)
765 return(DestroyKernel(kernel));
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];
775 case ManhattenKernel:
781 kernel->width = kernel->height = 3; /* default radius = 1 */
783 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
784 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
786 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
787 kernel->height*sizeof(double));
788 if (kernel->values == (double *) NULL)
789 return(DestroyKernel(kernel));
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];
799 case EuclideanKernel:
805 kernel->width = kernel->height = 3; /* default radius = 1 */
807 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
808 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
810 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
811 kernel->height*sizeof(double));
812 if (kernel->values == (double *) NULL)
813 return(DestroyKernel(kernel));
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];
823 /* Undefined Kernels */
824 case LaplacianKernel:
827 assert("Kernel Type has not been defined yet");
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;
839 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
847 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
851 % D e s t r o y K e r n e l %
855 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
857 % DestroyKernel() frees the memory used by a Convolution/Morphology kernel.
859 % The format of the DestroyKernel method is:
861 % MagickKernel *DestroyKernel(MagickKernel *kernel)
863 % A description of each parameter follows:
865 % o kernel: the Morphology/Convolution kernel to be destroyed
869 MagickExport MagickKernel *DestroyKernel(MagickKernel *kernel)
871 assert(kernel != (MagickKernel *) NULL);
872 kernel->values=(double *)RelinquishMagickMemory(kernel->values);
873 kernel=(MagickKernel *) RelinquishMagickMemory(kernel);
878 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
882 % K e r n e l N o r m a l i z e %
886 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
888 % KernelNormalize() normalize the kernel so its convolution output will
889 % be over a unit range.
891 % The format of the KernelNormalize method is:
893 % void KernelRotate (MagickKernel *kernel)
895 % A description of each parameter follows:
897 % o kernel: the Morphology/Convolution kernel
900 MagickExport void KernelNormalize(MagickKernel *kernel)
902 register unsigned long
905 for (i=0; i < kernel->width; i++)
906 kernel->values[i] /= (kernel->range_pos - kernel->range_neg);
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);
917 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
921 % K e r n e l P r i n t %
925 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
927 % KernelPrint() Print out the kernel details to standard error
929 % The format of the KernelNormalize method is:
931 % void KernelPrint (MagickKernel *kernel)
933 % A description of each parameter follows:
935 % o kernel: the Morphology/Convolution kernel
938 MagickExport void KernelPrint(MagickKernel *kernel)
944 "Kernel \"%s\" of size %lux%lu%+ld%+ld with value from %lg to %lg\n",
945 MagickOptionToMnemonic(MagickKernelOptions, 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");
961 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
965 % K e r n e l R o t a t e %
969 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
974 % The format of the KernelRotate method is:
976 % void KernelRotate (MagickKernel *kernel, double angle)
978 % A description of each parameter follows:
980 % o kernel: the Morphology/Convolution kernel
982 % o angle: angle to rotate in degrees
985 MagickExport void KernelRotate(MagickKernel *kernel, double angle)
987 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
989 ** TODO: expand beyond simple 90 degree rotates, flips and flops
992 /* Modulus the angle */
993 angle = fmod(angle, 360.0);
997 if ( 315.0 < angle || angle <= 45.0 )
998 return; /* no change! - At least at this time */
1000 switch (kernel->type) {
1001 /* These built-in kernels are cylindrical kernel, rotating is useless */
1002 case GaussianKernel:
1003 case LaplacianKernel:
1007 case ChebyshevKernel:
1008 case ManhattenKernel:
1009 case EuclideanKernel:
1012 /* These may be rotatable at non-90 angles in the future */
1013 /* but simply rotating them 90 degrees is useless */
1019 /* These only allows a +/-90 degree rotation (transpose) */
1021 case RectangleKernel:
1022 if ( 135.0 < angle && angle <= 225.0 )
1024 if ( 225.0 < angle && angle <= 315.0 )
1028 /* these are freely rotatable in 90 degree units */
1030 case UndefinedKernel:
1031 case UserDefinedKernel:
1035 if ( 135.0 < angle && angle <= 315.0 )
1037 /* Do a flop, this assumes kernel is horizontally symetrical. */
1038 /* Each kernel data row need to be reversed! */
1041 register unsigned long
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;
1049 kernel->offset_x = kernel->width - kernel->offset_x - 1;
1050 angle = fmod(angle+180.0, 360.0);
1052 if ( 45.0 < angle && angle <= 135.0 )
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. */
1059 kernel->width = kernel->height;
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);
1066 /* at this point angle should be between +45 and -45 (315) degrees */
1071 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1075 % M o r p h o l o g y I m a g e %
1079 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1081 % MorphologyImage() applies a user supplied kernel to the image according to
1082 % the given mophology method.
1084 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1085 % by the kernel generator.
1087 % The format of the MorphologyImage method is:
1089 % Image *MorphologyImage(const Image *image, const MorphologyMethod
1090 % method, const long iterations, const ChannelType channel,
1091 % const MagickKernel *kernel, ExceptionInfo *exception)
1093 % A description of each parameter follows:
1095 % o image: the image.
1097 % o method: the morphology method to be applied.
1099 % o iterations: apply the operation this many times (or no change).
1100 % A value of -1 means loop until no change found.
1101 % How this is applied may depend on the morphology method.
1102 % Typically this is a value of 1.
1104 % o channel: the channel type.
1106 % o kernel: An array of double representing the morphology kernel.
1107 % Warning: kernel may be normalized for a Convolve.
1109 % o exception: return any errors or warnings in this structure.
1112 % TODO: bias and auto-scale handling of the kernel for convolution
1113 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1114 % by the kernel generator.
1118 static inline double MagickMin(const MagickRealType x,const MagickRealType y)
1120 return( x < y ? x : y);
1122 static inline double MagickMax(const MagickRealType x,const MagickRealType y)
1124 return( x > y ? x : y);
1126 #define Minimize(assign,value) assign=MagickMin(assign,value)
1127 #define Maximize(assign,value) assign=MagickMax(assign,value)
1129 /* incr change if the value being assigned changed */
1130 #define Assign(channel,value) \
1131 { q->channel = ClampToQuantum(value); \
1132 if ( p[r].channel != q->channel ) changed++; \
1134 #define AssignIndex(value) \
1135 { q_indexes[x] = ClampToQuantum(value); \
1136 if ( p_indexes[r] != q_indexes[x] ) changed++; \
1139 /* Internal function
1140 * Apply the Morphology method with the given Kernel
1141 * And return the number of values changed.
1143 static unsigned long MorphologyApply(const Image *image, Image
1144 *result_image, const MorphologyMethod method, const ChannelType channel,
1145 const MagickKernel *kernel, ExceptionInfo *exception)
1147 #define MorphologyTag "Morphology/Image"
1167 Apply Morphology to Image.
1173 GetMagickPixelPacket(image,&bias);
1174 SetMagickPixelPacketBias(image,&bias);
1176 p_view=AcquireCacheView(image);
1177 q_view=AcquireCacheView(result_image);
1178 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1179 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1181 for (y=0; y < (long) image->rows; y++)
1186 register const PixelPacket
1189 register const IndexPacket
1190 *restrict p_indexes;
1192 register PixelPacket
1195 register IndexPacket
1196 *restrict q_indexes;
1204 if (status == MagickFalse)
1206 p=GetCacheViewVirtualPixels(p_view, -kernel->offset_x, y-kernel->offset_y,
1207 image->columns+kernel->width, kernel->height, exception);
1208 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1210 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1215 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1216 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1217 r = (image->columns+kernel->width)*kernel->offset_y+kernel->offset_x;
1218 for (x=0; x < (long) image->columns; x++)
1226 register const double
1229 register const PixelPacket
1232 register const IndexPacket
1233 *restrict k_indexes;
1238 /* Copy input to ouput image - removes need for 'cloning' new images */
1240 if (image->colorspace == CMYKColorspace)
1241 q_indexes[x] = p_indexes[r];
1245 case ConvolveMorphology:
1247 break; /* default result is the convolution bias */
1248 case DialateIntensityMorphology:
1249 case ErodeIntensityMorphology:
1250 /* result is the pixel as is */
1251 result.red = p[r].red;
1252 result.green = p[r].green;
1253 result.blue = p[r].blue;
1254 result.opacity = p[r].opacity;
1255 if ( image->colorspace == CMYKColorspace)
1256 result.index = p_indexes[r];
1259 /* most need to handle transparency as alpha */
1260 result.red = p[r].red;
1261 result.green = p[r].green;
1262 result.blue = p[r].blue;
1263 result.opacity = QuantumRange - p[r].opacity;
1264 if ( image->colorspace == CMYKColorspace)
1265 result.index = p_indexes[r];
1270 case ConvolveMorphology:
1271 /* Weighted Average of pixels */
1272 if (((channel & OpacityChannel) == 0) ||
1273 (image->matte == MagickFalse))
1275 /* Kernel Weighted Convolution (no transparency) */
1278 k_indexes = p_indexes;
1279 for (v=0; v < (long) kernel->height; v++) {
1280 for (u=0; u < (long) kernel->width; u++, k++) {
1281 if ( IsNan(*k) ) continue;
1282 result.red += (*k)*k_pixels[u].red;
1283 result.green += (*k)*k_pixels[u].green;
1284 result.blue += (*k)*k_pixels[u].blue;
1285 /* result.opacity += no involvment */
1286 if ( image->colorspace == CMYKColorspace)
1287 result.index += (*k)*k_indexes[u];
1289 k_pixels += image->columns+kernel->width;
1290 k_indexes += image->columns+kernel->width;
1292 if ((channel & RedChannel) != 0)
1293 Assign(red,result.red);
1294 if ((channel & GreenChannel) != 0)
1295 Assign(green,result.green);
1296 if ((channel & BlueChannel) != 0)
1297 Assign(blue,result.blue);
1298 /* no transparency involved */
1299 if ((channel & IndexChannel) != 0
1300 && image->colorspace == CMYKColorspace)
1301 AssignIndex(result.index);
1304 { /* Kernel & Alpha weighted Convolution */
1306 alpha, /* alpha value * kernel weighting */
1307 gamma; /* weighting divisor */
1312 k_indexes = p_indexes;
1313 for (v=0; v < (long) kernel->height; v++) {
1314 for (u=0; u < (long) kernel->width; u++, k++) {
1315 if ( IsNan(*k) ) continue;
1316 alpha=(*k)*(QuantumScale*(QuantumRange-
1317 k_pixels[u].opacity));
1319 result.red += alpha*k_pixels[u].red;
1320 result.green += alpha*k_pixels[u].green;
1321 result.blue += alpha*k_pixels[u].blue;
1322 result.opacity += (*k)*k_pixels[u].opacity;
1323 if ( image->colorspace == CMYKColorspace)
1324 result.index += alpha*k_indexes[u];
1326 k_pixels += image->columns+kernel->width;
1327 k_indexes += image->columns+kernel->width;
1329 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1330 if ((channel & RedChannel) != 0)
1331 Assign(red,gamma*result.red);
1332 if ((channel & GreenChannel) != 0)
1333 Assign(green,gamma*result.green);
1334 if ((channel & BlueChannel) != 0)
1335 Assign(blue,gamma*result.blue);
1336 if ((channel & OpacityChannel) != 0
1337 && image->matte == MagickTrue )
1338 Assign(opacity,result.opacity);
1339 if ((channel & IndexChannel) != 0
1340 && image->colorspace == CMYKColorspace)
1341 AssignIndex(gamma*result.index);
1345 case DialateMorphology:
1346 /* Maximize Value - Kernel should be boolean */
1349 k_indexes = p_indexes;
1350 for (v=0; v < (long) kernel->height; v++) {
1351 for (u=0; u < (long) kernel->width; u++, k++) {
1352 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1353 Maximize(result.red, k_pixels[u].red);
1354 Maximize(result.green, k_pixels[u].green);
1355 Maximize(result.blue, k_pixels[u].blue);
1356 Maximize(result.opacity, QuantumRange-k_pixels[u].opacity);
1357 if ( image->colorspace == CMYKColorspace)
1358 Maximize(result.index, k_indexes[u]);
1360 k_pixels += image->columns+kernel->width;
1361 k_indexes += image->columns+kernel->width;
1363 if ((channel & RedChannel) != 0)
1364 Assign(red,result.red);
1365 if ((channel & GreenChannel) != 0)
1366 Assign(green,result.green);
1367 if ((channel & BlueChannel) != 0)
1368 Assign(blue,result.blue);
1369 if ((channel & OpacityChannel) != 0
1370 && image->matte == MagickTrue )
1371 Assign(opacity,QuantumRange-result.opacity);
1372 if ((channel & IndexChannel) != 0
1373 && image->colorspace == CMYKColorspace)
1374 AssignIndex(result.index);
1377 case ErodeMorphology:
1378 /* Minimize Value - Kernel should be boolean */
1381 k_indexes = p_indexes;
1382 for (v=0; v < (long) kernel->height; v++) {
1383 for (u=0; u < (long) kernel->width; u++, k++) {
1384 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1385 Minimize(result.red, k_pixels[u].red);
1386 Minimize(result.green, k_pixels[u].green);
1387 Minimize(result.blue, k_pixels[u].blue);
1388 Minimize(result.opacity, QuantumRange-k_pixels[u].opacity);
1389 if ( image->colorspace == CMYKColorspace)
1390 Minimize(result.index, k_indexes[u]);
1392 k_pixels += image->columns+kernel->width;
1393 k_indexes += image->columns+kernel->width;
1395 if ((channel & RedChannel) != 0)
1396 Assign(red,result.red);
1397 if ((channel & GreenChannel) != 0)
1398 Assign(green,result.green);
1399 if ((channel & BlueChannel) != 0)
1400 Assign(blue,result.blue);
1401 if ((channel & OpacityChannel) != 0
1402 && image->matte == MagickTrue )
1403 Assign(opacity,QuantumRange-result.opacity);
1404 if ((channel & IndexChannel) != 0
1405 && image->colorspace == CMYKColorspace)
1406 AssignIndex(result.index);
1409 case DialateIntensityMorphology:
1410 /* Maximum Intensity Pixel - Kernel should be boolean */
1413 k_indexes = p_indexes;
1414 for (v=0; v < (long) kernel->height; v++) {
1415 for (u=0; u < (long) kernel->width; u++, k++) {
1416 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1417 if ( PixelIntensity(&p[r]) >
1418 PixelIntensity(&(k_pixels[u])) ) continue;
1419 result.red = k_pixels[u].red;
1420 result.green = k_pixels[u].green;
1421 result.blue = k_pixels[u].blue;
1422 result.opacity = k_pixels[u].opacity;
1423 if ( image->colorspace == CMYKColorspace)
1424 result.index = k_indexes[u];
1426 k_pixels += image->columns+kernel->width;
1427 k_indexes += image->columns+kernel->width;
1429 if ((channel & RedChannel) != 0)
1430 Assign(red,result.red);
1431 if ((channel & GreenChannel) != 0)
1432 Assign(green,result.green);
1433 if ((channel & BlueChannel) != 0)
1434 Assign(blue,result.blue);
1435 if ((channel & OpacityChannel) != 0
1436 && image->matte == MagickTrue )
1437 Assign(opacity,result.opacity);
1438 if ((channel & IndexChannel) != 0
1439 && image->colorspace == CMYKColorspace)
1440 AssignIndex(result.index);
1443 case ErodeIntensityMorphology:
1444 /* Minimum Intensity Pixel - Kernel should be boolean */
1447 k_indexes = p_indexes;
1448 for (v=0; v < (long) kernel->height; v++) {
1449 for (u=0; u < (long) kernel->width; u++, k++) {
1450 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1451 if ( PixelIntensity(&p[r]) <
1452 PixelIntensity(&(k_pixels[u])) ) continue;
1453 result.red = k_pixels[u].red;
1454 result.green = k_pixels[u].green;
1455 result.blue = k_pixels[u].blue;
1456 result.opacity = k_pixels[u].opacity;
1457 if ( image->colorspace == CMYKColorspace)
1458 result.index = k_indexes[u];
1460 k_pixels += image->columns+kernel->width;
1461 k_indexes += image->columns+kernel->width;
1463 if ((channel & RedChannel) != 0)
1464 Assign(red,result.red);
1465 if ((channel & GreenChannel) != 0)
1466 Assign(green,result.green);
1467 if ((channel & BlueChannel) != 0)
1468 Assign(blue,result.blue);
1469 if ((channel & OpacityChannel) != 0
1470 && image->matte == MagickTrue )
1471 Assign(opacity,result.opacity);
1472 if ((channel & IndexChannel) != 0
1473 && image->colorspace == CMYKColorspace)
1474 AssignIndex(result.index);
1477 case DistanceMorphology:
1479 /* No need to do distance morphology if all values are zero */
1480 /* Unfortunatally I have not been able to get this right! */
1481 if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1482 || ((channel & GreenChannel) == 0 && p[r].green == 0)
1483 || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1484 || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1485 || (( (channel & IndexChannel) == 0
1486 || image->colorspace != CMYKColorspace
1487 ) && p_indexes[x] ==0 )
1493 k_indexes = p_indexes;
1494 for (v=0; v < (long) kernel->height; v++) {
1495 for (u=0; u < (long) kernel->width; u++, k++) {
1496 if ( IsNan(*k) ) continue;
1497 Minimize(result.red, (*k)+k_pixels[u].red);
1498 Minimize(result.green, (*k)+k_pixels[u].green);
1499 Minimize(result.blue, (*k)+k_pixels[u].blue);
1500 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1501 if ( image->colorspace == CMYKColorspace)
1502 Minimize(result.index, (*k)+k_indexes[u]);
1504 k_pixels += image->columns+kernel->width;
1505 k_indexes += image->columns+kernel->width;
1508 if ((channel & RedChannel) != 0)
1509 Assign(red,result.red);
1510 if ((channel & GreenChannel) != 0)
1511 Assign(green,result.green);
1512 if ((channel & BlueChannel) != 0)
1513 Assign(blue,result.blue);
1514 if ((channel & OpacityChannel) != 0
1515 && image->matte == MagickTrue )
1516 Assign(opacity,QuantumRange-result.opacity);
1517 if ((channel & IndexChannel) != 0
1518 && image->colorspace == CMYKColorspace)
1519 AssignIndex(result.index);
1521 /* By returning the number of 'maximum' values still to process
1522 ** we can get the Distance iteration to finish faster.
1523 ** BUT this may cause an infinite loop on very large shapes,
1524 ** which may have a distance that reachs a maximum gradient.
1526 if ((channel & RedChannel) != 0)
1527 { q->red = ClampToQuantum(result.red);
1528 if ( q->red == QuantumRange ) changed++; /* more to do */
1530 if ((channel & GreenChannel) != 0)
1531 { q->green = ClampToQuantum(result.green);
1532 if ( q->green == QuantumRange ) changed++; /* more to do */
1534 if ((channel & BlueChannel) != 0)
1535 { q->blue = ClampToQuantum(result.blue);
1536 if ( q->blue == QuantumRange ) changed++; /* more to do */
1538 if ((channel & OpacityChannel) != 0)
1539 { q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1540 if ( q->opacity == 0 ) changed++; /* more to do */
1542 if (((channel & IndexChannel) != 0) &&
1543 (image->colorspace == CMYKColorspace))
1544 { q_indexes[x] = ClampToQuantum(result.index);
1545 if ( q_indexes[x] == QuantumRange ) changed++;
1550 case UndefinedMorphology:
1552 break; /* Do nothing */
1557 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1558 if (sync == MagickFalse)
1560 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1565 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1566 #pragma omp critical (MagickCore_MorphologyImage)
1568 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1569 if (proceed == MagickFalse)
1573 result_image->type=image->type;
1574 q_view=DestroyCacheView(q_view);
1575 p_view=DestroyCacheView(p_view);
1576 return(status ? changed : 0);
1579 MagickExport Image *MorphologyImage(const Image *image,
1580 const ChannelType channel, MorphologyMethod method, const long iterations,
1581 MagickKernel *kernel, ExceptionInfo *exception)
1592 assert(image != (Image *) NULL);
1593 assert(image->signature == MagickSignature);
1594 assert(exception != (ExceptionInfo *) NULL);
1595 assert(exception->signature == MagickSignature);
1597 if ( GetImageArtifact(image,"showkernel") != (const char *) NULL)
1598 KernelPrint(kernel);
1600 if ( iterations == 0 )
1601 return((Image *)NULL); /* null operation - nothing to do! */
1603 /* kernel must be valid at this point
1604 * (except maybe for posible future morphology methods like "Prune"
1606 assert(kernel != (MagickKernel *)NULL);
1610 if ( iterations < 0 )
1611 limit = image->columns > image->rows ? image->columns : image->rows;
1613 /* Special morphology cases */
1614 changed=MagickFalse;
1616 case CloseMorphology:
1617 new_image = MorphologyImage(image, DialateMorphology, iterations, channel,
1619 if (new_image == (Image *) NULL)
1620 return((Image *) NULL);
1621 method = ErodeMorphology;
1623 case OpenMorphology:
1624 new_image = MorphologyImage(image, ErodeMorphology, iterations, channel,
1626 if (new_image == (Image *) NULL)
1627 return((Image *) NULL);
1628 method = DialateMorphology;
1630 case CloseIntensityMorphology:
1631 new_image = MorphologyImage(image, DialateIntensityMorphology,
1632 iterations, channel, kernel, exception);
1633 if (new_image == (Image *) NULL)
1634 return((Image *) NULL);
1635 method = ErodeIntensityMorphology;
1637 case OpenIntensityMorphology:
1638 new_image = MorphologyImage(image, ErodeIntensityMorphology,
1639 iterations, channel, kernel, exception);
1640 if (new_image == (Image *) NULL)
1641 return((Image *) NULL);
1642 method = DialateIntensityMorphology;
1645 case ConvolveMorphology:
1646 KernelNormalize(kernel);
1649 /* Do a morphology just once at this point!
1650 This ensures a new_image has been generated, but allows us
1651 to skip the creation of 'old_image' if it isn't needed.
1653 new_image=CloneImage(image,0,0,MagickTrue,exception);
1654 if (new_image == (Image *) NULL)
1655 return((Image *) NULL);
1656 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1658 InheritException(exception,&new_image->exception);
1659 new_image=DestroyImage(new_image);
1660 return((Image *) NULL);
1662 changed = MorphologyApply(image,new_image,method,channel,kernel,
1665 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1666 fprintf(stderr, "Morphology %s:%lu => Changed %lu\n",
1667 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1671 /* Repeat the interative morphology until count or no change */
1672 if ( count < limit && changed > 0 ) {
1673 old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1674 if (old_image == (Image *) NULL)
1675 return(DestroyImage(new_image));
1676 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1678 InheritException(exception,&old_image->exception);
1679 old_image=DestroyImage(old_image);
1680 return(DestroyImage(new_image));
1682 while( count < limit && changed != 0 )
1684 Image *tmp = old_image;
1685 old_image = new_image;
1687 changed = MorphologyApply(old_image,new_image,method,channel,kernel,
1690 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1691 fprintf(stderr, "Morphology %s:%lu => Changed %lu\n",
1692 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1695 DestroyImage(old_image);