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 kernels, 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 kernel 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 freed using the DestroyKernelInfo() when you
146 % are finished with it. Do not free this memory yourself.
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 kernel 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 % You can define a 'list of kernels' which can be used by some morphology
173 % operators A list is defined as a semi-colon seperated list kernels.
175 % " kernel ; kernel ; kernel "
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 /* This was separated so that it could be used as a separate
189 ** array input handling function, such as for -color-matrix
191 static KernelInfo *ParseKernelArray(const char *kernel_string)
197 token[MaxTextExtent];
207 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
209 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
210 if (kernel == (KernelInfo *)NULL)
212 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
213 kernel->minimum = kernel->maximum = 0.0;
214 kernel->negative_range = kernel->positive_range = 0.0;
215 kernel->type = UserDefinedKernel;
216 kernel->next = (KernelInfo *) NULL;
217 kernel->signature = MagickSignature;
219 /* find end of this specific kernel definition string */
220 end = strchr(kernel_string, ';');
221 if ( end == (char *) NULL )
222 end = strchr(kernel_string, '\0');
224 /* Has a ':' in argument - New user kernel specification */
225 p = strchr(kernel_string, ':');
226 if ( p != (char *) NULL && p < end)
234 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
235 memcpy(token, kernel_string, (size_t) (p-kernel_string));
236 token[p-kernel_string] = '\0';
237 SetGeometryInfo(&args);
238 flags = ParseGeometry(token, &args);
240 /* Size handling and checks of geometry settings */
241 if ( (flags & WidthValue) == 0 ) /* if no width then */
242 args.rho = args.sigma; /* then width = height */
243 if ( args.rho < 1.0 ) /* if width too small */
244 args.rho = 1.0; /* then width = 1 */
245 if ( args.sigma < 1.0 ) /* if height too small */
246 args.sigma = args.rho; /* then height = width */
247 kernel->width = (unsigned long)args.rho;
248 kernel->height = (unsigned long)args.sigma;
250 /* Offset Handling and Checks */
251 if ( args.xi < 0.0 || args.psi < 0.0 )
252 return(DestroyKernelInfo(kernel));
253 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
254 : (long) (kernel->width-1)/2;
255 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
256 : (long) (kernel->height-1)/2;
257 if ( kernel->x >= (long) kernel->width ||
258 kernel->y >= (long) kernel->height )
259 return(DestroyKernelInfo(kernel));
261 p++; /* advance beyond the ':' */
264 { /* ELSE - Old old specification, forming odd-square kernel */
265 /* count up number of values given */
266 p=(const char *) kernel_string;
267 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
268 p++; /* ignore "'" chars for convolve filter usage - Cristy */
269 for (i=0; p < end; i++)
271 GetMagickToken(p,&p,token);
273 GetMagickToken(p,&p,token);
275 /* set the size of the kernel - old sized square */
276 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
277 kernel->x = kernel->y = (long) (kernel->width-1)/2;
278 p=(const char *) kernel_string;
279 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
280 p++; /* ignore "'" chars for convolve filter usage - Cristy */
283 /* Read in the kernel values from rest of input string argument */
284 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
285 kernel->height*sizeof(double));
286 if (kernel->values == (double *) NULL)
287 return(DestroyKernelInfo(kernel));
289 kernel->minimum = +MagickHuge;
290 kernel->maximum = -MagickHuge;
291 kernel->negative_range = kernel->positive_range = 0.0;
293 for (i=0; (i < (long) (kernel->width*kernel->height)) && (p < end); i++)
295 GetMagickToken(p,&p,token);
297 GetMagickToken(p,&p,token);
298 if ( LocaleCompare("nan",token) == 0
299 || LocaleCompare("-",token) == 0 ) {
300 kernel->values[i] = nan; /* do not include this value in kernel */
303 kernel->values[i] = StringToDouble(token);
304 ( kernel->values[i] < 0)
305 ? ( kernel->negative_range += kernel->values[i] )
306 : ( kernel->positive_range += kernel->values[i] );
307 Minimize(kernel->minimum, kernel->values[i]);
308 Maximize(kernel->maximum, kernel->values[i]);
312 /* sanity check -- no more values in kernel definition */
313 GetMagickToken(p,&p,token);
314 if ( *token != '\0' && *token != ';' && *token != '\'' )
315 return(DestroyKernelInfo(kernel));
318 /* this was the old method of handling a incomplete kernel */
319 if ( i < (long) (kernel->width*kernel->height) ) {
320 Minimize(kernel->minimum, kernel->values[i]);
321 Maximize(kernel->maximum, kernel->values[i]);
322 for ( ; i < (long) (kernel->width*kernel->height); i++)
323 kernel->values[i]=0.0;
326 /* Number of values for kernel was not enough - Report Error */
327 if ( i < (long) (kernel->width*kernel->height) )
328 return(DestroyKernelInfo(kernel));
331 /* check that we recieved at least one real (non-nan) value! */
332 if ( kernel->minimum == MagickHuge )
333 return(DestroyKernelInfo(kernel));
338 static KernelInfo *ParseNamedKernel(const char *kernel_string)
341 token[MaxTextExtent];
356 /* Parse special 'named' kernel */
357 GetMagickToken(kernel_string,&p,token);
358 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
359 if ( type < 0 || type == UserDefinedKernel )
360 return((KernelInfo *)NULL); /* not a valid named kernel */
362 while (((isspace((int) ((unsigned char) *p)) != 0) ||
363 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
366 end = strchr(p, ';'); /* end of this kernel defintion */
367 if ( end == (char *) NULL )
368 end = strchr(p, '\0');
370 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
371 memcpy(token, p, (size_t) (end-p));
373 SetGeometryInfo(&args);
374 flags = ParseGeometry(token, &args);
376 /* special handling of missing values in input string */
378 case RectangleKernel:
379 if ( (flags & WidthValue) == 0 ) /* if no width then */
380 args.rho = args.sigma; /* then width = height */
381 if ( args.rho < 1.0 ) /* if width too small */
382 args.rho = 3; /* then width = 3 */
383 if ( args.sigma < 1.0 ) /* if height too small */
384 args.sigma = args.rho; /* then height = width */
385 if ( (flags & XValue) == 0 ) /* center offset if not defined */
386 args.xi = (double)(((long)args.rho-1)/2);
387 if ( (flags & YValue) == 0 )
388 args.psi = (double)(((long)args.sigma-1)/2);
394 /* If no scale given (a 0 scale is valid! - set it to 1.0 */
395 if ( (flags & HeightValue) == 0 )
398 case ChebyshevKernel:
399 case ManhattenKernel:
400 case EuclideanKernel:
401 if ( (flags & HeightValue) == 0 )
402 args.sigma = 100.0; /* default distance scaling */
403 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
404 args.sigma = QuantumRange/args.sigma; /* maximum pixel distance */
405 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
406 args.sigma *= QuantumRange/100.0; /* percentage of color range */
412 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
415 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
422 token[MaxTextExtent];
427 /* If it does not start with an alpha - Its is a user defined kernel array */
428 GetMagickToken(kernel_string,NULL,token);
429 if (isalpha((int) ((unsigned char) *token)) == 0)
430 kernel = ParseKernelArray(kernel_string);
432 kernel = ParseNamedKernel(kernel_string);
434 /* was this the last (or only) kernel in the string */
435 next = strchr(kernel_string, ';');
436 if ( next == (char *) NULL )
439 /* Add the next kernel to the list */
440 kernel->next = AcquireKernelInfo( next+1 );
442 /* failed for some reason? */
443 if ( kernel->next == (KernelInfo *) NULL )
444 return(DestroyKernelInfo(kernel));
451 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
455 % A c q u i r e K e r n e l B u i l t I n %
459 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
461 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
462 % kernels used for special purposes such as gaussian blurring, skeleton
463 % pruning, and edge distance determination.
465 % They take a KernelType, and a set of geometry style arguments, which were
466 % typically decoded from a user supplied string, or from a more complex
467 % Morphology Method that was requested.
469 % The format of the AcquireKernalBuiltIn method is:
471 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
472 % const GeometryInfo args)
474 % A description of each parameter follows:
476 % o type: the pre-defined type of kernel wanted
478 % o args: arguments defining or modifying the kernel
480 % Convolution Kernels
482 % Gaussian "{radius},{sigma}"
483 % Generate a two-dimentional gaussian kernel, as used by -gaussian
484 % A sigma is required, (with the 'x'), due to historical reasons.
486 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
487 % the final size of the resulting kernel to a square 2*radius+1 in size.
488 % The radius should be at least 2 times that of the sigma value, or
489 % sever clipping and aliasing may result. If not given or set to 0 the
490 % radius will be determined so as to produce the best minimal error
491 % result, which is usally much larger than is normally needed.
493 % Blur "{radius},{sigma},{angle}"
494 % As per Gaussian, but generates a 1 dimensional or linear gaussian
495 % blur, at the angle given (current restricted to orthogonal angles).
496 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
498 % NOTE that two such blurs perpendicular to each other is equivelent to
499 % -blur and the previous gaussian, but is often 10 or more times faster.
501 % Comet "{width},{sigma},{angle}"
502 % Blur in one direction only, mush like how a bright object leaves
503 % a comet like trail. The Kernel is actually half a gaussian curve,
504 % Adding two such blurs in oppiste directions produces a Linear Blur.
506 % NOTE: that the first argument is the width of the kernel and not the
507 % radius of the kernel.
509 % # Still to be implemented...
511 % # Sharpen "{radius},{sigma}
512 % # Negated Gaussian (center zeroed and re-normalized),
513 % # with a 2 unit positive peak. -- Check On line documentation
515 % # Laplacian "{radius},{sigma}"
516 % # Laplacian (a mexican hat like) Function
518 % # LOG "{radius},{sigma1},{sigma2}
519 % # Laplacian of Gaussian
521 % # DOG "{radius},{sigma1},{sigma2}
522 % # Difference of two Gaussians
526 % # Set kernel values using a resize filter, and given scale (sigma)
527 % # Cylindrical or Linear. Is this posible with an image?
532 % Rectangle "{geometry}"
533 % Simply generate a rectangle of 1's with the size given. You can also
534 % specify the location of the 'control point', otherwise the closest
535 % pixel to the center of the rectangle is selected.
537 % Properly centered and odd sized rectangles work the best.
539 % Diamond "[{radius}[,{scale}]]"
540 % Generate a diamond shaped kernel with given radius to the points.
541 % Kernel size will again be radius*2+1 square and defaults to radius 1,
542 % generating a 3x3 kernel that is slightly larger than a square.
544 % Square "[{radius}[,{scale}]]"
545 % Generate a square shaped kernel of size radius*2+1, and defaulting
546 % to a 3x3 (radius 1).
548 % Note that using a larger radius for the "Square" or the "Diamond"
549 % is also equivelent to iterating the basic morphological method
550 % that many times. However However iterating with the smaller radius 1
551 % default is actually faster than using a larger kernel radius.
553 % Disk "[{radius}[,{scale}]]
554 % Generate a binary disk of the radius given, radius may be a float.
555 % Kernel size will be ceil(radius)*2+1 square.
556 % NOTE: Here are some disk shapes of specific interest
557 % "disk:1" => "diamond" or "cross:1"
558 % "disk:1.5" => "square"
559 % "disk:2" => "diamond:2"
560 % "disk:2.5" => a general disk shape of radius 2
561 % "disk:2.9" => "square:2"
562 % "disk:3.5" => default - octagonal/disk shape of radius 3
563 % "disk:4.2" => roughly octagonal shape of radius 4
564 % "disk:4.3" => a general disk shape of radius 4
565 % After this all the kernel shape becomes more and more circular.
567 % Because a "disk" is more circular when using a larger radius, using a
568 % larger radius is preferred over iterating the morphological operation.
570 % Plus "[{radius}[,{scale}]]"
571 % Generate a kernel in the shape of a 'plus' sign. The length of each
572 % arm is also the radius, which defaults to 2.
574 % This kernel is not a good general morphological kernel, but is used
575 % more for highlighting and marking any single pixels in an image using,
576 % a "Dilate" or "Erode" method as appropriate.
578 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
580 % Note that unlike other kernels iterating a plus does not produce the
581 % same result as using a larger radius for the cross.
583 % Distance Measuring Kernels
585 % Chebyshev "[{radius}][x{scale}[%!]]"
586 % Manhatten "[{radius}][x{scale}[%!]]"
587 % Euclidean "[{radius}][x{scale}[%!]]"
589 % Different types of distance measuring methods, which are used with the
590 % a 'Distance' morphology method for generating a gradient based on
591 % distance from an edge of a binary shape, though there is a technique
592 % for handling a anti-aliased shape.
594 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
595 % one to any neighbour, orthogonal or diagonal. One why of thinking of
596 % it is the number of squares a 'King' or 'Queen' in chess needs to
597 % traverse reach any other position on a chess board. It results in a
598 % 'square' like distance function, but one where diagonals are closer
601 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
602 % Cab metric), is the distance needed when you can only travel in
603 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
604 % in chess would travel. It results in a diamond like distances, where
605 % diagonals are further than expected.
607 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
608 % However by default the kernel size only has a radius of 1, which
609 % limits the distance to 'Knight' like moves, with only orthogonal and
610 % diagonal measurements being correct. As such for the default kernel
611 % you will get octagonal like distance function, which is reasonally
614 % However if you use a larger radius such as "Euclidean:4" you will
615 % get a much smoother distance gradient from the edge of the shape.
616 % Of course a larger kernel is slower to use, and generally not needed.
618 % To allow the use of fractional distances that you get with diagonals
619 % the actual distance is scaled by a fixed value which the user can
620 % provide. This is not actually nessary for either ""Chebyshev" or
621 % "Manhatten" distance kernels, but is done for all three distance
622 % kernels. If no scale is provided it is set to a value of 100,
623 % allowing for a maximum distance measurement of 655 pixels using a Q16
624 % version of IM, from any edge. However for small images this can
625 % result in quite a dark gradient.
627 % See the 'Distance' Morphological Method, for information of how it is
630 % # Hit-n-Miss Kernel-Lists -- Still to be implemented
632 % # specifically for Pruning, Thinning, Thickening
636 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
637 const GeometryInfo *args)
650 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
652 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
653 if (kernel == (KernelInfo *) NULL)
655 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
656 kernel->minimum = kernel->maximum = 0.0;
657 kernel->negative_range = kernel->positive_range = 0.0;
659 kernel->next = (KernelInfo *) NULL;
660 kernel->signature = MagickSignature;
663 /* Convolution Kernels */
666 sigma = fabs(args->sigma);
668 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
670 kernel->width = kernel->height =
671 GetOptimalKernelWidth2D(args->rho,sigma);
672 kernel->x = kernel->y = (long) (kernel->width-1)/2;
673 kernel->negative_range = kernel->positive_range = 0.0;
674 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
675 kernel->height*sizeof(double));
676 if (kernel->values == (double *) NULL)
677 return(DestroyKernelInfo(kernel));
679 sigma = 2.0*sigma*sigma; /* simplify the expression */
680 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
681 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
682 kernel->positive_range += (
684 exp(-((double)(u*u+v*v))/sigma)
685 /* / (MagickPI*sigma) */ );
687 kernel->maximum = kernel->values[
688 kernel->y*kernel->width+kernel->x ];
690 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
696 sigma = fabs(args->sigma);
698 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
700 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
701 kernel->x = (long) (kernel->width-1)/2;
704 kernel->negative_range = kernel->positive_range = 0.0;
705 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
706 kernel->height*sizeof(double));
707 if (kernel->values == (double *) NULL)
708 return(DestroyKernelInfo(kernel));
712 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
713 ** It generates a gaussian 3 times the width, and compresses it into
714 ** the expected range. This produces a closer normalization of the
715 ** resulting kernel, especially for very low sigma values.
716 ** As such while wierd it is prefered.
718 ** I am told this method originally came from Photoshop.
720 sigma *= KernelRank; /* simplify expanded curve */
721 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
722 (void) ResetMagickMemory(kernel->values,0, (size_t)
723 kernel->width*sizeof(double));
724 for ( u=-v; u <= v; u++) {
725 kernel->values[(u+v)/KernelRank] +=
726 exp(-((double)(u*u))/(2.0*sigma*sigma))
727 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
729 for (i=0; i < (long) kernel->width; i++)
730 kernel->positive_range += kernel->values[i];
732 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
733 kernel->positive_range += (
735 exp(-((double)(u*u))/(2.0*sigma*sigma))
736 /* / (MagickSQ2PI*sigma) */ );
739 kernel->maximum = kernel->values[ kernel->x ];
740 /* Note that neither methods above generate a normalized kernel,
741 ** though it gets close. The kernel may be 'clipped' by a user defined
742 ** radius, producing a smaller (darker) kernel. Also for very small
743 ** sigma's (> 0.1) the central value becomes larger than one, and thus
744 ** producing a very bright kernel.
747 /* Normalize the 1D Gaussian Kernel
749 ** Because of this the divisor in the above kernel generator is
750 ** not needed, so is not done above.
752 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
754 /* rotate the kernel by given angle */
755 RotateKernelInfo(kernel, args->xi);
760 sigma = fabs(args->sigma);
762 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
764 if ( args->rho < 1.0 )
765 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
767 kernel->width = (unsigned long)args->rho;
768 kernel->x = kernel->y = 0;
770 kernel->negative_range = kernel->positive_range = 0.0;
771 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
772 kernel->height*sizeof(double));
773 if (kernel->values == (double *) NULL)
774 return(DestroyKernelInfo(kernel));
776 /* A comet blur is half a gaussian curve, so that the object is
777 ** blurred in one direction only. This may not be quite the right
778 ** curve so may change in the future. The function must be normalised.
782 sigma *= KernelRank; /* simplify expanded curve */
783 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
784 (void) ResetMagickMemory(kernel->values,0, (size_t)
785 kernel->width*sizeof(double));
786 for ( u=0; u < v; u++) {
787 kernel->values[u/KernelRank] +=
788 exp(-((double)(u*u))/(2.0*sigma*sigma))
789 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
791 for (i=0; i < (long) kernel->width; i++)
792 kernel->positive_range += kernel->values[i];
794 for ( i=0; i < (long) kernel->width; i++)
795 kernel->positive_range += (
797 exp(-((double)(i*i))/(2.0*sigma*sigma))
798 /* / (MagickSQ2PI*sigma) */ );
801 kernel->maximum = kernel->values[0];
803 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
804 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
807 /* Boolean Kernels */
808 case RectangleKernel:
812 if ( type == SquareKernel )
815 kernel->width = kernel->height = 3; /* default radius = 1 */
817 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
818 kernel->x = kernel->y = (long) (kernel->width-1)/2;
822 /* NOTE: user defaults set in "AcquireKernelInfo()" */
823 if ( args->rho < 1.0 || args->sigma < 1.0 )
824 return(DestroyKernelInfo(kernel)); /* invalid args given */
825 kernel->width = (unsigned long)args->rho;
826 kernel->height = (unsigned long)args->sigma;
827 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
828 args->psi < 0.0 || args->psi > (double)kernel->height )
829 return(DestroyKernelInfo(kernel)); /* invalid args given */
830 kernel->x = (long) args->xi;
831 kernel->y = (long) args->psi;
834 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
835 kernel->height*sizeof(double));
836 if (kernel->values == (double *) NULL)
837 return(DestroyKernelInfo(kernel));
839 /* set all kernel values to 1.0 */
840 u=(long) kernel->width*kernel->height;
841 for ( i=0; i < u; i++)
842 kernel->values[i] = scale;
843 kernel->minimum = kernel->maximum = scale; /* a flat shape */
844 kernel->positive_range = scale*u;
850 kernel->width = kernel->height = 3; /* default radius = 1 */
852 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
853 kernel->x = kernel->y = (long) (kernel->width-1)/2;
855 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
856 kernel->height*sizeof(double));
857 if (kernel->values == (double *) NULL)
858 return(DestroyKernelInfo(kernel));
860 /* set all kernel values within diamond area to scale given */
861 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
862 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
863 if ((labs(u)+labs(v)) <= (long)kernel->x)
864 kernel->positive_range += kernel->values[i] = args->sigma;
866 kernel->values[i] = nan;
867 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
875 limit = (long)(args->rho*args->rho);
876 if (args->rho < 0.1) /* default radius approx 3.5 */
877 kernel->width = kernel->height = 7L, limit = 10L;
879 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
880 kernel->x = kernel->y = (long) (kernel->width-1)/2;
882 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
883 kernel->height*sizeof(double));
884 if (kernel->values == (double *) NULL)
885 return(DestroyKernelInfo(kernel));
887 /* set all kernel values within disk area to 1.0 */
888 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
889 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
890 if ((u*u+v*v) <= limit)
891 kernel->positive_range += kernel->values[i] = args->sigma;
893 kernel->values[i] = nan;
894 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
900 kernel->width = kernel->height = 5; /* default radius 2 */
902 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
903 kernel->x = kernel->y = (long) (kernel->width-1)/2;
905 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
906 kernel->height*sizeof(double));
907 if (kernel->values == (double *) NULL)
908 return(DestroyKernelInfo(kernel));
910 /* set all kernel values along axises to 1.0 */
911 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
912 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
913 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
914 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
915 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
918 /* Distance Measuring Kernels */
919 case ChebyshevKernel:
922 kernel->width = kernel->height = 3; /* default radius = 1 */
924 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
925 kernel->x = kernel->y = (long) (kernel->width-1)/2;
927 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
928 kernel->height*sizeof(double));
929 if (kernel->values == (double *) NULL)
930 return(DestroyKernelInfo(kernel));
932 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
933 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
934 kernel->positive_range += ( kernel->values[i] =
935 args->sigma*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
936 kernel->maximum = kernel->values[0];
939 case ManhattenKernel:
942 kernel->width = kernel->height = 3; /* default radius = 1 */
944 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
945 kernel->x = kernel->y = (long) (kernel->width-1)/2;
947 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
948 kernel->height*sizeof(double));
949 if (kernel->values == (double *) NULL)
950 return(DestroyKernelInfo(kernel));
952 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
953 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
954 kernel->positive_range += ( kernel->values[i] =
955 args->sigma*(labs(u)+labs(v)) );
956 kernel->maximum = kernel->values[0];
959 case EuclideanKernel:
962 kernel->width = kernel->height = 3; /* default radius = 1 */
964 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
965 kernel->x = kernel->y = (long) (kernel->width-1)/2;
967 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
968 kernel->height*sizeof(double));
969 if (kernel->values == (double *) NULL)
970 return(DestroyKernelInfo(kernel));
972 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
973 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
974 kernel->positive_range += ( kernel->values[i] =
975 args->sigma*sqrt((double)(u*u+v*v)) );
976 kernel->maximum = kernel->values[0];
979 /* Undefined Kernels */
980 case LaplacianKernel:
983 perror("Kernel Type has not been defined yet");
986 /* Generate a No-Op minimal kernel - 1x1 pixel */
987 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
988 if (kernel->values == (double *) NULL)
989 return(DestroyKernelInfo(kernel));
990 kernel->width = kernel->height = 1;
991 kernel->x = kernel->x = 0;
992 kernel->type = UndefinedKernel;
994 kernel->positive_range =
995 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
1003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1007 % C l o n e K e r n e l I n f o %
1011 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1013 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
1014 % can be modified without effecting the original. The cloned kernel should
1015 % be destroyed using DestoryKernelInfo() when no longer needed.
1017 % The format of the CloneKernelInfo method is:
1019 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1021 % A description of each parameter follows:
1023 % o kernel: the Morphology/Convolution kernel to be cloned
1026 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1034 assert(kernel != (KernelInfo *) NULL);
1035 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1036 if (new_kernel == (KernelInfo *) NULL)
1038 *new_kernel=(*kernel); /* copy values in structure */
1040 /* replace the values with a copy of the values */
1041 new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1042 kernel->height*sizeof(double));
1043 if (new_kernel->values == (double *) NULL)
1044 return(DestroyKernelInfo(new_kernel));
1045 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1046 new_kernel->values[i]=kernel->values[i];
1048 /* Also clone the next kernel in the kernel list */
1049 if ( kernel->next != (KernelInfo *) NULL ) {
1050 new_kernel->next = CloneKernelInfo(kernel->next);
1051 if ( new_kernel->next == (KernelInfo *) NULL )
1052 return(DestroyKernelInfo(new_kernel));
1059 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1063 % D e s t r o y K e r n e l I n f o %
1067 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1069 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
1072 % The format of the DestroyKernelInfo method is:
1074 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1076 % A description of each parameter follows:
1078 % o kernel: the Morphology/Convolution kernel to be destroyed
1082 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1084 assert(kernel != (KernelInfo *) NULL);
1086 if ( kernel->next != (KernelInfo *) NULL )
1087 kernel->next = DestroyKernelInfo(kernel->next);
1089 kernel->values = (double *)RelinquishMagickMemory(kernel->values);
1090 kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
1095 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1099 % M o r p h o l o g y I m a g e C h a n n e l %
1103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1105 % MorphologyImageChannel() applies a user supplied kernel to the image
1106 % according to the given mophology method.
1108 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1109 % by the kernel generator.
1111 % The format of the MorphologyImage method is:
1113 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
1114 % const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
1115 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1116 % channel,MorphologyMethod method,const long iterations,
1117 % KernelInfo *kernel,ExceptionInfo *exception)
1119 % A description of each parameter follows:
1121 % o image: the image.
1123 % o method: the morphology method to be applied.
1125 % o iterations: apply the operation this many times (or no change).
1126 % A value of -1 means loop until no change found.
1127 % How this is applied may depend on the morphology method.
1128 % Typically this is a value of 1.
1130 % o channel: the channel type.
1132 % o kernel: An array of double representing the morphology kernel.
1133 % Warning: kernel may be normalized for the Convolve method.
1135 % o exception: return any errors or warnings in this structure.
1138 % TODO: bias and auto-scale handling of the kernel for convolution
1139 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1140 % by the kernel generator.
1145 /* Internal function
1146 * Apply the low-level Morphology Method Primatives using the given Kernel
1147 * Returning the number of pixels that changed.
1148 * Two pre-created images must be provided, no image is created.
1150 static unsigned long MorphologyPrimative(const Image *image, Image
1151 *result_image, const MorphologyMethod method, const ChannelType channel,
1152 const KernelInfo *kernel, ExceptionInfo *exception)
1154 #define MorphologyTag "Morphology/Image"
1171 /* Only the most basic morphology is actually performed by this routine */
1174 Apply Basic Morphology to Image.
1180 GetMagickPixelPacket(image,&bias);
1181 SetMagickPixelPacketBias(image,&bias);
1182 /* Future: handle auto-bias from user, based on kernel input */
1184 p_view=AcquireCacheView(image);
1185 q_view=AcquireCacheView(result_image);
1187 /* Some methods (including convolve) needs use a reflected kernel.
1188 * Adjust 'origin' offsets for this reflected kernel.
1193 case ConvolveMorphology:
1194 case DilateMorphology:
1195 case DilateIntensityMorphology:
1196 case DistanceMorphology:
1197 /* kernel needs to used with reflection about origin */
1198 offx = (long) kernel->width-offx-1;
1199 offy = (long) kernel->height-offy-1;
1201 case ErodeMorphology:
1202 case ErodeIntensityMorphology:
1203 case HitAndMissMorphology:
1204 case ThinningMorphology:
1205 case ThickenMorphology:
1206 /* kernel is user as is, without reflection */
1209 perror("Not a low level Morphology Method");
1213 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1214 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1216 for (y=0; y < (long) image->rows; y++)
1221 register const PixelPacket
1224 register const IndexPacket
1225 *restrict p_indexes;
1227 register PixelPacket
1230 register IndexPacket
1231 *restrict q_indexes;
1239 if (status == MagickFalse)
1241 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1242 image->columns+kernel->width, kernel->height, exception);
1243 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1245 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1250 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1251 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1252 r = (image->columns+kernel->width)*offy+offx; /* constant */
1254 for (x=0; x < (long) image->columns; x++)
1262 register const double
1265 register const PixelPacket
1268 register const IndexPacket
1269 *restrict k_indexes;
1276 /* Copy input to ouput image for unused channels
1277 * This removes need for 'cloning' a new image every iteration
1280 if (image->colorspace == CMYKColorspace)
1281 q_indexes[x] = p_indexes[r];
1288 min.index = (MagickRealType) QuantumRange;
1293 max.index = (MagickRealType) 0;
1294 /* original pixel value */
1295 result.red = (MagickRealType) p[r].red;
1296 result.green = (MagickRealType) p[r].green;
1297 result.blue = (MagickRealType) p[r].blue;
1298 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1299 if ( image->colorspace == CMYKColorspace)
1300 result.index = (MagickRealType) p_indexes[r];
1303 case ConvolveMorphology:
1304 /* Set the user defined bias of the weighted average output
1306 ** FUTURE: provide some way for internal functions to disable
1307 ** user provided bias and scaling effects.
1311 case DilateIntensityMorphology:
1312 case ErodeIntensityMorphology:
1313 result.red = 0.0; /* flag indicating when first match found */
1320 case ConvolveMorphology:
1321 /* Weighted Average of pixels using reflected kernel
1323 ** NOTE for correct working of this operation for asymetrical
1324 ** kernels, the kernel needs to be applied in its reflected form.
1325 ** That is its values needs to be reversed.
1327 ** Correlation is actually the same as this but without reflecting
1328 ** the kernel, and thus 'lower-level' that Convolution. However
1329 ** as Convolution is the more common method used, and it does not
1330 ** really cost us much in terms of processing to use a reflected
1331 ** kernel, so it is Convolution that is implemented.
1333 ** Correlation will have its kernel reflected before calling
1334 ** this function to do a Convolve.
1336 ** For more details of Correlation vs Convolution see
1337 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
1339 if (((channel & SyncChannels) != 0 ) &&
1340 (image->matte == MagickTrue))
1341 { /* Channel has a 'Sync' Flag, and Alpha Channel enabled.
1342 ** Weight the color channels with Alpha Channel so that
1343 ** transparent pixels are not part of the results.
1346 alpha, /* color channel weighting : kernel*alpha */
1347 gamma; /* divisor, sum of weighting values */
1350 k = &kernel->values[ kernel->width*kernel->height-1 ];
1352 k_indexes = p_indexes;
1353 for (v=0; v < (long) kernel->height; v++) {
1354 for (u=0; u < (long) kernel->width; u++, k--) {
1355 if ( IsNan(*k) ) continue;
1356 alpha=(*k)*(QuantumScale*(QuantumRange-
1357 k_pixels[u].opacity));
1359 result.red += alpha*k_pixels[u].red;
1360 result.green += alpha*k_pixels[u].green;
1361 result.blue += alpha*k_pixels[u].blue;
1362 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1363 if ( image->colorspace == CMYKColorspace)
1364 result.index += alpha*k_indexes[u];
1366 k_pixels += image->columns+kernel->width;
1367 k_indexes += image->columns+kernel->width;
1369 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1370 result.red *= gamma;
1371 result.green *= gamma;
1372 result.blue *= gamma;
1373 result.opacity *= gamma;
1374 result.index *= gamma;
1378 /* No 'Sync' flag, or no Alpha involved.
1379 ** Convolution is simple individual channel weigthed sum.
1381 k = &kernel->values[ kernel->width*kernel->height-1 ];
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) ) continue;
1387 result.red += (*k)*k_pixels[u].red;
1388 result.green += (*k)*k_pixels[u].green;
1389 result.blue += (*k)*k_pixels[u].blue;
1390 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1391 if ( image->colorspace == CMYKColorspace)
1392 result.index += (*k)*k_indexes[u];
1394 k_pixels += image->columns+kernel->width;
1395 k_indexes += image->columns+kernel->width;
1400 case ErodeMorphology:
1401 /* Minimum Value within kernel neighbourhood
1403 ** NOTE that the kernel is not reflected for this operation!
1405 ** NOTE: in normal Greyscale Morphology, the kernel value should
1406 ** be added to the real value, this is currently not done, due to
1407 ** the nature of the boolean kernels being used.
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;
1415 Minimize(min.red, (double) k_pixels[u].red);
1416 Minimize(min.green, (double) k_pixels[u].green);
1417 Minimize(min.blue, (double) k_pixels[u].blue);
1418 Minimize(min.opacity,
1419 QuantumRange-(double) k_pixels[u].opacity);
1420 if ( image->colorspace == CMYKColorspace)
1421 Minimize(min.index, (double) k_indexes[u]);
1423 k_pixels += image->columns+kernel->width;
1424 k_indexes += image->columns+kernel->width;
1429 case DilateMorphology:
1430 /* Maximum Value within kernel neighbourhood
1432 ** NOTE for correct working of this operation for asymetrical
1433 ** kernels, the kernel needs to be applied in its reflected form.
1434 ** That is its values needs to be reversed.
1436 ** NOTE: in normal Greyscale Morphology, the kernel value should
1437 ** be added to the real value, this is currently not done, due to
1438 ** the nature of the boolean kernels being used.
1441 k = &kernel->values[ kernel->width*kernel->height-1 ];
1443 k_indexes = p_indexes;
1444 for (v=0; v < (long) kernel->height; v++) {
1445 for (u=0; u < (long) kernel->width; u++, k--) {
1446 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1447 Maximize(max.red, (double) k_pixels[u].red);
1448 Maximize(max.green, (double) k_pixels[u].green);
1449 Maximize(max.blue, (double) k_pixels[u].blue);
1450 Maximize(max.opacity,
1451 QuantumRange-(double) k_pixels[u].opacity);
1452 if ( image->colorspace == CMYKColorspace)
1453 Maximize(max.index, (double) k_indexes[u]);
1455 k_pixels += image->columns+kernel->width;
1456 k_indexes += image->columns+kernel->width;
1460 case HitAndMissMorphology:
1461 case ThinningMorphology:
1462 case ThickenMorphology:
1463 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
1465 ** NOTE that the kernel is not reflected for this operation,
1466 ** and consists of both foreground and background pixel
1467 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
1468 ** with either Nan or 0.5 values for don't care.
1470 ** Note that this can produce negative results, though really
1471 ** only a positive match has any real value.
1475 k_indexes = p_indexes;
1476 for (v=0; v < (long) kernel->height; v++) {
1477 for (u=0; u < (long) kernel->width; u++, k++) {
1478 if ( IsNan(*k) ) continue;
1480 { /* minimim of foreground pixels */
1481 Minimize(min.red, (double) k_pixels[u].red);
1482 Minimize(min.green, (double) k_pixels[u].green);
1483 Minimize(min.blue, (double) k_pixels[u].blue);
1484 Minimize(min.opacity,
1485 QuantumRange-(double) k_pixels[u].opacity);
1486 if ( image->colorspace == CMYKColorspace)
1487 Minimize(min.index, (double) k_indexes[u]);
1489 else if ( (*k) < 0.3 )
1490 { /* maximum of background pixels */
1491 Maximize(max.red, (double) k_pixels[u].red);
1492 Maximize(max.green, (double) k_pixels[u].green);
1493 Maximize(max.blue, (double) k_pixels[u].blue);
1494 Maximize(max.opacity,
1495 QuantumRange-(double) k_pixels[u].opacity);
1496 if ( image->colorspace == CMYKColorspace)
1497 Maximize(max.index, (double) k_indexes[u]);
1500 k_pixels += image->columns+kernel->width;
1501 k_indexes += image->columns+kernel->width;
1503 /* Pattern Match only if min fg larger than min bg pixels */
1504 min.red -= max.red; Maximize( min.red, 0.0 );
1505 min.green -= max.green; Maximize( min.green, 0.0 );
1506 min.blue -= max.blue; Maximize( min.blue, 0.0 );
1507 min.opacity -= max.opacity; Maximize( min.opacity, 0.0 );
1508 min.index -= max.index; Maximize( min.index, 0.0 );
1511 case ErodeIntensityMorphology:
1512 /* Select Pixel with Minimum Intensity within kernel neighbourhood
1514 ** WARNING: the intensity test fails for CMYK and does not
1515 ** take into account the moderating effect of teh alpha channel
1516 ** on the intensity.
1518 ** NOTE that the kernel is not reflected for this operation!
1522 k_indexes = p_indexes;
1523 for (v=0; v < (long) kernel->height; v++) {
1524 for (u=0; u < (long) kernel->width; u++, k++) {
1525 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1526 if ( result.red == 0.0 ||
1527 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1528 /* copy the whole pixel - no channel selection */
1530 if ( result.red > 0.0 ) changed++;
1534 k_pixels += image->columns+kernel->width;
1535 k_indexes += image->columns+kernel->width;
1539 case DilateIntensityMorphology:
1540 /* Select Pixel with Maximum Intensity within kernel neighbourhood
1542 ** WARNING: the intensity test fails for CMYK and does not
1543 ** take into account the moderating effect of teh alpha channel
1544 ** on the intensity.
1546 ** NOTE for correct working of this operation for asymetrical
1547 ** kernels, the kernel needs to be applied in its reflected form.
1548 ** That is its values needs to be reversed.
1550 k = &kernel->values[ kernel->width*kernel->height-1 ];
1552 k_indexes = p_indexes;
1553 for (v=0; v < (long) kernel->height; v++) {
1554 for (u=0; u < (long) kernel->width; u++, k--) {
1555 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1556 if ( result.red == 0.0 ||
1557 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1558 /* copy the whole pixel - no channel selection */
1560 if ( result.red > 0.0 ) changed++;
1564 k_pixels += image->columns+kernel->width;
1565 k_indexes += image->columns+kernel->width;
1570 case DistanceMorphology:
1571 /* Add kernel Value and select the minimum value found.
1572 ** The result is a iterative distance from edge of image shape.
1574 ** All Distance Kernels are symetrical, but that may not always
1575 ** be the case. For example how about a distance from left edges?
1576 ** To work correctly with asymetrical kernels the reflected kernel
1577 ** needs to be applied.
1579 ** Actually this is really a GreyErode with a negative kernel!
1582 k = &kernel->values[ kernel->width*kernel->height-1 ];
1584 k_indexes = p_indexes;
1585 for (v=0; v < (long) kernel->height; v++) {
1586 for (u=0; u < (long) kernel->width; u++, k--) {
1587 if ( IsNan(*k) ) continue;
1588 Minimize(result.red, (*k)+k_pixels[u].red);
1589 Minimize(result.green, (*k)+k_pixels[u].green);
1590 Minimize(result.blue, (*k)+k_pixels[u].blue);
1591 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1592 if ( image->colorspace == CMYKColorspace)
1593 Minimize(result.index, (*k)+k_indexes[u]);
1595 k_pixels += image->columns+kernel->width;
1596 k_indexes += image->columns+kernel->width;
1600 case UndefinedMorphology:
1602 break; /* Do nothing */
1604 /* Final mathematics of results (combine with original image?)
1606 ** NOTE: Difference Morphology operators Edge* and *Hat could also
1607 ** be done here but works better with iteration as a image difference
1608 ** in the controling function (below). Thicken and Thinning however
1609 ** should be done here so thay can be iterated correctly.
1612 case HitAndMissMorphology:
1613 case ErodeMorphology:
1614 result = min; /* minimum of neighbourhood */
1616 case DilateMorphology:
1617 result = max; /* maximum of neighbourhood */
1619 case ThinningMorphology:
1620 /* subtract pattern match from original */
1621 result.red -= min.red;
1622 result.green -= min.green;
1623 result.blue -= min.blue;
1624 result.opacity -= min.opacity;
1625 result.index -= min.index;
1627 case ThickenMorphology:
1628 /* Union with original image (maximize) - or should this be + */
1629 Maximize( result.red, min.red );
1630 Maximize( result.green, min.green );
1631 Maximize( result.blue, min.blue );
1632 Maximize( result.opacity, min.opacity );
1633 Maximize( result.index, min.index );
1636 /* result directly calculated or assigned */
1639 /* Assign the resulting pixel values - Clamping Result */
1641 case UndefinedMorphology:
1642 case DilateIntensityMorphology:
1643 case ErodeIntensityMorphology:
1644 break; /* full pixel was directly assigned - not a channel method */
1646 if ((channel & RedChannel) != 0)
1647 q->red = ClampToQuantum(result.red);
1648 if ((channel & GreenChannel) != 0)
1649 q->green = ClampToQuantum(result.green);
1650 if ((channel & BlueChannel) != 0)
1651 q->blue = ClampToQuantum(result.blue);
1652 if ((channel & OpacityChannel) != 0
1653 && image->matte == MagickTrue )
1654 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1655 if ((channel & IndexChannel) != 0
1656 && image->colorspace == CMYKColorspace)
1657 q_indexes[x] = ClampToQuantum(result.index);
1660 /* Count up changed pixels */
1661 if ( ( p[r].red != q->red )
1662 || ( p[r].green != q->green )
1663 || ( p[r].blue != q->blue )
1664 || ( p[r].opacity != q->opacity )
1665 || ( image->colorspace == CMYKColorspace &&
1666 p_indexes[r] != q_indexes[x] ) )
1667 changed++; /* The pixel had some value changed! */
1671 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1672 if (sync == MagickFalse)
1674 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1679 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1680 #pragma omp critical (MagickCore_MorphologyImage)
1682 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1683 if (proceed == MagickFalse)
1687 result_image->type=image->type;
1688 q_view=DestroyCacheView(q_view);
1689 p_view=DestroyCacheView(p_view);
1690 return(status ? (unsigned long) changed : 0);
1694 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1695 method, const long iterations,const KernelInfo *kernel, ExceptionInfo
1701 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1702 iterations,kernel,exception);
1703 return(morphology_image);
1707 MagickExport Image *MorphologyImageChannel(const Image *image,
1708 const ChannelType channel,const MorphologyMethod method,
1709 const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
1730 assert(image != (Image *) NULL);
1731 assert(image->signature == MagickSignature);
1732 assert(kernel != (KernelInfo *) NULL);
1733 assert(kernel->signature == MagickSignature);
1734 assert(exception != (ExceptionInfo *) NULL);
1735 assert(exception->signature == MagickSignature);
1737 if ( iterations == 0 )
1738 return((Image *)NULL); /* null operation - nothing to do! */
1740 /* kernel must be valid at this point
1741 * (except maybe for posible future morphology methods like "Prune"
1743 assert(kernel != (KernelInfo *)NULL);
1745 new_image = (Image *) NULL; /* for cleanup */
1746 old_image = (Image *) NULL;
1747 grad_image = (Image *) NULL;
1748 curr_kernel = (KernelInfo *) NULL;
1749 count = 0; /* interation count */
1750 changed = 1; /* was last run succesfull */
1751 curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
1752 curr_method = method; /* to be changed as nessary */
1754 limit = (unsigned long) iterations;
1755 if ( iterations < 0 )
1756 limit = image->columns > image->rows ? image->columns : image->rows;
1758 /* Third-level morphology methods */
1759 switch( curr_method ) {
1760 case EdgeMorphology:
1761 grad_image = MorphologyImageChannel(image, channel,
1762 DilateMorphology, iterations, curr_kernel, exception);
1763 if ( grad_image == (Image *) NULL )
1766 case EdgeInMorphology:
1767 curr_method = ErodeMorphology;
1769 case EdgeOutMorphology:
1770 curr_method = DilateMorphology;
1772 case TopHatMorphology:
1773 curr_method = OpenMorphology;
1775 case BottomHatMorphology:
1776 curr_method = CloseMorphology;
1779 break; /* not a third-level method */
1782 /* Second-level morphology methods */
1783 switch( curr_method ) {
1784 case OpenMorphology:
1785 /* Open is a Erode then a Dilate without reflection */
1786 new_image = MorphologyImageChannel(image, channel,
1787 ErodeMorphology, iterations, curr_kernel, exception);
1788 if (new_image == (Image *) NULL)
1790 curr_method = DilateMorphology;
1792 case OpenIntensityMorphology:
1793 new_image = MorphologyImageChannel(image, channel,
1794 ErodeIntensityMorphology, iterations, curr_kernel, exception);
1795 if (new_image == (Image *) NULL)
1797 curr_method = DilateIntensityMorphology;
1800 case CloseMorphology:
1801 /* Close is a Dilate then Erode using reflected kernel */
1802 /* A reflected kernel is needed for a Close */
1803 if ( curr_kernel == kernel )
1804 curr_kernel = CloneKernelInfo(kernel);
1805 if (curr_kernel == (KernelInfo *) NULL)
1807 RotateKernelInfo(curr_kernel,180);
1808 new_image = MorphologyImageChannel(image, channel,
1809 DilateMorphology, iterations, curr_kernel, exception);
1810 if (new_image == (Image *) NULL)
1812 curr_method = ErodeMorphology;
1814 case CloseIntensityMorphology:
1815 /* A reflected kernel is needed for a Close */
1816 if ( curr_kernel == kernel )
1817 curr_kernel = CloneKernelInfo(kernel);
1818 if (curr_kernel == (KernelInfo *) NULL)
1820 RotateKernelInfo(curr_kernel,180);
1821 new_image = MorphologyImageChannel(image, channel,
1822 DilateIntensityMorphology, iterations, curr_kernel, exception);
1823 if (new_image == (Image *) NULL)
1825 curr_method = ErodeIntensityMorphology;
1828 case CorrelateMorphology:
1829 /* A Correlation is actually a Convolution with a reflected kernel.
1830 ** However a Convolution is a weighted sum with a reflected kernel.
1831 ** It may seem stange to convert a Correlation into a Convolution
1832 ** as the Correleation is the simplier method, but Convolution is
1833 ** much more commonly used, and it makes sense to implement it directly
1834 ** so as to avoid the need to duplicate the kernel when it is not
1835 ** required (which is typically the default).
1837 if ( curr_kernel == kernel )
1838 curr_kernel = CloneKernelInfo(kernel);
1839 if (curr_kernel == (KernelInfo *) NULL)
1841 RotateKernelInfo(curr_kernel,180);
1842 curr_method = ConvolveMorphology;
1843 /* FALL-THRU into Correlate (weigthed sum without reflection) */
1845 case ConvolveMorphology:
1846 { /* Scale or Normalize kernel, according to user wishes
1847 ** before using it for the Convolve/Correlate method.
1849 ** FUTURE: provide some way for internal functions to disable
1850 ** user bias and scaling effects.
1853 *artifact = GetImageArtifact(image,"convolve:scale");
1854 if ( artifact != (char *)NULL ) {
1860 if ( curr_kernel == kernel )
1861 curr_kernel = CloneKernelInfo(kernel);
1862 if (curr_kernel == (KernelInfo *) NULL)
1865 flags = (GeometryFlags) ParseGeometry(artifact, &args);
1866 ScaleKernelInfo(curr_kernel, args.rho, flags);
1869 /* FALL-THRU to do the first, and typically the only iteration */
1872 /* Do a single iteration using the Low-Level Morphology method!
1873 ** This ensures a "new_image" has been generated, but allows us to skip
1874 ** the creation of 'old_image' if no more iterations are needed.
1876 ** The "curr_method" should also be set to a low-level method that is
1877 ** understood by the MorphologyPrimative() internal function.
1879 new_image=CloneImage(image,0,0,MagickTrue,exception);
1880 if (new_image == (Image *) NULL)
1882 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1884 InheritException(exception,&new_image->exception);
1887 count++; /* iteration count */
1888 changed = MorphologyPrimative(image,new_image,curr_method,channel,
1889 curr_kernel, exception);
1890 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1891 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
1892 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1893 count, 0L, changed);
1897 * + "curr_method" should be set to a low-level morphology method
1898 * + "count=1" if the first iteration of the first kernel has been done.
1899 * + "new_image" is defined and contains the current resulting image
1902 /* The Hit-And-Miss Morphology is not 'iterated' against the resulting
1903 ** image from the previous morphology step. It must always be applied
1904 ** to the original image, with the results collected into a union (maximimum
1905 ** or lighten) of all the results. Multiple kernels may be applied but
1906 ** an iteration of the morphology does nothing, so is ignored.
1908 ** The first kernel is guranteed to have been done, so lets continue the
1909 ** process, with the other kernels in the kernel list.
1911 if ( method == HitAndMissMorphology )
1913 if ( curr_kernel->next != (KernelInfo *) NULL ) {
1914 /* create a second working image */
1915 old_image = CloneImage(image,0,0,MagickTrue,exception);
1916 if (old_image == (Image *) NULL)
1918 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1920 InheritException(exception,&old_image->exception);
1924 /* loop through rest of the kernels */
1925 this_kernel=curr_kernel->next;
1927 while( this_kernel != (KernelInfo *) NULL )
1929 changed = MorphologyPrimative(image,old_image,curr_method,channel,
1930 this_kernel,exception);
1931 (void) CompositeImageChannel(new_image,
1932 (ChannelType) (channel & ~SyncChannels), LightenCompositeOp,
1934 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1935 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
1936 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1937 count, kernel_number, changed);
1938 this_kernel = this_kernel->next;
1941 old_image=DestroyImage(old_image);
1946 /* Repeat the low-level morphology over all kernels
1947 until iteration count limit or no change from any kernel is found */
1948 if ( ( count < limit && changed > 0 ) ||
1949 curr_kernel->next != (KernelInfo *) NULL ) {
1951 /* create a second working image */
1952 old_image = CloneImage(image,0,0,MagickTrue,exception);
1953 if (old_image == (Image *) NULL)
1955 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1957 InheritException(exception,&old_image->exception);
1961 /* reset variables for the first/next iteration, or next kernel) */
1963 this_kernel = curr_kernel;
1964 total_changed = count != 0 ? changed : 0;
1965 if ( count != 0 && this_kernel != (KernelInfo *) NULL ) {
1966 count = 0; /* first iteration is not yet finished! */
1967 this_kernel = curr_kernel->next;
1969 total_changed = changed;
1972 while ( count < limit ) {
1974 while ( this_kernel != (KernelInfo *) NULL ) {
1975 Image *tmp = old_image;
1976 old_image = new_image;
1978 changed = MorphologyPrimative(old_image,new_image,curr_method,channel,
1979 this_kernel,exception);
1980 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1981 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
1982 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1983 count, kernel_number, changed);
1984 total_changed += changed;
1985 this_kernel = this_kernel->next;
1988 if ( total_changed == 0 )
1989 break; /* no changes after processing all kernels - ABORT */
1990 /* prepare for next loop */
1993 this_kernel = curr_kernel;
1995 old_image=DestroyImage(old_image);
1998 /* finished with kernel - destary any copy that was made */
1999 if ( curr_kernel != kernel )
2000 curr_kernel=DestroyKernelInfo(curr_kernel);
2002 /* Third-level Subtractive methods post-processing
2004 ** The removal of any 'Sync' channel flag in the Image Compositon below
2005 ** ensures the compose method is applied in a purely mathematical way, only
2006 ** the selected channels, without any normal 'alpha blending' normally
2007 ** associated with the compose method.
2009 ** Note "method" here is the 'original' morphological method, and not the
2010 ** 'current' morphological method used above to generate "new_image".
2013 case EdgeOutMorphology:
2014 case EdgeInMorphology:
2015 case TopHatMorphology:
2016 case BottomHatMorphology:
2017 /* Get Difference relative to the original image */
2018 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2019 DifferenceCompositeOp, image, 0, 0);
2021 case EdgeMorphology:
2022 /* Difference the Eroded image from the saved Dilated image */
2023 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2024 DifferenceCompositeOp, grad_image, 0, 0);
2025 grad_image=DestroyImage(grad_image);
2033 /* Yes goto's are bad, but in this case it makes cleanup lot more efficient */
2035 if ( new_image != (Image *) NULL )
2036 DestroyImage(new_image);
2038 if ( old_image != (Image *) NULL )
2039 DestroyImage(old_image);
2040 if ( curr_kernel != (KernelInfo *) NULL && curr_kernel != kernel )
2041 curr_kernel=DestroyKernelInfo(curr_kernel);
2046 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2050 + R o t a t e K e r n e l I n f o %
2054 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2056 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
2057 % restricted to 90 degree angles, but this may be improved in the future.
2059 % The format of the RotateKernelInfo method is:
2061 % void RotateKernelInfo(KernelInfo *kernel, double angle)
2063 % A description of each parameter follows:
2065 % o kernel: the Morphology/Convolution kernel
2067 % o angle: angle to rotate in degrees
2069 % This function is only internel to this module, as it is not finalized,
2070 % especially with regard to non-orthogonal angles, and rotation of larger
2073 static void RotateKernelInfo(KernelInfo *kernel, double angle)
2075 /* angle the lower kernels first */
2076 if ( kernel->next != (KernelInfo *) NULL)
2077 RotateKernelInfo(kernel->next, angle);
2079 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
2081 ** TODO: expand beyond simple 90 degree rotates, flips and flops
2084 /* Modulus the angle */
2085 angle = fmod(angle, 360.0);
2089 if ( 315.0 < angle || angle <= 45.0 )
2090 return; /* no change! - At least at this time */
2092 switch (kernel->type) {
2093 /* These built-in kernels are cylindrical kernels, rotating is useless */
2094 case GaussianKernel:
2095 case LaplacianKernel:
2099 case ChebyshevKernel:
2100 case ManhattenKernel:
2101 case EuclideanKernel:
2104 /* These may be rotatable at non-90 angles in the future */
2105 /* but simply rotating them in multiples of 90 degrees is useless */
2111 /* These only allows a +/-90 degree rotation (by transpose) */
2112 /* A 180 degree rotation is useless */
2114 case RectangleKernel:
2115 if ( 135.0 < angle && angle <= 225.0 )
2117 if ( 225.0 < angle && angle <= 315.0 )
2121 /* these are freely rotatable in 90 degree units */
2123 case UndefinedKernel:
2124 case UserDefinedKernel:
2127 if ( 135.0 < angle && angle <= 225.0 )
2129 /* For a 180 degree rotation - also know as a reflection */
2130 /* This is actually a very very common operation! */
2131 /* Basically all that is needed is a reversal of the kernel data! */
2138 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
2139 t=k[i], k[i]=k[j], k[j]=t;
2141 kernel->x = (long) kernel->width - kernel->x - 1;
2142 kernel->y = (long) kernel->height - kernel->y - 1;
2143 angle = fmod(angle+180.0, 360.0);
2145 if ( 45.0 < angle && angle <= 135.0 )
2146 { /* Do a transpose and a flop, of the image, which results in a 90
2147 * degree rotation using two mirror operations.
2149 * WARNING: this assumes the original image was a 1 dimentional image
2150 * but currently that is the only built-ins it is applied to.
2154 t = (long) kernel->width;
2155 kernel->width = kernel->height;
2156 kernel->height = (unsigned long) t;
2158 kernel->x = kernel->y;
2160 angle = fmod(450.0 - angle, 360.0);
2162 /* At this point angle should be between -45 (315) and +45 degrees
2163 * In the future some form of non-orthogonal angled rotates could be
2164 * performed here, posibily with a linear kernel restriction.
2168 Not currently in use!
2169 { /* Do a flop, this assumes kernel is horizontally symetrical.
2170 * Each row of the kernel needs to be reversed!
2179 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2180 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2181 t=k[x], k[x]=k[r], k[r]=t;
2183 kernel->x = kernel->width - kernel->x - 1;
2184 angle = fmod(angle+180.0, 360.0);
2191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2195 % S c a l e K e r n e l I n f o %
2199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2201 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
2202 % without normalization of the sum of the kernel values (as per given flags).
2204 % By default (no flags given) the values within the kernel is scaled
2205 % directly using given scaling factor without change.
2207 % If any 'normalize_flags' are given the kernel will first be normalized and
2208 % then further scaled by the scaling factor value given. A 'PercentValue'
2209 % flag will cause the given scaling factor to be divided by one hundred
2212 % Kernel normalization ('normalize_flags' given) is designed to ensure that
2213 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
2214 % morphology methods will fall into -1.0 to +1.0 range. Note that for
2215 % non-HDRI versions of IM this may cause images to have any negative results
2216 % clipped, unless some 'bias' is used.
2218 % More specifically. Kernels which only contain positive values (such as a
2219 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
2220 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
2222 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
2223 % the kernel will be scaled by the absolute of the sum of kernel values, so
2224 % that it will generally fall within the +/- 1.0 range.
2226 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
2227 % will be scaled by just the sum of the postive values, so that its output
2228 % range will again fall into the +/- 1.0 range.
2230 % For special kernels designed for locating shapes using 'Correlate', (often
2231 % only containing +1 and -1 values, representing foreground/brackground
2232 % matching) a special normalization method is provided to scale the positive
2233 % values seperatally to those of the negative values, so the kernel will be
2234 % forced to become a zero-sum kernel better suited to such searches.
2236 % WARNING: Correct normalization of the kernel assumes that the '*_range'
2237 % attributes within the kernel structure have been correctly set during the
2240 % NOTE: The values used for 'normalize_flags' have been selected specifically
2241 % to match the use of geometry options, so that '!' means NormalizeValue,
2242 % '^' means CorrelateNormalizeValue, and '%' means PercentValue. All other
2243 % GeometryFlags values are ignored.
2245 % The format of the ScaleKernelInfo method is:
2247 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
2248 % const MagickStatusType normalize_flags )
2250 % A description of each parameter follows:
2252 % o kernel: the Morphology/Convolution kernel
2255 % multiply all values (after normalization) by this factor if not
2256 % zero. If the kernel is normalized regardless of any flags.
2258 % o normalize_flags:
2259 % GeometryFlags defining normalization method to use.
2260 % specifically: NormalizeValue, CorrelateNormalizeValue,
2261 % and/or PercentValue
2263 % This function is internal to this module only at this time, but can be
2264 % exported to other modules if needed.
2266 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
2267 const double scaling_factor,const GeometryFlags normalize_flags)
2276 /* scale the lower kernels first */
2277 if ( kernel->next != (KernelInfo *) NULL)
2278 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
2281 if ( (normalize_flags&NormalizeValue) != 0 ) {
2282 /* normalize kernel appropriately */
2283 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
2284 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
2286 pos_scale = kernel->positive_range; /* special zero-summing kernel */
2288 /* force kernel into being a normalized zero-summing kernel */
2289 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
2290 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
2291 ? kernel->positive_range : 1.0;
2292 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
2293 ? -kernel->negative_range : 1.0;
2296 neg_scale = pos_scale;
2298 /* finialize scaling_factor for positive and negative components */
2299 pos_scale = scaling_factor/pos_scale;
2300 neg_scale = scaling_factor/neg_scale;
2301 if ( (normalize_flags&PercentValue) != 0 ) {
2306 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2307 if ( ! IsNan(kernel->values[i]) )
2308 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
2310 /* convolution output range */
2311 kernel->positive_range *= pos_scale;
2312 kernel->negative_range *= neg_scale;
2313 /* maximum and minimum values in kernel */
2314 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
2315 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
2317 /* swap kernel settings if user scaling factor is negative */
2318 if ( scaling_factor < MagickEpsilon ) {
2320 t = kernel->positive_range;
2321 kernel->positive_range = kernel->negative_range;
2322 kernel->negative_range = t;
2323 t = kernel->maximum;
2324 kernel->maximum = kernel->minimum;
2325 kernel->minimum = 1;
2332 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2336 + S h o w K e r n e l I n f o %
2340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2342 % ShowKernelInfo() outputs the details of the given kernel defination to
2343 % standard error, generally due to a users 'showkernel' option request.
2345 % The format of the ShowKernel method is:
2347 % void ShowKernelInfo(KernelInfo *kernel)
2349 % A description of each parameter follows:
2351 % o kernel: the Morphology/Convolution kernel
2353 % This function is internal to this module only at this time. That may change
2356 MagickExport void ShowKernelInfo(KernelInfo *kernel)
2364 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
2366 fprintf(stderr, "Kernel ");
2367 if ( kernel->next != (KernelInfo *) NULL )
2368 fprintf(stderr, " #%ld", c );
2369 fprintf(stderr, " \"%s\" of size %lux%lu%+ld%+ld ",
2370 MagickOptionToMnemonic(MagickKernelOptions, k->type),
2371 kernel->width, k->height,
2372 kernel->x, kernel->y );
2374 " with values from %.*lg to %.*lg\n",
2375 GetMagickPrecision(), k->minimum,
2376 GetMagickPrecision(), k->maximum);
2377 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2378 GetMagickPrecision(), k->negative_range,
2379 GetMagickPrecision(), k->positive_range,
2380 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2381 for (i=v=0; v < (long) k->height; v++) {
2382 fprintf(stderr,"%2ld:",v);
2383 for (u=0; u < (long) k->width; u++, i++)
2384 if ( IsNan(k->values[i]) )
2385 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2387 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2388 GetMagickPrecision(), k->values[i]);
2389 fprintf(stderr,"\n");
2395 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2399 + Z e r o K e r n e l N a n s %
2403 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2405 % ZeroKernelNans() replaces any special 'nan' value that may be present in
2406 % the kernel with a zero value. This is typically done when the kernel will
2407 % be used in special hardware (GPU) convolution processors, to simply
2410 % The format of the ZeroKernelNans method is:
2412 % voidZeroKernelNans (KernelInfo *kernel)
2414 % A description of each parameter follows:
2416 % o kernel: the Morphology/Convolution kernel
2418 % FUTURE: return the information in a string for API usage.
2420 MagickExport void ZeroKernelNans(KernelInfo *kernel)
2425 /* scale the lower kernels first */
2426 if ( kernel->next != (KernelInfo *) NULL)
2427 ZeroKernelNans(kernel->next);
2429 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2430 if ( IsNan(kernel->values[i]) )
2431 kernel->values[i] = 0.0;