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 ExpandKernelInfo(KernelInfo *, double),
111 RotateKernelInfo(KernelInfo *, double);
114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
118 % A c q u i r e K e r n e l I n f o %
122 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
124 % AcquireKernelInfo() takes the given string (generally supplied by the
125 % user) and converts it into a Morphology/Convolution Kernel. This allows
126 % users to specify a kernel from a number of pre-defined kernels, or to fully
127 % specify their own kernel for a specific Convolution or Morphology
130 % The kernel so generated can be any rectangular array of floating point
131 % values (doubles) with the 'control point' or 'pixel being affected'
132 % anywhere within that array of values.
134 % Previously IM was restricted to a square of odd size using the exact
135 % center as origin, this is no longer the case, and any rectangular kernel
136 % with any value being declared the origin. This in turn allows the use of
137 % highly asymmetrical kernels.
139 % The floating point values in the kernel can also include a special value
140 % known as 'nan' or 'not a number' to indicate that this value is not part
141 % of the kernel array. This allows you to shaped the kernel within its
142 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
143 % shape. However at least one non-nan value must be provided for correct
144 % working of a kernel.
146 % The returned kernel should be freed using the DestroyKernelInfo() when you
147 % are finished with it. Do not free this memory yourself.
149 % Input kernel defintion strings can consist of any of three types.
152 % Select from one of the built in kernels, using the name and
153 % geometry arguments supplied. See AcquireKernelBuiltIn()
155 % "WxH[+X+Y]:num, num, num ..."
156 % a kernel of size W by H, with W*H floating point numbers following.
157 % the 'center' can be optionally be defined at +X+Y (such that +0+0
158 % is top left corner). If not defined the pixel in the center, for
159 % odd sizes, or to the immediate top or left of center for even sizes
160 % is automatically selected.
162 % "num, num, num, num, ..."
163 % list of floating point numbers defining an 'old style' odd sized
164 % square kernel. At least 9 values should be provided for a 3x3
165 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
166 % Values can be space or comma separated. This is not recommended.
168 % Note that 'name' kernels will start with an alphabetic character while the
169 % new kernel specification has a ':' character in its specification string.
170 % If neither is the case, it is assumed an old style of a simple list of
171 % numbers generating a odd-sized square kernel has been given.
173 % You can define a 'list of kernels' which can be used by some morphology
174 % operators A list is defined as a semi-colon seperated list kernels.
176 % " kernel ; kernel ; kernel ; "
178 % Extra ';' characters are simply ignored.
180 % The format of the AcquireKernal method is:
182 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
184 % A description of each parameter follows:
186 % o kernel_string: the Morphology/Convolution kernel wanted.
190 /* This was separated so that it could be used as a separate
191 ** array input handling function, such as for -color-matrix
193 static KernelInfo *ParseKernelArray(const char *kernel_string)
199 token[MaxTextExtent];
209 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
211 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
212 if (kernel == (KernelInfo *)NULL)
214 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
215 kernel->minimum = kernel->maximum = 0.0;
216 kernel->negative_range = kernel->positive_range = 0.0;
217 kernel->type = UserDefinedKernel;
218 kernel->next = (KernelInfo *) NULL;
219 kernel->signature = MagickSignature;
221 /* find end of this specific kernel definition string */
222 end = strchr(kernel_string, ';');
223 if ( end == (char *) NULL )
224 end = strchr(kernel_string, '\0');
226 /* Has a ':' in argument - New user kernel specification */
227 p = strchr(kernel_string, ':');
228 if ( p != (char *) NULL && p < end)
236 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
237 memcpy(token, kernel_string, (size_t) (p-kernel_string));
238 token[p-kernel_string] = '\0';
239 SetGeometryInfo(&args);
240 flags = ParseGeometry(token, &args);
242 /* Size handling and checks of geometry settings */
243 if ( (flags & WidthValue) == 0 ) /* if no width then */
244 args.rho = args.sigma; /* then width = height */
245 if ( args.rho < 1.0 ) /* if width too small */
246 args.rho = 1.0; /* then width = 1 */
247 if ( args.sigma < 1.0 ) /* if height too small */
248 args.sigma = args.rho; /* then height = width */
249 kernel->width = (unsigned long)args.rho;
250 kernel->height = (unsigned long)args.sigma;
252 /* Offset Handling and Checks */
253 if ( args.xi < 0.0 || args.psi < 0.0 )
254 return(DestroyKernelInfo(kernel));
255 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
256 : (long) (kernel->width-1)/2;
257 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
258 : (long) (kernel->height-1)/2;
259 if ( kernel->x >= (long) kernel->width ||
260 kernel->y >= (long) kernel->height )
261 return(DestroyKernelInfo(kernel));
263 p++; /* advance beyond the ':' */
266 { /* ELSE - Old old specification, forming odd-square kernel */
267 /* count up number of values given */
268 p=(const char *) kernel_string;
269 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
270 p++; /* ignore "'" chars for convolve filter usage - Cristy */
271 for (i=0; p < end; i++)
273 GetMagickToken(p,&p,token);
275 GetMagickToken(p,&p,token);
277 /* set the size of the kernel - old sized square */
278 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
279 kernel->x = kernel->y = (long) (kernel->width-1)/2;
280 p=(const char *) kernel_string;
281 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
282 p++; /* ignore "'" chars for convolve filter usage - Cristy */
285 /* Read in the kernel values from rest of input string argument */
286 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
287 kernel->height*sizeof(double));
288 if (kernel->values == (double *) NULL)
289 return(DestroyKernelInfo(kernel));
291 kernel->minimum = +MagickHuge;
292 kernel->maximum = -MagickHuge;
293 kernel->negative_range = kernel->positive_range = 0.0;
295 for (i=0; (i < (long) (kernel->width*kernel->height)) && (p < end); i++)
297 GetMagickToken(p,&p,token);
299 GetMagickToken(p,&p,token);
300 if ( LocaleCompare("nan",token) == 0
301 || LocaleCompare("-",token) == 0 ) {
302 kernel->values[i] = nan; /* do not include this value in kernel */
305 kernel->values[i] = StringToDouble(token);
306 ( kernel->values[i] < 0)
307 ? ( kernel->negative_range += kernel->values[i] )
308 : ( kernel->positive_range += kernel->values[i] );
309 Minimize(kernel->minimum, kernel->values[i]);
310 Maximize(kernel->maximum, kernel->values[i]);
314 /* sanity check -- no more values in kernel definition */
315 GetMagickToken(p,&p,token);
316 if ( *token != '\0' && *token != ';' && *token != '\'' )
317 return(DestroyKernelInfo(kernel));
320 /* this was the old method of handling a incomplete kernel */
321 if ( i < (long) (kernel->width*kernel->height) ) {
322 Minimize(kernel->minimum, kernel->values[i]);
323 Maximize(kernel->maximum, kernel->values[i]);
324 for ( ; i < (long) (kernel->width*kernel->height); i++)
325 kernel->values[i]=0.0;
328 /* Number of values for kernel was not enough - Report Error */
329 if ( i < (long) (kernel->width*kernel->height) )
330 return(DestroyKernelInfo(kernel));
333 /* check that we recieved at least one real (non-nan) value! */
334 if ( kernel->minimum == MagickHuge )
335 return(DestroyKernelInfo(kernel));
340 static KernelInfo *ParseNamedKernel(const char *kernel_string)
343 token[MaxTextExtent];
358 /* Parse special 'named' kernel */
359 GetMagickToken(kernel_string,&p,token);
360 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
361 if ( type < 0 || type == UserDefinedKernel )
362 return((KernelInfo *)NULL); /* not a valid named kernel */
364 while (((isspace((int) ((unsigned char) *p)) != 0) ||
365 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
368 end = strchr(p, ';'); /* end of this kernel defintion */
369 if ( end == (char *) NULL )
370 end = strchr(p, '\0');
372 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
373 memcpy(token, p, (size_t) (end-p));
375 SetGeometryInfo(&args);
376 flags = ParseGeometry(token, &args);
379 /* For Debugging Geometry Input */
380 fprintf(stderr, "Geometry = %04x : %lf x %lf %+lf %+lf\n",
381 flags, args.rho, args.sigma, args.xi, args.psi );
384 /* special handling of missing values in input string */
386 case RectangleKernel:
387 if ( (flags & WidthValue) == 0 ) /* if no width then */
388 args.rho = args.sigma; /* then width = height */
389 if ( args.rho < 1.0 ) /* if width too small */
390 args.rho = 3; /* then width = 3 */
391 if ( args.sigma < 1.0 ) /* if height too small */
392 args.sigma = args.rho; /* then height = width */
393 if ( (flags & XValue) == 0 ) /* center offset if not defined */
394 args.xi = (double)(((long)args.rho-1)/2);
395 if ( (flags & YValue) == 0 )
396 args.psi = (double)(((long)args.sigma-1)/2);
402 /* If no scale given (a 0 scale is valid! - set it to 1.0 */
403 if ( (flags & HeightValue) == 0 )
406 case ChebyshevKernel:
407 case ManhattenKernel:
408 case EuclideanKernel:
409 if ( (flags & HeightValue) == 0 )
410 args.sigma = 100.0; /* default distance scaling */
411 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
412 args.sigma = QuantumRange/args.sigma; /* maximum pixel distance */
413 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
414 args.sigma *= QuantumRange/100.0; /* percentage of color range */
420 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
423 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
432 token[MaxTextExtent];
441 kernel = last_kernel = NULL;
444 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
446 /* ignore multiple ';' following each other */
447 if ( *token != ';' ) {
449 /* tokens starting with alpha is a Named kernel */
450 if (isalpha((int) *token) == 0)
451 new_kernel = ParseKernelArray(p);
452 else /* otherwise a user defined kernel array */
453 new_kernel = ParseNamedKernel(p);
455 /* Error handling -- this is not proper error handling! */
456 if ( new_kernel == (KernelInfo *) NULL ) {
457 fprintf(stderr, "Failed to parse kernel number #%lu\n", kernel_number);
458 if ( kernel != (KernelInfo *) NULL )
459 kernel=DestroyKernelInfo(kernel);
460 return((KernelInfo *) NULL);
463 /* initialise or append the kernel list */
464 if ( last_kernel == (KernelInfo *) NULL )
465 kernel = last_kernel = new_kernel;
467 last_kernel->next = new_kernel;
468 last_kernel = new_kernel;
472 /* look for the next kernel in list */
474 if ( p == (char *) NULL )
484 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
488 % A c q u i r e K e r n e l B u i l t I n %
492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
494 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
495 % kernels used for special purposes such as gaussian blurring, skeleton
496 % pruning, and edge distance determination.
498 % They take a KernelType, and a set of geometry style arguments, which were
499 % typically decoded from a user supplied string, or from a more complex
500 % Morphology Method that was requested.
502 % The format of the AcquireKernalBuiltIn method is:
504 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
505 % const GeometryInfo args)
507 % A description of each parameter follows:
509 % o type: the pre-defined type of kernel wanted
511 % o args: arguments defining or modifying the kernel
513 % Convolution Kernels
515 % Gaussian:{radius},{sigma}
516 % Generate a two-dimentional gaussian kernel, as used by -gaussian.
517 % A sigma is required.
519 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
520 % the final size of the resulting kernel to a square 2*radius+1 in size.
521 % The radius should be at least 2 times that of the sigma value, or
522 % sever clipping and aliasing may result. If not given or set to 0 the
523 % radius will be determined so as to produce the best minimal error
524 % result, which is usally much larger than is normally needed.
526 % Blur:{radius},{sigma},{angle}
527 % As per Gaussian, but generates a 1 dimensional or linear gaussian
528 % blur, at the angle given (current restricted to orthogonal angles).
529 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
530 % Angle can be rotated in multiples of 90 degrees.
532 % Note that two such blurs perpendicular to each other is equivelent to
533 % the far large "Gaussian" kernel, but much faster to apply. This is how
534 % the -blur operator works.
536 % Comet:{width},{sigma},{angle}
537 % Blur in one direction only, much like how a bright object leaves
538 % a comet like trail. The Kernel is actually half a gaussian curve,
539 % Adding two such blurs in opposite directions produces a Blur Kernel.
540 % Angle can be rotated in multiples of 90 degrees.
542 % Note that the first argument is the width of the kernel and not the
543 % radius of the kernel.
545 % # Still to be implemented...
547 % # Sharpen:{radius},{sigma}
548 % # Negated Gaussian (center zeroed and re-normalized),
549 % # with a 2 unit positive peak. -- Check On line documentation
551 % # LOG:{radius},{sigma1},{sigma2}
552 % # Laplacian of Gaussian
554 % # DOG:{radius},{sigma1},{sigma2}
555 % # Difference of two Gaussians
559 % # Set kernel values using a resize filter, and given scale (sigma)
560 % # Cylindrical or Linear. Is this posible with an image?
563 % Named Constant Convolution Kernels
566 % The 3x3 sobel convolution kernel. Angle may be given in multiples
567 % of 45 degrees. Kernel is unscaled by default so some normalization
568 % may be required to ensure results are not clipped.
569 % Default kernel is -1,0,1
574 % Generate Lapacian kernel of the type specified. (1 is the default)
575 % Type 0 default square laplacian 3x3: all -1's with central 8
576 % Type 1 3x3: central 4 edge -1 corner 0
577 % Type 2 3x3: central 4 edge 1 corner -2
578 % Type 3 a 5x5 laplacian
579 % Type 5 a 7x7 laplacian
583 % Rectangle:{geometry}
584 % Simply generate a rectangle of 1's with the size given. You can also
585 % specify the location of the 'control point', otherwise the closest
586 % pixel to the center of the rectangle is selected.
588 % Properly centered and odd sized rectangles work the best.
590 % Diamond:[{radius}[,{scale}]]
591 % Generate a diamond shaped kernel with given radius to the points.
592 % Kernel size will again be radius*2+1 square and defaults to radius 1,
593 % generating a 3x3 kernel that is slightly larger than a square.
595 % Square:[{radius}[,{scale}]]
596 % Generate a square shaped kernel of size radius*2+1, and defaulting
597 % to a 3x3 (radius 1).
599 % Note that using a larger radius for the "Square" or the "Diamond"
600 % is also equivelent to iterating the basic morphological method
601 % that many times. However However iterating with the smaller radius 1
602 % default is actually faster than using a larger kernel radius.
604 % Disk:[{radius}[,{scale}]]
605 % Generate a binary disk of the radius given, radius may be a float.
606 % Kernel size will be ceil(radius)*2+1 square.
607 % NOTE: Here are some disk shapes of specific interest
608 % "disk:1" => "diamond" or "cross:1"
609 % "disk:1.5" => "square"
610 % "disk:2" => "diamond:2"
611 % "disk:2.5" => a general disk shape of radius 2
612 % "disk:2.9" => "square:2"
613 % "disk:3.5" => default - octagonal/disk shape of radius 3
614 % "disk:4.2" => roughly octagonal shape of radius 4
615 % "disk:4.3" => a general disk shape of radius 4
616 % After this all the kernel shape becomes more and more circular.
618 % Because a "disk" is more circular when using a larger radius, using a
619 % larger radius is preferred over iterating the morphological operation.
621 % Plus:[{radius}[,{scale}]]
622 % Generate a kernel in the shape of a 'plus' sign. The length of each
623 % arm is also the radius, which defaults to 2.
625 % This kernel is not a good general morphological kernel, but is used
626 % more for highlighting and marking any single pixels in an image using,
627 % a "Dilate" or "Erode" method as appropriate.
629 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
631 % Note that unlike other kernels iterating a plus does not produce the
632 % same result as using a larger radius for the cross.
634 % Distance Measuring Kernels
636 % Chebyshev "[{radius}][x{scale}[%!]]"
637 % Manhatten "[{radius}][x{scale}[%!]]"
638 % Euclidean "[{radius}][x{scale}[%!]]"
640 % Different types of distance measuring methods, which are used with the
641 % a 'Distance' morphology method for generating a gradient based on
642 % distance from an edge of a binary shape, though there is a technique
643 % for handling a anti-aliased shape.
645 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
646 % one to any neighbour, orthogonal or diagonal. One why of thinking of
647 % it is the number of squares a 'King' or 'Queen' in chess needs to
648 % traverse reach any other position on a chess board. It results in a
649 % 'square' like distance function, but one where diagonals are closer
652 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
653 % Cab metric), is the distance needed when you can only travel in
654 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
655 % in chess would travel. It results in a diamond like distances, where
656 % diagonals are further than expected.
658 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
659 % However by default the kernel size only has a radius of 1, which
660 % limits the distance to 'Knight' like moves, with only orthogonal and
661 % diagonal measurements being correct. As such for the default kernel
662 % you will get octagonal like distance function, which is reasonally
665 % However if you use a larger radius such as "Euclidean:4" you will
666 % get a much smoother distance gradient from the edge of the shape.
667 % Of course a larger kernel is slower to use, and generally not needed.
669 % To allow the use of fractional distances that you get with diagonals
670 % the actual distance is scaled by a fixed value which the user can
671 % provide. This is not actually nessary for either ""Chebyshev" or
672 % "Manhatten" distance kernels, but is done for all three distance
673 % kernels. If no scale is provided it is set to a value of 100,
674 % allowing for a maximum distance measurement of 655 pixels using a Q16
675 % version of IM, from any edge. However for small images this can
676 % result in quite a dark gradient.
678 % See the 'Distance' Morphological Method, for information of how it is
681 % # Hit-n-Miss Kernel-Lists -- Still to be implemented
683 % # specifically for Pruning, Thinning, Thickening
687 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
688 const GeometryInfo *args)
701 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
703 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
704 if (kernel == (KernelInfo *) NULL)
706 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
707 kernel->minimum = kernel->maximum = 0.0;
708 kernel->negative_range = kernel->positive_range = 0.0;
710 kernel->next = (KernelInfo *) NULL;
711 kernel->signature = MagickSignature;
714 /* Convolution Kernels */
717 sigma = fabs(args->sigma);
719 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
721 kernel->width = kernel->height =
722 GetOptimalKernelWidth2D(args->rho,sigma);
723 kernel->x = kernel->y = (long) (kernel->width-1)/2;
724 kernel->negative_range = kernel->positive_range = 0.0;
725 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
726 kernel->height*sizeof(double));
727 if (kernel->values == (double *) NULL)
728 return(DestroyKernelInfo(kernel));
730 sigma = 2.0*sigma*sigma; /* simplify the expression */
731 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
732 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
733 kernel->positive_range += (
735 exp(-((double)(u*u+v*v))/sigma)
736 /* / (MagickPI*sigma) */ );
738 kernel->maximum = kernel->values[
739 kernel->y*kernel->width+kernel->x ];
741 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
747 sigma = fabs(args->sigma);
749 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
751 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
752 kernel->x = (long) (kernel->width-1)/2;
755 kernel->negative_range = kernel->positive_range = 0.0;
756 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
757 kernel->height*sizeof(double));
758 if (kernel->values == (double *) NULL)
759 return(DestroyKernelInfo(kernel));
763 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
764 ** It generates a gaussian 3 times the width, and compresses it into
765 ** the expected range. This produces a closer normalization of the
766 ** resulting kernel, especially for very low sigma values.
767 ** As such while wierd it is prefered.
769 ** I am told this method originally came from Photoshop.
771 sigma *= KernelRank; /* simplify expanded curve */
772 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
773 (void) ResetMagickMemory(kernel->values,0, (size_t)
774 kernel->width*sizeof(double));
775 for ( u=-v; u <= v; u++) {
776 kernel->values[(u+v)/KernelRank] +=
777 exp(-((double)(u*u))/(2.0*sigma*sigma))
778 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
780 for (i=0; i < (long) kernel->width; i++)
781 kernel->positive_range += kernel->values[i];
783 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
784 kernel->positive_range += (
786 exp(-((double)(u*u))/(2.0*sigma*sigma))
787 /* / (MagickSQ2PI*sigma) */ );
790 kernel->maximum = kernel->values[ kernel->x ];
791 /* Note that neither methods above generate a normalized kernel,
792 ** though it gets close. The kernel may be 'clipped' by a user defined
793 ** radius, producing a smaller (darker) kernel. Also for very small
794 ** sigma's (> 0.1) the central value becomes larger than one, and thus
795 ** producing a very bright kernel.
798 /* Normalize the 1D Gaussian Kernel
800 ** Because of this the divisor in the above kernel generator is
801 ** not needed, so is not done above.
803 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
805 /* rotate the kernel by given angle */
806 RotateKernelInfo(kernel, args->xi);
811 sigma = fabs(args->sigma);
813 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
815 if ( args->rho < 1.0 )
816 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
818 kernel->width = (unsigned long)args->rho;
819 kernel->x = kernel->y = 0;
821 kernel->negative_range = kernel->positive_range = 0.0;
822 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
823 kernel->height*sizeof(double));
824 if (kernel->values == (double *) NULL)
825 return(DestroyKernelInfo(kernel));
827 /* A comet blur is half a gaussian curve, so that the object is
828 ** blurred in one direction only. This may not be quite the right
829 ** curve so may change in the future. The function must be normalised.
833 sigma *= KernelRank; /* simplify expanded curve */
834 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
835 (void) ResetMagickMemory(kernel->values,0, (size_t)
836 kernel->width*sizeof(double));
837 for ( u=0; u < v; u++) {
838 kernel->values[u/KernelRank] +=
839 exp(-((double)(u*u))/(2.0*sigma*sigma))
840 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
842 for (i=0; i < (long) kernel->width; i++)
843 kernel->positive_range += kernel->values[i];
845 for ( i=0; i < (long) kernel->width; i++)
846 kernel->positive_range += (
848 exp(-((double)(i*i))/(2.0*sigma*sigma))
849 /* / (MagickSQ2PI*sigma) */ );
852 kernel->maximum = kernel->values[0];
854 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
855 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
858 /* Convolution Kernels - Well Known Constants */
860 { kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
861 kernel=ParseKernelArray("3x3:-1,0,1 -2,0,2 -1,0,1");
862 if (kernel == (KernelInfo *) NULL)
865 RotateKernelInfo(kernel, args->rho); /* Rotate by angle */
868 case LaplacianKernel:
869 { kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
870 switch ( (int) args->rho ) {
872 kernel=ParseKernelArray("3x3: 0,-1,0 -1,4,-1 0,-1,0");
875 kernel=ParseKernelArray("3x3: -2,1,-2 1,4,1 -2,1,-2");
878 default: /* default laplacian 'edge' filter */
879 kernel=ParseKernelArray("3x3: -1,-1,-1 -1,8,-1 -1,-1,-1");
881 case 4: /* a 5x5 laplacian */
882 kernel=ParseKernelArray(
883 "5x5: -4,-1,0,-1,-4 -1,2,3,2,-1 0,3,4,3,0 -1,2,3,2,-1 -4,-1,0,-1,-4");
885 case 5: /* a 7x7 laplacian */
886 kernel=ParseKernelArray(
887 "7x7:-10,-5,-2,-1,-2,-5,-10 -5,0,3,4,3,0,-5 -2,3,6,7,6,3,-2 -1,4,7,8,7,4,-1 -2,3,6,7,6,3,-2 -5,0,3,4,3,0,-5 -10,-5,-2,-1,-2,-5,-10" );
890 if (kernel == (KernelInfo *) NULL)
895 /* Boolean Kernels */
896 case RectangleKernel:
900 if ( type == SquareKernel )
903 kernel->width = kernel->height = 3; /* default radius = 1 */
905 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
906 kernel->x = kernel->y = (long) (kernel->width-1)/2;
910 /* NOTE: user defaults set in "AcquireKernelInfo()" */
911 if ( args->rho < 1.0 || args->sigma < 1.0 )
912 return(DestroyKernelInfo(kernel)); /* invalid args given */
913 kernel->width = (unsigned long)args->rho;
914 kernel->height = (unsigned long)args->sigma;
915 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
916 args->psi < 0.0 || args->psi > (double)kernel->height )
917 return(DestroyKernelInfo(kernel)); /* invalid args given */
918 kernel->x = (long) args->xi;
919 kernel->y = (long) args->psi;
922 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
923 kernel->height*sizeof(double));
924 if (kernel->values == (double *) NULL)
925 return(DestroyKernelInfo(kernel));
927 /* set all kernel values to 1.0 */
928 u=(long) kernel->width*kernel->height;
929 for ( i=0; i < u; i++)
930 kernel->values[i] = scale;
931 kernel->minimum = kernel->maximum = scale; /* a flat shape */
932 kernel->positive_range = scale*u;
938 kernel->width = kernel->height = 3; /* default radius = 1 */
940 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
941 kernel->x = kernel->y = (long) (kernel->width-1)/2;
943 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
944 kernel->height*sizeof(double));
945 if (kernel->values == (double *) NULL)
946 return(DestroyKernelInfo(kernel));
948 /* set all kernel values within diamond area to scale given */
949 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
950 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
951 if ((labs(u)+labs(v)) <= (long)kernel->x)
952 kernel->positive_range += kernel->values[i] = args->sigma;
954 kernel->values[i] = nan;
955 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
963 limit = (long)(args->rho*args->rho);
964 if (args->rho < 0.1) /* default radius approx 3.5 */
965 kernel->width = kernel->height = 7L, limit = 10L;
967 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
968 kernel->x = kernel->y = (long) (kernel->width-1)/2;
970 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
971 kernel->height*sizeof(double));
972 if (kernel->values == (double *) NULL)
973 return(DestroyKernelInfo(kernel));
975 /* set all kernel values within disk area to 1.0 */
976 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
977 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
978 if ((u*u+v*v) <= limit)
979 kernel->positive_range += kernel->values[i] = args->sigma;
981 kernel->values[i] = nan;
982 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
988 kernel->width = kernel->height = 5; /* default radius 2 */
990 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
991 kernel->x = kernel->y = (long) (kernel->width-1)/2;
993 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
994 kernel->height*sizeof(double));
995 if (kernel->values == (double *) NULL)
996 return(DestroyKernelInfo(kernel));
998 /* set all kernel values along axises to 1.0 */
999 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1000 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1001 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1002 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1003 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1006 /* Distance Measuring Kernels */
1007 case ChebyshevKernel:
1009 if (args->rho < 1.0)
1010 kernel->width = kernel->height = 3; /* default radius = 1 */
1012 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
1013 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1015 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1016 kernel->height*sizeof(double));
1017 if (kernel->values == (double *) NULL)
1018 return(DestroyKernelInfo(kernel));
1020 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1021 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1022 kernel->positive_range += ( kernel->values[i] =
1023 args->sigma*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
1024 kernel->maximum = kernel->values[0];
1027 case ManhattenKernel:
1029 if (args->rho < 1.0)
1030 kernel->width = kernel->height = 3; /* default radius = 1 */
1032 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
1033 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1035 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1036 kernel->height*sizeof(double));
1037 if (kernel->values == (double *) NULL)
1038 return(DestroyKernelInfo(kernel));
1040 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1041 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1042 kernel->positive_range += ( kernel->values[i] =
1043 args->sigma*(labs(u)+labs(v)) );
1044 kernel->maximum = kernel->values[0];
1047 case EuclideanKernel:
1049 if (args->rho < 1.0)
1050 kernel->width = kernel->height = 3; /* default radius = 1 */
1052 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
1053 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1055 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1056 kernel->height*sizeof(double));
1057 if (kernel->values == (double *) NULL)
1058 return(DestroyKernelInfo(kernel));
1060 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1061 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1062 kernel->positive_range += ( kernel->values[i] =
1063 args->sigma*sqrt((double)(u*u+v*v)) );
1064 kernel->maximum = kernel->values[0];
1067 /* Undefined Kernels */
1070 perror("Kernel Type has not been defined yet\n");
1073 /* Generate a No-Op minimal kernel - 1x1 pixel */
1074 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
1075 if (kernel->values == (double *) NULL)
1076 return(DestroyKernelInfo(kernel));
1077 kernel->width = kernel->height = 1;
1078 kernel->x = kernel->x = 0;
1079 kernel->type = UndefinedKernel;
1081 kernel->positive_range =
1082 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
1090 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1094 % C l o n e K e r n e l I n f o %
1098 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1100 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
1101 % can be modified without effecting the original. The cloned kernel should
1102 % be destroyed using DestoryKernelInfo() when no longer needed.
1104 % The format of the CloneKernelInfo method is:
1106 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1108 % A description of each parameter follows:
1110 % o kernel: the Morphology/Convolution kernel to be cloned
1113 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1121 assert(kernel != (KernelInfo *) NULL);
1122 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1123 if (new_kernel == (KernelInfo *) NULL)
1125 *new_kernel=(*kernel); /* copy values in structure */
1127 /* replace the values with a copy of the values */
1128 new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1129 kernel->height*sizeof(double));
1130 if (new_kernel->values == (double *) NULL)
1131 return(DestroyKernelInfo(new_kernel));
1132 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1133 new_kernel->values[i]=kernel->values[i];
1135 /* Also clone the next kernel in the kernel list */
1136 if ( kernel->next != (KernelInfo *) NULL ) {
1137 new_kernel->next = CloneKernelInfo(kernel->next);
1138 if ( new_kernel->next == (KernelInfo *) NULL )
1139 return(DestroyKernelInfo(new_kernel));
1146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1150 % D e s t r o y K e r n e l I n f o %
1154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1156 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
1159 % The format of the DestroyKernelInfo method is:
1161 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1163 % A description of each parameter follows:
1165 % o kernel: the Morphology/Convolution kernel to be destroyed
1169 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1171 assert(kernel != (KernelInfo *) NULL);
1173 if ( kernel->next != (KernelInfo *) NULL )
1174 kernel->next = DestroyKernelInfo(kernel->next);
1176 kernel->values = (double *)RelinquishMagickMemory(kernel->values);
1177 kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
1182 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1186 % E x p a n d K e r n e l I n f o %
1190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1192 % ExpandKernelInfo() takes a single kernel, and expands it into a list
1193 % of kernels each incrementally rotated the angle given.
1195 % WARNING: 45 degree rotations only works for 3x3 kernels.
1196 % While 90 degree roatations only works for linear and square kernels
1198 % The format of the RotateKernelInfo method is:
1200 % void ExpandKernelInfo(KernelInfo *kernel, double angle)
1202 % A description of each parameter follows:
1204 % o kernel: the Morphology/Convolution kernel
1206 % o angle: angle to rotate in degrees
1208 % This function is only internel to this module, as it is not finalized,
1209 % especially with regard to non-orthogonal angles, and rotation of larger
1212 static void ExpandKernelInfo(KernelInfo *kernel, double angle)
1221 for (a=angle; a<355.0; a+=angle) {
1222 new = CloneKernelInfo(last);
1223 RotateKernelInfo(new, angle);
1230 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1234 % M o r p h o l o g y I m a g e C h a n n e l %
1238 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1240 % MorphologyImageChannel() applies a user supplied kernel to the image
1241 % according to the given mophology method.
1243 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1244 % by the kernel generator.
1246 % The format of the MorphologyImage method is:
1248 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
1249 % const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
1250 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1251 % channel,MorphologyMethod method,const long iterations,
1252 % KernelInfo *kernel,ExceptionInfo *exception)
1254 % A description of each parameter follows:
1256 % o image: the image.
1258 % o method: the morphology method to be applied.
1260 % o iterations: apply the operation this many times (or no change).
1261 % A value of -1 means loop until no change found.
1262 % How this is applied may depend on the morphology method.
1263 % Typically this is a value of 1.
1265 % o channel: the channel type.
1267 % o kernel: An array of double representing the morphology kernel.
1268 % Warning: kernel may be normalized for the Convolve method.
1270 % o exception: return any errors or warnings in this structure.
1273 % TODO: bias and auto-scale handling of the kernel for convolution
1274 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1275 % by the kernel generator.
1280 /* Internal function
1281 * Apply the low-level Morphology Method Primatives using the given Kernel
1282 * Returning the number of pixels that changed.
1283 * Two pre-created images must be provided, no image is created.
1285 static unsigned long MorphologyPrimative(const Image *image, Image
1286 *result_image, const MorphologyMethod method, const ChannelType channel,
1287 const KernelInfo *kernel, ExceptionInfo *exception)
1289 #define MorphologyTag "Morphology/Image"
1306 /* Only the most basic morphology is actually performed by this routine */
1309 Apply Basic Morphology to Image.
1315 GetMagickPixelPacket(image,&bias);
1316 SetMagickPixelPacketBias(image,&bias);
1317 /* Future: handle auto-bias from user, based on kernel input */
1319 p_view=AcquireCacheView(image);
1320 q_view=AcquireCacheView(result_image);
1322 /* Some methods (including convolve) needs use a reflected kernel.
1323 * Adjust 'origin' offsets for this reflected kernel.
1328 case ConvolveMorphology:
1329 case DilateMorphology:
1330 case DilateIntensityMorphology:
1331 case DistanceMorphology:
1332 /* kernel needs to used with reflection about origin */
1333 offx = (long) kernel->width-offx-1;
1334 offy = (long) kernel->height-offy-1;
1336 case ErodeMorphology:
1337 case ErodeIntensityMorphology:
1338 case HitAndMissMorphology:
1339 case ThinningMorphology:
1340 case ThickenMorphology:
1341 /* kernel is user as is, without reflection */
1344 perror("Not a low level Morphology Method");
1348 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1349 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1351 for (y=0; y < (long) image->rows; y++)
1356 register const PixelPacket
1359 register const IndexPacket
1360 *restrict p_indexes;
1362 register PixelPacket
1365 register IndexPacket
1366 *restrict q_indexes;
1374 if (status == MagickFalse)
1376 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1377 image->columns+kernel->width, kernel->height, exception);
1378 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1380 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1385 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1386 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1387 r = (image->columns+kernel->width)*offy+offx; /* constant */
1389 for (x=0; x < (long) image->columns; x++)
1397 register const double
1400 register const PixelPacket
1403 register const IndexPacket
1404 *restrict k_indexes;
1411 /* Copy input to ouput image for unused channels
1412 * This removes need for 'cloning' a new image every iteration
1415 if (image->colorspace == CMYKColorspace)
1416 q_indexes[x] = p_indexes[r];
1423 min.index = (MagickRealType) QuantumRange;
1428 max.index = (MagickRealType) 0;
1429 /* original pixel value */
1430 result.red = (MagickRealType) p[r].red;
1431 result.green = (MagickRealType) p[r].green;
1432 result.blue = (MagickRealType) p[r].blue;
1433 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1435 if ( image->colorspace == CMYKColorspace)
1436 result.index = (MagickRealType) p_indexes[r];
1439 case ConvolveMorphology:
1440 /* Set the user defined bias of the weighted average output
1442 ** FUTURE: provide some way for internal functions to disable
1443 ** user provided bias and scaling effects.
1447 case DilateIntensityMorphology:
1448 case ErodeIntensityMorphology:
1449 result.red = 0.0; /* flag indicating when first match found */
1456 case ConvolveMorphology:
1457 /* Weighted Average of pixels using reflected kernel
1459 ** NOTE for correct working of this operation for asymetrical
1460 ** kernels, the kernel needs to be applied in its reflected form.
1461 ** That is its values needs to be reversed.
1463 ** Correlation is actually the same as this but without reflecting
1464 ** the kernel, and thus 'lower-level' that Convolution. However
1465 ** as Convolution is the more common method used, and it does not
1466 ** really cost us much in terms of processing to use a reflected
1467 ** kernel, so it is Convolution that is implemented.
1469 ** Correlation will have its kernel reflected before calling
1470 ** this function to do a Convolve.
1472 ** For more details of Correlation vs Convolution see
1473 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
1475 if (((channel & SyncChannels) != 0 ) &&
1476 (image->matte == MagickTrue))
1477 { /* Channel has a 'Sync' Flag, and Alpha Channel enabled.
1478 ** Weight the color channels with Alpha Channel so that
1479 ** transparent pixels are not part of the results.
1482 alpha, /* color channel weighting : kernel*alpha */
1483 gamma; /* divisor, sum of weighting values */
1486 k = &kernel->values[ kernel->width*kernel->height-1 ];
1488 k_indexes = p_indexes;
1489 for (v=0; v < (long) kernel->height; v++) {
1490 for (u=0; u < (long) kernel->width; u++, k--) {
1491 if ( IsNan(*k) ) continue;
1492 alpha=(*k)*(QuantumScale*(QuantumRange-
1493 k_pixels[u].opacity));
1495 result.red += alpha*k_pixels[u].red;
1496 result.green += alpha*k_pixels[u].green;
1497 result.blue += alpha*k_pixels[u].blue;
1498 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1499 if ( image->colorspace == CMYKColorspace)
1500 result.index += alpha*k_indexes[u];
1502 k_pixels += image->columns+kernel->width;
1503 k_indexes += image->columns+kernel->width;
1505 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1506 result.red *= gamma;
1507 result.green *= gamma;
1508 result.blue *= gamma;
1509 result.opacity *= gamma;
1510 result.index *= gamma;
1514 /* No 'Sync' flag, or no Alpha involved.
1515 ** Convolution is simple individual channel weigthed sum.
1517 k = &kernel->values[ kernel->width*kernel->height-1 ];
1519 k_indexes = p_indexes;
1520 for (v=0; v < (long) kernel->height; v++) {
1521 for (u=0; u < (long) kernel->width; u++, k--) {
1522 if ( IsNan(*k) ) continue;
1523 result.red += (*k)*k_pixels[u].red;
1524 result.green += (*k)*k_pixels[u].green;
1525 result.blue += (*k)*k_pixels[u].blue;
1526 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1527 if ( image->colorspace == CMYKColorspace)
1528 result.index += (*k)*k_indexes[u];
1530 k_pixels += image->columns+kernel->width;
1531 k_indexes += image->columns+kernel->width;
1536 case ErodeMorphology:
1537 /* Minimum Value within kernel neighbourhood
1539 ** NOTE that the kernel is not reflected for this operation!
1541 ** NOTE: in normal Greyscale Morphology, the kernel value should
1542 ** be added to the real value, this is currently not done, due to
1543 ** the nature of the boolean kernels being used.
1547 k_indexes = p_indexes;
1548 for (v=0; v < (long) kernel->height; v++) {
1549 for (u=0; u < (long) kernel->width; u++, k++) {
1550 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1551 Minimize(min.red, (double) k_pixels[u].red);
1552 Minimize(min.green, (double) k_pixels[u].green);
1553 Minimize(min.blue, (double) k_pixels[u].blue);
1554 Minimize(min.opacity,
1555 QuantumRange-(double) k_pixels[u].opacity);
1556 if ( image->colorspace == CMYKColorspace)
1557 Minimize(min.index, (double) k_indexes[u]);
1559 k_pixels += image->columns+kernel->width;
1560 k_indexes += image->columns+kernel->width;
1565 case DilateMorphology:
1566 /* Maximum Value within kernel neighbourhood
1568 ** NOTE for correct working of this operation for asymetrical
1569 ** kernels, the kernel needs to be applied in its reflected form.
1570 ** That is its values needs to be reversed.
1572 ** NOTE: in normal Greyscale Morphology, the kernel value should
1573 ** be added to the real value, this is currently not done, due to
1574 ** the nature of the boolean kernels being used.
1577 k = &kernel->values[ kernel->width*kernel->height-1 ];
1579 k_indexes = p_indexes;
1580 for (v=0; v < (long) kernel->height; v++) {
1581 for (u=0; u < (long) kernel->width; u++, k--) {
1582 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1583 Maximize(max.red, (double) k_pixels[u].red);
1584 Maximize(max.green, (double) k_pixels[u].green);
1585 Maximize(max.blue, (double) k_pixels[u].blue);
1586 Maximize(max.opacity,
1587 QuantumRange-(double) k_pixels[u].opacity);
1588 if ( image->colorspace == CMYKColorspace)
1589 Maximize(max.index, (double) k_indexes[u]);
1591 k_pixels += image->columns+kernel->width;
1592 k_indexes += image->columns+kernel->width;
1596 case HitAndMissMorphology:
1597 case ThinningMorphology:
1598 case ThickenMorphology:
1599 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
1601 ** NOTE that the kernel is not reflected for this operation,
1602 ** and consists of both foreground and background pixel
1603 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
1604 ** with either Nan or 0.5 values for don't care.
1606 ** Note that this can produce negative results, though really
1607 ** only a positive match has any real value.
1611 k_indexes = p_indexes;
1612 for (v=0; v < (long) kernel->height; v++) {
1613 for (u=0; u < (long) kernel->width; u++, k++) {
1614 if ( IsNan(*k) ) continue;
1616 { /* minimim of foreground pixels */
1617 Minimize(min.red, (double) k_pixels[u].red);
1618 Minimize(min.green, (double) k_pixels[u].green);
1619 Minimize(min.blue, (double) k_pixels[u].blue);
1620 Minimize(min.opacity,
1621 QuantumRange-(double) k_pixels[u].opacity);
1622 if ( image->colorspace == CMYKColorspace)
1623 Minimize(min.index, (double) k_indexes[u]);
1625 else if ( (*k) < 0.3 )
1626 { /* maximum of background pixels */
1627 Maximize(max.red, (double) k_pixels[u].red);
1628 Maximize(max.green, (double) k_pixels[u].green);
1629 Maximize(max.blue, (double) k_pixels[u].blue);
1630 Maximize(max.opacity,
1631 QuantumRange-(double) k_pixels[u].opacity);
1632 if ( image->colorspace == CMYKColorspace)
1633 Maximize(max.index, (double) k_indexes[u]);
1636 k_pixels += image->columns+kernel->width;
1637 k_indexes += image->columns+kernel->width;
1639 /* Pattern Match only if min fg larger than min bg pixels */
1640 min.red -= max.red; Maximize( min.red, 0.0 );
1641 min.green -= max.green; Maximize( min.green, 0.0 );
1642 min.blue -= max.blue; Maximize( min.blue, 0.0 );
1643 min.opacity -= max.opacity; Maximize( min.opacity, 0.0 );
1644 min.index -= max.index; Maximize( min.index, 0.0 );
1647 case ErodeIntensityMorphology:
1648 /* Select Pixel with Minimum Intensity within kernel neighbourhood
1650 ** WARNING: the intensity test fails for CMYK and does not
1651 ** take into account the moderating effect of teh alpha channel
1652 ** on the intensity.
1654 ** NOTE that the kernel is not reflected for this operation!
1658 k_indexes = p_indexes;
1659 for (v=0; v < (long) kernel->height; v++) {
1660 for (u=0; u < (long) kernel->width; u++, k++) {
1661 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1662 if ( result.red == 0.0 ||
1663 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1664 /* copy the whole pixel - no channel selection */
1666 if ( result.red > 0.0 ) changed++;
1670 k_pixels += image->columns+kernel->width;
1671 k_indexes += image->columns+kernel->width;
1675 case DilateIntensityMorphology:
1676 /* Select Pixel with Maximum Intensity within kernel neighbourhood
1678 ** WARNING: the intensity test fails for CMYK and does not
1679 ** take into account the moderating effect of teh alpha channel
1680 ** on the intensity.
1682 ** NOTE for correct working of this operation for asymetrical
1683 ** kernels, the kernel needs to be applied in its reflected form.
1684 ** That is its values needs to be reversed.
1686 k = &kernel->values[ kernel->width*kernel->height-1 ];
1688 k_indexes = p_indexes;
1689 for (v=0; v < (long) kernel->height; v++) {
1690 for (u=0; u < (long) kernel->width; u++, k--) {
1691 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1692 if ( result.red == 0.0 ||
1693 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1694 /* copy the whole pixel - no channel selection */
1696 if ( result.red > 0.0 ) changed++;
1700 k_pixels += image->columns+kernel->width;
1701 k_indexes += image->columns+kernel->width;
1706 case DistanceMorphology:
1707 /* Add kernel Value and select the minimum value found.
1708 ** The result is a iterative distance from edge of image shape.
1710 ** All Distance Kernels are symetrical, but that may not always
1711 ** be the case. For example how about a distance from left edges?
1712 ** To work correctly with asymetrical kernels the reflected kernel
1713 ** needs to be applied.
1715 ** Actually this is really a GreyErode with a negative kernel!
1718 k = &kernel->values[ kernel->width*kernel->height-1 ];
1720 k_indexes = p_indexes;
1721 for (v=0; v < (long) kernel->height; v++) {
1722 for (u=0; u < (long) kernel->width; u++, k--) {
1723 if ( IsNan(*k) ) continue;
1724 Minimize(result.red, (*k)+k_pixels[u].red);
1725 Minimize(result.green, (*k)+k_pixels[u].green);
1726 Minimize(result.blue, (*k)+k_pixels[u].blue);
1727 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1728 if ( image->colorspace == CMYKColorspace)
1729 Minimize(result.index, (*k)+k_indexes[u]);
1731 k_pixels += image->columns+kernel->width;
1732 k_indexes += image->columns+kernel->width;
1736 case UndefinedMorphology:
1738 break; /* Do nothing */
1740 /* Final mathematics of results (combine with original image?)
1742 ** NOTE: Difference Morphology operators Edge* and *Hat could also
1743 ** be done here but works better with iteration as a image difference
1744 ** in the controling function (below). Thicken and Thinning however
1745 ** should be done here so thay can be iterated correctly.
1748 case HitAndMissMorphology:
1749 case ErodeMorphology:
1750 result = min; /* minimum of neighbourhood */
1752 case DilateMorphology:
1753 result = max; /* maximum of neighbourhood */
1755 case ThinningMorphology:
1756 /* subtract pattern match from original */
1757 result.red -= min.red;
1758 result.green -= min.green;
1759 result.blue -= min.blue;
1760 result.opacity -= min.opacity;
1761 result.index -= min.index;
1763 case ThickenMorphology:
1764 /* Union with original image (maximize) - or should this be + */
1765 Maximize( result.red, min.red );
1766 Maximize( result.green, min.green );
1767 Maximize( result.blue, min.blue );
1768 Maximize( result.opacity, min.opacity );
1769 Maximize( result.index, min.index );
1772 /* result directly calculated or assigned */
1775 /* Assign the resulting pixel values - Clamping Result */
1777 case UndefinedMorphology:
1778 case DilateIntensityMorphology:
1779 case ErodeIntensityMorphology:
1780 break; /* full pixel was directly assigned - not a channel method */
1782 if ((channel & RedChannel) != 0)
1783 q->red = ClampToQuantum(result.red);
1784 if ((channel & GreenChannel) != 0)
1785 q->green = ClampToQuantum(result.green);
1786 if ((channel & BlueChannel) != 0)
1787 q->blue = ClampToQuantum(result.blue);
1788 if ((channel & OpacityChannel) != 0
1789 && image->matte == MagickTrue )
1790 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1791 if ((channel & IndexChannel) != 0
1792 && image->colorspace == CMYKColorspace)
1793 q_indexes[x] = ClampToQuantum(result.index);
1796 /* Count up changed pixels */
1797 if ( ( p[r].red != q->red )
1798 || ( p[r].green != q->green )
1799 || ( p[r].blue != q->blue )
1800 || ( p[r].opacity != q->opacity )
1801 || ( image->colorspace == CMYKColorspace &&
1802 p_indexes[r] != q_indexes[x] ) )
1803 changed++; /* The pixel had some value changed! */
1807 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1808 if (sync == MagickFalse)
1810 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1815 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1816 #pragma omp critical (MagickCore_MorphologyImage)
1818 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1819 if (proceed == MagickFalse)
1823 result_image->type=image->type;
1824 q_view=DestroyCacheView(q_view);
1825 p_view=DestroyCacheView(p_view);
1826 return(status ? (unsigned long) changed : 0);
1830 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1831 method, const long iterations,const KernelInfo *kernel, ExceptionInfo
1837 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1838 iterations,kernel,exception);
1839 return(morphology_image);
1843 MagickExport Image *MorphologyImageChannel(const Image *image,
1844 const ChannelType channel,const MorphologyMethod method,
1845 const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
1866 assert(image != (Image *) NULL);
1867 assert(image->signature == MagickSignature);
1868 assert(kernel != (KernelInfo *) NULL);
1869 assert(kernel->signature == MagickSignature);
1870 assert(exception != (ExceptionInfo *) NULL);
1871 assert(exception->signature == MagickSignature);
1873 if ( iterations == 0 )
1874 return((Image *)NULL); /* null operation - nothing to do! */
1876 /* kernel must be valid at this point
1877 * (except maybe for posible future morphology methods like "Prune"
1879 assert(kernel != (KernelInfo *)NULL);
1881 new_image = (Image *) NULL; /* for cleanup */
1882 old_image = (Image *) NULL;
1883 grad_image = (Image *) NULL;
1884 curr_kernel = (KernelInfo *) NULL;
1885 count = 0; /* interation count */
1886 changed = 1; /* was last run succesfull */
1887 curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
1888 curr_method = method; /* to be changed as nessary */
1890 limit = (unsigned long) iterations;
1891 if ( iterations < 0 )
1892 limit = image->columns > image->rows ? image->columns : image->rows;
1894 /* Third-level morphology methods */
1895 switch( curr_method ) {
1896 case EdgeMorphology:
1897 grad_image = MorphologyImageChannel(image, channel,
1898 DilateMorphology, iterations, curr_kernel, exception);
1899 if ( grad_image == (Image *) NULL )
1902 case EdgeInMorphology:
1903 curr_method = ErodeMorphology;
1905 case EdgeOutMorphology:
1906 curr_method = DilateMorphology;
1908 case TopHatMorphology:
1909 curr_method = OpenMorphology;
1911 case BottomHatMorphology:
1912 curr_method = CloseMorphology;
1915 break; /* not a third-level method */
1918 /* Second-level morphology methods */
1919 switch( curr_method ) {
1920 case OpenMorphology:
1921 /* Open is a Erode then a Dilate without reflection */
1922 new_image = MorphologyImageChannel(image, channel,
1923 ErodeMorphology, iterations, curr_kernel, exception);
1924 if (new_image == (Image *) NULL)
1926 curr_method = DilateMorphology;
1928 case OpenIntensityMorphology:
1929 new_image = MorphologyImageChannel(image, channel,
1930 ErodeIntensityMorphology, iterations, curr_kernel, exception);
1931 if (new_image == (Image *) NULL)
1933 curr_method = DilateIntensityMorphology;
1936 case CloseMorphology:
1937 /* Close is a Dilate then Erode using reflected kernel */
1938 /* A reflected kernel is needed for a Close */
1939 if ( curr_kernel == kernel )
1940 curr_kernel = CloneKernelInfo(kernel);
1941 if (curr_kernel == (KernelInfo *) NULL)
1943 RotateKernelInfo(curr_kernel,180);
1944 new_image = MorphologyImageChannel(image, channel,
1945 DilateMorphology, iterations, curr_kernel, exception);
1946 if (new_image == (Image *) NULL)
1948 curr_method = ErodeMorphology;
1950 case CloseIntensityMorphology:
1951 /* A reflected kernel is needed for a Close */
1952 if ( curr_kernel == kernel )
1953 curr_kernel = CloneKernelInfo(kernel);
1954 if (curr_kernel == (KernelInfo *) NULL)
1956 RotateKernelInfo(curr_kernel,180);
1957 new_image = MorphologyImageChannel(image, channel,
1958 DilateIntensityMorphology, iterations, curr_kernel, exception);
1959 if (new_image == (Image *) NULL)
1961 curr_method = ErodeIntensityMorphology;
1964 case CorrelateMorphology:
1965 /* A Correlation is actually a Convolution with a reflected kernel.
1966 ** However a Convolution is a weighted sum with a reflected kernel.
1967 ** It may seem stange to convert a Correlation into a Convolution
1968 ** as the Correleation is the simplier method, but Convolution is
1969 ** much more commonly used, and it makes sense to implement it directly
1970 ** so as to avoid the need to duplicate the kernel when it is not
1971 ** required (which is typically the default).
1973 if ( curr_kernel == kernel )
1974 curr_kernel = CloneKernelInfo(kernel);
1975 if (curr_kernel == (KernelInfo *) NULL)
1977 RotateKernelInfo(curr_kernel,180);
1978 curr_method = ConvolveMorphology;
1979 /* FALL-THRU into Correlate (weigthed sum without reflection) */
1981 case ConvolveMorphology:
1982 { /* Scale or Normalize kernel, according to user wishes
1983 ** before using it for the Convolve/Correlate method.
1985 ** FUTURE: provide some way for internal functions to disable
1986 ** user bias and scaling effects.
1989 *artifact = GetImageArtifact(image,"convolve:scale");
1990 if ( artifact != (char *)NULL ) {
1996 if ( curr_kernel == kernel )
1997 curr_kernel = CloneKernelInfo(kernel);
1998 if (curr_kernel == (KernelInfo *) NULL)
2001 flags = (GeometryFlags) ParseGeometry(artifact, &args);
2002 ScaleKernelInfo(curr_kernel, args.rho, flags);
2005 /* FALL-THRU to do the first, and typically the only iteration */
2008 /* Do a single iteration using the Low-Level Morphology method!
2009 ** This ensures a "new_image" has been generated, but allows us to skip
2010 ** the creation of 'old_image' if no more iterations are needed.
2012 ** The "curr_method" should also be set to a low-level method that is
2013 ** understood by the MorphologyPrimative() internal function.
2015 new_image=CloneImage(image,0,0,MagickTrue,exception);
2016 if (new_image == (Image *) NULL)
2018 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
2020 InheritException(exception,&new_image->exception);
2023 count++; /* iteration count */
2024 changed = MorphologyPrimative(image,new_image,curr_method,channel,
2025 curr_kernel, exception);
2026 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2027 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
2028 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2029 count, 0L, changed);
2033 * + "curr_method" should be set to a low-level morphology method
2034 * + "count=1" if the first iteration of the first kernel has been done.
2035 * + "new_image" is defined and contains the current resulting image
2038 /* The Hit-And-Miss Morphology is not 'iterated' against the resulting
2039 ** image from the previous morphology step. It must always be applied
2040 ** to the original image, with the results collected into a union (maximimum
2041 ** or lighten) of all the results. Multiple kernels may be applied but
2042 ** an iteration of the morphology does nothing, so is ignored.
2044 ** The first kernel is guranteed to have been done, so lets continue the
2045 ** process, with the other kernels in the kernel list.
2047 if ( method == HitAndMissMorphology )
2049 if ( curr_kernel->next != (KernelInfo *) NULL ) {
2050 /* create a second working image */
2051 old_image = CloneImage(image,0,0,MagickTrue,exception);
2052 if (old_image == (Image *) NULL)
2054 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
2056 InheritException(exception,&old_image->exception);
2060 /* loop through rest of the kernels */
2061 this_kernel=curr_kernel->next;
2063 while( this_kernel != (KernelInfo *) NULL )
2065 changed = MorphologyPrimative(image,old_image,curr_method,channel,
2066 this_kernel,exception);
2067 (void) CompositeImageChannel(new_image,
2068 (ChannelType) (channel & ~SyncChannels), LightenCompositeOp,
2070 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2071 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
2072 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2073 count, kernel_number, changed);
2074 this_kernel = this_kernel->next;
2077 old_image=DestroyImage(old_image);
2082 /* Repeat the low-level morphology over all kernels
2083 until iteration count limit or no change from any kernel is found */
2084 if ( ( count < limit && changed > 0 ) ||
2085 curr_kernel->next != (KernelInfo *) NULL ) {
2087 /* create a second working image */
2088 old_image = CloneImage(image,0,0,MagickTrue,exception);
2089 if (old_image == (Image *) NULL)
2091 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
2093 InheritException(exception,&old_image->exception);
2097 /* reset variables for the first/next iteration, or next kernel) */
2099 this_kernel = curr_kernel;
2100 total_changed = count != 0 ? changed : 0;
2101 if ( count != 0 && this_kernel != (KernelInfo *) NULL ) {
2102 count = 0; /* first iteration is not yet finished! */
2103 this_kernel = curr_kernel->next;
2105 total_changed = changed;
2108 while ( count < limit ) {
2110 while ( this_kernel != (KernelInfo *) NULL ) {
2111 Image *tmp = old_image;
2112 old_image = new_image;
2114 changed = MorphologyPrimative(old_image,new_image,curr_method,channel,
2115 this_kernel,exception);
2116 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2117 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
2118 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2119 count, kernel_number, changed);
2120 total_changed += changed;
2121 this_kernel = this_kernel->next;
2124 if ( total_changed == 0 )
2125 break; /* no changes after processing all kernels - ABORT */
2126 /* prepare for next loop */
2129 this_kernel = curr_kernel;
2131 old_image=DestroyImage(old_image);
2134 /* finished with kernel - destary any copy that was made */
2135 if ( curr_kernel != kernel )
2136 curr_kernel=DestroyKernelInfo(curr_kernel);
2138 /* Third-level Subtractive methods post-processing
2140 ** The removal of any 'Sync' channel flag in the Image Compositon below
2141 ** ensures the compose method is applied in a purely mathematical way, only
2142 ** the selected channels, without any normal 'alpha blending' normally
2143 ** associated with the compose method.
2145 ** Note "method" here is the 'original' morphological method, and not the
2146 ** 'current' morphological method used above to generate "new_image".
2149 case EdgeOutMorphology:
2150 case EdgeInMorphology:
2151 case TopHatMorphology:
2152 case BottomHatMorphology:
2153 /* Get Difference relative to the original image */
2154 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2155 DifferenceCompositeOp, image, 0, 0);
2157 case EdgeMorphology:
2158 /* Difference the Eroded image from the saved Dilated image */
2159 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2160 DifferenceCompositeOp, grad_image, 0, 0);
2161 grad_image=DestroyImage(grad_image);
2169 /* Yes goto's are bad, but in this case it makes cleanup lot more efficient */
2171 if ( new_image != (Image *) NULL )
2172 DestroyImage(new_image);
2174 if ( old_image != (Image *) NULL )
2175 DestroyImage(old_image);
2176 if ( curr_kernel != (KernelInfo *) NULL && curr_kernel != kernel )
2177 curr_kernel=DestroyKernelInfo(curr_kernel);
2182 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2186 + R o t a t e K e r n e l I n f o %
2190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2192 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
2193 % restricted to 90 degree angles, but this may be improved in the future.
2195 % The format of the RotateKernelInfo method is:
2197 % void RotateKernelInfo(KernelInfo *kernel, double angle)
2199 % A description of each parameter follows:
2201 % o kernel: the Morphology/Convolution kernel
2203 % o angle: angle to rotate in degrees
2205 % This function is only internel to this module, as it is not finalized,
2206 % especially with regard to non-orthogonal angles, and rotation of larger
2209 static void RotateKernelInfo(KernelInfo *kernel, double angle)
2211 /* angle the lower kernels first */
2212 if ( kernel->next != (KernelInfo *) NULL)
2213 RotateKernelInfo(kernel->next, angle);
2215 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
2217 ** TODO: expand beyond simple 90 degree rotates, flips and flops
2220 /* Modulus the angle */
2221 angle = fmod(angle, 360.0);
2225 if ( 337.5 < angle || angle <= 22.5 )
2226 return; /* no change! - At least at this time */
2228 switch (kernel->type) {
2229 /* These built-in kernels are cylindrical kernels, rotating is useless */
2230 case GaussianKernel:
2231 case LaplacianKernel:
2235 case ChebyshevKernel:
2236 case ManhattenKernel:
2237 case EuclideanKernel:
2240 /* These may be rotatable at non-90 angles in the future */
2241 /* but simply rotating them in multiples of 90 degrees is useless */
2247 /* These only allows a +/-90 degree rotation (by transpose) */
2248 /* A 180 degree rotation is useless */
2250 case RectangleKernel:
2251 if ( 135.0 < angle && angle <= 225.0 )
2253 if ( 225.0 < angle && angle <= 315.0 )
2257 /* these are freely rotatable in 90 degree units */
2260 case UndefinedKernel:
2261 case UserDefinedKernel:
2264 /* Attempt rotations by 45 degrees */
2265 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
2267 if ( kernel->width == 3 && kernel->height == 3 )
2268 { /* Rotate a 3x3 square by 45 degree angle */
2269 MagickRealType t = kernel->values[0];
2270 kernel->values[0] = kernel->values[1];
2271 kernel->values[1] = kernel->values[2];
2272 kernel->values[2] = kernel->values[5];
2273 kernel->values[5] = kernel->values[8];
2274 kernel->values[8] = kernel->values[7];
2275 kernel->values[7] = kernel->values[6];
2276 kernel->values[6] = kernel->values[3];
2277 kernel->values[3] = t;
2278 angle = fmod(angle+45.0, 360.0);
2281 perror("Unable to rotate non-3x3 kernel by 45 degrees");
2283 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
2285 if ( kernel->width == 1 || kernel->height == 1 )
2286 { /* Do a transpose of the image, which results in a 90
2287 ** degree rotation of a 1 dimentional kernel
2291 t = (long) kernel->width;
2292 kernel->width = kernel->height;
2293 kernel->height = (unsigned long) t;
2295 kernel->x = kernel->y;
2297 if ( kernel->width == 1 )
2298 angle = fmod(angle+270.0, 360.0); /* angle -90 degrees */
2300 angle = fmod(angle+90.0, 360.0); /* angle +90 degrees */
2302 else if ( kernel->width == kernel->height )
2303 { /* Rotate a square array of values by 90 degrees */
2304 register unsigned long
2306 register MagickRealType
2309 for( i=0, x=kernel->width-1; i<=x; i++, x--)
2310 for( j=0, y=kernel->height-1; j<y; j++, y--)
2311 { t = k[i+j*kernel->width];
2312 k[i+j*kernel->width] = k[y+i*kernel->width];
2313 k[y+i*kernel->width] = k[x+y*kernel->width];
2314 k[x+y*kernel->width] = k[j+x*kernel->width];
2315 k[j+x*kernel->width] = t;
2317 angle = fmod(angle+90.0, 360.0); /* angle +90 degrees */
2320 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
2322 if ( 135.0 < angle && angle <= 225.0 )
2324 /* For a 180 degree rotation - also know as a reflection */
2325 /* This is actually a very very common operation! */
2326 /* Basically all that is needed is a reversal of the kernel data! */
2333 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
2334 t=k[i], k[i]=k[j], k[j]=t;
2336 kernel->x = (long) kernel->width - kernel->x - 1;
2337 kernel->y = (long) kernel->height - kernel->y - 1;
2338 angle = fmod(angle+180.0, 360.0); /* angle+180 degrees */
2340 /* At this point angle should at least between -45 (315) and +45 degrees
2341 * In the future some form of non-orthogonal angled rotates could be
2342 * performed here, posibily with a linear kernel restriction.
2346 { /* Do a Flop by reversing each row.
2355 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2356 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2357 t=k[x], k[x]=k[r], k[r]=t;
2359 kernel->x = kernel->width - kernel->x - 1;
2360 angle = fmod(angle+180.0, 360.0);
2367 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2371 % S c a l e K e r n e l I n f o %
2375 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2377 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
2378 % without normalization of the sum of the kernel values (as per given flags).
2380 % By default (no flags given) the values within the kernel is scaled
2381 % directly using given scaling factor without change.
2383 % If any 'normalize_flags' are given the kernel will first be normalized and
2384 % then further scaled by the scaling factor value given. A 'PercentValue'
2385 % flag will cause the given scaling factor to be divided by one hundred
2388 % Kernel normalization ('normalize_flags' given) is designed to ensure that
2389 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
2390 % morphology methods will fall into -1.0 to +1.0 range. Note that for
2391 % non-HDRI versions of IM this may cause images to have any negative results
2392 % clipped, unless some 'bias' is used.
2394 % More specifically. Kernels which only contain positive values (such as a
2395 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
2396 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
2398 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
2399 % the kernel will be scaled by the absolute of the sum of kernel values, so
2400 % that it will generally fall within the +/- 1.0 range.
2402 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
2403 % will be scaled by just the sum of the postive values, so that its output
2404 % range will again fall into the +/- 1.0 range.
2406 % For special kernels designed for locating shapes using 'Correlate', (often
2407 % only containing +1 and -1 values, representing foreground/brackground
2408 % matching) a special normalization method is provided to scale the positive
2409 % values seperatally to those of the negative values, so the kernel will be
2410 % forced to become a zero-sum kernel better suited to such searches.
2412 % WARNING: Correct normalization of the kernel assumes that the '*_range'
2413 % attributes within the kernel structure have been correctly set during the
2416 % NOTE: The values used for 'normalize_flags' have been selected specifically
2417 % to match the use of geometry options, so that '!' means NormalizeValue,
2418 % '^' means CorrelateNormalizeValue, and '%' means PercentValue. All other
2419 % GeometryFlags values are ignored.
2421 % The format of the ScaleKernelInfo method is:
2423 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
2424 % const MagickStatusType normalize_flags )
2426 % A description of each parameter follows:
2428 % o kernel: the Morphology/Convolution kernel
2431 % multiply all values (after normalization) by this factor if not
2432 % zero. If the kernel is normalized regardless of any flags.
2434 % o normalize_flags:
2435 % GeometryFlags defining normalization method to use.
2436 % specifically: NormalizeValue, CorrelateNormalizeValue,
2437 % and/or PercentValue
2439 % This function is internal to this module only at this time, but can be
2440 % exported to other modules if needed.
2442 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
2443 const double scaling_factor,const GeometryFlags normalize_flags)
2452 /* scale the lower kernels first */
2453 if ( kernel->next != (KernelInfo *) NULL)
2454 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
2457 if ( (normalize_flags&NormalizeValue) != 0 ) {
2458 /* normalize kernel appropriately */
2459 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
2460 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
2462 pos_scale = kernel->positive_range; /* special zero-summing kernel */
2464 /* force kernel into being a normalized zero-summing kernel */
2465 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
2466 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
2467 ? kernel->positive_range : 1.0;
2468 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
2469 ? -kernel->negative_range : 1.0;
2472 neg_scale = pos_scale;
2474 /* finialize scaling_factor for positive and negative components */
2475 pos_scale = scaling_factor/pos_scale;
2476 neg_scale = scaling_factor/neg_scale;
2477 if ( (normalize_flags&PercentValue) != 0 ) {
2482 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2483 if ( ! IsNan(kernel->values[i]) )
2484 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
2486 /* convolution output range */
2487 kernel->positive_range *= pos_scale;
2488 kernel->negative_range *= neg_scale;
2489 /* maximum and minimum values in kernel */
2490 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
2491 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
2493 /* swap kernel settings if user scaling factor is negative */
2494 if ( scaling_factor < MagickEpsilon ) {
2496 t = kernel->positive_range;
2497 kernel->positive_range = kernel->negative_range;
2498 kernel->negative_range = t;
2499 t = kernel->maximum;
2500 kernel->maximum = kernel->minimum;
2501 kernel->minimum = 1;
2508 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2512 + S h o w K e r n e l I n f o %
2516 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2518 % ShowKernelInfo() outputs the details of the given kernel defination to
2519 % standard error, generally due to a users 'showkernel' option request.
2521 % The format of the ShowKernel method is:
2523 % void ShowKernelInfo(KernelInfo *kernel)
2525 % A description of each parameter follows:
2527 % o kernel: the Morphology/Convolution kernel
2529 % This function is internal to this module only at this time. That may change
2532 MagickExport void ShowKernelInfo(KernelInfo *kernel)
2540 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
2542 fprintf(stderr, "Kernel ");
2543 if ( kernel->next != (KernelInfo *) NULL )
2544 fprintf(stderr, " #%ld", c );
2545 fprintf(stderr, " \"%s\" of size %lux%lu%+ld%+ld ",
2546 MagickOptionToMnemonic(MagickKernelOptions, k->type),
2547 kernel->width, k->height,
2548 kernel->x, kernel->y );
2550 " with values from %.*lg to %.*lg\n",
2551 GetMagickPrecision(), k->minimum,
2552 GetMagickPrecision(), k->maximum);
2553 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2554 GetMagickPrecision(), k->negative_range,
2555 GetMagickPrecision(), k->positive_range,
2556 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2557 for (i=v=0; v < (long) k->height; v++) {
2558 fprintf(stderr,"%2ld:",v);
2559 for (u=0; u < (long) k->width; u++, i++)
2560 if ( IsNan(k->values[i]) )
2561 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2563 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2564 GetMagickPrecision(), k->values[i]);
2565 fprintf(stderr,"\n");
2571 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2575 + Z e r o K e r n e l N a n s %
2579 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2581 % ZeroKernelNans() replaces any special 'nan' value that may be present in
2582 % the kernel with a zero value. This is typically done when the kernel will
2583 % be used in special hardware (GPU) convolution processors, to simply
2586 % The format of the ZeroKernelNans method is:
2588 % voidZeroKernelNans (KernelInfo *kernel)
2590 % A description of each parameter follows:
2592 % o kernel: the Morphology/Convolution kernel
2594 % FUTURE: return the information in a string for API usage.
2596 MagickExport void ZeroKernelNans(KernelInfo *kernel)
2601 /* scale the lower kernels first */
2602 if ( kernel->next != (KernelInfo *) NULL)
2603 ZeroKernelNans(kernel->next);
2605 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2606 if ( IsNan(kernel->values[i]) )
2607 kernel->values[i] = 0.0;