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)
109 /* Currently these are only internal to this module */
111 RotateKernelInfo(KernelInfo *, double),
112 ScaleKernelInfo(KernelInfo *, double);
115 *CloneKernelInfo(KernelInfo *);
119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123 % A c q u i r e K e r n e l I n f o %
127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
129 % AcquireKernelInfo() takes the given string (generally supplied by the
130 % user) and converts it into a Morphology/Convolution Kernel. This allows
131 % users to specify a kernel from a number of pre-defined kernels, or to fully
132 % specify their own kernel for a specific Convolution or Morphology
135 % The kernel so generated can be any rectangular array of floating point
136 % values (doubles) with the 'control point' or 'pixel being affected'
137 % anywhere within that array of values.
139 % Previously IM was restricted to a square of odd size using the exact
140 % center as origin, this is no longer the case, and any rectangular kernel
141 % with any value being declared the origin. This in turn allows the use of
142 % highly asymmetrical kernels.
144 % The floating point values in the kernel can also include a special value
145 % known as 'nan' or 'not a number' to indicate that this value is not part
146 % of the kernel array. This allows you to shaped the kernel within its
147 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
148 % shape. However at least one non-nan value must be provided for correct
149 % working of a kernel.
151 % The returned kernel should be free using the DestroyKernelInfo() when you
152 % are finished with it.
154 % Input kernel defintion strings can consist of any of three types.
157 % Select from one of the built in kernels, using the name and
158 % geometry arguments supplied. See AcquireKernelBuiltIn()
160 % "WxH[+X+Y]:num, num, num ..."
161 % a kernal of size W by H, with W*H floating point numbers following.
162 % the 'center' can be optionally be defined at +X+Y (such that +0+0
163 % is top left corner). If not defined the pixel in the center, for
164 % odd sizes, or to the immediate top or left of center for even sizes
165 % is automatically selected.
167 % "num, num, num, num, ..."
168 % list of floating point numbers defining an 'old style' odd sized
169 % square kernel. At least 9 values should be provided for a 3x3
170 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
171 % Values can be space or comma separated. This is not recommended.
173 % Note that 'name' kernels will start with an alphabetic character while the
174 % new kernel specification has a ':' character in its specification string.
175 % If neither is the case, it is assumed an old style of a simple list of
176 % numbers generating a odd-sized square kernel has been given.
178 % The format of the AcquireKernal method is:
180 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
182 % A description of each parameter follows:
184 % o kernel_string: the Morphology/Convolution kernel wanted.
188 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
194 token[MaxTextExtent];
209 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
211 assert(kernel_string != (const char *) NULL);
212 SetGeometryInfo(&args);
214 /* does it start with an alpha - Return a builtin kernel */
215 GetMagickToken(kernel_string,&p,token);
216 if ( isalpha((int)token[0]) )
221 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
222 if ( type < 0 || type == UserDefinedKernel )
223 return((KernelInfo *)NULL);
225 while (((isspace((int) ((unsigned char) *p)) != 0) ||
226 (*p == ',') || (*p == ':' )) && (*p != '\0'))
228 flags = ParseGeometry(p, &args);
230 /* special handling of missing values in input string */
232 case RectangleKernel:
233 if ( (flags & WidthValue) == 0 ) /* if no width then */
234 args.rho = args.sigma; /* then width = height */
235 if ( args.rho < 1.0 ) /* if width too small */
236 args.rho = 3; /* then width = 3 */
237 if ( args.sigma < 1.0 ) /* if height too small */
238 args.sigma = args.rho; /* then height = width */
239 if ( (flags & XValue) == 0 ) /* center offset if not defined */
240 args.xi = (double)(((long)args.rho-1)/2);
241 if ( (flags & YValue) == 0 )
242 args.psi = (double)(((long)args.sigma-1)/2);
248 if ( (flags & HeightValue) == 0 ) /* if no scale */
249 args.sigma = 1.0; /* then scale = 1.0 */
255 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
258 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
259 if (kernel == (KernelInfo *)NULL)
261 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
262 kernel->type = UserDefinedKernel;
263 kernel->signature = MagickSignature;
265 /* Has a ':' in argument - New user kernel specification */
266 p = strchr(kernel_string, ':');
267 if ( p != (char *) NULL)
269 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
270 memcpy(token, kernel_string, (size_t) (p-kernel_string));
271 token[p-kernel_string] = '\0';
272 flags = ParseGeometry(token, &args);
274 /* Size handling and checks of geometry settings */
275 if ( (flags & WidthValue) == 0 ) /* if no width then */
276 args.rho = args.sigma; /* then width = height */
277 if ( args.rho < 1.0 ) /* if width too small */
278 args.rho = 1.0; /* then width = 1 */
279 if ( args.sigma < 1.0 ) /* if height too small */
280 args.sigma = args.rho; /* then height = width */
281 kernel->width = (unsigned long)args.rho;
282 kernel->height = (unsigned long)args.sigma;
284 /* Offset Handling and Checks */
285 if ( args.xi < 0.0 || args.psi < 0.0 )
286 return(DestroyKernelInfo(kernel));
287 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
288 : (long) (kernel->width-1)/2;
289 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
290 : (long) (kernel->height-1)/2;
291 if ( kernel->x >= (long) kernel->width ||
292 kernel->y >= (long) kernel->height )
293 return(DestroyKernelInfo(kernel));
295 p++; /* advance beyond the ':' */
298 { /* ELSE - Old old kernel specification, forming odd-square kernel */
299 /* count up number of values given */
300 p=(const char *) kernel_string;
301 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
302 p++; /* ignore "'" chars for convolve filter usage - Cristy */
303 for (i=0; *p != '\0'; i++)
305 GetMagickToken(p,&p,token);
307 GetMagickToken(p,&p,token);
309 /* set the size of the kernel - old sized square */
310 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
311 kernel->x = kernel->y = (long) (kernel->width-1)/2;
312 p=(const char *) kernel_string;
313 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
314 p++; /* ignore "'" chars for convolve filter usage - Cristy */
317 /* Read in the kernel values from rest of input string argument */
318 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
319 kernel->height*sizeof(double));
320 if (kernel->values == (double *) NULL)
321 return(DestroyKernelInfo(kernel));
323 kernel->minimum = +MagickHuge;
324 kernel->maximum = -MagickHuge;
325 kernel->negative_range = kernel->positive_range = 0.0;
326 for (i=0; (i < (long) (kernel->width*kernel->height)) && (*p != '\0'); i++)
328 GetMagickToken(p,&p,token);
330 GetMagickToken(p,&p,token);
331 if ( LocaleCompare("nan",token) == 0
332 || LocaleCompare("-",token) == 0 ) {
333 kernel->values[i] = nan; /* do not include this value in kernel */
336 kernel->values[i] = StringToDouble(token);
337 ( kernel->values[i] < 0)
338 ? ( kernel->negative_range += kernel->values[i] )
339 : ( kernel->positive_range += kernel->values[i] );
340 Minimize(kernel->minimum, kernel->values[i]);
341 Maximize(kernel->maximum, kernel->values[i]);
344 /* check that we recieved at least one real (non-nan) value! */
345 if ( kernel->minimum == MagickHuge )
346 return(DestroyKernelInfo(kernel));
348 /* This should not be needed for a fully defined kernel
349 * Perhaps an error should be reported instead!
350 * Kept for backward compatibility.
352 if ( i < (long) (kernel->width*kernel->height) ) {
353 Minimize(kernel->minimum, kernel->values[i]);
354 Maximize(kernel->maximum, kernel->values[i]);
355 for ( ; i < (long) (kernel->width*kernel->height); i++)
356 kernel->values[i]=0.0;
363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
367 % A c q u i r e K e r n e l B u i l t I n %
371 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
373 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
374 % kernels used for special purposes such as gaussian blurring, skeleton
375 % pruning, and edge distance determination.
377 % They take a KernelType, and a set of geometry style arguments, which were
378 % typically decoded from a user supplied string, or from a more complex
379 % Morphology Method that was requested.
381 % The format of the AcquireKernalBuiltIn method is:
383 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
384 % const GeometryInfo args)
386 % A description of each parameter follows:
388 % o type: the pre-defined type of kernel wanted
390 % o args: arguments defining or modifying the kernel
392 % Convolution Kernels
394 % Gaussian "{radius},{sigma}"
395 % Generate a two-dimentional gaussian kernel, as used by -gaussian
396 % A sigma is required, (with the 'x'), due to historical reasons.
398 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
399 % the final size of the resulting kernel to a square 2*radius+1 in size.
400 % The radius should be at least 2 times that of the sigma value, or
401 % sever clipping and aliasing may result. If not given or set to 0 the
402 % radius will be determined so as to produce the best minimal error
403 % result, which is usally much larger than is normally needed.
405 % Blur "{radius},{sigma},{angle}"
406 % As per Gaussian, but generates a 1 dimensional or linear gaussian
407 % blur, at the angle given (current restricted to orthogonal angles).
408 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
410 % NOTE that two such blurs perpendicular to each other is equivelent to
411 % -blur and the previous gaussian, but is often 10 or more times faster.
413 % Comet "{width},{sigma},{angle}"
414 % Blur in one direction only, mush like how a bright object leaves
415 % a comet like trail. The Kernel is actually half a gaussian curve,
416 % Adding two such blurs in oppiste directions produces a Linear Blur.
418 % NOTE: that the first argument is the width of the kernel and not the
419 % radius of the kernel.
421 % # Still to be implemented...
423 % # Sharpen "{radius},{sigma}
424 % # Negated Gaussian (center zeroed and re-normalized),
425 % # with a 2 unit positive peak. -- Check On line documentation
427 % # Laplacian "{radius},{sigma}"
428 % # Laplacian (a mexican hat like) Function
430 % # LOG "{radius},{sigma1},{sigma2}
431 % # Laplacian of Gaussian
433 % # DOG "{radius},{sigma1},{sigma2}
434 % # Difference of two Gaussians
438 % # Set kernel values using a resize filter, and given scale (sigma)
439 % # Cylindrical or Linear. Is this posible with an image?
444 % Rectangle "{geometry}"
445 % Simply generate a rectangle of 1's with the size given. You can also
446 % specify the location of the 'control point', otherwise the closest
447 % pixel to the center of the rectangle is selected.
449 % Properly centered and odd sized rectangles work the best.
451 % Diamond "[{radius}[,{scale}]]"
452 % Generate a diamond shaped kernal with given radius to the points.
453 % Kernel size will again be radius*2+1 square and defaults to radius 1,
454 % generating a 3x3 kernel that is slightly larger than a square.
456 % Square "[{radius}[,{scale}]]"
457 % Generate a square shaped kernel of size radius*2+1, and defaulting
458 % to a 3x3 (radius 1).
460 % Note that using a larger radius for the "Square" or the "Diamond"
461 % is also equivelent to iterating the basic morphological method
462 % that many times. However However iterating with the smaller radius 1
463 % default is actually faster than using a larger kernel radius.
465 % Disk "[{radius}[,{scale}]]
466 % Generate a binary disk of the radius given, radius may be a float.
467 % Kernel size will be ceil(radius)*2+1 square.
468 % NOTE: Here are some disk shapes of specific interest
469 % "disk:1" => "diamond" or "cross:1"
470 % "disk:1.5" => "square"
471 % "disk:2" => "diamond:2"
472 % "disk:2.5" => a general disk shape of radius 2
473 % "disk:2.9" => "square:2"
474 % "disk:3.5" => default - octagonal/disk shape of radius 3
475 % "disk:4.2" => roughly octagonal shape of radius 4
476 % "disk:4.3" => a general disk shape of radius 4
477 % After this all the kernel shape becomes more and more circular.
479 % Because a "disk" is more circular when using a larger radius, using a
480 % larger radius is preferred over iterating the morphological operation.
482 % Plus "[{radius}[,{scale}]]"
483 % Generate a kernel in the shape of a 'plus' sign. The length of each
484 % arm is also the radius, which defaults to 2.
486 % This kernel is not a good general morphological kernel, but is used
487 % more for highlighting and marking any single pixels in an image using,
488 % a "Dilate" or "Erode" method as appropriate.
490 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
492 % Note that unlike other kernels iterating a plus does not produce the
493 % same result as using a larger radius for the cross.
495 % Distance Measuring Kernels
497 % Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1)
498 % Manhatten "[{radius}][x{scale}]" square grid distance (default r=1)
499 % Euclidean "[{radius}][x{scale}]" direct distance (default r=1)
501 % Different types of distance measuring methods, which are used with the
502 % a 'Distance' morphology method for generating a gradient based on
503 % distance from an edge of a binary shape, though there is a technique
504 % for handling a anti-aliased shape.
506 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
507 % one to any neighbour, orthogonal or diagonal. One why of thinking of
508 % it is the number of squares a 'King' or 'Queen' in chess needs to
509 % traverse reach any other position on a chess board. It results in a
510 % 'square' like distance function, but one where diagonals are closer
513 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
514 % Cab metric), is the distance needed when you can only travel in
515 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
516 % in chess would travel. It results in a diamond like distances, where
517 % diagonals are further than expected.
519 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
520 % However by default the kernel size only has a radius of 1, which
521 % limits the distance to 'Knight' like moves, with only orthogonal and
522 % diagonal measurements being correct. As such for the default kernel
523 % you will get octagonal like distance function, which is reasonally
526 % However if you use a larger radius such as "Euclidean:4" you will
527 % get a much smoother distance gradient from the edge of the shape.
528 % Of course a larger kernel is slower to use, and generally not needed.
530 % To allow the use of fractional distances that you get with diagonals
531 % the actual distance is scaled by a fixed value which the user can
532 % provide. This is not actually nessary for either ""Chebyshev" or
533 % "Manhatten" distance kernels, but is done for all three distance
534 % kernels. If no scale is provided it is set to a value of 100,
535 % allowing for a maximum distance measurement of 655 pixels using a Q16
536 % version of IM, from any edge. However for small images this can
537 % result in quite a dark gradient.
539 % See the 'Distance' Morphological Method, for information of how it is
542 % # Hit-n-Miss Kernel-Lists -- Still to be implemented
544 % # specifically for Pruning, Thinning, Thickening
548 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
549 const GeometryInfo *args)
562 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
564 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
565 if (kernel == (KernelInfo *) NULL)
567 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
568 kernel->minimum = kernel->maximum = 0.0;
569 kernel->negative_range = kernel->positive_range = 0.0;
571 kernel->signature = MagickSignature;
574 /* Convolution Kernels */
577 sigma = fabs(args->sigma);
579 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
581 kernel->width = kernel->height =
582 GetOptimalKernelWidth2D(args->rho,sigma);
583 kernel->x = kernel->y = (long) (kernel->width-1)/2;
584 kernel->negative_range = kernel->positive_range = 0.0;
585 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
586 kernel->height*sizeof(double));
587 if (kernel->values == (double *) NULL)
588 return(DestroyKernelInfo(kernel));
590 sigma = 2.0*sigma*sigma; /* simplify the expression */
591 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
592 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
593 kernel->positive_range += (
595 exp(-((double)(u*u+v*v))/sigma)
596 /* / (MagickPI*sigma) */ );
598 kernel->maximum = kernel->values[
599 kernel->y*kernel->width+kernel->x ];
601 ScaleKernelInfo(kernel, 0.0); /* Normalize Kernel Values */
607 sigma = fabs(args->sigma);
609 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
611 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
612 kernel->x = (long) (kernel->width-1)/2;
615 kernel->negative_range = kernel->positive_range = 0.0;
616 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
617 kernel->height*sizeof(double));
618 if (kernel->values == (double *) NULL)
619 return(DestroyKernelInfo(kernel));
623 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
624 ** It generates a gaussian 3 times the width, and compresses it into
625 ** the expected range. This produces a closer normalization of the
626 ** resulting kernel, especially for very low sigma values.
627 ** As such while wierd it is prefered.
629 ** I am told this method originally came from Photoshop.
631 sigma *= KernelRank; /* simplify expanded curve */
632 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
633 (void) ResetMagickMemory(kernel->values,0, (size_t)
634 kernel->width*sizeof(double));
635 for ( u=-v; u <= v; u++) {
636 kernel->values[(u+v)/KernelRank] +=
637 exp(-((double)(u*u))/(2.0*sigma*sigma))
638 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
640 for (i=0; i < (long) kernel->width; i++)
641 kernel->positive_range += kernel->values[i];
643 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
644 kernel->positive_range += (
646 exp(-((double)(u*u))/(2.0*sigma*sigma))
647 /* / (MagickSQ2PI*sigma) */ );
650 kernel->maximum = kernel->values[ kernel->x ];
651 /* Note that neither methods above generate a normalized kernel,
652 ** though it gets close. The kernel may be 'clipped' by a user defined
653 ** radius, producing a smaller (darker) kernel. Also for very small
654 ** sigma's (> 0.1) the central value becomes larger than one, and thus
655 ** producing a very bright kernel.
658 /* Normalize the 1D Gaussian Kernel
660 ** Because of this the divisor in the above kernel generator is
661 ** not needed, so is not done above.
663 ScaleKernelInfo(kernel, 0.0); /* Normalize Kernel Values */
665 /* rotate the kernel by given angle */
666 RotateKernelInfo(kernel, args->xi);
671 sigma = fabs(args->sigma);
673 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
675 if ( args->rho < 1.0 )
676 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
678 kernel->width = (unsigned long)args->rho;
679 kernel->x = kernel->y = 0;
681 kernel->negative_range = kernel->positive_range = 0.0;
682 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
683 kernel->height*sizeof(double));
684 if (kernel->values == (double *) NULL)
685 return(DestroyKernelInfo(kernel));
687 /* A comet blur is half a gaussian curve, so that the object is
688 ** blurred in one direction only. This may not be quite the right
689 ** curve so may change in the future. The function must be normalised.
693 sigma *= KernelRank; /* simplify expanded curve */
694 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
695 (void) ResetMagickMemory(kernel->values,0, (size_t)
696 kernel->width*sizeof(double));
697 for ( u=0; u < v; u++) {
698 kernel->values[u/KernelRank] +=
699 exp(-((double)(u*u))/(2.0*sigma*sigma))
700 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
702 for (i=0; i < (long) kernel->width; i++)
703 kernel->positive_range += kernel->values[i];
705 for ( i=0; i < (long) kernel->width; i++)
706 kernel->positive_range += (
708 exp(-((double)(i*i))/(2.0*sigma*sigma))
709 /* / (MagickSQ2PI*sigma) */ );
712 kernel->maximum = kernel->values[0];
714 ScaleKernelInfo(kernel, 0.0); /* Normalize Kernel Values */
715 RotateKernelInfo(kernel, args->xi);
718 /* Boolean Kernels */
719 case RectangleKernel:
723 if ( type == SquareKernel )
726 kernel->width = kernel->height = 3; /* default radius = 1 */
728 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
729 kernel->x = kernel->y = (long) (kernel->width-1)/2;
733 /* NOTE: user defaults set in "AcquireKernelInfo()" */
734 if ( args->rho < 1.0 || args->sigma < 1.0 )
735 return(DestroyKernelInfo(kernel)); /* invalid args given */
736 kernel->width = (unsigned long)args->rho;
737 kernel->height = (unsigned long)args->sigma;
738 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
739 args->psi < 0.0 || args->psi > (double)kernel->height )
740 return(DestroyKernelInfo(kernel)); /* invalid args given */
741 kernel->x = (long) args->xi;
742 kernel->y = (long) args->psi;
745 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
746 kernel->height*sizeof(double));
747 if (kernel->values == (double *) NULL)
748 return(DestroyKernelInfo(kernel));
750 /* set all kernel values to 1.0 */
751 u=(long) kernel->width*kernel->height;
752 for ( i=0; i < u; i++)
753 kernel->values[i] = scale;
754 kernel->minimum = kernel->maximum = scale; /* a flat shape */
755 kernel->positive_range = scale*u;
761 kernel->width = kernel->height = 3; /* default radius = 1 */
763 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
764 kernel->x = kernel->y = (long) (kernel->width-1)/2;
766 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
767 kernel->height*sizeof(double));
768 if (kernel->values == (double *) NULL)
769 return(DestroyKernelInfo(kernel));
771 /* set all kernel values within diamond area to scale given */
772 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
773 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
774 if ((labs(u)+labs(v)) <= (long)kernel->x)
775 kernel->positive_range += kernel->values[i] = args->sigma;
777 kernel->values[i] = nan;
778 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
786 limit = (long)(args->rho*args->rho);
787 if (args->rho < 0.1) /* default radius approx 3.5 */
788 kernel->width = kernel->height = 7L, limit = 10L;
790 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
791 kernel->x = kernel->y = (long) (kernel->width-1)/2;
793 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
794 kernel->height*sizeof(double));
795 if (kernel->values == (double *) NULL)
796 return(DestroyKernelInfo(kernel));
798 /* set all kernel values within disk area to 1.0 */
799 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
800 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
801 if ((u*u+v*v) <= limit)
802 kernel->positive_range += kernel->values[i] = args->sigma;
804 kernel->values[i] = nan;
805 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
811 kernel->width = kernel->height = 5; /* default radius 2 */
813 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
814 kernel->x = kernel->y = (long) (kernel->width-1)/2;
816 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
817 kernel->height*sizeof(double));
818 if (kernel->values == (double *) NULL)
819 return(DestroyKernelInfo(kernel));
821 /* set all kernel values along axises to 1.0 */
822 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
823 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
824 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
825 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
826 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
829 /* Distance Measuring Kernels */
830 case ChebyshevKernel:
836 kernel->width = kernel->height = 3; /* default radius = 1 */
838 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
839 kernel->x = kernel->y = (long) (kernel->width-1)/2;
841 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
842 kernel->height*sizeof(double));
843 if (kernel->values == (double *) NULL)
844 return(DestroyKernelInfo(kernel));
846 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
847 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
848 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
849 kernel->positive_range += ( kernel->values[i] =
850 scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
851 kernel->maximum = kernel->values[0];
854 case ManhattenKernel:
860 kernel->width = kernel->height = 3; /* default radius = 1 */
862 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
863 kernel->x = kernel->y = (long) (kernel->width-1)/2;
865 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
866 kernel->height*sizeof(double));
867 if (kernel->values == (double *) NULL)
868 return(DestroyKernelInfo(kernel));
870 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
871 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
872 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
873 kernel->positive_range += ( kernel->values[i] =
874 scale*(labs(u)+labs(v)) );
875 kernel->maximum = kernel->values[0];
878 case EuclideanKernel:
884 kernel->width = kernel->height = 3; /* default radius = 1 */
886 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
887 kernel->x = kernel->y = (long) (kernel->width-1)/2;
889 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
890 kernel->height*sizeof(double));
891 if (kernel->values == (double *) NULL)
892 return(DestroyKernelInfo(kernel));
894 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
895 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
896 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
897 kernel->positive_range += ( kernel->values[i] =
898 scale*sqrt((double)(u*u+v*v)) );
899 kernel->maximum = kernel->values[0];
902 /* Undefined Kernels */
903 case LaplacianKernel:
906 perror("Kernel Type has not been defined yet");
909 /* Generate a No-Op minimal kernel - 1x1 pixel */
910 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
911 if (kernel->values == (double *) NULL)
912 return(DestroyKernelInfo(kernel));
913 kernel->width = kernel->height = 1;
914 kernel->x = kernel->x = 0;
915 kernel->type = UndefinedKernel;
917 kernel->positive_range =
918 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
926 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
930 + C l o n e K e r n e l I n f o %
934 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
936 % CloneKernelInfo() creates a new clone of the given Kernel so that its can
937 % be modified without effecting the original. The cloned kernel should be
938 % destroyed using DestoryKernelInfo() when no longer needed.
940 % The format of the DestroyKernelInfo method is:
942 % KernelInfo *CloneKernelInfo(KernelInfo *kernel)
944 % A description of each parameter follows:
946 % o kernel: the Morphology/Convolution kernel to be cloned
950 static KernelInfo *CloneKernelInfo(KernelInfo *kernel)
958 assert(kernel != (KernelInfo *) NULL);
960 new=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
961 if (new == (KernelInfo *) NULL)
963 *new = *kernel; /* copy values in structure */
965 new->values=(double *) AcquireQuantumMemory(kernel->width,
966 kernel->height*sizeof(double));
967 if (new->values == (double *) NULL)
968 return(DestroyKernelInfo(new));
970 for (i=0; i < (long) (kernel->width*kernel->height); i++)
971 new->values[i] = kernel->values[i];
977 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
981 % D e s t r o y K e r n e l I n f o %
985 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
987 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
990 % The format of the DestroyKernelInfo method is:
992 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
994 % A description of each parameter follows:
996 % o kernel: the Morphology/Convolution kernel to be destroyed
1000 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1002 assert(kernel != (KernelInfo *) NULL);
1004 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1005 kernel->height*sizeof(double));
1006 kernel->values=(double *)RelinquishMagickMemory(kernel->values);
1007 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
1012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1016 % M o r p h o l o g y I m a g e C h a n n e l %
1020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1022 % MorphologyImageChannel() applies a user supplied kernel to the image
1023 % according to the given mophology method.
1025 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1026 % by the kernel generator.
1028 % The format of the MorphologyImage method is:
1030 % Image *MorphologyImage(const Image *image, MorphologyMethod method,
1031 % const long iterations, KernelInfo *kernel, ExceptionInfo *exception)
1032 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1033 % channel, MorphologyMethod method, const long iterations, KernelInfo
1034 % *kernel, ExceptionInfo *exception)
1036 % A description of each parameter follows:
1038 % o image: the image.
1040 % o method: the morphology method to be applied.
1042 % o iterations: apply the operation this many times (or no change).
1043 % A value of -1 means loop until no change found.
1044 % How this is applied may depend on the morphology method.
1045 % Typically this is a value of 1.
1047 % o channel: the channel type.
1049 % o kernel: An array of double representing the morphology kernel.
1050 % Warning: kernel may be normalized for the Convolve method.
1052 % o exception: return any errors or warnings in this structure.
1055 % TODO: bias and auto-scale handling of the kernel for convolution
1056 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1057 % by the kernel generator.
1061 /* Internal function
1062 * Apply the Morphology method with the given Kernel
1063 * And return the number of pixels that changed.
1065 static unsigned long MorphologyApply(const Image *image, Image
1066 *result_image, const MorphologyMethod method, const ChannelType channel,
1067 const KernelInfo *kernel, ExceptionInfo *exception)
1069 #define MorphologyTag "Morphology/Image"
1086 /* Only the most basic morphology is actually performed by this routine */
1087 assert( method <= DistanceMorphology );
1090 Apply Basic Morphology to Image.
1096 GetMagickPixelPacket(image,&bias);
1097 SetMagickPixelPacketBias(image,&bias);
1098 /* Future: handle auto-bias from user, based on kernel input */
1100 p_view=AcquireCacheView(image);
1101 q_view=AcquireCacheView(result_image);
1103 /* Some methods (including convolve) needs use a reflected kernel.
1104 * Adjust 'origin' offsets for this reflected kernel.
1109 case ErodeMorphology:
1110 case ErodeIntensityMorphology:
1111 /* kernel is not reflected */
1114 /* kernel needs to be reflected */
1115 offx = (long) kernel->width-offx-1;
1116 offy = (long) kernel->height-offy-1;
1120 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1121 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1123 for (y=0; y < (long) image->rows; y++)
1128 register const PixelPacket
1131 register const IndexPacket
1132 *restrict p_indexes;
1134 register PixelPacket
1137 register IndexPacket
1138 *restrict q_indexes;
1146 if (status == MagickFalse)
1148 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1149 image->columns+kernel->width, kernel->height, exception);
1150 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1152 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1157 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1158 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1159 r = (image->columns+kernel->width)*offy+offx; /* constant */
1161 for (x=0; x < (long) image->columns; x++)
1169 register const double
1172 register const PixelPacket
1175 register const IndexPacket
1176 *restrict k_indexes;
1181 /* Copy input to ouput image for unused channels
1182 * This removes need for 'cloning' a new image every iteration
1185 if (image->colorspace == CMYKColorspace)
1186 q_indexes[x] = p_indexes[r];
1188 result.index=(MagickRealType) 0; /* stop compiler warnings */
1190 case ConvolveMorphology:
1192 break; /* default result is the convolution bias */
1193 case DilateMorphology:
1198 result.index = -MagickHuge;
1200 case ErodeMorphology:
1205 result.index = +MagickHuge;
1207 case DilateIntensityMorphology:
1208 case ErodeIntensityMorphology:
1209 result.red = 0.0; /* flag indicating first match found */
1212 /* Otherwise just start with the original pixel value */
1213 result.red = (MagickRealType) p[r].red;
1214 result.green = (MagickRealType) p[r].green;
1215 result.blue = (MagickRealType) p[r].blue;
1216 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1217 if ( image->colorspace == CMYKColorspace)
1218 result.index = (MagickRealType) p_indexes[r];
1223 case ConvolveMorphology:
1224 /* Weighted Average of pixels
1226 * NOTE for correct working of this operation for asymetrical
1227 * kernels, the kernel needs to be applied in its reflected form.
1228 * That is its values needs to be reversed.
1230 if (((channel & OpacityChannel) == 0) ||
1231 (image->matte == MagickFalse))
1233 /* Convolution (no transparency) */
1234 k = &kernel->values[ kernel->width*kernel->height-1 ];
1236 k_indexes = p_indexes;
1237 for (v=0; v < (long) kernel->height; v++) {
1238 for (u=0; u < (long) kernel->width; u++, k--) {
1239 if ( IsNan(*k) ) continue;
1240 result.red += (*k)*k_pixels[u].red;
1241 result.green += (*k)*k_pixels[u].green;
1242 result.blue += (*k)*k_pixels[u].blue;
1243 /* result.opacity += not involved here */
1244 if ( image->colorspace == CMYKColorspace)
1245 result.index += (*k)*k_indexes[u];
1247 k_pixels += image->columns+kernel->width;
1248 k_indexes += image->columns+kernel->width;
1252 { /* Kernel & Alpha weighted Convolution */
1254 alpha, /* alpha value * kernel weighting */
1255 gamma; /* weighting divisor */
1258 k = &kernel->values[ kernel->width*kernel->height-1 ];
1260 k_indexes = p_indexes;
1261 for (v=0; v < (long) kernel->height; v++) {
1262 for (u=0; u < (long) kernel->width; u++, k--) {
1263 if ( IsNan(*k) ) continue;
1264 alpha=(*k)*(QuantumScale*(QuantumRange-
1265 k_pixels[u].opacity));
1267 result.red += alpha*k_pixels[u].red;
1268 result.green += alpha*k_pixels[u].green;
1269 result.blue += alpha*k_pixels[u].blue;
1270 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1271 if ( image->colorspace == CMYKColorspace)
1272 result.index += alpha*k_indexes[u];
1274 k_pixels += image->columns+kernel->width;
1275 k_indexes += image->columns+kernel->width;
1277 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1278 result.red *= gamma;
1279 result.green *= gamma;
1280 result.blue *= gamma;
1281 result.opacity *= gamma;
1282 result.index *= gamma;
1286 case ErodeMorphology:
1287 /* Minimize Value within kernel shape
1289 * NOTE that the kernel is not reflected for this operation!
1291 * NOTE: in normal Greyscale Morphology, the kernel value should
1292 * be added to the real value, this is currently not done, due to
1293 * the nature of the boolean kernels being used.
1297 k_indexes = p_indexes;
1298 for (v=0; v < (long) kernel->height; v++) {
1299 for (u=0; u < (long) kernel->width; u++, k++) {
1300 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1301 Minimize(result.red, (double) k_pixels[u].red);
1302 Minimize(result.green, (double) k_pixels[u].green);
1303 Minimize(result.blue, (double) k_pixels[u].blue);
1304 Minimize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1305 if ( image->colorspace == CMYKColorspace)
1306 Minimize(result.index, (double) k_indexes[u]);
1308 k_pixels += image->columns+kernel->width;
1309 k_indexes += image->columns+kernel->width;
1313 case DilateMorphology:
1314 /* Maximize Value within kernel shape
1316 * NOTE for correct working of this operation for asymetrical
1317 * kernels, the kernel needs to be applied in its reflected form.
1318 * That is its values needs to be reversed.
1320 * NOTE: in normal Greyscale Morphology, the kernel value should
1321 * be added to the real value, this is currently not done, due to
1322 * the nature of the boolean kernels being used.
1325 k = &kernel->values[ kernel->width*kernel->height-1 ];
1327 k_indexes = p_indexes;
1328 for (v=0; v < (long) kernel->height; v++) {
1329 for (u=0; u < (long) kernel->width; u++, k--) {
1330 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1331 Maximize(result.red, (double) k_pixels[u].red);
1332 Maximize(result.green, (double) k_pixels[u].green);
1333 Maximize(result.blue, (double) k_pixels[u].blue);
1334 Maximize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1335 if ( image->colorspace == CMYKColorspace)
1336 Maximize(result.index, (double) k_indexes[u]);
1338 k_pixels += image->columns+kernel->width;
1339 k_indexes += image->columns+kernel->width;
1343 case ErodeIntensityMorphology:
1344 /* Select pixel with mimimum intensity within kernel shape
1346 * WARNING: the intensity test fails for CMYK and does not
1347 * take into account the moderating effect of teh alpha channel
1350 * NOTE that the kernel is not reflected for this operation!
1354 k_indexes = p_indexes;
1355 for (v=0; v < (long) kernel->height; v++) {
1356 for (u=0; u < (long) kernel->width; u++, k++) {
1357 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1358 if ( result.red == 0.0 ||
1359 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1360 /* copy the whole pixel - no channel selection */
1362 if ( result.red > 0.0 ) changed++;
1366 k_pixels += image->columns+kernel->width;
1367 k_indexes += image->columns+kernel->width;
1371 case DilateIntensityMorphology:
1372 /* Select pixel with maximum intensity within kernel shape
1374 * WARNING: the intensity test fails for CMYK and does not
1375 * take into account the moderating effect of teh alpha channel
1378 * NOTE for correct working of this operation for asymetrical
1379 * kernels, the kernel needs to be applied in its reflected form.
1380 * That is its values needs to be reversed.
1382 k = &kernel->values[ kernel->width*kernel->height-1 ];
1384 k_indexes = p_indexes;
1385 for (v=0; v < (long) kernel->height; v++) {
1386 for (u=0; u < (long) kernel->width; u++, k--) {
1387 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1388 if ( result.red == 0.0 ||
1389 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1390 /* copy the whole pixel - no channel selection */
1392 if ( result.red > 0.0 ) changed++;
1396 k_pixels += image->columns+kernel->width;
1397 k_indexes += image->columns+kernel->width;
1401 case DistanceMorphology:
1402 /* Add kernel value and select the minimum value found.
1403 * The result is a iterative distance from edge function.
1405 * All Distance Kernels are symetrical, but that may not always
1406 * be the case. For example how about a distance from left edges?
1407 * To make it work correctly for asymetrical kernels the reflected
1408 * kernel needs to be applied.
1411 /* No need to do distance morphology if original value is zero
1412 * Unfortunatally I have not been able to get this right
1413 * when channel selection also becomes involved. -- Arrgghhh
1415 if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1416 || ((channel & GreenChannel) == 0 && p[r].green == 0)
1417 || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1418 || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1419 || (( (channel & IndexChannel) == 0
1420 || image->colorspace != CMYKColorspace
1421 ) && p_indexes[x] ==0 )
1425 k = &kernel->values[ kernel->width*kernel->height-1 ];
1427 k_indexes = p_indexes;
1428 for (v=0; v < (long) kernel->height; v++) {
1429 for (u=0; u < (long) kernel->width; u++, k--) {
1430 if ( IsNan(*k) ) continue;
1431 Minimize(result.red, (*k)+k_pixels[u].red);
1432 Minimize(result.green, (*k)+k_pixels[u].green);
1433 Minimize(result.blue, (*k)+k_pixels[u].blue);
1434 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1435 if ( image->colorspace == CMYKColorspace)
1436 Minimize(result.index, (*k)+k_indexes[u]);
1438 k_pixels += image->columns+kernel->width;
1439 k_indexes += image->columns+kernel->width;
1443 case UndefinedMorphology:
1445 break; /* Do nothing */
1448 case UndefinedMorphology:
1449 case DilateIntensityMorphology:
1450 case ErodeIntensityMorphology:
1451 break; /* full pixel was directly assigned */
1453 /* Assign the results */
1454 if ((channel & RedChannel) != 0)
1455 q->red = ClampToQuantum(result.red);
1456 if ((channel & GreenChannel) != 0)
1457 q->green = ClampToQuantum(result.green);
1458 if ((channel & BlueChannel) != 0)
1459 q->blue = ClampToQuantum(result.blue);
1460 if ((channel & OpacityChannel) != 0
1461 && image->matte == MagickTrue )
1462 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1463 if ((channel & IndexChannel) != 0
1464 && image->colorspace == CMYKColorspace)
1465 q_indexes[x] = ClampToQuantum(result.index);
1468 if ( ( p[r].red != q->red )
1469 || ( p[r].green != q->green )
1470 || ( p[r].blue != q->blue )
1471 || ( p[r].opacity != q->opacity )
1472 || ( image->colorspace == CMYKColorspace &&
1473 p_indexes[r] != q_indexes[x] ) )
1474 changed++; /* The pixel had some value changed! */
1478 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1479 if (sync == MagickFalse)
1481 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1486 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1487 #pragma omp critical (MagickCore_MorphologyImage)
1489 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1490 if (proceed == MagickFalse)
1494 result_image->type=image->type;
1495 q_view=DestroyCacheView(q_view);
1496 p_view=DestroyCacheView(p_view);
1497 return(status ? (unsigned long) changed : 0);
1501 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1502 method, const long iterations,KernelInfo *kernel, ExceptionInfo *exception)
1507 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1508 iterations,kernel,exception);
1509 return(morphology_image);
1513 MagickExport Image *MorphologyImageChannel(const Image *image,
1514 const ChannelType channel, const MorphologyMethod method,
1515 const long iterations, KernelInfo *kernel, ExceptionInfo *exception)
1538 assert(image != (Image *) NULL);
1539 assert(image->signature == MagickSignature);
1540 assert(kernel != (KernelInfo *) NULL);
1541 assert(kernel->signature == MagickSignature);
1542 assert(exception != (ExceptionInfo *) NULL);
1543 assert(exception->signature == MagickSignature);
1545 if ( iterations == 0 )
1546 return((Image *)NULL); /* null operation - nothing to do! */
1548 /* kernel must be valid at this point
1549 * (except maybe for posible future morphology methods like "Prune"
1551 assert(kernel != (KernelInfo *)NULL);
1553 count = 0; /* interation count */
1554 changed = 1; /* if compound method assume image was changed */
1555 curr_kernel = kernel; /* allow kernel to be cloned and modified */
1556 curr_method = method; /* and the method to be changed */
1558 limit = (unsigned long) iterations;
1559 if ( iterations < 0 )
1560 limit = image->columns > image->rows ? image->columns : image->rows;
1562 /* Third-level morphology methods */
1563 switch( curr_method ) {
1564 case EdgeMorphology:
1565 grad_image = MorphologyImageChannel(image, channel,
1566 DilateMorphology, iterations, curr_kernel, exception);
1568 case EdgeInMorphology:
1569 curr_method = ErodeMorphology;
1571 case EdgeOutMorphology:
1572 curr_method = DilateMorphology;
1574 case TopHatMorphology:
1575 curr_method = OpenMorphology;
1577 case BottomHatMorphology:
1578 curr_method = CloseMorphology;
1584 /* Second-level morphology methods */
1585 switch( curr_method ) {
1586 case OpenMorphology: /* Erode then Dilate without reflection */
1587 new_image = MorphologyImageChannel(image, channel,
1588 ErodeMorphology, iterations, curr_kernel, exception);
1589 if (new_image == (Image *) NULL)
1590 return((Image *) NULL);
1591 curr_method = DilateMorphology;
1593 case CloseMorphology: /* Dilate then Erode using reflected kernel */
1594 curr_kernel = CloneKernelInfo(kernel);
1595 RotateKernelInfo(curr_kernel,180);
1596 new_image = MorphologyImageChannel(image, channel,
1597 DilateMorphology, iterations, curr_kernel, exception);
1598 if (new_image == (Image *) NULL)
1599 return((Image *) NULL);
1600 curr_method = ErodeMorphology;
1602 case OpenIntensityMorphology:
1603 new_image = MorphologyImageChannel(image, channel,
1604 ErodeIntensityMorphology, iterations, curr_kernel, exception);
1605 if (new_image == (Image *) NULL)
1606 return((Image *) NULL);
1607 curr_method = DilateIntensityMorphology;
1609 case CloseIntensityMorphology:
1610 curr_kernel = CloneKernelInfo(kernel); /* a reflected kernel is needed */
1611 RotateKernelInfo(curr_kernel,180);
1612 new_image = MorphologyImageChannel(image, channel,
1613 DilateIntensityMorphology, iterations, curr_kernel, exception);
1614 if (new_image == (Image *) NULL)
1615 return((Image *) NULL);
1616 curr_method = ErodeIntensityMorphology;
1619 case ConvolveMorphology:
1620 /* Scale or Normalize kernel, according to user wishes
1621 ** before using it for the convolution method.
1623 artifact = GetImageArtifact(image,"convolve:scale");
1624 if ( artifact != (char *)NULL ) {
1625 curr_kernel = CloneKernelInfo(kernel);
1626 ScaleKernelInfo(curr_kernel, StringToDouble(artifact) );
1631 /* Do a iteration using a Basic Morphology method just once!
1632 ** This ensures a new_image has been generated, but allows us
1633 ** to skip the creation of 'old_image' if it isn't needed.
1635 new_image=CloneImage(image,0,0,MagickTrue,exception);
1636 if (new_image == (Image *) NULL)
1637 return((Image *) NULL);
1638 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1640 InheritException(exception,&new_image->exception);
1641 new_image=DestroyImage(new_image);
1642 return((Image *) NULL);
1644 changed = MorphologyApply(image,new_image,curr_method,channel,curr_kernel,
1647 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1648 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1649 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1653 /* Repeat a Basic morphology until count or no change reached */
1654 if ( count < (long) limit && changed > 0 ) {
1655 old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1656 if (old_image == (Image *) NULL)
1657 return(DestroyImage(new_image));
1658 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1660 InheritException(exception,&old_image->exception);
1661 old_image=DestroyImage(old_image);
1662 return(DestroyImage(new_image));
1664 while( count < (long) limit && changed != 0 )
1666 Image *tmp = old_image;
1667 old_image = new_image;
1669 changed = MorphologyApply(old_image,new_image,curr_method,channel,
1670 curr_kernel, exception);
1672 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1673 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1674 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1677 old_image=DestroyImage(old_image);
1679 if ( curr_kernel != kernel )
1680 curr_kernel=DestroyKernelInfo(curr_kernel);
1682 /* Subtractive morphology cases */
1684 case EdgeOutMorphology:
1685 case EdgeInMorphology:
1686 case TopHatMorphology:
1687 case BottomHatMorphology:
1688 CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1691 case EdgeMorphology: /* subtract erode from dialate ??? */
1692 CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1694 grad_image=DestroyImage(grad_image);
1704 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1708 + R o t a t e K e r n e l I n f o %
1712 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1714 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
1715 % restricted to 90 degree angles, but this may be improved in the future.
1717 % The format of the RotateKernelInfo method is:
1719 % void RotateKernelInfo(KernelInfo *kernel, double angle)
1721 % A description of each parameter follows:
1723 % o kernel: the Morphology/Convolution kernel
1725 % o angle: angle to rotate in degrees
1727 % This function is only internel to this module, as it is not finalized,
1728 % especially with regard to non-orthogonal angles, and rotation of larger
1731 static void RotateKernelInfo(KernelInfo *kernel, double angle)
1733 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
1735 ** TODO: expand beyond simple 90 degree rotates, flips and flops
1738 /* Modulus the angle */
1739 angle = fmod(angle, 360.0);
1743 if ( 315.0 < angle || angle <= 45.0 )
1744 return; /* no change! - At least at this time */
1746 switch (kernel->type) {
1747 /* These built-in kernels are cylindrical kernels, rotating is useless */
1748 case GaussianKernel:
1749 case LaplacianKernel:
1753 case ChebyshevKernel:
1754 case ManhattenKernel:
1755 case EuclideanKernel:
1758 /* These may be rotatable at non-90 angles in the future */
1759 /* but simply rotating them in multiples of 90 degrees is useless */
1765 /* These only allows a +/-90 degree rotation (by transpose) */
1766 /* A 180 degree rotation is useless */
1768 case RectangleKernel:
1769 if ( 135.0 < angle && angle <= 225.0 )
1771 if ( 225.0 < angle && angle <= 315.0 )
1775 /* these are freely rotatable in 90 degree units */
1777 case UndefinedKernel:
1778 case UserDefinedKernel:
1781 if ( 135.0 < angle && angle <= 225.0 )
1783 /* For a 180 degree rotation - also know as a reflection */
1784 /* This is actually a very very common operation! */
1785 /* Basically all that is needed is a reversal of the kernel data! */
1792 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
1793 t=k[i], k[i]=k[j], k[j]=t;
1795 kernel->x = (long) kernel->width - kernel->x - 1;
1796 kernel->y = (long) kernel->width - kernel->y - 1;
1797 angle = fmod(angle+180.0, 360.0);
1799 if ( 45.0 < angle && angle <= 135.0 )
1800 { /* Do a transpose and a flop, of the image, which results in a 90
1801 * degree rotation using two mirror operations.
1803 * WARNING: this assumes the original image was a 1 dimentional image
1804 * but currently that is the only built-ins it is applied to.
1808 t = (long) kernel->width;
1809 kernel->width = kernel->height;
1810 kernel->height = (unsigned long) t;
1812 kernel->x = kernel->y;
1814 angle = fmod(450.0 - angle, 360.0);
1816 /* At this point angle should be between -45 (315) and +45 degrees
1817 * In the future some form of non-orthogonal angled rotates could be
1818 * performed here, posibily with a linear kernel restriction.
1822 Not currently in use!
1823 { /* Do a flop, this assumes kernel is horizontally symetrical.
1824 * Each row of the kernel needs to be reversed!
1833 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
1834 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
1835 t=k[x], k[x]=k[r], k[r]=t;
1837 kernel->x = kernel->width - kernel->x - 1;
1838 angle = fmod(angle+180.0, 360.0);
1845 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1849 + S c a l e K e r n e l I n f o %
1853 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1855 % ScaleKernelInfo() scales the kernel by the given amount. Scaling by value
1856 % of zero will result in a normalization of the kernel.
1858 % For positive kernels normalization scales the kernel so the addition os all
1859 % values is 1.0. While for kernels where values add to zero it is scaled
1860 % so that the convolution output range covers 1.0. In such 'zero kernels'
1861 % it is generally recomended that the user also provides a 50% bias to the
1864 % Correct normalization assumes the 'range_*' attributes of the kernel
1865 % structure have been correctly set during the kernel creation.
1867 % The format of the ScaleKernelInfo method is:
1869 % void ScaleKernelInfo(KernelInfo *kernel)
1871 % A description of each parameter follows:
1873 % o kernel: the Morphology/Convolution kernel
1875 % o scale: multiple all values by this, if zero normalize instead.
1877 % This function is internal to this module only at this time, but can be
1878 % exported to other modules if needed.
1880 static void ScaleKernelInfo(KernelInfo *kernel, double scale)
1885 if ( fabs(scale) < MagickEpsilon ) {
1886 if ( fabs(kernel->positive_range + kernel->negative_range) < MagickEpsilon )
1887 scale = 1/(kernel->positive_range - kernel->negative_range); /* zero kernels */
1889 scale = 1/(kernel->positive_range + kernel->negative_range); /* non-zero kernel */
1892 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1893 if ( ! IsNan(kernel->values[i]) )
1894 kernel->values[i] *= scale;
1896 kernel->positive_range *= scale; /* convolution output range */
1897 kernel->negative_range *= scale;
1898 kernel->maximum *= scale; /* maximum and minimum values in kernel */
1899 kernel->minimum *= scale;
1905 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1909 + S h o w K e r n e l I n f o %
1913 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1915 % ShowKernelInfo() outputs the details of the given kernel defination to
1916 % standard error, generally due to a users 'showkernel' option request.
1918 % The format of the ShowKernel method is:
1920 % void ShowKernelInfo(KernelInfo *kernel)
1922 % A description of each parameter follows:
1924 % o kernel: the Morphology/Convolution kernel
1926 % This function is internal to this module only at this time. That may change
1929 MagickExport void ShowKernelInfo(KernelInfo *kernel)
1935 "Kernel \"%s\" of size %lux%lu%+ld%+ld with values from %.*lg to %.*lg\n",
1936 MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
1937 kernel->width, kernel->height,
1938 kernel->x, kernel->y,
1939 GetMagickPrecision(), kernel->minimum,
1940 GetMagickPrecision(), kernel->maximum);
1941 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
1942 GetMagickPrecision(), kernel->negative_range,
1943 GetMagickPrecision(), kernel->positive_range,
1944 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
1945 for (i=v=0; v < (long) kernel->height; v++) {
1946 fprintf(stderr,"%2ld:",v);
1947 for (u=0; u < (long) kernel->width; u++, i++)
1948 if ( IsNan(kernel->values[i]) )
1949 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
1951 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
1952 GetMagickPrecision(), kernel->values[i]);
1953 fprintf(stderr,"\n");
1958 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1962 + Z e r o K e r n e l N a n s %
1966 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1968 % ZeroKernelNans() replaces any special 'nan' value that may be present in
1969 % the kernel with a zero value. This is typically done when the kernel will
1970 % be used in special hardware (GPU) convolution processors, to simply
1973 % The format of the ZeroKernelNans method is:
1975 % voidZeroKernelNans (KernelInfo *kernel)
1977 % A description of each parameter follows:
1979 % o kernel: the Morphology/Convolution kernel
1981 % FUTURE: return the information in a string for API usage.
1983 MagickExport void ZeroKernelNans(KernelInfo *kernel)
1988 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1989 if ( IsNan(kernel->values[i]) )
1990 kernel->values[i] = 0.0;