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/magick.h"
65 #include "magick/memory_.h"
66 #include "magick/monitor-private.h"
67 #include "magick/morphology.h"
68 #include "magick/option.h"
69 #include "magick/pixel-private.h"
70 #include "magick/prepress.h"
71 #include "magick/quantize.h"
72 #include "magick/registry.h"
73 #include "magick/semaphore.h"
74 #include "magick/splay-tree.h"
75 #include "magick/statistic.h"
76 #include "magick/string_.h"
77 #include "magick/string-private.h"
78 #include "magick/token.h"
82 * The following test is for special floating point numbers of value NaN (not
83 * a number), that may be used within a Kernel Definition. NaN's are defined
84 * as part of the IEEE standard for floating point number representation.
86 * These are used a Kernel value of NaN means that that kernal position is not
87 * part of the normal convolution or morphology process, and thus allowing the
88 * use of 'shaped' kernels.
90 * Special Properities Two NaN's are never equal, even if they are from the
91 * same variable That is the IsNaN() macro is only true if the value is NaN.
93 #define IsNan(a) ((a)!=(a))
96 * Other global definitions used by module
98 static inline double MagickMin(const double x,const double y)
100 return( x < y ? x : y);
102 static inline double MagickMax(const double x,const double y)
104 return( x > y ? x : y);
106 #define Minimize(assign,value) assign=MagickMin(assign,value)
107 #define Maximize(assign,value) assign=MagickMax(assign,value)
111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
115 % A c q u i r e K e r n e l F r o m S t r i n g %
119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 % AcquireKernelInfo() takes the given string (generally supplied by the
122 % user) and converts it into a Morphology/Convolution Kernel. This allows
123 % users to specify a kernel from a number of pre-defined kernels, or to fully
124 % specify their own kernel for a specific Convolution or Morphology
127 % The kernel so generated can be any rectangular array of floating point
128 % values (doubles) with the 'control point' or 'pixel being affected'
129 % anywhere within that array of values.
131 % ASIDE: Previously IM was restricted to a square of odd size using the exact
134 % The floating point values in the kernel can also include a special value
135 % known as 'NaN' or 'not a number' to indicate that this value is not part
136 % of the kernel array. This allows you to specify a non-rectangular shaped
137 % kernel, for use in Morphological operators, without the need for some type
140 % The returned kernel should be freed using the DestroyKernel() when you are
143 % Input kernel defintion strings can consist of any of three types.
146 % Select from one of the built in kernels, using the name and
147 % geometry arguments supplied. See AcquireKernelBuiltIn()
149 % "WxH[+X+Y]:num, num, num ..."
150 % a kernal of size W by H, with W*H floating point numbers following.
151 % the 'center' can be optionally be defined at +X+Y (such that +0+0
152 % is top left corner). If not defined the pixel in the center, for
153 % odd sizes, or to the immediate top or left of center for even sizes
154 % is automatically selected.
156 % "num, num, num, num, ..."
157 % list of floating point numbers defining an 'old style' odd sized
158 % square kernel. At least 9 values should be provided for a 3x3
159 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
160 % Values can be space or comma separated. This is not recommended.
162 % Note that 'name' kernels will start with an alphabetic character
163 % while the new kernel specification has a ':' character in its
164 % specification string. If neither is the case, it assumes it is the
165 % old style of a simple list of numbers.
167 % The format of the AcquireKernal method is:
169 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
171 % A description of each parameter follows:
173 % o kernel_string: the Morphology/Convolution kernel wanted.
177 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
183 token[MaxTextExtent];
185 register unsigned long
198 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
200 assert(kernel_string != (const char *) NULL);
201 SetGeometryInfo(&args);
203 /* does it start with an alpha - Return a builtin kernel */
204 GetMagickToken(kernel_string,&p,token);
205 if ( isalpha((int)token[0]) )
210 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
211 if ( type < 0 || type == UserDefinedKernel )
212 return((KernelInfo *)NULL);
214 while (((isspace((int) ((unsigned char) *p)) != 0) ||
215 (*p == ',') || (*p == ':' )) && (*p != '\0'))
217 flags = ParseGeometry(p, &args);
219 /* special handling of missing values in input string */
220 if ( type == RectangleKernel ) {
221 if ( (flags & WidthValue) == 0 ) /* if no width then */
222 args.rho = args.sigma; /* then width = height */
223 if ( args.rho < 1.0 ) /* if width too small */
224 args.rho = 3; /* then width = 3 */
225 if ( args.sigma < 1.0 ) /* if height too small */
226 args.sigma = args.rho; /* then height = width */
227 if ( (flags & XValue) == 0 ) /* center offset if not defined */
228 args.xi = (double)(((long)args.rho-1)/2);
229 if ( (flags & YValue) == 0 )
230 args.psi = (double)(((long)args.sigma-1)/2);
233 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
236 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
237 if (kernel == (KernelInfo *)NULL)
239 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
240 kernel->type = UserDefinedKernel;
241 kernel->signature = MagickSignature;
243 /* Has a ':' in argument - New user kernel specification */
244 p = strchr(kernel_string, ':');
245 if ( p != (char *) NULL)
247 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
248 memcpy(token, kernel_string, p-kernel_string);
249 token[p-kernel_string] = '\0';
250 flags = ParseGeometry(token, &args);
252 /* Size handling and checks of geometry settings */
253 if ( (flags & WidthValue) == 0 ) /* if no width then */
254 args.rho = args.sigma; /* then width = height */
255 if ( args.rho < 1.0 ) /* if width too small */
256 args.rho = 1.0; /* then width = 1 */
257 if ( args.sigma < 1.0 ) /* if height too small */
258 args.sigma = args.rho; /* then height = width */
259 kernel->width = (unsigned long)args.rho;
260 kernel->height = (unsigned long)args.sigma;
262 /* Offset Handling and Checks */
263 if ( args.xi < 0.0 || args.psi < 0.0 )
264 return(DestroyKernel(kernel));
265 kernel->offset_x = ((flags & XValue)!=0) ? (unsigned long)args.xi
266 : (kernel->width-1)/2;
267 kernel->offset_y = ((flags & YValue)!=0) ? (unsigned long)args.psi
268 : (kernel->height-1)/2;
269 if ( kernel->offset_x >= kernel->width ||
270 kernel->offset_y >= kernel->height )
271 return(DestroyKernel(kernel));
273 p++; /* advance beyond the ':' */
276 { /* ELSE - Old old kernel specification, forming odd-square kernel */
277 /* count up number of values given */
278 p=(const char *) kernel_string;
279 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
280 p++; /* ignore "'" chars for convolve filter usage - Cristy */
281 for (i=0; *p != '\0'; i++)
283 GetMagickToken(p,&p,token);
285 GetMagickToken(p,&p,token);
287 /* set the size of the kernel - old sized square */
288 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
289 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
290 p=(const char *) kernel_string;
291 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
292 p++; /* ignore "'" chars for convolve filter usage - Cristy */
295 /* Read in the kernel values from rest of input string argument */
296 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
297 kernel->height*sizeof(double));
298 if (kernel->values == (double *) NULL)
299 return(DestroyKernel(kernel));
301 kernel->value_min = +MagickHuge;
302 kernel->value_max = -MagickHuge;
303 kernel->range_neg = kernel->range_pos = 0.0;
304 for (i=0; (i < kernel->width*kernel->height) && (*p != '\0'); i++)
306 GetMagickToken(p,&p,token);
308 GetMagickToken(p,&p,token);
309 if ( LocaleCompare("nan",token) == 0
310 || LocaleCompare("-",token) == 0 ) {
311 kernel->values[i] = nan; /* do not include this value in kernel */
314 kernel->values[i] = StringToDouble(token);
315 ( kernel->values[i] < 0)
316 ? ( kernel->range_neg += kernel->values[i] )
317 : ( kernel->range_pos += kernel->values[i] );
318 Minimize(kernel->value_min, kernel->values[i]);
319 Maximize(kernel->value_max, kernel->values[i]);
323 /* This should not be needed for a fully defined defined kernel
324 * Perhaps an error should be reported instead!
326 if ( i < kernel->width*kernel->height ) {
327 Minimize(kernel->value_min, kernel->values[i]);
328 Maximize(kernel->value_max, kernel->values[i]);
329 for ( ; i < kernel->width*kernel->height; i++)
330 kernel->values[i]=0.0;
337 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
341 % A c q u i r e K e r n e l B u i l t I n %
345 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
347 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
348 % kernels used for special purposes such as gaussian blurring, skeleton
349 % pruning, and edge distance determination.
351 % They take a KernelType, and a set of geometry style arguments, which were
352 % typically decoded from a user supplied string, or from a more complex
353 % Morphology Method that was requested.
355 % The format of the AcquireKernalBuiltIn method is:
357 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
358 % const GeometryInfo args)
360 % A description of each parameter follows:
362 % o type: the pre-defined type of kernel wanted
364 % o args: arguments defining or modifying the kernel
366 % Convolution Kernels
368 % Gaussian "[{radius}]x{sigma}"
369 % Generate a two-dimentional gaussian kernel, as used by -gaussian
370 % A sigma is required, (with the 'x'), due to historical reasons.
372 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
373 % the final size of the resulting kernel to a square 2*radius+1 in size.
374 % The radius should be at least 2 times that of the sigma value, or
375 % sever clipping and aliasing may result. If not given or set to 0 the
376 % radius will be determined so as to produce the best minimal error
377 % result, which is usally much larger than is normally needed.
379 % Blur "[{radius}]x{sigma}[+angle]"
380 % As per Gaussian, but generates a 1 dimensional or linear gaussian
381 % blur, at the angle given (current restricted to orthogonal angles).
382 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
384 % NOTE that two such blurs perpendicular to each other is equivelent to
385 % -blur and the previous gaussian, but is often 10 or more times faster.
387 % Comet "[{width}]x{sigma}[+angle]"
388 % Blur in one direction only, mush like how a bright object leaves
389 % a comet like trail. The Kernel is actually half a gaussian curve,
390 % Adding two such blurs in oppiste directions produces a Linear Blur.
392 % NOTE: that the first argument is the width of the kernel and not the
393 % radius of the kernel.
395 % # Still to be implemented...
397 % # Laplacian "{radius}x{sigma}"
398 % # Laplacian (a mexican hat like) Function
400 % # LOG "{radius},{sigma1},{sigma2}
401 % # Laplacian of Gaussian
403 % # DOG "{radius},{sigma1},{sigma2}
404 % # Difference of Gaussians
408 % Rectangle "{geometry}"
409 % Simply generate a rectangle of 1's with the size given. You can also
410 % specify the location of the 'control point', otherwise the closest
411 % pixel to the center of the rectangle is selected.
413 % Properly centered and odd sized rectangles work the best.
415 % Diamond "[{radius}]"
416 % Generate a diamond shaped kernal with given radius to the points.
417 % Kernel size will again be radius*2+1 square and defaults to radius 1,
418 % generating a 3x3 kernel that is slightly larger than a square.
420 % Square "[{radius}]"
421 % Generate a square shaped kernel of size radius*2+1, and defaulting
422 % to a 3x3 (radius 1).
424 % Note that using a larger radius for the "Square" or the "Diamond"
425 % is also equivelent to iterating the basic morphological method
426 % that many times. However However iterating with the smaller radius 1
427 % default is actually faster than using a larger kernel radius.
430 % Generate a binary disk of the radius given, radius may be a float.
431 % Kernel size will be ceil(radius)*2+1 square.
432 % NOTE: Here are some disk shapes of specific interest
433 % "disk:1" => "diamond" or "cross:1"
434 % "disk:1.5" => "square"
435 % "disk:2" => "diamond:2"
436 % "disk:2.5" => default - radius 2 disk shape
437 % "disk:2.9" => "square:2"
438 % "disk:3.5" => octagonal/disk shape of radius 3
439 % "disk:4.2" => roughly octagonal shape of radius 4
440 % "disk:4.3" => disk shape of radius 4
441 % After this all the kernel shape becomes more and more circular.
443 % Because a "disk" is more circular when using a larger radius, using a
444 % larger radius is preferred over iterating the morphological operation.
447 % Generate a kernel in the shape of a 'plus' sign. The length of each
448 % arm is also the radius, which defaults to 2.
450 % This kernel is not a good general morphological kernel, but is used
451 % more for highlighting and marking any single pixels in an image using,
452 % a "Dilate" or "Erode" method as appropriate.
454 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
456 % Note that unlike other kernels iterating a plus does not produce the
457 % same result as using a larger radius for the cross.
459 % Distance Measuring Kernels
461 % Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1)
462 % Manhatten "[{radius}][x{scale}]" square grid distance (default r=1)
463 % Euclidean "[{radius}][x{scale}]" direct distance (default r=1)
465 % Different types of distance measuring methods, which are used with the
466 % a 'Distance' morphology method for generating a gradient based on
467 % distance from an edge of a binary shape, though there is a technique
468 % for handling a anti-aliased shape.
470 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
471 % one to any neighbour, orthogonal or diagonal. One why of thinking of
472 % it is the number of squares a 'King' or 'Queen' in chess needs to
473 % traverse reach any other position on a chess board. It results in a
474 % 'square' like distance function, but one where diagonals are closer
477 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
478 % Cab metric), is the distance needed when you can only travel in
479 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
480 % in chess would travel. It results in a diamond like distances, where
481 % diagonals are further than expected.
483 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
484 % However by default the kernel size only has a radius of 1, which
485 % limits the distance to 'Knight' like moves, with only orthogonal and
486 % diagonal measurements being correct. As such for the default kernel
487 % you will get octagonal like distance function, which is reasonally
490 % However if you use a larger radius such as "Euclidean:4" you will
491 % get a much smoother distance gradient from the edge of the shape.
492 % Of course a larger kernel is slower to use, and generally not needed.
494 % To allow the use of fractional distances that you get with diagonals
495 % the actual distance is scaled by a fixed value which the user can
496 % provide. This is not actually nessary for either ""Chebyshev" or
497 % "Manhatten" distance kernels, but is done for all three distance
498 % kernels. If no scale is provided it is set to a value of 100,
499 % allowing for a maximum distance measurement of 655 pixels using a Q16
500 % version of IM, from any edge. However for small images this can
501 % result in quite a dark gradient.
503 % See the 'Distance' Morphological Method, for information of how it is
508 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
509 const GeometryInfo *args)
514 register unsigned long
522 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
524 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
525 if (kernel == (KernelInfo *) NULL)
527 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
528 kernel->value_min = kernel->value_max = 0.0;
529 kernel->range_neg = kernel->range_pos = 0.0;
531 kernel->signature = MagickSignature;
534 /* Convolution Kernels */
537 sigma = fabs(args->sigma);
539 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
541 kernel->width = kernel->height =
542 GetOptimalKernelWidth2D(args->rho,sigma);
543 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
544 kernel->range_neg = kernel->range_pos = 0.0;
545 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
546 kernel->height*sizeof(double));
547 if (kernel->values == (double *) NULL)
548 return(DestroyKernel(kernel));
550 sigma = 2.0*sigma*sigma; /* simplify the expression */
551 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
552 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
553 kernel->range_pos += (
555 exp(-((double)(u*u+v*v))/sigma)
556 /* / (MagickPI*sigma) */ );
557 kernel->value_min = 0;
558 kernel->value_max = kernel->values[
559 kernel->offset_y*kernel->width+kernel->offset_x ];
561 KernelNormalize(kernel);
567 sigma = fabs(args->sigma);
569 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
571 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
572 kernel->offset_x = (kernel->width-1)/2;
574 kernel->offset_y = 0;
575 kernel->range_neg = kernel->range_pos = 0.0;
576 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
577 kernel->height*sizeof(double));
578 if (kernel->values == (double *) NULL)
579 return(DestroyKernel(kernel));
583 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
584 ** It generates a gaussian 3 times the width, and compresses it into
585 ** the expected range. This produces a closer normalization of the
586 ** resulting kernel, especially for very low sigma values.
587 ** As such while wierd it is prefered.
589 ** I am told this method originally came from Photoshop.
591 sigma *= KernelRank; /* simplify expanded curve */
592 v = (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
593 (void) ResetMagickMemory(kernel->values,0, (size_t)
594 kernel->width*sizeof(double));
595 for ( u=-v; u <= v; u++) {
596 kernel->values[(u+v)/KernelRank] +=
597 exp(-((double)(u*u))/(2.0*sigma*sigma))
598 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
600 for (i=0; i < kernel->width; i++)
601 kernel->range_pos += kernel->values[i];
603 for ( i=0, u=-kernel->offset_x; i < kernel->width; i++, u++)
604 kernel->range_pos += (
606 exp(-((double)(u*u))/(2.0*sigma*sigma))
607 /* / (MagickSQ2PI*sigma) */ );
609 kernel->value_min = 0;
610 kernel->value_max = kernel->values[ kernel->offset_x ];
611 /* Note that both the above methods do not generate a normalized
612 ** kernel, though it gets close. The kernel may be 'clipped' by a user
613 ** defined radius, producing a smaller (darker) kernel. Also for very
614 ** small sigma's (> 0.1) the central value becomes larger than one,
615 ** and thus producing a very bright kernel.
618 /* Normalize the 1D Gaussian Kernel
620 ** Because of this the divisor in the above kernel generator is
621 ** not needed, so is not done above.
623 KernelNormalize(kernel);
625 /* rotate the kernel by given angle */
626 KernelRotate(kernel, args->xi);
631 sigma = fabs(args->sigma);
633 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
635 if ( args->rho < 1.0 )
636 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
638 kernel->width = (unsigned long)args->rho;
639 kernel->offset_x = kernel->offset_y = 0;
641 kernel->range_neg = kernel->range_pos = 0.0;
642 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
643 kernel->height*sizeof(double));
644 if (kernel->values == (double *) NULL)
645 return(DestroyKernel(kernel));
647 /* A comet blur is half a gaussian curve, so that the object is
648 ** blurred in one direction only. This may not be quite the right
649 ** curve so may change in the future. The function must be normalised.
653 sigma *= KernelRank; /* simplify expanded curve */
654 v = kernel->width*KernelRank; /* start/end points to fit range */
655 (void) ResetMagickMemory(kernel->values,0, (size_t)
656 kernel->width*sizeof(double));
657 for ( u=0; u < v; u++) {
658 kernel->values[u/KernelRank] +=
659 exp(-((double)(u*u))/(2.0*sigma*sigma))
660 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
662 for (i=0; i < kernel->width; i++)
663 kernel->range_pos += kernel->values[i];
665 for ( i=0; i < kernel->width; i++)
666 kernel->range_pos += (
668 exp(-((double)(i*i))/(2.0*sigma*sigma))
669 /* / (MagickSQ2PI*sigma) */ );
671 kernel->value_min = 0;
672 kernel->value_max = kernel->values[0];
674 KernelNormalize(kernel);
675 KernelRotate(kernel, args->xi);
678 /* Boolean Kernels */
679 case RectangleKernel:
682 if ( type == SquareKernel )
685 kernel->width = kernel->height = 3; /* default radius = 1 */
687 kernel->width = kernel->height = 2*(long)args->rho+1;
688 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
691 /* NOTE: user defaults set in "AcquireKernelInfo()" */
692 if ( args->rho < 1.0 || args->sigma < 1.0 )
693 return(DestroyKernel(kernel)); /* invalid args given */
694 kernel->width = (unsigned long)args->rho;
695 kernel->height = (unsigned long)args->sigma;
696 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
697 args->psi < 0.0 || args->psi > (double)kernel->height )
698 return(DestroyKernel(kernel)); /* invalid args given */
699 kernel->offset_x = (unsigned long)args->xi;
700 kernel->offset_y = (unsigned long)args->psi;
702 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
703 kernel->height*sizeof(double));
704 if (kernel->values == (double *) NULL)
705 return(DestroyKernel(kernel));
707 u=kernel->width*kernel->height;
708 for ( i=0; i < (unsigned long)u; i++)
709 kernel->values[i] = 1.0;
711 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
712 kernel->range_pos = (double) u;
717 kernel->width = kernel->height = 3; /* default radius = 1 */
719 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
720 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
722 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
723 kernel->height*sizeof(double));
724 if (kernel->values == (double *) NULL)
725 return(DestroyKernel(kernel));
727 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
728 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
729 if ((labs(u)+labs(v)) <= (long)kernel->offset_x)
730 kernel->range_pos += kernel->values[i] = 1.0;
732 kernel->values[i] = nan;
733 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
741 limit = (long)(args->rho*args->rho);
742 if (args->rho < 0.1) /* default radius approx 2.5 */
743 kernel->width = kernel->height = 5L, limit = 5L;
745 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
746 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
748 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
749 kernel->height*sizeof(double));
750 if (kernel->values == (double *) NULL)
751 return(DestroyKernel(kernel));
753 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
754 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
755 if ((u*u+v*v) <= limit)
756 kernel->range_pos += kernel->values[i] = 1.0;
758 kernel->values[i] = nan;
759 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
765 kernel->width = kernel->height = 5; /* default radius 2 */
767 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
768 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
770 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
771 kernel->height*sizeof(double));
772 if (kernel->values == (double *) NULL)
773 return(DestroyKernel(kernel));
775 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
776 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
777 kernel->values[i] = (u == 0 || v == 0) ? 1.0 : nan;
778 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
779 kernel->range_pos = kernel->width*2.0 - 1.0;
782 /* Distance Measuring Kernels */
783 case ChebyshevKernel:
789 kernel->width = kernel->height = 3; /* default radius = 1 */
791 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
792 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
794 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
795 kernel->height*sizeof(double));
796 if (kernel->values == (double *) NULL)
797 return(DestroyKernel(kernel));
799 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
800 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
801 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
802 kernel->range_pos += ( kernel->values[i] =
803 scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
804 kernel->value_max = kernel->values[0];
807 case ManhattenKernel:
813 kernel->width = kernel->height = 3; /* default radius = 1 */
815 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
816 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
818 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
819 kernel->height*sizeof(double));
820 if (kernel->values == (double *) NULL)
821 return(DestroyKernel(kernel));
823 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
824 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
825 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
826 kernel->range_pos += ( kernel->values[i] =
827 scale*(labs(u)+labs(v)) );
828 kernel->value_max = kernel->values[0];
831 case EuclideanKernel:
837 kernel->width = kernel->height = 3; /* default radius = 1 */
839 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
840 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
842 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
843 kernel->height*sizeof(double));
844 if (kernel->values == (double *) NULL)
845 return(DestroyKernel(kernel));
847 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
848 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
849 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
850 kernel->range_pos += ( kernel->values[i] =
851 scale*sqrt((double)(u*u+v*v)) );
852 kernel->value_max = kernel->values[0];
855 /* Undefined Kernels */
856 case LaplacianKernel:
859 assert("Kernel Type has not been defined yet");
862 /* Generate a No-Op minimal kernel - 1x1 pixel */
863 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
864 if (kernel->values == (double *) NULL)
865 return(DestroyKernel(kernel));
866 kernel->width = kernel->height = 1;
867 kernel->offset_x = kernel->offset_x = 0;
868 kernel->type = UndefinedKernel;
871 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
879 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
883 % D e s t r o y K e r n e l %
887 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
889 % DestroyKernel() frees the memory used by a Convolution/Morphology kernel.
891 % The format of the DestroyKernel method is:
893 % KernelInfo *DestroyKernel(KernelInfo *kernel)
895 % A description of each parameter follows:
897 % o kernel: the Morphology/Convolution kernel to be destroyed
901 MagickExport KernelInfo *DestroyKernel(KernelInfo *kernel)
903 assert(kernel != (KernelInfo *) NULL);
904 kernel->values=(double *)RelinquishMagickMemory(kernel->values);
905 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
910 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
914 % K e r n e l N o r m a l i z e %
918 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
920 % KernelNormalize() normalize the kernel so its convolution output will
921 % be over a unit range.
923 % The format of the KernelNormalize method is:
925 % void KernelRotate (KernelInfo *kernel)
927 % A description of each parameter follows:
929 % o kernel: the Morphology/Convolution kernel
932 MagickExport void KernelNormalize(KernelInfo *kernel)
934 register unsigned long
937 for (i=0; i < kernel->width*kernel->height; i++)
938 kernel->values[i] /= (kernel->range_pos - kernel->range_neg);
940 kernel->range_pos /= (kernel->range_pos - kernel->range_neg);
941 kernel->range_neg /= (kernel->range_pos - kernel->range_neg);
942 kernel->value_max /= (kernel->range_pos - kernel->range_neg);
943 kernel->value_min /= (kernel->range_pos - kernel->range_neg);
949 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
953 % K e r n e l P r i n t %
957 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
959 % KernelPrint() Print out the kernel details to standard error
961 % The format of the KernelNormalize method is:
963 % void KernelPrint (KernelInfo *kernel)
965 % A description of each parameter follows:
967 % o kernel: the Morphology/Convolution kernel
970 MagickExport void KernelPrint(KernelInfo *kernel)
976 "Kernel \"%s\" of size %lux%lu%+ld%+ld with value from %lg to %lg\n",
977 MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
978 kernel->width, kernel->height,
979 kernel->offset_x, kernel->offset_y,
980 kernel->value_min, kernel->value_max);
981 fprintf(stderr, "Forming convolution output range from %lg to %lg%s\n",
982 kernel->range_neg, kernel->range_pos,
983 kernel->normalized == MagickTrue ? " (normalized)" : "" );
984 for (i=v=0; v < kernel->height; v++) {
985 fprintf(stderr,"%2ld:",v);
986 for (u=0; u < kernel->width; u++, i++)
987 if ( IsNan(kernel->values[i]) )
988 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
990 fprintf(stderr," %.*lg", GetMagickPrecision(), kernel->values[i]);
991 fprintf(stderr,"\n");
996 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1000 % K e r n e l R o t a t e %
1004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1006 % KernelRotate() rotates the kernel by the angle given. Currently it is
1007 % restricted to 90 degree angles, but this may be improved in the future.
1009 % The format of the KernelRotate method is:
1011 % void KernelRotate (KernelInfo *kernel, double angle)
1013 % A description of each parameter follows:
1015 % o kernel: the Morphology/Convolution kernel
1017 % o angle: angle to rotate in degrees
1020 MagickExport void KernelRotate(KernelInfo *kernel, double angle)
1022 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
1024 ** TODO: expand beyond simple 90 degree rotates, flips and flops
1027 /* Modulus the angle */
1028 angle = fmod(angle, 360.0);
1032 if ( 315.0 < angle || angle <= 45.0 )
1033 return; /* no change! - At least at this time */
1035 switch (kernel->type) {
1036 /* These built-in kernels are cylindrical kernels, rotating is useless */
1037 case GaussianKernel:
1038 case LaplacianKernel:
1042 case ChebyshevKernel:
1043 case ManhattenKernel:
1044 case EuclideanKernel:
1047 /* These may be rotatable at non-90 angles in the future */
1048 /* but simply rotating them in multiples of 90 degrees is useless */
1054 /* These only allows a +/-90 degree rotation (by transpose) */
1055 /* A 180 degree rotation is useless */
1057 case RectangleKernel:
1058 if ( 135.0 < angle && angle <= 225.0 )
1060 if ( 225.0 < angle && angle <= 315.0 )
1064 /* these are freely rotatable in 90 degree units */
1066 case UndefinedKernel:
1067 case UserDefinedKernel:
1070 if ( 135.0 < angle && angle <= 225.0 )
1072 /* for a 180 degree rotation - also know as a reflection */
1073 /* This is actually a very very common operation! */
1074 /* basically this is a reverse of the kernel data! */
1081 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
1082 t=k[i], k[i]=k[j], k[j]=t;
1084 kernel->offset_x = kernel->width - kernel->offset_x - 1;
1085 kernel->offset_y = kernel->width - kernel->offset_y - 1;
1086 angle = fmod(angle+180.0, 360.0);
1088 if ( 45.0 < angle && angle <= 135.0 )
1089 { /* Do a transpose and a flop, of the image, which results in a 90
1090 * degree rotation using two mirror operations.
1092 * WARNING: this assumes the original image was a 1 dimentional image
1093 * but currently that is the only built-ins it is applied to.
1098 kernel->width = kernel->height;
1100 t = kernel->offset_x;
1101 kernel->offset_x = kernel->offset_y;
1102 kernel->offset_y = t;
1103 angle = fmod(450.0 - angle, 360.0);
1105 /* at this point angle should be between -45 (315) and +45 degrees */
1108 { /* Do a flop, this assumes kernel is horizontally symetrical.
1109 * Each kernel data row need to be reversed!
1110 * This is currently not used, but by be used in the future.
1114 register unsigned long
1119 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
1120 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
1121 t=k[x], k[x]=k[r], k[r]=t;
1123 kernel->offset_x = kernel->width - kernel->offset_x - 1;
1124 angle = fmod(angle+180.0, 360.0);
1131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1135 % M o r p h o l o g y I m a g e C h a n n e l %
1139 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1141 % MorphologyImageChannel() applies a user supplied kernel to the image
1142 % according to the given mophology method.
1144 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1145 % by the kernel generator.
1147 % The format of the MorphologyImage method is:
1149 % Image *MorphologyImage(const Image *image, MorphologyMethod method,
1150 % const long iterations, KernelInfo *kernel, ExceptionInfo *exception)
1151 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1152 % channel, MorphologyMethod method, const long iterations, KernelInfo
1153 % *kernel, ExceptionInfo *exception)
1155 % A description of each parameter follows:
1157 % o image: the image.
1159 % o method: the morphology method to be applied.
1161 % o iterations: apply the operation this many times (or no change).
1162 % A value of -1 means loop until no change found.
1163 % How this is applied may depend on the morphology method.
1164 % Typically this is a value of 1.
1166 % o channel: the channel type.
1168 % o kernel: An array of double representing the morphology kernel.
1169 % Warning: kernel may be normalized for the Convolve method.
1171 % o exception: return any errors or warnings in this structure.
1174 % TODO: bias and auto-scale handling of the kernel for convolution
1175 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1176 % by the kernel generator.
1180 /* incr change if the value being assigned changed */
1181 #define Assign(channel,value) \
1182 { q->channel = ClampToQuantum(value); \
1183 if ( p[r].channel != q->channel ) changed++; \
1185 #define AssignIndex(value) \
1186 { q_indexes[x] = ClampToQuantum(value); \
1187 if ( p_indexes[r] != q_indexes[x] ) changed++; \
1190 /* Internal function
1191 * Apply the Morphology method with the given Kernel
1192 * And return the number of values changed.
1194 static unsigned long MorphologyApply(const Image *image, Image
1195 *result_image, const MorphologyMethod method, const ChannelType channel,
1196 const KernelInfo *kernel, ExceptionInfo *exception)
1198 #define MorphologyTag "Morphology/Image"
1218 Apply Morphology to Image.
1224 GetMagickPixelPacket(image,&bias);
1225 SetMagickPixelPacketBias(image,&bias);
1227 p_view=AcquireCacheView(image);
1228 q_view=AcquireCacheView(result_image);
1230 /* some methods (including convolve) needs to apply a reflected kernel
1231 * the offset for getting the kernel view needs to be adjusted for this
1234 offx = kernel->offset_x;
1235 offy = kernel->offset_y;
1237 case ErodeMorphology:
1238 case ErodeIntensityMorphology:
1239 /* kernel is not reflected */
1242 /* kernel needs to be reflected */
1243 offx = kernel->width-offx-1;
1244 offy = kernel->height-offy-1;
1248 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1249 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1251 for (y=0; y < image->rows; y++)
1256 register const PixelPacket
1259 register const IndexPacket
1260 *restrict p_indexes;
1262 register PixelPacket
1265 register IndexPacket
1266 *restrict q_indexes;
1268 register unsigned long
1274 if (status == MagickFalse)
1276 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1277 image->columns+kernel->width, kernel->height, exception);
1278 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1280 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1285 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1286 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1287 r = (image->columns+kernel->width)*offy+offx; /* constant */
1289 for (x=0; x < image->columns; x++)
1294 register unsigned long
1297 register const double
1300 register const PixelPacket
1303 register const IndexPacket
1304 *restrict k_indexes;
1309 /* Copy input to ouput image for unused channels
1310 * This removes need for 'cloning' new images
1313 if (image->colorspace == CMYKColorspace)
1314 q_indexes[x] = p_indexes[r];
1316 result.index=0; /* stop compiler warnings */
1318 case ConvolveMorphology:
1320 break; /* default result is the convolution bias */
1322 removed as it caused problems with change
1323 Really we want to set each to the first value found
1324 but not up date change if that is the first value!
1325 case DialateMorphology:
1330 result.index = -MagickHuge;
1332 case ErodeMorphology:
1337 result.index = +MagickHuge;
1340 case DialateIntensityMorphology:
1341 case ErodeIntensityMorphology:
1342 result.red = 0.0; /* was initial pixel found ? */
1345 /* Otherwise just start with the original pixel value */
1346 result.red = p[r].red;
1347 result.green = p[r].green;
1348 result.blue = p[r].blue;
1349 result.opacity = QuantumRange - p[r].opacity;
1350 if ( image->colorspace == CMYKColorspace)
1351 result.index = p_indexes[r];
1356 case ConvolveMorphology:
1357 /* Weighted Average of pixels
1359 * NOTE for correct working of this operation for asymetrical
1360 * kernels, the kernel needs to be applied in its reflected form.
1361 * That is its values needs to be reversed.
1363 if (((channel & OpacityChannel) == 0) ||
1364 (image->matte == MagickFalse))
1366 /* Convolution (no transparency) */
1367 k = &kernel->values[ kernel->width*kernel->height-1 ];
1369 k_indexes = p_indexes;
1370 for (v=0; v < kernel->height; v++) {
1371 for (u=0; u < kernel->width; u++, k--) {
1372 if ( IsNan(*k) ) continue;
1373 result.red += (*k)*k_pixels[u].red;
1374 result.green += (*k)*k_pixels[u].green;
1375 result.blue += (*k)*k_pixels[u].blue;
1376 /* result.opacity += no involvment */
1377 if ( image->colorspace == CMYKColorspace)
1378 result.index += (*k)*k_indexes[u];
1380 k_pixels += image->columns+kernel->width;
1381 k_indexes += image->columns+kernel->width;
1383 if ((channel & RedChannel) != 0)
1384 Assign(red,result.red);
1385 if ((channel & GreenChannel) != 0)
1386 Assign(green,result.green);
1387 if ((channel & BlueChannel) != 0)
1388 Assign(blue,result.blue);
1389 /* no transparency involved */
1390 if ((channel & IndexChannel) != 0
1391 && image->colorspace == CMYKColorspace)
1392 AssignIndex(result.index);
1395 { /* Kernel & Alpha weighted Convolution */
1397 alpha, /* alpha value * kernel weighting */
1398 gamma; /* weighting divisor */
1401 k = &kernel->values[ kernel->width*kernel->height-1 ];
1403 k_indexes = p_indexes;
1404 for (v=0; v < kernel->height; v++) {
1405 for (u=0; u < kernel->width; u++, k--) {
1406 if ( IsNan(*k) ) continue;
1407 alpha=(*k)*(QuantumScale*(QuantumRange-
1408 k_pixels[u].opacity));
1410 result.red += alpha*k_pixels[u].red;
1411 result.green += alpha*k_pixels[u].green;
1412 result.blue += alpha*k_pixels[u].blue;
1413 result.opacity += (*k)*k_pixels[u].opacity;
1414 if ( image->colorspace == CMYKColorspace)
1415 result.index += alpha*k_indexes[u];
1417 k_pixels += image->columns+kernel->width;
1418 k_indexes += image->columns+kernel->width;
1420 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1421 if ((channel & RedChannel) != 0)
1422 Assign(red,gamma*result.red);
1423 if ((channel & GreenChannel) != 0)
1424 Assign(green,gamma*result.green);
1425 if ((channel & BlueChannel) != 0)
1426 Assign(blue,gamma*result.blue);
1427 if ((channel & OpacityChannel) != 0
1428 && image->matte == MagickTrue )
1429 Assign(opacity,result.opacity);
1430 if ((channel & IndexChannel) != 0
1431 && image->colorspace == CMYKColorspace)
1432 AssignIndex(gamma*result.index);
1436 case DialateMorphology:
1437 /* Maximize Value within kernel shape
1439 * NOTE for correct working of this operation for asymetrical
1440 * kernels, the kernel needs to be applied in its reflected form.
1441 * That is its values needs to be reversed.
1443 * NOTE: in normal Greyscale Morphology, the kernel value should
1444 * be added to the real value, this is currently not done, due to
1445 * the nature of the boolean kernels being used.
1448 k = &kernel->values[ kernel->width*kernel->height-1 ];
1450 k_indexes = p_indexes;
1451 for (v=0; v < kernel->height; v++) {
1452 for (u=0; u < kernel->width; u++, k--) {
1453 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1454 Maximize(result.red, k_pixels[u].red);
1455 Maximize(result.green, k_pixels[u].green);
1456 Maximize(result.blue, k_pixels[u].blue);
1457 Maximize(result.opacity, QuantumRange-k_pixels[u].opacity);
1458 if ( image->colorspace == CMYKColorspace)
1459 Maximize(result.index, k_indexes[u]);
1461 k_pixels += image->columns+kernel->width;
1462 k_indexes += image->columns+kernel->width;
1464 if ((channel & RedChannel) != 0)
1465 Assign(red,result.red);
1466 if ((channel & GreenChannel) != 0)
1467 Assign(green,result.green);
1468 if ((channel & BlueChannel) != 0)
1469 Assign(blue,result.blue);
1470 if ((channel & OpacityChannel) != 0
1471 && image->matte == MagickTrue )
1472 Assign(opacity,QuantumRange-result.opacity);
1473 if ((channel & IndexChannel) != 0
1474 && image->colorspace == CMYKColorspace)
1475 AssignIndex(result.index);
1478 case ErodeMorphology:
1479 /* Minimize Value within kernel shape
1481 * NOTE that the kernel is not reflected for this operation!
1483 * NOTE: in normal Greyscale Morphology, the kernel value should
1484 * be added to the real value, this is currently not done, due to
1485 * the nature of the boolean kernels being used.
1489 k_indexes = p_indexes;
1490 for (v=0; v < kernel->height; v++) {
1491 for (u=0; u < kernel->width; u++, k++) {
1492 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1493 Minimize(result.red, k_pixels[u].red);
1494 Minimize(result.green, k_pixels[u].green);
1495 Minimize(result.blue, k_pixels[u].blue);
1496 Minimize(result.opacity, QuantumRange-k_pixels[u].opacity);
1497 if ( image->colorspace == CMYKColorspace)
1498 Minimize(result.index, k_indexes[u]);
1500 k_pixels += image->columns+kernel->width;
1501 k_indexes += image->columns+kernel->width;
1503 if ((channel & RedChannel) != 0)
1504 Assign(red,result.red);
1505 if ((channel & GreenChannel) != 0)
1506 Assign(green,result.green);
1507 if ((channel & BlueChannel) != 0)
1508 Assign(blue,result.blue);
1509 if ((channel & OpacityChannel) != 0
1510 && image->matte == MagickTrue )
1511 Assign(opacity,QuantumRange-result.opacity);
1512 if ((channel & IndexChannel) != 0
1513 && image->colorspace == CMYKColorspace)
1514 AssignIndex(result.index);
1517 case DialateIntensityMorphology:
1518 /* Select pixel with maximum intensity within kernel shape
1520 * WARNING: the intensity test fails for CMYK and does not
1521 * take into account the moderating effect of teh alpha channel
1524 * NOTE for correct working of this operation for asymetrical
1525 * kernels, the kernel needs to be applied in its reflected form.
1526 * That is its values needs to be reversed.
1528 k = &kernel->values[ kernel->width*kernel->height-1 ];
1530 k_indexes = p_indexes;
1531 for (v=0; v < kernel->height; v++) {
1532 for (u=0; u < kernel->width; u++, k--) {
1533 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1534 if ( result.red == 0.0 ||
1535 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1536 /* copy the whole pixel - no channel selection */
1538 if ( result.red > 0.0 ) changed++;
1542 k_pixels += image->columns+kernel->width;
1543 k_indexes += image->columns+kernel->width;
1547 case ErodeIntensityMorphology:
1548 /* Select pixel with mimimum intensity within kernel shape
1550 * WARNING: the intensity test fails for CMYK and does not
1551 * take into account the moderating effect of teh alpha channel
1554 * NOTE that the kernel is not reflected for this operation!
1558 k_indexes = p_indexes;
1559 for (v=0; v < kernel->height; v++) {
1560 for (u=0; u < kernel->width; u++, k++) {
1561 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1562 if ( result.red == 0.0 ||
1563 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1564 /* copy the whole pixel - no channel selection */
1566 if ( result.red > 0.0 ) changed++;
1570 k_pixels += image->columns+kernel->width;
1571 k_indexes += image->columns+kernel->width;
1575 case DistanceMorphology:
1576 /* Add kernel value and select the minimum value found.
1577 * The result is a iterative distance from edge function.
1579 * All Distance Kernels are symetrical, but that may not always
1580 * be the case. For example how about a distance from left edges?
1581 * To make it work correctly for asymetrical kernels the reflected
1582 * kernel needs to be applied.
1585 /* No need to do distance morphology if all values are zero */
1586 /* Unfortunatally I have not been able to get this right! */
1587 if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1588 || ((channel & GreenChannel) == 0 && p[r].green == 0)
1589 || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1590 || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1591 || (( (channel & IndexChannel) == 0
1592 || image->colorspace != CMYKColorspace
1593 ) && p_indexes[x] ==0 )
1597 k = &kernel->values[ kernel->width*kernel->height-1 ];
1599 k_indexes = p_indexes;
1600 for (v=0; v < kernel->height; v++) {
1601 for (u=0; u < kernel->width; u++, k--) {
1602 if ( IsNan(*k) ) continue;
1603 Minimize(result.red, (*k)+k_pixels[u].red);
1604 Minimize(result.green, (*k)+k_pixels[u].green);
1605 Minimize(result.blue, (*k)+k_pixels[u].blue);
1606 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1607 if ( image->colorspace == CMYKColorspace)
1608 Minimize(result.index, (*k)+k_indexes[u]);
1610 k_pixels += image->columns+kernel->width;
1611 k_indexes += image->columns+kernel->width;
1614 if ((channel & RedChannel) != 0)
1615 Assign(red,result.red);
1616 if ((channel & GreenChannel) != 0)
1617 Assign(green,result.green);
1618 if ((channel & BlueChannel) != 0)
1619 Assign(blue,result.blue);
1620 if ((channel & OpacityChannel) != 0
1621 && image->matte == MagickTrue )
1622 Assign(opacity,QuantumRange-result.opacity);
1623 if ((channel & IndexChannel) != 0
1624 && image->colorspace == CMYKColorspace)
1625 AssignIndex(result.index);
1627 /* By returning the number of 'maximum' values still to process
1628 ** we can get the Distance iteration to finish faster.
1629 ** BUT this may cause an infinite loop on very large shapes,
1630 ** which may have a distance gradient that reachs max value!
1632 if ((channel & RedChannel) != 0)
1633 { q->red = ClampToQuantum(result.red);
1634 if ( q->red == QuantumRange ) changed++; /* more to do */
1636 if ((channel & GreenChannel) != 0)
1637 { q->green = ClampToQuantum(result.green);
1638 if ( q->green == QuantumRange ) changed++; /* more to do */
1640 if ((channel & BlueChannel) != 0)
1641 { q->blue = ClampToQuantum(result.blue);
1642 if ( q->blue == QuantumRange ) changed++; /* more to do */
1644 if ((channel & OpacityChannel) != 0)
1645 { q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1646 if ( q->opacity == 0 ) changed++; /* more to do */
1648 if (((channel & IndexChannel) != 0) &&
1649 (image->colorspace == CMYKColorspace))
1650 { q_indexes[x] = ClampToQuantum(result.index);
1651 if ( q_indexes[x] == QuantumRange ) changed++;
1656 case UndefinedMorphology:
1658 break; /* Do nothing */
1663 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1664 if (sync == MagickFalse)
1666 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1671 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1672 #pragma omp critical (MagickCore_MorphologyImage)
1674 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1675 if (proceed == MagickFalse)
1679 result_image->type=image->type;
1680 q_view=DestroyCacheView(q_view);
1681 p_view=DestroyCacheView(p_view);
1682 return(status ? changed : 0);
1685 MagickExport Image *MorphologyImage(const Image *image,MorphologyMethod method,
1686 const long iterations,KernelInfo *kernel, ExceptionInfo *exception)
1691 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1692 iterations,kernel,exception);
1693 return(morphology_image);
1696 MagickExport Image *MorphologyImageChannel(const Image *image,
1697 const ChannelType channel, MorphologyMethod method, const long iterations,
1698 KernelInfo *kernel, ExceptionInfo *exception)
1709 assert(image != (Image *) NULL);
1710 assert(image->signature == MagickSignature);
1711 assert(exception != (ExceptionInfo *) NULL);
1712 assert(exception->signature == MagickSignature);
1714 if ( GetImageArtifact(image,"showkernel") != (const char *) NULL)
1715 KernelPrint(kernel);
1717 if ( iterations == 0 )
1718 return((Image *)NULL); /* null operation - nothing to do! */
1720 /* kernel must be valid at this point
1721 * (except maybe for posible future morphology methods like "Prune"
1723 assert(kernel != (KernelInfo *)NULL);
1728 if ( iterations < 0 )
1729 limit = image->columns > image->rows ? image->columns : image->rows;
1731 /* Special morphology cases */
1733 case CloseMorphology:
1734 new_image = MorphologyImageChannel(image, channel, DialateMorphology,
1735 iterations, kernel, exception);
1736 if (new_image == (Image *) NULL)
1737 return((Image *) NULL);
1738 method = ErodeMorphology;
1740 case OpenMorphology:
1741 new_image = MorphologyImageChannel(image, channel, ErodeMorphology,
1742 iterations, kernel, exception);
1743 if (new_image == (Image *) NULL)
1744 return((Image *) NULL);
1745 method = DialateMorphology;
1747 case CloseIntensityMorphology:
1748 new_image = MorphologyImageChannel(image, channel, DialateIntensityMorphology,
1749 iterations, kernel, exception);
1750 if (new_image == (Image *) NULL)
1751 return((Image *) NULL);
1752 method = ErodeIntensityMorphology;
1754 case OpenIntensityMorphology:
1755 new_image = MorphologyImageChannel(image, channel, ErodeIntensityMorphology,
1756 iterations, kernel, exception);
1757 if (new_image == (Image *) NULL)
1758 return((Image *) NULL);
1759 method = DialateIntensityMorphology;
1762 case ConvolveMorphology:
1763 KernelNormalize(kernel);
1766 /* Do a morphology just once at this point!
1767 This ensures a new_image has been generated, but allows us
1768 to skip the creation of 'old_image' if it isn't needed.
1770 new_image=CloneImage(image,0,0,MagickTrue,exception);
1771 if (new_image == (Image *) NULL)
1772 return((Image *) NULL);
1773 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1775 InheritException(exception,&new_image->exception);
1776 new_image=DestroyImage(new_image);
1777 return((Image *) NULL);
1779 changed = MorphologyApply(image,new_image,method,channel,kernel,
1782 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1783 fprintf(stderr, "Morphology %s:%lu => Changed %lu\n",
1784 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1788 /* Repeat the interative morphology until count or no change */
1789 if ( count < limit && changed > 0 ) {
1790 old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1791 if (old_image == (Image *) NULL)
1792 return(DestroyImage(new_image));
1793 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1795 InheritException(exception,&old_image->exception);
1796 old_image=DestroyImage(old_image);
1797 return(DestroyImage(new_image));
1799 while( count < limit && changed != 0 )
1801 Image *tmp = old_image;
1802 old_image = new_image;
1804 changed = MorphologyApply(old_image,new_image,method,channel,kernel,
1807 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1808 fprintf(stderr, "Morphology %s:%lu => Changed %lu\n",
1809 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1812 DestroyImage(old_image);