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-2015 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 #define Minimize(assign,value) assign=MagickMin(assign,value)
93 #define Maximize(assign,value) assign=MagickMax(assign,value)
95 /* Integer Factorial Function - for a Binomial kernel */
97 static inline size_t fact(size_t n)
100 for(f=1, l=2; l <= n; f=f*l, l++);
103 #elif 1 /* glibc floating point alternatives */
104 #define fact(n) ((size_t)tgamma((double)n+1))
106 #define fact(n) ((size_t)lgamma((double)n+1))
110 /* Currently these are only internal to this module */
112 CalcKernelMetaData(KernelInfo *),
113 ExpandMirrorKernelInfo(KernelInfo *),
114 ExpandRotateKernelInfo(KernelInfo *, const double),
115 RotateKernelInfo(KernelInfo *, double);
118 /* Quick function to find last kernel in a kernel list */
119 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
121 while (kernel->next != (KernelInfo *) NULL)
127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
131 % A c q u i r e K e r n e l I n f o %
135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
137 % AcquireKernelInfo() takes the given string (generally supplied by the
138 % user) and converts it into a Morphology/Convolution Kernel. This allows
139 % users to specify a kernel from a number of pre-defined kernels, or to fully
140 % specify their own kernel for a specific Convolution or Morphology
143 % The kernel so generated can be any rectangular array of floating point
144 % values (doubles) with the 'control point' or 'pixel being affected'
145 % anywhere within that array of values.
147 % Previously IM was restricted to a square of odd size using the exact
148 % center as origin, this is no longer the case, and any rectangular kernel
149 % with any value being declared the origin. This in turn allows the use of
150 % highly asymmetrical kernels.
152 % The floating point values in the kernel can also include a special value
153 % known as 'nan' or 'not a number' to indicate that this value is not part
154 % of the kernel array. This allows you to shaped the kernel within its
155 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
156 % shape. However at least one non-nan value must be provided for correct
157 % working of a kernel.
159 % The returned kernel should be freed using the DestroyKernelInfo() when you
160 % are finished with it. Do not free this memory yourself.
162 % Input kernel defintion strings can consist of any of three types.
165 % Select from one of the built in kernels, using the name and
166 % geometry arguments supplied. See AcquireKernelBuiltIn()
168 % "WxH[+X+Y][@><]:num, num, num ..."
169 % a kernel of size W by H, with W*H floating point numbers following.
170 % the 'center' can be optionally be defined at +X+Y (such that +0+0
171 % is top left corner). If not defined the pixel in the center, for
172 % odd sizes, or to the immediate top or left of center for even sizes
173 % is automatically selected.
175 % "num, num, num, num, ..."
176 % list of floating point numbers defining an 'old style' odd sized
177 % square kernel. At least 9 values should be provided for a 3x3
178 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
179 % Values can be space or comma separated. This is not recommended.
181 % You can define a 'list of kernels' which can be used by some morphology
182 % operators A list is defined as a semi-colon separated list kernels.
184 % " kernel ; kernel ; kernel ; "
186 % Any extra ';' characters, at start, end or between kernel defintions are
189 % The special flags will expand a single kernel, into a list of rotated
190 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
191 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
192 % The '<' also exands using 90-degree rotates, but giving a 180-degree
193 % reflected kernel before the +/- 90-degree rotations, which can be important
194 % for Thinning operations.
196 % Note that 'name' kernels will start with an alphabetic character while the
197 % new kernel specification has a ':' character in its specification string.
198 % If neither is the case, it is assumed an old style of a simple list of
199 % numbers generating a odd-sized square kernel has been given.
201 % The format of the AcquireKernal method is:
203 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
205 % A description of each parameter follows:
207 % o kernel_string: the Morphology/Convolution kernel wanted.
211 /* This was separated so that it could be used as a separate
212 ** array input handling function, such as for -color-matrix
214 static KernelInfo *ParseKernelArray(const char *kernel_string)
220 token[MaxTextExtent];
230 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
238 kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
239 if (kernel == (KernelInfo *) NULL)
241 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
242 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
243 kernel->negative_range = kernel->positive_range = 0.0;
244 kernel->type = UserDefinedKernel;
245 kernel->next = (KernelInfo *) NULL;
246 kernel->signature = MagickSignature;
247 if (kernel_string == (const char *) NULL)
250 /* find end of this specific kernel definition string */
251 end = strchr(kernel_string, ';');
252 if ( end == (char *) NULL )
253 end = strchr(kernel_string, '\0');
255 /* clear flags - for Expanding kernel lists thorugh rotations */
258 /* Has a ':' in argument - New user kernel specification
259 FUTURE: this split on ':' could be done by StringToken()
261 p = strchr(kernel_string, ':');
262 if ( p != (char *) NULL && p < end)
264 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
265 memcpy(token, kernel_string, (size_t) (p-kernel_string));
266 token[p-kernel_string] = '\0';
267 SetGeometryInfo(&args);
268 flags = ParseGeometry(token, &args);
270 /* Size handling and checks of geometry settings */
271 if ( (flags & WidthValue) == 0 ) /* if no width then */
272 args.rho = args.sigma; /* then width = height */
273 if ( args.rho < 1.0 ) /* if width too small */
274 args.rho = 1.0; /* then width = 1 */
275 if ( args.sigma < 1.0 ) /* if height too small */
276 args.sigma = args.rho; /* then height = width */
277 kernel->width = (size_t)args.rho;
278 kernel->height = (size_t)args.sigma;
280 /* Offset Handling and Checks */
281 if ( args.xi < 0.0 || args.psi < 0.0 )
282 return(DestroyKernelInfo(kernel));
283 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
284 : (ssize_t) (kernel->width-1)/2;
285 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
286 : (ssize_t) (kernel->height-1)/2;
287 if ( kernel->x >= (ssize_t) kernel->width ||
288 kernel->y >= (ssize_t) kernel->height )
289 return(DestroyKernelInfo(kernel));
291 p++; /* advance beyond the ':' */
294 { /* ELSE - Old old specification, forming odd-square kernel */
295 /* count up number of values given */
296 p=(const char *) kernel_string;
297 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
298 p++; /* ignore "'" chars for convolve filter usage - Cristy */
299 for (i=0; p < end; i++)
301 GetMagickToken(p,&p,token);
303 GetMagickToken(p,&p,token);
305 /* set the size of the kernel - old sized square */
306 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
307 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
308 p=(const char *) kernel_string;
309 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
310 p++; /* ignore "'" chars for convolve filter usage - Cristy */
313 /* Read in the kernel values from rest of input string argument */
314 kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
315 kernel->width,kernel->height*sizeof(*kernel->values)));
316 if (kernel->values == (MagickRealType *) NULL)
317 return(DestroyKernelInfo(kernel));
318 kernel->minimum=MagickMaximumValue;
319 kernel->maximum=(-MagickMaximumValue);
320 kernel->negative_range = kernel->positive_range = 0.0;
321 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
323 GetMagickToken(p,&p,token);
325 GetMagickToken(p,&p,token);
326 if ( LocaleCompare("nan",token) == 0
327 || LocaleCompare("-",token) == 0 ) {
328 kernel->values[i] = nan; /* this value is not part of neighbourhood */
331 kernel->values[i] = StringToDouble(token,(char **) NULL);
332 ( kernel->values[i] < 0)
333 ? ( kernel->negative_range += kernel->values[i] )
334 : ( kernel->positive_range += kernel->values[i] );
335 Minimize(kernel->minimum, kernel->values[i]);
336 Maximize(kernel->maximum, kernel->values[i]);
340 /* sanity check -- no more values in kernel definition */
341 GetMagickToken(p,&p,token);
342 if ( *token != '\0' && *token != ';' && *token != '\'' )
343 return(DestroyKernelInfo(kernel));
346 /* this was the old method of handling a incomplete kernel */
347 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
348 Minimize(kernel->minimum, kernel->values[i]);
349 Maximize(kernel->maximum, kernel->values[i]);
350 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
351 kernel->values[i]=0.0;
354 /* Number of values for kernel was not enough - Report Error */
355 if ( i < (ssize_t) (kernel->width*kernel->height) )
356 return(DestroyKernelInfo(kernel));
359 /* check that we recieved at least one real (non-nan) value! */
360 if (kernel->minimum == MagickMaximumValue)
361 return(DestroyKernelInfo(kernel));
363 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
364 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
365 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
366 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
367 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
368 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
373 static KernelInfo *ParseKernelName(const char *kernel_string,
374 ExceptionInfo *exception)
377 token[MaxTextExtent];
395 /* Parse special 'named' kernel */
396 GetMagickToken(kernel_string,&p,token);
397 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
398 if ( type < 0 || type == UserDefinedKernel )
399 return((KernelInfo *) NULL); /* not a valid named kernel */
401 while (((isspace((int) ((unsigned char) *p)) != 0) ||
402 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
405 end = strchr(p, ';'); /* end of this kernel defintion */
406 if ( end == (char *) NULL )
407 end = strchr(p, '\0');
409 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
410 memcpy(token, p, (size_t) (end-p));
412 SetGeometryInfo(&args);
413 flags = ParseGeometry(token, &args);
416 /* For Debugging Geometry Input */
417 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
418 flags, args.rho, args.sigma, args.xi, args.psi );
421 /* special handling of missing values in input string */
423 /* Shape Kernel Defaults */
425 if ( (flags & WidthValue) == 0 )
426 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
434 if ( (flags & HeightValue) == 0 )
435 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
438 if ( (flags & XValue) == 0 )
439 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
441 case RectangleKernel: /* Rectangle - set size defaults */
442 if ( (flags & WidthValue) == 0 ) /* if no width then */
443 args.rho = args.sigma; /* then width = height */
444 if ( args.rho < 1.0 ) /* if width too small */
445 args.rho = 3; /* then width = 3 */
446 if ( args.sigma < 1.0 ) /* if height too small */
447 args.sigma = args.rho; /* then height = width */
448 if ( (flags & XValue) == 0 ) /* center offset if not defined */
449 args.xi = (double)(((ssize_t)args.rho-1)/2);
450 if ( (flags & YValue) == 0 )
451 args.psi = (double)(((ssize_t)args.sigma-1)/2);
453 /* Distance Kernel Defaults */
454 case ChebyshevKernel:
455 case ManhattanKernel:
456 case OctagonalKernel:
457 case EuclideanKernel:
458 if ( (flags & HeightValue) == 0 ) /* no distance scale */
459 args.sigma = 100.0; /* default distance scaling */
460 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
461 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
462 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
463 args.sigma *= QuantumRange/100.0; /* percentage of color range */
469 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args, exception);
470 if ( kernel == (KernelInfo *) NULL )
473 /* global expand to rotated kernel list - only for single kernels */
474 if ( kernel->next == (KernelInfo *) NULL ) {
475 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
476 ExpandRotateKernelInfo(kernel, 45.0);
477 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
478 ExpandRotateKernelInfo(kernel, 90.0);
479 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
480 ExpandMirrorKernelInfo(kernel);
486 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string,
487 ExceptionInfo *exception)
495 token[MaxTextExtent];
500 if (kernel_string == (const char *) NULL)
501 return(ParseKernelArray(kernel_string));
503 kernel_cache=(char *) NULL;
504 if (*kernel_string == '@')
506 kernel_cache=FileToString(kernel_string+1,~0UL,exception);
507 if (kernel_cache == (char *) NULL)
508 return((KernelInfo *) NULL);
509 p=(const char *) kernel_cache;
512 while (GetMagickToken(p,NULL,token), *token != '\0')
514 /* ignore extra or multiple ';' kernel separators */
517 /* tokens starting with alpha is a Named kernel */
518 if (isalpha((int) ((unsigned char) *token)) != 0)
519 new_kernel=ParseKernelName(p,exception);
520 else /* otherwise a user defined kernel array */
521 new_kernel=ParseKernelArray(p);
523 /* Error handling -- this is not proper error handling! */
524 if (new_kernel == (KernelInfo *) NULL)
526 if (kernel != (KernelInfo *) NULL)
527 kernel=DestroyKernelInfo(kernel);
528 return((KernelInfo *) NULL);
531 /* initialise or append the kernel list */
532 if (kernel == (KernelInfo *) NULL)
535 LastKernelInfo(kernel)->next=new_kernel;
538 /* look for the next kernel in list */
540 if (p == (char *) NULL)
544 if (kernel_cache != (char *) NULL)
545 kernel_cache=DestroyString(kernel_cache);
550 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
554 % A c q u i r e K e r n e l B u i l t I n %
558 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
560 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
561 % kernels used for special purposes such as gaussian blurring, skeleton
562 % pruning, and edge distance determination.
564 % They take a KernelType, and a set of geometry style arguments, which were
565 % typically decoded from a user supplied string, or from a more complex
566 % Morphology Method that was requested.
568 % The format of the AcquireKernalBuiltIn method is:
570 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
571 % const GeometryInfo args)
573 % A description of each parameter follows:
575 % o type: the pre-defined type of kernel wanted
577 % o args: arguments defining or modifying the kernel
579 % Convolution Kernels
582 % The a No-Op or Scaling single element kernel.
584 % Gaussian:{radius},{sigma}
585 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
586 % The sigma for the curve is required. The resulting kernel is
589 % If 'sigma' is zero, you get a single pixel on a field of zeros.
591 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
592 % the final size of the resulting kernel to a square 2*radius+1 in size.
593 % The radius should be at least 2 times that of the sigma value, or
594 % sever clipping and aliasing may result. If not given or set to 0 the
595 % radius will be determined so as to produce the best minimal error
596 % result, which is usally much larger than is normally needed.
598 % LoG:{radius},{sigma}
599 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
600 % The supposed ideal edge detection, zero-summing kernel.
602 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
603 % approx 1.6 (according to wikipedia).
605 % DoG:{radius},{sigma1},{sigma2}
606 % "Difference of Gaussians" Kernel.
607 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
608 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
609 % The result is a zero-summing kernel.
611 % Blur:{radius},{sigma}[,{angle}]
612 % Generates a 1 dimensional or linear gaussian blur, at the angle given
613 % (current restricted to orthogonal angles). If a 'radius' is given the
614 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
615 % by a 90 degree angle.
617 % If 'sigma' is zero, you get a single pixel on a field of zeros.
619 % Note that two convolutions with two "Blur" kernels perpendicular to
620 % each other, is equivalent to a far larger "Gaussian" kernel with the
621 % same sigma value, However it is much faster to apply. This is how the
622 % "-blur" operator actually works.
624 % Comet:{width},{sigma},{angle}
625 % Blur in one direction only, much like how a bright object leaves
626 % a comet like trail. The Kernel is actually half a gaussian curve,
627 % Adding two such blurs in opposite directions produces a Blur Kernel.
628 % Angle can be rotated in multiples of 90 degrees.
630 % Note that the first argument is the width of the kernel and not the
631 % radius of the kernel.
633 % Binomial:[{radius}]
634 % Generate a discrete kernel using a 2 dimentional Pascel's Triangle
635 % of values. Used for special forma of image filters.
637 % # Still to be implemented...
641 % # Set kernel values using a resize filter, and given scale (sigma)
642 % # Cylindrical or Linear. Is this possible with an image?
645 % Named Constant Convolution Kernels
647 % All these are unscaled, zero-summing kernels by default. As such for
648 % non-HDRI version of ImageMagick some form of normalization, user scaling,
649 % and biasing the results is recommended, to prevent the resulting image
652 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
653 % 45 degrees to generate the 8 angled varients of each of the kernels.
656 % Discrete Lapacian Kernels, (without normalization)
657 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
658 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
659 % Type 2 : 3x3 with center:4 edge:1 corner:-2
660 % Type 3 : 3x3 with center:4 edge:-2 corner:1
661 % Type 5 : 5x5 laplacian
662 % Type 7 : 7x7 laplacian
663 % Type 15 : 5x5 LoG (sigma approx 1.4)
664 % Type 19 : 9x9 LoG (sigma approx 1.4)
667 % Sobel 'Edge' convolution kernel (3x3)
673 % Roberts convolution kernel (3x3)
679 % Prewitt Edge convolution kernel (3x3)
685 % Prewitt's "Compass" convolution kernel (3x3)
691 % Kirsch's "Compass" convolution kernel (3x3)
697 % Frei-Chen Edge Detector is based on a kernel that is similar to
698 % the Sobel Kernel, but is designed to be isotropic. That is it takes
699 % into account the distance of the diagonal in the kernel.
702 % | sqrt(2), 0, -sqrt(2) |
705 % FreiChen:{type},{angle}
707 % Frei-Chen Pre-weighted kernels...
709 % Type 0: default un-nomalized version shown above.
711 % Type 1: Orthogonal Kernel (same as type 11 below)
713 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
716 % Type 2: Diagonal form of Kernel...
718 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
721 % However this kernel is als at the heart of the FreiChen Edge Detection
722 % Process which uses a set of 9 specially weighted kernel. These 9
723 % kernels not be normalized, but directly applied to the image. The
724 % results is then added together, to produce the intensity of an edge in
725 % a specific direction. The square root of the pixel value can then be
726 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
727 % from each other, both the direction and the strength of the edge can be
730 % Type 10: All 9 of the following pre-weighted kernels...
732 % Type 11: | 1, 0, -1 |
733 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
736 % Type 12: | 1, sqrt(2), 1 |
737 % | 0, 0, 0 | / 2*sqrt(2)
740 % Type 13: | sqrt(2), -1, 0 |
741 % | -1, 0, 1 | / 2*sqrt(2)
744 % Type 14: | 0, 1, -sqrt(2) |
745 % | -1, 0, 1 | / 2*sqrt(2)
748 % Type 15: | 0, -1, 0 |
752 % Type 16: | 1, 0, -1 |
756 % Type 17: | 1, -2, 1 |
760 % Type 18: | -2, 1, -2 |
764 % Type 19: | 1, 1, 1 |
768 % The first 4 are for edge detection, the next 4 are for line detection
769 % and the last is to add a average component to the results.
771 % Using a special type of '-1' will return all 9 pre-weighted kernels
772 % as a multi-kernel list, so that you can use them directly (without
773 % normalization) with the special "-set option:morphology:compose Plus"
774 % setting to apply the full FreiChen Edge Detection Technique.
776 % If 'type' is large it will be taken to be an actual rotation angle for
777 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
778 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
780 % WARNING: The above was layed out as per
781 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
782 % But rotated 90 degrees so direction is from left rather than the top.
783 % I have yet to find any secondary confirmation of the above. The only
784 % other source found was actual source code at
785 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
786 % Neigher paper defineds the kernels in a way that looks locical or
787 % correct when taken as a whole.
791 % Diamond:[{radius}[,{scale}]]
792 % Generate a diamond shaped kernel with given radius to the points.
793 % Kernel size will again be radius*2+1 square and defaults to radius 1,
794 % generating a 3x3 kernel that is slightly larger than a square.
796 % Square:[{radius}[,{scale}]]
797 % Generate a square shaped kernel of size radius*2+1, and defaulting
798 % to a 3x3 (radius 1).
800 % Octagon:[{radius}[,{scale}]]
801 % Generate octagonal shaped kernel of given radius and constant scale.
802 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
803 % in "Diamond" kernel.
805 % Disk:[{radius}[,{scale}]]
806 % Generate a binary disk, thresholded at the radius given, the radius
807 % may be a float-point value. Final Kernel size is floor(radius)*2+1
808 % square. A radius of 5.3 is the default.
810 % NOTE: That a low radii Disk kernels produce the same results as
811 % many of the previously defined kernels, but differ greatly at larger
812 % radii. Here is a table of equivalences...
813 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
814 % "Disk:1.5" => "Square"
815 % "Disk:2" => "Diamond:2"
816 % "Disk:2.5" => "Octagon"
817 % "Disk:2.9" => "Square:2"
818 % "Disk:3.5" => "Octagon:3"
819 % "Disk:4.5" => "Octagon:4"
820 % "Disk:5.4" => "Octagon:5"
821 % "Disk:6.4" => "Octagon:6"
822 % All other Disk shapes are unique to this kernel, but because a "Disk"
823 % is more circular when using a larger radius, using a larger radius is
824 % preferred over iterating the morphological operation.
826 % Rectangle:{geometry}
827 % Simply generate a rectangle of 1's with the size given. You can also
828 % specify the location of the 'control point', otherwise the closest
829 % pixel to the center of the rectangle is selected.
831 % Properly centered and odd sized rectangles work the best.
833 % Symbol Dilation Kernels
835 % These kernel is not a good general morphological kernel, but is used
836 % more for highlighting and marking any single pixels in an image using,
837 % a "Dilate" method as appropriate.
839 % For the same reasons iterating these kernels does not produce the
840 % same result as using a larger radius for the symbol.
842 % Plus:[{radius}[,{scale}]]
843 % Cross:[{radius}[,{scale}]]
844 % Generate a kernel in the shape of a 'plus' or a 'cross' with
845 % a each arm the length of the given radius (default 2).
847 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
849 % Ring:{radius1},{radius2}[,{scale}]
850 % A ring of the values given that falls between the two radii.
851 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
852 % This is the 'edge' pixels of the default "Disk" kernel,
853 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
855 % Hit and Miss Kernels
857 % Peak:radius1,radius2
858 % Find any peak larger than the pixels the fall between the two radii.
859 % The default ring of pixels is as per "Ring".
861 % Find flat orthogonal edges of a binary shape
863 % Find 90 degree corners of a binary shape
865 % A special kernel to thin the 'outside' of diagonals
867 % Find end points of lines (for pruning a skeletion)
868 % Two types of lines ends (default to both) can be searched for
869 % Type 0: All line ends
870 % Type 1: single kernel for 4-conneected line ends
871 % Type 2: single kernel for simple line ends
873 % Find three line junctions (within a skeletion)
874 % Type 0: all line junctions
875 % Type 1: Y Junction kernel
876 % Type 2: Diagonal T Junction kernel
877 % Type 3: Orthogonal T Junction kernel
878 % Type 4: Diagonal X Junction kernel
879 % Type 5: Orthogonal + Junction kernel
881 % Find single pixel ridges or thin lines
882 % Type 1: Fine single pixel thick lines and ridges
883 % Type 2: Find two pixel thick lines and ridges
885 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
887 % Traditional skeleton generating kernels.
888 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
889 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
890 % Type 3: Thinning skeleton based on a ressearch paper by
891 % Dan S. Bloomberg (Default Type)
893 % A huge variety of Thinning Kernels designed to preserve conectivity.
894 % many other kernel sets use these kernels as source definitions.
895 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
896 % the super and sub notations used in the source research paper.
898 % Distance Measuring Kernels
900 % Different types of distance measuring methods, which are used with the
901 % a 'Distance' morphology method for generating a gradient based on
902 % distance from an edge of a binary shape, though there is a technique
903 % for handling a anti-aliased shape.
905 % See the 'Distance' Morphological Method, for information of how it is
908 % Chebyshev:[{radius}][x{scale}[%!]]
909 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
910 % is a value of one to any neighbour, orthogonal or diagonal. One why
911 % of thinking of it is the number of squares a 'King' or 'Queen' in
912 % chess needs to traverse reach any other position on a chess board.
913 % It results in a 'square' like distance function, but one where
914 % diagonals are given a value that is closer than expected.
916 % Manhattan:[{radius}][x{scale}[%!]]
917 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
918 % Cab distance metric), it is the distance needed when you can only
919 % travel in horizontal or vertical directions only. It is the
920 % distance a 'Rook' in chess would have to travel, and results in a
921 % diamond like distances, where diagonals are further than expected.
923 % Octagonal:[{radius}][x{scale}[%!]]
924 % An interleving of Manhatten and Chebyshev metrics producing an
925 % increasing octagonally shaped distance. Distances matches those of
926 % the "Octagon" shaped kernel of the same radius. The minimum radius
927 % and default is 2, producing a 5x5 kernel.
929 % Euclidean:[{radius}][x{scale}[%!]]
930 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
931 % However by default the kernel size only has a radius of 1, which
932 % limits the distance to 'Knight' like moves, with only orthogonal and
933 % diagonal measurements being correct. As such for the default kernel
934 % you will get octagonal like distance function.
936 % However using a larger radius such as "Euclidean:4" you will get a
937 % much smoother distance gradient from the edge of the shape. Especially
938 % if the image is pre-processed to include any anti-aliasing pixels.
939 % Of course a larger kernel is slower to use, and not always needed.
941 % The first three Distance Measuring Kernels will only generate distances
942 % of exact multiples of {scale} in binary images. As such you can use a
943 % scale of 1 without loosing any information. However you also need some
944 % scaling when handling non-binary anti-aliased shapes.
946 % The "Euclidean" Distance Kernel however does generate a non-integer
947 % fractional results, and as such scaling is vital even for binary shapes.
951 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
952 const GeometryInfo *args,ExceptionInfo *exception)
965 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
967 /* Generate a new empty kernel if needed */
968 kernel=(KernelInfo *) NULL;
970 case UndefinedKernel: /* These should not call this function */
971 case UserDefinedKernel:
972 assert("Should not call this function" != (char *) NULL);
974 case LaplacianKernel: /* Named Descrete Convolution Kernels */
975 case SobelKernel: /* these are defined using other kernels */
981 case EdgesKernel: /* Hit and Miss kernels */
983 case DiagonalsKernel:
985 case LineJunctionsKernel:
987 case ConvexHullKernel:
990 break; /* A pre-generated kernel is not needed */
992 /* set to 1 to do a compile-time check that we haven't missed anything */
1002 case RectangleKernel:
1009 case ChebyshevKernel:
1010 case ManhattanKernel:
1011 case OctangonalKernel:
1012 case EuclideanKernel:
1016 /* Generate the base Kernel Structure */
1017 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1018 if (kernel == (KernelInfo *) NULL)
1020 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1021 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1022 kernel->negative_range = kernel->positive_range = 0.0;
1023 kernel->type = type;
1024 kernel->next = (KernelInfo *) NULL;
1025 kernel->signature = MagickSignature;
1035 kernel->height = kernel->width = (size_t) 1;
1036 kernel->x = kernel->y = (ssize_t) 0;
1037 kernel->values=(MagickRealType *) MagickAssumeAligned(
1038 AcquireAlignedMemory(1,sizeof(*kernel->values)));
1039 if (kernel->values == (MagickRealType *) NULL)
1040 return(DestroyKernelInfo(kernel));
1041 kernel->maximum = kernel->values[0] = args->rho;
1045 case GaussianKernel:
1049 sigma = fabs(args->sigma),
1050 sigma2 = fabs(args->xi),
1053 if ( args->rho >= 1.0 )
1054 kernel->width = (size_t)args->rho*2+1;
1055 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1056 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1058 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1059 kernel->height = kernel->width;
1060 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1061 kernel->values=(MagickRealType *) MagickAssumeAligned(
1062 AcquireAlignedMemory(kernel->width,kernel->height*
1063 sizeof(*kernel->values)));
1064 if (kernel->values == (MagickRealType *) NULL)
1065 return(DestroyKernelInfo(kernel));
1067 /* WARNING: The following generates a 'sampled gaussian' kernel.
1068 * What we really want is a 'discrete gaussian' kernel.
1070 * How to do this is I don't know, but appears to be basied on the
1071 * Error Function 'erf()' (intergral of a gaussian)
1074 if ( type == GaussianKernel || type == DoGKernel )
1075 { /* Calculate a Gaussian, OR positive half of a DoG */
1076 if ( sigma > MagickEpsilon )
1077 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1078 B = (double) (1.0/(Magick2PI*sigma*sigma));
1079 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1080 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1081 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1083 else /* limiting case - a unity (normalized Dirac) kernel */
1084 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1085 kernel->width*kernel->height*sizeof(*kernel->values));
1086 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1090 if ( type == DoGKernel )
1091 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1092 if ( sigma2 > MagickEpsilon )
1093 { sigma = sigma2; /* simplify loop expressions */
1094 A = 1.0/(2.0*sigma*sigma);
1095 B = (double) (1.0/(Magick2PI*sigma*sigma));
1096 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1097 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1098 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1100 else /* limiting case - a unity (normalized Dirac) kernel */
1101 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1104 if ( type == LoGKernel )
1105 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1106 if ( sigma > MagickEpsilon )
1107 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1108 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1109 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1110 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1111 { R = ((double)(u*u+v*v))*A;
1112 kernel->values[i] = (1-R)*exp(-R)*B;
1115 else /* special case - generate a unity kernel */
1116 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1117 kernel->width*kernel->height*sizeof(*kernel->values));
1118 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1122 /* Note the above kernels may have been 'clipped' by a user defined
1123 ** radius, producing a smaller (darker) kernel. Also for very small
1124 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1125 ** producing a very bright kernel.
1127 ** Normalization will still be needed.
1130 /* Normalize the 2D Gaussian Kernel
1132 ** NB: a CorrelateNormalize performs a normal Normalize if
1133 ** there are no negative values.
1135 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1136 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1142 sigma = fabs(args->sigma),
1145 if ( args->rho >= 1.0 )
1146 kernel->width = (size_t)args->rho*2+1;
1148 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1150 kernel->x = (ssize_t) (kernel->width-1)/2;
1152 kernel->negative_range = kernel->positive_range = 0.0;
1153 kernel->values=(MagickRealType *) MagickAssumeAligned(
1154 AcquireAlignedMemory(kernel->width,kernel->height*
1155 sizeof(*kernel->values)));
1156 if (kernel->values == (MagickRealType *) NULL)
1157 return(DestroyKernelInfo(kernel));
1160 #define KernelRank 3
1161 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1162 ** It generates a gaussian 3 times the width, and compresses it into
1163 ** the expected range. This produces a closer normalization of the
1164 ** resulting kernel, especially for very low sigma values.
1165 ** As such while wierd it is prefered.
1167 ** I am told this method originally came from Photoshop.
1169 ** A properly normalized curve is generated (apart from edge clipping)
1170 ** even though we later normalize the result (for edge clipping)
1171 ** to allow the correct generation of a "Difference of Blurs".
1175 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1176 (void) ResetMagickMemory(kernel->values,0, (size_t)
1177 kernel->width*kernel->height*sizeof(*kernel->values));
1178 /* Calculate a Positive 1D Gaussian */
1179 if ( sigma > MagickEpsilon )
1180 { sigma *= KernelRank; /* simplify loop expressions */
1181 alpha = 1.0/(2.0*sigma*sigma);
1182 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1183 for ( u=-v; u <= v; u++) {
1184 kernel->values[(u+v)/KernelRank] +=
1185 exp(-((double)(u*u))*alpha)*beta;
1188 else /* special case - generate a unity kernel */
1189 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1191 /* Direct calculation without curve averaging
1192 This is equivelent to a KernelRank of 1 */
1194 /* Calculate a Positive Gaussian */
1195 if ( sigma > MagickEpsilon )
1196 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1197 beta = 1.0/(MagickSQ2PI*sigma);
1198 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1199 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1201 else /* special case - generate a unity kernel */
1202 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1203 kernel->width*kernel->height*sizeof(*kernel->values));
1204 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1207 /* Note the above kernel may have been 'clipped' by a user defined
1208 ** radius, producing a smaller (darker) kernel. Also for very small
1209 ** sigma's (> 0.1) the central value becomes larger than one, as a
1210 ** result of not generating a actual 'discrete' kernel, and thus
1211 ** producing a very bright 'impulse'.
1213 ** Becuase of these two factors Normalization is required!
1216 /* Normalize the 1D Gaussian Kernel
1218 ** NB: a CorrelateNormalize performs a normal Normalize if
1219 ** there are no negative values.
1221 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1222 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1224 /* rotate the 1D kernel by given angle */
1225 RotateKernelInfo(kernel, args->xi );
1230 sigma = fabs(args->sigma),
1233 if ( args->rho < 1.0 )
1234 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1236 kernel->width = (size_t)args->rho;
1237 kernel->x = kernel->y = 0;
1239 kernel->negative_range = kernel->positive_range = 0.0;
1240 kernel->values=(MagickRealType *) MagickAssumeAligned(
1241 AcquireAlignedMemory(kernel->width,kernel->height*
1242 sizeof(*kernel->values)));
1243 if (kernel->values == (MagickRealType *) NULL)
1244 return(DestroyKernelInfo(kernel));
1246 /* A comet blur is half a 1D gaussian curve, so that the object is
1247 ** blurred in one direction only. This may not be quite the right
1248 ** curve to use so may change in the future. The function must be
1249 ** normalised after generation, which also resolves any clipping.
1251 ** As we are normalizing and not subtracting gaussians,
1252 ** there is no need for a divisor in the gaussian formula
1254 ** It is less comples
1256 if ( sigma > MagickEpsilon )
1259 #define KernelRank 3
1260 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1261 (void) ResetMagickMemory(kernel->values,0, (size_t)
1262 kernel->width*sizeof(*kernel->values));
1263 sigma *= KernelRank; /* simplify the loop expression */
1264 A = 1.0/(2.0*sigma*sigma);
1265 /* B = 1.0/(MagickSQ2PI*sigma); */
1266 for ( u=0; u < v; u++) {
1267 kernel->values[u/KernelRank] +=
1268 exp(-((double)(u*u))*A);
1269 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1271 for (i=0; i < (ssize_t) kernel->width; i++)
1272 kernel->positive_range += kernel->values[i];
1274 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1275 /* B = 1.0/(MagickSQ2PI*sigma); */
1276 for ( i=0; i < (ssize_t) kernel->width; i++)
1277 kernel->positive_range +=
1278 kernel->values[i] = exp(-((double)(i*i))*A);
1279 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1282 else /* special case - generate a unity kernel */
1283 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1284 kernel->width*kernel->height*sizeof(*kernel->values));
1285 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1286 kernel->positive_range = 1.0;
1289 kernel->minimum = 0.0;
1290 kernel->maximum = kernel->values[0];
1291 kernel->negative_range = 0.0;
1293 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1294 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1297 case BinomialKernel:
1302 if (args->rho < 1.0)
1303 kernel->width = kernel->height = 3; /* default radius = 1 */
1305 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1306 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1308 order_f = fact(kernel->width-1);
1310 kernel->values=(MagickRealType *) MagickAssumeAligned(
1311 AcquireAlignedMemory(kernel->width,kernel->height*
1312 sizeof(*kernel->values)));
1313 if (kernel->values == (MagickRealType *) NULL)
1314 return(DestroyKernelInfo(kernel));
1316 /* set all kernel values within diamond area to scale given */
1317 for ( i=0, v=0; v < (ssize_t)kernel->height; v++)
1319 alpha = order_f / ( fact((size_t) v) * fact(kernel->height-v-1) );
1320 for ( u=0; u < (ssize_t)kernel->width; u++, i++)
1321 kernel->positive_range += kernel->values[i] = (double)
1322 (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-u-1) ));
1324 kernel->minimum = 1.0;
1325 kernel->maximum = kernel->values[kernel->x+kernel->y*kernel->width];
1326 kernel->negative_range = 0.0;
1331 Convolution Kernels - Well Known Named Constant Kernels
1333 case LaplacianKernel:
1334 { switch ( (int) args->rho ) {
1336 default: /* laplacian square filter -- default */
1337 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1339 case 1: /* laplacian diamond filter */
1340 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1343 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1346 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1348 case 5: /* a 5x5 laplacian */
1349 kernel=ParseKernelArray(
1350 "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");
1352 case 7: /* a 7x7 laplacian */
1353 kernel=ParseKernelArray(
1354 "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" );
1356 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1357 kernel=ParseKernelArray(
1358 "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");
1360 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1361 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1362 kernel=ParseKernelArray(
1363 "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");
1366 if (kernel == (KernelInfo *) NULL)
1368 kernel->type = type;
1372 { /* Simple Sobel Kernel */
1373 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1374 if (kernel == (KernelInfo *) NULL)
1376 kernel->type = type;
1377 RotateKernelInfo(kernel, args->rho);
1382 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1383 if (kernel == (KernelInfo *) NULL)
1385 kernel->type = type;
1386 RotateKernelInfo(kernel, args->rho);
1391 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1392 if (kernel == (KernelInfo *) NULL)
1394 kernel->type = type;
1395 RotateKernelInfo(kernel, args->rho);
1400 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1401 if (kernel == (KernelInfo *) NULL)
1403 kernel->type = type;
1404 RotateKernelInfo(kernel, args->rho);
1409 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1410 if (kernel == (KernelInfo *) NULL)
1412 kernel->type = type;
1413 RotateKernelInfo(kernel, args->rho);
1416 case FreiChenKernel:
1417 /* Direction is set to be left to right positive */
1418 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1419 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1420 { switch ( (int) args->rho ) {
1423 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1424 if (kernel == (KernelInfo *) NULL)
1426 kernel->type = type;
1427 kernel->values[3] = +(MagickRealType) MagickSQ2;
1428 kernel->values[5] = -(MagickRealType) MagickSQ2;
1429 CalcKernelMetaData(kernel); /* recalculate meta-data */
1432 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1433 if (kernel == (KernelInfo *) NULL)
1435 kernel->type = type;
1436 kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2;
1437 kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2;
1438 CalcKernelMetaData(kernel); /* recalculate meta-data */
1439 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1443 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19",exception);
1444 if (kernel == (KernelInfo *) NULL)
1450 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1451 if (kernel == (KernelInfo *) NULL)
1453 kernel->type = type;
1454 kernel->values[3] = +(MagickRealType) MagickSQ2;
1455 kernel->values[5] = -(MagickRealType) MagickSQ2;
1456 CalcKernelMetaData(kernel); /* recalculate meta-data */
1457 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1460 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1461 if (kernel == (KernelInfo *) NULL)
1463 kernel->type = type;
1464 kernel->values[1] = +(MagickRealType) MagickSQ2;
1465 kernel->values[7] = +(MagickRealType) MagickSQ2;
1466 CalcKernelMetaData(kernel);
1467 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1470 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1471 if (kernel == (KernelInfo *) NULL)
1473 kernel->type = type;
1474 kernel->values[0] = +(MagickRealType) MagickSQ2;
1475 kernel->values[8] = -(MagickRealType) MagickSQ2;
1476 CalcKernelMetaData(kernel);
1477 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1480 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1481 if (kernel == (KernelInfo *) NULL)
1483 kernel->type = type;
1484 kernel->values[2] = -(MagickRealType) MagickSQ2;
1485 kernel->values[6] = +(MagickRealType) MagickSQ2;
1486 CalcKernelMetaData(kernel);
1487 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1490 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1491 if (kernel == (KernelInfo *) NULL)
1493 kernel->type = type;
1494 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1497 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1498 if (kernel == (KernelInfo *) NULL)
1500 kernel->type = type;
1501 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1504 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1505 if (kernel == (KernelInfo *) NULL)
1507 kernel->type = type;
1508 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1511 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1512 if (kernel == (KernelInfo *) NULL)
1514 kernel->type = type;
1515 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1518 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1519 if (kernel == (KernelInfo *) NULL)
1521 kernel->type = type;
1522 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1525 if ( fabs(args->sigma) >= MagickEpsilon )
1526 /* Rotate by correctly supplied 'angle' */
1527 RotateKernelInfo(kernel, args->sigma);
1528 else if ( args->rho > 30.0 || args->rho < -30.0 )
1529 /* Rotate by out of bounds 'type' */
1530 RotateKernelInfo(kernel, args->rho);
1535 Boolean or Shaped Kernels
1539 if (args->rho < 1.0)
1540 kernel->width = kernel->height = 3; /* default radius = 1 */
1542 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1543 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1545 kernel->values=(MagickRealType *) MagickAssumeAligned(
1546 AcquireAlignedMemory(kernel->width,kernel->height*
1547 sizeof(*kernel->values)));
1548 if (kernel->values == (MagickRealType *) NULL)
1549 return(DestroyKernelInfo(kernel));
1551 /* set all kernel values within diamond area to scale given */
1552 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1553 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1554 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1555 kernel->positive_range += kernel->values[i] = args->sigma;
1557 kernel->values[i] = nan;
1558 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1562 case RectangleKernel:
1565 if ( type == SquareKernel )
1567 if (args->rho < 1.0)
1568 kernel->width = kernel->height = 3; /* default radius = 1 */
1570 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1571 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1572 scale = args->sigma;
1575 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1576 if ( args->rho < 1.0 || args->sigma < 1.0 )
1577 return(DestroyKernelInfo(kernel)); /* invalid args given */
1578 kernel->width = (size_t)args->rho;
1579 kernel->height = (size_t)args->sigma;
1580 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1581 args->psi < 0.0 || args->psi > (double)kernel->height )
1582 return(DestroyKernelInfo(kernel)); /* invalid args given */
1583 kernel->x = (ssize_t) args->xi;
1584 kernel->y = (ssize_t) args->psi;
1587 kernel->values=(MagickRealType *) MagickAssumeAligned(
1588 AcquireAlignedMemory(kernel->width,kernel->height*
1589 sizeof(*kernel->values)));
1590 if (kernel->values == (MagickRealType *) NULL)
1591 return(DestroyKernelInfo(kernel));
1593 /* set all kernel values to scale given */
1594 u=(ssize_t) (kernel->width*kernel->height);
1595 for ( i=0; i < u; i++)
1596 kernel->values[i] = scale;
1597 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1598 kernel->positive_range = scale*u;
1603 if (args->rho < 1.0)
1604 kernel->width = kernel->height = 5; /* default radius = 2 */
1606 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1607 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1609 kernel->values=(MagickRealType *) MagickAssumeAligned(
1610 AcquireAlignedMemory(kernel->width,kernel->height*
1611 sizeof(*kernel->values)));
1612 if (kernel->values == (MagickRealType *) NULL)
1613 return(DestroyKernelInfo(kernel));
1615 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1616 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1617 if ( (labs((long) u)+labs((long) v)) <=
1618 ((long)kernel->x + (long)(kernel->x/2)) )
1619 kernel->positive_range += kernel->values[i] = args->sigma;
1621 kernel->values[i] = nan;
1622 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1628 limit = (ssize_t)(args->rho*args->rho);
1630 if (args->rho < 0.4) /* default radius approx 4.3 */
1631 kernel->width = kernel->height = 9L, limit = 18L;
1633 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1634 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1636 kernel->values=(MagickRealType *) MagickAssumeAligned(
1637 AcquireAlignedMemory(kernel->width,kernel->height*
1638 sizeof(*kernel->values)));
1639 if (kernel->values == (MagickRealType *) NULL)
1640 return(DestroyKernelInfo(kernel));
1642 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1643 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1644 if ((u*u+v*v) <= limit)
1645 kernel->positive_range += kernel->values[i] = args->sigma;
1647 kernel->values[i] = nan;
1648 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1653 if (args->rho < 1.0)
1654 kernel->width = kernel->height = 5; /* default radius 2 */
1656 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1657 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1659 kernel->values=(MagickRealType *) MagickAssumeAligned(
1660 AcquireAlignedMemory(kernel->width,kernel->height*
1661 sizeof(*kernel->values)));
1662 if (kernel->values == (MagickRealType *) NULL)
1663 return(DestroyKernelInfo(kernel));
1665 /* set all kernel values along axises to given scale */
1666 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1667 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1668 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1669 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1670 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1675 if (args->rho < 1.0)
1676 kernel->width = kernel->height = 5; /* default radius 2 */
1678 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1679 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1681 kernel->values=(MagickRealType *) MagickAssumeAligned(
1682 AcquireAlignedMemory(kernel->width,kernel->height*
1683 sizeof(*kernel->values)));
1684 if (kernel->values == (MagickRealType *) NULL)
1685 return(DestroyKernelInfo(kernel));
1687 /* set all kernel values along axises to given scale */
1688 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1689 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1690 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1691 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1692 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1706 if (args->rho < args->sigma)
1708 kernel->width = ((size_t)args->sigma)*2+1;
1709 limit1 = (ssize_t)(args->rho*args->rho);
1710 limit2 = (ssize_t)(args->sigma*args->sigma);
1714 kernel->width = ((size_t)args->rho)*2+1;
1715 limit1 = (ssize_t)(args->sigma*args->sigma);
1716 limit2 = (ssize_t)(args->rho*args->rho);
1719 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1721 kernel->height = kernel->width;
1722 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1723 kernel->values=(MagickRealType *) MagickAssumeAligned(
1724 AcquireAlignedMemory(kernel->width,kernel->height*
1725 sizeof(*kernel->values)));
1726 if (kernel->values == (MagickRealType *) NULL)
1727 return(DestroyKernelInfo(kernel));
1729 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1730 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1731 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1732 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1733 { ssize_t radius=u*u+v*v;
1734 if (limit1 < radius && radius <= limit2)
1735 kernel->positive_range += kernel->values[i] = (double) scale;
1737 kernel->values[i] = nan;
1739 kernel->minimum = kernel->maximum = (double) scale;
1740 if ( type == PeaksKernel ) {
1741 /* set the central point in the middle */
1742 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1743 kernel->positive_range = 1.0;
1744 kernel->maximum = 1.0;
1750 kernel=AcquireKernelInfo("ThinSE:482",exception);
1751 if (kernel == (KernelInfo *) NULL)
1753 kernel->type = type;
1754 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1759 kernel=AcquireKernelInfo("ThinSE:87",exception);
1760 if (kernel == (KernelInfo *) NULL)
1762 kernel->type = type;
1763 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1766 case DiagonalsKernel:
1768 switch ( (int) args->rho ) {
1773 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1774 if (kernel == (KernelInfo *) NULL)
1776 kernel->type = type;
1777 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1778 if (new_kernel == (KernelInfo *) NULL)
1779 return(DestroyKernelInfo(kernel));
1780 new_kernel->type = type;
1781 LastKernelInfo(kernel)->next = new_kernel;
1782 ExpandMirrorKernelInfo(kernel);
1786 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1789 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1792 if (kernel == (KernelInfo *) NULL)
1794 kernel->type = type;
1795 RotateKernelInfo(kernel, args->sigma);
1798 case LineEndsKernel:
1799 { /* Kernels for finding the end of thin lines */
1800 switch ( (int) args->rho ) {
1803 /* set of kernels to find all end of lines */
1804 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>",exception));
1806 /* kernel for 4-connected line ends - no rotation */
1807 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1810 /* kernel to add for 8-connected lines - no rotation */
1811 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1814 /* kernel to add for orthogonal line ends - does not find corners */
1815 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1818 /* traditional line end - fails on last T end */
1819 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1822 if (kernel == (KernelInfo *) NULL)
1824 kernel->type = type;
1825 RotateKernelInfo(kernel, args->sigma);
1828 case LineJunctionsKernel:
1829 { /* kernels for finding the junctions of multiple lines */
1830 switch ( (int) args->rho ) {
1833 /* set of kernels to find all line junctions */
1834 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>",exception));
1837 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1840 /* Diagonal T Junctions */
1841 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1844 /* Orthogonal T Junctions */
1845 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1848 /* Diagonal X Junctions */
1849 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1852 /* Orthogonal X Junctions - minimal diamond kernel */
1853 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1856 if (kernel == (KernelInfo *) NULL)
1858 kernel->type = type;
1859 RotateKernelInfo(kernel, args->sigma);
1863 { /* Ridges - Ridge finding kernels */
1866 switch ( (int) args->rho ) {
1869 kernel=ParseKernelArray("3x1:0,1,0");
1870 if (kernel == (KernelInfo *) NULL)
1872 kernel->type = type;
1873 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1876 kernel=ParseKernelArray("4x1:0,1,1,0");
1877 if (kernel == (KernelInfo *) NULL)
1879 kernel->type = type;
1880 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1882 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1883 /* Unfortunatally we can not yet rotate a non-square kernel */
1884 /* But then we can't flip a non-symetrical kernel either */
1885 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1886 if (new_kernel == (KernelInfo *) NULL)
1887 return(DestroyKernelInfo(kernel));
1888 new_kernel->type = type;
1889 LastKernelInfo(kernel)->next = new_kernel;
1890 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1891 if (new_kernel == (KernelInfo *) NULL)
1892 return(DestroyKernelInfo(kernel));
1893 new_kernel->type = type;
1894 LastKernelInfo(kernel)->next = new_kernel;
1895 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1896 if (new_kernel == (KernelInfo *) NULL)
1897 return(DestroyKernelInfo(kernel));
1898 new_kernel->type = type;
1899 LastKernelInfo(kernel)->next = new_kernel;
1900 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1901 if (new_kernel == (KernelInfo *) NULL)
1902 return(DestroyKernelInfo(kernel));
1903 new_kernel->type = type;
1904 LastKernelInfo(kernel)->next = new_kernel;
1905 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1906 if (new_kernel == (KernelInfo *) NULL)
1907 return(DestroyKernelInfo(kernel));
1908 new_kernel->type = type;
1909 LastKernelInfo(kernel)->next = new_kernel;
1910 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1911 if (new_kernel == (KernelInfo *) NULL)
1912 return(DestroyKernelInfo(kernel));
1913 new_kernel->type = type;
1914 LastKernelInfo(kernel)->next = new_kernel;
1915 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1916 if (new_kernel == (KernelInfo *) NULL)
1917 return(DestroyKernelInfo(kernel));
1918 new_kernel->type = type;
1919 LastKernelInfo(kernel)->next = new_kernel;
1920 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1921 if (new_kernel == (KernelInfo *) NULL)
1922 return(DestroyKernelInfo(kernel));
1923 new_kernel->type = type;
1924 LastKernelInfo(kernel)->next = new_kernel;
1929 case ConvexHullKernel:
1933 /* first set of 8 kernels */
1934 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1935 if (kernel == (KernelInfo *) NULL)
1937 kernel->type = type;
1938 ExpandRotateKernelInfo(kernel, 90.0);
1939 /* append the mirror versions too - no flip function yet */
1940 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1941 if (new_kernel == (KernelInfo *) NULL)
1942 return(DestroyKernelInfo(kernel));
1943 new_kernel->type = type;
1944 ExpandRotateKernelInfo(new_kernel, 90.0);
1945 LastKernelInfo(kernel)->next = new_kernel;
1948 case SkeletonKernel:
1950 switch ( (int) args->rho ) {
1953 /* Traditional Skeleton...
1954 ** A cyclically rotated single kernel
1956 kernel=AcquireKernelInfo("ThinSE:482",exception);
1957 if (kernel == (KernelInfo *) NULL)
1959 kernel->type = type;
1960 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1963 /* HIPR Variation of the cyclic skeleton
1964 ** Corners of the traditional method made more forgiving,
1965 ** but the retain the same cyclic order.
1967 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;",exception);
1968 if (kernel == (KernelInfo *) NULL)
1970 if (kernel->next == (KernelInfo *) NULL)
1971 return(DestroyKernelInfo(kernel));
1972 kernel->type = type;
1973 kernel->next->type = type;
1974 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1977 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1978 ** "Connectivity-Preserving Morphological Image Thransformations"
1979 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1980 ** http://www.leptonica.com/papers/conn.pdf
1982 kernel=AcquireKernelInfo("ThinSE:41; ThinSE:42; ThinSE:43",
1984 if (kernel == (KernelInfo *) NULL)
1986 kernel->type = type;
1987 kernel->next->type = type;
1988 kernel->next->next->type = type;
1989 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1995 { /* Special kernels for general thinning, while preserving connections
1996 ** "Connectivity-Preserving Morphological Image Thransformations"
1997 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1998 ** http://www.leptonica.com/papers/conn.pdf
2000 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2002 ** Note kernels do not specify the origin pixel, allowing them
2003 ** to be used for both thickening and thinning operations.
2005 switch ( (int) args->rho ) {
2006 /* SE for 4-connected thinning */
2007 case 41: /* SE_4_1 */
2008 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
2010 case 42: /* SE_4_2 */
2011 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
2013 case 43: /* SE_4_3 */
2014 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
2016 case 44: /* SE_4_4 */
2017 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
2019 case 45: /* SE_4_5 */
2020 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
2022 case 46: /* SE_4_6 */
2023 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
2025 case 47: /* SE_4_7 */
2026 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
2028 case 48: /* SE_4_8 */
2029 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
2031 case 49: /* SE_4_9 */
2032 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
2034 /* SE for 8-connected thinning - negatives of the above */
2035 case 81: /* SE_8_0 */
2036 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
2038 case 82: /* SE_8_2 */
2039 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
2041 case 83: /* SE_8_3 */
2042 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
2044 case 84: /* SE_8_4 */
2045 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
2047 case 85: /* SE_8_5 */
2048 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2050 case 86: /* SE_8_6 */
2051 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2053 case 87: /* SE_8_7 */
2054 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2056 case 88: /* SE_8_8 */
2057 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2059 case 89: /* SE_8_9 */
2060 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2062 /* Special combined SE kernels */
2063 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2064 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2066 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2067 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2069 case 481: /* SE_48_1 - General Connected Corner Kernel */
2070 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2073 case 482: /* SE_48_2 - General Edge Kernel */
2074 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2077 if (kernel == (KernelInfo *) NULL)
2079 kernel->type = type;
2080 RotateKernelInfo(kernel, args->sigma);
2084 Distance Measuring Kernels
2086 case ChebyshevKernel:
2088 if (args->rho < 1.0)
2089 kernel->width = kernel->height = 3; /* default radius = 1 */
2091 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2092 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2094 kernel->values=(MagickRealType *) MagickAssumeAligned(
2095 AcquireAlignedMemory(kernel->width,kernel->height*
2096 sizeof(*kernel->values)));
2097 if (kernel->values == (MagickRealType *) NULL)
2098 return(DestroyKernelInfo(kernel));
2100 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2101 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2102 kernel->positive_range += ( kernel->values[i] =
2103 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2104 kernel->maximum = kernel->values[0];
2107 case ManhattanKernel:
2109 if (args->rho < 1.0)
2110 kernel->width = kernel->height = 3; /* default radius = 1 */
2112 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2113 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2115 kernel->values=(MagickRealType *) MagickAssumeAligned(
2116 AcquireAlignedMemory(kernel->width,kernel->height*
2117 sizeof(*kernel->values)));
2118 if (kernel->values == (MagickRealType *) NULL)
2119 return(DestroyKernelInfo(kernel));
2121 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2122 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2123 kernel->positive_range += ( kernel->values[i] =
2124 args->sigma*(labs((long) u)+labs((long) v)) );
2125 kernel->maximum = kernel->values[0];
2128 case OctagonalKernel:
2130 if (args->rho < 2.0)
2131 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2133 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2134 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2136 kernel->values=(MagickRealType *) MagickAssumeAligned(
2137 AcquireAlignedMemory(kernel->width,kernel->height*
2138 sizeof(*kernel->values)));
2139 if (kernel->values == (MagickRealType *) NULL)
2140 return(DestroyKernelInfo(kernel));
2142 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2143 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2146 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2147 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2148 kernel->positive_range += kernel->values[i] =
2149 args->sigma*MagickMax(r1,r2);
2151 kernel->maximum = kernel->values[0];
2154 case EuclideanKernel:
2156 if (args->rho < 1.0)
2157 kernel->width = kernel->height = 3; /* default radius = 1 */
2159 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2160 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2162 kernel->values=(MagickRealType *) MagickAssumeAligned(
2163 AcquireAlignedMemory(kernel->width,kernel->height*
2164 sizeof(*kernel->values)));
2165 if (kernel->values == (MagickRealType *) NULL)
2166 return(DestroyKernelInfo(kernel));
2168 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2169 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2170 kernel->positive_range += ( kernel->values[i] =
2171 args->sigma*sqrt((double)(u*u+v*v)) );
2172 kernel->maximum = kernel->values[0];
2177 /* No-Op Kernel - Basically just a single pixel on its own */
2178 kernel=ParseKernelArray("1:1");
2179 if (kernel == (KernelInfo *) NULL)
2181 kernel->type = UndefinedKernel;
2190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2194 % C l o n e K e r n e l I n f o %
2198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2200 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2201 % can be modified without effecting the original. The cloned kernel should
2202 % be destroyed using DestoryKernelInfo() when no longer needed.
2204 % The format of the CloneKernelInfo method is:
2206 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2208 % A description of each parameter follows:
2210 % o kernel: the Morphology/Convolution kernel to be cloned
2213 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2221 assert(kernel != (KernelInfo *) NULL);
2222 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2223 if (new_kernel == (KernelInfo *) NULL)
2225 *new_kernel=(*kernel); /* copy values in structure */
2227 /* replace the values with a copy of the values */
2228 new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2229 AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2230 if (new_kernel->values == (MagickRealType *) NULL)
2231 return(DestroyKernelInfo(new_kernel));
2232 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2233 new_kernel->values[i]=kernel->values[i];
2235 /* Also clone the next kernel in the kernel list */
2236 if ( kernel->next != (KernelInfo *) NULL ) {
2237 new_kernel->next = CloneKernelInfo(kernel->next);
2238 if ( new_kernel->next == (KernelInfo *) NULL )
2239 return(DestroyKernelInfo(new_kernel));
2246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2250 % D e s t r o y K e r n e l I n f o %
2254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2256 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2259 % The format of the DestroyKernelInfo method is:
2261 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2263 % A description of each parameter follows:
2265 % o kernel: the Morphology/Convolution kernel to be destroyed
2268 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2270 assert(kernel != (KernelInfo *) NULL);
2271 if (kernel->next != (KernelInfo *) NULL)
2272 kernel->next=DestroyKernelInfo(kernel->next);
2273 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2274 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2279 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2283 + E x p a n d M i r r o r K e r n e l I n f o %
2287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2289 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2290 % sequence of 90-degree rotated kernels but providing a reflected 180
2291 % rotatation, before the -/+ 90-degree rotations.
2293 % This special rotation order produces a better, more symetrical thinning of
2296 % The format of the ExpandMirrorKernelInfo method is:
2298 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2300 % A description of each parameter follows:
2302 % o kernel: the Morphology/Convolution kernel
2304 % This function is only internel to this module, as it is not finalized,
2305 % especially with regard to non-orthogonal angles, and rotation of larger
2310 static void FlopKernelInfo(KernelInfo *kernel)
2311 { /* Do a Flop by reversing each row. */
2319 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2320 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2321 t=k[x], k[x]=k[r], k[r]=t;
2323 kernel->x = kernel->width - kernel->x - 1;
2324 angle = fmod(angle+180.0, 360.0);
2328 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2336 clone = CloneKernelInfo(last);
2337 RotateKernelInfo(clone, 180); /* flip */
2338 LastKernelInfo(last)->next = clone;
2341 clone = CloneKernelInfo(last);
2342 RotateKernelInfo(clone, 90); /* transpose */
2343 LastKernelInfo(last)->next = clone;
2346 clone = CloneKernelInfo(last);
2347 RotateKernelInfo(clone, 180); /* flop */
2348 LastKernelInfo(last)->next = clone;
2354 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2358 + E x p a n d R o t a t e K e r n e l I n f o %
2362 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2364 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2365 % incrementally by the angle given, until the kernel repeats.
2367 % WARNING: 45 degree rotations only works for 3x3 kernels.
2368 % While 90 degree roatations only works for linear and square kernels
2370 % The format of the ExpandRotateKernelInfo method is:
2372 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2374 % A description of each parameter follows:
2376 % o kernel: the Morphology/Convolution kernel
2378 % o angle: angle to rotate in degrees
2380 % This function is only internel to this module, as it is not finalized,
2381 % especially with regard to non-orthogonal angles, and rotation of larger
2385 /* Internal Routine - Return true if two kernels are the same */
2386 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2387 const KernelInfo *kernel2)
2392 /* check size and origin location */
2393 if ( kernel1->width != kernel2->width
2394 || kernel1->height != kernel2->height
2395 || kernel1->x != kernel2->x
2396 || kernel1->y != kernel2->y )
2399 /* check actual kernel values */
2400 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2401 /* Test for Nan equivalence */
2402 if ( IfNaN(kernel1->values[i]) && !IfNaN(kernel2->values[i]) )
2404 if ( IfNaN(kernel2->values[i]) && !IfNaN(kernel1->values[i]) )
2406 /* Test actual values are equivalent */
2407 if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2414 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2421 DisableMSCWarning(4127)
2424 clone = CloneKernelInfo(last);
2425 RotateKernelInfo(clone, angle);
2426 if ( SameKernelInfo(kernel, clone) != MagickFalse )
2428 LastKernelInfo(last)->next = clone;
2431 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2436 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2440 + C a l c M e t a K e r n a l I n f o %
2444 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2446 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2447 % using the kernel values. This should only ne used if it is not possible to
2448 % calculate that meta-data in some easier way.
2450 % It is important that the meta-data is correct before ScaleKernelInfo() is
2451 % used to perform kernel normalization.
2453 % The format of the CalcKernelMetaData method is:
2455 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2457 % A description of each parameter follows:
2459 % o kernel: the Morphology/Convolution kernel to modify
2461 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2462 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2463 % however is not true for flat-shaped morphological kernels.
2465 % WARNING: Only the specific kernel pointed to is modified, not a list of
2468 % This is an internal function and not expected to be useful outside this
2469 % module. This could change however.
2471 static void CalcKernelMetaData(KernelInfo *kernel)
2476 kernel->minimum = kernel->maximum = 0.0;
2477 kernel->negative_range = kernel->positive_range = 0.0;
2478 for (i=0; i < (kernel->width*kernel->height); i++)
2480 if ( fabs(kernel->values[i]) < MagickEpsilon )
2481 kernel->values[i] = 0.0;
2482 ( kernel->values[i] < 0)
2483 ? ( kernel->negative_range += kernel->values[i] )
2484 : ( kernel->positive_range += kernel->values[i] );
2485 Minimize(kernel->minimum, kernel->values[i]);
2486 Maximize(kernel->maximum, kernel->values[i]);
2493 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2497 % M o r p h o l o g y A p p l y %
2501 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2503 % MorphologyApply() applies a morphological method, multiple times using
2504 % a list of multiple kernels. This is the method that should be called by
2505 % other 'operators' that internally use morphology operations as part of
2508 % It is basically equivalent to as MorphologyImage() (see below) but without
2509 % any user controls. This allows internel programs to use this method to
2510 % perform a specific task without possible interference by any API user
2511 % supplied settings.
2513 % It is MorphologyImage() task to extract any such user controls, and
2514 % pass them to this function for processing.
2516 % More specifically all given kernels should already be scaled, normalised,
2517 % and blended appropriatally before being parred to this routine. The
2518 % appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2520 % The format of the MorphologyApply method is:
2522 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2523 % const ssize_t iterations,const KernelInfo *kernel,
2524 % const CompositeMethod compose,const double bias,
2525 % ExceptionInfo *exception)
2527 % A description of each parameter follows:
2529 % o image: the source image
2531 % o method: the morphology method to be applied.
2533 % o iterations: apply the operation this many times (or no change).
2534 % A value of -1 means loop until no change found.
2535 % How this is applied may depend on the morphology method.
2536 % Typically this is a value of 1.
2538 % o channel: the channel type.
2540 % o kernel: An array of double representing the morphology kernel.
2542 % o compose: How to handle or merge multi-kernel results.
2543 % If 'UndefinedCompositeOp' use default for the Morphology method.
2544 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2545 % Otherwise merge the results using the compose method given.
2547 % o bias: Convolution Output Bias.
2549 % o exception: return any errors or warnings in this structure.
2552 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2553 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2554 ExceptionInfo *exception)
2556 #define MorphologyTag "Morphology/Image"
2582 assert(image != (Image *) NULL);
2583 assert(image->signature == MagickSignature);
2584 assert(morphology_image != (Image *) NULL);
2585 assert(morphology_image->signature == MagickSignature);
2586 assert(kernel != (KernelInfo *) NULL);
2587 assert(kernel->signature == MagickSignature);
2588 assert(exception != (ExceptionInfo *) NULL);
2589 assert(exception->signature == MagickSignature);
2592 image_view=AcquireVirtualCacheView(image,exception);
2593 morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2594 width=image->columns+kernel->width-1;
2599 case ConvolveMorphology:
2600 case DilateMorphology:
2601 case DilateIntensityMorphology:
2602 case IterativeDistanceMorphology:
2605 Kernel needs to used with reflection about origin.
2607 offset.x=(ssize_t) kernel->width-kernel->x-1;
2608 offset.y=(ssize_t) kernel->height-kernel->y-1;
2611 case ErodeMorphology:
2612 case ErodeIntensityMorphology:
2613 case HitAndMissMorphology:
2614 case ThinningMorphology:
2615 case ThickenMorphology:
2623 assert("Not a Primitive Morphology Method" != (char *) NULL);
2628 changes=(size_t *) AcquireQuantumMemory(GetOpenMPMaximumThreads(),
2630 if (changes == (size_t *) NULL)
2631 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2632 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2634 if ((method == ConvolveMorphology) && (kernel->width == 1))
2640 Special handling (for speed) of vertical (blur) kernels. This performs
2641 its handling in columns rather than in rows. This is only done
2642 for convolve as it is the only method that generates very large 1-D
2643 vertical kernels (such as a 'BlurKernel')
2645 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2646 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2647 magick_threads(image,morphology_image,image->columns,1)
2649 for (x=0; x < (ssize_t) image->columns; x++)
2652 id = GetOpenMPThreadId();
2654 register const Quantum
2666 if (status == MagickFalse)
2668 p=GetCacheViewVirtualPixels(image_view,x,-offset.y,1,image->rows+
2669 kernel->height-1,exception);
2670 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2671 morphology_image->rows,exception);
2672 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2677 center=(ssize_t) GetPixelChannels(image)*offset.y;
2678 for (y=0; y < (ssize_t) image->rows; y++)
2683 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2697 register const MagickRealType
2700 register const Quantum
2712 channel=GetPixelChannelChannel(image,i);
2713 traits=GetPixelChannelTraits(image,channel);
2714 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2715 if ((traits == UndefinedPixelTrait) ||
2716 (morphology_traits == UndefinedPixelTrait))
2718 if (((traits & CopyPixelTrait) != 0) ||
2719 (GetPixelReadMask(image,p+center) == 0))
2721 SetPixelChannel(morphology_image,channel,p[center+i],q);
2724 k=(&kernel->values[kernel->width*kernel->height-1]);
2729 if ((morphology_traits & BlendPixelTrait) == 0)
2730 for (v=0; v < (ssize_t) kernel->height; v++)
2732 for (u=0; u < (ssize_t) kernel->width; u++)
2734 if (IfNaN(*k) == MagickFalse)
2736 pixel+=(*k)*pixels[i];
2741 pixels+=GetPixelChannels(image);
2745 for (v=0; v < (ssize_t) kernel->height; v++)
2747 for (u=0; u < (ssize_t) kernel->width; u++)
2749 if (IfNaN(*k) == MagickFalse)
2751 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2752 pixel+=alpha*(*k)*pixels[i];
2757 pixels+=GetPixelChannels(image);
2760 if (fabs(pixel-p[center+i]) > MagickEpsilon)
2762 gamma=PerceptibleReciprocal(gamma);
2764 gamma*=(double) kernel->height*kernel->width/count;
2765 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2768 p+=GetPixelChannels(image);
2769 q+=GetPixelChannels(morphology_image);
2771 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2773 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2778 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2779 #pragma omp critical (MagickCore_MorphologyPrimitive)
2781 proceed=SetImageProgress(image,MorphologyTag,progress++,
2783 if (proceed == MagickFalse)
2787 morphology_image->type=image->type;
2788 morphology_view=DestroyCacheView(morphology_view);
2789 image_view=DestroyCacheView(image_view);
2790 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2791 changed+=changes[i];
2792 changes=(size_t *) RelinquishMagickMemory(changes);
2793 return(status ? (ssize_t) changed : 0);
2796 Normal handling of horizontal or rectangular kernels (row by row).
2798 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2799 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2800 magick_threads(image,morphology_image,image->rows,1)
2802 for (y=0; y < (ssize_t) image->rows; y++)
2805 id = GetOpenMPThreadId();
2807 register const Quantum
2819 if (status == MagickFalse)
2821 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,
2822 kernel->height,exception);
2823 q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2825 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2830 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
2831 GetPixelChannels(image)*offset.x);
2832 for (x=0; x < (ssize_t) image->columns; x++)
2837 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2854 register const MagickRealType
2857 register const Quantum
2869 channel=GetPixelChannelChannel(image,i);
2870 traits=GetPixelChannelTraits(image,channel);
2871 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2872 if ((traits == UndefinedPixelTrait) ||
2873 (morphology_traits == UndefinedPixelTrait))
2875 if (((traits & CopyPixelTrait) != 0) ||
2876 (GetPixelReadMask(image,p+center) == 0))
2878 SetPixelChannel(morphology_image,channel,p[center+i],q);
2883 minimum=(double) QuantumRange;
2884 count=kernel->width*kernel->height;
2887 case ConvolveMorphology: pixel=bias; break;
2888 case HitAndMissMorphology: pixel=(double) QuantumRange; break;
2889 case ThinningMorphology: pixel=(double) QuantumRange; break;
2890 case ThickenMorphology: pixel=(double) QuantumRange; break;
2891 case ErodeMorphology: pixel=(double) QuantumRange; break;
2892 case DilateMorphology: pixel=0.0; break;
2893 case ErodeIntensityMorphology:
2894 case DilateIntensityMorphology:
2895 case IterativeDistanceMorphology:
2897 pixel=(double) p[center+i];
2900 default: pixel=0; break;
2905 case ConvolveMorphology:
2908 Weighted Average of pixels using reflected kernel
2910 For correct working of this operation for asymetrical kernels,
2911 the kernel needs to be applied in its reflected form. That is
2912 its values needs to be reversed.
2914 Correlation is actually the same as this but without reflecting
2915 the kernel, and thus 'lower-level' that Convolution. However as
2916 Convolution is the more common method used, and it does not
2917 really cost us much in terms of processing to use a reflected
2918 kernel, so it is Convolution that is implemented.
2920 Correlation will have its kernel reflected before calling this
2921 function to do a Convolve.
2923 For more details of Correlation vs Convolution see
2924 http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2926 k=(&kernel->values[kernel->width*kernel->height-1]);
2928 if ((morphology_traits & BlendPixelTrait) == 0)
2933 for (v=0; v < (ssize_t) kernel->height; v++)
2935 for (u=0; u < (ssize_t) kernel->width; u++)
2937 if (IfNaN(*k) == MagickFalse)
2939 pixel+=(*k)*pixels[i];
2943 pixels+=GetPixelChannels(image);
2945 pixels+=(image->columns-1)*GetPixelChannels(image);
2953 for (v=0; v < (ssize_t) kernel->height; v++)
2955 for (u=0; u < (ssize_t) kernel->width; u++)
2957 if (IfNaN(*k) == MagickFalse)
2959 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2960 pixel+=alpha*(*k)*pixels[i];
2965 pixels+=GetPixelChannels(image);
2967 pixels+=(image->columns-1)*GetPixelChannels(image);
2971 case ErodeMorphology:
2974 Minimum value within kernel neighbourhood.
2976 The kernel is not reflected for this operation. In normal
2977 Greyscale Morphology, the kernel value should be added
2978 to the real value, this is currently not done, due to the
2979 nature of the boolean kernels being used.
2982 for (v=0; v < (ssize_t) kernel->height; v++)
2984 for (u=0; u < (ssize_t) kernel->width; u++)
2986 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
2988 if ((double) pixels[i] < pixel)
2989 pixel=(double) pixels[i];
2992 pixels+=GetPixelChannels(image);
2994 pixels+=(image->columns-1)*GetPixelChannels(image);
2998 case DilateMorphology:
3001 Maximum value within kernel neighbourhood.
3003 For correct working of this operation for asymetrical kernels,
3004 the kernel needs to be applied in its reflected form. That is
3005 its values needs to be reversed.
3007 In normal Greyscale Morphology, the kernel value should be
3008 added to the real value, this is currently not done, due to the
3009 nature of the boolean kernels being used.
3012 k=(&kernel->values[kernel->width*kernel->height-1]);
3013 for (v=0; v < (ssize_t) kernel->height; v++)
3015 for (u=0; u < (ssize_t) kernel->width; u++)
3017 if ((IfNaN(*k) == MagickFalse) && (*k > 0.5))
3019 if ((double) pixels[i] > pixel)
3020 pixel=(double) pixels[i];
3023 pixels+=GetPixelChannels(image);
3025 pixels+=(image->columns-1)*GetPixelChannels(image);
3029 case HitAndMissMorphology:
3030 case ThinningMorphology:
3031 case ThickenMorphology:
3034 Minimum of foreground pixel minus maxumum of background pixels.
3036 The kernel is not reflected for this operation, and consists
3037 of both foreground and background pixel neighbourhoods, 0.0 for
3038 background, and 1.0 for foreground with either Nan or 0.5 values
3041 This never produces a meaningless negative result. Such results
3042 cause Thinning/Thicken to not work correctly when used against a
3047 for (v=0; v < (ssize_t) kernel->height; v++)
3049 for (u=0; u < (ssize_t) kernel->width; u++)
3051 if (IfNaN(*k) == MagickFalse)
3055 if ((double) pixels[i] < pixel)
3056 pixel=(double) pixels[i];
3061 if ((double) pixels[i] > maximum)
3062 maximum=(double) pixels[i];
3067 pixels+=GetPixelChannels(image);
3069 pixels+=(image->columns-1)*GetPixelChannels(image);
3074 if (method == ThinningMorphology)
3075 pixel=(double) p[center+i]-pixel;
3077 if (method == ThickenMorphology)
3078 pixel+=(double) p[center+i]+pixel;
3081 case ErodeIntensityMorphology:
3084 Select pixel with minimum intensity within kernel neighbourhood.
3086 The kernel is not reflected for this operation.
3090 for (v=0; v < (ssize_t) kernel->height; v++)
3092 for (u=0; u < (ssize_t) kernel->width; u++)
3094 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3096 intensity=(double) GetPixelIntensity(image,pixels);
3097 if (intensity < minimum)
3099 pixel=(double) pixels[i];
3105 pixels+=GetPixelChannels(image);
3107 pixels+=(image->columns-1)*GetPixelChannels(image);
3111 case DilateIntensityMorphology:
3114 Select pixel with maximum intensity within kernel neighbourhood.
3116 The kernel is not reflected for this operation.
3119 k=(&kernel->values[kernel->width*kernel->height-1]);
3120 for (v=0; v < (ssize_t) kernel->height; v++)
3122 for (u=0; u < (ssize_t) kernel->width; u++)
3124 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3126 intensity=(double) GetPixelIntensity(image,pixels);
3127 if (intensity > maximum)
3129 pixel=(double) pixels[i];
3135 pixels+=GetPixelChannels(image);
3137 pixels+=(image->columns-1)*GetPixelChannels(image);
3141 case IterativeDistanceMorphology:
3144 Compute th iterative distance from black edge of a white image
3145 shape. Essentually white values are decreased to the smallest
3146 'distance from edge' it can find.
3148 It works by adding kernel values to the neighbourhood, and and
3149 select the minimum value found. The kernel is rotated before
3150 use, so kernel distances match resulting distances, when a user
3151 provided asymmetric kernel is applied.
3153 This code is nearly identical to True GrayScale Morphology but
3156 GreyDilate Kernel values added, maximum value found Kernel is
3159 GrayErode: Kernel values subtracted and minimum value found No
3160 kernel rotation used.
3162 Note the the Iterative Distance method is essentially a
3163 GrayErode, but with negative kernel values, and kernel rotation
3167 k=(&kernel->values[kernel->width*kernel->height-1]);
3168 for (v=0; v < (ssize_t) kernel->height; v++)
3170 for (u=0; u < (ssize_t) kernel->width; u++)
3172 if (IfNaN(*k) == MagickFalse)
3174 if ((pixels[i]+(*k)) < pixel)
3175 pixel=(double) pixels[i]+(*k);
3179 pixels+=GetPixelChannels(image);
3181 pixels+=(image->columns-1)*GetPixelChannels(image);
3185 case UndefinedMorphology:
3189 if (fabs(pixel-p[center+i]) > MagickEpsilon)
3191 gamma=PerceptibleReciprocal(gamma);
3193 gamma*=(double) kernel->height*kernel->width/count;
3194 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*pixel),q);
3196 p+=GetPixelChannels(image);
3197 q+=GetPixelChannels(morphology_image);
3199 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3201 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3206 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3207 #pragma omp critical (MagickCore_MorphologyPrimitive)
3209 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3210 if (proceed == MagickFalse)
3214 morphology_view=DestroyCacheView(morphology_view);
3215 image_view=DestroyCacheView(image_view);
3216 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
3217 changed+=changes[i];
3218 changes=(size_t *) RelinquishMagickMemory(changes);
3219 return(status ? (ssize_t) changed : -1);
3223 This is almost identical to the MorphologyPrimative() function above, but
3224 applies the primitive directly to the actual image using two passes, once in
3225 each direction, with the results of the previous (and current) row being
3228 That is after each row is 'Sync'ed' into the image, the next row makes use of
3229 those values as part of the calculation of the next row. It repeats, but
3230 going in the oppisite (bottom-up) direction.
3232 Because of this 're-use of results' this function can not make use of multi-
3233 threaded, parellel processing.
3235 static ssize_t MorphologyPrimitiveDirect(Image *image,
3236 const MorphologyMethod method,const KernelInfo *kernel,
3237 ExceptionInfo *exception)
3259 assert(image != (Image *) NULL);
3260 assert(image->signature == MagickSignature);
3261 assert(kernel != (KernelInfo *) NULL);
3262 assert(kernel->signature == MagickSignature);
3263 assert(exception != (ExceptionInfo *) NULL);
3264 assert(exception->signature == MagickSignature);
3270 case DistanceMorphology:
3271 case VoronoiMorphology:
3274 Kernel reflected about origin.
3276 offset.x=(ssize_t) kernel->width-kernel->x-1;
3277 offset.y=(ssize_t) kernel->height-kernel->y-1;
3288 Two views into same image, do not thread.
3290 image_view=AcquireVirtualCacheView(image,exception);
3291 morphology_view=AcquireAuthenticCacheView(image,exception);
3292 width=image->columns+kernel->width-1;
3293 for (y=0; y < (ssize_t) image->rows; y++)
3295 register const Quantum
3308 Read virtual pixels, and authentic pixels, from the same image! We read
3309 using virtual to get virtual pixel handling, but write back into the same
3312 Only top half of kernel is processed as we do a single pass downward
3313 through the image iterating the distance function as we go.
3315 if (status == MagickFalse)
3317 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t)
3318 offset.y+1,exception);
3319 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3321 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3326 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
3327 GetPixelChannels(image)*offset.x);
3328 for (x=0; x < (ssize_t) image->columns; x++)
3333 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3341 register const MagickRealType
3344 register const Quantum
3353 traits=GetPixelChannelTraits(image,(PixelChannel) i);
3354 if (traits == UndefinedPixelTrait)
3356 if (((traits & CopyPixelTrait) != 0) ||
3357 (GetPixelReadMask(image,p+center) == 0))
3360 pixel=(double) QuantumRange;
3363 case DistanceMorphology:
3365 k=(&kernel->values[kernel->width*kernel->height-1]);
3366 for (v=0; v <= offset.y; v++)
3368 for (u=0; u < (ssize_t) kernel->width; u++)
3370 if (IfNaN(*k) == MagickFalse)
3372 if ((pixels[i]+(*k)) < pixel)
3373 pixel=(double) pixels[i]+(*k);
3376 pixels+=GetPixelChannels(image);
3378 pixels+=(image->columns-1)*GetPixelChannels(image);
3380 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3381 pixels=q-offset.x*GetPixelChannels(image);
3382 for (u=0; u < offset.x; u++)
3384 if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3386 if ((pixels[i]+(*k)) < pixel)
3387 pixel=(double) pixels[i]+(*k);
3390 pixels+=GetPixelChannels(image);
3394 case VoronoiMorphology:
3396 k=(&kernel->values[kernel->width*kernel->height-1]);
3397 for (v=0; v < offset.y; v++)
3399 for (u=0; u < (ssize_t) kernel->width; u++)
3401 if (IfNaN(*k) == MagickFalse)
3403 if ((pixels[i]+(*k)) < pixel)
3404 pixel=(double) pixels[i]+(*k);
3407 pixels+=GetPixelChannels(image);
3409 pixels+=(image->columns-1)*GetPixelChannels(image);
3411 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3412 pixels=q-offset.x*GetPixelChannels(image);
3413 for (u=0; u < offset.x; u++)
3415 if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3417 if ((pixels[i]+(*k)) < pixel)
3418 pixel=(double) pixels[i]+(*k);
3421 pixels+=GetPixelChannels(image);
3428 if (fabs(pixel-q[i]) > MagickEpsilon)
3430 q[i]=ClampToQuantum(pixel);
3432 p+=GetPixelChannels(image);
3433 q+=GetPixelChannels(image);
3435 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3437 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3442 proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3443 if (proceed == MagickFalse)
3447 morphology_view=DestroyCacheView(morphology_view);
3448 image_view=DestroyCacheView(image_view);
3450 Do the reverse pass through the image.
3452 image_view=AcquireVirtualCacheView(image,exception);
3453 morphology_view=AcquireAuthenticCacheView(image,exception);
3454 for (y=(ssize_t) image->rows-1; y >= 0; y--)
3456 register const Quantum
3469 Read virtual pixels, and authentic pixels, from the same image. We
3470 read using virtual to get virtual pixel handling, but write back
3471 into the same image.
3473 Only the bottom half of the kernel is processed as we up the image.
3475 if (status == MagickFalse)
3477 p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t)
3478 kernel->y+1,exception);
3479 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3481 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3486 p+=(image->columns-1)*GetPixelChannels(image);
3487 q+=(image->columns-1)*GetPixelChannels(image);
3488 center=(ssize_t) (offset.x*GetPixelChannels(image));
3489 for (x=(ssize_t) image->columns-1; x >= 0; x--)
3494 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3502 register const MagickRealType
3505 register const Quantum
3514 traits=GetPixelChannelTraits(image,(PixelChannel) i);
3515 if (traits == UndefinedPixelTrait)
3517 if (((traits & CopyPixelTrait) != 0) ||
3518 (GetPixelReadMask(image,p+center) == 0))
3521 pixel=(double) QuantumRange;
3524 case DistanceMorphology:
3526 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3527 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3529 for (u=0; u < (ssize_t) kernel->width; u++)
3531 if (IfNaN(*k) == MagickFalse)
3533 if ((pixels[i]+(*k)) < pixel)
3534 pixel=(double) pixels[i]+(*k);
3537 pixels+=GetPixelChannels(image);
3539 pixels+=(image->columns-1)*GetPixelChannels(image);
3541 k=(&kernel->values[kernel->width*kernel->y+kernel->x-1]);
3543 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3545 pixels+=GetPixelChannels(image);
3546 if ((IfNaN(*k) == MagickFalse) &&
3547 ((x+u-offset.x) < (ssize_t) image->columns))
3549 if ((pixels[i]+(*k)) < pixel)
3550 pixel=(double) pixels[i]+(*k);
3556 case VoronoiMorphology:
3558 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3559 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3561 for (u=0; u < (ssize_t) kernel->width; u++)
3563 if (IfNaN(*k) == MagickFalse)
3565 if ((pixels[i]+(*k)) < pixel)
3566 pixel=(double) pixels[i]+(*k);
3569 pixels+=GetPixelChannels(image);
3571 pixels+=(image->columns-1)*GetPixelChannels(image);
3573 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3575 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3577 pixels+=GetPixelChannels(image);
3578 if ((IfNaN(*k) == MagickFalse) &&
3579 ((x+u-offset.x) < (ssize_t) image->columns))
3581 if ((pixels[i]+(*k)) < pixel)
3582 pixel=(double) pixels[i]+(*k);
3591 if (fabs(pixel-q[i]) > MagickEpsilon)
3593 q[i]=ClampToQuantum(pixel);
3595 p-=GetPixelChannels(image);
3596 q-=GetPixelChannels(image);
3598 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3600 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3605 proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3606 if (proceed == MagickFalse)
3610 morphology_view=DestroyCacheView(morphology_view);
3611 image_view=DestroyCacheView(image_view);
3612 return(status ? (ssize_t) changed : -1);
3616 Apply a Morphology by calling one of the above low level primitive
3617 application functions. This function handles any iteration loops,
3618 composition or re-iteration of results, and compound morphology methods that
3619 is based on multiple low-level (staged) morphology methods.
3621 Basically this provides the complex glue between the requested morphology
3622 method and raw low-level implementation (above).
3624 MagickPrivate Image *MorphologyApply(const Image *image,
3625 const MorphologyMethod method, const ssize_t iterations,
3626 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3627 ExceptionInfo *exception)
3633 *curr_image, /* Image we are working with or iterating */
3634 *work_image, /* secondary image for primitive iteration */
3635 *save_image, /* saved image - for 'edge' method only */
3636 *rslt_image; /* resultant image - after multi-kernel handling */
3639 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3640 *norm_kernel, /* the current normal un-reflected kernel */
3641 *rflt_kernel, /* the current reflected kernel (if needed) */
3642 *this_kernel; /* the kernel being applied */
3645 primitive; /* the current morphology primitive being applied */
3648 rslt_compose; /* multi-kernel compose method for results to use */
3651 special, /* do we use a direct modify function? */
3652 verbose; /* verbose output of results */
3655 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3656 method_limit, /* maximum number of compound method iterations */
3657 kernel_number, /* Loop 2: the kernel number being applied */
3658 stage_loop, /* Loop 3: primitive loop for compound morphology */
3659 stage_limit, /* how many primitives are in this compound */
3660 kernel_loop, /* Loop 4: iterate the kernel over image */
3661 kernel_limit, /* number of times to iterate kernel */
3662 count, /* total count of primitive steps applied */
3663 kernel_changed, /* total count of changed using iterated kernel */
3664 method_changed; /* total count of changed over method iteration */
3667 changed; /* number pixels changed by last primitive operation */
3670 v_info[MaxTextExtent];
3672 assert(image != (Image *) NULL);
3673 assert(image->signature == MagickSignature);
3674 assert(kernel != (KernelInfo *) NULL);
3675 assert(kernel->signature == MagickSignature);
3676 assert(exception != (ExceptionInfo *) NULL);
3677 assert(exception->signature == MagickSignature);
3679 count = 0; /* number of low-level morphology primitives performed */
3680 if ( iterations == 0 )
3681 return((Image *) NULL); /* null operation - nothing to do! */
3683 kernel_limit = (size_t) iterations;
3684 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3685 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3687 verbose = IsStringTrue(GetImageArtifact(image,"debug"));
3689 /* initialise for cleanup */
3690 curr_image = (Image *) image;
3691 curr_compose = image->compose;
3692 (void) curr_compose;
3693 work_image = save_image = rslt_image = (Image *) NULL;
3694 reflected_kernel = (KernelInfo *) NULL;
3696 /* Initialize specific methods
3697 * + which loop should use the given iteratations
3698 * + how many primitives make up the compound morphology
3699 * + multi-kernel compose method to use (by default)
3701 method_limit = 1; /* just do method once, unless otherwise set */
3702 stage_limit = 1; /* assume method is not a compound */
3703 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3704 rslt_compose = compose; /* and we are composing multi-kernels as given */
3706 case SmoothMorphology: /* 4 primitive compound morphology */
3709 case OpenMorphology: /* 2 primitive compound morphology */
3710 case OpenIntensityMorphology:
3711 case TopHatMorphology:
3712 case CloseMorphology:
3713 case CloseIntensityMorphology:
3714 case BottomHatMorphology:
3715 case EdgeMorphology:
3718 case HitAndMissMorphology:
3719 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3721 case ThinningMorphology:
3722 case ThickenMorphology:
3723 method_limit = kernel_limit; /* iterate the whole method */
3724 kernel_limit = 1; /* do not do kernel iteration */
3726 case DistanceMorphology:
3727 case VoronoiMorphology:
3728 special = MagickTrue; /* use special direct primative */
3734 /* Apply special methods with special requirments
3735 ** For example, single run only, or post-processing requirements
3737 if ( special != MagickFalse )
3739 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3740 if (rslt_image == (Image *) NULL)
3742 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3745 changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception);
3747 if ( IfMagickTrue(verbose) )
3748 (void) (void) FormatLocaleFile(stderr,
3749 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3750 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3751 1.0,0.0,1.0, (double) changed);
3756 if ( method == VoronoiMorphology ) {
3757 /* Preserve the alpha channel of input image - but turned it off */
3758 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3760 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3761 MagickTrue,0,0,exception);
3762 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3768 /* Handle user (caller) specified multi-kernel composition method */
3769 if ( compose != UndefinedCompositeOp )
3770 rslt_compose = compose; /* override default composition for method */
3771 if ( rslt_compose == UndefinedCompositeOp )
3772 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3774 /* Some methods require a reflected kernel to use with primitives.
3775 * Create the reflected kernel for those methods. */
3777 case CorrelateMorphology:
3778 case CloseMorphology:
3779 case CloseIntensityMorphology:
3780 case BottomHatMorphology:
3781 case SmoothMorphology:
3782 reflected_kernel = CloneKernelInfo(kernel);
3783 if (reflected_kernel == (KernelInfo *) NULL)
3785 RotateKernelInfo(reflected_kernel,180);
3791 /* Loops around more primitive morpholgy methods
3792 ** erose, dilate, open, close, smooth, edge, etc...
3794 /* Loop 1: iterate the compound method */
3797 while ( method_loop < method_limit && method_changed > 0 ) {
3801 /* Loop 2: iterate over each kernel in a multi-kernel list */
3802 norm_kernel = (KernelInfo *) kernel;
3803 this_kernel = (KernelInfo *) kernel;
3804 rflt_kernel = reflected_kernel;
3807 while ( norm_kernel != NULL ) {
3809 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3810 stage_loop = 0; /* the compound morphology stage number */
3811 while ( stage_loop < stage_limit ) {
3812 stage_loop++; /* The stage of the compound morphology */
3814 /* Select primitive morphology for this stage of compound method */
3815 this_kernel = norm_kernel; /* default use unreflected kernel */
3816 primitive = method; /* Assume method is a primitive */
3818 case ErodeMorphology: /* just erode */
3819 case EdgeInMorphology: /* erode and image difference */
3820 primitive = ErodeMorphology;
3822 case DilateMorphology: /* just dilate */
3823 case EdgeOutMorphology: /* dilate and image difference */
3824 primitive = DilateMorphology;
3826 case OpenMorphology: /* erode then dialate */
3827 case TopHatMorphology: /* open and image difference */
3828 primitive = ErodeMorphology;
3829 if ( stage_loop == 2 )
3830 primitive = DilateMorphology;
3832 case OpenIntensityMorphology:
3833 primitive = ErodeIntensityMorphology;
3834 if ( stage_loop == 2 )
3835 primitive = DilateIntensityMorphology;
3837 case CloseMorphology: /* dilate, then erode */
3838 case BottomHatMorphology: /* close and image difference */
3839 this_kernel = rflt_kernel; /* use the reflected kernel */
3840 primitive = DilateMorphology;
3841 if ( stage_loop == 2 )
3842 primitive = ErodeMorphology;
3844 case CloseIntensityMorphology:
3845 this_kernel = rflt_kernel; /* use the reflected kernel */
3846 primitive = DilateIntensityMorphology;
3847 if ( stage_loop == 2 )
3848 primitive = ErodeIntensityMorphology;
3850 case SmoothMorphology: /* open, close */
3851 switch ( stage_loop ) {
3852 case 1: /* start an open method, which starts with Erode */
3853 primitive = ErodeMorphology;
3855 case 2: /* now Dilate the Erode */
3856 primitive = DilateMorphology;
3858 case 3: /* Reflect kernel a close */
3859 this_kernel = rflt_kernel; /* use the reflected kernel */
3860 primitive = DilateMorphology;
3862 case 4: /* Finish the Close */
3863 this_kernel = rflt_kernel; /* use the reflected kernel */
3864 primitive = ErodeMorphology;
3868 case EdgeMorphology: /* dilate and erode difference */
3869 primitive = DilateMorphology;
3870 if ( stage_loop == 2 ) {
3871 save_image = curr_image; /* save the image difference */
3872 curr_image = (Image *) image;
3873 primitive = ErodeMorphology;
3876 case CorrelateMorphology:
3877 /* A Correlation is a Convolution with a reflected kernel.
3878 ** However a Convolution is a weighted sum using a reflected
3879 ** kernel. It may seem stange to convert a Correlation into a
3880 ** Convolution as the Correlation is the simplier method, but
3881 ** Convolution is much more commonly used, and it makes sense to
3882 ** implement it directly so as to avoid the need to duplicate the
3883 ** kernel when it is not required (which is typically the
3886 this_kernel = rflt_kernel; /* use the reflected kernel */
3887 primitive = ConvolveMorphology;
3892 assert( this_kernel != (KernelInfo *) NULL );
3894 /* Extra information for debugging compound operations */
3895 if ( IfMagickTrue(verbose) ) {
3896 if ( stage_limit > 1 )
3897 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
3898 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
3899 method_loop,(double) stage_loop);
3900 else if ( primitive != method )
3901 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
3902 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
3908 /* Loop 4: Iterate the kernel with primitive */
3912 while ( kernel_loop < kernel_limit && changed > 0 ) {
3913 kernel_loop++; /* the iteration of this kernel */
3915 /* Create a clone as the destination image, if not yet defined */
3916 if ( work_image == (Image *) NULL )
3918 work_image=CloneImage(image,0,0,MagickTrue,exception);
3919 if (work_image == (Image *) NULL)
3921 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
3925 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3927 changed = MorphologyPrimitive(curr_image, work_image, primitive,
3928 this_kernel, bias, exception);
3929 if ( IfMagickTrue(verbose) ) {
3930 if ( kernel_loop > 1 )
3931 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
3932 (void) (void) FormatLocaleFile(stderr,
3933 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3934 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
3935 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
3936 (double) (method_loop+kernel_loop-1),(double) kernel_number,
3937 (double) count,(double) changed);
3941 kernel_changed += changed;
3942 method_changed += changed;
3944 /* prepare next loop */
3945 { Image *tmp = work_image; /* swap images for iteration */
3946 work_image = curr_image;
3949 if ( work_image == image )
3950 work_image = (Image *) NULL; /* replace input 'image' */
3952 } /* End Loop 4: Iterate the kernel with primitive */
3954 if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed )
3955 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
3956 if ( IfMagickTrue(verbose) && stage_loop < stage_limit )
3957 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
3960 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3961 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
3962 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
3963 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
3964 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
3967 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3969 /* Final Post-processing for some Compound Methods
3971 ** The removal of any 'Sync' channel flag in the Image Compositon
3972 ** below ensures the methematical compose method is applied in a
3973 ** purely mathematical way, and only to the selected channels.
3974 ** Turn off SVG composition 'alpha blending'.
3977 case EdgeOutMorphology:
3978 case EdgeInMorphology:
3979 case TopHatMorphology:
3980 case BottomHatMorphology:
3981 if ( IfMagickTrue(verbose) )
3982 (void) FormatLocaleFile(stderr,
3983 "\n%s: Difference with original image",CommandOptionToMnemonic(
3984 MagickMorphologyOptions, method) );
3985 (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
3986 MagickTrue,0,0,exception);
3988 case EdgeMorphology:
3989 if ( IfMagickTrue(verbose) )
3990 (void) FormatLocaleFile(stderr,
3991 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
3992 MagickMorphologyOptions, method) );
3993 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
3994 MagickTrue,0,0,exception);
3995 save_image = DestroyImage(save_image); /* finished with save image */
4001 /* multi-kernel handling: re-iterate, or compose results */
4002 if ( kernel->next == (KernelInfo *) NULL )
4003 rslt_image = curr_image; /* just return the resulting image */
4004 else if ( rslt_compose == NoCompositeOp )
4005 { if ( IfMagickTrue(verbose) ) {
4006 if ( this_kernel->next != (KernelInfo *) NULL )
4007 (void) FormatLocaleFile(stderr, " (re-iterate)");
4009 (void) FormatLocaleFile(stderr, " (done)");
4011 rslt_image = curr_image; /* return result, and re-iterate */
4013 else if ( rslt_image == (Image *) NULL)
4014 { if ( IfMagickTrue(verbose) )
4015 (void) FormatLocaleFile(stderr, " (save for compose)");
4016 rslt_image = curr_image;
4017 curr_image = (Image *) image; /* continue with original image */
4020 { /* Add the new 'current' result to the composition
4022 ** The removal of any 'Sync' channel flag in the Image Compositon
4023 ** below ensures the methematical compose method is applied in a
4024 ** purely mathematical way, and only to the selected channels.
4025 ** IE: Turn off SVG composition 'alpha blending'.
4027 if ( IfMagickTrue(verbose) )
4028 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4029 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4030 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4032 curr_image = DestroyImage(curr_image);
4033 curr_image = (Image *) image; /* continue with original image */
4035 if ( IfMagickTrue(verbose) )
4036 (void) FormatLocaleFile(stderr, "\n");
4038 /* loop to the next kernel in a multi-kernel list */
4039 norm_kernel = norm_kernel->next;
4040 if ( rflt_kernel != (KernelInfo *) NULL )
4041 rflt_kernel = rflt_kernel->next;
4043 } /* End Loop 2: Loop over each kernel */
4045 } /* End Loop 1: compound method interation */
4049 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4051 if ( curr_image == rslt_image )
4052 curr_image = (Image *) NULL;
4053 if ( rslt_image != (Image *) NULL )
4054 rslt_image = DestroyImage(rslt_image);
4056 if ( curr_image == rslt_image || curr_image == image )
4057 curr_image = (Image *) NULL;
4058 if ( curr_image != (Image *) NULL )
4059 curr_image = DestroyImage(curr_image);
4060 if ( work_image != (Image *) NULL )
4061 work_image = DestroyImage(work_image);
4062 if ( save_image != (Image *) NULL )
4063 save_image = DestroyImage(save_image);
4064 if ( reflected_kernel != (KernelInfo *) NULL )
4065 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4071 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4075 % M o r p h o l o g y I m a g e %
4079 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4081 % MorphologyImage() applies a user supplied kernel to the image according to
4082 % the given mophology method.
4084 % This function applies any and all user defined settings before calling
4085 % the above internal function MorphologyApply().
4087 % User defined settings include...
4088 % * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4089 % * Kernel Scale/normalize settings ("-define convolve:scale=??")
4090 % This can also includes the addition of a scaled unity kernel.
4091 % * Show Kernel being applied ("-define showkernel=1")
4093 % Other operators that do not want user supplied options interfering,
4094 % especially "convolve:bias" and "showkernel" should use MorphologyApply()
4097 % The format of the MorphologyImage method is:
4099 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4100 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4102 % A description of each parameter follows:
4104 % o image: the image.
4106 % o method: the morphology method to be applied.
4108 % o iterations: apply the operation this many times (or no change).
4109 % A value of -1 means loop until no change found.
4110 % How this is applied may depend on the morphology method.
4111 % Typically this is a value of 1.
4113 % o kernel: An array of double representing the morphology kernel.
4114 % Warning: kernel may be normalized for the Convolve method.
4116 % o exception: return any errors or warnings in this structure.
4119 MagickExport Image *MorphologyImage(const Image *image,
4120 const MorphologyMethod method,const ssize_t iterations,
4121 const KernelInfo *kernel,ExceptionInfo *exception)
4135 curr_kernel = (KernelInfo *) kernel;
4137 compose = UndefinedCompositeOp; /* use default for method */
4139 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4140 * This is done BEFORE the ShowKernelInfo() function is called so that
4141 * users can see the results of the 'option:convolve:scale' option.
4143 if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4147 /* Get the bias value as it will be needed */
4148 artifact = GetImageArtifact(image,"convolve:bias");
4149 if ( artifact != (const char *) NULL) {
4150 if (IfMagickFalse(IsGeometry(artifact)))
4151 (void) ThrowMagickException(exception,GetMagickModule(),
4152 OptionWarning,"InvalidSetting","'%s' '%s'",
4153 "convolve:bias",artifact);
4155 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4158 /* Scale kernel according to user wishes */
4159 artifact = GetImageArtifact(image,"convolve:scale");
4160 if ( artifact != (const char *) NULL ) {
4161 if (IfMagickFalse(IsGeometry(artifact)))
4162 (void) ThrowMagickException(exception,GetMagickModule(),
4163 OptionWarning,"InvalidSetting","'%s' '%s'",
4164 "convolve:scale",artifact);
4166 if ( curr_kernel == kernel )
4167 curr_kernel = CloneKernelInfo(kernel);
4168 if (curr_kernel == (KernelInfo *) NULL)
4169 return((Image *) NULL);
4170 ScaleGeometryKernelInfo(curr_kernel, artifact);
4175 /* display the (normalized) kernel via stderr */
4176 if ( IfStringTrue(GetImageArtifact(image,"showkernel"))
4177 || IfStringTrue(GetImageArtifact(image,"convolve:showkernel"))
4178 || IfStringTrue(GetImageArtifact(image,"morphology:showkernel")) )
4179 ShowKernelInfo(curr_kernel);
4181 /* Override the default handling of multi-kernel morphology results
4182 * If 'Undefined' use the default method
4183 * If 'None' (default for 'Convolve') re-iterate previous result
4184 * Otherwise merge resulting images using compose method given.
4185 * Default for 'HitAndMiss' is 'Lighten'.
4192 artifact = GetImageArtifact(image,"morphology:compose");
4193 if ( artifact != (const char *) NULL) {
4194 parse=ParseCommandOption(MagickComposeOptions,
4195 MagickFalse,artifact);
4197 (void) ThrowMagickException(exception,GetMagickModule(),
4198 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4199 "morphology:compose",artifact);
4201 compose=(CompositeOperator)parse;
4204 /* Apply the Morphology */
4205 morphology_image = MorphologyApply(image,method,iterations,
4206 curr_kernel,compose,bias,exception);
4208 /* Cleanup and Exit */
4209 if ( curr_kernel != kernel )
4210 curr_kernel=DestroyKernelInfo(curr_kernel);
4211 return(morphology_image);
4215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4219 + R o t a t e K e r n e l I n f o %
4223 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4225 % RotateKernelInfo() rotates the kernel by the angle given.
4227 % Currently it is restricted to 90 degree angles, of either 1D kernels
4228 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4229 % It will ignore usless rotations for specific 'named' built-in kernels.
4231 % The format of the RotateKernelInfo method is:
4233 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4235 % A description of each parameter follows:
4237 % o kernel: the Morphology/Convolution kernel
4239 % o angle: angle to rotate in degrees
4241 % This function is currently internal to this module only, but can be exported
4242 % to other modules if needed.
4244 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4246 /* angle the lower kernels first */
4247 if ( kernel->next != (KernelInfo *) NULL)
4248 RotateKernelInfo(kernel->next, angle);
4250 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4252 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4255 /* Modulus the angle */
4256 angle = fmod(angle, 360.0);
4260 if ( 337.5 < angle || angle <= 22.5 )
4261 return; /* Near zero angle - no change! - At least not at this time */
4263 /* Handle special cases */
4264 switch (kernel->type) {
4265 /* These built-in kernels are cylindrical kernels, rotating is useless */
4266 case GaussianKernel:
4271 case LaplacianKernel:
4272 case ChebyshevKernel:
4273 case ManhattanKernel:
4274 case EuclideanKernel:
4277 /* These may be rotatable at non-90 angles in the future */
4278 /* but simply rotating them in multiples of 90 degrees is useless */
4285 /* These only allows a +/-90 degree rotation (by transpose) */
4286 /* A 180 degree rotation is useless */
4288 if ( 135.0 < angle && angle <= 225.0 )
4290 if ( 225.0 < angle && angle <= 315.0 )
4297 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4298 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4300 if ( kernel->width == 3 && kernel->height == 3 )
4301 { /* Rotate a 3x3 square by 45 degree angle */
4302 double t = kernel->values[0];
4303 kernel->values[0] = kernel->values[3];
4304 kernel->values[3] = kernel->values[6];
4305 kernel->values[6] = kernel->values[7];
4306 kernel->values[7] = kernel->values[8];
4307 kernel->values[8] = kernel->values[5];
4308 kernel->values[5] = kernel->values[2];
4309 kernel->values[2] = kernel->values[1];
4310 kernel->values[1] = t;
4311 /* rotate non-centered origin */
4312 if ( kernel->x != 1 || kernel->y != 1 ) {
4314 x = (ssize_t) kernel->x-1;
4315 y = (ssize_t) kernel->y-1;
4316 if ( x == y ) x = 0;
4317 else if ( x == 0 ) x = -y;
4318 else if ( x == -y ) y = 0;
4319 else if ( y == 0 ) y = x;
4320 kernel->x = (ssize_t) x+1;
4321 kernel->y = (ssize_t) y+1;
4323 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4324 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4327 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4329 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4331 if ( kernel->width == 1 || kernel->height == 1 )
4332 { /* Do a transpose of a 1 dimensional kernel,
4333 ** which results in a fast 90 degree rotation of some type.
4337 t = (ssize_t) kernel->width;
4338 kernel->width = kernel->height;
4339 kernel->height = (size_t) t;
4341 kernel->x = kernel->y;
4343 if ( kernel->width == 1 ) {
4344 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4345 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4347 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4348 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4351 else if ( kernel->width == kernel->height )
4352 { /* Rotate a square array of values by 90 degrees */
4356 register MagickRealType
4360 for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--)
4361 for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--)
4362 { t = k[i+j*kernel->width];
4363 k[i+j*kernel->width] = k[j+x*kernel->width];
4364 k[j+x*kernel->width] = k[x+y*kernel->width];
4365 k[x+y*kernel->width] = k[y+i*kernel->width];
4366 k[y+i*kernel->width] = t;
4369 /* rotate the origin - relative to center of array */
4370 { register ssize_t x,y;
4371 x = (ssize_t) (kernel->x*2-kernel->width+1);
4372 y = (ssize_t) (kernel->y*2-kernel->height+1);
4373 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4374 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4376 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4377 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4380 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4382 if ( 135.0 < angle && angle <= 225.0 )
4384 /* For a 180 degree rotation - also know as a reflection
4385 * This is actually a very very common operation!
4386 * Basically all that is needed is a reversal of the kernel data!
4387 * And a reflection of the origon
4392 register MagickRealType
4400 j=(ssize_t) (kernel->width*kernel->height-1);
4401 for (i=0; i < j; i++, j--)
4402 t=k[i], k[i]=k[j], k[j]=t;
4404 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4405 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4406 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4407 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4409 /* At this point angle should at least between -45 (315) and +45 degrees
4410 * In the future some form of non-orthogonal angled rotates could be
4411 * performed here, posibily with a linear kernel restriction.
4418 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4422 % S c a l e G e o m e t r y K e r n e l I n f o %
4426 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4428 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4429 % provided as a "-set option:convolve:scale {geometry}" user setting,
4430 % and modifies the kernel according to the parsed arguments of that setting.
4432 % The first argument (and any normalization flags) are passed to
4433 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4434 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4435 % into the scaled/normalized kernel.
4437 % The format of the ScaleGeometryKernelInfo method is:
4439 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4440 % const double scaling_factor,const MagickStatusType normalize_flags)
4442 % A description of each parameter follows:
4444 % o kernel: the Morphology/Convolution kernel to modify
4447 % The geometry string to parse, typically from the user provided
4448 % "-set option:convolve:scale {geometry}" setting.
4451 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4452 const char *geometry)
4460 SetGeometryInfo(&args);
4461 flags = ParseGeometry(geometry, &args);
4464 /* For Debugging Geometry Input */
4465 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4466 flags, args.rho, args.sigma, args.xi, args.psi );
4469 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4470 args.rho *= 0.01, args.sigma *= 0.01;
4472 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4474 if ( (flags & SigmaValue) == 0 )
4477 /* Scale/Normalize the input kernel */
4478 ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4480 /* Add Unity Kernel, for blending with original */
4481 if ( (flags & SigmaValue) != 0 )
4482 UnityAddKernelInfo(kernel, args.sigma);
4487 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4491 % S c a l e K e r n e l I n f o %
4495 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4497 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4498 % without normalization of the sum of the kernel values (as per given flags).
4500 % By default (no flags given) the values within the kernel is scaled
4501 % directly using given scaling factor without change.
4503 % If either of the two 'normalize_flags' are given the kernel will first be
4504 % normalized and then further scaled by the scaling factor value given.
4506 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4507 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4508 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4509 % non-HDRI versions of IM this may cause images to have any negative results
4510 % clipped, unless some 'bias' is used.
4512 % More specifically. Kernels which only contain positive values (such as a
4513 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4514 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4516 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4517 % the kernel will be scaled by the absolute of the sum of kernel values, so
4518 % that it will generally fall within the +/- 1.0 range.
4520 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4521 % will be scaled by just the sum of the postive values, so that its output
4522 % range will again fall into the +/- 1.0 range.
4524 % For special kernels designed for locating shapes using 'Correlate', (often
4525 % only containing +1 and -1 values, representing foreground/brackground
4526 % matching) a special normalization method is provided to scale the positive
4527 % values separately to those of the negative values, so the kernel will be
4528 % forced to become a zero-sum kernel better suited to such searches.
4530 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4531 % attributes within the kernel structure have been correctly set during the
4534 % NOTE: The values used for 'normalize_flags' have been selected specifically
4535 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4536 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4538 % The format of the ScaleKernelInfo method is:
4540 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4541 % const MagickStatusType normalize_flags )
4543 % A description of each parameter follows:
4545 % o kernel: the Morphology/Convolution kernel
4548 % multiply all values (after normalization) by this factor if not
4549 % zero. If the kernel is normalized regardless of any flags.
4551 % o normalize_flags:
4552 % GeometryFlags defining normalization method to use.
4553 % specifically: NormalizeValue, CorrelateNormalizeValue,
4554 % and/or PercentValue
4557 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4558 const double scaling_factor,const GeometryFlags normalize_flags)
4567 /* do the other kernels in a multi-kernel list first */
4568 if ( kernel->next != (KernelInfo *) NULL)
4569 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4571 /* Normalization of Kernel */
4573 if ( (normalize_flags&NormalizeValue) != 0 ) {
4574 if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4575 /* non-zero-summing kernel (generally positive) */
4576 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4578 /* zero-summing kernel */
4579 pos_scale = kernel->positive_range;
4581 /* Force kernel into a normalized zero-summing kernel */
4582 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4583 pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4584 ? kernel->positive_range : 1.0;
4585 neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4586 ? -kernel->negative_range : 1.0;
4589 neg_scale = pos_scale;
4591 /* finialize scaling_factor for positive and negative components */
4592 pos_scale = scaling_factor/pos_scale;
4593 neg_scale = scaling_factor/neg_scale;
4595 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4596 if ( ! IfNaN(kernel->values[i]) )
4597 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4599 /* convolution output range */
4600 kernel->positive_range *= pos_scale;
4601 kernel->negative_range *= neg_scale;
4602 /* maximum and minimum values in kernel */
4603 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4604 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4606 /* swap kernel settings if user's scaling factor is negative */
4607 if ( scaling_factor < MagickEpsilon ) {
4609 t = kernel->positive_range;
4610 kernel->positive_range = kernel->negative_range;
4611 kernel->negative_range = t;
4612 t = kernel->maximum;
4613 kernel->maximum = kernel->minimum;
4614 kernel->minimum = 1;
4621 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4625 % S h o w K e r n e l I n f o %
4629 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4631 % ShowKernelInfo() outputs the details of the given kernel defination to
4632 % standard error, generally due to a users 'showkernel' option request.
4634 % The format of the ShowKernel method is:
4636 % void ShowKernelInfo(const KernelInfo *kernel)
4638 % A description of each parameter follows:
4640 % o kernel: the Morphology/Convolution kernel
4643 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4651 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4653 (void) FormatLocaleFile(stderr, "Kernel");
4654 if ( kernel->next != (KernelInfo *) NULL )
4655 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4656 (void) FormatLocaleFile(stderr, " \"%s",
4657 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4658 if ( fabs(k->angle) >= MagickEpsilon )
4659 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4660 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4661 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4662 (void) FormatLocaleFile(stderr,
4663 " with values from %.*lg to %.*lg\n",
4664 GetMagickPrecision(), k->minimum,
4665 GetMagickPrecision(), k->maximum);
4666 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4667 GetMagickPrecision(), k->negative_range,
4668 GetMagickPrecision(), k->positive_range);
4669 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4670 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4671 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4672 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4674 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4675 GetMagickPrecision(), k->positive_range+k->negative_range);
4676 for (i=v=0; v < k->height; v++) {
4677 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4678 for (u=0; u < k->width; u++, i++)
4679 if ( IfNaN(k->values[i]) )
4680 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4682 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4683 GetMagickPrecision(), (double) k->values[i]);
4684 (void) FormatLocaleFile(stderr,"\n");
4690 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4694 % U n i t y A d d K e r n a l I n f o %
4698 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4700 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4701 % to the given pre-scaled and normalized Kernel. This in effect adds that
4702 % amount of the original image into the resulting convolution kernel. This
4703 % value is usually provided by the user as a percentage value in the
4704 % 'convolve:scale' setting.
4706 % The resulting effect is to convert the defined kernels into blended
4707 % soft-blurs, unsharp kernels or into sharpening kernels.
4709 % The format of the UnityAdditionKernelInfo method is:
4711 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4713 % A description of each parameter follows:
4715 % o kernel: the Morphology/Convolution kernel
4718 % scaling factor for the unity kernel to be added to
4722 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4725 /* do the other kernels in a multi-kernel list first */
4726 if ( kernel->next != (KernelInfo *) NULL)
4727 UnityAddKernelInfo(kernel->next, scale);
4729 /* Add the scaled unity kernel to the existing kernel */
4730 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4731 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4737 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4741 % Z e r o K e r n e l N a n s %
4745 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4747 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4748 % the kernel with a zero value. This is typically done when the kernel will
4749 % be used in special hardware (GPU) convolution processors, to simply
4752 % The format of the ZeroKernelNans method is:
4754 % void ZeroKernelNans (KernelInfo *kernel)
4756 % A description of each parameter follows:
4758 % o kernel: the Morphology/Convolution kernel
4761 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4766 /* do the other kernels in a multi-kernel list first */
4767 if ( kernel->next != (KernelInfo *) NULL)
4768 ZeroKernelNans(kernel->next);
4770 for (i=0; i < (kernel->width*kernel->height); i++)
4771 if ( IfNaN(kernel->values[i]) )
4772 kernel->values[i] = 0.0;