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-2013 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/memory-private.h"
68 #include "MagickCore/monitor-private.h"
69 #include "MagickCore/morphology.h"
70 #include "MagickCore/morphology-private.h"
71 #include "MagickCore/option.h"
72 #include "MagickCore/pixel-accessor.h"
73 #include "MagickCore/pixel-private.h"
74 #include "MagickCore/prepress.h"
75 #include "MagickCore/quantize.h"
76 #include "MagickCore/resource_.h"
77 #include "MagickCore/registry.h"
78 #include "MagickCore/semaphore.h"
79 #include "MagickCore/splay-tree.h"
80 #include "MagickCore/statistic.h"
81 #include "MagickCore/string_.h"
82 #include "MagickCore/string-private.h"
83 #include "MagickCore/thread-private.h"
84 #include "MagickCore/token.h"
85 #include "MagickCore/utility.h"
86 #include "MagickCore/utility-private.h"
89 Other global definitions used by module.
91 static inline double MagickMin(const double x,const double y)
93 return( x < y ? x : y);
95 static inline double MagickMax(const double x,const double y)
97 return( x > y ? x : y);
99 #define Minimize(assign,value) assign=MagickMin(assign,value)
100 #define Maximize(assign,value) assign=MagickMax(assign,value)
102 /* Integer Factorial Function - for a Binomial kernel */
104 static inline size_t fact(size_t n)
107 for(f=1, l=2; l <= n; f=f*l, l++);
110 #elif 1 /* glibc floating point alternatives */
111 #define fact(n) ((size_t)tgamma((double)n+1))
113 #define fact(n) ((size_t)lgamma((double)n+1))
117 /* Currently these are only internal to this module */
119 CalcKernelMetaData(KernelInfo *),
120 ExpandMirrorKernelInfo(KernelInfo *),
121 ExpandRotateKernelInfo(KernelInfo *, const double),
122 RotateKernelInfo(KernelInfo *, double);
125 /* Quick function to find last kernel in a kernel list */
126 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
128 while (kernel->next != (KernelInfo *) NULL)
129 kernel = kernel->next;
134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
138 % A c q u i r e K e r n e l I n f o %
142 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
144 % AcquireKernelInfo() takes the given string (generally supplied by the
145 % user) and converts it into a Morphology/Convolution Kernel. This allows
146 % users to specify a kernel from a number of pre-defined kernels, or to fully
147 % specify their own kernel for a specific Convolution or Morphology
150 % The kernel so generated can be any rectangular array of floating point
151 % values (doubles) with the 'control point' or 'pixel being affected'
152 % anywhere within that array of values.
154 % Previously IM was restricted to a square of odd size using the exact
155 % center as origin, this is no longer the case, and any rectangular kernel
156 % with any value being declared the origin. This in turn allows the use of
157 % highly asymmetrical kernels.
159 % The floating point values in the kernel can also include a special value
160 % known as 'nan' or 'not a number' to indicate that this value is not part
161 % of the kernel array. This allows you to shaped the kernel within its
162 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
163 % shape. However at least one non-nan value must be provided for correct
164 % working of a kernel.
166 % The returned kernel should be freed using the DestroyKernelInfo() when you
167 % are finished with it. Do not free this memory yourself.
169 % Input kernel defintion strings can consist of any of three types.
172 % Select from one of the built in kernels, using the name and
173 % geometry arguments supplied. See AcquireKernelBuiltIn()
175 % "WxH[+X+Y][@><]:num, num, num ..."
176 % a kernel of size W by H, with W*H floating point numbers following.
177 % the 'center' can be optionally be defined at +X+Y (such that +0+0
178 % is top left corner). If not defined the pixel in the center, for
179 % odd sizes, or to the immediate top or left of center for even sizes
180 % is automatically selected.
182 % "num, num, num, num, ..."
183 % list of floating point numbers defining an 'old style' odd sized
184 % square kernel. At least 9 values should be provided for a 3x3
185 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
186 % Values can be space or comma separated. This is not recommended.
188 % You can define a 'list of kernels' which can be used by some morphology
189 % operators A list is defined as a semi-colon separated list kernels.
191 % " kernel ; kernel ; kernel ; "
193 % Any extra ';' characters, at start, end or between kernel defintions are
196 % The special flags will expand a single kernel, into a list of rotated
197 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
198 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
199 % The '<' also exands using 90-degree rotates, but giving a 180-degree
200 % reflected kernel before the +/- 90-degree rotations, which can be important
201 % for Thinning operations.
203 % Note that 'name' kernels will start with an alphabetic character while the
204 % new kernel specification has a ':' character in its specification string.
205 % If neither is the case, it is assumed an old style of a simple list of
206 % numbers generating a odd-sized square kernel has been given.
208 % The format of the AcquireKernal method is:
210 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
212 % A description of each parameter follows:
214 % o kernel_string: the Morphology/Convolution kernel wanted.
218 /* This was separated so that it could be used as a separate
219 ** array input handling function, such as for -color-matrix
221 static KernelInfo *ParseKernelArray(const char *kernel_string)
227 token[MaxTextExtent];
237 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
245 kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
246 if (kernel == (KernelInfo *)NULL)
248 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
249 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
250 kernel->negative_range = kernel->positive_range = 0.0;
251 kernel->type = UserDefinedKernel;
252 kernel->next = (KernelInfo *) NULL;
253 kernel->signature = MagickSignature;
254 if (kernel_string == (const char *) NULL)
257 /* find end of this specific kernel definition string */
258 end = strchr(kernel_string, ';');
259 if ( end == (char *) NULL )
260 end = strchr(kernel_string, '\0');
262 /* clear flags - for Expanding kernel lists thorugh rotations */
265 /* Has a ':' in argument - New user kernel specification
266 FUTURE: this split on ':' could be done by StringToken()
268 p = strchr(kernel_string, ':');
269 if ( p != (char *) NULL && p < end)
271 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
272 memcpy(token, kernel_string, (size_t) (p-kernel_string));
273 token[p-kernel_string] = '\0';
274 SetGeometryInfo(&args);
275 flags = ParseGeometry(token, &args);
277 /* Size handling and checks of geometry settings */
278 if ( (flags & WidthValue) == 0 ) /* if no width then */
279 args.rho = args.sigma; /* then width = height */
280 if ( args.rho < 1.0 ) /* if width too small */
281 args.rho = 1.0; /* then width = 1 */
282 if ( args.sigma < 1.0 ) /* if height too small */
283 args.sigma = args.rho; /* then height = width */
284 kernel->width = (size_t)args.rho;
285 kernel->height = (size_t)args.sigma;
287 /* Offset Handling and Checks */
288 if ( args.xi < 0.0 || args.psi < 0.0 )
289 return(DestroyKernelInfo(kernel));
290 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
291 : (ssize_t) (kernel->width-1)/2;
292 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
293 : (ssize_t) (kernel->height-1)/2;
294 if ( kernel->x >= (ssize_t) kernel->width ||
295 kernel->y >= (ssize_t) kernel->height )
296 return(DestroyKernelInfo(kernel));
298 p++; /* advance beyond the ':' */
301 { /* ELSE - Old old specification, forming odd-square kernel */
302 /* count up number of values given */
303 p=(const char *) kernel_string;
304 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
305 p++; /* ignore "'" chars for convolve filter usage - Cristy */
306 for (i=0; p < end; i++)
308 GetMagickToken(p,&p,token);
310 GetMagickToken(p,&p,token);
312 /* set the size of the kernel - old sized square */
313 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
314 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
315 p=(const char *) kernel_string;
316 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
317 p++; /* ignore "'" chars for convolve filter usage - Cristy */
320 /* Read in the kernel values from rest of input string argument */
321 kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
322 kernel->width,kernel->height*sizeof(*kernel->values)));
323 if (kernel->values == (MagickRealType *) NULL)
324 return(DestroyKernelInfo(kernel));
325 kernel->minimum = +MagickHuge;
326 kernel->maximum = -MagickHuge;
327 kernel->negative_range = kernel->positive_range = 0.0;
328 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
330 GetMagickToken(p,&p,token);
332 GetMagickToken(p,&p,token);
333 if ( LocaleCompare("nan",token) == 0
334 || LocaleCompare("-",token) == 0 ) {
335 kernel->values[i] = nan; /* this value is not part of neighbourhood */
338 kernel->values[i] = StringToDouble(token,(char **) NULL);
339 ( kernel->values[i] < 0)
340 ? ( kernel->negative_range += kernel->values[i] )
341 : ( kernel->positive_range += kernel->values[i] );
342 Minimize(kernel->minimum, kernel->values[i]);
343 Maximize(kernel->maximum, kernel->values[i]);
347 /* sanity check -- no more values in kernel definition */
348 GetMagickToken(p,&p,token);
349 if ( *token != '\0' && *token != ';' && *token != '\'' )
350 return(DestroyKernelInfo(kernel));
353 /* this was the old method of handling a incomplete kernel */
354 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
355 Minimize(kernel->minimum, kernel->values[i]);
356 Maximize(kernel->maximum, kernel->values[i]);
357 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
358 kernel->values[i]=0.0;
361 /* Number of values for kernel was not enough - Report Error */
362 if ( i < (ssize_t) (kernel->width*kernel->height) )
363 return(DestroyKernelInfo(kernel));
366 /* check that we recieved at least one real (non-nan) value! */
367 if ( kernel->minimum == MagickHuge )
368 return(DestroyKernelInfo(kernel));
370 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
371 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
372 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
373 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
374 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
375 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
380 static KernelInfo *ParseKernelName(const char *kernel_string)
383 token[MaxTextExtent];
401 /* Parse special 'named' kernel */
402 GetMagickToken(kernel_string,&p,token);
403 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
404 if ( type < 0 || type == UserDefinedKernel )
405 return((KernelInfo *)NULL); /* not a valid named kernel */
407 while (((isspace((int) ((unsigned char) *p)) != 0) ||
408 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
411 end = strchr(p, ';'); /* end of this kernel defintion */
412 if ( end == (char *) NULL )
413 end = strchr(p, '\0');
415 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
416 memcpy(token, p, (size_t) (end-p));
418 SetGeometryInfo(&args);
419 flags = ParseGeometry(token, &args);
422 /* For Debugging Geometry Input */
423 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
424 flags, args.rho, args.sigma, args.xi, args.psi );
427 /* special handling of missing values in input string */
429 /* Shape Kernel Defaults */
431 if ( (flags & WidthValue) == 0 )
432 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
440 if ( (flags & HeightValue) == 0 )
441 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
444 if ( (flags & XValue) == 0 )
445 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
447 case RectangleKernel: /* Rectangle - set size defaults */
448 if ( (flags & WidthValue) == 0 ) /* if no width then */
449 args.rho = args.sigma; /* then width = height */
450 if ( args.rho < 1.0 ) /* if width too small */
451 args.rho = 3; /* then width = 3 */
452 if ( args.sigma < 1.0 ) /* if height too small */
453 args.sigma = args.rho; /* then height = width */
454 if ( (flags & XValue) == 0 ) /* center offset if not defined */
455 args.xi = (double)(((ssize_t)args.rho-1)/2);
456 if ( (flags & YValue) == 0 )
457 args.psi = (double)(((ssize_t)args.sigma-1)/2);
459 /* Distance Kernel Defaults */
460 case ChebyshevKernel:
461 case ManhattanKernel:
462 case OctagonalKernel:
463 case EuclideanKernel:
464 if ( (flags & HeightValue) == 0 ) /* no distance scale */
465 args.sigma = 100.0; /* default distance scaling */
466 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
467 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
468 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
469 args.sigma *= QuantumRange/100.0; /* percentage of color range */
475 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
476 if ( kernel == (KernelInfo *) NULL )
479 /* global expand to rotated kernel list - only for single kernels */
480 if ( kernel->next == (KernelInfo *) NULL ) {
481 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
482 ExpandRotateKernelInfo(kernel, 45.0);
483 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
484 ExpandRotateKernelInfo(kernel, 90.0);
485 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
486 ExpandMirrorKernelInfo(kernel);
492 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
500 token[MaxTextExtent];
508 if (kernel_string == (const char *) NULL)
509 return(ParseKernelArray(kernel_string));
514 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
516 /* ignore extra or multiple ';' kernel separators */
517 if ( *token != ';' ) {
519 /* tokens starting with alpha is a Named kernel */
520 if (isalpha((int) *token) != 0)
521 new_kernel = ParseKernelName(p);
522 else /* otherwise a user defined kernel array */
523 new_kernel = ParseKernelArray(p);
525 /* Error handling -- this is not proper error handling! */
526 if ( new_kernel == (KernelInfo *) NULL ) {
527 (void) FormatLocaleFile(stderr,"Failed to parse kernel number #%.20g\n",
528 (double) kernel_number);
529 if ( kernel != (KernelInfo *) NULL )
530 kernel=DestroyKernelInfo(kernel);
531 return((KernelInfo *) NULL);
534 /* initialise or append the kernel list */
535 if ( kernel == (KernelInfo *) NULL )
538 LastKernelInfo(kernel)->next = new_kernel;
541 /* look for the next kernel in list */
543 if ( p == (char *) NULL )
553 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
557 % A c q u i r e K e r n e l B u i l t I n %
561 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
563 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
564 % kernels used for special purposes such as gaussian blurring, skeleton
565 % pruning, and edge distance determination.
567 % They take a KernelType, and a set of geometry style arguments, which were
568 % typically decoded from a user supplied string, or from a more complex
569 % Morphology Method that was requested.
571 % The format of the AcquireKernalBuiltIn method is:
573 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
574 % const GeometryInfo args)
576 % A description of each parameter follows:
578 % o type: the pre-defined type of kernel wanted
580 % o args: arguments defining or modifying the kernel
582 % Convolution Kernels
585 % The a No-Op or Scaling single element kernel.
587 % Gaussian:{radius},{sigma}
588 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
589 % The sigma for the curve is required. The resulting kernel is
592 % If 'sigma' is zero, you get a single pixel on a field of zeros.
594 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
595 % the final size of the resulting kernel to a square 2*radius+1 in size.
596 % The radius should be at least 2 times that of the sigma value, or
597 % sever clipping and aliasing may result. If not given or set to 0 the
598 % radius will be determined so as to produce the best minimal error
599 % result, which is usally much larger than is normally needed.
601 % LoG:{radius},{sigma}
602 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
603 % The supposed ideal edge detection, zero-summing kernel.
605 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
606 % approx 1.6 (according to wikipedia).
608 % DoG:{radius},{sigma1},{sigma2}
609 % "Difference of Gaussians" Kernel.
610 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
611 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
612 % The result is a zero-summing kernel.
614 % Blur:{radius},{sigma}[,{angle}]
615 % Generates a 1 dimensional or linear gaussian blur, at the angle given
616 % (current restricted to orthogonal angles). If a 'radius' is given the
617 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
618 % by a 90 degree angle.
620 % If 'sigma' is zero, you get a single pixel on a field of zeros.
622 % Note that two convolutions with two "Blur" kernels perpendicular to
623 % each other, is equivalent to a far larger "Gaussian" kernel with the
624 % same sigma value, However it is much faster to apply. This is how the
625 % "-blur" operator actually works.
627 % Comet:{width},{sigma},{angle}
628 % Blur in one direction only, much like how a bright object leaves
629 % a comet like trail. The Kernel is actually half a gaussian curve,
630 % Adding two such blurs in opposite directions produces a Blur Kernel.
631 % Angle can be rotated in multiples of 90 degrees.
633 % Note that the first argument is the width of the kernel and not the
634 % radius of the kernel.
636 % Binomial:[{radius}]
637 % Generate a discrete kernel using a 2 dimentional Pascel's Triangle
638 % of values. Used for special forma of image filters.
640 % # Still to be implemented...
644 % # Set kernel values using a resize filter, and given scale (sigma)
645 % # Cylindrical or Linear. Is this possible with an image?
648 % Named Constant Convolution Kernels
650 % All these are unscaled, zero-summing kernels by default. As such for
651 % non-HDRI version of ImageMagick some form of normalization, user scaling,
652 % and biasing the results is recommended, to prevent the resulting image
655 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
656 % 45 degrees to generate the 8 angled varients of each of the kernels.
659 % Discrete Lapacian Kernels, (without normalization)
660 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
661 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
662 % Type 2 : 3x3 with center:4 edge:1 corner:-2
663 % Type 3 : 3x3 with center:4 edge:-2 corner:1
664 % Type 5 : 5x5 laplacian
665 % Type 7 : 7x7 laplacian
666 % Type 15 : 5x5 LoG (sigma approx 1.4)
667 % Type 19 : 9x9 LoG (sigma approx 1.4)
670 % Sobel 'Edge' convolution kernel (3x3)
676 % Roberts convolution kernel (3x3)
682 % Prewitt Edge convolution kernel (3x3)
688 % Prewitt's "Compass" convolution kernel (3x3)
694 % Kirsch's "Compass" convolution kernel (3x3)
700 % Frei-Chen Edge Detector is based on a kernel that is similar to
701 % the Sobel Kernel, but is designed to be isotropic. That is it takes
702 % into account the distance of the diagonal in the kernel.
705 % | sqrt(2), 0, -sqrt(2) |
708 % FreiChen:{type},{angle}
710 % Frei-Chen Pre-weighted kernels...
712 % Type 0: default un-nomalized version shown above.
714 % Type 1: Orthogonal Kernel (same as type 11 below)
716 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
719 % Type 2: Diagonal form of Kernel...
721 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
724 % However this kernel is als at the heart of the FreiChen Edge Detection
725 % Process which uses a set of 9 specially weighted kernel. These 9
726 % kernels not be normalized, but directly applied to the image. The
727 % results is then added together, to produce the intensity of an edge in
728 % a specific direction. The square root of the pixel value can then be
729 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
730 % from each other, both the direction and the strength of the edge can be
733 % Type 10: All 9 of the following pre-weighted kernels...
735 % Type 11: | 1, 0, -1 |
736 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
739 % Type 12: | 1, sqrt(2), 1 |
740 % | 0, 0, 0 | / 2*sqrt(2)
743 % Type 13: | sqrt(2), -1, 0 |
744 % | -1, 0, 1 | / 2*sqrt(2)
747 % Type 14: | 0, 1, -sqrt(2) |
748 % | -1, 0, 1 | / 2*sqrt(2)
751 % Type 15: | 0, -1, 0 |
755 % Type 16: | 1, 0, -1 |
759 % Type 17: | 1, -2, 1 |
763 % Type 18: | -2, 1, -2 |
767 % Type 19: | 1, 1, 1 |
771 % The first 4 are for edge detection, the next 4 are for line detection
772 % and the last is to add a average component to the results.
774 % Using a special type of '-1' will return all 9 pre-weighted kernels
775 % as a multi-kernel list, so that you can use them directly (without
776 % normalization) with the special "-set option:morphology:compose Plus"
777 % setting to apply the full FreiChen Edge Detection Technique.
779 % If 'type' is large it will be taken to be an actual rotation angle for
780 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
781 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
783 % WARNING: The above was layed out as per
784 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
785 % But rotated 90 degrees so direction is from left rather than the top.
786 % I have yet to find any secondary confirmation of the above. The only
787 % other source found was actual source code at
788 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
789 % Neigher paper defineds the kernels in a way that looks locical or
790 % correct when taken as a whole.
794 % Diamond:[{radius}[,{scale}]]
795 % Generate a diamond shaped kernel with given radius to the points.
796 % Kernel size will again be radius*2+1 square and defaults to radius 1,
797 % generating a 3x3 kernel that is slightly larger than a square.
799 % Square:[{radius}[,{scale}]]
800 % Generate a square shaped kernel of size radius*2+1, and defaulting
801 % to a 3x3 (radius 1).
803 % Octagon:[{radius}[,{scale}]]
804 % Generate octagonal shaped kernel of given radius and constant scale.
805 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
806 % in "Diamond" kernel.
808 % Disk:[{radius}[,{scale}]]
809 % Generate a binary disk, thresholded at the radius given, the radius
810 % may be a float-point value. Final Kernel size is floor(radius)*2+1
811 % square. A radius of 5.3 is the default.
813 % NOTE: That a low radii Disk kernels produce the same results as
814 % many of the previously defined kernels, but differ greatly at larger
815 % radii. Here is a table of equivalences...
816 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
817 % "Disk:1.5" => "Square"
818 % "Disk:2" => "Diamond:2"
819 % "Disk:2.5" => "Octagon"
820 % "Disk:2.9" => "Square:2"
821 % "Disk:3.5" => "Octagon:3"
822 % "Disk:4.5" => "Octagon:4"
823 % "Disk:5.4" => "Octagon:5"
824 % "Disk:6.4" => "Octagon:6"
825 % All other Disk shapes are unique to this kernel, but because a "Disk"
826 % is more circular when using a larger radius, using a larger radius is
827 % preferred over iterating the morphological operation.
829 % Rectangle:{geometry}
830 % Simply generate a rectangle of 1's with the size given. You can also
831 % specify the location of the 'control point', otherwise the closest
832 % pixel to the center of the rectangle is selected.
834 % Properly centered and odd sized rectangles work the best.
836 % Symbol Dilation Kernels
838 % These kernel is not a good general morphological kernel, but is used
839 % more for highlighting and marking any single pixels in an image using,
840 % a "Dilate" method as appropriate.
842 % For the same reasons iterating these kernels does not produce the
843 % same result as using a larger radius for the symbol.
845 % Plus:[{radius}[,{scale}]]
846 % Cross:[{radius}[,{scale}]]
847 % Generate a kernel in the shape of a 'plus' or a 'cross' with
848 % a each arm the length of the given radius (default 2).
850 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
852 % Ring:{radius1},{radius2}[,{scale}]
853 % A ring of the values given that falls between the two radii.
854 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
855 % This is the 'edge' pixels of the default "Disk" kernel,
856 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
858 % Hit and Miss Kernels
860 % Peak:radius1,radius2
861 % Find any peak larger than the pixels the fall between the two radii.
862 % The default ring of pixels is as per "Ring".
864 % Find flat orthogonal edges of a binary shape
866 % Find 90 degree corners of a binary shape
868 % A special kernel to thin the 'outside' of diagonals
870 % Find end points of lines (for pruning a skeletion)
871 % Two types of lines ends (default to both) can be searched for
872 % Type 0: All line ends
873 % Type 1: single kernel for 4-conneected line ends
874 % Type 2: single kernel for simple line ends
876 % Find three line junctions (within a skeletion)
877 % Type 0: all line junctions
878 % Type 1: Y Junction kernel
879 % Type 2: Diagonal T Junction kernel
880 % Type 3: Orthogonal T Junction kernel
881 % Type 4: Diagonal X Junction kernel
882 % Type 5: Orthogonal + Junction kernel
884 % Find single pixel ridges or thin lines
885 % Type 1: Fine single pixel thick lines and ridges
886 % Type 2: Find two pixel thick lines and ridges
888 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
890 % Traditional skeleton generating kernels.
891 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
892 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
893 % Type 3: Thinning skeleton based on a ressearch paper by
894 % Dan S. Bloomberg (Default Type)
896 % A huge variety of Thinning Kernels designed to preserve conectivity.
897 % many other kernel sets use these kernels as source definitions.
898 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
899 % the super and sub notations used in the source research paper.
901 % Distance Measuring Kernels
903 % Different types of distance measuring methods, which are used with the
904 % a 'Distance' morphology method for generating a gradient based on
905 % distance from an edge of a binary shape, though there is a technique
906 % for handling a anti-aliased shape.
908 % See the 'Distance' Morphological Method, for information of how it is
911 % Chebyshev:[{radius}][x{scale}[%!]]
912 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
913 % is a value of one to any neighbour, orthogonal or diagonal. One why
914 % of thinking of it is the number of squares a 'King' or 'Queen' in
915 % chess needs to traverse reach any other position on a chess board.
916 % It results in a 'square' like distance function, but one where
917 % diagonals are given a value that is closer than expected.
919 % Manhattan:[{radius}][x{scale}[%!]]
920 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
921 % Cab distance metric), it is the distance needed when you can only
922 % travel in horizontal or vertical directions only. It is the
923 % distance a 'Rook' in chess would have to travel, and results in a
924 % diamond like distances, where diagonals are further than expected.
926 % Octagonal:[{radius}][x{scale}[%!]]
927 % An interleving of Manhatten and Chebyshev metrics producing an
928 % increasing octagonally shaped distance. Distances matches those of
929 % the "Octagon" shaped kernel of the same radius. The minimum radius
930 % and default is 2, producing a 5x5 kernel.
932 % Euclidean:[{radius}][x{scale}[%!]]
933 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
934 % However by default the kernel size only has a radius of 1, which
935 % limits the distance to 'Knight' like moves, with only orthogonal and
936 % diagonal measurements being correct. As such for the default kernel
937 % you will get octagonal like distance function.
939 % However using a larger radius such as "Euclidean:4" you will get a
940 % much smoother distance gradient from the edge of the shape. Especially
941 % if the image is pre-processed to include any anti-aliasing pixels.
942 % Of course a larger kernel is slower to use, and not always needed.
944 % The first three Distance Measuring Kernels will only generate distances
945 % of exact multiples of {scale} in binary images. As such you can use a
946 % scale of 1 without loosing any information. However you also need some
947 % scaling when handling non-binary anti-aliased shapes.
949 % The "Euclidean" Distance Kernel however does generate a non-integer
950 % fractional results, and as such scaling is vital even for binary shapes.
954 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
955 const GeometryInfo *args)
968 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
970 /* Generate a new empty kernel if needed */
971 kernel=(KernelInfo *) NULL;
973 case UndefinedKernel: /* These should not call this function */
974 case UserDefinedKernel:
975 assert("Should not call this function" != (char *)NULL);
977 case LaplacianKernel: /* Named Descrete Convolution Kernels */
978 case SobelKernel: /* these are defined using other kernels */
984 case EdgesKernel: /* Hit and Miss kernels */
986 case DiagonalsKernel:
988 case LineJunctionsKernel:
990 case ConvexHullKernel:
993 break; /* A pre-generated kernel is not needed */
995 /* set to 1 to do a compile-time check that we haven't missed anything */
1002 case BinomialKernel:
1005 case RectangleKernel:
1012 case ChebyshevKernel:
1013 case ManhattanKernel:
1014 case OctangonalKernel:
1015 case EuclideanKernel:
1019 /* Generate the base Kernel Structure */
1020 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1021 if (kernel == (KernelInfo *) NULL)
1023 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1024 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1025 kernel->negative_range = kernel->positive_range = 0.0;
1026 kernel->type = type;
1027 kernel->next = (KernelInfo *) NULL;
1028 kernel->signature = MagickSignature;
1038 kernel->height = kernel->width = (size_t) 1;
1039 kernel->x = kernel->y = (ssize_t) 0;
1040 kernel->values=(MagickRealType *) MagickAssumeAligned(
1041 AcquireAlignedMemory(1,sizeof(*kernel->values)));
1042 if (kernel->values == (MagickRealType *) NULL)
1043 return(DestroyKernelInfo(kernel));
1044 kernel->maximum = kernel->values[0] = args->rho;
1048 case GaussianKernel:
1052 sigma = fabs(args->sigma),
1053 sigma2 = fabs(args->xi),
1056 if ( args->rho >= 1.0 )
1057 kernel->width = (size_t)args->rho*2+1;
1058 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1059 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1061 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1062 kernel->height = kernel->width;
1063 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1064 kernel->values=(MagickRealType *) MagickAssumeAligned(
1065 AcquireAlignedMemory(kernel->width,kernel->height*
1066 sizeof(*kernel->values)));
1067 if (kernel->values == (MagickRealType *) NULL)
1068 return(DestroyKernelInfo(kernel));
1070 /* WARNING: The following generates a 'sampled gaussian' kernel.
1071 * What we really want is a 'discrete gaussian' kernel.
1073 * How to do this is I don't know, but appears to be basied on the
1074 * Error Function 'erf()' (intergral of a gaussian)
1077 if ( type == GaussianKernel || type == DoGKernel )
1078 { /* Calculate a Gaussian, OR positive half of a DoG */
1079 if ( sigma > MagickEpsilon )
1080 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1081 B = (double) (1.0/(Magick2PI*sigma*sigma));
1082 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1083 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1084 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1086 else /* limiting case - a unity (normalized Dirac) kernel */
1087 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1088 kernel->width*kernel->height*sizeof(*kernel->values));
1089 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1093 if ( type == DoGKernel )
1094 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1095 if ( sigma2 > MagickEpsilon )
1096 { sigma = sigma2; /* simplify loop expressions */
1097 A = 1.0/(2.0*sigma*sigma);
1098 B = (double) (1.0/(Magick2PI*sigma*sigma));
1099 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1100 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1101 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1103 else /* limiting case - a unity (normalized Dirac) kernel */
1104 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1107 if ( type == LoGKernel )
1108 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1109 if ( sigma > MagickEpsilon )
1110 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1111 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1112 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1113 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1114 { R = ((double)(u*u+v*v))*A;
1115 kernel->values[i] = (1-R)*exp(-R)*B;
1118 else /* special case - generate a unity kernel */
1119 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1120 kernel->width*kernel->height*sizeof(*kernel->values));
1121 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1125 /* Note the above kernels may have been 'clipped' by a user defined
1126 ** radius, producing a smaller (darker) kernel. Also for very small
1127 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1128 ** producing a very bright kernel.
1130 ** Normalization will still be needed.
1133 /* Normalize the 2D Gaussian Kernel
1135 ** NB: a CorrelateNormalize performs a normal Normalize if
1136 ** there are no negative values.
1138 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1139 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1145 sigma = fabs(args->sigma),
1148 if ( args->rho >= 1.0 )
1149 kernel->width = (size_t)args->rho*2+1;
1151 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1153 kernel->x = (ssize_t) (kernel->width-1)/2;
1155 kernel->negative_range = kernel->positive_range = 0.0;
1156 kernel->values=(MagickRealType *) MagickAssumeAligned(
1157 AcquireAlignedMemory(kernel->width,kernel->height*
1158 sizeof(*kernel->values)));
1159 if (kernel->values == (MagickRealType *) NULL)
1160 return(DestroyKernelInfo(kernel));
1163 #define KernelRank 3
1164 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1165 ** It generates a gaussian 3 times the width, and compresses it into
1166 ** the expected range. This produces a closer normalization of the
1167 ** resulting kernel, especially for very low sigma values.
1168 ** As such while wierd it is prefered.
1170 ** I am told this method originally came from Photoshop.
1172 ** A properly normalized curve is generated (apart from edge clipping)
1173 ** even though we later normalize the result (for edge clipping)
1174 ** to allow the correct generation of a "Difference of Blurs".
1178 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1179 (void) ResetMagickMemory(kernel->values,0, (size_t)
1180 kernel->width*kernel->height*sizeof(*kernel->values));
1181 /* Calculate a Positive 1D Gaussian */
1182 if ( sigma > MagickEpsilon )
1183 { sigma *= KernelRank; /* simplify loop expressions */
1184 alpha = 1.0/(2.0*sigma*sigma);
1185 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1186 for ( u=-v; u <= v; u++) {
1187 kernel->values[(u+v)/KernelRank] +=
1188 exp(-((double)(u*u))*alpha)*beta;
1191 else /* special case - generate a unity kernel */
1192 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1194 /* Direct calculation without curve averaging
1195 This is equivelent to a KernelRank of 1 */
1197 /* Calculate a Positive Gaussian */
1198 if ( sigma > MagickEpsilon )
1199 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1200 beta = 1.0/(MagickSQ2PI*sigma);
1201 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1202 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1204 else /* special case - generate a unity kernel */
1205 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1206 kernel->width*kernel->height*sizeof(*kernel->values));
1207 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1210 /* Note the above kernel may have been 'clipped' by a user defined
1211 ** radius, producing a smaller (darker) kernel. Also for very small
1212 ** sigma's (> 0.1) the central value becomes larger than one, as a
1213 ** result of not generating a actual 'discrete' kernel, and thus
1214 ** producing a very bright 'impulse'.
1216 ** Becuase of these two factors Normalization is required!
1219 /* Normalize the 1D Gaussian Kernel
1221 ** NB: a CorrelateNormalize performs a normal Normalize if
1222 ** there are no negative values.
1224 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1225 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1227 /* rotate the 1D kernel by given angle */
1228 RotateKernelInfo(kernel, args->xi );
1233 sigma = fabs(args->sigma),
1236 if ( args->rho < 1.0 )
1237 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1239 kernel->width = (size_t)args->rho;
1240 kernel->x = kernel->y = 0;
1242 kernel->negative_range = kernel->positive_range = 0.0;
1243 kernel->values=(MagickRealType *) MagickAssumeAligned(
1244 AcquireAlignedMemory(kernel->width,kernel->height*
1245 sizeof(*kernel->values)));
1246 if (kernel->values == (MagickRealType *) NULL)
1247 return(DestroyKernelInfo(kernel));
1249 /* A comet blur is half a 1D gaussian curve, so that the object is
1250 ** blurred in one direction only. This may not be quite the right
1251 ** curve to use so may change in the future. The function must be
1252 ** normalised after generation, which also resolves any clipping.
1254 ** As we are normalizing and not subtracting gaussians,
1255 ** there is no need for a divisor in the gaussian formula
1257 ** It is less comples
1259 if ( sigma > MagickEpsilon )
1262 #define KernelRank 3
1263 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1264 (void) ResetMagickMemory(kernel->values,0, (size_t)
1265 kernel->width*sizeof(*kernel->values));
1266 sigma *= KernelRank; /* simplify the loop expression */
1267 A = 1.0/(2.0*sigma*sigma);
1268 /* B = 1.0/(MagickSQ2PI*sigma); */
1269 for ( u=0; u < v; u++) {
1270 kernel->values[u/KernelRank] +=
1271 exp(-((double)(u*u))*A);
1272 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1274 for (i=0; i < (ssize_t) kernel->width; i++)
1275 kernel->positive_range += kernel->values[i];
1277 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1278 /* B = 1.0/(MagickSQ2PI*sigma); */
1279 for ( i=0; i < (ssize_t) kernel->width; i++)
1280 kernel->positive_range +=
1281 kernel->values[i] = exp(-((double)(i*i))*A);
1282 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1285 else /* special case - generate a unity kernel */
1286 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1287 kernel->width*kernel->height*sizeof(*kernel->values));
1288 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1289 kernel->positive_range = 1.0;
1292 kernel->minimum = 0.0;
1293 kernel->maximum = kernel->values[0];
1294 kernel->negative_range = 0.0;
1296 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1297 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1300 case BinomialKernel:
1305 if (args->rho < 1.0)
1306 kernel->width = kernel->height = 3; /* default radius = 1 */
1308 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1309 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1311 order_f = fact(kernel->width-1);
1313 kernel->values=(MagickRealType *) MagickAssumeAligned(
1314 AcquireAlignedMemory(kernel->width,kernel->height*
1315 sizeof(*kernel->values)));
1316 if (kernel->values == (MagickRealType *) NULL)
1317 return(DestroyKernelInfo(kernel));
1319 /* set all kernel values within diamond area to scale given */
1320 for ( i=0, v=0; v < (ssize_t)kernel->height; v++)
1322 alpha = order_f / ( fact((size_t) v) * fact(kernel->height-v-1) );
1323 for ( u=0; u < (ssize_t)kernel->width; u++, i++)
1324 kernel->positive_range += kernel->values[i] = (double)
1325 (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-u-1) ));
1327 kernel->minimum = 1.0;
1328 kernel->maximum = kernel->values[kernel->x+kernel->y*kernel->width];
1329 kernel->negative_range = 0.0;
1334 Convolution Kernels - Well Known Named Constant Kernels
1336 case LaplacianKernel:
1337 { switch ( (int) args->rho ) {
1339 default: /* laplacian square filter -- default */
1340 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1342 case 1: /* laplacian diamond filter */
1343 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1346 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1349 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1351 case 5: /* a 5x5 laplacian */
1352 kernel=ParseKernelArray(
1353 "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");
1355 case 7: /* a 7x7 laplacian */
1356 kernel=ParseKernelArray(
1357 "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" );
1359 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1360 kernel=ParseKernelArray(
1361 "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");
1363 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1364 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1365 kernel=ParseKernelArray(
1366 "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");
1369 if (kernel == (KernelInfo *) NULL)
1371 kernel->type = type;
1375 { /* Simple Sobel Kernel */
1376 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1377 if (kernel == (KernelInfo *) NULL)
1379 kernel->type = type;
1380 RotateKernelInfo(kernel, args->rho);
1385 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1386 if (kernel == (KernelInfo *) NULL)
1388 kernel->type = type;
1389 RotateKernelInfo(kernel, args->rho);
1394 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1395 if (kernel == (KernelInfo *) NULL)
1397 kernel->type = type;
1398 RotateKernelInfo(kernel, args->rho);
1403 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1404 if (kernel == (KernelInfo *) NULL)
1406 kernel->type = type;
1407 RotateKernelInfo(kernel, args->rho);
1412 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1413 if (kernel == (KernelInfo *) NULL)
1415 kernel->type = type;
1416 RotateKernelInfo(kernel, args->rho);
1419 case FreiChenKernel:
1420 /* Direction is set to be left to right positive */
1421 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1422 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1423 { switch ( (int) args->rho ) {
1426 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1427 if (kernel == (KernelInfo *) NULL)
1429 kernel->type = type;
1430 kernel->values[3] = +(MagickRealType) MagickSQ2;
1431 kernel->values[5] = -(MagickRealType) MagickSQ2;
1432 CalcKernelMetaData(kernel); /* recalculate meta-data */
1435 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1436 if (kernel == (KernelInfo *) NULL)
1438 kernel->type = type;
1439 kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2;
1440 kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2;
1441 CalcKernelMetaData(kernel); /* recalculate meta-data */
1442 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1445 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1446 if (kernel == (KernelInfo *) NULL)
1451 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1452 if (kernel == (KernelInfo *) NULL)
1454 kernel->type = type;
1455 kernel->values[3] = +(MagickRealType) MagickSQ2;
1456 kernel->values[5] = -(MagickRealType) MagickSQ2;
1457 CalcKernelMetaData(kernel); /* recalculate meta-data */
1458 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1461 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1462 if (kernel == (KernelInfo *) NULL)
1464 kernel->type = type;
1465 kernel->values[1] = +(MagickRealType) MagickSQ2;
1466 kernel->values[7] = +(MagickRealType) MagickSQ2;
1467 CalcKernelMetaData(kernel);
1468 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1471 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1472 if (kernel == (KernelInfo *) NULL)
1474 kernel->type = type;
1475 kernel->values[0] = +(MagickRealType) MagickSQ2;
1476 kernel->values[8] = -(MagickRealType) MagickSQ2;
1477 CalcKernelMetaData(kernel);
1478 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1481 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1482 if (kernel == (KernelInfo *) NULL)
1484 kernel->type = type;
1485 kernel->values[2] = -(MagickRealType) MagickSQ2;
1486 kernel->values[6] = +(MagickRealType) MagickSQ2;
1487 CalcKernelMetaData(kernel);
1488 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1491 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1492 if (kernel == (KernelInfo *) NULL)
1494 kernel->type = type;
1495 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1498 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1499 if (kernel == (KernelInfo *) NULL)
1501 kernel->type = type;
1502 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1505 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1506 if (kernel == (KernelInfo *) NULL)
1508 kernel->type = type;
1509 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1512 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1513 if (kernel == (KernelInfo *) NULL)
1515 kernel->type = type;
1516 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1519 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1520 if (kernel == (KernelInfo *) NULL)
1522 kernel->type = type;
1523 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1526 if ( fabs(args->sigma) >= MagickEpsilon )
1527 /* Rotate by correctly supplied 'angle' */
1528 RotateKernelInfo(kernel, args->sigma);
1529 else if ( args->rho > 30.0 || args->rho < -30.0 )
1530 /* Rotate by out of bounds 'type' */
1531 RotateKernelInfo(kernel, args->rho);
1536 Boolean or Shaped Kernels
1540 if (args->rho < 1.0)
1541 kernel->width = kernel->height = 3; /* default radius = 1 */
1543 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1544 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1546 kernel->values=(MagickRealType *) MagickAssumeAligned(
1547 AcquireAlignedMemory(kernel->width,kernel->height*
1548 sizeof(*kernel->values)));
1549 if (kernel->values == (MagickRealType *) NULL)
1550 return(DestroyKernelInfo(kernel));
1552 /* set all kernel values within diamond area to scale given */
1553 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1554 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1555 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1556 kernel->positive_range += kernel->values[i] = args->sigma;
1558 kernel->values[i] = nan;
1559 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1563 case RectangleKernel:
1566 if ( type == SquareKernel )
1568 if (args->rho < 1.0)
1569 kernel->width = kernel->height = 3; /* default radius = 1 */
1571 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1572 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1573 scale = args->sigma;
1576 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1577 if ( args->rho < 1.0 || args->sigma < 1.0 )
1578 return(DestroyKernelInfo(kernel)); /* invalid args given */
1579 kernel->width = (size_t)args->rho;
1580 kernel->height = (size_t)args->sigma;
1581 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1582 args->psi < 0.0 || args->psi > (double)kernel->height )
1583 return(DestroyKernelInfo(kernel)); /* invalid args given */
1584 kernel->x = (ssize_t) args->xi;
1585 kernel->y = (ssize_t) args->psi;
1588 kernel->values=(MagickRealType *) MagickAssumeAligned(
1589 AcquireAlignedMemory(kernel->width,kernel->height*
1590 sizeof(*kernel->values)));
1591 if (kernel->values == (MagickRealType *) NULL)
1592 return(DestroyKernelInfo(kernel));
1594 /* set all kernel values to scale given */
1595 u=(ssize_t) (kernel->width*kernel->height);
1596 for ( i=0; i < u; i++)
1597 kernel->values[i] = scale;
1598 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1599 kernel->positive_range = scale*u;
1604 if (args->rho < 1.0)
1605 kernel->width = kernel->height = 5; /* default radius = 2 */
1607 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1608 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1610 kernel->values=(MagickRealType *) MagickAssumeAligned(
1611 AcquireAlignedMemory(kernel->width,kernel->height*
1612 sizeof(*kernel->values)));
1613 if (kernel->values == (MagickRealType *) NULL)
1614 return(DestroyKernelInfo(kernel));
1616 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1617 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1618 if ( (labs((long) u)+labs((long) v)) <=
1619 ((long)kernel->x + (long)(kernel->x/2)) )
1620 kernel->positive_range += kernel->values[i] = args->sigma;
1622 kernel->values[i] = nan;
1623 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1629 limit = (ssize_t)(args->rho*args->rho);
1631 if (args->rho < 0.4) /* default radius approx 4.3 */
1632 kernel->width = kernel->height = 9L, limit = 18L;
1634 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1635 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1637 kernel->values=(MagickRealType *) MagickAssumeAligned(
1638 AcquireAlignedMemory(kernel->width,kernel->height*
1639 sizeof(*kernel->values)));
1640 if (kernel->values == (MagickRealType *) NULL)
1641 return(DestroyKernelInfo(kernel));
1643 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1644 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1645 if ((u*u+v*v) <= limit)
1646 kernel->positive_range += kernel->values[i] = args->sigma;
1648 kernel->values[i] = nan;
1649 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1654 if (args->rho < 1.0)
1655 kernel->width = kernel->height = 5; /* default radius 2 */
1657 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1658 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1660 kernel->values=(MagickRealType *) MagickAssumeAligned(
1661 AcquireAlignedMemory(kernel->width,kernel->height*
1662 sizeof(*kernel->values)));
1663 if (kernel->values == (MagickRealType *) NULL)
1664 return(DestroyKernelInfo(kernel));
1666 /* set all kernel values along axises to given scale */
1667 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1668 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1669 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1670 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1671 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1676 if (args->rho < 1.0)
1677 kernel->width = kernel->height = 5; /* default radius 2 */
1679 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1680 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1682 kernel->values=(MagickRealType *) MagickAssumeAligned(
1683 AcquireAlignedMemory(kernel->width,kernel->height*
1684 sizeof(*kernel->values)));
1685 if (kernel->values == (MagickRealType *) NULL)
1686 return(DestroyKernelInfo(kernel));
1688 /* set all kernel values along axises to given scale */
1689 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1690 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1691 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1692 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1693 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1707 if (args->rho < args->sigma)
1709 kernel->width = ((size_t)args->sigma)*2+1;
1710 limit1 = (ssize_t)(args->rho*args->rho);
1711 limit2 = (ssize_t)(args->sigma*args->sigma);
1715 kernel->width = ((size_t)args->rho)*2+1;
1716 limit1 = (ssize_t)(args->sigma*args->sigma);
1717 limit2 = (ssize_t)(args->rho*args->rho);
1720 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1722 kernel->height = kernel->width;
1723 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1724 kernel->values=(MagickRealType *) MagickAssumeAligned(
1725 AcquireAlignedMemory(kernel->width,kernel->height*
1726 sizeof(*kernel->values)));
1727 if (kernel->values == (MagickRealType *) NULL)
1728 return(DestroyKernelInfo(kernel));
1730 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1731 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1732 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1733 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1734 { ssize_t radius=u*u+v*v;
1735 if (limit1 < radius && radius <= limit2)
1736 kernel->positive_range += kernel->values[i] = (double) scale;
1738 kernel->values[i] = nan;
1740 kernel->minimum = kernel->maximum = (double) scale;
1741 if ( type == PeaksKernel ) {
1742 /* set the central point in the middle */
1743 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1744 kernel->positive_range = 1.0;
1745 kernel->maximum = 1.0;
1751 kernel=AcquireKernelInfo("ThinSE:482");
1752 if (kernel == (KernelInfo *) NULL)
1754 kernel->type = type;
1755 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1760 kernel=AcquireKernelInfo("ThinSE:87");
1761 if (kernel == (KernelInfo *) NULL)
1763 kernel->type = type;
1764 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1767 case DiagonalsKernel:
1769 switch ( (int) args->rho ) {
1774 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1775 if (kernel == (KernelInfo *) NULL)
1777 kernel->type = type;
1778 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1779 if (new_kernel == (KernelInfo *) NULL)
1780 return(DestroyKernelInfo(kernel));
1781 new_kernel->type = type;
1782 LastKernelInfo(kernel)->next = new_kernel;
1783 ExpandMirrorKernelInfo(kernel);
1787 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1790 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1793 if (kernel == (KernelInfo *) NULL)
1795 kernel->type = type;
1796 RotateKernelInfo(kernel, args->sigma);
1799 case LineEndsKernel:
1800 { /* Kernels for finding the end of thin lines */
1801 switch ( (int) args->rho ) {
1804 /* set of kernels to find all end of lines */
1805 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1807 /* kernel for 4-connected line ends - no rotation */
1808 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1811 /* kernel to add for 8-connected lines - no rotation */
1812 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1815 /* kernel to add for orthogonal line ends - does not find corners */
1816 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1819 /* traditional line end - fails on last T end */
1820 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1823 if (kernel == (KernelInfo *) NULL)
1825 kernel->type = type;
1826 RotateKernelInfo(kernel, args->sigma);
1829 case LineJunctionsKernel:
1830 { /* kernels for finding the junctions of multiple lines */
1831 switch ( (int) args->rho ) {
1834 /* set of kernels to find all line junctions */
1835 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1838 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1841 /* Diagonal T Junctions */
1842 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1845 /* Orthogonal T Junctions */
1846 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1849 /* Diagonal X Junctions */
1850 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1853 /* Orthogonal X Junctions - minimal diamond kernel */
1854 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1857 if (kernel == (KernelInfo *) NULL)
1859 kernel->type = type;
1860 RotateKernelInfo(kernel, args->sigma);
1864 { /* Ridges - Ridge finding kernels */
1867 switch ( (int) args->rho ) {
1870 kernel=ParseKernelArray("3x1:0,1,0");
1871 if (kernel == (KernelInfo *) NULL)
1873 kernel->type = type;
1874 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1877 kernel=ParseKernelArray("4x1:0,1,1,0");
1878 if (kernel == (KernelInfo *) NULL)
1880 kernel->type = type;
1881 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1883 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1884 /* Unfortunatally we can not yet rotate a non-square kernel */
1885 /* But then we can't flip a non-symetrical kernel either */
1886 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1887 if (new_kernel == (KernelInfo *) NULL)
1888 return(DestroyKernelInfo(kernel));
1889 new_kernel->type = type;
1890 LastKernelInfo(kernel)->next = new_kernel;
1891 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1892 if (new_kernel == (KernelInfo *) NULL)
1893 return(DestroyKernelInfo(kernel));
1894 new_kernel->type = type;
1895 LastKernelInfo(kernel)->next = new_kernel;
1896 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1897 if (new_kernel == (KernelInfo *) NULL)
1898 return(DestroyKernelInfo(kernel));
1899 new_kernel->type = type;
1900 LastKernelInfo(kernel)->next = new_kernel;
1901 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1902 if (new_kernel == (KernelInfo *) NULL)
1903 return(DestroyKernelInfo(kernel));
1904 new_kernel->type = type;
1905 LastKernelInfo(kernel)->next = new_kernel;
1906 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1907 if (new_kernel == (KernelInfo *) NULL)
1908 return(DestroyKernelInfo(kernel));
1909 new_kernel->type = type;
1910 LastKernelInfo(kernel)->next = new_kernel;
1911 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1912 if (new_kernel == (KernelInfo *) NULL)
1913 return(DestroyKernelInfo(kernel));
1914 new_kernel->type = type;
1915 LastKernelInfo(kernel)->next = new_kernel;
1916 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1917 if (new_kernel == (KernelInfo *) NULL)
1918 return(DestroyKernelInfo(kernel));
1919 new_kernel->type = type;
1920 LastKernelInfo(kernel)->next = new_kernel;
1921 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1922 if (new_kernel == (KernelInfo *) NULL)
1923 return(DestroyKernelInfo(kernel));
1924 new_kernel->type = type;
1925 LastKernelInfo(kernel)->next = new_kernel;
1930 case ConvexHullKernel:
1934 /* first set of 8 kernels */
1935 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1936 if (kernel == (KernelInfo *) NULL)
1938 kernel->type = type;
1939 ExpandRotateKernelInfo(kernel, 90.0);
1940 /* append the mirror versions too - no flip function yet */
1941 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1942 if (new_kernel == (KernelInfo *) NULL)
1943 return(DestroyKernelInfo(kernel));
1944 new_kernel->type = type;
1945 ExpandRotateKernelInfo(new_kernel, 90.0);
1946 LastKernelInfo(kernel)->next = new_kernel;
1949 case SkeletonKernel:
1951 switch ( (int) args->rho ) {
1954 /* Traditional Skeleton...
1955 ** A cyclically rotated single kernel
1957 kernel=AcquireKernelInfo("ThinSE:482");
1958 if (kernel == (KernelInfo *) NULL)
1960 kernel->type = type;
1961 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1964 /* HIPR Variation of the cyclic skeleton
1965 ** Corners of the traditional method made more forgiving,
1966 ** but the retain the same cyclic order.
1968 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1969 if (kernel == (KernelInfo *) NULL)
1971 if (kernel->next == (KernelInfo *) NULL)
1972 return(DestroyKernelInfo(kernel));
1973 kernel->type = type;
1974 kernel->next->type = type;
1975 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1978 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1979 ** "Connectivity-Preserving Morphological Image Thransformations"
1980 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1981 ** http://www.leptonica.com/papers/conn.pdf
1983 kernel=AcquireKernelInfo(
1984 "ThinSE:41; ThinSE:42; ThinSE:43");
1985 if (kernel == (KernelInfo *) NULL)
1987 kernel->type = type;
1988 kernel->next->type = type;
1989 kernel->next->next->type = type;
1990 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1996 { /* Special kernels for general thinning, while preserving connections
1997 ** "Connectivity-Preserving Morphological Image Thransformations"
1998 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1999 ** http://www.leptonica.com/papers/conn.pdf
2001 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2003 ** Note kernels do not specify the origin pixel, allowing them
2004 ** to be used for both thickening and thinning operations.
2006 switch ( (int) args->rho ) {
2007 /* SE for 4-connected thinning */
2008 case 41: /* SE_4_1 */
2009 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
2011 case 42: /* SE_4_2 */
2012 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
2014 case 43: /* SE_4_3 */
2015 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
2017 case 44: /* SE_4_4 */
2018 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
2020 case 45: /* SE_4_5 */
2021 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
2023 case 46: /* SE_4_6 */
2024 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
2026 case 47: /* SE_4_7 */
2027 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
2029 case 48: /* SE_4_8 */
2030 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
2032 case 49: /* SE_4_9 */
2033 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
2035 /* SE for 8-connected thinning - negatives of the above */
2036 case 81: /* SE_8_0 */
2037 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
2039 case 82: /* SE_8_2 */
2040 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
2042 case 83: /* SE_8_3 */
2043 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
2045 case 84: /* SE_8_4 */
2046 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
2048 case 85: /* SE_8_5 */
2049 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2051 case 86: /* SE_8_6 */
2052 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2054 case 87: /* SE_8_7 */
2055 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2057 case 88: /* SE_8_8 */
2058 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2060 case 89: /* SE_8_9 */
2061 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2063 /* Special combined SE kernels */
2064 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2065 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2067 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2068 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2070 case 481: /* SE_48_1 - General Connected Corner Kernel */
2071 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2074 case 482: /* SE_48_2 - General Edge Kernel */
2075 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2078 if (kernel == (KernelInfo *) NULL)
2080 kernel->type = type;
2081 RotateKernelInfo(kernel, args->sigma);
2085 Distance Measuring Kernels
2087 case ChebyshevKernel:
2089 if (args->rho < 1.0)
2090 kernel->width = kernel->height = 3; /* default radius = 1 */
2092 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2093 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2095 kernel->values=(MagickRealType *) MagickAssumeAligned(
2096 AcquireAlignedMemory(kernel->width,kernel->height*
2097 sizeof(*kernel->values)));
2098 if (kernel->values == (MagickRealType *) NULL)
2099 return(DestroyKernelInfo(kernel));
2101 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2102 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2103 kernel->positive_range += ( kernel->values[i] =
2104 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2105 kernel->maximum = kernel->values[0];
2108 case ManhattanKernel:
2110 if (args->rho < 1.0)
2111 kernel->width = kernel->height = 3; /* default radius = 1 */
2113 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2114 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2116 kernel->values=(MagickRealType *) MagickAssumeAligned(
2117 AcquireAlignedMemory(kernel->width,kernel->height*
2118 sizeof(*kernel->values)));
2119 if (kernel->values == (MagickRealType *) NULL)
2120 return(DestroyKernelInfo(kernel));
2122 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2123 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2124 kernel->positive_range += ( kernel->values[i] =
2125 args->sigma*(labs((long) u)+labs((long) v)) );
2126 kernel->maximum = kernel->values[0];
2129 case OctagonalKernel:
2131 if (args->rho < 2.0)
2132 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2134 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2135 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2137 kernel->values=(MagickRealType *) MagickAssumeAligned(
2138 AcquireAlignedMemory(kernel->width,kernel->height*
2139 sizeof(*kernel->values)));
2140 if (kernel->values == (MagickRealType *) NULL)
2141 return(DestroyKernelInfo(kernel));
2143 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2144 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2147 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2148 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2149 kernel->positive_range += kernel->values[i] =
2150 args->sigma*MagickMax(r1,r2);
2152 kernel->maximum = kernel->values[0];
2155 case EuclideanKernel:
2157 if (args->rho < 1.0)
2158 kernel->width = kernel->height = 3; /* default radius = 1 */
2160 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2161 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2163 kernel->values=(MagickRealType *) MagickAssumeAligned(
2164 AcquireAlignedMemory(kernel->width,kernel->height*
2165 sizeof(*kernel->values)));
2166 if (kernel->values == (MagickRealType *) NULL)
2167 return(DestroyKernelInfo(kernel));
2169 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2170 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2171 kernel->positive_range += ( kernel->values[i] =
2172 args->sigma*sqrt((double)(u*u+v*v)) );
2173 kernel->maximum = kernel->values[0];
2178 /* No-Op Kernel - Basically just a single pixel on its own */
2179 kernel=ParseKernelArray("1:1");
2180 if (kernel == (KernelInfo *) NULL)
2182 kernel->type = UndefinedKernel;
2191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2195 % C l o n e K e r n e l I n f o %
2199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2201 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2202 % can be modified without effecting the original. The cloned kernel should
2203 % be destroyed using DestoryKernelInfo() when no longer needed.
2205 % The format of the CloneKernelInfo method is:
2207 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2209 % A description of each parameter follows:
2211 % o kernel: the Morphology/Convolution kernel to be cloned
2214 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2222 assert(kernel != (KernelInfo *) NULL);
2223 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2224 if (new_kernel == (KernelInfo *) NULL)
2226 *new_kernel=(*kernel); /* copy values in structure */
2228 /* replace the values with a copy of the values */
2229 new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2230 AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2231 if (new_kernel->values == (MagickRealType *) NULL)
2232 return(DestroyKernelInfo(new_kernel));
2233 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2234 new_kernel->values[i]=kernel->values[i];
2236 /* Also clone the next kernel in the kernel list */
2237 if ( kernel->next != (KernelInfo *) NULL ) {
2238 new_kernel->next = CloneKernelInfo(kernel->next);
2239 if ( new_kernel->next == (KernelInfo *) NULL )
2240 return(DestroyKernelInfo(new_kernel));
2247 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2251 % D e s t r o y K e r n e l I n f o %
2255 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2257 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2260 % The format of the DestroyKernelInfo method is:
2262 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2264 % A description of each parameter follows:
2266 % o kernel: the Morphology/Convolution kernel to be destroyed
2269 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2271 assert(kernel != (KernelInfo *) NULL);
2272 if ( kernel->next != (KernelInfo *) NULL )
2273 kernel->next=DestroyKernelInfo(kernel->next);
2274 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2275 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2284 + E x p a n d M i r r o r K e r n e l I n f o %
2288 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2290 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2291 % sequence of 90-degree rotated kernels but providing a reflected 180
2292 % rotatation, before the -/+ 90-degree rotations.
2294 % This special rotation order produces a better, more symetrical thinning of
2297 % The format of the ExpandMirrorKernelInfo method is:
2299 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2301 % A description of each parameter follows:
2303 % o kernel: the Morphology/Convolution kernel
2305 % This function is only internel to this module, as it is not finalized,
2306 % especially with regard to non-orthogonal angles, and rotation of larger
2311 static void FlopKernelInfo(KernelInfo *kernel)
2312 { /* Do a Flop by reversing each row. */
2320 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2321 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2322 t=k[x], k[x]=k[r], k[r]=t;
2324 kernel->x = kernel->width - kernel->x - 1;
2325 angle = fmod(angle+180.0, 360.0);
2329 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2337 clone = CloneKernelInfo(last);
2338 RotateKernelInfo(clone, 180); /* flip */
2339 LastKernelInfo(last)->next = clone;
2342 clone = CloneKernelInfo(last);
2343 RotateKernelInfo(clone, 90); /* transpose */
2344 LastKernelInfo(last)->next = clone;
2347 clone = CloneKernelInfo(last);
2348 RotateKernelInfo(clone, 180); /* flop */
2349 LastKernelInfo(last)->next = clone;
2355 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2359 + E x p a n d R o t a t e K e r n e l I n f o %
2363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2365 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2366 % incrementally by the angle given, until the kernel repeats.
2368 % WARNING: 45 degree rotations only works for 3x3 kernels.
2369 % While 90 degree roatations only works for linear and square kernels
2371 % The format of the ExpandRotateKernelInfo method is:
2373 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2375 % A description of each parameter follows:
2377 % o kernel: the Morphology/Convolution kernel
2379 % o angle: angle to rotate in degrees
2381 % This function is only internel to this module, as it is not finalized,
2382 % especially with regard to non-orthogonal angles, and rotation of larger
2386 /* Internal Routine - Return true if two kernels are the same */
2387 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2388 const KernelInfo *kernel2)
2393 /* check size and origin location */
2394 if ( kernel1->width != kernel2->width
2395 || kernel1->height != kernel2->height
2396 || kernel1->x != kernel2->x
2397 || kernel1->y != kernel2->y )
2400 /* check actual kernel values */
2401 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2402 /* Test for Nan equivalence */
2403 if ( IsNaN(kernel1->values[i]) && !IsNaN(kernel2->values[i]) )
2405 if ( IsNaN(kernel2->values[i]) && !IsNaN(kernel1->values[i]) )
2407 /* Test actual values are equivalent */
2408 if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2415 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2423 clone = CloneKernelInfo(last);
2424 RotateKernelInfo(clone, angle);
2425 if ( SameKernelInfo(kernel, clone) == MagickTrue )
2427 LastKernelInfo(last)->next = clone;
2430 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2439 + C a l c M e t a K e r n a l I n f o %
2443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2445 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2446 % using the kernel values. This should only ne used if it is not possible to
2447 % calculate that meta-data in some easier way.
2449 % It is important that the meta-data is correct before ScaleKernelInfo() is
2450 % used to perform kernel normalization.
2452 % The format of the CalcKernelMetaData method is:
2454 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2456 % A description of each parameter follows:
2458 % o kernel: the Morphology/Convolution kernel to modify
2460 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2461 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2462 % however is not true for flat-shaped morphological kernels.
2464 % WARNING: Only the specific kernel pointed to is modified, not a list of
2467 % This is an internal function and not expected to be useful outside this
2468 % module. This could change however.
2470 static void CalcKernelMetaData(KernelInfo *kernel)
2475 kernel->minimum = kernel->maximum = 0.0;
2476 kernel->negative_range = kernel->positive_range = 0.0;
2477 for (i=0; i < (kernel->width*kernel->height); i++)
2479 if ( fabs(kernel->values[i]) < MagickEpsilon )
2480 kernel->values[i] = 0.0;
2481 ( kernel->values[i] < 0)
2482 ? ( kernel->negative_range += kernel->values[i] )
2483 : ( kernel->positive_range += kernel->values[i] );
2484 Minimize(kernel->minimum, kernel->values[i]);
2485 Maximize(kernel->maximum, kernel->values[i]);
2492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2496 % M o r p h o l o g y A p p l y %
2500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2502 % MorphologyApply() applies a morphological method, multiple times using
2503 % a list of multiple kernels. This is the method that should be called by
2504 % other 'operators' that internally use morphology operations as part of
2507 % It is basically equivalent to as MorphologyImage() (see below) but
2508 % without any user controls. This allows internel programs to use this
2509 % function, to actually perform a specific task without possible interference
2510 % by any API user supplied settings.
2512 % It is MorphologyImage() task to extract any such user controls, and
2513 % pass them to this function for processing.
2515 % More specifically all given kernels should already be scaled, normalised,
2516 % and blended appropriatally before being parred to this routine. The
2517 % appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2519 % The format of the MorphologyApply method is:
2521 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2522 % const ssize_t iterations,const KernelInfo *kernel,
2523 % const CompositeMethod compose,const double bias,
2524 % ExceptionInfo *exception)
2526 % A description of each parameter follows:
2528 % o image: the source image
2530 % o method: the morphology method to be applied.
2532 % o iterations: apply the operation this many times (or no change).
2533 % A value of -1 means loop until no change found.
2534 % How this is applied may depend on the morphology method.
2535 % Typically this is a value of 1.
2537 % o channel: the channel type.
2539 % o kernel: An array of double representing the morphology kernel.
2541 % o compose: How to handle or merge multi-kernel results.
2542 % If 'UndefinedCompositeOp' use default for the Morphology method.
2543 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2544 % Otherwise merge the results using the compose method given.
2546 % o bias: Convolution Output Bias.
2548 % o exception: return any errors or warnings in this structure.
2552 /* Apply a Morphology Primative to an image using the given kernel.
2553 ** Two pre-created images must be provided, and no image is created.
2554 ** It returns the number of pixels that changed between the images
2555 ** for result convergence determination.
2557 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2558 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2559 ExceptionInfo *exception)
2561 #define MorphologyTag "Morphology/Image"
2580 assert(image != (Image *) NULL);
2581 assert(image->signature == MagickSignature);
2582 assert(morphology_image != (Image *) NULL);
2583 assert(morphology_image->signature == MagickSignature);
2584 assert(kernel != (KernelInfo *) NULL);
2585 assert(kernel->signature == MagickSignature);
2586 assert(exception != (ExceptionInfo *) NULL);
2587 assert(exception->signature == MagickSignature);
2593 image_view=AcquireVirtualCacheView(image,exception);
2594 morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2595 virt_width=image->columns+kernel->width-1;
2597 /* Some methods (including convolve) needs use a reflected kernel.
2598 * Adjust 'origin' offsets to loop though kernel as a reflection.
2603 case ConvolveMorphology:
2604 case DilateMorphology:
2605 case DilateIntensityMorphology:
2606 case IterativeDistanceMorphology:
2607 /* kernel needs to used with reflection about origin */
2608 offx = (ssize_t) kernel->width-offx-1;
2609 offy = (ssize_t) kernel->height-offy-1;
2611 case ErodeMorphology:
2612 case ErodeIntensityMorphology:
2613 case HitAndMissMorphology:
2614 case ThinningMorphology:
2615 case ThickenMorphology:
2616 /* kernel is used as is, without reflection */
2619 assert("Not a Primitive Morphology Method" != (char *) NULL);
2623 if (method == ConvolveMorphology && kernel->width == 1)
2629 Anthony Thyssen, 14 June 2010
2631 Special handling (for speed) of vertical (blur) kernels. This
2632 performs its handling in columns rather than in rows. This is
2633 only done for convolve as it is the only method that generates very
2634 large 1-D vertical kernels (such as a 'BlurKernel')
2636 Timing tests (on single CPU laptop). Using a vertical 1-d Blue with
2637 normal row-by-row (below):
2638 time convert logo: -morphology Convolve Blur:0x10+90 null: 0.807u
2639 Using this column method
2640 time convert logo: -morphology Convolve Blur:0x10+90 null: 0.620u
2642 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2643 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2644 magick_threads(image,morphology_image,image->columns,1)
2646 for (x=0; x < (ssize_t) image->columns; x++)
2648 register const Quantum
2660 if (status == MagickFalse)
2662 p=GetCacheViewVirtualPixels(image_view,x,-offy,1,image->rows+
2663 kernel->height-1,exception);
2664 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2665 morphology_image->rows,exception);
2666 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2671 center=(ssize_t) GetPixelChannels(image)*offy;
2672 for (y=0; y < (ssize_t) image->rows; y++)
2677 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2691 register const MagickRealType
2694 register const Quantum
2703 channel=GetPixelChannelChannel(image,i);
2704 traits=GetPixelChannelTraits(image,channel);
2705 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2706 if ((traits == UndefinedPixelTrait) ||
2707 (morphology_traits == UndefinedPixelTrait))
2709 if (((morphology_traits & CopyPixelTrait) != 0) ||
2710 (GetPixelMask(image,p) != 0))
2712 SetPixelChannel(morphology_image,channel,p[center+i],q);
2715 k=(&kernel->values[kernel->height-1]);
2719 if ((morphology_traits & BlendPixelTrait) == 0)
2724 for (v=0; v < (ssize_t) kernel->height; v++)
2726 for (u=0; u < (ssize_t) kernel->width; u++)
2728 if (IsNaN(*k) != MagickFalse)
2730 pixel+=(*k)*pixels[i];
2733 pixels+=GetPixelChannels(image);
2736 gamma=PerceptibleReciprocal(gamma);
2738 if (fabs(pixel-p[center+i]) > MagickEpsilon)
2740 SetPixelChannel(morphology_image,channel,ClampToQuantum(pixel),
2747 for (v=0; v < (ssize_t) kernel->width; v++)
2749 for (u=0; u < (ssize_t) kernel->width; u++)
2751 if (IsNaN(*k) != MagickFalse)
2753 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2754 pixel+=(*k)*alpha*pixels[i];
2757 pixels+=GetPixelChannels(image);
2760 gamma=PerceptibleReciprocal(gamma);
2762 if (fabs(pixel-p[center+i]) > MagickEpsilon)
2764 SetPixelChannel(morphology_image,channel,ClampToQuantum(pixel),q);
2766 p+=GetPixelChannels(image);
2767 q+=GetPixelChannels(morphology_image);
2769 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2771 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2776 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2777 #pragma omp critical (MagickCore_MorphologyImage)
2779 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2780 if (proceed == MagickFalse)
2784 morphology_image->type=image->type;
2785 morphology_view=DestroyCacheView(morphology_view);
2786 image_view=DestroyCacheView(image_view);
2787 return(status ? (ssize_t) changed : 0);
2790 Normal handling of horizontal or rectangular kernels (row by row).
2792 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2793 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2794 magick_threads(image,morphology_image,image->rows,1)
2796 for (y=0; y < (ssize_t) image->rows; y++)
2798 register const Quantum
2810 if (status == MagickFalse)
2812 p=GetCacheViewVirtualPixels(image_view,-offx,y-offy,virt_width,
2813 kernel->height,exception);
2814 q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2816 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2821 /* offset to origin in 'p'. while 'q' points to it directly */
2822 r = GetPixelChannels(image)*virt_width*offy + GetPixelChannels(image)*offx;
2824 for (x=0; x < (ssize_t) image->columns; x++)
2831 register const MagickRealType
2834 register const Quantum
2843 /* Copy input image to the output image for unused channels
2844 * This removes need for 'cloning' a new image every iteration
2846 SetPixelRed(morphology_image,GetPixelRed(image,p+r),q);
2847 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r),q);
2848 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r),q);
2849 if (image->colorspace == CMYKColorspace)
2850 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r),q);
2857 min.black = (double) QuantumRange;
2862 max.black = (double) 0;
2863 /* default result is the original pixel value */
2864 result.red = (double) GetPixelRed(image,p+r);
2865 result.green = (double) GetPixelGreen(image,p+r);
2866 result.blue = (double) GetPixelBlue(image,p+r);
2868 if (image->colorspace == CMYKColorspace)
2869 result.black = (double) GetPixelBlack(image,p+r);
2870 result.alpha=(double) GetPixelAlpha(image,p+r);
2873 case ConvolveMorphology:
2874 /* Set the bias of the weighted average output */
2879 result.black = bias;
2881 case DilateIntensityMorphology:
2882 case ErodeIntensityMorphology:
2883 /* use a boolean flag indicating when first match found */
2884 result.red = 0.0; /* result is not used otherwise */
2891 case ConvolveMorphology:
2892 /* Weighted Average of pixels using reflected kernel
2894 ** NOTE for correct working of this operation for asymetrical
2895 ** kernels, the kernel needs to be applied in its reflected form.
2896 ** That is its values needs to be reversed.
2898 ** Correlation is actually the same as this but without reflecting
2899 ** the kernel, and thus 'lower-level' that Convolution. However
2900 ** as Convolution is the more common method used, and it does not
2901 ** really cost us much in terms of processing to use a reflected
2902 ** kernel, so it is Convolution that is implemented.
2904 ** Correlation will have its kernel reflected before calling
2905 ** this function to do a Convolve.
2907 ** For more details of Correlation vs Convolution see
2908 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2910 k = &kernel->values[ kernel->width*kernel->height-1 ];
2912 if ( (image->channel_mask != DefaultChannels) ||
2913 (image->alpha_trait != BlendPixelTrait) )
2914 { /* No 'Sync' involved.
2915 ** Convolution is simple greyscale channel operation
2917 for (v=0; v < (ssize_t) kernel->height; v++) {
2918 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2919 if ( IsNaN(*k) ) continue;
2921 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2922 result.green += (*k)*
2923 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2924 result.blue += (*k)*
2925 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2926 if (image->colorspace == CMYKColorspace)
2927 result.black += (*k)*
2928 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2929 result.alpha += (*k)*
2930 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2932 k_pixels += virt_width*GetPixelChannels(image);
2934 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2935 SetPixelRed(morphology_image,ClampToQuantum(result.red),
2937 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2938 SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2940 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2941 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2943 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2944 (image->colorspace == CMYKColorspace))
2945 SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2947 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2948 (image->alpha_trait == BlendPixelTrait))
2949 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
2953 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2954 ** Weight the color channels with Alpha Channel so that
2955 ** transparent pixels are not part of the results.
2958 gamma; /* divisor, sum of color alpha weighting */
2961 alpha; /* alpha weighting for colors : alpha */
2964 count; /* alpha valus collected, number kernel values */
2968 for (v=0; v < (ssize_t) kernel->height; v++) {
2969 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2970 if ( IsNaN(*k) ) continue;
2971 alpha=QuantumScale*GetPixelAlpha(image,
2972 k_pixels+u*GetPixelChannels(image));
2973 gamma += alpha; /* normalize alpha weights only */
2974 count++; /* number of alpha values collected */
2975 alpha=alpha*(*k); /* include kernel weighting now */
2976 result.red += alpha*
2977 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2978 result.green += alpha*
2979 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2980 result.blue += alpha*
2981 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2982 if (image->colorspace == CMYKColorspace)
2983 result.black += alpha*
2984 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2985 result.alpha += (*k)*
2986 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2988 k_pixels += virt_width*GetPixelChannels(image);
2990 /* Sync'ed channels, all channels are modified */
2991 gamma=(double)count/(fabs((double) gamma) < MagickEpsilon ? MagickEpsilon : gamma);
2992 SetPixelRed(morphology_image,
2993 ClampToQuantum(gamma*result.red),q);
2994 SetPixelGreen(morphology_image,
2995 ClampToQuantum(gamma*result.green),q);
2996 SetPixelBlue(morphology_image,
2997 ClampToQuantum(gamma*result.blue),q);
2998 if (image->colorspace == CMYKColorspace)
2999 SetPixelBlack(morphology_image,
3000 ClampToQuantum(gamma*result.black),q);
3001 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3005 case ErodeMorphology:
3006 /* Minimum Value within kernel neighbourhood
3008 ** NOTE that the kernel is not reflected for this operation!
3010 ** NOTE: in normal Greyscale Morphology, the kernel value should
3011 ** be added to the real value, this is currently not done, due to
3012 ** the nature of the boolean kernels being used.
3016 for (v=0; v < (ssize_t) kernel->height; v++) {
3017 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3018 if ( IsNaN(*k) || (*k) < 0.5 ) continue;
3019 Minimize(min.red, (double)
3020 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3021 Minimize(min.green, (double)
3022 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3023 Minimize(min.blue, (double)
3024 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3025 Minimize(min.alpha, (double)
3026 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3027 if (image->colorspace == CMYKColorspace)
3028 Minimize(min.black, (double)
3029 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3031 k_pixels += virt_width*GetPixelChannels(image);
3035 case DilateMorphology:
3036 /* Maximum Value within kernel neighbourhood
3038 ** NOTE for correct working of this operation for asymetrical
3039 ** kernels, the kernel needs to be applied in its reflected form.
3040 ** That is its values needs to be reversed.
3042 ** NOTE: in normal Greyscale Morphology, the kernel value should
3043 ** be added to the real value, this is currently not done, due to
3044 ** the nature of the boolean kernels being used.
3047 k = &kernel->values[ kernel->width*kernel->height-1 ];
3049 for (v=0; v < (ssize_t) kernel->height; v++) {
3050 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3051 if ( IsNaN(*k) || (*k) < 0.5 ) continue;
3052 Maximize(max.red, (double)
3053 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3054 Maximize(max.green, (double)
3055 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3056 Maximize(max.blue, (double)
3057 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3058 Maximize(max.alpha, (double)
3059 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3060 if (image->colorspace == CMYKColorspace)
3061 Maximize(max.black, (double)
3062 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3064 k_pixels += virt_width*GetPixelChannels(image);
3068 case HitAndMissMorphology:
3069 case ThinningMorphology:
3070 case ThickenMorphology:
3071 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3073 ** NOTE that the kernel is not reflected for this operation,
3074 ** and consists of both foreground and background pixel
3075 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3076 ** with either Nan or 0.5 values for don't care.
3078 ** Note that this will never produce a meaningless negative
3079 ** result. Such results can cause Thinning/Thicken to not work
3080 ** correctly when used against a greyscale image.
3084 for (v=0; v < (ssize_t) kernel->height; v++) {
3085 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3086 if ( IsNaN(*k) ) continue;
3088 { /* minimim of foreground pixels */
3089 Minimize(min.red, (double)
3090 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3091 Minimize(min.green, (double)
3092 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3093 Minimize(min.blue, (double)
3094 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3095 Minimize(min.alpha,(double)
3096 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3097 if ( image->colorspace == CMYKColorspace)
3098 Minimize(min.black,(double)
3099 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3101 else if ( (*k) < 0.3 )
3102 { /* maximum of background pixels */
3103 Maximize(max.red, (double)
3104 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3105 Maximize(max.green, (double)
3106 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3107 Maximize(max.blue, (double)
3108 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3109 Maximize(max.alpha,(double)
3110 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3111 if (image->colorspace == CMYKColorspace)
3112 Maximize(max.black, (double)
3113 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3116 k_pixels += virt_width*GetPixelChannels(image);
3118 /* Pattern Match if difference is positive */
3119 min.red -= max.red; Maximize( min.red, 0.0 );
3120 min.green -= max.green; Maximize( min.green, 0.0 );
3121 min.blue -= max.blue; Maximize( min.blue, 0.0 );
3122 min.black -= max.black; Maximize( min.black, 0.0 );
3123 min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3126 case ErodeIntensityMorphology:
3127 /* Select Pixel with Minimum Intensity within kernel neighbourhood
3129 ** WARNING: the intensity test fails for CMYK and does not
3130 ** take into account the moderating effect of the alpha channel
3131 ** on the intensity.
3133 ** NOTE that the kernel is not reflected for this operation!
3137 for (v=0; v < (ssize_t) kernel->height; v++) {
3138 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3139 if ( IsNaN(*k) || (*k) < 0.5 ) continue;
3140 if ( result.red == 0.0 ||
3141 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) < GetPixelIntensity(morphology_image,q) ) {
3142 /* copy the whole pixel - no channel selection */
3143 SetPixelRed(morphology_image,GetPixelRed(image,
3144 k_pixels+u*GetPixelChannels(image)),q);
3145 SetPixelGreen(morphology_image,GetPixelGreen(image,
3146 k_pixels+u*GetPixelChannels(image)),q);
3147 SetPixelBlue(morphology_image,GetPixelBlue(image,
3148 k_pixels+u*GetPixelChannels(image)),q);
3149 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3150 k_pixels+u*GetPixelChannels(image)),q);
3151 if ( result.red > 0.0 ) changed++;
3155 k_pixels += virt_width*GetPixelChannels(image);
3159 case DilateIntensityMorphology:
3160 /* Select Pixel with Maximum Intensity within kernel neighbourhood
3162 ** WARNING: the intensity test fails for CMYK and does not
3163 ** take into account the moderating effect of the alpha channel
3164 ** on the intensity (yet).
3166 ** NOTE for correct working of this operation for asymetrical
3167 ** kernels, the kernel needs to be applied in its reflected form.
3168 ** That is its values needs to be reversed.
3170 k = &kernel->values[ kernel->width*kernel->height-1 ];
3172 for (v=0; v < (ssize_t) kernel->height; v++) {
3173 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3174 if ( IsNaN(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3175 if ( result.red == 0.0 ||
3176 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) > GetPixelIntensity(morphology_image,q) ) {
3177 /* copy the whole pixel - no channel selection */
3178 SetPixelRed(morphology_image,GetPixelRed(image,
3179 k_pixels+u*GetPixelChannels(image)),q);
3180 SetPixelGreen(morphology_image,GetPixelGreen(image,
3181 k_pixels+u*GetPixelChannels(image)),q);
3182 SetPixelBlue(morphology_image,GetPixelBlue(image,
3183 k_pixels+u*GetPixelChannels(image)),q);
3184 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3185 k_pixels+u*GetPixelChannels(image)),q);
3186 if ( result.red > 0.0 ) changed++;
3190 k_pixels += virt_width*GetPixelChannels(image);
3194 case IterativeDistanceMorphology:
3195 /* Work out an iterative distance from black edge of a white image
3196 ** shape. Essentually white values are decreased to the smallest
3197 ** 'distance from edge' it can find.
3199 ** It works by adding kernel values to the neighbourhood, and and
3200 ** select the minimum value found. The kernel is rotated before
3201 ** use, so kernel distances match resulting distances, when a user
3202 ** provided asymmetric kernel is applied.
3205 ** This code is almost identical to True GrayScale Morphology But
3208 ** GreyDilate Kernel values added, maximum value found Kernel is
3209 ** rotated before use.
3211 ** GrayErode: Kernel values subtracted and minimum value found No
3212 ** kernel rotation used.
3214 ** Note the the Iterative Distance method is essentially a
3215 ** GrayErode, but with negative kernel values, and kernel
3216 ** rotation applied.
3218 k = &kernel->values[ kernel->width*kernel->height-1 ];
3220 for (v=0; v < (ssize_t) kernel->height; v++) {
3221 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3222 if ( IsNaN(*k) ) continue;
3223 Minimize(result.red, (*k)+(double)
3224 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3225 Minimize(result.green, (*k)+(double)
3226 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3227 Minimize(result.blue, (*k)+(double)
3228 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3229 Minimize(result.alpha, (*k)+(double)
3230 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3231 if ( image->colorspace == CMYKColorspace)
3232 Maximize(result.black, (*k)+(double)
3233 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3235 k_pixels += virt_width*GetPixelChannels(image);
3239 case UndefinedMorphology:
3241 break; /* Do nothing */
3243 /* Final mathematics of results (combine with original image?)
3245 ** NOTE: Difference Morphology operators Edge* and *Hat could also
3246 ** be done here but works better with iteration as a image difference
3247 ** in the controling function (below). Thicken and Thinning however
3248 ** should be done here so thay can be iterated correctly.
3251 case HitAndMissMorphology:
3252 case ErodeMorphology:
3253 result = min; /* minimum of neighbourhood */
3255 case DilateMorphology:
3256 result = max; /* maximum of neighbourhood */
3258 case ThinningMorphology:
3259 /* subtract pattern match from original */
3260 result.red -= min.red;
3261 result.green -= min.green;
3262 result.blue -= min.blue;
3263 result.black -= min.black;
3264 result.alpha -= min.alpha;
3266 case ThickenMorphology:
3267 /* Add the pattern matchs to the original */
3268 result.red += min.red;
3269 result.green += min.green;
3270 result.blue += min.blue;
3271 result.black += min.black;
3272 result.alpha += min.alpha;
3275 /* result directly calculated or assigned */
3278 /* Assign the resulting pixel values - Clamping Result */
3280 case UndefinedMorphology:
3281 case ConvolveMorphology:
3282 case DilateIntensityMorphology:
3283 case ErodeIntensityMorphology:
3284 break; /* full pixel was directly assigned - not a channel method */
3286 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3287 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
3288 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3289 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
3290 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3291 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
3292 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3293 (image->colorspace == CMYKColorspace))
3294 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
3295 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3296 (image->alpha_trait == BlendPixelTrait))
3297 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3300 /* Count up changed pixels */
3301 if ((GetPixelRed(image,p+r) != GetPixelRed(morphology_image,q)) ||
3302 (GetPixelGreen(image,p+r) != GetPixelGreen(morphology_image,q)) ||
3303 (GetPixelBlue(image,p+r) != GetPixelBlue(morphology_image,q)) ||
3304 (GetPixelAlpha(image,p+r) != GetPixelAlpha(morphology_image,q)) ||
3305 ((image->colorspace == CMYKColorspace) &&
3306 (GetPixelBlack(image,p+r) != GetPixelBlack(morphology_image,q))))
3307 changed++; /* The pixel was changed in some way! */
3308 p+=GetPixelChannels(image);
3309 q+=GetPixelChannels(morphology_image);
3311 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3313 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3318 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3319 #pragma omp critical (MagickCore_MorphologyImage)
3321 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3322 if (proceed == MagickFalse)
3326 morphology_view=DestroyCacheView(morphology_view);
3327 image_view=DestroyCacheView(image_view);
3328 return(status ? (ssize_t)changed : -1);
3331 /* This is almost identical to the MorphologyPrimative() function above,
3332 ** but will apply the primitive directly to the actual image using two
3333 ** passes, once in each direction, with the results of the previous (and
3334 ** current) row being re-used.
3336 ** That is after each row is 'Sync'ed' into the image, the next row will
3337 ** make use of those values as part of the calculation of the next row.
3338 ** It then repeats, but going in the oppisite (bottom-up) direction.
3340 ** Because of this 're-use of results' this function can not make use
3341 ** of multi-threaded, parellel processing.
3343 static ssize_t MorphologyPrimitiveDirect(Image *image,
3344 const MorphologyMethod method,const KernelInfo *kernel,
3345 ExceptionInfo *exception)
3368 assert(image != (Image *) NULL);
3369 assert(image->signature == MagickSignature);
3370 assert(kernel != (KernelInfo *) NULL);
3371 assert(kernel->signature == MagickSignature);
3372 assert(exception != (ExceptionInfo *) NULL);
3373 assert(exception->signature == MagickSignature);
3375 /* Some methods (including convolve) needs use a reflected kernel.
3376 * Adjust 'origin' offsets to loop though kernel as a reflection.
3381 case DistanceMorphology:
3382 case VoronoiMorphology:
3383 /* kernel needs to used with reflection about origin */
3384 offx = (ssize_t) kernel->width-offx-1;
3385 offy = (ssize_t) kernel->height-offy-1;
3388 case ?????Morphology:
3389 /* kernel is used as is, without reflection */
3393 assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
3397 /* DO NOT THREAD THIS CODE! */
3398 /* two views into same image (virtual, and actual) */
3399 virt_view=AcquireVirtualCacheView(image,exception);
3400 auth_view=AcquireAuthenticCacheView(image,exception);
3401 virt_width=image->columns+kernel->width-1;
3403 for (y=0; y < (ssize_t) image->rows; y++)
3405 register const Quantum
3417 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3418 ** we read using virtual to get virtual pixel handling, but write back
3419 ** into the same image.
3421 ** Only top half of kernel is processed as we do a single pass downward
3422 ** through the image iterating the distance function as we go.
3424 if (status == MagickFalse)
3426 p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
3428 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3430 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3432 if (status == MagickFalse)
3435 /* offset to origin in 'p'. while 'q' points to it directly */
3436 r = (ssize_t) GetPixelChannels(image)*virt_width*offy + GetPixelChannels(image)*offx;
3438 for (x=0; x < (ssize_t) image->columns; x++)
3443 register const MagickRealType
3446 register const Quantum
3455 /* Starting Defaults */
3456 GetPixelInfo(image,&result);
3457 GetPixelInfoPixel(image,q,&result);
3458 if ( method != VoronoiMorphology )
3459 result.alpha = QuantumRange - result.alpha;
3462 case DistanceMorphology:
3463 /* Add kernel Value and select the minimum value found. */
3464 k = &kernel->values[ kernel->width*kernel->height-1 ];
3466 for (v=0; v <= (ssize_t) offy; v++) {
3467 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3468 if ( IsNaN(*k) ) continue;
3469 Minimize(result.red, (*k)+
3470 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3471 Minimize(result.green, (*k)+
3472 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3473 Minimize(result.blue, (*k)+
3474 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3475 if (image->colorspace == CMYKColorspace)
3476 Minimize(result.black,(*k)+
3477 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3478 Minimize(result.alpha, (*k)+
3479 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3481 k_pixels += virt_width*GetPixelChannels(image);
3483 /* repeat with the just processed pixels of this row */
3484 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3485 k_pixels = q-offx*GetPixelChannels(image);
3486 for (u=0; u < (ssize_t) offx; u++, k--) {
3487 if ( x+u-offx < 0 ) continue; /* off the edge! */
3488 if ( IsNaN(*k) ) continue;
3489 Minimize(result.red, (*k)+
3490 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3491 Minimize(result.green, (*k)+
3492 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3493 Minimize(result.blue, (*k)+
3494 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3495 if (image->colorspace == CMYKColorspace)
3496 Minimize(result.black,(*k)+
3497 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3498 Minimize(result.alpha,(*k)+
3499 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3502 case VoronoiMorphology:
3503 /* Apply Distance to 'Matte' channel, while coping the color
3504 ** values of the closest pixel.
3506 ** This is experimental, and realy the 'alpha' component should
3507 ** be completely separate 'masking' channel so that alpha can
3508 ** also be used as part of the results.
3510 k = &kernel->values[ kernel->width*kernel->height-1 ];
3512 for (v=0; v <= (ssize_t) offy; v++) {
3513 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3514 if ( IsNaN(*k) ) continue;
3515 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3517 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3522 k_pixels += virt_width*GetPixelChannels(image);
3524 /* repeat with the just processed pixels of this row */
3525 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3526 k_pixels = q-offx*GetPixelChannels(image);
3527 for (u=0; u < (ssize_t) offx; u++, k--) {
3528 if ( x+u-offx < 0 ) continue; /* off the edge! */
3529 if ( IsNaN(*k) ) continue;
3530 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3532 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3539 /* result directly calculated or assigned */
3542 /* Assign the resulting pixel values - Clamping Result */
3544 case VoronoiMorphology:
3545 SetPixelInfoPixel(image,&result,q);
3548 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3549 SetPixelRed(image,ClampToQuantum(result.red),q);
3550 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3551 SetPixelGreen(image,ClampToQuantum(result.green),q);
3552 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3553 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3554 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3555 (image->colorspace == CMYKColorspace))
3556 SetPixelBlack(image,ClampToQuantum(result.black),q);
3557 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3558 (image->alpha_trait == BlendPixelTrait))
3559 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3562 /* Count up changed pixels */
3563 if ((GetPixelRed(image,p+r) != GetPixelRed(image,q)) ||
3564 (GetPixelGreen(image,p+r) != GetPixelGreen(image,q)) ||
3565 (GetPixelBlue(image,p+r) != GetPixelBlue(image,q)) ||
3566 (GetPixelAlpha(image,p+r) != GetPixelAlpha(image,q)) ||
3567 ((image->colorspace == CMYKColorspace) &&
3568 (GetPixelBlack(image,p+r) != GetPixelBlack(image,q))))
3569 changed++; /* The pixel was changed in some way! */
3571 p+=GetPixelChannels(image); /* increment pixel buffers */
3572 q+=GetPixelChannels(image);
3575 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3577 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3578 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3584 /* Do the reversed pass through the image */
3585 for (y=(ssize_t)image->rows-1; y >= 0; y--)
3587 register const Quantum
3599 if (status == MagickFalse)
3601 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3602 ** we read using virtual to get virtual pixel handling, but write back
3603 ** into the same image.
3605 ** Only the bottom half of the kernel will be processes as we
3608 p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
3609 kernel->y+1,exception);
3610 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3612 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3614 if (status == MagickFalse)
3617 /* adjust positions to end of row */
3618 p += (image->columns-1)*GetPixelChannels(image);
3619 q += (image->columns-1)*GetPixelChannels(image);
3621 /* offset to origin in 'p'. while 'q' points to it directly */
3622 r = GetPixelChannels(image)*offx;
3624 for (x=(ssize_t)image->columns-1; x >= 0; x--)
3629 register const MagickRealType
3632 register const Quantum
3641 /* Default - previously modified pixel */
3642 GetPixelInfo(image,&result);
3643 GetPixelInfoPixel(image,q,&result);
3644 if ( method != VoronoiMorphology )
3645 result.alpha = QuantumRange - result.alpha;
3648 case DistanceMorphology:
3649 /* Add kernel Value and select the minimum value found. */
3650 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3652 for (v=offy; v < (ssize_t) kernel->height; v++) {
3653 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3654 if ( IsNaN(*k) ) continue;
3655 Minimize(result.red, (*k)+
3656 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3657 Minimize(result.green, (*k)+
3658 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3659 Minimize(result.blue, (*k)+
3660 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3661 if ( image->colorspace == CMYKColorspace)
3662 Minimize(result.black,(*k)+
3663 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3664 Minimize(result.alpha, (*k)+
3665 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3667 k_pixels += virt_width*GetPixelChannels(image);
3669 /* repeat with the just processed pixels of this row */
3670 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3671 k_pixels = q-offx*GetPixelChannels(image);
3672 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3673 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3674 if ( IsNaN(*k) ) continue;
3675 Minimize(result.red, (*k)+
3676 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3677 Minimize(result.green, (*k)+
3678 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3679 Minimize(result.blue, (*k)+
3680 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3681 if ( image->colorspace == CMYKColorspace)
3682 Minimize(result.black, (*k)+
3683 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3684 Minimize(result.alpha, (*k)+
3685 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3688 case VoronoiMorphology:
3689 /* Apply Distance to 'Matte' channel, coping the closest color.
3691 ** This is experimental, and realy the 'alpha' component should
3692 ** be completely separate 'masking' channel.
3694 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3696 for (v=offy; v < (ssize_t) kernel->height; v++) {
3697 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3698 if ( IsNaN(*k) ) continue;
3699 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3701 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3706 k_pixels += virt_width*GetPixelChannels(image);
3708 /* repeat with the just processed pixels of this row */
3709 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3710 k_pixels = q-offx*GetPixelChannels(image);
3711 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3712 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3713 if ( IsNaN(*k) ) continue;
3714 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3716 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3723 /* result directly calculated or assigned */
3726 /* Assign the resulting pixel values - Clamping Result */
3728 case VoronoiMorphology:
3729 SetPixelInfoPixel(image,&result,q);
3732 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3733 SetPixelRed(image,ClampToQuantum(result.red),q);
3734 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3735 SetPixelGreen(image,ClampToQuantum(result.green),q);
3736 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3737 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3738 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3739 (image->colorspace == CMYKColorspace))
3740 SetPixelBlack(image,ClampToQuantum(result.black),q);
3741 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3742 (image->alpha_trait == BlendPixelTrait))
3743 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3746 /* Count up changed pixels */
3747 if ( (GetPixelRed(image,p+r) != GetPixelRed(image,q))
3748 || (GetPixelGreen(image,p+r) != GetPixelGreen(image,q))
3749 || (GetPixelBlue(image,p+r) != GetPixelBlue(image,q))
3750 || (GetPixelAlpha(image,p+r) != GetPixelAlpha(image,q))
3751 || ((image->colorspace == CMYKColorspace) &&
3752 (GetPixelBlack(image,p+r) != GetPixelBlack(image,q))))
3753 changed++; /* The pixel was changed in some way! */
3755 p-=GetPixelChannels(image); /* go backward through pixel buffers */
3756 q-=GetPixelChannels(image);
3758 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3760 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3761 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3767 auth_view=DestroyCacheView(auth_view);
3768 virt_view=DestroyCacheView(virt_view);
3769 return(status ? (ssize_t) changed : -1);
3772 /* Apply a Morphology by calling one of the above low level primitive
3773 ** application functions. This function handles any iteration loops,
3774 ** composition or re-iteration of results, and compound morphology methods
3775 ** that is based on multiple low-level (staged) morphology methods.
3777 ** Basically this provides the complex glue between the requested morphology
3778 ** method and raw low-level implementation (above).
3780 MagickPrivate Image *MorphologyApply(const Image *image,
3781 const MorphologyMethod method, const ssize_t iterations,
3782 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3783 ExceptionInfo *exception)
3789 *curr_image, /* Image we are working with or iterating */
3790 *work_image, /* secondary image for primitive iteration */
3791 *save_image, /* saved image - for 'edge' method only */
3792 *rslt_image; /* resultant image - after multi-kernel handling */
3795 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3796 *norm_kernel, /* the current normal un-reflected kernel */
3797 *rflt_kernel, /* the current reflected kernel (if needed) */
3798 *this_kernel; /* the kernel being applied */
3801 primitive; /* the current morphology primitive being applied */
3804 rslt_compose; /* multi-kernel compose method for results to use */
3807 special, /* do we use a direct modify function? */
3808 verbose; /* verbose output of results */
3811 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3812 method_limit, /* maximum number of compound method iterations */
3813 kernel_number, /* Loop 2: the kernel number being applied */
3814 stage_loop, /* Loop 3: primitive loop for compound morphology */
3815 stage_limit, /* how many primitives are in this compound */
3816 kernel_loop, /* Loop 4: iterate the kernel over image */
3817 kernel_limit, /* number of times to iterate kernel */
3818 count, /* total count of primitive steps applied */
3819 kernel_changed, /* total count of changed using iterated kernel */
3820 method_changed; /* total count of changed over method iteration */
3823 changed; /* number pixels changed by last primitive operation */
3828 assert(image != (Image *) NULL);
3829 assert(image->signature == MagickSignature);
3830 assert(kernel != (KernelInfo *) NULL);
3831 assert(kernel->signature == MagickSignature);
3832 assert(exception != (ExceptionInfo *) NULL);
3833 assert(exception->signature == MagickSignature);
3835 count = 0; /* number of low-level morphology primitives performed */
3836 if ( iterations == 0 )
3837 return((Image *)NULL); /* null operation - nothing to do! */
3839 kernel_limit = (size_t) iterations;
3840 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3841 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3843 verbose = IsStringTrue(GetImageArtifact(image,"verbose"));
3845 /* initialise for cleanup */
3846 curr_image = (Image *) image;
3847 curr_compose = image->compose;
3848 (void) curr_compose;
3849 work_image = save_image = rslt_image = (Image *) NULL;
3850 reflected_kernel = (KernelInfo *) NULL;
3852 /* Initialize specific methods
3853 * + which loop should use the given iteratations
3854 * + how many primitives make up the compound morphology
3855 * + multi-kernel compose method to use (by default)
3857 method_limit = 1; /* just do method once, unless otherwise set */
3858 stage_limit = 1; /* assume method is not a compound */
3859 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3860 rslt_compose = compose; /* and we are composing multi-kernels as given */
3862 case SmoothMorphology: /* 4 primitive compound morphology */
3865 case OpenMorphology: /* 2 primitive compound morphology */
3866 case OpenIntensityMorphology:
3867 case TopHatMorphology:
3868 case CloseMorphology:
3869 case CloseIntensityMorphology:
3870 case BottomHatMorphology:
3871 case EdgeMorphology:
3874 case HitAndMissMorphology:
3875 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3877 case ThinningMorphology:
3878 case ThickenMorphology:
3879 method_limit = kernel_limit; /* iterate the whole method */
3880 kernel_limit = 1; /* do not do kernel iteration */
3882 case DistanceMorphology:
3883 case VoronoiMorphology:
3884 special = MagickTrue; /* use special direct primative */
3890 /* Apply special methods with special requirments
3891 ** For example, single run only, or post-processing requirements
3893 if ( special == MagickTrue )
3895 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3896 if (rslt_image == (Image *) NULL)
3898 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3901 changed = MorphologyPrimitiveDirect(rslt_image, method,
3904 if ( IfMagickTrue(verbose) )
3905 (void) (void) FormatLocaleFile(stderr,
3906 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3907 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3908 1.0,0.0,1.0, (double) changed);
3913 if ( method == VoronoiMorphology ) {
3914 /* Preserve the alpha channel of input image - but turned off */
3915 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3917 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3918 MagickTrue,0,0,exception);
3919 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3925 /* Handle user (caller) specified multi-kernel composition method */
3926 if ( compose != UndefinedCompositeOp )
3927 rslt_compose = compose; /* override default composition for method */
3928 if ( rslt_compose == UndefinedCompositeOp )
3929 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3931 /* Some methods require a reflected kernel to use with primitives.
3932 * Create the reflected kernel for those methods. */
3934 case CorrelateMorphology:
3935 case CloseMorphology:
3936 case CloseIntensityMorphology:
3937 case BottomHatMorphology:
3938 case SmoothMorphology:
3939 reflected_kernel = CloneKernelInfo(kernel);
3940 if (reflected_kernel == (KernelInfo *) NULL)
3942 RotateKernelInfo(reflected_kernel,180);
3948 /* Loops around more primitive morpholgy methods
3949 ** erose, dilate, open, close, smooth, edge, etc...
3951 /* Loop 1: iterate the compound method */
3954 while ( method_loop < method_limit && method_changed > 0 ) {
3958 /* Loop 2: iterate over each kernel in a multi-kernel list */
3959 norm_kernel = (KernelInfo *) kernel;
3960 this_kernel = (KernelInfo *) kernel;
3961 rflt_kernel = reflected_kernel;
3964 while ( norm_kernel != NULL ) {
3966 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3967 stage_loop = 0; /* the compound morphology stage number */
3968 while ( stage_loop < stage_limit ) {
3969 stage_loop++; /* The stage of the compound morphology */
3971 /* Select primitive morphology for this stage of compound method */
3972 this_kernel = norm_kernel; /* default use unreflected kernel */
3973 primitive = method; /* Assume method is a primitive */
3975 case ErodeMorphology: /* just erode */
3976 case EdgeInMorphology: /* erode and image difference */
3977 primitive = ErodeMorphology;
3979 case DilateMorphology: /* just dilate */
3980 case EdgeOutMorphology: /* dilate and image difference */
3981 primitive = DilateMorphology;
3983 case OpenMorphology: /* erode then dialate */
3984 case TopHatMorphology: /* open and image difference */
3985 primitive = ErodeMorphology;
3986 if ( stage_loop == 2 )
3987 primitive = DilateMorphology;
3989 case OpenIntensityMorphology:
3990 primitive = ErodeIntensityMorphology;
3991 if ( stage_loop == 2 )
3992 primitive = DilateIntensityMorphology;
3994 case CloseMorphology: /* dilate, then erode */
3995 case BottomHatMorphology: /* close and image difference */
3996 this_kernel = rflt_kernel; /* use the reflected kernel */
3997 primitive = DilateMorphology;
3998 if ( stage_loop == 2 )
3999 primitive = ErodeMorphology;
4001 case CloseIntensityMorphology:
4002 this_kernel = rflt_kernel; /* use the reflected kernel */
4003 primitive = DilateIntensityMorphology;
4004 if ( stage_loop == 2 )
4005 primitive = ErodeIntensityMorphology;
4007 case SmoothMorphology: /* open, close */
4008 switch ( stage_loop ) {
4009 case 1: /* start an open method, which starts with Erode */
4010 primitive = ErodeMorphology;
4012 case 2: /* now Dilate the Erode */
4013 primitive = DilateMorphology;
4015 case 3: /* Reflect kernel a close */
4016 this_kernel = rflt_kernel; /* use the reflected kernel */
4017 primitive = DilateMorphology;
4019 case 4: /* Finish the Close */
4020 this_kernel = rflt_kernel; /* use the reflected kernel */
4021 primitive = ErodeMorphology;
4025 case EdgeMorphology: /* dilate and erode difference */
4026 primitive = DilateMorphology;
4027 if ( stage_loop == 2 ) {
4028 save_image = curr_image; /* save the image difference */
4029 curr_image = (Image *) image;
4030 primitive = ErodeMorphology;
4033 case CorrelateMorphology:
4034 /* A Correlation is a Convolution with a reflected kernel.
4035 ** However a Convolution is a weighted sum using a reflected
4036 ** kernel. It may seem stange to convert a Correlation into a
4037 ** Convolution as the Correlation is the simplier method, but
4038 ** Convolution is much more commonly used, and it makes sense to
4039 ** implement it directly so as to avoid the need to duplicate the
4040 ** kernel when it is not required (which is typically the
4043 this_kernel = rflt_kernel; /* use the reflected kernel */
4044 primitive = ConvolveMorphology;
4049 assert( this_kernel != (KernelInfo *) NULL );
4051 /* Extra information for debugging compound operations */
4052 if ( IfMagickTrue(verbose) ) {
4053 if ( stage_limit > 1 )
4054 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
4055 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
4056 method_loop,(double) stage_loop);
4057 else if ( primitive != method )
4058 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
4059 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
4065 /* Loop 4: Iterate the kernel with primitive */
4069 while ( kernel_loop < kernel_limit && changed > 0 ) {
4070 kernel_loop++; /* the iteration of this kernel */
4072 /* Create a clone as the destination image, if not yet defined */
4073 if ( work_image == (Image *) NULL )
4075 work_image=CloneImage(image,0,0,MagickTrue,exception);
4076 if (work_image == (Image *) NULL)
4078 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
4080 /* work_image->type=image->type; ??? */
4083 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4085 changed = MorphologyPrimitive(curr_image, work_image, primitive,
4086 this_kernel, bias, exception);
4088 if ( IfMagickTrue(verbose) ) {
4089 if ( kernel_loop > 1 )
4090 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4091 (void) (void) FormatLocaleFile(stderr,
4092 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4093 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4094 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4095 (double) (method_loop+kernel_loop-1),(double) kernel_number,
4096 (double) count,(double) changed);
4100 kernel_changed += changed;
4101 method_changed += changed;
4103 /* prepare next loop */
4104 { Image *tmp = work_image; /* swap images for iteration */
4105 work_image = curr_image;
4108 if ( work_image == image )
4109 work_image = (Image *) NULL; /* replace input 'image' */
4111 } /* End Loop 4: Iterate the kernel with primitive */
4113 if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed )
4114 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
4115 if ( IfMagickTrue(verbose) && stage_loop < stage_limit )
4116 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4119 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4120 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
4121 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
4122 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
4123 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
4126 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4128 /* Final Post-processing for some Compound Methods
4130 ** The removal of any 'Sync' channel flag in the Image Compositon
4131 ** below ensures the methematical compose method is applied in a
4132 ** purely mathematical way, and only to the selected channels.
4133 ** Turn off SVG composition 'alpha blending'.
4136 case EdgeOutMorphology:
4137 case EdgeInMorphology:
4138 case TopHatMorphology:
4139 case BottomHatMorphology:
4140 if ( IfMagickTrue(verbose) )
4141 (void) FormatLocaleFile(stderr,
4142 "\n%s: Difference with original image",CommandOptionToMnemonic(
4143 MagickMorphologyOptions, method) );
4144 (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
4145 MagickTrue,0,0,exception);
4147 case EdgeMorphology:
4148 if ( IfMagickTrue(verbose) )
4149 (void) FormatLocaleFile(stderr,
4150 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
4151 MagickMorphologyOptions, method) );
4152 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
4153 MagickTrue,0,0,exception);
4154 save_image = DestroyImage(save_image); /* finished with save image */
4160 /* multi-kernel handling: re-iterate, or compose results */
4161 if ( kernel->next == (KernelInfo *) NULL )
4162 rslt_image = curr_image; /* just return the resulting image */
4163 else if ( rslt_compose == NoCompositeOp )
4164 { if ( IfMagickTrue(verbose) ) {
4165 if ( this_kernel->next != (KernelInfo *) NULL )
4166 (void) FormatLocaleFile(stderr, " (re-iterate)");
4168 (void) FormatLocaleFile(stderr, " (done)");
4170 rslt_image = curr_image; /* return result, and re-iterate */
4172 else if ( rslt_image == (Image *) NULL)
4173 { if ( IfMagickTrue(verbose) )
4174 (void) FormatLocaleFile(stderr, " (save for compose)");
4175 rslt_image = curr_image;
4176 curr_image = (Image *) image; /* continue with original image */
4179 { /* Add the new 'current' result to the composition
4181 ** The removal of any 'Sync' channel flag in the Image Compositon
4182 ** below ensures the methematical compose method is applied in a
4183 ** purely mathematical way, and only to the selected channels.
4184 ** IE: Turn off SVG composition 'alpha blending'.
4186 if ( IfMagickTrue(verbose) )
4187 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4188 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4189 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4191 curr_image = DestroyImage(curr_image);
4192 curr_image = (Image *) image; /* continue with original image */
4194 if ( IfMagickTrue(verbose) )
4195 (void) FormatLocaleFile(stderr, "\n");
4197 /* loop to the next kernel in a multi-kernel list */
4198 norm_kernel = norm_kernel->next;
4199 if ( rflt_kernel != (KernelInfo *) NULL )
4200 rflt_kernel = rflt_kernel->next;
4202 } /* End Loop 2: Loop over each kernel */
4204 } /* End Loop 1: compound method interation */
4208 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4210 if ( curr_image == rslt_image )
4211 curr_image = (Image *) NULL;
4212 if ( rslt_image != (Image *) NULL )
4213 rslt_image = DestroyImage(rslt_image);
4215 if ( curr_image == rslt_image || curr_image == image )
4216 curr_image = (Image *) NULL;
4217 if ( curr_image != (Image *) NULL )
4218 curr_image = DestroyImage(curr_image);
4219 if ( work_image != (Image *) NULL )
4220 work_image = DestroyImage(work_image);
4221 if ( save_image != (Image *) NULL )
4222 save_image = DestroyImage(save_image);
4223 if ( reflected_kernel != (KernelInfo *) NULL )
4224 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4230 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4234 % M o r p h o l o g y I m a g e %
4238 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4240 % MorphologyImage() applies a user supplied kernel to the image
4241 % according to the given mophology method.
4243 % This function applies any and all user defined settings before calling
4244 % the above internal function MorphologyApply().
4246 % User defined settings include...
4247 % * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4248 % * Kernel Scale/normalize settings ("-define convolve:scale=??")
4249 % This can also includes the addition of a scaled unity kernel.
4250 % * Show Kernel being applied ("-define showkernel=1")
4252 % Other operators that do not want user supplied options interfering,
4253 % especially "convolve:bias" and "showkernel" should use MorphologyApply()
4256 % The format of the MorphologyImage method is:
4258 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4259 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4261 % A description of each parameter follows:
4263 % o image: the image.
4265 % o method: the morphology method to be applied.
4267 % o iterations: apply the operation this many times (or no change).
4268 % A value of -1 means loop until no change found.
4269 % How this is applied may depend on the morphology method.
4270 % Typically this is a value of 1.
4272 % o kernel: An array of double representing the morphology kernel.
4273 % Warning: kernel may be normalized for the Convolve method.
4275 % o exception: return any errors or warnings in this structure.
4278 MagickExport Image *MorphologyImage(const Image *image,
4279 const MorphologyMethod method,const ssize_t iterations,
4280 const KernelInfo *kernel,ExceptionInfo *exception)
4294 curr_kernel = (KernelInfo *) kernel;
4296 compose = UndefinedCompositeOp; /* use default for method */
4298 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4299 * This is done BEFORE the ShowKernelInfo() function is called so that
4300 * users can see the results of the 'option:convolve:scale' option.
4302 if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4306 /* Get the bias value as it will be needed */
4307 artifact = GetImageArtifact(image,"convolve:bias");
4308 if ( artifact != (const char *) NULL) {
4309 if (IfMagickFalse(IsGeometry(artifact)))
4310 (void) ThrowMagickException(exception,GetMagickModule(),
4311 OptionWarning,"InvalidSetting","'%s' '%s'",
4312 "convolve:bias",artifact);
4314 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4317 /* Scale kernel according to user wishes */
4318 artifact = GetImageArtifact(image,"convolve:scale");
4319 if ( artifact != (const char *)NULL ) {
4320 if (IfMagickFalse(IsGeometry(artifact)))
4321 (void) ThrowMagickException(exception,GetMagickModule(),
4322 OptionWarning,"InvalidSetting","'%s' '%s'",
4323 "convolve:scale",artifact);
4325 if ( curr_kernel == kernel )
4326 curr_kernel = CloneKernelInfo(kernel);
4327 if (curr_kernel == (KernelInfo *) NULL)
4328 return((Image *) NULL);
4329 ScaleGeometryKernelInfo(curr_kernel, artifact);
4334 /* display the (normalized) kernel via stderr */
4335 if ( IfStringTrue(GetImageArtifact(image,"showkernel"))
4336 || IfStringTrue(GetImageArtifact(image,"convolve:showkernel"))
4337 || IfStringTrue(GetImageArtifact(image,"morphology:showkernel")) )
4338 ShowKernelInfo(curr_kernel);
4340 /* Override the default handling of multi-kernel morphology results
4341 * If 'Undefined' use the default method
4342 * If 'None' (default for 'Convolve') re-iterate previous result
4343 * Otherwise merge resulting images using compose method given.
4344 * Default for 'HitAndMiss' is 'Lighten'.
4351 artifact = GetImageArtifact(image,"morphology:compose");
4352 if ( artifact != (const char *) NULL) {
4353 parse=ParseCommandOption(MagickComposeOptions,
4354 MagickFalse,artifact);
4356 (void) ThrowMagickException(exception,GetMagickModule(),
4357 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4358 "morphology:compose",artifact);
4360 compose=(CompositeOperator)parse;
4363 /* Apply the Morphology */
4364 morphology_image = MorphologyApply(image,method,iterations,
4365 curr_kernel,compose,bias,exception);
4367 /* Cleanup and Exit */
4368 if ( curr_kernel != kernel )
4369 curr_kernel=DestroyKernelInfo(curr_kernel);
4370 return(morphology_image);
4374 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4378 + R o t a t e K e r n e l I n f o %
4382 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4384 % RotateKernelInfo() rotates the kernel by the angle given.
4386 % Currently it is restricted to 90 degree angles, of either 1D kernels
4387 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4388 % It will ignore usless rotations for specific 'named' built-in kernels.
4390 % The format of the RotateKernelInfo method is:
4392 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4394 % A description of each parameter follows:
4396 % o kernel: the Morphology/Convolution kernel
4398 % o angle: angle to rotate in degrees
4400 % This function is currently internal to this module only, but can be exported
4401 % to other modules if needed.
4403 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4405 /* angle the lower kernels first */
4406 if ( kernel->next != (KernelInfo *) NULL)
4407 RotateKernelInfo(kernel->next, angle);
4409 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4411 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4414 /* Modulus the angle */
4415 angle = fmod(angle, 360.0);
4419 if ( 337.5 < angle || angle <= 22.5 )
4420 return; /* Near zero angle - no change! - At least not at this time */
4422 /* Handle special cases */
4423 switch (kernel->type) {
4424 /* These built-in kernels are cylindrical kernels, rotating is useless */
4425 case GaussianKernel:
4430 case LaplacianKernel:
4431 case ChebyshevKernel:
4432 case ManhattanKernel:
4433 case EuclideanKernel:
4436 /* These may be rotatable at non-90 angles in the future */
4437 /* but simply rotating them in multiples of 90 degrees is useless */
4444 /* These only allows a +/-90 degree rotation (by transpose) */
4445 /* A 180 degree rotation is useless */
4447 if ( 135.0 < angle && angle <= 225.0 )
4449 if ( 225.0 < angle && angle <= 315.0 )
4456 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4457 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4459 if ( kernel->width == 3 && kernel->height == 3 )
4460 { /* Rotate a 3x3 square by 45 degree angle */
4461 double t = kernel->values[0];
4462 kernel->values[0] = kernel->values[3];
4463 kernel->values[3] = kernel->values[6];
4464 kernel->values[6] = kernel->values[7];
4465 kernel->values[7] = kernel->values[8];
4466 kernel->values[8] = kernel->values[5];
4467 kernel->values[5] = kernel->values[2];
4468 kernel->values[2] = kernel->values[1];
4469 kernel->values[1] = t;
4470 /* rotate non-centered origin */
4471 if ( kernel->x != 1 || kernel->y != 1 ) {
4473 x = (ssize_t) kernel->x-1;
4474 y = (ssize_t) kernel->y-1;
4475 if ( x == y ) x = 0;
4476 else if ( x == 0 ) x = -y;
4477 else if ( x == -y ) y = 0;
4478 else if ( y == 0 ) y = x;
4479 kernel->x = (ssize_t) x+1;
4480 kernel->y = (ssize_t) y+1;
4482 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4483 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4486 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4488 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4490 if ( kernel->width == 1 || kernel->height == 1 )
4491 { /* Do a transpose of a 1 dimensional kernel,
4492 ** which results in a fast 90 degree rotation of some type.
4496 t = (ssize_t) kernel->width;
4497 kernel->width = kernel->height;
4498 kernel->height = (size_t) t;
4500 kernel->x = kernel->y;
4502 if ( kernel->width == 1 ) {
4503 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4504 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4506 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4507 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4510 else if ( kernel->width == kernel->height )
4511 { /* Rotate a square array of values by 90 degrees */
4515 register MagickRealType
4519 for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--)
4520 for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--)
4521 { t = k[i+j*kernel->width];
4522 k[i+j*kernel->width] = k[j+x*kernel->width];
4523 k[j+x*kernel->width] = k[x+y*kernel->width];
4524 k[x+y*kernel->width] = k[y+i*kernel->width];
4525 k[y+i*kernel->width] = t;
4528 /* rotate the origin - relative to center of array */
4529 { register ssize_t x,y;
4530 x = (ssize_t) (kernel->x*2-kernel->width+1);
4531 y = (ssize_t) (kernel->y*2-kernel->height+1);
4532 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4533 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4535 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4536 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4539 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4541 if ( 135.0 < angle && angle <= 225.0 )
4543 /* For a 180 degree rotation - also know as a reflection
4544 * This is actually a very very common operation!
4545 * Basically all that is needed is a reversal of the kernel data!
4546 * And a reflection of the origon
4551 register MagickRealType
4559 j=(ssize_t) (kernel->width*kernel->height-1);
4560 for (i=0; i < j; i++, j--)
4561 t=k[i], k[i]=k[j], k[j]=t;
4563 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4564 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4565 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4566 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4568 /* At this point angle should at least between -45 (315) and +45 degrees
4569 * In the future some form of non-orthogonal angled rotates could be
4570 * performed here, posibily with a linear kernel restriction.
4577 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4581 % S c a l e G e o m e t r y K e r n e l I n f o %
4585 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4587 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4588 % provided as a "-set option:convolve:scale {geometry}" user setting,
4589 % and modifies the kernel according to the parsed arguments of that setting.
4591 % The first argument (and any normalization flags) are passed to
4592 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4593 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4594 % into the scaled/normalized kernel.
4596 % The format of the ScaleGeometryKernelInfo method is:
4598 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4599 % const double scaling_factor,const MagickStatusType normalize_flags)
4601 % A description of each parameter follows:
4603 % o kernel: the Morphology/Convolution kernel to modify
4606 % The geometry string to parse, typically from the user provided
4607 % "-set option:convolve:scale {geometry}" setting.
4610 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4611 const char *geometry)
4619 SetGeometryInfo(&args);
4620 flags = ParseGeometry(geometry, &args);
4623 /* For Debugging Geometry Input */
4624 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4625 flags, args.rho, args.sigma, args.xi, args.psi );
4628 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4629 args.rho *= 0.01, args.sigma *= 0.01;
4631 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4633 if ( (flags & SigmaValue) == 0 )
4636 /* Scale/Normalize the input kernel */
4637 ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4639 /* Add Unity Kernel, for blending with original */
4640 if ( (flags & SigmaValue) != 0 )
4641 UnityAddKernelInfo(kernel, args.sigma);
4646 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4650 % S c a l e K e r n e l I n f o %
4654 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4656 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4657 % without normalization of the sum of the kernel values (as per given flags).
4659 % By default (no flags given) the values within the kernel is scaled
4660 % directly using given scaling factor without change.
4662 % If either of the two 'normalize_flags' are given the kernel will first be
4663 % normalized and then further scaled by the scaling factor value given.
4665 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4666 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4667 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4668 % non-HDRI versions of IM this may cause images to have any negative results
4669 % clipped, unless some 'bias' is used.
4671 % More specifically. Kernels which only contain positive values (such as a
4672 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4673 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4675 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4676 % the kernel will be scaled by the absolute of the sum of kernel values, so
4677 % that it will generally fall within the +/- 1.0 range.
4679 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4680 % will be scaled by just the sum of the postive values, so that its output
4681 % range will again fall into the +/- 1.0 range.
4683 % For special kernels designed for locating shapes using 'Correlate', (often
4684 % only containing +1 and -1 values, representing foreground/brackground
4685 % matching) a special normalization method is provided to scale the positive
4686 % values separately to those of the negative values, so the kernel will be
4687 % forced to become a zero-sum kernel better suited to such searches.
4689 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4690 % attributes within the kernel structure have been correctly set during the
4693 % NOTE: The values used for 'normalize_flags' have been selected specifically
4694 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4695 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4697 % The format of the ScaleKernelInfo method is:
4699 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4700 % const MagickStatusType normalize_flags )
4702 % A description of each parameter follows:
4704 % o kernel: the Morphology/Convolution kernel
4707 % multiply all values (after normalization) by this factor if not
4708 % zero. If the kernel is normalized regardless of any flags.
4710 % o normalize_flags:
4711 % GeometryFlags defining normalization method to use.
4712 % specifically: NormalizeValue, CorrelateNormalizeValue,
4713 % and/or PercentValue
4716 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4717 const double scaling_factor,const GeometryFlags normalize_flags)
4726 /* do the other kernels in a multi-kernel list first */
4727 if ( kernel->next != (KernelInfo *) NULL)
4728 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4730 /* Normalization of Kernel */
4732 if ( (normalize_flags&NormalizeValue) != 0 ) {
4733 if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4734 /* non-zero-summing kernel (generally positive) */
4735 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4737 /* zero-summing kernel */
4738 pos_scale = kernel->positive_range;
4740 /* Force kernel into a normalized zero-summing kernel */
4741 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4742 pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4743 ? kernel->positive_range : 1.0;
4744 neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4745 ? -kernel->negative_range : 1.0;
4748 neg_scale = pos_scale;
4750 /* finialize scaling_factor for positive and negative components */
4751 pos_scale = scaling_factor/pos_scale;
4752 neg_scale = scaling_factor/neg_scale;
4754 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4755 if ( ! IsNaN(kernel->values[i]) )
4756 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4758 /* convolution output range */
4759 kernel->positive_range *= pos_scale;
4760 kernel->negative_range *= neg_scale;
4761 /* maximum and minimum values in kernel */
4762 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4763 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4765 /* swap kernel settings if user's scaling factor is negative */
4766 if ( scaling_factor < MagickEpsilon ) {
4768 t = kernel->positive_range;
4769 kernel->positive_range = kernel->negative_range;
4770 kernel->negative_range = t;
4771 t = kernel->maximum;
4772 kernel->maximum = kernel->minimum;
4773 kernel->minimum = 1;
4780 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4784 % S h o w K e r n e l I n f o %
4788 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4790 % ShowKernelInfo() outputs the details of the given kernel defination to
4791 % standard error, generally due to a users 'showkernel' option request.
4793 % The format of the ShowKernel method is:
4795 % void ShowKernelInfo(const KernelInfo *kernel)
4797 % A description of each parameter follows:
4799 % o kernel: the Morphology/Convolution kernel
4802 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4810 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4812 (void) FormatLocaleFile(stderr, "Kernel");
4813 if ( kernel->next != (KernelInfo *) NULL )
4814 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4815 (void) FormatLocaleFile(stderr, " \"%s",
4816 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4817 if ( fabs(k->angle) >= MagickEpsilon )
4818 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4819 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4820 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4821 (void) FormatLocaleFile(stderr,
4822 " with values from %.*lg to %.*lg\n",
4823 GetMagickPrecision(), k->minimum,
4824 GetMagickPrecision(), k->maximum);
4825 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4826 GetMagickPrecision(), k->negative_range,
4827 GetMagickPrecision(), k->positive_range);
4828 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4829 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4830 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4831 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4833 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4834 GetMagickPrecision(), k->positive_range+k->negative_range);
4835 for (i=v=0; v < k->height; v++) {
4836 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4837 for (u=0; u < k->width; u++, i++)
4838 if ( IsNaN(k->values[i]) )
4839 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4841 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4842 GetMagickPrecision(), (double) k->values[i]);
4843 (void) FormatLocaleFile(stderr,"\n");
4849 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4853 % U n i t y A d d K e r n a l I n f o %
4857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4859 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4860 % to the given pre-scaled and normalized Kernel. This in effect adds that
4861 % amount of the original image into the resulting convolution kernel. This
4862 % value is usually provided by the user as a percentage value in the
4863 % 'convolve:scale' setting.
4865 % The resulting effect is to convert the defined kernels into blended
4866 % soft-blurs, unsharp kernels or into sharpening kernels.
4868 % The format of the UnityAdditionKernelInfo method is:
4870 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4872 % A description of each parameter follows:
4874 % o kernel: the Morphology/Convolution kernel
4877 % scaling factor for the unity kernel to be added to
4881 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4884 /* do the other kernels in a multi-kernel list first */
4885 if ( kernel->next != (KernelInfo *) NULL)
4886 UnityAddKernelInfo(kernel->next, scale);
4888 /* Add the scaled unity kernel to the existing kernel */
4889 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4890 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4896 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4900 % Z e r o K e r n e l N a n s %
4904 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4906 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4907 % the kernel with a zero value. This is typically done when the kernel will
4908 % be used in special hardware (GPU) convolution processors, to simply
4911 % The format of the ZeroKernelNans method is:
4913 % void ZeroKernelNans (KernelInfo *kernel)
4915 % A description of each parameter follows:
4917 % o kernel: the Morphology/Convolution kernel
4920 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4925 /* do the other kernels in a multi-kernel list first */
4926 if ( kernel->next != (KernelInfo *) NULL)
4927 ZeroKernelNans(kernel->next);
4929 for (i=0; i < (kernel->width*kernel->height); i++)
4930 if ( IsNaN(kernel->values[i]) )
4931 kernel->values[i] = 0.0;