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-2012 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % http://www.imagemagick.org/script/license.php %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 % Morpology is the the application of various kernels, of any size and even
37 % shape, to a image in various ways (typically binary, but not always).
39 % Convolution (weighted sum or average) is just one specific type of
40 % morphology. Just one that is very common for image bluring and sharpening
41 % effects. Not only 2D Gaussian blurring, but also 2-pass 1D Blurring.
43 % This module provides not only a general morphology function, and the ability
44 % to apply more advanced or iterative morphologies, but also functions for the
45 % generation of many different types of kernel arrays from user supplied
46 % arguments. Prehaps even the generation of a kernel from a small image.
52 #include "MagickCore/studio.h"
53 #include "MagickCore/artifact.h"
54 #include "MagickCore/cache-view.h"
55 #include "MagickCore/color-private.h"
56 #include "MagickCore/enhance.h"
57 #include "MagickCore/exception.h"
58 #include "MagickCore/exception-private.h"
59 #include "MagickCore/gem.h"
60 #include "MagickCore/gem-private.h"
61 #include "MagickCore/hashmap.h"
62 #include "MagickCore/image.h"
63 #include "MagickCore/image-private.h"
64 #include "MagickCore/list.h"
65 #include "MagickCore/magick.h"
66 #include "MagickCore/memory_.h"
67 #include "MagickCore/memory-private.h"
68 #include "MagickCore/monitor-private.h"
69 #include "MagickCore/morphology.h"
70 #include "MagickCore/morphology-private.h"
71 #include "MagickCore/option.h"
72 #include "MagickCore/pixel-accessor.h"
73 #include "MagickCore/prepress.h"
74 #include "MagickCore/quantize.h"
75 #include "MagickCore/resource_.h"
76 #include "MagickCore/registry.h"
77 #include "MagickCore/semaphore.h"
78 #include "MagickCore/splay-tree.h"
79 #include "MagickCore/statistic.h"
80 #include "MagickCore/string_.h"
81 #include "MagickCore/string-private.h"
82 #include "MagickCore/thread-private.h"
83 #include "MagickCore/token.h"
84 #include "MagickCore/utility.h"
85 #include "MagickCore/utility-private.h"
89 ** The following test is for special floating point numbers of value NaN (not
90 ** a number), that may be used within a Kernel Definition. NaN's are defined
91 ** as part of the IEEE standard for floating point number representation.
93 ** These are used as a Kernel value to mean that this kernel position is not
94 ** part of the kernel neighbourhood for convolution or morphology processing,
95 ** and thus should be ignored. This allows the use of 'shaped' kernels.
97 ** The special property that two NaN's are never equal, even if they are from
98 ** the same variable allow you to test if a value is special NaN value.
100 ** This macro IsNaN() is thus is only true if the value given is NaN.
102 #define IsNan(a) ((a)!=(a))
105 Other global definitions used by module.
107 static inline double MagickMin(const double x,const double y)
109 return( x < y ? x : y);
111 static inline double MagickMax(const double x,const double y)
113 return( x > y ? x : y);
115 #define Minimize(assign,value) assign=MagickMin(assign,value)
116 #define Maximize(assign,value) assign=MagickMax(assign,value)
118 /* Integer Factorial Function - for a Binomial kernel */
120 static inline size_t fact(size_t n)
123 for(f=1, l=2; l <= n; f=f*l, l++);
126 #elif 1 /* glibc floating point alternatives */
127 #define fact(n) ((size_t)tgamma((double)n+1))
129 #define fact(n) ((size_t)lgamma((double)n+1))
133 /* Currently these are only internal to this module */
135 CalcKernelMetaData(KernelInfo *),
136 ExpandMirrorKernelInfo(KernelInfo *),
137 ExpandRotateKernelInfo(KernelInfo *, const double),
138 RotateKernelInfo(KernelInfo *, double);
141 /* Quick function to find last kernel in a kernel list */
142 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
144 while (kernel->next != (KernelInfo *) NULL)
145 kernel = kernel->next;
150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
154 % A c q u i r e K e r n e l I n f o %
158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
160 % AcquireKernelInfo() takes the given string (generally supplied by the
161 % user) and converts it into a Morphology/Convolution Kernel. This allows
162 % users to specify a kernel from a number of pre-defined kernels, or to fully
163 % specify their own kernel for a specific Convolution or Morphology
166 % The kernel so generated can be any rectangular array of floating point
167 % values (doubles) with the 'control point' or 'pixel being affected'
168 % anywhere within that array of values.
170 % Previously IM was restricted to a square of odd size using the exact
171 % center as origin, this is no longer the case, and any rectangular kernel
172 % with any value being declared the origin. This in turn allows the use of
173 % highly asymmetrical kernels.
175 % The floating point values in the kernel can also include a special value
176 % known as 'nan' or 'not a number' to indicate that this value is not part
177 % of the kernel array. This allows you to shaped the kernel within its
178 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
179 % shape. However at least one non-nan value must be provided for correct
180 % working of a kernel.
182 % The returned kernel should be freed using the DestroyKernelInfo() when you
183 % are finished with it. Do not free this memory yourself.
185 % Input kernel defintion strings can consist of any of three types.
188 % Select from one of the built in kernels, using the name and
189 % geometry arguments supplied. See AcquireKernelBuiltIn()
191 % "WxH[+X+Y][@><]:num, num, num ..."
192 % a kernel of size W by H, with W*H floating point numbers following.
193 % the 'center' can be optionally be defined at +X+Y (such that +0+0
194 % is top left corner). If not defined the pixel in the center, for
195 % odd sizes, or to the immediate top or left of center for even sizes
196 % is automatically selected.
198 % "num, num, num, num, ..."
199 % list of floating point numbers defining an 'old style' odd sized
200 % square kernel. At least 9 values should be provided for a 3x3
201 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
202 % Values can be space or comma separated. This is not recommended.
204 % You can define a 'list of kernels' which can be used by some morphology
205 % operators A list is defined as a semi-colon separated list kernels.
207 % " kernel ; kernel ; kernel ; "
209 % Any extra ';' characters, at start, end or between kernel defintions are
212 % The special flags will expand a single kernel, into a list of rotated
213 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
214 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
215 % The '<' also exands using 90-degree rotates, but giving a 180-degree
216 % reflected kernel before the +/- 90-degree rotations, which can be important
217 % for Thinning operations.
219 % Note that 'name' kernels will start with an alphabetic character while the
220 % new kernel specification has a ':' character in its specification string.
221 % If neither is the case, it is assumed an old style of a simple list of
222 % numbers generating a odd-sized square kernel has been given.
224 % The format of the AcquireKernal method is:
226 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
228 % A description of each parameter follows:
230 % o kernel_string: the Morphology/Convolution kernel wanted.
234 /* This was separated so that it could be used as a separate
235 ** array input handling function, such as for -color-matrix
237 static KernelInfo *ParseKernelArray(const char *kernel_string)
243 token[MaxTextExtent];
253 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
261 kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
262 if (kernel == (KernelInfo *)NULL)
264 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
265 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
266 kernel->negative_range = kernel->positive_range = 0.0;
267 kernel->type = UserDefinedKernel;
268 kernel->next = (KernelInfo *) NULL;
269 kernel->signature = MagickSignature;
270 if (kernel_string == (const char *) NULL)
273 /* find end of this specific kernel definition string */
274 end = strchr(kernel_string, ';');
275 if ( end == (char *) NULL )
276 end = strchr(kernel_string, '\0');
278 /* clear flags - for Expanding kernel lists thorugh rotations */
281 /* Has a ':' in argument - New user kernel specification
282 FUTURE: this split on ':' could be done by StringToken()
284 p = strchr(kernel_string, ':');
285 if ( p != (char *) NULL && p < end)
287 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
288 memcpy(token, kernel_string, (size_t) (p-kernel_string));
289 token[p-kernel_string] = '\0';
290 SetGeometryInfo(&args);
291 flags = ParseGeometry(token, &args);
293 /* Size handling and checks of geometry settings */
294 if ( (flags & WidthValue) == 0 ) /* if no width then */
295 args.rho = args.sigma; /* then width = height */
296 if ( args.rho < 1.0 ) /* if width too small */
297 args.rho = 1.0; /* then width = 1 */
298 if ( args.sigma < 1.0 ) /* if height too small */
299 args.sigma = args.rho; /* then height = width */
300 kernel->width = (size_t)args.rho;
301 kernel->height = (size_t)args.sigma;
303 /* Offset Handling and Checks */
304 if ( args.xi < 0.0 || args.psi < 0.0 )
305 return(DestroyKernelInfo(kernel));
306 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
307 : (ssize_t) (kernel->width-1)/2;
308 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
309 : (ssize_t) (kernel->height-1)/2;
310 if ( kernel->x >= (ssize_t) kernel->width ||
311 kernel->y >= (ssize_t) kernel->height )
312 return(DestroyKernelInfo(kernel));
314 p++; /* advance beyond the ':' */
317 { /* ELSE - Old old specification, forming odd-square kernel */
318 /* count up number of values given */
319 p=(const char *) kernel_string;
320 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
321 p++; /* ignore "'" chars for convolve filter usage - Cristy */
322 for (i=0; p < end; i++)
324 GetMagickToken(p,&p,token);
326 GetMagickToken(p,&p,token);
328 /* set the size of the kernel - old sized square */
329 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
330 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
331 p=(const char *) kernel_string;
332 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
333 p++; /* ignore "'" chars for convolve filter usage - Cristy */
336 /* Read in the kernel values from rest of input string argument */
337 kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
338 kernel->width,kernel->height*sizeof(*kernel->values)));
339 if (kernel->values == (MagickRealType *) NULL)
340 return(DestroyKernelInfo(kernel));
341 kernel->minimum = +MagickHuge;
342 kernel->maximum = -MagickHuge;
343 kernel->negative_range = kernel->positive_range = 0.0;
344 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
346 GetMagickToken(p,&p,token);
348 GetMagickToken(p,&p,token);
349 if ( LocaleCompare("nan",token) == 0
350 || LocaleCompare("-",token) == 0 ) {
351 kernel->values[i] = nan; /* this value is not part of neighbourhood */
354 kernel->values[i] = StringToDouble(token,(char **) NULL);
355 ( kernel->values[i] < 0)
356 ? ( kernel->negative_range += kernel->values[i] )
357 : ( kernel->positive_range += kernel->values[i] );
358 Minimize(kernel->minimum, kernel->values[i]);
359 Maximize(kernel->maximum, kernel->values[i]);
363 /* sanity check -- no more values in kernel definition */
364 GetMagickToken(p,&p,token);
365 if ( *token != '\0' && *token != ';' && *token != '\'' )
366 return(DestroyKernelInfo(kernel));
369 /* this was the old method of handling a incomplete kernel */
370 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
371 Minimize(kernel->minimum, kernel->values[i]);
372 Maximize(kernel->maximum, kernel->values[i]);
373 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
374 kernel->values[i]=0.0;
377 /* Number of values for kernel was not enough - Report Error */
378 if ( i < (ssize_t) (kernel->width*kernel->height) )
379 return(DestroyKernelInfo(kernel));
382 /* check that we recieved at least one real (non-nan) value! */
383 if ( kernel->minimum == MagickHuge )
384 return(DestroyKernelInfo(kernel));
386 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
387 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
388 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
389 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
390 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
391 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
396 static KernelInfo *ParseKernelName(const char *kernel_string)
399 token[MaxTextExtent];
417 /* Parse special 'named' kernel */
418 GetMagickToken(kernel_string,&p,token);
419 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
420 if ( type < 0 || type == UserDefinedKernel )
421 return((KernelInfo *)NULL); /* not a valid named kernel */
423 while (((isspace((int) ((unsigned char) *p)) != 0) ||
424 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
427 end = strchr(p, ';'); /* end of this kernel defintion */
428 if ( end == (char *) NULL )
429 end = strchr(p, '\0');
431 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
432 memcpy(token, p, (size_t) (end-p));
434 SetGeometryInfo(&args);
435 flags = ParseGeometry(token, &args);
438 /* For Debugging Geometry Input */
439 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
440 flags, args.rho, args.sigma, args.xi, args.psi );
443 /* special handling of missing values in input string */
445 /* Shape Kernel Defaults */
447 if ( (flags & WidthValue) == 0 )
448 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
456 if ( (flags & HeightValue) == 0 )
457 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
460 if ( (flags & XValue) == 0 )
461 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
463 case RectangleKernel: /* Rectangle - set size defaults */
464 if ( (flags & WidthValue) == 0 ) /* if no width then */
465 args.rho = args.sigma; /* then width = height */
466 if ( args.rho < 1.0 ) /* if width too small */
467 args.rho = 3; /* then width = 3 */
468 if ( args.sigma < 1.0 ) /* if height too small */
469 args.sigma = args.rho; /* then height = width */
470 if ( (flags & XValue) == 0 ) /* center offset if not defined */
471 args.xi = (double)(((ssize_t)args.rho-1)/2);
472 if ( (flags & YValue) == 0 )
473 args.psi = (double)(((ssize_t)args.sigma-1)/2);
475 /* Distance Kernel Defaults */
476 case ChebyshevKernel:
477 case ManhattanKernel:
478 case OctagonalKernel:
479 case EuclideanKernel:
480 if ( (flags & HeightValue) == 0 ) /* no distance scale */
481 args.sigma = 100.0; /* default distance scaling */
482 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
483 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
484 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
485 args.sigma *= QuantumRange/100.0; /* percentage of color range */
491 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
492 if ( kernel == (KernelInfo *) NULL )
495 /* global expand to rotated kernel list - only for single kernels */
496 if ( kernel->next == (KernelInfo *) NULL ) {
497 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
498 ExpandRotateKernelInfo(kernel, 45.0);
499 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
500 ExpandRotateKernelInfo(kernel, 90.0);
501 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
502 ExpandMirrorKernelInfo(kernel);
508 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
516 token[MaxTextExtent];
524 if (kernel_string == (const char *) NULL)
525 return(ParseKernelArray(kernel_string));
530 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
532 /* ignore extra or multiple ';' kernel separators */
533 if ( *token != ';' ) {
535 /* tokens starting with alpha is a Named kernel */
536 if (isalpha((int) *token) != 0)
537 new_kernel = ParseKernelName(p);
538 else /* otherwise a user defined kernel array */
539 new_kernel = ParseKernelArray(p);
541 /* Error handling -- this is not proper error handling! */
542 if ( new_kernel == (KernelInfo *) NULL ) {
543 (void) FormatLocaleFile(stderr, "Failed to parse kernel number #%.20g\n",
544 (double) kernel_number);
545 if ( kernel != (KernelInfo *) NULL )
546 kernel=DestroyKernelInfo(kernel);
547 return((KernelInfo *) NULL);
550 /* initialise or append the kernel list */
551 if ( kernel == (KernelInfo *) NULL )
554 LastKernelInfo(kernel)->next = new_kernel;
557 /* look for the next kernel in list */
559 if ( p == (char *) NULL )
569 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
573 % A c q u i r e K e r n e l B u i l t I n %
577 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
579 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
580 % kernels used for special purposes such as gaussian blurring, skeleton
581 % pruning, and edge distance determination.
583 % They take a KernelType, and a set of geometry style arguments, which were
584 % typically decoded from a user supplied string, or from a more complex
585 % Morphology Method that was requested.
587 % The format of the AcquireKernalBuiltIn method is:
589 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
590 % const GeometryInfo args)
592 % A description of each parameter follows:
594 % o type: the pre-defined type of kernel wanted
596 % o args: arguments defining or modifying the kernel
598 % Convolution Kernels
601 % The a No-Op or Scaling single element kernel.
603 % Gaussian:{radius},{sigma}
604 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
605 % The sigma for the curve is required. The resulting kernel is
608 % If 'sigma' is zero, you get a single pixel on a field of zeros.
610 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
611 % the final size of the resulting kernel to a square 2*radius+1 in size.
612 % The radius should be at least 2 times that of the sigma value, or
613 % sever clipping and aliasing may result. If not given or set to 0 the
614 % radius will be determined so as to produce the best minimal error
615 % result, which is usally much larger than is normally needed.
617 % LoG:{radius},{sigma}
618 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
619 % The supposed ideal edge detection, zero-summing kernel.
621 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
622 % approx 1.6 (according to wikipedia).
624 % DoG:{radius},{sigma1},{sigma2}
625 % "Difference of Gaussians" Kernel.
626 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
627 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
628 % The result is a zero-summing kernel.
630 % Blur:{radius},{sigma}[,{angle}]
631 % Generates a 1 dimensional or linear gaussian blur, at the angle given
632 % (current restricted to orthogonal angles). If a 'radius' is given the
633 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
634 % by a 90 degree angle.
636 % If 'sigma' is zero, you get a single pixel on a field of zeros.
638 % Note that two convolutions with two "Blur" kernels perpendicular to
639 % each other, is equivalent to a far larger "Gaussian" kernel with the
640 % same sigma value, However it is much faster to apply. This is how the
641 % "-blur" operator actually works.
643 % Comet:{width},{sigma},{angle}
644 % Blur in one direction only, much like how a bright object leaves
645 % a comet like trail. The Kernel is actually half a gaussian curve,
646 % Adding two such blurs in opposite directions produces a Blur Kernel.
647 % Angle can be rotated in multiples of 90 degrees.
649 % Note that the first argument is the width of the kernel and not the
650 % radius of the kernel.
652 % Binomial:[{radius}]
653 % Generate a discrete kernel using a 2 dimentional Pascel's Triangle
654 % of values. Used for special forma of image filters.
656 % # Still to be implemented...
660 % # Set kernel values using a resize filter, and given scale (sigma)
661 % # Cylindrical or Linear. Is this possible with an image?
664 % Named Constant Convolution Kernels
666 % All these are unscaled, zero-summing kernels by default. As such for
667 % non-HDRI version of ImageMagick some form of normalization, user scaling,
668 % and biasing the results is recommended, to prevent the resulting image
671 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
672 % 45 degrees to generate the 8 angled varients of each of the kernels.
675 % Discrete Lapacian Kernels, (without normalization)
676 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
677 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
678 % Type 2 : 3x3 with center:4 edge:1 corner:-2
679 % Type 3 : 3x3 with center:4 edge:-2 corner:1
680 % Type 5 : 5x5 laplacian
681 % Type 7 : 7x7 laplacian
682 % Type 15 : 5x5 LoG (sigma approx 1.4)
683 % Type 19 : 9x9 LoG (sigma approx 1.4)
686 % Sobel 'Edge' convolution kernel (3x3)
692 % Roberts convolution kernel (3x3)
698 % Prewitt Edge convolution kernel (3x3)
704 % Prewitt's "Compass" convolution kernel (3x3)
710 % Kirsch's "Compass" convolution kernel (3x3)
716 % Frei-Chen Edge Detector is based on a kernel that is similar to
717 % the Sobel Kernel, but is designed to be isotropic. That is it takes
718 % into account the distance of the diagonal in the kernel.
721 % | sqrt(2), 0, -sqrt(2) |
724 % FreiChen:{type},{angle}
726 % Frei-Chen Pre-weighted kernels...
728 % Type 0: default un-nomalized version shown above.
730 % Type 1: Orthogonal Kernel (same as type 11 below)
732 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
735 % Type 2: Diagonal form of Kernel...
737 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
740 % However this kernel is als at the heart of the FreiChen Edge Detection
741 % Process which uses a set of 9 specially weighted kernel. These 9
742 % kernels not be normalized, but directly applied to the image. The
743 % results is then added together, to produce the intensity of an edge in
744 % a specific direction. The square root of the pixel value can then be
745 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
746 % from each other, both the direction and the strength of the edge can be
749 % Type 10: All 9 of the following pre-weighted kernels...
751 % Type 11: | 1, 0, -1 |
752 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
755 % Type 12: | 1, sqrt(2), 1 |
756 % | 0, 0, 0 | / 2*sqrt(2)
759 % Type 13: | sqrt(2), -1, 0 |
760 % | -1, 0, 1 | / 2*sqrt(2)
763 % Type 14: | 0, 1, -sqrt(2) |
764 % | -1, 0, 1 | / 2*sqrt(2)
767 % Type 15: | 0, -1, 0 |
771 % Type 16: | 1, 0, -1 |
775 % Type 17: | 1, -2, 1 |
779 % Type 18: | -2, 1, -2 |
783 % Type 19: | 1, 1, 1 |
787 % The first 4 are for edge detection, the next 4 are for line detection
788 % and the last is to add a average component to the results.
790 % Using a special type of '-1' will return all 9 pre-weighted kernels
791 % as a multi-kernel list, so that you can use them directly (without
792 % normalization) with the special "-set option:morphology:compose Plus"
793 % setting to apply the full FreiChen Edge Detection Technique.
795 % If 'type' is large it will be taken to be an actual rotation angle for
796 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
797 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
799 % WARNING: The above was layed out as per
800 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
801 % But rotated 90 degrees so direction is from left rather than the top.
802 % I have yet to find any secondary confirmation of the above. The only
803 % other source found was actual source code at
804 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
805 % Neigher paper defineds the kernels in a way that looks locical or
806 % correct when taken as a whole.
810 % Diamond:[{radius}[,{scale}]]
811 % Generate a diamond shaped kernel with given radius to the points.
812 % Kernel size will again be radius*2+1 square and defaults to radius 1,
813 % generating a 3x3 kernel that is slightly larger than a square.
815 % Square:[{radius}[,{scale}]]
816 % Generate a square shaped kernel of size radius*2+1, and defaulting
817 % to a 3x3 (radius 1).
819 % Octagon:[{radius}[,{scale}]]
820 % Generate octagonal shaped kernel of given radius and constant scale.
821 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
822 % in "Diamond" kernel.
824 % Disk:[{radius}[,{scale}]]
825 % Generate a binary disk, thresholded at the radius given, the radius
826 % may be a float-point value. Final Kernel size is floor(radius)*2+1
827 % square. A radius of 5.3 is the default.
829 % NOTE: That a low radii Disk kernels produce the same results as
830 % many of the previously defined kernels, but differ greatly at larger
831 % radii. Here is a table of equivalences...
832 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
833 % "Disk:1.5" => "Square"
834 % "Disk:2" => "Diamond:2"
835 % "Disk:2.5" => "Octagon"
836 % "Disk:2.9" => "Square:2"
837 % "Disk:3.5" => "Octagon:3"
838 % "Disk:4.5" => "Octagon:4"
839 % "Disk:5.4" => "Octagon:5"
840 % "Disk:6.4" => "Octagon:6"
841 % All other Disk shapes are unique to this kernel, but because a "Disk"
842 % is more circular when using a larger radius, using a larger radius is
843 % preferred over iterating the morphological operation.
845 % Rectangle:{geometry}
846 % Simply generate a rectangle of 1's with the size given. You can also
847 % specify the location of the 'control point', otherwise the closest
848 % pixel to the center of the rectangle is selected.
850 % Properly centered and odd sized rectangles work the best.
852 % Symbol Dilation Kernels
854 % These kernel is not a good general morphological kernel, but is used
855 % more for highlighting and marking any single pixels in an image using,
856 % a "Dilate" method as appropriate.
858 % For the same reasons iterating these kernels does not produce the
859 % same result as using a larger radius for the symbol.
861 % Plus:[{radius}[,{scale}]]
862 % Cross:[{radius}[,{scale}]]
863 % Generate a kernel in the shape of a 'plus' or a 'cross' with
864 % a each arm the length of the given radius (default 2).
866 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
868 % Ring:{radius1},{radius2}[,{scale}]
869 % A ring of the values given that falls between the two radii.
870 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
871 % This is the 'edge' pixels of the default "Disk" kernel,
872 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
874 % Hit and Miss Kernels
876 % Peak:radius1,radius2
877 % Find any peak larger than the pixels the fall between the two radii.
878 % The default ring of pixels is as per "Ring".
880 % Find flat orthogonal edges of a binary shape
882 % Find 90 degree corners of a binary shape
884 % A special kernel to thin the 'outside' of diagonals
886 % Find end points of lines (for pruning a skeletion)
887 % Two types of lines ends (default to both) can be searched for
888 % Type 0: All line ends
889 % Type 1: single kernel for 4-conneected line ends
890 % Type 2: single kernel for simple line ends
892 % Find three line junctions (within a skeletion)
893 % Type 0: all line junctions
894 % Type 1: Y Junction kernel
895 % Type 2: Diagonal T Junction kernel
896 % Type 3: Orthogonal T Junction kernel
897 % Type 4: Diagonal X Junction kernel
898 % Type 5: Orthogonal + Junction kernel
900 % Find single pixel ridges or thin lines
901 % Type 1: Fine single pixel thick lines and ridges
902 % Type 2: Find two pixel thick lines and ridges
904 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
906 % Traditional skeleton generating kernels.
907 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
908 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
909 % Type 3: Thinning skeleton based on a ressearch paper by
910 % Dan S. Bloomberg (Default Type)
912 % A huge variety of Thinning Kernels designed to preserve conectivity.
913 % many other kernel sets use these kernels as source definitions.
914 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
915 % the super and sub notations used in the source research paper.
917 % Distance Measuring Kernels
919 % Different types of distance measuring methods, which are used with the
920 % a 'Distance' morphology method for generating a gradient based on
921 % distance from an edge of a binary shape, though there is a technique
922 % for handling a anti-aliased shape.
924 % See the 'Distance' Morphological Method, for information of how it is
927 % Chebyshev:[{radius}][x{scale}[%!]]
928 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
929 % is a value of one to any neighbour, orthogonal or diagonal. One why
930 % of thinking of it is the number of squares a 'King' or 'Queen' in
931 % chess needs to traverse reach any other position on a chess board.
932 % It results in a 'square' like distance function, but one where
933 % diagonals are given a value that is closer than expected.
935 % Manhattan:[{radius}][x{scale}[%!]]
936 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
937 % Cab distance metric), it is the distance needed when you can only
938 % travel in horizontal or vertical directions only. It is the
939 % distance a 'Rook' in chess would have to travel, and results in a
940 % diamond like distances, where diagonals are further than expected.
942 % Octagonal:[{radius}][x{scale}[%!]]
943 % An interleving of Manhatten and Chebyshev metrics producing an
944 % increasing octagonally shaped distance. Distances matches those of
945 % the "Octagon" shaped kernel of the same radius. The minimum radius
946 % and default is 2, producing a 5x5 kernel.
948 % Euclidean:[{radius}][x{scale}[%!]]
949 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
950 % However by default the kernel size only has a radius of 1, which
951 % limits the distance to 'Knight' like moves, with only orthogonal and
952 % diagonal measurements being correct. As such for the default kernel
953 % you will get octagonal like distance function.
955 % However using a larger radius such as "Euclidean:4" you will get a
956 % much smoother distance gradient from the edge of the shape. Especially
957 % if the image is pre-processed to include any anti-aliasing pixels.
958 % Of course a larger kernel is slower to use, and not always needed.
960 % The first three Distance Measuring Kernels will only generate distances
961 % of exact multiples of {scale} in binary images. As such you can use a
962 % scale of 1 without loosing any information. However you also need some
963 % scaling when handling non-binary anti-aliased shapes.
965 % The "Euclidean" Distance Kernel however does generate a non-integer
966 % fractional results, and as such scaling is vital even for binary shapes.
970 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
971 const GeometryInfo *args)
984 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
986 /* Generate a new empty kernel if needed */
987 kernel=(KernelInfo *) NULL;
989 case UndefinedKernel: /* These should not call this function */
990 case UserDefinedKernel:
991 assert("Should not call this function" != (char *)NULL);
993 case LaplacianKernel: /* Named Descrete Convolution Kernels */
994 case SobelKernel: /* these are defined using other kernels */
1000 case EdgesKernel: /* Hit and Miss kernels */
1002 case DiagonalsKernel:
1003 case LineEndsKernel:
1004 case LineJunctionsKernel:
1006 case ConvexHullKernel:
1007 case SkeletonKernel:
1009 break; /* A pre-generated kernel is not needed */
1011 /* set to 1 to do a compile-time check that we haven't missed anything */
1013 case GaussianKernel:
1018 case BinomialKernel:
1021 case RectangleKernel:
1028 case ChebyshevKernel:
1029 case ManhattanKernel:
1030 case OctangonalKernel:
1031 case EuclideanKernel:
1035 /* Generate the base Kernel Structure */
1036 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1037 if (kernel == (KernelInfo *) NULL)
1039 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1040 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1041 kernel->negative_range = kernel->positive_range = 0.0;
1042 kernel->type = type;
1043 kernel->next = (KernelInfo *) NULL;
1044 kernel->signature = MagickSignature;
1054 kernel->height = kernel->width = (size_t) 1;
1055 kernel->x = kernel->y = (ssize_t) 0;
1056 kernel->values=(MagickRealType *) MagickAssumeAligned(
1057 AcquireAlignedMemory(1,sizeof(*kernel->values)));
1058 if (kernel->values == (MagickRealType *) NULL)
1059 return(DestroyKernelInfo(kernel));
1060 kernel->maximum = kernel->values[0] = args->rho;
1064 case GaussianKernel:
1068 sigma = fabs(args->sigma),
1069 sigma2 = fabs(args->xi),
1072 if ( args->rho >= 1.0 )
1073 kernel->width = (size_t)args->rho*2+1;
1074 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1075 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1077 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1078 kernel->height = kernel->width;
1079 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1080 kernel->values=(MagickRealType *) MagickAssumeAligned(
1081 AcquireAlignedMemory(kernel->width,kernel->height*
1082 sizeof(*kernel->values)));
1083 if (kernel->values == (MagickRealType *) NULL)
1084 return(DestroyKernelInfo(kernel));
1086 /* WARNING: The following generates a 'sampled gaussian' kernel.
1087 * What we really want is a 'discrete gaussian' kernel.
1089 * How to do this is I don't know, but appears to be basied on the
1090 * Error Function 'erf()' (intergral of a gaussian)
1093 if ( type == GaussianKernel || type == DoGKernel )
1094 { /* Calculate a Gaussian, OR positive half of a DoG */
1095 if ( sigma > MagickEpsilon )
1096 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1097 B = (double) (1.0/(Magick2PI*sigma*sigma));
1098 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1099 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1100 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1102 else /* limiting case - a unity (normalized Dirac) kernel */
1103 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1104 kernel->width*kernel->height*sizeof(double));
1105 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1109 if ( type == DoGKernel )
1110 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1111 if ( sigma2 > MagickEpsilon )
1112 { sigma = sigma2; /* simplify loop expressions */
1113 A = 1.0/(2.0*sigma*sigma);
1114 B = (double) (1.0/(Magick2PI*sigma*sigma));
1115 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1116 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1117 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1119 else /* limiting case - a unity (normalized Dirac) kernel */
1120 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1123 if ( type == LoGKernel )
1124 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1125 if ( sigma > MagickEpsilon )
1126 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1127 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1128 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1129 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1130 { R = ((double)(u*u+v*v))*A;
1131 kernel->values[i] = (1-R)*exp(-R)*B;
1134 else /* special case - generate a unity kernel */
1135 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1136 kernel->width*kernel->height*sizeof(double));
1137 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1141 /* Note the above kernels may have been 'clipped' by a user defined
1142 ** radius, producing a smaller (darker) kernel. Also for very small
1143 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1144 ** producing a very bright kernel.
1146 ** Normalization will still be needed.
1149 /* Normalize the 2D Gaussian Kernel
1151 ** NB: a CorrelateNormalize performs a normal Normalize if
1152 ** there are no negative values.
1154 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1155 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1161 sigma = fabs(args->sigma),
1164 if ( args->rho >= 1.0 )
1165 kernel->width = (size_t)args->rho*2+1;
1167 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1169 kernel->x = (ssize_t) (kernel->width-1)/2;
1171 kernel->negative_range = kernel->positive_range = 0.0;
1172 kernel->values=(MagickRealType *) MagickAssumeAligned(
1173 AcquireAlignedMemory(kernel->width,kernel->height*
1174 sizeof(*kernel->values)));
1175 if (kernel->values == (MagickRealType *) NULL)
1176 return(DestroyKernelInfo(kernel));
1179 #define KernelRank 3
1180 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1181 ** It generates a gaussian 3 times the width, and compresses it into
1182 ** the expected range. This produces a closer normalization of the
1183 ** resulting kernel, especially for very low sigma values.
1184 ** As such while wierd it is prefered.
1186 ** I am told this method originally came from Photoshop.
1188 ** A properly normalized curve is generated (apart from edge clipping)
1189 ** even though we later normalize the result (for edge clipping)
1190 ** to allow the correct generation of a "Difference of Blurs".
1194 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1195 (void) ResetMagickMemory(kernel->values,0, (size_t)
1196 kernel->width*kernel->height*sizeof(double));
1197 /* Calculate a Positive 1D Gaussian */
1198 if ( sigma > MagickEpsilon )
1199 { sigma *= KernelRank; /* simplify loop expressions */
1200 alpha = 1.0/(2.0*sigma*sigma);
1201 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1202 for ( u=-v; u <= v; u++) {
1203 kernel->values[(u+v)/KernelRank] +=
1204 exp(-((double)(u*u))*alpha)*beta;
1207 else /* special case - generate a unity kernel */
1208 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1210 /* Direct calculation without curve averaging
1211 This is equivelent to a KernelRank of 1 */
1213 /* Calculate a Positive Gaussian */
1214 if ( sigma > MagickEpsilon )
1215 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1216 beta = 1.0/(MagickSQ2PI*sigma);
1217 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1218 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1220 else /* special case - generate a unity kernel */
1221 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1222 kernel->width*kernel->height*sizeof(double));
1223 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1226 /* Note the above kernel may have been 'clipped' by a user defined
1227 ** radius, producing a smaller (darker) kernel. Also for very small
1228 ** sigma's (> 0.1) the central value becomes larger than one, as a
1229 ** result of not generating a actual 'discrete' kernel, and thus
1230 ** producing a very bright 'impulse'.
1232 ** Becuase of these two factors Normalization is required!
1235 /* Normalize the 1D Gaussian Kernel
1237 ** NB: a CorrelateNormalize performs a normal Normalize if
1238 ** there are no negative values.
1240 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1241 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1243 /* rotate the 1D kernel by given angle */
1244 RotateKernelInfo(kernel, args->xi );
1249 sigma = fabs(args->sigma),
1252 if ( args->rho < 1.0 )
1253 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1255 kernel->width = (size_t)args->rho;
1256 kernel->x = kernel->y = 0;
1258 kernel->negative_range = kernel->positive_range = 0.0;
1259 kernel->values=(MagickRealType *) MagickAssumeAligned(
1260 AcquireAlignedMemory(kernel->width,kernel->height*
1261 sizeof(*kernel->values)));
1262 if (kernel->values == (MagickRealType *) NULL)
1263 return(DestroyKernelInfo(kernel));
1265 /* A comet blur is half a 1D gaussian curve, so that the object is
1266 ** blurred in one direction only. This may not be quite the right
1267 ** curve to use so may change in the future. The function must be
1268 ** normalised after generation, which also resolves any clipping.
1270 ** As we are normalizing and not subtracting gaussians,
1271 ** there is no need for a divisor in the gaussian formula
1273 ** It is less comples
1275 if ( sigma > MagickEpsilon )
1278 #define KernelRank 3
1279 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1280 (void) ResetMagickMemory(kernel->values,0, (size_t)
1281 kernel->width*sizeof(double));
1282 sigma *= KernelRank; /* simplify the loop expression */
1283 A = 1.0/(2.0*sigma*sigma);
1284 /* B = 1.0/(MagickSQ2PI*sigma); */
1285 for ( u=0; u < v; u++) {
1286 kernel->values[u/KernelRank] +=
1287 exp(-((double)(u*u))*A);
1288 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1290 for (i=0; i < (ssize_t) kernel->width; i++)
1291 kernel->positive_range += kernel->values[i];
1293 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1294 /* B = 1.0/(MagickSQ2PI*sigma); */
1295 for ( i=0; i < (ssize_t) kernel->width; i++)
1296 kernel->positive_range +=
1297 kernel->values[i] = exp(-((double)(i*i))*A);
1298 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1301 else /* special case - generate a unity kernel */
1302 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1303 kernel->width*kernel->height*sizeof(double));
1304 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1305 kernel->positive_range = 1.0;
1308 kernel->minimum = 0.0;
1309 kernel->maximum = kernel->values[0];
1310 kernel->negative_range = 0.0;
1312 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1313 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1316 case BinomialKernel:
1321 if (args->rho < 1.0)
1322 kernel->width = kernel->height = 3; /* default radius = 1 */
1324 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1325 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1327 order_f = fact(kernel->width-1);
1329 kernel->values=(MagickRealType *) MagickAssumeAligned(
1330 AcquireAlignedMemory(kernel->width,kernel->height*
1331 sizeof(*kernel->values)));
1332 if (kernel->values == (MagickRealType *) NULL)
1333 return(DestroyKernelInfo(kernel));
1335 /* set all kernel values within diamond area to scale given */
1336 for ( i=0, v=0; v < (ssize_t)kernel->height; v++)
1338 alpha = order_f / ( fact((size_t) v) * fact(kernel->height-v-1) );
1339 for ( u=0; u < (ssize_t)kernel->width; u++, i++)
1340 kernel->positive_range += kernel->values[i] = (double)
1341 (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-u-1) ));
1343 kernel->minimum = 1.0;
1344 kernel->maximum = kernel->values[kernel->x+kernel->y*kernel->width];
1345 kernel->negative_range = 0.0;
1350 Convolution Kernels - Well Known Named Constant Kernels
1352 case LaplacianKernel:
1353 { switch ( (int) args->rho ) {
1355 default: /* laplacian square filter -- default */
1356 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1358 case 1: /* laplacian diamond filter */
1359 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1362 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1365 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1367 case 5: /* a 5x5 laplacian */
1368 kernel=ParseKernelArray(
1369 "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");
1371 case 7: /* a 7x7 laplacian */
1372 kernel=ParseKernelArray(
1373 "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" );
1375 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1376 kernel=ParseKernelArray(
1377 "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");
1379 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1380 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1381 kernel=ParseKernelArray(
1382 "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");
1385 if (kernel == (KernelInfo *) NULL)
1387 kernel->type = type;
1391 { /* Simple Sobel Kernel */
1392 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1393 if (kernel == (KernelInfo *) NULL)
1395 kernel->type = type;
1396 RotateKernelInfo(kernel, args->rho);
1401 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1402 if (kernel == (KernelInfo *) NULL)
1404 kernel->type = type;
1405 RotateKernelInfo(kernel, args->rho);
1410 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1411 if (kernel == (KernelInfo *) NULL)
1413 kernel->type = type;
1414 RotateKernelInfo(kernel, args->rho);
1419 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1420 if (kernel == (KernelInfo *) NULL)
1422 kernel->type = type;
1423 RotateKernelInfo(kernel, args->rho);
1428 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1429 if (kernel == (KernelInfo *) NULL)
1431 kernel->type = type;
1432 RotateKernelInfo(kernel, args->rho);
1435 case FreiChenKernel:
1436 /* Direction is set to be left to right positive */
1437 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1438 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1439 { switch ( (int) args->rho ) {
1442 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1443 if (kernel == (KernelInfo *) NULL)
1445 kernel->type = type;
1446 kernel->values[3] = +(MagickRealType) MagickSQ2;
1447 kernel->values[5] = -(MagickRealType) MagickSQ2;
1448 CalcKernelMetaData(kernel); /* recalculate meta-data */
1451 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1452 if (kernel == (KernelInfo *) NULL)
1454 kernel->type = type;
1455 kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2;
1456 kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2;
1457 CalcKernelMetaData(kernel); /* recalculate meta-data */
1458 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1461 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1462 if (kernel == (KernelInfo *) NULL)
1467 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1468 if (kernel == (KernelInfo *) NULL)
1470 kernel->type = type;
1471 kernel->values[3] = +(MagickRealType) MagickSQ2;
1472 kernel->values[5] = -(MagickRealType) MagickSQ2;
1473 CalcKernelMetaData(kernel); /* recalculate meta-data */
1474 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1477 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1478 if (kernel == (KernelInfo *) NULL)
1480 kernel->type = type;
1481 kernel->values[1] = +(double) MagickSQ2;
1482 kernel->values[7] = +(double) MagickSQ2;
1483 CalcKernelMetaData(kernel);
1484 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1487 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1488 if (kernel == (KernelInfo *) NULL)
1490 kernel->type = type;
1491 kernel->values[0] = +(MagickRealType) MagickSQ2;
1492 kernel->values[8] = -(MagickRealType) MagickSQ2;
1493 CalcKernelMetaData(kernel);
1494 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1497 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1498 if (kernel == (KernelInfo *) NULL)
1500 kernel->type = type;
1501 kernel->values[2] = -(MagickRealType) MagickSQ2;
1502 kernel->values[6] = +(MagickRealType) MagickSQ2;
1503 CalcKernelMetaData(kernel);
1504 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1507 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1508 if (kernel == (KernelInfo *) NULL)
1510 kernel->type = type;
1511 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1514 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1515 if (kernel == (KernelInfo *) NULL)
1517 kernel->type = type;
1518 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1521 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1522 if (kernel == (KernelInfo *) NULL)
1524 kernel->type = type;
1525 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1528 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1529 if (kernel == (KernelInfo *) NULL)
1531 kernel->type = type;
1532 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1535 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1536 if (kernel == (KernelInfo *) NULL)
1538 kernel->type = type;
1539 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1542 if ( fabs(args->sigma) >= MagickEpsilon )
1543 /* Rotate by correctly supplied 'angle' */
1544 RotateKernelInfo(kernel, args->sigma);
1545 else if ( args->rho > 30.0 || args->rho < -30.0 )
1546 /* Rotate by out of bounds 'type' */
1547 RotateKernelInfo(kernel, args->rho);
1552 Boolean or Shaped Kernels
1556 if (args->rho < 1.0)
1557 kernel->width = kernel->height = 3; /* default radius = 1 */
1559 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1560 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1562 kernel->values=(MagickRealType *) MagickAssumeAligned(
1563 AcquireAlignedMemory(kernel->width,kernel->height*
1564 sizeof(*kernel->values)));
1565 if (kernel->values == (MagickRealType *) NULL)
1566 return(DestroyKernelInfo(kernel));
1568 /* set all kernel values within diamond area to scale given */
1569 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1570 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1571 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1572 kernel->positive_range += kernel->values[i] = args->sigma;
1574 kernel->values[i] = nan;
1575 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1579 case RectangleKernel:
1582 if ( type == SquareKernel )
1584 if (args->rho < 1.0)
1585 kernel->width = kernel->height = 3; /* default radius = 1 */
1587 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1588 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1589 scale = args->sigma;
1592 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1593 if ( args->rho < 1.0 || args->sigma < 1.0 )
1594 return(DestroyKernelInfo(kernel)); /* invalid args given */
1595 kernel->width = (size_t)args->rho;
1596 kernel->height = (size_t)args->sigma;
1597 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1598 args->psi < 0.0 || args->psi > (double)kernel->height )
1599 return(DestroyKernelInfo(kernel)); /* invalid args given */
1600 kernel->x = (ssize_t) args->xi;
1601 kernel->y = (ssize_t) args->psi;
1604 kernel->values=(MagickRealType *) MagickAssumeAligned(
1605 AcquireAlignedMemory(kernel->width,kernel->height*
1606 sizeof(*kernel->values)));
1607 if (kernel->values == (MagickRealType *) NULL)
1608 return(DestroyKernelInfo(kernel));
1610 /* set all kernel values to scale given */
1611 u=(ssize_t) (kernel->width*kernel->height);
1612 for ( i=0; i < u; i++)
1613 kernel->values[i] = scale;
1614 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1615 kernel->positive_range = scale*u;
1620 if (args->rho < 1.0)
1621 kernel->width = kernel->height = 5; /* default radius = 2 */
1623 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1624 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1626 kernel->values=(MagickRealType *) MagickAssumeAligned(
1627 AcquireAlignedMemory(kernel->width,kernel->height*
1628 sizeof(*kernel->values)));
1629 if (kernel->values == (MagickRealType *) NULL)
1630 return(DestroyKernelInfo(kernel));
1632 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1633 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1634 if ( (labs((long) u)+labs((long) v)) <=
1635 ((long)kernel->x + (long)(kernel->x/2)) )
1636 kernel->positive_range += kernel->values[i] = args->sigma;
1638 kernel->values[i] = nan;
1639 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1645 limit = (ssize_t)(args->rho*args->rho);
1647 if (args->rho < 0.4) /* default radius approx 4.3 */
1648 kernel->width = kernel->height = 9L, limit = 18L;
1650 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1651 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1653 kernel->values=(MagickRealType *) MagickAssumeAligned(
1654 AcquireAlignedMemory(kernel->width,kernel->height*
1655 sizeof(*kernel->values)));
1656 if (kernel->values == (MagickRealType *) NULL)
1657 return(DestroyKernelInfo(kernel));
1659 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1660 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1661 if ((u*u+v*v) <= limit)
1662 kernel->positive_range += kernel->values[i] = args->sigma;
1664 kernel->values[i] = nan;
1665 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1670 if (args->rho < 1.0)
1671 kernel->width = kernel->height = 5; /* default radius 2 */
1673 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1674 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1676 kernel->values=(MagickRealType *) MagickAssumeAligned(
1677 AcquireAlignedMemory(kernel->width,kernel->height*
1678 sizeof(*kernel->values)));
1679 if (kernel->values == (MagickRealType *) NULL)
1680 return(DestroyKernelInfo(kernel));
1682 /* set all kernel values along axises to given scale */
1683 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1684 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1685 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1686 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1687 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1692 if (args->rho < 1.0)
1693 kernel->width = kernel->height = 5; /* default radius 2 */
1695 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1696 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1698 kernel->values=(MagickRealType *) MagickAssumeAligned(
1699 AcquireAlignedMemory(kernel->width,kernel->height*
1700 sizeof(*kernel->values)));
1701 if (kernel->values == (MagickRealType *) NULL)
1702 return(DestroyKernelInfo(kernel));
1704 /* set all kernel values along axises to given scale */
1705 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1706 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1707 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1708 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1709 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1723 if (args->rho < args->sigma)
1725 kernel->width = ((size_t)args->sigma)*2+1;
1726 limit1 = (ssize_t)(args->rho*args->rho);
1727 limit2 = (ssize_t)(args->sigma*args->sigma);
1731 kernel->width = ((size_t)args->rho)*2+1;
1732 limit1 = (ssize_t)(args->sigma*args->sigma);
1733 limit2 = (ssize_t)(args->rho*args->rho);
1736 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1738 kernel->height = kernel->width;
1739 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1740 kernel->values=(MagickRealType *) MagickAssumeAligned(
1741 AcquireAlignedMemory(kernel->width,kernel->height*
1742 sizeof(*kernel->values)));
1743 if (kernel->values == (MagickRealType *) NULL)
1744 return(DestroyKernelInfo(kernel));
1746 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1747 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1748 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1749 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1750 { ssize_t radius=u*u+v*v;
1751 if (limit1 < radius && radius <= limit2)
1752 kernel->positive_range += kernel->values[i] = (double) scale;
1754 kernel->values[i] = nan;
1756 kernel->minimum = kernel->maximum = (double) scale;
1757 if ( type == PeaksKernel ) {
1758 /* set the central point in the middle */
1759 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1760 kernel->positive_range = 1.0;
1761 kernel->maximum = 1.0;
1767 kernel=AcquireKernelInfo("ThinSE:482");
1768 if (kernel == (KernelInfo *) NULL)
1770 kernel->type = type;
1771 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1776 kernel=AcquireKernelInfo("ThinSE:87");
1777 if (kernel == (KernelInfo *) NULL)
1779 kernel->type = type;
1780 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1783 case DiagonalsKernel:
1785 switch ( (int) args->rho ) {
1790 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1791 if (kernel == (KernelInfo *) NULL)
1793 kernel->type = type;
1794 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1795 if (new_kernel == (KernelInfo *) NULL)
1796 return(DestroyKernelInfo(kernel));
1797 new_kernel->type = type;
1798 LastKernelInfo(kernel)->next = new_kernel;
1799 ExpandMirrorKernelInfo(kernel);
1803 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1806 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1809 if (kernel == (KernelInfo *) NULL)
1811 kernel->type = type;
1812 RotateKernelInfo(kernel, args->sigma);
1815 case LineEndsKernel:
1816 { /* Kernels for finding the end of thin lines */
1817 switch ( (int) args->rho ) {
1820 /* set of kernels to find all end of lines */
1821 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1823 /* kernel for 4-connected line ends - no rotation */
1824 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1827 /* kernel to add for 8-connected lines - no rotation */
1828 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1831 /* kernel to add for orthogonal line ends - does not find corners */
1832 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1835 /* traditional line end - fails on last T end */
1836 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1839 if (kernel == (KernelInfo *) NULL)
1841 kernel->type = type;
1842 RotateKernelInfo(kernel, args->sigma);
1845 case LineJunctionsKernel:
1846 { /* kernels for finding the junctions of multiple lines */
1847 switch ( (int) args->rho ) {
1850 /* set of kernels to find all line junctions */
1851 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1854 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1857 /* Diagonal T Junctions */
1858 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1861 /* Orthogonal T Junctions */
1862 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1865 /* Diagonal X Junctions */
1866 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1869 /* Orthogonal X Junctions - minimal diamond kernel */
1870 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1873 if (kernel == (KernelInfo *) NULL)
1875 kernel->type = type;
1876 RotateKernelInfo(kernel, args->sigma);
1880 { /* Ridges - Ridge finding kernels */
1883 switch ( (int) args->rho ) {
1886 kernel=ParseKernelArray("3x1:0,1,0");
1887 if (kernel == (KernelInfo *) NULL)
1889 kernel->type = type;
1890 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1893 kernel=ParseKernelArray("4x1:0,1,1,0");
1894 if (kernel == (KernelInfo *) NULL)
1896 kernel->type = type;
1897 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1899 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1900 /* Unfortunatally we can not yet rotate a non-square kernel */
1901 /* But then we can't flip a non-symetrical kernel either */
1902 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1903 if (new_kernel == (KernelInfo *) NULL)
1904 return(DestroyKernelInfo(kernel));
1905 new_kernel->type = type;
1906 LastKernelInfo(kernel)->next = new_kernel;
1907 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1908 if (new_kernel == (KernelInfo *) NULL)
1909 return(DestroyKernelInfo(kernel));
1910 new_kernel->type = type;
1911 LastKernelInfo(kernel)->next = new_kernel;
1912 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1913 if (new_kernel == (KernelInfo *) NULL)
1914 return(DestroyKernelInfo(kernel));
1915 new_kernel->type = type;
1916 LastKernelInfo(kernel)->next = new_kernel;
1917 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1918 if (new_kernel == (KernelInfo *) NULL)
1919 return(DestroyKernelInfo(kernel));
1920 new_kernel->type = type;
1921 LastKernelInfo(kernel)->next = new_kernel;
1922 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1923 if (new_kernel == (KernelInfo *) NULL)
1924 return(DestroyKernelInfo(kernel));
1925 new_kernel->type = type;
1926 LastKernelInfo(kernel)->next = new_kernel;
1927 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1928 if (new_kernel == (KernelInfo *) NULL)
1929 return(DestroyKernelInfo(kernel));
1930 new_kernel->type = type;
1931 LastKernelInfo(kernel)->next = new_kernel;
1932 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1933 if (new_kernel == (KernelInfo *) NULL)
1934 return(DestroyKernelInfo(kernel));
1935 new_kernel->type = type;
1936 LastKernelInfo(kernel)->next = new_kernel;
1937 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1938 if (new_kernel == (KernelInfo *) NULL)
1939 return(DestroyKernelInfo(kernel));
1940 new_kernel->type = type;
1941 LastKernelInfo(kernel)->next = new_kernel;
1946 case ConvexHullKernel:
1950 /* first set of 8 kernels */
1951 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1952 if (kernel == (KernelInfo *) NULL)
1954 kernel->type = type;
1955 ExpandRotateKernelInfo(kernel, 90.0);
1956 /* append the mirror versions too - no flip function yet */
1957 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1958 if (new_kernel == (KernelInfo *) NULL)
1959 return(DestroyKernelInfo(kernel));
1960 new_kernel->type = type;
1961 ExpandRotateKernelInfo(new_kernel, 90.0);
1962 LastKernelInfo(kernel)->next = new_kernel;
1965 case SkeletonKernel:
1967 switch ( (int) args->rho ) {
1970 /* Traditional Skeleton...
1971 ** A cyclically rotated single kernel
1973 kernel=AcquireKernelInfo("ThinSE:482");
1974 if (kernel == (KernelInfo *) NULL)
1976 kernel->type = type;
1977 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1980 /* HIPR Variation of the cyclic skeleton
1981 ** Corners of the traditional method made more forgiving,
1982 ** but the retain the same cyclic order.
1984 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1985 if (kernel == (KernelInfo *) NULL)
1987 if (kernel->next == (KernelInfo *) NULL)
1988 return(DestroyKernelInfo(kernel));
1989 kernel->type = type;
1990 kernel->next->type = type;
1991 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1994 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1995 ** "Connectivity-Preserving Morphological Image Thransformations"
1996 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1997 ** http://www.leptonica.com/papers/conn.pdf
1999 kernel=AcquireKernelInfo(
2000 "ThinSE:41; ThinSE:42; ThinSE:43");
2001 if (kernel == (KernelInfo *) NULL)
2003 kernel->type = type;
2004 kernel->next->type = type;
2005 kernel->next->next->type = type;
2006 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
2012 { /* Special kernels for general thinning, while preserving connections
2013 ** "Connectivity-Preserving Morphological Image Thransformations"
2014 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
2015 ** http://www.leptonica.com/papers/conn.pdf
2017 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2019 ** Note kernels do not specify the origin pixel, allowing them
2020 ** to be used for both thickening and thinning operations.
2022 switch ( (int) args->rho ) {
2023 /* SE for 4-connected thinning */
2024 case 41: /* SE_4_1 */
2025 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
2027 case 42: /* SE_4_2 */
2028 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
2030 case 43: /* SE_4_3 */
2031 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
2033 case 44: /* SE_4_4 */
2034 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
2036 case 45: /* SE_4_5 */
2037 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
2039 case 46: /* SE_4_6 */
2040 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
2042 case 47: /* SE_4_7 */
2043 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
2045 case 48: /* SE_4_8 */
2046 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
2048 case 49: /* SE_4_9 */
2049 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
2051 /* SE for 8-connected thinning - negatives of the above */
2052 case 81: /* SE_8_0 */
2053 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
2055 case 82: /* SE_8_2 */
2056 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
2058 case 83: /* SE_8_3 */
2059 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
2061 case 84: /* SE_8_4 */
2062 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
2064 case 85: /* SE_8_5 */
2065 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2067 case 86: /* SE_8_6 */
2068 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2070 case 87: /* SE_8_7 */
2071 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2073 case 88: /* SE_8_8 */
2074 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2076 case 89: /* SE_8_9 */
2077 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2079 /* Special combined SE kernels */
2080 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2081 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2083 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2084 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2086 case 481: /* SE_48_1 - General Connected Corner Kernel */
2087 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2090 case 482: /* SE_48_2 - General Edge Kernel */
2091 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2094 if (kernel == (KernelInfo *) NULL)
2096 kernel->type = type;
2097 RotateKernelInfo(kernel, args->sigma);
2101 Distance Measuring Kernels
2103 case ChebyshevKernel:
2105 if (args->rho < 1.0)
2106 kernel->width = kernel->height = 3; /* default radius = 1 */
2108 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2109 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2111 kernel->values=(MagickRealType *) MagickAssumeAligned(
2112 AcquireAlignedMemory(kernel->width,kernel->height*
2113 sizeof(*kernel->values)));
2114 if (kernel->values == (MagickRealType *) NULL)
2115 return(DestroyKernelInfo(kernel));
2117 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2118 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2119 kernel->positive_range += ( kernel->values[i] =
2120 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2121 kernel->maximum = kernel->values[0];
2124 case ManhattanKernel:
2126 if (args->rho < 1.0)
2127 kernel->width = kernel->height = 3; /* default radius = 1 */
2129 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2130 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2132 kernel->values=(MagickRealType *) MagickAssumeAligned(
2133 AcquireAlignedMemory(kernel->width,kernel->height*
2134 sizeof(*kernel->values)));
2135 if (kernel->values == (MagickRealType *) NULL)
2136 return(DestroyKernelInfo(kernel));
2138 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2139 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2140 kernel->positive_range += ( kernel->values[i] =
2141 args->sigma*(labs((long) u)+labs((long) v)) );
2142 kernel->maximum = kernel->values[0];
2145 case OctagonalKernel:
2147 if (args->rho < 2.0)
2148 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2150 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2151 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2153 kernel->values=(MagickRealType *) MagickAssumeAligned(
2154 AcquireAlignedMemory(kernel->width,kernel->height*
2155 sizeof(*kernel->values)));
2156 if (kernel->values == (MagickRealType *) NULL)
2157 return(DestroyKernelInfo(kernel));
2159 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2160 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2163 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2164 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2165 kernel->positive_range += kernel->values[i] =
2166 args->sigma*MagickMax(r1,r2);
2168 kernel->maximum = kernel->values[0];
2171 case EuclideanKernel:
2173 if (args->rho < 1.0)
2174 kernel->width = kernel->height = 3; /* default radius = 1 */
2176 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2177 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2179 kernel->values=(MagickRealType *) MagickAssumeAligned(
2180 AcquireAlignedMemory(kernel->width,kernel->height*
2181 sizeof(*kernel->values)));
2182 if (kernel->values == (MagickRealType *) NULL)
2183 return(DestroyKernelInfo(kernel));
2185 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2186 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2187 kernel->positive_range += ( kernel->values[i] =
2188 args->sigma*sqrt((double)(u*u+v*v)) );
2189 kernel->maximum = kernel->values[0];
2194 /* No-Op Kernel - Basically just a single pixel on its own */
2195 kernel=ParseKernelArray("1:1");
2196 if (kernel == (KernelInfo *) NULL)
2198 kernel->type = UndefinedKernel;
2207 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2211 % C l o n e K e r n e l I n f o %
2215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2217 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2218 % can be modified without effecting the original. The cloned kernel should
2219 % be destroyed using DestoryKernelInfo() when no longer needed.
2221 % The format of the CloneKernelInfo method is:
2223 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2225 % A description of each parameter follows:
2227 % o kernel: the Morphology/Convolution kernel to be cloned
2230 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2238 assert(kernel != (KernelInfo *) NULL);
2239 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2240 if (new_kernel == (KernelInfo *) NULL)
2242 *new_kernel=(*kernel); /* copy values in structure */
2244 /* replace the values with a copy of the values */
2245 new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2246 AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2247 if (new_kernel->values == (MagickRealType *) NULL)
2248 return(DestroyKernelInfo(new_kernel));
2249 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2250 new_kernel->values[i]=kernel->values[i];
2252 /* Also clone the next kernel in the kernel list */
2253 if ( kernel->next != (KernelInfo *) NULL ) {
2254 new_kernel->next = CloneKernelInfo(kernel->next);
2255 if ( new_kernel->next == (KernelInfo *) NULL )
2256 return(DestroyKernelInfo(new_kernel));
2263 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2267 % D e s t r o y K e r n e l I n f o %
2271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2273 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2276 % The format of the DestroyKernelInfo method is:
2278 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2280 % A description of each parameter follows:
2282 % o kernel: the Morphology/Convolution kernel to be destroyed
2285 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2287 assert(kernel != (KernelInfo *) NULL);
2288 if ( kernel->next != (KernelInfo *) NULL )
2289 kernel->next=DestroyKernelInfo(kernel->next);
2290 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2291 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2296 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2300 + E x p a n d M i r r o r K e r n e l I n f o %
2304 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2306 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2307 % sequence of 90-degree rotated kernels but providing a reflected 180
2308 % rotatation, before the -/+ 90-degree rotations.
2310 % This special rotation order produces a better, more symetrical thinning of
2313 % The format of the ExpandMirrorKernelInfo method is:
2315 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2317 % A description of each parameter follows:
2319 % o kernel: the Morphology/Convolution kernel
2321 % This function is only internel to this module, as it is not finalized,
2322 % especially with regard to non-orthogonal angles, and rotation of larger
2327 static void FlopKernelInfo(KernelInfo *kernel)
2328 { /* Do a Flop by reversing each row. */
2336 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2337 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2338 t=k[x], k[x]=k[r], k[r]=t;
2340 kernel->x = kernel->width - kernel->x - 1;
2341 angle = fmod(angle+180.0, 360.0);
2345 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2353 clone = CloneKernelInfo(last);
2354 RotateKernelInfo(clone, 180); /* flip */
2355 LastKernelInfo(last)->next = clone;
2358 clone = CloneKernelInfo(last);
2359 RotateKernelInfo(clone, 90); /* transpose */
2360 LastKernelInfo(last)->next = clone;
2363 clone = CloneKernelInfo(last);
2364 RotateKernelInfo(clone, 180); /* flop */
2365 LastKernelInfo(last)->next = clone;
2371 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2375 + E x p a n d R o t a t e K e r n e l I n f o %
2379 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2381 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2382 % incrementally by the angle given, until the kernel repeats.
2384 % WARNING: 45 degree rotations only works for 3x3 kernels.
2385 % While 90 degree roatations only works for linear and square kernels
2387 % The format of the ExpandRotateKernelInfo method is:
2389 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2391 % A description of each parameter follows:
2393 % o kernel: the Morphology/Convolution kernel
2395 % o angle: angle to rotate in degrees
2397 % This function is only internel to this module, as it is not finalized,
2398 % especially with regard to non-orthogonal angles, and rotation of larger
2402 /* Internal Routine - Return true if two kernels are the same */
2403 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2404 const KernelInfo *kernel2)
2409 /* check size and origin location */
2410 if ( kernel1->width != kernel2->width
2411 || kernel1->height != kernel2->height
2412 || kernel1->x != kernel2->x
2413 || kernel1->y != kernel2->y )
2416 /* check actual kernel values */
2417 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2418 /* Test for Nan equivalence */
2419 if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2421 if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2423 /* Test actual values are equivalent */
2424 if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2431 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2439 clone = CloneKernelInfo(last);
2440 RotateKernelInfo(clone, angle);
2441 if ( SameKernelInfo(kernel, clone) == MagickTrue )
2443 LastKernelInfo(last)->next = clone;
2446 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2451 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2455 + C a l c M e t a K e r n a l I n f o %
2459 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2461 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2462 % using the kernel values. This should only ne used if it is not possible to
2463 % calculate that meta-data in some easier way.
2465 % It is important that the meta-data is correct before ScaleKernelInfo() is
2466 % used to perform kernel normalization.
2468 % The format of the CalcKernelMetaData method is:
2470 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2472 % A description of each parameter follows:
2474 % o kernel: the Morphology/Convolution kernel to modify
2476 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2477 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2478 % however is not true for flat-shaped morphological kernels.
2480 % WARNING: Only the specific kernel pointed to is modified, not a list of
2483 % This is an internal function and not expected to be useful outside this
2484 % module. This could change however.
2486 static void CalcKernelMetaData(KernelInfo *kernel)
2491 kernel->minimum = kernel->maximum = 0.0;
2492 kernel->negative_range = kernel->positive_range = 0.0;
2493 for (i=0; i < (kernel->width*kernel->height); i++)
2495 if ( fabs(kernel->values[i]) < MagickEpsilon )
2496 kernel->values[i] = 0.0;
2497 ( kernel->values[i] < 0)
2498 ? ( kernel->negative_range += kernel->values[i] )
2499 : ( kernel->positive_range += kernel->values[i] );
2500 Minimize(kernel->minimum, kernel->values[i]);
2501 Maximize(kernel->maximum, kernel->values[i]);
2508 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2512 % M o r p h o l o g y A p p l y %
2516 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2518 % MorphologyApply() applies a morphological method, multiple times using
2519 % a list of multiple kernels. This is the method that should be called by
2520 % other 'operators' that internally use morphology operations as part of
2523 % It is basically equivalent to as MorphologyImage() (see below) but
2524 % without any user controls. This allows internel programs to use this
2525 % function, to actually perform a specific task without possible interference
2526 % by any API user supplied settings.
2528 % It is MorphologyImage() task to extract any such user controls, and
2529 % pass them to this function for processing.
2531 % More specifically all given kernels should already be scaled, normalised,
2532 % and blended appropriatally before being parred to this routine. The
2533 % appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2535 % The format of the MorphologyApply method is:
2537 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2538 % const ssize_t iterations,const KernelInfo *kernel,
2539 % const CompositeMethod compose,const double bias,
2540 % ExceptionInfo *exception)
2542 % A description of each parameter follows:
2544 % o image: the source image
2546 % o method: the morphology method to be applied.
2548 % o iterations: apply the operation this many times (or no change).
2549 % A value of -1 means loop until no change found.
2550 % How this is applied may depend on the morphology method.
2551 % Typically this is a value of 1.
2553 % o channel: the channel type.
2555 % o kernel: An array of double representing the morphology kernel.
2557 % o compose: How to handle or merge multi-kernel results.
2558 % If 'UndefinedCompositeOp' use default for the Morphology method.
2559 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2560 % Otherwise merge the results using the compose method given.
2562 % o bias: Convolution Output Bias.
2564 % o exception: return any errors or warnings in this structure.
2568 /* Apply a Morphology Primative to an image using the given kernel.
2569 ** Two pre-created images must be provided, and no image is created.
2570 ** It returns the number of pixels that changed between the images
2571 ** for result convergence determination.
2573 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2574 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2575 ExceptionInfo *exception)
2577 #define MorphologyTag "Morphology/Image"
2596 assert(image != (Image *) NULL);
2597 assert(image->signature == MagickSignature);
2598 assert(morphology_image != (Image *) NULL);
2599 assert(morphology_image->signature == MagickSignature);
2600 assert(kernel != (KernelInfo *) NULL);
2601 assert(kernel->signature == MagickSignature);
2602 assert(exception != (ExceptionInfo *) NULL);
2603 assert(exception->signature == MagickSignature);
2609 image_view=AcquireVirtualCacheView(image,exception);
2610 morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2611 virt_width=image->columns+kernel->width-1;
2613 /* Some methods (including convolve) needs use a reflected kernel.
2614 * Adjust 'origin' offsets to loop though kernel as a reflection.
2619 case ConvolveMorphology:
2620 case DilateMorphology:
2621 case DilateIntensityMorphology:
2622 case IterativeDistanceMorphology:
2623 /* kernel needs to used with reflection about origin */
2624 offx = (ssize_t) kernel->width-offx-1;
2625 offy = (ssize_t) kernel->height-offy-1;
2627 case ErodeMorphology:
2628 case ErodeIntensityMorphology:
2629 case HitAndMissMorphology:
2630 case ThinningMorphology:
2631 case ThickenMorphology:
2632 /* kernel is used as is, without reflection */
2635 assert("Not a Primitive Morphology Method" != (char *) NULL);
2639 if ( method == ConvolveMorphology && kernel->width == 1 )
2640 { /* Special handling (for speed) of vertical (blur) kernels.
2641 ** This performs its handling in columns rather than in rows.
2642 ** This is only done for convolve as it is the only method that
2643 ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
2645 ** Timing tests (on single CPU laptop)
2646 ** Using a vertical 1-d Blue with normal row-by-row (below)
2647 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2649 ** Using this column method
2650 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2653 ** Anthony Thyssen, 14 June 2010
2658 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2659 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2660 dynamic_number_threads(image,image->columns,image->rows,1)
2662 for (x=0; x < (ssize_t) image->columns; x++)
2664 register const Quantum
2676 if (status == MagickFalse)
2678 p=GetCacheViewVirtualPixels(image_view,x,-offy,1,image->rows+
2679 kernel->height-1,exception);
2680 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2681 morphology_image->rows,exception);
2682 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2687 /* offset to origin in 'p'. while 'q' points to it directly */
2690 for (y=0; y < (ssize_t) image->rows; y++)
2695 register const MagickRealType
2698 register const Quantum
2704 /* Copy input image to the output image for unused channels
2705 * This removes need for 'cloning' a new image every iteration
2707 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2708 GetPixelChannels(image)),q);
2709 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2710 GetPixelChannels(image)),q);
2711 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2712 GetPixelChannels(image)),q);
2713 if (image->colorspace == CMYKColorspace)
2714 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2715 GetPixelChannels(image)),q);
2717 /* Set the bias of the weighted average output */
2722 result.black = bias;
2725 /* Weighted Average of pixels using reflected kernel
2727 ** NOTE for correct working of this operation for asymetrical
2728 ** kernels, the kernel needs to be applied in its reflected form.
2729 ** That is its values needs to be reversed.
2731 k = &kernel->values[ kernel->height-1 ];
2733 if ( (image->channel_mask != DefaultChannels) ||
2734 (image->alpha_trait != BlendPixelTrait) )
2735 { /* No 'Sync' involved.
2736 ** Convolution is just a simple greyscale channel operation
2738 for (v=0; v < (ssize_t) kernel->height; v++) {
2739 if ( IsNan(*k) ) continue;
2740 result.red += (*k)*GetPixelRed(image,k_pixels);
2741 result.green += (*k)*GetPixelGreen(image,k_pixels);
2742 result.blue += (*k)*GetPixelBlue(image,k_pixels);
2743 if (image->colorspace == CMYKColorspace)
2744 result.black+=(*k)*GetPixelBlack(image,k_pixels);
2745 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2747 k_pixels+=GetPixelChannels(image);
2749 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2750 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
2751 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2752 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
2753 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2754 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
2755 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2756 (image->colorspace == CMYKColorspace))
2757 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
2758 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2759 (image->alpha_trait == BlendPixelTrait))
2760 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2763 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2764 ** Weight the color channels with Alpha Channel so that
2765 ** transparent pixels are not part of the results.
2768 alpha, /* alpha weighting for colors : alpha */
2769 gamma; /* divisor, sum of color alpha weighting */
2771 count; /* alpha valus collected, number kernel values */
2775 for (v=0; v < (ssize_t) kernel->height; v++) {
2776 if ( IsNan(*k) ) continue;
2777 alpha=QuantumScale*GetPixelAlpha(image,k_pixels);
2778 gamma += alpha; /* normalize alpha weights only */
2779 count++; /* number of alpha values collected */
2780 alpha*=(*k); /* include kernel weighting now */
2781 result.red += alpha*GetPixelRed(image,k_pixels);
2782 result.green += alpha*GetPixelGreen(image,k_pixels);
2783 result.blue += alpha*GetPixelBlue(image,k_pixels);
2784 if (image->colorspace == CMYKColorspace)
2785 result.black += alpha*GetPixelBlack(image,k_pixels);
2786 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2788 k_pixels+=GetPixelChannels(image);
2790 /* Sync'ed channels, all channels are modified */
2791 gamma=(double)count/(fabs((double) gamma) < MagickEpsilon ? MagickEpsilon : gamma);
2792 SetPixelRed(morphology_image,ClampToQuantum(gamma*result.red),q);
2793 SetPixelGreen(morphology_image,ClampToQuantum(gamma*result.green),q);
2794 SetPixelBlue(morphology_image,ClampToQuantum(gamma*result.blue),q);
2795 if (image->colorspace == CMYKColorspace)
2796 SetPixelBlack(morphology_image,ClampToQuantum(gamma*result.black),q);
2797 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2800 /* Count up changed pixels */
2801 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q))
2802 || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q))
2803 || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q))
2804 || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q))
2805 || ((image->colorspace == CMYKColorspace) &&
2806 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
2807 changed++; /* The pixel was changed in some way! */
2808 p+=GetPixelChannels(image);
2809 q+=GetPixelChannels(morphology_image);
2811 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2813 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2818 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2819 #pragma omp critical (MagickCore_MorphologyImage)
2821 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2822 if (proceed == MagickFalse)
2826 morphology_image->type=image->type;
2827 morphology_view=DestroyCacheView(morphology_view);
2828 image_view=DestroyCacheView(image_view);
2829 return(status ? (ssize_t) changed : 0);
2833 ** Normal handling of horizontal or rectangular kernels (row by row)
2835 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2836 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2837 dynamic_number_threads(image,image->columns,image->rows,1)
2839 for (y=0; y < (ssize_t) image->rows; y++)
2841 register const Quantum
2853 if (status == MagickFalse)
2855 p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
2856 kernel->height, exception);
2857 q=GetCacheViewAuthenticPixels(morphology_view,0,y,
2858 morphology_image->columns,1,exception);
2859 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2864 /* offset to origin in 'p'. while 'q' points to it directly */
2865 r = virt_width*offy + offx;
2867 for (x=0; x < (ssize_t) image->columns; x++)
2874 register const MagickRealType
2877 register const Quantum
2886 /* Copy input image to the output image for unused channels
2887 * This removes need for 'cloning' a new image every iteration
2889 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2890 GetPixelChannels(image)),q);
2891 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2892 GetPixelChannels(image)),q);
2893 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2894 GetPixelChannels(image)),q);
2895 if (image->colorspace == CMYKColorspace)
2896 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2897 GetPixelChannels(image)),q);
2904 min.black = (double) QuantumRange;
2909 max.black = (double) 0;
2910 /* default result is the original pixel value */
2911 result.red = (double) GetPixelRed(image,p+r*GetPixelChannels(image));
2912 result.green = (double) GetPixelGreen(image,p+r*GetPixelChannels(image));
2913 result.blue = (double) GetPixelBlue(image,p+r*GetPixelChannels(image));
2915 if (image->colorspace == CMYKColorspace)
2916 result.black = (double) GetPixelBlack(image,p+r*GetPixelChannels(image));
2917 result.alpha=(double) GetPixelAlpha(image,p+r*GetPixelChannels(image));
2920 case ConvolveMorphology:
2921 /* Set the bias of the weighted average output */
2926 result.black = bias;
2928 case DilateIntensityMorphology:
2929 case ErodeIntensityMorphology:
2930 /* use a boolean flag indicating when first match found */
2931 result.red = 0.0; /* result is not used otherwise */
2938 case ConvolveMorphology:
2939 /* Weighted Average of pixels using reflected kernel
2941 ** NOTE for correct working of this operation for asymetrical
2942 ** kernels, the kernel needs to be applied in its reflected form.
2943 ** That is its values needs to be reversed.
2945 ** Correlation is actually the same as this but without reflecting
2946 ** the kernel, and thus 'lower-level' that Convolution. However
2947 ** as Convolution is the more common method used, and it does not
2948 ** really cost us much in terms of processing to use a reflected
2949 ** kernel, so it is Convolution that is implemented.
2951 ** Correlation will have its kernel reflected before calling
2952 ** this function to do a Convolve.
2954 ** For more details of Correlation vs Convolution see
2955 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2957 k = &kernel->values[ kernel->width*kernel->height-1 ];
2959 if ( (image->channel_mask != DefaultChannels) ||
2960 (image->alpha_trait != BlendPixelTrait) )
2961 { /* No 'Sync' involved.
2962 ** Convolution is simple greyscale channel operation
2964 for (v=0; v < (ssize_t) kernel->height; v++) {
2965 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2966 if ( IsNan(*k) ) continue;
2968 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2969 result.green += (*k)*
2970 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2971 result.blue += (*k)*
2972 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2973 if (image->colorspace == CMYKColorspace)
2974 result.black += (*k)*
2975 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2976 result.alpha += (*k)*
2977 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2979 k_pixels += virt_width*GetPixelChannels(image);
2981 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2982 SetPixelRed(morphology_image,ClampToQuantum(result.red),
2984 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2985 SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2987 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2988 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2990 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2991 (image->colorspace == CMYKColorspace))
2992 SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2994 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2995 (image->alpha_trait == BlendPixelTrait))
2996 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
3000 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
3001 ** Weight the color channels with Alpha Channel so that
3002 ** transparent pixels are not part of the results.
3005 alpha, /* alpha weighting for colors : alpha */
3006 gamma; /* divisor, sum of color alpha weighting */
3008 count; /* alpha valus collected, number kernel values */
3012 for (v=0; v < (ssize_t) kernel->height; v++) {
3013 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3014 if ( IsNan(*k) ) continue;
3015 alpha=QuantumScale*GetPixelAlpha(image,
3016 k_pixels+u*GetPixelChannels(image));
3017 gamma += alpha; /* normalize alpha weights only */
3018 count++; /* number of alpha values collected */
3019 alpha=alpha*(*k); /* include kernel weighting now */
3020 result.red += alpha*
3021 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
3022 result.green += alpha*
3023 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
3024 result.blue += alpha*
3025 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
3026 if (image->colorspace == CMYKColorspace)
3027 result.black += alpha*
3028 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
3029 result.alpha += (*k)*
3030 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
3032 k_pixels += virt_width*GetPixelChannels(image);
3034 /* Sync'ed channels, all channels are modified */
3035 gamma=(double)count/(fabs((double) gamma) < MagickEpsilon ? MagickEpsilon : gamma);
3036 SetPixelRed(morphology_image,
3037 ClampToQuantum(gamma*result.red),q);
3038 SetPixelGreen(morphology_image,
3039 ClampToQuantum(gamma*result.green),q);
3040 SetPixelBlue(morphology_image,
3041 ClampToQuantum(gamma*result.blue),q);
3042 if (image->colorspace == CMYKColorspace)
3043 SetPixelBlack(morphology_image,
3044 ClampToQuantum(gamma*result.black),q);
3045 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3049 case ErodeMorphology:
3050 /* Minimum Value within kernel neighbourhood
3052 ** NOTE that the kernel is not reflected for this operation!
3054 ** NOTE: in normal Greyscale Morphology, the kernel value should
3055 ** be added to the real value, this is currently not done, due to
3056 ** the nature of the boolean kernels being used.
3060 for (v=0; v < (ssize_t) kernel->height; v++) {
3061 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3062 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3063 Minimize(min.red, (double)
3064 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3065 Minimize(min.green, (double)
3066 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3067 Minimize(min.blue, (double)
3068 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3069 Minimize(min.alpha, (double)
3070 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3071 if (image->colorspace == CMYKColorspace)
3072 Minimize(min.black, (double)
3073 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3075 k_pixels += virt_width*GetPixelChannels(image);
3079 case DilateMorphology:
3080 /* Maximum Value within kernel neighbourhood
3082 ** NOTE for correct working of this operation for asymetrical
3083 ** kernels, the kernel needs to be applied in its reflected form.
3084 ** That is its values needs to be reversed.
3086 ** NOTE: in normal Greyscale Morphology, the kernel value should
3087 ** be added to the real value, this is currently not done, due to
3088 ** the nature of the boolean kernels being used.
3091 k = &kernel->values[ kernel->width*kernel->height-1 ];
3093 for (v=0; v < (ssize_t) kernel->height; v++) {
3094 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3095 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3096 Maximize(max.red, (double)
3097 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3098 Maximize(max.green, (double)
3099 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3100 Maximize(max.blue, (double)
3101 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3102 Maximize(max.alpha, (double)
3103 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3104 if (image->colorspace == CMYKColorspace)
3105 Maximize(max.black, (double)
3106 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3108 k_pixels += virt_width*GetPixelChannels(image);
3112 case HitAndMissMorphology:
3113 case ThinningMorphology:
3114 case ThickenMorphology:
3115 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3117 ** NOTE that the kernel is not reflected for this operation,
3118 ** and consists of both foreground and background pixel
3119 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3120 ** with either Nan or 0.5 values for don't care.
3122 ** Note that this will never produce a meaningless negative
3123 ** result. Such results can cause Thinning/Thicken to not work
3124 ** correctly when used against a greyscale image.
3128 for (v=0; v < (ssize_t) kernel->height; v++) {
3129 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3130 if ( IsNan(*k) ) continue;
3132 { /* minimim of foreground pixels */
3133 Minimize(min.red, (double)
3134 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3135 Minimize(min.green, (double)
3136 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3137 Minimize(min.blue, (double)
3138 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3139 Minimize(min.alpha,(double)
3140 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3141 if ( image->colorspace == CMYKColorspace)
3142 Minimize(min.black,(double)
3143 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3145 else if ( (*k) < 0.3 )
3146 { /* maximum of background pixels */
3147 Maximize(max.red, (double)
3148 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3149 Maximize(max.green, (double)
3150 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3151 Maximize(max.blue, (double)
3152 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3153 Maximize(max.alpha,(double)
3154 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3155 if (image->colorspace == CMYKColorspace)
3156 Maximize(max.black, (double)
3157 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3160 k_pixels += virt_width*GetPixelChannels(image);
3162 /* Pattern Match if difference is positive */
3163 min.red -= max.red; Maximize( min.red, 0.0 );
3164 min.green -= max.green; Maximize( min.green, 0.0 );
3165 min.blue -= max.blue; Maximize( min.blue, 0.0 );
3166 min.black -= max.black; Maximize( min.black, 0.0 );
3167 min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3170 case ErodeIntensityMorphology:
3171 /* Select Pixel with Minimum Intensity within kernel neighbourhood
3173 ** WARNING: the intensity test fails for CMYK and does not
3174 ** take into account the moderating effect of the alpha channel
3175 ** on the intensity.
3177 ** NOTE that the kernel is not reflected for this operation!
3181 for (v=0; v < (ssize_t) kernel->height; v++) {
3182 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3183 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3184 if ( result.red == 0.0 ||
3185 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) < GetPixelIntensity(morphology_image,q) ) {
3186 /* copy the whole pixel - no channel selection */
3187 SetPixelRed(morphology_image,GetPixelRed(image,
3188 k_pixels+u*GetPixelChannels(image)),q);
3189 SetPixelGreen(morphology_image,GetPixelGreen(image,
3190 k_pixels+u*GetPixelChannels(image)),q);
3191 SetPixelBlue(morphology_image,GetPixelBlue(image,
3192 k_pixels+u*GetPixelChannels(image)),q);
3193 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3194 k_pixels+u*GetPixelChannels(image)),q);
3195 if ( result.red > 0.0 ) changed++;
3199 k_pixels += virt_width*GetPixelChannels(image);
3203 case DilateIntensityMorphology:
3204 /* Select Pixel with Maximum Intensity within kernel neighbourhood
3206 ** WARNING: the intensity test fails for CMYK and does not
3207 ** take into account the moderating effect of the alpha channel
3208 ** on the intensity (yet).
3210 ** NOTE for correct working of this operation for asymetrical
3211 ** kernels, the kernel needs to be applied in its reflected form.
3212 ** That is its values needs to be reversed.
3214 k = &kernel->values[ kernel->width*kernel->height-1 ];
3216 for (v=0; v < (ssize_t) kernel->height; v++) {
3217 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3218 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3219 if ( result.red == 0.0 ||
3220 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) > GetPixelIntensity(morphology_image,q) ) {
3221 /* copy the whole pixel - no channel selection */
3222 SetPixelRed(morphology_image,GetPixelRed(image,
3223 k_pixels+u*GetPixelChannels(image)),q);
3224 SetPixelGreen(morphology_image,GetPixelGreen(image,
3225 k_pixels+u*GetPixelChannels(image)),q);
3226 SetPixelBlue(morphology_image,GetPixelBlue(image,
3227 k_pixels+u*GetPixelChannels(image)),q);
3228 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3229 k_pixels+u*GetPixelChannels(image)),q);
3230 if ( result.red > 0.0 ) changed++;
3234 k_pixels += virt_width*GetPixelChannels(image);
3238 case IterativeDistanceMorphology:
3239 /* Work out an iterative distance from black edge of a white image
3240 ** shape. Essentually white values are decreased to the smallest
3241 ** 'distance from edge' it can find.
3243 ** It works by adding kernel values to the neighbourhood, and and
3244 ** select the minimum value found. The kernel is rotated before
3245 ** use, so kernel distances match resulting distances, when a user
3246 ** provided asymmetric kernel is applied.
3249 ** This code is almost identical to True GrayScale Morphology But
3252 ** GreyDilate Kernel values added, maximum value found Kernel is
3253 ** rotated before use.
3255 ** GrayErode: Kernel values subtracted and minimum value found No
3256 ** kernel rotation used.
3258 ** Note the the Iterative Distance method is essentially a
3259 ** GrayErode, but with negative kernel values, and kernel
3260 ** rotation applied.
3262 k = &kernel->values[ kernel->width*kernel->height-1 ];
3264 for (v=0; v < (ssize_t) kernel->height; v++) {
3265 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3266 if ( IsNan(*k) ) continue;
3267 Minimize(result.red, (*k)+(double)
3268 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3269 Minimize(result.green, (*k)+(double)
3270 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3271 Minimize(result.blue, (*k)+(double)
3272 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3273 Minimize(result.alpha, (*k)+(double)
3274 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3275 if ( image->colorspace == CMYKColorspace)
3276 Maximize(result.black, (*k)+(double)
3277 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3279 k_pixels += virt_width*GetPixelChannels(image);
3283 case UndefinedMorphology:
3285 break; /* Do nothing */
3287 /* Final mathematics of results (combine with original image?)
3289 ** NOTE: Difference Morphology operators Edge* and *Hat could also
3290 ** be done here but works better with iteration as a image difference
3291 ** in the controling function (below). Thicken and Thinning however
3292 ** should be done here so thay can be iterated correctly.
3295 case HitAndMissMorphology:
3296 case ErodeMorphology:
3297 result = min; /* minimum of neighbourhood */
3299 case DilateMorphology:
3300 result = max; /* maximum of neighbourhood */
3302 case ThinningMorphology:
3303 /* subtract pattern match from original */
3304 result.red -= min.red;
3305 result.green -= min.green;
3306 result.blue -= min.blue;
3307 result.black -= min.black;
3308 result.alpha -= min.alpha;
3310 case ThickenMorphology:
3311 /* Add the pattern matchs to the original */
3312 result.red += min.red;
3313 result.green += min.green;
3314 result.blue += min.blue;
3315 result.black += min.black;
3316 result.alpha += min.alpha;
3319 /* result directly calculated or assigned */
3322 /* Assign the resulting pixel values - Clamping Result */
3324 case UndefinedMorphology:
3325 case ConvolveMorphology:
3326 case DilateIntensityMorphology:
3327 case ErodeIntensityMorphology:
3328 break; /* full pixel was directly assigned - not a channel method */
3330 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3331 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
3332 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3333 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
3334 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3335 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
3336 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3337 (image->colorspace == CMYKColorspace))
3338 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
3339 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3340 (image->alpha_trait == BlendPixelTrait))
3341 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3344 /* Count up changed pixels */
3345 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q)) ||
3346 (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q)) ||
3347 (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q)) ||
3348 (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q)) ||
3349 ((image->colorspace == CMYKColorspace) &&
3350 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
3351 changed++; /* The pixel was changed in some way! */
3352 p+=GetPixelChannels(image);
3353 q+=GetPixelChannels(morphology_image);
3355 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3357 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3362 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3363 #pragma omp critical (MagickCore_MorphologyImage)
3365 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3366 if (proceed == MagickFalse)
3370 morphology_view=DestroyCacheView(morphology_view);
3371 image_view=DestroyCacheView(image_view);
3372 return(status ? (ssize_t)changed : -1);
3375 /* This is almost identical to the MorphologyPrimative() function above,
3376 ** but will apply the primitive directly to the actual image using two
3377 ** passes, once in each direction, with the results of the previous (and
3378 ** current) row being re-used.
3380 ** That is after each row is 'Sync'ed' into the image, the next row will
3381 ** make use of those values as part of the calculation of the next row.
3382 ** It then repeats, but going in the oppisite (bottom-up) direction.
3384 ** Because of this 're-use of results' this function can not make use
3385 ** of multi-threaded, parellel processing.
3387 static ssize_t MorphologyPrimitiveDirect(Image *image,
3388 const MorphologyMethod method,const KernelInfo *kernel,
3389 ExceptionInfo *exception)
3412 assert(image != (Image *) NULL);
3413 assert(image->signature == MagickSignature);
3414 assert(kernel != (KernelInfo *) NULL);
3415 assert(kernel->signature == MagickSignature);
3416 assert(exception != (ExceptionInfo *) NULL);
3417 assert(exception->signature == MagickSignature);
3419 /* Some methods (including convolve) needs use a reflected kernel.
3420 * Adjust 'origin' offsets to loop though kernel as a reflection.
3425 case DistanceMorphology:
3426 case VoronoiMorphology:
3427 /* kernel needs to used with reflection about origin */
3428 offx = (ssize_t) kernel->width-offx-1;
3429 offy = (ssize_t) kernel->height-offy-1;
3432 case ?????Morphology:
3433 /* kernel is used as is, without reflection */
3437 assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
3441 /* DO NOT THREAD THIS CODE! */
3442 /* two views into same image (virtual, and actual) */
3443 virt_view=AcquireVirtualCacheView(image,exception);
3444 auth_view=AcquireAuthenticCacheView(image,exception);
3445 virt_width=image->columns+kernel->width-1;
3447 for (y=0; y < (ssize_t) image->rows; y++)
3449 register const Quantum
3461 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3462 ** we read using virtual to get virtual pixel handling, but write back
3463 ** into the same image.
3465 ** Only top half of kernel is processed as we do a single pass downward
3466 ** through the image iterating the distance function as we go.
3468 if (status == MagickFalse)
3470 p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
3472 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3474 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3476 if (status == MagickFalse)
3479 /* offset to origin in 'p'. while 'q' points to it directly */
3480 r = (ssize_t) virt_width*offy + offx;
3482 for (x=0; x < (ssize_t) image->columns; x++)
3487 register const MagickRealType
3490 register const Quantum
3499 /* Starting Defaults */
3500 GetPixelInfo(image,&result);
3501 GetPixelInfoPixel(image,q,&result);
3502 if ( method != VoronoiMorphology )
3503 result.alpha = QuantumRange - result.alpha;
3506 case DistanceMorphology:
3507 /* Add kernel Value and select the minimum value found. */
3508 k = &kernel->values[ kernel->width*kernel->height-1 ];
3510 for (v=0; v <= (ssize_t) offy; v++) {
3511 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3512 if ( IsNan(*k) ) continue;
3513 Minimize(result.red, (*k)+
3514 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3515 Minimize(result.green, (*k)+
3516 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3517 Minimize(result.blue, (*k)+
3518 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3519 if (image->colorspace == CMYKColorspace)
3520 Minimize(result.black,(*k)+
3521 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3522 Minimize(result.alpha, (*k)+
3523 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3525 k_pixels += virt_width*GetPixelChannels(image);
3527 /* repeat with the just processed pixels of this row */
3528 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3529 k_pixels = q-offx*GetPixelChannels(image);
3530 for (u=0; u < (ssize_t) offx; u++, k--) {
3531 if ( x+u-offx < 0 ) continue; /* off the edge! */
3532 if ( IsNan(*k) ) continue;
3533 Minimize(result.red, (*k)+
3534 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3535 Minimize(result.green, (*k)+
3536 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3537 Minimize(result.blue, (*k)+
3538 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3539 if (image->colorspace == CMYKColorspace)
3540 Minimize(result.black,(*k)+
3541 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3542 Minimize(result.alpha,(*k)+
3543 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3546 case VoronoiMorphology:
3547 /* Apply Distance to 'Matte' channel, while coping the color
3548 ** values of the closest pixel.
3550 ** This is experimental, and realy the 'alpha' component should
3551 ** be completely separate 'masking' channel so that alpha can
3552 ** also be used as part of the results.
3554 k = &kernel->values[ kernel->width*kernel->height-1 ];
3556 for (v=0; v <= (ssize_t) offy; v++) {
3557 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3558 if ( IsNan(*k) ) continue;
3559 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3561 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3566 k_pixels += virt_width*GetPixelChannels(image);
3568 /* repeat with the just processed pixels of this row */
3569 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3570 k_pixels = q-offx*GetPixelChannels(image);
3571 for (u=0; u < (ssize_t) offx; u++, k--) {
3572 if ( x+u-offx < 0 ) continue; /* off the edge! */
3573 if ( IsNan(*k) ) continue;
3574 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3576 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3583 /* result directly calculated or assigned */
3586 /* Assign the resulting pixel values - Clamping Result */
3588 case VoronoiMorphology:
3589 SetPixelInfoPixel(image,&result,q);
3592 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3593 SetPixelRed(image,ClampToQuantum(result.red),q);
3594 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3595 SetPixelGreen(image,ClampToQuantum(result.green),q);
3596 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3597 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3598 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3599 (image->colorspace == CMYKColorspace))
3600 SetPixelBlack(image,ClampToQuantum(result.black),q);
3601 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3602 (image->alpha_trait == BlendPixelTrait))
3603 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3606 /* Count up changed pixels */
3607 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q)) ||
3608 (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q)) ||
3609 (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q)) ||
3610 (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q)) ||
3611 ((image->colorspace == CMYKColorspace) &&
3612 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3613 changed++; /* The pixel was changed in some way! */
3615 p+=GetPixelChannels(image); /* increment pixel buffers */
3616 q+=GetPixelChannels(image);
3619 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3621 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3622 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3628 /* Do the reversed pass through the image */
3629 for (y=(ssize_t)image->rows-1; y >= 0; y--)
3631 register const Quantum
3643 if (status == MagickFalse)
3645 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3646 ** we read using virtual to get virtual pixel handling, but write back
3647 ** into the same image.
3649 ** Only the bottom half of the kernel will be processes as we
3652 p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
3653 kernel->y+1,exception);
3654 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3656 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3658 if (status == MagickFalse)
3661 /* adjust positions to end of row */
3662 p += (image->columns-1)*GetPixelChannels(image);
3663 q += (image->columns-1)*GetPixelChannels(image);
3665 /* offset to origin in 'p'. while 'q' points to it directly */
3668 for (x=(ssize_t)image->columns-1; x >= 0; x--)
3673 register const MagickRealType
3676 register const Quantum
3685 /* Default - previously modified pixel */
3686 GetPixelInfo(image,&result);
3687 GetPixelInfoPixel(image,q,&result);
3688 if ( method != VoronoiMorphology )
3689 result.alpha = QuantumRange - result.alpha;
3692 case DistanceMorphology:
3693 /* Add kernel Value and select the minimum value found. */
3694 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3696 for (v=offy; v < (ssize_t) kernel->height; v++) {
3697 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3698 if ( IsNan(*k) ) continue;
3699 Minimize(result.red, (*k)+
3700 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3701 Minimize(result.green, (*k)+
3702 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3703 Minimize(result.blue, (*k)+
3704 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3705 if ( image->colorspace == CMYKColorspace)
3706 Minimize(result.black,(*k)+
3707 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3708 Minimize(result.alpha, (*k)+
3709 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3711 k_pixels += virt_width*GetPixelChannels(image);
3713 /* repeat with the just processed pixels of this row */
3714 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3715 k_pixels = q-offx*GetPixelChannels(image);
3716 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3717 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3718 if ( IsNan(*k) ) continue;
3719 Minimize(result.red, (*k)+
3720 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3721 Minimize(result.green, (*k)+
3722 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3723 Minimize(result.blue, (*k)+
3724 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3725 if ( image->colorspace == CMYKColorspace)
3726 Minimize(result.black, (*k)+
3727 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3728 Minimize(result.alpha, (*k)+
3729 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3732 case VoronoiMorphology:
3733 /* Apply Distance to 'Matte' channel, coping the closest color.
3735 ** This is experimental, and realy the 'alpha' component should
3736 ** be completely separate 'masking' channel.
3738 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3740 for (v=offy; v < (ssize_t) kernel->height; v++) {
3741 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3742 if ( IsNan(*k) ) continue;
3743 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3745 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3750 k_pixels += virt_width*GetPixelChannels(image);
3752 /* repeat with the just processed pixels of this row */
3753 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3754 k_pixels = q-offx*GetPixelChannels(image);
3755 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3756 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3757 if ( IsNan(*k) ) continue;
3758 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3760 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3767 /* result directly calculated or assigned */
3770 /* Assign the resulting pixel values - Clamping Result */
3772 case VoronoiMorphology:
3773 SetPixelInfoPixel(image,&result,q);
3776 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3777 SetPixelRed(image,ClampToQuantum(result.red),q);
3778 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3779 SetPixelGreen(image,ClampToQuantum(result.green),q);
3780 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3781 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3782 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3783 (image->colorspace == CMYKColorspace))
3784 SetPixelBlack(image,ClampToQuantum(result.black),q);
3785 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3786 (image->alpha_trait == BlendPixelTrait))
3787 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3790 /* Count up changed pixels */
3791 if ( (GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q))
3792 || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q))
3793 || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q))
3794 || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q))
3795 || ((image->colorspace == CMYKColorspace) &&
3796 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3797 changed++; /* The pixel was changed in some way! */
3799 p-=GetPixelChannels(image); /* go backward through pixel buffers */
3800 q-=GetPixelChannels(image);
3802 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3804 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3805 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3811 auth_view=DestroyCacheView(auth_view);
3812 virt_view=DestroyCacheView(virt_view);
3813 return(status ? (ssize_t) changed : -1);
3816 /* Apply a Morphology by calling one of the above low level primitive
3817 ** application functions. This function handles any iteration loops,
3818 ** composition or re-iteration of results, and compound morphology methods
3819 ** that is based on multiple low-level (staged) morphology methods.
3821 ** Basically this provides the complex glue between the requested morphology
3822 ** method and raw low-level implementation (above).
3824 MagickPrivate Image *MorphologyApply(const Image *image,
3825 const MorphologyMethod method, const ssize_t iterations,
3826 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3827 ExceptionInfo *exception)
3833 *curr_image, /* Image we are working with or iterating */
3834 *work_image, /* secondary image for primitive iteration */
3835 *save_image, /* saved image - for 'edge' method only */
3836 *rslt_image; /* resultant image - after multi-kernel handling */
3839 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3840 *norm_kernel, /* the current normal un-reflected kernel */
3841 *rflt_kernel, /* the current reflected kernel (if needed) */
3842 *this_kernel; /* the kernel being applied */
3845 primitive; /* the current morphology primitive being applied */
3848 rslt_compose; /* multi-kernel compose method for results to use */
3851 special, /* do we use a direct modify function? */
3852 verbose; /* verbose output of results */
3855 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3856 method_limit, /* maximum number of compound method iterations */
3857 kernel_number, /* Loop 2: the kernel number being applied */
3858 stage_loop, /* Loop 3: primitive loop for compound morphology */
3859 stage_limit, /* how many primitives are in this compound */
3860 kernel_loop, /* Loop 4: iterate the kernel over image */
3861 kernel_limit, /* number of times to iterate kernel */
3862 count, /* total count of primitive steps applied */
3863 kernel_changed, /* total count of changed using iterated kernel */
3864 method_changed; /* total count of changed over method iteration */
3867 changed; /* number pixels changed by last primitive operation */
3872 assert(image != (Image *) NULL);
3873 assert(image->signature == MagickSignature);
3874 assert(kernel != (KernelInfo *) NULL);
3875 assert(kernel->signature == MagickSignature);
3876 assert(exception != (ExceptionInfo *) NULL);
3877 assert(exception->signature == MagickSignature);
3879 count = 0; /* number of low-level morphology primitives performed */
3880 if ( iterations == 0 )
3881 return((Image *)NULL); /* null operation - nothing to do! */
3883 kernel_limit = (size_t) iterations;
3884 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3885 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3887 verbose = IsStringTrue(GetImageArtifact(image,"verbose"));
3889 /* initialise for cleanup */
3890 curr_image = (Image *) image;
3891 curr_compose = image->compose;
3892 (void) curr_compose;
3893 work_image = save_image = rslt_image = (Image *) NULL;
3894 reflected_kernel = (KernelInfo *) NULL;
3896 /* Initialize specific methods
3897 * + which loop should use the given iteratations
3898 * + how many primitives make up the compound morphology
3899 * + multi-kernel compose method to use (by default)
3901 method_limit = 1; /* just do method once, unless otherwise set */
3902 stage_limit = 1; /* assume method is not a compound */
3903 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3904 rslt_compose = compose; /* and we are composing multi-kernels as given */
3906 case SmoothMorphology: /* 4 primitive compound morphology */
3909 case OpenMorphology: /* 2 primitive compound morphology */
3910 case OpenIntensityMorphology:
3911 case TopHatMorphology:
3912 case CloseMorphology:
3913 case CloseIntensityMorphology:
3914 case BottomHatMorphology:
3915 case EdgeMorphology:
3918 case HitAndMissMorphology:
3919 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3921 case ThinningMorphology:
3922 case ThickenMorphology:
3923 method_limit = kernel_limit; /* iterate the whole method */
3924 kernel_limit = 1; /* do not do kernel iteration */
3926 case DistanceMorphology:
3927 case VoronoiMorphology:
3928 special = MagickTrue; /* use special direct primative */
3934 /* Apply special methods with special requirments
3935 ** For example, single run only, or post-processing requirements
3937 if ( special == MagickTrue )
3939 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3940 if (rslt_image == (Image *) NULL)
3942 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3945 changed = MorphologyPrimitiveDirect(rslt_image, method,
3948 if ( IfMagickTrue(verbose) )
3949 (void) (void) FormatLocaleFile(stderr,
3950 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3951 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3952 1.0,0.0,1.0, (double) changed);
3957 if ( method == VoronoiMorphology ) {
3958 /* Preserve the alpha channel of input image - but turned off */
3959 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3961 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3962 MagickTrue,0,0,exception);
3963 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3969 /* Handle user (caller) specified multi-kernel composition method */
3970 if ( compose != UndefinedCompositeOp )
3971 rslt_compose = compose; /* override default composition for method */
3972 if ( rslt_compose == UndefinedCompositeOp )
3973 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3975 /* Some methods require a reflected kernel to use with primitives.
3976 * Create the reflected kernel for those methods. */
3978 case CorrelateMorphology:
3979 case CloseMorphology:
3980 case CloseIntensityMorphology:
3981 case BottomHatMorphology:
3982 case SmoothMorphology:
3983 reflected_kernel = CloneKernelInfo(kernel);
3984 if (reflected_kernel == (KernelInfo *) NULL)
3986 RotateKernelInfo(reflected_kernel,180);
3992 /* Loops around more primitive morpholgy methods
3993 ** erose, dilate, open, close, smooth, edge, etc...
3995 /* Loop 1: iterate the compound method */
3998 while ( method_loop < method_limit && method_changed > 0 ) {
4002 /* Loop 2: iterate over each kernel in a multi-kernel list */
4003 norm_kernel = (KernelInfo *) kernel;
4004 this_kernel = (KernelInfo *) kernel;
4005 rflt_kernel = reflected_kernel;
4008 while ( norm_kernel != NULL ) {
4010 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
4011 stage_loop = 0; /* the compound morphology stage number */
4012 while ( stage_loop < stage_limit ) {
4013 stage_loop++; /* The stage of the compound morphology */
4015 /* Select primitive morphology for this stage of compound method */
4016 this_kernel = norm_kernel; /* default use unreflected kernel */
4017 primitive = method; /* Assume method is a primitive */
4019 case ErodeMorphology: /* just erode */
4020 case EdgeInMorphology: /* erode and image difference */
4021 primitive = ErodeMorphology;
4023 case DilateMorphology: /* just dilate */
4024 case EdgeOutMorphology: /* dilate and image difference */
4025 primitive = DilateMorphology;
4027 case OpenMorphology: /* erode then dialate */
4028 case TopHatMorphology: /* open and image difference */
4029 primitive = ErodeMorphology;
4030 if ( stage_loop == 2 )
4031 primitive = DilateMorphology;
4033 case OpenIntensityMorphology:
4034 primitive = ErodeIntensityMorphology;
4035 if ( stage_loop == 2 )
4036 primitive = DilateIntensityMorphology;
4038 case CloseMorphology: /* dilate, then erode */
4039 case BottomHatMorphology: /* close and image difference */
4040 this_kernel = rflt_kernel; /* use the reflected kernel */
4041 primitive = DilateMorphology;
4042 if ( stage_loop == 2 )
4043 primitive = ErodeMorphology;
4045 case CloseIntensityMorphology:
4046 this_kernel = rflt_kernel; /* use the reflected kernel */
4047 primitive = DilateIntensityMorphology;
4048 if ( stage_loop == 2 )
4049 primitive = ErodeIntensityMorphology;
4051 case SmoothMorphology: /* open, close */
4052 switch ( stage_loop ) {
4053 case 1: /* start an open method, which starts with Erode */
4054 primitive = ErodeMorphology;
4056 case 2: /* now Dilate the Erode */
4057 primitive = DilateMorphology;
4059 case 3: /* Reflect kernel a close */
4060 this_kernel = rflt_kernel; /* use the reflected kernel */
4061 primitive = DilateMorphology;
4063 case 4: /* Finish the Close */
4064 this_kernel = rflt_kernel; /* use the reflected kernel */
4065 primitive = ErodeMorphology;
4069 case EdgeMorphology: /* dilate and erode difference */
4070 primitive = DilateMorphology;
4071 if ( stage_loop == 2 ) {
4072 save_image = curr_image; /* save the image difference */
4073 curr_image = (Image *) image;
4074 primitive = ErodeMorphology;
4077 case CorrelateMorphology:
4078 /* A Correlation is a Convolution with a reflected kernel.
4079 ** However a Convolution is a weighted sum using a reflected
4080 ** kernel. It may seem stange to convert a Correlation into a
4081 ** Convolution as the Correlation is the simplier method, but
4082 ** Convolution is much more commonly used, and it makes sense to
4083 ** implement it directly so as to avoid the need to duplicate the
4084 ** kernel when it is not required (which is typically the
4087 this_kernel = rflt_kernel; /* use the reflected kernel */
4088 primitive = ConvolveMorphology;
4093 assert( this_kernel != (KernelInfo *) NULL );
4095 /* Extra information for debugging compound operations */
4096 if ( IfMagickTrue(verbose) ) {
4097 if ( stage_limit > 1 )
4098 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
4099 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
4100 method_loop,(double) stage_loop);
4101 else if ( primitive != method )
4102 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
4103 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
4109 /* Loop 4: Iterate the kernel with primitive */
4113 while ( kernel_loop < kernel_limit && changed > 0 ) {
4114 kernel_loop++; /* the iteration of this kernel */
4116 /* Create a clone as the destination image, if not yet defined */
4117 if ( work_image == (Image *) NULL )
4119 work_image=CloneImage(image,0,0,MagickTrue,exception);
4120 if (work_image == (Image *) NULL)
4122 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
4124 /* work_image->type=image->type; ??? */
4127 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4129 changed = MorphologyPrimitive(curr_image, work_image, primitive,
4130 this_kernel, bias, exception);
4132 if ( IfMagickTrue(verbose) ) {
4133 if ( kernel_loop > 1 )
4134 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4135 (void) (void) FormatLocaleFile(stderr,
4136 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4137 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4138 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4139 (double) (method_loop+kernel_loop-1),(double) kernel_number,
4140 (double) count,(double) changed);
4144 kernel_changed += changed;
4145 method_changed += changed;
4147 /* prepare next loop */
4148 { Image *tmp = work_image; /* swap images for iteration */
4149 work_image = curr_image;
4152 if ( work_image == image )
4153 work_image = (Image *) NULL; /* replace input 'image' */
4155 } /* End Loop 4: Iterate the kernel with primitive */
4157 if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed )
4158 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
4159 if ( IfMagickTrue(verbose) && stage_loop < stage_limit )
4160 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4163 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4164 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
4165 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
4166 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
4167 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
4170 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4172 /* Final Post-processing for some Compound Methods
4174 ** The removal of any 'Sync' channel flag in the Image Compositon
4175 ** below ensures the methematical compose method is applied in a
4176 ** purely mathematical way, and only to the selected channels.
4177 ** Turn off SVG composition 'alpha blending'.
4180 case EdgeOutMorphology:
4181 case EdgeInMorphology:
4182 case TopHatMorphology:
4183 case BottomHatMorphology:
4184 if ( IfMagickTrue(verbose) )
4185 (void) FormatLocaleFile(stderr,
4186 "\n%s: Difference with original image",CommandOptionToMnemonic(
4187 MagickMorphologyOptions, method) );
4188 (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
4189 MagickTrue,0,0,exception);
4191 case EdgeMorphology:
4192 if ( IfMagickTrue(verbose) )
4193 (void) FormatLocaleFile(stderr,
4194 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
4195 MagickMorphologyOptions, method) );
4196 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
4197 MagickTrue,0,0,exception);
4198 save_image = DestroyImage(save_image); /* finished with save image */
4204 /* multi-kernel handling: re-iterate, or compose results */
4205 if ( kernel->next == (KernelInfo *) NULL )
4206 rslt_image = curr_image; /* just return the resulting image */
4207 else if ( rslt_compose == NoCompositeOp )
4208 { if ( IfMagickTrue(verbose) ) {
4209 if ( this_kernel->next != (KernelInfo *) NULL )
4210 (void) FormatLocaleFile(stderr, " (re-iterate)");
4212 (void) FormatLocaleFile(stderr, " (done)");
4214 rslt_image = curr_image; /* return result, and re-iterate */
4216 else if ( rslt_image == (Image *) NULL)
4217 { if ( IfMagickTrue(verbose) )
4218 (void) FormatLocaleFile(stderr, " (save for compose)");
4219 rslt_image = curr_image;
4220 curr_image = (Image *) image; /* continue with original image */
4223 { /* Add the new 'current' result to the composition
4225 ** The removal of any 'Sync' channel flag in the Image Compositon
4226 ** below ensures the methematical compose method is applied in a
4227 ** purely mathematical way, and only to the selected channels.
4228 ** IE: Turn off SVG composition 'alpha blending'.
4230 if ( IfMagickTrue(verbose) )
4231 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4232 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4233 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4235 curr_image = DestroyImage(curr_image);
4236 curr_image = (Image *) image; /* continue with original image */
4238 if ( IfMagickTrue(verbose) )
4239 (void) FormatLocaleFile(stderr, "\n");
4241 /* loop to the next kernel in a multi-kernel list */
4242 norm_kernel = norm_kernel->next;
4243 if ( rflt_kernel != (KernelInfo *) NULL )
4244 rflt_kernel = rflt_kernel->next;
4246 } /* End Loop 2: Loop over each kernel */
4248 } /* End Loop 1: compound method interation */
4252 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4254 if ( curr_image == rslt_image )
4255 curr_image = (Image *) NULL;
4256 if ( rslt_image != (Image *) NULL )
4257 rslt_image = DestroyImage(rslt_image);
4259 if ( curr_image == rslt_image || curr_image == image )
4260 curr_image = (Image *) NULL;
4261 if ( curr_image != (Image *) NULL )
4262 curr_image = DestroyImage(curr_image);
4263 if ( work_image != (Image *) NULL )
4264 work_image = DestroyImage(work_image);
4265 if ( save_image != (Image *) NULL )
4266 save_image = DestroyImage(save_image);
4267 if ( reflected_kernel != (KernelInfo *) NULL )
4268 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4274 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4278 % M o r p h o l o g y I m a g e %
4282 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4284 % MorphologyImage() applies a user supplied kernel to the image
4285 % according to the given mophology method.
4287 % This function applies any and all user defined settings before calling
4288 % the above internal function MorphologyApply().
4290 % User defined settings include...
4291 % * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4292 % * Kernel Scale/normalize settings ("-define convolve:scale=??")
4293 % This can also includes the addition of a scaled unity kernel.
4294 % * Show Kernel being applied ("-define showkernel=1")
4296 % Other operators that do not want user supplied options interfering,
4297 % especially "convolve:bias" and "showkernel" should use MorphologyApply()
4300 % The format of the MorphologyImage method is:
4302 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4303 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4305 % A description of each parameter follows:
4307 % o image: the image.
4309 % o method: the morphology method to be applied.
4311 % o iterations: apply the operation this many times (or no change).
4312 % A value of -1 means loop until no change found.
4313 % How this is applied may depend on the morphology method.
4314 % Typically this is a value of 1.
4316 % o kernel: An array of double representing the morphology kernel.
4317 % Warning: kernel may be normalized for the Convolve method.
4319 % o exception: return any errors or warnings in this structure.
4322 MagickExport Image *MorphologyImage(const Image *image,
4323 const MorphologyMethod method,const ssize_t iterations,
4324 const KernelInfo *kernel,ExceptionInfo *exception)
4338 curr_kernel = (KernelInfo *) kernel;
4340 compose = UndefinedCompositeOp; /* use default for method */
4342 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4343 * This is done BEFORE the ShowKernelInfo() function is called so that
4344 * users can see the results of the 'option:convolve:scale' option.
4346 if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4350 /* Get the bias value as it will be needed */
4351 artifact = GetImageArtifact(image,"convolve:bias");
4352 if ( artifact != (const char *) NULL) {
4353 if (IfMagickFalse(IsGeometry(artifact)))
4354 (void) ThrowMagickException(exception,GetMagickModule(),
4355 OptionWarning,"InvalidSetting","'%s' '%s'",
4356 "convolve:bias",artifact);
4358 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4361 /* Scale kernel according to user wishes */
4362 artifact = GetImageArtifact(image,"convolve:scale");
4363 if ( artifact != (const char *)NULL ) {
4364 if (IfMagickFalse(IsGeometry(artifact)))
4365 (void) ThrowMagickException(exception,GetMagickModule(),
4366 OptionWarning,"InvalidSetting","'%s' '%s'",
4367 "convolve:scale",artifact);
4369 if ( curr_kernel == kernel )
4370 curr_kernel = CloneKernelInfo(kernel);
4371 if (curr_kernel == (KernelInfo *) NULL)
4372 return((Image *) NULL);
4373 ScaleGeometryKernelInfo(curr_kernel, artifact);
4378 /* display the (normalized) kernel via stderr */
4379 if ( IfStringTrue(GetImageArtifact(image,"showkernel"))
4380 || IfStringTrue(GetImageArtifact(image,"convolve:showkernel"))
4381 || IfStringTrue(GetImageArtifact(image,"morphology:showkernel")) )
4382 ShowKernelInfo(curr_kernel);
4384 /* Override the default handling of multi-kernel morphology results
4385 * If 'Undefined' use the default method
4386 * If 'None' (default for 'Convolve') re-iterate previous result
4387 * Otherwise merge resulting images using compose method given.
4388 * Default for 'HitAndMiss' is 'Lighten'.
4395 artifact = GetImageArtifact(image,"morphology:compose");
4396 if ( artifact != (const char *) NULL) {
4397 parse=ParseCommandOption(MagickComposeOptions,
4398 MagickFalse,artifact);
4400 (void) ThrowMagickException(exception,GetMagickModule(),
4401 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4402 "morphology:compose",artifact);
4404 compose=(CompositeOperator)parse;
4407 /* Apply the Morphology */
4408 morphology_image = MorphologyApply(image,method,iterations,
4409 curr_kernel,compose,bias,exception);
4411 /* Cleanup and Exit */
4412 if ( curr_kernel != kernel )
4413 curr_kernel=DestroyKernelInfo(curr_kernel);
4414 return(morphology_image);
4418 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4422 + R o t a t e K e r n e l I n f o %
4426 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4428 % RotateKernelInfo() rotates the kernel by the angle given.
4430 % Currently it is restricted to 90 degree angles, of either 1D kernels
4431 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4432 % It will ignore usless rotations for specific 'named' built-in kernels.
4434 % The format of the RotateKernelInfo method is:
4436 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4438 % A description of each parameter follows:
4440 % o kernel: the Morphology/Convolution kernel
4442 % o angle: angle to rotate in degrees
4444 % This function is currently internal to this module only, but can be exported
4445 % to other modules if needed.
4447 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4449 /* angle the lower kernels first */
4450 if ( kernel->next != (KernelInfo *) NULL)
4451 RotateKernelInfo(kernel->next, angle);
4453 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4455 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4458 /* Modulus the angle */
4459 angle = fmod(angle, 360.0);
4463 if ( 337.5 < angle || angle <= 22.5 )
4464 return; /* Near zero angle - no change! - At least not at this time */
4466 /* Handle special cases */
4467 switch (kernel->type) {
4468 /* These built-in kernels are cylindrical kernels, rotating is useless */
4469 case GaussianKernel:
4474 case LaplacianKernel:
4475 case ChebyshevKernel:
4476 case ManhattanKernel:
4477 case EuclideanKernel:
4480 /* These may be rotatable at non-90 angles in the future */
4481 /* but simply rotating them in multiples of 90 degrees is useless */
4488 /* These only allows a +/-90 degree rotation (by transpose) */
4489 /* A 180 degree rotation is useless */
4491 if ( 135.0 < angle && angle <= 225.0 )
4493 if ( 225.0 < angle && angle <= 315.0 )
4500 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4501 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4503 if ( kernel->width == 3 && kernel->height == 3 )
4504 { /* Rotate a 3x3 square by 45 degree angle */
4505 double t = kernel->values[0];
4506 kernel->values[0] = kernel->values[3];
4507 kernel->values[3] = kernel->values[6];
4508 kernel->values[6] = kernel->values[7];
4509 kernel->values[7] = kernel->values[8];
4510 kernel->values[8] = kernel->values[5];
4511 kernel->values[5] = kernel->values[2];
4512 kernel->values[2] = kernel->values[1];
4513 kernel->values[1] = t;
4514 /* rotate non-centered origin */
4515 if ( kernel->x != 1 || kernel->y != 1 ) {
4517 x = (ssize_t) kernel->x-1;
4518 y = (ssize_t) kernel->y-1;
4519 if ( x == y ) x = 0;
4520 else if ( x == 0 ) x = -y;
4521 else if ( x == -y ) y = 0;
4522 else if ( y == 0 ) y = x;
4523 kernel->x = (ssize_t) x+1;
4524 kernel->y = (ssize_t) y+1;
4526 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4527 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4530 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4532 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4534 if ( kernel->width == 1 || kernel->height == 1 )
4535 { /* Do a transpose of a 1 dimensional kernel,
4536 ** which results in a fast 90 degree rotation of some type.
4540 t = (ssize_t) kernel->width;
4541 kernel->width = kernel->height;
4542 kernel->height = (size_t) t;
4544 kernel->x = kernel->y;
4546 if ( kernel->width == 1 ) {
4547 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4548 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4550 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4551 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4554 else if ( kernel->width == kernel->height )
4555 { /* Rotate a square array of values by 90 degrees */
4559 register MagickRealType
4563 for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--)
4564 for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--)
4565 { t = k[i+j*kernel->width];
4566 k[i+j*kernel->width] = k[j+x*kernel->width];
4567 k[j+x*kernel->width] = k[x+y*kernel->width];
4568 k[x+y*kernel->width] = k[y+i*kernel->width];
4569 k[y+i*kernel->width] = t;
4572 /* rotate the origin - relative to center of array */
4573 { register ssize_t x,y;
4574 x = (ssize_t) (kernel->x*2-kernel->width+1);
4575 y = (ssize_t) (kernel->y*2-kernel->height+1);
4576 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4577 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4579 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4580 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4583 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4585 if ( 135.0 < angle && angle <= 225.0 )
4587 /* For a 180 degree rotation - also know as a reflection
4588 * This is actually a very very common operation!
4589 * Basically all that is needed is a reversal of the kernel data!
4590 * And a reflection of the origon
4595 register MagickRealType
4603 j=(ssize_t) (kernel->width*kernel->height-1);
4604 for (i=0; i < j; i++, j--)
4605 t=k[i], k[i]=k[j], k[j]=t;
4607 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4608 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4609 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4610 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4612 /* At this point angle should at least between -45 (315) and +45 degrees
4613 * In the future some form of non-orthogonal angled rotates could be
4614 * performed here, posibily with a linear kernel restriction.
4621 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4625 % S c a l e G e o m e t r y K e r n e l I n f o %
4629 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4631 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4632 % provided as a "-set option:convolve:scale {geometry}" user setting,
4633 % and modifies the kernel according to the parsed arguments of that setting.
4635 % The first argument (and any normalization flags) are passed to
4636 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4637 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4638 % into the scaled/normalized kernel.
4640 % The format of the ScaleGeometryKernelInfo method is:
4642 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4643 % const double scaling_factor,const MagickStatusType normalize_flags)
4645 % A description of each parameter follows:
4647 % o kernel: the Morphology/Convolution kernel to modify
4650 % The geometry string to parse, typically from the user provided
4651 % "-set option:convolve:scale {geometry}" setting.
4654 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4655 const char *geometry)
4664 SetGeometryInfo(&args);
4665 flags = ParseGeometry(geometry, &args);
4668 /* For Debugging Geometry Input */
4669 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4670 flags, args.rho, args.sigma, args.xi, args.psi );
4673 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4674 args.rho *= 0.01, args.sigma *= 0.01;
4676 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4678 if ( (flags & SigmaValue) == 0 )
4681 /* Scale/Normalize the input kernel */
4682 ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4684 /* Add Unity Kernel, for blending with original */
4685 if ( (flags & SigmaValue) != 0 )
4686 UnityAddKernelInfo(kernel, args.sigma);
4691 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4695 % S c a l e K e r n e l I n f o %
4699 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4701 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4702 % without normalization of the sum of the kernel values (as per given flags).
4704 % By default (no flags given) the values within the kernel is scaled
4705 % directly using given scaling factor without change.
4707 % If either of the two 'normalize_flags' are given the kernel will first be
4708 % normalized and then further scaled by the scaling factor value given.
4710 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4711 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4712 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4713 % non-HDRI versions of IM this may cause images to have any negative results
4714 % clipped, unless some 'bias' is used.
4716 % More specifically. Kernels which only contain positive values (such as a
4717 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4718 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4720 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4721 % the kernel will be scaled by the absolute of the sum of kernel values, so
4722 % that it will generally fall within the +/- 1.0 range.
4724 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4725 % will be scaled by just the sum of the postive values, so that its output
4726 % range will again fall into the +/- 1.0 range.
4728 % For special kernels designed for locating shapes using 'Correlate', (often
4729 % only containing +1 and -1 values, representing foreground/brackground
4730 % matching) a special normalization method is provided to scale the positive
4731 % values separately to those of the negative values, so the kernel will be
4732 % forced to become a zero-sum kernel better suited to such searches.
4734 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4735 % attributes within the kernel structure have been correctly set during the
4738 % NOTE: The values used for 'normalize_flags' have been selected specifically
4739 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4740 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4742 % The format of the ScaleKernelInfo method is:
4744 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4745 % const MagickStatusType normalize_flags )
4747 % A description of each parameter follows:
4749 % o kernel: the Morphology/Convolution kernel
4752 % multiply all values (after normalization) by this factor if not
4753 % zero. If the kernel is normalized regardless of any flags.
4755 % o normalize_flags:
4756 % GeometryFlags defining normalization method to use.
4757 % specifically: NormalizeValue, CorrelateNormalizeValue,
4758 % and/or PercentValue
4761 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4762 const double scaling_factor,const GeometryFlags normalize_flags)
4771 /* do the other kernels in a multi-kernel list first */
4772 if ( kernel->next != (KernelInfo *) NULL)
4773 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4775 /* Normalization of Kernel */
4777 if ( (normalize_flags&NormalizeValue) != 0 ) {
4778 if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4779 /* non-zero-summing kernel (generally positive) */
4780 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4782 /* zero-summing kernel */
4783 pos_scale = kernel->positive_range;
4785 /* Force kernel into a normalized zero-summing kernel */
4786 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4787 pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4788 ? kernel->positive_range : 1.0;
4789 neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4790 ? -kernel->negative_range : 1.0;
4793 neg_scale = pos_scale;
4795 /* finialize scaling_factor for positive and negative components */
4796 pos_scale = scaling_factor/pos_scale;
4797 neg_scale = scaling_factor/neg_scale;
4799 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4800 if ( ! IsNan(kernel->values[i]) )
4801 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4803 /* convolution output range */
4804 kernel->positive_range *= pos_scale;
4805 kernel->negative_range *= neg_scale;
4806 /* maximum and minimum values in kernel */
4807 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4808 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4810 /* swap kernel settings if user's scaling factor is negative */
4811 if ( scaling_factor < MagickEpsilon ) {
4813 t = kernel->positive_range;
4814 kernel->positive_range = kernel->negative_range;
4815 kernel->negative_range = t;
4816 t = kernel->maximum;
4817 kernel->maximum = kernel->minimum;
4818 kernel->minimum = 1;
4825 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4829 % S h o w K e r n e l I n f o %
4833 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4835 % ShowKernelInfo() outputs the details of the given kernel defination to
4836 % standard error, generally due to a users 'showkernel' option request.
4838 % The format of the ShowKernel method is:
4840 % void ShowKernelInfo(const KernelInfo *kernel)
4842 % A description of each parameter follows:
4844 % o kernel: the Morphology/Convolution kernel
4847 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4855 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4857 (void) FormatLocaleFile(stderr, "Kernel");
4858 if ( kernel->next != (KernelInfo *) NULL )
4859 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4860 (void) FormatLocaleFile(stderr, " \"%s",
4861 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4862 if ( fabs(k->angle) >= MagickEpsilon )
4863 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4864 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4865 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4866 (void) FormatLocaleFile(stderr,
4867 " with values from %.*lg to %.*lg\n",
4868 GetMagickPrecision(), k->minimum,
4869 GetMagickPrecision(), k->maximum);
4870 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4871 GetMagickPrecision(), k->negative_range,
4872 GetMagickPrecision(), k->positive_range);
4873 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4874 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4875 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4876 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4878 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4879 GetMagickPrecision(), k->positive_range+k->negative_range);
4880 for (i=v=0; v < k->height; v++) {
4881 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4882 for (u=0; u < k->width; u++, i++)
4883 if ( IsNan(k->values[i]) )
4884 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4886 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4887 GetMagickPrecision(), (double) k->values[i]);
4888 (void) FormatLocaleFile(stderr,"\n");
4894 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4898 % U n i t y A d d K e r n a l I n f o %
4902 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4904 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4905 % to the given pre-scaled and normalized Kernel. This in effect adds that
4906 % amount of the original image into the resulting convolution kernel. This
4907 % value is usually provided by the user as a percentage value in the
4908 % 'convolve:scale' setting.
4910 % The resulting effect is to convert the defined kernels into blended
4911 % soft-blurs, unsharp kernels or into sharpening kernels.
4913 % The format of the UnityAdditionKernelInfo method is:
4915 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4917 % A description of each parameter follows:
4919 % o kernel: the Morphology/Convolution kernel
4922 % scaling factor for the unity kernel to be added to
4926 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4929 /* do the other kernels in a multi-kernel list first */
4930 if ( kernel->next != (KernelInfo *) NULL)
4931 UnityAddKernelInfo(kernel->next, scale);
4933 /* Add the scaled unity kernel to the existing kernel */
4934 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4935 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4941 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4945 % Z e r o K e r n e l N a n s %
4949 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4951 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4952 % the kernel with a zero value. This is typically done when the kernel will
4953 % be used in special hardware (GPU) convolution processors, to simply
4956 % The format of the ZeroKernelNans method is:
4958 % void ZeroKernelNans (KernelInfo *kernel)
4960 % A description of each parameter follows:
4962 % o kernel: the Morphology/Convolution kernel
4965 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4970 /* do the other kernels in a multi-kernel list first */
4971 if ( kernel->next != (KernelInfo *) NULL)
4972 ZeroKernelNans(kernel->next);
4974 for (i=0; i < (kernel->width*kernel->height); i++)
4975 if ( IsNan(kernel->values[i]) )
4976 kernel->values[i] = 0.0;