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 internal to this module
110 * Eventually these may become 'private' to library
111 * OR may become externally available to API's
113 static MagickExport void
114 NormalizeKernel(KernelInfo *),
115 RotateKernel(KernelInfo *, double),
116 ShowKernel(KernelInfo *);
120 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
124 % A c q u i r e K e r n e l I n f o %
128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
130 % AcquireKernelInfo() takes the given string (generally supplied by the
131 % user) and converts it into a Morphology/Convolution Kernel. This allows
132 % users to specify a kernel from a number of pre-defined kernels, or to fully
133 % specify their own kernel for a specific Convolution or Morphology
136 % The kernel so generated can be any rectangular array of floating point
137 % values (doubles) with the 'control point' or 'pixel being affected'
138 % anywhere within that array of values.
140 % Previously IM was restricted to a square of odd size using the exact
141 % center as origin, this is no longer the case, and any rectangular kernel
142 % with any value being declared the origin. This in turn allows the use of
143 % highly asymmetrical kernels.
145 % The floating point values in the kernel can also include a special value
146 % known as 'nan' or 'not a number' to indicate that this value is not part
147 % of the kernel array. This allows you to shaped the kernel within its
148 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
149 % shape. However at least one non-nan value must be provided for correct
150 % working of a kernel.
152 % The returned kernel should be free using the DestroyKernelInfo() when you
153 % are finished with it.
155 % Input kernel defintion strings can consist of any of three types.
158 % Select from one of the built in kernels, using the name and
159 % geometry arguments supplied. See AcquireKernelBuiltIn()
161 % "WxH[+X+Y]:num, num, num ..."
162 % a kernal of size W by H, with W*H floating point numbers following.
163 % the 'center' can be optionally be defined at +X+Y (such that +0+0
164 % is top left corner). If not defined the pixel in the center, for
165 % odd sizes, or to the immediate top or left of center for even sizes
166 % is automatically selected.
168 % "num, num, num, num, ..."
169 % list of floating point numbers defining an 'old style' odd sized
170 % square kernel. At least 9 values should be provided for a 3x3
171 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
172 % Values can be space or comma separated. This is not recommended.
174 % Note that 'name' kernels will start with an alphabetic character while the
175 % new kernel specification has a ':' character in its specification string.
176 % If neither is the case, it is assumed an old style of a simple list of
177 % numbers generating a odd-sized square kernel has been given.
179 % The format of the AcquireKernal method is:
181 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
183 % A description of each parameter follows:
185 % o kernel_string: the Morphology/Convolution kernel wanted.
189 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
195 token[MaxTextExtent];
197 register unsigned long
210 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
212 assert(kernel_string != (const char *) NULL);
213 SetGeometryInfo(&args);
215 /* does it start with an alpha - Return a builtin kernel */
216 GetMagickToken(kernel_string,&p,token);
217 if ( isalpha((int)token[0]) )
222 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
223 if ( type < 0 || type == UserDefinedKernel )
224 return((KernelInfo *)NULL);
226 while (((isspace((int) ((unsigned char) *p)) != 0) ||
227 (*p == ',') || (*p == ':' )) && (*p != '\0'))
229 flags = ParseGeometry(p, &args);
231 /* special handling of missing values in input string */
232 if ( type == 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);
245 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
248 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
249 if (kernel == (KernelInfo *)NULL)
251 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
252 kernel->type = UserDefinedKernel;
253 kernel->signature = MagickSignature;
255 /* Has a ':' in argument - New user kernel specification */
256 p = strchr(kernel_string, ':');
257 if ( p != (char *) NULL)
259 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
260 memcpy(token, kernel_string, p-kernel_string);
261 token[p-kernel_string] = '\0';
262 flags = ParseGeometry(token, &args);
264 /* Size handling and checks of geometry settings */
265 if ( (flags & WidthValue) == 0 ) /* if no width then */
266 args.rho = args.sigma; /* then width = height */
267 if ( args.rho < 1.0 ) /* if width too small */
268 args.rho = 1.0; /* then width = 1 */
269 if ( args.sigma < 1.0 ) /* if height too small */
270 args.sigma = args.rho; /* then height = width */
271 kernel->width = (unsigned long)args.rho;
272 kernel->height = (unsigned long)args.sigma;
274 /* Offset Handling and Checks */
275 if ( args.xi < 0.0 || args.psi < 0.0 )
276 return(DestroyKernelInfo(kernel));
277 kernel->offset_x = ((flags & XValue)!=0) ? (unsigned long)args.xi
278 : (kernel->width-1)/2;
279 kernel->offset_y = ((flags & YValue)!=0) ? (unsigned long)args.psi
280 : (kernel->height-1)/2;
281 if ( kernel->offset_x >= kernel->width ||
282 kernel->offset_y >= kernel->height )
283 return(DestroyKernelInfo(kernel));
285 p++; /* advance beyond the ':' */
288 { /* ELSE - Old old kernel specification, forming odd-square kernel */
289 /* count up number of values given */
290 p=(const char *) kernel_string;
291 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
292 p++; /* ignore "'" chars for convolve filter usage - Cristy */
293 for (i=0; *p != '\0'; i++)
295 GetMagickToken(p,&p,token);
297 GetMagickToken(p,&p,token);
299 /* set the size of the kernel - old sized square */
300 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
301 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
302 p=(const char *) kernel_string;
303 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
304 p++; /* ignore "'" chars for convolve filter usage - Cristy */
307 /* Read in the kernel values from rest of input string argument */
308 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
309 kernel->height*sizeof(double));
310 if (kernel->values == (double *) NULL)
311 return(DestroyKernelInfo(kernel));
313 kernel->value_min = +MagickHuge;
314 kernel->value_max = -MagickHuge;
315 kernel->range_neg = kernel->range_pos = 0.0;
316 for (i=0; (i < kernel->width*kernel->height) && (*p != '\0'); i++)
318 GetMagickToken(p,&p,token);
320 GetMagickToken(p,&p,token);
321 if ( LocaleCompare("nan",token) == 0
322 || LocaleCompare("-",token) == 0 ) {
323 kernel->values[i] = nan; /* do not include this value in kernel */
326 kernel->values[i] = StringToDouble(token);
327 ( kernel->values[i] < 0)
328 ? ( kernel->range_neg += kernel->values[i] )
329 : ( kernel->range_pos += kernel->values[i] );
330 Minimize(kernel->value_min, kernel->values[i]);
331 Maximize(kernel->value_max, kernel->values[i]);
335 /* This should not be needed for a fully defined defined kernel
336 * Perhaps an error should be reported instead!
338 if ( i < kernel->width*kernel->height ) {
339 Minimize(kernel->value_min, kernel->values[i]);
340 Maximize(kernel->value_max, kernel->values[i]);
341 for ( ; i < kernel->width*kernel->height; i++)
342 kernel->values[i]=0.0;
349 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
353 % A c q u i r e K e r n e l B u i l t I n %
357 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
359 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
360 % kernels used for special purposes such as gaussian blurring, skeleton
361 % pruning, and edge distance determination.
363 % They take a KernelType, and a set of geometry style arguments, which were
364 % typically decoded from a user supplied string, or from a more complex
365 % Morphology Method that was requested.
367 % The format of the AcquireKernalBuiltIn method is:
369 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
370 % const GeometryInfo args)
372 % A description of each parameter follows:
374 % o type: the pre-defined type of kernel wanted
376 % o args: arguments defining or modifying the kernel
378 % Convolution Kernels
380 % Gaussian "[{radius}]x{sigma}"
381 % Generate a two-dimentional gaussian kernel, as used by -gaussian
382 % A sigma is required, (with the 'x'), due to historical reasons.
384 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
385 % the final size of the resulting kernel to a square 2*radius+1 in size.
386 % The radius should be at least 2 times that of the sigma value, or
387 % sever clipping and aliasing may result. If not given or set to 0 the
388 % radius will be determined so as to produce the best minimal error
389 % result, which is usally much larger than is normally needed.
391 % Blur "[{radius}]x{sigma}[+angle]"
392 % As per Gaussian, but generates a 1 dimensional or linear gaussian
393 % blur, at the angle given (current restricted to orthogonal angles).
394 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
396 % NOTE that two such blurs perpendicular to each other is equivelent to
397 % -blur and the previous gaussian, but is often 10 or more times faster.
399 % Comet "[{width}]x{sigma}[+angle]"
400 % Blur in one direction only, mush like how a bright object leaves
401 % a comet like trail. The Kernel is actually half a gaussian curve,
402 % Adding two such blurs in oppiste directions produces a Linear Blur.
404 % NOTE: that the first argument is the width of the kernel and not the
405 % radius of the kernel.
407 % # Still to be implemented...
409 % # Laplacian "{radius}x{sigma}"
410 % # Laplacian (a mexican hat like) Function
412 % # LOG "{radius},{sigma1},{sigma2}
413 % # Laplacian of Gaussian
415 % # DOG "{radius},{sigma1},{sigma2}
416 % # Difference of Gaussians
420 % Rectangle "{geometry}"
421 % Simply generate a rectangle of 1's with the size given. You can also
422 % specify the location of the 'control point', otherwise the closest
423 % pixel to the center of the rectangle is selected.
425 % Properly centered and odd sized rectangles work the best.
427 % Diamond "[{radius}]"
428 % Generate a diamond shaped kernal with given radius to the points.
429 % Kernel size will again be radius*2+1 square and defaults to radius 1,
430 % generating a 3x3 kernel that is slightly larger than a square.
432 % Square "[{radius}]"
433 % Generate a square shaped kernel of size radius*2+1, and defaulting
434 % to a 3x3 (radius 1).
436 % Note that using a larger radius for the "Square" or the "Diamond"
437 % is also equivelent to iterating the basic morphological method
438 % that many times. However However iterating with the smaller radius 1
439 % default is actually faster than using a larger kernel radius.
442 % Generate a binary disk of the radius given, radius may be a float.
443 % Kernel size will be ceil(radius)*2+1 square.
444 % NOTE: Here are some disk shapes of specific interest
445 % "disk:1" => "diamond" or "cross:1"
446 % "disk:1.5" => "square"
447 % "disk:2" => "diamond:2"
448 % "disk:2.5" => a general disk shape of radius 2
449 % "disk:2.9" => "square:2"
450 % "disk:3.5" => default - octagonal/disk shape of radius 3
451 % "disk:4.2" => roughly octagonal shape of radius 4
452 % "disk:4.3" => a general disk shape of radius 4
453 % After this all the kernel shape becomes more and more circular.
455 % Because a "disk" is more circular when using a larger radius, using a
456 % larger radius is preferred over iterating the morphological operation.
459 % Generate a kernel in the shape of a 'plus' sign. The length of each
460 % arm is also the radius, which defaults to 2.
462 % This kernel is not a good general morphological kernel, but is used
463 % more for highlighting and marking any single pixels in an image using,
464 % a "Dilate" or "Erode" method as appropriate.
466 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
468 % Note that unlike other kernels iterating a plus does not produce the
469 % same result as using a larger radius for the cross.
471 % Distance Measuring Kernels
473 % Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1)
474 % Manhatten "[{radius}][x{scale}]" square grid distance (default r=1)
475 % Euclidean "[{radius}][x{scale}]" direct distance (default r=1)
477 % Different types of distance measuring methods, which are used with the
478 % a 'Distance' morphology method for generating a gradient based on
479 % distance from an edge of a binary shape, though there is a technique
480 % for handling a anti-aliased shape.
482 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
483 % one to any neighbour, orthogonal or diagonal. One why of thinking of
484 % it is the number of squares a 'King' or 'Queen' in chess needs to
485 % traverse reach any other position on a chess board. It results in a
486 % 'square' like distance function, but one where diagonals are closer
489 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
490 % Cab metric), is the distance needed when you can only travel in
491 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
492 % in chess would travel. It results in a diamond like distances, where
493 % diagonals are further than expected.
495 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
496 % However by default the kernel size only has a radius of 1, which
497 % limits the distance to 'Knight' like moves, with only orthogonal and
498 % diagonal measurements being correct. As such for the default kernel
499 % you will get octagonal like distance function, which is reasonally
502 % However if you use a larger radius such as "Euclidean:4" you will
503 % get a much smoother distance gradient from the edge of the shape.
504 % Of course a larger kernel is slower to use, and generally not needed.
506 % To allow the use of fractional distances that you get with diagonals
507 % the actual distance is scaled by a fixed value which the user can
508 % provide. This is not actually nessary for either ""Chebyshev" or
509 % "Manhatten" distance kernels, but is done for all three distance
510 % kernels. If no scale is provided it is set to a value of 100,
511 % allowing for a maximum distance measurement of 655 pixels using a Q16
512 % version of IM, from any edge. However for small images this can
513 % result in quite a dark gradient.
515 % See the 'Distance' Morphological Method, for information of how it is
520 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
521 const GeometryInfo *args)
526 register unsigned long
534 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
536 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
537 if (kernel == (KernelInfo *) NULL)
539 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
540 kernel->value_min = kernel->value_max = 0.0;
541 kernel->range_neg = kernel->range_pos = 0.0;
543 kernel->signature = MagickSignature;
546 /* Convolution Kernels */
549 sigma = fabs(args->sigma);
551 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
553 kernel->width = kernel->height =
554 GetOptimalKernelWidth2D(args->rho,sigma);
555 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
556 kernel->range_neg = kernel->range_pos = 0.0;
557 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
558 kernel->height*sizeof(double));
559 if (kernel->values == (double *) NULL)
560 return(DestroyKernelInfo(kernel));
562 sigma = 2.0*sigma*sigma; /* simplify the expression */
563 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
564 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
565 kernel->range_pos += (
567 exp(-((double)(u*u+v*v))/sigma)
568 /* / (MagickPI*sigma) */ );
569 kernel->value_min = 0;
570 kernel->value_max = kernel->values[
571 kernel->offset_y*kernel->width+kernel->offset_x ];
573 NormalizeKernel(kernel);
579 sigma = fabs(args->sigma);
581 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
583 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
584 kernel->offset_x = (kernel->width-1)/2;
586 kernel->offset_y = 0;
587 kernel->range_neg = kernel->range_pos = 0.0;
588 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
589 kernel->height*sizeof(double));
590 if (kernel->values == (double *) NULL)
591 return(DestroyKernelInfo(kernel));
595 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
596 ** It generates a gaussian 3 times the width, and compresses it into
597 ** the expected range. This produces a closer normalization of the
598 ** resulting kernel, especially for very low sigma values.
599 ** As such while wierd it is prefered.
601 ** I am told this method originally came from Photoshop.
603 sigma *= KernelRank; /* simplify expanded curve */
604 v = (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
605 (void) ResetMagickMemory(kernel->values,0, (size_t)
606 kernel->width*sizeof(double));
607 for ( u=-v; u <= v; u++) {
608 kernel->values[(u+v)/KernelRank] +=
609 exp(-((double)(u*u))/(2.0*sigma*sigma))
610 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
612 for (i=0; i < kernel->width; i++)
613 kernel->range_pos += kernel->values[i];
615 for ( i=0, u=-kernel->offset_x; i < kernel->width; i++, u++)
616 kernel->range_pos += (
618 exp(-((double)(u*u))/(2.0*sigma*sigma))
619 /* / (MagickSQ2PI*sigma) */ );
621 kernel->value_min = 0;
622 kernel->value_max = kernel->values[ kernel->offset_x ];
623 /* Note that both the above methods do not generate a normalized
624 ** kernel, though it gets close. The kernel may be 'clipped' by a user
625 ** defined radius, producing a smaller (darker) kernel. Also for very
626 ** small sigma's (> 0.1) the central value becomes larger than one,
627 ** and thus producing a very bright kernel.
630 /* Normalize the 1D Gaussian Kernel
632 ** Because of this the divisor in the above kernel generator is
633 ** not needed, so is not done above.
635 NormalizeKernel(kernel);
637 /* rotate the kernel by given angle */
638 RotateKernel(kernel, args->xi);
643 sigma = fabs(args->sigma);
645 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
647 if ( args->rho < 1.0 )
648 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
650 kernel->width = (unsigned long)args->rho;
651 kernel->offset_x = kernel->offset_y = 0;
653 kernel->range_neg = kernel->range_pos = 0.0;
654 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
655 kernel->height*sizeof(double));
656 if (kernel->values == (double *) NULL)
657 return(DestroyKernelInfo(kernel));
659 /* A comet blur is half a gaussian curve, so that the object is
660 ** blurred in one direction only. This may not be quite the right
661 ** curve so may change in the future. The function must be normalised.
665 sigma *= KernelRank; /* simplify expanded curve */
666 v = kernel->width*KernelRank; /* start/end points to fit range */
667 (void) ResetMagickMemory(kernel->values,0, (size_t)
668 kernel->width*sizeof(double));
669 for ( u=0; u < v; u++) {
670 kernel->values[u/KernelRank] +=
671 exp(-((double)(u*u))/(2.0*sigma*sigma))
672 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
674 for (i=0; i < kernel->width; i++)
675 kernel->range_pos += kernel->values[i];
677 for ( i=0; i < kernel->width; i++)
678 kernel->range_pos += (
680 exp(-((double)(i*i))/(2.0*sigma*sigma))
681 /* / (MagickSQ2PI*sigma) */ );
683 kernel->value_min = 0;
684 kernel->value_max = kernel->values[0];
686 NormalizeKernel(kernel);
687 RotateKernel(kernel, args->xi);
690 /* Boolean Kernels */
691 case RectangleKernel:
694 if ( type == SquareKernel )
697 kernel->width = kernel->height = 3; /* default radius = 1 */
699 kernel->width = kernel->height = 2*(long)args->rho+1;
700 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
703 /* NOTE: user defaults set in "AcquireKernelInfo()" */
704 if ( args->rho < 1.0 || args->sigma < 1.0 )
705 return(DestroyKernelInfo(kernel)); /* invalid args given */
706 kernel->width = (unsigned long)args->rho;
707 kernel->height = (unsigned long)args->sigma;
708 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
709 args->psi < 0.0 || args->psi > (double)kernel->height )
710 return(DestroyKernelInfo(kernel)); /* invalid args given */
711 kernel->offset_x = (unsigned long)args->xi;
712 kernel->offset_y = (unsigned long)args->psi;
714 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
715 kernel->height*sizeof(double));
716 if (kernel->values == (double *) NULL)
717 return(DestroyKernelInfo(kernel));
719 u=kernel->width*kernel->height;
720 for ( i=0; i < (unsigned long)u; i++)
721 kernel->values[i] = 1.0;
723 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
724 kernel->range_pos = (double) u;
729 kernel->width = kernel->height = 3; /* default radius = 1 */
731 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
732 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
734 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
735 kernel->height*sizeof(double));
736 if (kernel->values == (double *) NULL)
737 return(DestroyKernelInfo(kernel));
739 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
740 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
741 if ((labs(u)+labs(v)) <= (long)kernel->offset_x)
742 kernel->range_pos += kernel->values[i] = 1.0;
744 kernel->values[i] = nan;
745 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
753 limit = (long)(args->rho*args->rho);
754 if (args->rho < 0.1) /* default radius approx 3.5 */
755 kernel->width = kernel->height = 7L, limit = 10L;
757 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
758 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
760 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
761 kernel->height*sizeof(double));
762 if (kernel->values == (double *) NULL)
763 return(DestroyKernelInfo(kernel));
765 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
766 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
767 if ((u*u+v*v) <= limit)
768 kernel->range_pos += kernel->values[i] = 1.0;
770 kernel->values[i] = nan;
771 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
777 kernel->width = kernel->height = 5; /* default radius 2 */
779 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
780 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
782 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
783 kernel->height*sizeof(double));
784 if (kernel->values == (double *) NULL)
785 return(DestroyKernelInfo(kernel));
787 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
788 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
789 kernel->values[i] = (u == 0 || v == 0) ? 1.0 : nan;
790 kernel->value_min = kernel->value_max = 1.0; /* a flat kernel */
791 kernel->range_pos = kernel->width*2.0 - 1.0;
794 /* Distance Measuring Kernels */
795 case ChebyshevKernel:
801 kernel->width = kernel->height = 3; /* default radius = 1 */
803 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
804 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
806 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
807 kernel->height*sizeof(double));
808 if (kernel->values == (double *) NULL)
809 return(DestroyKernelInfo(kernel));
811 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
812 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
813 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
814 kernel->range_pos += ( kernel->values[i] =
815 scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
816 kernel->value_max = kernel->values[0];
819 case ManhattenKernel:
825 kernel->width = kernel->height = 3; /* default radius = 1 */
827 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
828 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
830 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
831 kernel->height*sizeof(double));
832 if (kernel->values == (double *) NULL)
833 return(DestroyKernelInfo(kernel));
835 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
836 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
837 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
838 kernel->range_pos += ( kernel->values[i] =
839 scale*(labs(u)+labs(v)) );
840 kernel->value_max = kernel->values[0];
843 case EuclideanKernel:
849 kernel->width = kernel->height = 3; /* default radius = 1 */
851 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
852 kernel->offset_x = kernel->offset_y = (kernel->width-1)/2;
854 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
855 kernel->height*sizeof(double));
856 if (kernel->values == (double *) NULL)
857 return(DestroyKernelInfo(kernel));
859 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
860 for ( i=0, v=-kernel->offset_y; v <= (long)kernel->offset_y; v++)
861 for ( u=-kernel->offset_x; u <= (long)kernel->offset_x; u++, i++)
862 kernel->range_pos += ( kernel->values[i] =
863 scale*sqrt((double)(u*u+v*v)) );
864 kernel->value_max = kernel->values[0];
867 /* Undefined Kernels */
868 case LaplacianKernel:
871 assert("Kernel Type has not been defined yet");
874 /* Generate a No-Op minimal kernel - 1x1 pixel */
875 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
876 if (kernel->values == (double *) NULL)
877 return(DestroyKernelInfo(kernel));
878 kernel->width = kernel->height = 1;
879 kernel->offset_x = kernel->offset_x = 0;
880 kernel->type = UndefinedKernel;
883 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
891 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
895 % D e s t r o y K e r n e l I n f o %
899 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
901 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
904 % The format of the DestroyKernelInfo method is:
906 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
908 % A description of each parameter follows:
910 % o kernel: the Morphology/Convolution kernel to be destroyed
914 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
916 assert(kernel != (KernelInfo *) NULL);
917 kernel->values=(double *)RelinquishMagickMemory(kernel->values);
918 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
923 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
927 % M o r p h o l o g y I m a g e C h a n n e l %
931 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
933 % MorphologyImageChannel() applies a user supplied kernel to the image
934 % according to the given mophology method.
936 % The given kernel is assumed to have been pre-scaled appropriatally, usally
937 % by the kernel generator.
939 % The format of the MorphologyImage method is:
941 % Image *MorphologyImage(const Image *image, MorphologyMethod method,
942 % const long iterations, KernelInfo *kernel, ExceptionInfo *exception)
943 % Image *MorphologyImageChannel(const Image *image, const ChannelType
944 % channel, MorphologyMethod method, const long iterations, KernelInfo
945 % *kernel, ExceptionInfo *exception)
947 % A description of each parameter follows:
949 % o image: the image.
951 % o method: the morphology method to be applied.
953 % o iterations: apply the operation this many times (or no change).
954 % A value of -1 means loop until no change found.
955 % How this is applied may depend on the morphology method.
956 % Typically this is a value of 1.
958 % o channel: the channel type.
960 % o kernel: An array of double representing the morphology kernel.
961 % Warning: kernel may be normalized for the Convolve method.
963 % o exception: return any errors or warnings in this structure.
966 % TODO: bias and auto-scale handling of the kernel for convolution
967 % The given kernel is assumed to have been pre-scaled appropriatally, usally
968 % by the kernel generator.
973 * Apply the Morphology method with the given Kernel
974 * And return the number of pixels that changed.
976 static unsigned long MorphologyApply(const Image *image, Image
977 *result_image, const MorphologyMethod method, const ChannelType channel,
978 const KernelInfo *kernel, ExceptionInfo *exception)
980 #define MorphologyTag "Morphology/Image"
1000 Apply Morphology to Image.
1006 GetMagickPixelPacket(image,&bias);
1007 SetMagickPixelPacketBias(image,&bias);
1009 p_view=AcquireCacheView(image);
1010 q_view=AcquireCacheView(result_image);
1012 /* some methods (including convolve) needs to apply a reflected kernel
1013 * the offset for getting the kernel view needs to be adjusted for this
1016 offx = kernel->offset_x;
1017 offy = kernel->offset_y;
1019 case ErodeMorphology:
1020 case ErodeIntensityMorphology:
1021 /* kernel is not reflected */
1024 /* kernel needs to be reflected */
1025 offx = kernel->width-offx-1;
1026 offy = kernel->height-offy-1;
1030 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1031 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1033 for (y=0; y < image->rows; y++)
1038 register const PixelPacket
1041 register const IndexPacket
1042 *restrict p_indexes;
1044 register PixelPacket
1047 register IndexPacket
1048 *restrict q_indexes;
1050 register unsigned long
1056 if (status == MagickFalse)
1058 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1059 image->columns+kernel->width, kernel->height, exception);
1060 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1062 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1067 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1068 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1069 r = (image->columns+kernel->width)*offy+offx; /* constant */
1071 for (x=0; x < image->columns; x++)
1076 register unsigned long
1079 register const double
1082 register const PixelPacket
1085 register const IndexPacket
1086 *restrict k_indexes;
1091 /* Copy input to ouput image for unused channels
1092 * This removes need for 'cloning' a new image every iteration
1095 if (image->colorspace == CMYKColorspace)
1096 q_indexes[x] = p_indexes[r];
1098 result.index=0; /* stop compiler warnings */
1100 case ConvolveMorphology:
1102 break; /* default result is the convolution bias */
1104 case DilateMorphology:
1109 result.index = -MagickHuge;
1111 case ErodeMorphology:
1116 result.index = +MagickHuge;
1120 /* Otherwise just start with the original pixel value */
1121 result.red = p[r].red;
1122 result.green = p[r].green;
1123 result.blue = p[r].blue;
1124 result.opacity = QuantumRange - p[r].opacity;
1125 if ( image->colorspace == CMYKColorspace)
1126 result.index = p_indexes[r];
1131 case ConvolveMorphology:
1132 /* Weighted Average of pixels
1134 * NOTE for correct working of this operation for asymetrical
1135 * kernels, the kernel needs to be applied in its reflected form.
1136 * That is its values needs to be reversed.
1138 if (((channel & OpacityChannel) == 0) ||
1139 (image->matte == MagickFalse))
1141 /* Convolution (no transparency) */
1142 k = &kernel->values[ kernel->width*kernel->height-1 ];
1144 k_indexes = p_indexes;
1145 for (v=0; v < kernel->height; v++) {
1146 for (u=0; u < kernel->width; u++, k--) {
1147 if ( IsNan(*k) ) continue;
1148 result.red += (*k)*k_pixels[u].red;
1149 result.green += (*k)*k_pixels[u].green;
1150 result.blue += (*k)*k_pixels[u].blue;
1151 /* result.opacity += not involved here */
1152 if ( image->colorspace == CMYKColorspace)
1153 result.index += (*k)*k_indexes[u];
1155 k_pixels += image->columns+kernel->width;
1156 k_indexes += image->columns+kernel->width;
1160 { /* Kernel & Alpha weighted Convolution */
1162 alpha, /* alpha value * kernel weighting */
1163 gamma; /* weighting divisor */
1166 k = &kernel->values[ kernel->width*kernel->height-1 ];
1168 k_indexes = p_indexes;
1169 for (v=0; v < kernel->height; v++) {
1170 for (u=0; u < kernel->width; u++, k--) {
1171 if ( IsNan(*k) ) continue;
1172 alpha=(*k)*(QuantumScale*(QuantumRange-
1173 k_pixels[u].opacity));
1175 result.red += alpha*k_pixels[u].red;
1176 result.green += alpha*k_pixels[u].green;
1177 result.blue += alpha*k_pixels[u].blue;
1178 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1179 if ( image->colorspace == CMYKColorspace)
1180 result.index += alpha*k_indexes[u];
1182 k_pixels += image->columns+kernel->width;
1183 k_indexes += image->columns+kernel->width;
1185 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1186 result.red *= gamma;
1187 result.green *= gamma;
1188 result.blue *= gamma;
1189 result.opacity *= gamma;
1190 result.index *= gamma;
1194 case DilateMorphology:
1195 /* Maximize Value within kernel shape
1197 * NOTE for correct working of this operation for asymetrical
1198 * kernels, the kernel needs to be applied in its reflected form.
1199 * That is its values needs to be reversed.
1201 * NOTE: in normal Greyscale Morphology, the kernel value should
1202 * be added to the real value, this is currently not done, due to
1203 * the nature of the boolean kernels being used.
1206 k = &kernel->values[ kernel->width*kernel->height-1 ];
1208 k_indexes = p_indexes;
1209 for (v=0; v < kernel->height; v++) {
1210 for (u=0; u < kernel->width; u++, k--) {
1211 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1212 Maximize(result.red, k_pixels[u].red);
1213 Maximize(result.green, k_pixels[u].green);
1214 Maximize(result.blue, k_pixels[u].blue);
1215 Maximize(result.opacity, QuantumRange-k_pixels[u].opacity);
1216 if ( image->colorspace == CMYKColorspace)
1217 Maximize(result.index, k_indexes[u]);
1219 k_pixels += image->columns+kernel->width;
1220 k_indexes += image->columns+kernel->width;
1224 case ErodeMorphology:
1225 /* Minimize Value within kernel shape
1227 * NOTE that the kernel is not reflected for this operation!
1229 * NOTE: in normal Greyscale Morphology, the kernel value should
1230 * be added to the real value, this is currently not done, due to
1231 * the nature of the boolean kernels being used.
1235 k_indexes = p_indexes;
1236 for (v=0; v < kernel->height; v++) {
1237 for (u=0; u < kernel->width; u++, k++) {
1238 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1239 Minimize(result.red, k_pixels[u].red);
1240 Minimize(result.green, k_pixels[u].green);
1241 Minimize(result.blue, k_pixels[u].blue);
1242 Minimize(result.opacity, QuantumRange-k_pixels[u].opacity);
1243 if ( image->colorspace == CMYKColorspace)
1244 Minimize(result.index, k_indexes[u]);
1246 k_pixels += image->columns+kernel->width;
1247 k_indexes += image->columns+kernel->width;
1251 case DilateIntensityMorphology:
1252 /* Select pixel with maximum intensity within kernel shape
1254 * WARNING: the intensity test fails for CMYK and does not
1255 * take into account the moderating effect of teh alpha channel
1258 * NOTE for correct working of this operation for asymetrical
1259 * kernels, the kernel needs to be applied in its reflected form.
1260 * That is its values needs to be reversed.
1262 k = &kernel->values[ kernel->width*kernel->height-1 ];
1264 k_indexes = p_indexes;
1265 for (v=0; v < kernel->height; v++) {
1266 for (u=0; u < kernel->width; u++, k--) {
1267 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1268 if ( result.red == 0.0 ||
1269 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1270 /* copy the whole pixel - no channel selection */
1272 if ( result.red > 0.0 ) changed++;
1276 k_pixels += image->columns+kernel->width;
1277 k_indexes += image->columns+kernel->width;
1281 case ErodeIntensityMorphology:
1282 /* Select pixel with mimimum intensity within kernel shape
1284 * WARNING: the intensity test fails for CMYK and does not
1285 * take into account the moderating effect of teh alpha channel
1288 * NOTE that the kernel is not reflected for this operation!
1292 k_indexes = p_indexes;
1293 for (v=0; v < kernel->height; v++) {
1294 for (u=0; u < kernel->width; u++, k++) {
1295 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1296 if ( result.red == 0.0 ||
1297 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1298 /* copy the whole pixel - no channel selection */
1300 if ( result.red > 0.0 ) changed++;
1304 k_pixels += image->columns+kernel->width;
1305 k_indexes += image->columns+kernel->width;
1309 case DistanceMorphology:
1310 /* Add kernel value and select the minimum value found.
1311 * The result is a iterative distance from edge function.
1313 * All Distance Kernels are symetrical, but that may not always
1314 * be the case. For example how about a distance from left edges?
1315 * To make it work correctly for asymetrical kernels the reflected
1316 * kernel needs to be applied.
1319 /* No need to do distance morphology if original value is zero
1320 * Unfortunatally I have not been able to get this right
1321 * when channel selection also becomes involved. -- Arrgghhh
1323 if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1324 || ((channel & GreenChannel) == 0 && p[r].green == 0)
1325 || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1326 || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1327 || (( (channel & IndexChannel) == 0
1328 || image->colorspace != CMYKColorspace
1329 ) && p_indexes[x] ==0 )
1333 k = &kernel->values[ kernel->width*kernel->height-1 ];
1335 k_indexes = p_indexes;
1336 for (v=0; v < kernel->height; v++) {
1337 for (u=0; u < kernel->width; u++, k--) {
1338 if ( IsNan(*k) ) continue;
1339 Minimize(result.red, (*k)+k_pixels[u].red);
1340 Minimize(result.green, (*k)+k_pixels[u].green);
1341 Minimize(result.blue, (*k)+k_pixels[u].blue);
1342 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1343 if ( image->colorspace == CMYKColorspace)
1344 Minimize(result.index, (*k)+k_indexes[u]);
1346 k_pixels += image->columns+kernel->width;
1347 k_indexes += image->columns+kernel->width;
1351 case UndefinedMorphology:
1353 break; /* Do nothing */
1356 case UndefinedMorphology:
1357 case DilateIntensityMorphology:
1358 case ErodeIntensityMorphology:
1359 break; /* full pixel was directly assigned */
1361 /* Assign the results */
1362 if ((channel & RedChannel) != 0)
1363 q->red = ClampToQuantum(result.red);
1364 if ((channel & GreenChannel) != 0)
1365 q->green = ClampToQuantum(result.green);
1366 if ((channel & BlueChannel) != 0)
1367 q->blue = ClampToQuantum(result.blue);
1368 if ((channel & OpacityChannel) != 0
1369 && image->matte == MagickTrue )
1370 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1371 if ((channel & IndexChannel) != 0
1372 && image->colorspace == CMYKColorspace)
1373 q_indexes[x] = ClampToQuantum(result.index);
1376 if ( ( p[r].red != q->red )
1377 || ( p[r].green != q->green )
1378 || ( p[r].blue != q->blue )
1379 || ( p[r].opacity != q->opacity )
1380 || ( image->colorspace == CMYKColorspace &&
1381 p_indexes[r] != q_indexes[x] ) )
1382 changed++; /* The pixel had some value changed! */
1386 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1387 if (sync == MagickFalse)
1389 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1394 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1395 #pragma omp critical (MagickCore_MorphologyImage)
1397 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1398 if (proceed == MagickFalse)
1402 result_image->type=image->type;
1403 q_view=DestroyCacheView(q_view);
1404 p_view=DestroyCacheView(p_view);
1405 return(status ? changed : 0);
1408 MagickExport Image *MorphologyImage(const Image *image,MorphologyMethod method,
1409 const long iterations,KernelInfo *kernel, ExceptionInfo *exception)
1414 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1415 iterations,kernel,exception);
1416 return(morphology_image);
1419 MagickExport Image *MorphologyImageChannel(const Image *image,
1420 const ChannelType channel, MorphologyMethod method, const long iterations,
1421 KernelInfo *kernel, ExceptionInfo *exception)
1432 assert(image != (Image *) NULL);
1433 assert(image->signature == MagickSignature);
1434 assert(exception != (ExceptionInfo *) NULL);
1435 assert(exception->signature == MagickSignature);
1437 if ( GetImageArtifact(image,"showkernel") != (const char *) NULL)
1438 ShowKernel(kernel); /* request to display the kernel to stderr */
1440 if ( iterations == 0 )
1441 return((Image *)NULL); /* null operation - nothing to do! */
1443 /* kernel must be valid at this point
1444 * (except maybe for posible future morphology methods like "Prune"
1446 assert(kernel != (KernelInfo *)NULL);
1451 if ( iterations < 0 )
1452 limit = image->columns > image->rows ? image->columns : image->rows;
1454 /* Special morphology cases */
1456 case CloseMorphology:
1457 new_image = MorphologyImageChannel(image, channel, DilateMorphology,
1458 iterations, kernel, exception);
1459 if (new_image == (Image *) NULL)
1460 return((Image *) NULL);
1461 method = ErodeMorphology;
1463 case OpenMorphology:
1464 new_image = MorphologyImageChannel(image, channel, ErodeMorphology,
1465 iterations, kernel, exception);
1466 if (new_image == (Image *) NULL)
1467 return((Image *) NULL);
1468 method = DilateMorphology;
1470 case CloseIntensityMorphology:
1471 new_image = MorphologyImageChannel(image, channel, DilateIntensityMorphology,
1472 iterations, kernel, exception);
1473 if (new_image == (Image *) NULL)
1474 return((Image *) NULL);
1475 method = ErodeIntensityMorphology;
1477 case OpenIntensityMorphology:
1478 new_image = MorphologyImageChannel(image, channel, ErodeIntensityMorphology,
1479 iterations, kernel, exception);
1480 if (new_image == (Image *) NULL)
1481 return((Image *) NULL);
1482 method = DilateIntensityMorphology;
1485 case ConvolveMorphology:
1486 /* NormalizeKernel(kernel); removed, by default use kernel as defined */
1487 /* TODO: auto-bias, auto-scaling and user-scaling of kernel according
1488 * to expert user settings */
1491 /* Do a morphology just once at this point!
1492 This ensures a new_image has been generated, but allows us
1493 to skip the creation of 'old_image' if it isn't needed.
1495 new_image=CloneImage(image,0,0,MagickTrue,exception);
1496 if (new_image == (Image *) NULL)
1497 return((Image *) NULL);
1498 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1500 InheritException(exception,&new_image->exception);
1501 new_image=DestroyImage(new_image);
1502 return((Image *) NULL);
1504 changed = MorphologyApply(image,new_image,method,channel,kernel,
1507 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1508 fprintf(stderr, "Morphology %s:%lu => Changed %lu\n",
1509 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1513 /* Repeat the interative morphology until count or no change */
1514 if ( count < limit && changed > 0 ) {
1515 old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1516 if (old_image == (Image *) NULL)
1517 return(DestroyImage(new_image));
1518 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1520 InheritException(exception,&old_image->exception);
1521 old_image=DestroyImage(old_image);
1522 return(DestroyImage(new_image));
1524 while( count < limit && changed != 0 )
1526 Image *tmp = old_image;
1527 old_image = new_image;
1529 changed = MorphologyApply(old_image,new_image,method,channel,kernel,
1532 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1533 fprintf(stderr, "Morphology %s:%lu => Changed %lu\n",
1534 MagickOptionToMnemonic(MagickMorphologyOptions, method),
1537 DestroyImage(old_image);
1544 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1548 + N o r m a l i z e K e r n e l %
1552 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1554 % NormalizeKernel() normalize the kernel so its convolution output will
1555 % be over a unit range.
1557 % Assumes the 'range_*' attributes of the kernel structure have been
1558 % correctly set during the kernel creation.
1560 % The format of the NormalizeKernel method is:
1562 % void NormalizeKernel(KernelInfo *kernel)
1564 % A description of each parameter follows:
1566 % o kernel: the Morphology/Convolution kernel
1569 static void NormalizeKernel(KernelInfo *kernel)
1571 register unsigned long
1574 for (i=0; i < kernel->width*kernel->height; i++)
1575 kernel->values[i] /= (kernel->range_pos - kernel->range_neg);
1577 kernel->range_pos /= (kernel->range_pos - kernel->range_neg);
1578 kernel->range_neg /= (kernel->range_pos - kernel->range_neg);
1579 kernel->value_max /= (kernel->range_pos - kernel->range_neg);
1580 kernel->value_min /= (kernel->range_pos - kernel->range_neg);
1586 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1590 % R o t a t e K e r n e l %
1594 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1596 % RotateKernel() rotates the kernel by the angle given. Currently it is
1597 % restricted to 90 degree angles, but this may be improved in the future.
1599 % The format of the RotateKernel method is:
1601 % void RotateKernel(KernelInfo *kernel, double angle)
1603 % A description of each parameter follows:
1605 % o kernel: the Morphology/Convolution kernel
1607 % o angle: angle to rotate in degrees
1610 static void RotateKernel(KernelInfo *kernel, double angle)
1612 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
1614 ** TODO: expand beyond simple 90 degree rotates, flips and flops
1617 /* Modulus the angle */
1618 angle = fmod(angle, 360.0);
1622 if ( 315.0 < angle || angle <= 45.0 )
1623 return; /* no change! - At least at this time */
1625 switch (kernel->type) {
1626 /* These built-in kernels are cylindrical kernels, rotating is useless */
1627 case GaussianKernel:
1628 case LaplacianKernel:
1632 case ChebyshevKernel:
1633 case ManhattenKernel:
1634 case EuclideanKernel:
1637 /* These may be rotatable at non-90 angles in the future */
1638 /* but simply rotating them in multiples of 90 degrees is useless */
1644 /* These only allows a +/-90 degree rotation (by transpose) */
1645 /* A 180 degree rotation is useless */
1647 case RectangleKernel:
1648 if ( 135.0 < angle && angle <= 225.0 )
1650 if ( 225.0 < angle && angle <= 315.0 )
1654 /* these are freely rotatable in 90 degree units */
1656 case UndefinedKernel:
1657 case UserDefinedKernel:
1660 if ( 135.0 < angle && angle <= 225.0 )
1662 /* For a 180 degree rotation - also know as a reflection */
1663 /* This is actually a very very common operation! */
1664 /* Basically all that is needed is a reversal of the kernel data! */
1671 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
1672 t=k[i], k[i]=k[j], k[j]=t;
1674 kernel->offset_x = kernel->width - kernel->offset_x - 1;
1675 kernel->offset_y = kernel->width - kernel->offset_y - 1;
1676 angle = fmod(angle+180.0, 360.0);
1678 if ( 45.0 < angle && angle <= 135.0 )
1679 { /* Do a transpose and a flop, of the image, which results in a 90
1680 * degree rotation using two mirror operations.
1682 * WARNING: this assumes the original image was a 1 dimentional image
1683 * but currently that is the only built-ins it is applied to.
1688 kernel->width = kernel->height;
1690 t = kernel->offset_x;
1691 kernel->offset_x = kernel->offset_y;
1692 kernel->offset_y = t;
1693 angle = fmod(450.0 - angle, 360.0);
1695 /* At this point angle should be between -45 (315) and +45 degrees
1696 * In the future some form of non-orthogonal angled rotates could be
1697 * performed here, posibily with a linear kernel restriction.
1701 Not currently in use!
1702 { /* Do a flop, this assumes kernel is horizontally symetrical.
1703 * Each row of the kernel needs to be reversed!
1707 register unsigned long
1712 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
1713 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
1714 t=k[x], k[x]=k[r], k[r]=t;
1716 kernel->offset_x = kernel->width - kernel->offset_x - 1;
1717 angle = fmod(angle+180.0, 360.0);
1724 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1728 % S h o w K e r n e l %
1732 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1734 % ShowKernel() Output the details of the given kernel defination to
1735 % standard error, as per a users 'showkernel' option request.
1737 % The format of the ShowKernel method is:
1739 % void KernelPrint (KernelInfo *kernel)
1741 % A description of each parameter follows:
1743 % o kernel: the Morphology/Convolution kernel
1745 % FUTURE: return the information in a string for API usage.
1747 static void ShowKernel(KernelInfo *kernel)
1753 "Kernel \"%s\" of size %lux%lu%+ld%+ld with values from %lg to %lg\n",
1754 MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
1755 kernel->width, kernel->height,
1756 kernel->offset_x, kernel->offset_y,
1757 kernel->value_min, kernel->value_max);
1758 fprintf(stderr, "Forming convolution output range from %lg to %lg%s\n",
1759 kernel->range_neg, kernel->range_pos,
1760 kernel->normalized == MagickTrue ? " (normalized)" : "" );
1761 for (i=v=0; v < kernel->height; v++) {
1762 fprintf(stderr,"%2ld:",v);
1763 for (u=0; u < kernel->width; u++, i++)
1764 if ( IsNan(kernel->values[i]) )
1765 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
1767 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
1768 GetMagickPrecision(), kernel->values[i]);
1769 fprintf(stderr,"\n");