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/gem-private.h"
61 #include "MagickCore/hashmap.h"
62 #include "MagickCore/image.h"
63 #include "MagickCore/image-private.h"
64 #include "MagickCore/list.h"
65 #include "MagickCore/magick.h"
66 #include "MagickCore/memory_.h"
67 #include "MagickCore/monitor-private.h"
68 #include "MagickCore/morphology.h"
69 #include "MagickCore/morphology-private.h"
70 #include "MagickCore/option.h"
71 #include "MagickCore/pixel-accessor.h"
72 #include "MagickCore/prepress.h"
73 #include "MagickCore/quantize.h"
74 #include "MagickCore/registry.h"
75 #include "MagickCore/semaphore.h"
76 #include "MagickCore/splay-tree.h"
77 #include "MagickCore/statistic.h"
78 #include "MagickCore/string_.h"
79 #include "MagickCore/string-private.h"
80 #include "MagickCore/token.h"
81 #include "MagickCore/utility.h"
82 #include "MagickCore/utility-private.h"
86 ** The following test is for special floating point numbers of value NaN (not
87 ** a number), that may be used within a Kernel Definition. NaN's are defined
88 ** as part of the IEEE standard for floating point number representation.
90 ** These are used as a Kernel value to mean that this kernel position is not
91 ** part of the kernel neighbourhood for convolution or morphology processing,
92 ** and thus should be ignored. This allows the use of 'shaped' kernels.
94 ** The special properity that two NaN's are never equal, even if they are from
95 ** the same variable allow you to test if a value is special NaN value.
97 ** This macro IsNaN() is thus is only true if the value given is NaN.
99 #define IsNan(a) ((a)!=(a))
102 Other global definitions used by module.
104 static inline double MagickMin(const double x,const double y)
106 return( x < y ? x : y);
108 static inline double MagickMax(const double x,const double y)
110 return( x > y ? x : y);
112 #define Minimize(assign,value) assign=MagickMin(assign,value)
113 #define Maximize(assign,value) assign=MagickMax(assign,value)
115 /* Currently these are only internal to this module */
117 CalcKernelMetaData(KernelInfo *),
118 ExpandMirrorKernelInfo(KernelInfo *),
119 ExpandRotateKernelInfo(KernelInfo *, const double),
120 RotateKernelInfo(KernelInfo *, double);
123 /* Quick function to find last kernel in a kernel list */
124 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
126 while (kernel->next != (KernelInfo *) NULL)
127 kernel = kernel->next;
132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
136 % A c q u i r e K e r n e l I n f o %
140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
142 % AcquireKernelInfo() takes the given string (generally supplied by the
143 % user) and converts it into a Morphology/Convolution Kernel. This allows
144 % users to specify a kernel from a number of pre-defined kernels, or to fully
145 % specify their own kernel for a specific Convolution or Morphology
148 % The kernel so generated can be any rectangular array of floating point
149 % values (doubles) with the 'control point' or 'pixel being affected'
150 % anywhere within that array of values.
152 % Previously IM was restricted to a square of odd size using the exact
153 % center as origin, this is no longer the case, and any rectangular kernel
154 % with any value being declared the origin. This in turn allows the use of
155 % highly asymmetrical kernels.
157 % The floating point values in the kernel can also include a special value
158 % known as 'nan' or 'not a number' to indicate that this value is not part
159 % of the kernel array. This allows you to shaped the kernel within its
160 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
161 % shape. However at least one non-nan value must be provided for correct
162 % working of a kernel.
164 % The returned kernel should be freed using the DestroyKernelInfo() when you
165 % are finished with it. Do not free this memory yourself.
167 % Input kernel defintion strings can consist of any of three types.
170 % Select from one of the built in kernels, using the name and
171 % geometry arguments supplied. See AcquireKernelBuiltIn()
173 % "WxH[+X+Y][@><]:num, num, num ..."
174 % a kernel of size W by H, with W*H floating point numbers following.
175 % the 'center' can be optionally be defined at +X+Y (such that +0+0
176 % is top left corner). If not defined the pixel in the center, for
177 % odd sizes, or to the immediate top or left of center for even sizes
178 % is automatically selected.
180 % "num, num, num, num, ..."
181 % list of floating point numbers defining an 'old style' odd sized
182 % square kernel. At least 9 values should be provided for a 3x3
183 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
184 % Values can be space or comma separated. This is not recommended.
186 % You can define a 'list of kernels' which can be used by some morphology
187 % operators A list is defined as a semi-colon separated list kernels.
189 % " kernel ; kernel ; kernel ; "
191 % Any extra ';' characters, at start, end or between kernel defintions are
194 % The special flags will expand a single kernel, into a list of rotated
195 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
196 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
197 % The '<' also exands using 90-degree rotates, but giving a 180-degree
198 % reflected kernel before the +/- 90-degree rotations, which can be important
199 % for Thinning operations.
201 % Note that 'name' kernels will start with an alphabetic character while the
202 % new kernel specification has a ':' character in its specification string.
203 % If neither is the case, it is assumed an old style of a simple list of
204 % numbers generating a odd-sized square kernel has been given.
206 % The format of the AcquireKernal method is:
208 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
210 % A description of each parameter follows:
212 % o kernel_string: the Morphology/Convolution kernel wanted.
216 /* This was separated so that it could be used as a separate
217 ** array input handling function, such as for -color-matrix
219 static KernelInfo *ParseKernelArray(const char *kernel_string)
225 token[MaxTextExtent];
235 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
243 kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
244 if (kernel == (KernelInfo *)NULL)
246 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
247 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
248 kernel->negative_range = kernel->positive_range = 0.0;
249 kernel->type = UserDefinedKernel;
250 kernel->next = (KernelInfo *) NULL;
251 kernel->signature = MagickSignature;
252 if (kernel_string == (const char *) NULL)
255 /* find end of this specific kernel definition string */
256 end = strchr(kernel_string, ';');
257 if ( end == (char *) NULL )
258 end = strchr(kernel_string, '\0');
260 /* clear flags - for Expanding kernel lists thorugh rotations */
263 /* Has a ':' in argument - New user kernel specification */
264 p = strchr(kernel_string, ':');
265 if ( p != (char *) NULL && p < end)
267 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
268 memcpy(token, kernel_string, (size_t) (p-kernel_string));
269 token[p-kernel_string] = '\0';
270 SetGeometryInfo(&args);
271 flags = ParseGeometry(token, &args);
273 /* Size handling and checks of geometry settings */
274 if ( (flags & WidthValue) == 0 ) /* if no width then */
275 args.rho = args.sigma; /* then width = height */
276 if ( args.rho < 1.0 ) /* if width too small */
277 args.rho = 1.0; /* then width = 1 */
278 if ( args.sigma < 1.0 ) /* if height too small */
279 args.sigma = args.rho; /* then height = width */
280 kernel->width = (size_t)args.rho;
281 kernel->height = (size_t)args.sigma;
283 /* Offset Handling and Checks */
284 if ( args.xi < 0.0 || args.psi < 0.0 )
285 return(DestroyKernelInfo(kernel));
286 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
287 : (ssize_t) (kernel->width-1)/2;
288 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
289 : (ssize_t) (kernel->height-1)/2;
290 if ( kernel->x >= (ssize_t) kernel->width ||
291 kernel->y >= (ssize_t) kernel->height )
292 return(DestroyKernelInfo(kernel));
294 p++; /* advance beyond the ':' */
297 { /* ELSE - Old old specification, forming odd-square kernel */
298 /* count up number of values given */
299 p=(const char *) kernel_string;
300 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
301 p++; /* ignore "'" chars for convolve filter usage - Cristy */
302 for (i=0; p < end; i++)
304 GetMagickToken(p,&p,token);
306 GetMagickToken(p,&p,token);
308 /* set the size of the kernel - old sized square */
309 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
310 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
311 p=(const char *) kernel_string;
312 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
313 p++; /* ignore "'" chars for convolve filter usage - Cristy */
316 /* Read in the kernel values from rest of input string argument */
317 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
318 kernel->height*sizeof(*kernel->values));
319 if (kernel->values == (double *) NULL)
320 return(DestroyKernelInfo(kernel));
321 kernel->minimum = +MagickHuge;
322 kernel->maximum = -MagickHuge;
323 kernel->negative_range = kernel->positive_range = 0.0;
324 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
326 GetMagickToken(p,&p,token);
328 GetMagickToken(p,&p,token);
329 if ( LocaleCompare("nan",token) == 0
330 || LocaleCompare("-",token) == 0 ) {
331 kernel->values[i] = nan; /* do not include this value in kernel */
334 kernel->values[i] = StringToDouble(token,(char **) NULL);
335 ( kernel->values[i] < 0)
336 ? ( kernel->negative_range += kernel->values[i] )
337 : ( kernel->positive_range += kernel->values[i] );
338 Minimize(kernel->minimum, kernel->values[i]);
339 Maximize(kernel->maximum, kernel->values[i]);
343 /* sanity check -- no more values in kernel definition */
344 GetMagickToken(p,&p,token);
345 if ( *token != '\0' && *token != ';' && *token != '\'' )
346 return(DestroyKernelInfo(kernel));
349 /* this was the old method of handling a incomplete kernel */
350 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
351 Minimize(kernel->minimum, kernel->values[i]);
352 Maximize(kernel->maximum, kernel->values[i]);
353 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
354 kernel->values[i]=0.0;
357 /* Number of values for kernel was not enough - Report Error */
358 if ( i < (ssize_t) (kernel->width*kernel->height) )
359 return(DestroyKernelInfo(kernel));
362 /* check that we recieved at least one real (non-nan) value! */
363 if ( kernel->minimum == MagickHuge )
364 return(DestroyKernelInfo(kernel));
366 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
367 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
368 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
369 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
370 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
371 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
376 static KernelInfo *ParseKernelName(const char *kernel_string)
379 token[MaxTextExtent];
397 /* Parse special 'named' kernel */
398 GetMagickToken(kernel_string,&p,token);
399 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
400 if ( type < 0 || type == UserDefinedKernel )
401 return((KernelInfo *)NULL); /* not a valid named kernel */
403 while (((isspace((int) ((unsigned char) *p)) != 0) ||
404 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
407 end = strchr(p, ';'); /* end of this kernel defintion */
408 if ( end == (char *) NULL )
409 end = strchr(p, '\0');
411 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
412 memcpy(token, p, (size_t) (end-p));
414 SetGeometryInfo(&args);
415 flags = ParseGeometry(token, &args);
418 /* For Debugging Geometry Input */
419 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
420 flags, args.rho, args.sigma, args.xi, args.psi );
423 /* special handling of missing values in input string */
425 /* Shape Kernel Defaults */
427 if ( (flags & WidthValue) == 0 )
428 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
436 if ( (flags & HeightValue) == 0 )
437 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
440 if ( (flags & XValue) == 0 )
441 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
443 case RectangleKernel: /* Rectangle - set size defaults */
444 if ( (flags & WidthValue) == 0 ) /* if no width then */
445 args.rho = args.sigma; /* then width = height */
446 if ( args.rho < 1.0 ) /* if width too small */
447 args.rho = 3; /* then width = 3 */
448 if ( args.sigma < 1.0 ) /* if height too small */
449 args.sigma = args.rho; /* then height = width */
450 if ( (flags & XValue) == 0 ) /* center offset if not defined */
451 args.xi = (double)(((ssize_t)args.rho-1)/2);
452 if ( (flags & YValue) == 0 )
453 args.psi = (double)(((ssize_t)args.sigma-1)/2);
455 /* Distance Kernel Defaults */
456 case ChebyshevKernel:
457 case ManhattanKernel:
458 case OctagonalKernel:
459 case EuclideanKernel:
460 if ( (flags & HeightValue) == 0 ) /* no distance scale */
461 args.sigma = 100.0; /* default distance scaling */
462 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
463 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
464 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
465 args.sigma *= QuantumRange/100.0; /* percentage of color range */
471 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
472 if ( kernel == (KernelInfo *) NULL )
475 /* global expand to rotated kernel list - only for single kernels */
476 if ( kernel->next == (KernelInfo *) NULL ) {
477 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
478 ExpandRotateKernelInfo(kernel, 45.0);
479 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
480 ExpandRotateKernelInfo(kernel, 90.0);
481 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
482 ExpandMirrorKernelInfo(kernel);
488 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
496 token[MaxTextExtent];
504 if (kernel_string == (const char *) NULL)
505 return(ParseKernelArray(kernel_string));
510 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
512 /* ignore extra or multiple ';' kernel separators */
513 if ( *token != ';' ) {
515 /* tokens starting with alpha is a Named kernel */
516 if (isalpha((int) *token) != 0)
517 new_kernel = ParseKernelName(p);
518 else /* otherwise a user defined kernel array */
519 new_kernel = ParseKernelArray(p);
521 /* Error handling -- this is not proper error handling! */
522 if ( new_kernel == (KernelInfo *) NULL ) {
523 (void) FormatLocaleFile(stderr, "Failed to parse kernel number #%.20g\n",
524 (double) kernel_number);
525 if ( kernel != (KernelInfo *) NULL )
526 kernel=DestroyKernelInfo(kernel);
527 return((KernelInfo *) NULL);
530 /* initialise or append the kernel list */
531 if ( kernel == (KernelInfo *) NULL )
534 LastKernelInfo(kernel)->next = new_kernel;
537 /* look for the next kernel in list */
539 if ( p == (char *) NULL )
549 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
553 % A c q u i r e K e r n e l B u i l t I n %
557 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
559 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
560 % kernels used for special purposes such as gaussian blurring, skeleton
561 % pruning, and edge distance determination.
563 % They take a KernelType, and a set of geometry style arguments, which were
564 % typically decoded from a user supplied string, or from a more complex
565 % Morphology Method that was requested.
567 % The format of the AcquireKernalBuiltIn method is:
569 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
570 % const GeometryInfo args)
572 % A description of each parameter follows:
574 % o type: the pre-defined type of kernel wanted
576 % o args: arguments defining or modifying the kernel
578 % Convolution Kernels
581 % The a No-Op or Scaling single element kernel.
583 % Gaussian:{radius},{sigma}
584 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
585 % The sigma for the curve is required. The resulting kernel is
588 % If 'sigma' is zero, you get a single pixel on a field of zeros.
590 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
591 % the final size of the resulting kernel to a square 2*radius+1 in size.
592 % The radius should be at least 2 times that of the sigma value, or
593 % sever clipping and aliasing may result. If not given or set to 0 the
594 % radius will be determined so as to produce the best minimal error
595 % result, which is usally much larger than is normally needed.
597 % LoG:{radius},{sigma}
598 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
599 % The supposed ideal edge detection, zero-summing kernel.
601 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
602 % approx 1.6 (according to wikipedia).
604 % DoG:{radius},{sigma1},{sigma2}
605 % "Difference of Gaussians" Kernel.
606 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
607 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
608 % The result is a zero-summing kernel.
610 % Blur:{radius},{sigma}[,{angle}]
611 % Generates a 1 dimensional or linear gaussian blur, at the angle given
612 % (current restricted to orthogonal angles). If a 'radius' is given the
613 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
614 % by a 90 degree angle.
616 % If 'sigma' is zero, you get a single pixel on a field of zeros.
618 % Note that two convolutions with two "Blur" kernels perpendicular to
619 % each other, is equivalent to a far larger "Gaussian" kernel with the
620 % same sigma value, However it is much faster to apply. This is how the
621 % "-blur" operator actually works.
623 % Comet:{width},{sigma},{angle}
624 % Blur in one direction only, much like how a bright object leaves
625 % a comet like trail. The Kernel is actually half a gaussian curve,
626 % Adding two such blurs in opposite directions produces a Blur Kernel.
627 % Angle can be rotated in multiples of 90 degrees.
629 % Note that the first argument is the width of the kernel and not the
630 % radius of the kernel.
632 % # Still to be implemented...
636 % # Set kernel values using a resize filter, and given scale (sigma)
637 % # Cylindrical or Linear. Is this possible with an image?
640 % Named Constant Convolution Kernels
642 % All these are unscaled, zero-summing kernels by default. As such for
643 % non-HDRI version of ImageMagick some form of normalization, user scaling,
644 % and biasing the results is recommended, to prevent the resulting image
647 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
648 % 45 degrees to generate the 8 angled varients of each of the kernels.
651 % Discrete Lapacian Kernels, (without normalization)
652 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
653 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
654 % Type 2 : 3x3 with center:4 edge:1 corner:-2
655 % Type 3 : 3x3 with center:4 edge:-2 corner:1
656 % Type 5 : 5x5 laplacian
657 % Type 7 : 7x7 laplacian
658 % Type 15 : 5x5 LoG (sigma approx 1.4)
659 % Type 19 : 9x9 LoG (sigma approx 1.4)
662 % Sobel 'Edge' convolution kernel (3x3)
668 % Roberts convolution kernel (3x3)
674 % Prewitt Edge convolution kernel (3x3)
680 % Prewitt's "Compass" convolution kernel (3x3)
686 % Kirsch's "Compass" convolution kernel (3x3)
692 % Frei-Chen Edge Detector is based on a kernel that is similar to
693 % the Sobel Kernel, but is designed to be isotropic. That is it takes
694 % into account the distance of the diagonal in the kernel.
697 % | sqrt(2), 0, -sqrt(2) |
700 % FreiChen:{type},{angle}
702 % Frei-Chen Pre-weighted kernels...
704 % Type 0: default un-nomalized version shown above.
706 % Type 1: Orthogonal Kernel (same as type 11 below)
708 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
711 % Type 2: Diagonal form of Kernel...
713 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
716 % However this kernel is als at the heart of the FreiChen Edge Detection
717 % Process which uses a set of 9 specially weighted kernel. These 9
718 % kernels not be normalized, but directly applied to the image. The
719 % results is then added together, to produce the intensity of an edge in
720 % a specific direction. The square root of the pixel value can then be
721 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
722 % from each other, both the direction and the strength of the edge can be
725 % Type 10: All 9 of the following pre-weighted kernels...
727 % Type 11: | 1, 0, -1 |
728 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
731 % Type 12: | 1, sqrt(2), 1 |
732 % | 0, 0, 0 | / 2*sqrt(2)
735 % Type 13: | sqrt(2), -1, 0 |
736 % | -1, 0, 1 | / 2*sqrt(2)
739 % Type 14: | 0, 1, -sqrt(2) |
740 % | -1, 0, 1 | / 2*sqrt(2)
743 % Type 15: | 0, -1, 0 |
747 % Type 16: | 1, 0, -1 |
751 % Type 17: | 1, -2, 1 |
755 % Type 18: | -2, 1, -2 |
759 % Type 19: | 1, 1, 1 |
763 % The first 4 are for edge detection, the next 4 are for line detection
764 % and the last is to add a average component to the results.
766 % Using a special type of '-1' will return all 9 pre-weighted kernels
767 % as a multi-kernel list, so that you can use them directly (without
768 % normalization) with the special "-set option:morphology:compose Plus"
769 % setting to apply the full FreiChen Edge Detection Technique.
771 % If 'type' is large it will be taken to be an actual rotation angle for
772 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
773 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
775 % WARNING: The above was layed out as per
776 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
777 % But rotated 90 degrees so direction is from left rather than the top.
778 % I have yet to find any secondary confirmation of the above. The only
779 % other source found was actual source code at
780 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
781 % Neigher paper defineds the kernels in a way that looks locical or
782 % correct when taken as a whole.
786 % Diamond:[{radius}[,{scale}]]
787 % Generate a diamond shaped kernel with given radius to the points.
788 % Kernel size will again be radius*2+1 square and defaults to radius 1,
789 % generating a 3x3 kernel that is slightly larger than a square.
791 % Square:[{radius}[,{scale}]]
792 % Generate a square shaped kernel of size radius*2+1, and defaulting
793 % to a 3x3 (radius 1).
795 % Octagon:[{radius}[,{scale}]]
796 % Generate octagonal shaped kernel of given radius and constant scale.
797 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
798 % in "Diamond" kernel.
800 % Disk:[{radius}[,{scale}]]
801 % Generate a binary disk, thresholded at the radius given, the radius
802 % may be a float-point value. Final Kernel size is floor(radius)*2+1
803 % square. A radius of 5.3 is the default.
805 % NOTE: That a low radii Disk kernels produce the same results as
806 % many of the previously defined kernels, but differ greatly at larger
807 % radii. Here is a table of equivalences...
808 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
809 % "Disk:1.5" => "Square"
810 % "Disk:2" => "Diamond:2"
811 % "Disk:2.5" => "Octagon"
812 % "Disk:2.9" => "Square:2"
813 % "Disk:3.5" => "Octagon:3"
814 % "Disk:4.5" => "Octagon:4"
815 % "Disk:5.4" => "Octagon:5"
816 % "Disk:6.4" => "Octagon:6"
817 % All other Disk shapes are unique to this kernel, but because a "Disk"
818 % is more circular when using a larger radius, using a larger radius is
819 % preferred over iterating the morphological operation.
821 % Rectangle:{geometry}
822 % Simply generate a rectangle of 1's with the size given. You can also
823 % specify the location of the 'control point', otherwise the closest
824 % pixel to the center of the rectangle is selected.
826 % Properly centered and odd sized rectangles work the best.
828 % Symbol Dilation Kernels
830 % These kernel is not a good general morphological kernel, but is used
831 % more for highlighting and marking any single pixels in an image using,
832 % a "Dilate" method as appropriate.
834 % For the same reasons iterating these kernels does not produce the
835 % same result as using a larger radius for the symbol.
837 % Plus:[{radius}[,{scale}]]
838 % Cross:[{radius}[,{scale}]]
839 % Generate a kernel in the shape of a 'plus' or a 'cross' with
840 % a each arm the length of the given radius (default 2).
842 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
844 % Ring:{radius1},{radius2}[,{scale}]
845 % A ring of the values given that falls between the two radii.
846 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
847 % This is the 'edge' pixels of the default "Disk" kernel,
848 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
850 % Hit and Miss Kernels
852 % Peak:radius1,radius2
853 % Find any peak larger than the pixels the fall between the two radii.
854 % The default ring of pixels is as per "Ring".
856 % Find flat orthogonal edges of a binary shape
858 % Find 90 degree corners of a binary shape
860 % A special kernel to thin the 'outside' of diagonals
862 % Find end points of lines (for pruning a skeletion)
863 % Two types of lines ends (default to both) can be searched for
864 % Type 0: All line ends
865 % Type 1: single kernel for 4-conneected line ends
866 % Type 2: single kernel for simple line ends
868 % Find three line junctions (within a skeletion)
869 % Type 0: all line junctions
870 % Type 1: Y Junction kernel
871 % Type 2: Diagonal T Junction kernel
872 % Type 3: Orthogonal T Junction kernel
873 % Type 4: Diagonal X Junction kernel
874 % Type 5: Orthogonal + Junction kernel
876 % Find single pixel ridges or thin lines
877 % Type 1: Fine single pixel thick lines and ridges
878 % Type 2: Find two pixel thick lines and ridges
880 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
882 % Traditional skeleton generating kernels.
883 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
884 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
885 % Type 3: Thinning skeleton based on a ressearch paper by
886 % Dan S. Bloomberg (Default Type)
888 % A huge variety of Thinning Kernels designed to preserve conectivity.
889 % many other kernel sets use these kernels as source definitions.
890 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
891 % the super and sub notations used in the source research paper.
893 % Distance Measuring Kernels
895 % Different types of distance measuring methods, which are used with the
896 % a 'Distance' morphology method for generating a gradient based on
897 % distance from an edge of a binary shape, though there is a technique
898 % for handling a anti-aliased shape.
900 % See the 'Distance' Morphological Method, for information of how it is
903 % Chebyshev:[{radius}][x{scale}[%!]]
904 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
905 % is a value of one to any neighbour, orthogonal or diagonal. One why
906 % of thinking of it is the number of squares a 'King' or 'Queen' in
907 % chess needs to traverse reach any other position on a chess board.
908 % It results in a 'square' like distance function, but one where
909 % diagonals are given a value that is closer than expected.
911 % Manhattan:[{radius}][x{scale}[%!]]
912 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
913 % Cab distance metric), it is the distance needed when you can only
914 % travel in horizontal or vertical directions only. It is the
915 % distance a 'Rook' in chess would have to travel, and results in a
916 % diamond like distances, where diagonals are further than expected.
918 % Octagonal:[{radius}][x{scale}[%!]]
919 % An interleving of Manhatten and Chebyshev metrics producing an
920 % increasing octagonally shaped distance. Distances matches those of
921 % the "Octagon" shaped kernel of the same radius. The minimum radius
922 % and default is 2, producing a 5x5 kernel.
924 % Euclidean:[{radius}][x{scale}[%!]]
925 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
926 % However by default the kernel size only has a radius of 1, which
927 % limits the distance to 'Knight' like moves, with only orthogonal and
928 % diagonal measurements being correct. As such for the default kernel
929 % you will get octagonal like distance function.
931 % However using a larger radius such as "Euclidean:4" you will get a
932 % much smoother distance gradient from the edge of the shape. Especially
933 % if the image is pre-processed to include any anti-aliasing pixels.
934 % Of course a larger kernel is slower to use, and not always needed.
936 % The first three Distance Measuring Kernels will only generate distances
937 % of exact multiples of {scale} in binary images. As such you can use a
938 % scale of 1 without loosing any information. However you also need some
939 % scaling when handling non-binary anti-aliased shapes.
941 % The "Euclidean" Distance Kernel however does generate a non-integer
942 % fractional results, and as such scaling is vital even for binary shapes.
946 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
947 const GeometryInfo *args)
960 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
962 /* Generate a new empty kernel if needed */
963 kernel=(KernelInfo *) NULL;
965 case UndefinedKernel: /* These should not call this function */
966 case UserDefinedKernel:
967 assert("Should not call this function" != (char *)NULL);
969 case LaplacianKernel: /* Named Descrete Convolution Kernels */
970 case SobelKernel: /* these are defined using other kernels */
976 case EdgesKernel: /* Hit and Miss kernels */
978 case DiagonalsKernel:
980 case LineJunctionsKernel:
982 case ConvexHullKernel:
985 break; /* A pre-generated kernel is not needed */
987 /* set to 1 to do a compile-time check that we haven't missed anything */
996 case RectangleKernel:
1003 case ChebyshevKernel:
1004 case ManhattanKernel:
1005 case OctangonalKernel:
1006 case EuclideanKernel:
1010 /* Generate the base Kernel Structure */
1011 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1012 if (kernel == (KernelInfo *) NULL)
1014 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1015 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1016 kernel->negative_range = kernel->positive_range = 0.0;
1017 kernel->type = type;
1018 kernel->next = (KernelInfo *) NULL;
1019 kernel->signature = MagickSignature;
1029 kernel->height = kernel->width = (size_t) 1;
1030 kernel->x = kernel->y = (ssize_t) 0;
1031 kernel->values=(double *) AcquireAlignedMemory(1,
1032 sizeof(*kernel->values));
1033 if (kernel->values == (double *) NULL)
1034 return(DestroyKernelInfo(kernel));
1035 kernel->maximum = kernel->values[0] = args->rho;
1039 case GaussianKernel:
1043 sigma = fabs(args->sigma),
1044 sigma2 = fabs(args->xi),
1047 if ( args->rho >= 1.0 )
1048 kernel->width = (size_t)args->rho*2+1;
1049 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1050 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1052 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1053 kernel->height = kernel->width;
1054 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1055 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1056 kernel->height*sizeof(*kernel->values));
1057 if (kernel->values == (double *) NULL)
1058 return(DestroyKernelInfo(kernel));
1060 /* WARNING: The following generates a 'sampled gaussian' kernel.
1061 * What we really want is a 'discrete gaussian' kernel.
1063 * How to do this is I don't know, but appears to be basied on the
1064 * Error Function 'erf()' (intergral of a gaussian)
1067 if ( type == GaussianKernel || type == DoGKernel )
1068 { /* Calculate a Gaussian, OR positive half of a DoG */
1069 if ( sigma > MagickEpsilon )
1070 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1071 B = (double) (1.0/(Magick2PI*sigma*sigma));
1072 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1073 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1074 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1076 else /* limiting case - a unity (normalized Dirac) kernel */
1077 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1078 kernel->width*kernel->height*sizeof(double));
1079 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1083 if ( type == DoGKernel )
1084 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1085 if ( sigma2 > MagickEpsilon )
1086 { sigma = sigma2; /* simplify loop expressions */
1087 A = 1.0/(2.0*sigma*sigma);
1088 B = (double) (1.0/(Magick2PI*sigma*sigma));
1089 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1090 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1091 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1093 else /* limiting case - a unity (normalized Dirac) kernel */
1094 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1097 if ( type == LoGKernel )
1098 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1099 if ( sigma > MagickEpsilon )
1100 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1101 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1102 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1103 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1104 { R = ((double)(u*u+v*v))*A;
1105 kernel->values[i] = (1-R)*exp(-R)*B;
1108 else /* special case - generate a unity kernel */
1109 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1110 kernel->width*kernel->height*sizeof(double));
1111 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1115 /* Note the above kernels may have been 'clipped' by a user defined
1116 ** radius, producing a smaller (darker) kernel. Also for very small
1117 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1118 ** producing a very bright kernel.
1120 ** Normalization will still be needed.
1123 /* Normalize the 2D Gaussian Kernel
1125 ** NB: a CorrelateNormalize performs a normal Normalize if
1126 ** there are no negative values.
1128 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1129 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1135 sigma = fabs(args->sigma),
1138 if ( args->rho >= 1.0 )
1139 kernel->width = (size_t)args->rho*2+1;
1141 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1143 kernel->x = (ssize_t) (kernel->width-1)/2;
1145 kernel->negative_range = kernel->positive_range = 0.0;
1146 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1147 kernel->height*sizeof(*kernel->values));
1148 if (kernel->values == (double *) NULL)
1149 return(DestroyKernelInfo(kernel));
1152 #define KernelRank 3
1153 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1154 ** It generates a gaussian 3 times the width, and compresses it into
1155 ** the expected range. This produces a closer normalization of the
1156 ** resulting kernel, especially for very low sigma values.
1157 ** As such while wierd it is prefered.
1159 ** I am told this method originally came from Photoshop.
1161 ** A properly normalized curve is generated (apart from edge clipping)
1162 ** even though we later normalize the result (for edge clipping)
1163 ** to allow the correct generation of a "Difference of Blurs".
1167 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1168 (void) ResetMagickMemory(kernel->values,0, (size_t)
1169 kernel->width*kernel->height*sizeof(double));
1170 /* Calculate a Positive 1D Gaussian */
1171 if ( sigma > MagickEpsilon )
1172 { sigma *= KernelRank; /* simplify loop expressions */
1173 alpha = 1.0/(2.0*sigma*sigma);
1174 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1175 for ( u=-v; u <= v; u++) {
1176 kernel->values[(u+v)/KernelRank] +=
1177 exp(-((double)(u*u))*alpha)*beta;
1180 else /* special case - generate a unity kernel */
1181 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1183 /* Direct calculation without curve averaging */
1185 /* Calculate a Positive Gaussian */
1186 if ( sigma > MagickEpsilon )
1187 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1188 beta = 1.0/(MagickSQ2PI*sigma);
1189 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1190 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1192 else /* special case - generate a unity kernel */
1193 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1194 kernel->width*kernel->height*sizeof(double));
1195 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1198 /* Note the above kernel may have been 'clipped' by a user defined
1199 ** radius, producing a smaller (darker) kernel. Also for very small
1200 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1201 ** producing a very bright kernel.
1203 ** Normalization will still be needed.
1206 /* Normalize the 1D Gaussian Kernel
1208 ** NB: a CorrelateNormalize performs a normal Normalize if
1209 ** there are no negative values.
1211 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1212 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1214 /* rotate the 1D kernel by given angle */
1215 RotateKernelInfo(kernel, args->xi );
1220 sigma = fabs(args->sigma),
1223 if ( args->rho < 1.0 )
1224 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1226 kernel->width = (size_t)args->rho;
1227 kernel->x = kernel->y = 0;
1229 kernel->negative_range = kernel->positive_range = 0.0;
1230 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1231 kernel->height*sizeof(*kernel->values));
1232 if (kernel->values == (double *) NULL)
1233 return(DestroyKernelInfo(kernel));
1235 /* A comet blur is half a 1D gaussian curve, so that the object is
1236 ** blurred in one direction only. This may not be quite the right
1237 ** curve to use so may change in the future. The function must be
1238 ** normalised after generation, which also resolves any clipping.
1240 ** As we are normalizing and not subtracting gaussians,
1241 ** there is no need for a divisor in the gaussian formula
1243 ** It is less comples
1245 if ( sigma > MagickEpsilon )
1248 #define KernelRank 3
1249 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1250 (void) ResetMagickMemory(kernel->values,0, (size_t)
1251 kernel->width*sizeof(double));
1252 sigma *= KernelRank; /* simplify the loop expression */
1253 A = 1.0/(2.0*sigma*sigma);
1254 /* B = 1.0/(MagickSQ2PI*sigma); */
1255 for ( u=0; u < v; u++) {
1256 kernel->values[u/KernelRank] +=
1257 exp(-((double)(u*u))*A);
1258 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1260 for (i=0; i < (ssize_t) kernel->width; i++)
1261 kernel->positive_range += kernel->values[i];
1263 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1264 /* B = 1.0/(MagickSQ2PI*sigma); */
1265 for ( i=0; i < (ssize_t) kernel->width; i++)
1266 kernel->positive_range +=
1268 exp(-((double)(i*i))*A);
1269 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1272 else /* special case - generate a unity kernel */
1273 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1274 kernel->width*kernel->height*sizeof(double));
1275 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1276 kernel->positive_range = 1.0;
1279 kernel->minimum = 0.0;
1280 kernel->maximum = kernel->values[0];
1281 kernel->negative_range = 0.0;
1283 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1284 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1289 Convolution Kernels - Well Known Named Constant Kernels
1291 case LaplacianKernel:
1292 { switch ( (int) args->rho ) {
1294 default: /* laplacian square filter -- default */
1295 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1297 case 1: /* laplacian diamond filter */
1298 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1301 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1304 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1306 case 5: /* a 5x5 laplacian */
1307 kernel=ParseKernelArray(
1308 "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");
1310 case 7: /* a 7x7 laplacian */
1311 kernel=ParseKernelArray(
1312 "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" );
1314 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1315 kernel=ParseKernelArray(
1316 "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");
1318 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1319 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1320 kernel=ParseKernelArray(
1321 "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");
1324 if (kernel == (KernelInfo *) NULL)
1326 kernel->type = type;
1330 { /* Simple Sobel Kernel */
1331 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1332 if (kernel == (KernelInfo *) NULL)
1334 kernel->type = type;
1335 RotateKernelInfo(kernel, args->rho);
1340 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1341 if (kernel == (KernelInfo *) NULL)
1343 kernel->type = type;
1344 RotateKernelInfo(kernel, args->rho);
1349 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1350 if (kernel == (KernelInfo *) NULL)
1352 kernel->type = type;
1353 RotateKernelInfo(kernel, args->rho);
1358 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1359 if (kernel == (KernelInfo *) NULL)
1361 kernel->type = type;
1362 RotateKernelInfo(kernel, args->rho);
1367 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1368 if (kernel == (KernelInfo *) NULL)
1370 kernel->type = type;
1371 RotateKernelInfo(kernel, args->rho);
1374 case FreiChenKernel:
1375 /* Direction is set to be left to right positive */
1376 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1377 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1378 { switch ( (int) args->rho ) {
1381 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1382 if (kernel == (KernelInfo *) NULL)
1384 kernel->type = type;
1385 kernel->values[3] = +MagickSQ2;
1386 kernel->values[5] = -MagickSQ2;
1387 CalcKernelMetaData(kernel); /* recalculate meta-data */
1390 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1391 if (kernel == (KernelInfo *) NULL)
1393 kernel->type = type;
1394 kernel->values[1] = kernel->values[3] = +MagickSQ2;
1395 kernel->values[5] = kernel->values[7] = -MagickSQ2;
1396 CalcKernelMetaData(kernel); /* recalculate meta-data */
1397 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1400 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1401 if (kernel == (KernelInfo *) NULL)
1406 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1407 if (kernel == (KernelInfo *) NULL)
1409 kernel->type = type;
1410 kernel->values[3] = +MagickSQ2;
1411 kernel->values[5] = -MagickSQ2;
1412 CalcKernelMetaData(kernel); /* recalculate meta-data */
1413 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1416 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1417 if (kernel == (KernelInfo *) NULL)
1419 kernel->type = type;
1420 kernel->values[1] = +MagickSQ2;
1421 kernel->values[7] = +MagickSQ2;
1422 CalcKernelMetaData(kernel);
1423 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1426 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1427 if (kernel == (KernelInfo *) NULL)
1429 kernel->type = type;
1430 kernel->values[0] = +MagickSQ2;
1431 kernel->values[8] = -MagickSQ2;
1432 CalcKernelMetaData(kernel);
1433 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1436 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1437 if (kernel == (KernelInfo *) NULL)
1439 kernel->type = type;
1440 kernel->values[2] = -MagickSQ2;
1441 kernel->values[6] = +MagickSQ2;
1442 CalcKernelMetaData(kernel);
1443 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1446 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1447 if (kernel == (KernelInfo *) NULL)
1449 kernel->type = type;
1450 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1453 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1454 if (kernel == (KernelInfo *) NULL)
1456 kernel->type = type;
1457 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1460 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1461 if (kernel == (KernelInfo *) NULL)
1463 kernel->type = type;
1464 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1467 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1468 if (kernel == (KernelInfo *) NULL)
1470 kernel->type = type;
1471 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1474 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1475 if (kernel == (KernelInfo *) NULL)
1477 kernel->type = type;
1478 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1481 if ( fabs(args->sigma) > MagickEpsilon )
1482 /* Rotate by correctly supplied 'angle' */
1483 RotateKernelInfo(kernel, args->sigma);
1484 else if ( args->rho > 30.0 || args->rho < -30.0 )
1485 /* Rotate by out of bounds 'type' */
1486 RotateKernelInfo(kernel, args->rho);
1491 Boolean or Shaped Kernels
1495 if (args->rho < 1.0)
1496 kernel->width = kernel->height = 3; /* default radius = 1 */
1498 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1499 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1501 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1502 kernel->height*sizeof(*kernel->values));
1503 if (kernel->values == (double *) NULL)
1504 return(DestroyKernelInfo(kernel));
1506 /* set all kernel values within diamond area to scale given */
1507 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1508 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1509 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1510 kernel->positive_range += kernel->values[i] = args->sigma;
1512 kernel->values[i] = nan;
1513 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1517 case RectangleKernel:
1520 if ( type == SquareKernel )
1522 if (args->rho < 1.0)
1523 kernel->width = kernel->height = 3; /* default radius = 1 */
1525 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1526 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1527 scale = args->sigma;
1530 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1531 if ( args->rho < 1.0 || args->sigma < 1.0 )
1532 return(DestroyKernelInfo(kernel)); /* invalid args given */
1533 kernel->width = (size_t)args->rho;
1534 kernel->height = (size_t)args->sigma;
1535 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1536 args->psi < 0.0 || args->psi > (double)kernel->height )
1537 return(DestroyKernelInfo(kernel)); /* invalid args given */
1538 kernel->x = (ssize_t) args->xi;
1539 kernel->y = (ssize_t) args->psi;
1542 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1543 kernel->height*sizeof(*kernel->values));
1544 if (kernel->values == (double *) NULL)
1545 return(DestroyKernelInfo(kernel));
1547 /* set all kernel values to scale given */
1548 u=(ssize_t) (kernel->width*kernel->height);
1549 for ( i=0; i < u; i++)
1550 kernel->values[i] = scale;
1551 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1552 kernel->positive_range = scale*u;
1557 if (args->rho < 1.0)
1558 kernel->width = kernel->height = 5; /* default radius = 2 */
1560 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1561 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1563 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1564 kernel->height*sizeof(*kernel->values));
1565 if (kernel->values == (double *) NULL)
1566 return(DestroyKernelInfo(kernel));
1568 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1569 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1570 if ( (labs((long) u)+labs((long) v)) <=
1571 ((long)kernel->x + (long)(kernel->x/2)) )
1572 kernel->positive_range += kernel->values[i] = args->sigma;
1574 kernel->values[i] = nan;
1575 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1581 limit = (ssize_t)(args->rho*args->rho);
1583 if (args->rho < 0.4) /* default radius approx 4.3 */
1584 kernel->width = kernel->height = 9L, limit = 18L;
1586 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1587 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1589 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1590 kernel->height*sizeof(*kernel->values));
1591 if (kernel->values == (double *) NULL)
1592 return(DestroyKernelInfo(kernel));
1594 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1595 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1596 if ((u*u+v*v) <= limit)
1597 kernel->positive_range += kernel->values[i] = args->sigma;
1599 kernel->values[i] = nan;
1600 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1605 if (args->rho < 1.0)
1606 kernel->width = kernel->height = 5; /* default radius 2 */
1608 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1609 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1611 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1612 kernel->height*sizeof(*kernel->values));
1613 if (kernel->values == (double *) NULL)
1614 return(DestroyKernelInfo(kernel));
1616 /* set all kernel values along axises to given scale */
1617 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1618 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1619 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1620 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1621 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1626 if (args->rho < 1.0)
1627 kernel->width = kernel->height = 5; /* default radius 2 */
1629 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1630 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1632 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1633 kernel->height*sizeof(*kernel->values));
1634 if (kernel->values == (double *) NULL)
1635 return(DestroyKernelInfo(kernel));
1637 /* set all kernel values along axises to given scale */
1638 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1639 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1640 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1641 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1642 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1656 if (args->rho < args->sigma)
1658 kernel->width = ((size_t)args->sigma)*2+1;
1659 limit1 = (ssize_t)(args->rho*args->rho);
1660 limit2 = (ssize_t)(args->sigma*args->sigma);
1664 kernel->width = ((size_t)args->rho)*2+1;
1665 limit1 = (ssize_t)(args->sigma*args->sigma);
1666 limit2 = (ssize_t)(args->rho*args->rho);
1669 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1671 kernel->height = kernel->width;
1672 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1673 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1674 kernel->height*sizeof(*kernel->values));
1675 if (kernel->values == (double *) NULL)
1676 return(DestroyKernelInfo(kernel));
1678 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1679 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1680 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1681 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1682 { ssize_t radius=u*u+v*v;
1683 if (limit1 < radius && radius <= limit2)
1684 kernel->positive_range += kernel->values[i] = (double) scale;
1686 kernel->values[i] = nan;
1688 kernel->minimum = kernel->maximum = (double) scale;
1689 if ( type == PeaksKernel ) {
1690 /* set the central point in the middle */
1691 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1692 kernel->positive_range = 1.0;
1693 kernel->maximum = 1.0;
1699 kernel=AcquireKernelInfo("ThinSE:482");
1700 if (kernel == (KernelInfo *) NULL)
1702 kernel->type = type;
1703 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1708 kernel=AcquireKernelInfo("ThinSE:87");
1709 if (kernel == (KernelInfo *) NULL)
1711 kernel->type = type;
1712 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1715 case DiagonalsKernel:
1717 switch ( (int) args->rho ) {
1722 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1723 if (kernel == (KernelInfo *) NULL)
1725 kernel->type = type;
1726 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1727 if (new_kernel == (KernelInfo *) NULL)
1728 return(DestroyKernelInfo(kernel));
1729 new_kernel->type = type;
1730 LastKernelInfo(kernel)->next = new_kernel;
1731 ExpandMirrorKernelInfo(kernel);
1735 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1738 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1741 if (kernel == (KernelInfo *) NULL)
1743 kernel->type = type;
1744 RotateKernelInfo(kernel, args->sigma);
1747 case LineEndsKernel:
1748 { /* Kernels for finding the end of thin lines */
1749 switch ( (int) args->rho ) {
1752 /* set of kernels to find all end of lines */
1753 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1755 /* kernel for 4-connected line ends - no rotation */
1756 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1759 /* kernel to add for 8-connected lines - no rotation */
1760 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1763 /* kernel to add for orthogonal line ends - does not find corners */
1764 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1767 /* traditional line end - fails on last T end */
1768 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1771 if (kernel == (KernelInfo *) NULL)
1773 kernel->type = type;
1774 RotateKernelInfo(kernel, args->sigma);
1777 case LineJunctionsKernel:
1778 { /* kernels for finding the junctions of multiple lines */
1779 switch ( (int) args->rho ) {
1782 /* set of kernels to find all line junctions */
1783 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1786 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1789 /* Diagonal T Junctions */
1790 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1793 /* Orthogonal T Junctions */
1794 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1797 /* Diagonal X Junctions */
1798 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1801 /* Orthogonal X Junctions - minimal diamond kernel */
1802 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1805 if (kernel == (KernelInfo *) NULL)
1807 kernel->type = type;
1808 RotateKernelInfo(kernel, args->sigma);
1812 { /* Ridges - Ridge finding kernels */
1815 switch ( (int) args->rho ) {
1818 kernel=ParseKernelArray("3x1:0,1,0");
1819 if (kernel == (KernelInfo *) NULL)
1821 kernel->type = type;
1822 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1825 kernel=ParseKernelArray("4x1:0,1,1,0");
1826 if (kernel == (KernelInfo *) NULL)
1828 kernel->type = type;
1829 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1831 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1832 /* Unfortunatally we can not yet rotate a non-square kernel */
1833 /* But then we can't flip a non-symetrical kernel either */
1834 new_kernel=ParseKernelArray("4x3+1+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+2+1:0,1,1,- -,1,1,- -,1,1,0");
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+1+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("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
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+1: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+2: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+1:-,-,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;
1869 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1870 if (new_kernel == (KernelInfo *) NULL)
1871 return(DestroyKernelInfo(kernel));
1872 new_kernel->type = type;
1873 LastKernelInfo(kernel)->next = new_kernel;
1878 case ConvexHullKernel:
1882 /* first set of 8 kernels */
1883 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1884 if (kernel == (KernelInfo *) NULL)
1886 kernel->type = type;
1887 ExpandRotateKernelInfo(kernel, 90.0);
1888 /* append the mirror versions too - no flip function yet */
1889 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1890 if (new_kernel == (KernelInfo *) NULL)
1891 return(DestroyKernelInfo(kernel));
1892 new_kernel->type = type;
1893 ExpandRotateKernelInfo(new_kernel, 90.0);
1894 LastKernelInfo(kernel)->next = new_kernel;
1897 case SkeletonKernel:
1899 switch ( (int) args->rho ) {
1902 /* Traditional Skeleton...
1903 ** A cyclically rotated single kernel
1905 kernel=AcquireKernelInfo("ThinSE:482");
1906 if (kernel == (KernelInfo *) NULL)
1908 kernel->type = type;
1909 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1912 /* HIPR Variation of the cyclic skeleton
1913 ** Corners of the traditional method made more forgiving,
1914 ** but the retain the same cyclic order.
1916 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1917 if (kernel == (KernelInfo *) NULL)
1919 if (kernel->next == (KernelInfo *) NULL)
1920 return(DestroyKernelInfo(kernel));
1921 kernel->type = type;
1922 kernel->next->type = type;
1923 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1926 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1927 ** "Connectivity-Preserving Morphological Image Thransformations"
1928 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1929 ** http://www.leptonica.com/papers/conn.pdf
1931 kernel=AcquireKernelInfo(
1932 "ThinSE:41; ThinSE:42; ThinSE:43");
1933 if (kernel == (KernelInfo *) NULL)
1935 kernel->type = type;
1936 kernel->next->type = type;
1937 kernel->next->next->type = type;
1938 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1944 { /* Special kernels for general thinning, while preserving connections
1945 ** "Connectivity-Preserving Morphological Image Thransformations"
1946 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1947 ** http://www.leptonica.com/papers/conn.pdf
1949 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
1951 ** Note kernels do not specify the origin pixel, allowing them
1952 ** to be used for both thickening and thinning operations.
1954 switch ( (int) args->rho ) {
1955 /* SE for 4-connected thinning */
1956 case 41: /* SE_4_1 */
1957 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
1959 case 42: /* SE_4_2 */
1960 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
1962 case 43: /* SE_4_3 */
1963 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
1965 case 44: /* SE_4_4 */
1966 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
1968 case 45: /* SE_4_5 */
1969 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
1971 case 46: /* SE_4_6 */
1972 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
1974 case 47: /* SE_4_7 */
1975 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
1977 case 48: /* SE_4_8 */
1978 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
1980 case 49: /* SE_4_9 */
1981 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
1983 /* SE for 8-connected thinning - negatives of the above */
1984 case 81: /* SE_8_0 */
1985 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
1987 case 82: /* SE_8_2 */
1988 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
1990 case 83: /* SE_8_3 */
1991 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
1993 case 84: /* SE_8_4 */
1994 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
1996 case 85: /* SE_8_5 */
1997 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
1999 case 86: /* SE_8_6 */
2000 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2002 case 87: /* SE_8_7 */
2003 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2005 case 88: /* SE_8_8 */
2006 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2008 case 89: /* SE_8_9 */
2009 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2011 /* Special combined SE kernels */
2012 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2013 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2015 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2016 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2018 case 481: /* SE_48_1 - General Connected Corner Kernel */
2019 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2022 case 482: /* SE_48_2 - General Edge Kernel */
2023 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2026 if (kernel == (KernelInfo *) NULL)
2028 kernel->type = type;
2029 RotateKernelInfo(kernel, args->sigma);
2033 Distance Measuring Kernels
2035 case ChebyshevKernel:
2037 if (args->rho < 1.0)
2038 kernel->width = kernel->height = 3; /* default radius = 1 */
2040 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2041 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2043 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2044 kernel->height*sizeof(*kernel->values));
2045 if (kernel->values == (double *) NULL)
2046 return(DestroyKernelInfo(kernel));
2048 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2049 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2050 kernel->positive_range += ( kernel->values[i] =
2051 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2052 kernel->maximum = kernel->values[0];
2055 case ManhattanKernel:
2057 if (args->rho < 1.0)
2058 kernel->width = kernel->height = 3; /* default radius = 1 */
2060 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2061 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2063 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2064 kernel->height*sizeof(*kernel->values));
2065 if (kernel->values == (double *) NULL)
2066 return(DestroyKernelInfo(kernel));
2068 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2069 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2070 kernel->positive_range += ( kernel->values[i] =
2071 args->sigma*(labs((long) u)+labs((long) v)) );
2072 kernel->maximum = kernel->values[0];
2075 case OctagonalKernel:
2077 if (args->rho < 2.0)
2078 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2080 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2081 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2083 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2084 kernel->height*sizeof(*kernel->values));
2085 if (kernel->values == (double *) NULL)
2086 return(DestroyKernelInfo(kernel));
2088 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2089 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2092 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2093 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2094 kernel->positive_range += kernel->values[i] =
2095 args->sigma*MagickMax(r1,r2);
2097 kernel->maximum = kernel->values[0];
2100 case EuclideanKernel:
2102 if (args->rho < 1.0)
2103 kernel->width = kernel->height = 3; /* default radius = 1 */
2105 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2106 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2108 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2109 kernel->height*sizeof(*kernel->values));
2110 if (kernel->values == (double *) NULL)
2111 return(DestroyKernelInfo(kernel));
2113 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2114 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2115 kernel->positive_range += ( kernel->values[i] =
2116 args->sigma*sqrt((double)(u*u+v*v)) );
2117 kernel->maximum = kernel->values[0];
2122 /* No-Op Kernel - Basically just a single pixel on its own */
2123 kernel=ParseKernelArray("1:1");
2124 if (kernel == (KernelInfo *) NULL)
2126 kernel->type = UndefinedKernel;
2135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2139 % C l o n e K e r n e l I n f o %
2143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2145 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2146 % can be modified without effecting the original. The cloned kernel should
2147 % be destroyed using DestroyKernelInfo() when no longer needed.
2149 % The format of the CloneKernelInfo method is:
2151 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2153 % A description of each parameter follows:
2155 % o kernel: the Morphology/Convolution kernel to be cloned
2158 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2166 assert(kernel != (KernelInfo *) NULL);
2167 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2168 if (new_kernel == (KernelInfo *) NULL)
2170 *new_kernel=(*kernel); /* copy values in structure */
2172 /* replace the values with a copy of the values */
2173 new_kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2174 kernel->height*sizeof(*kernel->values));
2175 if (new_kernel->values == (double *) NULL)
2176 return(DestroyKernelInfo(new_kernel));
2177 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2178 new_kernel->values[i]=kernel->values[i];
2180 /* Also clone the next kernel in the kernel list */
2181 if ( kernel->next != (KernelInfo *) NULL ) {
2182 new_kernel->next = CloneKernelInfo(kernel->next);
2183 if ( new_kernel->next == (KernelInfo *) NULL )
2184 return(DestroyKernelInfo(new_kernel));
2191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2195 % D e s t r o y K e r n e l I n f o %
2199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2201 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2204 % The format of the DestroyKernelInfo method is:
2206 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2208 % A description of each parameter follows:
2210 % o kernel: the Morphology/Convolution kernel to be destroyed
2213 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2215 assert(kernel != (KernelInfo *) NULL);
2216 if ( kernel->next != (KernelInfo *) NULL )
2217 kernel->next=DestroyKernelInfo(kernel->next);
2218 kernel->values=(double *) RelinquishAlignedMemory(kernel->values);
2219 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2224 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2228 + E x p a n d M i r r o r K e r n e l I n f o %
2232 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2234 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2235 % sequence of 90-degree rotated kernels but providing a reflected 180
2236 % rotatation, before the -/+ 90-degree rotations.
2238 % This special rotation order produces a better, more symetrical thinning of
2241 % The format of the ExpandMirrorKernelInfo method is:
2243 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2245 % A description of each parameter follows:
2247 % o kernel: the Morphology/Convolution kernel
2249 % This function is only internel to this module, as it is not finalized,
2250 % especially with regard to non-orthogonal angles, and rotation of larger
2255 static void FlopKernelInfo(KernelInfo *kernel)
2256 { /* Do a Flop by reversing each row. */
2264 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2265 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2266 t=k[x], k[x]=k[r], k[r]=t;
2268 kernel->x = kernel->width - kernel->x - 1;
2269 angle = fmod(angle+180.0, 360.0);
2273 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2281 clone = CloneKernelInfo(last);
2282 RotateKernelInfo(clone, 180); /* flip */
2283 LastKernelInfo(last)->next = clone;
2286 clone = CloneKernelInfo(last);
2287 RotateKernelInfo(clone, 90); /* transpose */
2288 LastKernelInfo(last)->next = clone;
2291 clone = CloneKernelInfo(last);
2292 RotateKernelInfo(clone, 180); /* flop */
2293 LastKernelInfo(last)->next = clone;
2299 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2303 + E x p a n d R o t a t e K e r n e l I n f o %
2307 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2309 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2310 % incrementally by the angle given, until the kernel repeats.
2312 % WARNING: 45 degree rotations only works for 3x3 kernels.
2313 % While 90 degree roatations only works for linear and square kernels
2315 % The format of the ExpandRotateKernelInfo method is:
2317 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2319 % A description of each parameter follows:
2321 % o kernel: the Morphology/Convolution kernel
2323 % o angle: angle to rotate in degrees
2325 % This function is only internel to this module, as it is not finalized,
2326 % especially with regard to non-orthogonal angles, and rotation of larger
2330 /* Internal Routine - Return true if two kernels are the same */
2331 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2332 const KernelInfo *kernel2)
2337 /* check size and origin location */
2338 if ( kernel1->width != kernel2->width
2339 || kernel1->height != kernel2->height
2340 || kernel1->x != kernel2->x
2341 || kernel1->y != kernel2->y )
2344 /* check actual kernel values */
2345 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2346 /* Test for Nan equivalence */
2347 if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2349 if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2351 /* Test actual values are equivalent */
2352 if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
2359 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2367 clone = CloneKernelInfo(last);
2368 RotateKernelInfo(clone, angle);
2369 if ( SameKernelInfo(kernel, clone) == MagickTrue )
2371 LastKernelInfo(last)->next = clone;
2374 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2379 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2383 + C a l c M e t a K e r n a l I n f o %
2387 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2389 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2390 % using the kernel values. This should only ne used if it is not possible to
2391 % calculate that meta-data in some easier way.
2393 % It is important that the meta-data is correct before ScaleKernelInfo() is
2394 % used to perform kernel normalization.
2396 % The format of the CalcKernelMetaData method is:
2398 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2400 % A description of each parameter follows:
2402 % o kernel: the Morphology/Convolution kernel to modify
2404 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2405 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2406 % however is not true for flat-shaped morphological kernels.
2408 % WARNING: Only the specific kernel pointed to is modified, not a list of
2411 % This is an internal function and not expected to be useful outside this
2412 % module. This could change however.
2414 static void CalcKernelMetaData(KernelInfo *kernel)
2419 kernel->minimum = kernel->maximum = 0.0;
2420 kernel->negative_range = kernel->positive_range = 0.0;
2421 for (i=0; i < (kernel->width*kernel->height); i++)
2423 if ( fabs(kernel->values[i]) < MagickEpsilon )
2424 kernel->values[i] = 0.0;
2425 ( kernel->values[i] < 0)
2426 ? ( kernel->negative_range += kernel->values[i] )
2427 : ( kernel->positive_range += kernel->values[i] );
2428 Minimize(kernel->minimum, kernel->values[i]);
2429 Maximize(kernel->maximum, kernel->values[i]);
2436 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2440 % M o r p h o l o g y A p p l y %
2444 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2446 % MorphologyApply() applies a morphological method, multiple times using
2447 % a list of multiple kernels.
2449 % It is basically equivalent to as MorphologyImage() (see below) but
2450 % without any user controls. This allows internel programs to use this
2451 % function, to actually perform a specific task without possible interference
2452 % by any API user supplied settings.
2454 % It is MorphologyImage() task to extract any such user controls, and
2455 % pass them to this function for processing.
2457 % More specifically kernels are not normalized/scaled/blended by the
2458 % 'convolve:scale' Image Artifact (setting), nor is the convolve bias
2459 % (-bias setting or image->bias) loooked at, but must be supplied from the
2460 % function arguments.
2462 % The format of the MorphologyApply method is:
2464 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2465 % const ssize_t iterations,const KernelInfo *kernel,
2466 % const CompositeMethod compose,const double bias,
2467 % ExceptionInfo *exception)
2469 % A description of each parameter follows:
2471 % o image: the source image
2473 % o method: the morphology method to be applied.
2475 % o iterations: apply the operation this many times (or no change).
2476 % A value of -1 means loop until no change found.
2477 % How this is applied may depend on the morphology method.
2478 % Typically this is a value of 1.
2480 % o channel: the channel type.
2482 % o kernel: An array of double representing the morphology kernel.
2484 % o compose: How to handle or merge multi-kernel results.
2485 % If 'UndefinedCompositeOp' use default for the Morphology method.
2486 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2487 % Otherwise merge the results using the compose method given.
2489 % o bias: Convolution Output Bias.
2491 % o exception: return any errors or warnings in this structure.
2495 /* Apply a Morphology Primative to an image using the given kernel.
2496 ** Two pre-created images must be provided, and no image is created.
2497 ** It returns the number of pixels that changed between the images
2498 ** for result convergence determination.
2500 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2501 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2502 ExceptionInfo *exception)
2504 #define MorphologyTag "Morphology/Image"
2523 assert(image != (Image *) NULL);
2524 assert(image->signature == MagickSignature);
2525 assert(morphology_image != (Image *) NULL);
2526 assert(morphology_image->signature == MagickSignature);
2527 assert(kernel != (KernelInfo *) NULL);
2528 assert(kernel->signature == MagickSignature);
2529 assert(exception != (ExceptionInfo *) NULL);
2530 assert(exception->signature == MagickSignature);
2536 image_view=AcquireCacheView(image);
2537 morphology_view=AcquireCacheView(morphology_image);
2538 virt_width=image->columns+kernel->width-1;
2540 /* Some methods (including convolve) needs use a reflected kernel.
2541 * Adjust 'origin' offsets to loop though kernel as a reflection.
2546 case ConvolveMorphology:
2547 case DilateMorphology:
2548 case DilateIntensityMorphology:
2549 /*case DistanceMorphology:*/
2550 /* kernel needs to used with reflection about origin */
2551 offx = (ssize_t) kernel->width-offx-1;
2552 offy = (ssize_t) kernel->height-offy-1;
2554 case ErodeMorphology:
2555 case ErodeIntensityMorphology:
2556 case HitAndMissMorphology:
2557 case ThinningMorphology:
2558 case ThickenMorphology:
2559 /* kernel is used as is, without reflection */
2562 assert("Not a Primitive Morphology Method" != (char *) NULL);
2566 if ( method == ConvolveMorphology && kernel->width == 1 )
2567 { /* Special handling (for speed) of vertical (blur) kernels.
2568 ** This performs its handling in columns rather than in rows.
2569 ** This is only done for convolve as it is the only method that
2570 ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
2572 ** Timing tests (on single CPU laptop)
2573 ** Using a vertical 1-d Blue with normal row-by-row (below)
2574 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2576 ** Using this column method
2577 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2580 ** Anthony Thyssen, 14 June 2010
2585 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2586 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2588 for (x=0; x < (ssize_t) image->columns; x++)
2590 register const Quantum
2602 if (status == MagickFalse)
2604 p=GetCacheViewVirtualPixels(image_view,x,-offy,1,image->rows+
2605 kernel->height-1,exception);
2606 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2607 morphology_image->rows,exception);
2608 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2613 /* offset to origin in 'p'. while 'q' points to it directly */
2616 for (y=0; y < (ssize_t) image->rows; y++)
2624 register const double
2627 register const Quantum
2630 /* Copy input image to the output image for unused channels
2631 * This removes need for 'cloning' a new image every iteration
2633 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2634 GetPixelChannels(image)),q);
2635 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2636 GetPixelChannels(image)),q);
2637 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2638 GetPixelChannels(image)),q);
2639 if (image->colorspace == CMYKColorspace)
2640 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2641 GetPixelChannels(image)),q);
2643 /* Set the bias of the weighted average output */
2648 result.black = bias;
2651 /* Weighted Average of pixels using reflected kernel
2653 ** NOTE for correct working of this operation for asymetrical
2654 ** kernels, the kernel needs to be applied in its reflected form.
2655 ** That is its values needs to be reversed.
2657 k = &kernel->values[ kernel->height-1 ];
2659 if ( (image->channel_mask != DefaultChannels) || (image->matte == MagickFalse) )
2660 { /* No 'Sync' involved.
2661 ** Convolution is simple greyscale channel operation
2663 for (v=0; v < (ssize_t) kernel->height; v++) {
2664 if ( IsNan(*k) ) continue;
2665 result.red += (*k)*GetPixelRed(image,k_pixels);
2666 result.green += (*k)*GetPixelGreen(image,k_pixels);
2667 result.blue += (*k)*GetPixelBlue(image,k_pixels);
2668 if (image->colorspace == CMYKColorspace)
2669 result.black+=(*k)*GetPixelBlack(image,k_pixels);
2670 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2672 k_pixels+=GetPixelChannels(image);
2674 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2675 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
2676 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2677 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
2678 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2679 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
2680 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2681 (image->colorspace == CMYKColorspace))
2682 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
2683 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2684 (image->matte == MagickTrue))
2685 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2688 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2689 ** Weight the color channels with Alpha Channel so that
2690 ** transparent pixels are not part of the results.
2693 alpha, /* alpha weighting of colors : kernel*alpha */
2694 gamma; /* divisor, sum of color weighting values */
2697 for (v=0; v < (ssize_t) kernel->height; v++) {
2698 if ( IsNan(*k) ) continue;
2699 alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels));
2701 result.red += alpha*GetPixelRed(image,k_pixels);
2702 result.green += alpha*GetPixelGreen(image,k_pixels);
2703 result.blue += alpha*GetPixelBlue(image,k_pixels);
2704 if (image->colorspace == CMYKColorspace)
2705 result.black += alpha*GetPixelBlack(image,k_pixels);
2706 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2708 k_pixels+=GetPixelChannels(image);
2710 /* Sync'ed channels, all channels are modified */
2711 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2712 SetPixelRed(morphology_image,
2713 ClampToQuantum(gamma*result.red),q);
2714 SetPixelGreen(morphology_image,
2715 ClampToQuantum(gamma*result.green),q);
2716 SetPixelBlue(morphology_image,
2717 ClampToQuantum(gamma*result.blue),q);
2718 if (image->colorspace == CMYKColorspace)
2719 SetPixelBlack(morphology_image,
2720 ClampToQuantum(gamma*result.black),q);
2721 SetPixelAlpha(morphology_image,
2722 ClampToQuantum(result.alpha),q);
2725 /* Count up changed pixels */
2726 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q))
2727 || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q))
2728 || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q))
2729 || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q))
2730 || ((image->colorspace == CMYKColorspace) &&
2731 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
2732 changed++; /* The pixel was changed in some way! */
2733 p+=GetPixelChannels(image);
2734 q+=GetPixelChannels(morphology_image);
2736 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2738 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2743 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2744 #pragma omp critical (MagickCore_MorphologyImage)
2746 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2747 if (proceed == MagickFalse)
2751 morphology_image->type=image->type;
2752 morphology_view=DestroyCacheView(morphology_view);
2753 image_view=DestroyCacheView(image_view);
2754 return(status ? (ssize_t) changed : 0);
2758 ** Normal handling of horizontal or rectangular kernels (row by row)
2760 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2761 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2763 for (y=0; y < (ssize_t) image->rows; y++)
2765 register const Quantum
2777 if (status == MagickFalse)
2779 p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
2780 kernel->height, exception);
2781 q=GetCacheViewAuthenticPixels(morphology_view,0,y,
2782 morphology_image->columns,1,exception);
2783 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2788 /* offset to origin in 'p'. while 'q' points to it directly */
2789 r = virt_width*offy + offx;
2791 for (x=0; x < (ssize_t) image->columns; x++)
2799 register const double
2802 register const Quantum
2810 /* Copy input image to the output image for unused channels
2811 * This removes need for 'cloning' a new image every iteration
2813 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2814 GetPixelChannels(image)),q);
2815 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2816 GetPixelChannels(image)),q);
2817 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2818 GetPixelChannels(image)),q);
2819 if (image->colorspace == CMYKColorspace)
2820 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2821 GetPixelChannels(image)),q);
2828 min.black = (MagickRealType) QuantumRange;
2833 max.black = (MagickRealType) 0;
2834 /* default result is the original pixel value */
2835 result.red = (MagickRealType) GetPixelRed(image,p+r*GetPixelChannels(image));
2836 result.green = (MagickRealType) GetPixelGreen(image,p+r*GetPixelChannels(image));
2837 result.blue = (MagickRealType) GetPixelBlue(image,p+r*GetPixelChannels(image));
2839 if (image->colorspace == CMYKColorspace)
2840 result.black = (MagickRealType) GetPixelBlack(image,p+r*GetPixelChannels(image));
2841 result.alpha=(MagickRealType) GetPixelAlpha(image,p+r*GetPixelChannels(image));
2844 case ConvolveMorphology:
2845 /* Set the bias of the weighted average output */
2850 result.black = bias;
2852 case DilateIntensityMorphology:
2853 case ErodeIntensityMorphology:
2854 /* use a boolean flag indicating when first match found */
2855 result.red = 0.0; /* result is not used otherwise */
2862 case ConvolveMorphology:
2863 /* Weighted Average of pixels using reflected kernel
2865 ** NOTE for correct working of this operation for asymetrical
2866 ** kernels, the kernel needs to be applied in its reflected form.
2867 ** That is its values needs to be reversed.
2869 ** Correlation is actually the same as this but without reflecting
2870 ** the kernel, and thus 'lower-level' that Convolution. However
2871 ** as Convolution is the more common method used, and it does not
2872 ** really cost us much in terms of processing to use a reflected
2873 ** kernel, so it is Convolution that is implemented.
2875 ** Correlation will have its kernel reflected before calling
2876 ** this function to do a Convolve.
2878 ** For more details of Correlation vs Convolution see
2879 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2881 k = &kernel->values[ kernel->width*kernel->height-1 ];
2883 if ( (image->channel_mask != DefaultChannels) ||
2884 (image->matte == MagickFalse) )
2885 { /* No 'Sync' involved.
2886 ** Convolution is simple greyscale channel operation
2888 for (v=0; v < (ssize_t) kernel->height; v++) {
2889 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2890 if ( IsNan(*k) ) continue;
2892 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2893 result.green += (*k)*
2894 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2895 result.blue += (*k)*
2896 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2897 if (image->colorspace == CMYKColorspace)
2898 result.black += (*k)*
2899 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2900 result.alpha += (*k)*
2901 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2903 k_pixels += virt_width*GetPixelChannels(image);
2905 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2906 SetPixelRed(morphology_image,ClampToQuantum(result.red),
2908 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2909 SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2911 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2912 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2914 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2915 (image->colorspace == CMYKColorspace))
2916 SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2918 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2919 (image->matte == MagickTrue))
2920 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
2924 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2925 ** Weight the color channels with Alpha Channel so that
2926 ** transparent pixels are not part of the results.
2929 alpha, /* alpha weighting of colors : kernel*alpha */
2930 gamma; /* divisor, sum of color weighting values */
2933 for (v=0; v < (ssize_t) kernel->height; v++) {
2934 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2935 if ( IsNan(*k) ) continue;
2936 alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels+u*
2937 GetPixelChannels(image)));
2939 result.red += alpha*
2940 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2941 result.green += alpha*
2942 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2943 result.blue += alpha*
2944 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2945 if (image->colorspace == CMYKColorspace)
2946 result.black+=alpha*
2947 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2948 result.alpha += (*k)*
2949 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2951 k_pixels += virt_width*GetPixelChannels(image);
2953 /* Sync'ed channels, all channels are modified */
2954 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2955 SetPixelRed(morphology_image,
2956 ClampToQuantum(gamma*result.red),q);
2957 SetPixelGreen(morphology_image,
2958 ClampToQuantum(gamma*result.green),q);
2959 SetPixelBlue(morphology_image,
2960 ClampToQuantum(gamma*result.blue),q);
2961 if (image->colorspace == CMYKColorspace)
2962 SetPixelBlack(morphology_image,
2963 ClampToQuantum(gamma*result.black),q);
2964 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2968 case ErodeMorphology:
2969 /* Minimum Value within kernel neighbourhood
2971 ** NOTE that the kernel is not reflected for this operation!
2973 ** NOTE: in normal Greyscale Morphology, the kernel value should
2974 ** be added to the real value, this is currently not done, due to
2975 ** the nature of the boolean kernels being used.
2979 for (v=0; v < (ssize_t) kernel->height; v++) {
2980 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2981 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2982 Minimize(min.red, (double)
2983 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
2984 Minimize(min.green, (double)
2985 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
2986 Minimize(min.blue, (double)
2987 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
2988 if (image->colorspace == CMYKColorspace)
2989 Minimize(min.black,(double)
2990 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
2991 Minimize(min.alpha,(double)
2992 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
2994 k_pixels += virt_width*GetPixelChannels(image);
2998 case DilateMorphology:
2999 /* Maximum Value within kernel neighbourhood
3001 ** NOTE for correct working of this operation for asymetrical
3002 ** kernels, the kernel needs to be applied in its reflected form.
3003 ** That is its values needs to be reversed.
3005 ** NOTE: in normal Greyscale Morphology, the kernel value should
3006 ** be added to the real value, this is currently not done, due to
3007 ** the nature of the boolean kernels being used.
3010 k = &kernel->values[ kernel->width*kernel->height-1 ];
3012 for (v=0; v < (ssize_t) kernel->height; v++) {
3013 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3014 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3015 Maximize(max.red, (double)
3016 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3017 Maximize(max.green, (double)
3018 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3019 Maximize(max.blue, (double)
3020 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3021 if (image->colorspace == CMYKColorspace)
3022 Maximize(max.black, (double)
3023 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3024 Maximize(max.alpha,(double)
3025 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3027 k_pixels += virt_width*GetPixelChannels(image);
3031 case HitAndMissMorphology:
3032 case ThinningMorphology:
3033 case ThickenMorphology:
3034 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3036 ** NOTE that the kernel is not reflected for this operation,
3037 ** and consists of both foreground and background pixel
3038 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3039 ** with either Nan or 0.5 values for don't care.
3041 ** Note that this will never produce a meaningless negative
3042 ** result. Such results can cause Thinning/Thicken to not work
3043 ** correctly when used against a greyscale image.
3047 for (v=0; v < (ssize_t) kernel->height; v++) {
3048 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3049 if ( IsNan(*k) ) continue;
3051 { /* minimim of foreground pixels */
3052 Minimize(min.red, (double)
3053 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3054 Minimize(min.green, (double)
3055 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3056 Minimize(min.blue, (double)
3057 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3058 if ( image->colorspace == CMYKColorspace)
3059 Minimize(min.black,(double)
3060 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3061 Minimize(min.alpha,(double)
3062 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3064 else if ( (*k) < 0.3 )
3065 { /* maximum of background pixels */
3066 Maximize(max.red, (double)
3067 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3068 Maximize(max.green, (double)
3069 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3070 Maximize(max.blue, (double)
3071 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3072 if (image->colorspace == CMYKColorspace)
3073 Maximize(max.black, (double)
3074 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3075 Maximize(max.alpha,(double)
3076 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3079 k_pixels += virt_width*GetPixelChannels(image);
3081 /* Pattern Match if difference is positive */
3082 min.red -= max.red; Maximize( min.red, 0.0 );
3083 min.green -= max.green; Maximize( min.green, 0.0 );
3084 min.blue -= max.blue; Maximize( min.blue, 0.0 );
3085 min.black -= max.black; Maximize( min.black, 0.0 );
3086 min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3089 case ErodeIntensityMorphology:
3090 /* Select Pixel with Minimum Intensity within kernel neighbourhood
3092 ** WARNING: the intensity test fails for CMYK and does not
3093 ** take into account the moderating effect of the alpha channel
3094 ** on the intensity.
3096 ** NOTE that the kernel is not reflected for this operation!
3100 for (v=0; v < (ssize_t) kernel->height; v++) {
3101 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3102 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3103 if ( result.red == 0.0 ||
3104 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) < GetPixelIntensity(morphology_image,q) ) {
3105 /* copy the whole pixel - no channel selection */
3106 SetPixelRed(morphology_image,GetPixelRed(image,
3107 k_pixels+u*GetPixelChannels(image)),q);
3108 SetPixelGreen(morphology_image,GetPixelGreen(image,
3109 k_pixels+u*GetPixelChannels(image)),q);
3110 SetPixelBlue(morphology_image,GetPixelBlue(image,
3111 k_pixels+u*GetPixelChannels(image)),q);
3112 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3113 k_pixels+u*GetPixelChannels(image)),q);
3114 if ( result.red > 0.0 ) changed++;
3118 k_pixels += virt_width*GetPixelChannels(image);
3122 case DilateIntensityMorphology:
3123 /* Select Pixel with Maximum Intensity within kernel neighbourhood
3125 ** WARNING: the intensity test fails for CMYK and does not
3126 ** take into account the moderating effect of the alpha channel
3127 ** on the intensity (yet).
3129 ** NOTE for correct working of this operation for asymetrical
3130 ** kernels, the kernel needs to be applied in its reflected form.
3131 ** That is its values needs to be reversed.
3133 k = &kernel->values[ kernel->width*kernel->height-1 ];
3135 for (v=0; v < (ssize_t) kernel->height; v++) {
3136 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3137 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3138 if ( result.red == 0.0 ||
3139 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) > GetPixelIntensity(morphology_image,q) ) {
3140 /* copy the whole pixel - no channel selection */
3141 SetPixelRed(morphology_image,GetPixelRed(image,
3142 k_pixels+u*GetPixelChannels(image)),q);
3143 SetPixelGreen(morphology_image,GetPixelGreen(image,
3144 k_pixels+u*GetPixelChannels(image)),q);
3145 SetPixelBlue(morphology_image,GetPixelBlue(image,
3146 k_pixels+u*GetPixelChannels(image)),q);
3147 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3148 k_pixels+u*GetPixelChannels(image)),q);
3149 if ( result.red > 0.0 ) changed++;
3153 k_pixels += virt_width*GetPixelChannels(image);
3157 This code has been obsoleted by the MorphologyPrimitiveDirect() function.
3158 However it is still (almost) correct coding for Grayscale Morphology.
3161 GrayErode is equivalent but with kernel values subtracted from pixels
3162 without the kernel rotation
3163 GreyDilate is equivalent but using Maximum() instead of Minimum()
3164 using kernel rotation
3166 It has thus been preserved for future implementation of those methods.
3168 case DistanceMorphology:
3169 /* Add kernel Value and select the minimum value found.
3170 ** The result is a iterative distance from edge of image shape.
3172 ** All Distance Kernels are symetrical, but that may not always
3173 ** be the case. For example how about a distance from left edges?
3174 ** To work correctly with asymetrical kernels the reflected kernel
3175 ** needs to be applied.
3177 k = &kernel->values[ kernel->width*kernel->height-1 ];
3179 for (v=0; v < (ssize_t) kernel->height; v++) {
3180 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3181 if ( IsNan(*k) ) continue;
3182 Minimize(result.red, (*k)+k_pixels[u].red);
3183 Minimize(result.green, (*k)+k_pixels[u].green);
3184 Minimize(result.blue, (*k)+k_pixels[u].blue);
3185 Minimize(result.alpha, (*k)+k_pixels[u].alpha);
3186 if ( image->colorspace == CMYKColorspace)
3187 Minimize(result.black,(*k)+GetPixelBlack(p_image,k_indexes+u));
3189 k_pixels += virt_width*GetPixelChannels(image);
3193 case UndefinedMorphology:
3195 break; /* Do nothing */
3197 /* Final mathematics of results (combine with original image?)
3199 ** NOTE: Difference Morphology operators Edge* and *Hat could also
3200 ** be done here but works better with iteration as a image difference
3201 ** in the controling function (below). Thicken and Thinning however
3202 ** should be done here so thay can be iterated correctly.
3205 case HitAndMissMorphology:
3206 case ErodeMorphology:
3207 result = min; /* minimum of neighbourhood */
3209 case DilateMorphology:
3210 result = max; /* maximum of neighbourhood */
3212 case ThinningMorphology:
3213 /* subtract pattern match from original */
3214 result.red -= min.red;
3215 result.green -= min.green;
3216 result.blue -= min.blue;
3217 result.black -= min.black;
3218 result.alpha -= min.alpha;
3220 case ThickenMorphology:
3221 /* Add the pattern matchs to the original */
3222 result.red += min.red;
3223 result.green += min.green;
3224 result.blue += min.blue;
3225 result.black += min.black;
3226 result.alpha += min.alpha;
3229 /* result directly calculated or assigned */
3232 /* Assign the resulting pixel values - Clamping Result */
3234 case UndefinedMorphology:
3235 case ConvolveMorphology:
3236 case DilateIntensityMorphology:
3237 case ErodeIntensityMorphology:
3238 break; /* full pixel was directly assigned - not a channel method */
3240 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3241 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
3242 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3243 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
3244 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3245 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
3246 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3247 (image->colorspace == CMYKColorspace))
3248 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
3249 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3250 (image->matte == MagickTrue))
3251 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3254 /* Count up changed pixels */
3255 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q)) ||
3256 (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q)) ||
3257 (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q)) ||
3258 (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q)) ||
3259 ((image->colorspace == CMYKColorspace) &&
3260 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
3261 changed++; /* The pixel was changed in some way! */
3262 p+=GetPixelChannels(image);
3263 q+=GetPixelChannels(morphology_image);
3265 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3267 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3272 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3273 #pragma omp critical (MagickCore_MorphologyImage)
3275 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3276 if (proceed == MagickFalse)
3280 morphology_view=DestroyCacheView(morphology_view);
3281 image_view=DestroyCacheView(image_view);
3282 return(status ? (ssize_t)changed : -1);
3285 /* This is almost identical to the MorphologyPrimative() function above,
3286 ** but will apply the primitive directly to the image in two passes.
3288 ** That is after each row is 'Sync'ed' into the image, the next row will
3289 ** make use of those values as part of the calculation of the next row.
3290 ** It then repeats, but going in the oppisite (bottom-up) direction.
3292 ** Because of this 'iterative' handling this function can not make use
3293 ** of multi-threaded, parellel processing.
3295 static ssize_t MorphologyPrimitiveDirect(Image *image,
3296 const MorphologyMethod method,const KernelInfo *kernel,
3297 ExceptionInfo *exception)
3320 assert(image != (Image *) NULL);
3321 assert(image->signature == MagickSignature);
3322 assert(kernel != (KernelInfo *) NULL);
3323 assert(kernel->signature == MagickSignature);
3324 assert(exception != (ExceptionInfo *) NULL);
3325 assert(exception->signature == MagickSignature);
3327 /* Some methods (including convolve) needs use a reflected kernel.
3328 * Adjust 'origin' offsets to loop though kernel as a reflection.
3333 case DistanceMorphology:
3334 case VoronoiMorphology:
3335 /* kernel needs to used with reflection about origin */
3336 offx = (ssize_t) kernel->width-offx-1;
3337 offy = (ssize_t) kernel->height-offy-1;
3340 case ?????Morphology:
3341 /* kernel is used as is, without reflection */
3345 assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
3349 /* DO NOT THREAD THIS CODE! */
3350 /* two views into same image (virtual, and actual) */
3351 virt_view=AcquireCacheView(image);
3352 auth_view=AcquireCacheView(image);
3353 virt_width=image->columns+kernel->width-1;
3355 for (y=0; y < (ssize_t) image->rows; y++)
3357 register const Quantum
3369 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3370 ** we read using virtual to get virtual pixel handling, but write back
3371 ** into the same image.
3373 ** Only top half of kernel is processed as we do a single pass downward
3374 ** through the image iterating the distance function as we go.
3376 if (status == MagickFalse)
3378 p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
3380 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3382 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3384 if (status == MagickFalse)
3387 /* offset to origin in 'p'. while 'q' points to it directly */
3388 r = (ssize_t) virt_width*offy + offx;
3390 for (x=0; x < (ssize_t) image->columns; x++)
3398 register const double
3401 register const Quantum
3407 /* Starting Defaults */
3408 GetPixelInfo(image,&result);
3409 SetPixelInfo(image,q,&result);
3410 if ( method != VoronoiMorphology )
3411 result.alpha = QuantumRange - result.alpha;
3414 case DistanceMorphology:
3415 /* Add kernel Value and select the minimum value found. */
3416 k = &kernel->values[ kernel->width*kernel->height-1 ];
3418 for (v=0; v <= (ssize_t) offy; v++) {
3419 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3420 if ( IsNan(*k) ) continue;
3421 Minimize(result.red, (*k)+
3422 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3423 Minimize(result.green, (*k)+
3424 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3425 Minimize(result.blue, (*k)+
3426 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3427 if (image->colorspace == CMYKColorspace)
3428 Minimize(result.black,(*k)+
3429 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3430 Minimize(result.alpha, (*k)+
3431 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3433 k_pixels += virt_width*GetPixelChannels(image);
3435 /* repeat with the just processed pixels of this row */
3436 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3437 k_pixels = q-offx*GetPixelChannels(image);
3438 for (u=0; u < (ssize_t) offx; u++, k--) {
3439 if ( x+u-offx < 0 ) continue; /* off the edge! */
3440 if ( IsNan(*k) ) continue;
3441 Minimize(result.red, (*k)+
3442 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3443 Minimize(result.green, (*k)+
3444 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3445 Minimize(result.blue, (*k)+
3446 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3447 if (image->colorspace == CMYKColorspace)
3448 Minimize(result.black,(*k)+
3449 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3450 Minimize(result.alpha,(*k)+
3451 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3454 case VoronoiMorphology:
3455 /* Apply Distance to 'Matte' channel, coping the closest color.
3457 ** This is experimental, and realy the 'alpha' component should
3458 ** be completely separate 'masking' channel.
3460 k = &kernel->values[ kernel->width*kernel->height-1 ];
3462 for (v=0; v <= (ssize_t) offy; v++) {
3463 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3464 if ( IsNan(*k) ) continue;
3465 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3467 SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3472 k_pixels += virt_width*GetPixelChannels(image);
3474 /* repeat with the just processed pixels of this row */
3475 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3476 k_pixels = q-offx*GetPixelChannels(image);
3477 for (u=0; u < (ssize_t) offx; u++, k--) {
3478 if ( x+u-offx < 0 ) continue; /* off the edge! */
3479 if ( IsNan(*k) ) continue;
3480 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3482 SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3489 /* result directly calculated or assigned */
3492 /* Assign the resulting pixel values - Clamping Result */
3494 case VoronoiMorphology:
3495 SetPixelPixelInfo(image,&result,q);
3498 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3499 SetPixelRed(image,ClampToQuantum(result.red),q);
3500 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3501 SetPixelGreen(image,ClampToQuantum(result.green),q);
3502 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3503 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3504 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3505 (image->colorspace == CMYKColorspace))
3506 SetPixelBlack(image,ClampToQuantum(result.black),q);
3507 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3508 (image->matte == MagickTrue))
3509 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3512 /* Count up changed pixels */
3513 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q)) ||
3514 (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q)) ||
3515 (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q)) ||
3516 (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q)) ||
3517 ((image->colorspace == CMYKColorspace) &&
3518 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3519 changed++; /* The pixel was changed in some way! */
3521 p+=GetPixelChannels(image); /* increment pixel buffers */
3522 q+=GetPixelChannels(image);
3525 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3527 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3528 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3534 /* Do the reversed pass through the image */
3535 for (y=(ssize_t)image->rows-1; y >= 0; y--)
3537 register const Quantum
3549 if (status == MagickFalse)
3551 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3552 ** we read using virtual to get virtual pixel handling, but write back
3553 ** into the same image.
3555 ** Only the bottom half of the kernel will be processes as we
3558 p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
3559 kernel->y+1,exception);
3560 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3562 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3564 if (status == MagickFalse)
3567 /* adjust positions to end of row */
3568 p += (image->columns-1)*GetPixelChannels(image);
3569 q += (image->columns-1)*GetPixelChannels(image);
3571 /* offset to origin in 'p'. while 'q' points to it directly */
3574 for (x=(ssize_t)image->columns-1; x >= 0; x--)
3582 register const double
3585 register const Quantum
3591 /* Default - previously modified pixel */
3592 GetPixelInfo(image,&result);
3593 SetPixelInfo(image,q,&result);
3594 if ( method != VoronoiMorphology )
3595 result.alpha = QuantumRange - result.alpha;
3598 case DistanceMorphology:
3599 /* Add kernel Value and select the minimum value found. */
3600 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3602 for (v=offy; v < (ssize_t) kernel->height; v++) {
3603 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3604 if ( IsNan(*k) ) continue;
3605 Minimize(result.red, (*k)+
3606 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3607 Minimize(result.green, (*k)+
3608 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3609 Minimize(result.blue, (*k)+
3610 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3611 if ( image->colorspace == CMYKColorspace)
3612 Minimize(result.black,(*k)+
3613 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3614 Minimize(result.alpha, (*k)+
3615 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3617 k_pixels += virt_width*GetPixelChannels(image);
3619 /* repeat with the just processed pixels of this row */
3620 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3621 k_pixels = q-offx*GetPixelChannels(image);
3622 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3623 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3624 if ( IsNan(*k) ) continue;
3625 Minimize(result.red, (*k)+
3626 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3627 Minimize(result.green, (*k)+
3628 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3629 Minimize(result.blue, (*k)+
3630 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3631 if ( image->colorspace == CMYKColorspace)
3632 Minimize(result.black, (*k)+
3633 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3634 Minimize(result.alpha, (*k)+
3635 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3638 case VoronoiMorphology:
3639 /* Apply Distance to 'Matte' channel, coping the closest color.
3641 ** This is experimental, and realy the 'alpha' component should
3642 ** be completely separate 'masking' channel.
3644 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3646 for (v=offy; v < (ssize_t) kernel->height; v++) {
3647 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3648 if ( IsNan(*k) ) continue;
3649 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3651 SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3656 k_pixels += virt_width*GetPixelChannels(image);
3658 /* repeat with the just processed pixels of this row */
3659 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3660 k_pixels = q-offx*GetPixelChannels(image);
3661 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3662 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3663 if ( IsNan(*k) ) continue;
3664 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3666 SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3673 /* result directly calculated or assigned */
3676 /* Assign the resulting pixel values - Clamping Result */
3678 case VoronoiMorphology:
3679 SetPixelPixelInfo(image,&result,q);
3682 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3683 SetPixelRed(image,ClampToQuantum(result.red),q);
3684 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3685 SetPixelGreen(image,ClampToQuantum(result.green),q);
3686 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3687 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3688 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3689 (image->colorspace == CMYKColorspace))
3690 SetPixelBlack(image,ClampToQuantum(result.black),q);
3691 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3692 (image->matte == MagickTrue))
3693 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3696 /* Count up changed pixels */
3697 if ( (GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q))
3698 || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q))
3699 || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q))
3700 || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q))
3701 || ((image->colorspace == CMYKColorspace) &&
3702 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3703 changed++; /* The pixel was changed in some way! */
3705 p-=GetPixelChannels(image); /* go backward through pixel buffers */
3706 q-=GetPixelChannels(image);
3708 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3710 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3711 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3717 auth_view=DestroyCacheView(auth_view);
3718 virt_view=DestroyCacheView(virt_view);
3719 return(status ? (ssize_t) changed : -1);
3722 /* Apply a Morphology by calling theabove low level primitive application
3723 ** functions. This function handles any iteration loops, composition or
3724 ** re-iteration of results, and compound morphology methods that is based
3725 ** on multiple low-level (staged) morphology methods.
3727 ** Basically this provides the complex grue between the requested morphology
3728 ** method and raw low-level implementation (above).
3730 MagickPrivate Image *MorphologyApply(const Image *image,
3731 const MorphologyMethod method, const ssize_t iterations,
3732 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3733 ExceptionInfo *exception)
3739 *curr_image, /* Image we are working with or iterating */
3740 *work_image, /* secondary image for primitive iteration */
3741 *save_image, /* saved image - for 'edge' method only */
3742 *rslt_image; /* resultant image - after multi-kernel handling */
3745 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3746 *norm_kernel, /* the current normal un-reflected kernel */
3747 *rflt_kernel, /* the current reflected kernel (if needed) */
3748 *this_kernel; /* the kernel being applied */
3751 primitive; /* the current morphology primitive being applied */
3754 rslt_compose; /* multi-kernel compose method for results to use */
3757 special, /* do we use a direct modify function? */
3758 verbose; /* verbose output of results */
3761 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3762 method_limit, /* maximum number of compound method iterations */
3763 kernel_number, /* Loop 2: the kernel number being applied */
3764 stage_loop, /* Loop 3: primitive loop for compound morphology */
3765 stage_limit, /* how many primitives are in this compound */
3766 kernel_loop, /* Loop 4: iterate the kernel over image */
3767 kernel_limit, /* number of times to iterate kernel */
3768 count, /* total count of primitive steps applied */
3769 kernel_changed, /* total count of changed using iterated kernel */
3770 method_changed; /* total count of changed over method iteration */
3773 changed; /* number pixels changed by last primitive operation */
3778 assert(image != (Image *) NULL);
3779 assert(image->signature == MagickSignature);
3780 assert(kernel != (KernelInfo *) NULL);
3781 assert(kernel->signature == MagickSignature);
3782 assert(exception != (ExceptionInfo *) NULL);
3783 assert(exception->signature == MagickSignature);
3785 count = 0; /* number of low-level morphology primitives performed */
3786 if ( iterations == 0 )
3787 return((Image *)NULL); /* null operation - nothing to do! */
3789 kernel_limit = (size_t) iterations;
3790 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3791 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3793 verbose = IsMagickTrue(GetImageArtifact(image,"verbose"));
3795 /* initialise for cleanup */
3796 curr_image = (Image *) image;
3797 curr_compose = image->compose;
3798 (void) curr_compose;
3799 work_image = save_image = rslt_image = (Image *) NULL;
3800 reflected_kernel = (KernelInfo *) NULL;
3802 /* Initialize specific methods
3803 * + which loop should use the given iteratations
3804 * + how many primitives make up the compound morphology
3805 * + multi-kernel compose method to use (by default)
3807 method_limit = 1; /* just do method once, unless otherwise set */
3808 stage_limit = 1; /* assume method is not a compound */
3809 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3810 rslt_compose = compose; /* and we are composing multi-kernels as given */
3812 case SmoothMorphology: /* 4 primitive compound morphology */
3815 case OpenMorphology: /* 2 primitive compound morphology */
3816 case OpenIntensityMorphology:
3817 case TopHatMorphology:
3818 case CloseMorphology:
3819 case CloseIntensityMorphology:
3820 case BottomHatMorphology:
3821 case EdgeMorphology:
3824 case HitAndMissMorphology:
3825 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3827 case ThinningMorphology:
3828 case ThickenMorphology:
3829 method_limit = kernel_limit; /* iterate the whole method */
3830 kernel_limit = 1; /* do not do kernel iteration */
3832 case DistanceMorphology:
3833 case VoronoiMorphology:
3834 special = MagickTrue;
3840 /* Apply special methods with special requirments
3841 ** For example, single run only, or post-processing requirements
3843 if ( special == MagickTrue )
3845 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3846 if (rslt_image == (Image *) NULL)
3848 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3851 changed = MorphologyPrimitiveDirect(rslt_image, method,
3854 if ( verbose == MagickTrue )
3855 (void) (void) FormatLocaleFile(stderr,
3856 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3857 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3858 1.0,0.0,1.0, (double) changed);
3863 if ( method == VoronoiMorphology ) {
3864 /* Preserve the alpha channel of input image - but turned off */
3865 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3867 (void) CompositeImage(rslt_image, CopyOpacityCompositeOp, image, 0, 0,
3869 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3875 /* Handle user (caller) specified multi-kernel composition method */
3876 if ( compose != UndefinedCompositeOp )
3877 rslt_compose = compose; /* override default composition for method */
3878 if ( rslt_compose == UndefinedCompositeOp )
3879 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3881 /* Some methods require a reflected kernel to use with primitives.
3882 * Create the reflected kernel for those methods. */
3884 case CorrelateMorphology:
3885 case CloseMorphology:
3886 case CloseIntensityMorphology:
3887 case BottomHatMorphology:
3888 case SmoothMorphology:
3889 reflected_kernel = CloneKernelInfo(kernel);
3890 if (reflected_kernel == (KernelInfo *) NULL)
3892 RotateKernelInfo(reflected_kernel,180);
3898 /* Loops around more primitive morpholgy methods
3899 ** erose, dilate, open, close, smooth, edge, etc...
3901 /* Loop 1: iterate the compound method */
3904 while ( method_loop < method_limit && method_changed > 0 ) {
3908 /* Loop 2: iterate over each kernel in a multi-kernel list */
3909 norm_kernel = (KernelInfo *) kernel;
3910 this_kernel = (KernelInfo *) kernel;
3911 rflt_kernel = reflected_kernel;
3914 while ( norm_kernel != NULL ) {
3916 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3917 stage_loop = 0; /* the compound morphology stage number */
3918 while ( stage_loop < stage_limit ) {
3919 stage_loop++; /* The stage of the compound morphology */
3921 /* Select primitive morphology for this stage of compound method */
3922 this_kernel = norm_kernel; /* default use unreflected kernel */
3923 primitive = method; /* Assume method is a primitive */
3925 case ErodeMorphology: /* just erode */
3926 case EdgeInMorphology: /* erode and image difference */
3927 primitive = ErodeMorphology;
3929 case DilateMorphology: /* just dilate */
3930 case EdgeOutMorphology: /* dilate and image difference */
3931 primitive = DilateMorphology;
3933 case OpenMorphology: /* erode then dialate */
3934 case TopHatMorphology: /* open and image difference */
3935 primitive = ErodeMorphology;
3936 if ( stage_loop == 2 )
3937 primitive = DilateMorphology;
3939 case OpenIntensityMorphology:
3940 primitive = ErodeIntensityMorphology;
3941 if ( stage_loop == 2 )
3942 primitive = DilateIntensityMorphology;
3944 case CloseMorphology: /* dilate, then erode */
3945 case BottomHatMorphology: /* close and image difference */
3946 this_kernel = rflt_kernel; /* use the reflected kernel */
3947 primitive = DilateMorphology;
3948 if ( stage_loop == 2 )
3949 primitive = ErodeMorphology;
3951 case CloseIntensityMorphology:
3952 this_kernel = rflt_kernel; /* use the reflected kernel */
3953 primitive = DilateIntensityMorphology;
3954 if ( stage_loop == 2 )
3955 primitive = ErodeIntensityMorphology;
3957 case SmoothMorphology: /* open, close */
3958 switch ( stage_loop ) {
3959 case 1: /* start an open method, which starts with Erode */
3960 primitive = ErodeMorphology;
3962 case 2: /* now Dilate the Erode */
3963 primitive = DilateMorphology;
3965 case 3: /* Reflect kernel a close */
3966 this_kernel = rflt_kernel; /* use the reflected kernel */
3967 primitive = DilateMorphology;
3969 case 4: /* Finish the Close */
3970 this_kernel = rflt_kernel; /* use the reflected kernel */
3971 primitive = ErodeMorphology;
3975 case EdgeMorphology: /* dilate and erode difference */
3976 primitive = DilateMorphology;
3977 if ( stage_loop == 2 ) {
3978 save_image = curr_image; /* save the image difference */
3979 curr_image = (Image *) image;
3980 primitive = ErodeMorphology;
3983 case CorrelateMorphology:
3984 /* A Correlation is a Convolution with a reflected kernel.
3985 ** However a Convolution is a weighted sum using a reflected
3986 ** kernel. It may seem stange to convert a Correlation into a
3987 ** Convolution as the Correlation is the simplier method, but
3988 ** Convolution is much more commonly used, and it makes sense to
3989 ** implement it directly so as to avoid the need to duplicate the
3990 ** kernel when it is not required (which is typically the
3993 this_kernel = rflt_kernel; /* use the reflected kernel */
3994 primitive = ConvolveMorphology;
3999 assert( this_kernel != (KernelInfo *) NULL );
4001 /* Extra information for debugging compound operations */
4002 if ( verbose == MagickTrue ) {
4003 if ( stage_limit > 1 )
4004 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
4005 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
4006 method_loop,(double) stage_loop);
4007 else if ( primitive != method )
4008 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
4009 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
4015 /* Loop 4: Iterate the kernel with primitive */
4019 while ( kernel_loop < kernel_limit && changed > 0 ) {
4020 kernel_loop++; /* the iteration of this kernel */
4022 /* Create a clone as the destination image, if not yet defined */
4023 if ( work_image == (Image *) NULL )
4025 work_image=CloneImage(image,0,0,MagickTrue,exception);
4026 if (work_image == (Image *) NULL)
4028 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
4030 /* work_image->type=image->type; ??? */
4033 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4035 changed = MorphologyPrimitive(curr_image, work_image, primitive,
4036 this_kernel, bias, exception);
4038 if ( verbose == MagickTrue ) {
4039 if ( kernel_loop > 1 )
4040 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4041 (void) (void) FormatLocaleFile(stderr,
4042 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4043 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4044 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4045 (double) (method_loop+kernel_loop-1),(double) kernel_number,
4046 (double) count,(double) changed);
4050 kernel_changed += changed;
4051 method_changed += changed;
4053 /* prepare next loop */
4054 { Image *tmp = work_image; /* swap images for iteration */
4055 work_image = curr_image;
4058 if ( work_image == image )
4059 work_image = (Image *) NULL; /* replace input 'image' */
4061 } /* End Loop 4: Iterate the kernel with primitive */
4063 if ( verbose == MagickTrue && kernel_changed != (size_t)changed )
4064 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
4065 if ( verbose == MagickTrue && stage_loop < stage_limit )
4066 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4069 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4070 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
4071 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
4072 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
4073 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
4076 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4078 /* Final Post-processing for some Compound Methods
4080 ** The removal of any 'Sync' channel flag in the Image Compositon
4081 ** below ensures the methematical compose method is applied in a
4082 ** purely mathematical way, and only to the selected channels.
4083 ** Turn off SVG composition 'alpha blending'.
4086 case EdgeOutMorphology:
4087 case EdgeInMorphology:
4088 case TopHatMorphology:
4089 case BottomHatMorphology:
4090 if ( verbose == MagickTrue )
4091 (void) FormatLocaleFile(stderr,
4092 "\n%s: Difference with original image",CommandOptionToMnemonic(
4093 MagickMorphologyOptions, method) );
4094 (void) CompositeImage(curr_image,DifferenceCompositeOp,image,0,0,
4097 case EdgeMorphology:
4098 if ( verbose == MagickTrue )
4099 (void) FormatLocaleFile(stderr,
4100 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
4101 MagickMorphologyOptions, method) );
4102 (void) CompositeImage(curr_image,DifferenceCompositeOp,save_image,0,
4104 save_image = DestroyImage(save_image); /* finished with save image */
4110 /* multi-kernel handling: re-iterate, or compose results */
4111 if ( kernel->next == (KernelInfo *) NULL )
4112 rslt_image = curr_image; /* just return the resulting image */
4113 else if ( rslt_compose == NoCompositeOp )
4114 { if ( verbose == MagickTrue ) {
4115 if ( this_kernel->next != (KernelInfo *) NULL )
4116 (void) FormatLocaleFile(stderr, " (re-iterate)");
4118 (void) FormatLocaleFile(stderr, " (done)");
4120 rslt_image = curr_image; /* return result, and re-iterate */
4122 else if ( rslt_image == (Image *) NULL)
4123 { if ( verbose == MagickTrue )
4124 (void) FormatLocaleFile(stderr, " (save for compose)");
4125 rslt_image = curr_image;
4126 curr_image = (Image *) image; /* continue with original image */
4129 { /* Add the new 'current' result to the composition
4131 ** The removal of any 'Sync' channel flag in the Image Compositon
4132 ** below ensures the methematical compose method is applied in a
4133 ** purely mathematical way, and only to the selected channels.
4134 ** IE: Turn off SVG composition 'alpha blending'.
4136 if ( verbose == MagickTrue )
4137 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4138 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4139 (void) CompositeImage(rslt_image, rslt_compose, curr_image, 0, 0,
4141 curr_image = DestroyImage(curr_image);
4142 curr_image = (Image *) image; /* continue with original image */
4144 if ( verbose == MagickTrue )
4145 (void) FormatLocaleFile(stderr, "\n");
4147 /* loop to the next kernel in a multi-kernel list */
4148 norm_kernel = norm_kernel->next;
4149 if ( rflt_kernel != (KernelInfo *) NULL )
4150 rflt_kernel = rflt_kernel->next;
4152 } /* End Loop 2: Loop over each kernel */
4154 } /* End Loop 1: compound method interation */
4158 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4160 if ( curr_image == rslt_image )
4161 curr_image = (Image *) NULL;
4162 if ( rslt_image != (Image *) NULL )
4163 rslt_image = DestroyImage(rslt_image);
4165 if ( curr_image == rslt_image || curr_image == image )
4166 curr_image = (Image *) NULL;
4167 if ( curr_image != (Image *) NULL )
4168 curr_image = DestroyImage(curr_image);
4169 if ( work_image != (Image *) NULL )
4170 work_image = DestroyImage(work_image);
4171 if ( save_image != (Image *) NULL )
4172 save_image = DestroyImage(save_image);
4173 if ( reflected_kernel != (KernelInfo *) NULL )
4174 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4180 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4184 % M o r p h o l o g y I m a g e %
4188 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4190 % MorphologyImage() applies a user supplied kernel to the image
4191 % according to the given mophology method.
4193 % This function applies any and all user defined settings before calling
4194 % the above internal function MorphologyApply().
4196 % User defined settings include...
4197 % * Output Bias for Convolution and correlation ("-bias")
4198 % * Kernel Scale/normalize settings ("-set 'option:convolve:scale'")
4199 % This can also includes the addition of a scaled unity kernel.
4200 % * Show Kernel being applied ("-set option:showkernel 1")
4202 % The format of the MorphologyImage method is:
4204 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4205 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4207 % Image *MorphologyImage(const Image *image, const ChannelType
4208 % channel,MorphologyMethod method,const ssize_t iterations,
4209 % KernelInfo *kernel,ExceptionInfo *exception)
4211 % A description of each parameter follows:
4213 % o image: the image.
4215 % o method: the morphology method to be applied.
4217 % o iterations: apply the operation this many times (or no change).
4218 % A value of -1 means loop until no change found.
4219 % How this is applied may depend on the morphology method.
4220 % Typically this is a value of 1.
4222 % o kernel: An array of double representing the morphology kernel.
4223 % Warning: kernel may be normalized for the Convolve method.
4225 % o exception: return any errors or warnings in this structure.
4228 MagickExport Image *MorphologyImage(const Image *image,
4229 const MorphologyMethod method,const ssize_t iterations,
4230 const KernelInfo *kernel,ExceptionInfo *exception)
4242 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4243 * This is done BEFORE the ShowKernelInfo() function is called so that
4244 * users can see the results of the 'option:convolve:scale' option.
4246 curr_kernel = (KernelInfo *) kernel;
4247 if ( method == ConvolveMorphology || method == CorrelateMorphology )
4251 artifact = GetImageArtifact(image,"convolve:scale");
4252 if ( artifact != (const char *)NULL ) {
4253 if ( curr_kernel == kernel )
4254 curr_kernel = CloneKernelInfo(kernel);
4255 if (curr_kernel == (KernelInfo *) NULL) {
4256 curr_kernel=DestroyKernelInfo(curr_kernel);
4257 return((Image *) NULL);
4259 ScaleGeometryKernelInfo(curr_kernel, artifact);
4263 /* display the (normalized) kernel via stderr */
4264 if ( IsMagickTrue(GetImageArtifact(image,"showkernel"))
4265 || IsMagickTrue(GetImageArtifact(image,"convolve:showkernel"))
4266 || IsMagickTrue(GetImageArtifact(image,"morphology:showkernel")) )
4267 ShowKernelInfo(curr_kernel);
4269 /* Override the default handling of multi-kernel morphology results
4270 * If 'Undefined' use the default method
4271 * If 'None' (default for 'Convolve') re-iterate previous result
4272 * Otherwise merge resulting images using compose method given.
4273 * Default for 'HitAndMiss' is 'Lighten'.
4277 artifact = GetImageArtifact(image,"morphology:compose");
4278 compose = UndefinedCompositeOp; /* use default for method */
4279 if ( artifact != (const char *) NULL)
4280 compose=(CompositeOperator) ParseCommandOption(MagickComposeOptions,
4281 MagickFalse,artifact);
4283 /* Apply the Morphology */
4284 morphology_image = MorphologyApply(image, method, iterations,
4285 curr_kernel, compose, image->bias, exception);
4287 /* Cleanup and Exit */
4288 if ( curr_kernel != kernel )
4289 curr_kernel=DestroyKernelInfo(curr_kernel);
4290 return(morphology_image);
4294 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4298 + R o t a t e K e r n e l I n f o %
4302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4304 % RotateKernelInfo() rotates the kernel by the angle given.
4306 % Currently it is restricted to 90 degree angles, of either 1D kernels
4307 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4308 % It will ignore usless rotations for specific 'named' built-in kernels.
4310 % The format of the RotateKernelInfo method is:
4312 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4314 % A description of each parameter follows:
4316 % o kernel: the Morphology/Convolution kernel
4318 % o angle: angle to rotate in degrees
4320 % This function is currently internal to this module only, but can be exported
4321 % to other modules if needed.
4323 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4325 /* angle the lower kernels first */
4326 if ( kernel->next != (KernelInfo *) NULL)
4327 RotateKernelInfo(kernel->next, angle);
4329 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4331 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4334 /* Modulus the angle */
4335 angle = fmod(angle, 360.0);
4339 if ( 337.5 < angle || angle <= 22.5 )
4340 return; /* Near zero angle - no change! - At least not at this time */
4342 /* Handle special cases */
4343 switch (kernel->type) {
4344 /* These built-in kernels are cylindrical kernels, rotating is useless */
4345 case GaussianKernel:
4350 case LaplacianKernel:
4351 case ChebyshevKernel:
4352 case ManhattanKernel:
4353 case EuclideanKernel:
4356 /* These may be rotatable at non-90 angles in the future */
4357 /* but simply rotating them in multiples of 90 degrees is useless */
4364 /* These only allows a +/-90 degree rotation (by transpose) */
4365 /* A 180 degree rotation is useless */
4367 if ( 135.0 < angle && angle <= 225.0 )
4369 if ( 225.0 < angle && angle <= 315.0 )
4376 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4377 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4379 if ( kernel->width == 3 && kernel->height == 3 )
4380 { /* Rotate a 3x3 square by 45 degree angle */
4381 MagickRealType t = kernel->values[0];
4382 kernel->values[0] = kernel->values[3];
4383 kernel->values[3] = kernel->values[6];
4384 kernel->values[6] = kernel->values[7];
4385 kernel->values[7] = kernel->values[8];
4386 kernel->values[8] = kernel->values[5];
4387 kernel->values[5] = kernel->values[2];
4388 kernel->values[2] = kernel->values[1];
4389 kernel->values[1] = t;
4390 /* rotate non-centered origin */
4391 if ( kernel->x != 1 || kernel->y != 1 ) {
4393 x = (ssize_t) kernel->x-1;
4394 y = (ssize_t) kernel->y-1;
4395 if ( x == y ) x = 0;
4396 else if ( x == 0 ) x = -y;
4397 else if ( x == -y ) y = 0;
4398 else if ( y == 0 ) y = x;
4399 kernel->x = (ssize_t) x+1;
4400 kernel->y = (ssize_t) y+1;
4402 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4403 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4406 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4408 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4410 if ( kernel->width == 1 || kernel->height == 1 )
4411 { /* Do a transpose of a 1 dimensional kernel,
4412 ** which results in a fast 90 degree rotation of some type.
4416 t = (ssize_t) kernel->width;
4417 kernel->width = kernel->height;
4418 kernel->height = (size_t) t;
4420 kernel->x = kernel->y;
4422 if ( kernel->width == 1 ) {
4423 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4424 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4426 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4427 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4430 else if ( kernel->width == kernel->height )
4431 { /* Rotate a square array of values by 90 degrees */
4434 register MagickRealType
4437 for( i=0, x=kernel->width-1; i<=x; i++, x--)
4438 for( j=0, y=kernel->height-1; j<y; j++, y--)
4439 { t = k[i+j*kernel->width];
4440 k[i+j*kernel->width] = k[j+x*kernel->width];
4441 k[j+x*kernel->width] = k[x+y*kernel->width];
4442 k[x+y*kernel->width] = k[y+i*kernel->width];
4443 k[y+i*kernel->width] = t;
4446 /* rotate the origin - relative to center of array */
4447 { register ssize_t x,y;
4448 x = (ssize_t) (kernel->x*2-kernel->width+1);
4449 y = (ssize_t) (kernel->y*2-kernel->height+1);
4450 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4451 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4453 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4454 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4457 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4459 if ( 135.0 < angle && angle <= 225.0 )
4461 /* For a 180 degree rotation - also know as a reflection
4462 * This is actually a very very common operation!
4463 * Basically all that is needed is a reversal of the kernel data!
4464 * And a reflection of the origon
4472 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
4473 t=k[i], k[i]=k[j], k[j]=t;
4475 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4476 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4477 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4478 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4480 /* At this point angle should at least between -45 (315) and +45 degrees
4481 * In the future some form of non-orthogonal angled rotates could be
4482 * performed here, posibily with a linear kernel restriction.
4489 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4493 % S c a l e G e o m e t r y K e r n e l I n f o %
4497 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4499 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4500 % provided as a "-set option:convolve:scale {geometry}" user setting,
4501 % and modifies the kernel according to the parsed arguments of that setting.
4503 % The first argument (and any normalization flags) are passed to
4504 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4505 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4506 % into the scaled/normalized kernel.
4508 % The format of the ScaleGeometryKernelInfo method is:
4510 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4511 % const double scaling_factor,const MagickStatusType normalize_flags)
4513 % A description of each parameter follows:
4515 % o kernel: the Morphology/Convolution kernel to modify
4518 % The geometry string to parse, typically from the user provided
4519 % "-set option:convolve:scale {geometry}" setting.
4522 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4523 const char *geometry)
4530 SetGeometryInfo(&args);
4531 flags = (GeometryFlags) ParseGeometry(geometry, &args);
4534 /* For Debugging Geometry Input */
4535 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4536 flags, args.rho, args.sigma, args.xi, args.psi );
4539 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4540 args.rho *= 0.01, args.sigma *= 0.01;
4542 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4544 if ( (flags & SigmaValue) == 0 )
4547 /* Scale/Normalize the input kernel */
4548 ScaleKernelInfo(kernel, args.rho, flags);
4550 /* Add Unity Kernel, for blending with original */
4551 if ( (flags & SigmaValue) != 0 )
4552 UnityAddKernelInfo(kernel, args.sigma);
4557 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4561 % S c a l e K e r n e l I n f o %
4565 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4567 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4568 % without normalization of the sum of the kernel values (as per given flags).
4570 % By default (no flags given) the values within the kernel is scaled
4571 % directly using given scaling factor without change.
4573 % If either of the two 'normalize_flags' are given the kernel will first be
4574 % normalized and then further scaled by the scaling factor value given.
4576 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4577 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4578 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4579 % non-HDRI versions of IM this may cause images to have any negative results
4580 % clipped, unless some 'bias' is used.
4582 % More specifically. Kernels which only contain positive values (such as a
4583 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4584 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4586 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4587 % the kernel will be scaled by the absolute of the sum of kernel values, so
4588 % that it will generally fall within the +/- 1.0 range.
4590 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4591 % will be scaled by just the sum of the postive values, so that its output
4592 % range will again fall into the +/- 1.0 range.
4594 % For special kernels designed for locating shapes using 'Correlate', (often
4595 % only containing +1 and -1 values, representing foreground/brackground
4596 % matching) a special normalization method is provided to scale the positive
4597 % values separately to those of the negative values, so the kernel will be
4598 % forced to become a zero-sum kernel better suited to such searches.
4600 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4601 % attributes within the kernel structure have been correctly set during the
4604 % NOTE: The values used for 'normalize_flags' have been selected specifically
4605 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4606 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4608 % The format of the ScaleKernelInfo method is:
4610 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4611 % const MagickStatusType normalize_flags )
4613 % A description of each parameter follows:
4615 % o kernel: the Morphology/Convolution kernel
4618 % multiply all values (after normalization) by this factor if not
4619 % zero. If the kernel is normalized regardless of any flags.
4621 % o normalize_flags:
4622 % GeometryFlags defining normalization method to use.
4623 % specifically: NormalizeValue, CorrelateNormalizeValue,
4624 % and/or PercentValue
4627 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4628 const double scaling_factor,const GeometryFlags normalize_flags)
4637 /* do the other kernels in a multi-kernel list first */
4638 if ( kernel->next != (KernelInfo *) NULL)
4639 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4641 /* Normalization of Kernel */
4643 if ( (normalize_flags&NormalizeValue) != 0 ) {
4644 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4645 /* non-zero-summing kernel (generally positive) */
4646 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4648 /* zero-summing kernel */
4649 pos_scale = kernel->positive_range;
4651 /* Force kernel into a normalized zero-summing kernel */
4652 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4653 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4654 ? kernel->positive_range : 1.0;
4655 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
4656 ? -kernel->negative_range : 1.0;
4659 neg_scale = pos_scale;
4661 /* finialize scaling_factor for positive and negative components */
4662 pos_scale = scaling_factor/pos_scale;
4663 neg_scale = scaling_factor/neg_scale;
4665 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4666 if ( ! IsNan(kernel->values[i]) )
4667 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4669 /* convolution output range */
4670 kernel->positive_range *= pos_scale;
4671 kernel->negative_range *= neg_scale;
4672 /* maximum and minimum values in kernel */
4673 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4674 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4676 /* swap kernel settings if user's scaling factor is negative */
4677 if ( scaling_factor < MagickEpsilon ) {
4679 t = kernel->positive_range;
4680 kernel->positive_range = kernel->negative_range;
4681 kernel->negative_range = t;
4682 t = kernel->maximum;
4683 kernel->maximum = kernel->minimum;
4684 kernel->minimum = 1;
4691 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4695 % S h o w K e r n e l I n f o %
4699 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4701 % ShowKernelInfo() outputs the details of the given kernel defination to
4702 % standard error, generally due to a users 'showkernel' option request.
4704 % The format of the ShowKernel method is:
4706 % void ShowKernelInfo(const KernelInfo *kernel)
4708 % A description of each parameter follows:
4710 % o kernel: the Morphology/Convolution kernel
4713 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4721 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4723 (void) FormatLocaleFile(stderr, "Kernel");
4724 if ( kernel->next != (KernelInfo *) NULL )
4725 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4726 (void) FormatLocaleFile(stderr, " \"%s",
4727 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4728 if ( fabs(k->angle) > MagickEpsilon )
4729 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4730 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4731 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4732 (void) FormatLocaleFile(stderr,
4733 " with values from %.*lg to %.*lg\n",
4734 GetMagickPrecision(), k->minimum,
4735 GetMagickPrecision(), k->maximum);
4736 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4737 GetMagickPrecision(), k->negative_range,
4738 GetMagickPrecision(), k->positive_range);
4739 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4740 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4741 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4742 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4744 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4745 GetMagickPrecision(), k->positive_range+k->negative_range);
4746 for (i=v=0; v < k->height; v++) {
4747 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4748 for (u=0; u < k->width; u++, i++)
4749 if ( IsNan(k->values[i]) )
4750 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4752 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4753 GetMagickPrecision(), k->values[i]);
4754 (void) FormatLocaleFile(stderr,"\n");
4760 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4764 % U n i t y A d d K e r n a l I n f o %
4768 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4770 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4771 % to the given pre-scaled and normalized Kernel. This in effect adds that
4772 % amount of the original image into the resulting convolution kernel. This
4773 % value is usually provided by the user as a percentage value in the
4774 % 'convolve:scale' setting.
4776 % The resulting effect is to convert the defined kernels into blended
4777 % soft-blurs, unsharp kernels or into sharpening kernels.
4779 % The format of the UnityAdditionKernelInfo method is:
4781 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4783 % A description of each parameter follows:
4785 % o kernel: the Morphology/Convolution kernel
4788 % scaling factor for the unity kernel to be added to
4792 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4795 /* do the other kernels in a multi-kernel list first */
4796 if ( kernel->next != (KernelInfo *) NULL)
4797 UnityAddKernelInfo(kernel->next, scale);
4799 /* Add the scaled unity kernel to the existing kernel */
4800 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4801 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4807 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4811 % Z e r o K e r n e l N a n s %
4815 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4817 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4818 % the kernel with a zero value. This is typically done when the kernel will
4819 % be used in special hardware (GPU) convolution processors, to simply
4822 % The format of the ZeroKernelNans method is:
4824 % void ZeroKernelNans (KernelInfo *kernel)
4826 % A description of each parameter follows:
4828 % o kernel: the Morphology/Convolution kernel
4831 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4836 /* do the other kernels in a multi-kernel list first */
4837 if ( kernel->next != (KernelInfo *) NULL)
4838 ZeroKernelNans(kernel->next);
4840 for (i=0; i < (kernel->width*kernel->height); i++)
4841 if ( IsNan(kernel->values[i]) )
4842 kernel->values[i] = 0.0;