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-2014 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/channel.h"
56 #include "MagickCore/color-private.h"
57 #include "MagickCore/enhance.h"
58 #include "MagickCore/exception.h"
59 #include "MagickCore/exception-private.h"
60 #include "MagickCore/gem.h"
61 #include "MagickCore/gem-private.h"
62 #include "MagickCore/hashmap.h"
63 #include "MagickCore/image.h"
64 #include "MagickCore/image-private.h"
65 #include "MagickCore/list.h"
66 #include "MagickCore/magick.h"
67 #include "MagickCore/memory_.h"
68 #include "MagickCore/memory-private.h"
69 #include "MagickCore/monitor-private.h"
70 #include "MagickCore/morphology.h"
71 #include "MagickCore/morphology-private.h"
72 #include "MagickCore/option.h"
73 #include "MagickCore/pixel-accessor.h"
74 #include "MagickCore/pixel-private.h"
75 #include "MagickCore/prepress.h"
76 #include "MagickCore/quantize.h"
77 #include "MagickCore/resource_.h"
78 #include "MagickCore/registry.h"
79 #include "MagickCore/semaphore.h"
80 #include "MagickCore/splay-tree.h"
81 #include "MagickCore/statistic.h"
82 #include "MagickCore/string_.h"
83 #include "MagickCore/string-private.h"
84 #include "MagickCore/thread-private.h"
85 #include "MagickCore/token.h"
86 #include "MagickCore/utility.h"
87 #include "MagickCore/utility-private.h"
90 Other global definitions used by module.
92 static inline double MagickMin(const double x,const double y)
94 return( x < y ? x : y);
96 static inline double MagickMax(const double x,const double y)
98 return( x > y ? x : y);
100 #define Minimize(assign,value) assign=MagickMin(assign,value)
101 #define Maximize(assign,value) assign=MagickMax(assign,value)
103 /* Integer Factorial Function - for a Binomial kernel */
105 static inline size_t fact(size_t n)
108 for(f=1, l=2; l <= n; f=f*l, l++);
111 #elif 1 /* glibc floating point alternatives */
112 #define fact(n) ((size_t)tgamma((double)n+1))
114 #define fact(n) ((size_t)lgamma((double)n+1))
118 /* Currently these are only internal to this module */
120 CalcKernelMetaData(KernelInfo *),
121 ExpandMirrorKernelInfo(KernelInfo *),
122 ExpandRotateKernelInfo(KernelInfo *, const double),
123 RotateKernelInfo(KernelInfo *, double);
126 /* Quick function to find last kernel in a kernel list */
127 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
129 while (kernel->next != (KernelInfo *) NULL)
130 kernel = kernel->next;
135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139 % A c q u i r e K e r n e l I n f o %
143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145 % AcquireKernelInfo() takes the given string (generally supplied by the
146 % user) and converts it into a Morphology/Convolution Kernel. This allows
147 % users to specify a kernel from a number of pre-defined kernels, or to fully
148 % specify their own kernel for a specific Convolution or Morphology
151 % The kernel so generated can be any rectangular array of floating point
152 % values (doubles) with the 'control point' or 'pixel being affected'
153 % anywhere within that array of values.
155 % Previously IM was restricted to a square of odd size using the exact
156 % center as origin, this is no longer the case, and any rectangular kernel
157 % with any value being declared the origin. This in turn allows the use of
158 % highly asymmetrical kernels.
160 % The floating point values in the kernel can also include a special value
161 % known as 'nan' or 'not a number' to indicate that this value is not part
162 % of the kernel array. This allows you to shaped the kernel within its
163 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
164 % shape. However at least one non-nan value must be provided for correct
165 % working of a kernel.
167 % The returned kernel should be freed using the DestroyKernelInfo() when you
168 % are finished with it. Do not free this memory yourself.
170 % Input kernel defintion strings can consist of any of three types.
173 % Select from one of the built in kernels, using the name and
174 % geometry arguments supplied. See AcquireKernelBuiltIn()
176 % "WxH[+X+Y][@><]:num, num, num ..."
177 % a kernel of size W by H, with W*H floating point numbers following.
178 % the 'center' can be optionally be defined at +X+Y (such that +0+0
179 % is top left corner). If not defined the pixel in the center, for
180 % odd sizes, or to the immediate top or left of center for even sizes
181 % is automatically selected.
183 % "num, num, num, num, ..."
184 % list of floating point numbers defining an 'old style' odd sized
185 % square kernel. At least 9 values should be provided for a 3x3
186 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
187 % Values can be space or comma separated. This is not recommended.
189 % You can define a 'list of kernels' which can be used by some morphology
190 % operators A list is defined as a semi-colon separated list kernels.
192 % " kernel ; kernel ; kernel ; "
194 % Any extra ';' characters, at start, end or between kernel defintions are
197 % The special flags will expand a single kernel, into a list of rotated
198 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
199 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
200 % The '<' also exands using 90-degree rotates, but giving a 180-degree
201 % reflected kernel before the +/- 90-degree rotations, which can be important
202 % for Thinning operations.
204 % Note that 'name' kernels will start with an alphabetic character while the
205 % new kernel specification has a ':' character in its specification string.
206 % If neither is the case, it is assumed an old style of a simple list of
207 % numbers generating a odd-sized square kernel has been given.
209 % The format of the AcquireKernal method is:
211 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
213 % A description of each parameter follows:
215 % o kernel_string: the Morphology/Convolution kernel wanted.
219 /* This was separated so that it could be used as a separate
220 ** array input handling function, such as for -color-matrix
222 static KernelInfo *ParseKernelArray(const char *kernel_string)
228 token[MaxTextExtent];
238 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
246 kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
247 if (kernel == (KernelInfo *)NULL)
249 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
250 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
251 kernel->negative_range = kernel->positive_range = 0.0;
252 kernel->type = UserDefinedKernel;
253 kernel->next = (KernelInfo *) NULL;
254 kernel->signature = MagickSignature;
255 if (kernel_string == (const char *) NULL)
258 /* find end of this specific kernel definition string */
259 end = strchr(kernel_string, ';');
260 if ( end == (char *) NULL )
261 end = strchr(kernel_string, '\0');
263 /* clear flags - for Expanding kernel lists thorugh rotations */
266 /* Has a ':' in argument - New user kernel specification
267 FUTURE: this split on ':' could be done by StringToken()
269 p = strchr(kernel_string, ':');
270 if ( p != (char *) NULL && p < end)
272 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
273 memcpy(token, kernel_string, (size_t) (p-kernel_string));
274 token[p-kernel_string] = '\0';
275 SetGeometryInfo(&args);
276 flags = ParseGeometry(token, &args);
278 /* Size handling and checks of geometry settings */
279 if ( (flags & WidthValue) == 0 ) /* if no width then */
280 args.rho = args.sigma; /* then width = height */
281 if ( args.rho < 1.0 ) /* if width too small */
282 args.rho = 1.0; /* then width = 1 */
283 if ( args.sigma < 1.0 ) /* if height too small */
284 args.sigma = args.rho; /* then height = width */
285 kernel->width = (size_t)args.rho;
286 kernel->height = (size_t)args.sigma;
288 /* Offset Handling and Checks */
289 if ( args.xi < 0.0 || args.psi < 0.0 )
290 return(DestroyKernelInfo(kernel));
291 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
292 : (ssize_t) (kernel->width-1)/2;
293 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
294 : (ssize_t) (kernel->height-1)/2;
295 if ( kernel->x >= (ssize_t) kernel->width ||
296 kernel->y >= (ssize_t) kernel->height )
297 return(DestroyKernelInfo(kernel));
299 p++; /* advance beyond the ':' */
302 { /* ELSE - Old old specification, forming odd-square kernel */
303 /* count up number of values given */
304 p=(const char *) kernel_string;
305 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
306 p++; /* ignore "'" chars for convolve filter usage - Cristy */
307 for (i=0; p < end; i++)
309 GetMagickToken(p,&p,token);
311 GetMagickToken(p,&p,token);
313 /* set the size of the kernel - old sized square */
314 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
315 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
316 p=(const char *) kernel_string;
317 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
318 p++; /* ignore "'" chars for convolve filter usage - Cristy */
321 /* Read in the kernel values from rest of input string argument */
322 kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
323 kernel->width,kernel->height*sizeof(*kernel->values)));
324 if (kernel->values == (MagickRealType *) NULL)
325 return(DestroyKernelInfo(kernel));
326 kernel->minimum = +MagickMaximumValue;
327 kernel->maximum = -MagickMaximumValue;
328 kernel->negative_range = kernel->positive_range = 0.0;
329 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
331 GetMagickToken(p,&p,token);
333 GetMagickToken(p,&p,token);
334 if ( LocaleCompare("nan",token) == 0
335 || LocaleCompare("-",token) == 0 ) {
336 kernel->values[i] = nan; /* this value is not part of neighbourhood */
339 kernel->values[i] = StringToDouble(token,(char **) NULL);
340 ( kernel->values[i] < 0)
341 ? ( kernel->negative_range += kernel->values[i] )
342 : ( kernel->positive_range += kernel->values[i] );
343 Minimize(kernel->minimum, kernel->values[i]);
344 Maximize(kernel->maximum, kernel->values[i]);
348 /* sanity check -- no more values in kernel definition */
349 GetMagickToken(p,&p,token);
350 if ( *token != '\0' && *token != ';' && *token != '\'' )
351 return(DestroyKernelInfo(kernel));
354 /* this was the old method of handling a incomplete kernel */
355 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
356 Minimize(kernel->minimum, kernel->values[i]);
357 Maximize(kernel->maximum, kernel->values[i]);
358 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
359 kernel->values[i]=0.0;
362 /* Number of values for kernel was not enough - Report Error */
363 if ( i < (ssize_t) (kernel->width*kernel->height) )
364 return(DestroyKernelInfo(kernel));
367 /* check that we recieved at least one real (non-nan) value! */
368 if ( kernel->minimum == MagickMaximumValue )
369 return(DestroyKernelInfo(kernel));
371 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
372 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
373 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
374 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
375 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
376 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
381 static KernelInfo *ParseKernelName(const char *kernel_string)
384 token[MaxTextExtent];
402 /* Parse special 'named' kernel */
403 GetMagickToken(kernel_string,&p,token);
404 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
405 if ( type < 0 || type == UserDefinedKernel )
406 return((KernelInfo *)NULL); /* not a valid named kernel */
408 while (((isspace((int) ((unsigned char) *p)) != 0) ||
409 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
412 end = strchr(p, ';'); /* end of this kernel defintion */
413 if ( end == (char *) NULL )
414 end = strchr(p, '\0');
416 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
417 memcpy(token, p, (size_t) (end-p));
419 SetGeometryInfo(&args);
420 flags = ParseGeometry(token, &args);
423 /* For Debugging Geometry Input */
424 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
425 flags, args.rho, args.sigma, args.xi, args.psi );
428 /* special handling of missing values in input string */
430 /* Shape Kernel Defaults */
432 if ( (flags & WidthValue) == 0 )
433 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
441 if ( (flags & HeightValue) == 0 )
442 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
445 if ( (flags & XValue) == 0 )
446 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
448 case RectangleKernel: /* Rectangle - set size defaults */
449 if ( (flags & WidthValue) == 0 ) /* if no width then */
450 args.rho = args.sigma; /* then width = height */
451 if ( args.rho < 1.0 ) /* if width too small */
452 args.rho = 3; /* then width = 3 */
453 if ( args.sigma < 1.0 ) /* if height too small */
454 args.sigma = args.rho; /* then height = width */
455 if ( (flags & XValue) == 0 ) /* center offset if not defined */
456 args.xi = (double)(((ssize_t)args.rho-1)/2);
457 if ( (flags & YValue) == 0 )
458 args.psi = (double)(((ssize_t)args.sigma-1)/2);
460 /* Distance Kernel Defaults */
461 case ChebyshevKernel:
462 case ManhattanKernel:
463 case OctagonalKernel:
464 case EuclideanKernel:
465 if ( (flags & HeightValue) == 0 ) /* no distance scale */
466 args.sigma = 100.0; /* default distance scaling */
467 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
468 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
469 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
470 args.sigma *= QuantumRange/100.0; /* percentage of color range */
476 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
477 if ( kernel == (KernelInfo *) NULL )
480 /* global expand to rotated kernel list - only for single kernels */
481 if ( kernel->next == (KernelInfo *) NULL ) {
482 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
483 ExpandRotateKernelInfo(kernel, 45.0);
484 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
485 ExpandRotateKernelInfo(kernel, 90.0);
486 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
487 ExpandMirrorKernelInfo(kernel);
493 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
501 token[MaxTextExtent];
509 if (kernel_string == (const char *) NULL)
510 return(ParseKernelArray(kernel_string));
515 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
517 /* ignore extra or multiple ';' kernel separators */
518 if ( *token != ';' ) {
520 /* tokens starting with alpha is a Named kernel */
521 if (isalpha((int) ((unsigned char) *token)) != 0)
522 new_kernel = ParseKernelName(p);
523 else /* otherwise a user defined kernel array */
524 new_kernel = ParseKernelArray(p);
526 /* Error handling -- this is not proper error handling! */
527 if ( new_kernel == (KernelInfo *) NULL ) {
528 (void) FormatLocaleFile(stderr,"Failed to parse kernel number #%.20g\n",
529 (double) kernel_number);
530 if ( kernel != (KernelInfo *) NULL )
531 kernel=DestroyKernelInfo(kernel);
532 return((KernelInfo *) NULL);
535 /* initialise or append the kernel list */
536 if ( kernel == (KernelInfo *) NULL )
539 LastKernelInfo(kernel)->next = new_kernel;
542 /* look for the next kernel in list */
544 if ( p == (char *) NULL )
554 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
558 % A c q u i r e K e r n e l B u i l t I n %
562 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
564 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
565 % kernels used for special purposes such as gaussian blurring, skeleton
566 % pruning, and edge distance determination.
568 % They take a KernelType, and a set of geometry style arguments, which were
569 % typically decoded from a user supplied string, or from a more complex
570 % Morphology Method that was requested.
572 % The format of the AcquireKernalBuiltIn method is:
574 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
575 % const GeometryInfo args)
577 % A description of each parameter follows:
579 % o type: the pre-defined type of kernel wanted
581 % o args: arguments defining or modifying the kernel
583 % Convolution Kernels
586 % The a No-Op or Scaling single element kernel.
588 % Gaussian:{radius},{sigma}
589 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
590 % The sigma for the curve is required. The resulting kernel is
593 % If 'sigma' is zero, you get a single pixel on a field of zeros.
595 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
596 % the final size of the resulting kernel to a square 2*radius+1 in size.
597 % The radius should be at least 2 times that of the sigma value, or
598 % sever clipping and aliasing may result. If not given or set to 0 the
599 % radius will be determined so as to produce the best minimal error
600 % result, which is usally much larger than is normally needed.
602 % LoG:{radius},{sigma}
603 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
604 % The supposed ideal edge detection, zero-summing kernel.
606 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
607 % approx 1.6 (according to wikipedia).
609 % DoG:{radius},{sigma1},{sigma2}
610 % "Difference of Gaussians" Kernel.
611 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
612 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
613 % The result is a zero-summing kernel.
615 % Blur:{radius},{sigma}[,{angle}]
616 % Generates a 1 dimensional or linear gaussian blur, at the angle given
617 % (current restricted to orthogonal angles). If a 'radius' is given the
618 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
619 % by a 90 degree angle.
621 % If 'sigma' is zero, you get a single pixel on a field of zeros.
623 % Note that two convolutions with two "Blur" kernels perpendicular to
624 % each other, is equivalent to a far larger "Gaussian" kernel with the
625 % same sigma value, However it is much faster to apply. This is how the
626 % "-blur" operator actually works.
628 % Comet:{width},{sigma},{angle}
629 % Blur in one direction only, much like how a bright object leaves
630 % a comet like trail. The Kernel is actually half a gaussian curve,
631 % Adding two such blurs in opposite directions produces a Blur Kernel.
632 % Angle can be rotated in multiples of 90 degrees.
634 % Note that the first argument is the width of the kernel and not the
635 % radius of the kernel.
637 % Binomial:[{radius}]
638 % Generate a discrete kernel using a 2 dimentional Pascel's Triangle
639 % of values. Used for special forma of image filters.
641 % # Still to be implemented...
645 % # Set kernel values using a resize filter, and given scale (sigma)
646 % # Cylindrical or Linear. Is this possible with an image?
649 % Named Constant Convolution Kernels
651 % All these are unscaled, zero-summing kernels by default. As such for
652 % non-HDRI version of ImageMagick some form of normalization, user scaling,
653 % and biasing the results is recommended, to prevent the resulting image
656 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
657 % 45 degrees to generate the 8 angled varients of each of the kernels.
660 % Discrete Lapacian Kernels, (without normalization)
661 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
662 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
663 % Type 2 : 3x3 with center:4 edge:1 corner:-2
664 % Type 3 : 3x3 with center:4 edge:-2 corner:1
665 % Type 5 : 5x5 laplacian
666 % Type 7 : 7x7 laplacian
667 % Type 15 : 5x5 LoG (sigma approx 1.4)
668 % Type 19 : 9x9 LoG (sigma approx 1.4)
671 % Sobel 'Edge' convolution kernel (3x3)
677 % Roberts convolution kernel (3x3)
683 % Prewitt Edge convolution kernel (3x3)
689 % Prewitt's "Compass" convolution kernel (3x3)
695 % Kirsch's "Compass" convolution kernel (3x3)
701 % Frei-Chen Edge Detector is based on a kernel that is similar to
702 % the Sobel Kernel, but is designed to be isotropic. That is it takes
703 % into account the distance of the diagonal in the kernel.
706 % | sqrt(2), 0, -sqrt(2) |
709 % FreiChen:{type},{angle}
711 % Frei-Chen Pre-weighted kernels...
713 % Type 0: default un-nomalized version shown above.
715 % Type 1: Orthogonal Kernel (same as type 11 below)
717 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
720 % Type 2: Diagonal form of Kernel...
722 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
725 % However this kernel is als at the heart of the FreiChen Edge Detection
726 % Process which uses a set of 9 specially weighted kernel. These 9
727 % kernels not be normalized, but directly applied to the image. The
728 % results is then added together, to produce the intensity of an edge in
729 % a specific direction. The square root of the pixel value can then be
730 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
731 % from each other, both the direction and the strength of the edge can be
734 % Type 10: All 9 of the following pre-weighted kernels...
736 % Type 11: | 1, 0, -1 |
737 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
740 % Type 12: | 1, sqrt(2), 1 |
741 % | 0, 0, 0 | / 2*sqrt(2)
744 % Type 13: | sqrt(2), -1, 0 |
745 % | -1, 0, 1 | / 2*sqrt(2)
748 % Type 14: | 0, 1, -sqrt(2) |
749 % | -1, 0, 1 | / 2*sqrt(2)
752 % Type 15: | 0, -1, 0 |
756 % Type 16: | 1, 0, -1 |
760 % Type 17: | 1, -2, 1 |
764 % Type 18: | -2, 1, -2 |
768 % Type 19: | 1, 1, 1 |
772 % The first 4 are for edge detection, the next 4 are for line detection
773 % and the last is to add a average component to the results.
775 % Using a special type of '-1' will return all 9 pre-weighted kernels
776 % as a multi-kernel list, so that you can use them directly (without
777 % normalization) with the special "-set option:morphology:compose Plus"
778 % setting to apply the full FreiChen Edge Detection Technique.
780 % If 'type' is large it will be taken to be an actual rotation angle for
781 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
782 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
784 % WARNING: The above was layed out as per
785 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
786 % But rotated 90 degrees so direction is from left rather than the top.
787 % I have yet to find any secondary confirmation of the above. The only
788 % other source found was actual source code at
789 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
790 % Neigher paper defineds the kernels in a way that looks locical or
791 % correct when taken as a whole.
795 % Diamond:[{radius}[,{scale}]]
796 % Generate a diamond shaped kernel with given radius to the points.
797 % Kernel size will again be radius*2+1 square and defaults to radius 1,
798 % generating a 3x3 kernel that is slightly larger than a square.
800 % Square:[{radius}[,{scale}]]
801 % Generate a square shaped kernel of size radius*2+1, and defaulting
802 % to a 3x3 (radius 1).
804 % Octagon:[{radius}[,{scale}]]
805 % Generate octagonal shaped kernel of given radius and constant scale.
806 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
807 % in "Diamond" kernel.
809 % Disk:[{radius}[,{scale}]]
810 % Generate a binary disk, thresholded at the radius given, the radius
811 % may be a float-point value. Final Kernel size is floor(radius)*2+1
812 % square. A radius of 5.3 is the default.
814 % NOTE: That a low radii Disk kernels produce the same results as
815 % many of the previously defined kernels, but differ greatly at larger
816 % radii. Here is a table of equivalences...
817 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
818 % "Disk:1.5" => "Square"
819 % "Disk:2" => "Diamond:2"
820 % "Disk:2.5" => "Octagon"
821 % "Disk:2.9" => "Square:2"
822 % "Disk:3.5" => "Octagon:3"
823 % "Disk:4.5" => "Octagon:4"
824 % "Disk:5.4" => "Octagon:5"
825 % "Disk:6.4" => "Octagon:6"
826 % All other Disk shapes are unique to this kernel, but because a "Disk"
827 % is more circular when using a larger radius, using a larger radius is
828 % preferred over iterating the morphological operation.
830 % Rectangle:{geometry}
831 % Simply generate a rectangle of 1's with the size given. You can also
832 % specify the location of the 'control point', otherwise the closest
833 % pixel to the center of the rectangle is selected.
835 % Properly centered and odd sized rectangles work the best.
837 % Symbol Dilation Kernels
839 % These kernel is not a good general morphological kernel, but is used
840 % more for highlighting and marking any single pixels in an image using,
841 % a "Dilate" method as appropriate.
843 % For the same reasons iterating these kernels does not produce the
844 % same result as using a larger radius for the symbol.
846 % Plus:[{radius}[,{scale}]]
847 % Cross:[{radius}[,{scale}]]
848 % Generate a kernel in the shape of a 'plus' or a 'cross' with
849 % a each arm the length of the given radius (default 2).
851 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
853 % Ring:{radius1},{radius2}[,{scale}]
854 % A ring of the values given that falls between the two radii.
855 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
856 % This is the 'edge' pixels of the default "Disk" kernel,
857 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
859 % Hit and Miss Kernels
861 % Peak:radius1,radius2
862 % Find any peak larger than the pixels the fall between the two radii.
863 % The default ring of pixels is as per "Ring".
865 % Find flat orthogonal edges of a binary shape
867 % Find 90 degree corners of a binary shape
869 % A special kernel to thin the 'outside' of diagonals
871 % Find end points of lines (for pruning a skeletion)
872 % Two types of lines ends (default to both) can be searched for
873 % Type 0: All line ends
874 % Type 1: single kernel for 4-conneected line ends
875 % Type 2: single kernel for simple line ends
877 % Find three line junctions (within a skeletion)
878 % Type 0: all line junctions
879 % Type 1: Y Junction kernel
880 % Type 2: Diagonal T Junction kernel
881 % Type 3: Orthogonal T Junction kernel
882 % Type 4: Diagonal X Junction kernel
883 % Type 5: Orthogonal + Junction kernel
885 % Find single pixel ridges or thin lines
886 % Type 1: Fine single pixel thick lines and ridges
887 % Type 2: Find two pixel thick lines and ridges
889 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
891 % Traditional skeleton generating kernels.
892 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
893 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
894 % Type 3: Thinning skeleton based on a ressearch paper by
895 % Dan S. Bloomberg (Default Type)
897 % A huge variety of Thinning Kernels designed to preserve conectivity.
898 % many other kernel sets use these kernels as source definitions.
899 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
900 % the super and sub notations used in the source research paper.
902 % Distance Measuring Kernels
904 % Different types of distance measuring methods, which are used with the
905 % a 'Distance' morphology method for generating a gradient based on
906 % distance from an edge of a binary shape, though there is a technique
907 % for handling a anti-aliased shape.
909 % See the 'Distance' Morphological Method, for information of how it is
912 % Chebyshev:[{radius}][x{scale}[%!]]
913 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
914 % is a value of one to any neighbour, orthogonal or diagonal. One why
915 % of thinking of it is the number of squares a 'King' or 'Queen' in
916 % chess needs to traverse reach any other position on a chess board.
917 % It results in a 'square' like distance function, but one where
918 % diagonals are given a value that is closer than expected.
920 % Manhattan:[{radius}][x{scale}[%!]]
921 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
922 % Cab distance metric), it is the distance needed when you can only
923 % travel in horizontal or vertical directions only. It is the
924 % distance a 'Rook' in chess would have to travel, and results in a
925 % diamond like distances, where diagonals are further than expected.
927 % Octagonal:[{radius}][x{scale}[%!]]
928 % An interleving of Manhatten and Chebyshev metrics producing an
929 % increasing octagonally shaped distance. Distances matches those of
930 % the "Octagon" shaped kernel of the same radius. The minimum radius
931 % and default is 2, producing a 5x5 kernel.
933 % Euclidean:[{radius}][x{scale}[%!]]
934 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
935 % However by default the kernel size only has a radius of 1, which
936 % limits the distance to 'Knight' like moves, with only orthogonal and
937 % diagonal measurements being correct. As such for the default kernel
938 % you will get octagonal like distance function.
940 % However using a larger radius such as "Euclidean:4" you will get a
941 % much smoother distance gradient from the edge of the shape. Especially
942 % if the image is pre-processed to include any anti-aliasing pixels.
943 % Of course a larger kernel is slower to use, and not always needed.
945 % The first three Distance Measuring Kernels will only generate distances
946 % of exact multiples of {scale} in binary images. As such you can use a
947 % scale of 1 without loosing any information. However you also need some
948 % scaling when handling non-binary anti-aliased shapes.
950 % The "Euclidean" Distance Kernel however does generate a non-integer
951 % fractional results, and as such scaling is vital even for binary shapes.
955 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
956 const GeometryInfo *args)
969 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
971 /* Generate a new empty kernel if needed */
972 kernel=(KernelInfo *) NULL;
974 case UndefinedKernel: /* These should not call this function */
975 case UserDefinedKernel:
976 assert("Should not call this function" != (char *)NULL);
978 case LaplacianKernel: /* Named Descrete Convolution Kernels */
979 case SobelKernel: /* these are defined using other kernels */
985 case EdgesKernel: /* Hit and Miss kernels */
987 case DiagonalsKernel:
989 case LineJunctionsKernel:
991 case ConvexHullKernel:
994 break; /* A pre-generated kernel is not needed */
996 /* set to 1 to do a compile-time check that we haven't missed anything */
1003 case BinomialKernel:
1006 case RectangleKernel:
1013 case ChebyshevKernel:
1014 case ManhattanKernel:
1015 case OctangonalKernel:
1016 case EuclideanKernel:
1020 /* Generate the base Kernel Structure */
1021 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1022 if (kernel == (KernelInfo *) NULL)
1024 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1025 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1026 kernel->negative_range = kernel->positive_range = 0.0;
1027 kernel->type = type;
1028 kernel->next = (KernelInfo *) NULL;
1029 kernel->signature = MagickSignature;
1039 kernel->height = kernel->width = (size_t) 1;
1040 kernel->x = kernel->y = (ssize_t) 0;
1041 kernel->values=(MagickRealType *) MagickAssumeAligned(
1042 AcquireAlignedMemory(1,sizeof(*kernel->values)));
1043 if (kernel->values == (MagickRealType *) NULL)
1044 return(DestroyKernelInfo(kernel));
1045 kernel->maximum = kernel->values[0] = args->rho;
1049 case GaussianKernel:
1053 sigma = fabs(args->sigma),
1054 sigma2 = fabs(args->xi),
1057 if ( args->rho >= 1.0 )
1058 kernel->width = (size_t)args->rho*2+1;
1059 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1060 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1062 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1063 kernel->height = kernel->width;
1064 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1065 kernel->values=(MagickRealType *) MagickAssumeAligned(
1066 AcquireAlignedMemory(kernel->width,kernel->height*
1067 sizeof(*kernel->values)));
1068 if (kernel->values == (MagickRealType *) NULL)
1069 return(DestroyKernelInfo(kernel));
1071 /* WARNING: The following generates a 'sampled gaussian' kernel.
1072 * What we really want is a 'discrete gaussian' kernel.
1074 * How to do this is I don't know, but appears to be basied on the
1075 * Error Function 'erf()' (intergral of a gaussian)
1078 if ( type == GaussianKernel || type == DoGKernel )
1079 { /* Calculate a Gaussian, OR positive half of a DoG */
1080 if ( sigma > MagickEpsilon )
1081 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1082 B = (double) (1.0/(Magick2PI*sigma*sigma));
1083 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1084 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1085 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1087 else /* limiting case - a unity (normalized Dirac) kernel */
1088 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1089 kernel->width*kernel->height*sizeof(*kernel->values));
1090 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1094 if ( type == DoGKernel )
1095 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1096 if ( sigma2 > MagickEpsilon )
1097 { sigma = sigma2; /* simplify loop expressions */
1098 A = 1.0/(2.0*sigma*sigma);
1099 B = (double) (1.0/(Magick2PI*sigma*sigma));
1100 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1101 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1102 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1104 else /* limiting case - a unity (normalized Dirac) kernel */
1105 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1108 if ( type == LoGKernel )
1109 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1110 if ( sigma > MagickEpsilon )
1111 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1112 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1113 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1114 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1115 { R = ((double)(u*u+v*v))*A;
1116 kernel->values[i] = (1-R)*exp(-R)*B;
1119 else /* special case - generate a unity kernel */
1120 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1121 kernel->width*kernel->height*sizeof(*kernel->values));
1122 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1126 /* Note the above kernels may have been 'clipped' by a user defined
1127 ** radius, producing a smaller (darker) kernel. Also for very small
1128 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1129 ** producing a very bright kernel.
1131 ** Normalization will still be needed.
1134 /* Normalize the 2D Gaussian Kernel
1136 ** NB: a CorrelateNormalize performs a normal Normalize if
1137 ** there are no negative values.
1139 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1140 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1146 sigma = fabs(args->sigma),
1149 if ( args->rho >= 1.0 )
1150 kernel->width = (size_t)args->rho*2+1;
1152 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1154 kernel->x = (ssize_t) (kernel->width-1)/2;
1156 kernel->negative_range = kernel->positive_range = 0.0;
1157 kernel->values=(MagickRealType *) MagickAssumeAligned(
1158 AcquireAlignedMemory(kernel->width,kernel->height*
1159 sizeof(*kernel->values)));
1160 if (kernel->values == (MagickRealType *) NULL)
1161 return(DestroyKernelInfo(kernel));
1164 #define KernelRank 3
1165 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1166 ** It generates a gaussian 3 times the width, and compresses it into
1167 ** the expected range. This produces a closer normalization of the
1168 ** resulting kernel, especially for very low sigma values.
1169 ** As such while wierd it is prefered.
1171 ** I am told this method originally came from Photoshop.
1173 ** A properly normalized curve is generated (apart from edge clipping)
1174 ** even though we later normalize the result (for edge clipping)
1175 ** to allow the correct generation of a "Difference of Blurs".
1179 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1180 (void) ResetMagickMemory(kernel->values,0, (size_t)
1181 kernel->width*kernel->height*sizeof(*kernel->values));
1182 /* Calculate a Positive 1D Gaussian */
1183 if ( sigma > MagickEpsilon )
1184 { sigma *= KernelRank; /* simplify loop expressions */
1185 alpha = 1.0/(2.0*sigma*sigma);
1186 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1187 for ( u=-v; u <= v; u++) {
1188 kernel->values[(u+v)/KernelRank] +=
1189 exp(-((double)(u*u))*alpha)*beta;
1192 else /* special case - generate a unity kernel */
1193 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1195 /* Direct calculation without curve averaging
1196 This is equivelent to a KernelRank of 1 */
1198 /* Calculate a Positive Gaussian */
1199 if ( sigma > MagickEpsilon )
1200 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1201 beta = 1.0/(MagickSQ2PI*sigma);
1202 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1203 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1205 else /* special case - generate a unity kernel */
1206 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1207 kernel->width*kernel->height*sizeof(*kernel->values));
1208 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1211 /* Note the above kernel may have been 'clipped' by a user defined
1212 ** radius, producing a smaller (darker) kernel. Also for very small
1213 ** sigma's (> 0.1) the central value becomes larger than one, as a
1214 ** result of not generating a actual 'discrete' kernel, and thus
1215 ** producing a very bright 'impulse'.
1217 ** Becuase of these two factors Normalization is required!
1220 /* Normalize the 1D Gaussian Kernel
1222 ** NB: a CorrelateNormalize performs a normal Normalize if
1223 ** there are no negative values.
1225 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1226 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1228 /* rotate the 1D kernel by given angle */
1229 RotateKernelInfo(kernel, args->xi );
1234 sigma = fabs(args->sigma),
1237 if ( args->rho < 1.0 )
1238 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1240 kernel->width = (size_t)args->rho;
1241 kernel->x = kernel->y = 0;
1243 kernel->negative_range = kernel->positive_range = 0.0;
1244 kernel->values=(MagickRealType *) MagickAssumeAligned(
1245 AcquireAlignedMemory(kernel->width,kernel->height*
1246 sizeof(*kernel->values)));
1247 if (kernel->values == (MagickRealType *) NULL)
1248 return(DestroyKernelInfo(kernel));
1250 /* A comet blur is half a 1D gaussian curve, so that the object is
1251 ** blurred in one direction only. This may not be quite the right
1252 ** curve to use so may change in the future. The function must be
1253 ** normalised after generation, which also resolves any clipping.
1255 ** As we are normalizing and not subtracting gaussians,
1256 ** there is no need for a divisor in the gaussian formula
1258 ** It is less comples
1260 if ( sigma > MagickEpsilon )
1263 #define KernelRank 3
1264 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1265 (void) ResetMagickMemory(kernel->values,0, (size_t)
1266 kernel->width*sizeof(*kernel->values));
1267 sigma *= KernelRank; /* simplify the loop expression */
1268 A = 1.0/(2.0*sigma*sigma);
1269 /* B = 1.0/(MagickSQ2PI*sigma); */
1270 for ( u=0; u < v; u++) {
1271 kernel->values[u/KernelRank] +=
1272 exp(-((double)(u*u))*A);
1273 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1275 for (i=0; i < (ssize_t) kernel->width; i++)
1276 kernel->positive_range += kernel->values[i];
1278 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1279 /* B = 1.0/(MagickSQ2PI*sigma); */
1280 for ( i=0; i < (ssize_t) kernel->width; i++)
1281 kernel->positive_range +=
1282 kernel->values[i] = exp(-((double)(i*i))*A);
1283 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1286 else /* special case - generate a unity kernel */
1287 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1288 kernel->width*kernel->height*sizeof(*kernel->values));
1289 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1290 kernel->positive_range = 1.0;
1293 kernel->minimum = 0.0;
1294 kernel->maximum = kernel->values[0];
1295 kernel->negative_range = 0.0;
1297 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1298 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1301 case BinomialKernel:
1306 if (args->rho < 1.0)
1307 kernel->width = kernel->height = 3; /* default radius = 1 */
1309 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1310 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1312 order_f = fact(kernel->width-1);
1314 kernel->values=(MagickRealType *) MagickAssumeAligned(
1315 AcquireAlignedMemory(kernel->width,kernel->height*
1316 sizeof(*kernel->values)));
1317 if (kernel->values == (MagickRealType *) NULL)
1318 return(DestroyKernelInfo(kernel));
1320 /* set all kernel values within diamond area to scale given */
1321 for ( i=0, v=0; v < (ssize_t)kernel->height; v++)
1323 alpha = order_f / ( fact((size_t) v) * fact(kernel->height-v-1) );
1324 for ( u=0; u < (ssize_t)kernel->width; u++, i++)
1325 kernel->positive_range += kernel->values[i] = (double)
1326 (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-u-1) ));
1328 kernel->minimum = 1.0;
1329 kernel->maximum = kernel->values[kernel->x+kernel->y*kernel->width];
1330 kernel->negative_range = 0.0;
1335 Convolution Kernels - Well Known Named Constant Kernels
1337 case LaplacianKernel:
1338 { switch ( (int) args->rho ) {
1340 default: /* laplacian square filter -- default */
1341 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1343 case 1: /* laplacian diamond filter */
1344 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1347 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1350 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1352 case 5: /* a 5x5 laplacian */
1353 kernel=ParseKernelArray(
1354 "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");
1356 case 7: /* a 7x7 laplacian */
1357 kernel=ParseKernelArray(
1358 "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" );
1360 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1361 kernel=ParseKernelArray(
1362 "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");
1364 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1365 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1366 kernel=ParseKernelArray(
1367 "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");
1370 if (kernel == (KernelInfo *) NULL)
1372 kernel->type = type;
1376 { /* Simple Sobel Kernel */
1377 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1378 if (kernel == (KernelInfo *) NULL)
1380 kernel->type = type;
1381 RotateKernelInfo(kernel, args->rho);
1386 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1387 if (kernel == (KernelInfo *) NULL)
1389 kernel->type = type;
1390 RotateKernelInfo(kernel, args->rho);
1395 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1396 if (kernel == (KernelInfo *) NULL)
1398 kernel->type = type;
1399 RotateKernelInfo(kernel, args->rho);
1404 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1405 if (kernel == (KernelInfo *) NULL)
1407 kernel->type = type;
1408 RotateKernelInfo(kernel, args->rho);
1413 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1414 if (kernel == (KernelInfo *) NULL)
1416 kernel->type = type;
1417 RotateKernelInfo(kernel, args->rho);
1420 case FreiChenKernel:
1421 /* Direction is set to be left to right positive */
1422 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1423 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1424 { switch ( (int) args->rho ) {
1427 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1428 if (kernel == (KernelInfo *) NULL)
1430 kernel->type = type;
1431 kernel->values[3] = +(MagickRealType) MagickSQ2;
1432 kernel->values[5] = -(MagickRealType) MagickSQ2;
1433 CalcKernelMetaData(kernel); /* recalculate meta-data */
1436 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1437 if (kernel == (KernelInfo *) NULL)
1439 kernel->type = type;
1440 kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2;
1441 kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2;
1442 CalcKernelMetaData(kernel); /* recalculate meta-data */
1443 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1446 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1447 if (kernel == (KernelInfo *) NULL)
1452 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1453 if (kernel == (KernelInfo *) NULL)
1455 kernel->type = type;
1456 kernel->values[3] = +(MagickRealType) MagickSQ2;
1457 kernel->values[5] = -(MagickRealType) MagickSQ2;
1458 CalcKernelMetaData(kernel); /* recalculate meta-data */
1459 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1462 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1463 if (kernel == (KernelInfo *) NULL)
1465 kernel->type = type;
1466 kernel->values[1] = +(MagickRealType) MagickSQ2;
1467 kernel->values[7] = +(MagickRealType) MagickSQ2;
1468 CalcKernelMetaData(kernel);
1469 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1472 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1473 if (kernel == (KernelInfo *) NULL)
1475 kernel->type = type;
1476 kernel->values[0] = +(MagickRealType) MagickSQ2;
1477 kernel->values[8] = -(MagickRealType) MagickSQ2;
1478 CalcKernelMetaData(kernel);
1479 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1482 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1483 if (kernel == (KernelInfo *) NULL)
1485 kernel->type = type;
1486 kernel->values[2] = -(MagickRealType) MagickSQ2;
1487 kernel->values[6] = +(MagickRealType) MagickSQ2;
1488 CalcKernelMetaData(kernel);
1489 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1492 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1493 if (kernel == (KernelInfo *) NULL)
1495 kernel->type = type;
1496 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1499 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1500 if (kernel == (KernelInfo *) NULL)
1502 kernel->type = type;
1503 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1506 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1507 if (kernel == (KernelInfo *) NULL)
1509 kernel->type = type;
1510 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1513 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1514 if (kernel == (KernelInfo *) NULL)
1516 kernel->type = type;
1517 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1520 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1521 if (kernel == (KernelInfo *) NULL)
1523 kernel->type = type;
1524 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1527 if ( fabs(args->sigma) >= MagickEpsilon )
1528 /* Rotate by correctly supplied 'angle' */
1529 RotateKernelInfo(kernel, args->sigma);
1530 else if ( args->rho > 30.0 || args->rho < -30.0 )
1531 /* Rotate by out of bounds 'type' */
1532 RotateKernelInfo(kernel, args->rho);
1537 Boolean or Shaped Kernels
1541 if (args->rho < 1.0)
1542 kernel->width = kernel->height = 3; /* default radius = 1 */
1544 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1545 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1547 kernel->values=(MagickRealType *) MagickAssumeAligned(
1548 AcquireAlignedMemory(kernel->width,kernel->height*
1549 sizeof(*kernel->values)));
1550 if (kernel->values == (MagickRealType *) NULL)
1551 return(DestroyKernelInfo(kernel));
1553 /* set all kernel values within diamond area to scale given */
1554 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1555 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1556 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1557 kernel->positive_range += kernel->values[i] = args->sigma;
1559 kernel->values[i] = nan;
1560 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1564 case RectangleKernel:
1567 if ( type == SquareKernel )
1569 if (args->rho < 1.0)
1570 kernel->width = kernel->height = 3; /* default radius = 1 */
1572 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1573 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1574 scale = args->sigma;
1577 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1578 if ( args->rho < 1.0 || args->sigma < 1.0 )
1579 return(DestroyKernelInfo(kernel)); /* invalid args given */
1580 kernel->width = (size_t)args->rho;
1581 kernel->height = (size_t)args->sigma;
1582 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1583 args->psi < 0.0 || args->psi > (double)kernel->height )
1584 return(DestroyKernelInfo(kernel)); /* invalid args given */
1585 kernel->x = (ssize_t) args->xi;
1586 kernel->y = (ssize_t) args->psi;
1589 kernel->values=(MagickRealType *) MagickAssumeAligned(
1590 AcquireAlignedMemory(kernel->width,kernel->height*
1591 sizeof(*kernel->values)));
1592 if (kernel->values == (MagickRealType *) NULL)
1593 return(DestroyKernelInfo(kernel));
1595 /* set all kernel values to scale given */
1596 u=(ssize_t) (kernel->width*kernel->height);
1597 for ( i=0; i < u; i++)
1598 kernel->values[i] = scale;
1599 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1600 kernel->positive_range = scale*u;
1605 if (args->rho < 1.0)
1606 kernel->width = kernel->height = 5; /* default radius = 2 */
1608 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1609 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1611 kernel->values=(MagickRealType *) MagickAssumeAligned(
1612 AcquireAlignedMemory(kernel->width,kernel->height*
1613 sizeof(*kernel->values)));
1614 if (kernel->values == (MagickRealType *) NULL)
1615 return(DestroyKernelInfo(kernel));
1617 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1618 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1619 if ( (labs((long) u)+labs((long) v)) <=
1620 ((long)kernel->x + (long)(kernel->x/2)) )
1621 kernel->positive_range += kernel->values[i] = args->sigma;
1623 kernel->values[i] = nan;
1624 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1630 limit = (ssize_t)(args->rho*args->rho);
1632 if (args->rho < 0.4) /* default radius approx 4.3 */
1633 kernel->width = kernel->height = 9L, limit = 18L;
1635 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1636 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1638 kernel->values=(MagickRealType *) MagickAssumeAligned(
1639 AcquireAlignedMemory(kernel->width,kernel->height*
1640 sizeof(*kernel->values)));
1641 if (kernel->values == (MagickRealType *) NULL)
1642 return(DestroyKernelInfo(kernel));
1644 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1645 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1646 if ((u*u+v*v) <= limit)
1647 kernel->positive_range += kernel->values[i] = args->sigma;
1649 kernel->values[i] = nan;
1650 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1655 if (args->rho < 1.0)
1656 kernel->width = kernel->height = 5; /* default radius 2 */
1658 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1659 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1661 kernel->values=(MagickRealType *) MagickAssumeAligned(
1662 AcquireAlignedMemory(kernel->width,kernel->height*
1663 sizeof(*kernel->values)));
1664 if (kernel->values == (MagickRealType *) NULL)
1665 return(DestroyKernelInfo(kernel));
1667 /* set all kernel values along axises to given scale */
1668 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1669 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1670 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1671 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1672 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1677 if (args->rho < 1.0)
1678 kernel->width = kernel->height = 5; /* default radius 2 */
1680 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1681 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1683 kernel->values=(MagickRealType *) MagickAssumeAligned(
1684 AcquireAlignedMemory(kernel->width,kernel->height*
1685 sizeof(*kernel->values)));
1686 if (kernel->values == (MagickRealType *) NULL)
1687 return(DestroyKernelInfo(kernel));
1689 /* set all kernel values along axises to given scale */
1690 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1691 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1692 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1693 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1694 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1708 if (args->rho < args->sigma)
1710 kernel->width = ((size_t)args->sigma)*2+1;
1711 limit1 = (ssize_t)(args->rho*args->rho);
1712 limit2 = (ssize_t)(args->sigma*args->sigma);
1716 kernel->width = ((size_t)args->rho)*2+1;
1717 limit1 = (ssize_t)(args->sigma*args->sigma);
1718 limit2 = (ssize_t)(args->rho*args->rho);
1721 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1723 kernel->height = kernel->width;
1724 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1725 kernel->values=(MagickRealType *) MagickAssumeAligned(
1726 AcquireAlignedMemory(kernel->width,kernel->height*
1727 sizeof(*kernel->values)));
1728 if (kernel->values == (MagickRealType *) NULL)
1729 return(DestroyKernelInfo(kernel));
1731 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1732 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1733 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1734 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1735 { ssize_t radius=u*u+v*v;
1736 if (limit1 < radius && radius <= limit2)
1737 kernel->positive_range += kernel->values[i] = (double) scale;
1739 kernel->values[i] = nan;
1741 kernel->minimum = kernel->maximum = (double) scale;
1742 if ( type == PeaksKernel ) {
1743 /* set the central point in the middle */
1744 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1745 kernel->positive_range = 1.0;
1746 kernel->maximum = 1.0;
1752 kernel=AcquireKernelInfo("ThinSE:482");
1753 if (kernel == (KernelInfo *) NULL)
1755 kernel->type = type;
1756 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1761 kernel=AcquireKernelInfo("ThinSE:87");
1762 if (kernel == (KernelInfo *) NULL)
1764 kernel->type = type;
1765 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1768 case DiagonalsKernel:
1770 switch ( (int) args->rho ) {
1775 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1776 if (kernel == (KernelInfo *) NULL)
1778 kernel->type = type;
1779 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1780 if (new_kernel == (KernelInfo *) NULL)
1781 return(DestroyKernelInfo(kernel));
1782 new_kernel->type = type;
1783 LastKernelInfo(kernel)->next = new_kernel;
1784 ExpandMirrorKernelInfo(kernel);
1788 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1791 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1794 if (kernel == (KernelInfo *) NULL)
1796 kernel->type = type;
1797 RotateKernelInfo(kernel, args->sigma);
1800 case LineEndsKernel:
1801 { /* Kernels for finding the end of thin lines */
1802 switch ( (int) args->rho ) {
1805 /* set of kernels to find all end of lines */
1806 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1808 /* kernel for 4-connected line ends - no rotation */
1809 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1812 /* kernel to add for 8-connected lines - no rotation */
1813 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1816 /* kernel to add for orthogonal line ends - does not find corners */
1817 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1820 /* traditional line end - fails on last T end */
1821 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1824 if (kernel == (KernelInfo *) NULL)
1826 kernel->type = type;
1827 RotateKernelInfo(kernel, args->sigma);
1830 case LineJunctionsKernel:
1831 { /* kernels for finding the junctions of multiple lines */
1832 switch ( (int) args->rho ) {
1835 /* set of kernels to find all line junctions */
1836 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1839 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1842 /* Diagonal T Junctions */
1843 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1846 /* Orthogonal T Junctions */
1847 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1850 /* Diagonal X Junctions */
1851 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1854 /* Orthogonal X Junctions - minimal diamond kernel */
1855 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1858 if (kernel == (KernelInfo *) NULL)
1860 kernel->type = type;
1861 RotateKernelInfo(kernel, args->sigma);
1865 { /* Ridges - Ridge finding kernels */
1868 switch ( (int) args->rho ) {
1871 kernel=ParseKernelArray("3x1:0,1,0");
1872 if (kernel == (KernelInfo *) NULL)
1874 kernel->type = type;
1875 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1878 kernel=ParseKernelArray("4x1:0,1,1,0");
1879 if (kernel == (KernelInfo *) NULL)
1881 kernel->type = type;
1882 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1884 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1885 /* Unfortunatally we can not yet rotate a non-square kernel */
1886 /* But then we can't flip a non-symetrical kernel either */
1887 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1888 if (new_kernel == (KernelInfo *) NULL)
1889 return(DestroyKernelInfo(kernel));
1890 new_kernel->type = type;
1891 LastKernelInfo(kernel)->next = new_kernel;
1892 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1893 if (new_kernel == (KernelInfo *) NULL)
1894 return(DestroyKernelInfo(kernel));
1895 new_kernel->type = type;
1896 LastKernelInfo(kernel)->next = new_kernel;
1897 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1898 if (new_kernel == (KernelInfo *) NULL)
1899 return(DestroyKernelInfo(kernel));
1900 new_kernel->type = type;
1901 LastKernelInfo(kernel)->next = new_kernel;
1902 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1903 if (new_kernel == (KernelInfo *) NULL)
1904 return(DestroyKernelInfo(kernel));
1905 new_kernel->type = type;
1906 LastKernelInfo(kernel)->next = new_kernel;
1907 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1908 if (new_kernel == (KernelInfo *) NULL)
1909 return(DestroyKernelInfo(kernel));
1910 new_kernel->type = type;
1911 LastKernelInfo(kernel)->next = new_kernel;
1912 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1913 if (new_kernel == (KernelInfo *) NULL)
1914 return(DestroyKernelInfo(kernel));
1915 new_kernel->type = type;
1916 LastKernelInfo(kernel)->next = new_kernel;
1917 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1918 if (new_kernel == (KernelInfo *) NULL)
1919 return(DestroyKernelInfo(kernel));
1920 new_kernel->type = type;
1921 LastKernelInfo(kernel)->next = new_kernel;
1922 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1923 if (new_kernel == (KernelInfo *) NULL)
1924 return(DestroyKernelInfo(kernel));
1925 new_kernel->type = type;
1926 LastKernelInfo(kernel)->next = new_kernel;
1931 case ConvexHullKernel:
1935 /* first set of 8 kernels */
1936 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1937 if (kernel == (KernelInfo *) NULL)
1939 kernel->type = type;
1940 ExpandRotateKernelInfo(kernel, 90.0);
1941 /* append the mirror versions too - no flip function yet */
1942 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1943 if (new_kernel == (KernelInfo *) NULL)
1944 return(DestroyKernelInfo(kernel));
1945 new_kernel->type = type;
1946 ExpandRotateKernelInfo(new_kernel, 90.0);
1947 LastKernelInfo(kernel)->next = new_kernel;
1950 case SkeletonKernel:
1952 switch ( (int) args->rho ) {
1955 /* Traditional Skeleton...
1956 ** A cyclically rotated single kernel
1958 kernel=AcquireKernelInfo("ThinSE:482");
1959 if (kernel == (KernelInfo *) NULL)
1961 kernel->type = type;
1962 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1965 /* HIPR Variation of the cyclic skeleton
1966 ** Corners of the traditional method made more forgiving,
1967 ** but the retain the same cyclic order.
1969 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1970 if (kernel == (KernelInfo *) NULL)
1972 if (kernel->next == (KernelInfo *) NULL)
1973 return(DestroyKernelInfo(kernel));
1974 kernel->type = type;
1975 kernel->next->type = type;
1976 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1979 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1980 ** "Connectivity-Preserving Morphological Image Thransformations"
1981 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1982 ** http://www.leptonica.com/papers/conn.pdf
1984 kernel=AcquireKernelInfo(
1985 "ThinSE:41; ThinSE:42; ThinSE:43");
1986 if (kernel == (KernelInfo *) NULL)
1988 kernel->type = type;
1989 kernel->next->type = type;
1990 kernel->next->next->type = type;
1991 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1997 { /* Special kernels for general thinning, while preserving connections
1998 ** "Connectivity-Preserving Morphological Image Thransformations"
1999 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
2000 ** http://www.leptonica.com/papers/conn.pdf
2002 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2004 ** Note kernels do not specify the origin pixel, allowing them
2005 ** to be used for both thickening and thinning operations.
2007 switch ( (int) args->rho ) {
2008 /* SE for 4-connected thinning */
2009 case 41: /* SE_4_1 */
2010 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
2012 case 42: /* SE_4_2 */
2013 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
2015 case 43: /* SE_4_3 */
2016 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
2018 case 44: /* SE_4_4 */
2019 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
2021 case 45: /* SE_4_5 */
2022 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
2024 case 46: /* SE_4_6 */
2025 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
2027 case 47: /* SE_4_7 */
2028 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
2030 case 48: /* SE_4_8 */
2031 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
2033 case 49: /* SE_4_9 */
2034 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
2036 /* SE for 8-connected thinning - negatives of the above */
2037 case 81: /* SE_8_0 */
2038 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
2040 case 82: /* SE_8_2 */
2041 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
2043 case 83: /* SE_8_3 */
2044 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
2046 case 84: /* SE_8_4 */
2047 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
2049 case 85: /* SE_8_5 */
2050 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2052 case 86: /* SE_8_6 */
2053 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2055 case 87: /* SE_8_7 */
2056 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2058 case 88: /* SE_8_8 */
2059 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2061 case 89: /* SE_8_9 */
2062 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2064 /* Special combined SE kernels */
2065 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2066 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2068 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2069 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2071 case 481: /* SE_48_1 - General Connected Corner Kernel */
2072 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2075 case 482: /* SE_48_2 - General Edge Kernel */
2076 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2079 if (kernel == (KernelInfo *) NULL)
2081 kernel->type = type;
2082 RotateKernelInfo(kernel, args->sigma);
2086 Distance Measuring Kernels
2088 case ChebyshevKernel:
2090 if (args->rho < 1.0)
2091 kernel->width = kernel->height = 3; /* default radius = 1 */
2093 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2094 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2096 kernel->values=(MagickRealType *) MagickAssumeAligned(
2097 AcquireAlignedMemory(kernel->width,kernel->height*
2098 sizeof(*kernel->values)));
2099 if (kernel->values == (MagickRealType *) NULL)
2100 return(DestroyKernelInfo(kernel));
2102 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2103 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2104 kernel->positive_range += ( kernel->values[i] =
2105 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2106 kernel->maximum = kernel->values[0];
2109 case ManhattanKernel:
2111 if (args->rho < 1.0)
2112 kernel->width = kernel->height = 3; /* default radius = 1 */
2114 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2115 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2117 kernel->values=(MagickRealType *) MagickAssumeAligned(
2118 AcquireAlignedMemory(kernel->width,kernel->height*
2119 sizeof(*kernel->values)));
2120 if (kernel->values == (MagickRealType *) NULL)
2121 return(DestroyKernelInfo(kernel));
2123 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2124 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2125 kernel->positive_range += ( kernel->values[i] =
2126 args->sigma*(labs((long) u)+labs((long) v)) );
2127 kernel->maximum = kernel->values[0];
2130 case OctagonalKernel:
2132 if (args->rho < 2.0)
2133 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2135 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2136 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2138 kernel->values=(MagickRealType *) MagickAssumeAligned(
2139 AcquireAlignedMemory(kernel->width,kernel->height*
2140 sizeof(*kernel->values)));
2141 if (kernel->values == (MagickRealType *) NULL)
2142 return(DestroyKernelInfo(kernel));
2144 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2145 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2148 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2149 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2150 kernel->positive_range += kernel->values[i] =
2151 args->sigma*MagickMax(r1,r2);
2153 kernel->maximum = kernel->values[0];
2156 case EuclideanKernel:
2158 if (args->rho < 1.0)
2159 kernel->width = kernel->height = 3; /* default radius = 1 */
2161 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2162 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2164 kernel->values=(MagickRealType *) MagickAssumeAligned(
2165 AcquireAlignedMemory(kernel->width,kernel->height*
2166 sizeof(*kernel->values)));
2167 if (kernel->values == (MagickRealType *) NULL)
2168 return(DestroyKernelInfo(kernel));
2170 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2171 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2172 kernel->positive_range += ( kernel->values[i] =
2173 args->sigma*sqrt((double)(u*u+v*v)) );
2174 kernel->maximum = kernel->values[0];
2179 /* No-Op Kernel - Basically just a single pixel on its own */
2180 kernel=ParseKernelArray("1:1");
2181 if (kernel == (KernelInfo *) NULL)
2183 kernel->type = UndefinedKernel;
2192 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2196 % C l o n e K e r n e l I n f o %
2200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2202 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2203 % can be modified without effecting the original. The cloned kernel should
2204 % be destroyed using DestoryKernelInfo() when no longer needed.
2206 % The format of the CloneKernelInfo method is:
2208 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2210 % A description of each parameter follows:
2212 % o kernel: the Morphology/Convolution kernel to be cloned
2215 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2223 assert(kernel != (KernelInfo *) NULL);
2224 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2225 if (new_kernel == (KernelInfo *) NULL)
2227 *new_kernel=(*kernel); /* copy values in structure */
2229 /* replace the values with a copy of the values */
2230 new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2231 AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2232 if (new_kernel->values == (MagickRealType *) NULL)
2233 return(DestroyKernelInfo(new_kernel));
2234 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2235 new_kernel->values[i]=kernel->values[i];
2237 /* Also clone the next kernel in the kernel list */
2238 if ( kernel->next != (KernelInfo *) NULL ) {
2239 new_kernel->next = CloneKernelInfo(kernel->next);
2240 if ( new_kernel->next == (KernelInfo *) NULL )
2241 return(DestroyKernelInfo(new_kernel));
2248 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2252 % D e s t r o y K e r n e l I n f o %
2256 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2258 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2261 % The format of the DestroyKernelInfo method is:
2263 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2265 % A description of each parameter follows:
2267 % o kernel: the Morphology/Convolution kernel to be destroyed
2270 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2272 assert(kernel != (KernelInfo *) NULL);
2273 if ( kernel->next != (KernelInfo *) NULL )
2274 kernel->next=DestroyKernelInfo(kernel->next);
2275 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2276 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2281 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2285 + E x p a n d M i r r o r K e r n e l I n f o %
2289 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2291 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2292 % sequence of 90-degree rotated kernels but providing a reflected 180
2293 % rotatation, before the -/+ 90-degree rotations.
2295 % This special rotation order produces a better, more symetrical thinning of
2298 % The format of the ExpandMirrorKernelInfo method is:
2300 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2302 % A description of each parameter follows:
2304 % o kernel: the Morphology/Convolution kernel
2306 % This function is only internel to this module, as it is not finalized,
2307 % especially with regard to non-orthogonal angles, and rotation of larger
2312 static void FlopKernelInfo(KernelInfo *kernel)
2313 { /* Do a Flop by reversing each row. */
2321 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2322 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2323 t=k[x], k[x]=k[r], k[r]=t;
2325 kernel->x = kernel->width - kernel->x - 1;
2326 angle = fmod(angle+180.0, 360.0);
2330 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2338 clone = CloneKernelInfo(last);
2339 RotateKernelInfo(clone, 180); /* flip */
2340 LastKernelInfo(last)->next = clone;
2343 clone = CloneKernelInfo(last);
2344 RotateKernelInfo(clone, 90); /* transpose */
2345 LastKernelInfo(last)->next = clone;
2348 clone = CloneKernelInfo(last);
2349 RotateKernelInfo(clone, 180); /* flop */
2350 LastKernelInfo(last)->next = clone;
2356 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2360 + E x p a n d R o t a t e K e r n e l I n f o %
2364 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2366 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2367 % incrementally by the angle given, until the kernel repeats.
2369 % WARNING: 45 degree rotations only works for 3x3 kernels.
2370 % While 90 degree roatations only works for linear and square kernels
2372 % The format of the ExpandRotateKernelInfo method is:
2374 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2376 % A description of each parameter follows:
2378 % o kernel: the Morphology/Convolution kernel
2380 % o angle: angle to rotate in degrees
2382 % This function is only internel to this module, as it is not finalized,
2383 % especially with regard to non-orthogonal angles, and rotation of larger
2387 /* Internal Routine - Return true if two kernels are the same */
2388 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2389 const KernelInfo *kernel2)
2394 /* check size and origin location */
2395 if ( kernel1->width != kernel2->width
2396 || kernel1->height != kernel2->height
2397 || kernel1->x != kernel2->x
2398 || kernel1->y != kernel2->y )
2401 /* check actual kernel values */
2402 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2403 /* Test for Nan equivalence */
2404 if ( IfNaN(kernel1->values[i]) && !IfNaN(kernel2->values[i]) )
2406 if ( IfNaN(kernel2->values[i]) && !IfNaN(kernel1->values[i]) )
2408 /* Test actual values are equivalent */
2409 if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2416 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2423 DisableMSCWarning(4127)
2426 clone = CloneKernelInfo(last);
2427 RotateKernelInfo(clone, angle);
2428 if ( SameKernelInfo(kernel, clone) != MagickFalse )
2430 LastKernelInfo(last)->next = clone;
2433 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2438 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2442 + C a l c M e t a K e r n a l I n f o %
2446 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2448 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2449 % using the kernel values. This should only ne used if it is not possible to
2450 % calculate that meta-data in some easier way.
2452 % It is important that the meta-data is correct before ScaleKernelInfo() is
2453 % used to perform kernel normalization.
2455 % The format of the CalcKernelMetaData method is:
2457 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2459 % A description of each parameter follows:
2461 % o kernel: the Morphology/Convolution kernel to modify
2463 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2464 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2465 % however is not true for flat-shaped morphological kernels.
2467 % WARNING: Only the specific kernel pointed to is modified, not a list of
2470 % This is an internal function and not expected to be useful outside this
2471 % module. This could change however.
2473 static void CalcKernelMetaData(KernelInfo *kernel)
2478 kernel->minimum = kernel->maximum = 0.0;
2479 kernel->negative_range = kernel->positive_range = 0.0;
2480 for (i=0; i < (kernel->width*kernel->height); i++)
2482 if ( fabs(kernel->values[i]) < MagickEpsilon )
2483 kernel->values[i] = 0.0;
2484 ( kernel->values[i] < 0)
2485 ? ( kernel->negative_range += kernel->values[i] )
2486 : ( kernel->positive_range += kernel->values[i] );
2487 Minimize(kernel->minimum, kernel->values[i]);
2488 Maximize(kernel->maximum, kernel->values[i]);
2495 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2499 % M o r p h o l o g y A p p l y %
2503 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2505 % MorphologyApply() applies a morphological method, multiple times using
2506 % a list of multiple kernels. This is the method that should be called by
2507 % other 'operators' that internally use morphology operations as part of
2510 % It is basically equivalent to as MorphologyImage() (see below) but without
2511 % any user controls. This allows internel programs to use this method to
2512 % perform a specific task without possible interference by any API user
2513 % supplied settings.
2515 % It is MorphologyImage() task to extract any such user controls, and
2516 % pass them to this function for processing.
2518 % More specifically all given kernels should already be scaled, normalised,
2519 % and blended appropriatally before being parred to this routine. The
2520 % appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2522 % The format of the MorphologyApply method is:
2524 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2525 % const ssize_t iterations,const KernelInfo *kernel,
2526 % const CompositeMethod compose,const double bias,
2527 % ExceptionInfo *exception)
2529 % A description of each parameter follows:
2531 % o image: the source image
2533 % o method: the morphology method to be applied.
2535 % o iterations: apply the operation this many times (or no change).
2536 % A value of -1 means loop until no change found.
2537 % How this is applied may depend on the morphology method.
2538 % Typically this is a value of 1.
2540 % o channel: the channel type.
2542 % o kernel: An array of double representing the morphology kernel.
2544 % o compose: How to handle or merge multi-kernel results.
2545 % If 'UndefinedCompositeOp' use default for the Morphology method.
2546 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2547 % Otherwise merge the results using the compose method given.
2549 % o bias: Convolution Output Bias.
2551 % o exception: return any errors or warnings in this structure.
2554 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2555 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2556 ExceptionInfo *exception)
2558 #define MorphologyTag "Morphology/Image"
2584 assert(image != (Image *) NULL);
2585 assert(image->signature == MagickSignature);
2586 assert(morphology_image != (Image *) NULL);
2587 assert(morphology_image->signature == MagickSignature);
2588 assert(kernel != (KernelInfo *) NULL);
2589 assert(kernel->signature == MagickSignature);
2590 assert(exception != (ExceptionInfo *) NULL);
2591 assert(exception->signature == MagickSignature);
2594 image_view=AcquireVirtualCacheView(image,exception);
2595 morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2596 width=image->columns+kernel->width-1;
2601 case ConvolveMorphology:
2602 case DilateMorphology:
2603 case DilateIntensityMorphology:
2604 case IterativeDistanceMorphology:
2607 Kernel needs to used with reflection about origin.
2609 offset.x=(ssize_t) kernel->width-kernel->x-1;
2610 offset.y=(ssize_t) kernel->height-kernel->y-1;
2613 case ErodeMorphology:
2614 case ErodeIntensityMorphology:
2615 case HitAndMissMorphology:
2616 case ThinningMorphology:
2617 case ThickenMorphology:
2625 assert("Not a Primitive Morphology Method" != (char *) NULL);
2630 changes=(size_t *) AcquireQuantumMemory(GetOpenMPMaximumThreads(),
2632 if (changes == (size_t *) NULL)
2633 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2634 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2636 if ((method == ConvolveMorphology) && (kernel->width == 1))
2639 id = GetOpenMPThreadId();
2645 Special handling (for speed) of vertical (blur) kernels. This performs
2646 its handling in columns rather than in rows. This is only done
2647 for convolve as it is the only method that generates very large 1-D
2648 vertical kernels (such as a 'BlurKernel')
2650 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2651 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2652 magick_threads(image,morphology_image,image->columns,1)
2654 for (x=0; x < (ssize_t) image->columns; x++)
2656 register const Quantum
2668 if (status == MagickFalse)
2670 p=GetCacheViewVirtualPixels(image_view,x,-offset.y,1,image->rows+
2671 kernel->height-1,exception);
2672 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2673 morphology_image->rows,exception);
2674 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2679 center=(ssize_t) GetPixelChannels(image)*offset.y;
2680 for (y=0; y < (ssize_t) image->rows; y++)
2685 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2699 register const MagickRealType
2702 register const Quantum
2714 channel=GetPixelChannelChannel(image,i);
2715 traits=GetPixelChannelTraits(image,channel);
2716 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2717 if ((traits == UndefinedPixelTrait) ||
2718 (morphology_traits == UndefinedPixelTrait))
2720 if (((morphology_traits & CopyPixelTrait) != 0) ||
2721 (GetPixelReadMask(image,p+center) == 0))
2723 SetPixelChannel(morphology_image,channel,p[center+i],q);
2726 k=(&kernel->values[kernel->width*kernel->height-1]);
2730 if ((morphology_traits & BlendPixelTrait) == 0)
2735 for (v=0; v < (ssize_t) kernel->height; v++)
2737 for (u=0; u < (ssize_t) kernel->width; u++)
2739 if (IfNaN(*k) == MagickFalse)
2741 pixel+=(*k)*pixels[i];
2745 pixels+=GetPixelChannels(image);
2748 if (fabs(pixel-p[center+i]) > MagickEpsilon)
2750 gamma=(double) kernel->height*kernel->width/count;
2751 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2759 for (v=0; v < (ssize_t) kernel->height; v++)
2761 for (u=0; u < (ssize_t) kernel->width; u++)
2763 if (IfNaN(*k) == MagickFalse)
2765 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2766 pixel+=(*k)*alpha*pixels[i];
2771 pixels+=GetPixelChannels(image);
2774 if (fabs(pixel-p[center+i]) > MagickEpsilon)
2776 gamma=PerceptibleReciprocal(gamma);
2777 gamma*=(double) kernel->height*kernel->width/count;
2778 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2781 p+=GetPixelChannels(image);
2782 q+=GetPixelChannels(morphology_image);
2784 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2786 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2791 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2792 #pragma omp critical (MagickCore_MorphologyPrimitive)
2794 proceed=SetImageProgress(image,MorphologyTag,progress++,
2796 if (proceed == MagickFalse)
2800 morphology_image->type=image->type;
2801 morphology_view=DestroyCacheView(morphology_view);
2802 image_view=DestroyCacheView(image_view);
2803 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2804 changed+=changes[i];
2805 changes=(size_t *) RelinquishMagickMemory(changes);
2806 return(status ? (ssize_t) changed : 0);
2809 Normal handling of horizontal or rectangular kernels (row by row).
2811 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2812 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2813 magick_threads(image,morphology_image,image->rows,1)
2815 for (y=0; y < (ssize_t) image->rows; y++)
2818 id = GetOpenMPThreadId();
2820 register const Quantum
2832 if (status == MagickFalse)
2834 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,
2835 kernel->height,exception);
2836 q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2838 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2843 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
2844 GetPixelChannels(image)*offset.x);
2845 for (x=0; x < (ssize_t) image->columns; x++)
2850 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2866 register const MagickRealType
2869 register const Quantum
2881 channel=GetPixelChannelChannel(image,i);
2882 traits=GetPixelChannelTraits(image,channel);
2883 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2884 if ((traits == UndefinedPixelTrait) ||
2885 (morphology_traits == UndefinedPixelTrait))
2887 if (((morphology_traits & CopyPixelTrait) != 0) ||
2888 (GetPixelReadMask(image,p+center) == 0))
2890 SetPixelChannel(morphology_image,channel,p[center+i],q);
2895 minimum=(double) QuantumRange;
2896 count=kernel->width*kernel->height;
2899 case ConvolveMorphology: pixel=bias; break;
2900 case HitAndMissMorphology: pixel=(double) QuantumRange; break;
2901 case ThinningMorphology: pixel=(double) QuantumRange; break;
2902 case ThickenMorphology: pixel=(double) QuantumRange; break;
2903 case ErodeMorphology: pixel=(double) QuantumRange; break;
2904 case DilateMorphology: pixel=0.0; break;
2905 case ErodeIntensityMorphology:
2906 case DilateIntensityMorphology:
2907 case IterativeDistanceMorphology:
2909 pixel=(double) p[center+i];
2912 default: pixel=0; break;
2917 case ConvolveMorphology:
2920 Weighted Average of pixels using reflected kernel
2922 For correct working of this operation for asymetrical
2923 kernels, the kernel needs to be applied in its reflected form.
2924 That is its values needs to be reversed.
2926 Correlation is actually the same as this but without reflecting
2927 the kernel, and thus 'lower-level' that Convolution. However
2928 as Convolution is the more common method used, and it does not
2929 really cost us much in terms of processing to use a reflected
2930 kernel, so it is Convolution that is implemented.
2932 Correlation will have its kernel reflected before calling
2933 this function to do a Convolve.
2935 For more details of Correlation vs Convolution see
2936 http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2938 k=(&kernel->values[kernel->width*kernel->height-1]);
2940 if ((morphology_traits & BlendPixelTrait) == 0)
2945 for (v=0; v < (ssize_t) kernel->height; v++)
2947 for (u=0; u < (ssize_t) kernel->width; u++)
2949 if (IfNaN(*k) == MagickFalse)
2951 pixel+=(*k)*pixels[i];
2955 pixels+=GetPixelChannels(image);
2957 pixels+=(image->columns-1)*GetPixelChannels(image);
2964 for (v=0; v < (ssize_t) kernel->height; v++)
2966 for (u=0; u < (ssize_t) kernel->width; u++)
2968 if (IfNaN(*k) == MagickFalse)
2970 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2971 pixel+=(*k)*alpha*pixels[i];
2976 pixels+=GetPixelChannels(image);
2978 pixels+=(image->columns-1)*GetPixelChannels(image);
2982 case ErodeMorphology:
2985 Minimum value within kernel neighbourhood.
2987 The kernel is not reflected for this operation. In normal
2988 Greyscale Morphology, the kernel value should be added
2989 to the real value, this is currently not done, due to the
2990 nature of the boolean kernels being used.
2993 for (v=0; v < (ssize_t) kernel->height; v++)
2995 for (u=0; u < (ssize_t) kernel->width; u++)
2997 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
2999 if ((double) pixels[i] < pixel)
3000 pixel=(double) pixels[i];
3003 pixels+=GetPixelChannels(image);
3005 pixels+=(image->columns-1)*GetPixelChannels(image);
3009 case DilateMorphology:
3012 Maximum value within kernel neighbourhood.
3014 For correct working of this operation for asymetrical kernels,
3015 the kernel needs to be applied in its reflected form. That is
3016 its values needs to be reversed.
3018 In normal Greyscale Morphology, the kernel value should be
3019 added to the real value, this is currently not done, due to the
3020 nature of the boolean kernels being used.
3023 k=(&kernel->values[kernel->width*kernel->height-1]);
3024 for (v=0; v < (ssize_t) kernel->height; v++)
3026 for (u=0; u < (ssize_t) kernel->width; u++)
3028 if ((IfNaN(*k) == MagickFalse) && (*k > 0.5))
3030 if ((double) pixels[i] > pixel)
3031 pixel=(double) pixels[i];
3035 pixels+=GetPixelChannels(image);
3037 pixels+=(image->columns-1)*GetPixelChannels(image);
3041 case HitAndMissMorphology:
3042 case ThinningMorphology:
3043 case ThickenMorphology:
3046 Minimum of foreground pixel minus maxumum of background pixels.
3048 The kernel is not reflected for this operation, and consists
3049 of both foreground and background pixel neighbourhoods, 0.0 for
3050 background, and 1.0 for foreground with either Nan or 0.5 values
3053 This never produces a meaningless negative result. Such results
3054 cause Thinning/Thicken to not work correctly when used against a
3059 for (v=0; v < (ssize_t) kernel->height; v++)
3061 for (u=0; u < (ssize_t) kernel->width; u++)
3063 if (IfNaN(*k) == MagickFalse)
3067 if ((double) pixels[i] < pixel)
3068 pixel=(double) pixels[i];
3073 if ((double) pixels[i] > maximum)
3074 maximum=(double) pixels[i];
3079 pixels+=GetPixelChannels(image);
3081 pixels+=(image->columns-1)*GetPixelChannels(image);
3086 if (method == ThinningMorphology)
3087 pixel=(double) p[center+i]-pixel;
3089 if (method == ThickenMorphology)
3090 pixel+=(double) p[center+i]+pixel;
3093 case ErodeIntensityMorphology:
3096 Select pixel with minimum intensity within kernel neighbourhood.
3098 The kernel is not reflected for this operation.
3102 for (v=0; v < (ssize_t) kernel->height; v++)
3104 for (u=0; u < (ssize_t) kernel->width; u++)
3106 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3108 if (GetPixelIntensity(image,pixels) < minimum)
3110 pixel=(double) pixels[i];
3111 minimum=GetPixelIntensity(image,pixels);
3116 pixels+=GetPixelChannels(image);
3118 pixels+=(image->columns-1)*GetPixelChannels(image);
3122 case DilateIntensityMorphology:
3125 Select pixel with maximum intensity within kernel neighbourhood.
3127 The kernel is not reflected for this operation.
3130 k=(&kernel->values[kernel->width*kernel->height-1]);
3131 for (v=0; v < (ssize_t) kernel->height; v++)
3133 for (u=0; u < (ssize_t) kernel->width; u++)
3135 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3137 if (GetPixelIntensity(image,pixels) > maximum)
3139 pixel=(double) pixels[i];
3140 maximum=GetPixelIntensity(image,pixels);
3145 pixels+=GetPixelChannels(image);
3147 pixels+=(image->columns-1)*GetPixelChannels(image);
3151 case IterativeDistanceMorphology:
3154 Compute th iterative distance from black edge of a white image
3155 shape. Essentually white values are decreased to the smallest
3156 'distance from edge' it can find.
3158 It works by adding kernel values to the neighbourhood, and and
3159 select the minimum value found. The kernel is rotated before
3160 use, so kernel distances match resulting distances, when a user
3161 provided asymmetric kernel is applied.
3163 This code is nearly identical to True GrayScale Morphology but
3166 GreyDilate Kernel values added, maximum value found Kernel is
3169 GrayErode: Kernel values subtracted and minimum value found No
3170 kernel rotation used.
3172 Note the the Iterative Distance method is essentially a
3173 GrayErode, but with negative kernel values, and kernel rotation
3177 k=(&kernel->values[kernel->width*kernel->height-1]);
3178 for (v=0; v < (ssize_t) kernel->height; v++)
3180 for (u=0; u < (ssize_t) kernel->width; u++)
3182 if (IfNaN(*k) == MagickFalse)
3184 if ((pixels[i]+(*k)) < pixel)
3185 pixel=(double) pixels[i]+(*k);
3189 pixels+=GetPixelChannels(image);
3191 pixels+=(image->columns-1)*GetPixelChannels(image);
3195 case UndefinedMorphology:
3199 if (fabs(pixel-p[center+i]) > MagickEpsilon)
3201 gamma=PerceptibleReciprocal(gamma);
3202 gamma*=(double) kernel->height*kernel->width/count;
3203 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*pixel),q);
3205 p+=GetPixelChannels(image);
3206 q+=GetPixelChannels(morphology_image);
3208 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3210 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3215 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3216 #pragma omp critical (MagickCore_MorphologyPrimitive)
3218 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3219 if (proceed == MagickFalse)
3223 morphology_view=DestroyCacheView(morphology_view);
3224 image_view=DestroyCacheView(image_view);
3225 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
3226 changed+=changes[i];
3227 changes=(size_t *) RelinquishMagickMemory(changes);
3228 return(status ? (ssize_t) changed : -1);
3232 This is almost identical to the MorphologyPrimative() function above, but
3233 applies the primitive directly to the actual image using two passes, once in
3234 each direction, with the results of the previous (and current) row being
3237 That is after each row is 'Sync'ed' into the image, the next row makes use of
3238 those values as part of the calculation of the next row. It repeats, but
3239 going in the oppisite (bottom-up) direction.
3241 Because of this 're-use of results' this function can not make use of multi-
3242 threaded, parellel processing.
3244 static ssize_t MorphologyPrimitiveDirect(Image *image,
3245 const MorphologyMethod method,const KernelInfo *kernel,
3246 ExceptionInfo *exception)
3268 assert(image != (Image *) NULL);
3269 assert(image->signature == MagickSignature);
3270 assert(kernel != (KernelInfo *) NULL);
3271 assert(kernel->signature == MagickSignature);
3272 assert(exception != (ExceptionInfo *) NULL);
3273 assert(exception->signature == MagickSignature);
3279 case DistanceMorphology:
3280 case VoronoiMorphology:
3283 Kernel reflected about origin.
3285 offset.x=(ssize_t) kernel->width-kernel->x-1;
3286 offset.y=(ssize_t) kernel->height-kernel->y-1;
3297 Two views into same image, do not thread.
3299 image_view=AcquireVirtualCacheView(image,exception);
3300 morphology_view=AcquireAuthenticCacheView(image,exception);
3301 width=image->columns+kernel->width-1;
3302 for (y=0; y < (ssize_t) image->rows; y++)
3304 register const Quantum
3317 Read virtual pixels, and authentic pixels, from the same image! We read
3318 using virtual to get virtual pixel handling, but write back into the same
3321 Only top half of kernel is processed as we do a single pass downward
3322 through the image iterating the distance function as we go.
3324 if (status == MagickFalse)
3326 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t)
3327 offset.y+1,exception);
3328 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3330 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3335 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
3336 GetPixelChannels(image)*offset.x);
3337 for (x=0; x < (ssize_t) image->columns; x++)
3342 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3350 register const MagickRealType
3353 register const Quantum
3362 traits=GetPixelChannelTraits(image,(PixelChannel) i);
3363 if (traits == UndefinedPixelTrait)
3365 if (((traits & CopyPixelTrait) != 0) ||
3366 (GetPixelReadMask(image,p+center) == 0))
3369 pixel=(double) QuantumRange;
3372 case DistanceMorphology:
3374 k=(&kernel->values[kernel->width*kernel->height-1]);
3375 for (v=0; v <= offset.y; v++)
3377 for (u=0; u < (ssize_t) kernel->width; u++)
3379 if (IfNaN(*k) == MagickFalse)
3381 if ((pixels[i]+(*k)) < pixel)
3382 pixel=(double) pixels[i]+(*k);
3385 pixels+=GetPixelChannels(image);
3387 pixels+=(image->columns-1)*GetPixelChannels(image);
3389 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3390 pixels=q-offset.x*GetPixelChannels(image);
3391 for (u=0; u < offset.x; u++)
3393 if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3395 if ((pixels[i]+(*k)) < pixel)
3396 pixel=(double) pixels[i]+(*k);
3399 pixels+=GetPixelChannels(image);
3403 case VoronoiMorphology:
3405 k=(&kernel->values[kernel->width*kernel->height-1]);
3406 for (v=0; v < offset.y; v++)
3408 for (u=0; u < (ssize_t) kernel->width; u++)
3410 if (IfNaN(*k) == MagickFalse)
3412 if ((pixels[i]+(*k)) < pixel)
3413 pixel=(double) pixels[i]+(*k);
3416 pixels+=GetPixelChannels(image);
3418 pixels+=(image->columns-1)*GetPixelChannels(image);
3420 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3421 pixels=q-offset.x*GetPixelChannels(image);
3422 for (u=0; u < offset.x; u++)
3424 if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3426 if ((pixels[i]+(*k)) < pixel)
3427 pixel=(double) pixels[i]+(*k);
3430 pixels+=GetPixelChannels(image);
3437 if (fabs(pixel-q[i]) > MagickEpsilon)
3439 q[i]=ClampToQuantum(pixel);
3441 p+=GetPixelChannels(image);
3442 q+=GetPixelChannels(image);
3444 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3446 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3451 proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3452 if (proceed == MagickFalse)
3456 morphology_view=DestroyCacheView(morphology_view);
3457 image_view=DestroyCacheView(image_view);
3459 Do the reverse pass through the image.
3461 image_view=AcquireVirtualCacheView(image,exception);
3462 morphology_view=AcquireAuthenticCacheView(image,exception);
3463 for (y=(ssize_t) image->rows-1; y >= 0; y--)
3465 register const Quantum
3478 Read virtual pixels, and authentic pixels, from the same image. We
3479 read using virtual to get virtual pixel handling, but write back
3480 into the same image.
3482 Only the bottom half of the kernel is processed as we up the image.
3484 if (status == MagickFalse)
3486 p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t)
3487 kernel->y+1,exception);
3488 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3490 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3495 p+=(image->columns-1)*GetPixelChannels(image);
3496 q+=(image->columns-1)*GetPixelChannels(image);
3497 center=(ssize_t) (offset.x*GetPixelChannels(image));
3498 for (x=(ssize_t) image->columns-1; x >= 0; x--)
3503 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3511 register const MagickRealType
3514 register const Quantum
3523 traits=GetPixelChannelTraits(image,(PixelChannel) i);
3524 if (traits == UndefinedPixelTrait)
3526 if (((traits & CopyPixelTrait) != 0) ||
3527 (GetPixelReadMask(image,p+center) == 0))
3530 pixel=(double) QuantumRange;
3533 case DistanceMorphology:
3535 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3536 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3538 for (u=0; u < (ssize_t) kernel->width; u++)
3540 if (IfNaN(*k) == MagickFalse)
3542 if ((pixels[i]+(*k)) < pixel)
3543 pixel=(double) pixels[i]+(*k);
3546 pixels+=GetPixelChannels(image);
3548 pixels+=(image->columns-1)*GetPixelChannels(image);
3550 k=(&kernel->values[kernel->width*kernel->y+kernel->x-1]);
3551 pixels=q-offset.x*GetPixelChannels(image);
3552 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3554 if ((IfNaN(*k) == MagickFalse) &&
3555 ((x+u-offset.x) < (ssize_t) image->columns))
3557 if ((pixels[i]+(*k)) < pixel)
3558 pixel=(double) pixels[i]+(*k);
3561 pixels+=GetPixelChannels(image);
3565 case VoronoiMorphology:
3567 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3568 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3570 for (u=0; u < (ssize_t) kernel->width; u++)
3572 if (IfNaN(*k) == MagickFalse)
3574 if ((pixels[i]+(*k)) < pixel)
3575 pixel=(double) pixels[i]+(*k);
3578 pixels+=GetPixelChannels(image);
3580 pixels+=(image->columns-1)*GetPixelChannels(image);
3582 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3583 pixels=q-offset.x*GetPixelChannels(image);
3584 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3586 if ((IfNaN(*k) == MagickFalse) &&
3587 ((x+u-offset.x) < (ssize_t) image->columns))
3589 if ((pixels[i]+(*k)) < pixel)
3590 pixel=(double) pixels[i]+(*k);
3593 pixels+=GetPixelChannels(image);
3600 if (fabs(pixel-q[i]) > MagickEpsilon)
3602 q[i]=ClampToQuantum(pixel);
3604 p-=GetPixelChannels(image);
3605 q-=GetPixelChannels(image);
3607 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3609 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3614 proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3615 if (proceed == MagickFalse)
3619 morphology_view=DestroyCacheView(morphology_view);
3620 image_view=DestroyCacheView(image_view);
3621 return(status ? (ssize_t) changed : -1);
3625 Apply a Morphology by calling one of the above low level primitive
3626 application functions. This function handles any iteration loops,
3627 composition or re-iteration of results, and compound morphology methods that
3628 is based on multiple low-level (staged) morphology methods.
3630 Basically this provides the complex glue between the requested morphology
3631 method and raw low-level implementation (above).
3633 MagickPrivate Image *MorphologyApply(const Image *image,
3634 const MorphologyMethod method, const ssize_t iterations,
3635 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3636 ExceptionInfo *exception)
3642 *curr_image, /* Image we are working with or iterating */
3643 *work_image, /* secondary image for primitive iteration */
3644 *save_image, /* saved image - for 'edge' method only */
3645 *rslt_image; /* resultant image - after multi-kernel handling */
3648 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3649 *norm_kernel, /* the current normal un-reflected kernel */
3650 *rflt_kernel, /* the current reflected kernel (if needed) */
3651 *this_kernel; /* the kernel being applied */
3654 primitive; /* the current morphology primitive being applied */
3657 rslt_compose; /* multi-kernel compose method for results to use */
3660 special, /* do we use a direct modify function? */
3661 verbose; /* verbose output of results */
3664 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3665 method_limit, /* maximum number of compound method iterations */
3666 kernel_number, /* Loop 2: the kernel number being applied */
3667 stage_loop, /* Loop 3: primitive loop for compound morphology */
3668 stage_limit, /* how many primitives are in this compound */
3669 kernel_loop, /* Loop 4: iterate the kernel over image */
3670 kernel_limit, /* number of times to iterate kernel */
3671 count, /* total count of primitive steps applied */
3672 kernel_changed, /* total count of changed using iterated kernel */
3673 method_changed; /* total count of changed over method iteration */
3676 changed; /* number pixels changed by last primitive operation */
3681 assert(image != (Image *) NULL);
3682 assert(image->signature == MagickSignature);
3683 assert(kernel != (KernelInfo *) NULL);
3684 assert(kernel->signature == MagickSignature);
3685 assert(exception != (ExceptionInfo *) NULL);
3686 assert(exception->signature == MagickSignature);
3688 count = 0; /* number of low-level morphology primitives performed */
3689 if ( iterations == 0 )
3690 return((Image *)NULL); /* null operation - nothing to do! */
3692 kernel_limit = (size_t) iterations;
3693 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3694 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3696 verbose = IsStringTrue(GetImageArtifact(image,"debug"));
3698 /* initialise for cleanup */
3699 curr_image = (Image *) image;
3700 curr_compose = image->compose;
3701 (void) curr_compose;
3702 work_image = save_image = rslt_image = (Image *) NULL;
3703 reflected_kernel = (KernelInfo *) NULL;
3705 /* Initialize specific methods
3706 * + which loop should use the given iteratations
3707 * + how many primitives make up the compound morphology
3708 * + multi-kernel compose method to use (by default)
3710 method_limit = 1; /* just do method once, unless otherwise set */
3711 stage_limit = 1; /* assume method is not a compound */
3712 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3713 rslt_compose = compose; /* and we are composing multi-kernels as given */
3715 case SmoothMorphology: /* 4 primitive compound morphology */
3718 case OpenMorphology: /* 2 primitive compound morphology */
3719 case OpenIntensityMorphology:
3720 case TopHatMorphology:
3721 case CloseMorphology:
3722 case CloseIntensityMorphology:
3723 case BottomHatMorphology:
3724 case EdgeMorphology:
3727 case HitAndMissMorphology:
3728 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3730 case ThinningMorphology:
3731 case ThickenMorphology:
3732 method_limit = kernel_limit; /* iterate the whole method */
3733 kernel_limit = 1; /* do not do kernel iteration */
3735 case DistanceMorphology:
3736 case VoronoiMorphology:
3737 special = MagickTrue; /* use special direct primative */
3743 /* Apply special methods with special requirments
3744 ** For example, single run only, or post-processing requirements
3746 if ( special != MagickFalse )
3748 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3749 if (rslt_image == (Image *) NULL)
3751 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3754 changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception);
3756 if ( IfMagickTrue(verbose) )
3757 (void) (void) FormatLocaleFile(stderr,
3758 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3759 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3760 1.0,0.0,1.0, (double) changed);
3765 if ( method == VoronoiMorphology ) {
3766 /* Preserve the alpha channel of input image - but turned it off */
3767 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3769 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3770 MagickTrue,0,0,exception);
3771 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3777 /* Handle user (caller) specified multi-kernel composition method */
3778 if ( compose != UndefinedCompositeOp )
3779 rslt_compose = compose; /* override default composition for method */
3780 if ( rslt_compose == UndefinedCompositeOp )
3781 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3783 /* Some methods require a reflected kernel to use with primitives.
3784 * Create the reflected kernel for those methods. */
3786 case CorrelateMorphology:
3787 case CloseMorphology:
3788 case CloseIntensityMorphology:
3789 case BottomHatMorphology:
3790 case SmoothMorphology:
3791 reflected_kernel = CloneKernelInfo(kernel);
3792 if (reflected_kernel == (KernelInfo *) NULL)
3794 RotateKernelInfo(reflected_kernel,180);
3800 /* Loops around more primitive morpholgy methods
3801 ** erose, dilate, open, close, smooth, edge, etc...
3803 /* Loop 1: iterate the compound method */
3806 while ( method_loop < method_limit && method_changed > 0 ) {
3810 /* Loop 2: iterate over each kernel in a multi-kernel list */
3811 norm_kernel = (KernelInfo *) kernel;
3812 this_kernel = (KernelInfo *) kernel;
3813 rflt_kernel = reflected_kernel;
3816 while ( norm_kernel != NULL ) {
3818 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3819 stage_loop = 0; /* the compound morphology stage number */
3820 while ( stage_loop < stage_limit ) {
3821 stage_loop++; /* The stage of the compound morphology */
3823 /* Select primitive morphology for this stage of compound method */
3824 this_kernel = norm_kernel; /* default use unreflected kernel */
3825 primitive = method; /* Assume method is a primitive */
3827 case ErodeMorphology: /* just erode */
3828 case EdgeInMorphology: /* erode and image difference */
3829 primitive = ErodeMorphology;
3831 case DilateMorphology: /* just dilate */
3832 case EdgeOutMorphology: /* dilate and image difference */
3833 primitive = DilateMorphology;
3835 case OpenMorphology: /* erode then dialate */
3836 case TopHatMorphology: /* open and image difference */
3837 primitive = ErodeMorphology;
3838 if ( stage_loop == 2 )
3839 primitive = DilateMorphology;
3841 case OpenIntensityMorphology:
3842 primitive = ErodeIntensityMorphology;
3843 if ( stage_loop == 2 )
3844 primitive = DilateIntensityMorphology;
3846 case CloseMorphology: /* dilate, then erode */
3847 case BottomHatMorphology: /* close and image difference */
3848 this_kernel = rflt_kernel; /* use the reflected kernel */
3849 primitive = DilateMorphology;
3850 if ( stage_loop == 2 )
3851 primitive = ErodeMorphology;
3853 case CloseIntensityMorphology:
3854 this_kernel = rflt_kernel; /* use the reflected kernel */
3855 primitive = DilateIntensityMorphology;
3856 if ( stage_loop == 2 )
3857 primitive = ErodeIntensityMorphology;
3859 case SmoothMorphology: /* open, close */
3860 switch ( stage_loop ) {
3861 case 1: /* start an open method, which starts with Erode */
3862 primitive = ErodeMorphology;
3864 case 2: /* now Dilate the Erode */
3865 primitive = DilateMorphology;
3867 case 3: /* Reflect kernel a close */
3868 this_kernel = rflt_kernel; /* use the reflected kernel */
3869 primitive = DilateMorphology;
3871 case 4: /* Finish the Close */
3872 this_kernel = rflt_kernel; /* use the reflected kernel */
3873 primitive = ErodeMorphology;
3877 case EdgeMorphology: /* dilate and erode difference */
3878 primitive = DilateMorphology;
3879 if ( stage_loop == 2 ) {
3880 save_image = curr_image; /* save the image difference */
3881 curr_image = (Image *) image;
3882 primitive = ErodeMorphology;
3885 case CorrelateMorphology:
3886 /* A Correlation is a Convolution with a reflected kernel.
3887 ** However a Convolution is a weighted sum using a reflected
3888 ** kernel. It may seem stange to convert a Correlation into a
3889 ** Convolution as the Correlation is the simplier method, but
3890 ** Convolution is much more commonly used, and it makes sense to
3891 ** implement it directly so as to avoid the need to duplicate the
3892 ** kernel when it is not required (which is typically the
3895 this_kernel = rflt_kernel; /* use the reflected kernel */
3896 primitive = ConvolveMorphology;
3901 assert( this_kernel != (KernelInfo *) NULL );
3903 /* Extra information for debugging compound operations */
3904 if ( IfMagickTrue(verbose) ) {
3905 if ( stage_limit > 1 )
3906 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
3907 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
3908 method_loop,(double) stage_loop);
3909 else if ( primitive != method )
3910 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
3911 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
3917 /* Loop 4: Iterate the kernel with primitive */
3921 while ( kernel_loop < kernel_limit && changed > 0 ) {
3922 kernel_loop++; /* the iteration of this kernel */
3924 /* Create a clone as the destination image, if not yet defined */
3925 if ( work_image == (Image *) NULL )
3927 work_image=CloneImage(image,0,0,MagickTrue,exception);
3928 if (work_image == (Image *) NULL)
3930 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
3934 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3936 changed = MorphologyPrimitive(curr_image, work_image, primitive,
3937 this_kernel, bias, exception);
3939 if ( IfMagickTrue(verbose) ) {
3940 if ( kernel_loop > 1 )
3941 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
3942 (void) (void) FormatLocaleFile(stderr,
3943 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3944 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
3945 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
3946 (double) (method_loop+kernel_loop-1),(double) kernel_number,
3947 (double) count,(double) changed);
3951 kernel_changed += changed;
3952 method_changed += changed;
3954 /* prepare next loop */
3955 { Image *tmp = work_image; /* swap images for iteration */
3956 work_image = curr_image;
3959 if ( work_image == image )
3960 work_image = (Image *) NULL; /* replace input 'image' */
3962 } /* End Loop 4: Iterate the kernel with primitive */
3964 if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed )
3965 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
3966 if ( IfMagickTrue(verbose) && stage_loop < stage_limit )
3967 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
3970 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3971 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
3972 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
3973 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
3974 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
3977 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3979 /* Final Post-processing for some Compound Methods
3981 ** The removal of any 'Sync' channel flag in the Image Compositon
3982 ** below ensures the methematical compose method is applied in a
3983 ** purely mathematical way, and only to the selected channels.
3984 ** Turn off SVG composition 'alpha blending'.
3987 case EdgeOutMorphology:
3988 case EdgeInMorphology:
3989 case TopHatMorphology:
3990 case BottomHatMorphology:
3991 if ( IfMagickTrue(verbose) )
3992 (void) FormatLocaleFile(stderr,
3993 "\n%s: Difference with original image",CommandOptionToMnemonic(
3994 MagickMorphologyOptions, method) );
3995 (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
3996 MagickTrue,0,0,exception);
3998 case EdgeMorphology:
3999 if ( IfMagickTrue(verbose) )
4000 (void) FormatLocaleFile(stderr,
4001 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
4002 MagickMorphologyOptions, method) );
4003 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
4004 MagickTrue,0,0,exception);
4005 save_image = DestroyImage(save_image); /* finished with save image */
4011 /* multi-kernel handling: re-iterate, or compose results */
4012 if ( kernel->next == (KernelInfo *) NULL )
4013 rslt_image = curr_image; /* just return the resulting image */
4014 else if ( rslt_compose == NoCompositeOp )
4015 { if ( IfMagickTrue(verbose) ) {
4016 if ( this_kernel->next != (KernelInfo *) NULL )
4017 (void) FormatLocaleFile(stderr, " (re-iterate)");
4019 (void) FormatLocaleFile(stderr, " (done)");
4021 rslt_image = curr_image; /* return result, and re-iterate */
4023 else if ( rslt_image == (Image *) NULL)
4024 { if ( IfMagickTrue(verbose) )
4025 (void) FormatLocaleFile(stderr, " (save for compose)");
4026 rslt_image = curr_image;
4027 curr_image = (Image *) image; /* continue with original image */
4030 { /* Add the new 'current' result to the composition
4032 ** The removal of any 'Sync' channel flag in the Image Compositon
4033 ** below ensures the methematical compose method is applied in a
4034 ** purely mathematical way, and only to the selected channels.
4035 ** IE: Turn off SVG composition 'alpha blending'.
4037 if ( IfMagickTrue(verbose) )
4038 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4039 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4040 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4042 curr_image = DestroyImage(curr_image);
4043 curr_image = (Image *) image; /* continue with original image */
4045 if ( IfMagickTrue(verbose) )
4046 (void) FormatLocaleFile(stderr, "\n");
4048 /* loop to the next kernel in a multi-kernel list */
4049 norm_kernel = norm_kernel->next;
4050 if ( rflt_kernel != (KernelInfo *) NULL )
4051 rflt_kernel = rflt_kernel->next;
4053 } /* End Loop 2: Loop over each kernel */
4055 } /* End Loop 1: compound method interation */
4059 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4061 if ( curr_image == rslt_image )
4062 curr_image = (Image *) NULL;
4063 if ( rslt_image != (Image *) NULL )
4064 rslt_image = DestroyImage(rslt_image);
4066 if ( curr_image == rslt_image || curr_image == image )
4067 curr_image = (Image *) NULL;
4068 if ( curr_image != (Image *) NULL )
4069 curr_image = DestroyImage(curr_image);
4070 if ( work_image != (Image *) NULL )
4071 work_image = DestroyImage(work_image);
4072 if ( save_image != (Image *) NULL )
4073 save_image = DestroyImage(save_image);
4074 if ( reflected_kernel != (KernelInfo *) NULL )
4075 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4081 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4085 % M o r p h o l o g y I m a g e %
4089 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4091 % MorphologyImage() applies a user supplied kernel to the image according to
4092 % the given mophology method.
4094 % This function applies any and all user defined settings before calling
4095 % the above internal function MorphologyApply().
4097 % User defined settings include...
4098 % * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4099 % * Kernel Scale/normalize settings ("-define convolve:scale=??")
4100 % This can also includes the addition of a scaled unity kernel.
4101 % * Show Kernel being applied ("-define showkernel=1")
4103 % Other operators that do not want user supplied options interfering,
4104 % especially "convolve:bias" and "showkernel" should use MorphologyApply()
4107 % The format of the MorphologyImage method is:
4109 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4110 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4112 % A description of each parameter follows:
4114 % o image: the image.
4116 % o method: the morphology method to be applied.
4118 % o iterations: apply the operation this many times (or no change).
4119 % A value of -1 means loop until no change found.
4120 % How this is applied may depend on the morphology method.
4121 % Typically this is a value of 1.
4123 % o kernel: An array of double representing the morphology kernel.
4124 % Warning: kernel may be normalized for the Convolve method.
4126 % o exception: return any errors or warnings in this structure.
4129 MagickExport Image *MorphologyImage(const Image *image,
4130 const MorphologyMethod method,const ssize_t iterations,
4131 const KernelInfo *kernel,ExceptionInfo *exception)
4145 curr_kernel = (KernelInfo *) kernel;
4147 compose = UndefinedCompositeOp; /* use default for method */
4149 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4150 * This is done BEFORE the ShowKernelInfo() function is called so that
4151 * users can see the results of the 'option:convolve:scale' option.
4153 if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4157 /* Get the bias value as it will be needed */
4158 artifact = GetImageArtifact(image,"convolve:bias");
4159 if ( artifact != (const char *) NULL) {
4160 if (IfMagickFalse(IsGeometry(artifact)))
4161 (void) ThrowMagickException(exception,GetMagickModule(),
4162 OptionWarning,"InvalidSetting","'%s' '%s'",
4163 "convolve:bias",artifact);
4165 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4168 /* Scale kernel according to user wishes */
4169 artifact = GetImageArtifact(image,"convolve:scale");
4170 if ( artifact != (const char *)NULL ) {
4171 if (IfMagickFalse(IsGeometry(artifact)))
4172 (void) ThrowMagickException(exception,GetMagickModule(),
4173 OptionWarning,"InvalidSetting","'%s' '%s'",
4174 "convolve:scale",artifact);
4176 if ( curr_kernel == kernel )
4177 curr_kernel = CloneKernelInfo(kernel);
4178 if (curr_kernel == (KernelInfo *) NULL)
4179 return((Image *) NULL);
4180 ScaleGeometryKernelInfo(curr_kernel, artifact);
4185 /* display the (normalized) kernel via stderr */
4186 if ( IfStringTrue(GetImageArtifact(image,"showkernel"))
4187 || IfStringTrue(GetImageArtifact(image,"convolve:showkernel"))
4188 || IfStringTrue(GetImageArtifact(image,"morphology:showkernel")) )
4189 ShowKernelInfo(curr_kernel);
4191 /* Override the default handling of multi-kernel morphology results
4192 * If 'Undefined' use the default method
4193 * If 'None' (default for 'Convolve') re-iterate previous result
4194 * Otherwise merge resulting images using compose method given.
4195 * Default for 'HitAndMiss' is 'Lighten'.
4202 artifact = GetImageArtifact(image,"morphology:compose");
4203 if ( artifact != (const char *) NULL) {
4204 parse=ParseCommandOption(MagickComposeOptions,
4205 MagickFalse,artifact);
4207 (void) ThrowMagickException(exception,GetMagickModule(),
4208 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4209 "morphology:compose",artifact);
4211 compose=(CompositeOperator)parse;
4214 /* Apply the Morphology */
4215 morphology_image = MorphologyApply(image,method,iterations,
4216 curr_kernel,compose,bias,exception);
4218 /* Cleanup and Exit */
4219 if ( curr_kernel != kernel )
4220 curr_kernel=DestroyKernelInfo(curr_kernel);
4221 return(morphology_image);
4225 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4229 + R o t a t e K e r n e l I n f o %
4233 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4235 % RotateKernelInfo() rotates the kernel by the angle given.
4237 % Currently it is restricted to 90 degree angles, of either 1D kernels
4238 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4239 % It will ignore usless rotations for specific 'named' built-in kernels.
4241 % The format of the RotateKernelInfo method is:
4243 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4245 % A description of each parameter follows:
4247 % o kernel: the Morphology/Convolution kernel
4249 % o angle: angle to rotate in degrees
4251 % This function is currently internal to this module only, but can be exported
4252 % to other modules if needed.
4254 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4256 /* angle the lower kernels first */
4257 if ( kernel->next != (KernelInfo *) NULL)
4258 RotateKernelInfo(kernel->next, angle);
4260 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4262 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4265 /* Modulus the angle */
4266 angle = fmod(angle, 360.0);
4270 if ( 337.5 < angle || angle <= 22.5 )
4271 return; /* Near zero angle - no change! - At least not at this time */
4273 /* Handle special cases */
4274 switch (kernel->type) {
4275 /* These built-in kernels are cylindrical kernels, rotating is useless */
4276 case GaussianKernel:
4281 case LaplacianKernel:
4282 case ChebyshevKernel:
4283 case ManhattanKernel:
4284 case EuclideanKernel:
4287 /* These may be rotatable at non-90 angles in the future */
4288 /* but simply rotating them in multiples of 90 degrees is useless */
4295 /* These only allows a +/-90 degree rotation (by transpose) */
4296 /* A 180 degree rotation is useless */
4298 if ( 135.0 < angle && angle <= 225.0 )
4300 if ( 225.0 < angle && angle <= 315.0 )
4307 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4308 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4310 if ( kernel->width == 3 && kernel->height == 3 )
4311 { /* Rotate a 3x3 square by 45 degree angle */
4312 double t = kernel->values[0];
4313 kernel->values[0] = kernel->values[3];
4314 kernel->values[3] = kernel->values[6];
4315 kernel->values[6] = kernel->values[7];
4316 kernel->values[7] = kernel->values[8];
4317 kernel->values[8] = kernel->values[5];
4318 kernel->values[5] = kernel->values[2];
4319 kernel->values[2] = kernel->values[1];
4320 kernel->values[1] = t;
4321 /* rotate non-centered origin */
4322 if ( kernel->x != 1 || kernel->y != 1 ) {
4324 x = (ssize_t) kernel->x-1;
4325 y = (ssize_t) kernel->y-1;
4326 if ( x == y ) x = 0;
4327 else if ( x == 0 ) x = -y;
4328 else if ( x == -y ) y = 0;
4329 else if ( y == 0 ) y = x;
4330 kernel->x = (ssize_t) x+1;
4331 kernel->y = (ssize_t) y+1;
4333 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4334 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4337 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4339 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4341 if ( kernel->width == 1 || kernel->height == 1 )
4342 { /* Do a transpose of a 1 dimensional kernel,
4343 ** which results in a fast 90 degree rotation of some type.
4347 t = (ssize_t) kernel->width;
4348 kernel->width = kernel->height;
4349 kernel->height = (size_t) t;
4351 kernel->x = kernel->y;
4353 if ( kernel->width == 1 ) {
4354 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4355 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4357 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4358 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4361 else if ( kernel->width == kernel->height )
4362 { /* Rotate a square array of values by 90 degrees */
4366 register MagickRealType
4370 for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--)
4371 for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--)
4372 { t = k[i+j*kernel->width];
4373 k[i+j*kernel->width] = k[j+x*kernel->width];
4374 k[j+x*kernel->width] = k[x+y*kernel->width];
4375 k[x+y*kernel->width] = k[y+i*kernel->width];
4376 k[y+i*kernel->width] = t;
4379 /* rotate the origin - relative to center of array */
4380 { register ssize_t x,y;
4381 x = (ssize_t) (kernel->x*2-kernel->width+1);
4382 y = (ssize_t) (kernel->y*2-kernel->height+1);
4383 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4384 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4386 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4387 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4390 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4392 if ( 135.0 < angle && angle <= 225.0 )
4394 /* For a 180 degree rotation - also know as a reflection
4395 * This is actually a very very common operation!
4396 * Basically all that is needed is a reversal of the kernel data!
4397 * And a reflection of the origon
4402 register MagickRealType
4410 j=(ssize_t) (kernel->width*kernel->height-1);
4411 for (i=0; i < j; i++, j--)
4412 t=k[i], k[i]=k[j], k[j]=t;
4414 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4415 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4416 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4417 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4419 /* At this point angle should at least between -45 (315) and +45 degrees
4420 * In the future some form of non-orthogonal angled rotates could be
4421 * performed here, posibily with a linear kernel restriction.
4428 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4432 % S c a l e G e o m e t r y K e r n e l I n f o %
4436 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4438 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4439 % provided as a "-set option:convolve:scale {geometry}" user setting,
4440 % and modifies the kernel according to the parsed arguments of that setting.
4442 % The first argument (and any normalization flags) are passed to
4443 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4444 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4445 % into the scaled/normalized kernel.
4447 % The format of the ScaleGeometryKernelInfo method is:
4449 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4450 % const double scaling_factor,const MagickStatusType normalize_flags)
4452 % A description of each parameter follows:
4454 % o kernel: the Morphology/Convolution kernel to modify
4457 % The geometry string to parse, typically from the user provided
4458 % "-set option:convolve:scale {geometry}" setting.
4461 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4462 const char *geometry)
4470 SetGeometryInfo(&args);
4471 flags = ParseGeometry(geometry, &args);
4474 /* For Debugging Geometry Input */
4475 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4476 flags, args.rho, args.sigma, args.xi, args.psi );
4479 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4480 args.rho *= 0.01, args.sigma *= 0.01;
4482 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4484 if ( (flags & SigmaValue) == 0 )
4487 /* Scale/Normalize the input kernel */
4488 ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4490 /* Add Unity Kernel, for blending with original */
4491 if ( (flags & SigmaValue) != 0 )
4492 UnityAddKernelInfo(kernel, args.sigma);
4497 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4501 % S c a l e K e r n e l I n f o %
4505 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4507 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4508 % without normalization of the sum of the kernel values (as per given flags).
4510 % By default (no flags given) the values within the kernel is scaled
4511 % directly using given scaling factor without change.
4513 % If either of the two 'normalize_flags' are given the kernel will first be
4514 % normalized and then further scaled by the scaling factor value given.
4516 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4517 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4518 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4519 % non-HDRI versions of IM this may cause images to have any negative results
4520 % clipped, unless some 'bias' is used.
4522 % More specifically. Kernels which only contain positive values (such as a
4523 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4524 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4526 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4527 % the kernel will be scaled by the absolute of the sum of kernel values, so
4528 % that it will generally fall within the +/- 1.0 range.
4530 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4531 % will be scaled by just the sum of the postive values, so that its output
4532 % range will again fall into the +/- 1.0 range.
4534 % For special kernels designed for locating shapes using 'Correlate', (often
4535 % only containing +1 and -1 values, representing foreground/brackground
4536 % matching) a special normalization method is provided to scale the positive
4537 % values separately to those of the negative values, so the kernel will be
4538 % forced to become a zero-sum kernel better suited to such searches.
4540 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4541 % attributes within the kernel structure have been correctly set during the
4544 % NOTE: The values used for 'normalize_flags' have been selected specifically
4545 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4546 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4548 % The format of the ScaleKernelInfo method is:
4550 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4551 % const MagickStatusType normalize_flags )
4553 % A description of each parameter follows:
4555 % o kernel: the Morphology/Convolution kernel
4558 % multiply all values (after normalization) by this factor if not
4559 % zero. If the kernel is normalized regardless of any flags.
4561 % o normalize_flags:
4562 % GeometryFlags defining normalization method to use.
4563 % specifically: NormalizeValue, CorrelateNormalizeValue,
4564 % and/or PercentValue
4567 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4568 const double scaling_factor,const GeometryFlags normalize_flags)
4577 /* do the other kernels in a multi-kernel list first */
4578 if ( kernel->next != (KernelInfo *) NULL)
4579 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4581 /* Normalization of Kernel */
4583 if ( (normalize_flags&NormalizeValue) != 0 ) {
4584 if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4585 /* non-zero-summing kernel (generally positive) */
4586 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4588 /* zero-summing kernel */
4589 pos_scale = kernel->positive_range;
4591 /* Force kernel into a normalized zero-summing kernel */
4592 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4593 pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4594 ? kernel->positive_range : 1.0;
4595 neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4596 ? -kernel->negative_range : 1.0;
4599 neg_scale = pos_scale;
4601 /* finialize scaling_factor for positive and negative components */
4602 pos_scale = scaling_factor/pos_scale;
4603 neg_scale = scaling_factor/neg_scale;
4605 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4606 if ( ! IfNaN(kernel->values[i]) )
4607 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4609 /* convolution output range */
4610 kernel->positive_range *= pos_scale;
4611 kernel->negative_range *= neg_scale;
4612 /* maximum and minimum values in kernel */
4613 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4614 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4616 /* swap kernel settings if user's scaling factor is negative */
4617 if ( scaling_factor < MagickEpsilon ) {
4619 t = kernel->positive_range;
4620 kernel->positive_range = kernel->negative_range;
4621 kernel->negative_range = t;
4622 t = kernel->maximum;
4623 kernel->maximum = kernel->minimum;
4624 kernel->minimum = 1;
4631 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4635 % S h o w K e r n e l I n f o %
4639 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4641 % ShowKernelInfo() outputs the details of the given kernel defination to
4642 % standard error, generally due to a users 'showkernel' option request.
4644 % The format of the ShowKernel method is:
4646 % void ShowKernelInfo(const KernelInfo *kernel)
4648 % A description of each parameter follows:
4650 % o kernel: the Morphology/Convolution kernel
4653 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4661 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4663 (void) FormatLocaleFile(stderr, "Kernel");
4664 if ( kernel->next != (KernelInfo *) NULL )
4665 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4666 (void) FormatLocaleFile(stderr, " \"%s",
4667 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4668 if ( fabs(k->angle) >= MagickEpsilon )
4669 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4670 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4671 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4672 (void) FormatLocaleFile(stderr,
4673 " with values from %.*lg to %.*lg\n",
4674 GetMagickPrecision(), k->minimum,
4675 GetMagickPrecision(), k->maximum);
4676 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4677 GetMagickPrecision(), k->negative_range,
4678 GetMagickPrecision(), k->positive_range);
4679 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4680 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4681 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4682 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4684 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4685 GetMagickPrecision(), k->positive_range+k->negative_range);
4686 for (i=v=0; v < k->height; v++) {
4687 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4688 for (u=0; u < k->width; u++, i++)
4689 if ( IfNaN(k->values[i]) )
4690 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4692 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4693 GetMagickPrecision(), (double) k->values[i]);
4694 (void) FormatLocaleFile(stderr,"\n");
4700 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4704 % U n i t y A d d K e r n a l I n f o %
4708 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4710 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4711 % to the given pre-scaled and normalized Kernel. This in effect adds that
4712 % amount of the original image into the resulting convolution kernel. This
4713 % value is usually provided by the user as a percentage value in the
4714 % 'convolve:scale' setting.
4716 % The resulting effect is to convert the defined kernels into blended
4717 % soft-blurs, unsharp kernels or into sharpening kernels.
4719 % The format of the UnityAdditionKernelInfo method is:
4721 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4723 % A description of each parameter follows:
4725 % o kernel: the Morphology/Convolution kernel
4728 % scaling factor for the unity kernel to be added to
4732 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4735 /* do the other kernels in a multi-kernel list first */
4736 if ( kernel->next != (KernelInfo *) NULL)
4737 UnityAddKernelInfo(kernel->next, scale);
4739 /* Add the scaled unity kernel to the existing kernel */
4740 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4741 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4747 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4751 % Z e r o K e r n e l N a n s %
4755 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4757 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4758 % the kernel with a zero value. This is typically done when the kernel will
4759 % be used in special hardware (GPU) convolution processors, to simply
4762 % The format of the ZeroKernelNans method is:
4764 % void ZeroKernelNans (KernelInfo *kernel)
4766 % A description of each parameter follows:
4768 % o kernel: the Morphology/Convolution kernel
4771 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4776 /* do the other kernels in a multi-kernel list first */
4777 if ( kernel->next != (KernelInfo *) NULL)
4778 ZeroKernelNans(kernel->next);
4780 for (i=0; i < (kernel->width*kernel->height); i++)
4781 if ( IfNaN(kernel->values[i]) )
4782 kernel->values[i] = 0.0;