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 RotateKernel(KernelInfo *, double),
112 ScaleKernel(KernelInfo *, double),
113 ShowKernel(KernelInfo *);
117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 % A c q u i r e K e r n e l I n f o %
125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127 % AcquireKernelInfo() takes the given string (generally supplied by the
128 % user) and converts it into a Morphology/Convolution Kernel. This allows
129 % users to specify a kernel from a number of pre-defined kernels, or to fully
130 % specify their own kernel for a specific Convolution or Morphology
133 % The kernel so generated can be any rectangular array of floating point
134 % values (doubles) with the 'control point' or 'pixel being affected'
135 % anywhere within that array of values.
137 % Previously IM was restricted to a square of odd size using the exact
138 % center as origin, this is no longer the case, and any rectangular kernel
139 % with any value being declared the origin. This in turn allows the use of
140 % highly asymmetrical kernels.
142 % The floating point values in the kernel can also include a special value
143 % known as 'nan' or 'not a number' to indicate that this value is not part
144 % of the kernel array. This allows you to shaped the kernel within its
145 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
146 % shape. However at least one non-nan value must be provided for correct
147 % working of a kernel.
149 % The returned kernel should be free using the DestroyKernelInfo() when you
150 % are finished with it.
152 % Input kernel defintion strings can consist of any of three types.
155 % Select from one of the built in kernels, using the name and
156 % geometry arguments supplied. See AcquireKernelBuiltIn()
158 % "WxH[+X+Y]:num, num, num ..."
159 % a kernal of size W by H, with W*H floating point numbers following.
160 % the 'center' can be optionally be defined at +X+Y (such that +0+0
161 % is top left corner). If not defined the pixel in the center, for
162 % odd sizes, or to the immediate top or left of center for even sizes
163 % is automatically selected.
165 % "num, num, num, num, ..."
166 % list of floating point numbers defining an 'old style' odd sized
167 % square kernel. At least 9 values should be provided for a 3x3
168 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
169 % Values can be space or comma separated. This is not recommended.
171 % Note that 'name' kernels will start with an alphabetic character while the
172 % new kernel specification has a ':' character in its specification string.
173 % If neither is the case, it is assumed an old style of a simple list of
174 % numbers generating a odd-sized square kernel has been given.
176 % The format of the AcquireKernal method is:
178 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
180 % A description of each parameter follows:
182 % o kernel_string: the Morphology/Convolution kernel wanted.
186 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
192 token[MaxTextExtent];
207 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
209 assert(kernel_string != (const char *) NULL);
210 SetGeometryInfo(&args);
212 /* does it start with an alpha - Return a builtin kernel */
213 GetMagickToken(kernel_string,&p,token);
214 if ( isalpha((int)token[0]) )
219 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
220 if ( type < 0 || type == UserDefinedKernel )
221 return((KernelInfo *)NULL);
223 while (((isspace((int) ((unsigned char) *p)) != 0) ||
224 (*p == ',') || (*p == ':' )) && (*p != '\0'))
226 flags = ParseGeometry(p, &args);
228 /* special handling of missing values in input string */
229 if ( type == RectangleKernel ) {
230 if ( (flags & WidthValue) == 0 ) /* if no width then */
231 args.rho = args.sigma; /* then width = height */
232 if ( args.rho < 1.0 ) /* if width too small */
233 args.rho = 3; /* then width = 3 */
234 if ( args.sigma < 1.0 ) /* if height too small */
235 args.sigma = args.rho; /* then height = width */
236 if ( (flags & XValue) == 0 ) /* center offset if not defined */
237 args.xi = (double)(((long)args.rho-1)/2);
238 if ( (flags & YValue) == 0 )
239 args.psi = (double)(((long)args.sigma-1)/2);
242 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
245 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
246 if (kernel == (KernelInfo *)NULL)
248 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
249 kernel->type = UserDefinedKernel;
250 kernel->signature = MagickSignature;
252 /* Has a ':' in argument - New user kernel specification */
253 p = strchr(kernel_string, ':');
254 if ( p != (char *) NULL)
256 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
257 memcpy(token, kernel_string, (size_t) (p-kernel_string));
258 token[p-kernel_string] = '\0';
259 flags = ParseGeometry(token, &args);
261 /* Size handling and checks of geometry settings */
262 if ( (flags & WidthValue) == 0 ) /* if no width then */
263 args.rho = args.sigma; /* then width = height */
264 if ( args.rho < 1.0 ) /* if width too small */
265 args.rho = 1.0; /* then width = 1 */
266 if ( args.sigma < 1.0 ) /* if height too small */
267 args.sigma = args.rho; /* then height = width */
268 kernel->width = (unsigned long)args.rho;
269 kernel->height = (unsigned long)args.sigma;
271 /* Offset Handling and Checks */
272 if ( args.xi < 0.0 || args.psi < 0.0 )
273 return(DestroyKernelInfo(kernel));
274 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
275 : (long) (kernel->width-1)/2;
276 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
277 : (long) (kernel->height-1)/2;
278 if ( kernel->x >= (long) kernel->width ||
279 kernel->y >= (long) kernel->height )
280 return(DestroyKernelInfo(kernel));
282 p++; /* advance beyond the ':' */
285 { /* ELSE - Old old kernel specification, forming odd-square kernel */
286 /* count up number of values given */
287 p=(const char *) kernel_string;
288 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
289 p++; /* ignore "'" chars for convolve filter usage - Cristy */
290 for (i=0; *p != '\0'; i++)
292 GetMagickToken(p,&p,token);
294 GetMagickToken(p,&p,token);
296 /* set the size of the kernel - old sized square */
297 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
298 kernel->x = kernel->y = (long) (kernel->width-1)/2;
299 p=(const char *) kernel_string;
300 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
301 p++; /* ignore "'" chars for convolve filter usage - Cristy */
304 /* Read in the kernel values from rest of input string argument */
305 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
306 kernel->height*sizeof(double));
307 if (kernel->values == (double *) NULL)
308 return(DestroyKernelInfo(kernel));
310 kernel->minimum = +MagickHuge;
311 kernel->maximum = -MagickHuge;
312 kernel->negative_range = kernel->positive_range = 0.0;
313 for (i=0; (i < (long) (kernel->width*kernel->height)) && (*p != '\0'); i++)
315 GetMagickToken(p,&p,token);
317 GetMagickToken(p,&p,token);
318 if ( LocaleCompare("nan",token) == 0
319 || LocaleCompare("-",token) == 0 ) {
320 kernel->values[i] = nan; /* do not include this value in kernel */
323 kernel->values[i] = StringToDouble(token);
324 ( kernel->values[i] < 0)
325 ? ( kernel->negative_range += kernel->values[i] )
326 : ( kernel->positive_range += kernel->values[i] );
327 Minimize(kernel->minimum, kernel->values[i]);
328 Maximize(kernel->maximum, kernel->values[i]);
331 /* check that we recieved at least one real (non-nan) value! */
332 if ( kernel->minimum == MagickHuge )
333 return(DestroyKernelInfo(kernel));
335 /* This should not be needed for a fully defined kernel
336 * Perhaps an error should be reported instead!
337 * Kept for backward compatibility.
339 if ( i < (long) (kernel->width*kernel->height) ) {
340 Minimize(kernel->minimum, kernel->values[i]);
341 Maximize(kernel->maximum, kernel->values[i]);
342 for ( ; i < (long) (kernel->width*kernel->height); i++)
343 kernel->values[i]=0.0;
350 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
354 % A c q u i r e K e r n e l B u i l t I n %
358 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
360 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
361 % kernels used for special purposes such as gaussian blurring, skeleton
362 % pruning, and edge distance determination.
364 % They take a KernelType, and a set of geometry style arguments, which were
365 % typically decoded from a user supplied string, or from a more complex
366 % Morphology Method that was requested.
368 % The format of the AcquireKernalBuiltIn method is:
370 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
371 % const GeometryInfo args)
373 % A description of each parameter follows:
375 % o type: the pre-defined type of kernel wanted
377 % o args: arguments defining or modifying the kernel
379 % Convolution Kernels
381 % Gaussian "[{radius}]x{sigma}"
382 % Generate a two-dimentional gaussian kernel, as used by -gaussian
383 % A sigma is required, (with the 'x'), due to historical reasons.
385 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
386 % the final size of the resulting kernel to a square 2*radius+1 in size.
387 % The radius should be at least 2 times that of the sigma value, or
388 % sever clipping and aliasing may result. If not given or set to 0 the
389 % radius will be determined so as to produce the best minimal error
390 % result, which is usally much larger than is normally needed.
392 % Blur "[{radius}]x{sigma}[+angle]"
393 % As per Gaussian, but generates a 1 dimensional or linear gaussian
394 % blur, at the angle given (current restricted to orthogonal angles).
395 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
397 % NOTE that two such blurs perpendicular to each other is equivelent to
398 % -blur and the previous gaussian, but is often 10 or more times faster.
400 % Comet "[{width}]x{sigma}[+angle]"
401 % Blur in one direction only, mush like how a bright object leaves
402 % a comet like trail. The Kernel is actually half a gaussian curve,
403 % Adding two such blurs in oppiste directions produces a Linear Blur.
405 % NOTE: that the first argument is the width of the kernel and not the
406 % radius of the kernel.
408 % # Still to be implemented...
410 % # Laplacian "{radius}x{sigma}"
411 % # Laplacian (a mexican hat like) Function
413 % # LOG "{radius},{sigma1},{sigma2}
414 % # Laplacian of Gaussian
416 % # DOG "{radius},{sigma1},{sigma2}
417 % # Difference of Gaussians
421 % Rectangle "{geometry}"
422 % Simply generate a rectangle of 1's with the size given. You can also
423 % specify the location of the 'control point', otherwise the closest
424 % pixel to the center of the rectangle is selected.
426 % Properly centered and odd sized rectangles work the best.
428 % Diamond "[{radius}]"
429 % Generate a diamond shaped kernal with given radius to the points.
430 % Kernel size will again be radius*2+1 square and defaults to radius 1,
431 % generating a 3x3 kernel that is slightly larger than a square.
433 % Square "[{radius}]"
434 % Generate a square shaped kernel of size radius*2+1, and defaulting
435 % to a 3x3 (radius 1).
437 % Note that using a larger radius for the "Square" or the "Diamond"
438 % is also equivelent to iterating the basic morphological method
439 % that many times. However However iterating with the smaller radius 1
440 % default is actually faster than using a larger kernel radius.
443 % Generate a binary disk of the radius given, radius may be a float.
444 % Kernel size will be ceil(radius)*2+1 square.
445 % NOTE: Here are some disk shapes of specific interest
446 % "disk:1" => "diamond" or "cross:1"
447 % "disk:1.5" => "square"
448 % "disk:2" => "diamond:2"
449 % "disk:2.5" => a general disk shape of radius 2
450 % "disk:2.9" => "square:2"
451 % "disk:3.5" => default - octagonal/disk shape of radius 3
452 % "disk:4.2" => roughly octagonal shape of radius 4
453 % "disk:4.3" => a general disk shape of radius 4
454 % After this all the kernel shape becomes more and more circular.
456 % Because a "disk" is more circular when using a larger radius, using a
457 % larger radius is preferred over iterating the morphological operation.
460 % Generate a kernel in the shape of a 'plus' sign. The length of each
461 % arm is also the radius, which defaults to 2.
463 % This kernel is not a good general morphological kernel, but is used
464 % more for highlighting and marking any single pixels in an image using,
465 % a "Dilate" or "Erode" method as appropriate.
467 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
469 % Note that unlike other kernels iterating a plus does not produce the
470 % same result as using a larger radius for the cross.
472 % Distance Measuring Kernels
474 % Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1)
475 % Manhatten "[{radius}][x{scale}]" square grid distance (default r=1)
476 % Euclidean "[{radius}][x{scale}]" direct distance (default r=1)
478 % Different types of distance measuring methods, which are used with the
479 % a 'Distance' morphology method for generating a gradient based on
480 % distance from an edge of a binary shape, though there is a technique
481 % for handling a anti-aliased shape.
483 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
484 % one to any neighbour, orthogonal or diagonal. One why of thinking of
485 % it is the number of squares a 'King' or 'Queen' in chess needs to
486 % traverse reach any other position on a chess board. It results in a
487 % 'square' like distance function, but one where diagonals are closer
490 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
491 % Cab metric), is the distance needed when you can only travel in
492 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
493 % in chess would travel. It results in a diamond like distances, where
494 % diagonals are further than expected.
496 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
497 % However by default the kernel size only has a radius of 1, which
498 % limits the distance to 'Knight' like moves, with only orthogonal and
499 % diagonal measurements being correct. As such for the default kernel
500 % you will get octagonal like distance function, which is reasonally
503 % However if you use a larger radius such as "Euclidean:4" you will
504 % get a much smoother distance gradient from the edge of the shape.
505 % Of course a larger kernel is slower to use, and generally not needed.
507 % To allow the use of fractional distances that you get with diagonals
508 % the actual distance is scaled by a fixed value which the user can
509 % provide. This is not actually nessary for either ""Chebyshev" or
510 % "Manhatten" distance kernels, but is done for all three distance
511 % kernels. If no scale is provided it is set to a value of 100,
512 % allowing for a maximum distance measurement of 655 pixels using a Q16
513 % version of IM, from any edge. However for small images this can
514 % result in quite a dark gradient.
516 % See the 'Distance' Morphological Method, for information of how it is
521 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
522 const GeometryInfo *args)
535 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
537 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
538 if (kernel == (KernelInfo *) NULL)
540 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
541 kernel->minimum = kernel->maximum = 0.0;
542 kernel->negative_range = kernel->positive_range = 0.0;
544 kernel->signature = MagickSignature;
547 /* Convolution Kernels */
550 sigma = fabs(args->sigma);
552 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
554 kernel->width = kernel->height =
555 GetOptimalKernelWidth2D(args->rho,sigma);
556 kernel->x = kernel->y = (long) (kernel->width-1)/2;
557 kernel->negative_range = kernel->positive_range = 0.0;
558 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
559 kernel->height*sizeof(double));
560 if (kernel->values == (double *) NULL)
561 return(DestroyKernelInfo(kernel));
563 sigma = 2.0*sigma*sigma; /* simplify the expression */
564 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
565 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
566 kernel->positive_range += (
568 exp(-((double)(u*u+v*v))/sigma)
569 /* / (MagickPI*sigma) */ );
571 kernel->maximum = kernel->values[
572 kernel->y*kernel->width+kernel->x ];
574 ScaleKernel(kernel, 0.0); /* Normalize Kernel Values */
580 sigma = fabs(args->sigma);
582 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
584 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
585 kernel->x = (long) (kernel->width-1)/2;
588 kernel->negative_range = kernel->positive_range = 0.0;
589 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
590 kernel->height*sizeof(double));
591 if (kernel->values == (double *) NULL)
592 return(DestroyKernelInfo(kernel));
596 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
597 ** It generates a gaussian 3 times the width, and compresses it into
598 ** the expected range. This produces a closer normalization of the
599 ** resulting kernel, especially for very low sigma values.
600 ** As such while wierd it is prefered.
602 ** I am told this method originally came from Photoshop.
604 sigma *= KernelRank; /* simplify expanded curve */
605 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
606 (void) ResetMagickMemory(kernel->values,0, (size_t)
607 kernel->width*sizeof(double));
608 for ( u=-v; u <= v; u++) {
609 kernel->values[(u+v)/KernelRank] +=
610 exp(-((double)(u*u))/(2.0*sigma*sigma))
611 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
613 for (i=0; i < (long) kernel->width; i++)
614 kernel->positive_range += kernel->values[i];
616 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
617 kernel->positive_range += (
619 exp(-((double)(u*u))/(2.0*sigma*sigma))
620 /* / (MagickSQ2PI*sigma) */ );
623 kernel->maximum = kernel->values[ kernel->x ];
624 /* Note that neither methods above generate a normalized kernel,
625 ** though it gets close. The kernel may be 'clipped' by a user defined
626 ** radius, producing a smaller (darker) kernel. Also for very small
627 ** sigma's (> 0.1) the central value becomes larger than one, and thus
628 ** producing a very bright kernel.
631 /* Normalize the 1D Gaussian Kernel
633 ** Because of this the divisor in the above kernel generator is
634 ** not needed, so is not done above.
636 ScaleKernel(kernel, 0.0); /* Normalize Kernel Values */
638 /* rotate the kernel by given angle */
639 RotateKernel(kernel, args->xi);
644 sigma = fabs(args->sigma);
646 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
648 if ( args->rho < 1.0 )
649 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
651 kernel->width = (unsigned long)args->rho;
652 kernel->x = kernel->y = 0;
654 kernel->negative_range = kernel->positive_range = 0.0;
655 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
656 kernel->height*sizeof(double));
657 if (kernel->values == (double *) NULL)
658 return(DestroyKernelInfo(kernel));
660 /* A comet blur is half a gaussian curve, so that the object is
661 ** blurred in one direction only. This may not be quite the right
662 ** curve so may change in the future. The function must be normalised.
666 sigma *= KernelRank; /* simplify expanded curve */
667 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
668 (void) ResetMagickMemory(kernel->values,0, (size_t)
669 kernel->width*sizeof(double));
670 for ( u=0; u < v; u++) {
671 kernel->values[u/KernelRank] +=
672 exp(-((double)(u*u))/(2.0*sigma*sigma))
673 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
675 for (i=0; i < (long) kernel->width; i++)
676 kernel->positive_range += kernel->values[i];
678 for ( i=0; i < (long) kernel->width; i++)
679 kernel->positive_range += (
681 exp(-((double)(i*i))/(2.0*sigma*sigma))
682 /* / (MagickSQ2PI*sigma) */ );
685 kernel->maximum = kernel->values[0];
687 ScaleKernel(kernel, 0.0); /* Normalize Kernel Values */
688 RotateKernel(kernel, args->xi);
691 /* Boolean Kernels */
692 case RectangleKernel:
695 if ( type == SquareKernel )
698 kernel->width = kernel->height = 3; /* default radius = 1 */
700 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
701 kernel->x = kernel->y = (long) (kernel->width-1)/2;
704 /* NOTE: user defaults set in "AcquireKernelInfo()" */
705 if ( args->rho < 1.0 || args->sigma < 1.0 )
706 return(DestroyKernelInfo(kernel)); /* invalid args given */
707 kernel->width = (unsigned long)args->rho;
708 kernel->height = (unsigned long)args->sigma;
709 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
710 args->psi < 0.0 || args->psi > (double)kernel->height )
711 return(DestroyKernelInfo(kernel)); /* invalid args given */
712 kernel->x = (long) args->xi;
713 kernel->y = (long) args->psi;
715 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
716 kernel->height*sizeof(double));
717 if (kernel->values == (double *) NULL)
718 return(DestroyKernelInfo(kernel));
720 /* set all kernel values to 1.0 */
721 u=(long) kernel->width*kernel->height;
722 for ( i=0; i < u; i++)
723 kernel->values[i] = 1.0;
724 kernel->minimum = kernel->maximum = 1.0; /* a flat shape */
725 kernel->positive_range = (double) u;
731 kernel->width = kernel->height = 3; /* default radius = 1 */
733 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
734 kernel->x = kernel->y = (long) (kernel->width-1)/2;
736 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
737 kernel->height*sizeof(double));
738 if (kernel->values == (double *) NULL)
739 return(DestroyKernelInfo(kernel));
741 /* set all kernel values within diamond area to 1.0 */
742 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
743 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
744 if ((labs(u)+labs(v)) <= (long)kernel->x)
745 kernel->positive_range += kernel->values[i] = 1.0;
747 kernel->values[i] = nan;
748 kernel->minimum = kernel->maximum = 1.0; /* a flat shape */
756 limit = (long)(args->rho*args->rho);
757 if (args->rho < 0.1) /* default radius approx 3.5 */
758 kernel->width = kernel->height = 7L, limit = 10L;
760 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
761 kernel->x = kernel->y = (long) (kernel->width-1)/2;
763 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
764 kernel->height*sizeof(double));
765 if (kernel->values == (double *) NULL)
766 return(DestroyKernelInfo(kernel));
768 /* set all kernel values within disk area to 1.0 */
769 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
770 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
771 if ((u*u+v*v) <= limit)
772 kernel->positive_range += kernel->values[i] = 1.0;
774 kernel->values[i] = nan;
775 kernel->minimum = kernel->maximum = 1.0; /* a flat shape */
781 kernel->width = kernel->height = 5; /* default radius 2 */
783 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
784 kernel->x = kernel->y = (long) (kernel->width-1)/2;
786 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
787 kernel->height*sizeof(double));
788 if (kernel->values == (double *) NULL)
789 return(DestroyKernelInfo(kernel));
791 /* set all kernel values along axises to 1.0 */
792 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
793 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
794 kernel->values[i] = (u == 0 || v == 0) ? 1.0 : nan;
795 kernel->minimum = kernel->maximum = 1.0; /* a flat shape */
796 kernel->positive_range = (double) kernel->width*2.0 - 1.0;
799 /* Distance Measuring Kernels */
800 case ChebyshevKernel:
806 kernel->width = kernel->height = 3; /* default radius = 1 */
808 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
809 kernel->x = kernel->y = (long) (kernel->width-1)/2;
811 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
812 kernel->height*sizeof(double));
813 if (kernel->values == (double *) NULL)
814 return(DestroyKernelInfo(kernel));
816 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
817 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
818 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
819 kernel->positive_range += ( kernel->values[i] =
820 scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
821 kernel->maximum = kernel->values[0];
824 case ManhattenKernel:
830 kernel->width = kernel->height = 3; /* default radius = 1 */
832 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
833 kernel->x = kernel->y = (long) (kernel->width-1)/2;
835 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
836 kernel->height*sizeof(double));
837 if (kernel->values == (double *) NULL)
838 return(DestroyKernelInfo(kernel));
840 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
841 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
842 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
843 kernel->positive_range += ( kernel->values[i] =
844 scale*(labs(u)+labs(v)) );
845 kernel->maximum = kernel->values[0];
848 case EuclideanKernel:
854 kernel->width = kernel->height = 3; /* default radius = 1 */
856 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
857 kernel->x = kernel->y = (long) (kernel->width-1)/2;
859 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
860 kernel->height*sizeof(double));
861 if (kernel->values == (double *) NULL)
862 return(DestroyKernelInfo(kernel));
864 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
865 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
866 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
867 kernel->positive_range += ( kernel->values[i] =
868 scale*sqrt((double)(u*u+v*v)) );
869 kernel->maximum = kernel->values[0];
872 /* Undefined Kernels */
873 case LaplacianKernel:
876 perror("Kernel Type has not been defined yet");
879 /* Generate a No-Op minimal kernel - 1x1 pixel */
880 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
881 if (kernel->values == (double *) NULL)
882 return(DestroyKernelInfo(kernel));
883 kernel->width = kernel->height = 1;
884 kernel->x = kernel->x = 0;
885 kernel->type = UndefinedKernel;
887 kernel->positive_range =
888 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
896 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
900 % D e s t r o y K e r n e l I n f o %
904 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
906 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
909 % The format of the DestroyKernelInfo method is:
911 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
913 % A description of each parameter follows:
915 % o kernel: the Morphology/Convolution kernel to be destroyed
919 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
921 assert(kernel != (KernelInfo *) NULL);
922 kernel->values=(double *)RelinquishMagickMemory(kernel->values);
923 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
928 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
932 % M o r p h o l o g y I m a g e C h a n n e l %
936 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
938 % MorphologyImageChannel() applies a user supplied kernel to the image
939 % according to the given mophology method.
941 % The given kernel is assumed to have been pre-scaled appropriatally, usally
942 % by the kernel generator.
944 % The format of the MorphologyImage method is:
946 % Image *MorphologyImage(const Image *image, MorphologyMethod method,
947 % const long iterations, KernelInfo *kernel, ExceptionInfo *exception)
948 % Image *MorphologyImageChannel(const Image *image, const ChannelType
949 % channel, MorphologyMethod method, const long iterations, KernelInfo
950 % *kernel, ExceptionInfo *exception)
952 % A description of each parameter follows:
954 % o image: the image.
956 % o method: the morphology method to be applied.
958 % o iterations: apply the operation this many times (or no change).
959 % A value of -1 means loop until no change found.
960 % How this is applied may depend on the morphology method.
961 % Typically this is a value of 1.
963 % o channel: the channel type.
965 % o kernel: An array of double representing the morphology kernel.
966 % Warning: kernel may be normalized for the Convolve method.
968 % o exception: return any errors or warnings in this structure.
971 % TODO: bias and auto-scale handling of the kernel for convolution
972 % The given kernel is assumed to have been pre-scaled appropriatally, usally
973 % by the kernel generator.
978 * Apply the Morphology method with the given Kernel
979 * And return the number of pixels that changed.
981 static unsigned long MorphologyApply(const Image *image, Image
982 *result_image, const MorphologyMethod method, const ChannelType channel,
983 const KernelInfo *kernel, ExceptionInfo *exception)
985 #define MorphologyTag "Morphology/Image"
1003 Apply Morphology to Image.
1009 GetMagickPixelPacket(image,&bias);
1010 SetMagickPixelPacketBias(image,&bias);
1011 /* Future: handle auto-bias from user, based on kernel input */
1013 p_view=AcquireCacheView(image);
1014 q_view=AcquireCacheView(result_image);
1016 /* Some methods (including convolve) needs use a reflected kernel.
1017 * Adjust 'origin' offsets for this reflected kernel.
1022 case ErodeMorphology:
1023 case ErodeIntensityMorphology:
1024 /* kernel is not reflected */
1027 /* kernel needs to be reflected */
1028 offx = (long) kernel->width-offx-1;
1029 offy = (long) kernel->height-offy-1;
1033 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1034 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1036 for (y=0; y < (long) image->rows; y++)
1041 register const PixelPacket
1044 register const IndexPacket
1045 *restrict p_indexes;
1047 register PixelPacket
1050 register IndexPacket
1051 *restrict q_indexes;
1059 if (status == MagickFalse)
1061 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1062 image->columns+kernel->width, kernel->height, exception);
1063 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1065 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1070 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1071 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1072 r = (image->columns+kernel->width)*offy+offx; /* constant */
1074 for (x=0; x < (long) image->columns; x++)
1082 register const double
1085 register const PixelPacket
1088 register const IndexPacket
1089 *restrict k_indexes;
1094 /* Copy input to ouput image for unused channels
1095 * This removes need for 'cloning' a new image every iteration
1098 if (image->colorspace == CMYKColorspace)
1099 q_indexes[x] = p_indexes[r];
1101 result.index=(MagickRealType) 0; /* stop compiler warnings */
1103 case ConvolveMorphology:
1105 break; /* default result is the convolution bias */
1107 case DilateMorphology:
1112 result.index = -MagickHuge;
1114 case ErodeMorphology:
1119 result.index = +MagickHuge;
1123 /* Otherwise just start with the original pixel value */
1124 result.red = (MagickRealType) p[r].red;
1125 result.green = (MagickRealType) p[r].green;
1126 result.blue = (MagickRealType) p[r].blue;
1127 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1128 if ( image->colorspace == CMYKColorspace)
1129 result.index = (MagickRealType) p_indexes[r];
1134 case ConvolveMorphology:
1135 /* Weighted Average of pixels
1137 * NOTE for correct working of this operation for asymetrical
1138 * kernels, the kernel needs to be applied in its reflected form.
1139 * That is its values needs to be reversed.
1141 if (((channel & OpacityChannel) == 0) ||
1142 (image->matte == MagickFalse))
1144 /* Convolution (no transparency) */
1145 k = &kernel->values[ kernel->width*kernel->height-1 ];
1147 k_indexes = p_indexes;
1148 for (v=0; v < (long) kernel->height; v++) {
1149 for (u=0; u < (long) kernel->width; u++, k--) {
1150 if ( IsNan(*k) ) continue;
1151 result.red += (*k)*k_pixels[u].red;
1152 result.green += (*k)*k_pixels[u].green;
1153 result.blue += (*k)*k_pixels[u].blue;
1154 /* result.opacity += not involved here */
1155 if ( image->colorspace == CMYKColorspace)
1156 result.index += (*k)*k_indexes[u];
1158 k_pixels += image->columns+kernel->width;
1159 k_indexes += image->columns+kernel->width;
1163 { /* Kernel & Alpha weighted Convolution */
1165 alpha, /* alpha value * kernel weighting */
1166 gamma; /* weighting divisor */
1169 k = &kernel->values[ kernel->width*kernel->height-1 ];
1171 k_indexes = p_indexes;
1172 for (v=0; v < (long) kernel->height; v++) {
1173 for (u=0; u < (long) kernel->width; u++, k--) {
1174 if ( IsNan(*k) ) continue;
1175 alpha=(*k)*(QuantumScale*(QuantumRange-
1176 k_pixels[u].opacity));
1178 result.red += alpha*k_pixels[u].red;
1179 result.green += alpha*k_pixels[u].green;
1180 result.blue += alpha*k_pixels[u].blue;
1181 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1182 if ( image->colorspace == CMYKColorspace)
1183 result.index += alpha*k_indexes[u];
1185 k_pixels += image->columns+kernel->width;
1186 k_indexes += image->columns+kernel->width;
1188 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1189 result.red *= gamma;
1190 result.green *= gamma;
1191 result.blue *= gamma;
1192 result.opacity *= gamma;
1193 result.index *= gamma;
1197 case DilateMorphology:
1198 /* Maximize Value within kernel shape
1200 * NOTE for correct working of this operation for asymetrical
1201 * kernels, the kernel needs to be applied in its reflected form.
1202 * That is its values needs to be reversed.
1204 * NOTE: in normal Greyscale Morphology, the kernel value should
1205 * be added to the real value, this is currently not done, due to
1206 * the nature of the boolean kernels being used.
1209 k = &kernel->values[ kernel->width*kernel->height-1 ];
1211 k_indexes = p_indexes;
1212 for (v=0; v < (long) kernel->height; v++) {
1213 for (u=0; u < (long) kernel->width; u++, k--) {
1214 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1215 Maximize(result.red, (double) k_pixels[u].red);
1216 Maximize(result.green, (double) k_pixels[u].green);
1217 Maximize(result.blue, (double) k_pixels[u].blue);
1218 Maximize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1219 if ( image->colorspace == CMYKColorspace)
1220 Maximize(result.index, (double) k_indexes[u]);
1222 k_pixels += image->columns+kernel->width;
1223 k_indexes += image->columns+kernel->width;
1227 case ErodeMorphology:
1228 /* Minimize Value within kernel shape
1230 * NOTE that the kernel is not reflected for this operation!
1232 * NOTE: in normal Greyscale Morphology, the kernel value should
1233 * be added to the real value, this is currently not done, due to
1234 * the nature of the boolean kernels being used.
1238 k_indexes = p_indexes;
1239 for (v=0; v < (long) kernel->height; v++) {
1240 for (u=0; u < (long) kernel->width; u++, k++) {
1241 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1242 Minimize(result.red, (double) k_pixels[u].red);
1243 Minimize(result.green, (double) k_pixels[u].green);
1244 Minimize(result.blue, (double) k_pixels[u].blue);
1245 Minimize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1246 if ( image->colorspace == CMYKColorspace)
1247 Minimize(result.index, (double) k_indexes[u]);
1249 k_pixels += image->columns+kernel->width;
1250 k_indexes += image->columns+kernel->width;
1254 case DilateIntensityMorphology:
1255 /* Select pixel with maximum intensity within kernel shape
1257 * WARNING: the intensity test fails for CMYK and does not
1258 * take into account the moderating effect of teh alpha channel
1261 * NOTE for correct working of this operation for asymetrical
1262 * kernels, the kernel needs to be applied in its reflected form.
1263 * That is its values needs to be reversed.
1265 k = &kernel->values[ kernel->width*kernel->height-1 ];
1267 k_indexes = p_indexes;
1268 for (v=0; v < (long) kernel->height; v++) {
1269 for (u=0; u < (long) kernel->width; u++, k--) {
1270 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1271 if ( result.red == 0.0 ||
1272 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1273 /* copy the whole pixel - no channel selection */
1275 if ( result.red > 0.0 ) changed++;
1279 k_pixels += image->columns+kernel->width;
1280 k_indexes += image->columns+kernel->width;
1284 case ErodeIntensityMorphology:
1285 /* Select pixel with mimimum intensity within kernel shape
1287 * WARNING: the intensity test fails for CMYK and does not
1288 * take into account the moderating effect of teh alpha channel
1291 * NOTE that the kernel is not reflected for this operation!
1295 k_indexes = p_indexes;
1296 for (v=0; v < (long) kernel->height; v++) {
1297 for (u=0; u < (long) kernel->width; u++, k++) {
1298 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1299 if ( result.red == 0.0 ||
1300 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1301 /* copy the whole pixel - no channel selection */
1303 if ( result.red > 0.0 ) changed++;
1307 k_pixels += image->columns+kernel->width;
1308 k_indexes += image->columns+kernel->width;
1312 case DistanceMorphology:
1313 /* Add kernel value and select the minimum value found.
1314 * The result is a iterative distance from edge function.
1316 * All Distance Kernels are symetrical, but that may not always
1317 * be the case. For example how about a distance from left edges?
1318 * To make it work correctly for asymetrical kernels the reflected
1319 * kernel needs to be applied.
1322 /* No need to do distance morphology if original value is zero
1323 * Unfortunatally I have not been able to get this right
1324 * when channel selection also becomes involved. -- Arrgghhh
1326 if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1327 || ((channel & GreenChannel) == 0 && p[r].green == 0)
1328 || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1329 || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1330 || (( (channel & IndexChannel) == 0
1331 || image->colorspace != CMYKColorspace
1332 ) && p_indexes[x] ==0 )
1336 k = &kernel->values[ kernel->width*kernel->height-1 ];
1338 k_indexes = p_indexes;
1339 for (v=0; v < (long) kernel->height; v++) {
1340 for (u=0; u < (long) kernel->width; u++, k--) {
1341 if ( IsNan(*k) ) continue;
1342 Minimize(result.red, (*k)+k_pixels[u].red);
1343 Minimize(result.green, (*k)+k_pixels[u].green);
1344 Minimize(result.blue, (*k)+k_pixels[u].blue);
1345 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1346 if ( image->colorspace == CMYKColorspace)
1347 Minimize(result.index, (*k)+k_indexes[u]);
1349 k_pixels += image->columns+kernel->width;
1350 k_indexes += image->columns+kernel->width;
1354 case UndefinedMorphology:
1356 break; /* Do nothing */
1359 case UndefinedMorphology:
1360 case DilateIntensityMorphology:
1361 case ErodeIntensityMorphology:
1362 break; /* full pixel was directly assigned */
1364 /* Assign the results */
1365 if ((channel & RedChannel) != 0)
1366 q->red = ClampToQuantum(result.red);
1367 if ((channel & GreenChannel) != 0)
1368 q->green = ClampToQuantum(result.green);
1369 if ((channel & BlueChannel) != 0)
1370 q->blue = ClampToQuantum(result.blue);
1371 if ((channel & OpacityChannel) != 0
1372 && image->matte == MagickTrue )
1373 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1374 if ((channel & IndexChannel) != 0
1375 && image->colorspace == CMYKColorspace)
1376 q_indexes[x] = ClampToQuantum(result.index);
1379 if ( ( p[r].red != q->red )
1380 || ( p[r].green != q->green )
1381 || ( p[r].blue != q->blue )
1382 || ( p[r].opacity != q->opacity )
1383 || ( image->colorspace == CMYKColorspace &&
1384 p_indexes[r] != q_indexes[x] ) )
1385 changed++; /* The pixel had some value changed! */
1389 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1390 if (sync == MagickFalse)
1392 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1397 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1398 #pragma omp critical (MagickCore_MorphologyImage)
1400 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1401 if (proceed == MagickFalse)
1405 result_image->type=image->type;
1406 q_view=DestroyCacheView(q_view);
1407 p_view=DestroyCacheView(p_view);
1408 return(status ? (unsigned long) changed : 0);
1411 MagickExport Image *MorphologyImage(const Image *image,MorphologyMethod method,
1412 const long iterations,KernelInfo *kernel, ExceptionInfo *exception)
1417 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1418 iterations,kernel,exception);
1419 return(morphology_image);
1422 MagickExport Image *MorphologyImageChannel(const Image *image,
1423 const ChannelType channel, MorphologyMethod method, const long iterations,
1424 KernelInfo *kernel, ExceptionInfo *exception)
1440 assert(image != (Image *) NULL);
1441 assert(image->signature == MagickSignature);
1442 assert(exception != (ExceptionInfo *) NULL);
1443 assert(exception->signature == MagickSignature);
1445 if ( GetImageArtifact(image,"showkernel") != (const char *) NULL)
1446 ShowKernel(kernel); /* request to display the kernel to stderr */
1448 if ( iterations == 0 )
1449 return((Image *)NULL); /* null operation - nothing to do! */
1451 /* kernel must be valid at this point
1452 * (except maybe for posible future morphology methods like "Prune"
1454 assert(kernel != (KernelInfo *)NULL);
1458 limit = (unsigned long) iterations;
1459 if ( iterations < 0 )
1460 limit = image->columns > image->rows ? image->columns : image->rows;
1462 /* Special morphology cases */
1464 case CloseMorphology:
1465 new_image = MorphologyImageChannel(image, channel, DilateMorphology,
1466 iterations, kernel, exception);
1467 if (new_image == (Image *) NULL)
1468 return((Image *) NULL);
1469 method = ErodeMorphology;
1471 case OpenMorphology:
1472 new_image = MorphologyImageChannel(image, channel, ErodeMorphology,
1473 iterations, kernel, exception);
1474 if (new_image == (Image *) NULL)
1475 return((Image *) NULL);
1476 method = DilateMorphology;
1478 case CloseIntensityMorphology:
1479 new_image = MorphologyImageChannel(image, channel, DilateIntensityMorphology,
1480 iterations, kernel, exception);
1481 if (new_image == (Image *) NULL)
1482 return((Image *) NULL);
1483 method = ErodeIntensityMorphology;
1485 case OpenIntensityMorphology:
1486 new_image = MorphologyImageChannel(image, channel, ErodeIntensityMorphology,
1487 iterations, kernel, exception);
1488 if (new_image == (Image *) NULL)
1489 return((Image *) NULL);
1490 method = DilateIntensityMorphology;
1493 case ConvolveMorphology:
1494 /* Scale or Normalize kernel according to user wishes
1495 ** WARNING: this directly modifies the kernel
1496 ** which probably should not be done.
1498 artifact = GetImageArtifact(image,"convolve:scale");
1499 if ( artifact != (char *)NULL )
1500 ScaleKernel(kernel, StringToDouble(artifact) );
1503 /* Do a morphology iteration just once at this point!
1504 This ensures a new_image has been generated, but allows us
1505 to skip the creation of 'old_image' if it isn't needed.
1507 new_image=CloneImage(image,0,0,MagickTrue,exception);
1508 if (new_image == (Image *) NULL)
1509 return((Image *) NULL);
1510 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1512 InheritException(exception,&new_image->exception);
1513 new_image=DestroyImage(new_image);
1514 return((Image *) NULL);
1516 changed = MorphologyApply(image,new_image,method,channel,kernel,
1519 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1520 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1521 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1525 /* Repeat an interative morphology until count or no change reached */
1526 if ( count < (long) limit && changed > 0 ) {
1527 old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1528 if (old_image == (Image *) NULL)
1529 return(DestroyImage(new_image));
1530 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1532 InheritException(exception,&old_image->exception);
1533 old_image=DestroyImage(old_image);
1534 return(DestroyImage(new_image));
1536 while( count < (long) limit && changed != 0 )
1538 Image *tmp = old_image;
1539 old_image = new_image;
1541 changed = MorphologyApply(old_image,new_image,method,channel,kernel,
1544 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1545 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1546 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1549 old_image=DestroyImage(old_image);
1556 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1560 + R o t a t e K e r n e l %
1564 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1566 % RotateKernel() rotates the kernel by the angle given. Currently it is
1567 % restricted to 90 degree angles, but this may be improved in the future.
1569 % The format of the RotateKernel method is:
1571 % void RotateKernel(KernelInfo *kernel, double angle)
1573 % A description of each parameter follows:
1575 % o kernel: the Morphology/Convolution kernel
1577 % o angle: angle to rotate in degrees
1579 % This function is only internel to this module, as it is not finalized,
1580 % especially with regard to non-orthogonal angles, and rotation of larger
1583 static void RotateKernel(KernelInfo *kernel, double angle)
1585 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
1587 ** TODO: expand beyond simple 90 degree rotates, flips and flops
1590 /* Modulus the angle */
1591 angle = fmod(angle, 360.0);
1595 if ( 315.0 < angle || angle <= 45.0 )
1596 return; /* no change! - At least at this time */
1598 switch (kernel->type) {
1599 /* These built-in kernels are cylindrical kernels, rotating is useless */
1600 case GaussianKernel:
1601 case LaplacianKernel:
1605 case ChebyshevKernel:
1606 case ManhattenKernel:
1607 case EuclideanKernel:
1610 /* These may be rotatable at non-90 angles in the future */
1611 /* but simply rotating them in multiples of 90 degrees is useless */
1617 /* These only allows a +/-90 degree rotation (by transpose) */
1618 /* A 180 degree rotation is useless */
1620 case RectangleKernel:
1621 if ( 135.0 < angle && angle <= 225.0 )
1623 if ( 225.0 < angle && angle <= 315.0 )
1627 /* these are freely rotatable in 90 degree units */
1629 case UndefinedKernel:
1630 case UserDefinedKernel:
1633 if ( 135.0 < angle && angle <= 225.0 )
1635 /* For a 180 degree rotation - also know as a reflection */
1636 /* This is actually a very very common operation! */
1637 /* Basically all that is needed is a reversal of the kernel data! */
1644 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
1645 t=k[i], k[i]=k[j], k[j]=t;
1647 kernel->x = (long) kernel->width - kernel->x - 1;
1648 kernel->y = (long) kernel->width - kernel->y - 1;
1649 angle = fmod(angle+180.0, 360.0);
1651 if ( 45.0 < angle && angle <= 135.0 )
1652 { /* Do a transpose and a flop, of the image, which results in a 90
1653 * degree rotation using two mirror operations.
1655 * WARNING: this assumes the original image was a 1 dimentional image
1656 * but currently that is the only built-ins it is applied to.
1660 t = (long) kernel->width;
1661 kernel->width = kernel->height;
1662 kernel->height = (unsigned long) t;
1664 kernel->x = kernel->y;
1666 angle = fmod(450.0 - angle, 360.0);
1668 /* At this point angle should be between -45 (315) and +45 degrees
1669 * In the future some form of non-orthogonal angled rotates could be
1670 * performed here, posibily with a linear kernel restriction.
1674 Not currently in use!
1675 { /* Do a flop, this assumes kernel is horizontally symetrical.
1676 * Each row of the kernel needs to be reversed!
1685 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
1686 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
1687 t=k[x], k[x]=k[r], k[r]=t;
1689 kernel->x = kernel->width - kernel->x - 1;
1690 angle = fmod(angle+180.0, 360.0);
1697 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1701 + S c a l e K e r n e l %
1705 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1707 % ScaleKernel() scales the kernel by the given amount. Scaling by value of
1708 % zero will result in a normalization of the kernel.
1710 % For positive kernels normalization scales the kernel so the addition os all
1711 % values is 1.0. While for kernels where values add to zero it is scaled
1712 % so that the convolution output range covers 1.0. In such 'zero kernels'
1713 % it is generally recomended that the user also provides a 50% bias to the
1716 % Correct normalization assumes the 'range_*' attributes of the kernel
1717 % structure have been correctly set during the kernel creation.
1719 % The format of the ScaleKernel method is:
1721 % void ScaleKernel(KernelInfo *kernel)
1723 % A description of each parameter follows:
1725 % o kernel: the Morphology/Convolution kernel
1727 % o scale: multiple all values by this, if zero normalize instead.
1729 % This function is internal to this module only at this time, but can be
1730 % exported to other modules if needed.
1732 static void ScaleKernel(KernelInfo *kernel, double scale)
1737 if ( fabs(scale) < MagickEpsilon ) {
1738 if ( fabs(kernel->positive_range + kernel->negative_range) < MagickEpsilon )
1739 scale = 1/(kernel->positive_range - kernel->negative_range); /* zero kernels */
1741 scale = 1/(kernel->positive_range + kernel->negative_range); /* non-zero kernel */
1744 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1745 if ( ! IsNan(kernel->values[i]) )
1746 kernel->values[i] *= scale;
1748 kernel->positive_range *= scale; /* convolution output range */
1749 kernel->negative_range *= scale;
1750 kernel->maximum *= scale; /* maximum and minimum values in kernel */
1751 kernel->minimum *= scale;
1757 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1761 + S h o w K e r n e l %
1765 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1767 % ShowKernel() Output the details of the given kernel defination to
1768 % standard error, as per a users 'showkernel' option request.
1770 % The format of the ShowKernel method is:
1772 % void ShowKernel(KernelInfo *kernel)
1774 % A description of each parameter follows:
1776 % o kernel: the Morphology/Convolution kernel
1778 % This function is internal to this module only at this time. That may change
1781 static void ShowKernel(KernelInfo *kernel)
1787 "Kernel \"%s\" of size %lux%lu%+ld%+ld with values from %.*lg to %.*lg\n",
1788 MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
1789 kernel->width, kernel->height,
1790 kernel->x, kernel->y,
1791 GetMagickPrecision(), kernel->minimum,
1792 GetMagickPrecision(), kernel->maximum);
1793 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
1794 GetMagickPrecision(), kernel->negative_range,
1795 GetMagickPrecision(), kernel->positive_range,
1796 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
1797 for (i=v=0; v < (long) kernel->height; v++) {
1798 fprintf(stderr,"%2ld:",v);
1799 for (u=0; u < (long) kernel->width; u++, i++)
1800 if ( IsNan(kernel->values[i]) )
1801 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
1803 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
1804 GetMagickPrecision(), kernel->values[i]);
1805 fprintf(stderr,"\n");
1810 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1814 + Z e r o K e r n e l N a n s %
1818 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1820 % ZeroKernelNans() replaces any special 'nan' value that may be present in
1821 % the kernel with a zero value. This is typically done when the kernel will
1822 % be used in special hardware (GPU) convolution processors, to simply
1825 % The format of the ZeroKernelNans method is:
1827 % voidZeroKernelNans (KernelInfo *kernel)
1829 % A description of each parameter follows:
1831 % o kernel: the Morphology/Convolution kernel
1833 % FUTURE: return the information in a string for API usage.
1835 MagickExport void ZeroKernelNans(KernelInfo *kernel)
1840 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1841 if ( IsNan(kernel->values[i]) )
1842 kernel->values[i] = 0.0;