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/morphology-private.h"
69 #include "magick/option.h"
70 #include "magick/pixel-private.h"
71 #include "magick/prepress.h"
72 #include "magick/quantize.h"
73 #include "magick/registry.h"
74 #include "magick/semaphore.h"
75 #include "magick/splay-tree.h"
76 #include "magick/statistic.h"
77 #include "magick/string_.h"
78 #include "magick/string-private.h"
79 #include "magick/token.h"
83 ** The following test is for special floating point numbers of value NaN (not
84 ** a number), that may be used within a Kernel Definition. NaN's are defined
85 ** as part of the IEEE standard for floating point number representation.
87 ** These are used as a Kernel value to mean that this kernel position is not
88 ** part of the kernel neighbourhood for convolution or morphology processing,
89 ** and thus should be ignored. This allows the use of 'shaped' kernels.
91 ** The special properity that two NaN's are never equal, even if they are from
92 ** the same variable allow you to test if a value is special NaN value.
94 ** This macro IsNaN() is thus is only true if the value given is NaN.
96 #define IsNan(a) ((a)!=(a))
99 Other global definitions used by module.
101 static inline double MagickMin(const double x,const double y)
103 return( x < y ? x : y);
105 static inline double MagickMax(const double x,const double y)
107 return( x > y ? x : y);
109 #define Minimize(assign,value) assign=MagickMin(assign,value)
110 #define Maximize(assign,value) assign=MagickMax(assign,value)
112 /* Currently these are only internal to this module */
114 CalcKernelMetaData(KernelInfo *),
115 ExpandMirrorKernelInfo(KernelInfo *),
116 ExpandRotateKernelInfo(KernelInfo *, const double),
117 RotateKernelInfo(KernelInfo *, double);
120 /* Quick function to find last kernel in a kernel list */
121 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
123 while (kernel->next != (KernelInfo *) NULL)
124 kernel = kernel->next;
130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134 % A c q u i r e K e r n e l I n f o %
138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
140 % AcquireKernelInfo() takes the given string (generally supplied by the
141 % user) and converts it into a Morphology/Convolution Kernel. This allows
142 % users to specify a kernel from a number of pre-defined kernels, or to fully
143 % specify their own kernel for a specific Convolution or Morphology
146 % The kernel so generated can be any rectangular array of floating point
147 % values (doubles) with the 'control point' or 'pixel being affected'
148 % anywhere within that array of values.
150 % Previously IM was restricted to a square of odd size using the exact
151 % center as origin, this is no longer the case, and any rectangular kernel
152 % with any value being declared the origin. This in turn allows the use of
153 % highly asymmetrical kernels.
155 % The floating point values in the kernel can also include a special value
156 % known as 'nan' or 'not a number' to indicate that this value is not part
157 % of the kernel array. This allows you to shaped the kernel within its
158 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
159 % shape. However at least one non-nan value must be provided for correct
160 % working of a kernel.
162 % The returned kernel should be freed using the DestroyKernelInfo() when you
163 % are finished with it. Do not free this memory yourself.
165 % Input kernel defintion strings can consist of any of three types.
168 % Select from one of the built in kernels, using the name and
169 % geometry arguments supplied. See AcquireKernelBuiltIn()
171 % "WxH[+X+Y][@><]:num, num, num ..."
172 % a kernel of size W by H, with W*H floating point numbers following.
173 % the 'center' can be optionally be defined at +X+Y (such that +0+0
174 % is top left corner). If not defined the pixel in the center, for
175 % odd sizes, or to the immediate top or left of center for even sizes
176 % is automatically selected.
178 % "num, num, num, num, ..."
179 % list of floating point numbers defining an 'old style' odd sized
180 % square kernel. At least 9 values should be provided for a 3x3
181 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
182 % Values can be space or comma separated. This is not recommended.
184 % You can define a 'list of kernels' which can be used by some morphology
185 % operators A list is defined as a semi-colon seperated list kernels.
187 % " kernel ; kernel ; kernel ; "
189 % Any extra ';' characters, at start, end or between kernel defintions are
192 % The special flags will expand a single kernel, into a list of rotated
193 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
194 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
195 % The '<' also exands using 90-degree rotates, but giving a 180-degree
196 % reflected kernel before the +/- 90-degree rotations, which can be important
197 % for Thinning operations.
199 % Note that 'name' kernels will start with an alphabetic character while the
200 % new kernel specification has a ':' character in its specification string.
201 % If neither is the case, it is assumed an old style of a simple list of
202 % numbers generating a odd-sized square kernel has been given.
204 % The format of the AcquireKernal method is:
206 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
208 % A description of each parameter follows:
210 % o kernel_string: the Morphology/Convolution kernel wanted.
214 /* This was separated so that it could be used as a separate
215 ** array input handling function, such as for -color-matrix
217 static KernelInfo *ParseKernelArray(const char *kernel_string)
223 token[MaxTextExtent];
233 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
241 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
242 if (kernel == (KernelInfo *)NULL)
244 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
245 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
246 kernel->negative_range = kernel->positive_range = 0.0;
247 kernel->type = UserDefinedKernel;
248 kernel->next = (KernelInfo *) NULL;
249 kernel->signature = MagickSignature;
251 /* find end of this specific kernel definition string */
252 end = strchr(kernel_string, ';');
253 if ( end == (char *) NULL )
254 end = strchr(kernel_string, '\0');
256 /* clear flags - for Expanding kernal lists thorugh rotations */
259 /* Has a ':' in argument - New user kernel specification */
260 p = strchr(kernel_string, ':');
261 if ( p != (char *) NULL && p < end)
263 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
264 memcpy(token, kernel_string, (size_t) (p-kernel_string));
265 token[p-kernel_string] = '\0';
266 SetGeometryInfo(&args);
267 flags = ParseGeometry(token, &args);
269 /* Size handling and checks of geometry settings */
270 if ( (flags & WidthValue) == 0 ) /* if no width then */
271 args.rho = args.sigma; /* then width = height */
272 if ( args.rho < 1.0 ) /* if width too small */
273 args.rho = 1.0; /* then width = 1 */
274 if ( args.sigma < 1.0 ) /* if height too small */
275 args.sigma = args.rho; /* then height = width */
276 kernel->width = (size_t)args.rho;
277 kernel->height = (size_t)args.sigma;
279 /* Offset Handling and Checks */
280 if ( args.xi < 0.0 || args.psi < 0.0 )
281 return(DestroyKernelInfo(kernel));
282 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
283 : (ssize_t) (kernel->width-1)/2;
284 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
285 : (ssize_t) (kernel->height-1)/2;
286 if ( kernel->x >= (ssize_t) kernel->width ||
287 kernel->y >= (ssize_t) kernel->height )
288 return(DestroyKernelInfo(kernel));
290 p++; /* advance beyond the ':' */
293 { /* ELSE - Old old specification, forming odd-square kernel */
294 /* count up number of values given */
295 p=(const char *) kernel_string;
296 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
297 p++; /* ignore "'" chars for convolve filter usage - Cristy */
298 for (i=0; p < end; i++)
300 GetMagickToken(p,&p,token);
302 GetMagickToken(p,&p,token);
304 /* set the size of the kernel - old sized square */
305 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
306 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
307 p=(const char *) kernel_string;
308 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
309 p++; /* ignore "'" chars for convolve filter usage - Cristy */
312 /* Read in the kernel values from rest of input string argument */
313 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
314 kernel->height*sizeof(double));
315 if (kernel->values == (double *) NULL)
316 return(DestroyKernelInfo(kernel));
318 kernel->minimum = +MagickHuge;
319 kernel->maximum = -MagickHuge;
320 kernel->negative_range = kernel->positive_range = 0.0;
322 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
324 GetMagickToken(p,&p,token);
326 GetMagickToken(p,&p,token);
327 if ( LocaleCompare("nan",token) == 0
328 || LocaleCompare("-",token) == 0 ) {
329 kernel->values[i] = nan; /* do not include this value in kernel */
332 kernel->values[i] = StringToDouble(token);
333 ( kernel->values[i] < 0)
334 ? ( kernel->negative_range += kernel->values[i] )
335 : ( kernel->positive_range += kernel->values[i] );
336 Minimize(kernel->minimum, kernel->values[i]);
337 Maximize(kernel->maximum, kernel->values[i]);
341 /* sanity check -- no more values in kernel definition */
342 GetMagickToken(p,&p,token);
343 if ( *token != '\0' && *token != ';' && *token != '\'' )
344 return(DestroyKernelInfo(kernel));
347 /* this was the old method of handling a incomplete kernel */
348 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
349 Minimize(kernel->minimum, kernel->values[i]);
350 Maximize(kernel->maximum, kernel->values[i]);
351 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
352 kernel->values[i]=0.0;
355 /* Number of values for kernel was not enough - Report Error */
356 if ( i < (ssize_t) (kernel->width*kernel->height) )
357 return(DestroyKernelInfo(kernel));
360 /* check that we recieved at least one real (non-nan) value! */
361 if ( kernel->minimum == MagickHuge )
362 return(DestroyKernelInfo(kernel));
364 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
365 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
366 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
367 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
368 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
369 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
374 static KernelInfo *ParseKernelName(const char *kernel_string)
380 token[MaxTextExtent];
395 /* Parse special 'named' kernel */
396 GetMagickToken(kernel_string,&p,token);
397 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
398 if ( type < 0 || type == UserDefinedKernel )
399 return((KernelInfo *)NULL); /* not a valid named kernel */
401 while (((isspace((int) ((unsigned char) *p)) != 0) ||
402 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
405 end = strchr(p, ';'); /* end of this kernel defintion */
406 if ( end == (char *) NULL )
407 end = strchr(p, '\0');
409 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
410 memcpy(token, p, (size_t) (end-p));
412 SetGeometryInfo(&args);
413 flags = ParseGeometry(token, &args);
416 /* For Debugging Geometry Input */
417 fprintf(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
418 flags, args.rho, args.sigma, args.xi, args.psi );
421 /* special handling of missing values in input string */
423 case RectangleKernel:
424 if ( (flags & WidthValue) == 0 ) /* if no width then */
425 args.rho = args.sigma; /* then width = height */
426 if ( args.rho < 1.0 ) /* if width too small */
427 args.rho = 3; /* then width = 3 */
428 if ( args.sigma < 1.0 ) /* if height too small */
429 args.sigma = args.rho; /* then height = width */
430 if ( (flags & XValue) == 0 ) /* center offset if not defined */
431 args.xi = (double)(((ssize_t)args.rho-1)/2);
432 if ( (flags & YValue) == 0 )
433 args.psi = (double)(((ssize_t)args.sigma-1)/2);
440 /* If no scale given (a 0 scale is valid! - set it to 1.0 */
441 if ( (flags & HeightValue) == 0 )
445 if ( (flags & XValue) == 0 )
448 case ChebyshevKernel:
449 case ManhattanKernel:
450 case EuclideanKernel:
451 if ( (flags & HeightValue) == 0 ) /* no distance scale */
452 args.sigma = 100.0; /* default distance scaling */
453 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
454 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
455 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
456 args.sigma *= QuantumRange/100.0; /* percentage of color range */
462 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
464 /* global expand to rotated kernel list - only for single kernels */
465 if ( kernel->next == (KernelInfo *) NULL ) {
466 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
467 ExpandRotateKernelInfo(kernel, 45.0);
468 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
469 ExpandRotateKernelInfo(kernel, 90.0);
470 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
471 ExpandMirrorKernelInfo(kernel);
477 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
485 token[MaxTextExtent];
497 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
499 /* ignore extra or multiple ';' kernel seperators */
500 if ( *token != ';' ) {
502 /* tokens starting with alpha is a Named kernel */
503 if (isalpha((int) *token) != 0)
504 new_kernel = ParseKernelName(p);
505 else /* otherwise a user defined kernel array */
506 new_kernel = ParseKernelArray(p);
508 /* Error handling -- this is not proper error handling! */
509 if ( new_kernel == (KernelInfo *) NULL ) {
510 fprintf(stderr, "Failed to parse kernel number #%.20g\n",(double)
512 if ( kernel != (KernelInfo *) NULL )
513 kernel=DestroyKernelInfo(kernel);
514 return((KernelInfo *) NULL);
517 /* initialise or append the kernel list */
518 if ( kernel == (KernelInfo *) NULL )
521 LastKernelInfo(kernel)->next = new_kernel;
524 /* look for the next kernel in list */
526 if ( p == (char *) NULL )
536 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
540 % A c q u i r e K e r n e l B u i l t I n %
544 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
546 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
547 % kernels used for special purposes such as gaussian blurring, skeleton
548 % pruning, and edge distance determination.
550 % They take a KernelType, and a set of geometry style arguments, which were
551 % typically decoded from a user supplied string, or from a more complex
552 % Morphology Method that was requested.
554 % The format of the AcquireKernalBuiltIn method is:
556 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
557 % const GeometryInfo args)
559 % A description of each parameter follows:
561 % o type: the pre-defined type of kernel wanted
563 % o args: arguments defining or modifying the kernel
565 % Convolution Kernels
568 % the No-Op kernel, also requivelent to Gaussian of sigma zero.
569 % Basically a 3x3 kernel of a 1 surrounded by zeros.
571 % Gaussian:{radius},{sigma}
572 % Generate a two-dimentional gaussian kernel, as used by -gaussian.
573 % The sigma for the curve is required. The resulting kernel is
576 % If 'sigma' is zero, you get a single pixel on a field of zeros.
578 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
579 % the final size of the resulting kernel to a square 2*radius+1 in size.
580 % The radius should be at least 2 times that of the sigma value, or
581 % sever clipping and aliasing may result. If not given or set to 0 the
582 % radius will be determined so as to produce the best minimal error
583 % result, which is usally much larger than is normally needed.
585 % LoG:{radius},{sigma}
586 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
587 % The supposed ideal edge detection, zero-summing kernel.
589 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
590 % approx 1.6 (according to wikipedia).
592 % DoG:{radius},{sigma1},{sigma2}
593 % "Difference of Gaussians" Kernel.
594 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
595 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
596 % The result is a zero-summing kernel.
598 % Blur:{radius},{sigma}[,{angle}]
599 % Generates a 1 dimensional or linear gaussian blur, at the angle given
600 % (current restricted to orthogonal angles). If a 'radius' is given the
601 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
602 % by a 90 degree angle.
604 % If 'sigma' is zero, you get a single pixel on a field of zeros.
606 % Note that two convolutions with two "Blur" kernels perpendicular to
607 % each other, is equivelent to a far larger "Gaussian" kernel with the
608 % same sigma value, However it is much faster to apply. This is how the
609 % "-blur" operator actually works.
611 % Comet:{width},{sigma},{angle}
612 % Blur in one direction only, much like how a bright object leaves
613 % a comet like trail. The Kernel is actually half a gaussian curve,
614 % Adding two such blurs in opposite directions produces a Blur Kernel.
615 % Angle can be rotated in multiples of 90 degrees.
617 % Note that the first argument is the width of the kernel and not the
618 % radius of the kernel.
620 % # Still to be implemented...
624 % # Set kernel values using a resize filter, and given scale (sigma)
625 % # Cylindrical or Linear. Is this posible with an image?
628 % Named Constant Convolution Kernels
630 % All these are unscaled, zero-summing kernels by default. As such for
631 % non-HDRI version of ImageMagick some form of normalization, user scaling,
632 % and biasing the results is recommended, to prevent the resulting image
635 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
636 % 45 degrees to generate the 8 angled varients of each of the kernels.
639 % Discrete Lapacian Kernels, (without normalization)
640 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
641 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
642 % Type 2 : 3x3 with center:4 edge:1 corner:-2
643 % Type 3 : 3x3 with center:4 edge:-2 corner:1
644 % Type 5 : 5x5 laplacian
645 % Type 7 : 7x7 laplacian
646 % Type 15 : 5x5 LoG (sigma approx 1.4)
647 % Type 19 : 9x9 LoG (sigma approx 1.4)
650 % Sobel 'Edge' convolution kernel (3x3)
655 % Sobel:{type},{angle}
656 % Type 0: default un-nomalized version shown above.
658 % Type 1: As default but pre-normalized
663 % Type 2: Diagonal version with same normalization as 1
669 % Roberts convolution kernel (3x3)
675 % Prewitt Edge convolution kernel (3x3)
681 % Prewitt's "Compass" convolution kernel (3x3)
687 % Kirsch's "Compass" convolution kernel (3x3)
693 % Frei-Chen Edge Detector is based on a kernel that is similar to
694 % the Sobel Kernel, but is designed to be isotropic. That is it takes
695 % into account the distance of the diagonal in the kernel.
698 % | sqrt(2), 0, -sqrt(2) |
701 % FreiChen:{type},{angle}
703 % Frei-Chen Pre-weighted kernels...
705 % Type 0: default un-nomalized version shown above.
707 % Type 1: Orthogonal Kernel (same as type 11 below)
709 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
712 % Type 2: Diagonal form of Kernel...
714 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
717 % However this kernel is als at the heart of the FreiChen Edge Detection
718 % Process which uses a set of 9 specially weighted kernel. These 9
719 % kernels not be normalized, but directly applied to the image. The
720 % results is then added together, to produce the intensity of an edge in
721 % a specific direction. The square root of the pixel value can then be
722 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
723 % from each other, both the direction and the strength of the edge can be
726 % Type 10: All 9 of the following pre-weighted kernels...
728 % Type 11: | 1, 0, -1 |
729 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
732 % Type 12: | 1, sqrt(2), 1 |
733 % | 0, 0, 0 | / 2*sqrt(2)
736 % Type 13: | sqrt(2), -1, 0 |
737 % | -1, 0, 1 | / 2*sqrt(2)
740 % Type 14: | 0, 1, -sqrt(2) |
741 % | -1, 0, 1 | / 2*sqrt(2)
744 % Type 15: | 0, -1, 0 |
748 % Type 16: | 1, 0, -1 |
752 % Type 17: | 1, -2, 1 |
756 % Type 18: | -2, 1, -2 |
760 % Type 19: | 1, 1, 1 |
764 % The first 4 are for edge detection, the next 4 are for line detection
765 % and the last is to add a average component to the results.
767 % Using a special type of '-1' will return all 9 pre-weighted kernels
768 % as a multi-kernel list, so that you can use them directly (without
769 % normalization) with the special "-set option:morphology:compose Plus"
770 % setting to apply the full FreiChen Edge Detection Technique.
772 % If 'type' is large it will be taken to be an actual rotation angle for
773 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
774 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
776 % WARNING: The above was layed out as per
777 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
778 % But rotated 90 degrees so direction is from left rather than the top.
779 % I have yet to find any secondary confirmation of the above. The only
780 % other source found was actual source code at
781 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
782 % Neigher paper defineds the kernels in a way that looks locical or
783 % correct when taken as a whole.
787 % Diamond:[{radius}[,{scale}]]
788 % Generate a diamond shaped kernel with given radius to the points.
789 % Kernel size will again be radius*2+1 square and defaults to radius 1,
790 % generating a 3x3 kernel that is slightly larger than a square.
792 % Square:[{radius}[,{scale}]]
793 % Generate a square shaped kernel of size radius*2+1, and defaulting
794 % to a 3x3 (radius 1).
796 % Note that using a larger radius for the "Square" or the "Diamond" is
797 % also equivelent to iterating the basic morphological method that many
798 % times. However iterating with the smaller radius is actually faster
799 % than using a larger kernel radius.
801 % Rectangle:{geometry}
802 % Simply generate a rectangle of 1's with the size given. You can also
803 % specify the location of the 'control point', otherwise the closest
804 % pixel to the center of the rectangle is selected.
806 % Properly centered and odd sized rectangles work the best.
808 % Disk:[{radius}[,{scale}]]
809 % Generate a binary disk of the radius given, radius may be a float.
810 % Kernel size will be ceil(radius)*2+1 square.
811 % NOTE: Here are some disk shapes of specific interest
812 % "Disk:1" => "diamond" or "cross:1"
813 % "Disk:1.5" => "square"
814 % "Disk:2" => "diamond:2"
815 % "Disk:2.5" => a general disk shape of radius 2
816 % "Disk:2.9" => "square:2"
817 % "Disk:3.5" => default - octagonal/disk shape of radius 3
818 % "Disk:4.2" => roughly octagonal shape of radius 4
819 % "Disk:4.3" => a general disk shape of radius 4
820 % After this all the kernel shape becomes more and more circular.
822 % Because a "disk" is more circular when using a larger radius, using a
823 % larger radius is preferred over iterating the morphological operation.
825 % Symbol Dilation Kernels
827 % These kernel is not a good general morphological kernel, but is used
828 % more for highlighting and marking any single pixels in an image using,
829 % a "Dilate" method as appropriate.
831 % For the same reasons iterating these kernels does not produce the
832 % same result as using a larger radius for the symbol.
834 % Plus:[{radius}[,{scale}]]
835 % Cross:[{radius}[,{scale}]]
836 % Generate a kernel in the shape of a 'plus' or a 'cross' with
837 % a each arm the length of the given radius (default 2).
839 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
841 % Ring:{radius1},{radius2}[,{scale}]
842 % A ring of the values given that falls between the two radii.
843 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
844 % This is the 'edge' pixels of the default "Disk" kernel,
845 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
847 % Hit and Miss Kernels
849 % Peak:radius1,radius2
850 % Find any peak larger than the pixels the fall between the two radii.
851 % The default ring of pixels is as per "Ring".
853 % Find flat orthogonal edges of a binary shape
855 % Find 90 degree corners of a binary shape
857 % Find end points of lines (for pruning a skeletion)
858 % Two types of lines ends (default to both) can be searched for
859 % Type 0: All line ends
860 % Type 1: single kernel for 4-conneected line ends
861 % Type 2: single kernel for simple line ends
863 % Find three line junctions (within a skeletion)
864 % Type 0: all line junctions
865 % Type 1: Y Junction kernel
866 % Type 2: Diagonal T Junction kernel
867 % Type 3: Orthogonal T Junction kernel
868 % Type 4: Diagonal X Junction kernel
869 % Type 5: Orthogonal + Junction kernel
871 % Find single pixel ridges or thin lines
872 % Type 1: Fine single pixel thick lines and ridges
873 % Type 2: Find two pixel thick lines and ridges
875 % Octagonal thicken kernel, to generate convex hulls of 45 degrees
877 % Traditional skeleton generating kernels.
878 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
879 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
880 % Type 3: Experimental Variation to try to present left-right symmetry
881 % Type 4: Experimental Variation to preserve left-right symmetry
883 % Distance Measuring Kernels
885 % Different types of distance measuring methods, which are used with the
886 % a 'Distance' morphology method for generating a gradient based on
887 % distance from an edge of a binary shape, though there is a technique
888 % for handling a anti-aliased shape.
890 % See the 'Distance' Morphological Method, for information of how it is
893 % Chebyshev:[{radius}][x{scale}[%!]]
894 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
895 % one to any neighbour, orthogonal or diagonal. One why of thinking of
896 % it is the number of squares a 'King' or 'Queen' in chess needs to
897 % traverse reach any other position on a chess board. It results in a
898 % 'square' like distance function, but one where diagonals are closer
901 % Manhattan:[{radius}][x{scale}[%!]]
902 % Manhattan Distance (also known as Rectilinear Distance, or the Taxi
903 % Cab metric), is the distance needed when you can only travel in
904 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
905 % in chess would travel. It results in a diamond like distances, where
906 % diagonals are further than expected.
908 % Euclidean:[{radius}][x{scale}[%!]]
909 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
910 % However by default the kernel size only has a radius of 1, which
911 % limits the distance to 'Knight' like moves, with only orthogonal and
912 % diagonal measurements being correct. As such for the default kernel
913 % you will get octagonal like distance function, which is reasonally
916 % However if you use a larger radius such as "Euclidean:4" you will
917 % get a much smoother distance gradient from the edge of the shape.
918 % Of course a larger kernel is slower to use, and generally not needed.
920 % To allow the use of fractional distances that you get with diagonals
921 % the actual distance is scaled by a fixed value which the user can
922 % provide. This is not actually nessary for either ""Chebyshev" or
923 % "Manhattan" distance kernels, but is done for all three distance
924 % kernels. If no scale is provided it is set to a value of 100,
925 % allowing for a maximum distance measurement of 655 pixels using a Q16
926 % version of IM, from any edge. However for small images this can
927 % result in quite a dark gradient.
931 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
932 const GeometryInfo *args)
945 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
947 /* Generate a new empty kernel if needed */
948 kernel=(KernelInfo *) NULL;
950 case UndefinedKernel: /* These should not call this function */
951 case UserDefinedKernel:
953 case UnityKernel: /* Named Descrete Convolution Kernels */
954 case LaplacianKernel:
961 case EdgesKernel: /* Hit and Miss kernels */
963 case ThinDiagonalsKernel:
965 case LineJunctionsKernel:
967 case ConvexHullKernel:
969 break; /* A pre-generated kernel is not needed */
971 /* set to 1 to do a compile-time check that we haven't missed anything */
979 case RectangleKernel:
985 case ChebyshevKernel:
986 case ManhattanKernel:
987 case EuclideanKernel:
991 /* Generate the base Kernel Structure */
992 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
993 if (kernel == (KernelInfo *) NULL)
995 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
996 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
997 kernel->negative_range = kernel->positive_range = 0.0;
999 kernel->next = (KernelInfo *) NULL;
1000 kernel->signature = MagickSignature;
1005 /* Convolution Kernels */
1006 case GaussianKernel:
1010 sigma = fabs(args->sigma),
1011 sigma2 = fabs(args->xi),
1014 if ( args->rho >= 1.0 )
1015 kernel->width = (size_t)args->rho*2+1;
1016 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1017 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1019 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1020 kernel->height = kernel->width;
1021 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1022 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1023 kernel->height*sizeof(double));
1024 if (kernel->values == (double *) NULL)
1025 return(DestroyKernelInfo(kernel));
1027 /* WARNING: The following generates a 'sampled gaussian' kernel.
1028 * What we really want is a 'discrete gaussian' kernel.
1030 * How to do this is currently not known, but appears to be
1031 * basied on the Error Function 'erf()' (intergral of a gaussian)
1034 if ( type == GaussianKernel || type == DoGKernel )
1035 { /* Calculate a Gaussian, OR positive half of a DoG */
1036 if ( sigma > MagickEpsilon )
1037 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1038 B = 1.0/(Magick2PI*sigma*sigma);
1039 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1040 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1041 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1043 else /* limiting case - a unity (normalized Dirac) kernel */
1044 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1045 kernel->width*kernel->height*sizeof(double));
1046 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1050 if ( type == DoGKernel )
1051 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1052 if ( sigma2 > MagickEpsilon )
1053 { sigma = sigma2; /* simplify loop expressions */
1054 A = 1.0/(2.0*sigma*sigma);
1055 B = 1.0/(Magick2PI*sigma*sigma);
1056 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1057 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1058 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1060 else /* limiting case - a unity (normalized Dirac) kernel */
1061 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1064 if ( type == LoGKernel )
1065 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1066 if ( sigma > MagickEpsilon )
1067 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1068 B = 1.0/(MagickPI*sigma*sigma*sigma*sigma);
1069 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1070 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1071 { R = ((double)(u*u+v*v))*A;
1072 kernel->values[i] = (1-R)*exp(-R)*B;
1075 else /* special case - generate a unity kernel */
1076 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1077 kernel->width*kernel->height*sizeof(double));
1078 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1082 /* Note the above kernels may have been 'clipped' by a user defined
1083 ** radius, producing a smaller (darker) kernel. Also for very small
1084 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1085 ** producing a very bright kernel.
1087 ** Normalization will still be needed.
1090 /* Normalize the 2D Gaussian Kernel
1092 ** NB: a CorrelateNormalize performs a normal Normalize if
1093 ** there are no negative values.
1095 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1096 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1102 sigma = fabs(args->sigma),
1105 if ( args->rho >= 1.0 )
1106 kernel->width = (size_t)args->rho*2+1;
1108 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1110 kernel->x = (ssize_t) (kernel->width-1)/2;
1112 kernel->negative_range = kernel->positive_range = 0.0;
1113 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1114 kernel->height*sizeof(double));
1115 if (kernel->values == (double *) NULL)
1116 return(DestroyKernelInfo(kernel));
1119 #define KernelRank 3
1120 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1121 ** It generates a gaussian 3 times the width, and compresses it into
1122 ** the expected range. This produces a closer normalization of the
1123 ** resulting kernel, especially for very low sigma values.
1124 ** As such while wierd it is prefered.
1126 ** I am told this method originally came from Photoshop.
1128 ** A properly normalized curve is generated (apart from edge clipping)
1129 ** even though we later normalize the result (for edge clipping)
1130 ** to allow the correct generation of a "Difference of Blurs".
1134 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1135 (void) ResetMagickMemory(kernel->values,0, (size_t)
1136 kernel->width*kernel->height*sizeof(double));
1137 /* Calculate a Positive 1D Gaussian */
1138 if ( sigma > MagickEpsilon )
1139 { sigma *= KernelRank; /* simplify loop expressions */
1140 alpha = 1.0/(2.0*sigma*sigma);
1141 beta= 1.0/(MagickSQ2PI*sigma );
1142 for ( u=-v; u <= v; u++) {
1143 kernel->values[(u+v)/KernelRank] +=
1144 exp(-((double)(u*u))*alpha)*beta;
1147 else /* special case - generate a unity kernel */
1148 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1150 /* Direct calculation without curve averaging */
1152 /* Calculate a Positive Gaussian */
1153 if ( sigma > MagickEpsilon )
1154 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1155 beta = 1.0/(MagickSQ2PI*sigma);
1156 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1157 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1159 else /* special case - generate a unity kernel */
1160 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1161 kernel->width*kernel->height*sizeof(double));
1162 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1165 /* Note the above kernel may have been 'clipped' by a user defined
1166 ** radius, producing a smaller (darker) kernel. Also for very small
1167 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1168 ** producing a very bright kernel.
1170 ** Normalization will still be needed.
1173 /* Normalize the 1D Gaussian Kernel
1175 ** NB: a CorrelateNormalize performs a normal Normalize if
1176 ** there are no negative values.
1178 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1179 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1181 /* rotate the 1D kernel by given angle */
1182 RotateKernelInfo(kernel, args->xi );
1187 sigma = fabs(args->sigma),
1190 if ( args->rho < 1.0 )
1191 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1193 kernel->width = (size_t)args->rho;
1194 kernel->x = kernel->y = 0;
1196 kernel->negative_range = kernel->positive_range = 0.0;
1197 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1198 kernel->height*sizeof(double));
1199 if (kernel->values == (double *) NULL)
1200 return(DestroyKernelInfo(kernel));
1202 /* A comet blur is half a 1D gaussian curve, so that the object is
1203 ** blurred in one direction only. This may not be quite the right
1204 ** curve to use so may change in the future. The function must be
1205 ** normalised after generation, which also resolves any clipping.
1207 ** As we are normalizing and not subtracting gaussians,
1208 ** there is no need for a divisor in the gaussian formula
1210 ** It is less comples
1212 if ( sigma > MagickEpsilon )
1215 #define KernelRank 3
1216 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1217 (void) ResetMagickMemory(kernel->values,0, (size_t)
1218 kernel->width*sizeof(double));
1219 sigma *= KernelRank; /* simplify the loop expression */
1220 A = 1.0/(2.0*sigma*sigma);
1221 /* B = 1.0/(MagickSQ2PI*sigma); */
1222 for ( u=0; u < v; u++) {
1223 kernel->values[u/KernelRank] +=
1224 exp(-((double)(u*u))*A);
1225 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1227 for (i=0; i < (ssize_t) kernel->width; i++)
1228 kernel->positive_range += kernel->values[i];
1230 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1231 /* B = 1.0/(MagickSQ2PI*sigma); */
1232 for ( i=0; i < (ssize_t) kernel->width; i++)
1233 kernel->positive_range +=
1235 exp(-((double)(i*i))*A);
1236 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1239 else /* special case - generate a unity kernel */
1240 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1241 kernel->width*kernel->height*sizeof(double));
1242 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1243 kernel->positive_range = 1.0;
1246 kernel->minimum = 0.0;
1247 kernel->maximum = kernel->values[0];
1248 kernel->negative_range = 0.0;
1250 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1251 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1255 /* Convolution Kernels - Well Known Constants */
1256 case LaplacianKernel:
1257 { switch ( (int) args->rho ) {
1259 default: /* laplacian square filter -- default */
1260 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1262 case 1: /* laplacian diamond filter */
1263 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1266 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1269 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1271 case 5: /* a 5x5 laplacian */
1272 kernel=ParseKernelArray(
1273 "5: -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");
1275 case 7: /* a 7x7 laplacian */
1276 kernel=ParseKernelArray(
1277 "7:-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" );
1279 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1280 kernel=ParseKernelArray(
1281 "5: 0,0,-1,0,0 0,-1,-2,-1,0 -1,-2,16,-2,-1 0,-1,-2,-1,0 0,0,-1,0,0");
1283 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1284 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1285 kernel=ParseKernelArray(
1286 "9: 0,-1,-1,-2,-2,-2,-1,-1,0 -1,-2,-4,-5,-5,-5,-4,-2,-1 -1,-4,-5,-3,-0,-3,-5,-4,-1 -2,-5,-3,12,24,12,-3,-5,-2 -2,-5,-0,24,40,24,-0,-5,-2 -2,-5,-3,12,24,12,-3,-5,-2 -1,-4,-5,-3,-0,-3,-5,-4,-1 -1,-2,-4,-5,-5,-5,-4,-2,-1 0,-1,-1,-2,-2,-2,-1,-1,0");
1289 if (kernel == (KernelInfo *) NULL)
1291 kernel->type = type;
1296 { /* Sobel with optional 'sub-types' */
1297 switch ( (int) args->rho ) {
1300 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1301 if (kernel == (KernelInfo *) NULL)
1303 kernel->type = type;
1306 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1307 if (kernel == (KernelInfo *) NULL)
1309 kernel->type = type;
1310 ScaleKernelInfo(kernel, 0.25, NoValue);
1313 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1314 if (kernel == (KernelInfo *) NULL)
1316 kernel->type = type;
1317 ScaleKernelInfo(kernel, 0.25, NoValue);
1320 if ( fabs(args->sigma) > MagickEpsilon )
1321 /* Rotate by correctly supplied 'angle' */
1322 RotateKernelInfo(kernel, args->sigma);
1323 else if ( args->rho > 30.0 || args->rho < -30.0 )
1324 /* Rotate by out of bounds 'type' */
1325 RotateKernelInfo(kernel, args->rho);
1329 { /* Simple Sobel Kernel */
1330 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1331 if (kernel == (KernelInfo *) NULL)
1333 kernel->type = type;
1334 RotateKernelInfo(kernel, args->rho);
1340 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1341 if (kernel == (KernelInfo *) NULL)
1343 kernel->type = type;
1344 RotateKernelInfo(kernel, args->rho);
1349 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1350 if (kernel == (KernelInfo *) NULL)
1352 kernel->type = type;
1353 RotateKernelInfo(kernel, args->rho);
1358 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1359 if (kernel == (KernelInfo *) NULL)
1361 kernel->type = type;
1362 RotateKernelInfo(kernel, args->rho);
1367 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1368 if (kernel == (KernelInfo *) NULL)
1370 kernel->type = type;
1371 RotateKernelInfo(kernel, args->rho);
1374 case FreiChenKernel:
1375 /* Direction is set to be left to right positive */
1376 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1377 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1378 { switch ( (int) args->rho ) {
1381 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1382 if (kernel == (KernelInfo *) NULL)
1384 kernel->type = type;
1385 kernel->values[3] = +MagickSQ2;
1386 kernel->values[5] = -MagickSQ2;
1387 CalcKernelMetaData(kernel); /* recalculate meta-data */
1390 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1391 if (kernel == (KernelInfo *) NULL)
1393 kernel->type = type;
1394 kernel->values[1] = kernel->values[3] = +MagickSQ2;
1395 kernel->values[5] = kernel->values[7] = -MagickSQ2;
1396 CalcKernelMetaData(kernel); /* recalculate meta-data */
1397 ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1400 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1401 if (kernel == (KernelInfo *) NULL)
1406 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1407 if (kernel == (KernelInfo *) NULL)
1409 kernel->type = type;
1410 kernel->values[3] = +MagickSQ2;
1411 kernel->values[5] = -MagickSQ2;
1412 CalcKernelMetaData(kernel); /* recalculate meta-data */
1413 ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1416 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1417 if (kernel == (KernelInfo *) NULL)
1419 kernel->type = type;
1420 kernel->values[1] = +MagickSQ2;
1421 kernel->values[7] = +MagickSQ2;
1422 CalcKernelMetaData(kernel);
1423 ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1426 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1427 if (kernel == (KernelInfo *) NULL)
1429 kernel->type = type;
1430 kernel->values[0] = +MagickSQ2;
1431 kernel->values[8] = -MagickSQ2;
1432 CalcKernelMetaData(kernel);
1433 ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1436 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1437 if (kernel == (KernelInfo *) NULL)
1439 kernel->type = type;
1440 kernel->values[2] = -MagickSQ2;
1441 kernel->values[6] = +MagickSQ2;
1442 CalcKernelMetaData(kernel);
1443 ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
1446 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1447 if (kernel == (KernelInfo *) NULL)
1449 kernel->type = type;
1450 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1453 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1454 if (kernel == (KernelInfo *) NULL)
1456 kernel->type = type;
1457 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1460 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1461 if (kernel == (KernelInfo *) NULL)
1463 kernel->type = type;
1464 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1467 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1468 if (kernel == (KernelInfo *) NULL)
1470 kernel->type = type;
1471 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1474 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1475 if (kernel == (KernelInfo *) NULL)
1477 kernel->type = type;
1478 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1481 if ( fabs(args->sigma) > MagickEpsilon )
1482 /* Rotate by correctly supplied 'angle' */
1483 RotateKernelInfo(kernel, args->sigma);
1484 else if ( args->rho > 30.0 || args->rho < -30.0 )
1485 /* Rotate by out of bounds 'type' */
1486 RotateKernelInfo(kernel, args->rho);
1490 /* Boolean Kernels */
1493 if (args->rho < 1.0)
1494 kernel->width = kernel->height = 3; /* default radius = 1 */
1496 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1497 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1499 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1500 kernel->height*sizeof(double));
1501 if (kernel->values == (double *) NULL)
1502 return(DestroyKernelInfo(kernel));
1504 /* set all kernel values within diamond area to scale given */
1505 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1506 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1507 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1508 kernel->positive_range += kernel->values[i] = args->sigma;
1510 kernel->values[i] = nan;
1511 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1515 case RectangleKernel:
1518 if ( type == SquareKernel )
1520 if (args->rho < 1.0)
1521 kernel->width = kernel->height = 3; /* default radius = 1 */
1523 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1524 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1525 scale = args->sigma;
1528 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1529 if ( args->rho < 1.0 || args->sigma < 1.0 )
1530 return(DestroyKernelInfo(kernel)); /* invalid args given */
1531 kernel->width = (size_t)args->rho;
1532 kernel->height = (size_t)args->sigma;
1533 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1534 args->psi < 0.0 || args->psi > (double)kernel->height )
1535 return(DestroyKernelInfo(kernel)); /* invalid args given */
1536 kernel->x = (ssize_t) args->xi;
1537 kernel->y = (ssize_t) args->psi;
1540 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1541 kernel->height*sizeof(double));
1542 if (kernel->values == (double *) NULL)
1543 return(DestroyKernelInfo(kernel));
1545 /* set all kernel values to scale given */
1546 u=(ssize_t) (kernel->width*kernel->height);
1547 for ( i=0; i < u; i++)
1548 kernel->values[i] = scale;
1549 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1550 kernel->positive_range = scale*u;
1556 limit = (ssize_t)(args->rho*args->rho);
1558 if (args->rho < 0.4) /* default radius approx 3.5 */
1559 kernel->width = kernel->height = 7L, limit = 10L;
1561 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1562 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1564 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1565 kernel->height*sizeof(double));
1566 if (kernel->values == (double *) NULL)
1567 return(DestroyKernelInfo(kernel));
1569 /* set all kernel values within disk area to scale given */
1570 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1571 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1572 if ((u*u+v*v) <= limit)
1573 kernel->positive_range += kernel->values[i] = args->sigma;
1575 kernel->values[i] = nan;
1576 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1581 if (args->rho < 1.0)
1582 kernel->width = kernel->height = 5; /* default radius 2 */
1584 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1585 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1587 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1588 kernel->height*sizeof(double));
1589 if (kernel->values == (double *) NULL)
1590 return(DestroyKernelInfo(kernel));
1592 /* set all kernel values along axises to given scale */
1593 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1594 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1595 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1596 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1597 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1602 if (args->rho < 1.0)
1603 kernel->width = kernel->height = 5; /* default radius 2 */
1605 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1606 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1608 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1609 kernel->height*sizeof(double));
1610 if (kernel->values == (double *) NULL)
1611 return(DestroyKernelInfo(kernel));
1613 /* set all kernel values along axises to given scale */
1614 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1615 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1616 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1617 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1618 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1621 /* HitAndMiss Kernels */
1630 if (args->rho < args->sigma)
1632 kernel->width = ((size_t)args->sigma)*2+1;
1633 limit1 = (ssize_t)(args->rho*args->rho);
1634 limit2 = (ssize_t)(args->sigma*args->sigma);
1638 kernel->width = ((size_t)args->rho)*2+1;
1639 limit1 = (ssize_t)(args->sigma*args->sigma);
1640 limit2 = (ssize_t)(args->rho*args->rho);
1643 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1645 kernel->height = kernel->width;
1646 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1647 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1648 kernel->height*sizeof(double));
1649 if (kernel->values == (double *) NULL)
1650 return(DestroyKernelInfo(kernel));
1652 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1653 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1654 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1655 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1656 { ssize_t radius=u*u+v*v;
1657 if (limit1 < radius && radius <= limit2)
1658 kernel->positive_range += kernel->values[i] = (double) scale;
1660 kernel->values[i] = nan;
1662 kernel->minimum = kernel->maximum = (double) scale;
1663 if ( type == PeaksKernel ) {
1664 /* set the central point in the middle */
1665 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1666 kernel->positive_range = 1.0;
1667 kernel->maximum = 1.0;
1673 kernel=ParseKernelArray("3: 0,0,0 -,1,- 1,1,1");
1674 if (kernel == (KernelInfo *) NULL)
1676 kernel->type = type;
1677 ExpandMirrorKernelInfo(kernel); /* mirror expansion of other kernels */
1682 kernel=ParseKernelArray("3: 0,0,- 0,1,1 -,1,-");
1683 if (kernel == (KernelInfo *) NULL)
1685 kernel->type = type;
1686 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1689 case ThinDiagonalsKernel:
1691 switch ( (int) args->rho ) {
1696 kernel=ParseKernelArray("3: 0,0,0 0,1,1 1,1,-");
1697 if (kernel == (KernelInfo *) NULL)
1699 kernel->type = type;
1700 new_kernel=ParseKernelArray("3: 0,0,1 0,1,1 0,1,-");
1701 if (new_kernel == (KernelInfo *) NULL)
1702 return(DestroyKernelInfo(kernel));
1703 new_kernel->type = type;
1704 LastKernelInfo(kernel)->next = new_kernel;
1705 ExpandMirrorKernelInfo(kernel);
1709 kernel=ParseKernelArray("3: 0,0,0 0,1,1 1,1,-");
1710 if (kernel == (KernelInfo *) NULL)
1712 kernel->type = type;
1713 RotateKernelInfo(kernel, args->sigma);
1716 kernel=ParseKernelArray("3: 0,0,1 0,1,1 0,1,-");
1717 if (kernel == (KernelInfo *) NULL)
1719 kernel->type = type;
1720 RotateKernelInfo(kernel, args->sigma);
1725 case LineEndsKernel:
1726 { /* Kernels for finding the end of thin lines */
1727 switch ( (int) args->rho ) {
1730 /* set of kernels to find all end of lines */
1731 kernel=AcquireKernelInfo("LineEnds:1>;LineEnds:2>");
1732 if (kernel == (KernelInfo *) NULL)
1736 /* kernel for 4-connected line ends - no rotation */
1737 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1738 if (kernel == (KernelInfo *) NULL)
1740 kernel->type = type;
1741 RotateKernelInfo(kernel, args->sigma);
1744 /* kernel to add for 8-connected lines - no rotation */
1745 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1746 if (kernel == (KernelInfo *) NULL)
1748 kernel->type = type;
1749 RotateKernelInfo(kernel, args->sigma);
1752 /* kernel to add for orthogonal line ends - does not find corners */
1753 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1754 if (kernel == (KernelInfo *) NULL)
1756 kernel->type = type;
1757 RotateKernelInfo(kernel, args->sigma);
1760 /* traditional line end - fails on last T end */
1761 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1762 if (kernel == (KernelInfo *) NULL)
1764 kernel->type = type;
1765 RotateKernelInfo(kernel, args->sigma);
1770 case LineJunctionsKernel:
1771 { /* kernels for finding the junctions of multiple lines */
1772 switch ( (int) args->rho ) {
1775 /* set of kernels to find all line junctions */
1776 kernel=AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>");
1777 if (kernel == (KernelInfo *) NULL)
1782 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1783 if (kernel == (KernelInfo *) NULL)
1785 kernel->type = type;
1786 RotateKernelInfo(kernel, args->sigma);
1789 /* Diagonal T Junctions */
1790 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1791 if (kernel == (KernelInfo *) NULL)
1793 kernel->type = type;
1794 RotateKernelInfo(kernel, args->sigma);
1797 /* Orthogonal T Junctions */
1798 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1799 if (kernel == (KernelInfo *) NULL)
1801 kernel->type = type;
1802 RotateKernelInfo(kernel, args->sigma);
1805 /* Diagonal X Junctions */
1806 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1807 if (kernel == (KernelInfo *) NULL)
1809 kernel->type = type;
1810 RotateKernelInfo(kernel, args->sigma);
1813 /* Orthogonal X Junctions - minimal diamond kernel */
1814 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1815 if (kernel == (KernelInfo *) NULL)
1817 kernel->type = type;
1818 RotateKernelInfo(kernel, args->sigma);
1824 { /* Ridges - Ridge finding kernels */
1827 switch ( (int) args->rho ) {
1830 kernel=ParseKernelArray("3x1:0,1,0");
1831 if (kernel == (KernelInfo *) NULL)
1833 kernel->type = type;
1834 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1837 kernel=ParseKernelArray("4x1:0,1,1,0");
1838 if (kernel == (KernelInfo *) NULL)
1840 kernel->type = type;
1841 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1843 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1844 /* Unfortunatally we can not yet rotate a non-square kernel */
1845 /* But then we can't flip a non-symetrical kernel either */
1846 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1847 if (new_kernel == (KernelInfo *) NULL)
1848 return(DestroyKernelInfo(kernel));
1849 new_kernel->type = type;
1850 LastKernelInfo(kernel)->next = new_kernel;
1851 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1852 if (new_kernel == (KernelInfo *) NULL)
1853 return(DestroyKernelInfo(kernel));
1854 new_kernel->type = type;
1855 LastKernelInfo(kernel)->next = new_kernel;
1856 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1857 if (new_kernel == (KernelInfo *) NULL)
1858 return(DestroyKernelInfo(kernel));
1859 new_kernel->type = type;
1860 LastKernelInfo(kernel)->next = new_kernel;
1861 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1862 if (new_kernel == (KernelInfo *) NULL)
1863 return(DestroyKernelInfo(kernel));
1864 new_kernel->type = type;
1865 LastKernelInfo(kernel)->next = new_kernel;
1866 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1867 if (new_kernel == (KernelInfo *) NULL)
1868 return(DestroyKernelInfo(kernel));
1869 new_kernel->type = type;
1870 LastKernelInfo(kernel)->next = new_kernel;
1871 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1872 if (new_kernel == (KernelInfo *) NULL)
1873 return(DestroyKernelInfo(kernel));
1874 new_kernel->type = type;
1875 LastKernelInfo(kernel)->next = new_kernel;
1876 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1877 if (new_kernel == (KernelInfo *) NULL)
1878 return(DestroyKernelInfo(kernel));
1879 new_kernel->type = type;
1880 LastKernelInfo(kernel)->next = new_kernel;
1881 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1882 if (new_kernel == (KernelInfo *) NULL)
1883 return(DestroyKernelInfo(kernel));
1884 new_kernel->type = type;
1885 LastKernelInfo(kernel)->next = new_kernel;
1890 case ConvexHullKernel:
1894 /* first set of 8 kernels */
1895 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1896 if (kernel == (KernelInfo *) NULL)
1898 kernel->type = type;
1899 ExpandRotateKernelInfo(kernel, 90.0);
1900 /* append the mirror versions too - no flip function yet */
1901 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1902 if (new_kernel == (KernelInfo *) NULL)
1903 return(DestroyKernelInfo(kernel));
1904 new_kernel->type = type;
1905 ExpandRotateKernelInfo(new_kernel, 90.0);
1906 LastKernelInfo(kernel)->next = new_kernel;
1909 case SkeletonKernel:
1913 switch ( (int) args->rho ) {
1916 /* Traditional Skeleton...
1917 ** A cyclically rotated single kernel
1919 kernel=ParseKernelArray("3: 0,0,0 -,1,- 1,1,1");
1920 if (kernel == (KernelInfo *) NULL)
1922 kernel->type = type;
1923 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1926 /* HIPR Variation of the cyclic skeleton
1927 ** Corners of the traditional method made more forgiving,
1928 ** but the retain the same cyclic order.
1930 kernel=ParseKernelArray("3: 0,0,0 -,1,- 1,1,1");
1931 if (kernel == (KernelInfo *) NULL)
1933 kernel->type = type;
1934 new_kernel=ParseKernelArray("3: -,0,0 1,1,0 -,1,-");
1935 if (new_kernel == (KernelInfo *) NULL)
1937 new_kernel->type = type;
1938 LastKernelInfo(kernel)->next = new_kernel;
1939 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1944 /* Distance Measuring Kernels */
1945 case ChebyshevKernel:
1947 if (args->rho < 1.0)
1948 kernel->width = kernel->height = 3; /* default radius = 1 */
1950 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1951 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1953 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1954 kernel->height*sizeof(double));
1955 if (kernel->values == (double *) NULL)
1956 return(DestroyKernelInfo(kernel));
1958 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1959 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1960 kernel->positive_range += ( kernel->values[i] =
1961 args->sigma*((labs((long) u)>labs((long) v)) ? labs((long) u) : labs((long) v)) );
1962 kernel->maximum = kernel->values[0];
1965 case ManhattanKernel:
1967 if (args->rho < 1.0)
1968 kernel->width = kernel->height = 3; /* default radius = 1 */
1970 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1971 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1973 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1974 kernel->height*sizeof(double));
1975 if (kernel->values == (double *) NULL)
1976 return(DestroyKernelInfo(kernel));
1978 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1979 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1980 kernel->positive_range += ( kernel->values[i] =
1981 args->sigma*(labs((long) u)+labs((long) v)) );
1982 kernel->maximum = kernel->values[0];
1985 case EuclideanKernel:
1987 if (args->rho < 1.0)
1988 kernel->width = kernel->height = 3; /* default radius = 1 */
1990 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1991 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1993 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1994 kernel->height*sizeof(double));
1995 if (kernel->values == (double *) NULL)
1996 return(DestroyKernelInfo(kernel));
1998 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1999 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2000 kernel->positive_range += ( kernel->values[i] =
2001 args->sigma*sqrt((double)(u*u+v*v)) );
2002 kernel->maximum = kernel->values[0];
2008 /* Unity or No-Op Kernel - Basically just a single pixel on its own */
2009 kernel=ParseKernelArray("1:1");
2010 if (kernel == (KernelInfo *) NULL)
2012 kernel->type = ( type == UnityKernel ) ? UnityKernel : UndefinedKernel;
2022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2026 % C l o n e K e r n e l I n f o %
2030 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2032 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2033 % can be modified without effecting the original. The cloned kernel should
2034 % be destroyed using DestoryKernelInfo() when no longer needed.
2036 % The format of the CloneKernelInfo method is:
2038 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2040 % A description of each parameter follows:
2042 % o kernel: the Morphology/Convolution kernel to be cloned
2045 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2053 assert(kernel != (KernelInfo *) NULL);
2054 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2055 if (new_kernel == (KernelInfo *) NULL)
2057 *new_kernel=(*kernel); /* copy values in structure */
2059 /* replace the values with a copy of the values */
2060 new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2061 kernel->height*sizeof(double));
2062 if (new_kernel->values == (double *) NULL)
2063 return(DestroyKernelInfo(new_kernel));
2064 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2065 new_kernel->values[i]=kernel->values[i];
2067 /* Also clone the next kernel in the kernel list */
2068 if ( kernel->next != (KernelInfo *) NULL ) {
2069 new_kernel->next = CloneKernelInfo(kernel->next);
2070 if ( new_kernel->next == (KernelInfo *) NULL )
2071 return(DestroyKernelInfo(new_kernel));
2078 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2082 % D e s t r o y K e r n e l I n f o %
2086 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2088 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2091 % The format of the DestroyKernelInfo method is:
2093 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2095 % A description of each parameter follows:
2097 % o kernel: the Morphology/Convolution kernel to be destroyed
2100 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2102 assert(kernel != (KernelInfo *) NULL);
2104 if ( kernel->next != (KernelInfo *) NULL )
2105 kernel->next = DestroyKernelInfo(kernel->next);
2107 kernel->values = (double *)RelinquishMagickMemory(kernel->values);
2108 kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
2113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2117 + E x p a n d M i r r o r K e r n e l I n f o %
2121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2123 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2124 % sequence of 90-degree rotated kernels but providing a reflected 180
2125 % rotatation, before the -/+ 90-degree rotations.
2127 % This special rotation order produces a better, more symetrical thinning of
2130 % The format of the ExpandMirrorKernelInfo method is:
2132 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2134 % A description of each parameter follows:
2136 % o kernel: the Morphology/Convolution kernel
2138 % This function is only internel to this module, as it is not finalized,
2139 % especially with regard to non-orthogonal angles, and rotation of larger
2144 static void FlopKernelInfo(KernelInfo *kernel)
2145 { /* Do a Flop by reversing each row. */
2153 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2154 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2155 t=k[x], k[x]=k[r], k[r]=t;
2157 kernel->x = kernel->width - kernel->x - 1;
2158 angle = fmod(angle+180.0, 360.0);
2162 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2170 clone = CloneKernelInfo(last);
2171 RotateKernelInfo(clone, 180); /* flip */
2172 LastKernelInfo(last)->next = clone;
2175 clone = CloneKernelInfo(last);
2176 RotateKernelInfo(clone, 90); /* transpose */
2177 LastKernelInfo(last)->next = clone;
2180 clone = CloneKernelInfo(last);
2181 RotateKernelInfo(clone, 180); /* flop */
2182 LastKernelInfo(last)->next = clone;
2188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2192 + E x p a n d R o t a t e K e r n e l I n f o %
2196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2198 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2199 % incrementally by the angle given, until the first kernel repeats.
2201 % WARNING: 45 degree rotations only works for 3x3 kernels.
2202 % While 90 degree roatations only works for linear and square kernels
2204 % The format of the ExpandRotateKernelInfo method is:
2206 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2208 % A description of each parameter follows:
2210 % o kernel: the Morphology/Convolution kernel
2212 % o angle: angle to rotate in degrees
2214 % This function is only internel to this module, as it is not finalized,
2215 % especially with regard to non-orthogonal angles, and rotation of larger
2219 /* Internal Routine - Return true if two kernels are the same */
2220 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2221 const KernelInfo *kernel2)
2226 /* check size and origin location */
2227 if ( kernel1->width != kernel2->width
2228 || kernel1->height != kernel2->height
2229 || kernel1->x != kernel2->x
2230 || kernel1->y != kernel2->y )
2233 /* check actual kernel values */
2234 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2235 /* Test for Nan equivelence */
2236 if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2238 if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2240 /* Test actual values are equivelent */
2241 if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
2248 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2256 clone = CloneKernelInfo(last);
2257 RotateKernelInfo(clone, angle);
2258 if ( SameKernelInfo(kernel, clone) == MagickTrue )
2260 LastKernelInfo(last)->next = clone;
2263 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2268 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2272 + C a l c M e t a K e r n a l I n f o %
2276 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2278 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2279 % using the kernel values. This should only ne used if it is not posible to
2280 % calculate that meta-data in some easier way.
2282 % It is important that the meta-data is correct before ScaleKernelInfo() is
2283 % used to perform kernel normalization.
2285 % The format of the CalcKernelMetaData method is:
2287 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2289 % A description of each parameter follows:
2291 % o kernel: the Morphology/Convolution kernel to modify
2293 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2294 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2295 % however is not true for flat-shaped morphological kernels.
2297 % WARNING: Only the specific kernel pointed to is modified, not a list of
2300 % This is an internal function and not expected to be useful outside this
2301 % module. This could change however.
2303 static void CalcKernelMetaData(KernelInfo *kernel)
2308 kernel->minimum = kernel->maximum = 0.0;
2309 kernel->negative_range = kernel->positive_range = 0.0;
2310 for (i=0; i < (kernel->width*kernel->height); i++)
2312 if ( fabs(kernel->values[i]) < MagickEpsilon )
2313 kernel->values[i] = 0.0;
2314 ( kernel->values[i] < 0)
2315 ? ( kernel->negative_range += kernel->values[i] )
2316 : ( kernel->positive_range += kernel->values[i] );
2317 Minimize(kernel->minimum, kernel->values[i]);
2318 Maximize(kernel->maximum, kernel->values[i]);
2325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2329 % M o r p h o l o g y A p p l y %
2333 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2335 % MorphologyApply() applies a morphological method, multiple times using
2336 % a list of multiple kernels.
2338 % It is basically equivelent to as MorphologyImageChannel() (see below) but
2339 % without any user controls. This allows internel programs to use this
2340 % function, to actually perform a specific task without posible interference
2341 % by any API user supplied settings.
2343 % It is MorphologyImageChannel() task to extract any such user controls, and
2344 % pass them to this function for processing.
2346 % More specifically kernels are not normalized/scaled/blended by the
2347 % 'convolve:scale' Image Artifact (setting), nor is the convolve bias
2348 % (-bias setting or image->bias) loooked at, but must be supplied from the
2349 % function arguments.
2351 % The format of the MorphologyApply method is:
2353 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2354 % const ssize_t iterations,const KernelInfo *kernel,
2355 % const CompositeMethod compose, const double bias,
2356 % ExceptionInfo *exception)
2358 % A description of each parameter follows:
2360 % o image: the source image
2362 % o method: the morphology method to be applied.
2364 % o iterations: apply the operation this many times (or no change).
2365 % A value of -1 means loop until no change found.
2366 % How this is applied may depend on the morphology method.
2367 % Typically this is a value of 1.
2369 % o channel: the channel type.
2371 % o kernel: An array of double representing the morphology kernel.
2373 % o compose: How to handle or merge multi-kernel results.
2374 % If 'UndefinedCompositeOp' use default for the Morphology method.
2375 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2376 % Otherwise merge the results using the compose method given.
2378 % o bias: Convolution Output Bias.
2380 % o exception: return any errors or warnings in this structure.
2385 /* Apply a Morphology Primative to an image using the given kernel.
2386 ** Two pre-created images must be provided, no image is created.
2387 ** It returns the number of pixels that changed betwene the images
2388 ** for convergence determination.
2390 static size_t MorphologyPrimitive(const Image *image, Image
2391 *result_image, const MorphologyMethod method, const ChannelType channel,
2392 const KernelInfo *kernel,const double bias,ExceptionInfo *exception)
2394 #define MorphologyTag "Morphology/Image"
2410 assert(image != (Image *) NULL);
2411 assert(image->signature == MagickSignature);
2412 assert(result_image != (Image *) NULL);
2413 assert(result_image->signature == MagickSignature);
2414 assert(kernel != (KernelInfo *) NULL);
2415 assert(kernel->signature == MagickSignature);
2416 assert(exception != (ExceptionInfo *) NULL);
2417 assert(exception->signature == MagickSignature);
2423 p_view=AcquireCacheView(image);
2424 q_view=AcquireCacheView(result_image);
2426 /* Some methods (including convolve) needs use a reflected kernel.
2427 * Adjust 'origin' offsets to loop though kernel as a reflection.
2432 case ConvolveMorphology:
2433 case DilateMorphology:
2434 case DilateIntensityMorphology:
2435 case DistanceMorphology:
2436 /* kernel needs to used with reflection about origin */
2437 offx = (ssize_t) kernel->width-offx-1;
2438 offy = (ssize_t) kernel->height-offy-1;
2440 case ErodeMorphology:
2441 case ErodeIntensityMorphology:
2442 case HitAndMissMorphology:
2443 case ThinningMorphology:
2444 case ThickenMorphology:
2445 /* kernel is used as is, without reflection */
2448 assert("Not a Primitive Morphology Method" != (char *) NULL);
2453 if ( method == ConvolveMorphology && kernel->width == 1 )
2454 { /* Special handling (for speed) of vertical (blur) kernels.
2455 ** This performs its handling in columns rather than in rows.
2456 ** This is only done fo convolve as it is the only method that
2457 ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
2459 ** Timing tests (on single CPU laptop)
2460 ** Using a vertical 1-d Blue with normal row-by-row (below)
2461 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2463 ** Using this column method
2464 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2467 ** Anthony Thyssen, 14 June 2010
2472 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2473 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2475 for (x=0; x < (ssize_t) image->columns; x++)
2477 register const PixelPacket
2480 register const IndexPacket
2481 *restrict p_indexes;
2483 register PixelPacket
2486 register IndexPacket
2487 *restrict q_indexes;
2495 if (status == MagickFalse)
2497 p=GetCacheViewVirtualPixels(p_view, x, -offy,1,
2498 image->rows+kernel->height, exception);
2499 q=GetCacheViewAuthenticPixels(q_view,x,0,1,result_image->rows,exception);
2500 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2505 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
2506 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
2507 r = offy; /* offset to the origin pixel in 'p' */
2509 for (y=0; y < (ssize_t) image->rows; y++)
2514 register const double
2517 register const PixelPacket
2520 register const IndexPacket
2521 *restrict k_indexes;
2526 /* Copy input image to the output image for unused channels
2527 * This removes need for 'cloning' a new image every iteration
2530 if (image->colorspace == CMYKColorspace)
2531 q_indexes[y] = p_indexes[r];
2533 /* Set the bias of the weighted average output */
2538 result.index = bias;
2541 /* Weighted Average of pixels using reflected kernel
2543 ** NOTE for correct working of this operation for asymetrical
2544 ** kernels, the kernel needs to be applied in its reflected form.
2545 ** That is its values needs to be reversed.
2547 k = &kernel->values[ kernel->height-1 ];
2549 k_indexes = p_indexes;
2550 if ( ((channel & SyncChannels) == 0 ) ||
2551 (image->matte == MagickFalse) )
2552 { /* No 'Sync' involved.
2553 ** Convolution is simple greyscale channel operation
2555 for (v=0; v < (ssize_t) kernel->height; v++) {
2556 if ( IsNan(*k) ) continue;
2557 result.red += (*k)*k_pixels->red;
2558 result.green += (*k)*k_pixels->green;
2559 result.blue += (*k)*k_pixels->blue;
2560 result.opacity += (*k)*k_pixels->opacity;
2561 if ( image->colorspace == CMYKColorspace)
2562 result.index += (*k)*(*k_indexes);
2567 if ((channel & RedChannel) != 0)
2568 q->red = ClampToQuantum(result.red);
2569 if ((channel & GreenChannel) != 0)
2570 q->green = ClampToQuantum(result.green);
2571 if ((channel & BlueChannel) != 0)
2572 q->blue = ClampToQuantum(result.blue);
2573 if ((channel & OpacityChannel) != 0
2574 && image->matte == MagickTrue )
2575 q->opacity = ClampToQuantum(result.opacity);
2576 if ((channel & IndexChannel) != 0
2577 && image->colorspace == CMYKColorspace)
2578 q_indexes[x] = ClampToQuantum(result.index);
2581 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2582 ** Weight the color channels with Alpha Channel so that
2583 ** transparent pixels are not part of the results.
2586 alpha, /* alpha weighting of colors : kernel*alpha */
2587 gamma; /* divisor, sum of color weighting values */
2590 for (v=0; v < (ssize_t) kernel->height; v++) {
2591 if ( IsNan(*k) ) continue;
2592 alpha=(*k)*(QuantumScale*(QuantumRange-k_pixels->opacity));
2594 result.red += alpha*k_pixels->red;
2595 result.green += alpha*k_pixels->green;
2596 result.blue += alpha*k_pixels->blue;
2597 result.opacity += (*k)*k_pixels->opacity;
2598 if ( image->colorspace == CMYKColorspace)
2599 result.index += alpha*(*k_indexes);
2604 /* Sync'ed channels, all channels are modified */
2605 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2606 q->red = ClampToQuantum(gamma*result.red);
2607 q->green = ClampToQuantum(gamma*result.green);
2608 q->blue = ClampToQuantum(gamma*result.blue);
2609 q->opacity = ClampToQuantum(result.opacity);
2610 if (image->colorspace == CMYKColorspace)
2611 q_indexes[x] = ClampToQuantum(gamma*result.index);
2614 /* Count up changed pixels */
2615 if ( ( p[r].red != q->red )
2616 || ( p[r].green != q->green )
2617 || ( p[r].blue != q->blue )
2618 || ( p[r].opacity != q->opacity )
2619 || ( image->colorspace == CMYKColorspace &&
2620 p_indexes[r] != q_indexes[x] ) )
2621 changed++; /* The pixel was changed in some way! */
2625 if ( SyncCacheViewAuthenticPixels(q_view,exception) == MagickFalse)
2627 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2632 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2633 #pragma omp critical (MagickCore_MorphologyImage)
2635 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2636 if (proceed == MagickFalse)
2640 result_image->type=image->type;
2641 q_view=DestroyCacheView(q_view);
2642 p_view=DestroyCacheView(p_view);
2643 return(status ? (size_t) changed : 0);
2647 ** Normal handling of horizontal or rectangular kernels (row by row)
2649 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2650 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2652 for (y=0; y < (ssize_t) image->rows; y++)
2654 register const PixelPacket
2657 register const IndexPacket
2658 *restrict p_indexes;
2660 register PixelPacket
2663 register IndexPacket
2664 *restrict q_indexes;
2672 if (status == MagickFalse)
2674 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
2675 image->columns+kernel->width, kernel->height, exception);
2676 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
2678 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2683 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
2684 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
2685 r = (image->columns+kernel->width)*offy+offx; /* offset to origin in 'p' */
2687 for (x=0; x < (ssize_t) image->columns; x++)
2695 register const double
2698 register const PixelPacket
2701 register const IndexPacket
2702 *restrict k_indexes;
2709 /* Copy input image to the output image for unused channels
2710 * This removes need for 'cloning' a new image every iteration
2713 if (image->colorspace == CMYKColorspace)
2714 q_indexes[x] = p_indexes[r];
2721 min.index = (MagickRealType) QuantumRange;
2726 max.index = (MagickRealType) 0;
2727 /* default result is the original pixel value */
2728 result.red = (MagickRealType) p[r].red;
2729 result.green = (MagickRealType) p[r].green;
2730 result.blue = (MagickRealType) p[r].blue;
2731 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
2733 if ( image->colorspace == CMYKColorspace)
2734 result.index = (MagickRealType) p_indexes[r];
2737 case ConvolveMorphology:
2738 /* Set the bias of the weighted average output */
2743 result.index = bias;
2745 case DilateIntensityMorphology:
2746 case ErodeIntensityMorphology:
2747 /* use a boolean flag indicating when first match found */
2748 result.red = 0.0; /* result is not used otherwise */
2755 case ConvolveMorphology:
2756 /* Weighted Average of pixels using reflected kernel
2758 ** NOTE for correct working of this operation for asymetrical
2759 ** kernels, the kernel needs to be applied in its reflected form.
2760 ** That is its values needs to be reversed.
2762 ** Correlation is actually the same as this but without reflecting
2763 ** the kernel, and thus 'lower-level' that Convolution. However
2764 ** as Convolution is the more common method used, and it does not
2765 ** really cost us much in terms of processing to use a reflected
2766 ** kernel, so it is Convolution that is implemented.
2768 ** Correlation will have its kernel reflected before calling
2769 ** this function to do a Convolve.
2771 ** For more details of Correlation vs Convolution see
2772 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2774 k = &kernel->values[ kernel->width*kernel->height-1 ];
2776 k_indexes = p_indexes;
2777 if ( ((channel & SyncChannels) == 0 ) ||
2778 (image->matte == MagickFalse) )
2779 { /* No 'Sync' involved.
2780 ** Convolution is simple greyscale channel operation
2782 for (v=0; v < (ssize_t) kernel->height; v++) {
2783 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2784 if ( IsNan(*k) ) continue;
2785 result.red += (*k)*k_pixels[u].red;
2786 result.green += (*k)*k_pixels[u].green;
2787 result.blue += (*k)*k_pixels[u].blue;
2788 result.opacity += (*k)*k_pixels[u].opacity;
2789 if ( image->colorspace == CMYKColorspace)
2790 result.index += (*k)*k_indexes[u];
2792 k_pixels += image->columns+kernel->width;
2793 k_indexes += image->columns+kernel->width;
2795 if ((channel & RedChannel) != 0)
2796 q->red = ClampToQuantum(result.red);
2797 if ((channel & GreenChannel) != 0)
2798 q->green = ClampToQuantum(result.green);
2799 if ((channel & BlueChannel) != 0)
2800 q->blue = ClampToQuantum(result.blue);
2801 if ((channel & OpacityChannel) != 0
2802 && image->matte == MagickTrue )
2803 q->opacity = ClampToQuantum(result.opacity);
2804 if ((channel & IndexChannel) != 0
2805 && image->colorspace == CMYKColorspace)
2806 q_indexes[x] = ClampToQuantum(result.index);
2809 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2810 ** Weight the color channels with Alpha Channel so that
2811 ** transparent pixels are not part of the results.
2814 alpha, /* alpha weighting of colors : kernel*alpha */
2815 gamma; /* divisor, sum of color weighting values */
2818 for (v=0; v < (ssize_t) kernel->height; v++) {
2819 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2820 if ( IsNan(*k) ) continue;
2821 alpha=(*k)*(QuantumScale*(QuantumRange-
2822 k_pixels[u].opacity));
2824 result.red += alpha*k_pixels[u].red;
2825 result.green += alpha*k_pixels[u].green;
2826 result.blue += alpha*k_pixels[u].blue;
2827 result.opacity += (*k)*k_pixels[u].opacity;
2828 if ( image->colorspace == CMYKColorspace)
2829 result.index += alpha*k_indexes[u];
2831 k_pixels += image->columns+kernel->width;
2832 k_indexes += image->columns+kernel->width;
2834 /* Sync'ed channels, all channels are modified */
2835 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2836 q->red = ClampToQuantum(gamma*result.red);
2837 q->green = ClampToQuantum(gamma*result.green);
2838 q->blue = ClampToQuantum(gamma*result.blue);
2839 q->opacity = ClampToQuantum(result.opacity);
2840 if (image->colorspace == CMYKColorspace)
2841 q_indexes[x] = ClampToQuantum(gamma*result.index);
2845 case ErodeMorphology:
2846 /* Minimum Value within kernel neighbourhood
2848 ** NOTE that the kernel is not reflected for this operation!
2850 ** NOTE: in normal Greyscale Morphology, the kernel value should
2851 ** be added to the real value, this is currently not done, due to
2852 ** the nature of the boolean kernels being used.
2856 k_indexes = p_indexes;
2857 for (v=0; v < (ssize_t) kernel->height; v++) {
2858 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2859 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2860 Minimize(min.red, (double) k_pixels[u].red);
2861 Minimize(min.green, (double) k_pixels[u].green);
2862 Minimize(min.blue, (double) k_pixels[u].blue);
2863 Minimize(min.opacity,
2864 QuantumRange-(double) k_pixels[u].opacity);
2865 if ( image->colorspace == CMYKColorspace)
2866 Minimize(min.index, (double) k_indexes[u]);
2868 k_pixels += image->columns+kernel->width;
2869 k_indexes += image->columns+kernel->width;
2873 case DilateMorphology:
2874 /* Maximum Value within kernel neighbourhood
2876 ** NOTE for correct working of this operation for asymetrical
2877 ** kernels, the kernel needs to be applied in its reflected form.
2878 ** That is its values needs to be reversed.
2880 ** NOTE: in normal Greyscale Morphology, the kernel value should
2881 ** be added to the real value, this is currently not done, due to
2882 ** the nature of the boolean kernels being used.
2885 k = &kernel->values[ kernel->width*kernel->height-1 ];
2887 k_indexes = p_indexes;
2888 for (v=0; v < (ssize_t) kernel->height; v++) {
2889 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2890 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2891 Maximize(max.red, (double) k_pixels[u].red);
2892 Maximize(max.green, (double) k_pixels[u].green);
2893 Maximize(max.blue, (double) k_pixels[u].blue);
2894 Maximize(max.opacity,
2895 QuantumRange-(double) k_pixels[u].opacity);
2896 if ( image->colorspace == CMYKColorspace)
2897 Maximize(max.index, (double) k_indexes[u]);
2899 k_pixels += image->columns+kernel->width;
2900 k_indexes += image->columns+kernel->width;
2904 case HitAndMissMorphology:
2905 case ThinningMorphology:
2906 case ThickenMorphology:
2907 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
2909 ** NOTE that the kernel is not reflected for this operation,
2910 ** and consists of both foreground and background pixel
2911 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
2912 ** with either Nan or 0.5 values for don't care.
2914 ** Note that this will never produce a meaningless negative
2915 ** result. Such results can cause Thinning/Thicken to not work
2916 ** correctly when used against a greyscale image.
2920 k_indexes = p_indexes;
2921 for (v=0; v < (ssize_t) kernel->height; v++) {
2922 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2923 if ( IsNan(*k) ) continue;
2925 { /* minimim of foreground pixels */
2926 Minimize(min.red, (double) k_pixels[u].red);
2927 Minimize(min.green, (double) k_pixels[u].green);
2928 Minimize(min.blue, (double) k_pixels[u].blue);
2929 Minimize(min.opacity,
2930 QuantumRange-(double) k_pixels[u].opacity);
2931 if ( image->colorspace == CMYKColorspace)
2932 Minimize(min.index, (double) k_indexes[u]);
2934 else if ( (*k) < 0.3 )
2935 { /* maximum of background pixels */
2936 Maximize(max.red, (double) k_pixels[u].red);
2937 Maximize(max.green, (double) k_pixels[u].green);
2938 Maximize(max.blue, (double) k_pixels[u].blue);
2939 Maximize(max.opacity,
2940 QuantumRange-(double) k_pixels[u].opacity);
2941 if ( image->colorspace == CMYKColorspace)
2942 Maximize(max.index, (double) k_indexes[u]);
2945 k_pixels += image->columns+kernel->width;
2946 k_indexes += image->columns+kernel->width;
2948 /* Pattern Match if difference is positive */
2949 min.red -= max.red; Maximize( min.red, 0.0 );
2950 min.green -= max.green; Maximize( min.green, 0.0 );
2951 min.blue -= max.blue; Maximize( min.blue, 0.0 );
2952 min.opacity -= max.opacity; Maximize( min.opacity, 0.0 );
2953 min.index -= max.index; Maximize( min.index, 0.0 );
2956 case ErodeIntensityMorphology:
2957 /* Select Pixel with Minimum Intensity within kernel neighbourhood
2959 ** WARNING: the intensity test fails for CMYK and does not
2960 ** take into account the moderating effect of the alpha channel
2961 ** on the intensity.
2963 ** NOTE that the kernel is not reflected for this operation!
2967 k_indexes = p_indexes;
2968 for (v=0; v < (ssize_t) kernel->height; v++) {
2969 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2970 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2971 if ( result.red == 0.0 ||
2972 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
2973 /* copy the whole pixel - no channel selection */
2975 if ( result.red > 0.0 ) changed++;
2979 k_pixels += image->columns+kernel->width;
2980 k_indexes += image->columns+kernel->width;
2984 case DilateIntensityMorphology:
2985 /* Select Pixel with Maximum Intensity within kernel neighbourhood
2987 ** WARNING: the intensity test fails for CMYK and does not
2988 ** take into account the moderating effect of the alpha channel
2989 ** on the intensity (yet).
2991 ** NOTE for correct working of this operation for asymetrical
2992 ** kernels, the kernel needs to be applied in its reflected form.
2993 ** That is its values needs to be reversed.
2995 k = &kernel->values[ kernel->width*kernel->height-1 ];
2997 k_indexes = p_indexes;
2998 for (v=0; v < (ssize_t) kernel->height; v++) {
2999 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3000 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3001 if ( result.red == 0.0 ||
3002 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
3003 /* copy the whole pixel - no channel selection */
3005 if ( result.red > 0.0 ) changed++;
3009 k_pixels += image->columns+kernel->width;
3010 k_indexes += image->columns+kernel->width;
3015 case DistanceMorphology:
3016 /* Add kernel Value and select the minimum value found.
3017 ** The result is a iterative distance from edge of image shape.
3019 ** All Distance Kernels are symetrical, but that may not always
3020 ** be the case. For example how about a distance from left edges?
3021 ** To work correctly with asymetrical kernels the reflected kernel
3022 ** needs to be applied.
3024 ** Actually this is really a GreyErode with a negative kernel!
3027 k = &kernel->values[ kernel->width*kernel->height-1 ];
3029 k_indexes = p_indexes;
3030 for (v=0; v < (ssize_t) kernel->height; v++) {
3031 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3032 if ( IsNan(*k) ) continue;
3033 Minimize(result.red, (*k)+k_pixels[u].red);
3034 Minimize(result.green, (*k)+k_pixels[u].green);
3035 Minimize(result.blue, (*k)+k_pixels[u].blue);
3036 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
3037 if ( image->colorspace == CMYKColorspace)
3038 Minimize(result.index, (*k)+k_indexes[u]);
3040 k_pixels += image->columns+kernel->width;
3041 k_indexes += image->columns+kernel->width;
3045 case UndefinedMorphology:
3047 break; /* Do nothing */
3049 /* Final mathematics of results (combine with original image?)
3051 ** NOTE: Difference Morphology operators Edge* and *Hat could also
3052 ** be done here but works better with iteration as a image difference
3053 ** in the controling function (below). Thicken and Thinning however
3054 ** should be done here so thay can be iterated correctly.
3057 case HitAndMissMorphology:
3058 case ErodeMorphology:
3059 result = min; /* minimum of neighbourhood */
3061 case DilateMorphology:
3062 result = max; /* maximum of neighbourhood */
3064 case ThinningMorphology:
3065 /* subtract pattern match from original */
3066 result.red -= min.red;
3067 result.green -= min.green;
3068 result.blue -= min.blue;
3069 result.opacity -= min.opacity;
3070 result.index -= min.index;
3072 case ThickenMorphology:
3073 /* Add the pattern matchs to the original */
3074 result.red += min.red;
3075 result.green += min.green;
3076 result.blue += min.blue;
3077 result.opacity += min.opacity;
3078 result.index += min.index;
3081 /* result directly calculated or assigned */
3084 /* Assign the resulting pixel values - Clamping Result */
3086 case UndefinedMorphology:
3087 case ConvolveMorphology:
3088 case DilateIntensityMorphology:
3089 case ErodeIntensityMorphology:
3090 break; /* full pixel was directly assigned - not a channel method */
3092 if ((channel & RedChannel) != 0)
3093 q->red = ClampToQuantum(result.red);
3094 if ((channel & GreenChannel) != 0)
3095 q->green = ClampToQuantum(result.green);
3096 if ((channel & BlueChannel) != 0)
3097 q->blue = ClampToQuantum(result.blue);
3098 if ((channel & OpacityChannel) != 0
3099 && image->matte == MagickTrue )
3100 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
3101 if ((channel & IndexChannel) != 0
3102 && image->colorspace == CMYKColorspace)
3103 q_indexes[x] = ClampToQuantum(result.index);
3106 /* Count up changed pixels */
3107 if ( ( p[r].red != q->red )
3108 || ( p[r].green != q->green )
3109 || ( p[r].blue != q->blue )
3110 || ( p[r].opacity != q->opacity )
3111 || ( image->colorspace == CMYKColorspace &&
3112 p_indexes[r] != q_indexes[x] ) )
3113 changed++; /* The pixel was changed in some way! */
3117 if ( SyncCacheViewAuthenticPixels(q_view,exception) == MagickFalse)
3119 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3124 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3125 #pragma omp critical (MagickCore_MorphologyImage)
3127 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3128 if (proceed == MagickFalse)
3132 result_image->type=image->type;
3133 q_view=DestroyCacheView(q_view);
3134 p_view=DestroyCacheView(p_view);
3135 return(status ? (size_t) changed : 0);
3139 MagickExport Image *MorphologyApply(const Image *image, const ChannelType
3140 channel,const MorphologyMethod method, const ssize_t iterations,
3141 const KernelInfo *kernel, const CompositeOperator compose,
3142 const double bias, ExceptionInfo *exception)
3145 *curr_image, /* Image we are working with or iterating */
3146 *work_image, /* secondary image for primative iteration */
3147 *save_image, /* saved image - for 'edge' method only */
3148 *rslt_image; /* resultant image - after multi-kernel handling */
3151 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3152 *norm_kernel, /* the current normal un-reflected kernel */
3153 *rflt_kernel, /* the current reflected kernel (if needed) */
3154 *this_kernel; /* the kernel being applied */
3157 primative; /* the current morphology primative being applied */
3160 rslt_compose; /* multi-kernel compose method for results to use */
3163 verbose; /* verbose output of results */
3166 method_loop, /* Loop 1: number of compound method iterations */
3167 method_limit, /* maximum number of compound method iterations */
3168 kernel_number, /* Loop 2: the kernel number being applied */
3169 stage_loop, /* Loop 3: primative loop for compound morphology */
3170 stage_limit, /* how many primatives in this compound */
3171 kernel_loop, /* Loop 4: iterate the kernel (basic morphology) */
3172 kernel_limit, /* number of times to iterate kernel */
3173 count, /* total count of primative steps applied */
3174 changed, /* number pixels changed by last primative operation */
3175 kernel_changed, /* total count of changed using iterated kernel */
3176 method_changed; /* total count of changed over method iteration */
3181 assert(image != (Image *) NULL);
3182 assert(image->signature == MagickSignature);
3183 assert(kernel != (KernelInfo *) NULL);
3184 assert(kernel->signature == MagickSignature);
3185 assert(exception != (ExceptionInfo *) NULL);
3186 assert(exception->signature == MagickSignature);
3188 count = 0; /* number of low-level morphology primatives performed */
3189 if ( iterations == 0 )
3190 return((Image *)NULL); /* null operation - nothing to do! */
3192 kernel_limit = (size_t) iterations;
3193 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3194 kernel_limit = image->columns > image->rows ? image->columns : image->rows;
3196 verbose = ( GetImageArtifact(image,"verbose") != (const char *) NULL ) ?
3197 MagickTrue : MagickFalse;
3199 /* initialise for cleanup */
3200 curr_image = (Image *) image;
3201 work_image = save_image = rslt_image = (Image *) NULL;
3202 reflected_kernel = (KernelInfo *) NULL;
3204 /* Initialize specific methods
3205 * + which loop should use the given iteratations
3206 * + how many primatives make up the compound morphology
3207 * + multi-kernel compose method to use (by default)
3209 method_limit = 1; /* just do method once, unless otherwise set */
3210 stage_limit = 1; /* assume method is not a compount */
3211 rslt_compose = compose; /* and we are composing multi-kernels as given */
3213 case SmoothMorphology: /* 4 primative compound morphology */
3216 case OpenMorphology: /* 2 primative compound morphology */
3217 case OpenIntensityMorphology:
3218 case TopHatMorphology:
3219 case CloseMorphology:
3220 case CloseIntensityMorphology:
3221 case BottomHatMorphology:
3222 case EdgeMorphology:
3225 case HitAndMissMorphology:
3226 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3228 case ThinningMorphology:
3229 case ThickenMorphology:
3230 method_limit = kernel_limit; /* iterate the whole method */
3231 kernel_limit = 1; /* do not do kernel iteration */
3237 /* Handle user (caller) specified multi-kernel composition method */
3238 if ( compose != UndefinedCompositeOp )
3239 rslt_compose = compose; /* override default composition for method */
3240 if ( rslt_compose == UndefinedCompositeOp )
3241 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3243 /* Some methods require a reflected kernel to use with primatives.
3244 * Create the reflected kernel for those methods. */
3246 case CorrelateMorphology:
3247 case CloseMorphology:
3248 case CloseIntensityMorphology:
3249 case BottomHatMorphology:
3250 case SmoothMorphology:
3251 reflected_kernel = CloneKernelInfo(kernel);
3252 if (reflected_kernel == (KernelInfo *) NULL)
3254 RotateKernelInfo(reflected_kernel,180);
3260 /* Loop 1: iterate the compound method */
3263 while ( method_loop < method_limit && method_changed > 0 ) {
3267 /* Loop 2: iterate over each kernel in a multi-kernel list */
3268 norm_kernel = (KernelInfo *) kernel;
3269 this_kernel = (KernelInfo *) kernel;
3270 rflt_kernel = reflected_kernel;
3273 while ( norm_kernel != NULL ) {
3275 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3276 stage_loop = 0; /* the compound morphology stage number */
3277 while ( stage_loop < stage_limit ) {
3278 stage_loop++; /* The stage of the compound morphology */
3280 /* Select primative morphology for this stage of compound method */
3281 this_kernel = norm_kernel; /* default use unreflected kernel */
3282 primative = method; /* Assume method is a primative */
3284 case ErodeMorphology: /* just erode */
3285 case EdgeInMorphology: /* erode and image difference */
3286 primative = ErodeMorphology;
3288 case DilateMorphology: /* just dilate */
3289 case EdgeOutMorphology: /* dilate and image difference */
3290 primative = DilateMorphology;
3292 case OpenMorphology: /* erode then dialate */
3293 case TopHatMorphology: /* open and image difference */
3294 primative = ErodeMorphology;
3295 if ( stage_loop == 2 )
3296 primative = DilateMorphology;
3298 case OpenIntensityMorphology:
3299 primative = ErodeIntensityMorphology;
3300 if ( stage_loop == 2 )
3301 primative = DilateIntensityMorphology;
3303 case CloseMorphology: /* dilate, then erode */
3304 case BottomHatMorphology: /* close and image difference */
3305 this_kernel = rflt_kernel; /* use the reflected kernel */
3306 primative = DilateMorphology;
3307 if ( stage_loop == 2 )
3308 primative = ErodeMorphology;
3310 case CloseIntensityMorphology:
3311 this_kernel = rflt_kernel; /* use the reflected kernel */
3312 primative = DilateIntensityMorphology;
3313 if ( stage_loop == 2 )
3314 primative = ErodeIntensityMorphology;
3316 case SmoothMorphology: /* open, close */
3317 switch ( stage_loop ) {
3318 case 1: /* start an open method, which starts with Erode */
3319 primative = ErodeMorphology;
3321 case 2: /* now Dilate the Erode */
3322 primative = DilateMorphology;
3324 case 3: /* Reflect kernel a close */
3325 this_kernel = rflt_kernel; /* use the reflected kernel */
3326 primative = DilateMorphology;
3328 case 4: /* Finish the Close */
3329 this_kernel = rflt_kernel; /* use the reflected kernel */
3330 primative = ErodeMorphology;
3334 case EdgeMorphology: /* dilate and erode difference */
3335 primative = DilateMorphology;
3336 if ( stage_loop == 2 ) {
3337 save_image = curr_image; /* save the image difference */
3338 curr_image = (Image *) image;
3339 primative = ErodeMorphology;
3342 case CorrelateMorphology:
3343 /* A Correlation is a Convolution with a reflected kernel.
3344 ** However a Convolution is a weighted sum using a reflected
3345 ** kernel. It may seem stange to convert a Correlation into a
3346 ** Convolution as the Correlation is the simplier method, but
3347 ** Convolution is much more commonly used, and it makes sense to
3348 ** implement it directly so as to avoid the need to duplicate the
3349 ** kernel when it is not required (which is typically the
3352 this_kernel = rflt_kernel; /* use the reflected kernel */
3353 primative = ConvolveMorphology;
3358 assert( this_kernel != (KernelInfo *) NULL );
3360 /* Extra information for debugging compound operations */
3361 if ( verbose == MagickTrue ) {
3362 if ( stage_limit > 1 )
3363 (void) FormatMagickString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
3364 MagickOptionToMnemonic(MagickMorphologyOptions,method),(double)
3365 method_loop,(double) stage_loop);
3366 else if ( primative != method )
3367 (void) FormatMagickString(v_info, MaxTextExtent, "%s:%.20g -> ",
3368 MagickOptionToMnemonic(MagickMorphologyOptions, method),(double)
3374 /* Loop 4: Iterate the kernel with primative */
3378 while ( kernel_loop < kernel_limit && changed > 0 ) {
3379 kernel_loop++; /* the iteration of this kernel */
3381 /* Create a destination image, if not yet defined */
3382 if ( work_image == (Image *) NULL )
3384 work_image=CloneImage(image,0,0,MagickTrue,exception);
3385 if (work_image == (Image *) NULL)
3387 if (SetImageStorageClass(work_image,DirectClass) == MagickFalse)
3389 InheritException(exception,&work_image->exception);
3394 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3396 changed = MorphologyPrimitive(curr_image, work_image, primative,
3397 channel, this_kernel, bias, exception);
3398 kernel_changed += changed;
3399 method_changed += changed;
3401 if ( verbose == MagickTrue ) {
3402 if ( kernel_loop > 1 )
3403 fprintf(stderr, "\n"); /* add end-of-line from previous */
3404 (void) fprintf(stderr, "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3405 v_info,MagickOptionToMnemonic(MagickMorphologyOptions,
3406 primative),(this_kernel == rflt_kernel ) ? "*" : "",
3407 (double) (method_loop+kernel_loop-1),(double) kernel_number,
3408 (double) count,(double) changed);
3410 /* prepare next loop */
3411 { Image *tmp = work_image; /* swap images for iteration */
3412 work_image = curr_image;
3415 if ( work_image == image )
3416 work_image = (Image *) NULL; /* replace input 'image' */
3418 } /* End Loop 4: Iterate the kernel with primative */
3420 if ( verbose == MagickTrue && kernel_changed != changed )
3421 fprintf(stderr, " Total %.20g",(double) kernel_changed);
3422 if ( verbose == MagickTrue && stage_loop < stage_limit )
3423 fprintf(stderr, "\n"); /* add end-of-line before looping */
3426 fprintf(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3427 fprintf(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
3428 fprintf(stderr, " work =0x%lx\n", (unsigned long)work_image);
3429 fprintf(stderr, " save =0x%lx\n", (unsigned long)save_image);
3430 fprintf(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
3433 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3435 /* Final Post-processing for some Compound Methods
3437 ** The removal of any 'Sync' channel flag in the Image Compositon
3438 ** below ensures the methematical compose method is applied in a
3439 ** purely mathematical way, and only to the selected channels.
3440 ** Turn off SVG composition 'alpha blending'.
3443 case EdgeOutMorphology:
3444 case EdgeInMorphology:
3445 case TopHatMorphology:
3446 case BottomHatMorphology:
3447 if ( verbose == MagickTrue )
3448 fprintf(stderr, "\n%s: Difference with original image",
3449 MagickOptionToMnemonic(MagickMorphologyOptions, method) );
3450 (void) CompositeImageChannel(curr_image,
3451 (ChannelType) (channel & ~SyncChannels),
3452 DifferenceCompositeOp, image, 0, 0);
3454 case EdgeMorphology:
3455 if ( verbose == MagickTrue )
3456 fprintf(stderr, "\n%s: Difference of Dilate and Erode",
3457 MagickOptionToMnemonic(MagickMorphologyOptions, method) );
3458 (void) CompositeImageChannel(curr_image,
3459 (ChannelType) (channel & ~SyncChannels),
3460 DifferenceCompositeOp, save_image, 0, 0);
3461 save_image = DestroyImage(save_image); /* finished with save image */
3467 /* multi-kernel handling: re-iterate, or compose results */
3468 if ( kernel->next == (KernelInfo *) NULL )
3469 rslt_image = curr_image; /* just return the resulting image */
3470 else if ( rslt_compose == NoCompositeOp )
3471 { if ( verbose == MagickTrue ) {
3472 if ( this_kernel->next != (KernelInfo *) NULL )
3473 fprintf(stderr, " (re-iterate)");
3475 fprintf(stderr, " (done)");
3477 rslt_image = curr_image; /* return result, and re-iterate */
3479 else if ( rslt_image == (Image *) NULL)
3480 { if ( verbose == MagickTrue )
3481 fprintf(stderr, " (save for compose)");
3482 rslt_image = curr_image;
3483 curr_image = (Image *) image; /* continue with original image */
3486 { /* add the new 'current' result to the composition
3488 ** The removal of any 'Sync' channel flag in the Image Compositon
3489 ** below ensures the methematical compose method is applied in a
3490 ** purely mathematical way, and only to the selected channels.
3491 ** Turn off SVG composition 'alpha blending'.
3493 ** The compose image order is specifically so that the new image can
3494 ** be subtarcted 'Minus' from the collected result, to allow you to
3495 ** convert a HitAndMiss methd into a Thinning method.
3497 if ( verbose == MagickTrue )
3498 fprintf(stderr, " (compose \"%s\")",
3499 MagickOptionToMnemonic(MagickComposeOptions, rslt_compose) );
3500 (void) CompositeImageChannel(curr_image,
3501 (ChannelType) (channel & ~SyncChannels), rslt_compose,
3503 rslt_image = DestroyImage(rslt_image);
3504 rslt_image = curr_image;
3505 curr_image = (Image *) image; /* continue with original image */
3507 if ( verbose == MagickTrue )
3508 fprintf(stderr, "\n");
3510 /* loop to the next kernel in a multi-kernel list */
3511 norm_kernel = norm_kernel->next;
3512 if ( rflt_kernel != (KernelInfo *) NULL )
3513 rflt_kernel = rflt_kernel->next;
3515 } /* End Loop 2: Loop over each kernel */
3517 } /* End Loop 1: compound method interation */
3521 /* Yes goto's are bad, but it makes cleanup lot more efficient */
3523 if ( curr_image != (Image *) NULL &&
3524 curr_image != rslt_image &&
3525 curr_image != image )
3526 curr_image = DestroyImage(curr_image);
3527 if ( rslt_image != (Image *) NULL )
3528 rslt_image = DestroyImage(rslt_image);
3530 if ( curr_image != (Image *) NULL &&
3531 curr_image != rslt_image &&
3532 curr_image != image )
3533 curr_image = DestroyImage(curr_image);
3534 if ( work_image != (Image *) NULL )
3535 work_image = DestroyImage(work_image);
3536 if ( save_image != (Image *) NULL )
3537 save_image = DestroyImage(save_image);
3538 if ( reflected_kernel != (KernelInfo *) NULL )
3539 reflected_kernel = DestroyKernelInfo(reflected_kernel);
3544 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3548 % M o r p h o l o g y I m a g e C h a n n e l %
3552 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3554 % MorphologyImageChannel() applies a user supplied kernel to the image
3555 % according to the given mophology method.
3557 % This function applies any and all user defined settings before calling
3558 % the above internal function MorphologyApply().
3560 % User defined settings include...
3561 % * Output Bias for Convolution and correlation ("-bias")
3562 % * Kernel Scale/normalize settings ("-set 'option:convolve:scale'")
3563 % This can also includes the addition of a scaled unity kernel.
3564 % * Show Kernel being applied ("-set option:showkernel 1")
3566 % The format of the MorphologyImage method is:
3568 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
3569 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
3571 % Image *MorphologyImageChannel(const Image *image, const ChannelType
3572 % channel,MorphologyMethod method,const ssize_t iterations,
3573 % KernelInfo *kernel,ExceptionInfo *exception)
3575 % A description of each parameter follows:
3577 % o image: the image.
3579 % o method: the morphology method to be applied.
3581 % o iterations: apply the operation this many times (or no change).
3582 % A value of -1 means loop until no change found.
3583 % How this is applied may depend on the morphology method.
3584 % Typically this is a value of 1.
3586 % o channel: the channel type.
3588 % o kernel: An array of double representing the morphology kernel.
3589 % Warning: kernel may be normalized for the Convolve method.
3591 % o exception: return any errors or warnings in this structure.
3595 MagickExport Image *MorphologyImageChannel(const Image *image,
3596 const ChannelType channel,const MorphologyMethod method,
3597 const ssize_t iterations,const KernelInfo *kernel,ExceptionInfo *exception)
3612 /* Apply Convolve/Correlate Normalization and Scaling Factors.
3613 * This is done BEFORE the ShowKernelInfo() function is called so that
3614 * users can see the results of the 'option:convolve:scale' option.
3616 curr_kernel = (KernelInfo *) kernel;
3617 if ( method == ConvolveMorphology || method == CorrelateMorphology )
3619 artifact = GetImageArtifact(image,"convolve:scale");
3620 if ( artifact != (const char *)NULL ) {
3621 if ( curr_kernel == kernel )
3622 curr_kernel = CloneKernelInfo(kernel);
3623 if (curr_kernel == (KernelInfo *) NULL) {
3624 curr_kernel=DestroyKernelInfo(curr_kernel);
3625 return((Image *) NULL);
3627 ScaleGeometryKernelInfo(curr_kernel, artifact);
3631 /* display the (normalized) kernel via stderr */
3632 artifact = GetImageArtifact(image,"showkernel");
3633 if ( artifact == (const char *) NULL)
3634 artifact = GetImageArtifact(image,"convolve:showkernel");
3635 if ( artifact == (const char *) NULL)
3636 artifact = GetImageArtifact(image,"morphology:showkernel");
3637 if ( artifact != (const char *) NULL)
3638 ShowKernelInfo(curr_kernel);
3640 /* Override the default handling of multi-kernel morphology results
3641 * If 'Undefined' use the default method
3642 * If 'None' (default for 'Convolve') re-iterate previous result
3643 * Otherwise merge resulting images using compose method given.
3644 * Default for 'HitAndMiss' is 'Lighten'.
3646 compose = UndefinedCompositeOp; /* use default for method */
3647 artifact = GetImageArtifact(image,"morphology:compose");
3648 if ( artifact != (const char *) NULL)
3649 compose = (CompositeOperator) ParseMagickOption(
3650 MagickComposeOptions,MagickFalse,artifact);
3652 /* Apply the Morphology */
3653 morphology_image = MorphologyApply(image, channel, method, iterations,
3654 curr_kernel, compose, image->bias, exception);
3656 /* Cleanup and Exit */
3657 if ( curr_kernel != kernel )
3658 curr_kernel=DestroyKernelInfo(curr_kernel);
3659 return(morphology_image);
3662 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
3663 method, const ssize_t iterations,const KernelInfo *kernel, ExceptionInfo
3669 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
3670 iterations,kernel,exception);
3671 return(morphology_image);
3675 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3679 + R o t a t e K e r n e l I n f o %
3683 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3685 % RotateKernelInfo() rotates the kernel by the angle given.
3687 % Currently it is restricted to 90 degree angles, of either 1D kernels
3688 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
3689 % It will ignore usless rotations for specific 'named' built-in kernels.
3691 % The format of the RotateKernelInfo method is:
3693 % void RotateKernelInfo(KernelInfo *kernel, double angle)
3695 % A description of each parameter follows:
3697 % o kernel: the Morphology/Convolution kernel
3699 % o angle: angle to rotate in degrees
3701 % This function is currently internal to this module only, but can be exported
3702 % to other modules if needed.
3704 static void RotateKernelInfo(KernelInfo *kernel, double angle)
3706 /* angle the lower kernels first */
3707 if ( kernel->next != (KernelInfo *) NULL)
3708 RotateKernelInfo(kernel->next, angle);
3710 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
3712 ** TODO: expand beyond simple 90 degree rotates, flips and flops
3715 /* Modulus the angle */
3716 angle = fmod(angle, 360.0);
3720 if ( 337.5 < angle || angle <= 22.5 )
3721 return; /* Near zero angle - no change! - At least not at this time */
3723 /* Handle special cases */
3724 switch (kernel->type) {
3725 /* These built-in kernels are cylindrical kernels, rotating is useless */
3726 case GaussianKernel:
3731 case LaplacianKernel:
3732 case ChebyshevKernel:
3733 case ManhattanKernel:
3734 case EuclideanKernel:
3737 /* These may be rotatable at non-90 angles in the future */
3738 /* but simply rotating them in multiples of 90 degrees is useless */
3745 /* These only allows a +/-90 degree rotation (by transpose) */
3746 /* A 180 degree rotation is useless */
3748 case RectangleKernel:
3749 if ( 135.0 < angle && angle <= 225.0 )
3751 if ( 225.0 < angle && angle <= 315.0 )
3758 /* Attempt rotations by 45 degrees */
3759 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
3761 if ( kernel->width == 3 && kernel->height == 3 )
3762 { /* Rotate a 3x3 square by 45 degree angle */
3763 MagickRealType t = kernel->values[0];
3764 kernel->values[0] = kernel->values[3];
3765 kernel->values[3] = kernel->values[6];
3766 kernel->values[6] = kernel->values[7];
3767 kernel->values[7] = kernel->values[8];
3768 kernel->values[8] = kernel->values[5];
3769 kernel->values[5] = kernel->values[2];
3770 kernel->values[2] = kernel->values[1];
3771 kernel->values[1] = t;
3772 /* rotate non-centered origin */
3773 if ( kernel->x != 1 || kernel->y != 1 ) {
3775 x = (ssize_t) kernel->x-1;
3776 y = (ssize_t) kernel->y-1;
3777 if ( x == y ) x = 0;
3778 else if ( x == 0 ) x = -y;
3779 else if ( x == -y ) y = 0;
3780 else if ( y == 0 ) y = x;
3781 kernel->x = (ssize_t) x+1;
3782 kernel->y = (ssize_t) y+1;
3784 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
3785 kernel->angle = fmod(kernel->angle+45.0, 360.0);
3788 perror("Unable to rotate non-3x3 kernel by 45 degrees");
3790 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
3792 if ( kernel->width == 1 || kernel->height == 1 )
3793 { /* Do a transpose of a 1 dimentional kernel,
3794 ** which results in a fast 90 degree rotation of some type.
3798 t = (ssize_t) kernel->width;
3799 kernel->width = kernel->height;
3800 kernel->height = (size_t) t;
3802 kernel->x = kernel->y;
3804 if ( kernel->width == 1 ) {
3805 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
3806 kernel->angle = fmod(kernel->angle+90.0, 360.0);
3808 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
3809 kernel->angle = fmod(kernel->angle+270.0, 360.0);
3812 else if ( kernel->width == kernel->height )
3813 { /* Rotate a square array of values by 90 degrees */
3816 register MagickRealType
3819 for( i=0, x=kernel->width-1; i<=x; i++, x--)
3820 for( j=0, y=kernel->height-1; j<y; j++, y--)
3821 { t = k[i+j*kernel->width];
3822 k[i+j*kernel->width] = k[j+x*kernel->width];
3823 k[j+x*kernel->width] = k[x+y*kernel->width];
3824 k[x+y*kernel->width] = k[y+i*kernel->width];
3825 k[y+i*kernel->width] = t;
3828 /* rotate the origin - relative to center of array */
3829 { register ssize_t x,y;
3830 x = (ssize_t) (kernel->x*2-kernel->width+1);
3831 y = (ssize_t) (kernel->y*2-kernel->height+1);
3832 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
3833 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
3835 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
3836 kernel->angle = fmod(kernel->angle+90.0, 360.0);
3839 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
3841 if ( 135.0 < angle && angle <= 225.0 )
3843 /* For a 180 degree rotation - also know as a reflection
3844 * This is actually a very very common operation!
3845 * Basically all that is needed is a reversal of the kernel data!
3846 * And a reflection of the origon
3854 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
3855 t=k[i], k[i]=k[j], k[j]=t;
3857 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
3858 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
3859 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
3860 kernel->angle = fmod(kernel->angle+180.0, 360.0);
3862 /* At this point angle should at least between -45 (315) and +45 degrees
3863 * In the future some form of non-orthogonal angled rotates could be
3864 * performed here, posibily with a linear kernel restriction.
3871 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3875 % S c a l e G e o m e t r y K e r n e l I n f o %
3879 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3881 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
3882 % provided as a "-set option:convolve:scale {geometry}" user setting,
3883 % and modifies the kernel according to the parsed arguments of that setting.
3885 % The first argument (and any normalization flags) are passed to
3886 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
3887 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
3888 % into the scaled/normalized kernel.
3890 % The format of the ScaleKernelInfo method is:
3892 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
3893 % const MagickStatusType normalize_flags )
3895 % A description of each parameter follows:
3897 % o kernel: the Morphology/Convolution kernel to modify
3900 % The geometry string to parse, typically from the user provided
3901 % "-set option:convolve:scale {geometry}" setting.
3904 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
3905 const char *geometry)
3912 SetGeometryInfo(&args);
3913 flags = (GeometryFlags) ParseGeometry(geometry, &args);
3916 /* For Debugging Geometry Input */
3917 fprintf(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
3918 flags, args.rho, args.sigma, args.xi, args.psi );
3921 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
3922 args.rho *= 0.01, args.sigma *= 0.01;
3924 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
3926 if ( (flags & SigmaValue) == 0 )
3929 /* Scale/Normalize the input kernel */
3930 ScaleKernelInfo(kernel, args.rho, flags);
3932 /* Add Unity Kernel, for blending with original */
3933 if ( (flags & SigmaValue) != 0 )
3934 UnityAddKernelInfo(kernel, args.sigma);
3939 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3943 % S c a l e K e r n e l I n f o %
3947 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3949 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
3950 % without normalization of the sum of the kernel values (as per given flags).
3952 % By default (no flags given) the values within the kernel is scaled
3953 % directly using given scaling factor without change.
3955 % If either of the two 'normalize_flags' are given the kernel will first be
3956 % normalized and then further scaled by the scaling factor value given.
3958 % Kernel normalization ('normalize_flags' given) is designed to ensure that
3959 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
3960 % morphology methods will fall into -1.0 to +1.0 range. Note that for
3961 % non-HDRI versions of IM this may cause images to have any negative results
3962 % clipped, unless some 'bias' is used.
3964 % More specifically. Kernels which only contain positive values (such as a
3965 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
3966 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
3968 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
3969 % the kernel will be scaled by the absolute of the sum of kernel values, so
3970 % that it will generally fall within the +/- 1.0 range.
3972 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
3973 % will be scaled by just the sum of the postive values, so that its output
3974 % range will again fall into the +/- 1.0 range.
3976 % For special kernels designed for locating shapes using 'Correlate', (often
3977 % only containing +1 and -1 values, representing foreground/brackground
3978 % matching) a special normalization method is provided to scale the positive
3979 % values seperatally to those of the negative values, so the kernel will be
3980 % forced to become a zero-sum kernel better suited to such searches.
3982 % WARNING: Correct normalization of the kernel assumes that the '*_range'
3983 % attributes within the kernel structure have been correctly set during the
3986 % NOTE: The values used for 'normalize_flags' have been selected specifically
3987 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
3988 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
3990 % The format of the ScaleKernelInfo method is:
3992 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
3993 % const MagickStatusType normalize_flags )
3995 % A description of each parameter follows:
3997 % o kernel: the Morphology/Convolution kernel
4000 % multiply all values (after normalization) by this factor if not
4001 % zero. If the kernel is normalized regardless of any flags.
4003 % o normalize_flags:
4004 % GeometryFlags defining normalization method to use.
4005 % specifically: NormalizeValue, CorrelateNormalizeValue,
4006 % and/or PercentValue
4009 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4010 const double scaling_factor,const GeometryFlags normalize_flags)
4019 /* do the other kernels in a multi-kernel list first */
4020 if ( kernel->next != (KernelInfo *) NULL)
4021 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4023 /* Normalization of Kernel */
4025 if ( (normalize_flags&NormalizeValue) != 0 ) {
4026 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4027 /* non-zero-summing kernel (generally positive) */
4028 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4030 /* zero-summing kernel */
4031 pos_scale = kernel->positive_range;
4033 /* Force kernel into a normalized zero-summing kernel */
4034 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4035 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4036 ? kernel->positive_range : 1.0;
4037 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
4038 ? -kernel->negative_range : 1.0;
4041 neg_scale = pos_scale;
4043 /* finialize scaling_factor for positive and negative components */
4044 pos_scale = scaling_factor/pos_scale;
4045 neg_scale = scaling_factor/neg_scale;
4047 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4048 if ( ! IsNan(kernel->values[i]) )
4049 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4051 /* convolution output range */
4052 kernel->positive_range *= pos_scale;
4053 kernel->negative_range *= neg_scale;
4054 /* maximum and minimum values in kernel */
4055 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4056 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4058 /* swap kernel settings if user's scaling factor is negative */
4059 if ( scaling_factor < MagickEpsilon ) {
4061 t = kernel->positive_range;
4062 kernel->positive_range = kernel->negative_range;
4063 kernel->negative_range = t;
4064 t = kernel->maximum;
4065 kernel->maximum = kernel->minimum;
4066 kernel->minimum = 1;
4073 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4077 % S h o w K e r n e l I n f o %
4081 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4083 % ShowKernelInfo() outputs the details of the given kernel defination to
4084 % standard error, generally due to a users 'showkernel' option request.
4086 % The format of the ShowKernel method is:
4088 % void ShowKernelInfo(KernelInfo *kernel)
4090 % A description of each parameter follows:
4092 % o kernel: the Morphology/Convolution kernel
4095 MagickExport void ShowKernelInfo(KernelInfo *kernel)
4103 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4105 fprintf(stderr, "Kernel");
4106 if ( kernel->next != (KernelInfo *) NULL )
4107 fprintf(stderr, " #%lu", (unsigned long) c );
4108 fprintf(stderr, " \"%s",
4109 MagickOptionToMnemonic(MagickKernelOptions, k->type) );
4110 if ( fabs(k->angle) > MagickEpsilon )
4111 fprintf(stderr, "@%lg", k->angle);
4112 fprintf(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long) k->width,
4113 (unsigned long) k->height,(long) k->x,(long) k->y);
4115 " with values from %.*lg to %.*lg\n",
4116 GetMagickPrecision(), k->minimum,
4117 GetMagickPrecision(), k->maximum);
4118 fprintf(stderr, "Forming a output range from %.*lg to %.*lg",
4119 GetMagickPrecision(), k->negative_range,
4120 GetMagickPrecision(), k->positive_range);
4121 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4122 fprintf(stderr, " (Zero-Summing)\n");
4123 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4124 fprintf(stderr, " (Normalized)\n");
4126 fprintf(stderr, " (Sum %.*lg)\n",
4127 GetMagickPrecision(), k->positive_range+k->negative_range);
4128 for (i=v=0; v < k->height; v++) {
4129 fprintf(stderr, "%2lu:", (unsigned long) v );
4130 for (u=0; u < k->width; u++, i++)
4131 if ( IsNan(k->values[i]) )
4132 fprintf(stderr," %*s", GetMagickPrecision()+3, "nan");
4134 fprintf(stderr," %*.*lg", GetMagickPrecision()+3,
4135 GetMagickPrecision(), k->values[i]);
4136 fprintf(stderr,"\n");
4142 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4146 % U n i t y A d d K e r n a l I n f o %
4150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4152 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4153 % to the given pre-scaled and normalized Kernel. This in effect adds that
4154 % amount of the original image into the resulting convolution kernel. This
4155 % value is usually provided by the user as a percentage value in the
4156 % 'convolve:scale' setting.
4158 % The resulting effect is to convert the defined kernels into blended
4159 % soft-blurs, unsharp kernels or into sharpening kernels.
4161 % The format of the UnityAdditionKernelInfo method is:
4163 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4165 % A description of each parameter follows:
4167 % o kernel: the Morphology/Convolution kernel
4170 % scaling factor for the unity kernel to be added to
4174 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4177 /* do the other kernels in a multi-kernel list first */
4178 if ( kernel->next != (KernelInfo *) NULL)
4179 UnityAddKernelInfo(kernel->next, scale);
4181 /* Add the scaled unity kernel to the existing kernel */
4182 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4183 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4193 % Z e r o K e r n e l N a n s %
4197 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4199 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4200 % the kernel with a zero value. This is typically done when the kernel will
4201 % be used in special hardware (GPU) convolution processors, to simply
4204 % The format of the ZeroKernelNans method is:
4206 % void ZeroKernelNans (KernelInfo *kernel)
4208 % A description of each parameter follows:
4210 % o kernel: the Morphology/Convolution kernel
4213 MagickExport void ZeroKernelNans(KernelInfo *kernel)
4218 /* do the other kernels in a multi-kernel list first */
4219 if ( kernel->next != (KernelInfo *) NULL)
4220 ZeroKernelNans(kernel->next);
4222 for (i=0; i < (kernel->width*kernel->height); i++)
4223 if ( IsNan(kernel->values[i]) )
4224 kernel->values[i] = 0.0;