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 MorphologyImage() (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 MorphologyImage() 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 ssize_t iterations,const KernelInfo *kernel,
2463 % const CompositeMethod compose,const double bias,
2464 % 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 iterations: apply the operation this many times (or no change).
2473 % A value of -1 means loop until no change found.
2474 % How this is applied may depend on the morphology method.
2475 % Typically this is a value of 1.
2477 % o channel: the channel type.
2479 % o kernel: An array of double representing the morphology kernel.
2481 % o compose: How to handle or merge multi-kernel results.
2482 % If 'UndefinedCompositeOp' use default for the Morphology method.
2483 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2484 % Otherwise merge the results using the compose method given.
2486 % o bias: Convolution Output Bias.
2488 % o exception: return any errors or warnings in this structure.
2492 /* Apply a Morphology Primative to an image using the given kernel.
2493 ** Two pre-created images must be provided, and no image is created.
2494 ** It returns the number of pixels that changed between the images
2495 ** for result convergence determination.
2497 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2498 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2499 ExceptionInfo *exception)
2501 #define MorphologyTag "Morphology/Image"
2520 assert(image != (Image *) NULL);
2521 assert(image->signature == MagickSignature);
2522 assert(morphology_image != (Image *) NULL);
2523 assert(morphology_image->signature == MagickSignature);
2524 assert(kernel != (KernelInfo *) NULL);
2525 assert(kernel->signature == MagickSignature);
2526 assert(exception != (ExceptionInfo *) NULL);
2527 assert(exception->signature == MagickSignature);
2533 image_view=AcquireCacheView(image);
2534 morphology_view=AcquireCacheView(morphology_image);
2535 virt_width=image->columns+kernel->width-1;
2537 /* Some methods (including convolve) needs use a reflected kernel.
2538 * Adjust 'origin' offsets to loop though kernel as a reflection.
2543 case ConvolveMorphology:
2544 case DilateMorphology:
2545 case DilateIntensityMorphology:
2546 /*case DistanceMorphology:*/
2547 /* kernel needs to used with reflection about origin */
2548 offx = (ssize_t) kernel->width-offx-1;
2549 offy = (ssize_t) kernel->height-offy-1;
2551 case ErodeMorphology:
2552 case ErodeIntensityMorphology:
2553 case HitAndMissMorphology:
2554 case ThinningMorphology:
2555 case ThickenMorphology:
2556 /* kernel is used as is, without reflection */
2559 assert("Not a Primitive Morphology Method" != (char *) NULL);
2564 if ( method == ConvolveMorphology && kernel->width == 1 )
2565 { /* Special handling (for speed) of vertical (blur) kernels.
2566 ** This performs its handling in columns rather than in rows.
2567 ** This is only done for convolve as it is the only method that
2568 ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
2570 ** Timing tests (on single CPU laptop)
2571 ** Using a vertical 1-d Blue with normal row-by-row (below)
2572 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2574 ** Using this column method
2575 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2578 ** Anthony Thyssen, 14 June 2010
2583 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2584 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2586 for (x=0; x < (ssize_t) image->columns; x++)
2588 register const Quantum
2600 if (status == MagickFalse)
2602 p=GetCacheViewVirtualPixels(image_view,x,-offy,1,image->rows+
2603 kernel->height-1,exception);
2604 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2605 morphology_image->rows,exception);
2606 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2611 /* offset to origin in 'p'. while 'q' points to it directly */
2614 for (y=0; y < (ssize_t) image->rows; y++)
2619 register const double
2622 register const Quantum
2628 /* Copy input image to the output image for unused channels
2629 * This removes need for 'cloning' a new image every iteration
2631 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2632 GetPixelComponents(image)),q);
2633 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2634 GetPixelComponents(image)),q);
2635 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2636 GetPixelComponents(image)),q);
2637 if (image->colorspace == CMYKColorspace)
2638 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2639 GetPixelComponents(image)),q);
2641 /* Set the bias of the weighted average output */
2646 result.black = bias;
2649 /* Weighted Average of pixels using reflected kernel
2651 ** NOTE for correct working of this operation for asymetrical
2652 ** kernels, the kernel needs to be applied in its reflected form.
2653 ** That is its values needs to be reversed.
2655 k = &kernel->values[ kernel->height-1 ];
2657 if ( (image->sync == MagickFalse) || (image->matte == MagickFalse) )
2658 { /* No 'Sync' involved.
2659 ** Convolution is simple greyscale channel operation
2661 for (v=0; v < (ssize_t) kernel->height; v++) {
2662 if ( IsNan(*k) ) continue;
2663 result.red += (*k)*GetPixelRed(image,k_pixels);
2664 result.green += (*k)*GetPixelGreen(image,k_pixels);
2665 result.blue += (*k)*GetPixelBlue(image,k_pixels);
2666 if (image->colorspace == CMYKColorspace)
2667 result.black+=(*k)*GetPixelBlack(image,k_pixels);
2668 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2672 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
2673 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
2674 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
2675 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
2676 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
2677 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
2678 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
2679 (image->colorspace == CMYKColorspace))
2680 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
2681 if (((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0) &&
2682 (image->matte == MagickTrue))
2683 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2686 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2687 ** Weight the color channels with Alpha Channel so that
2688 ** transparent pixels are not part of the results.
2691 alpha, /* alpha weighting of colors : kernel*alpha */
2692 gamma; /* divisor, sum of color weighting values */
2695 for (v=0; v < (ssize_t) kernel->height; v++) {
2696 if ( IsNan(*k) ) continue;
2697 alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels));
2699 result.red += alpha*GetPixelRed(image,k_pixels);
2700 result.green += alpha*GetPixelGreen(image,k_pixels);
2701 result.blue += alpha*GetPixelBlue(image,k_pixels);
2702 if (image->colorspace == CMYKColorspace)
2703 result.black += alpha*GetPixelBlack(image,k_pixels);
2704 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2708 /* Sync'ed channels, all channels are modified */
2709 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2710 SetPixelRed(morphology_image,
2711 ClampToQuantum(gamma*result.red),q);
2712 SetPixelGreen(morphology_image,
2713 ClampToQuantum(gamma*result.green),q);
2714 SetPixelBlue(morphology_image,
2715 ClampToQuantum(gamma*result.blue),q);
2716 if (image->colorspace == CMYKColorspace)
2717 SetPixelBlack(morphology_image,
2718 ClampToQuantum(gamma*result.black),q);
2719 SetPixelAlpha(morphology_image,
2720 ClampToQuantum(result.alpha),q);
2723 /* Count up changed pixels */
2724 if ((GetPixelRed(image,p+r*GetPixelComponents(image)) != GetPixelRed(morphology_image,q))
2725 || (GetPixelGreen(image,p+r*GetPixelComponents(image)) != GetPixelGreen(morphology_image,q))
2726 || (GetPixelBlue(image,p+r*GetPixelComponents(image)) != GetPixelBlue(morphology_image,q))
2727 || (GetPixelAlpha(image,p+r*GetPixelComponents(image)) != GetPixelAlpha(morphology_image,q))
2728 || ((image->colorspace == CMYKColorspace) &&
2729 (GetPixelBlack(image,p+r*GetPixelComponents(image)) != GetPixelBlack(morphology_image,q))))
2730 changed++; /* The pixel was changed in some way! */
2731 p+=GetPixelComponents(image);
2732 q+=GetPixelComponents(morphology_image);
2734 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2736 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2741 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2742 #pragma omp critical (MagickCore_MorphologyImage)
2744 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2745 if (proceed == MagickFalse)
2749 morphology_image->type=image->type;
2750 morphology_view=DestroyCacheView(morphology_view);
2751 image_view=DestroyCacheView(image_view);
2752 return(status ? (ssize_t) changed : 0);
2756 ** Normal handling of horizontal or rectangular kernels (row by row)
2758 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2759 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2761 for (y=0; y < (ssize_t) image->rows; y++)
2763 register const Quantum
2775 if (status == MagickFalse)
2777 p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
2778 kernel->height, exception);
2779 q=GetCacheViewAuthenticPixels(morphology_view,0,y,
2780 morphology_image->columns,1,exception);
2781 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2786 /* offset to origin in 'p'. while 'q' points to it directly */
2787 r = virt_width*offy + offx;
2789 for (x=0; x < (ssize_t) image->columns; x++)
2797 register const double
2800 register const Quantum
2808 /* Copy input image to the output image for unused channels
2809 * This removes need for 'cloning' a new image every iteration
2811 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2812 GetPixelComponents(image)),q);
2813 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2814 GetPixelComponents(image)),q);
2815 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2816 GetPixelComponents(image)),q);
2817 if (image->colorspace == CMYKColorspace)
2818 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2819 GetPixelComponents(image)),q);
2826 min.black = (MagickRealType) QuantumRange;
2831 max.black = (MagickRealType) 0;
2832 /* default result is the original pixel value */
2833 result.red = (MagickRealType) GetPixelRed(image,p+r*GetPixelComponents(image));
2834 result.green = (MagickRealType) GetPixelGreen(image,p+r*GetPixelComponents(image));
2835 result.blue = (MagickRealType) GetPixelBlue(image,p+r*GetPixelComponents(image));
2837 if (image->colorspace == CMYKColorspace)
2838 result.black = (MagickRealType) GetPixelBlack(image,p+r*GetPixelComponents(image));
2839 result.alpha=(MagickRealType) GetPixelAlpha(image,p+r*GetPixelComponents(image));
2842 case ConvolveMorphology:
2843 /* Set the bias of the weighted average output */
2848 result.black = bias;
2850 case DilateIntensityMorphology:
2851 case ErodeIntensityMorphology:
2852 /* use a boolean flag indicating when first match found */
2853 result.red = 0.0; /* result is not used otherwise */
2860 case ConvolveMorphology:
2861 /* Weighted Average of pixels using reflected kernel
2863 ** NOTE for correct working of this operation for asymetrical
2864 ** kernels, the kernel needs to be applied in its reflected form.
2865 ** That is its values needs to be reversed.
2867 ** Correlation is actually the same as this but without reflecting
2868 ** the kernel, and thus 'lower-level' that Convolution. However
2869 ** as Convolution is the more common method used, and it does not
2870 ** really cost us much in terms of processing to use a reflected
2871 ** kernel, so it is Convolution that is implemented.
2873 ** Correlation will have its kernel reflected before calling
2874 ** this function to do a Convolve.
2876 ** For more details of Correlation vs Convolution see
2877 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2879 k = &kernel->values[ kernel->width*kernel->height-1 ];
2881 if ( (image->sync == MagickFalse) ||
2882 (image->matte == MagickFalse) )
2883 { /* No 'Sync' involved.
2884 ** Convolution is simple greyscale channel operation
2886 for (v=0; v < (ssize_t) kernel->height; v++) {
2887 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2888 if ( IsNan(*k) ) continue;
2890 GetPixelRed(image,k_pixels+u*GetPixelComponents(image));
2891 result.green += (*k)*
2892 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image));
2893 result.blue += (*k)*
2894 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image));
2895 if (image->colorspace == CMYKColorspace)
2896 result.black += (*k)*
2897 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image));
2898 result.alpha += (*k)*
2899 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image));
2901 k_pixels += virt_width*GetPixelComponents(image);
2903 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
2904 SetPixelRed(morphology_image,ClampToQuantum(result.red),
2906 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
2907 SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2909 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
2910 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2912 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
2913 (image->colorspace == CMYKColorspace))
2914 SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2916 if (((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0) &&
2917 (image->matte == MagickTrue))
2918 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
2922 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2923 ** Weight the color channels with Alpha Channel so that
2924 ** transparent pixels are not part of the results.
2927 alpha, /* alpha weighting of colors : kernel*alpha */
2928 gamma; /* divisor, sum of color weighting values */
2931 for (v=0; v < (ssize_t) kernel->height; v++) {
2932 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2933 if ( IsNan(*k) ) continue;
2934 alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels+u*
2935 GetPixelComponents(image)));
2937 result.red += alpha*
2938 GetPixelRed(image,k_pixels+u*GetPixelComponents(image));
2939 result.green += alpha*
2940 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image));
2941 result.blue += alpha*
2942 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image));
2943 if (image->colorspace == CMYKColorspace)
2944 result.black+=alpha*
2945 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image));
2946 result.alpha += (*k)*
2947 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image));
2949 k_pixels += virt_width*GetPixelComponents(image);
2951 /* Sync'ed channels, all channels are modified */
2952 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2953 SetPixelRed(morphology_image,
2954 ClampToQuantum(gamma*result.red),q);
2955 SetPixelGreen(morphology_image,
2956 ClampToQuantum(gamma*result.green),q);
2957 SetPixelBlue(morphology_image,
2958 ClampToQuantum(gamma*result.blue),q);
2959 if (image->colorspace == CMYKColorspace)
2960 SetPixelBlack(morphology_image,
2961 ClampToQuantum(gamma*result.black),q);
2962 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2966 case ErodeMorphology:
2967 /* Minimum Value within kernel neighbourhood
2969 ** NOTE that the kernel is not reflected for this operation!
2971 ** NOTE: in normal Greyscale Morphology, the kernel value should
2972 ** be added to the real value, this is currently not done, due to
2973 ** the nature of the boolean kernels being used.
2977 for (v=0; v < (ssize_t) kernel->height; v++) {
2978 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2979 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2980 Minimize(min.red, (double)
2981 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
2982 Minimize(min.green, (double)
2983 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
2984 Minimize(min.blue, (double)
2985 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
2986 if (image->colorspace == CMYKColorspace)
2987 Minimize(min.black,(double)
2988 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
2989 Minimize(min.alpha,(double)
2990 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
2992 k_pixels += virt_width*GetPixelComponents(image);
2996 case DilateMorphology:
2997 /* Maximum Value within kernel neighbourhood
2999 ** NOTE for correct working of this operation for asymetrical
3000 ** kernels, the kernel needs to be applied in its reflected form.
3001 ** That is its values needs to be reversed.
3003 ** NOTE: in normal Greyscale Morphology, the kernel value should
3004 ** be added to the real value, this is currently not done, due to
3005 ** the nature of the boolean kernels being used.
3008 k = &kernel->values[ kernel->width*kernel->height-1 ];
3010 for (v=0; v < (ssize_t) kernel->height; v++) {
3011 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3012 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3013 Maximize(max.red, (double)
3014 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3015 Maximize(max.green, (double)
3016 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3017 Maximize(max.blue, (double)
3018 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3019 if (image->colorspace == CMYKColorspace)
3020 Maximize(max.black, (double)
3021 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3022 Maximize(max.alpha,(double)
3023 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3025 k_pixels += virt_width*GetPixelComponents(image);
3029 case HitAndMissMorphology:
3030 case ThinningMorphology:
3031 case ThickenMorphology:
3032 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3034 ** NOTE that the kernel is not reflected for this operation,
3035 ** and consists of both foreground and background pixel
3036 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3037 ** with either Nan or 0.5 values for don't care.
3039 ** Note that this will never produce a meaningless negative
3040 ** result. Such results can cause Thinning/Thicken to not work
3041 ** correctly when used against a greyscale image.
3045 for (v=0; v < (ssize_t) kernel->height; v++) {
3046 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3047 if ( IsNan(*k) ) continue;
3049 { /* minimim of foreground pixels */
3050 Minimize(min.red, (double)
3051 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3052 Minimize(min.green, (double)
3053 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3054 Minimize(min.blue, (double)
3055 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3056 if ( image->colorspace == CMYKColorspace)
3057 Minimize(min.black,(double)
3058 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3059 Minimize(min.alpha,(double)
3060 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3062 else if ( (*k) < 0.3 )
3063 { /* maximum of background pixels */
3064 Maximize(max.red, (double)
3065 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3066 Maximize(max.green, (double)
3067 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3068 Maximize(max.blue, (double)
3069 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3070 if (image->colorspace == CMYKColorspace)
3071 Maximize(max.black, (double)
3072 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3073 Maximize(max.alpha,(double)
3074 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3077 k_pixels += virt_width*GetPixelComponents(image);
3079 /* Pattern Match if difference is positive */
3080 min.red -= max.red; Maximize( min.red, 0.0 );
3081 min.green -= max.green; Maximize( min.green, 0.0 );
3082 min.blue -= max.blue; Maximize( min.blue, 0.0 );
3083 min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3084 min.black -= max.black; Maximize( min.black, 0.0 );
3087 case ErodeIntensityMorphology:
3088 /* Select Pixel with Minimum Intensity within kernel neighbourhood
3090 ** WARNING: the intensity test fails for CMYK and does not
3091 ** take into account the moderating effect of the alpha channel
3092 ** on the intensity.
3094 ** NOTE that the kernel is not reflected for this operation!
3098 for (v=0; v < (ssize_t) kernel->height; v++) {
3099 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3100 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3101 if ( result.red == 0.0 ||
3102 GetPixelIntensity(image,&(k_pixels[u])) < GetPixelIntensity(morphology_image,q) ) {
3103 /* copy the whole pixel - no channel selection */
3104 SetPixelRed(morphology_image,GetPixelRed(image,
3105 k_pixels+u*GetPixelComponents(image)),q);
3106 SetPixelGreen(morphology_image,GetPixelGreen(image,
3107 k_pixels+u*GetPixelComponents(image)),q);
3108 SetPixelBlue(morphology_image,GetPixelBlue(image,
3109 k_pixels+u*GetPixelComponents(image)),q);
3110 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3111 k_pixels+u*GetPixelComponents(image)),q);
3112 if ( result.red > 0.0 ) changed++;
3116 k_pixels += virt_width*GetPixelComponents(image);
3120 case DilateIntensityMorphology:
3121 /* Select Pixel with Maximum Intensity within kernel neighbourhood
3123 ** WARNING: the intensity test fails for CMYK and does not
3124 ** take into account the moderating effect of the alpha channel
3125 ** on the intensity (yet).
3127 ** NOTE for correct working of this operation for asymetrical
3128 ** kernels, the kernel needs to be applied in its reflected form.
3129 ** That is its values needs to be reversed.
3131 k = &kernel->values[ kernel->width*kernel->height-1 ];
3133 for (v=0; v < (ssize_t) kernel->height; v++) {
3134 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3135 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3136 if ( result.red == 0.0 ||
3137 GetPixelIntensity(image,&(k_pixels[u])) > GetPixelIntensity(morphology_image,q) ) {
3138 /* copy the whole pixel - no channel selection */
3139 SetPixelRed(morphology_image,GetPixelRed(image,
3140 k_pixels+u*GetPixelComponents(image)),q);
3141 SetPixelGreen(morphology_image,GetPixelGreen(image,
3142 k_pixels+u*GetPixelComponents(image)),q);
3143 SetPixelBlue(morphology_image,GetPixelBlue(image,
3144 k_pixels+u*GetPixelComponents(image)),q);
3145 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3146 k_pixels+u*GetPixelComponents(image)),q);
3147 if ( result.red > 0.0 ) changed++;
3151 k_pixels += virt_width*GetPixelComponents(image);
3155 This code has been obsoleted by the MorphologyPrimitiveDirect() function.
3156 However it is still (almost) correct coding for Grayscale Morphology.
3159 GrayErode is equivalent but with kernel values subtracted from pixels
3160 without the kernel rotation
3161 GreyDilate is equivalent but using Maximum() instead of Minimum()
3162 useing kernel rotation
3164 case DistanceMorphology:
3165 /* Add kernel Value and select the minimum value found.
3166 ** The result is a iterative distance from edge of image shape.
3168 ** All Distance Kernels are symetrical, but that may not always
3169 ** be the case. For example how about a distance from left edges?
3170 ** To work correctly with asymetrical kernels the reflected kernel
3171 ** needs to be applied.
3173 k = &kernel->values[ kernel->width*kernel->height-1 ];
3175 for (v=0; v < (ssize_t) kernel->height; v++) {
3176 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3177 if ( IsNan(*k) ) continue;
3178 Minimize(result.red, (*k)+k_pixels[u].red);
3179 Minimize(result.green, (*k)+k_pixels[u].green);
3180 Minimize(result.blue, (*k)+k_pixels[u].blue);
3181 Minimize(result.alpha, (*k)+k_pixels[u].alpha);
3182 if ( image->colorspace == CMYKColorspace)
3183 Minimize(result.black,(*k)+GetPixelBlack(p_image,k_indexes+u));
3185 k_pixels += virt_width*GetPixelComponents(image);
3189 case UndefinedMorphology:
3191 break; /* Do nothing */
3193 /* Final mathematics of results (combine with original image?)
3195 ** NOTE: Difference Morphology operators Edge* and *Hat could also
3196 ** be done here but works better with iteration as a image difference
3197 ** in the controling function (below). Thicken and Thinning however
3198 ** should be done here so thay can be iterated correctly.
3201 case HitAndMissMorphology:
3202 case ErodeMorphology:
3203 result = min; /* minimum of neighbourhood */
3205 case DilateMorphology:
3206 result = max; /* maximum of neighbourhood */
3208 case ThinningMorphology:
3209 /* subtract pattern match from original */
3210 result.red -= min.red;
3211 result.green -= min.green;
3212 result.blue -= min.blue;
3213 result.alpha -= min.alpha;
3214 result.black -= min.black;
3216 case ThickenMorphology:
3217 /* Add the pattern matchs to the original */
3218 result.red += min.red;
3219 result.green += min.green;
3220 result.blue += min.blue;
3221 result.alpha += min.alpha;
3222 result.black += min.black;
3225 /* result directly calculated or assigned */
3228 /* Assign the resulting pixel values - Clamping Result */
3230 case UndefinedMorphology:
3231 case ConvolveMorphology:
3232 case DilateIntensityMorphology:
3233 case ErodeIntensityMorphology:
3234 break; /* full pixel was directly assigned - not a channel method */
3236 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
3237 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
3238 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
3239 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
3240 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
3241 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
3242 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
3243 (image->colorspace == CMYKColorspace))
3244 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
3245 if (((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0) &&
3246 (image->matte == MagickTrue))
3247 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3250 /* Count up changed pixels */
3251 if ((GetPixelRed(image,p+r*GetPixelComponents(image)) != GetPixelRed(morphology_image,q)) ||
3252 (GetPixelGreen(image,p+r*GetPixelComponents(image)) != GetPixelGreen(morphology_image,q)) ||
3253 (GetPixelBlue(image,p+r*GetPixelComponents(image)) != GetPixelBlue(morphology_image,q)) ||
3254 (GetPixelAlpha(image,p+r*GetPixelComponents(image)) != GetPixelAlpha(morphology_image,q)) ||
3255 ((image->colorspace == CMYKColorspace) &&
3256 (GetPixelBlack(image,p+r*GetPixelComponents(image)) != GetPixelBlack(morphology_image,q))))
3257 changed++; /* The pixel was changed in some way! */
3258 p+=GetPixelComponents(image);
3259 q+=GetPixelComponents(morphology_image);
3261 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3263 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3268 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3269 #pragma omp critical (MagickCore_MorphologyImage)
3271 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3272 if (proceed == MagickFalse)
3276 morphology_view=DestroyCacheView(morphology_view);
3277 image_view=DestroyCacheView(image_view);
3278 return(status ? (ssize_t)changed : -1);
3281 /* This is almost identical to the MorphologyPrimative() function above,
3282 ** but will apply the primitive directly to the image in two passes.
3284 ** That is after each row is 'Sync'ed' into the image, the next row will
3285 ** make use of those values as part of the calculation of the next row.
3286 ** It then repeats, but going in the oppisite (bottom-up) direction.
3288 ** Because of this 'iterative' handling this function can not make use
3289 ** of multi-threaded, parellel processing.
3291 static ssize_t MorphologyPrimitiveDirect(Image *image,
3292 const MorphologyMethod method,const KernelInfo *kernel,
3293 ExceptionInfo *exception)
3316 assert(image != (Image *) NULL);
3317 assert(image->signature == MagickSignature);
3318 assert(kernel != (KernelInfo *) NULL);
3319 assert(kernel->signature == MagickSignature);
3320 assert(exception != (ExceptionInfo *) NULL);
3321 assert(exception->signature == MagickSignature);
3323 /* Some methods (including convolve) needs use a reflected kernel.
3324 * Adjust 'origin' offsets to loop though kernel as a reflection.
3329 case DistanceMorphology:
3330 case VoronoiMorphology:
3331 /* kernel needs to used with reflection about origin */
3332 offx = (ssize_t) kernel->width-offx-1;
3333 offy = (ssize_t) kernel->height-offy-1;
3336 case ?????Morphology:
3337 /* kernel is used as is, without reflection */
3341 assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
3345 /* DO NOT THREAD THIS CODE! */
3346 /* two views into same image (virtual, and actual) */
3347 virt_view=AcquireCacheView(image);
3348 auth_view=AcquireCacheView(image);
3349 virt_width=image->columns+kernel->width-1;
3351 for (y=0; y < (ssize_t) image->rows; y++)
3353 register const Quantum
3365 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3366 ** we read using virtual to get virtual pixel handling, but write back
3367 ** into the same image.
3369 ** Only top half of kernel is processed as we do a single pass downward
3370 ** through the image iterating the distance function as we go.
3372 if (status == MagickFalse)
3374 p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
3376 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3378 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3380 if (status == MagickFalse)
3383 /* offset to origin in 'p'. while 'q' points to it directly */
3384 r = (ssize_t) virt_width*offy + offx;
3386 for (x=0; x < (ssize_t) image->columns; x++)
3394 register const double
3397 register const Quantum
3403 /* Starting Defaults */
3404 GetPixelInfo(image,&result);
3405 SetPixelInfo(image,q,&result);
3406 if ( method != VoronoiMorphology )
3407 result.alpha = QuantumRange - result.alpha;
3410 case DistanceMorphology:
3411 /* Add kernel Value and select the minimum value found. */
3412 k = &kernel->values[ kernel->width*kernel->height-1 ];
3414 for (v=0; v <= (ssize_t) offy; v++) {
3415 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3416 if ( IsNan(*k) ) continue;
3417 Minimize(result.red, (*k)+
3418 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3419 Minimize(result.green, (*k)+
3420 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3421 Minimize(result.blue, (*k)+
3422 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3423 if (image->colorspace == CMYKColorspace)
3424 Minimize(result.black,(*k)+
3425 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3426 Minimize(result.alpha, (*k)+
3427 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3429 k_pixels += virt_width*GetPixelComponents(image);
3431 /* repeat with the just processed pixels of this row */
3432 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3434 for (u=0; u < (ssize_t) offx; u++, k--) {
3435 if ( x+u-offx < 0 ) continue; /* off the edge! */
3436 if ( IsNan(*k) ) continue;
3437 Minimize(result.red, (*k)+
3438 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3439 Minimize(result.green, (*k)+
3440 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3441 Minimize(result.blue, (*k)+
3442 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3443 if (image->colorspace == CMYKColorspace)
3444 Minimize(result.black,(*k)+
3445 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3446 Minimize(result.alpha,(*k)+
3447 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3450 case VoronoiMorphology:
3451 /* Apply Distance to 'Matte' channel, coping the closest color.
3453 ** This is experimental, and realy the 'alpha' component should
3454 ** be completely separate 'masking' channel.
3456 k = &kernel->values[ kernel->width*kernel->height-1 ];
3458 for (v=0; v <= (ssize_t) offy; v++) {
3459 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3460 if ( IsNan(*k) ) continue;
3461 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)) )
3463 SetPixelInfo(image,k_pixels+u*GetPixelComponents(image),
3468 k_pixels += virt_width*GetPixelComponents(image);
3470 /* repeat with the just processed pixels of this row */
3471 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3473 for (u=0; u < (ssize_t) offx; u++, k--) {
3474 if ( x+u-offx < 0 ) continue; /* off the edge! */
3475 if ( IsNan(*k) ) continue;
3476 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)) )
3478 SetPixelInfo(image,k_pixels+u*GetPixelComponents(image),
3485 /* result directly calculated or assigned */
3488 /* Assign the resulting pixel values - Clamping Result */
3490 case VoronoiMorphology:
3491 SetPixelPixelInfo(image,&result,q);
3494 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
3495 SetPixelRed(image,ClampToQuantum(result.red),q);
3496 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
3497 SetPixelGreen(image,ClampToQuantum(result.green),q);
3498 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
3499 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3500 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
3501 (image->colorspace == CMYKColorspace))
3502 SetPixelBlack(image,ClampToQuantum(result.black),q);
3503 if ((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0 &&
3504 (image->matte == MagickTrue))
3505 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3508 /* Count up changed pixels */
3509 if ((GetPixelRed(image,p+r*GetPixelComponents(image)) != GetPixelRed(image,q)) ||
3510 (GetPixelGreen(image,p+r*GetPixelComponents(image)) != GetPixelGreen(image,q)) ||
3511 (GetPixelBlue(image,p+r*GetPixelComponents(image)) != GetPixelBlue(image,q)) ||
3512 (GetPixelAlpha(image,p+r*GetPixelComponents(image)) != GetPixelAlpha(image,q)) ||
3513 ((image->colorspace == CMYKColorspace) &&
3514 (GetPixelBlack(image,p+r*GetPixelComponents(image)) != GetPixelBlack(image,q))))
3515 changed++; /* The pixel was changed in some way! */
3517 p+=GetPixelComponents(image); /* increment pixel buffers */
3518 q+=GetPixelComponents(image);
3521 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3523 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3524 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3530 /* Do the reversed pass through the image */
3531 for (y=(ssize_t)image->rows-1; y >= 0; y--)
3533 register const Quantum
3545 if (status == MagickFalse)
3547 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3548 ** we read using virtual to get virtual pixel handling, but write back
3549 ** into the same image.
3551 ** Only the bottom half of the kernel will be processes as we
3554 p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
3555 kernel->y+1,exception);
3556 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3558 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3560 if (status == MagickFalse)
3563 /* adjust positions to end of row */
3564 p += (image->columns-1)*GetPixelComponents(image);
3565 q += (image->columns-1)*GetPixelComponents(image);
3567 /* offset to origin in 'p'. while 'q' points to it directly */
3570 for (x=(ssize_t)image->columns-1; x >= 0; x--)
3578 register const double
3581 register const Quantum
3587 /* Default - previously modified pixel */
3588 GetPixelInfo(image,&result);
3589 SetPixelInfo(image,q,&result);
3590 if ( method != VoronoiMorphology )
3591 result.alpha = QuantumRange - result.alpha;
3594 case DistanceMorphology:
3595 /* Add kernel Value and select the minimum value found. */
3596 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3598 for (v=offy; v < (ssize_t) kernel->height; v++) {
3599 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3600 if ( IsNan(*k) ) continue;
3601 Minimize(result.red, (*k)+
3602 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3603 Minimize(result.green, (*k)+
3604 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3605 Minimize(result.blue, (*k)+
3606 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3607 if ( image->colorspace == CMYKColorspace)
3608 Minimize(result.black,(*k)+
3609 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3610 Minimize(result.alpha, (*k)+
3611 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3613 k_pixels += virt_width*GetPixelComponents(image);
3615 /* repeat with the just processed pixels of this row */
3616 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3617 k_pixels = q-offx*GetPixelComponents(image);
3618 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3619 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3620 if ( IsNan(*k) ) continue;
3621 Minimize(result.red, (*k)+
3622 GetPixelRed(image,k_pixels+u*GetPixelComponents(image)));
3623 Minimize(result.green, (*k)+
3624 GetPixelGreen(image,k_pixels+u*GetPixelComponents(image)));
3625 Minimize(result.blue, (*k)+
3626 GetPixelBlue(image,k_pixels+u*GetPixelComponents(image)));
3627 if ( image->colorspace == CMYKColorspace)
3628 Minimize(result.black, (*k)+
3629 GetPixelBlack(image,k_pixels+u*GetPixelComponents(image)));
3630 Minimize(result.alpha, (*k)+
3631 GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)));
3634 case VoronoiMorphology:
3635 /* Apply Distance to 'Matte' channel, coping the closest color.
3637 ** This is experimental, and realy the 'alpha' component should
3638 ** be completely separate 'masking' channel.
3640 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3642 for (v=offy; v < (ssize_t) kernel->height; v++) {
3643 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3644 if ( IsNan(*k) ) continue;
3645 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)) )
3647 SetPixelInfo(image,k_pixels+u*GetPixelComponents(image),
3652 k_pixels += virt_width*GetPixelComponents(image);
3654 /* repeat with the just processed pixels of this row */
3655 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3656 k_pixels = q-offx*GetPixelComponents(image);
3657 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3658 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3659 if ( IsNan(*k) ) continue;
3660 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelComponents(image)) )
3662 SetPixelInfo(image,k_pixels+u*GetPixelComponents(image),
3669 /* result directly calculated or assigned */
3672 /* Assign the resulting pixel values - Clamping Result */
3674 case VoronoiMorphology:
3675 SetPixelPixelInfo(image,&result,q);
3678 if ((GetPixelRedTraits(image) & ActivePixelTrait) != 0)
3679 SetPixelRed(image,ClampToQuantum(result.red),q);
3680 if ((GetPixelGreenTraits(image) & ActivePixelTrait) != 0)
3681 SetPixelGreen(image,ClampToQuantum(result.green),q);
3682 if ((GetPixelBlueTraits(image) & ActivePixelTrait) != 0)
3683 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3684 if (((GetPixelBlackTraits(image) & ActivePixelTrait) != 0) &&
3685 (image->colorspace == CMYKColorspace))
3686 SetPixelBlack(image,ClampToQuantum(result.black),q);
3687 if ((GetPixelAlphaTraits(image) & ActivePixelTrait) != 0 &&
3688 (image->matte == MagickTrue))
3689 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3692 /* Count up changed pixels */
3693 if ( (GetPixelRed(image,p+r*GetPixelComponents(image)) != GetPixelRed(image,q))
3694 || (GetPixelGreen(image,p+r*GetPixelComponents(image)) != GetPixelGreen(image,q))
3695 || (GetPixelBlue(image,p+r*GetPixelComponents(image)) != GetPixelBlue(image,q))
3696 || (GetPixelAlpha(image,p+r*GetPixelComponents(image)) != GetPixelAlpha(image,q))
3697 || ((image->colorspace == CMYKColorspace) &&
3698 (GetPixelBlack(image,p+r*GetPixelComponents(image)) != GetPixelBlack(image,q))))
3699 changed++; /* The pixel was changed in some way! */
3701 p-=GetPixelComponents(image); /* go backward through pixel buffers */
3702 q-=GetPixelComponents(image);
3704 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3706 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3707 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3713 auth_view=DestroyCacheView(auth_view);
3714 virt_view=DestroyCacheView(virt_view);
3715 return(status ? (ssize_t) changed : -1);
3718 /* Apply a Morphology by calling theabove low level primitive application
3719 ** functions. This function handles any iteration loops, composition or
3720 ** re-iteration of results, and compound morphology methods that is based
3721 ** on multiple low-level (staged) morphology methods.
3723 ** Basically this provides the complex grue between the requested morphology
3724 ** method and raw low-level implementation (above).
3726 MagickExport Image *MorphologyApply(const Image *image,
3727 const MorphologyMethod method, const ssize_t iterations,
3728 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3729 ExceptionInfo *exception)
3735 *curr_image, /* Image we are working with or iterating */
3736 *work_image, /* secondary image for primitive iteration */
3737 *save_image, /* saved image - for 'edge' method only */
3738 *rslt_image; /* resultant image - after multi-kernel handling */
3741 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3742 *norm_kernel, /* the current normal un-reflected kernel */
3743 *rflt_kernel, /* the current reflected kernel (if needed) */
3744 *this_kernel; /* the kernel being applied */
3747 primitive; /* the current morphology primitive being applied */
3750 rslt_compose; /* multi-kernel compose method for results to use */
3753 special, /* do we use a direct modify function? */
3754 verbose; /* verbose output of results */
3757 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3758 method_limit, /* maximum number of compound method iterations */
3759 kernel_number, /* Loop 2: the kernel number being applied */
3760 stage_loop, /* Loop 3: primitive loop for compound morphology */
3761 stage_limit, /* how many primitives are in this compound */
3762 kernel_loop, /* Loop 4: iterate the kernel over image */
3763 kernel_limit, /* number of times to iterate kernel */
3764 count, /* total count of primitive steps applied */
3765 kernel_changed, /* total count of changed using iterated kernel */
3766 method_changed; /* total count of changed over method iteration */
3769 changed; /* number pixels changed by last primitive operation */
3774 assert(image != (Image *) NULL);
3775 assert(image->signature == MagickSignature);
3776 assert(kernel != (KernelInfo *) NULL);
3777 assert(kernel->signature == MagickSignature);
3778 assert(exception != (ExceptionInfo *) NULL);
3779 assert(exception->signature == MagickSignature);
3781 count = 0; /* number of low-level morphology primitives performed */
3782 if ( iterations == 0 )
3783 return((Image *)NULL); /* null operation - nothing to do! */
3785 kernel_limit = (size_t) iterations;
3786 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3787 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3789 verbose = IsMagickTrue(GetImageArtifact(image,"verbose"));
3791 /* initialise for cleanup */
3792 curr_image = (Image *) image;
3793 curr_compose = image->compose;
3794 (void) curr_compose;
3795 work_image = save_image = rslt_image = (Image *) NULL;
3796 reflected_kernel = (KernelInfo *) NULL;
3798 /* Initialize specific methods
3799 * + which loop should use the given iteratations
3800 * + how many primitives make up the compound morphology
3801 * + multi-kernel compose method to use (by default)
3803 method_limit = 1; /* just do method once, unless otherwise set */
3804 stage_limit = 1; /* assume method is not a compound */
3805 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3806 rslt_compose = compose; /* and we are composing multi-kernels as given */
3808 case SmoothMorphology: /* 4 primitive compound morphology */
3811 case OpenMorphology: /* 2 primitive compound morphology */
3812 case OpenIntensityMorphology:
3813 case TopHatMorphology:
3814 case CloseMorphology:
3815 case CloseIntensityMorphology:
3816 case BottomHatMorphology:
3817 case EdgeMorphology:
3820 case HitAndMissMorphology:
3821 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3823 case ThinningMorphology:
3824 case ThickenMorphology:
3825 method_limit = kernel_limit; /* iterate the whole method */
3826 kernel_limit = 1; /* do not do kernel iteration */
3828 case DistanceMorphology:
3829 case VoronoiMorphology:
3830 special = MagickTrue;
3836 /* Apply special methods with special requirments
3837 ** For example, single run only, or post-processing requirements
3839 if ( special == MagickTrue )
3841 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3842 if (rslt_image == (Image *) NULL)
3844 if (SetImageStorageClass(rslt_image,DirectClass) == MagickFalse)
3846 InheritException(exception,&rslt_image->exception);
3850 changed = MorphologyPrimitiveDirect(rslt_image, method,
3853 if ( verbose == MagickTrue )
3854 (void) (void) FormatLocaleFile(stderr,
3855 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3856 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3857 1.0,0.0,1.0, (double) changed);
3862 if ( method == VoronoiMorphology ) {
3863 /* Preserve the alpha channel of input image - but turned off */
3864 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel);
3865 (void) CompositeImage(rslt_image, CopyOpacityCompositeOp, image, 0, 0);
3866 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel);
3871 /* Handle user (caller) specified multi-kernel composition method */
3872 if ( compose != UndefinedCompositeOp )
3873 rslt_compose = compose; /* override default composition for method */
3874 if ( rslt_compose == UndefinedCompositeOp )
3875 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3877 /* Some methods require a reflected kernel to use with primitives.
3878 * Create the reflected kernel for those methods. */
3880 case CorrelateMorphology:
3881 case CloseMorphology:
3882 case CloseIntensityMorphology:
3883 case BottomHatMorphology:
3884 case SmoothMorphology:
3885 reflected_kernel = CloneKernelInfo(kernel);
3886 if (reflected_kernel == (KernelInfo *) NULL)
3888 RotateKernelInfo(reflected_kernel,180);
3894 /* Loops around more primitive morpholgy methods
3895 ** erose, dilate, open, close, smooth, edge, etc...
3897 /* Loop 1: iterate the compound method */
3900 while ( method_loop < method_limit && method_changed > 0 ) {
3904 /* Loop 2: iterate over each kernel in a multi-kernel list */
3905 norm_kernel = (KernelInfo *) kernel;
3906 this_kernel = (KernelInfo *) kernel;
3907 rflt_kernel = reflected_kernel;
3910 while ( norm_kernel != NULL ) {
3912 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3913 stage_loop = 0; /* the compound morphology stage number */
3914 while ( stage_loop < stage_limit ) {
3915 stage_loop++; /* The stage of the compound morphology */
3917 /* Select primitive morphology for this stage of compound method */
3918 this_kernel = norm_kernel; /* default use unreflected kernel */
3919 primitive = method; /* Assume method is a primitive */
3921 case ErodeMorphology: /* just erode */
3922 case EdgeInMorphology: /* erode and image difference */
3923 primitive = ErodeMorphology;
3925 case DilateMorphology: /* just dilate */
3926 case EdgeOutMorphology: /* dilate and image difference */
3927 primitive = DilateMorphology;
3929 case OpenMorphology: /* erode then dialate */
3930 case TopHatMorphology: /* open and image difference */
3931 primitive = ErodeMorphology;
3932 if ( stage_loop == 2 )
3933 primitive = DilateMorphology;
3935 case OpenIntensityMorphology:
3936 primitive = ErodeIntensityMorphology;
3937 if ( stage_loop == 2 )
3938 primitive = DilateIntensityMorphology;
3940 case CloseMorphology: /* dilate, then erode */
3941 case BottomHatMorphology: /* close and image difference */
3942 this_kernel = rflt_kernel; /* use the reflected kernel */
3943 primitive = DilateMorphology;
3944 if ( stage_loop == 2 )
3945 primitive = ErodeMorphology;
3947 case CloseIntensityMorphology:
3948 this_kernel = rflt_kernel; /* use the reflected kernel */
3949 primitive = DilateIntensityMorphology;
3950 if ( stage_loop == 2 )
3951 primitive = ErodeIntensityMorphology;
3953 case SmoothMorphology: /* open, close */
3954 switch ( stage_loop ) {
3955 case 1: /* start an open method, which starts with Erode */
3956 primitive = ErodeMorphology;
3958 case 2: /* now Dilate the Erode */
3959 primitive = DilateMorphology;
3961 case 3: /* Reflect kernel a close */
3962 this_kernel = rflt_kernel; /* use the reflected kernel */
3963 primitive = DilateMorphology;
3965 case 4: /* Finish the Close */
3966 this_kernel = rflt_kernel; /* use the reflected kernel */
3967 primitive = ErodeMorphology;
3971 case EdgeMorphology: /* dilate and erode difference */
3972 primitive = DilateMorphology;
3973 if ( stage_loop == 2 ) {
3974 save_image = curr_image; /* save the image difference */
3975 curr_image = (Image *) image;
3976 primitive = ErodeMorphology;
3979 case CorrelateMorphology:
3980 /* A Correlation is a Convolution with a reflected kernel.
3981 ** However a Convolution is a weighted sum using a reflected
3982 ** kernel. It may seem stange to convert a Correlation into a
3983 ** Convolution as the Correlation is the simplier method, but
3984 ** Convolution is much more commonly used, and it makes sense to
3985 ** implement it directly so as to avoid the need to duplicate the
3986 ** kernel when it is not required (which is typically the
3989 this_kernel = rflt_kernel; /* use the reflected kernel */
3990 primitive = ConvolveMorphology;
3995 assert( this_kernel != (KernelInfo *) NULL );
3997 /* Extra information for debugging compound operations */
3998 if ( verbose == MagickTrue ) {
3999 if ( stage_limit > 1 )
4000 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
4001 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
4002 method_loop,(double) stage_loop);
4003 else if ( primitive != method )
4004 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
4005 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
4011 /* Loop 4: Iterate the kernel with primitive */
4015 while ( kernel_loop < kernel_limit && changed > 0 ) {
4016 kernel_loop++; /* the iteration of this kernel */
4018 /* Create a clone as the destination image, if not yet defined */
4019 if ( work_image == (Image *) NULL )
4021 work_image=CloneImage(image,0,0,MagickTrue,exception);
4022 if (work_image == (Image *) NULL)
4024 if (SetImageStorageClass(work_image,DirectClass) == MagickFalse)
4026 InheritException(exception,&work_image->exception);
4029 /* work_image->type=image->type; ??? */
4032 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4034 changed = MorphologyPrimitive(curr_image, work_image, primitive,
4035 this_kernel, bias, exception);
4037 if ( verbose == MagickTrue ) {
4038 if ( kernel_loop > 1 )
4039 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4040 (void) (void) FormatLocaleFile(stderr,
4041 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4042 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4043 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4044 (double) (method_loop+kernel_loop-1),(double) kernel_number,
4045 (double) count,(double) changed);
4049 kernel_changed += changed;
4050 method_changed += changed;
4052 /* prepare next loop */
4053 { Image *tmp = work_image; /* swap images for iteration */
4054 work_image = curr_image;
4057 if ( work_image == image )
4058 work_image = (Image *) NULL; /* replace input 'image' */
4060 } /* End Loop 4: Iterate the kernel with primitive */
4062 if ( verbose == MagickTrue && kernel_changed != (size_t)changed )
4063 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
4064 if ( verbose == MagickTrue && stage_loop < stage_limit )
4065 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4068 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4069 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
4070 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
4071 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
4072 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
4075 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4077 /* Final Post-processing for some Compound Methods
4079 ** The removal of any 'Sync' channel flag in the Image Compositon
4080 ** below ensures the methematical compose method is applied in a
4081 ** purely mathematical way, and only to the selected channels.
4082 ** Turn off SVG composition 'alpha blending'.
4085 case EdgeOutMorphology:
4086 case EdgeInMorphology:
4087 case TopHatMorphology:
4088 case BottomHatMorphology:
4089 if ( verbose == MagickTrue )
4090 (void) FormatLocaleFile(stderr, "\n%s: Difference with original image",
4091 CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4092 curr_image->sync=MagickFalse;
4093 (void) CompositeImage(curr_image,DifferenceCompositeOp,image,0,0);
4095 case EdgeMorphology:
4096 if ( verbose == MagickTrue )
4097 (void) FormatLocaleFile(stderr, "\n%s: Difference of Dilate and Erode",
4098 CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4099 curr_image->sync=MagickFalse;
4100 (void) CompositeImage(curr_image,DifferenceCompositeOp,save_image,0,
4102 save_image = DestroyImage(save_image); /* finished with save image */
4108 /* multi-kernel handling: re-iterate, or compose results */
4109 if ( kernel->next == (KernelInfo *) NULL )
4110 rslt_image = curr_image; /* just return the resulting image */
4111 else if ( rslt_compose == NoCompositeOp )
4112 { if ( verbose == MagickTrue ) {
4113 if ( this_kernel->next != (KernelInfo *) NULL )
4114 (void) FormatLocaleFile(stderr, " (re-iterate)");
4116 (void) FormatLocaleFile(stderr, " (done)");
4118 rslt_image = curr_image; /* return result, and re-iterate */
4120 else if ( rslt_image == (Image *) NULL)
4121 { if ( verbose == MagickTrue )
4122 (void) FormatLocaleFile(stderr, " (save for compose)");
4123 rslt_image = curr_image;
4124 curr_image = (Image *) image; /* continue with original image */
4127 { /* Add the new 'current' result to the composition
4129 ** The removal of any 'Sync' channel flag in the Image Compositon
4130 ** below ensures the methematical compose method is applied in a
4131 ** purely mathematical way, and only to the selected channels.
4132 ** IE: Turn off SVG composition 'alpha blending'.
4134 if ( verbose == MagickTrue )
4135 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4136 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4137 rslt_image->sync=MagickFalse;
4138 (void) CompositeImage(rslt_image, rslt_compose, curr_image, 0, 0);
4139 curr_image = DestroyImage(curr_image);
4140 curr_image = (Image *) image; /* continue with original image */
4142 if ( verbose == MagickTrue )
4143 (void) FormatLocaleFile(stderr, "\n");
4145 /* loop to the next kernel in a multi-kernel list */
4146 norm_kernel = norm_kernel->next;
4147 if ( rflt_kernel != (KernelInfo *) NULL )
4148 rflt_kernel = rflt_kernel->next;
4150 } /* End Loop 2: Loop over each kernel */
4152 } /* End Loop 1: compound method interation */
4156 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4158 if ( curr_image == rslt_image )
4159 curr_image = (Image *) NULL;
4160 if ( rslt_image != (Image *) NULL )
4161 rslt_image = DestroyImage(rslt_image);
4163 if ( curr_image == rslt_image || curr_image == image )
4164 curr_image = (Image *) NULL;
4165 if ( curr_image != (Image *) NULL )
4166 curr_image = DestroyImage(curr_image);
4167 if ( work_image != (Image *) NULL )
4168 work_image = DestroyImage(work_image);
4169 if ( save_image != (Image *) NULL )
4170 save_image = DestroyImage(save_image);
4171 if ( reflected_kernel != (KernelInfo *) NULL )
4172 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4178 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4182 % M o r p h o l o g y I m a g e %
4186 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4188 % MorphologyImage() applies a user supplied kernel to the image
4189 % according to the given mophology method.
4191 % This function applies any and all user defined settings before calling
4192 % the above internal function MorphologyApply().
4194 % User defined settings include...
4195 % * Output Bias for Convolution and correlation ("-bias")
4196 % * Kernel Scale/normalize settings ("-set 'option:convolve:scale'")
4197 % This can also includes the addition of a scaled unity kernel.
4198 % * Show Kernel being applied ("-set option:showkernel 1")
4200 % The format of the MorphologyImage method is:
4202 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4203 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4205 % Image *MorphologyImage(const Image *image, const ChannelType
4206 % channel,MorphologyMethod method,const ssize_t iterations,
4207 % KernelInfo *kernel,ExceptionInfo *exception)
4209 % A description of each parameter follows:
4211 % o image: the image.
4213 % o method: the morphology method to be applied.
4215 % o iterations: apply the operation this many times (or no change).
4216 % A value of -1 means loop until no change found.
4217 % How this is applied may depend on the morphology method.
4218 % Typically this is a value of 1.
4220 % o kernel: An array of double representing the morphology kernel.
4221 % Warning: kernel may be normalized for the Convolve method.
4223 % o exception: return any errors or warnings in this structure.
4226 MagickExport Image *MorphologyImage(const Image *image,
4227 const MorphologyMethod method,const ssize_t iterations,
4228 const KernelInfo *kernel,ExceptionInfo *exception)
4240 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4241 * This is done BEFORE the ShowKernelInfo() function is called so that
4242 * users can see the results of the 'option:convolve:scale' option.
4244 curr_kernel = (KernelInfo *) kernel;
4245 if ( method == ConvolveMorphology || method == CorrelateMorphology )
4249 artifact = GetImageArtifact(image,"convolve:scale");
4250 if ( artifact != (const char *)NULL ) {
4251 if ( curr_kernel == kernel )
4252 curr_kernel = CloneKernelInfo(kernel);
4253 if (curr_kernel == (KernelInfo *) NULL) {
4254 curr_kernel=DestroyKernelInfo(curr_kernel);
4255 return((Image *) NULL);
4257 ScaleGeometryKernelInfo(curr_kernel, artifact);
4261 /* display the (normalized) kernel via stderr */
4262 if ( IsMagickTrue(GetImageArtifact(image,"showkernel"))
4263 || IsMagickTrue(GetImageArtifact(image,"convolve:showkernel"))
4264 || IsMagickTrue(GetImageArtifact(image,"morphology:showkernel")) )
4265 ShowKernelInfo(curr_kernel);
4267 /* Override the default handling of multi-kernel morphology results
4268 * If 'Undefined' use the default method
4269 * If 'None' (default for 'Convolve') re-iterate previous result
4270 * Otherwise merge resulting images using compose method given.
4271 * Default for 'HitAndMiss' is 'Lighten'.
4275 artifact = GetImageArtifact(image,"morphology:compose");
4276 compose = UndefinedCompositeOp; /* use default for method */
4277 if ( artifact != (const char *) NULL)
4278 compose = (CompositeOperator) ParseCommandOption(
4279 MagickComposeOptions,MagickFalse,artifact);
4281 /* Apply the Morphology */
4282 morphology_image = MorphologyApply(image, method, iterations,
4283 curr_kernel, compose, image->bias, exception);
4285 /* Cleanup and Exit */
4286 if ( curr_kernel != kernel )
4287 curr_kernel=DestroyKernelInfo(curr_kernel);
4288 return(morphology_image);
4292 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4296 + R o t a t e K e r n e l I n f o %
4300 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4302 % RotateKernelInfo() rotates the kernel by the angle given.
4304 % Currently it is restricted to 90 degree angles, of either 1D kernels
4305 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4306 % It will ignore usless rotations for specific 'named' built-in kernels.
4308 % The format of the RotateKernelInfo method is:
4310 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4312 % A description of each parameter follows:
4314 % o kernel: the Morphology/Convolution kernel
4316 % o angle: angle to rotate in degrees
4318 % This function is currently internal to this module only, but can be exported
4319 % to other modules if needed.
4321 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4323 /* angle the lower kernels first */
4324 if ( kernel->next != (KernelInfo *) NULL)
4325 RotateKernelInfo(kernel->next, angle);
4327 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4329 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4332 /* Modulus the angle */
4333 angle = fmod(angle, 360.0);
4337 if ( 337.5 < angle || angle <= 22.5 )
4338 return; /* Near zero angle - no change! - At least not at this time */
4340 /* Handle special cases */
4341 switch (kernel->type) {
4342 /* These built-in kernels are cylindrical kernels, rotating is useless */
4343 case GaussianKernel:
4348 case LaplacianKernel:
4349 case ChebyshevKernel:
4350 case ManhattanKernel:
4351 case EuclideanKernel:
4354 /* These may be rotatable at non-90 angles in the future */
4355 /* but simply rotating them in multiples of 90 degrees is useless */
4362 /* These only allows a +/-90 degree rotation (by transpose) */
4363 /* A 180 degree rotation is useless */
4365 case RectangleKernel:
4366 if ( 135.0 < angle && angle <= 225.0 )
4368 if ( 225.0 < angle && angle <= 315.0 )
4375 /* Attempt rotations by 45 degrees */
4376 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4378 if ( kernel->width == 3 && kernel->height == 3 )
4379 { /* Rotate a 3x3 square by 45 degree angle */
4380 MagickRealType t = kernel->values[0];
4381 kernel->values[0] = kernel->values[3];
4382 kernel->values[3] = kernel->values[6];
4383 kernel->values[6] = kernel->values[7];
4384 kernel->values[7] = kernel->values[8];
4385 kernel->values[8] = kernel->values[5];
4386 kernel->values[5] = kernel->values[2];
4387 kernel->values[2] = kernel->values[1];
4388 kernel->values[1] = t;
4389 /* rotate non-centered origin */
4390 if ( kernel->x != 1 || kernel->y != 1 ) {
4392 x = (ssize_t) kernel->x-1;
4393 y = (ssize_t) kernel->y-1;
4394 if ( x == y ) x = 0;
4395 else if ( x == 0 ) x = -y;
4396 else if ( x == -y ) y = 0;
4397 else if ( y == 0 ) y = x;
4398 kernel->x = (ssize_t) x+1;
4399 kernel->y = (ssize_t) y+1;
4401 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4402 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4405 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4407 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4409 if ( kernel->width == 1 || kernel->height == 1 )
4410 { /* Do a transpose of a 1 dimensional kernel,
4411 ** which results in a fast 90 degree rotation of some type.
4415 t = (ssize_t) kernel->width;
4416 kernel->width = kernel->height;
4417 kernel->height = (size_t) t;
4419 kernel->x = kernel->y;
4421 if ( kernel->width == 1 ) {
4422 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4423 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4425 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4426 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4429 else if ( kernel->width == kernel->height )
4430 { /* Rotate a square array of values by 90 degrees */
4433 register MagickRealType
4436 for( i=0, x=kernel->width-1; i<=x; i++, x--)
4437 for( j=0, y=kernel->height-1; j<y; j++, y--)
4438 { t = k[i+j*kernel->width];
4439 k[i+j*kernel->width] = k[j+x*kernel->width];
4440 k[j+x*kernel->width] = k[x+y*kernel->width];
4441 k[x+y*kernel->width] = k[y+i*kernel->width];
4442 k[y+i*kernel->width] = t;
4445 /* rotate the origin - relative to center of array */
4446 { register ssize_t x,y;
4447 x = (ssize_t) (kernel->x*2-kernel->width+1);
4448 y = (ssize_t) (kernel->y*2-kernel->height+1);
4449 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4450 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4452 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4453 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4456 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4458 if ( 135.0 < angle && angle <= 225.0 )
4460 /* For a 180 degree rotation - also know as a reflection
4461 * This is actually a very very common operation!
4462 * Basically all that is needed is a reversal of the kernel data!
4463 * And a reflection of the origon
4471 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
4472 t=k[i], k[i]=k[j], k[j]=t;
4474 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4475 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4476 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4477 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4479 /* At this point angle should at least between -45 (315) and +45 degrees
4480 * In the future some form of non-orthogonal angled rotates could be
4481 * performed here, posibily with a linear kernel restriction.
4488 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4492 % S c a l e G e o m e t r y K e r n e l I n f o %
4496 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4498 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4499 % provided as a "-set option:convolve:scale {geometry}" user setting,
4500 % and modifies the kernel according to the parsed arguments of that setting.
4502 % The first argument (and any normalization flags) are passed to
4503 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4504 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4505 % into the scaled/normalized kernel.
4507 % The format of the ScaleGeometryKernelInfo method is:
4509 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4510 % const double scaling_factor,const MagickStatusType normalize_flags)
4512 % A description of each parameter follows:
4514 % o kernel: the Morphology/Convolution kernel to modify
4517 % The geometry string to parse, typically from the user provided
4518 % "-set option:convolve:scale {geometry}" setting.
4521 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4522 const char *geometry)
4529 SetGeometryInfo(&args);
4530 flags = (GeometryFlags) ParseGeometry(geometry, &args);
4533 /* For Debugging Geometry Input */
4534 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4535 flags, args.rho, args.sigma, args.xi, args.psi );
4538 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4539 args.rho *= 0.01, args.sigma *= 0.01;
4541 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4543 if ( (flags & SigmaValue) == 0 )
4546 /* Scale/Normalize the input kernel */
4547 ScaleKernelInfo(kernel, args.rho, flags);
4549 /* Add Unity Kernel, for blending with original */
4550 if ( (flags & SigmaValue) != 0 )
4551 UnityAddKernelInfo(kernel, args.sigma);
4556 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4560 % S c a l e K e r n e l I n f o %
4564 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4566 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4567 % without normalization of the sum of the kernel values (as per given flags).
4569 % By default (no flags given) the values within the kernel is scaled
4570 % directly using given scaling factor without change.
4572 % If either of the two 'normalize_flags' are given the kernel will first be
4573 % normalized and then further scaled by the scaling factor value given.
4575 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4576 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4577 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4578 % non-HDRI versions of IM this may cause images to have any negative results
4579 % clipped, unless some 'bias' is used.
4581 % More specifically. Kernels which only contain positive values (such as a
4582 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4583 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4585 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4586 % the kernel will be scaled by the absolute of the sum of kernel values, so
4587 % that it will generally fall within the +/- 1.0 range.
4589 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4590 % will be scaled by just the sum of the postive values, so that its output
4591 % range will again fall into the +/- 1.0 range.
4593 % For special kernels designed for locating shapes using 'Correlate', (often
4594 % only containing +1 and -1 values, representing foreground/brackground
4595 % matching) a special normalization method is provided to scale the positive
4596 % values separately to those of the negative values, so the kernel will be
4597 % forced to become a zero-sum kernel better suited to such searches.
4599 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4600 % attributes within the kernel structure have been correctly set during the
4603 % NOTE: The values used for 'normalize_flags' have been selected specifically
4604 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4605 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4607 % The format of the ScaleKernelInfo method is:
4609 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4610 % const MagickStatusType normalize_flags )
4612 % A description of each parameter follows:
4614 % o kernel: the Morphology/Convolution kernel
4617 % multiply all values (after normalization) by this factor if not
4618 % zero. If the kernel is normalized regardless of any flags.
4620 % o normalize_flags:
4621 % GeometryFlags defining normalization method to use.
4622 % specifically: NormalizeValue, CorrelateNormalizeValue,
4623 % and/or PercentValue
4626 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4627 const double scaling_factor,const GeometryFlags normalize_flags)
4636 /* do the other kernels in a multi-kernel list first */
4637 if ( kernel->next != (KernelInfo *) NULL)
4638 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4640 /* Normalization of Kernel */
4642 if ( (normalize_flags&NormalizeValue) != 0 ) {
4643 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4644 /* non-zero-summing kernel (generally positive) */
4645 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4647 /* zero-summing kernel */
4648 pos_scale = kernel->positive_range;
4650 /* Force kernel into a normalized zero-summing kernel */
4651 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4652 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4653 ? kernel->positive_range : 1.0;
4654 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
4655 ? -kernel->negative_range : 1.0;
4658 neg_scale = pos_scale;
4660 /* finialize scaling_factor for positive and negative components */
4661 pos_scale = scaling_factor/pos_scale;
4662 neg_scale = scaling_factor/neg_scale;
4664 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4665 if ( ! IsNan(kernel->values[i]) )
4666 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4668 /* convolution output range */
4669 kernel->positive_range *= pos_scale;
4670 kernel->negative_range *= neg_scale;
4671 /* maximum and minimum values in kernel */
4672 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4673 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4675 /* swap kernel settings if user's scaling factor is negative */
4676 if ( scaling_factor < MagickEpsilon ) {
4678 t = kernel->positive_range;
4679 kernel->positive_range = kernel->negative_range;
4680 kernel->negative_range = t;
4681 t = kernel->maximum;
4682 kernel->maximum = kernel->minimum;
4683 kernel->minimum = 1;
4690 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4694 % S h o w K e r n e l I n f o %
4698 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4700 % ShowKernelInfo() outputs the details of the given kernel defination to
4701 % standard error, generally due to a users 'showkernel' option request.
4703 % The format of the ShowKernel method is:
4705 % void ShowKernelInfo(KernelInfo *kernel)
4707 % A description of each parameter follows:
4709 % o kernel: the Morphology/Convolution kernel
4712 MagickExport void ShowKernelInfo(KernelInfo *kernel)
4720 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4722 (void) FormatLocaleFile(stderr, "Kernel");
4723 if ( kernel->next != (KernelInfo *) NULL )
4724 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4725 (void) FormatLocaleFile(stderr, " \"%s",
4726 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4727 if ( fabs(k->angle) > MagickEpsilon )
4728 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4729 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4730 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4731 (void) FormatLocaleFile(stderr,
4732 " with values from %.*lg to %.*lg\n",
4733 GetMagickPrecision(), k->minimum,
4734 GetMagickPrecision(), k->maximum);
4735 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4736 GetMagickPrecision(), k->negative_range,
4737 GetMagickPrecision(), k->positive_range);
4738 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4739 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4740 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4741 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4743 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4744 GetMagickPrecision(), k->positive_range+k->negative_range);
4745 for (i=v=0; v < k->height; v++) {
4746 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4747 for (u=0; u < k->width; u++, i++)
4748 if ( IsNan(k->values[i]) )
4749 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4751 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4752 GetMagickPrecision(), k->values[i]);
4753 (void) FormatLocaleFile(stderr,"\n");
4759 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4763 % U n i t y A d d K e r n a l I n f o %
4767 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4769 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4770 % to the given pre-scaled and normalized Kernel. This in effect adds that
4771 % amount of the original image into the resulting convolution kernel. This
4772 % value is usually provided by the user as a percentage value in the
4773 % 'convolve:scale' setting.
4775 % The resulting effect is to convert the defined kernels into blended
4776 % soft-blurs, unsharp kernels or into sharpening kernels.
4778 % The format of the UnityAdditionKernelInfo method is:
4780 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4782 % A description of each parameter follows:
4784 % o kernel: the Morphology/Convolution kernel
4787 % scaling factor for the unity kernel to be added to
4791 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4794 /* do the other kernels in a multi-kernel list first */
4795 if ( kernel->next != (KernelInfo *) NULL)
4796 UnityAddKernelInfo(kernel->next, scale);
4798 /* Add the scaled unity kernel to the existing kernel */
4799 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4800 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4806 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4810 % Z e r o K e r n e l N a n s %
4814 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4816 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4817 % the kernel with a zero value. This is typically done when the kernel will
4818 % be used in special hardware (GPU) convolution processors, to simply
4821 % The format of the ZeroKernelNans method is:
4823 % void ZeroKernelNans (KernelInfo *kernel)
4825 % A description of each parameter follows:
4827 % o kernel: the Morphology/Convolution kernel
4830 MagickExport void ZeroKernelNans(KernelInfo *kernel)
4835 /* do the other kernels in a multi-kernel list first */
4836 if ( kernel->next != (KernelInfo *) NULL)
4837 ZeroKernelNans(kernel->next);
4839 for (i=0; i < (kernel->width*kernel->height); i++)
4840 if ( IsNan(kernel->values[i]) )
4841 kernel->values[i] = 0.0;