2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 % M M OOO RRRR PPPP H H OOO L OOO GGGG Y Y %
7 % MM MM O O R R P P H H O O L O O G Y Y %
8 % M M M O O RRRR PPPP HHHHH O O L O O G GGG Y %
9 % M M O O R R P H H O O L O O G G Y %
10 % M M OOO R R P H H OOO LLLLL OOO GGG Y %
13 % MagickCore Morphology Methods %
20 % Copyright 1999-2010 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % http://www.imagemagick.org/script/license.php %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 % Morpology is the the application of various kernals, of any size and even
37 % shape, to a image in various ways (typically binary, but not always).
39 % Convolution (weighted sum or average) is just one specific type of
40 % morphology. Just one that is very common for image bluring and sharpening
41 % effects. Not only 2D Gaussian blurring, but also 2-pass 1D Blurring.
43 % This module provides not only a general morphology function, and the ability
44 % to apply more advanced or iterative morphologies, but also functions for the
45 % generation of many different types of kernel arrays from user supplied
46 % arguments. Prehaps even the generation of a kernel from a small image.
52 #include "magick/studio.h"
53 #include "magick/artifact.h"
54 #include "magick/cache-view.h"
55 #include "magick/color-private.h"
56 #include "magick/enhance.h"
57 #include "magick/exception.h"
58 #include "magick/exception-private.h"
59 #include "magick/gem.h"
60 #include "magick/hashmap.h"
61 #include "magick/image.h"
62 #include "magick/image-private.h"
63 #include "magick/list.h"
64 #include "magick/magick.h"
65 #include "magick/memory_.h"
66 #include "magick/monitor-private.h"
67 #include "magick/morphology.h"
68 #include "magick/option.h"
69 #include "magick/pixel-private.h"
70 #include "magick/prepress.h"
71 #include "magick/quantize.h"
72 #include "magick/registry.h"
73 #include "magick/semaphore.h"
74 #include "magick/splay-tree.h"
75 #include "magick/statistic.h"
76 #include "magick/string_.h"
77 #include "magick/string-private.h"
78 #include "magick/token.h"
82 * The following test is for special floating point numbers of value NaN (not
83 * a number), that may be used within a Kernel Definition. NaN's are defined
84 * as part of the IEEE standard for floating point number representation.
86 * These are used a Kernel value of NaN means that that kernal position is not
87 * part of the normal convolution or morphology process, and thus allowing the
88 * use of 'shaped' kernels.
90 * Special Properities Two NaN's are never equal, even if they are from the
91 * same variable That is the IsNaN() macro is only true if the value is NaN.
93 #define IsNan(a) ((a)!=(a))
96 * Other global definitions used by module
98 static inline double MagickMin(const double x,const double y)
100 return( x < y ? x : y);
102 static inline double MagickMax(const double x,const double y)
104 return( x > y ? x : y);
106 #define Minimize(assign,value) assign=MagickMin(assign,value)
107 #define Maximize(assign,value) assign=MagickMax(assign,value)
109 /* Currently these are only internal to this module */
111 RotateKernelInfo(KernelInfo *, double),
112 ScaleKernelInfo(KernelInfo *, const double, const MagickStatusType);
115 *CloneKernelInfo(const KernelInfo *);
119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123 % A c q u i r e K e r n e l I n f o %
127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
129 % AcquireKernelInfo() takes the given string (generally supplied by the
130 % user) and converts it into a Morphology/Convolution Kernel. This allows
131 % users to specify a kernel from a number of pre-defined kernels, or to fully
132 % specify their own kernel for a specific Convolution or Morphology
135 % The kernel so generated can be any rectangular array of floating point
136 % values (doubles) with the 'control point' or 'pixel being affected'
137 % anywhere within that array of values.
139 % Previously IM was restricted to a square of odd size using the exact
140 % center as origin, this is no longer the case, and any rectangular kernel
141 % with any value being declared the origin. This in turn allows the use of
142 % highly asymmetrical kernels.
144 % The floating point values in the kernel can also include a special value
145 % known as 'nan' or 'not a number' to indicate that this value is not part
146 % of the kernel array. This allows you to shaped the kernel within its
147 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
148 % shape. However at least one non-nan value must be provided for correct
149 % working of a kernel.
151 % The returned kernel should be free using the DestroyKernelInfo() when you
152 % are finished with it.
154 % Input kernel defintion strings can consist of any of three types.
157 % Select from one of the built in kernels, using the name and
158 % geometry arguments supplied. See AcquireKernelBuiltIn()
160 % "WxH[+X+Y]:num, num, num ..."
161 % a kernal of size W by H, with W*H floating point numbers following.
162 % the 'center' can be optionally be defined at +X+Y (such that +0+0
163 % is top left corner). If not defined the pixel in the center, for
164 % odd sizes, or to the immediate top or left of center for even sizes
165 % is automatically selected.
167 % "num, num, num, num, ..."
168 % list of floating point numbers defining an 'old style' odd sized
169 % square kernel. At least 9 values should be provided for a 3x3
170 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
171 % Values can be space or comma separated. This is not recommended.
173 % Note that 'name' kernels will start with an alphabetic character while the
174 % new kernel specification has a ':' character in its specification string.
175 % If neither is the case, it is assumed an old style of a simple list of
176 % numbers generating a odd-sized square kernel has been given.
178 % The format of the AcquireKernal method is:
180 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
182 % A description of each parameter follows:
184 % o kernel_string: the Morphology/Convolution kernel wanted.
188 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
194 token[MaxTextExtent];
209 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
211 assert(kernel_string != (const char *) NULL);
212 SetGeometryInfo(&args);
214 /* does it start with an alpha - Return a builtin kernel */
215 GetMagickToken(kernel_string,&p,token);
216 if ( isalpha((int)token[0]) )
221 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
222 if ( type < 0 || type == UserDefinedKernel )
223 return((KernelInfo *)NULL);
225 while (((isspace((int) ((unsigned char) *p)) != 0) ||
226 (*p == ',') || (*p == ':' )) && (*p != '\0'))
228 flags = ParseGeometry(p, &args);
230 /* special handling of missing values in input string */
232 case RectangleKernel:
233 if ( (flags & WidthValue) == 0 ) /* if no width then */
234 args.rho = args.sigma; /* then width = height */
235 if ( args.rho < 1.0 ) /* if width too small */
236 args.rho = 3; /* then width = 3 */
237 if ( args.sigma < 1.0 ) /* if height too small */
238 args.sigma = args.rho; /* then height = width */
239 if ( (flags & XValue) == 0 ) /* center offset if not defined */
240 args.xi = (double)(((long)args.rho-1)/2);
241 if ( (flags & YValue) == 0 )
242 args.psi = (double)(((long)args.sigma-1)/2);
248 if ( (flags & HeightValue) == 0 ) /* if no scale */
249 args.sigma = 1.0; /* then scale = 1.0 */
255 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
258 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
259 if (kernel == (KernelInfo *)NULL)
261 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
262 kernel->type = UserDefinedKernel;
263 kernel->signature = MagickSignature;
265 /* Has a ':' in argument - New user kernel specification */
266 p = strchr(kernel_string, ':');
267 if ( p != (char *) NULL)
269 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
270 memcpy(token, kernel_string, (size_t) (p-kernel_string));
271 token[p-kernel_string] = '\0';
272 flags = ParseGeometry(token, &args);
274 /* Size handling and checks of geometry settings */
275 if ( (flags & WidthValue) == 0 ) /* if no width then */
276 args.rho = args.sigma; /* then width = height */
277 if ( args.rho < 1.0 ) /* if width too small */
278 args.rho = 1.0; /* then width = 1 */
279 if ( args.sigma < 1.0 ) /* if height too small */
280 args.sigma = args.rho; /* then height = width */
281 kernel->width = (unsigned long)args.rho;
282 kernel->height = (unsigned long)args.sigma;
284 /* Offset Handling and Checks */
285 if ( args.xi < 0.0 || args.psi < 0.0 )
286 return(DestroyKernelInfo(kernel));
287 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
288 : (long) (kernel->width-1)/2;
289 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
290 : (long) (kernel->height-1)/2;
291 if ( kernel->x >= (long) kernel->width ||
292 kernel->y >= (long) kernel->height )
293 return(DestroyKernelInfo(kernel));
295 p++; /* advance beyond the ':' */
298 { /* ELSE - Old old kernel specification, forming odd-square kernel */
299 /* count up number of values given */
300 p=(const char *) kernel_string;
301 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
302 p++; /* ignore "'" chars for convolve filter usage - Cristy */
303 for (i=0; *p != '\0'; i++)
305 GetMagickToken(p,&p,token);
307 GetMagickToken(p,&p,token);
309 /* set the size of the kernel - old sized square */
310 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
311 kernel->x = kernel->y = (long) (kernel->width-1)/2;
312 p=(const char *) kernel_string;
313 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
314 p++; /* ignore "'" chars for convolve filter usage - Cristy */
317 /* Read in the kernel values from rest of input string argument */
318 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
319 kernel->height*sizeof(double));
320 if (kernel->values == (double *) NULL)
321 return(DestroyKernelInfo(kernel));
323 kernel->minimum = +MagickHuge;
324 kernel->maximum = -MagickHuge;
325 kernel->negative_range = kernel->positive_range = 0.0;
326 for (i=0; (i < (long) (kernel->width*kernel->height)) && (*p != '\0'); i++)
328 GetMagickToken(p,&p,token);
330 GetMagickToken(p,&p,token);
331 if ( LocaleCompare("nan",token) == 0
332 || LocaleCompare("-",token) == 0 ) {
333 kernel->values[i] = nan; /* do not include this value in kernel */
336 kernel->values[i] = StringToDouble(token);
337 ( kernel->values[i] < 0)
338 ? ( kernel->negative_range += kernel->values[i] )
339 : ( kernel->positive_range += kernel->values[i] );
340 Minimize(kernel->minimum, kernel->values[i]);
341 Maximize(kernel->maximum, kernel->values[i]);
344 /* check that we recieved at least one real (non-nan) value! */
345 if ( kernel->minimum == MagickHuge )
346 return(DestroyKernelInfo(kernel));
348 /* This should not be needed for a fully defined kernel
349 * Perhaps an error should be reported instead!
350 * Kept for backward compatibility.
352 if ( i < (long) (kernel->width*kernel->height) ) {
353 Minimize(kernel->minimum, kernel->values[i]);
354 Maximize(kernel->maximum, kernel->values[i]);
355 for ( ; i < (long) (kernel->width*kernel->height); i++)
356 kernel->values[i]=0.0;
363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
367 % A c q u i r e K e r n e l B u i l t I n %
371 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
373 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
374 % kernels used for special purposes such as gaussian blurring, skeleton
375 % pruning, and edge distance determination.
377 % They take a KernelType, and a set of geometry style arguments, which were
378 % typically decoded from a user supplied string, or from a more complex
379 % Morphology Method that was requested.
381 % The format of the AcquireKernalBuiltIn method is:
383 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
384 % const GeometryInfo args)
386 % A description of each parameter follows:
388 % o type: the pre-defined type of kernel wanted
390 % o args: arguments defining or modifying the kernel
392 % Convolution Kernels
394 % Gaussian "{radius},{sigma}"
395 % Generate a two-dimentional gaussian kernel, as used by -gaussian
396 % A sigma is required, (with the 'x'), due to historical reasons.
398 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
399 % the final size of the resulting kernel to a square 2*radius+1 in size.
400 % The radius should be at least 2 times that of the sigma value, or
401 % sever clipping and aliasing may result. If not given or set to 0 the
402 % radius will be determined so as to produce the best minimal error
403 % result, which is usally much larger than is normally needed.
405 % Blur "{radius},{sigma},{angle}"
406 % As per Gaussian, but generates a 1 dimensional or linear gaussian
407 % blur, at the angle given (current restricted to orthogonal angles).
408 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
410 % NOTE that two such blurs perpendicular to each other is equivelent to
411 % -blur and the previous gaussian, but is often 10 or more times faster.
413 % Comet "{width},{sigma},{angle}"
414 % Blur in one direction only, mush like how a bright object leaves
415 % a comet like trail. The Kernel is actually half a gaussian curve,
416 % Adding two such blurs in oppiste directions produces a Linear Blur.
418 % NOTE: that the first argument is the width of the kernel and not the
419 % radius of the kernel.
421 % # Still to be implemented...
423 % # Sharpen "{radius},{sigma}
424 % # Negated Gaussian (center zeroed and re-normalized),
425 % # with a 2 unit positive peak. -- Check On line documentation
427 % # Laplacian "{radius},{sigma}"
428 % # Laplacian (a mexican hat like) Function
430 % # LOG "{radius},{sigma1},{sigma2}
431 % # Laplacian of Gaussian
433 % # DOG "{radius},{sigma1},{sigma2}
434 % # Difference of two Gaussians
438 % # Set kernel values using a resize filter, and given scale (sigma)
439 % # Cylindrical or Linear. Is this posible with an image?
444 % Rectangle "{geometry}"
445 % Simply generate a rectangle of 1's with the size given. You can also
446 % specify the location of the 'control point', otherwise the closest
447 % pixel to the center of the rectangle is selected.
449 % Properly centered and odd sized rectangles work the best.
451 % Diamond "[{radius}[,{scale}]]"
452 % Generate a diamond shaped kernal with given radius to the points.
453 % Kernel size will again be radius*2+1 square and defaults to radius 1,
454 % generating a 3x3 kernel that is slightly larger than a square.
456 % Square "[{radius}[,{scale}]]"
457 % Generate a square shaped kernel of size radius*2+1, and defaulting
458 % to a 3x3 (radius 1).
460 % Note that using a larger radius for the "Square" or the "Diamond"
461 % is also equivelent to iterating the basic morphological method
462 % that many times. However However iterating with the smaller radius 1
463 % default is actually faster than using a larger kernel radius.
465 % Disk "[{radius}[,{scale}]]
466 % Generate a binary disk of the radius given, radius may be a float.
467 % Kernel size will be ceil(radius)*2+1 square.
468 % NOTE: Here are some disk shapes of specific interest
469 % "disk:1" => "diamond" or "cross:1"
470 % "disk:1.5" => "square"
471 % "disk:2" => "diamond:2"
472 % "disk:2.5" => a general disk shape of radius 2
473 % "disk:2.9" => "square:2"
474 % "disk:3.5" => default - octagonal/disk shape of radius 3
475 % "disk:4.2" => roughly octagonal shape of radius 4
476 % "disk:4.3" => a general disk shape of radius 4
477 % After this all the kernel shape becomes more and more circular.
479 % Because a "disk" is more circular when using a larger radius, using a
480 % larger radius is preferred over iterating the morphological operation.
482 % Plus "[{radius}[,{scale}]]"
483 % Generate a kernel in the shape of a 'plus' sign. The length of each
484 % arm is also the radius, which defaults to 2.
486 % This kernel is not a good general morphological kernel, but is used
487 % more for highlighting and marking any single pixels in an image using,
488 % a "Dilate" or "Erode" method as appropriate.
490 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
492 % Note that unlike other kernels iterating a plus does not produce the
493 % same result as using a larger radius for the cross.
495 % Distance Measuring Kernels
497 % Chebyshev "[{radius}][x{scale}]" largest x or y distance (default r=1)
498 % Manhatten "[{radius}][x{scale}]" square grid distance (default r=1)
499 % Euclidean "[{radius}][x{scale}]" direct distance (default r=1)
501 % Different types of distance measuring methods, which are used with the
502 % a 'Distance' morphology method for generating a gradient based on
503 % distance from an edge of a binary shape, though there is a technique
504 % for handling a anti-aliased shape.
506 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
507 % one to any neighbour, orthogonal or diagonal. One why of thinking of
508 % it is the number of squares a 'King' or 'Queen' in chess needs to
509 % traverse reach any other position on a chess board. It results in a
510 % 'square' like distance function, but one where diagonals are closer
513 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
514 % Cab metric), is the distance needed when you can only travel in
515 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
516 % in chess would travel. It results in a diamond like distances, where
517 % diagonals are further than expected.
519 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
520 % However by default the kernel size only has a radius of 1, which
521 % limits the distance to 'Knight' like moves, with only orthogonal and
522 % diagonal measurements being correct. As such for the default kernel
523 % you will get octagonal like distance function, which is reasonally
526 % However if you use a larger radius such as "Euclidean:4" you will
527 % get a much smoother distance gradient from the edge of the shape.
528 % Of course a larger kernel is slower to use, and generally not needed.
530 % To allow the use of fractional distances that you get with diagonals
531 % the actual distance is scaled by a fixed value which the user can
532 % provide. This is not actually nessary for either ""Chebyshev" or
533 % "Manhatten" distance kernels, but is done for all three distance
534 % kernels. If no scale is provided it is set to a value of 100,
535 % allowing for a maximum distance measurement of 655 pixels using a Q16
536 % version of IM, from any edge. However for small images this can
537 % result in quite a dark gradient.
539 % See the 'Distance' Morphological Method, for information of how it is
542 % # Hit-n-Miss Kernel-Lists -- Still to be implemented
544 % # specifically for Pruning, Thinning, Thickening
548 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
549 const GeometryInfo *args)
562 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
564 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
565 if (kernel == (KernelInfo *) NULL)
567 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
568 kernel->minimum = kernel->maximum = 0.0;
569 kernel->negative_range = kernel->positive_range = 0.0;
571 kernel->signature = MagickSignature;
574 /* Convolution Kernels */
577 sigma = fabs(args->sigma);
579 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
581 kernel->width = kernel->height =
582 GetOptimalKernelWidth2D(args->rho,sigma);
583 kernel->x = kernel->y = (long) (kernel->width-1)/2;
584 kernel->negative_range = kernel->positive_range = 0.0;
585 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
586 kernel->height*sizeof(double));
587 if (kernel->values == (double *) NULL)
588 return(DestroyKernelInfo(kernel));
590 sigma = 2.0*sigma*sigma; /* simplify the expression */
591 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
592 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
593 kernel->positive_range += (
595 exp(-((double)(u*u+v*v))/sigma)
596 /* / (MagickPI*sigma) */ );
598 kernel->maximum = kernel->values[
599 kernel->y*kernel->width+kernel->x ];
601 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
607 sigma = fabs(args->sigma);
609 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
611 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
612 kernel->x = (long) (kernel->width-1)/2;
615 kernel->negative_range = kernel->positive_range = 0.0;
616 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
617 kernel->height*sizeof(double));
618 if (kernel->values == (double *) NULL)
619 return(DestroyKernelInfo(kernel));
623 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
624 ** It generates a gaussian 3 times the width, and compresses it into
625 ** the expected range. This produces a closer normalization of the
626 ** resulting kernel, especially for very low sigma values.
627 ** As such while wierd it is prefered.
629 ** I am told this method originally came from Photoshop.
631 sigma *= KernelRank; /* simplify expanded curve */
632 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
633 (void) ResetMagickMemory(kernel->values,0, (size_t)
634 kernel->width*sizeof(double));
635 for ( u=-v; u <= v; u++) {
636 kernel->values[(u+v)/KernelRank] +=
637 exp(-((double)(u*u))/(2.0*sigma*sigma))
638 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
640 for (i=0; i < (long) kernel->width; i++)
641 kernel->positive_range += kernel->values[i];
643 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
644 kernel->positive_range += (
646 exp(-((double)(u*u))/(2.0*sigma*sigma))
647 /* / (MagickSQ2PI*sigma) */ );
650 kernel->maximum = kernel->values[ kernel->x ];
651 /* Note that neither methods above generate a normalized kernel,
652 ** though it gets close. The kernel may be 'clipped' by a user defined
653 ** radius, producing a smaller (darker) kernel. Also for very small
654 ** sigma's (> 0.1) the central value becomes larger than one, and thus
655 ** producing a very bright kernel.
658 /* Normalize the 1D Gaussian Kernel
660 ** Because of this the divisor in the above kernel generator is
661 ** not needed, so is not done above.
663 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
665 /* rotate the kernel by given angle */
666 RotateKernelInfo(kernel, args->xi);
671 sigma = fabs(args->sigma);
673 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
675 if ( args->rho < 1.0 )
676 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
678 kernel->width = (unsigned long)args->rho;
679 kernel->x = kernel->y = 0;
681 kernel->negative_range = kernel->positive_range = 0.0;
682 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
683 kernel->height*sizeof(double));
684 if (kernel->values == (double *) NULL)
685 return(DestroyKernelInfo(kernel));
687 /* A comet blur is half a gaussian curve, so that the object is
688 ** blurred in one direction only. This may not be quite the right
689 ** curve so may change in the future. The function must be normalised.
693 sigma *= KernelRank; /* simplify expanded curve */
694 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
695 (void) ResetMagickMemory(kernel->values,0, (size_t)
696 kernel->width*sizeof(double));
697 for ( u=0; u < v; u++) {
698 kernel->values[u/KernelRank] +=
699 exp(-((double)(u*u))/(2.0*sigma*sigma))
700 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
702 for (i=0; i < (long) kernel->width; i++)
703 kernel->positive_range += kernel->values[i];
705 for ( i=0; i < (long) kernel->width; i++)
706 kernel->positive_range += (
708 exp(-((double)(i*i))/(2.0*sigma*sigma))
709 /* / (MagickSQ2PI*sigma) */ );
712 kernel->maximum = kernel->values[0];
714 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
715 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
718 /* Boolean Kernels */
719 case RectangleKernel:
723 if ( type == SquareKernel )
726 kernel->width = kernel->height = 3; /* default radius = 1 */
728 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
729 kernel->x = kernel->y = (long) (kernel->width-1)/2;
733 /* NOTE: user defaults set in "AcquireKernelInfo()" */
734 if ( args->rho < 1.0 || args->sigma < 1.0 )
735 return(DestroyKernelInfo(kernel)); /* invalid args given */
736 kernel->width = (unsigned long)args->rho;
737 kernel->height = (unsigned long)args->sigma;
738 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
739 args->psi < 0.0 || args->psi > (double)kernel->height )
740 return(DestroyKernelInfo(kernel)); /* invalid args given */
741 kernel->x = (long) args->xi;
742 kernel->y = (long) args->psi;
745 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
746 kernel->height*sizeof(double));
747 if (kernel->values == (double *) NULL)
748 return(DestroyKernelInfo(kernel));
750 /* set all kernel values to 1.0 */
751 u=(long) kernel->width*kernel->height;
752 for ( i=0; i < u; i++)
753 kernel->values[i] = scale;
754 kernel->minimum = kernel->maximum = scale; /* a flat shape */
755 kernel->positive_range = scale*u;
761 kernel->width = kernel->height = 3; /* default radius = 1 */
763 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
764 kernel->x = kernel->y = (long) (kernel->width-1)/2;
766 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
767 kernel->height*sizeof(double));
768 if (kernel->values == (double *) NULL)
769 return(DestroyKernelInfo(kernel));
771 /* set all kernel values within diamond area to scale given */
772 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
773 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
774 if ((labs(u)+labs(v)) <= (long)kernel->x)
775 kernel->positive_range += kernel->values[i] = args->sigma;
777 kernel->values[i] = nan;
778 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
786 limit = (long)(args->rho*args->rho);
787 if (args->rho < 0.1) /* default radius approx 3.5 */
788 kernel->width = kernel->height = 7L, limit = 10L;
790 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
791 kernel->x = kernel->y = (long) (kernel->width-1)/2;
793 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
794 kernel->height*sizeof(double));
795 if (kernel->values == (double *) NULL)
796 return(DestroyKernelInfo(kernel));
798 /* set all kernel values within disk area to 1.0 */
799 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
800 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
801 if ((u*u+v*v) <= limit)
802 kernel->positive_range += kernel->values[i] = args->sigma;
804 kernel->values[i] = nan;
805 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
811 kernel->width = kernel->height = 5; /* default radius 2 */
813 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
814 kernel->x = kernel->y = (long) (kernel->width-1)/2;
816 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
817 kernel->height*sizeof(double));
818 if (kernel->values == (double *) NULL)
819 return(DestroyKernelInfo(kernel));
821 /* set all kernel values along axises to 1.0 */
822 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
823 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
824 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
825 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
826 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
829 /* Distance Measuring Kernels */
830 case ChebyshevKernel:
836 kernel->width = kernel->height = 3; /* default radius = 1 */
838 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
839 kernel->x = kernel->y = (long) (kernel->width-1)/2;
841 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
842 kernel->height*sizeof(double));
843 if (kernel->values == (double *) NULL)
844 return(DestroyKernelInfo(kernel));
846 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
847 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
848 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
849 kernel->positive_range += ( kernel->values[i] =
850 scale*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
851 kernel->maximum = kernel->values[0];
854 case ManhattenKernel:
860 kernel->width = kernel->height = 3; /* default radius = 1 */
862 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
863 kernel->x = kernel->y = (long) (kernel->width-1)/2;
865 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
866 kernel->height*sizeof(double));
867 if (kernel->values == (double *) NULL)
868 return(DestroyKernelInfo(kernel));
870 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
871 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
872 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
873 kernel->positive_range += ( kernel->values[i] =
874 scale*(labs(u)+labs(v)) );
875 kernel->maximum = kernel->values[0];
878 case EuclideanKernel:
884 kernel->width = kernel->height = 3; /* default radius = 1 */
886 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
887 kernel->x = kernel->y = (long) (kernel->width-1)/2;
889 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
890 kernel->height*sizeof(double));
891 if (kernel->values == (double *) NULL)
892 return(DestroyKernelInfo(kernel));
894 scale = (args->sigma < 1.0) ? 100.0 : args->sigma;
895 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
896 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
897 kernel->positive_range += ( kernel->values[i] =
898 scale*sqrt((double)(u*u+v*v)) );
899 kernel->maximum = kernel->values[0];
902 /* Undefined Kernels */
903 case LaplacianKernel:
906 perror("Kernel Type has not been defined yet");
909 /* Generate a No-Op minimal kernel - 1x1 pixel */
910 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
911 if (kernel->values == (double *) NULL)
912 return(DestroyKernelInfo(kernel));
913 kernel->width = kernel->height = 1;
914 kernel->x = kernel->x = 0;
915 kernel->type = UndefinedKernel;
917 kernel->positive_range =
918 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
926 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
930 + C l o n e K e r n e l I n f o %
934 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
936 % CloneKernelInfo() creates a new clone of the given Kernel so that its can
937 % be modified without effecting the original. The cloned kernel should be
938 % destroyed using DestoryKernelInfo() when no longer needed.
940 % The format of the DestroyKernelInfo method is:
942 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
944 % A description of each parameter follows:
946 % o kernel: the Morphology/Convolution kernel to be cloned
950 static KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
958 assert(kernel != (KernelInfo *) NULL);
960 new=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
961 if (new == (KernelInfo *) NULL)
963 *new = *kernel; /* copy values in structure */
965 new->values=(double *) AcquireQuantumMemory(kernel->width,
966 kernel->height*sizeof(double));
967 if (new->values == (double *) NULL)
968 return(DestroyKernelInfo(new));
970 for (i=0; i < (long) (kernel->width*kernel->height); i++)
971 new->values[i] = kernel->values[i];
977 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
981 % D e s t r o y K e r n e l I n f o %
985 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
987 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
990 % The format of the DestroyKernelInfo method is:
992 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
994 % A description of each parameter follows:
996 % o kernel: the Morphology/Convolution kernel to be destroyed
1000 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1002 assert(kernel != (KernelInfo *) NULL);
1004 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1005 kernel->height*sizeof(double));
1006 kernel->values=(double *)RelinquishMagickMemory(kernel->values);
1007 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
1012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1016 % M o r p h o l o g y I m a g e C h a n n e l %
1020 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1022 % MorphologyImageChannel() applies a user supplied kernel to the image
1023 % according to the given mophology method.
1025 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1026 % by the kernel generator.
1028 % The format of the MorphologyImage method is:
1030 % Image *MorphologyImage(const Image *image, MorphologyMethod method,
1031 % const long iterations, KernelInfo *kernel, ExceptionInfo *exception)
1032 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1033 % channel, MorphologyMethod method, const long iterations, KernelInfo
1034 % *kernel, ExceptionInfo *exception)
1036 % A description of each parameter follows:
1038 % o image: the image.
1040 % o method: the morphology method to be applied.
1042 % o iterations: apply the operation this many times (or no change).
1043 % A value of -1 means loop until no change found.
1044 % How this is applied may depend on the morphology method.
1045 % Typically this is a value of 1.
1047 % o channel: the channel type.
1049 % o kernel: An array of double representing the morphology kernel.
1050 % Warning: kernel may be normalized for the Convolve method.
1052 % o exception: return any errors or warnings in this structure.
1055 % TODO: bias and auto-scale handling of the kernel for convolution
1056 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1057 % by the kernel generator.
1062 /* Internal function
1063 * Apply the Low-Level Morphology Method using the given Kernel
1064 * Returning the number of pixels that changed.
1065 * Two pre-created images must be provided, no image is created.
1067 static unsigned long MorphologyApply(const Image *image, Image
1068 *result_image, const MorphologyMethod method, const ChannelType channel,
1069 const KernelInfo *kernel, ExceptionInfo *exception)
1071 #define MorphologyTag "Morphology/Image"
1088 /* Only the most basic morphology is actually performed by this routine */
1091 Apply Basic Morphology to Image.
1097 GetMagickPixelPacket(image,&bias);
1098 SetMagickPixelPacketBias(image,&bias);
1099 /* Future: handle auto-bias from user, based on kernel input */
1101 p_view=AcquireCacheView(image);
1102 q_view=AcquireCacheView(result_image);
1104 /* Some methods (including convolve) needs use a reflected kernel.
1105 * Adjust 'origin' offsets for this reflected kernel.
1110 case ErodeMorphology:
1111 case ErodeIntensityMorphology:
1112 /* kernel is user as is, without reflection */
1114 case ConvolveMorphology:
1115 case DilateMorphology:
1116 case DilateIntensityMorphology:
1117 case DistanceMorphology:
1118 /* kernel needs to used with reflection */
1119 offx = (long) kernel->width-offx-1;
1120 offy = (long) kernel->height-offy-1;
1123 perror("Not a low level Morpholgy Method");
1127 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1128 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1130 for (y=0; y < (long) image->rows; y++)
1135 register const PixelPacket
1138 register const IndexPacket
1139 *restrict p_indexes;
1141 register PixelPacket
1144 register IndexPacket
1145 *restrict q_indexes;
1153 if (status == MagickFalse)
1155 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1156 image->columns+kernel->width, kernel->height, exception);
1157 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1159 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1164 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1165 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1166 r = (image->columns+kernel->width)*offy+offx; /* constant */
1168 for (x=0; x < (long) image->columns; x++)
1176 register const double
1179 register const PixelPacket
1182 register const IndexPacket
1183 *restrict k_indexes;
1188 /* Copy input to ouput image for unused channels
1189 * This removes need for 'cloning' a new image every iteration
1192 if (image->colorspace == CMYKColorspace)
1193 q_indexes[x] = p_indexes[r];
1195 result.green=(MagickRealType) 0;
1196 result.blue=(MagickRealType) 0;
1197 result.opacity=(MagickRealType) 0;
1198 result.index=(MagickRealType) 0;
1200 case ConvolveMorphology:
1201 /* Set the user defined bias of the weighted average output
1203 ** FUTURE: provide some way for internal functions to disable
1204 ** user defined bias and scaling effects.
1208 case DilateMorphology:
1213 result.index = -MagickHuge;
1215 case ErodeMorphology:
1220 result.index = +MagickHuge;
1222 case DilateIntensityMorphology:
1223 case ErodeIntensityMorphology:
1224 result.red = 0.0; /* flag indicating first match found */
1227 /* Otherwise just start with the original pixel value */
1228 result.red = (MagickRealType) p[r].red;
1229 result.green = (MagickRealType) p[r].green;
1230 result.blue = (MagickRealType) p[r].blue;
1231 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1232 if ( image->colorspace == CMYKColorspace)
1233 result.index = (MagickRealType) p_indexes[r];
1238 case ConvolveMorphology:
1239 /* Weighted Average of pixels using reflected kernel
1241 ** NOTE for correct working of this operation for asymetrical
1242 ** kernels, the kernel needs to be applied in its reflected form.
1243 ** That is its values needs to be reversed.
1245 ** Correlation is actually the same as this but without reflecting
1246 ** the kernel, and thus 'lower-level' that Convolution. However
1247 ** as Convolution is the more common method used, and it does not
1248 ** really cost us much in terms of processing to use a reflected
1249 ** kernel it is Convolution that is implemented.
1251 ** Correlation will have its kernel reflected before calling
1252 ** this function to do a Convolve.
1254 ** For more details of Correlation vs Convolution see
1255 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
1257 if (((channel & OpacityChannel) == 0) ||
1258 (image->matte == MagickFalse))
1260 /* Convolution without transparency effects */
1261 k = &kernel->values[ kernel->width*kernel->height-1 ];
1263 k_indexes = p_indexes;
1264 for (v=0; v < (long) kernel->height; v++) {
1265 for (u=0; u < (long) kernel->width; u++, k--) {
1266 if ( IsNan(*k) ) continue;
1267 result.red += (*k)*k_pixels[u].red;
1268 result.green += (*k)*k_pixels[u].green;
1269 result.blue += (*k)*k_pixels[u].blue;
1270 /* result.opacity += not involved here */
1271 if ( image->colorspace == CMYKColorspace)
1272 result.index += (*k)*k_indexes[u];
1274 k_pixels += image->columns+kernel->width;
1275 k_indexes += image->columns+kernel->width;
1279 { /* Kernel & Alpha weighted Convolution */
1281 alpha, /* alpha value * kernel weighting */
1282 gamma; /* weighting divisor */
1285 k = &kernel->values[ kernel->width*kernel->height-1 ];
1287 k_indexes = p_indexes;
1288 for (v=0; v < (long) kernel->height; v++) {
1289 for (u=0; u < (long) kernel->width; u++, k--) {
1290 if ( IsNan(*k) ) continue;
1291 alpha=(*k)*(QuantumScale*(QuantumRange-
1292 k_pixels[u].opacity));
1294 result.red += alpha*k_pixels[u].red;
1295 result.green += alpha*k_pixels[u].green;
1296 result.blue += alpha*k_pixels[u].blue;
1297 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1298 if ( image->colorspace == CMYKColorspace)
1299 result.index += alpha*k_indexes[u];
1301 k_pixels += image->columns+kernel->width;
1302 k_indexes += image->columns+kernel->width;
1304 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1305 result.red *= gamma;
1306 result.green *= gamma;
1307 result.blue *= gamma;
1308 result.opacity *= gamma;
1309 result.index *= gamma;
1313 case ErodeMorphology:
1314 /* Minimize Value within kernel neighbourhood
1316 ** NOTE that the kernel is not reflected for this operation!
1318 ** NOTE: in normal Greyscale Morphology, the kernel value should
1319 ** be added to the real value, this is currently not done, due to
1320 ** the nature of the boolean kernels being used.
1324 k_indexes = p_indexes;
1325 for (v=0; v < (long) kernel->height; v++) {
1326 for (u=0; u < (long) kernel->width; u++, k++) {
1327 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1328 Minimize(result.red, (double) k_pixels[u].red);
1329 Minimize(result.green, (double) k_pixels[u].green);
1330 Minimize(result.blue, (double) k_pixels[u].blue);
1331 Minimize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1332 if ( image->colorspace == CMYKColorspace)
1333 Minimize(result.index, (double) k_indexes[u]);
1335 k_pixels += image->columns+kernel->width;
1336 k_indexes += image->columns+kernel->width;
1340 case DilateMorphology:
1341 /* Maximize Value within kernel neighbourhood
1343 ** NOTE for correct working of this operation for asymetrical
1344 ** kernels, the kernel needs to be applied in its reflected form.
1345 ** That is its values needs to be reversed.
1347 ** NOTE: in normal Greyscale Morphology, the kernel value should
1348 ** be added to the real value, this is currently not done, due to
1349 ** the nature of the boolean kernels being used.
1352 k = &kernel->values[ kernel->width*kernel->height-1 ];
1354 k_indexes = p_indexes;
1355 for (v=0; v < (long) kernel->height; v++) {
1356 for (u=0; u < (long) kernel->width; u++, k--) {
1357 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1358 Maximize(result.red, (double) k_pixels[u].red);
1359 Maximize(result.green, (double) k_pixels[u].green);
1360 Maximize(result.blue, (double) k_pixels[u].blue);
1361 Maximize(result.opacity, QuantumRange-(double) k_pixels[u].opacity);
1362 if ( image->colorspace == CMYKColorspace)
1363 Maximize(result.index, (double) k_indexes[u]);
1365 k_pixels += image->columns+kernel->width;
1366 k_indexes += image->columns+kernel->width;
1370 case ErodeIntensityMorphology:
1371 /* Select Pixel with Minimum Intensity within kernel neighbourhood
1373 ** WARNING: the intensity test fails for CMYK and does not
1374 ** take into account the moderating effect of teh alpha channel
1375 ** on the intensity.
1377 ** NOTE that the kernel is not reflected for this operation!
1381 k_indexes = p_indexes;
1382 for (v=0; v < (long) kernel->height; v++) {
1383 for (u=0; u < (long) kernel->width; u++, k++) {
1384 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1385 if ( result.red == 0.0 ||
1386 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1387 /* copy the whole pixel - no channel selection */
1389 if ( result.red > 0.0 ) changed++;
1393 k_pixels += image->columns+kernel->width;
1394 k_indexes += image->columns+kernel->width;
1398 case DilateIntensityMorphology:
1399 /* Select Pixel with Maximum Intensity within kernel neighbourhood
1401 ** WARNING: the intensity test fails for CMYK and does not
1402 ** take into account the moderating effect of teh alpha channel
1403 ** on the intensity.
1405 ** NOTE for correct working of this operation for asymetrical
1406 ** kernels, the kernel needs to be applied in its reflected form.
1407 ** That is its values needs to be reversed.
1409 k = &kernel->values[ kernel->width*kernel->height-1 ];
1411 k_indexes = p_indexes;
1412 for (v=0; v < (long) kernel->height; v++) {
1413 for (u=0; u < (long) kernel->width; u++, k--) {
1414 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1415 if ( result.red == 0.0 ||
1416 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1417 /* copy the whole pixel - no channel selection */
1419 if ( result.red > 0.0 ) changed++;
1423 k_pixels += image->columns+kernel->width;
1424 k_indexes += image->columns+kernel->width;
1428 case DistanceMorphology:
1429 /* Add kernel Value and select the minimum value found.
1430 ** The result is a iterative distance from edge of image shape.
1432 ** All Distance Kernels are symetrical, but that may not always
1433 ** be the case. For example how about a distance from left edges?
1434 ** To work correctly with asymetrical kernels the reflected kernel
1435 ** needs to be applied.
1438 /* No need to do distance morphology if original value is zero
1439 ** Unfortunatally I have not been able to get this right
1440 ** when channel selection also becomes involved. -- Arrgghhh
1442 if ( ((channel & RedChannel) == 0 && p[r].red == 0)
1443 || ((channel & GreenChannel) == 0 && p[r].green == 0)
1444 || ((channel & BlueChannel) == 0 && p[r].blue == 0)
1445 || ((channel & OpacityChannel) == 0 && p[r].opacity == 0)
1446 || (( (channel & IndexChannel) == 0
1447 || image->colorspace != CMYKColorspace
1448 ) && p_indexes[x] ==0 )
1452 k = &kernel->values[ kernel->width*kernel->height-1 ];
1454 k_indexes = p_indexes;
1455 for (v=0; v < (long) kernel->height; v++) {
1456 for (u=0; u < (long) kernel->width; u++, k--) {
1457 if ( IsNan(*k) ) continue;
1458 Minimize(result.red, (*k)+k_pixels[u].red);
1459 Minimize(result.green, (*k)+k_pixels[u].green);
1460 Minimize(result.blue, (*k)+k_pixels[u].blue);
1461 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1462 if ( image->colorspace == CMYKColorspace)
1463 Minimize(result.index, (*k)+k_indexes[u]);
1465 k_pixels += image->columns+kernel->width;
1466 k_indexes += image->columns+kernel->width;
1470 case UndefinedMorphology:
1472 break; /* Do nothing */
1475 case UndefinedMorphology:
1476 case DilateIntensityMorphology:
1477 case ErodeIntensityMorphology:
1478 break; /* full pixel was directly assigned - not a channel method */
1480 /* Assign the results */
1481 if ((channel & RedChannel) != 0)
1482 q->red = ClampToQuantum(result.red);
1483 if ((channel & GreenChannel) != 0)
1484 q->green = ClampToQuantum(result.green);
1485 if ((channel & BlueChannel) != 0)
1486 q->blue = ClampToQuantum(result.blue);
1487 if ((channel & OpacityChannel) != 0
1488 && image->matte == MagickTrue )
1489 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1490 if ((channel & IndexChannel) != 0
1491 && image->colorspace == CMYKColorspace)
1492 q_indexes[x] = ClampToQuantum(result.index);
1495 if ( ( p[r].red != q->red )
1496 || ( p[r].green != q->green )
1497 || ( p[r].blue != q->blue )
1498 || ( p[r].opacity != q->opacity )
1499 || ( image->colorspace == CMYKColorspace &&
1500 p_indexes[r] != q_indexes[x] ) )
1501 changed++; /* The pixel had some value changed! */
1505 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1506 if (sync == MagickFalse)
1508 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1513 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1514 #pragma omp critical (MagickCore_MorphologyImage)
1516 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1517 if (proceed == MagickFalse)
1521 result_image->type=image->type;
1522 q_view=DestroyCacheView(q_view);
1523 p_view=DestroyCacheView(p_view);
1524 return(status ? (unsigned long) changed : 0);
1528 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1529 method, const long iterations,const KernelInfo *kernel, ExceptionInfo
1535 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1536 iterations,kernel,exception);
1537 return(morphology_image);
1541 MagickExport Image *MorphologyImageChannel(const Image *image, const
1542 ChannelType channel, const MorphologyMethod method, const long
1543 iterations, const KernelInfo *kernel, ExceptionInfo *exception)
1566 assert(image != (Image *) NULL);
1567 assert(image->signature == MagickSignature);
1568 assert(kernel != (KernelInfo *) NULL);
1569 assert(kernel->signature == MagickSignature);
1570 assert(exception != (ExceptionInfo *) NULL);
1571 assert(exception->signature == MagickSignature);
1573 if ( iterations == 0 )
1574 return((Image *)NULL); /* null operation - nothing to do! */
1576 /* kernel must be valid at this point
1577 * (except maybe for posible future morphology methods like "Prune"
1579 assert(kernel != (KernelInfo *)NULL);
1581 count = 0; /* interation count */
1582 changed = 1; /* if compound method assume image was changed */
1583 curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
1584 curr_method = method; /* to be changed as nessary */
1586 limit = (unsigned long) iterations;
1587 if ( iterations < 0 )
1588 limit = image->columns > image->rows ? image->columns : image->rows;
1590 /* Third-level morphology methods */
1591 grad_image=(Image *) NULL;
1592 switch( curr_method ) {
1593 case EdgeMorphology:
1594 grad_image = MorphologyImageChannel(image, channel,
1595 DilateMorphology, iterations, curr_kernel, exception);
1597 case EdgeInMorphology:
1598 curr_method = ErodeMorphology;
1600 case EdgeOutMorphology:
1601 curr_method = DilateMorphology;
1603 case TopHatMorphology:
1604 curr_method = OpenMorphology;
1606 case BottomHatMorphology:
1607 curr_method = CloseMorphology;
1610 break; /* not a third-level method */
1613 /* Second-level morphology methods */
1614 switch( curr_method ) {
1615 case OpenMorphology:
1616 /* Open is a Erode then a Dilate without reflection */
1617 new_image = MorphologyImageChannel(image, channel,
1618 ErodeMorphology, iterations, curr_kernel, exception);
1619 if (new_image == (Image *) NULL)
1620 return((Image *) NULL);
1621 curr_method = DilateMorphology;
1623 case OpenIntensityMorphology:
1624 new_image = MorphologyImageChannel(image, channel,
1625 ErodeIntensityMorphology, iterations, curr_kernel, exception);
1626 if (new_image == (Image *) NULL)
1627 return((Image *) NULL);
1628 curr_method = DilateIntensityMorphology;
1631 case CloseMorphology:
1632 /* Close is a Dilate then Erode using reflected kernel */
1633 /* A reflected kernel is needed for a Close */
1634 if ( curr_kernel == kernel )
1635 curr_kernel = CloneKernelInfo(kernel);
1636 RotateKernelInfo(curr_kernel,180);
1637 new_image = MorphologyImageChannel(image, channel,
1638 DilateMorphology, iterations, curr_kernel, exception);
1639 if (new_image == (Image *) NULL)
1640 return((Image *) NULL);
1641 curr_method = ErodeMorphology;
1643 case CloseIntensityMorphology:
1644 /* A reflected kernel is needed for a Close */
1645 if ( curr_kernel == kernel )
1646 curr_kernel = CloneKernelInfo(kernel);
1647 RotateKernelInfo(curr_kernel,180);
1648 new_image = MorphologyImageChannel(image, channel,
1649 DilateIntensityMorphology, iterations, curr_kernel, exception);
1650 if (new_image == (Image *) NULL)
1651 return((Image *) NULL);
1652 curr_method = ErodeIntensityMorphology;
1655 case CorrelateMorphology:
1656 /* A Correlation is actually a Convolution with a reflected kernel.
1657 ** However a Convolution is a weighted sum with a reflected kernel.
1658 ** It may seem stange to convert a Correlation into a Convolution
1659 ** as the Correleation is the simplier method, but Convolution is
1660 ** much more commonly used, and it makes sense to implement it directly
1661 ** so as to avoid the need to duplicate the kernel when it is not
1662 ** required (which is typically the default).
1664 if ( curr_kernel == kernel )
1665 curr_kernel = CloneKernelInfo(kernel);
1666 RotateKernelInfo(curr_kernel,180);
1667 curr_method = ConvolveMorphology;
1668 /* FALL-THRU into Correlate (weigthed sum without reflection) */
1670 case ConvolveMorphology:
1671 /* Scale or Normalize kernel, according to user wishes
1672 ** before using it for the Convolve/Correlate method.
1674 ** FUTURE: provide some way for internal functions to disable
1675 ** user bias and scaling effects.
1677 artifact = GetImageArtifact(image,"convolve:scale");
1678 if ( artifact != (char *)NULL ) {
1684 if ( curr_kernel == kernel )
1685 curr_kernel = CloneKernelInfo(kernel);
1688 flags = ParseGeometry(artifact, &args);
1689 ScaleKernelInfo(curr_kernel, args.rho, flags);
1691 /* FALL-THRU to do the first, and typically the only iteration */
1694 /* Do a single iteration using the Low-Level Morphology method!
1695 ** This ensures a "new_image" has been generated, but allows us to skip
1696 ** the creation of 'old_image' if no more iterations are needed.
1698 ** The "curr_method" should also be set to a low-level method that is
1699 ** understood by the MorphologyApply() internal function.
1701 new_image=CloneImage(image,0,0,MagickTrue,exception);
1702 if (new_image == (Image *) NULL)
1703 return((Image *) NULL);
1704 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1706 InheritException(exception,&new_image->exception);
1707 new_image=DestroyImage(new_image);
1708 return((Image *) NULL);
1710 changed = MorphologyApply(image,new_image,curr_method,channel,curr_kernel,
1713 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1714 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1715 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1720 /* At this point the "curr_method" should not only be set to a low-level
1721 ** method that is understood by the MorphologyApply() internal function,
1722 ** but "new_image" should now be defined, as the image to apply the
1723 ** "curr_method" to.
1726 /* Repeat the low-level morphology until count or no change reached */
1727 if ( count < (long) limit && changed > 0 ) {
1728 old_image = CloneImage(new_image,0,0,MagickTrue,exception);
1729 if (old_image == (Image *) NULL)
1730 return(DestroyImage(new_image));
1731 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1733 InheritException(exception,&old_image->exception);
1734 old_image=DestroyImage(old_image);
1735 return(DestroyImage(new_image));
1737 while( count < (long) limit && changed != 0 )
1739 Image *tmp = old_image;
1740 old_image = new_image;
1742 changed = MorphologyApply(old_image,new_image,curr_method,channel,
1743 curr_kernel, exception);
1745 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1746 fprintf(stderr, "Morphology %s:%ld => Changed %lu\n",
1747 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1750 old_image=DestroyImage(old_image);
1753 /* We are finished with kernel - destroy it if we made a clone */
1754 if ( curr_kernel != kernel )
1755 curr_kernel=DestroyKernelInfo(curr_kernel);
1757 /* Third-level Subtractive methods post-processing */
1759 case EdgeOutMorphology:
1760 case EdgeInMorphology:
1761 case TopHatMorphology:
1762 case BottomHatMorphology:
1763 /* Get Difference relative to the original image */
1764 CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1767 case EdgeMorphology: /* subtract the Erode from a Dilate */
1768 CompositeImageChannel(new_image, channel, DifferenceCompositeOp,
1770 grad_image=DestroyImage(grad_image);
1780 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1784 + R o t a t e K e r n e l I n f o %
1788 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1790 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
1791 % restricted to 90 degree angles, but this may be improved in the future.
1793 % The format of the RotateKernelInfo method is:
1795 % void RotateKernelInfo(KernelInfo *kernel, double angle)
1797 % A description of each parameter follows:
1799 % o kernel: the Morphology/Convolution kernel
1801 % o angle: angle to rotate in degrees
1803 % This function is only internel to this module, as it is not finalized,
1804 % especially with regard to non-orthogonal angles, and rotation of larger
1807 static void RotateKernelInfo(KernelInfo *kernel, double angle)
1809 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
1811 ** TODO: expand beyond simple 90 degree rotates, flips and flops
1814 /* Modulus the angle */
1815 angle = fmod(angle, 360.0);
1819 if ( 315.0 < angle || angle <= 45.0 )
1820 return; /* no change! - At least at this time */
1822 switch (kernel->type) {
1823 /* These built-in kernels are cylindrical kernels, rotating is useless */
1824 case GaussianKernel:
1825 case LaplacianKernel:
1829 case ChebyshevKernel:
1830 case ManhattenKernel:
1831 case EuclideanKernel:
1834 /* These may be rotatable at non-90 angles in the future */
1835 /* but simply rotating them in multiples of 90 degrees is useless */
1841 /* These only allows a +/-90 degree rotation (by transpose) */
1842 /* A 180 degree rotation is useless */
1844 case RectangleKernel:
1845 if ( 135.0 < angle && angle <= 225.0 )
1847 if ( 225.0 < angle && angle <= 315.0 )
1851 /* these are freely rotatable in 90 degree units */
1853 case UndefinedKernel:
1854 case UserDefinedKernel:
1857 if ( 135.0 < angle && angle <= 225.0 )
1859 /* For a 180 degree rotation - also know as a reflection */
1860 /* This is actually a very very common operation! */
1861 /* Basically all that is needed is a reversal of the kernel data! */
1868 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
1869 t=k[i], k[i]=k[j], k[j]=t;
1871 kernel->x = (long) kernel->width - kernel->x - 1;
1872 kernel->y = (long) kernel->height - kernel->y - 1;
1873 angle = fmod(angle+180.0, 360.0);
1875 if ( 45.0 < angle && angle <= 135.0 )
1876 { /* Do a transpose and a flop, of the image, which results in a 90
1877 * degree rotation using two mirror operations.
1879 * WARNING: this assumes the original image was a 1 dimentional image
1880 * but currently that is the only built-ins it is applied to.
1884 t = (long) kernel->width;
1885 kernel->width = kernel->height;
1886 kernel->height = (unsigned long) t;
1888 kernel->x = kernel->y;
1890 angle = fmod(450.0 - angle, 360.0);
1892 /* At this point angle should be between -45 (315) and +45 degrees
1893 * In the future some form of non-orthogonal angled rotates could be
1894 * performed here, posibily with a linear kernel restriction.
1898 Not currently in use!
1899 { /* Do a flop, this assumes kernel is horizontally symetrical.
1900 * Each row of the kernel needs to be reversed!
1909 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
1910 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
1911 t=k[x], k[x]=k[r], k[r]=t;
1913 kernel->x = kernel->width - kernel->x - 1;
1914 angle = fmod(angle+180.0, 360.0);
1921 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1925 + S c a l e K e r n e l I n f o %
1929 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1931 % ScaleKernelInfo() scales the kernel by the given amount, with or without
1932 % normalization of the sum of the kernel values.
1934 % By default (no flags given) the values within the kernel is scaled
1935 % according the given scaling factor.
1937 % If any 'normalize_flags' are given the kernel will be normalized and then
1938 % further scaled by the scaleing factor value given. A 'PercentValue' flag
1939 % will cause the given scaling factor to be divided by one hundred percent.
1941 % Kernel normalization ('normalize_flags' given) is designed to ensure that
1942 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
1943 % morphology methods will fall into -1.0 to +1.0 range. Note however that
1944 % for non-HDRI versions of IM this may cause images to have any negative
1945 % results clipped, unless some 'clip' any negative output from 'Convolve'
1946 % with the use of some kernels.
1948 % More specifically. Kernels which only contain positive values (such as a
1949 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
1950 % ensuring a 0.0 to +1.0 convolution output range for non-HDRI images.
1952 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
1953 % the kernel will be scaled by the absolute of the sum of kernel values, so
1954 % that it will generally fall within the +/- 1.0 range.
1956 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
1957 % will be scaled by just the sum of the postive values, so that its output
1958 % range will again fall into the +/- 1.0 range.
1960 % For special kernels designed for locating shapes using 'Correlate', (often
1961 % only containing +1 and -1 values, representing foreground/brackground
1962 % matching) a special normalization method is provided to scale the positive
1963 % values seperatally to those of the negative values, so the kernel will be
1964 % forced to become a zero-sum kernel better suited to such searches.
1966 % WARNING: Correct normalization of the kernal assumes that the '*_range'
1967 % attributes within the kernel structure have been correctly set during the
1970 % NOTE: The values used for 'normalize_flags' have been selected specifically
1971 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
1972 % means CorrelateNormalizeValue, and '%' means PercentValue. All other
1973 % GeometryFlags values are ignored.
1975 % The format of the ScaleKernelInfo method is:
1977 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
1978 % const MagickStatusType normalize_flags )
1980 % A description of each parameter follows:
1982 % o kernel: the Morphology/Convolution kernel
1985 % multiply all values (after normalization) by this factor if not
1986 % zero. If the kernel is normalized regardless of any flags.
1988 % o normalize_flags:
1989 % GeometryFlags defining normalization method to use.
1990 % specifically: NormalizeValue, CorrelateNormalizeValue,
1991 % and/or PercentValue
1993 % This function is internal to this module only at this time, but can be
1994 % exported to other modules if needed.
1996 static void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
1997 const GeometryFlags normalize_flags)
2007 if ( (normalize_flags&NormalizeValue) != 0 ) {
2008 /* normalize kernel appropriately */
2009 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
2010 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
2012 pos_scale = kernel->positive_range; /* special zero-summing kernel */
2014 /* force kernel into being a normalized zero-summing kernel */
2015 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
2016 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
2017 ? kernel->positive_range : 1.0;
2018 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
2019 ? -kernel->negative_range : 1.0;
2022 neg_scale = pos_scale;
2024 /* finialize scaling_factor for positive and negative components */
2025 pos_scale = scaling_factor/pos_scale;
2026 neg_scale = scaling_factor/neg_scale;
2027 if ( (normalize_flags&PercentValue) != 0 ) {
2032 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2033 if ( ! IsNan(kernel->values[i]) )
2034 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
2036 /* convolution output range */
2037 kernel->positive_range *= pos_scale;
2038 kernel->negative_range *= neg_scale;
2039 /* maximum and minimum values in kernel */
2040 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
2041 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
2043 /* swap kernel settings if user scaling factor is negative */
2044 if ( scaling_factor < MagickEpsilon ) {
2046 t = kernel->positive_range;
2047 kernel->positive_range = kernel->negative_range;
2048 kernel->negative_range = t;
2049 t = kernel->maximum;
2050 kernel->maximum = kernel->minimum;
2051 kernel->minimum = 1;
2058 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2062 + S h o w K e r n e l I n f o %
2066 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2068 % ShowKernelInfo() outputs the details of the given kernel defination to
2069 % standard error, generally due to a users 'showkernel' option request.
2071 % The format of the ShowKernel method is:
2073 % void ShowKernelInfo(KernelInfo *kernel)
2075 % A description of each parameter follows:
2077 % o kernel: the Morphology/Convolution kernel
2079 % This function is internal to this module only at this time. That may change
2082 MagickExport void ShowKernelInfo(KernelInfo *kernel)
2088 "Kernel \"%s\" of size %lux%lu%+ld%+ld with values from %.*lg to %.*lg\n",
2089 MagickOptionToMnemonic(MagickKernelOptions, kernel->type),
2090 kernel->width, kernel->height,
2091 kernel->x, kernel->y,
2092 GetMagickPrecision(), kernel->minimum,
2093 GetMagickPrecision(), kernel->maximum);
2094 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2095 GetMagickPrecision(), kernel->negative_range,
2096 GetMagickPrecision(), kernel->positive_range,
2097 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2098 for (i=v=0; v < (long) kernel->height; v++) {
2099 fprintf(stderr,"%2ld:",v);
2100 for (u=0; u < (long) kernel->width; u++, i++)
2101 if ( IsNan(kernel->values[i]) )
2102 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2104 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2105 GetMagickPrecision(), kernel->values[i]);
2106 fprintf(stderr,"\n");
2111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2115 + Z e r o K e r n e l N a n s %
2119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2121 % ZeroKernelNans() replaces any special 'nan' value that may be present in
2122 % the kernel with a zero value. This is typically done when the kernel will
2123 % be used in special hardware (GPU) convolution processors, to simply
2126 % The format of the ZeroKernelNans method is:
2128 % voidZeroKernelNans (KernelInfo *kernel)
2130 % A description of each parameter follows:
2132 % o kernel: the Morphology/Convolution kernel
2134 % FUTURE: return the information in a string for API usage.
2136 MagickExport void ZeroKernelNans(KernelInfo *kernel)
2141 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2142 if ( IsNan(kernel->values[i]) )
2143 kernel->values[i] = 0.0;