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-2011 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 "MagickCore/studio.h"
53 #include "MagickCore/artifact.h"
54 #include "MagickCore/cache-view.h"
55 #include "MagickCore/color-private.h"
56 #include "MagickCore/enhance.h"
57 #include "MagickCore/exception.h"
58 #include "MagickCore/exception-private.h"
59 #include "MagickCore/gem.h"
60 #include "MagickCore/hashmap.h"
61 #include "MagickCore/image.h"
62 #include "MagickCore/image-private.h"
63 #include "MagickCore/list.h"
64 #include "MagickCore/magick.h"
65 #include "MagickCore/memory_.h"
66 #include "MagickCore/monitor-private.h"
67 #include "MagickCore/morphology.h"
68 #include "MagickCore/morphology-private.h"
69 #include "MagickCore/option.h"
70 #include "MagickCore/pixel-accessor.h"
71 #include "MagickCore/prepress.h"
72 #include "MagickCore/quantize.h"
73 #include "MagickCore/registry.h"
74 #include "MagickCore/semaphore.h"
75 #include "MagickCore/splay-tree.h"
76 #include "MagickCore/statistic.h"
77 #include "MagickCore/string_.h"
78 #include "MagickCore/string-private.h"
79 #include "MagickCore/token.h"
80 #include "MagickCore/utility.h"
84 ** The following test is for special floating point numbers of value NaN (not
85 ** a number), that may be used within a Kernel Definition. NaN's are defined
86 ** as part of the IEEE standard for floating point number representation.
88 ** These are used as a Kernel value to mean that this kernel position is not
89 ** part of the kernel neighbourhood for convolution or morphology processing,
90 ** and thus should be ignored. This allows the use of 'shaped' kernels.
92 ** The special properity that two NaN's are never equal, even if they are from
93 ** the same variable allow you to test if a value is special NaN value.
95 ** This macro IsNaN() is thus is only true if the value given is NaN.
97 #define IsNan(a) ((a)!=(a))
100 Other global definitions used by module.
102 static inline double MagickMin(const double x,const double y)
104 return( x < y ? x : y);
106 static inline double MagickMax(const double x,const double y)
108 return( x > y ? x : y);
110 #define Minimize(assign,value) assign=MagickMin(assign,value)
111 #define Maximize(assign,value) assign=MagickMax(assign,value)
113 /* Currently these are only internal to this module */
115 CalcKernelMetaData(KernelInfo *),
116 ExpandMirrorKernelInfo(KernelInfo *),
117 ExpandRotateKernelInfo(KernelInfo *, const double),
118 RotateKernelInfo(KernelInfo *, double);
121 /* Quick function to find last kernel in a kernel list */
122 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
124 while (kernel->next != (KernelInfo *) NULL)
125 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 separated 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 kernel 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] = InterpretLocaleValue(token,(char **) NULL);
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)
377 token[MaxTextExtent];
395 /* Parse special 'named' kernel */
396 GetMagickToken(kernel_string,&p,token);
397 type=ParseCommandOption(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 (void) FormatLocaleFile(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 /* Shape Kernel Defaults */
425 if ( (flags & WidthValue) == 0 )
426 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
434 if ( (flags & HeightValue) == 0 )
435 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
438 if ( (flags & XValue) == 0 )
439 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
441 case RectangleKernel: /* Rectangle - set size defaults */
442 if ( (flags & WidthValue) == 0 ) /* if no width then */
443 args.rho = args.sigma; /* then width = height */
444 if ( args.rho < 1.0 ) /* if width too small */
445 args.rho = 3; /* then width = 3 */
446 if ( args.sigma < 1.0 ) /* if height too small */
447 args.sigma = args.rho; /* then height = width */
448 if ( (flags & XValue) == 0 ) /* center offset if not defined */
449 args.xi = (double)(((ssize_t)args.rho-1)/2);
450 if ( (flags & YValue) == 0 )
451 args.psi = (double)(((ssize_t)args.sigma-1)/2);
453 /* Distance Kernel Defaults */
454 case ChebyshevKernel:
455 case ManhattanKernel:
456 case OctagonalKernel:
457 case EuclideanKernel:
458 if ( (flags & HeightValue) == 0 ) /* no distance scale */
459 args.sigma = 100.0; /* default distance scaling */
460 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
461 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
462 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
463 args.sigma *= QuantumRange/100.0; /* percentage of color range */
469 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
470 if ( kernel == (KernelInfo *) NULL )
473 /* global expand to rotated kernel list - only for single kernels */
474 if ( kernel->next == (KernelInfo *) NULL ) {
475 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
476 ExpandRotateKernelInfo(kernel, 45.0);
477 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
478 ExpandRotateKernelInfo(kernel, 90.0);
479 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
480 ExpandMirrorKernelInfo(kernel);
486 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
494 token[MaxTextExtent];
506 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
508 /* ignore extra or multiple ';' kernel separators */
509 if ( *token != ';' ) {
511 /* tokens starting with alpha is a Named kernel */
512 if (isalpha((int) *token) != 0)
513 new_kernel = ParseKernelName(p);
514 else /* otherwise a user defined kernel array */
515 new_kernel = ParseKernelArray(p);
517 /* Error handling -- this is not proper error handling! */
518 if ( new_kernel == (KernelInfo *) NULL ) {
519 (void) FormatLocaleFile(stderr, "Failed to parse kernel number #%.20g\n",
520 (double) kernel_number);
521 if ( kernel != (KernelInfo *) NULL )
522 kernel=DestroyKernelInfo(kernel);
523 return((KernelInfo *) NULL);
526 /* initialise or append the kernel list */
527 if ( kernel == (KernelInfo *) NULL )
530 LastKernelInfo(kernel)->next = new_kernel;
533 /* look for the next kernel in list */
535 if ( p == (char *) NULL )
545 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
549 % A c q u i r e K e r n e l B u i l t I n %
553 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
555 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
556 % kernels used for special purposes such as gaussian blurring, skeleton
557 % pruning, and edge distance determination.
559 % They take a KernelType, and a set of geometry style arguments, which were
560 % typically decoded from a user supplied string, or from a more complex
561 % Morphology Method that was requested.
563 % The format of the AcquireKernalBuiltIn method is:
565 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
566 % const GeometryInfo args)
568 % A description of each parameter follows:
570 % o type: the pre-defined type of kernel wanted
572 % o args: arguments defining or modifying the kernel
574 % Convolution Kernels
577 % The a No-Op or Scaling single element kernel.
579 % Gaussian:{radius},{sigma}
580 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
581 % The sigma for the curve is required. The resulting kernel is
584 % If 'sigma' is zero, you get a single pixel on a field of zeros.
586 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
587 % the final size of the resulting kernel to a square 2*radius+1 in size.
588 % The radius should be at least 2 times that of the sigma value, or
589 % sever clipping and aliasing may result. If not given or set to 0 the
590 % radius will be determined so as to produce the best minimal error
591 % result, which is usally much larger than is normally needed.
593 % LoG:{radius},{sigma}
594 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
595 % The supposed ideal edge detection, zero-summing kernel.
597 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
598 % approx 1.6 (according to wikipedia).
600 % DoG:{radius},{sigma1},{sigma2}
601 % "Difference of Gaussians" Kernel.
602 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
603 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
604 % The result is a zero-summing kernel.
606 % Blur:{radius},{sigma}[,{angle}]
607 % Generates a 1 dimensional or linear gaussian blur, at the angle given
608 % (current restricted to orthogonal angles). If a 'radius' is given the
609 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
610 % by a 90 degree angle.
612 % If 'sigma' is zero, you get a single pixel on a field of zeros.
614 % Note that two convolutions with two "Blur" kernels perpendicular to
615 % each other, is equivalent to a far larger "Gaussian" kernel with the
616 % same sigma value, However it is much faster to apply. This is how the
617 % "-blur" operator actually works.
619 % Comet:{width},{sigma},{angle}
620 % Blur in one direction only, much like how a bright object leaves
621 % a comet like trail. The Kernel is actually half a gaussian curve,
622 % Adding two such blurs in opposite directions produces a Blur Kernel.
623 % Angle can be rotated in multiples of 90 degrees.
625 % Note that the first argument is the width of the kernel and not the
626 % radius of the kernel.
628 % # Still to be implemented...
632 % # Set kernel values using a resize filter, and given scale (sigma)
633 % # Cylindrical or Linear. Is this possible with an image?
636 % Named Constant Convolution Kernels
638 % All these are unscaled, zero-summing kernels by default. As such for
639 % non-HDRI version of ImageMagick some form of normalization, user scaling,
640 % and biasing the results is recommended, to prevent the resulting image
643 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
644 % 45 degrees to generate the 8 angled varients of each of the kernels.
647 % Discrete Lapacian Kernels, (without normalization)
648 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
649 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
650 % Type 2 : 3x3 with center:4 edge:1 corner:-2
651 % Type 3 : 3x3 with center:4 edge:-2 corner:1
652 % Type 5 : 5x5 laplacian
653 % Type 7 : 7x7 laplacian
654 % Type 15 : 5x5 LoG (sigma approx 1.4)
655 % Type 19 : 9x9 LoG (sigma approx 1.4)
658 % Sobel 'Edge' convolution kernel (3x3)
664 % Roberts convolution kernel (3x3)
670 % Prewitt Edge convolution kernel (3x3)
676 % Prewitt's "Compass" convolution kernel (3x3)
682 % Kirsch's "Compass" convolution kernel (3x3)
688 % Frei-Chen Edge Detector is based on a kernel that is similar to
689 % the Sobel Kernel, but is designed to be isotropic. That is it takes
690 % into account the distance of the diagonal in the kernel.
693 % | sqrt(2), 0, -sqrt(2) |
696 % FreiChen:{type},{angle}
698 % Frei-Chen Pre-weighted kernels...
700 % Type 0: default un-nomalized version shown above.
702 % Type 1: Orthogonal Kernel (same as type 11 below)
704 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
707 % Type 2: Diagonal form of Kernel...
709 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
712 % However this kernel is als at the heart of the FreiChen Edge Detection
713 % Process which uses a set of 9 specially weighted kernel. These 9
714 % kernels not be normalized, but directly applied to the image. The
715 % results is then added together, to produce the intensity of an edge in
716 % a specific direction. The square root of the pixel value can then be
717 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
718 % from each other, both the direction and the strength of the edge can be
721 % Type 10: All 9 of the following pre-weighted kernels...
723 % Type 11: | 1, 0, -1 |
724 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
727 % Type 12: | 1, sqrt(2), 1 |
728 % | 0, 0, 0 | / 2*sqrt(2)
731 % Type 13: | sqrt(2), -1, 0 |
732 % | -1, 0, 1 | / 2*sqrt(2)
735 % Type 14: | 0, 1, -sqrt(2) |
736 % | -1, 0, 1 | / 2*sqrt(2)
739 % Type 15: | 0, -1, 0 |
743 % Type 16: | 1, 0, -1 |
747 % Type 17: | 1, -2, 1 |
751 % Type 18: | -2, 1, -2 |
755 % Type 19: | 1, 1, 1 |
759 % The first 4 are for edge detection, the next 4 are for line detection
760 % and the last is to add a average component to the results.
762 % Using a special type of '-1' will return all 9 pre-weighted kernels
763 % as a multi-kernel list, so that you can use them directly (without
764 % normalization) with the special "-set option:morphology:compose Plus"
765 % setting to apply the full FreiChen Edge Detection Technique.
767 % If 'type' is large it will be taken to be an actual rotation angle for
768 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
769 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
771 % WARNING: The above was layed out as per
772 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
773 % But rotated 90 degrees so direction is from left rather than the top.
774 % I have yet to find any secondary confirmation of the above. The only
775 % other source found was actual source code at
776 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
777 % Neigher paper defineds the kernels in a way that looks locical or
778 % correct when taken as a whole.
782 % Diamond:[{radius}[,{scale}]]
783 % Generate a diamond shaped kernel with given radius to the points.
784 % Kernel size will again be radius*2+1 square and defaults to radius 1,
785 % generating a 3x3 kernel that is slightly larger than a square.
787 % Square:[{radius}[,{scale}]]
788 % Generate a square shaped kernel of size radius*2+1, and defaulting
789 % to a 3x3 (radius 1).
791 % Octagon:[{radius}[,{scale}]]
792 % Generate octagonal shaped kernel of given radius and constant scale.
793 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
794 % in "Diamond" kernel.
796 % Disk:[{radius}[,{scale}]]
797 % Generate a binary disk, thresholded at the radius given, the radius
798 % may be a float-point value. Final Kernel size is floor(radius)*2+1
799 % square. A radius of 5.3 is the default.
801 % NOTE: That a low radii Disk kernels produce the same results as
802 % many of the previously defined kernels, but differ greatly at larger
803 % radii. Here is a table of equivalences...
804 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
805 % "Disk:1.5" => "Square"
806 % "Disk:2" => "Diamond:2"
807 % "Disk:2.5" => "Octagon"
808 % "Disk:2.9" => "Square:2"
809 % "Disk:3.5" => "Octagon:3"
810 % "Disk:4.5" => "Octagon:4"
811 % "Disk:5.4" => "Octagon:5"
812 % "Disk:6.4" => "Octagon:6"
813 % All other Disk shapes are unique to this kernel, but because a "Disk"
814 % is more circular when using a larger radius, using a larger radius is
815 % preferred over iterating the morphological operation.
817 % Rectangle:{geometry}
818 % Simply generate a rectangle of 1's with the size given. You can also
819 % specify the location of the 'control point', otherwise the closest
820 % pixel to the center of the rectangle is selected.
822 % Properly centered and odd sized rectangles work the best.
824 % Symbol Dilation Kernels
826 % These kernel is not a good general morphological kernel, but is used
827 % more for highlighting and marking any single pixels in an image using,
828 % a "Dilate" method as appropriate.
830 % For the same reasons iterating these kernels does not produce the
831 % same result as using a larger radius for the symbol.
833 % Plus:[{radius}[,{scale}]]
834 % Cross:[{radius}[,{scale}]]
835 % Generate a kernel in the shape of a 'plus' or a 'cross' with
836 % a each arm the length of the given radius (default 2).
838 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
840 % Ring:{radius1},{radius2}[,{scale}]
841 % A ring of the values given that falls between the two radii.
842 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
843 % This is the 'edge' pixels of the default "Disk" kernel,
844 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
846 % Hit and Miss Kernels
848 % Peak:radius1,radius2
849 % Find any peak larger than the pixels the fall between the two radii.
850 % The default ring of pixels is as per "Ring".
852 % Find flat orthogonal edges of a binary shape
854 % Find 90 degree corners of a binary shape
856 % A special kernel to thin the 'outside' of diagonals
858 % Find end points of lines (for pruning a skeletion)
859 % Two types of lines ends (default to both) can be searched for
860 % Type 0: All line ends
861 % Type 1: single kernel for 4-conneected line ends
862 % Type 2: single kernel for simple line ends
864 % Find three line junctions (within a skeletion)
865 % Type 0: all line junctions
866 % Type 1: Y Junction kernel
867 % Type 2: Diagonal T Junction kernel
868 % Type 3: Orthogonal T Junction kernel
869 % Type 4: Diagonal X Junction kernel
870 % Type 5: Orthogonal + Junction kernel
872 % Find single pixel ridges or thin lines
873 % Type 1: Fine single pixel thick lines and ridges
874 % Type 2: Find two pixel thick lines and ridges
876 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
878 % Traditional skeleton generating kernels.
879 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
880 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
881 % Type 3: Thinning skeleton based on a ressearch paper by
882 % Dan S. Bloomberg (Default Type)
884 % A huge variety of Thinning Kernels designed to preserve conectivity.
885 % many other kernel sets use these kernels as source definitions.
886 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
887 % the super and sub notations used in the source research paper.
889 % Distance Measuring Kernels
891 % Different types of distance measuring methods, which are used with the
892 % a 'Distance' morphology method for generating a gradient based on
893 % distance from an edge of a binary shape, though there is a technique
894 % for handling a anti-aliased shape.
896 % See the 'Distance' Morphological Method, for information of how it is
899 % Chebyshev:[{radius}][x{scale}[%!]]
900 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
901 % is a value of one to any neighbour, orthogonal or diagonal. One why
902 % of thinking of it is the number of squares a 'King' or 'Queen' in
903 % chess needs to traverse reach any other position on a chess board.
904 % It results in a 'square' like distance function, but one where
905 % diagonals are given a value that is closer than expected.
907 % Manhattan:[{radius}][x{scale}[%!]]
908 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
909 % Cab distance metric), it is the distance needed when you can only
910 % travel in horizontal or vertical directions only. It is the
911 % distance a 'Rook' in chess would have to travel, and results in a
912 % diamond like distances, where diagonals are further than expected.
914 % Octagonal:[{radius}][x{scale}[%!]]
915 % An interleving of Manhatten and Chebyshev metrics producing an
916 % increasing octagonally shaped distance. Distances matches those of
917 % the "Octagon" shaped kernel of the same radius. The minimum radius
918 % and default is 2, producing a 5x5 kernel.
920 % Euclidean:[{radius}][x{scale}[%!]]
921 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
922 % However by default the kernel size only has a radius of 1, which
923 % limits the distance to 'Knight' like moves, with only orthogonal and
924 % diagonal measurements being correct. As such for the default kernel
925 % you will get octagonal like distance function.
927 % However using a larger radius such as "Euclidean:4" you will get a
928 % much smoother distance gradient from the edge of the shape. Especially
929 % if the image is pre-processed to include any anti-aliasing pixels.
930 % Of course a larger kernel is slower to use, and not always needed.
932 % The first three Distance Measuring Kernels will only generate distances
933 % of exact multiples of {scale} in binary images. As such you can use a
934 % scale of 1 without loosing any information. However you also need some
935 % scaling when handling non-binary anti-aliased shapes.
937 % The "Euclidean" Distance Kernel however does generate a non-integer
938 % fractional results, and as such scaling is vital even for binary shapes.
942 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
943 const GeometryInfo *args)
956 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
958 /* Generate a new empty kernel if needed */
959 kernel=(KernelInfo *) NULL;
961 case UndefinedKernel: /* These should not call this function */
962 case UserDefinedKernel:
963 assert("Should not call this function" != (char *)NULL);
965 case LaplacianKernel: /* Named Descrete Convolution Kernels */
966 case SobelKernel: /* these are defined using other kernels */
972 case EdgesKernel: /* Hit and Miss kernels */
974 case DiagonalsKernel:
976 case LineJunctionsKernel:
978 case ConvexHullKernel:
981 break; /* A pre-generated kernel is not needed */
983 /* set to 1 to do a compile-time check that we haven't missed anything */
992 case RectangleKernel:
999 case ChebyshevKernel:
1000 case ManhattanKernel:
1001 case OctangonalKernel:
1002 case EuclideanKernel:
1006 /* Generate the base Kernel Structure */
1007 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1008 if (kernel == (KernelInfo *) NULL)
1010 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1011 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1012 kernel->negative_range = kernel->positive_range = 0.0;
1013 kernel->type = type;
1014 kernel->next = (KernelInfo *) NULL;
1015 kernel->signature = MagickSignature;
1025 kernel->height = kernel->width = (size_t) 1;
1026 kernel->x = kernel->y = (ssize_t) 0;
1027 kernel->values=(double *) AcquireQuantumMemory(1,sizeof(double));
1028 if (kernel->values == (double *) NULL)
1029 return(DestroyKernelInfo(kernel));
1030 kernel->maximum = kernel->values[0] = args->rho;
1034 case GaussianKernel:
1038 sigma = fabs(args->sigma),
1039 sigma2 = fabs(args->xi),
1042 if ( args->rho >= 1.0 )
1043 kernel->width = (size_t)args->rho*2+1;
1044 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1045 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1047 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1048 kernel->height = kernel->width;
1049 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1050 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1051 kernel->height*sizeof(double));
1052 if (kernel->values == (double *) NULL)
1053 return(DestroyKernelInfo(kernel));
1055 /* WARNING: The following generates a 'sampled gaussian' kernel.
1056 * What we really want is a 'discrete gaussian' kernel.
1058 * How to do this is I don't know, but appears to be basied on the
1059 * Error Function 'erf()' (intergral of a gaussian)
1062 if ( type == GaussianKernel || type == DoGKernel )
1063 { /* Calculate a Gaussian, OR positive half of a DoG */
1064 if ( sigma > MagickEpsilon )
1065 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1066 B = (double) (1.0/(Magick2PI*sigma*sigma));
1067 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1068 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1069 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1071 else /* limiting case - a unity (normalized Dirac) kernel */
1072 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1073 kernel->width*kernel->height*sizeof(double));
1074 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1078 if ( type == DoGKernel )
1079 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1080 if ( sigma2 > MagickEpsilon )
1081 { sigma = sigma2; /* simplify loop expressions */
1082 A = 1.0/(2.0*sigma*sigma);
1083 B = (double) (1.0/(Magick2PI*sigma*sigma));
1084 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1085 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1086 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1088 else /* limiting case - a unity (normalized Dirac) kernel */
1089 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1092 if ( type == LoGKernel )
1093 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1094 if ( sigma > MagickEpsilon )
1095 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1096 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1097 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1098 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1099 { R = ((double)(u*u+v*v))*A;
1100 kernel->values[i] = (1-R)*exp(-R)*B;
1103 else /* special case - generate a unity kernel */
1104 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1105 kernel->width*kernel->height*sizeof(double));
1106 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1110 /* Note the above kernels may have been 'clipped' by a user defined
1111 ** radius, producing a smaller (darker) kernel. Also for very small
1112 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1113 ** producing a very bright kernel.
1115 ** Normalization will still be needed.
1118 /* Normalize the 2D Gaussian Kernel
1120 ** NB: a CorrelateNormalize performs a normal Normalize if
1121 ** there are no negative values.
1123 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1124 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1130 sigma = fabs(args->sigma),
1133 if ( args->rho >= 1.0 )
1134 kernel->width = (size_t)args->rho*2+1;
1136 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1138 kernel->x = (ssize_t) (kernel->width-1)/2;
1140 kernel->negative_range = kernel->positive_range = 0.0;
1141 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1142 kernel->height*sizeof(double));
1143 if (kernel->values == (double *) NULL)
1144 return(DestroyKernelInfo(kernel));
1147 #define KernelRank 3
1148 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1149 ** It generates a gaussian 3 times the width, and compresses it into
1150 ** the expected range. This produces a closer normalization of the
1151 ** resulting kernel, especially for very low sigma values.
1152 ** As such while wierd it is prefered.
1154 ** I am told this method originally came from Photoshop.
1156 ** A properly normalized curve is generated (apart from edge clipping)
1157 ** even though we later normalize the result (for edge clipping)
1158 ** to allow the correct generation of a "Difference of Blurs".
1162 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1163 (void) ResetMagickMemory(kernel->values,0, (size_t)
1164 kernel->width*kernel->height*sizeof(double));
1165 /* Calculate a Positive 1D Gaussian */
1166 if ( sigma > MagickEpsilon )
1167 { sigma *= KernelRank; /* simplify loop expressions */
1168 alpha = 1.0/(2.0*sigma*sigma);
1169 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1170 for ( u=-v; u <= v; u++) {
1171 kernel->values[(u+v)/KernelRank] +=
1172 exp(-((double)(u*u))*alpha)*beta;
1175 else /* special case - generate a unity kernel */
1176 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1178 /* Direct calculation without curve averaging */
1180 /* Calculate a Positive Gaussian */
1181 if ( sigma > MagickEpsilon )
1182 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1183 beta = 1.0/(MagickSQ2PI*sigma);
1184 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1185 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1187 else /* special case - generate a unity kernel */
1188 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1189 kernel->width*kernel->height*sizeof(double));
1190 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1193 /* Note the above kernel may have been 'clipped' by a user defined
1194 ** radius, producing a smaller (darker) kernel. Also for very small
1195 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1196 ** producing a very bright kernel.
1198 ** Normalization will still be needed.
1201 /* Normalize the 1D Gaussian Kernel
1203 ** NB: a CorrelateNormalize performs a normal Normalize if
1204 ** there are no negative values.
1206 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1207 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1209 /* rotate the 1D kernel by given angle */
1210 RotateKernelInfo(kernel, args->xi );
1215 sigma = fabs(args->sigma),
1218 if ( args->rho < 1.0 )
1219 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1221 kernel->width = (size_t)args->rho;
1222 kernel->x = kernel->y = 0;
1224 kernel->negative_range = kernel->positive_range = 0.0;
1225 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1226 kernel->height*sizeof(double));
1227 if (kernel->values == (double *) NULL)
1228 return(DestroyKernelInfo(kernel));
1230 /* A comet blur is half a 1D gaussian curve, so that the object is
1231 ** blurred in one direction only. This may not be quite the right
1232 ** curve to use so may change in the future. The function must be
1233 ** normalised after generation, which also resolves any clipping.
1235 ** As we are normalizing and not subtracting gaussians,
1236 ** there is no need for a divisor in the gaussian formula
1238 ** It is less comples
1240 if ( sigma > MagickEpsilon )
1243 #define KernelRank 3
1244 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1245 (void) ResetMagickMemory(kernel->values,0, (size_t)
1246 kernel->width*sizeof(double));
1247 sigma *= KernelRank; /* simplify the loop expression */
1248 A = 1.0/(2.0*sigma*sigma);
1249 /* B = 1.0/(MagickSQ2PI*sigma); */
1250 for ( u=0; u < v; u++) {
1251 kernel->values[u/KernelRank] +=
1252 exp(-((double)(u*u))*A);
1253 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1255 for (i=0; i < (ssize_t) kernel->width; i++)
1256 kernel->positive_range += kernel->values[i];
1258 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1259 /* B = 1.0/(MagickSQ2PI*sigma); */
1260 for ( i=0; i < (ssize_t) kernel->width; i++)
1261 kernel->positive_range +=
1263 exp(-((double)(i*i))*A);
1264 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1267 else /* special case - generate a unity kernel */
1268 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1269 kernel->width*kernel->height*sizeof(double));
1270 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1271 kernel->positive_range = 1.0;
1274 kernel->minimum = 0.0;
1275 kernel->maximum = kernel->values[0];
1276 kernel->negative_range = 0.0;
1278 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1279 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1284 Convolution Kernels - Well Known Named Constant Kernels
1286 case LaplacianKernel:
1287 { switch ( (int) args->rho ) {
1289 default: /* laplacian square filter -- default */
1290 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1292 case 1: /* laplacian diamond filter */
1293 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1296 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1299 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1301 case 5: /* a 5x5 laplacian */
1302 kernel=ParseKernelArray(
1303 "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");
1305 case 7: /* a 7x7 laplacian */
1306 kernel=ParseKernelArray(
1307 "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" );
1309 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1310 kernel=ParseKernelArray(
1311 "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");
1313 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1314 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1315 kernel=ParseKernelArray(
1316 "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");
1319 if (kernel == (KernelInfo *) NULL)
1321 kernel->type = type;
1325 { /* Simple Sobel Kernel */
1326 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1327 if (kernel == (KernelInfo *) NULL)
1329 kernel->type = type;
1330 RotateKernelInfo(kernel, args->rho);
1335 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1336 if (kernel == (KernelInfo *) NULL)
1338 kernel->type = type;
1339 RotateKernelInfo(kernel, args->rho);
1344 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1345 if (kernel == (KernelInfo *) NULL)
1347 kernel->type = type;
1348 RotateKernelInfo(kernel, args->rho);
1353 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1354 if (kernel == (KernelInfo *) NULL)
1356 kernel->type = type;
1357 RotateKernelInfo(kernel, args->rho);
1362 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1363 if (kernel == (KernelInfo *) NULL)
1365 kernel->type = type;
1366 RotateKernelInfo(kernel, args->rho);
1369 case FreiChenKernel:
1370 /* Direction is set to be left to right positive */
1371 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1372 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1373 { switch ( (int) args->rho ) {
1376 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1377 if (kernel == (KernelInfo *) NULL)
1379 kernel->type = type;
1380 kernel->values[3] = +MagickSQ2;
1381 kernel->values[5] = -MagickSQ2;
1382 CalcKernelMetaData(kernel); /* recalculate meta-data */
1385 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1386 if (kernel == (KernelInfo *) NULL)
1388 kernel->type = type;
1389 kernel->values[1] = kernel->values[3] = +MagickSQ2;
1390 kernel->values[5] = kernel->values[7] = -MagickSQ2;
1391 CalcKernelMetaData(kernel); /* recalculate meta-data */
1392 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1395 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1396 if (kernel == (KernelInfo *) NULL)
1401 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1402 if (kernel == (KernelInfo *) NULL)
1404 kernel->type = type;
1405 kernel->values[3] = +MagickSQ2;
1406 kernel->values[5] = -MagickSQ2;
1407 CalcKernelMetaData(kernel); /* recalculate meta-data */
1408 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1411 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1412 if (kernel == (KernelInfo *) NULL)
1414 kernel->type = type;
1415 kernel->values[1] = +MagickSQ2;
1416 kernel->values[7] = +MagickSQ2;
1417 CalcKernelMetaData(kernel);
1418 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1421 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1422 if (kernel == (KernelInfo *) NULL)
1424 kernel->type = type;
1425 kernel->values[0] = +MagickSQ2;
1426 kernel->values[8] = -MagickSQ2;
1427 CalcKernelMetaData(kernel);
1428 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1431 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1432 if (kernel == (KernelInfo *) NULL)
1434 kernel->type = type;
1435 kernel->values[2] = -MagickSQ2;
1436 kernel->values[6] = +MagickSQ2;
1437 CalcKernelMetaData(kernel);
1438 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1441 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1442 if (kernel == (KernelInfo *) NULL)
1444 kernel->type = type;
1445 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1448 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1449 if (kernel == (KernelInfo *) NULL)
1451 kernel->type = type;
1452 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1455 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1456 if (kernel == (KernelInfo *) NULL)
1458 kernel->type = type;
1459 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1462 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1463 if (kernel == (KernelInfo *) NULL)
1465 kernel->type = type;
1466 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1469 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1470 if (kernel == (KernelInfo *) NULL)
1472 kernel->type = type;
1473 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1476 if ( fabs(args->sigma) > MagickEpsilon )
1477 /* Rotate by correctly supplied 'angle' */
1478 RotateKernelInfo(kernel, args->sigma);
1479 else if ( args->rho > 30.0 || args->rho < -30.0 )
1480 /* Rotate by out of bounds 'type' */
1481 RotateKernelInfo(kernel, args->rho);
1486 Boolean or Shaped Kernels
1490 if (args->rho < 1.0)
1491 kernel->width = kernel->height = 3; /* default radius = 1 */
1493 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1494 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1496 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1497 kernel->height*sizeof(double));
1498 if (kernel->values == (double *) NULL)
1499 return(DestroyKernelInfo(kernel));
1501 /* set all kernel values within diamond area to scale given */
1502 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1503 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1504 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1505 kernel->positive_range += kernel->values[i] = args->sigma;
1507 kernel->values[i] = nan;
1508 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1512 case RectangleKernel:
1515 if ( type == SquareKernel )
1517 if (args->rho < 1.0)
1518 kernel->width = kernel->height = 3; /* default radius = 1 */
1520 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1521 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1522 scale = args->sigma;
1525 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1526 if ( args->rho < 1.0 || args->sigma < 1.0 )
1527 return(DestroyKernelInfo(kernel)); /* invalid args given */
1528 kernel->width = (size_t)args->rho;
1529 kernel->height = (size_t)args->sigma;
1530 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1531 args->psi < 0.0 || args->psi > (double)kernel->height )
1532 return(DestroyKernelInfo(kernel)); /* invalid args given */
1533 kernel->x = (ssize_t) args->xi;
1534 kernel->y = (ssize_t) args->psi;
1537 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1538 kernel->height*sizeof(double));
1539 if (kernel->values == (double *) NULL)
1540 return(DestroyKernelInfo(kernel));
1542 /* set all kernel values to scale given */
1543 u=(ssize_t) (kernel->width*kernel->height);
1544 for ( i=0; i < u; i++)
1545 kernel->values[i] = scale;
1546 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1547 kernel->positive_range = scale*u;
1552 if (args->rho < 1.0)
1553 kernel->width = kernel->height = 5; /* default radius = 2 */
1555 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1556 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1558 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1559 kernel->height*sizeof(double));
1560 if (kernel->values == (double *) NULL)
1561 return(DestroyKernelInfo(kernel));
1563 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1564 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1565 if ( (labs((long) u)+labs((long) v)) <=
1566 ((long)kernel->x + (long)(kernel->x/2)) )
1567 kernel->positive_range += kernel->values[i] = args->sigma;
1569 kernel->values[i] = nan;
1570 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1576 limit = (ssize_t)(args->rho*args->rho);
1578 if (args->rho < 0.4) /* default radius approx 4.3 */
1579 kernel->width = kernel->height = 9L, limit = 18L;
1581 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1582 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1584 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1585 kernel->height*sizeof(double));
1586 if (kernel->values == (double *) NULL)
1587 return(DestroyKernelInfo(kernel));
1589 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1590 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1591 if ((u*u+v*v) <= limit)
1592 kernel->positive_range += kernel->values[i] = args->sigma;
1594 kernel->values[i] = nan;
1595 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1600 if (args->rho < 1.0)
1601 kernel->width = kernel->height = 5; /* default radius 2 */
1603 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1604 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1606 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1607 kernel->height*sizeof(double));
1608 if (kernel->values == (double *) NULL)
1609 return(DestroyKernelInfo(kernel));
1611 /* set all kernel values along axises to given scale */
1612 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1613 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1614 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1615 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1616 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1621 if (args->rho < 1.0)
1622 kernel->width = kernel->height = 5; /* default radius 2 */
1624 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1625 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1627 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1628 kernel->height*sizeof(double));
1629 if (kernel->values == (double *) NULL)
1630 return(DestroyKernelInfo(kernel));
1632 /* set all kernel values along axises to given scale */
1633 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1634 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1635 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1636 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1637 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1651 if (args->rho < args->sigma)
1653 kernel->width = ((size_t)args->sigma)*2+1;
1654 limit1 = (ssize_t)(args->rho*args->rho);
1655 limit2 = (ssize_t)(args->sigma*args->sigma);
1659 kernel->width = ((size_t)args->rho)*2+1;
1660 limit1 = (ssize_t)(args->sigma*args->sigma);
1661 limit2 = (ssize_t)(args->rho*args->rho);
1664 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1666 kernel->height = kernel->width;
1667 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1668 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1669 kernel->height*sizeof(double));
1670 if (kernel->values == (double *) NULL)
1671 return(DestroyKernelInfo(kernel));
1673 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1674 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1675 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1676 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1677 { ssize_t radius=u*u+v*v;
1678 if (limit1 < radius && radius <= limit2)
1679 kernel->positive_range += kernel->values[i] = (double) scale;
1681 kernel->values[i] = nan;
1683 kernel->minimum = kernel->maximum = (double) scale;
1684 if ( type == PeaksKernel ) {
1685 /* set the central point in the middle */
1686 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1687 kernel->positive_range = 1.0;
1688 kernel->maximum = 1.0;
1694 kernel=AcquireKernelInfo("ThinSE:482");
1695 if (kernel == (KernelInfo *) NULL)
1697 kernel->type = type;
1698 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1703 kernel=AcquireKernelInfo("ThinSE:87");
1704 if (kernel == (KernelInfo *) NULL)
1706 kernel->type = type;
1707 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1710 case DiagonalsKernel:
1712 switch ( (int) args->rho ) {
1717 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1718 if (kernel == (KernelInfo *) NULL)
1720 kernel->type = type;
1721 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1722 if (new_kernel == (KernelInfo *) NULL)
1723 return(DestroyKernelInfo(kernel));
1724 new_kernel->type = type;
1725 LastKernelInfo(kernel)->next = new_kernel;
1726 ExpandMirrorKernelInfo(kernel);
1730 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1733 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1736 if (kernel == (KernelInfo *) NULL)
1738 kernel->type = type;
1739 RotateKernelInfo(kernel, args->sigma);
1742 case LineEndsKernel:
1743 { /* Kernels for finding the end of thin lines */
1744 switch ( (int) args->rho ) {
1747 /* set of kernels to find all end of lines */
1748 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1750 /* kernel for 4-connected line ends - no rotation */
1751 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1754 /* kernel to add for 8-connected lines - no rotation */
1755 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1758 /* kernel to add for orthogonal line ends - does not find corners */
1759 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1762 /* traditional line end - fails on last T end */
1763 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1766 if (kernel == (KernelInfo *) NULL)
1768 kernel->type = type;
1769 RotateKernelInfo(kernel, args->sigma);
1772 case LineJunctionsKernel:
1773 { /* kernels for finding the junctions of multiple lines */
1774 switch ( (int) args->rho ) {
1777 /* set of kernels to find all line junctions */
1778 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1781 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1784 /* Diagonal T Junctions */
1785 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1788 /* Orthogonal T Junctions */
1789 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1792 /* Diagonal X Junctions */
1793 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1796 /* Orthogonal X Junctions - minimal diamond kernel */
1797 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1800 if (kernel == (KernelInfo *) NULL)
1802 kernel->type = type;
1803 RotateKernelInfo(kernel, args->sigma);
1807 { /* Ridges - Ridge finding kernels */
1810 switch ( (int) args->rho ) {
1813 kernel=ParseKernelArray("3x1:0,1,0");
1814 if (kernel == (KernelInfo *) NULL)
1816 kernel->type = type;
1817 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1820 kernel=ParseKernelArray("4x1:0,1,1,0");
1821 if (kernel == (KernelInfo *) NULL)
1823 kernel->type = type;
1824 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1826 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1827 /* Unfortunatally we can not yet rotate a non-square kernel */
1828 /* But then we can't flip a non-symetrical kernel either */
1829 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1830 if (new_kernel == (KernelInfo *) NULL)
1831 return(DestroyKernelInfo(kernel));
1832 new_kernel->type = type;
1833 LastKernelInfo(kernel)->next = new_kernel;
1834 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1835 if (new_kernel == (KernelInfo *) NULL)
1836 return(DestroyKernelInfo(kernel));
1837 new_kernel->type = type;
1838 LastKernelInfo(kernel)->next = new_kernel;
1839 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1840 if (new_kernel == (KernelInfo *) NULL)
1841 return(DestroyKernelInfo(kernel));
1842 new_kernel->type = type;
1843 LastKernelInfo(kernel)->next = new_kernel;
1844 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1845 if (new_kernel == (KernelInfo *) NULL)
1846 return(DestroyKernelInfo(kernel));
1847 new_kernel->type = type;
1848 LastKernelInfo(kernel)->next = new_kernel;
1849 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1850 if (new_kernel == (KernelInfo *) NULL)
1851 return(DestroyKernelInfo(kernel));
1852 new_kernel->type = type;
1853 LastKernelInfo(kernel)->next = new_kernel;
1854 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1855 if (new_kernel == (KernelInfo *) NULL)
1856 return(DestroyKernelInfo(kernel));
1857 new_kernel->type = type;
1858 LastKernelInfo(kernel)->next = new_kernel;
1859 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1860 if (new_kernel == (KernelInfo *) NULL)
1861 return(DestroyKernelInfo(kernel));
1862 new_kernel->type = type;
1863 LastKernelInfo(kernel)->next = new_kernel;
1864 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1865 if (new_kernel == (KernelInfo *) NULL)
1866 return(DestroyKernelInfo(kernel));
1867 new_kernel->type = type;
1868 LastKernelInfo(kernel)->next = new_kernel;
1873 case ConvexHullKernel:
1877 /* first set of 8 kernels */
1878 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1879 if (kernel == (KernelInfo *) NULL)
1881 kernel->type = type;
1882 ExpandRotateKernelInfo(kernel, 90.0);
1883 /* append the mirror versions too - no flip function yet */
1884 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1885 if (new_kernel == (KernelInfo *) NULL)
1886 return(DestroyKernelInfo(kernel));
1887 new_kernel->type = type;
1888 ExpandRotateKernelInfo(new_kernel, 90.0);
1889 LastKernelInfo(kernel)->next = new_kernel;
1892 case SkeletonKernel:
1894 switch ( (int) args->rho ) {
1897 /* Traditional Skeleton...
1898 ** A cyclically rotated single kernel
1900 kernel=AcquireKernelInfo("ThinSE:482");
1901 if (kernel == (KernelInfo *) NULL)
1903 kernel->type = type;
1904 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1907 /* HIPR Variation of the cyclic skeleton
1908 ** Corners of the traditional method made more forgiving,
1909 ** but the retain the same cyclic order.
1911 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1912 if (kernel == (KernelInfo *) NULL)
1914 if (kernel->next == (KernelInfo *) NULL)
1915 return(DestroyKernelInfo(kernel));
1916 kernel->type = type;
1917 kernel->next->type = type;
1918 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1921 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1922 ** "Connectivity-Preserving Morphological Image Thransformations"
1923 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1924 ** http://www.leptonica.com/papers/conn.pdf
1926 kernel=AcquireKernelInfo(
1927 "ThinSE:41; ThinSE:42; ThinSE:43");
1928 if (kernel == (KernelInfo *) NULL)
1930 kernel->type = type;
1931 kernel->next->type = type;
1932 kernel->next->next->type = type;
1933 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1939 { /* Special kernels for general thinning, while preserving connections
1940 ** "Connectivity-Preserving Morphological Image Thransformations"
1941 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1942 ** http://www.leptonica.com/papers/conn.pdf
1944 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
1946 ** Note kernels do not specify the origin pixel, allowing them
1947 ** to be used for both thickening and thinning operations.
1949 switch ( (int) args->rho ) {
1950 /* SE for 4-connected thinning */
1951 case 41: /* SE_4_1 */
1952 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
1954 case 42: /* SE_4_2 */
1955 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
1957 case 43: /* SE_4_3 */
1958 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
1960 case 44: /* SE_4_4 */
1961 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
1963 case 45: /* SE_4_5 */
1964 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
1966 case 46: /* SE_4_6 */
1967 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
1969 case 47: /* SE_4_7 */
1970 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
1972 case 48: /* SE_4_8 */
1973 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
1975 case 49: /* SE_4_9 */
1976 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
1978 /* SE for 8-connected thinning - negatives of the above */
1979 case 81: /* SE_8_0 */
1980 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
1982 case 82: /* SE_8_2 */
1983 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
1985 case 83: /* SE_8_3 */
1986 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
1988 case 84: /* SE_8_4 */
1989 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
1991 case 85: /* SE_8_5 */
1992 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
1994 case 86: /* SE_8_6 */
1995 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
1997 case 87: /* SE_8_7 */
1998 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2000 case 88: /* SE_8_8 */
2001 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2003 case 89: /* SE_8_9 */
2004 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2006 /* Special combined SE kernels */
2007 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2008 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2010 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2011 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2013 case 481: /* SE_48_1 - General Connected Corner Kernel */
2014 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2017 case 482: /* SE_48_2 - General Edge Kernel */
2018 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2021 if (kernel == (KernelInfo *) NULL)
2023 kernel->type = type;
2024 RotateKernelInfo(kernel, args->sigma);
2028 Distance Measuring Kernels
2030 case ChebyshevKernel:
2032 if (args->rho < 1.0)
2033 kernel->width = kernel->height = 3; /* default radius = 1 */
2035 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2036 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2038 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2039 kernel->height*sizeof(double));
2040 if (kernel->values == (double *) NULL)
2041 return(DestroyKernelInfo(kernel));
2043 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2044 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2045 kernel->positive_range += ( kernel->values[i] =
2046 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2047 kernel->maximum = kernel->values[0];
2050 case ManhattanKernel:
2052 if (args->rho < 1.0)
2053 kernel->width = kernel->height = 3; /* default radius = 1 */
2055 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2056 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2058 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2059 kernel->height*sizeof(double));
2060 if (kernel->values == (double *) NULL)
2061 return(DestroyKernelInfo(kernel));
2063 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2064 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2065 kernel->positive_range += ( kernel->values[i] =
2066 args->sigma*(labs((long) u)+labs((long) v)) );
2067 kernel->maximum = kernel->values[0];
2070 case OctagonalKernel:
2072 if (args->rho < 2.0)
2073 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2075 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2076 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2078 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2079 kernel->height*sizeof(double));
2080 if (kernel->values == (double *) NULL)
2081 return(DestroyKernelInfo(kernel));
2083 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2084 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2087 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2088 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2089 kernel->positive_range += kernel->values[i] =
2090 args->sigma*MagickMax(r1,r2);
2092 kernel->maximum = kernel->values[0];
2095 case EuclideanKernel:
2097 if (args->rho < 1.0)
2098 kernel->width = kernel->height = 3; /* default radius = 1 */
2100 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2101 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2103 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2104 kernel->height*sizeof(double));
2105 if (kernel->values == (double *) NULL)
2106 return(DestroyKernelInfo(kernel));
2108 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2109 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2110 kernel->positive_range += ( kernel->values[i] =
2111 args->sigma*sqrt((double)(u*u+v*v)) );
2112 kernel->maximum = kernel->values[0];
2117 /* No-Op Kernel - Basically just a single pixel on its own */
2118 kernel=ParseKernelArray("1:1");
2119 if (kernel == (KernelInfo *) NULL)
2121 kernel->type = UndefinedKernel;
2130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2134 % C l o n e K e r n e l I n f o %
2138 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2140 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2141 % can be modified without effecting the original. The cloned kernel should
2142 % be destroyed using DestoryKernelInfo() when no longer needed.
2144 % The format of the CloneKernelInfo method is:
2146 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2148 % A description of each parameter follows:
2150 % o kernel: the Morphology/Convolution kernel to be cloned
2153 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2161 assert(kernel != (KernelInfo *) NULL);
2162 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2163 if (new_kernel == (KernelInfo *) NULL)
2165 *new_kernel=(*kernel); /* copy values in structure */
2167 /* replace the values with a copy of the values */
2168 new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
2169 kernel->height*sizeof(double));
2170 if (new_kernel->values == (double *) NULL)
2171 return(DestroyKernelInfo(new_kernel));
2172 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2173 new_kernel->values[i]=kernel->values[i];
2175 /* Also clone the next kernel in the kernel list */
2176 if ( kernel->next != (KernelInfo *) NULL ) {
2177 new_kernel->next = CloneKernelInfo(kernel->next);
2178 if ( new_kernel->next == (KernelInfo *) NULL )
2179 return(DestroyKernelInfo(new_kernel));
2186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2190 % D e s t r o y K e r n e l I n f o %
2194 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2196 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2199 % The format of the DestroyKernelInfo method is:
2201 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2203 % A description of each parameter follows:
2205 % o kernel: the Morphology/Convolution kernel to be destroyed
2208 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2210 assert(kernel != (KernelInfo *) NULL);
2212 if ( kernel->next != (KernelInfo *) NULL )
2213 kernel->next = DestroyKernelInfo(kernel->next);
2215 kernel->values = (double *)RelinquishMagickMemory(kernel->values);
2216 kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
2221 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2225 + E x p a n d M i r r o r K e r n e l I n f o %
2229 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2231 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2232 % sequence of 90-degree rotated kernels but providing a reflected 180
2233 % rotatation, before the -/+ 90-degree rotations.
2235 % This special rotation order produces a better, more symetrical thinning of
2238 % The format of the ExpandMirrorKernelInfo method is:
2240 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2242 % A description of each parameter follows:
2244 % o kernel: the Morphology/Convolution kernel
2246 % This function is only internel to this module, as it is not finalized,
2247 % especially with regard to non-orthogonal angles, and rotation of larger
2252 static void FlopKernelInfo(KernelInfo *kernel)
2253 { /* Do a Flop by reversing each row. */
2261 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2262 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2263 t=k[x], k[x]=k[r], k[r]=t;
2265 kernel->x = kernel->width - kernel->x - 1;
2266 angle = fmod(angle+180.0, 360.0);
2270 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2278 clone = CloneKernelInfo(last);
2279 RotateKernelInfo(clone, 180); /* flip */
2280 LastKernelInfo(last)->next = clone;
2283 clone = CloneKernelInfo(last);
2284 RotateKernelInfo(clone, 90); /* transpose */
2285 LastKernelInfo(last)->next = clone;
2288 clone = CloneKernelInfo(last);
2289 RotateKernelInfo(clone, 180); /* flop */
2290 LastKernelInfo(last)->next = clone;
2296 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2300 + E x p a n d R o t a t e K e r n e l I n f o %
2304 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2306 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2307 % incrementally by the angle given, until the kernel repeats.
2309 % WARNING: 45 degree rotations only works for 3x3 kernels.
2310 % While 90 degree roatations only works for linear and square kernels
2312 % The format of the ExpandRotateKernelInfo method is:
2314 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2316 % A description of each parameter follows:
2318 % o kernel: the Morphology/Convolution kernel
2320 % o angle: angle to rotate in degrees
2322 % This function is only internel to this module, as it is not finalized,
2323 % especially with regard to non-orthogonal angles, and rotation of larger
2327 /* Internal Routine - Return true if two kernels are the same */
2328 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2329 const KernelInfo *kernel2)
2334 /* check size and origin location */
2335 if ( kernel1->width != kernel2->width
2336 || kernel1->height != kernel2->height
2337 || kernel1->x != kernel2->x
2338 || kernel1->y != kernel2->y )
2341 /* check actual kernel values */
2342 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2343 /* Test for Nan equivalence */
2344 if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2346 if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2348 /* Test actual values are equivalent */
2349 if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
2356 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2364 clone = CloneKernelInfo(last);
2365 RotateKernelInfo(clone, angle);
2366 if ( SameKernelInfo(kernel, clone) == MagickTrue )
2368 LastKernelInfo(last)->next = clone;
2371 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2376 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2380 + C a l c M e t a K e r n a l I n f o %
2384 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2386 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2387 % using the kernel values. This should only ne used if it is not possible to
2388 % calculate that meta-data in some easier way.
2390 % It is important that the meta-data is correct before ScaleKernelInfo() is
2391 % used to perform kernel normalization.
2393 % The format of the CalcKernelMetaData method is:
2395 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2397 % A description of each parameter follows:
2399 % o kernel: the Morphology/Convolution kernel to modify
2401 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2402 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2403 % however is not true for flat-shaped morphological kernels.
2405 % WARNING: Only the specific kernel pointed to is modified, not a list of
2408 % This is an internal function and not expected to be useful outside this
2409 % module. This could change however.
2411 static void CalcKernelMetaData(KernelInfo *kernel)
2416 kernel->minimum = kernel->maximum = 0.0;
2417 kernel->negative_range = kernel->positive_range = 0.0;
2418 for (i=0; i < (kernel->width*kernel->height); i++)
2420 if ( fabs(kernel->values[i]) < MagickEpsilon )
2421 kernel->values[i] = 0.0;
2422 ( kernel->values[i] < 0)
2423 ? ( kernel->negative_range += kernel->values[i] )
2424 : ( kernel->positive_range += kernel->values[i] );
2425 Minimize(kernel->minimum, kernel->values[i]);
2426 Maximize(kernel->maximum, kernel->values[i]);
2433 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2437 % M o r p h o l o g y A p p l y %
2441 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2443 % MorphologyApply() applies a morphological method, multiple times using
2444 % a list of multiple kernels.
2446 % It is basically equivalent to as MorphologyImageChannel() (see below) but
2447 % without any user controls. This allows internel programs to use this
2448 % function, to actually perform a specific task without possible interference
2449 % by any API user supplied settings.
2451 % It is MorphologyImageChannel() task to extract any such user controls, and
2452 % pass them to this function for processing.
2454 % More specifically kernels are not normalized/scaled/blended by the
2455 % 'convolve:scale' Image Artifact (setting), nor is the convolve bias
2456 % (-bias setting or image->bias) loooked at, but must be supplied from the
2457 % function arguments.
2459 % The format of the MorphologyApply method is:
2461 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2462 % const ChannelType channel, const ssize_t iterations,
2463 % const KernelInfo *kernel, const CompositeMethod compose,
2464 % const double bias, ExceptionInfo *exception)
2466 % A description of each parameter follows:
2468 % o image: the source image
2470 % o method: the morphology method to be applied.
2472 % o channel: the channels to which the operations are applied
2473 % The channel 'sync' flag determines if 'alpha weighting' is
2474 % applied for convolution style operations.
2476 % o iterations: apply the operation this many times (or no change).
2477 % A value of -1 means loop until no change found.
2478 % How this is applied may depend on the morphology method.
2479 % Typically this is a value of 1.
2481 % o channel: the channel type.
2483 % o kernel: An array of double representing the morphology kernel.
2485 % o compose: How to handle or merge multi-kernel results.
2486 % If 'UndefinedCompositeOp' use default for the Morphology method.
2487 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2488 % Otherwise merge the results using the compose method given.
2490 % o bias: Convolution Output Bias.
2492 % o exception: return any errors or warnings in this structure.
2496 /* Apply a Morphology Primative to an image using the given kernel.
2497 ** Two pre-created images must be provided, and no image is created.
2498 ** It returns the number of pixels that changed between the images
2499 ** for result convergence determination.
2501 static ssize_t MorphologyPrimitive(const Image *image, Image *morphology_image,
2502 const MorphologyMethod method, const ChannelType channel,
2503 const KernelInfo *kernel,const double bias,ExceptionInfo *exception)
2505 #define MorphologyTag "Morphology/Image"
2524 assert(image != (Image *) NULL);
2525 assert(image->signature == MagickSignature);
2526 assert(morphology_image != (Image *) NULL);
2527 assert(morphology_image->signature == MagickSignature);
2528 assert(kernel != (KernelInfo *) NULL);
2529 assert(kernel->signature == MagickSignature);
2530 assert(exception != (ExceptionInfo *) NULL);
2531 assert(exception->signature == MagickSignature);
2537 image_view=AcquireCacheView(image);
2538 morphology_view=AcquireCacheView(morphology_image);
2539 virt_width=image->columns+kernel->width-1;
2541 /* Some methods (including convolve) needs use a reflected kernel.
2542 * Adjust 'origin' offsets to loop though kernel as a reflection.
2547 case ConvolveMorphology:
2548 case DilateMorphology:
2549 case DilateIntensityMorphology:
2550 /*case DistanceMorphology:*/
2551 /* kernel needs to used with reflection about origin */
2552 offx = (ssize_t) kernel->width-offx-1;
2553 offy = (ssize_t) kernel->height-offy-1;
2555 case ErodeMorphology:
2556 case ErodeIntensityMorphology:
2557 case HitAndMissMorphology:
2558 case ThinningMorphology:
2559 case ThickenMorphology:
2560 /* kernel is used as is, without reflection */
2563 assert("Not a Primitive Morphology Method" != (char *) NULL);
2568 if ( method == ConvolveMorphology && kernel->width == 1 )
2569 { /* Special handling (for speed) of vertical (blur) kernels.
2570 ** This performs its handling in columns rather than in rows.
2571 ** This is only done for convolve as it is the only method that
2572 ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
2574 ** Timing tests (on single CPU laptop)
2575 ** Using a vertical 1-d Blue with normal row-by-row (below)
2576 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2578 ** Using this column method
2579 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2582 ** Anthony Thyssen, 14 June 2010
2587 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2588 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2590 for (x=0; x < (ssize_t) image->columns; x++)
2592 register const Quantum
2604 if (status == MagickFalse)
2606 p=GetCacheViewVirtualPixels(image_view, x, -offy,1,
2607 image->rows+kernel->height-1, exception);
2608 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,morphology_image->rows,exception);
2609 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2614 /* offset to origin in 'p'. while 'q' points to it directly */
2617 for (y=0; y < (ssize_t) image->rows; y++)
2622 register const double
2625 register const Quantum
2631 /* Copy input image to the output image for unused channels
2632 * This removes need for 'cloning' a new image every iteration
2634 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2635 GetPixelComponents(image)),q);
2636 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2637 GetPixelComponents(image)),q);
2638 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2639 GetPixelComponents(image)),q);
2640 if (image->colorspace == CMYKColorspace)
2641 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2642 GetPixelComponents(image)),q);
2644 /* Set the bias of the weighted average output */
2649 result.black = bias;
2652 /* Weighted Average of pixels using reflected kernel
2654 ** NOTE for correct working of this operation for asymetrical
2655 ** kernels, the kernel needs to be applied in its reflected form.
2656 ** That is its values needs to be reversed.
2658 k = &kernel->values[ kernel->height-1 ];
2660 if ( ((channel & SyncChannels) == 0 ) ||
2661 (image->matte == MagickFalse) )
2662 { /* No 'Sync' involved.
2663 ** Convolution is simple greyscale channel operation
2665 for (v=0; v < (ssize_t) kernel->height; v++) {
2666 if ( IsNan(*k) ) continue;
2667 result.red += (*k)*GetPixelRed(image,k_pixels);
2668 result.green += (*k)*GetPixelGreen(image,k_pixels);
2669 result.blue += (*k)*GetPixelBlue(image,k_pixels);
2670 if (image->colorspace == CMYKColorspace)
2671 result.black+=(*k)*GetPixelBlack(image,k_pixels);
2672 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2676 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
2677 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
2678 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
2679 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
2680 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
2681 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
2682 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
2683 (image->colorspace == CMYKColorspace))
2684 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
2685 if (((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0) &&
2686 (image->matte == MagickTrue))
2687 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2690 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2691 ** Weight the color channels with Alpha Channel so that
2692 ** transparent pixels are not part of the results.
2695 alpha, /* alpha weighting of colors : kernel*alpha */
2696 gamma; /* divisor, sum of color weighting values */
2699 for (v=0; v < (ssize_t) kernel->height; v++) {
2700 if ( IsNan(*k) ) continue;
2701 alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels));
2703 result.red += alpha*GetPixelRed(image,k_pixels);
2704 result.green += alpha*GetPixelGreen(image,k_pixels);
2705 result.blue += alpha*GetPixelBlue(image,k_pixels);
2706 if (image->colorspace == CMYKColorspace)
2707 result.black += alpha*GetPixelBlack(image,k_pixels);
2708 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2712 /* Sync'ed channels, all channels are modified */
2713 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2714 SetPixelRed(morphology_image,
2715 ClampToQuantum(gamma*result.red),q);
2716 SetPixelGreen(morphology_image,
2717 ClampToQuantum(gamma*result.green),q);
2718 SetPixelBlue(morphology_image,
2719 ClampToQuantum(gamma*result.blue),q);
2720 if (image->colorspace == CMYKColorspace)
2721 SetPixelBlack(morphology_image,
2722 ClampToQuantum(gamma*result.black),q);
2723 SetPixelAlpha(morphology_image,
2724 ClampToQuantum(result.alpha),q);
2727 /* Count up changed pixels */
2728 if ((GetPixelRed(image,p+r) != GetPixelRed(morphology_image,q))
2729 || (GetPixelGreen(image,p+r) != GetPixelGreen(morphology_image,q))
2730 || (GetPixelBlue(image,p+r) != GetPixelBlue(morphology_image,q))
2731 || (GetPixelAlpha(image,p+r) != GetPixelAlpha(morphology_image,q))
2732 || ((image->colorspace == CMYKColorspace) &&
2733 (GetPixelBlack(image,p+r) != GetPixelBlack(morphology_image,q))))
2734 changed++; /* The pixel was changed in some way! */
2735 p+=GetPixelComponents(image);
2736 q+=GetPixelComponents(morphology_image);
2738 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2740 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2745 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2746 #pragma omp critical (MagickCore_MorphologyImage)
2748 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2749 if (proceed == MagickFalse)
2753 morphology_image->type=image->type;
2754 morphology_view=DestroyCacheView(morphology_view);
2755 image_view=DestroyCacheView(image_view);
2756 return(status ? (ssize_t) changed : 0);
2760 ** Normal handling of horizontal or rectangular kernels (row by row)
2762 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2763 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2765 for (y=0; y < (ssize_t) image->rows; y++)
2767 register const Quantum
2779 if (status == MagickFalse)
2781 p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
2782 kernel->height, exception);
2783 q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,1,
2785 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2790 /* offset to origin in 'p'. while 'q' points to it directly */
2791 r = virt_width*offy + offx;
2793 for (x=0; x < (ssize_t) image->columns; x++)
2801 register const double
2804 register const Quantum
2812 /* Copy input image to the output image for unused channels
2813 * This removes need for 'cloning' a new image every iteration
2815 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2816 GetPixelComponents(image)),q);
2817 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2818 GetPixelComponents(image)),q);
2819 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2820 GetPixelComponents(image)),q);
2821 if (image->colorspace == CMYKColorspace)
2822 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2823 GetPixelComponents(image)),q);
2830 min.black = (MagickRealType) QuantumRange;
2835 max.black = (MagickRealType) 0;
2836 /* default result is the original pixel value */
2837 result.red = (MagickRealType) GetPixelRed(image,p+r);
2838 result.green = (MagickRealType) GetPixelGreen(image,p+r);
2839 result.blue = (MagickRealType) GetPixelBlue(image,p+r);
2841 if (image->colorspace == CMYKColorspace)
2842 result.black = (MagickRealType) GetPixelBlack(image,p+r);
2843 result.alpha=(MagickRealType) GetPixelAlpha(image,p+r);
2846 case ConvolveMorphology:
2847 /* Set the bias of the weighted average output */
2852 result.black = bias;
2854 case DilateIntensityMorphology:
2855 case ErodeIntensityMorphology:
2856 /* use a boolean flag indicating when first match found */
2857 result.red = 0.0; /* result is not used otherwise */
2864 case ConvolveMorphology:
2865 /* Weighted Average of pixels using reflected kernel
2867 ** NOTE for correct working of this operation for asymetrical
2868 ** kernels, the kernel needs to be applied in its reflected form.
2869 ** That is its values needs to be reversed.
2871 ** Correlation is actually the same as this but without reflecting
2872 ** the kernel, and thus 'lower-level' that Convolution. However
2873 ** as Convolution is the more common method used, and it does not
2874 ** really cost us much in terms of processing to use a reflected
2875 ** kernel, so it is Convolution that is implemented.
2877 ** Correlation will have its kernel reflected before calling
2878 ** this function to do a Convolve.
2880 ** For more details of Correlation vs Convolution see
2881 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2883 k = &kernel->values[ kernel->width*kernel->height-1 ];
2885 if ( ((channel & SyncChannels) == 0 ) ||
2886 (image->matte == MagickFalse) )
2887 { /* No 'Sync' involved.
2888 ** Convolution is simple greyscale channel operation
2890 for (v=0; v < (ssize_t) kernel->height; v++) {
2891 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2892 if ( IsNan(*k) ) continue;
2894 GetPixelRed(image,k_pixels+u*GetPixelComponents(image));
2895 result.green += (*k)*
2896 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image));
2897 result.blue += (*k)*
2898 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image));
2899 if (image->colorspace == CMYKColorspace)
2900 result.black += (*k)*
2901 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image));
2902 result.alpha += (*k)*
2903 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image));
2905 k_pixels += virt_width*GetPixelComponents(image);
2907 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
2908 SetPixelRed(morphology_image,ClampToQuantum(result.red),
2910 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
2911 SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2913 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
2914 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2916 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
2917 (image->colorspace == CMYKColorspace))
2918 SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2920 if (((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0) &&
2921 (image->matte == MagickTrue))
2922 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
2926 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2927 ** Weight the color channels with Alpha Channel so that
2928 ** transparent pixels are not part of the results.
2931 alpha, /* alpha weighting of colors : kernel*alpha */
2932 gamma; /* divisor, sum of color weighting values */
2935 for (v=0; v < (ssize_t) kernel->height; v++) {
2936 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2937 if ( IsNan(*k) ) continue;
2938 alpha=(*k)*(QuantumScale*
2939 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
2941 result.red += alpha*
2942 GetPixelRed(image,k_pixels+u*GetPixelComponents(image));
2943 result.green += alpha*
2944 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image));
2945 result.blue += alpha*
2946 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image));
2947 if (image->colorspace == CMYKColorspace)
2948 result.black+=alpha*
2949 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image));
2950 result.alpha += (*k)*
2951 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image));
2953 k_pixels += virt_width*GetPixelComponents(image);
2955 /* Sync'ed channels, all channels are modified */
2956 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2957 SetPixelRed(morphology_image,
2958 ClampToQuantum(gamma*result.red),q);
2959 SetPixelGreen(morphology_image,
2960 ClampToQuantum(gamma*result.green),q);
2961 SetPixelBlue(morphology_image,
2962 ClampToQuantum(gamma*result.blue),q);
2963 if (image->colorspace == CMYKColorspace)
2964 SetPixelBlack(morphology_image,
2965 ClampToQuantum(gamma*result.black),q);
2966 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2970 case ErodeMorphology:
2971 /* Minimum Value within kernel neighbourhood
2973 ** NOTE that the kernel is not reflected for this operation!
2975 ** NOTE: in normal Greyscale Morphology, the kernel value should
2976 ** be added to the real value, this is currently not done, due to
2977 ** the nature of the boolean kernels being used.
2981 for (v=0; v < (ssize_t) kernel->height; v++) {
2982 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2983 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2984 Minimize(min.red, (double)
2985 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
2986 Minimize(min.green, (double)
2987 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
2988 Minimize(min.blue, (double)
2989 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
2990 if (image->colorspace == CMYKColorspace)
2991 Minimize(min.black,(double)
2992 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
2993 Minimize(min.alpha,(double)
2994 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
2996 k_pixels += virt_width*GetPixelComponents(image);
3000 case DilateMorphology:
3001 /* Maximum Value within kernel neighbourhood
3003 ** NOTE for correct working of this operation for asymetrical
3004 ** kernels, the kernel needs to be applied in its reflected form.
3005 ** That is its values needs to be reversed.
3007 ** NOTE: in normal Greyscale Morphology, the kernel value should
3008 ** be added to the real value, this is currently not done, due to
3009 ** the nature of the boolean kernels being used.
3012 k = &kernel->values[ kernel->width*kernel->height-1 ];
3014 for (v=0; v < (ssize_t) kernel->height; v++) {
3015 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3016 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3017 Maximize(max.red, (double)
3018 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3019 Maximize(max.green, (double)
3020 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3021 Maximize(max.blue, (double)
3022 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3023 if (image->colorspace == CMYKColorspace)
3024 Maximize(max.black, (double)
3025 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3026 Maximize(max.alpha,(double)
3027 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3029 k_pixels += virt_width*GetPixelComponents(image);
3033 case HitAndMissMorphology:
3034 case ThinningMorphology:
3035 case ThickenMorphology:
3036 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3038 ** NOTE that the kernel is not reflected for this operation,
3039 ** and consists of both foreground and background pixel
3040 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3041 ** with either Nan or 0.5 values for don't care.
3043 ** Note that this will never produce a meaningless negative
3044 ** result. Such results can cause Thinning/Thicken to not work
3045 ** correctly when used against a greyscale image.
3049 for (v=0; v < (ssize_t) kernel->height; v++) {
3050 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3051 if ( IsNan(*k) ) continue;
3053 { /* minimim of foreground pixels */
3054 Minimize(min.red, (double)
3055 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3056 Minimize(min.green, (double)
3057 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3058 Minimize(min.blue, (double)
3059 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3060 if ( image->colorspace == CMYKColorspace)
3061 Minimize(min.black,(double)
3062 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3063 Minimize(min.alpha,(double)
3064 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3066 else if ( (*k) < 0.3 )
3067 { /* maximum of background pixels */
3068 Maximize(max.red, (double)
3069 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3070 Maximize(max.green, (double)
3071 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3072 Maximize(max.blue, (double)
3073 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3074 if (image->colorspace == CMYKColorspace)
3075 Maximize(max.black, (double)
3076 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3077 Maximize(max.alpha,(double)
3078 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3081 k_pixels += virt_width*GetPixelComponents(image);
3083 /* Pattern Match if difference is positive */
3084 min.red -= max.red; Maximize( min.red, 0.0 );
3085 min.green -= max.green; Maximize( min.green, 0.0 );
3086 min.blue -= max.blue; Maximize( min.blue, 0.0 );
3087 min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3088 min.black -= max.black; Maximize( min.black, 0.0 );
3091 case ErodeIntensityMorphology:
3092 /* Select Pixel with Minimum Intensity within kernel neighbourhood
3094 ** WARNING: the intensity test fails for CMYK and does not
3095 ** take into account the moderating effect of the alpha channel
3096 ** on the intensity.
3098 ** NOTE that the kernel is not reflected for this operation!
3102 for (v=0; v < (ssize_t) kernel->height; v++) {
3103 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3104 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3105 if ( result.red == 0.0 ||
3106 GetPixelIntensity(image,&(k_pixels[u])) < GetPixelIntensity(morphology_image,q) ) {
3107 /* copy the whole pixel - no channel selection */
3109 if ( result.red > 0.0 ) changed++;
3113 k_pixels += virt_width*GetPixelComponents(image);
3117 case DilateIntensityMorphology:
3118 /* Select Pixel with Maximum Intensity within kernel neighbourhood
3120 ** WARNING: the intensity test fails for CMYK and does not
3121 ** take into account the moderating effect of the alpha channel
3122 ** on the intensity (yet).
3124 ** NOTE for correct working of this operation for asymetrical
3125 ** kernels, the kernel needs to be applied in its reflected form.
3126 ** That is its values needs to be reversed.
3128 k = &kernel->values[ kernel->width*kernel->height-1 ];
3130 for (v=0; v < (ssize_t) kernel->height; v++) {
3131 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3132 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3133 if ( result.red == 0.0 ||
3134 GetPixelIntensity(image,&(k_pixels[u])) > GetPixelIntensity(morphology_image,q) ) {
3135 /* copy the whole pixel - no channel selection */
3137 if ( result.red > 0.0 ) changed++;
3141 k_pixels += virt_width*GetPixelComponents(image);
3145 This code has been obsoleted by the MorphologyPrimitiveDirect() function.
3146 However it is still (almost) correct coding for Grayscale Morphology.
3149 GrayErode is equivalent but with kernel values subtracted from pixels
3150 without the kernel rotation
3151 GreyDilate is equivalent but using Maximum() instead of Minimum()
3152 useing kernel rotation
3154 case DistanceMorphology:
3155 /* Add kernel Value and select the minimum value found.
3156 ** The result is a iterative distance from edge of image shape.
3158 ** All Distance Kernels are symetrical, but that may not always
3159 ** be the case. For example how about a distance from left edges?
3160 ** To work correctly with asymetrical kernels the reflected kernel
3161 ** needs to be applied.
3163 k = &kernel->values[ kernel->width*kernel->height-1 ];
3165 for (v=0; v < (ssize_t) kernel->height; v++) {
3166 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3167 if ( IsNan(*k) ) continue;
3168 Minimize(result.red, (*k)+k_pixels[u].red);
3169 Minimize(result.green, (*k)+k_pixels[u].green);
3170 Minimize(result.blue, (*k)+k_pixels[u].blue);
3171 Minimize(result.alpha, (*k)+k_pixels[u].alpha);
3172 if ( image->colorspace == CMYKColorspace)
3173 Minimize(result.black,(*k)+GetPixelBlack(p_image,k_indexes+u));
3175 k_pixels += virt_width*GetPixelComponents(image);
3179 case UndefinedMorphology:
3181 break; /* Do nothing */
3183 /* Final mathematics of results (combine with original image?)
3185 ** NOTE: Difference Morphology operators Edge* and *Hat could also
3186 ** be done here but works better with iteration as a image difference
3187 ** in the controling function (below). Thicken and Thinning however
3188 ** should be done here so thay can be iterated correctly.
3191 case HitAndMissMorphology:
3192 case ErodeMorphology:
3193 result = min; /* minimum of neighbourhood */
3195 case DilateMorphology:
3196 result = max; /* maximum of neighbourhood */
3198 case ThinningMorphology:
3199 /* subtract pattern match from original */
3200 result.red -= min.red;
3201 result.green -= min.green;
3202 result.blue -= min.blue;
3203 result.alpha -= min.alpha;
3204 result.black -= min.black;
3206 case ThickenMorphology:
3207 /* Add the pattern matchs to the original */
3208 result.red += min.red;
3209 result.green += min.green;
3210 result.blue += min.blue;
3211 result.alpha += min.alpha;
3212 result.black += min.black;
3215 /* result directly calculated or assigned */
3218 /* Assign the resulting pixel values - Clamping Result */
3220 case UndefinedMorphology:
3221 case ConvolveMorphology:
3222 case DilateIntensityMorphology:
3223 case ErodeIntensityMorphology:
3224 break; /* full pixel was directly assigned - not a channel method */
3226 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
3227 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
3228 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
3229 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
3230 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
3231 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
3232 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
3233 (image->colorspace == CMYKColorspace))
3234 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
3235 if (((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0) &&
3236 (image->matte == MagickTrue))
3237 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3240 /* Count up changed pixels */
3241 if ((GetPixelRed(image,p+r) != GetPixelRed(morphology_image,q)) ||
3242 (GetPixelGreen(image,p+r) != GetPixelGreen(morphology_image,q)) ||
3243 (GetPixelBlue(image,p+r) != GetPixelBlue(morphology_image,q)) ||
3244 (GetPixelAlpha(image,p+r) != GetPixelAlpha(morphology_image,q)) ||
3245 ((image->colorspace == CMYKColorspace) &&
3246 (GetPixelBlack(image,p+r) != GetPixelBlack(morphology_image,q))))
3247 changed++; /* The pixel was changed in some way! */
3248 p+=GetPixelComponents(image);
3249 q+=GetPixelComponents(morphology_image);
3251 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3253 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3258 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3259 #pragma omp critical (MagickCore_MorphologyImage)
3261 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3262 if (proceed == MagickFalse)
3266 morphology_view=DestroyCacheView(morphology_view);
3267 image_view=DestroyCacheView(image_view);
3268 return(status ? (ssize_t)changed : -1);
3271 /* This is almost identical to the MorphologyPrimative() function above,
3272 ** but will apply the primitive directly to the image in two passes.
3274 ** That is after each row is 'Sync'ed' into the image, the next row will
3275 ** make use of those values as part of the calculation of the next row.
3276 ** It then repeats, but going in the oppisite (bottom-up) direction.
3278 ** Because of this 'iterative' handling this function can not make use
3279 ** of multi-threaded, parellel processing.
3281 static ssize_t MorphologyPrimitiveDirect(Image *image,
3282 const MorphologyMethod method, const ChannelType channel,
3283 const KernelInfo *kernel,ExceptionInfo *exception)
3306 assert(image != (Image *) NULL);
3307 assert(image->signature == MagickSignature);
3308 assert(kernel != (KernelInfo *) NULL);
3309 assert(kernel->signature == MagickSignature);
3310 assert(exception != (ExceptionInfo *) NULL);
3311 assert(exception->signature == MagickSignature);
3313 /* Some methods (including convolve) needs use a reflected kernel.
3314 * Adjust 'origin' offsets to loop though kernel as a reflection.
3319 case DistanceMorphology:
3320 case VoronoiMorphology:
3321 /* kernel needs to used with reflection about origin */
3322 offx = (ssize_t) kernel->width-offx-1;
3323 offy = (ssize_t) kernel->height-offy-1;
3326 case ?????Morphology:
3327 /* kernel is used as is, without reflection */
3331 assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
3335 /* DO NOT THREAD THIS CODE! */
3336 /* two views into same image (virtual, and actual) */
3337 virt_view=AcquireCacheView(image);
3338 auth_view=AcquireCacheView(image);
3339 virt_width=image->columns+kernel->width-1;
3341 for (y=0; y < (ssize_t) image->rows; y++)
3343 register const Quantum
3355 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3356 ** we read using virtual to get virtual pixel handling, but write back
3357 ** into the same image.
3359 ** Only top half of kernel is processed as we do a single pass downward
3360 ** through the image iterating the distance function as we go.
3362 if (status == MagickFalse)
3364 p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
3366 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3368 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3370 if (status == MagickFalse)
3373 /* offset to origin in 'p'. while 'q' points to it directly */
3374 r = (ssize_t) virt_width*offy + offx;
3376 for (x=0; x < (ssize_t) image->columns; x++)
3384 register const double
3387 register const Quantum
3393 /* Starting Defaults */
3394 GetPixelInfo(image,&result);
3395 SetPixelInfo(image,q,&result);
3396 if ( method != VoronoiMorphology )
3397 result.alpha = QuantumRange - result.alpha;
3400 case DistanceMorphology:
3401 /* Add kernel Value and select the minimum value found. */
3402 k = &kernel->values[ kernel->width*kernel->height-1 ];
3404 for (v=0; v <= (ssize_t) offy; v++) {
3405 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3406 if ( IsNan(*k) ) continue;
3407 Minimize(result.red, (*k)+
3408 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3409 Minimize(result.green, (*k)+
3410 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3411 Minimize(result.blue, (*k)+
3412 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3413 if (image->colorspace == CMYKColorspace)
3414 Minimize(result.black,(*k)+
3415 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3416 Minimize(result.alpha, (*k)+
3417 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3419 k_pixels += virt_width*GetPixelComponents(image);
3421 /* repeat with the just processed pixels of this row */
3422 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3424 for (u=0; u < (ssize_t) offx; u++, k--) {
3425 if ( x+u-offx < 0 ) continue; /* off the edge! */
3426 if ( IsNan(*k) ) continue;
3427 Minimize(result.red, (*k)+
3428 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3429 Minimize(result.green, (*k)+
3430 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3431 Minimize(result.blue, (*k)+
3432 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3433 if (image->colorspace == CMYKColorspace)
3434 Minimize(result.black,(*k)+
3435 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3436 Minimize(result.alpha,(*k)+
3437 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3440 case VoronoiMorphology:
3441 /* Apply Distance to 'Matte' channel, coping the closest color.
3443 ** This is experimental, and realy the 'alpha' component should
3444 ** be completely separate 'masking' channel.
3446 k = &kernel->values[ kernel->width*kernel->height-1 ];
3448 for (v=0; v <= (ssize_t) offy; v++) {
3449 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3450 if ( IsNan(*k) ) continue;
3451 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)) )
3453 SetPixelInfo(image,k_pixels+u*GetPixelComponents(image),
3458 k_pixels += virt_width*GetPixelComponents(image);
3460 /* repeat with the just processed pixels of this row */
3461 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3463 for (u=0; u < (ssize_t) offx; u++, k--) {
3464 if ( x+u-offx < 0 ) continue; /* off the edge! */
3465 if ( IsNan(*k) ) continue;
3466 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)) )
3468 SetPixelInfo(image,k_pixels+u*GetPixelComponents(image),
3475 /* result directly calculated or assigned */
3478 /* Assign the resulting pixel values - Clamping Result */
3480 case VoronoiMorphology:
3481 SetPixelPixelInfo(image,&result,q);
3484 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
3485 SetPixelRed(image,ClampToQuantum(result.red),q);
3486 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
3487 SetPixelGreen(image,ClampToQuantum(result.green),q);
3488 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
3489 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3490 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
3491 (image->colorspace == CMYKColorspace))
3492 SetPixelBlack(image,ClampToQuantum(result.black),q);
3493 if ((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0 &&
3494 (image->matte == MagickTrue))
3495 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3498 /* Count up changed pixels */
3499 if ((GetPixelRed(image,p+r) != GetPixelRed(image,q)) ||
3500 (GetPixelGreen(image,p+r) != GetPixelGreen(image,q)) ||
3501 (GetPixelBlue(image,p+r) != GetPixelBlue(image,q)) ||
3502 (GetPixelAlpha(image,p+r) != GetPixelAlpha(image,q)) ||
3503 ((image->colorspace == CMYKColorspace) &&
3504 (GetPixelBlack(image,p+r) != GetPixelBlack(image,q))))
3505 changed++; /* The pixel was changed in some way! */
3507 p+=GetPixelComponents(image); /* increment pixel buffers */
3508 q+=GetPixelComponents(image);
3511 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3513 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3514 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3520 /* Do the reversed pass through the image */
3521 for (y=(ssize_t)image->rows-1; y >= 0; y--)
3523 register const Quantum
3535 if (status == MagickFalse)
3537 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3538 ** we read using virtual to get virtual pixel handling, but write back
3539 ** into the same image.
3541 ** Only the bottom half of the kernel will be processes as we
3544 p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
3545 kernel->y+1,exception);
3546 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3548 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3550 if (status == MagickFalse)
3553 /* adjust positions to end of row */
3554 p += image->columns-1;
3555 q += image->columns-1;
3557 /* offset to origin in 'p'. while 'q' points to it directly */
3560 for (x=(ssize_t)image->columns-1; x >= 0; x--)
3568 register const double
3571 register const Quantum
3577 /* Default - previously modified pixel */
3578 GetPixelInfo(image,&result);
3579 SetPixelInfo(image,q,&result);
3580 if ( method != VoronoiMorphology )
3581 result.alpha = QuantumRange - result.alpha;
3584 case DistanceMorphology:
3585 /* Add kernel Value and select the minimum value found. */
3586 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3588 for (v=offy; v < (ssize_t) kernel->height; v++) {
3589 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3590 if ( IsNan(*k) ) continue;
3591 Minimize(result.red, (*k)+
3592 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3593 Minimize(result.green, (*k)+
3594 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3595 Minimize(result.blue, (*k)+
3596 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3597 if ( image->colorspace == CMYKColorspace)
3598 Minimize(result.black,(*k)+
3599 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3600 Minimize(result.alpha, (*k)+
3601 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3603 k_pixels += virt_width*GetPixelComponents(image);
3605 /* repeat with the just processed pixels of this row */
3606 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3608 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3609 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3610 if ( IsNan(*k) ) continue;
3611 Minimize(result.red, (*k)+
3612 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3613 Minimize(result.green, (*k)+
3614 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3615 Minimize(result.blue, (*k)+
3616 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3617 if ( image->colorspace == CMYKColorspace)
3618 Minimize(result.black, (*k)+
3619 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3620 Minimize(result.alpha, (*k)+
3621 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3624 case VoronoiMorphology:
3625 /* Apply Distance to 'Matte' channel, coping the closest color.
3627 ** This is experimental, and realy the 'alpha' component should
3628 ** be completely separate 'masking' channel.
3630 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3632 for (v=offy; v < (ssize_t) kernel->height; v++) {
3633 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3634 if ( IsNan(*k) ) continue;
3635 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)) )
3637 SetPixelInfo(image,k_pixels+u*GetPixelComponents(image),
3642 k_pixels += virt_width*GetPixelComponents(image);
3644 /* repeat with the just processed pixels of this row */
3645 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3647 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3648 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3649 if ( IsNan(*k) ) continue;
3650 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)) )
3652 SetPixelInfo(image,k_pixels+u*GetPixelComponents(image),
3659 /* result directly calculated or assigned */
3662 /* Assign the resulting pixel values - Clamping Result */
3664 case VoronoiMorphology:
3665 SetPixelPixelInfo(image,&result,q);
3668 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
3669 SetPixelRed(image,ClampToQuantum(result.red),q);
3670 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
3671 SetPixelGreen(image,ClampToQuantum(result.green),q);
3672 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
3673 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3674 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
3675 (image->colorspace == CMYKColorspace))
3676 SetPixelBlack(image,ClampToQuantum(result.black),q);
3677 if ((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0 &&
3678 (image->matte == MagickTrue))
3679 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3682 /* Count up changed pixels */
3683 if ( (GetPixelRed(image,p+r) != GetPixelRed(image,q))
3684 || (GetPixelGreen(image,p+r) != GetPixelGreen(image,q))
3685 || (GetPixelBlue(image,p+r) != GetPixelBlue(image,q))
3686 || (GetPixelAlpha(image,p+r) != GetPixelAlpha(image,q))
3687 || ((image->colorspace == CMYKColorspace) &&
3688 (GetPixelBlack(image,p+r) != GetPixelBlack(image,q))))
3689 changed++; /* The pixel was changed in some way! */
3691 p--; /* go backward through pixel buffers */
3694 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3696 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3697 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3703 auth_view=DestroyCacheView(auth_view);
3704 virt_view=DestroyCacheView(virt_view);
3705 return(status ? (ssize_t) changed : -1);
3708 /* Apply a Morphology by calling theabove low level primitive application
3709 ** functions. This function handles any iteration loops, composition or
3710 ** re-iteration of results, and compound morphology methods that is based
3711 ** on multiple low-level (staged) morphology methods.
3713 ** Basically this provides the complex grue between the requested morphology
3714 ** method and raw low-level implementation (above).
3716 MagickExport Image *MorphologyApply(const Image *image, const ChannelType
3717 channel,const MorphologyMethod method, const ssize_t iterations,
3718 const KernelInfo *kernel, const CompositeOperator compose,
3719 const double bias, ExceptionInfo *exception)
3725 *curr_image, /* Image we are working with or iterating */
3726 *work_image, /* secondary image for primitive iteration */
3727 *save_image, /* saved image - for 'edge' method only */
3728 *rslt_image; /* resultant image - after multi-kernel handling */
3731 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3732 *norm_kernel, /* the current normal un-reflected kernel */
3733 *rflt_kernel, /* the current reflected kernel (if needed) */
3734 *this_kernel; /* the kernel being applied */
3737 primitive; /* the current morphology primitive being applied */
3740 rslt_compose; /* multi-kernel compose method for results to use */
3743 special, /* do we use a direct modify function? */
3744 verbose; /* verbose output of results */
3747 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3748 method_limit, /* maximum number of compound method iterations */
3749 kernel_number, /* Loop 2: the kernel number being applied */
3750 stage_loop, /* Loop 3: primitive loop for compound morphology */
3751 stage_limit, /* how many primitives are in this compound */
3752 kernel_loop, /* Loop 4: iterate the kernel over image */
3753 kernel_limit, /* number of times to iterate kernel */
3754 count, /* total count of primitive steps applied */
3755 kernel_changed, /* total count of changed using iterated kernel */
3756 method_changed; /* total count of changed over method iteration */
3759 changed; /* number pixels changed by last primitive operation */
3764 assert(image != (Image *) NULL);
3765 assert(image->signature == MagickSignature);
3766 assert(kernel != (KernelInfo *) NULL);
3767 assert(kernel->signature == MagickSignature);
3768 assert(exception != (ExceptionInfo *) NULL);
3769 assert(exception->signature == MagickSignature);
3771 count = 0; /* number of low-level morphology primitives performed */
3772 if ( iterations == 0 )
3773 return((Image *)NULL); /* null operation - nothing to do! */
3775 kernel_limit = (size_t) iterations;
3776 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3777 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3779 verbose = IsMagickTrue(GetImageArtifact(image,"verbose"));
3781 /* initialise for cleanup */
3782 curr_image = (Image *) image;
3783 curr_compose = image->compose;
3784 (void) curr_compose;
3785 work_image = save_image = rslt_image = (Image *) NULL;
3786 reflected_kernel = (KernelInfo *) NULL;
3788 /* Initialize specific methods
3789 * + which loop should use the given iteratations
3790 * + how many primitives make up the compound morphology
3791 * + multi-kernel compose method to use (by default)
3793 method_limit = 1; /* just do method once, unless otherwise set */
3794 stage_limit = 1; /* assume method is not a compound */
3795 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3796 rslt_compose = compose; /* and we are composing multi-kernels as given */
3798 case SmoothMorphology: /* 4 primitive compound morphology */
3801 case OpenMorphology: /* 2 primitive compound morphology */
3802 case OpenIntensityMorphology:
3803 case TopHatMorphology:
3804 case CloseMorphology:
3805 case CloseIntensityMorphology:
3806 case BottomHatMorphology:
3807 case EdgeMorphology:
3810 case HitAndMissMorphology:
3811 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3813 case ThinningMorphology:
3814 case ThickenMorphology:
3815 method_limit = kernel_limit; /* iterate the whole method */
3816 kernel_limit = 1; /* do not do kernel iteration */
3818 case DistanceMorphology:
3819 case VoronoiMorphology:
3820 special = MagickTrue;
3826 /* Apply special methods with special requirments
3827 ** For example, single run only, or post-processing requirements
3829 if ( special == MagickTrue )
3831 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3832 if (rslt_image == (Image *) NULL)
3834 if (SetImageStorageClass(rslt_image,DirectClass) == MagickFalse)
3836 InheritException(exception,&rslt_image->exception);
3840 changed = MorphologyPrimitiveDirect(rslt_image, method,
3841 channel, kernel, exception);
3843 if ( verbose == MagickTrue )
3844 (void) (void) FormatLocaleFile(stderr,
3845 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3846 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3847 1.0,0.0,1.0, (double) changed);
3852 if ( method == VoronoiMorphology ) {
3853 /* Preserve the alpha channel of input image - but turned off */
3854 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel);
3855 (void) CompositeImageChannel(rslt_image, DefaultChannels,
3856 CopyOpacityCompositeOp, image, 0, 0);
3857 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel);
3862 /* Handle user (caller) specified multi-kernel composition method */
3863 if ( compose != UndefinedCompositeOp )
3864 rslt_compose = compose; /* override default composition for method */
3865 if ( rslt_compose == UndefinedCompositeOp )
3866 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3868 /* Some methods require a reflected kernel to use with primitives.
3869 * Create the reflected kernel for those methods. */
3871 case CorrelateMorphology:
3872 case CloseMorphology:
3873 case CloseIntensityMorphology:
3874 case BottomHatMorphology:
3875 case SmoothMorphology:
3876 reflected_kernel = CloneKernelInfo(kernel);
3877 if (reflected_kernel == (KernelInfo *) NULL)
3879 RotateKernelInfo(reflected_kernel,180);
3885 /* Loops around more primitive morpholgy methods
3886 ** erose, dilate, open, close, smooth, edge, etc...
3888 /* Loop 1: iterate the compound method */
3891 while ( method_loop < method_limit && method_changed > 0 ) {
3895 /* Loop 2: iterate over each kernel in a multi-kernel list */
3896 norm_kernel = (KernelInfo *) kernel;
3897 this_kernel = (KernelInfo *) kernel;
3898 rflt_kernel = reflected_kernel;
3901 while ( norm_kernel != NULL ) {
3903 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3904 stage_loop = 0; /* the compound morphology stage number */
3905 while ( stage_loop < stage_limit ) {
3906 stage_loop++; /* The stage of the compound morphology */
3908 /* Select primitive morphology for this stage of compound method */
3909 this_kernel = norm_kernel; /* default use unreflected kernel */
3910 primitive = method; /* Assume method is a primitive */
3912 case ErodeMorphology: /* just erode */
3913 case EdgeInMorphology: /* erode and image difference */
3914 primitive = ErodeMorphology;
3916 case DilateMorphology: /* just dilate */
3917 case EdgeOutMorphology: /* dilate and image difference */
3918 primitive = DilateMorphology;
3920 case OpenMorphology: /* erode then dialate */
3921 case TopHatMorphology: /* open and image difference */
3922 primitive = ErodeMorphology;
3923 if ( stage_loop == 2 )
3924 primitive = DilateMorphology;
3926 case OpenIntensityMorphology:
3927 primitive = ErodeIntensityMorphology;
3928 if ( stage_loop == 2 )
3929 primitive = DilateIntensityMorphology;
3931 case CloseMorphology: /* dilate, then erode */
3932 case BottomHatMorphology: /* close and image difference */
3933 this_kernel = rflt_kernel; /* use the reflected kernel */
3934 primitive = DilateMorphology;
3935 if ( stage_loop == 2 )
3936 primitive = ErodeMorphology;
3938 case CloseIntensityMorphology:
3939 this_kernel = rflt_kernel; /* use the reflected kernel */
3940 primitive = DilateIntensityMorphology;
3941 if ( stage_loop == 2 )
3942 primitive = ErodeIntensityMorphology;
3944 case SmoothMorphology: /* open, close */
3945 switch ( stage_loop ) {
3946 case 1: /* start an open method, which starts with Erode */
3947 primitive = ErodeMorphology;
3949 case 2: /* now Dilate the Erode */
3950 primitive = DilateMorphology;
3952 case 3: /* Reflect kernel a close */
3953 this_kernel = rflt_kernel; /* use the reflected kernel */
3954 primitive = DilateMorphology;
3956 case 4: /* Finish the Close */
3957 this_kernel = rflt_kernel; /* use the reflected kernel */
3958 primitive = ErodeMorphology;
3962 case EdgeMorphology: /* dilate and erode difference */
3963 primitive = DilateMorphology;
3964 if ( stage_loop == 2 ) {
3965 save_image = curr_image; /* save the image difference */
3966 curr_image = (Image *) image;
3967 primitive = ErodeMorphology;
3970 case CorrelateMorphology:
3971 /* A Correlation is a Convolution with a reflected kernel.
3972 ** However a Convolution is a weighted sum using a reflected
3973 ** kernel. It may seem stange to convert a Correlation into a
3974 ** Convolution as the Correlation is the simplier method, but
3975 ** Convolution is much more commonly used, and it makes sense to
3976 ** implement it directly so as to avoid the need to duplicate the
3977 ** kernel when it is not required (which is typically the
3980 this_kernel = rflt_kernel; /* use the reflected kernel */
3981 primitive = ConvolveMorphology;
3986 assert( this_kernel != (KernelInfo *) NULL );
3988 /* Extra information for debugging compound operations */
3989 if ( verbose == MagickTrue ) {
3990 if ( stage_limit > 1 )
3991 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
3992 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
3993 method_loop,(double) stage_loop);
3994 else if ( primitive != method )
3995 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
3996 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
4002 /* Loop 4: Iterate the kernel with primitive */
4006 while ( kernel_loop < kernel_limit && changed > 0 ) {
4007 kernel_loop++; /* the iteration of this kernel */
4009 /* Create a clone as the destination image, if not yet defined */
4010 if ( work_image == (Image *) NULL )
4012 work_image=CloneImage(image,0,0,MagickTrue,exception);
4013 if (work_image == (Image *) NULL)
4015 if (SetImageStorageClass(work_image,DirectClass) == MagickFalse)
4017 InheritException(exception,&work_image->exception);
4020 /* work_image->type=image->type; ??? */
4023 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4025 changed = MorphologyPrimitive(curr_image, work_image, primitive,
4026 channel, this_kernel, bias, exception);
4028 if ( verbose == MagickTrue ) {
4029 if ( kernel_loop > 1 )
4030 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4031 (void) (void) FormatLocaleFile(stderr,
4032 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4033 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4034 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4035 (double) (method_loop+kernel_loop-1),(double) kernel_number,
4036 (double) count,(double) changed);
4040 kernel_changed += changed;
4041 method_changed += changed;
4043 /* prepare next loop */
4044 { Image *tmp = work_image; /* swap images for iteration */
4045 work_image = curr_image;
4048 if ( work_image == image )
4049 work_image = (Image *) NULL; /* replace input 'image' */
4051 } /* End Loop 4: Iterate the kernel with primitive */
4053 if ( verbose == MagickTrue && kernel_changed != (size_t)changed )
4054 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
4055 if ( verbose == MagickTrue && stage_loop < stage_limit )
4056 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4059 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4060 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
4061 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
4062 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
4063 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
4066 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4068 /* Final Post-processing for some Compound Methods
4070 ** The removal of any 'Sync' channel flag in the Image Compositon
4071 ** below ensures the methematical compose method is applied in a
4072 ** purely mathematical way, and only to the selected channels.
4073 ** Turn off SVG composition 'alpha blending'.
4076 case EdgeOutMorphology:
4077 case EdgeInMorphology:
4078 case TopHatMorphology:
4079 case BottomHatMorphology:
4080 if ( verbose == MagickTrue )
4081 (void) FormatLocaleFile(stderr, "\n%s: Difference with original image",
4082 CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4083 (void) CompositeImageChannel(curr_image,
4084 (ChannelType) (channel & ~SyncChannels),
4085 DifferenceCompositeOp, image, 0, 0);
4087 case EdgeMorphology:
4088 if ( verbose == MagickTrue )
4089 (void) FormatLocaleFile(stderr, "\n%s: Difference of Dilate and Erode",
4090 CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4091 (void) CompositeImageChannel(curr_image,
4092 (ChannelType) (channel & ~SyncChannels),
4093 DifferenceCompositeOp, save_image, 0, 0);
4094 save_image = DestroyImage(save_image); /* finished with save image */
4100 /* multi-kernel handling: re-iterate, or compose results */
4101 if ( kernel->next == (KernelInfo *) NULL )
4102 rslt_image = curr_image; /* just return the resulting image */
4103 else if ( rslt_compose == NoCompositeOp )
4104 { if ( verbose == MagickTrue ) {
4105 if ( this_kernel->next != (KernelInfo *) NULL )
4106 (void) FormatLocaleFile(stderr, " (re-iterate)");
4108 (void) FormatLocaleFile(stderr, " (done)");
4110 rslt_image = curr_image; /* return result, and re-iterate */
4112 else if ( rslt_image == (Image *) NULL)
4113 { if ( verbose == MagickTrue )
4114 (void) FormatLocaleFile(stderr, " (save for compose)");
4115 rslt_image = curr_image;
4116 curr_image = (Image *) image; /* continue with original image */
4119 { /* Add the new 'current' result to the composition
4121 ** The removal of any 'Sync' channel flag in the Image Compositon
4122 ** below ensures the methematical compose method is applied in a
4123 ** purely mathematical way, and only to the selected channels.
4124 ** IE: Turn off SVG composition 'alpha blending'.
4126 if ( verbose == MagickTrue )
4127 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4128 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4129 (void) CompositeImageChannel(rslt_image,
4130 (ChannelType) (channel & ~SyncChannels), rslt_compose,
4132 curr_image = DestroyImage(curr_image);
4133 curr_image = (Image *) image; /* continue with original image */
4135 if ( verbose == MagickTrue )
4136 (void) FormatLocaleFile(stderr, "\n");
4138 /* loop to the next kernel in a multi-kernel list */
4139 norm_kernel = norm_kernel->next;
4140 if ( rflt_kernel != (KernelInfo *) NULL )
4141 rflt_kernel = rflt_kernel->next;
4143 } /* End Loop 2: Loop over each kernel */
4145 } /* End Loop 1: compound method interation */
4149 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4151 if ( curr_image == rslt_image )
4152 curr_image = (Image *) NULL;
4153 if ( rslt_image != (Image *) NULL )
4154 rslt_image = DestroyImage(rslt_image);
4156 if ( curr_image == rslt_image || curr_image == image )
4157 curr_image = (Image *) NULL;
4158 if ( curr_image != (Image *) NULL )
4159 curr_image = DestroyImage(curr_image);
4160 if ( work_image != (Image *) NULL )
4161 work_image = DestroyImage(work_image);
4162 if ( save_image != (Image *) NULL )
4163 save_image = DestroyImage(save_image);
4164 if ( reflected_kernel != (KernelInfo *) NULL )
4165 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4175 % M o r p h o l o g y I m a g e C h a n n e l %
4179 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4181 % MorphologyImageChannel() applies a user supplied kernel to the image
4182 % according to the given mophology method.
4184 % This function applies any and all user defined settings before calling
4185 % the above internal function MorphologyApply().
4187 % User defined settings include...
4188 % * Output Bias for Convolution and correlation ("-bias")
4189 % * Kernel Scale/normalize settings ("-set 'option:convolve:scale'")
4190 % This can also includes the addition of a scaled unity kernel.
4191 % * Show Kernel being applied ("-set option:showkernel 1")
4193 % The format of the MorphologyImage method is:
4195 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4196 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4198 % Image *MorphologyImageChannel(const Image *image, const ChannelType
4199 % channel,MorphologyMethod method,const ssize_t iterations,
4200 % KernelInfo *kernel,ExceptionInfo *exception)
4202 % A description of each parameter follows:
4204 % o image: the image.
4206 % o method: the morphology method to be applied.
4208 % o iterations: apply the operation this many times (or no change).
4209 % A value of -1 means loop until no change found.
4210 % How this is applied may depend on the morphology method.
4211 % Typically this is a value of 1.
4213 % o channel: the channel type.
4215 % o kernel: An array of double representing the morphology kernel.
4216 % Warning: kernel may be normalized for the Convolve method.
4218 % o exception: return any errors or warnings in this structure.
4222 MagickExport Image *MorphologyImageChannel(const Image *image,
4223 const ChannelType channel,const MorphologyMethod method,
4224 const ssize_t iterations,const KernelInfo *kernel,ExceptionInfo *exception)
4236 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4237 * This is done BEFORE the ShowKernelInfo() function is called so that
4238 * users can see the results of the 'option:convolve:scale' option.
4240 curr_kernel = (KernelInfo *) kernel;
4241 if ( method == ConvolveMorphology || method == CorrelateMorphology )
4245 artifact = GetImageArtifact(image,"convolve:scale");
4246 if ( artifact != (const char *)NULL ) {
4247 if ( curr_kernel == kernel )
4248 curr_kernel = CloneKernelInfo(kernel);
4249 if (curr_kernel == (KernelInfo *) NULL) {
4250 curr_kernel=DestroyKernelInfo(curr_kernel);
4251 return((Image *) NULL);
4253 ScaleGeometryKernelInfo(curr_kernel, artifact);
4257 /* display the (normalized) kernel via stderr */
4258 if ( IsMagickTrue(GetImageArtifact(image,"showkernel"))
4259 || IsMagickTrue(GetImageArtifact(image,"convolve:showkernel"))
4260 || IsMagickTrue(GetImageArtifact(image,"morphology:showkernel")) )
4261 ShowKernelInfo(curr_kernel);
4263 /* Override the default handling of multi-kernel morphology results
4264 * If 'Undefined' use the default method
4265 * If 'None' (default for 'Convolve') re-iterate previous result
4266 * Otherwise merge resulting images using compose method given.
4267 * Default for 'HitAndMiss' is 'Lighten'.
4271 artifact = GetImageArtifact(image,"morphology:compose");
4272 compose = UndefinedCompositeOp; /* use default for method */
4273 if ( artifact != (const char *) NULL)
4274 compose = (CompositeOperator) ParseCommandOption(
4275 MagickComposeOptions,MagickFalse,artifact);
4277 /* Apply the Morphology */
4278 morphology_image = MorphologyApply(image, channel, method, iterations,
4279 curr_kernel, compose, image->bias, exception);
4281 /* Cleanup and Exit */
4282 if ( curr_kernel != kernel )
4283 curr_kernel=DestroyKernelInfo(curr_kernel);
4284 return(morphology_image);
4287 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
4288 method, const ssize_t iterations,const KernelInfo *kernel, ExceptionInfo
4294 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
4295 iterations,kernel,exception);
4296 return(morphology_image);
4300 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4304 + R o t a t e K e r n e l I n f o %
4308 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4310 % RotateKernelInfo() rotates the kernel by the angle given.
4312 % Currently it is restricted to 90 degree angles, of either 1D kernels
4313 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4314 % It will ignore usless rotations for specific 'named' built-in kernels.
4316 % The format of the RotateKernelInfo method is:
4318 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4320 % A description of each parameter follows:
4322 % o kernel: the Morphology/Convolution kernel
4324 % o angle: angle to rotate in degrees
4326 % This function is currently internal to this module only, but can be exported
4327 % to other modules if needed.
4329 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4331 /* angle the lower kernels first */
4332 if ( kernel->next != (KernelInfo *) NULL)
4333 RotateKernelInfo(kernel->next, angle);
4335 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4337 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4340 /* Modulus the angle */
4341 angle = fmod(angle, 360.0);
4345 if ( 337.5 < angle || angle <= 22.5 )
4346 return; /* Near zero angle - no change! - At least not at this time */
4348 /* Handle special cases */
4349 switch (kernel->type) {
4350 /* These built-in kernels are cylindrical kernels, rotating is useless */
4351 case GaussianKernel:
4356 case LaplacianKernel:
4357 case ChebyshevKernel:
4358 case ManhattanKernel:
4359 case EuclideanKernel:
4362 /* These may be rotatable at non-90 angles in the future */
4363 /* but simply rotating them in multiples of 90 degrees is useless */
4370 /* These only allows a +/-90 degree rotation (by transpose) */
4371 /* A 180 degree rotation is useless */
4373 case RectangleKernel:
4374 if ( 135.0 < angle && angle <= 225.0 )
4376 if ( 225.0 < angle && angle <= 315.0 )
4383 /* Attempt rotations by 45 degrees */
4384 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4386 if ( kernel->width == 3 && kernel->height == 3 )
4387 { /* Rotate a 3x3 square by 45 degree angle */
4388 MagickRealType t = kernel->values[0];
4389 kernel->values[0] = kernel->values[3];
4390 kernel->values[3] = kernel->values[6];
4391 kernel->values[6] = kernel->values[7];
4392 kernel->values[7] = kernel->values[8];
4393 kernel->values[8] = kernel->values[5];
4394 kernel->values[5] = kernel->values[2];
4395 kernel->values[2] = kernel->values[1];
4396 kernel->values[1] = t;
4397 /* rotate non-centered origin */
4398 if ( kernel->x != 1 || kernel->y != 1 ) {
4400 x = (ssize_t) kernel->x-1;
4401 y = (ssize_t) kernel->y-1;
4402 if ( x == y ) x = 0;
4403 else if ( x == 0 ) x = -y;
4404 else if ( x == -y ) y = 0;
4405 else if ( y == 0 ) y = x;
4406 kernel->x = (ssize_t) x+1;
4407 kernel->y = (ssize_t) y+1;
4409 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4410 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4413 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4415 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4417 if ( kernel->width == 1 || kernel->height == 1 )
4418 { /* Do a transpose of a 1 dimensional kernel,
4419 ** which results in a fast 90 degree rotation of some type.
4423 t = (ssize_t) kernel->width;
4424 kernel->width = kernel->height;
4425 kernel->height = (size_t) t;
4427 kernel->x = kernel->y;
4429 if ( kernel->width == 1 ) {
4430 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4431 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4433 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4434 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4437 else if ( kernel->width == kernel->height )
4438 { /* Rotate a square array of values by 90 degrees */
4441 register MagickRealType
4444 for( i=0, x=kernel->width-1; i<=x; i++, x--)
4445 for( j=0, y=kernel->height-1; j<y; j++, y--)
4446 { t = k[i+j*kernel->width];
4447 k[i+j*kernel->width] = k[j+x*kernel->width];
4448 k[j+x*kernel->width] = k[x+y*kernel->width];
4449 k[x+y*kernel->width] = k[y+i*kernel->width];
4450 k[y+i*kernel->width] = t;
4453 /* rotate the origin - relative to center of array */
4454 { register ssize_t x,y;
4455 x = (ssize_t) (kernel->x*2-kernel->width+1);
4456 y = (ssize_t) (kernel->y*2-kernel->height+1);
4457 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4458 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4460 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4461 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4464 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4466 if ( 135.0 < angle && angle <= 225.0 )
4468 /* For a 180 degree rotation - also know as a reflection
4469 * This is actually a very very common operation!
4470 * Basically all that is needed is a reversal of the kernel data!
4471 * And a reflection of the origon
4479 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
4480 t=k[i], k[i]=k[j], k[j]=t;
4482 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4483 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4484 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4485 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4487 /* At this point angle should at least between -45 (315) and +45 degrees
4488 * In the future some form of non-orthogonal angled rotates could be
4489 * performed here, posibily with a linear kernel restriction.
4496 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4500 % S c a l e G e o m e t r y K e r n e l I n f o %
4504 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4506 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4507 % provided as a "-set option:convolve:scale {geometry}" user setting,
4508 % and modifies the kernel according to the parsed arguments of that setting.
4510 % The first argument (and any normalization flags) are passed to
4511 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4512 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4513 % into the scaled/normalized kernel.
4515 % The format of the ScaleGeometryKernelInfo method is:
4517 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4518 % const double scaling_factor,const MagickStatusType normalize_flags)
4520 % A description of each parameter follows:
4522 % o kernel: the Morphology/Convolution kernel to modify
4525 % The geometry string to parse, typically from the user provided
4526 % "-set option:convolve:scale {geometry}" setting.
4529 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4530 const char *geometry)
4537 SetGeometryInfo(&args);
4538 flags = (GeometryFlags) ParseGeometry(geometry, &args);
4541 /* For Debugging Geometry Input */
4542 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4543 flags, args.rho, args.sigma, args.xi, args.psi );
4546 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4547 args.rho *= 0.01, args.sigma *= 0.01;
4549 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4551 if ( (flags & SigmaValue) == 0 )
4554 /* Scale/Normalize the input kernel */
4555 ScaleKernelInfo(kernel, args.rho, flags);
4557 /* Add Unity Kernel, for blending with original */
4558 if ( (flags & SigmaValue) != 0 )
4559 UnityAddKernelInfo(kernel, args.sigma);
4564 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4568 % S c a l e K e r n e l I n f o %
4572 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4574 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4575 % without normalization of the sum of the kernel values (as per given flags).
4577 % By default (no flags given) the values within the kernel is scaled
4578 % directly using given scaling factor without change.
4580 % If either of the two 'normalize_flags' are given the kernel will first be
4581 % normalized and then further scaled by the scaling factor value given.
4583 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4584 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4585 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4586 % non-HDRI versions of IM this may cause images to have any negative results
4587 % clipped, unless some 'bias' is used.
4589 % More specifically. Kernels which only contain positive values (such as a
4590 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4591 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4593 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4594 % the kernel will be scaled by the absolute of the sum of kernel values, so
4595 % that it will generally fall within the +/- 1.0 range.
4597 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4598 % will be scaled by just the sum of the postive values, so that its output
4599 % range will again fall into the +/- 1.0 range.
4601 % For special kernels designed for locating shapes using 'Correlate', (often
4602 % only containing +1 and -1 values, representing foreground/brackground
4603 % matching) a special normalization method is provided to scale the positive
4604 % values separately to those of the negative values, so the kernel will be
4605 % forced to become a zero-sum kernel better suited to such searches.
4607 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4608 % attributes within the kernel structure have been correctly set during the
4611 % NOTE: The values used for 'normalize_flags' have been selected specifically
4612 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4613 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4615 % The format of the ScaleKernelInfo method is:
4617 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4618 % const MagickStatusType normalize_flags )
4620 % A description of each parameter follows:
4622 % o kernel: the Morphology/Convolution kernel
4625 % multiply all values (after normalization) by this factor if not
4626 % zero. If the kernel is normalized regardless of any flags.
4628 % o normalize_flags:
4629 % GeometryFlags defining normalization method to use.
4630 % specifically: NormalizeValue, CorrelateNormalizeValue,
4631 % and/or PercentValue
4634 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4635 const double scaling_factor,const GeometryFlags normalize_flags)
4644 /* do the other kernels in a multi-kernel list first */
4645 if ( kernel->next != (KernelInfo *) NULL)
4646 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4648 /* Normalization of Kernel */
4650 if ( (normalize_flags&NormalizeValue) != 0 ) {
4651 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4652 /* non-zero-summing kernel (generally positive) */
4653 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4655 /* zero-summing kernel */
4656 pos_scale = kernel->positive_range;
4658 /* Force kernel into a normalized zero-summing kernel */
4659 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4660 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4661 ? kernel->positive_range : 1.0;
4662 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
4663 ? -kernel->negative_range : 1.0;
4666 neg_scale = pos_scale;
4668 /* finialize scaling_factor for positive and negative components */
4669 pos_scale = scaling_factor/pos_scale;
4670 neg_scale = scaling_factor/neg_scale;
4672 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4673 if ( ! IsNan(kernel->values[i]) )
4674 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4676 /* convolution output range */
4677 kernel->positive_range *= pos_scale;
4678 kernel->negative_range *= neg_scale;
4679 /* maximum and minimum values in kernel */
4680 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4681 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4683 /* swap kernel settings if user's scaling factor is negative */
4684 if ( scaling_factor < MagickEpsilon ) {
4686 t = kernel->positive_range;
4687 kernel->positive_range = kernel->negative_range;
4688 kernel->negative_range = t;
4689 t = kernel->maximum;
4690 kernel->maximum = kernel->minimum;
4691 kernel->minimum = 1;
4698 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4702 % S h o w K e r n e l I n f o %
4706 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4708 % ShowKernelInfo() outputs the details of the given kernel defination to
4709 % standard error, generally due to a users 'showkernel' option request.
4711 % The format of the ShowKernel method is:
4713 % void ShowKernelInfo(KernelInfo *kernel)
4715 % A description of each parameter follows:
4717 % o kernel: the Morphology/Convolution kernel
4720 MagickExport void ShowKernelInfo(KernelInfo *kernel)
4728 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4730 (void) FormatLocaleFile(stderr, "Kernel");
4731 if ( kernel->next != (KernelInfo *) NULL )
4732 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4733 (void) FormatLocaleFile(stderr, " \"%s",
4734 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4735 if ( fabs(k->angle) > MagickEpsilon )
4736 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4737 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4738 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4739 (void) FormatLocaleFile(stderr,
4740 " with values from %.*lg to %.*lg\n",
4741 GetMagickPrecision(), k->minimum,
4742 GetMagickPrecision(), k->maximum);
4743 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4744 GetMagickPrecision(), k->negative_range,
4745 GetMagickPrecision(), k->positive_range);
4746 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4747 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4748 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4749 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4751 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4752 GetMagickPrecision(), k->positive_range+k->negative_range);
4753 for (i=v=0; v < k->height; v++) {
4754 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4755 for (u=0; u < k->width; u++, i++)
4756 if ( IsNan(k->values[i]) )
4757 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4759 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4760 GetMagickPrecision(), k->values[i]);
4761 (void) FormatLocaleFile(stderr,"\n");
4767 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4771 % U n i t y A d d K e r n a l I n f o %
4775 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4777 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4778 % to the given pre-scaled and normalized Kernel. This in effect adds that
4779 % amount of the original image into the resulting convolution kernel. This
4780 % value is usually provided by the user as a percentage value in the
4781 % 'convolve:scale' setting.
4783 % The resulting effect is to convert the defined kernels into blended
4784 % soft-blurs, unsharp kernels or into sharpening kernels.
4786 % The format of the UnityAdditionKernelInfo method is:
4788 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4790 % A description of each parameter follows:
4792 % o kernel: the Morphology/Convolution kernel
4795 % scaling factor for the unity kernel to be added to
4799 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4802 /* do the other kernels in a multi-kernel list first */
4803 if ( kernel->next != (KernelInfo *) NULL)
4804 UnityAddKernelInfo(kernel->next, scale);
4806 /* Add the scaled unity kernel to the existing kernel */
4807 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4808 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4814 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4818 % Z e r o K e r n e l N a n s %
4822 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4824 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4825 % the kernel with a zero value. This is typically done when the kernel will
4826 % be used in special hardware (GPU) convolution processors, to simply
4829 % The format of the ZeroKernelNans method is:
4831 % void ZeroKernelNans (KernelInfo *kernel)
4833 % A description of each parameter follows:
4835 % o kernel: the Morphology/Convolution kernel
4838 MagickExport void ZeroKernelNans(KernelInfo *kernel)
4843 /* do the other kernels in a multi-kernel list first */
4844 if ( kernel->next != (KernelInfo *) NULL)
4845 ZeroKernelNans(kernel->next);
4847 for (i=0; i < (kernel->width*kernel->height); i++)
4848 if ( IsNan(kernel->values[i]) )
4849 kernel->values[i] = 0.0;