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 if (kernel_string == (const char *) NULL)
207 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
208 if (kernel == (KernelInfo *)NULL)
210 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
211 kernel->type=UserDefinedKernel;
212 kernel->signature=MagickSignature;
215 SetGeometryInfo(&args);
217 /* does it start with an alpha - Return a builtin kernel */
218 GetMagickToken(kernel_string,&p,token);
219 if (isalpha((int) ((unsigned char) *token)) != 0)
224 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
225 if ( type < 0 || type == UserDefinedKernel )
226 return((KernelInfo *)NULL);
228 while (((isspace((int) ((unsigned char) *p)) != 0) ||
229 (*p == ',') || (*p == ':' )) && (*p != '\0'))
231 flags = ParseGeometry(p, &args);
233 /* special handling of missing values in input string */
235 case RectangleKernel:
236 if ( (flags & WidthValue) == 0 ) /* if no width then */
237 args.rho = args.sigma; /* then width = height */
238 if ( args.rho < 1.0 ) /* if width too small */
239 args.rho = 3; /* then width = 3 */
240 if ( args.sigma < 1.0 ) /* if height too small */
241 args.sigma = args.rho; /* then height = width */
242 if ( (flags & XValue) == 0 ) /* center offset if not defined */
243 args.xi = (double)(((long)args.rho-1)/2);
244 if ( (flags & YValue) == 0 )
245 args.psi = (double)(((long)args.sigma-1)/2);
251 if ( (flags & HeightValue) == 0 ) /* if no scale */
252 args.sigma = 1.0; /* then scale = 1.0 */
258 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
261 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
262 if (kernel == (KernelInfo *)NULL)
264 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
265 kernel->type = UserDefinedKernel;
266 kernel->signature = MagickSignature;
268 /* Has a ':' in argument - New user kernel specification */
269 p = strchr(kernel_string, ':');
270 if ( p != (char *) NULL)
272 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
273 memcpy(token, kernel_string, (size_t) (p-kernel_string));
274 token[p-kernel_string] = '\0';
275 flags = ParseGeometry(token, &args);
277 /* Size handling and checks of geometry settings */
278 if ( (flags & WidthValue) == 0 ) /* if no width then */
279 args.rho = args.sigma; /* then width = height */
280 if ( args.rho < 1.0 ) /* if width too small */
281 args.rho = 1.0; /* then width = 1 */
282 if ( args.sigma < 1.0 ) /* if height too small */
283 args.sigma = args.rho; /* then height = width */
284 kernel->width = (unsigned long)args.rho;
285 kernel->height = (unsigned long)args.sigma;
287 /* Offset Handling and Checks */
288 if ( args.xi < 0.0 || args.psi < 0.0 )
289 return(DestroyKernelInfo(kernel));
290 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
291 : (long) (kernel->width-1)/2;
292 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
293 : (long) (kernel->height-1)/2;
294 if ( kernel->x >= (long) kernel->width ||
295 kernel->y >= (long) kernel->height )
296 return(DestroyKernelInfo(kernel));
298 p++; /* advance beyond the ':' */
301 { /* ELSE - Old old kernel specification, forming odd-square kernel */
302 /* count up number of values given */
303 p=(const char *) kernel_string;
304 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
305 p++; /* ignore "'" chars for convolve filter usage - Cristy */
306 for (i=0; *p != '\0'; i++)
308 GetMagickToken(p,&p,token);
310 GetMagickToken(p,&p,token);
312 /* set the size of the kernel - old sized square */
313 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
314 kernel->x = kernel->y = (long) (kernel->width-1)/2;
315 p=(const char *) kernel_string;
316 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
317 p++; /* ignore "'" chars for convolve filter usage - Cristy */
320 /* Read in the kernel values from rest of input string argument */
321 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
322 kernel->height*sizeof(double));
323 if (kernel->values == (double *) NULL)
324 return(DestroyKernelInfo(kernel));
326 kernel->minimum = +MagickHuge;
327 kernel->maximum = -MagickHuge;
328 kernel->negative_range = kernel->positive_range = 0.0;
329 for (i=0; (i < (long) (kernel->width*kernel->height)) && (*p != '\0'); i++)
331 GetMagickToken(p,&p,token);
333 GetMagickToken(p,&p,token);
334 if ( LocaleCompare("nan",token) == 0
335 || LocaleCompare("-",token) == 0 ) {
336 kernel->values[i] = nan; /* do not include this value in kernel */
339 kernel->values[i] = StringToDouble(token);
340 ( kernel->values[i] < 0)
341 ? ( kernel->negative_range += kernel->values[i] )
342 : ( kernel->positive_range += kernel->values[i] );
343 Minimize(kernel->minimum, kernel->values[i]);
344 Maximize(kernel->maximum, kernel->values[i]);
347 /* check that we recieved at least one real (non-nan) value! */
348 if ( kernel->minimum == MagickHuge )
349 return(DestroyKernelInfo(kernel));
351 /* This should not be needed for a fully defined kernel
352 * Perhaps an error should be reported instead!
353 * Kept for backward compatibility.
355 if ( i < (long) (kernel->width*kernel->height) ) {
356 Minimize(kernel->minimum, kernel->values[i]);
357 Maximize(kernel->maximum, kernel->values[i]);
358 for ( ; i < (long) (kernel->width*kernel->height); i++)
359 kernel->values[i]=0.0;
366 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
370 % A c q u i r e K e r n e l B u i l t I n %
374 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
376 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
377 % kernels used for special purposes such as gaussian blurring, skeleton
378 % pruning, and edge distance determination.
380 % They take a KernelType, and a set of geometry style arguments, which were
381 % typically decoded from a user supplied string, or from a more complex
382 % Morphology Method that was requested.
384 % The format of the AcquireKernalBuiltIn method is:
386 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
387 % const GeometryInfo args)
389 % A description of each parameter follows:
391 % o type: the pre-defined type of kernel wanted
393 % o args: arguments defining or modifying the kernel
395 % Convolution Kernels
397 % Gaussian "{radius},{sigma}"
398 % Generate a two-dimentional gaussian kernel, as used by -gaussian
399 % A sigma is required, (with the 'x'), due to historical reasons.
401 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
402 % the final size of the resulting kernel to a square 2*radius+1 in size.
403 % The radius should be at least 2 times that of the sigma value, or
404 % sever clipping and aliasing may result. If not given or set to 0 the
405 % radius will be determined so as to produce the best minimal error
406 % result, which is usally much larger than is normally needed.
408 % Blur "{radius},{sigma},{angle}"
409 % As per Gaussian, but generates a 1 dimensional or linear gaussian
410 % blur, at the angle given (current restricted to orthogonal angles).
411 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
413 % NOTE that two such blurs perpendicular to each other is equivelent to
414 % -blur and the previous gaussian, but is often 10 or more times faster.
416 % Comet "{width},{sigma},{angle}"
417 % Blur in one direction only, mush like how a bright object leaves
418 % a comet like trail. The Kernel is actually half a gaussian curve,
419 % Adding two such blurs in oppiste directions produces a Linear Blur.
421 % NOTE: that the first argument is the width of the kernel and not the
422 % radius of the kernel.
424 % # Still to be implemented...
426 % # Sharpen "{radius},{sigma}
427 % # Negated Gaussian (center zeroed and re-normalized),
428 % # with a 2 unit positive peak. -- Check On line documentation
430 % # Laplacian "{radius},{sigma}"
431 % # Laplacian (a mexican hat like) Function
433 % # LOG "{radius},{sigma1},{sigma2}
434 % # Laplacian of Gaussian
436 % # DOG "{radius},{sigma1},{sigma2}
437 % # Difference of two Gaussians
441 % # Set kernel values using a resize filter, and given scale (sigma)
442 % # Cylindrical or Linear. Is this posible with an image?
447 % Rectangle "{geometry}"
448 % Simply generate a rectangle of 1's with the size given. You can also
449 % specify the location of the 'control point', otherwise the closest
450 % pixel to the center of the rectangle is selected.
452 % Properly centered and odd sized rectangles work the best.
454 % Diamond "[{radius}[,{scale}]]"
455 % Generate a diamond shaped kernal with given radius to the points.
456 % Kernel size will again be radius*2+1 square and defaults to radius 1,
457 % generating a 3x3 kernel that is slightly larger than a square.
459 % Square "[{radius}[,{scale}]]"
460 % Generate a square shaped kernel of size radius*2+1, and defaulting
461 % to a 3x3 (radius 1).
463 % Note that using a larger radius for the "Square" or the "Diamond"
464 % is also equivelent to iterating the basic morphological method
465 % that many times. However However iterating with the smaller radius 1
466 % default is actually faster than using a larger kernel radius.
468 % Disk "[{radius}[,{scale}]]
469 % Generate a binary disk of the radius given, radius may be a float.
470 % Kernel size will be ceil(radius)*2+1 square.
471 % NOTE: Here are some disk shapes of specific interest
472 % "disk:1" => "diamond" or "cross:1"
473 % "disk:1.5" => "square"
474 % "disk:2" => "diamond:2"
475 % "disk:2.5" => a general disk shape of radius 2
476 % "disk:2.9" => "square:2"
477 % "disk:3.5" => default - octagonal/disk shape of radius 3
478 % "disk:4.2" => roughly octagonal shape of radius 4
479 % "disk:4.3" => a general disk shape of radius 4
480 % After this all the kernel shape becomes more and more circular.
482 % Because a "disk" is more circular when using a larger radius, using a
483 % larger radius is preferred over iterating the morphological operation.
485 % Plus "[{radius}[,{scale}]]"
486 % Generate a kernel in the shape of a 'plus' sign. The length of each
487 % arm is also the radius, which defaults to 2.
489 % This kernel is not a good general morphological kernel, but is used
490 % more for highlighting and marking any single pixels in an image using,
491 % a "Dilate" or "Erode" method as appropriate.
493 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
495 % Note that unlike other kernels iterating a plus does not produce the
496 % same result as using a larger radius for the cross.
498 % Distance Measuring Kernels
500 % Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1)
501 % Manhatten "[{radius}][x{scale}]" square grid distance (default r=1)
502 % Euclidean "[{radius}][x{scale}]" direct distance (default r=1)
504 % Different types of distance measuring methods, which are used with the
505 % a 'Distance' morphology method for generating a gradient based on
506 % distance from an edge of a binary shape, though there is a technique
507 % for handling a anti-aliased shape.
509 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
510 % one to any neighbour, orthogonal or diagonal. One why of thinking of
511 % it is the number of squares a 'King' or 'Queen' in chess needs to
512 % traverse reach any other position on a chess board. It results in a
513 % 'square' like distance function, but one where diagonals are closer
516 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
517 % Cab metric), is the distance needed when you can only travel in
518 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
519 % in chess would travel. It results in a diamond like distances, where
520 % diagonals are further than expected.
522 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
523 % However by default the kernel size only has a radius of 1, which
524 % limits the distance to 'Knight' like moves, with only orthogonal and
525 % diagonal measurements being correct. As such for the default kernel
526 % you will get octagonal like distance function, which is reasonally
529 % However if you use a larger radius such as "Euclidean:4" you will
530 % get a much smoother distance gradient from the edge of the shape.
531 % Of course a larger kernel is slower to use, and generally not needed.
533 % To allow the use of fractional distances that you get with diagonals
534 % the actual distance is scaled by a fixed value which the user can
535 % provide. This is not actually nessary for either ""Chebyshev" or
536 % "Manhatten" distance kernels, but is done for all three distance
537 % kernels. If no scale is provided it is set to a value of 100,
538 % allowing for a maximum distance measurement of 655 pixels using a Q16
539 % version of IM, from any edge. However for small images this can
540 % result in quite a dark gradient.
542 % See the 'Distance' Morphological Method, for information of how it is
545 % # Hit-n-Miss Kernel-Lists -- Still to be implemented
547 % # specifically for Pruning, Thinning, Thickening
551 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
552 const GeometryInfo *args)
565 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
567 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
568 if (kernel == (KernelInfo *) NULL)
570 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
571 kernel->minimum = kernel->maximum = 0.0;
572 kernel->negative_range = kernel->positive_range = 0.0;
574 kernel->signature = MagickSignature;
577 /* Convolution Kernels */
580 sigma = fabs(args->sigma);
582 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
584 kernel->width = kernel->height =
585 GetOptimalKernelWidth2D(args->rho,sigma);
586 kernel->x = kernel->y = (long) (kernel->width-1)/2;
587 kernel->negative_range = kernel->positive_range = 0.0;
588 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
589 kernel->height*sizeof(double));
590 if (kernel->values == (double *) NULL)
591 return(DestroyKernelInfo(kernel));
593 sigma = 2.0*sigma*sigma; /* simplify the expression */
594 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
595 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
596 kernel->positive_range += (
598 exp(-((double)(u*u+v*v))/sigma)
599 /* / (MagickPI*sigma) */ );
601 kernel->maximum = kernel->values[
602 kernel->y*kernel->width+kernel->x ];
604 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
610 sigma = fabs(args->sigma);
612 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
614 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
615 kernel->x = (long) (kernel->width-1)/2;
618 kernel->negative_range = kernel->positive_range = 0.0;
619 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
620 kernel->height*sizeof(double));
621 if (kernel->values == (double *) NULL)
622 return(DestroyKernelInfo(kernel));
626 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
627 ** It generates a gaussian 3 times the width, and compresses it into
628 ** the expected range. This produces a closer normalization of the
629 ** resulting kernel, especially for very low sigma values.
630 ** As such while wierd it is prefered.
632 ** I am told this method originally came from Photoshop.
634 sigma *= KernelRank; /* simplify expanded curve */
635 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
636 (void) ResetMagickMemory(kernel->values,0, (size_t)
637 kernel->width*sizeof(double));
638 for ( u=-v; u <= v; u++) {
639 kernel->values[(u+v)/KernelRank] +=
640 exp(-((double)(u*u))/(2.0*sigma*sigma))
641 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
643 for (i=0; i < (long) kernel->width; i++)
644 kernel->positive_range += kernel->values[i];
646 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
647 kernel->positive_range += (
649 exp(-((double)(u*u))/(2.0*sigma*sigma))
650 /* / (MagickSQ2PI*sigma) */ );
653 kernel->maximum = kernel->values[ kernel->x ];
654 /* Note that neither methods above generate a normalized kernel,
655 ** though it gets close. The kernel may be 'clipped' by a user defined
656 ** radius, producing a smaller (darker) kernel. Also for very small
657 ** sigma's (> 0.1) the central value becomes larger than one, and thus
658 ** producing a very bright kernel.
661 /* Normalize the 1D Gaussian Kernel
663 ** Because of this the divisor in the above kernel generator is
664 ** not needed, so is not done above.
666 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
668 /* rotate the kernel by given angle */
669 RotateKernelInfo(kernel, args->xi);
674 sigma = fabs(args->sigma);
676 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
678 if ( args->rho < 1.0 )
679 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
681 kernel->width = (unsigned long)args->rho;
682 kernel->x = kernel->y = 0;
684 kernel->negative_range = kernel->positive_range = 0.0;
685 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
686 kernel->height*sizeof(double));
687 if (kernel->values == (double *) NULL)
688 return(DestroyKernelInfo(kernel));
690 /* A comet blur is half a gaussian curve, so that the object is
691 ** blurred in one direction only. This may not be quite the right
692 ** curve so may change in the future. The function must be normalised.
696 sigma *= KernelRank; /* simplify expanded curve */
697 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
698 (void) ResetMagickMemory(kernel->values,0, (size_t)
699 kernel->width*sizeof(double));
700 for ( u=0; u < v; u++) {
701 kernel->values[u/KernelRank] +=
702 exp(-((double)(u*u))/(2.0*sigma*sigma))
703 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
705 for (i=0; i < (long) kernel->width; i++)
706 kernel->positive_range += kernel->values[i];
708 for ( i=0; i < (long) kernel->width; i++)
709 kernel->positive_range += (
711 exp(-((double)(i*i))/(2.0*sigma*sigma))
712 /* / (MagickSQ2PI*sigma) */ );
715 kernel->maximum = kernel->values[0];
717 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
718 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
721 /* Boolean Kernels */
722 case RectangleKernel:
726 if ( type == SquareKernel )
729 kernel->width = kernel->height = 3; /* default radius = 1 */
731 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
732 kernel->x = kernel->y = (long) (kernel->width-1)/2;
736 /* NOTE: user defaults set in "AcquireKernelInfo()" */
737 if ( args->rho < 1.0 || args->sigma < 1.0 )
738 return(DestroyKernelInfo(kernel)); /* invalid args given */
739 kernel->width = (unsigned long)args->rho;
740 kernel->height = (unsigned long)args->sigma;
741 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
742 args->psi < 0.0 || args->psi > (double)kernel->height )
743 return(DestroyKernelInfo(kernel)); /* invalid args given */
744 kernel->x = (long) args->xi;
745 kernel->y = (long) args->psi;
748 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
749 kernel->height*sizeof(double));
750 if (kernel->values == (double *) NULL)
751 return(DestroyKernelInfo(kernel));
753 /* set all kernel values to 1.0 */
754 u=(long) kernel->width*kernel->height;
755 for ( i=0; i < u; i++)
756 kernel->values[i] = scale;
757 kernel->minimum = kernel->maximum = scale; /* a flat shape */
758 kernel->positive_range = scale*u;
764 kernel->width = kernel->height = 3; /* default radius = 1 */
766 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
767 kernel->x = kernel->y = (long) (kernel->width-1)/2;
769 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
770 kernel->height*sizeof(double));
771 if (kernel->values == (double *) NULL)
772 return(DestroyKernelInfo(kernel));
774 /* set all kernel values within diamond area to scale given */
775 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
776 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
777 if ((labs(u)+labs(v)) <= (long)kernel->x)
778 kernel->positive_range += kernel->values[i] = args->sigma;
780 kernel->values[i] = nan;
781 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
789 limit = (long)(args->rho*args->rho);
790 if (args->rho < 0.1) /* default radius approx 3.5 */
791 kernel->width = kernel->height = 7L, limit = 10L;
793 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
794 kernel->x = kernel->y = (long) (kernel->width-1)/2;
796 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
797 kernel->height*sizeof(double));
798 if (kernel->values == (double *) NULL)
799 return(DestroyKernelInfo(kernel));
801 /* set all kernel values within disk area to 1.0 */
802 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
803 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
804 if ((u*u+v*v) <= limit)
805 kernel->positive_range += kernel->values[i] = args->sigma;
807 kernel->values[i] = nan;
808 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
814 kernel->width = kernel->height = 5; /* default radius 2 */
816 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
817 kernel->x = kernel->y = (long) (kernel->width-1)/2;
819 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
820 kernel->height*sizeof(double));
821 if (kernel->values == (double *) NULL)
822 return(DestroyKernelInfo(kernel));
824 /* set all kernel values along axises to 1.0 */
825 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
826 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
827 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
828 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
829 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
832 /* Distance Measuring Kernels */
833 case ChebyshevKernel:
839 kernel->width = kernel->height = 3; /* default radius = 1 */
841 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
842 kernel->x = kernel->y = (long) (kernel->width-1)/2;
844 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
845 kernel->height*sizeof(double));
846 if (kernel->values == (double *) NULL)
847 return(DestroyKernelInfo(kernel));
849 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
850 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
851 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
852 kernel->positive_range += ( kernel->values[i] =
853 scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
854 kernel->maximum = kernel->values[0];
857 case ManhattenKernel:
863 kernel->width = kernel->height = 3; /* default radius = 1 */
865 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
866 kernel->x = kernel->y = (long) (kernel->width-1)/2;
868 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
869 kernel->height*sizeof(double));
870 if (kernel->values == (double *) NULL)
871 return(DestroyKernelInfo(kernel));
873 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
874 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
875 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
876 kernel->positive_range += ( kernel->values[i] =
877 scale*(labs(u)+labs(v)) );
878 kernel->maximum = kernel->values[0];
881 case EuclideanKernel:
887 kernel->width = kernel->height = 3; /* default radius = 1 */
889 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
890 kernel->x = kernel->y = (long) (kernel->width-1)/2;
892 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
893 kernel->height*sizeof(double));
894 if (kernel->values == (double *) NULL)
895 return(DestroyKernelInfo(kernel));
897 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
898 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
899 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
900 kernel->positive_range += ( kernel->values[i] =
901 scale*sqrt((double)(u*u+v*v)) );
902 kernel->maximum = kernel->values[0];
905 /* Undefined Kernels */
906 case LaplacianKernel:
909 perror("Kernel Type has not been defined yet");
912 /* Generate a No-Op minimal kernel - 1x1 pixel */
913 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
914 if (kernel->values == (double *) NULL)
915 return(DestroyKernelInfo(kernel));
916 kernel->width = kernel->height = 1;
917 kernel->x = kernel->x = 0;
918 kernel->type = UndefinedKernel;
920 kernel->positive_range =
921 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
929 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
933 % C l o n e K e r n e l I n f o %
937 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
939 % CloneKernelInfo() creates a new clone of the given Kernel so that its can
940 % be modified without effecting the original. The cloned kernel should be
941 % destroyed using DestoryKernelInfo() when no longer needed.
943 % The format of the CloneKernelInfo method is:
945 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
947 % A description of each parameter follows:
949 % o kernel: the Morphology/Convolution kernel to be cloned
952 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
960 assert(kernel != (KernelInfo *) NULL);
962 new=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
963 if (new == (KernelInfo *) NULL)
965 *new = *kernel; /* copy values in structure */
967 new->values=(double *) AcquireQuantumMemory(kernel->width,
968 kernel->height*sizeof(double));
969 if (new->values == (double *) NULL)
970 return(DestroyKernelInfo(new));
972 for (i=0; i < (long) (kernel->width*kernel->height); i++)
973 new->values[i] = kernel->values[i];
979 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
983 % D e s t r o y K e r n e l I n f o %
987 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
989 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
992 % The format of the DestroyKernelInfo method is:
994 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
996 % A description of each parameter follows:
998 % o kernel: the Morphology/Convolution kernel to be destroyed
1002 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1004 assert(kernel != (KernelInfo *) NULL);
1006 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1007 kernel->height*sizeof(double));
1008 kernel->values=(double *)RelinquishMagickMemory(kernel->values);
1009 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
1014 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1018 % M o r p h o l o g y I m a g e C h a n n e l %
1022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1024 % MorphologyImageChannel() applies a user supplied kernel to the image
1025 % according to the given mophology method.
1027 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1028 % by the kernel generator.
1030 % The format of the MorphologyImage method is:
1032 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
1033 % const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
1034 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1035 % channel,MorphologyMethod method,const long iterations,
1036 % KernelInfo *kernel,ExceptionInfo *exception)
1038 % A description of each parameter follows:
1040 % o image: the image.
1042 % o method: the morphology method to be applied.
1044 % o iterations: apply the operation this many times (or no change).
1045 % A value of -1 means loop until no change found.
1046 % How this is applied may depend on the morphology method.
1047 % Typically this is a value of 1.
1049 % o channel: the channel type.
1051 % o kernel: An array of double representing the morphology kernel.
1052 % Warning: kernel may be normalized for the Convolve method.
1054 % o exception: return any errors or warnings in this structure.
1057 % TODO: bias and auto-scale handling of the kernel for convolution
1058 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1059 % by the kernel generator.
1064 /* Internal function
1065 * Apply the Low-Level Morphology Method using the given Kernel
1066 * Returning the number of pixels that changed.
1067 * Two pre-created images must be provided, no image is created.
1069 static unsigned long MorphologyApply(const Image *image, Image
1070 *result_image, const MorphologyMethod method, const ChannelType channel,
1071 const KernelInfo *kernel, ExceptionInfo *exception)
1073 #define MorphologyTag "Morphology/Image"
1090 /* Only the most basic morphology is actually performed by this routine */
1093 Apply Basic Morphology to Image.
1099 GetMagickPixelPacket(image,&bias);
1100 SetMagickPixelPacketBias(image,&bias);
1101 /* Future: handle auto-bias from user, based on kernel input */
1103 p_view=AcquireCacheView(image);
1104 q_view=AcquireCacheView(result_image);
1106 /* Some methods (including convolve) needs use a reflected kernel.
1107 * Adjust 'origin' offsets for this reflected kernel.
1112 case ErodeMorphology:
1113 case ErodeIntensityMorphology:
1114 /* kernel is user as is, without reflection */
1116 case ConvolveMorphology:
1117 case DilateMorphology:
1118 case DilateIntensityMorphology:
1119 case DistanceMorphology:
1120 /* kernel needs to used with reflection */
1121 offx = (long) kernel->width-offx-1;
1122 offy = (long) kernel->height-offy-1;
1125 perror("Not a low level Morpholgy Method");
1129 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1130 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1132 for (y=0; y < (long) image->rows; y++)
1137 register const PixelPacket
1140 register const IndexPacket
1141 *restrict p_indexes;
1143 register PixelPacket
1146 register IndexPacket
1147 *restrict q_indexes;
1155 if (status == MagickFalse)
1157 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1158 image->columns+kernel->width, kernel->height, exception);
1159 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1161 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1166 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1167 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1168 r = (image->columns+kernel->width)*offy+offx; /* constant */
1170 for (x=0; x < (long) image->columns; x++)
1178 register const double
1181 register const PixelPacket
1184 register const IndexPacket
1185 *restrict k_indexes;
1190 /* Copy input to ouput image for unused channels
1191 * This removes need for 'cloning' a new image every iteration
1194 if (image->colorspace == CMYKColorspace)
1195 q_indexes[x] = p_indexes[r];
1197 result.green=(MagickRealType) 0;
1198 result.blue=(MagickRealType) 0;
1199 result.opacity=(MagickRealType) 0;
1200 result.index=(MagickRealType) 0;
1202 case ConvolveMorphology:
1203 /* Set the user defined bias of the weighted average output
1205 ** FUTURE: provide some way for internal functions to disable
1206 ** user defined bias and scaling effects.
1210 case DilateMorphology:
1215 result.index = -MagickHuge;
1217 case ErodeMorphology:
1222 result.index = +MagickHuge;
1224 case DilateIntensityMorphology:
1225 case ErodeIntensityMorphology:
1226 result.red = 0.0; /* flag indicating first match found */
1229 /* Otherwise just start with the original pixel value */
1230 result.red = (MagickRealType) p[r].red;
1231 result.green = (MagickRealType) p[r].green;
1232 result.blue = (MagickRealType) p[r].blue;
1233 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1234 if ( image->colorspace == CMYKColorspace)
1235 result.index = (MagickRealType) p_indexes[r];
1240 case ConvolveMorphology:
1241 /* Weighted Average of pixels using reflected kernel
1243 ** NOTE for correct working of this operation for asymetrical
1244 ** kernels, the kernel needs to be applied in its reflected form.
1245 ** That is its values needs to be reversed.
1247 ** Correlation is actually the same as this but without reflecting
1248 ** the kernel, and thus 'lower-level' that Convolution. However
1249 ** as Convolution is the more common method used, and it does not
1250 ** really cost us much in terms of processing to use a reflected
1251 ** kernel it is Convolution that is implemented.
1253 ** Correlation will have its kernel reflected before calling
1254 ** this function to do a Convolve.
1256 ** For more details of Correlation vs Convolution see
1257 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
1259 if (((channel & OpacityChannel) == 0) ||
1260 (image->matte == MagickFalse))
1262 /* Convolution without transparency effects */
1263 k = &kernel->values[ kernel->width*kernel->height-1 ];
1265 k_indexes = p_indexes;
1266 for (v=0; v < (long) kernel->height; v++) {
1267 for (u=0; u < (long) kernel->width; u++, k--) {
1268 if ( IsNan(*k) ) continue;
1269 result.red += (*k)*k_pixels[u].red;
1270 result.green += (*k)*k_pixels[u].green;
1271 result.blue += (*k)*k_pixels[u].blue;
1272 /* result.opacity += not involved here */
1273 if ( image->colorspace == CMYKColorspace)
1274 result.index += (*k)*k_indexes[u];
1276 k_pixels += image->columns+kernel->width;
1277 k_indexes += image->columns+kernel->width;
1281 { /* Kernel & Alpha weighted Convolution */
1283 alpha, /* alpha value * kernel weighting */
1284 gamma; /* weighting divisor */
1287 k = &kernel->values[ kernel->width*kernel->height-1 ];
1289 k_indexes = p_indexes;
1290 for (v=0; v < (long) kernel->height; v++) {
1291 for (u=0; u < (long) kernel->width; u++, k--) {
1292 if ( IsNan(*k) ) continue;
1293 alpha=(*k)*(QuantumScale*(QuantumRange-
1294 k_pixels[u].opacity));
1296 result.red += alpha*k_pixels[u].red;
1297 result.green += alpha*k_pixels[u].green;
1298 result.blue += alpha*k_pixels[u].blue;
1299 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1300 if ( image->colorspace == CMYKColorspace)
1301 result.index += alpha*k_indexes[u];
1303 k_pixels += image->columns+kernel->width;
1304 k_indexes += image->columns+kernel->width;
1306 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1307 result.red *= gamma;
1308 result.green *= gamma;
1309 result.blue *= gamma;
1310 result.opacity *= gamma;
1311 result.index *= gamma;
1315 case ErodeMorphology:
1316 /* Minimize Value within kernel neighbourhood
1318 ** NOTE that the kernel is not reflected for this operation!
1320 ** NOTE: in normal Greyscale Morphology, the kernel value should
1321 ** be added to the real value, this is currently not done, due to
1322 ** the nature of the boolean kernels being used.
1326 k_indexes = p_indexes;
1327 for (v=0; v < (long) kernel->height; v++) {
1328 for (u=0; u < (long) kernel->width; u++, k++) {
1329 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1330 Minimize(result.red, (double) k_pixels[u].red);
1331 Minimize(result.green, (double) k_pixels[u].green);
1332 Minimize(result.blue, (double) k_pixels[u].blue);
1333 Minimize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1334 if ( image->colorspace == CMYKColorspace)
1335 Minimize(result.index, (double) k_indexes[u]);
1337 k_pixels += image->columns+kernel->width;
1338 k_indexes += image->columns+kernel->width;
1342 case DilateMorphology:
1343 /* Maximize Value within kernel neighbourhood
1345 ** NOTE for correct working of this operation for asymetrical
1346 ** kernels, the kernel needs to be applied in its reflected form.
1347 ** That is its values needs to be reversed.
1349 ** NOTE: in normal Greyscale Morphology, the kernel value should
1350 ** be added to the real value, this is currently not done, due to
1351 ** the nature of the boolean kernels being used.
1354 k = &kernel->values[ kernel->width*kernel->height-1 ];
1356 k_indexes = p_indexes;
1357 for (v=0; v < (long) kernel->height; v++) {
1358 for (u=0; u < (long) kernel->width; u++, k--) {
1359 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1360 Maximize(result.red, (double) k_pixels[u].red);
1361 Maximize(result.green, (double) k_pixels[u].green);
1362 Maximize(result.blue, (double) k_pixels[u].blue);
1363 Maximize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1364 if ( image->colorspace == CMYKColorspace)
1365 Maximize(result.index, (double) k_indexes[u]);
1367 k_pixels += image->columns+kernel->width;
1368 k_indexes += image->columns+kernel->width;
1372 case ErodeIntensityMorphology:
1373 /* Select Pixel with Minimum Intensity within kernel neighbourhood
1375 ** WARNING: the intensity test fails for CMYK and does not
1376 ** take into account the moderating effect of teh alpha channel
1377 ** on the intensity.
1379 ** NOTE that the kernel is not reflected for this operation!
1383 k_indexes = p_indexes;
1384 for (v=0; v < (long) kernel->height; v++) {
1385 for (u=0; u < (long) kernel->width; u++, k++) {
1386 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1387 if ( result.red == 0.0 ||
1388 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1389 /* copy the whole pixel - no channel selection */
1391 if ( result.red > 0.0 ) changed++;
1395 k_pixels += image->columns+kernel->width;
1396 k_indexes += image->columns+kernel->width;
1400 case DilateIntensityMorphology:
1401 /* Select Pixel with Maximum Intensity within kernel neighbourhood
1403 ** WARNING: the intensity test fails for CMYK and does not
1404 ** take into account the moderating effect of teh alpha channel
1405 ** on the intensity.
1407 ** NOTE for correct working of this operation for asymetrical
1408 ** kernels, the kernel needs to be applied in its reflected form.
1409 ** That is its values needs to be reversed.
1411 k = &kernel->values[ kernel->width*kernel->height-1 ];
1413 k_indexes = p_indexes;
1414 for (v=0; v < (long) kernel->height; v++) {
1415 for (u=0; u < (long) kernel->width; u++, k--) {
1416 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1417 if ( result.red == 0.0 ||
1418 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1419 /* copy the whole pixel - no channel selection */
1421 if ( result.red > 0.0 ) changed++;
1425 k_pixels += image->columns+kernel->width;
1426 k_indexes += image->columns+kernel->width;
1430 case DistanceMorphology:
1431 /* Add kernel Value and select the minimum value found.
1432 ** The result is a iterative distance from edge of image shape.
1434 ** All Distance Kernels are symetrical, but that may not always
1435 ** be the case. For example how about a distance from left edges?
1436 ** To work correctly with asymetrical kernels the reflected kernel
1437 ** needs to be applied.
1440 /* No need to do distance morphology if original value is zero
1441 ** Unfortunatally I have not been able to get this right
1442 ** when channel selection also becomes involved. -- Arrgghhh
1444 if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1445 || ((channel & GreenChannel) == 0 && p[r].green == 0)
1446 || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1447 || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1448 || (( (channel & IndexChannel) == 0
1449 || image->colorspace != CMYKColorspace
1450 ) && p_indexes[x] ==0 )
1454 k = &kernel->values[ kernel->width*kernel->height-1 ];
1456 k_indexes = p_indexes;
1457 for (v=0; v < (long) kernel->height; v++) {
1458 for (u=0; u < (long) kernel->width; u++, k--) {
1459 if ( IsNan(*k) ) continue;
1460 Minimize(result.red, (*k)+k_pixels[u].red);
1461 Minimize(result.green, (*k)+k_pixels[u].green);
1462 Minimize(result.blue, (*k)+k_pixels[u].blue);
1463 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1464 if ( image->colorspace == CMYKColorspace)
1465 Minimize(result.index, (*k)+k_indexes[u]);
1467 k_pixels += image->columns+kernel->width;
1468 k_indexes += image->columns+kernel->width;
1472 case UndefinedMorphology:
1474 break; /* Do nothing */
1477 case UndefinedMorphology:
1478 case DilateIntensityMorphology:
1479 case ErodeIntensityMorphology:
1480 break; /* full pixel was directly assigned - not a channel method */
1482 /* Assign the results */
1483 if ((channel & RedChannel) != 0)
1484 q->red = ClampToQuantum(result.red);
1485 if ((channel & GreenChannel) != 0)
1486 q->green = ClampToQuantum(result.green);
1487 if ((channel & BlueChannel) != 0)
1488 q->blue = ClampToQuantum(result.blue);
1489 if ((channel & OpacityChannel) != 0
1490 && image->matte == MagickTrue )
1491 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1492 if ((channel & IndexChannel) != 0
1493 && image->colorspace == CMYKColorspace)
1494 q_indexes[x] = ClampToQuantum(result.index);
1497 if ( ( p[r].red != q->red )
1498 || ( p[r].green != q->green )
1499 || ( p[r].blue != q->blue )
1500 || ( p[r].opacity != q->opacity )
1501 || ( image->colorspace == CMYKColorspace &&
1502 p_indexes[r] != q_indexes[x] ) )
1503 changed++; /* The pixel had some value changed! */
1507 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1508 if (sync == MagickFalse)
1510 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1515 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1516 #pragma omp critical (MagickCore_MorphologyImage)
1518 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1519 if (proceed == MagickFalse)
1523 result_image->type=image->type;
1524 q_view=DestroyCacheView(q_view);
1525 p_view=DestroyCacheView(p_view);
1526 return(status ? (unsigned long) changed : 0);
1530 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1531 method, const long iterations,const KernelInfo *kernel, ExceptionInfo
1537 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1538 iterations,kernel,exception);
1539 return(morphology_image);
1543 MagickExport Image *MorphologyImageChannel(const Image *image,
1544 const ChannelType channel,const MorphologyMethod method,
1545 const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
1568 assert(image != (Image *) NULL);
1569 assert(image->signature == MagickSignature);
1570 assert(kernel != (KernelInfo *) NULL);
1571 assert(kernel->signature == MagickSignature);
1572 assert(exception != (ExceptionInfo *) NULL);
1573 assert(exception->signature == MagickSignature);
1575 if ( iterations == 0 )
1576 return((Image *)NULL); /* null operation - nothing to do! */
1578 /* kernel must be valid at this point
1579 * (except maybe for posible future morphology methods like "Prune"
1581 assert(kernel != (KernelInfo *)NULL);
1583 count = 0; /* interation count */
1584 changed = 1; /* if compound method assume image was changed */
1585 curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
1586 curr_method = method; /* to be changed as nessary */
1588 limit = (unsigned long) iterations;
1589 if ( iterations < 0 )
1590 limit = image->columns > image->rows ? image->columns : image->rows;
1592 /* Third-level morphology methods */
1593 grad_image=(Image *) NULL;
1594 switch( curr_method ) {
1595 case EdgeMorphology:
1596 grad_image = MorphologyImageChannel(image, channel,
1597 DilateMorphology, iterations, curr_kernel, exception);
1599 case EdgeInMorphology:
1600 curr_method = ErodeMorphology;
1602 case EdgeOutMorphology:
1603 curr_method = DilateMorphology;
1605 case TopHatMorphology:
1606 curr_method = OpenMorphology;
1608 case BottomHatMorphology:
1609 curr_method = CloseMorphology;
1612 break; /* not a third-level method */
1615 /* Second-level morphology methods */
1616 switch( curr_method ) {
1617 case OpenMorphology:
1618 /* Open is a Erode then a Dilate without reflection */
1619 new_image = MorphologyImageChannel(image, channel,
1620 ErodeMorphology, iterations, curr_kernel, exception);
1621 if (new_image == (Image *) NULL)
1622 return((Image *) NULL);
1623 curr_method = DilateMorphology;
1625 case OpenIntensityMorphology:
1626 new_image = MorphologyImageChannel(image, channel,
1627 ErodeIntensityMorphology, iterations, curr_kernel, exception);
1628 if (new_image == (Image *) NULL)
1629 return((Image *) NULL);
1630 curr_method = DilateIntensityMorphology;
1633 case CloseMorphology:
1634 /* Close is a Dilate then Erode using reflected kernel */
1635 /* A reflected kernel is needed for a Close */
1636 if ( curr_kernel == kernel )
1637 curr_kernel = CloneKernelInfo(kernel);
1638 RotateKernelInfo(curr_kernel,180);
1639 new_image = MorphologyImageChannel(image, channel,
1640 DilateMorphology, iterations, curr_kernel, exception);
1641 if (new_image == (Image *) NULL)
1642 return((Image *) NULL);
1643 curr_method = ErodeMorphology;
1645 case CloseIntensityMorphology:
1646 /* A reflected kernel is needed for a Close */
1647 if ( curr_kernel == kernel )
1648 curr_kernel = CloneKernelInfo(kernel);
1649 RotateKernelInfo(curr_kernel,180);
1650 new_image = MorphologyImageChannel(image, channel,
1651 DilateIntensityMorphology, iterations, curr_kernel, exception);
1652 if (new_image == (Image *) NULL)
1653 return((Image *) NULL);
1654 curr_method = ErodeIntensityMorphology;
1657 case CorrelateMorphology:
1658 /* A Correlation is actually a Convolution with a reflected kernel.
1659 ** However a Convolution is a weighted sum with a reflected kernel.
1660 ** It may seem stange to convert a Correlation into a Convolution
1661 ** as the Correleation is the simplier method, but Convolution is
1662 ** much more commonly used, and it makes sense to implement it directly
1663 ** so as to avoid the need to duplicate the kernel when it is not
1664 ** required (which is typically the default).
1666 if ( curr_kernel == kernel )
1667 curr_kernel = CloneKernelInfo(kernel);
1668 RotateKernelInfo(curr_kernel,180);
1669 curr_method = ConvolveMorphology;
1670 /* FALL-THRU into Correlate (weigthed sum without reflection) */
1672 case ConvolveMorphology:
1673 /* Scale or Normalize kernel, according to user wishes
1674 ** before using it for the Convolve/Correlate method.
1676 ** FUTURE: provide some way for internal functions to disable
1677 ** user bias and scaling effects.
1679 artifact = GetImageArtifact(image,"convolve:scale");
1680 if ( artifact != (char *)NULL ) {
1686 if ( curr_kernel == kernel )
1687 curr_kernel = CloneKernelInfo(kernel);
1690 flags = ParseGeometry(artifact, &args);
1691 ScaleKernelInfo(curr_kernel, args.rho, flags);
1693 /* FALL-THRU to do the first, and typically the only iteration */
1696 /* Do a single iteration using the Low-Level Morphology method!
1697 ** This ensures a "new_image" has been generated, but allows us to skip
1698 ** the creation of 'old_image' if no more iterations are needed.
1700 ** The "curr_method" should also be set to a low-level method that is
1701 ** understood by the MorphologyApply() internal function.
1703 new_image=CloneImage(image,0,0,MagickTrue,exception);
1704 if (new_image == (Image *) NULL)
1705 return((Image *) NULL);
1706 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1708 InheritException(exception,&new_image->exception);
1709 new_image=DestroyImage(new_image);
1710 return((Image *) NULL);
1712 changed = MorphologyApply(image,new_image,curr_method,channel,curr_kernel,
1715 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1716 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1717 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1722 /* At this point the "curr_method" should not only be set to a low-level
1723 ** method that is understood by the MorphologyApply() internal function,
1724 ** but "new_image" should now be defined, as the image to apply the
1725 ** "curr_method" to.
1728 /* Repeat the low-level morphology until count or no change reached */
1729 if ( count < (long) limit && changed > 0 ) {
1730 old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1731 if (old_image == (Image *) NULL)
1732 return(DestroyImage(new_image));
1733 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1735 InheritException(exception,&old_image->exception);
1736 old_image=DestroyImage(old_image);
1737 return(DestroyImage(new_image));
1739 while( count < (long) limit && changed != 0 )
1741 Image *tmp = old_image;
1742 old_image = new_image;
1744 changed = MorphologyApply(old_image,new_image,curr_method,channel,
1745 curr_kernel, exception);
1747 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1748 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1749 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1752 old_image=DestroyImage(old_image);
1755 /* We are finished with kernel - destroy it if we made a clone */
1756 if ( curr_kernel != kernel )
1757 curr_kernel=DestroyKernelInfo(curr_kernel);
1759 /* Third-level Subtractive methods post-processing */
1761 case EdgeOutMorphology:
1762 case EdgeInMorphology:
1763 case TopHatMorphology:
1764 case BottomHatMorphology:
1765 /* Get Difference relative to the original image */
1766 (void) CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1769 case EdgeMorphology: /* subtract the Erode from a Dilate */
1770 (void) CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1772 grad_image=DestroyImage(grad_image);
1782 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1786 + R o t a t e K e r n e l I n f o %
1790 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1792 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
1793 % restricted to 90 degree angles, but this may be improved in the future.
1795 % The format of the RotateKernelInfo method is:
1797 % void RotateKernelInfo(KernelInfo *kernel, double angle)
1799 % A description of each parameter follows:
1801 % o kernel: the Morphology/Convolution kernel
1803 % o angle: angle to rotate in degrees
1805 % This function is only internel to this module, as it is not finalized,
1806 % especially with regard to non-orthogonal angles, and rotation of larger
1809 static void RotateKernelInfo(KernelInfo *kernel, double angle)
1811 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
1813 ** TODO: expand beyond simple 90 degree rotates, flips and flops
1816 /* Modulus the angle */
1817 angle = fmod(angle, 360.0);
1821 if ( 315.0 < angle || angle <= 45.0 )
1822 return; /* no change! - At least at this time */
1824 switch (kernel->type) {
1825 /* These built-in kernels are cylindrical kernels, rotating is useless */
1826 case GaussianKernel:
1827 case LaplacianKernel:
1831 case ChebyshevKernel:
1832 case ManhattenKernel:
1833 case EuclideanKernel:
1836 /* These may be rotatable at non-90 angles in the future */
1837 /* but simply rotating them in multiples of 90 degrees is useless */
1843 /* These only allows a +/-90 degree rotation (by transpose) */
1844 /* A 180 degree rotation is useless */
1846 case RectangleKernel:
1847 if ( 135.0 < angle && angle <= 225.0 )
1849 if ( 225.0 < angle && angle <= 315.0 )
1853 /* these are freely rotatable in 90 degree units */
1855 case UndefinedKernel:
1856 case UserDefinedKernel:
1859 if ( 135.0 < angle && angle <= 225.0 )
1861 /* For a 180 degree rotation - also know as a reflection */
1862 /* This is actually a very very common operation! */
1863 /* Basically all that is needed is a reversal of the kernel data! */
1870 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
1871 t=k[i], k[i]=k[j], k[j]=t;
1873 kernel->x = (long) kernel->width - kernel->x - 1;
1874 kernel->y = (long) kernel->height - kernel->y - 1;
1875 angle = fmod(angle+180.0, 360.0);
1877 if ( 45.0 < angle && angle <= 135.0 )
1878 { /* Do a transpose and a flop, of the image, which results in a 90
1879 * degree rotation using two mirror operations.
1881 * WARNING: this assumes the original image was a 1 dimentional image
1882 * but currently that is the only built-ins it is applied to.
1886 t = (long) kernel->width;
1887 kernel->width = kernel->height;
1888 kernel->height = (unsigned long) t;
1890 kernel->x = kernel->y;
1892 angle = fmod(450.0 - angle, 360.0);
1894 /* At this point angle should be between -45 (315) and +45 degrees
1895 * In the future some form of non-orthogonal angled rotates could be
1896 * performed here, posibily with a linear kernel restriction.
1900 Not currently in use!
1901 { /* Do a flop, this assumes kernel is horizontally symetrical.
1902 * Each row of the kernel needs to be reversed!
1911 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
1912 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
1913 t=k[x], k[x]=k[r], k[r]=t;
1915 kernel->x = kernel->width - kernel->x - 1;
1916 angle = fmod(angle+180.0, 360.0);
1923 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1927 % S c a l e K e r n e l I n f o %
1931 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1933 % ScaleKernelInfo() scales the kernel by the given amount, with or without
1934 % normalization of the sum of the kernel values.
1936 % By default (no flags given) the values within the kernel is scaled
1937 % according the given scaling factor.
1939 % If any 'normalize_flags' are given the kernel will be normalized and then
1940 % further scaled by the scaleing factor value given. A 'PercentValue' flag
1941 % will cause the given scaling factor to be divided by one hundred percent.
1943 % Kernel normalization ('normalize_flags' given) is designed to ensure that
1944 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
1945 % morphology methods will fall into -1.0 to +1.0 range. Note however that
1946 % for non-HDRI versions of IM this may cause images to have any negative
1947 % results clipped, unless some 'clip' any negative output from 'Convolve'
1948 % with the use of some kernels.
1950 % More specifically. Kernels which only contain positive values (such as a
1951 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
1952 % ensuring a 0.0 to +1.0 convolution output range for non-HDRI images.
1954 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
1955 % the kernel will be scaled by the absolute of the sum of kernel values, so
1956 % that it will generally fall within the +/- 1.0 range.
1958 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
1959 % will be scaled by just the sum of the postive values, so that its output
1960 % range will again fall into the +/- 1.0 range.
1962 % For special kernels designed for locating shapes using 'Correlate', (often
1963 % only containing +1 and -1 values, representing foreground/brackground
1964 % matching) a special normalization method is provided to scale the positive
1965 % values seperatally to those of the negative values, so the kernel will be
1966 % forced to become a zero-sum kernel better suited to such searches.
1968 % WARNING: Correct normalization of the kernal assumes that the '*_range'
1969 % attributes within the kernel structure have been correctly set during the
1972 % NOTE: The values used for 'normalize_flags' have been selected specifically
1973 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
1974 % means CorrelateNormalizeValue, and '%' means PercentValue. All other
1975 % GeometryFlags values are ignored.
1977 % The format of the ScaleKernelInfo method is:
1979 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
1980 % const MagickStatusType normalize_flags )
1982 % A description of each parameter follows:
1984 % o kernel: the Morphology/Convolution kernel
1987 % multiply all values (after normalization) by this factor if not
1988 % zero. If the kernel is normalized regardless of any flags.
1990 % o normalize_flags:
1991 % GeometryFlags defining normalization method to use.
1992 % specifically: NormalizeValue, CorrelateNormalizeValue,
1993 % and/or PercentValue
1995 % This function is internal to this module only at this time, but can be
1996 % exported to other modules if needed.
1998 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
1999 const double scaling_factor,const GeometryFlags normalize_flags)
2009 if ( (normalize_flags&NormalizeValue) != 0 ) {
2010 /* normalize kernel appropriately */
2011 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
2012 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
2014 pos_scale = kernel->positive_range; /* special zero-summing kernel */
2016 /* force kernel into being a normalized zero-summing kernel */
2017 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
2018 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
2019 ? kernel->positive_range : 1.0;
2020 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
2021 ? -kernel->negative_range : 1.0;
2024 neg_scale = pos_scale;
2026 /* finialize scaling_factor for positive and negative components */
2027 pos_scale = scaling_factor/pos_scale;
2028 neg_scale = scaling_factor/neg_scale;
2029 if ( (normalize_flags&PercentValue) != 0 ) {
2034 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2035 if ( ! IsNan(kernel->values[i]) )
2036 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
2038 /* convolution output range */
2039 kernel->positive_range *= pos_scale;
2040 kernel->negative_range *= neg_scale;
2041 /* maximum and minimum values in kernel */
2042 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
2043 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
2045 /* swap kernel settings if user scaling factor is negative */
2046 if ( scaling_factor < MagickEpsilon ) {
2048 t = kernel->positive_range;
2049 kernel->positive_range = kernel->negative_range;
2050 kernel->negative_range = t;
2051 t = kernel->maximum;
2052 kernel->maximum = kernel->minimum;
2053 kernel->minimum = 1;
2060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2064 + S h o w K e r n e l I n f o %
2068 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2070 % ShowKernelInfo() outputs the details of the given kernel defination to
2071 % standard error, generally due to a users 'showkernel' option request.
2073 % The format of the ShowKernel method is:
2075 % void ShowKernelInfo(KernelInfo *kernel)
2077 % A description of each parameter follows:
2079 % o kernel: the Morphology/Convolution kernel
2081 % This function is internal to this module only at this time. That may change
2084 MagickExport void ShowKernelInfo(KernelInfo *kernel)
2090 "Kernel \"%s\" of size %lux%lu%+ld%+ld with values from %.*lg to %.*lg\n",
2091 MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
2092 kernel->width, kernel->height,
2093 kernel->x, kernel->y,
2094 GetMagickPrecision(), kernel->minimum,
2095 GetMagickPrecision(), kernel->maximum);
2096 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2097 GetMagickPrecision(), kernel->negative_range,
2098 GetMagickPrecision(), kernel->positive_range,
2099 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2100 for (i=v=0; v < (long) kernel->height; v++) {
2101 fprintf(stderr,"%2ld:",v);
2102 for (u=0; u < (long) kernel->width; u++, i++)
2103 if ( IsNan(kernel->values[i]) )
2104 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2106 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2107 GetMagickPrecision(), kernel->values[i]);
2108 fprintf(stderr,"\n");
2113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2117 + Z e r o K e r n e l N a n s %
2121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2123 % ZeroKernelNans() replaces any special 'nan' value that may be present in
2124 % the kernel with a zero value. This is typically done when the kernel will
2125 % be used in special hardware (GPU) convolution processors, to simply
2128 % The format of the ZeroKernelNans method is:
2130 % voidZeroKernelNans (KernelInfo *kernel)
2132 % A description of each parameter follows:
2134 % o kernel: the Morphology/Convolution kernel
2136 % FUTURE: return the information in a string for API usage.
2138 MagickExport void ZeroKernelNans(KernelInfo *kernel)
2143 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2144 if ( IsNan(kernel->values[i]) )
2145 kernel->values[i] = 0.0;