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"
81 The following test is for special floating point numbers of value NaN (not
82 a number), that may be used within a Kernel Definition. NaN's are defined
83 as part of the IEEE standard for floating point number representation.
85 These are used a Kernel value of NaN means that that kernal position is not
86 part of the normal convolution or morphology process, and thus allowing the
87 use of 'shaped' kernels.
89 Special properities two NaN's are never equal, even if they are from the
90 same variable That is the IsNaN() macro is only true if the value is NaN.
92 #define IsNan(a) ((a)!=(a))
95 Other global definitions used by module.
97 static inline double MagickMin(const double x,const double y)
99 return( x < y ? x : y);
101 static inline double MagickMax(const double x,const double y)
103 return( x > y ? x : y);
105 #define Minimize(assign,value) assign=MagickMin(assign,value)
106 #define Maximize(assign,value) assign=MagickMax(assign,value)
108 /* Currently these are only internal to this module */
110 RotateKernelInfo(KernelInfo *, double);
113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117 % A c q u i r e K e r n e l I n f o %
121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123 % AcquireKernelInfo() takes the given string (generally supplied by the
124 % user) and converts it into a Morphology/Convolution Kernel. This allows
125 % users to specify a kernel from a number of pre-defined kernels, or to fully
126 % specify their own kernel for a specific Convolution or Morphology
129 % The kernel so generated can be any rectangular array of floating point
130 % values (doubles) with the 'control point' or 'pixel being affected'
131 % anywhere within that array of values.
133 % Previously IM was restricted to a square of odd size using the exact
134 % center as origin, this is no longer the case, and any rectangular kernel
135 % with any value being declared the origin. This in turn allows the use of
136 % highly asymmetrical kernels.
138 % The floating point values in the kernel can also include a special value
139 % known as 'nan' or 'not a number' to indicate that this value is not part
140 % of the kernel array. This allows you to shaped the kernel within its
141 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
142 % shape. However at least one non-nan value must be provided for correct
143 % working of a kernel.
145 % The returned kernel should be free using the DestroyKernelInfo() when you
146 % are finished with it.
148 % Input kernel defintion strings can consist of any of three types.
151 % Select from one of the built in kernels, using the name and
152 % geometry arguments supplied. See AcquireKernelBuiltIn()
154 % "WxH[+X+Y]:num, num, num ..."
155 % a kernal of size W by H, with W*H floating point numbers following.
156 % the 'center' can be optionally be defined at +X+Y (such that +0+0
157 % is top left corner). If not defined the pixel in the center, for
158 % odd sizes, or to the immediate top or left of center for even sizes
159 % is automatically selected.
161 % "num, num, num, num, ..."
162 % list of floating point numbers defining an 'old style' odd sized
163 % square kernel. At least 9 values should be provided for a 3x3
164 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
165 % Values can be space or comma separated. This is not recommended.
167 % Note that 'name' kernels will start with an alphabetic character while the
168 % new kernel specification has a ':' character in its specification string.
169 % If neither is the case, it is assumed an old style of a simple list of
170 % numbers generating a odd-sized square kernel has been given.
172 % The format of the AcquireKernal method is:
174 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
176 % A description of each parameter follows:
178 % o kernel_string: the Morphology/Convolution kernel wanted.
182 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
188 token[MaxTextExtent];
203 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
205 assert(kernel_string != (const char *) NULL);
206 SetGeometryInfo(&args);
208 /* does it start with an alpha - Return a builtin kernel */
209 GetMagickToken(kernel_string,&p,token);
210 if ( isalpha((int)token[0]) )
215 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
216 if ( type < 0 || type == UserDefinedKernel )
217 return((KernelInfo *)NULL);
219 while (((isspace((int) ((unsigned char) *p)) != 0) ||
220 (*p == ',') || (*p == ':' )) && (*p != '\0'))
222 flags = ParseGeometry(p, &args);
224 /* special handling of missing values in input string */
226 case RectangleKernel:
227 if ( (flags & WidthValue) == 0 ) /* if no width then */
228 args.rho = args.sigma; /* then width = height */
229 if ( args.rho < 1.0 ) /* if width too small */
230 args.rho = 3; /* then width = 3 */
231 if ( args.sigma < 1.0 ) /* if height too small */
232 args.sigma = args.rho; /* then height = width */
233 if ( (flags & XValue) == 0 ) /* center offset if not defined */
234 args.xi = (double)(((long)args.rho-1)/2);
235 if ( (flags & YValue) == 0 )
236 args.psi = (double)(((long)args.sigma-1)/2);
242 if ( (flags & HeightValue) == 0 ) /* if no scale */
243 args.sigma = 1.0; /* then scale = 1.0 */
249 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
252 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
253 if (kernel == (KernelInfo *)NULL)
255 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
256 kernel->type = UserDefinedKernel;
257 kernel->signature = MagickSignature;
259 /* Has a ':' in argument - New user kernel specification */
260 p = strchr(kernel_string, ':');
261 if ( p != (char *) NULL)
263 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
264 memcpy(token, kernel_string, (size_t) (p-kernel_string));
265 token[p-kernel_string] = '\0';
266 flags = ParseGeometry(token, &args);
268 /* Size handling and checks of geometry settings */
269 if ( (flags & WidthValue) == 0 ) /* if no width then */
270 args.rho = args.sigma; /* then width = height */
271 if ( args.rho < 1.0 ) /* if width too small */
272 args.rho = 1.0; /* then width = 1 */
273 if ( args.sigma < 1.0 ) /* if height too small */
274 args.sigma = args.rho; /* then height = width */
275 kernel->width = (unsigned long)args.rho;
276 kernel->height = (unsigned long)args.sigma;
278 /* Offset Handling and Checks */
279 if ( args.xi < 0.0 || args.psi < 0.0 )
280 return(DestroyKernelInfo(kernel));
281 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
282 : (long) (kernel->width-1)/2;
283 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
284 : (long) (kernel->height-1)/2;
285 if ( kernel->x >= (long) kernel->width ||
286 kernel->y >= (long) kernel->height )
287 return(DestroyKernelInfo(kernel));
289 p++; /* advance beyond the ':' */
292 { /* ELSE - Old old kernel specification, forming odd-square kernel */
293 /* count up number of values given */
294 p=(const char *) kernel_string;
295 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
296 p++; /* ignore "'" chars for convolve filter usage - Cristy */
297 for (i=0; *p != '\0'; i++)
299 GetMagickToken(p,&p,token);
301 GetMagickToken(p,&p,token);
303 /* set the size of the kernel - old sized square */
304 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
305 kernel->x = kernel->y = (long) (kernel->width-1)/2;
306 p=(const char *) kernel_string;
307 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
308 p++; /* ignore "'" chars for convolve filter usage - Cristy */
311 /* Read in the kernel values from rest of input string argument */
312 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
313 kernel->height*sizeof(double));
314 if (kernel->values == (double *) NULL)
315 return(DestroyKernelInfo(kernel));
317 kernel->minimum = +MagickHuge;
318 kernel->maximum = -MagickHuge;
319 kernel->negative_range = kernel->positive_range = 0.0;
320 for (i=0; (i < (long) (kernel->width*kernel->height)) && (*p != '\0'); i++)
322 GetMagickToken(p,&p,token);
324 GetMagickToken(p,&p,token);
325 if ( LocaleCompare("nan",token) == 0
326 || LocaleCompare("-",token) == 0 ) {
327 kernel->values[i] = nan; /* do not include this value in kernel */
330 kernel->values[i] = StringToDouble(token);
331 ( kernel->values[i] < 0)
332 ? ( kernel->negative_range += kernel->values[i] )
333 : ( kernel->positive_range += kernel->values[i] );
334 Minimize(kernel->minimum, kernel->values[i]);
335 Maximize(kernel->maximum, kernel->values[i]);
338 /* check that we recieved at least one real (non-nan) value! */
339 if ( kernel->minimum == MagickHuge )
340 return(DestroyKernelInfo(kernel));
342 /* This should not be needed for a fully defined kernel
343 * Perhaps an error should be reported instead!
344 * Kept for backward compatibility.
346 if ( i < (long) (kernel->width*kernel->height) ) {
347 Minimize(kernel->minimum, kernel->values[i]);
348 Maximize(kernel->maximum, kernel->values[i]);
349 for ( ; i < (long) (kernel->width*kernel->height); i++)
350 kernel->values[i]=0.0;
357 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
361 % A c q u i r e K e r n e l B u i l t I n %
365 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
367 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
368 % kernels used for special purposes such as gaussian blurring, skeleton
369 % pruning, and edge distance determination.
371 % They take a KernelType, and a set of geometry style arguments, which were
372 % typically decoded from a user supplied string, or from a more complex
373 % Morphology Method that was requested.
375 % The format of the AcquireKernalBuiltIn method is:
377 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
378 % const GeometryInfo args)
380 % A description of each parameter follows:
382 % o type: the pre-defined type of kernel wanted
384 % o args: arguments defining or modifying the kernel
386 % Convolution Kernels
388 % Gaussian "{radius},{sigma}"
389 % Generate a two-dimentional gaussian kernel, as used by -gaussian
390 % A sigma is required, (with the 'x'), due to historical reasons.
392 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
393 % the final size of the resulting kernel to a square 2*radius+1 in size.
394 % The radius should be at least 2 times that of the sigma value, or
395 % sever clipping and aliasing may result. If not given or set to 0 the
396 % radius will be determined so as to produce the best minimal error
397 % result, which is usally much larger than is normally needed.
399 % Blur "{radius},{sigma},{angle}"
400 % As per Gaussian, but generates a 1 dimensional or linear gaussian
401 % blur, at the angle given (current restricted to orthogonal angles).
402 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
404 % NOTE that two such blurs perpendicular to each other is equivelent to
405 % -blur and the previous gaussian, but is often 10 or more times faster.
407 % Comet "{width},{sigma},{angle}"
408 % Blur in one direction only, mush like how a bright object leaves
409 % a comet like trail. The Kernel is actually half a gaussian curve,
410 % Adding two such blurs in oppiste directions produces a Linear Blur.
412 % NOTE: that the first argument is the width of the kernel and not the
413 % radius of the kernel.
415 % # Still to be implemented...
417 % # Sharpen "{radius},{sigma}
418 % # Negated Gaussian (center zeroed and re-normalized),
419 % # with a 2 unit positive peak. -- Check On line documentation
421 % # Laplacian "{radius},{sigma}"
422 % # Laplacian (a mexican hat like) Function
424 % # LOG "{radius},{sigma1},{sigma2}
425 % # Laplacian of Gaussian
427 % # DOG "{radius},{sigma1},{sigma2}
428 % # Difference of two Gaussians
432 % # Set kernel values using a resize filter, and given scale (sigma)
433 % # Cylindrical or Linear. Is this posible with an image?
438 % Rectangle "{geometry}"
439 % Simply generate a rectangle of 1's with the size given. You can also
440 % specify the location of the 'control point', otherwise the closest
441 % pixel to the center of the rectangle is selected.
443 % Properly centered and odd sized rectangles work the best.
445 % Diamond "[{radius}[,{scale}]]"
446 % Generate a diamond shaped kernal with given radius to the points.
447 % Kernel size will again be radius*2+1 square and defaults to radius 1,
448 % generating a 3x3 kernel that is slightly larger than a square.
450 % Square "[{radius}[,{scale}]]"
451 % Generate a square shaped kernel of size radius*2+1, and defaulting
452 % to a 3x3 (radius 1).
454 % Note that using a larger radius for the "Square" or the "Diamond"
455 % is also equivelent to iterating the basic morphological method
456 % that many times. However However iterating with the smaller radius 1
457 % default is actually faster than using a larger kernel radius.
459 % Disk "[{radius}[,{scale}]]
460 % Generate a binary disk of the radius given, radius may be a float.
461 % Kernel size will be ceil(radius)*2+1 square.
462 % NOTE: Here are some disk shapes of specific interest
463 % "disk:1" => "diamond" or "cross:1"
464 % "disk:1.5" => "square"
465 % "disk:2" => "diamond:2"
466 % "disk:2.5" => a general disk shape of radius 2
467 % "disk:2.9" => "square:2"
468 % "disk:3.5" => default - octagonal/disk shape of radius 3
469 % "disk:4.2" => roughly octagonal shape of radius 4
470 % "disk:4.3" => a general disk shape of radius 4
471 % After this all the kernel shape becomes more and more circular.
473 % Because a "disk" is more circular when using a larger radius, using a
474 % larger radius is preferred over iterating the morphological operation.
476 % Plus "[{radius}[,{scale}]]"
477 % Generate a kernel in the shape of a 'plus' sign. The length of each
478 % arm is also the radius, which defaults to 2.
480 % This kernel is not a good general morphological kernel, but is used
481 % more for highlighting and marking any single pixels in an image using,
482 % a "Dilate" or "Erode" method as appropriate.
484 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
486 % Note that unlike other kernels iterating a plus does not produce the
487 % same result as using a larger radius for the cross.
489 % Distance Measuring Kernels
491 % Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1)
492 % Manhatten "[{radius}][x{scale}]" square grid distance (default r=1)
493 % Euclidean "[{radius}][x{scale}]" direct distance (default r=1)
495 % Different types of distance measuring methods, which are used with the
496 % a 'Distance' morphology method for generating a gradient based on
497 % distance from an edge of a binary shape, though there is a technique
498 % for handling a anti-aliased shape.
500 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
501 % one to any neighbour, orthogonal or diagonal. One why of thinking of
502 % it is the number of squares a 'King' or 'Queen' in chess needs to
503 % traverse reach any other position on a chess board. It results in a
504 % 'square' like distance function, but one where diagonals are closer
507 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
508 % Cab metric), is the distance needed when you can only travel in
509 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
510 % in chess would travel. It results in a diamond like distances, where
511 % diagonals are further than expected.
513 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
514 % However by default the kernel size only has a radius of 1, which
515 % limits the distance to 'Knight' like moves, with only orthogonal and
516 % diagonal measurements being correct. As such for the default kernel
517 % you will get octagonal like distance function, which is reasonally
520 % However if you use a larger radius such as "Euclidean:4" you will
521 % get a much smoother distance gradient from the edge of the shape.
522 % Of course a larger kernel is slower to use, and generally not needed.
524 % To allow the use of fractional distances that you get with diagonals
525 % the actual distance is scaled by a fixed value which the user can
526 % provide. This is not actually nessary for either ""Chebyshev" or
527 % "Manhatten" distance kernels, but is done for all three distance
528 % kernels. If no scale is provided it is set to a value of 100,
529 % allowing for a maximum distance measurement of 655 pixels using a Q16
530 % version of IM, from any edge. However for small images this can
531 % result in quite a dark gradient.
533 % See the 'Distance' Morphological Method, for information of how it is
536 % # Hit-n-Miss Kernel-Lists -- Still to be implemented
538 % # specifically for Pruning, Thinning, Thickening
542 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
543 const GeometryInfo *args)
556 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
558 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
559 if (kernel == (KernelInfo *) NULL)
561 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
562 kernel->minimum = kernel->maximum = 0.0;
563 kernel->negative_range = kernel->positive_range = 0.0;
565 kernel->signature = MagickSignature;
568 /* Convolution Kernels */
571 sigma = fabs(args->sigma);
573 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
575 kernel->width = kernel->height =
576 GetOptimalKernelWidth2D(args->rho,sigma);
577 kernel->x = kernel->y = (long) (kernel->width-1)/2;
578 kernel->negative_range = kernel->positive_range = 0.0;
579 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
580 kernel->height*sizeof(double));
581 if (kernel->values == (double *) NULL)
582 return(DestroyKernelInfo(kernel));
584 sigma = 2.0*sigma*sigma; /* simplify the expression */
585 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
586 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
587 kernel->positive_range += (
589 exp(-((double)(u*u+v*v))/sigma)
590 /* / (MagickPI*sigma) */ );
592 kernel->maximum = kernel->values[
593 kernel->y*kernel->width+kernel->x ];
595 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
601 sigma = fabs(args->sigma);
603 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
605 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
606 kernel->x = (long) (kernel->width-1)/2;
609 kernel->negative_range = kernel->positive_range = 0.0;
610 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
611 kernel->height*sizeof(double));
612 if (kernel->values == (double *) NULL)
613 return(DestroyKernelInfo(kernel));
617 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
618 ** It generates a gaussian 3 times the width, and compresses it into
619 ** the expected range. This produces a closer normalization of the
620 ** resulting kernel, especially for very low sigma values.
621 ** As such while wierd it is prefered.
623 ** I am told this method originally came from Photoshop.
625 sigma *= KernelRank; /* simplify expanded curve */
626 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
627 (void) ResetMagickMemory(kernel->values,0, (size_t)
628 kernel->width*sizeof(double));
629 for ( u=-v; u <= v; u++) {
630 kernel->values[(u+v)/KernelRank] +=
631 exp(-((double)(u*u))/(2.0*sigma*sigma))
632 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
634 for (i=0; i < (long) kernel->width; i++)
635 kernel->positive_range += kernel->values[i];
637 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
638 kernel->positive_range += (
640 exp(-((double)(u*u))/(2.0*sigma*sigma))
641 /* / (MagickSQ2PI*sigma) */ );
644 kernel->maximum = kernel->values[ kernel->x ];
645 /* Note that neither methods above generate a normalized kernel,
646 ** though it gets close. The kernel may be 'clipped' by a user defined
647 ** radius, producing a smaller (darker) kernel. Also for very small
648 ** sigma's (> 0.1) the central value becomes larger than one, and thus
649 ** producing a very bright kernel.
652 /* Normalize the 1D Gaussian Kernel
654 ** Because of this the divisor in the above kernel generator is
655 ** not needed, so is not done above.
657 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
659 /* rotate the kernel by given angle */
660 RotateKernelInfo(kernel, args->xi);
665 sigma = fabs(args->sigma);
667 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
669 if ( args->rho < 1.0 )
670 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
672 kernel->width = (unsigned long)args->rho;
673 kernel->x = kernel->y = 0;
675 kernel->negative_range = kernel->positive_range = 0.0;
676 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
677 kernel->height*sizeof(double));
678 if (kernel->values == (double *) NULL)
679 return(DestroyKernelInfo(kernel));
681 /* A comet blur is half a gaussian curve, so that the object is
682 ** blurred in one direction only. This may not be quite the right
683 ** curve so may change in the future. The function must be normalised.
687 sigma *= KernelRank; /* simplify expanded curve */
688 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
689 (void) ResetMagickMemory(kernel->values,0, (size_t)
690 kernel->width*sizeof(double));
691 for ( u=0; u < v; u++) {
692 kernel->values[u/KernelRank] +=
693 exp(-((double)(u*u))/(2.0*sigma*sigma))
694 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
696 for (i=0; i < (long) kernel->width; i++)
697 kernel->positive_range += kernel->values[i];
699 for ( i=0; i < (long) kernel->width; i++)
700 kernel->positive_range += (
702 exp(-((double)(i*i))/(2.0*sigma*sigma))
703 /* / (MagickSQ2PI*sigma) */ );
706 kernel->maximum = kernel->values[0];
708 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
709 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
712 /* Boolean Kernels */
713 case RectangleKernel:
717 if ( type == SquareKernel )
720 kernel->width = kernel->height = 3; /* default radius = 1 */
722 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
723 kernel->x = kernel->y = (long) (kernel->width-1)/2;
727 /* NOTE: user defaults set in "AcquireKernelInfo()" */
728 if ( args->rho < 1.0 || args->sigma < 1.0 )
729 return(DestroyKernelInfo(kernel)); /* invalid args given */
730 kernel->width = (unsigned long)args->rho;
731 kernel->height = (unsigned long)args->sigma;
732 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
733 args->psi < 0.0 || args->psi > (double)kernel->height )
734 return(DestroyKernelInfo(kernel)); /* invalid args given */
735 kernel->x = (long) args->xi;
736 kernel->y = (long) args->psi;
739 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
740 kernel->height*sizeof(double));
741 if (kernel->values == (double *) NULL)
742 return(DestroyKernelInfo(kernel));
744 /* set all kernel values to 1.0 */
745 u=(long) kernel->width*kernel->height;
746 for ( i=0; i < u; i++)
747 kernel->values[i] = scale;
748 kernel->minimum = kernel->maximum = scale; /* a flat shape */
749 kernel->positive_range = scale*u;
755 kernel->width = kernel->height = 3; /* default radius = 1 */
757 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
758 kernel->x = kernel->y = (long) (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 /* set all kernel values within diamond area to scale given */
766 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
767 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
768 if ((labs(u)+labs(v)) <= (long)kernel->x)
769 kernel->positive_range += kernel->values[i] = args->sigma;
771 kernel->values[i] = nan;
772 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
780 limit = (long)(args->rho*args->rho);
781 if (args->rho < 0.1) /* default radius approx 3.5 */
782 kernel->width = kernel->height = 7L, limit = 10L;
784 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
785 kernel->x = kernel->y = (long) (kernel->width-1)/2;
787 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
788 kernel->height*sizeof(double));
789 if (kernel->values == (double *) NULL)
790 return(DestroyKernelInfo(kernel));
792 /* set all kernel values within disk area to 1.0 */
793 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
794 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
795 if ((u*u+v*v) <= limit)
796 kernel->positive_range += kernel->values[i] = args->sigma;
798 kernel->values[i] = nan;
799 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
805 kernel->width = kernel->height = 5; /* default radius 2 */
807 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
808 kernel->x = kernel->y = (long) (kernel->width-1)/2;
810 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
811 kernel->height*sizeof(double));
812 if (kernel->values == (double *) NULL)
813 return(DestroyKernelInfo(kernel));
815 /* set all kernel values along axises to 1.0 */
816 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
817 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
818 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
819 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
820 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
823 /* Distance Measuring Kernels */
824 case ChebyshevKernel:
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)) ? labs(u) : labs(v)) );
845 kernel->maximum = kernel->values[0];
848 case ManhattenKernel:
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*(labs(u)+labs(v)) );
869 kernel->maximum = kernel->values[0];
872 case EuclideanKernel:
878 kernel->width = kernel->height = 3; /* default radius = 1 */
880 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
881 kernel->x = kernel->y = (long) (kernel->width-1)/2;
883 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
884 kernel->height*sizeof(double));
885 if (kernel->values == (double *) NULL)
886 return(DestroyKernelInfo(kernel));
888 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
889 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
890 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
891 kernel->positive_range += ( kernel->values[i] =
892 scale*sqrt((double)(u*u+v*v)) );
893 kernel->maximum = kernel->values[0];
896 /* Undefined Kernels */
897 case LaplacianKernel:
900 perror("Kernel Type has not been defined yet");
903 /* Generate a No-Op minimal kernel - 1x1 pixel */
904 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
905 if (kernel->values == (double *) NULL)
906 return(DestroyKernelInfo(kernel));
907 kernel->width = kernel->height = 1;
908 kernel->x = kernel->x = 0;
909 kernel->type = UndefinedKernel;
911 kernel->positive_range =
912 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
920 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
924 % C l o n e K e r n e l I n f o %
928 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
930 % CloneKernelInfo() creates a new clone of the given Kernel so that its can
931 % be modified without effecting the original. The cloned kernel should be
932 % destroyed using DestoryKernelInfo() when no longer needed.
934 % The format of the DestroyKernelInfo method is:
936 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
938 % A description of each parameter follows:
940 % o kernel: the Morphology/Convolution kernel to be cloned
943 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
951 assert(kernel != (KernelInfo *) NULL);
953 new=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
954 if (new == (KernelInfo *) NULL)
956 *new = *kernel; /* copy values in structure */
958 new->values=(double *) AcquireQuantumMemory(kernel->width,
959 kernel->height*sizeof(double));
960 if (new->values == (double *) NULL)
961 return(DestroyKernelInfo(new));
963 for (i=0; i < (long) (kernel->width*kernel->height); i++)
964 new->values[i] = kernel->values[i];
970 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
974 % D e s t r o y K e r n e l I n f o %
978 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
980 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
983 % The format of the DestroyKernelInfo method is:
985 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
987 % A description of each parameter follows:
989 % o kernel: the Morphology/Convolution kernel to be destroyed
993 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
995 assert(kernel != (KernelInfo *) NULL);
997 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
998 kernel->height*sizeof(double));
999 kernel->values=(double *)RelinquishMagickMemory(kernel->values);
1000 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
1005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1009 % M o r p h o l o g y I m a g e C h a n n e l %
1013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1015 % MorphologyImageChannel() applies a user supplied kernel to the image
1016 % according to the given mophology method.
1018 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1019 % by the kernel generator.
1021 % The format of the MorphologyImage method is:
1023 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
1024 % const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
1025 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1026 % channel,MorphologyMethod method,const long iterations,
1027 % KernelInfo *kernel,ExceptionInfo *exception)
1029 % A description of each parameter follows:
1031 % o image: the image.
1033 % o method: the morphology method to be applied.
1035 % o iterations: apply the operation this many times (or no change).
1036 % A value of -1 means loop until no change found.
1037 % How this is applied may depend on the morphology method.
1038 % Typically this is a value of 1.
1040 % o channel: the channel type.
1042 % o kernel: An array of double representing the morphology kernel.
1043 % Warning: kernel may be normalized for the Convolve method.
1045 % o exception: return any errors or warnings in this structure.
1048 % TODO: bias and auto-scale handling of the kernel for convolution
1049 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1050 % by the kernel generator.
1055 /* Internal function
1056 * Apply the Low-Level Morphology Method using the given Kernel
1057 * Returning the number of pixels that changed.
1058 * Two pre-created images must be provided, no image is created.
1060 static unsigned long MorphologyApply(const Image *image, Image
1061 *result_image, const MorphologyMethod method, const ChannelType channel,
1062 const KernelInfo *kernel, ExceptionInfo *exception)
1064 #define MorphologyTag "Morphology/Image"
1081 /* Only the most basic morphology is actually performed by this routine */
1084 Apply Basic Morphology to Image.
1090 GetMagickPixelPacket(image,&bias);
1091 SetMagickPixelPacketBias(image,&bias);
1092 /* Future: handle auto-bias from user, based on kernel input */
1094 p_view=AcquireCacheView(image);
1095 q_view=AcquireCacheView(result_image);
1097 /* Some methods (including convolve) needs use a reflected kernel.
1098 * Adjust 'origin' offsets for this reflected kernel.
1103 case ErodeMorphology:
1104 case ErodeIntensityMorphology:
1105 /* kernel is user as is, without reflection */
1107 case ConvolveMorphology:
1108 case DilateMorphology:
1109 case DilateIntensityMorphology:
1110 case DistanceMorphology:
1111 /* kernel needs to used with reflection */
1112 offx = (long) kernel->width-offx-1;
1113 offy = (long) kernel->height-offy-1;
1116 perror("Not a low level Morpholgy Method");
1120 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1121 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1123 for (y=0; y < (long) image->rows; y++)
1128 register const PixelPacket
1131 register const IndexPacket
1132 *restrict p_indexes;
1134 register PixelPacket
1137 register IndexPacket
1138 *restrict q_indexes;
1146 if (status == MagickFalse)
1148 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1149 image->columns+kernel->width, kernel->height, exception);
1150 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1152 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1157 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1158 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1159 r = (image->columns+kernel->width)*offy+offx; /* constant */
1161 for (x=0; x < (long) image->columns; x++)
1169 register const double
1172 register const PixelPacket
1175 register const IndexPacket
1176 *restrict k_indexes;
1181 /* Copy input to ouput image for unused channels
1182 * This removes need for 'cloning' a new image every iteration
1185 if (image->colorspace == CMYKColorspace)
1186 q_indexes[x] = p_indexes[r];
1188 result.green=(MagickRealType) 0;
1189 result.blue=(MagickRealType) 0;
1190 result.opacity=(MagickRealType) 0;
1191 result.index=(MagickRealType) 0;
1193 case ConvolveMorphology:
1194 /* Set the user defined bias of the weighted average output
1196 ** FUTURE: provide some way for internal functions to disable
1197 ** user defined bias and scaling effects.
1201 case DilateMorphology:
1206 result.index = -MagickHuge;
1208 case ErodeMorphology:
1213 result.index = +MagickHuge;
1215 case DilateIntensityMorphology:
1216 case ErodeIntensityMorphology:
1217 result.red = 0.0; /* flag indicating first match found */
1220 /* Otherwise just start with the original pixel value */
1221 result.red = (MagickRealType) p[r].red;
1222 result.green = (MagickRealType) p[r].green;
1223 result.blue = (MagickRealType) p[r].blue;
1224 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1225 if ( image->colorspace == CMYKColorspace)
1226 result.index = (MagickRealType) p_indexes[r];
1231 case ConvolveMorphology:
1232 /* Weighted Average of pixels using reflected kernel
1234 ** NOTE for correct working of this operation for asymetrical
1235 ** kernels, the kernel needs to be applied in its reflected form.
1236 ** That is its values needs to be reversed.
1238 ** Correlation is actually the same as this but without reflecting
1239 ** the kernel, and thus 'lower-level' that Convolution. However
1240 ** as Convolution is the more common method used, and it does not
1241 ** really cost us much in terms of processing to use a reflected
1242 ** kernel it is Convolution that is implemented.
1244 ** Correlation will have its kernel reflected before calling
1245 ** this function to do a Convolve.
1247 ** For more details of Correlation vs Convolution see
1248 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
1250 if (((channel & OpacityChannel) == 0) ||
1251 (image->matte == MagickFalse))
1253 /* Convolution without transparency effects */
1254 k = &kernel->values[ kernel->width*kernel->height-1 ];
1256 k_indexes = p_indexes;
1257 for (v=0; v < (long) kernel->height; v++) {
1258 for (u=0; u < (long) kernel->width; u++, k--) {
1259 if ( IsNan(*k) ) continue;
1260 result.red += (*k)*k_pixels[u].red;
1261 result.green += (*k)*k_pixels[u].green;
1262 result.blue += (*k)*k_pixels[u].blue;
1263 /* result.opacity += not involved here */
1264 if ( image->colorspace == CMYKColorspace)
1265 result.index += (*k)*k_indexes[u];
1267 k_pixels += image->columns+kernel->width;
1268 k_indexes += image->columns+kernel->width;
1272 { /* Kernel & Alpha weighted Convolution */
1274 alpha, /* alpha value * kernel weighting */
1275 gamma; /* weighting divisor */
1278 k = &kernel->values[ kernel->width*kernel->height-1 ];
1280 k_indexes = p_indexes;
1281 for (v=0; v < (long) kernel->height; v++) {
1282 for (u=0; u < (long) kernel->width; u++, k--) {
1283 if ( IsNan(*k) ) continue;
1284 alpha=(*k)*(QuantumScale*(QuantumRange-
1285 k_pixels[u].opacity));
1287 result.red += alpha*k_pixels[u].red;
1288 result.green += alpha*k_pixels[u].green;
1289 result.blue += alpha*k_pixels[u].blue;
1290 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1291 if ( image->colorspace == CMYKColorspace)
1292 result.index += alpha*k_indexes[u];
1294 k_pixels += image->columns+kernel->width;
1295 k_indexes += image->columns+kernel->width;
1297 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1298 result.red *= gamma;
1299 result.green *= gamma;
1300 result.blue *= gamma;
1301 result.opacity *= gamma;
1302 result.index *= gamma;
1306 case ErodeMorphology:
1307 /* Minimize Value within kernel neighbourhood
1309 ** NOTE that the kernel is not reflected for this operation!
1311 ** NOTE: in normal Greyscale Morphology, the kernel value should
1312 ** be added to the real value, this is currently not done, due to
1313 ** the nature of the boolean kernels being used.
1317 k_indexes = p_indexes;
1318 for (v=0; v < (long) kernel->height; v++) {
1319 for (u=0; u < (long) kernel->width; u++, k++) {
1320 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1321 Minimize(result.red, (double) k_pixels[u].red);
1322 Minimize(result.green, (double) k_pixels[u].green);
1323 Minimize(result.blue, (double) k_pixels[u].blue);
1324 Minimize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1325 if ( image->colorspace == CMYKColorspace)
1326 Minimize(result.index, (double) k_indexes[u]);
1328 k_pixels += image->columns+kernel->width;
1329 k_indexes += image->columns+kernel->width;
1333 case DilateMorphology:
1334 /* Maximize Value within kernel neighbourhood
1336 ** NOTE for correct working of this operation for asymetrical
1337 ** kernels, the kernel needs to be applied in its reflected form.
1338 ** That is its values needs to be reversed.
1340 ** NOTE: in normal Greyscale Morphology, the kernel value should
1341 ** be added to the real value, this is currently not done, due to
1342 ** the nature of the boolean kernels being used.
1345 k = &kernel->values[ kernel->width*kernel->height-1 ];
1347 k_indexes = p_indexes;
1348 for (v=0; v < (long) kernel->height; v++) {
1349 for (u=0; u < (long) kernel->width; u++, k--) {
1350 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1351 Maximize(result.red, (double) k_pixels[u].red);
1352 Maximize(result.green, (double) k_pixels[u].green);
1353 Maximize(result.blue, (double) k_pixels[u].blue);
1354 Maximize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1355 if ( image->colorspace == CMYKColorspace)
1356 Maximize(result.index, (double) k_indexes[u]);
1358 k_pixels += image->columns+kernel->width;
1359 k_indexes += image->columns+kernel->width;
1363 case ErodeIntensityMorphology:
1364 /* Select Pixel with Minimum Intensity within kernel neighbourhood
1366 ** WARNING: the intensity test fails for CMYK and does not
1367 ** take into account the moderating effect of teh alpha channel
1368 ** on the intensity.
1370 ** NOTE that the kernel is not reflected for this operation!
1374 k_indexes = p_indexes;
1375 for (v=0; v < (long) kernel->height; v++) {
1376 for (u=0; u < (long) kernel->width; u++, k++) {
1377 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1378 if ( result.red == 0.0 ||
1379 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1380 /* copy the whole pixel - no channel selection */
1382 if ( result.red > 0.0 ) changed++;
1386 k_pixels += image->columns+kernel->width;
1387 k_indexes += image->columns+kernel->width;
1391 case DilateIntensityMorphology:
1392 /* Select Pixel with Maximum Intensity within kernel neighbourhood
1394 ** WARNING: the intensity test fails for CMYK and does not
1395 ** take into account the moderating effect of teh alpha channel
1396 ** on the intensity.
1398 ** NOTE for correct working of this operation for asymetrical
1399 ** kernels, the kernel needs to be applied in its reflected form.
1400 ** That is its values needs to be reversed.
1402 k = &kernel->values[ kernel->width*kernel->height-1 ];
1404 k_indexes = p_indexes;
1405 for (v=0; v < (long) kernel->height; v++) {
1406 for (u=0; u < (long) kernel->width; u++, k--) {
1407 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1408 if ( result.red == 0.0 ||
1409 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1410 /* copy the whole pixel - no channel selection */
1412 if ( result.red > 0.0 ) changed++;
1416 k_pixels += image->columns+kernel->width;
1417 k_indexes += image->columns+kernel->width;
1421 case DistanceMorphology:
1422 /* Add kernel Value and select the minimum value found.
1423 ** The result is a iterative distance from edge of image shape.
1425 ** All Distance Kernels are symetrical, but that may not always
1426 ** be the case. For example how about a distance from left edges?
1427 ** To work correctly with asymetrical kernels the reflected kernel
1428 ** needs to be applied.
1431 /* No need to do distance morphology if original value is zero
1432 ** Unfortunatally I have not been able to get this right
1433 ** when channel selection also becomes involved. -- Arrgghhh
1435 if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1436 || ((channel & GreenChannel) == 0 && p[r].green == 0)
1437 || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1438 || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1439 || (( (channel & IndexChannel) == 0
1440 || image->colorspace != CMYKColorspace
1441 ) && p_indexes[x] ==0 )
1445 k = &kernel->values[ kernel->width*kernel->height-1 ];
1447 k_indexes = p_indexes;
1448 for (v=0; v < (long) kernel->height; v++) {
1449 for (u=0; u < (long) kernel->width; u++, k--) {
1450 if ( IsNan(*k) ) continue;
1451 Minimize(result.red, (*k)+k_pixels[u].red);
1452 Minimize(result.green, (*k)+k_pixels[u].green);
1453 Minimize(result.blue, (*k)+k_pixels[u].blue);
1454 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1455 if ( image->colorspace == CMYKColorspace)
1456 Minimize(result.index, (*k)+k_indexes[u]);
1458 k_pixels += image->columns+kernel->width;
1459 k_indexes += image->columns+kernel->width;
1463 case UndefinedMorphology:
1465 break; /* Do nothing */
1468 case UndefinedMorphology:
1469 case DilateIntensityMorphology:
1470 case ErodeIntensityMorphology:
1471 break; /* full pixel was directly assigned - not a channel method */
1473 /* Assign the results */
1474 if ((channel & RedChannel) != 0)
1475 q->red = ClampToQuantum(result.red);
1476 if ((channel & GreenChannel) != 0)
1477 q->green = ClampToQuantum(result.green);
1478 if ((channel & BlueChannel) != 0)
1479 q->blue = ClampToQuantum(result.blue);
1480 if ((channel & OpacityChannel) != 0
1481 && image->matte == MagickTrue )
1482 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1483 if ((channel & IndexChannel) != 0
1484 && image->colorspace == CMYKColorspace)
1485 q_indexes[x] = ClampToQuantum(result.index);
1488 if ( ( p[r].red != q->red )
1489 || ( p[r].green != q->green )
1490 || ( p[r].blue != q->blue )
1491 || ( p[r].opacity != q->opacity )
1492 || ( image->colorspace == CMYKColorspace &&
1493 p_indexes[r] != q_indexes[x] ) )
1494 changed++; /* The pixel had some value changed! */
1498 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1499 if (sync == MagickFalse)
1501 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1506 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1507 #pragma omp critical (MagickCore_MorphologyImage)
1509 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1510 if (proceed == MagickFalse)
1514 result_image->type=image->type;
1515 q_view=DestroyCacheView(q_view);
1516 p_view=DestroyCacheView(p_view);
1517 return(status ? (unsigned long) changed : 0);
1521 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1522 method, const long iterations,const KernelInfo *kernel, ExceptionInfo
1528 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1529 iterations,kernel,exception);
1530 return(morphology_image);
1534 MagickExport Image *MorphologyImageChannel(const Image *image,
1535 const ChannelType channel,const MorphologyMethod method,
1536 const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
1559 assert(image != (Image *) NULL);
1560 assert(image->signature == MagickSignature);
1561 assert(kernel != (KernelInfo *) NULL);
1562 assert(kernel->signature == MagickSignature);
1563 assert(exception != (ExceptionInfo *) NULL);
1564 assert(exception->signature == MagickSignature);
1566 if ( iterations == 0 )
1567 return((Image *)NULL); /* null operation - nothing to do! */
1569 /* kernel must be valid at this point
1570 * (except maybe for posible future morphology methods like "Prune"
1572 assert(kernel != (KernelInfo *)NULL);
1574 count = 0; /* interation count */
1575 changed = 1; /* if compound method assume image was changed */
1576 curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
1577 curr_method = method; /* to be changed as nessary */
1579 limit = (unsigned long) iterations;
1580 if ( iterations < 0 )
1581 limit = image->columns > image->rows ? image->columns : image->rows;
1583 /* Third-level morphology methods */
1584 grad_image=(Image *) NULL;
1585 switch( curr_method ) {
1586 case EdgeMorphology:
1587 grad_image = MorphologyImageChannel(image, channel,
1588 DilateMorphology, iterations, curr_kernel, exception);
1590 case EdgeInMorphology:
1591 curr_method = ErodeMorphology;
1593 case EdgeOutMorphology:
1594 curr_method = DilateMorphology;
1596 case TopHatMorphology:
1597 curr_method = OpenMorphology;
1599 case BottomHatMorphology:
1600 curr_method = CloseMorphology;
1603 break; /* not a third-level method */
1606 /* Second-level morphology methods */
1607 switch( curr_method ) {
1608 case OpenMorphology:
1609 /* Open is a Erode then a Dilate without reflection */
1610 new_image = MorphologyImageChannel(image, channel,
1611 ErodeMorphology, iterations, curr_kernel, exception);
1612 if (new_image == (Image *) NULL)
1613 return((Image *) NULL);
1614 curr_method = DilateMorphology;
1616 case OpenIntensityMorphology:
1617 new_image = MorphologyImageChannel(image, channel,
1618 ErodeIntensityMorphology, iterations, curr_kernel, exception);
1619 if (new_image == (Image *) NULL)
1620 return((Image *) NULL);
1621 curr_method = DilateIntensityMorphology;
1624 case CloseMorphology:
1625 /* Close is a Dilate then Erode using reflected kernel */
1626 /* A reflected kernel is needed for a Close */
1627 if ( curr_kernel == kernel )
1628 curr_kernel = CloneKernelInfo(kernel);
1629 RotateKernelInfo(curr_kernel,180);
1630 new_image = MorphologyImageChannel(image, channel,
1631 DilateMorphology, iterations, curr_kernel, exception);
1632 if (new_image == (Image *) NULL)
1633 return((Image *) NULL);
1634 curr_method = ErodeMorphology;
1636 case CloseIntensityMorphology:
1637 /* A reflected kernel is needed for a Close */
1638 if ( curr_kernel == kernel )
1639 curr_kernel = CloneKernelInfo(kernel);
1640 RotateKernelInfo(curr_kernel,180);
1641 new_image = MorphologyImageChannel(image, channel,
1642 DilateIntensityMorphology, iterations, curr_kernel, exception);
1643 if (new_image == (Image *) NULL)
1644 return((Image *) NULL);
1645 curr_method = ErodeIntensityMorphology;
1648 case CorrelateMorphology:
1649 /* A Correlation is actually a Convolution with a reflected kernel.
1650 ** However a Convolution is a weighted sum with a reflected kernel.
1651 ** It may seem stange to convert a Correlation into a Convolution
1652 ** as the Correleation is the simplier method, but Convolution is
1653 ** much more commonly used, and it makes sense to implement it directly
1654 ** so as to avoid the need to duplicate the kernel when it is not
1655 ** required (which is typically the default).
1657 if ( curr_kernel == kernel )
1658 curr_kernel = CloneKernelInfo(kernel);
1659 RotateKernelInfo(curr_kernel,180);
1660 curr_method = ConvolveMorphology;
1661 /* FALL-THRU into Correlate (weigthed sum without reflection) */
1663 case ConvolveMorphology:
1664 /* Scale or Normalize kernel, according to user wishes
1665 ** before using it for the Convolve/Correlate method.
1667 ** FUTURE: provide some way for internal functions to disable
1668 ** user bias and scaling effects.
1670 artifact = GetImageArtifact(image,"convolve:scale");
1671 if ( artifact != (char *)NULL ) {
1677 if ( curr_kernel == kernel )
1678 curr_kernel = CloneKernelInfo(kernel);
1681 flags = ParseGeometry(artifact, &args);
1682 ScaleKernelInfo(curr_kernel, args.rho, flags);
1684 /* FALL-THRU to do the first, and typically the only iteration */
1687 /* Do a single iteration using the Low-Level Morphology method!
1688 ** This ensures a "new_image" has been generated, but allows us to skip
1689 ** the creation of 'old_image' if no more iterations are needed.
1691 ** The "curr_method" should also be set to a low-level method that is
1692 ** understood by the MorphologyApply() internal function.
1694 new_image=CloneImage(image,0,0,MagickTrue,exception);
1695 if (new_image == (Image *) NULL)
1696 return((Image *) NULL);
1697 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1699 InheritException(exception,&new_image->exception);
1700 new_image=DestroyImage(new_image);
1701 return((Image *) NULL);
1703 changed = MorphologyApply(image,new_image,curr_method,channel,curr_kernel,
1706 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1707 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1708 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1713 /* At this point the "curr_method" should not only be set to a low-level
1714 ** method that is understood by the MorphologyApply() internal function,
1715 ** but "new_image" should now be defined, as the image to apply the
1716 ** "curr_method" to.
1719 /* Repeat the low-level morphology until count or no change reached */
1720 if ( count < (long) limit && changed > 0 ) {
1721 old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1722 if (old_image == (Image *) NULL)
1723 return(DestroyImage(new_image));
1724 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1726 InheritException(exception,&old_image->exception);
1727 old_image=DestroyImage(old_image);
1728 return(DestroyImage(new_image));
1730 while( count < (long) limit && changed != 0 )
1732 Image *tmp = old_image;
1733 old_image = new_image;
1735 changed = MorphologyApply(old_image,new_image,curr_method,channel,
1736 curr_kernel, exception);
1738 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1739 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1740 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1743 old_image=DestroyImage(old_image);
1746 /* We are finished with kernel - destroy it if we made a clone */
1747 if ( curr_kernel != kernel )
1748 curr_kernel=DestroyKernelInfo(curr_kernel);
1750 /* Third-level Subtractive methods post-processing */
1752 case EdgeOutMorphology:
1753 case EdgeInMorphology:
1754 case TopHatMorphology:
1755 case BottomHatMorphology:
1756 /* Get Difference relative to the original image */
1757 (void) CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1760 case EdgeMorphology: /* subtract the Erode from a Dilate */
1761 (void) CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1763 grad_image=DestroyImage(grad_image);
1773 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1777 + R o t a t e K e r n e l I n f o %
1781 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1783 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
1784 % restricted to 90 degree angles, but this may be improved in the future.
1786 % The format of the RotateKernelInfo method is:
1788 % void RotateKernelInfo(KernelInfo *kernel, double angle)
1790 % A description of each parameter follows:
1792 % o kernel: the Morphology/Convolution kernel
1794 % o angle: angle to rotate in degrees
1796 % This function is only internel to this module, as it is not finalized,
1797 % especially with regard to non-orthogonal angles, and rotation of larger
1800 static void RotateKernelInfo(KernelInfo *kernel, double angle)
1802 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
1804 ** TODO: expand beyond simple 90 degree rotates, flips and flops
1807 /* Modulus the angle */
1808 angle = fmod(angle, 360.0);
1812 if ( 315.0 < angle || angle <= 45.0 )
1813 return; /* no change! - At least at this time */
1815 switch (kernel->type) {
1816 /* These built-in kernels are cylindrical kernels, rotating is useless */
1817 case GaussianKernel:
1818 case LaplacianKernel:
1822 case ChebyshevKernel:
1823 case ManhattenKernel:
1824 case EuclideanKernel:
1827 /* These may be rotatable at non-90 angles in the future */
1828 /* but simply rotating them in multiples of 90 degrees is useless */
1834 /* These only allows a +/-90 degree rotation (by transpose) */
1835 /* A 180 degree rotation is useless */
1837 case RectangleKernel:
1838 if ( 135.0 < angle && angle <= 225.0 )
1840 if ( 225.0 < angle && angle <= 315.0 )
1844 /* these are freely rotatable in 90 degree units */
1846 case UndefinedKernel:
1847 case UserDefinedKernel:
1850 if ( 135.0 < angle && angle <= 225.0 )
1852 /* For a 180 degree rotation - also know as a reflection */
1853 /* This is actually a very very common operation! */
1854 /* Basically all that is needed is a reversal of the kernel data! */
1861 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
1862 t=k[i], k[i]=k[j], k[j]=t;
1864 kernel->x = (long) kernel->width - kernel->x - 1;
1865 kernel->y = (long) kernel->height - kernel->y - 1;
1866 angle = fmod(angle+180.0, 360.0);
1868 if ( 45.0 < angle && angle <= 135.0 )
1869 { /* Do a transpose and a flop, of the image, which results in a 90
1870 * degree rotation using two mirror operations.
1872 * WARNING: this assumes the original image was a 1 dimentional image
1873 * but currently that is the only built-ins it is applied to.
1877 t = (long) kernel->width;
1878 kernel->width = kernel->height;
1879 kernel->height = (unsigned long) t;
1881 kernel->x = kernel->y;
1883 angle = fmod(450.0 - angle, 360.0);
1885 /* At this point angle should be between -45 (315) and +45 degrees
1886 * In the future some form of non-orthogonal angled rotates could be
1887 * performed here, posibily with a linear kernel restriction.
1891 Not currently in use!
1892 { /* Do a flop, this assumes kernel is horizontally symetrical.
1893 * Each row of the kernel needs to be reversed!
1902 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
1903 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
1904 t=k[x], k[x]=k[r], k[r]=t;
1906 kernel->x = kernel->width - kernel->x - 1;
1907 angle = fmod(angle+180.0, 360.0);
1914 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1918 % S c a l e K e r n e l I n f o %
1922 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1924 % ScaleKernelInfo() scales the kernel by the given amount, with or without
1925 % normalization of the sum of the kernel values.
1927 % By default (no flags given) the values within the kernel is scaled
1928 % according the given scaling factor.
1930 % If any 'normalize_flags' are given the kernel will be normalized and then
1931 % further scaled by the scaleing factor value given. A 'PercentValue' flag
1932 % will cause the given scaling factor to be divided by one hundred percent.
1934 % Kernel normalization ('normalize_flags' given) is designed to ensure that
1935 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
1936 % morphology methods will fall into -1.0 to +1.0 range. Note however that
1937 % for non-HDRI versions of IM this may cause images to have any negative
1938 % results clipped, unless some 'clip' any negative output from 'Convolve'
1939 % with the use of some kernels.
1941 % More specifically. Kernels which only contain positive values (such as a
1942 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
1943 % ensuring a 0.0 to +1.0 convolution output range for non-HDRI images.
1945 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
1946 % the kernel will be scaled by the absolute of the sum of kernel values, so
1947 % that it will generally fall within the +/- 1.0 range.
1949 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
1950 % will be scaled by just the sum of the postive values, so that its output
1951 % range will again fall into the +/- 1.0 range.
1953 % For special kernels designed for locating shapes using 'Correlate', (often
1954 % only containing +1 and -1 values, representing foreground/brackground
1955 % matching) a special normalization method is provided to scale the positive
1956 % values seperatally to those of the negative values, so the kernel will be
1957 % forced to become a zero-sum kernel better suited to such searches.
1959 % WARNING: Correct normalization of the kernal assumes that the '*_range'
1960 % attributes within the kernel structure have been correctly set during the
1963 % NOTE: The values used for 'normalize_flags' have been selected specifically
1964 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
1965 % means CorrelateNormalizeValue, and '%' means PercentValue. All other
1966 % GeometryFlags values are ignored.
1968 % The format of the ScaleKernelInfo method is:
1970 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
1971 % const MagickStatusType normalize_flags )
1973 % A description of each parameter follows:
1975 % o kernel: the Morphology/Convolution kernel
1978 % multiply all values (after normalization) by this factor if not
1979 % zero. If the kernel is normalized regardless of any flags.
1981 % o normalize_flags:
1982 % GeometryFlags defining normalization method to use.
1983 % specifically: NormalizeValue, CorrelateNormalizeValue,
1984 % and/or PercentValue
1986 % This function is internal to this module only at this time, but can be
1987 % exported to other modules if needed.
1989 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
1990 const double scaling_factor,const GeometryFlags normalize_flags)
2000 if ( (normalize_flags&NormalizeValue) != 0 ) {
2001 /* normalize kernel appropriately */
2002 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
2003 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
2005 pos_scale = kernel->positive_range; /* special zero-summing kernel */
2007 /* force kernel into being a normalized zero-summing kernel */
2008 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
2009 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
2010 ? kernel->positive_range : 1.0;
2011 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
2012 ? -kernel->negative_range : 1.0;
2015 neg_scale = pos_scale;
2017 /* finialize scaling_factor for positive and negative components */
2018 pos_scale = scaling_factor/pos_scale;
2019 neg_scale = scaling_factor/neg_scale;
2020 if ( (normalize_flags&PercentValue) != 0 ) {
2025 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2026 if ( ! IsNan(kernel->values[i]) )
2027 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
2029 /* convolution output range */
2030 kernel->positive_range *= pos_scale;
2031 kernel->negative_range *= neg_scale;
2032 /* maximum and minimum values in kernel */
2033 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
2034 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
2036 /* swap kernel settings if user scaling factor is negative */
2037 if ( scaling_factor < MagickEpsilon ) {
2039 t = kernel->positive_range;
2040 kernel->positive_range = kernel->negative_range;
2041 kernel->negative_range = t;
2042 t = kernel->maximum;
2043 kernel->maximum = kernel->minimum;
2044 kernel->minimum = 1;
2051 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2055 + S h o w K e r n e l I n f o %
2059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2061 % ShowKernelInfo() outputs the details of the given kernel defination to
2062 % standard error, generally due to a users 'showkernel' option request.
2064 % The format of the ShowKernel method is:
2066 % void ShowKernelInfo(KernelInfo *kernel)
2068 % A description of each parameter follows:
2070 % o kernel: the Morphology/Convolution kernel
2072 % This function is internal to this module only at this time. That may change
2075 MagickExport void ShowKernelInfo(KernelInfo *kernel)
2081 "Kernel \"%s\" of size %lux%lu%+ld%+ld with values from %.*lg to %.*lg\n",
2082 MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
2083 kernel->width, kernel->height,
2084 kernel->x, kernel->y,
2085 GetMagickPrecision(), kernel->minimum,
2086 GetMagickPrecision(), kernel->maximum);
2087 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2088 GetMagickPrecision(), kernel->negative_range,
2089 GetMagickPrecision(), kernel->positive_range,
2090 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2091 for (i=v=0; v < (long) kernel->height; v++) {
2092 fprintf(stderr,"%2ld:",v);
2093 for (u=0; u < (long) kernel->width; u++, i++)
2094 if ( IsNan(kernel->values[i]) )
2095 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2097 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2098 GetMagickPrecision(), kernel->values[i]);
2099 fprintf(stderr,"\n");
2104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2108 + Z e r o K e r n e l N a n s %
2112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2114 % ZeroKernelNans() replaces any special 'nan' value that may be present in
2115 % the kernel with a zero value. This is typically done when the kernel will
2116 % be used in special hardware (GPU) convolution processors, to simply
2119 % The format of the ZeroKernelNans method is:
2121 % voidZeroKernelNans (KernelInfo *kernel)
2123 % A description of each parameter follows:
2125 % o kernel: the Morphology/Convolution kernel
2127 % FUTURE: return the information in a string for API usage.
2129 MagickExport void ZeroKernelNans(KernelInfo *kernel)
2134 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2135 if ( IsNan(kernel->values[i]) )
2136 kernel->values[i] = 0.0;