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/monitor-private.h"
68 #include "MagickCore/morphology.h"
69 #include "MagickCore/morphology-private.h"
70 #include "MagickCore/option.h"
71 #include "MagickCore/pixel-accessor.h"
72 #include "MagickCore/prepress.h"
73 #include "MagickCore/quantize.h"
74 #include "MagickCore/registry.h"
75 #include "MagickCore/semaphore.h"
76 #include "MagickCore/splay-tree.h"
77 #include "MagickCore/statistic.h"
78 #include "MagickCore/string_.h"
79 #include "MagickCore/string-private.h"
80 #include "MagickCore/token.h"
81 #include "MagickCore/utility.h"
82 #include "MagickCore/utility-private.h"
86 ** The following test is for special floating point numbers of value NaN (not
87 ** a number), that may be used within a Kernel Definition. NaN's are defined
88 ** as part of the IEEE standard for floating point number representation.
90 ** These are used as a Kernel value to mean that this kernel position is not
91 ** part of the kernel neighbourhood for convolution or morphology processing,
92 ** and thus should be ignored. This allows the use of 'shaped' kernels.
94 ** The special properity that two NaN's are never equal, even if they are from
95 ** the same variable allow you to test if a value is special NaN value.
97 ** This macro IsNaN() is thus is only true if the value given is NaN.
99 #define IsNan(a) ((a)!=(a))
102 Other global definitions used by module.
104 static inline double MagickMin(const double x,const double y)
106 return( x < y ? x : y);
108 static inline double MagickMax(const double x,const double y)
110 return( x > y ? x : y);
112 #define Minimize(assign,value) assign=MagickMin(assign,value)
113 #define Maximize(assign,value) assign=MagickMax(assign,value)
115 /* Currently these are only internal to this module */
117 CalcKernelMetaData(KernelInfo *),
118 ExpandMirrorKernelInfo(KernelInfo *),
119 ExpandRotateKernelInfo(KernelInfo *, const double),
120 RotateKernelInfo(KernelInfo *, double);
123 /* Quick function to find last kernel in a kernel list */
124 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
126 while (kernel->next != (KernelInfo *) NULL)
127 kernel = kernel->next;
132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
136 % A c q u i r e K e r n e l I n f o %
140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
142 % AcquireKernelInfo() takes the given string (generally supplied by the
143 % user) and converts it into a Morphology/Convolution Kernel. This allows
144 % users to specify a kernel from a number of pre-defined kernels, or to fully
145 % specify their own kernel for a specific Convolution or Morphology
148 % The kernel so generated can be any rectangular array of floating point
149 % values (doubles) with the 'control point' or 'pixel being affected'
150 % anywhere within that array of values.
152 % Previously IM was restricted to a square of odd size using the exact
153 % center as origin, this is no longer the case, and any rectangular kernel
154 % with any value being declared the origin. This in turn allows the use of
155 % highly asymmetrical kernels.
157 % The floating point values in the kernel can also include a special value
158 % known as 'nan' or 'not a number' to indicate that this value is not part
159 % of the kernel array. This allows you to shaped the kernel within its
160 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
161 % shape. However at least one non-nan value must be provided for correct
162 % working of a kernel.
164 % The returned kernel should be freed using the DestroyKernelInfo() when you
165 % are finished with it. Do not free this memory yourself.
167 % Input kernel defintion strings can consist of any of three types.
170 % Select from one of the built in kernels, using the name and
171 % geometry arguments supplied. See AcquireKernelBuiltIn()
173 % "WxH[+X+Y][@><]:num, num, num ..."
174 % a kernel of size W by H, with W*H floating point numbers following.
175 % the 'center' can be optionally be defined at +X+Y (such that +0+0
176 % is top left corner). If not defined the pixel in the center, for
177 % odd sizes, or to the immediate top or left of center for even sizes
178 % is automatically selected.
180 % "num, num, num, num, ..."
181 % list of floating point numbers defining an 'old style' odd sized
182 % square kernel. At least 9 values should be provided for a 3x3
183 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
184 % Values can be space or comma separated. This is not recommended.
186 % You can define a 'list of kernels' which can be used by some morphology
187 % operators A list is defined as a semi-colon separated list kernels.
189 % " kernel ; kernel ; kernel ; "
191 % Any extra ';' characters, at start, end or between kernel defintions are
194 % The special flags will expand a single kernel, into a list of rotated
195 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
196 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
197 % The '<' also exands using 90-degree rotates, but giving a 180-degree
198 % reflected kernel before the +/- 90-degree rotations, which can be important
199 % for Thinning operations.
201 % Note that 'name' kernels will start with an alphabetic character while the
202 % new kernel specification has a ':' character in its specification string.
203 % If neither is the case, it is assumed an old style of a simple list of
204 % numbers generating a odd-sized square kernel has been given.
206 % The format of the AcquireKernal method is:
208 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
210 % A description of each parameter follows:
212 % o kernel_string: the Morphology/Convolution kernel wanted.
216 /* This was separated so that it could be used as a separate
217 ** array input handling function, such as for -color-matrix
219 static KernelInfo *ParseKernelArray(const char *kernel_string)
225 token[MaxTextExtent];
235 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
243 kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
244 if (kernel == (KernelInfo *)NULL)
246 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
247 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
248 kernel->negative_range = kernel->positive_range = 0.0;
249 kernel->type = UserDefinedKernel;
250 kernel->next = (KernelInfo *) NULL;
251 kernel->signature = MagickSignature;
252 if (kernel_string == (const char *) NULL)
255 /* find end of this specific kernel definition string */
256 end = strchr(kernel_string, ';');
257 if ( end == (char *) NULL )
258 end = strchr(kernel_string, '\0');
260 /* clear flags - for Expanding kernel lists thorugh rotations */
263 /* Has a ':' in argument - New user kernel specification
264 FUTURE: this split on ':' could be done by StringToken()
266 p = strchr(kernel_string, ':');
267 if ( p != (char *) NULL && p < end)
269 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
270 memcpy(token, kernel_string, (size_t) (p-kernel_string));
271 token[p-kernel_string] = '\0';
272 SetGeometryInfo(&args);
273 flags = ParseGeometry(token, &args);
275 /* Size handling and checks of geometry settings */
276 if ( (flags & WidthValue) == 0 ) /* if no width then */
277 args.rho = args.sigma; /* then width = height */
278 if ( args.rho < 1.0 ) /* if width too small */
279 args.rho = 1.0; /* then width = 1 */
280 if ( args.sigma < 1.0 ) /* if height too small */
281 args.sigma = args.rho; /* then height = width */
282 kernel->width = (size_t)args.rho;
283 kernel->height = (size_t)args.sigma;
285 /* Offset Handling and Checks */
286 if ( args.xi < 0.0 || args.psi < 0.0 )
287 return(DestroyKernelInfo(kernel));
288 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
289 : (ssize_t) (kernel->width-1)/2;
290 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
291 : (ssize_t) (kernel->height-1)/2;
292 if ( kernel->x >= (ssize_t) kernel->width ||
293 kernel->y >= (ssize_t) kernel->height )
294 return(DestroyKernelInfo(kernel));
296 p++; /* advance beyond the ':' */
299 { /* ELSE - Old old specification, forming odd-square kernel */
300 /* count up number of values given */
301 p=(const char *) kernel_string;
302 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
303 p++; /* ignore "'" chars for convolve filter usage - Cristy */
304 for (i=0; p < end; i++)
306 GetMagickToken(p,&p,token);
308 GetMagickToken(p,&p,token);
310 /* set the size of the kernel - old sized square */
311 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
312 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
313 p=(const char *) kernel_string;
314 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
315 p++; /* ignore "'" chars for convolve filter usage - Cristy */
318 /* Read in the kernel values from rest of input string argument */
319 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
320 kernel->height*sizeof(*kernel->values));
321 if (kernel->values == (MagickRealType *) NULL)
322 return(DestroyKernelInfo(kernel));
323 kernel->minimum = +MagickHuge;
324 kernel->maximum = -MagickHuge;
325 kernel->negative_range = kernel->positive_range = 0.0;
326 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
328 GetMagickToken(p,&p,token);
330 GetMagickToken(p,&p,token);
331 if ( LocaleCompare("nan",token) == 0
332 || LocaleCompare("-",token) == 0 ) {
333 kernel->values[i] = nan; /* this value is not part of neighbourhood */
336 kernel->values[i] = StringToDouble(token,(char **) NULL);
337 ( kernel->values[i] < 0)
338 ? ( kernel->negative_range += kernel->values[i] )
339 : ( kernel->positive_range += kernel->values[i] );
340 Minimize(kernel->minimum, kernel->values[i]);
341 Maximize(kernel->maximum, kernel->values[i]);
345 /* sanity check -- no more values in kernel definition */
346 GetMagickToken(p,&p,token);
347 if ( *token != '\0' && *token != ';' && *token != '\'' )
348 return(DestroyKernelInfo(kernel));
351 /* this was the old method of handling a incomplete kernel */
352 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
353 Minimize(kernel->minimum, kernel->values[i]);
354 Maximize(kernel->maximum, kernel->values[i]);
355 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
356 kernel->values[i]=0.0;
359 /* Number of values for kernel was not enough - Report Error */
360 if ( i < (ssize_t) (kernel->width*kernel->height) )
361 return(DestroyKernelInfo(kernel));
364 /* check that we recieved at least one real (non-nan) value! */
365 if ( kernel->minimum == MagickHuge )
366 return(DestroyKernelInfo(kernel));
368 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
369 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
370 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
371 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
372 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
373 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
378 static KernelInfo *ParseKernelName(const char *kernel_string)
381 token[MaxTextExtent];
399 /* Parse special 'named' kernel */
400 GetMagickToken(kernel_string,&p,token);
401 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
402 if ( type < 0 || type == UserDefinedKernel )
403 return((KernelInfo *)NULL); /* not a valid named kernel */
405 while (((isspace((int) ((unsigned char) *p)) != 0) ||
406 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
409 end = strchr(p, ';'); /* end of this kernel defintion */
410 if ( end == (char *) NULL )
411 end = strchr(p, '\0');
413 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
414 memcpy(token, p, (size_t) (end-p));
416 SetGeometryInfo(&args);
417 flags = ParseGeometry(token, &args);
420 /* For Debugging Geometry Input */
421 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
422 flags, args.rho, args.sigma, args.xi, args.psi );
425 /* special handling of missing values in input string */
427 /* Shape Kernel Defaults */
429 if ( (flags & WidthValue) == 0 )
430 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
438 if ( (flags & HeightValue) == 0 )
439 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
442 if ( (flags & XValue) == 0 )
443 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
445 case RectangleKernel: /* Rectangle - set size defaults */
446 if ( (flags & WidthValue) == 0 ) /* if no width then */
447 args.rho = args.sigma; /* then width = height */
448 if ( args.rho < 1.0 ) /* if width too small */
449 args.rho = 3; /* then width = 3 */
450 if ( args.sigma < 1.0 ) /* if height too small */
451 args.sigma = args.rho; /* then height = width */
452 if ( (flags & XValue) == 0 ) /* center offset if not defined */
453 args.xi = (double)(((ssize_t)args.rho-1)/2);
454 if ( (flags & YValue) == 0 )
455 args.psi = (double)(((ssize_t)args.sigma-1)/2);
457 /* Distance Kernel Defaults */
458 case ChebyshevKernel:
459 case ManhattanKernel:
460 case OctagonalKernel:
461 case EuclideanKernel:
462 if ( (flags & HeightValue) == 0 ) /* no distance scale */
463 args.sigma = 100.0; /* default distance scaling */
464 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
465 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
466 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
467 args.sigma *= QuantumRange/100.0; /* percentage of color range */
473 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
474 if ( kernel == (KernelInfo *) NULL )
477 /* global expand to rotated kernel list - only for single kernels */
478 if ( kernel->next == (KernelInfo *) NULL ) {
479 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
480 ExpandRotateKernelInfo(kernel, 45.0);
481 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
482 ExpandRotateKernelInfo(kernel, 90.0);
483 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
484 ExpandMirrorKernelInfo(kernel);
490 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
498 token[MaxTextExtent];
506 if (kernel_string == (const char *) NULL)
507 return(ParseKernelArray(kernel_string));
512 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
514 /* ignore extra or multiple ';' kernel separators */
515 if ( *token != ';' ) {
517 /* tokens starting with alpha is a Named kernel */
518 if (isalpha((int) *token) != 0)
519 new_kernel = ParseKernelName(p);
520 else /* otherwise a user defined kernel array */
521 new_kernel = ParseKernelArray(p);
523 /* Error handling -- this is not proper error handling! */
524 if ( new_kernel == (KernelInfo *) NULL ) {
525 (void) FormatLocaleFile(stderr, "Failed to parse kernel number #%.20g\n",
526 (double) kernel_number);
527 if ( kernel != (KernelInfo *) NULL )
528 kernel=DestroyKernelInfo(kernel);
529 return((KernelInfo *) NULL);
532 /* initialise or append the kernel list */
533 if ( kernel == (KernelInfo *) NULL )
536 LastKernelInfo(kernel)->next = new_kernel;
539 /* look for the next kernel in list */
541 if ( p == (char *) NULL )
551 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
555 % A c q u i r e K e r n e l B u i l t I n %
559 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
561 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
562 % kernels used for special purposes such as gaussian blurring, skeleton
563 % pruning, and edge distance determination.
565 % They take a KernelType, and a set of geometry style arguments, which were
566 % typically decoded from a user supplied string, or from a more complex
567 % Morphology Method that was requested.
569 % The format of the AcquireKernalBuiltIn method is:
571 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
572 % const GeometryInfo args)
574 % A description of each parameter follows:
576 % o type: the pre-defined type of kernel wanted
578 % o args: arguments defining or modifying the kernel
580 % Convolution Kernels
583 % The a No-Op or Scaling single element kernel.
585 % Gaussian:{radius},{sigma}
586 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
587 % The sigma for the curve is required. The resulting kernel is
590 % If 'sigma' is zero, you get a single pixel on a field of zeros.
592 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
593 % the final size of the resulting kernel to a square 2*radius+1 in size.
594 % The radius should be at least 2 times that of the sigma value, or
595 % sever clipping and aliasing may result. If not given or set to 0 the
596 % radius will be determined so as to produce the best minimal error
597 % result, which is usally much larger than is normally needed.
599 % LoG:{radius},{sigma}
600 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
601 % The supposed ideal edge detection, zero-summing kernel.
603 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
604 % approx 1.6 (according to wikipedia).
606 % DoG:{radius},{sigma1},{sigma2}
607 % "Difference of Gaussians" Kernel.
608 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
609 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
610 % The result is a zero-summing kernel.
612 % Blur:{radius},{sigma}[,{angle}]
613 % Generates a 1 dimensional or linear gaussian blur, at the angle given
614 % (current restricted to orthogonal angles). If a 'radius' is given the
615 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
616 % by a 90 degree angle.
618 % If 'sigma' is zero, you get a single pixel on a field of zeros.
620 % Note that two convolutions with two "Blur" kernels perpendicular to
621 % each other, is equivalent to a far larger "Gaussian" kernel with the
622 % same sigma value, However it is much faster to apply. This is how the
623 % "-blur" operator actually works.
625 % Comet:{width},{sigma},{angle}
626 % Blur in one direction only, much like how a bright object leaves
627 % a comet like trail. The Kernel is actually half a gaussian curve,
628 % Adding two such blurs in opposite directions produces a Blur Kernel.
629 % Angle can be rotated in multiples of 90 degrees.
631 % Note that the first argument is the width of the kernel and not the
632 % radius of the kernel.
634 % # Still to be implemented...
638 % # Set kernel values using a resize filter, and given scale (sigma)
639 % # Cylindrical or Linear. Is this possible with an image?
642 % Named Constant Convolution Kernels
644 % All these are unscaled, zero-summing kernels by default. As such for
645 % non-HDRI version of ImageMagick some form of normalization, user scaling,
646 % and biasing the results is recommended, to prevent the resulting image
649 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
650 % 45 degrees to generate the 8 angled varients of each of the kernels.
653 % Discrete Lapacian Kernels, (without normalization)
654 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
655 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
656 % Type 2 : 3x3 with center:4 edge:1 corner:-2
657 % Type 3 : 3x3 with center:4 edge:-2 corner:1
658 % Type 5 : 5x5 laplacian
659 % Type 7 : 7x7 laplacian
660 % Type 15 : 5x5 LoG (sigma approx 1.4)
661 % Type 19 : 9x9 LoG (sigma approx 1.4)
664 % Sobel 'Edge' convolution kernel (3x3)
670 % Roberts convolution kernel (3x3)
676 % Prewitt Edge convolution kernel (3x3)
682 % Prewitt's "Compass" convolution kernel (3x3)
688 % Kirsch's "Compass" convolution kernel (3x3)
694 % Frei-Chen Edge Detector is based on a kernel that is similar to
695 % the Sobel Kernel, but is designed to be isotropic. That is it takes
696 % into account the distance of the diagonal in the kernel.
699 % | sqrt(2), 0, -sqrt(2) |
702 % FreiChen:{type},{angle}
704 % Frei-Chen Pre-weighted kernels...
706 % Type 0: default un-nomalized version shown above.
708 % Type 1: Orthogonal Kernel (same as type 11 below)
710 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
713 % Type 2: Diagonal form of Kernel...
715 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
718 % However this kernel is als at the heart of the FreiChen Edge Detection
719 % Process which uses a set of 9 specially weighted kernel. These 9
720 % kernels not be normalized, but directly applied to the image. The
721 % results is then added together, to produce the intensity of an edge in
722 % a specific direction. The square root of the pixel value can then be
723 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
724 % from each other, both the direction and the strength of the edge can be
727 % Type 10: All 9 of the following pre-weighted kernels...
729 % Type 11: | 1, 0, -1 |
730 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
733 % Type 12: | 1, sqrt(2), 1 |
734 % | 0, 0, 0 | / 2*sqrt(2)
737 % Type 13: | sqrt(2), -1, 0 |
738 % | -1, 0, 1 | / 2*sqrt(2)
741 % Type 14: | 0, 1, -sqrt(2) |
742 % | -1, 0, 1 | / 2*sqrt(2)
745 % Type 15: | 0, -1, 0 |
749 % Type 16: | 1, 0, -1 |
753 % Type 17: | 1, -2, 1 |
757 % Type 18: | -2, 1, -2 |
761 % Type 19: | 1, 1, 1 |
765 % The first 4 are for edge detection, the next 4 are for line detection
766 % and the last is to add a average component to the results.
768 % Using a special type of '-1' will return all 9 pre-weighted kernels
769 % as a multi-kernel list, so that you can use them directly (without
770 % normalization) with the special "-set option:morphology:compose Plus"
771 % setting to apply the full FreiChen Edge Detection Technique.
773 % If 'type' is large it will be taken to be an actual rotation angle for
774 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
775 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
777 % WARNING: The above was layed out as per
778 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
779 % But rotated 90 degrees so direction is from left rather than the top.
780 % I have yet to find any secondary confirmation of the above. The only
781 % other source found was actual source code at
782 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
783 % Neigher paper defineds the kernels in a way that looks locical or
784 % correct when taken as a whole.
788 % Diamond:[{radius}[,{scale}]]
789 % Generate a diamond shaped kernel with given radius to the points.
790 % Kernel size will again be radius*2+1 square and defaults to radius 1,
791 % generating a 3x3 kernel that is slightly larger than a square.
793 % Square:[{radius}[,{scale}]]
794 % Generate a square shaped kernel of size radius*2+1, and defaulting
795 % to a 3x3 (radius 1).
797 % Octagon:[{radius}[,{scale}]]
798 % Generate octagonal shaped kernel of given radius and constant scale.
799 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
800 % in "Diamond" kernel.
802 % Disk:[{radius}[,{scale}]]
803 % Generate a binary disk, thresholded at the radius given, the radius
804 % may be a float-point value. Final Kernel size is floor(radius)*2+1
805 % square. A radius of 5.3 is the default.
807 % NOTE: That a low radii Disk kernels produce the same results as
808 % many of the previously defined kernels, but differ greatly at larger
809 % radii. Here is a table of equivalences...
810 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
811 % "Disk:1.5" => "Square"
812 % "Disk:2" => "Diamond:2"
813 % "Disk:2.5" => "Octagon"
814 % "Disk:2.9" => "Square:2"
815 % "Disk:3.5" => "Octagon:3"
816 % "Disk:4.5" => "Octagon:4"
817 % "Disk:5.4" => "Octagon:5"
818 % "Disk:6.4" => "Octagon:6"
819 % All other Disk shapes are unique to this kernel, but because a "Disk"
820 % is more circular when using a larger radius, using a larger radius is
821 % preferred over iterating the morphological operation.
823 % Rectangle:{geometry}
824 % Simply generate a rectangle of 1's with the size given. You can also
825 % specify the location of the 'control point', otherwise the closest
826 % pixel to the center of the rectangle is selected.
828 % Properly centered and odd sized rectangles work the best.
830 % Symbol Dilation Kernels
832 % These kernel is not a good general morphological kernel, but is used
833 % more for highlighting and marking any single pixels in an image using,
834 % a "Dilate" method as appropriate.
836 % For the same reasons iterating these kernels does not produce the
837 % same result as using a larger radius for the symbol.
839 % Plus:[{radius}[,{scale}]]
840 % Cross:[{radius}[,{scale}]]
841 % Generate a kernel in the shape of a 'plus' or a 'cross' with
842 % a each arm the length of the given radius (default 2).
844 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
846 % Ring:{radius1},{radius2}[,{scale}]
847 % A ring of the values given that falls between the two radii.
848 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
849 % This is the 'edge' pixels of the default "Disk" kernel,
850 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
852 % Hit and Miss Kernels
854 % Peak:radius1,radius2
855 % Find any peak larger than the pixels the fall between the two radii.
856 % The default ring of pixels is as per "Ring".
858 % Find flat orthogonal edges of a binary shape
860 % Find 90 degree corners of a binary shape
862 % A special kernel to thin the 'outside' of diagonals
864 % Find end points of lines (for pruning a skeletion)
865 % Two types of lines ends (default to both) can be searched for
866 % Type 0: All line ends
867 % Type 1: single kernel for 4-conneected line ends
868 % Type 2: single kernel for simple line ends
870 % Find three line junctions (within a skeletion)
871 % Type 0: all line junctions
872 % Type 1: Y Junction kernel
873 % Type 2: Diagonal T Junction kernel
874 % Type 3: Orthogonal T Junction kernel
875 % Type 4: Diagonal X Junction kernel
876 % Type 5: Orthogonal + Junction kernel
878 % Find single pixel ridges or thin lines
879 % Type 1: Fine single pixel thick lines and ridges
880 % Type 2: Find two pixel thick lines and ridges
882 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
884 % Traditional skeleton generating kernels.
885 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
886 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
887 % Type 3: Thinning skeleton based on a ressearch paper by
888 % Dan S. Bloomberg (Default Type)
890 % A huge variety of Thinning Kernels designed to preserve conectivity.
891 % many other kernel sets use these kernels as source definitions.
892 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
893 % the super and sub notations used in the source research paper.
895 % Distance Measuring Kernels
897 % Different types of distance measuring methods, which are used with the
898 % a 'Distance' morphology method for generating a gradient based on
899 % distance from an edge of a binary shape, though there is a technique
900 % for handling a anti-aliased shape.
902 % See the 'Distance' Morphological Method, for information of how it is
905 % Chebyshev:[{radius}][x{scale}[%!]]
906 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
907 % is a value of one to any neighbour, orthogonal or diagonal. One why
908 % of thinking of it is the number of squares a 'King' or 'Queen' in
909 % chess needs to traverse reach any other position on a chess board.
910 % It results in a 'square' like distance function, but one where
911 % diagonals are given a value that is closer than expected.
913 % Manhattan:[{radius}][x{scale}[%!]]
914 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
915 % Cab distance metric), it is the distance needed when you can only
916 % travel in horizontal or vertical directions only. It is the
917 % distance a 'Rook' in chess would have to travel, and results in a
918 % diamond like distances, where diagonals are further than expected.
920 % Octagonal:[{radius}][x{scale}[%!]]
921 % An interleving of Manhatten and Chebyshev metrics producing an
922 % increasing octagonally shaped distance. Distances matches those of
923 % the "Octagon" shaped kernel of the same radius. The minimum radius
924 % and default is 2, producing a 5x5 kernel.
926 % Euclidean:[{radius}][x{scale}[%!]]
927 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
928 % However by default the kernel size only has a radius of 1, which
929 % limits the distance to 'Knight' like moves, with only orthogonal and
930 % diagonal measurements being correct. As such for the default kernel
931 % you will get octagonal like distance function.
933 % However using a larger radius such as "Euclidean:4" you will get a
934 % much smoother distance gradient from the edge of the shape. Especially
935 % if the image is pre-processed to include any anti-aliasing pixels.
936 % Of course a larger kernel is slower to use, and not always needed.
938 % The first three Distance Measuring Kernels will only generate distances
939 % of exact multiples of {scale} in binary images. As such you can use a
940 % scale of 1 without loosing any information. However you also need some
941 % scaling when handling non-binary anti-aliased shapes.
943 % The "Euclidean" Distance Kernel however does generate a non-integer
944 % fractional results, and as such scaling is vital even for binary shapes.
948 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
949 const GeometryInfo *args)
962 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
964 /* Generate a new empty kernel if needed */
965 kernel=(KernelInfo *) NULL;
967 case UndefinedKernel: /* These should not call this function */
968 case UserDefinedKernel:
969 assert("Should not call this function" != (char *)NULL);
971 case LaplacianKernel: /* Named Descrete Convolution Kernels */
972 case SobelKernel: /* these are defined using other kernels */
978 case EdgesKernel: /* Hit and Miss kernels */
980 case DiagonalsKernel:
982 case LineJunctionsKernel:
984 case ConvexHullKernel:
987 break; /* A pre-generated kernel is not needed */
989 /* set to 1 to do a compile-time check that we haven't missed anything */
998 case RectangleKernel:
1005 case ChebyshevKernel:
1006 case ManhattanKernel:
1007 case OctangonalKernel:
1008 case EuclideanKernel:
1012 /* Generate the base Kernel Structure */
1013 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1014 if (kernel == (KernelInfo *) NULL)
1016 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1017 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1018 kernel->negative_range = kernel->positive_range = 0.0;
1019 kernel->type = type;
1020 kernel->next = (KernelInfo *) NULL;
1021 kernel->signature = MagickSignature;
1031 kernel->height = kernel->width = (size_t) 1;
1032 kernel->x = kernel->y = (ssize_t) 0;
1033 kernel->values=(MagickRealType *) AcquireAlignedMemory(1,
1034 sizeof(*kernel->values));
1035 if (kernel->values == (MagickRealType *) NULL)
1036 return(DestroyKernelInfo(kernel));
1037 kernel->maximum = kernel->values[0] = args->rho;
1041 case GaussianKernel:
1045 sigma = fabs(args->sigma),
1046 sigma2 = fabs(args->xi),
1049 if ( args->rho >= 1.0 )
1050 kernel->width = (size_t)args->rho*2+1;
1051 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1052 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1054 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1055 kernel->height = kernel->width;
1056 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1057 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1058 kernel->height*sizeof(*kernel->values));
1059 if (kernel->values == (MagickRealType *) NULL)
1060 return(DestroyKernelInfo(kernel));
1062 /* WARNING: The following generates a 'sampled gaussian' kernel.
1063 * What we really want is a 'discrete gaussian' kernel.
1065 * How to do this is I don't know, but appears to be basied on the
1066 * Error Function 'erf()' (intergral of a gaussian)
1069 if ( type == GaussianKernel || type == DoGKernel )
1070 { /* Calculate a Gaussian, OR positive half of a DoG */
1071 if ( sigma > MagickEpsilon )
1072 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1073 B = (double) (1.0/(Magick2PI*sigma*sigma));
1074 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1075 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1076 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1078 else /* limiting case - a unity (normalized Dirac) kernel */
1079 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1080 kernel->width*kernel->height*sizeof(double));
1081 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1085 if ( type == DoGKernel )
1086 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1087 if ( sigma2 > MagickEpsilon )
1088 { sigma = sigma2; /* simplify loop expressions */
1089 A = 1.0/(2.0*sigma*sigma);
1090 B = (double) (1.0/(Magick2PI*sigma*sigma));
1091 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1092 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1093 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1095 else /* limiting case - a unity (normalized Dirac) kernel */
1096 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1099 if ( type == LoGKernel )
1100 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1101 if ( sigma > MagickEpsilon )
1102 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1103 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1104 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1105 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1106 { R = ((double)(u*u+v*v))*A;
1107 kernel->values[i] = (1-R)*exp(-R)*B;
1110 else /* special case - generate a unity kernel */
1111 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1112 kernel->width*kernel->height*sizeof(double));
1113 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1117 /* Note the above kernels may have been 'clipped' by a user defined
1118 ** radius, producing a smaller (darker) kernel. Also for very small
1119 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1120 ** producing a very bright kernel.
1122 ** Normalization will still be needed.
1125 /* Normalize the 2D Gaussian Kernel
1127 ** NB: a CorrelateNormalize performs a normal Normalize if
1128 ** there are no negative values.
1130 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1131 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1137 sigma = fabs(args->sigma),
1140 if ( args->rho >= 1.0 )
1141 kernel->width = (size_t)args->rho*2+1;
1143 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1145 kernel->x = (ssize_t) (kernel->width-1)/2;
1147 kernel->negative_range = kernel->positive_range = 0.0;
1148 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1149 kernel->height*sizeof(*kernel->values));
1150 if (kernel->values == (MagickRealType *) NULL)
1151 return(DestroyKernelInfo(kernel));
1154 #define KernelRank 3
1155 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1156 ** It generates a gaussian 3 times the width, and compresses it into
1157 ** the expected range. This produces a closer normalization of the
1158 ** resulting kernel, especially for very low sigma values.
1159 ** As such while wierd it is prefered.
1161 ** I am told this method originally came from Photoshop.
1163 ** A properly normalized curve is generated (apart from edge clipping)
1164 ** even though we later normalize the result (for edge clipping)
1165 ** to allow the correct generation of a "Difference of Blurs".
1169 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1170 (void) ResetMagickMemory(kernel->values,0, (size_t)
1171 kernel->width*kernel->height*sizeof(double));
1172 /* Calculate a Positive 1D Gaussian */
1173 if ( sigma > MagickEpsilon )
1174 { sigma *= KernelRank; /* simplify loop expressions */
1175 alpha = 1.0/(2.0*sigma*sigma);
1176 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1177 for ( u=-v; u <= v; u++) {
1178 kernel->values[(u+v)/KernelRank] +=
1179 exp(-((double)(u*u))*alpha)*beta;
1182 else /* special case - generate a unity kernel */
1183 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1185 /* Direct calculation without curve averaging */
1187 /* Calculate a Positive Gaussian */
1188 if ( sigma > MagickEpsilon )
1189 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1190 beta = 1.0/(MagickSQ2PI*sigma);
1191 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1192 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1194 else /* special case - generate a unity kernel */
1195 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1196 kernel->width*kernel->height*sizeof(double));
1197 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1200 /* Note the above kernel may have been 'clipped' by a user defined
1201 ** radius, producing a smaller (darker) kernel. Also for very small
1202 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1203 ** producing a very bright kernel.
1205 ** Normalization will still be needed.
1208 /* Normalize the 1D Gaussian Kernel
1210 ** NB: a CorrelateNormalize performs a normal Normalize if
1211 ** there are no negative values.
1213 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1214 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1216 /* rotate the 1D kernel by given angle */
1217 RotateKernelInfo(kernel, args->xi );
1222 sigma = fabs(args->sigma),
1225 if ( args->rho < 1.0 )
1226 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1228 kernel->width = (size_t)args->rho;
1229 kernel->x = kernel->y = 0;
1231 kernel->negative_range = kernel->positive_range = 0.0;
1232 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1233 kernel->height*sizeof(*kernel->values));
1234 if (kernel->values == (MagickRealType *) NULL)
1235 return(DestroyKernelInfo(kernel));
1237 /* A comet blur is half a 1D gaussian curve, so that the object is
1238 ** blurred in one direction only. This may not be quite the right
1239 ** curve to use so may change in the future. The function must be
1240 ** normalised after generation, which also resolves any clipping.
1242 ** As we are normalizing and not subtracting gaussians,
1243 ** there is no need for a divisor in the gaussian formula
1245 ** It is less comples
1247 if ( sigma > MagickEpsilon )
1250 #define KernelRank 3
1251 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1252 (void) ResetMagickMemory(kernel->values,0, (size_t)
1253 kernel->width*sizeof(double));
1254 sigma *= KernelRank; /* simplify the loop expression */
1255 A = 1.0/(2.0*sigma*sigma);
1256 /* B = 1.0/(MagickSQ2PI*sigma); */
1257 for ( u=0; u < v; u++) {
1258 kernel->values[u/KernelRank] +=
1259 exp(-((double)(u*u))*A);
1260 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1262 for (i=0; i < (ssize_t) kernel->width; i++)
1263 kernel->positive_range += kernel->values[i];
1265 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1266 /* B = 1.0/(MagickSQ2PI*sigma); */
1267 for ( i=0; i < (ssize_t) kernel->width; i++)
1268 kernel->positive_range +=
1270 exp(-((double)(i*i))*A);
1271 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1274 else /* special case - generate a unity kernel */
1275 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1276 kernel->width*kernel->height*sizeof(double));
1277 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1278 kernel->positive_range = 1.0;
1281 kernel->minimum = 0.0;
1282 kernel->maximum = kernel->values[0];
1283 kernel->negative_range = 0.0;
1285 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1286 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1291 Convolution Kernels - Well Known Named Constant Kernels
1293 case LaplacianKernel:
1294 { switch ( (int) args->rho ) {
1296 default: /* laplacian square filter -- default */
1297 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1299 case 1: /* laplacian diamond filter */
1300 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1303 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1306 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1308 case 5: /* a 5x5 laplacian */
1309 kernel=ParseKernelArray(
1310 "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");
1312 case 7: /* a 7x7 laplacian */
1313 kernel=ParseKernelArray(
1314 "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" );
1316 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1317 kernel=ParseKernelArray(
1318 "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");
1320 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1321 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1322 kernel=ParseKernelArray(
1323 "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");
1326 if (kernel == (KernelInfo *) NULL)
1328 kernel->type = type;
1332 { /* Simple Sobel Kernel */
1333 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1334 if (kernel == (KernelInfo *) NULL)
1336 kernel->type = type;
1337 RotateKernelInfo(kernel, args->rho);
1342 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1343 if (kernel == (KernelInfo *) NULL)
1345 kernel->type = type;
1346 RotateKernelInfo(kernel, args->rho);
1351 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1352 if (kernel == (KernelInfo *) NULL)
1354 kernel->type = type;
1355 RotateKernelInfo(kernel, args->rho);
1360 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1361 if (kernel == (KernelInfo *) NULL)
1363 kernel->type = type;
1364 RotateKernelInfo(kernel, args->rho);
1369 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1370 if (kernel == (KernelInfo *) NULL)
1372 kernel->type = type;
1373 RotateKernelInfo(kernel, args->rho);
1376 case FreiChenKernel:
1377 /* Direction is set to be left to right positive */
1378 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1379 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1380 { switch ( (int) args->rho ) {
1383 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1384 if (kernel == (KernelInfo *) NULL)
1386 kernel->type = type;
1387 kernel->values[3]+=(MagickRealType) MagickSQ2;
1388 kernel->values[5]-=(MagickRealType) MagickSQ2;
1389 CalcKernelMetaData(kernel); /* recalculate meta-data */
1392 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1393 if (kernel == (KernelInfo *) NULL)
1395 kernel->type = type;
1396 kernel->values[1] = kernel->values[3]+=(MagickRealType) MagickSQ2;
1397 kernel->values[5] = kernel->values[7]-=(MagickRealType) MagickSQ2;
1398 CalcKernelMetaData(kernel); /* recalculate meta-data */
1399 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1402 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1403 if (kernel == (KernelInfo *) NULL)
1408 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1409 if (kernel == (KernelInfo *) NULL)
1411 kernel->type = type;
1412 kernel->values[3]+=(MagickRealType) MagickSQ2;
1413 kernel->values[5]-=(MagickRealType) MagickSQ2;
1414 CalcKernelMetaData(kernel); /* recalculate meta-data */
1415 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1418 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1419 if (kernel == (KernelInfo *) NULL)
1421 kernel->type = type;
1422 kernel->values[1]+=(MagickRealType) MagickSQ2;
1423 kernel->values[7]+=(MagickRealType) MagickSQ2;
1424 CalcKernelMetaData(kernel);
1425 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1428 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1429 if (kernel == (KernelInfo *) NULL)
1431 kernel->type = type;
1432 kernel->values[0]+=(MagickRealType) MagickSQ2;
1433 kernel->values[8]-=(MagickRealType) MagickSQ2;
1434 CalcKernelMetaData(kernel);
1435 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1438 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1439 if (kernel == (KernelInfo *) NULL)
1441 kernel->type = type;
1442 kernel->values[2]-=(MagickRealType) MagickSQ2;
1443 kernel->values[6]+=(MagickRealType) MagickSQ2;
1444 CalcKernelMetaData(kernel);
1445 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1448 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1449 if (kernel == (KernelInfo *) NULL)
1451 kernel->type = type;
1452 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1455 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1456 if (kernel == (KernelInfo *) NULL)
1458 kernel->type = type;
1459 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1462 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1463 if (kernel == (KernelInfo *) NULL)
1465 kernel->type = type;
1466 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1469 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1470 if (kernel == (KernelInfo *) NULL)
1472 kernel->type = type;
1473 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1476 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1477 if (kernel == (KernelInfo *) NULL)
1479 kernel->type = type;
1480 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1483 if ( fabs(args->sigma) > MagickEpsilon )
1484 /* Rotate by correctly supplied 'angle' */
1485 RotateKernelInfo(kernel, args->sigma);
1486 else if ( args->rho > 30.0 || args->rho < -30.0 )
1487 /* Rotate by out of bounds 'type' */
1488 RotateKernelInfo(kernel, args->rho);
1493 Boolean or Shaped Kernels
1497 if (args->rho < 1.0)
1498 kernel->width = kernel->height = 3; /* default radius = 1 */
1500 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1501 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1503 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1504 kernel->height*sizeof(*kernel->values));
1505 if (kernel->values == (MagickRealType *) NULL)
1506 return(DestroyKernelInfo(kernel));
1508 /* set all kernel values within diamond area to scale given */
1509 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1510 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1511 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1512 kernel->positive_range += kernel->values[i] = args->sigma;
1514 kernel->values[i] = nan;
1515 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1519 case RectangleKernel:
1522 if ( type == SquareKernel )
1524 if (args->rho < 1.0)
1525 kernel->width = kernel->height = 3; /* default radius = 1 */
1527 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1528 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1529 scale = args->sigma;
1532 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1533 if ( args->rho < 1.0 || args->sigma < 1.0 )
1534 return(DestroyKernelInfo(kernel)); /* invalid args given */
1535 kernel->width = (size_t)args->rho;
1536 kernel->height = (size_t)args->sigma;
1537 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1538 args->psi < 0.0 || args->psi > (double)kernel->height )
1539 return(DestroyKernelInfo(kernel)); /* invalid args given */
1540 kernel->x = (ssize_t) args->xi;
1541 kernel->y = (ssize_t) args->psi;
1544 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1545 kernel->height*sizeof(*kernel->values));
1546 if (kernel->values == (MagickRealType *) NULL)
1547 return(DestroyKernelInfo(kernel));
1549 /* set all kernel values to scale given */
1550 u=(ssize_t) (kernel->width*kernel->height);
1551 for ( i=0; i < u; i++)
1552 kernel->values[i] = scale;
1553 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1554 kernel->positive_range = scale*u;
1559 if (args->rho < 1.0)
1560 kernel->width = kernel->height = 5; /* default radius = 2 */
1562 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1563 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1565 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1566 kernel->height*sizeof(*kernel->values));
1567 if (kernel->values == (MagickRealType *) NULL)
1568 return(DestroyKernelInfo(kernel));
1570 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1571 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1572 if ( (labs((long) u)+labs((long) v)) <=
1573 ((long)kernel->x + (long)(kernel->x/2)) )
1574 kernel->positive_range += kernel->values[i] = args->sigma;
1576 kernel->values[i] = nan;
1577 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1583 limit = (ssize_t)(args->rho*args->rho);
1585 if (args->rho < 0.4) /* default radius approx 4.3 */
1586 kernel->width = kernel->height = 9L, limit = 18L;
1588 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1589 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1591 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1592 kernel->height*sizeof(*kernel->values));
1593 if (kernel->values == (MagickRealType *) NULL)
1594 return(DestroyKernelInfo(kernel));
1596 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1597 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1598 if ((u*u+v*v) <= limit)
1599 kernel->positive_range += kernel->values[i] = args->sigma;
1601 kernel->values[i] = nan;
1602 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1607 if (args->rho < 1.0)
1608 kernel->width = kernel->height = 5; /* default radius 2 */
1610 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1611 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1613 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1614 kernel->height*sizeof(*kernel->values));
1615 if (kernel->values == (MagickRealType *) NULL)
1616 return(DestroyKernelInfo(kernel));
1618 /* set all kernel values along axises to given scale */
1619 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1620 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1621 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1622 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1623 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1628 if (args->rho < 1.0)
1629 kernel->width = kernel->height = 5; /* default radius 2 */
1631 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1632 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1634 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1635 kernel->height*sizeof(*kernel->values));
1636 if (kernel->values == (MagickRealType *) NULL)
1637 return(DestroyKernelInfo(kernel));
1639 /* set all kernel values along axises to given scale */
1640 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1641 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1642 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1643 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1644 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1658 if (args->rho < args->sigma)
1660 kernel->width = ((size_t)args->sigma)*2+1;
1661 limit1 = (ssize_t)(args->rho*args->rho);
1662 limit2 = (ssize_t)(args->sigma*args->sigma);
1666 kernel->width = ((size_t)args->rho)*2+1;
1667 limit1 = (ssize_t)(args->sigma*args->sigma);
1668 limit2 = (ssize_t)(args->rho*args->rho);
1671 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1673 kernel->height = kernel->width;
1674 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1675 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
1676 kernel->height*sizeof(*kernel->values));
1677 if (kernel->values == (MagickRealType *) NULL)
1678 return(DestroyKernelInfo(kernel));
1680 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1681 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1682 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1683 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1684 { ssize_t radius=u*u+v*v;
1685 if (limit1 < radius && radius <= limit2)
1686 kernel->positive_range += kernel->values[i] = (double) scale;
1688 kernel->values[i] = nan;
1690 kernel->minimum = kernel->maximum = (double) scale;
1691 if ( type == PeaksKernel ) {
1692 /* set the central point in the middle */
1693 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1694 kernel->positive_range = 1.0;
1695 kernel->maximum = 1.0;
1701 kernel=AcquireKernelInfo("ThinSE:482");
1702 if (kernel == (KernelInfo *) NULL)
1704 kernel->type = type;
1705 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1710 kernel=AcquireKernelInfo("ThinSE:87");
1711 if (kernel == (KernelInfo *) NULL)
1713 kernel->type = type;
1714 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1717 case DiagonalsKernel:
1719 switch ( (int) args->rho ) {
1724 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1725 if (kernel == (KernelInfo *) NULL)
1727 kernel->type = type;
1728 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1729 if (new_kernel == (KernelInfo *) NULL)
1730 return(DestroyKernelInfo(kernel));
1731 new_kernel->type = type;
1732 LastKernelInfo(kernel)->next = new_kernel;
1733 ExpandMirrorKernelInfo(kernel);
1737 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1740 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1743 if (kernel == (KernelInfo *) NULL)
1745 kernel->type = type;
1746 RotateKernelInfo(kernel, args->sigma);
1749 case LineEndsKernel:
1750 { /* Kernels for finding the end of thin lines */
1751 switch ( (int) args->rho ) {
1754 /* set of kernels to find all end of lines */
1755 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1757 /* kernel for 4-connected line ends - no rotation */
1758 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1761 /* kernel to add for 8-connected lines - no rotation */
1762 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1765 /* kernel to add for orthogonal line ends - does not find corners */
1766 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1769 /* traditional line end - fails on last T end */
1770 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1773 if (kernel == (KernelInfo *) NULL)
1775 kernel->type = type;
1776 RotateKernelInfo(kernel, args->sigma);
1779 case LineJunctionsKernel:
1780 { /* kernels for finding the junctions of multiple lines */
1781 switch ( (int) args->rho ) {
1784 /* set of kernels to find all line junctions */
1785 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1788 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1791 /* Diagonal T Junctions */
1792 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1795 /* Orthogonal T Junctions */
1796 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1799 /* Diagonal X Junctions */
1800 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1803 /* Orthogonal X Junctions - minimal diamond kernel */
1804 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1807 if (kernel == (KernelInfo *) NULL)
1809 kernel->type = type;
1810 RotateKernelInfo(kernel, args->sigma);
1814 { /* Ridges - Ridge finding kernels */
1817 switch ( (int) args->rho ) {
1820 kernel=ParseKernelArray("3x1:0,1,0");
1821 if (kernel == (KernelInfo *) NULL)
1823 kernel->type = type;
1824 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1827 kernel=ParseKernelArray("4x1:0,1,1,0");
1828 if (kernel == (KernelInfo *) NULL)
1830 kernel->type = type;
1831 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1833 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1834 /* Unfortunatally we can not yet rotate a non-square kernel */
1835 /* But then we can't flip a non-symetrical kernel either */
1836 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1837 if (new_kernel == (KernelInfo *) NULL)
1838 return(DestroyKernelInfo(kernel));
1839 new_kernel->type = type;
1840 LastKernelInfo(kernel)->next = new_kernel;
1841 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1842 if (new_kernel == (KernelInfo *) NULL)
1843 return(DestroyKernelInfo(kernel));
1844 new_kernel->type = type;
1845 LastKernelInfo(kernel)->next = new_kernel;
1846 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1847 if (new_kernel == (KernelInfo *) NULL)
1848 return(DestroyKernelInfo(kernel));
1849 new_kernel->type = type;
1850 LastKernelInfo(kernel)->next = new_kernel;
1851 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1852 if (new_kernel == (KernelInfo *) NULL)
1853 return(DestroyKernelInfo(kernel));
1854 new_kernel->type = type;
1855 LastKernelInfo(kernel)->next = new_kernel;
1856 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1857 if (new_kernel == (KernelInfo *) NULL)
1858 return(DestroyKernelInfo(kernel));
1859 new_kernel->type = type;
1860 LastKernelInfo(kernel)->next = new_kernel;
1861 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1862 if (new_kernel == (KernelInfo *) NULL)
1863 return(DestroyKernelInfo(kernel));
1864 new_kernel->type = type;
1865 LastKernelInfo(kernel)->next = new_kernel;
1866 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1867 if (new_kernel == (KernelInfo *) NULL)
1868 return(DestroyKernelInfo(kernel));
1869 new_kernel->type = type;
1870 LastKernelInfo(kernel)->next = new_kernel;
1871 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1872 if (new_kernel == (KernelInfo *) NULL)
1873 return(DestroyKernelInfo(kernel));
1874 new_kernel->type = type;
1875 LastKernelInfo(kernel)->next = new_kernel;
1880 case ConvexHullKernel:
1884 /* first set of 8 kernels */
1885 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1886 if (kernel == (KernelInfo *) NULL)
1888 kernel->type = type;
1889 ExpandRotateKernelInfo(kernel, 90.0);
1890 /* append the mirror versions too - no flip function yet */
1891 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1892 if (new_kernel == (KernelInfo *) NULL)
1893 return(DestroyKernelInfo(kernel));
1894 new_kernel->type = type;
1895 ExpandRotateKernelInfo(new_kernel, 90.0);
1896 LastKernelInfo(kernel)->next = new_kernel;
1899 case SkeletonKernel:
1901 switch ( (int) args->rho ) {
1904 /* Traditional Skeleton...
1905 ** A cyclically rotated single kernel
1907 kernel=AcquireKernelInfo("ThinSE:482");
1908 if (kernel == (KernelInfo *) NULL)
1910 kernel->type = type;
1911 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1914 /* HIPR Variation of the cyclic skeleton
1915 ** Corners of the traditional method made more forgiving,
1916 ** but the retain the same cyclic order.
1918 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1919 if (kernel == (KernelInfo *) NULL)
1921 if (kernel->next == (KernelInfo *) NULL)
1922 return(DestroyKernelInfo(kernel));
1923 kernel->type = type;
1924 kernel->next->type = type;
1925 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1928 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1929 ** "Connectivity-Preserving Morphological Image Thransformations"
1930 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1931 ** http://www.leptonica.com/papers/conn.pdf
1933 kernel=AcquireKernelInfo(
1934 "ThinSE:41; ThinSE:42; ThinSE:43");
1935 if (kernel == (KernelInfo *) NULL)
1937 kernel->type = type;
1938 kernel->next->type = type;
1939 kernel->next->next->type = type;
1940 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1946 { /* Special kernels for general thinning, while preserving connections
1947 ** "Connectivity-Preserving Morphological Image Thransformations"
1948 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1949 ** http://www.leptonica.com/papers/conn.pdf
1951 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
1953 ** Note kernels do not specify the origin pixel, allowing them
1954 ** to be used for both thickening and thinning operations.
1956 switch ( (int) args->rho ) {
1957 /* SE for 4-connected thinning */
1958 case 41: /* SE_4_1 */
1959 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
1961 case 42: /* SE_4_2 */
1962 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
1964 case 43: /* SE_4_3 */
1965 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
1967 case 44: /* SE_4_4 */
1968 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
1970 case 45: /* SE_4_5 */
1971 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
1973 case 46: /* SE_4_6 */
1974 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
1976 case 47: /* SE_4_7 */
1977 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
1979 case 48: /* SE_4_8 */
1980 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
1982 case 49: /* SE_4_9 */
1983 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
1985 /* SE for 8-connected thinning - negatives of the above */
1986 case 81: /* SE_8_0 */
1987 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
1989 case 82: /* SE_8_2 */
1990 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
1992 case 83: /* SE_8_3 */
1993 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
1995 case 84: /* SE_8_4 */
1996 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
1998 case 85: /* SE_8_5 */
1999 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2001 case 86: /* SE_8_6 */
2002 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2004 case 87: /* SE_8_7 */
2005 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2007 case 88: /* SE_8_8 */
2008 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2010 case 89: /* SE_8_9 */
2011 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2013 /* Special combined SE kernels */
2014 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2015 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2017 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2018 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2020 case 481: /* SE_48_1 - General Connected Corner Kernel */
2021 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2024 case 482: /* SE_48_2 - General Edge Kernel */
2025 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2028 if (kernel == (KernelInfo *) NULL)
2030 kernel->type = type;
2031 RotateKernelInfo(kernel, args->sigma);
2035 Distance Measuring Kernels
2037 case ChebyshevKernel:
2039 if (args->rho < 1.0)
2040 kernel->width = kernel->height = 3; /* default radius = 1 */
2042 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2043 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2045 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
2046 kernel->height*sizeof(*kernel->values));
2047 if (kernel->values == (MagickRealType *) NULL)
2048 return(DestroyKernelInfo(kernel));
2050 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2051 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2052 kernel->positive_range += ( kernel->values[i] =
2053 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2054 kernel->maximum = kernel->values[0];
2057 case ManhattanKernel:
2059 if (args->rho < 1.0)
2060 kernel->width = kernel->height = 3; /* default radius = 1 */
2062 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2063 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2065 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
2066 kernel->height*sizeof(*kernel->values));
2067 if (kernel->values == (MagickRealType *) NULL)
2068 return(DestroyKernelInfo(kernel));
2070 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2071 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2072 kernel->positive_range += ( kernel->values[i] =
2073 args->sigma*(labs((long) u)+labs((long) v)) );
2074 kernel->maximum = kernel->values[0];
2077 case OctagonalKernel:
2079 if (args->rho < 2.0)
2080 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2082 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2083 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2085 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
2086 kernel->height*sizeof(*kernel->values));
2087 if (kernel->values == (MagickRealType *) NULL)
2088 return(DestroyKernelInfo(kernel));
2090 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2091 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2094 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2095 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2096 kernel->positive_range += kernel->values[i] =
2097 args->sigma*MagickMax(r1,r2);
2099 kernel->maximum = kernel->values[0];
2102 case EuclideanKernel:
2104 if (args->rho < 1.0)
2105 kernel->width = kernel->height = 3; /* default radius = 1 */
2107 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2108 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2110 kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
2111 kernel->height*sizeof(*kernel->values));
2112 if (kernel->values == (MagickRealType *) NULL)
2113 return(DestroyKernelInfo(kernel));
2115 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2116 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2117 kernel->positive_range += ( kernel->values[i] =
2118 args->sigma*sqrt((double)(u*u+v*v)) );
2119 kernel->maximum = kernel->values[0];
2124 /* No-Op Kernel - Basically just a single pixel on its own */
2125 kernel=ParseKernelArray("1:1");
2126 if (kernel == (KernelInfo *) NULL)
2128 kernel->type = UndefinedKernel;
2137 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2141 % C l o n e K e r n e l I n f o %
2145 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2147 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2148 % can be modified without effecting the original. The cloned kernel should
2149 % be destroyed using DestoryKernelInfo() when no longer needed.
2151 % The format of the CloneKernelInfo method is:
2153 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2155 % A description of each parameter follows:
2157 % o kernel: the Morphology/Convolution kernel to be cloned
2160 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2168 assert(kernel != (KernelInfo *) NULL);
2169 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2170 if (new_kernel == (KernelInfo *) NULL)
2172 *new_kernel=(*kernel); /* copy values in structure */
2174 /* replace the values with a copy of the values */
2175 new_kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
2176 kernel->height*sizeof(*kernel->values));
2177 if (new_kernel->values == (MagickRealType *) NULL)
2178 return(DestroyKernelInfo(new_kernel));
2179 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2180 new_kernel->values[i]=kernel->values[i];
2182 /* Also clone the next kernel in the kernel list */
2183 if ( kernel->next != (KernelInfo *) NULL ) {
2184 new_kernel->next = CloneKernelInfo(kernel->next);
2185 if ( new_kernel->next == (KernelInfo *) NULL )
2186 return(DestroyKernelInfo(new_kernel));
2193 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2197 % D e s t r o y K e r n e l I n f o %
2201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2203 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2206 % The format of the DestroyKernelInfo method is:
2208 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2210 % A description of each parameter follows:
2212 % o kernel: the Morphology/Convolution kernel to be destroyed
2215 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2217 assert(kernel != (KernelInfo *) NULL);
2218 if ( kernel->next != (KernelInfo *) NULL )
2219 kernel->next=DestroyKernelInfo(kernel->next);
2220 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2221 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2226 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2230 + E x p a n d M i r r o r K e r n e l I n f o %
2234 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2236 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2237 % sequence of 90-degree rotated kernels but providing a reflected 180
2238 % rotatation, before the -/+ 90-degree rotations.
2240 % This special rotation order produces a better, more symetrical thinning of
2243 % The format of the ExpandMirrorKernelInfo method is:
2245 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2247 % A description of each parameter follows:
2249 % o kernel: the Morphology/Convolution kernel
2251 % This function is only internel to this module, as it is not finalized,
2252 % especially with regard to non-orthogonal angles, and rotation of larger
2257 static void FlopKernelInfo(KernelInfo *kernel)
2258 { /* Do a Flop by reversing each row. */
2266 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2267 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2268 t=k[x], k[x]=k[r], k[r]=t;
2270 kernel->x = kernel->width - kernel->x - 1;
2271 angle = fmod(angle+180.0, 360.0);
2275 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2283 clone = CloneKernelInfo(last);
2284 RotateKernelInfo(clone, 180); /* flip */
2285 LastKernelInfo(last)->next = clone;
2288 clone = CloneKernelInfo(last);
2289 RotateKernelInfo(clone, 90); /* transpose */
2290 LastKernelInfo(last)->next = clone;
2293 clone = CloneKernelInfo(last);
2294 RotateKernelInfo(clone, 180); /* flop */
2295 LastKernelInfo(last)->next = clone;
2301 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2305 + E x p a n d R o t a t e K e r n e l I n f o %
2309 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2311 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2312 % incrementally by the angle given, until the kernel repeats.
2314 % WARNING: 45 degree rotations only works for 3x3 kernels.
2315 % While 90 degree roatations only works for linear and square kernels
2317 % The format of the ExpandRotateKernelInfo method is:
2319 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2321 % A description of each parameter follows:
2323 % o kernel: the Morphology/Convolution kernel
2325 % o angle: angle to rotate in degrees
2327 % This function is only internel to this module, as it is not finalized,
2328 % especially with regard to non-orthogonal angles, and rotation of larger
2332 /* Internal Routine - Return true if two kernels are the same */
2333 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2334 const KernelInfo *kernel2)
2339 /* check size and origin location */
2340 if ( kernel1->width != kernel2->width
2341 || kernel1->height != kernel2->height
2342 || kernel1->x != kernel2->x
2343 || kernel1->y != kernel2->y )
2346 /* check actual kernel values */
2347 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2348 /* Test for Nan equivalence */
2349 if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2351 if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2353 /* Test actual values are equivalent */
2354 if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
2361 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2369 clone = CloneKernelInfo(last);
2370 RotateKernelInfo(clone, angle);
2371 if ( SameKernelInfo(kernel, clone) == MagickTrue )
2373 LastKernelInfo(last)->next = clone;
2376 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2381 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2385 + C a l c M e t a K e r n a l I n f o %
2389 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2391 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2392 % using the kernel values. This should only ne used if it is not possible to
2393 % calculate that meta-data in some easier way.
2395 % It is important that the meta-data is correct before ScaleKernelInfo() is
2396 % used to perform kernel normalization.
2398 % The format of the CalcKernelMetaData method is:
2400 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2402 % A description of each parameter follows:
2404 % o kernel: the Morphology/Convolution kernel to modify
2406 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2407 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2408 % however is not true for flat-shaped morphological kernels.
2410 % WARNING: Only the specific kernel pointed to is modified, not a list of
2413 % This is an internal function and not expected to be useful outside this
2414 % module. This could change however.
2416 static void CalcKernelMetaData(KernelInfo *kernel)
2421 kernel->minimum = kernel->maximum = 0.0;
2422 kernel->negative_range = kernel->positive_range = 0.0;
2423 for (i=0; i < (kernel->width*kernel->height); i++)
2425 if ( fabs(kernel->values[i]) < MagickEpsilon )
2426 kernel->values[i] = 0.0;
2427 ( kernel->values[i] < 0)
2428 ? ( kernel->negative_range += kernel->values[i] )
2429 : ( kernel->positive_range += kernel->values[i] );
2430 Minimize(kernel->minimum, kernel->values[i]);
2431 Maximize(kernel->maximum, kernel->values[i]);
2438 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2442 % M o r p h o l o g y A p p l y %
2446 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2448 % MorphologyApply() applies a morphological method, multiple times using
2449 % a list of multiple kernels. This is the method that should be called by
2450 % other 'operators' that internally use morphology operations as part of
2453 % It is basically equivalent to as MorphologyImage() (see below) but
2454 % without any user controls. This allows internel programs to use this
2455 % function, to actually perform a specific task without possible interference
2456 % by any API user supplied settings.
2458 % It is MorphologyImage() task to extract any such user controls, and
2459 % pass them to this function for processing.
2461 % More specifically all given kernels should already be scaled, normalised,
2462 % and blended appropriatally before being parred to this routine. The
2463 % appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2465 % The format of the MorphologyApply method is:
2467 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2468 % const ssize_t iterations,const KernelInfo *kernel,
2469 % const CompositeMethod compose,const double bias,
2470 % ExceptionInfo *exception)
2472 % A description of each parameter follows:
2474 % o image: the source image
2476 % o method: the morphology method to be applied.
2478 % o iterations: apply the operation this many times (or no change).
2479 % A value of -1 means loop until no change found.
2480 % How this is applied may depend on the morphology method.
2481 % Typically this is a value of 1.
2483 % o channel: the channel type.
2485 % o kernel: An array of double representing the morphology kernel.
2487 % o compose: How to handle or merge multi-kernel results.
2488 % If 'UndefinedCompositeOp' use default for the Morphology method.
2489 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2490 % Otherwise merge the results using the compose method given.
2492 % o bias: Convolution Output Bias.
2494 % o exception: return any errors or warnings in this structure.
2498 /* Apply a Morphology Primative to an image using the given kernel.
2499 ** Two pre-created images must be provided, and no image is created.
2500 ** It returns the number of pixels that changed between the images
2501 ** for result convergence determination.
2503 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2504 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2505 ExceptionInfo *exception)
2507 #define MorphologyTag "Morphology/Image"
2526 assert(image != (Image *) NULL);
2527 assert(image->signature == MagickSignature);
2528 assert(morphology_image != (Image *) NULL);
2529 assert(morphology_image->signature == MagickSignature);
2530 assert(kernel != (KernelInfo *) NULL);
2531 assert(kernel->signature == MagickSignature);
2532 assert(exception != (ExceptionInfo *) NULL);
2533 assert(exception->signature == MagickSignature);
2539 image_view=AcquireCacheView(image);
2540 morphology_view=AcquireCacheView(morphology_image);
2541 virt_width=image->columns+kernel->width-1;
2543 /* Some methods (including convolve) needs use a reflected kernel.
2544 * Adjust 'origin' offsets to loop though kernel as a reflection.
2549 case ConvolveMorphology:
2550 case DilateMorphology:
2551 case DilateIntensityMorphology:
2552 case IterativeDistanceMorphology:
2553 /* kernel needs to used with reflection about origin */
2554 offx = (ssize_t) kernel->width-offx-1;
2555 offy = (ssize_t) kernel->height-offy-1;
2557 case ErodeMorphology:
2558 case ErodeIntensityMorphology:
2559 case HitAndMissMorphology:
2560 case ThinningMorphology:
2561 case ThickenMorphology:
2562 /* kernel is used as is, without reflection */
2565 assert("Not a Primitive Morphology Method" != (char *) NULL);
2569 if ( method == ConvolveMorphology && kernel->width == 1 )
2570 { /* Special handling (for speed) of vertical (blur) kernels.
2571 ** This performs its handling in columns rather than in rows.
2572 ** This is only done for convolve as it is the only method that
2573 ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
2575 ** Timing tests (on single CPU laptop)
2576 ** Using a vertical 1-d Blue with normal row-by-row (below)
2577 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2579 ** Using this column method
2580 ** time convert logo: -morphology Convolve Blur:0x10+90 null:
2583 ** Anthony Thyssen, 14 June 2010
2588 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2589 #pragma omp parallel for schedule(static,4) shared(progress,status)
2591 for (x=0; x < (ssize_t) image->columns; x++)
2593 register const Quantum
2605 if (status == MagickFalse)
2607 p=GetCacheViewVirtualPixels(image_view,x,-offy,1,image->rows+
2608 kernel->height-1,exception);
2609 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2610 morphology_image->rows,exception);
2611 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2616 /* offset to origin in 'p'. while 'q' points to it directly */
2619 for (y=0; y < (ssize_t) image->rows; y++)
2627 register const double
2630 register const Quantum
2633 /* Copy input image to the output image for unused channels
2634 * This removes need for 'cloning' a new image every iteration
2636 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2637 GetPixelChannels(image)),q);
2638 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2639 GetPixelChannels(image)),q);
2640 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2641 GetPixelChannels(image)),q);
2642 if (image->colorspace == CMYKColorspace)
2643 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2644 GetPixelChannels(image)),q);
2646 /* Set the bias of the weighted average output */
2651 result.black = bias;
2654 /* Weighted Average of pixels using reflected kernel
2656 ** NOTE for correct working of this operation for asymetrical
2657 ** kernels, the kernel needs to be applied in its reflected form.
2658 ** That is its values needs to be reversed.
2660 k = &kernel->values[ kernel->height-1 ];
2662 if ( (image->channel_mask != DefaultChannels) ||
2663 (image->matte == MagickFalse) )
2664 { /* No 'Sync' involved.
2665 ** Convolution is just a simple greyscale channel operation
2667 for (v=0; v < (ssize_t) kernel->height; v++) {
2668 if ( IsNan(*k) ) continue;
2669 result.red += (*k)*GetPixelRed(image,k_pixels);
2670 result.green += (*k)*GetPixelGreen(image,k_pixels);
2671 result.blue += (*k)*GetPixelBlue(image,k_pixels);
2672 if (image->colorspace == CMYKColorspace)
2673 result.black+=(*k)*GetPixelBlack(image,k_pixels);
2674 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2676 k_pixels+=GetPixelChannels(image);
2678 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2679 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
2680 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2681 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
2682 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2683 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
2684 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2685 (image->colorspace == CMYKColorspace))
2686 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
2687 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2688 (image->matte == MagickTrue))
2689 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2692 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2693 ** Weight the color channels with Alpha Channel so that
2694 ** transparent pixels are not part of the results.
2697 alpha, /* alpha weighting for colors : alpha */
2698 gamma; /* divisor, sum of color alpha weighting */
2700 count; /* alpha valus collected, number kernel values */
2704 for (v=0; v < (ssize_t) kernel->height; v++) {
2705 if ( IsNan(*k) ) continue;
2706 alpha=QuantumScale*GetPixelAlpha(image,k_pixels);
2707 gamma += alpha; /* normalize alpha weights only */
2708 count++; /* number of alpha values collected */
2709 alpha*=(*k); /* include kernel weighting now */
2710 result.red += alpha*GetPixelRed(image,k_pixels);
2711 result.green += alpha*GetPixelGreen(image,k_pixels);
2712 result.blue += alpha*GetPixelBlue(image,k_pixels);
2713 if (image->colorspace == CMYKColorspace)
2714 result.black += alpha*GetPixelBlack(image,k_pixels);
2715 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2717 k_pixels+=GetPixelChannels(image);
2719 /* Sync'ed channels, all channels are modified */
2720 gamma=(double)count/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2721 SetPixelRed(morphology_image,ClampToQuantum(gamma*result.red),q);
2722 SetPixelGreen(morphology_image,ClampToQuantum(gamma*result.green),q);
2723 SetPixelBlue(morphology_image,ClampToQuantum(gamma*result.blue),q);
2724 if (image->colorspace == CMYKColorspace)
2725 SetPixelBlack(morphology_image,ClampToQuantum(gamma*result.black),q);
2726 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2729 /* Count up changed pixels */
2730 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q))
2731 || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q))
2732 || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q))
2733 || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q))
2734 || ((image->colorspace == CMYKColorspace) &&
2735 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
2736 changed++; /* The pixel was changed in some way! */
2737 p+=GetPixelChannels(image);
2738 q+=GetPixelChannels(morphology_image);
2740 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2742 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2747 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2748 #pragma omp critical (MagickCore_MorphologyImage)
2750 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2751 if (proceed == MagickFalse)
2755 morphology_image->type=image->type;
2756 morphology_view=DestroyCacheView(morphology_view);
2757 image_view=DestroyCacheView(image_view);
2758 return(status ? (ssize_t) changed : 0);
2762 ** Normal handling of horizontal or rectangular kernels (row by row)
2764 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2765 #pragma omp parallel for schedule(static,4) shared(progress,status)
2767 for (y=0; y < (ssize_t) image->rows; y++)
2769 register const Quantum
2781 if (status == MagickFalse)
2783 p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
2784 kernel->height, exception);
2785 q=GetCacheViewAuthenticPixels(morphology_view,0,y,
2786 morphology_image->columns,1,exception);
2787 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2792 /* offset to origin in 'p'. while 'q' points to it directly */
2793 r = virt_width*offy + offx;
2795 for (x=0; x < (ssize_t) image->columns; x++)
2803 register const double
2806 register const Quantum
2814 /* Copy input image to the output image for unused channels
2815 * This removes need for 'cloning' a new image every iteration
2817 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2818 GetPixelChannels(image)),q);
2819 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2820 GetPixelChannels(image)),q);
2821 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2822 GetPixelChannels(image)),q);
2823 if (image->colorspace == CMYKColorspace)
2824 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2825 GetPixelChannels(image)),q);
2832 min.black = (MagickRealType) QuantumRange;
2837 max.black = (MagickRealType) 0;
2838 /* default result is the original pixel value */
2839 result.red = (MagickRealType) GetPixelRed(image,p+r*GetPixelChannels(image));
2840 result.green = (MagickRealType) GetPixelGreen(image,p+r*GetPixelChannels(image));
2841 result.blue = (MagickRealType) GetPixelBlue(image,p+r*GetPixelChannels(image));
2843 if (image->colorspace == CMYKColorspace)
2844 result.black = (MagickRealType) GetPixelBlack(image,p+r*GetPixelChannels(image));
2845 result.alpha=(MagickRealType) GetPixelAlpha(image,p+r*GetPixelChannels(image));
2848 case ConvolveMorphology:
2849 /* Set the bias of the weighted average output */
2854 result.black = bias;
2856 case DilateIntensityMorphology:
2857 case ErodeIntensityMorphology:
2858 /* use a boolean flag indicating when first match found */
2859 result.red = 0.0; /* result is not used otherwise */
2866 case ConvolveMorphology:
2867 /* Weighted Average of pixels using reflected kernel
2869 ** NOTE for correct working of this operation for asymetrical
2870 ** kernels, the kernel needs to be applied in its reflected form.
2871 ** That is its values needs to be reversed.
2873 ** Correlation is actually the same as this but without reflecting
2874 ** the kernel, and thus 'lower-level' that Convolution. However
2875 ** as Convolution is the more common method used, and it does not
2876 ** really cost us much in terms of processing to use a reflected
2877 ** kernel, so it is Convolution that is implemented.
2879 ** Correlation will have its kernel reflected before calling
2880 ** this function to do a Convolve.
2882 ** For more details of Correlation vs Convolution see
2883 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2885 k = &kernel->values[ kernel->width*kernel->height-1 ];
2887 if ( (image->channel_mask != DefaultChannels) ||
2888 (image->matte == MagickFalse) )
2889 { /* No 'Sync' involved.
2890 ** Convolution is simple greyscale channel operation
2892 for (v=0; v < (ssize_t) kernel->height; v++) {
2893 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2894 if ( IsNan(*k) ) continue;
2896 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2897 result.green += (*k)*
2898 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2899 result.blue += (*k)*
2900 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2901 if (image->colorspace == CMYKColorspace)
2902 result.black += (*k)*
2903 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2904 result.alpha += (*k)*
2905 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2907 k_pixels += virt_width*GetPixelChannels(image);
2909 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2910 SetPixelRed(morphology_image,ClampToQuantum(result.red),
2912 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2913 SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2915 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2916 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2918 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2919 (image->colorspace == CMYKColorspace))
2920 SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2922 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2923 (image->matte == MagickTrue))
2924 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
2928 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2929 ** Weight the color channels with Alpha Channel so that
2930 ** transparent pixels are not part of the results.
2933 alpha, /* alpha weighting for colors : alpha */
2934 gamma; /* divisor, sum of color alpha weighting */
2936 count; /* alpha valus collected, number kernel values */
2940 for (v=0; v < (ssize_t) kernel->height; v++) {
2941 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2942 if ( IsNan(*k) ) continue;
2943 alpha=QuantumScale*GetPixelAlpha(image,
2944 k_pixels+u*GetPixelChannels(image));
2945 gamma += alpha; /* normalize alpha weights only */
2946 count++; /* number of alpha values collected */
2947 alpha=alpha*(*k); /* include kernel weighting now */
2948 result.red += alpha*
2949 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2950 result.green += alpha*
2951 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2952 result.blue += alpha*
2953 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2954 if (image->colorspace == CMYKColorspace)
2955 result.black += alpha*
2956 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2957 result.alpha += (*k)*
2958 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2960 k_pixels += virt_width*GetPixelChannels(image);
2962 /* Sync'ed channels, all channels are modified */
2963 gamma=(double)count/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2964 SetPixelRed(morphology_image,
2965 ClampToQuantum(gamma*result.red),q);
2966 SetPixelGreen(morphology_image,
2967 ClampToQuantum(gamma*result.green),q);
2968 SetPixelBlue(morphology_image,
2969 ClampToQuantum(gamma*result.blue),q);
2970 if (image->colorspace == CMYKColorspace)
2971 SetPixelBlack(morphology_image,
2972 ClampToQuantum(gamma*result.black),q);
2973 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2977 case ErodeMorphology:
2978 /* Minimum Value within kernel neighbourhood
2980 ** NOTE that the kernel is not reflected for this operation!
2982 ** NOTE: in normal Greyscale Morphology, the kernel value should
2983 ** be added to the real value, this is currently not done, due to
2984 ** the nature of the boolean kernels being used.
2988 for (v=0; v < (ssize_t) kernel->height; v++) {
2989 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2990 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2991 Minimize(min.red, (double)
2992 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
2993 Minimize(min.green, (double)
2994 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
2995 Minimize(min.blue, (double)
2996 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
2997 Minimize(min.alpha, (double)
2998 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
2999 if (image->colorspace == CMYKColorspace)
3000 Minimize(min.black, (double)
3001 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3003 k_pixels += virt_width*GetPixelChannels(image);
3007 case DilateMorphology:
3008 /* Maximum Value within kernel neighbourhood
3010 ** NOTE for correct working of this operation for asymetrical
3011 ** kernels, the kernel needs to be applied in its reflected form.
3012 ** That is its values needs to be reversed.
3014 ** NOTE: in normal Greyscale Morphology, the kernel value should
3015 ** be added to the real value, this is currently not done, due to
3016 ** the nature of the boolean kernels being used.
3019 k = &kernel->values[ kernel->width*kernel->height-1 ];
3021 for (v=0; v < (ssize_t) kernel->height; v++) {
3022 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3023 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3024 Maximize(max.red, (double)
3025 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3026 Maximize(max.green, (double)
3027 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3028 Maximize(max.blue, (double)
3029 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3030 Maximize(max.alpha, (double)
3031 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3032 if (image->colorspace == CMYKColorspace)
3033 Maximize(max.black, (double)
3034 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3036 k_pixels += virt_width*GetPixelChannels(image);
3040 case HitAndMissMorphology:
3041 case ThinningMorphology:
3042 case ThickenMorphology:
3043 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3045 ** NOTE that the kernel is not reflected for this operation,
3046 ** and consists of both foreground and background pixel
3047 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3048 ** with either Nan or 0.5 values for don't care.
3050 ** Note that this will never produce a meaningless negative
3051 ** result. Such results can cause Thinning/Thicken to not work
3052 ** correctly when used against a greyscale image.
3056 for (v=0; v < (ssize_t) kernel->height; v++) {
3057 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3058 if ( IsNan(*k) ) continue;
3060 { /* minimim of foreground pixels */
3061 Minimize(min.red, (double)
3062 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3063 Minimize(min.green, (double)
3064 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3065 Minimize(min.blue, (double)
3066 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3067 Minimize(min.alpha,(double)
3068 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3069 if ( image->colorspace == CMYKColorspace)
3070 Minimize(min.black,(double)
3071 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3073 else if ( (*k) < 0.3 )
3074 { /* maximum of background pixels */
3075 Maximize(max.red, (double)
3076 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3077 Maximize(max.green, (double)
3078 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3079 Maximize(max.blue, (double)
3080 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3081 Maximize(max.alpha,(double)
3082 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3083 if (image->colorspace == CMYKColorspace)
3084 Maximize(max.black, (double)
3085 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3088 k_pixels += virt_width*GetPixelChannels(image);
3090 /* Pattern Match if difference is positive */
3091 min.red -= max.red; Maximize( min.red, 0.0 );
3092 min.green -= max.green; Maximize( min.green, 0.0 );
3093 min.blue -= max.blue; Maximize( min.blue, 0.0 );
3094 min.black -= max.black; Maximize( min.black, 0.0 );
3095 min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3098 case ErodeIntensityMorphology:
3099 /* Select Pixel with Minimum Intensity within kernel neighbourhood
3101 ** WARNING: the intensity test fails for CMYK and does not
3102 ** take into account the moderating effect of the alpha channel
3103 ** on the intensity.
3105 ** NOTE that the kernel is not reflected for this operation!
3109 for (v=0; v < (ssize_t) kernel->height; v++) {
3110 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3111 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3112 if ( result.red == 0.0 ||
3113 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) < GetPixelIntensity(morphology_image,q) ) {
3114 /* copy the whole pixel - no channel selection */
3115 SetPixelRed(morphology_image,GetPixelRed(image,
3116 k_pixels+u*GetPixelChannels(image)),q);
3117 SetPixelGreen(morphology_image,GetPixelGreen(image,
3118 k_pixels+u*GetPixelChannels(image)),q);
3119 SetPixelBlue(morphology_image,GetPixelBlue(image,
3120 k_pixels+u*GetPixelChannels(image)),q);
3121 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3122 k_pixels+u*GetPixelChannels(image)),q);
3123 if ( result.red > 0.0 ) changed++;
3127 k_pixels += virt_width*GetPixelChannels(image);
3131 case DilateIntensityMorphology:
3132 /* Select Pixel with Maximum Intensity within kernel neighbourhood
3134 ** WARNING: the intensity test fails for CMYK and does not
3135 ** take into account the moderating effect of the alpha channel
3136 ** on the intensity (yet).
3138 ** NOTE for correct working of this operation for asymetrical
3139 ** kernels, the kernel needs to be applied in its reflected form.
3140 ** That is its values needs to be reversed.
3142 k = &kernel->values[ kernel->width*kernel->height-1 ];
3144 for (v=0; v < (ssize_t) kernel->height; v++) {
3145 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3146 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3147 if ( result.red == 0.0 ||
3148 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) > GetPixelIntensity(morphology_image,q) ) {
3149 /* copy the whole pixel - no channel selection */
3150 SetPixelRed(morphology_image,GetPixelRed(image,
3151 k_pixels+u*GetPixelChannels(image)),q);
3152 SetPixelGreen(morphology_image,GetPixelGreen(image,
3153 k_pixels+u*GetPixelChannels(image)),q);
3154 SetPixelBlue(morphology_image,GetPixelBlue(image,
3155 k_pixels+u*GetPixelChannels(image)),q);
3156 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3157 k_pixels+u*GetPixelChannels(image)),q);
3158 if ( result.red > 0.0 ) changed++;
3162 k_pixels += virt_width*GetPixelChannels(image);
3166 case IterativeDistanceMorphology:
3167 /* Work out an iterative distance from black edge of a white image
3168 ** shape. Essentually white values are decreased to the smallest
3169 ** 'distance from edge' it can find.
3171 ** It works by adding kernel values to the neighbourhood, and and
3172 ** select the minimum value found. The kernel is rotated before
3173 ** use, so kernel distances match resulting distances, when a user
3174 ** provided asymmetric kernel is applied.
3177 ** This code is almost identical to True GrayScale Morphology But
3180 ** GreyDilate Kernel values added, maximum value found Kernel is
3181 ** rotated before use.
3183 ** GrayErode: Kernel values subtracted and minimum value found No
3184 ** kernel rotation used.
3186 ** Note the the Iterative Distance method is essentially a
3187 ** GrayErode, but with negative kernel values, and kernel
3188 ** rotation applied.
3190 k = &kernel->values[ kernel->width*kernel->height-1 ];
3192 for (v=0; v < (ssize_t) kernel->height; v++) {
3193 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3194 if ( IsNan(*k) ) continue;
3195 Minimize(result.red, (*k)+(double)
3196 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3197 Minimize(result.green, (*k)+(double)
3198 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3199 Minimize(result.blue, (*k)+(double)
3200 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3201 Minimize(result.alpha, (*k)+(double)
3202 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3203 if ( image->colorspace == CMYKColorspace)
3204 Maximize(result.black, (*k)+(double)
3205 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3207 k_pixels += virt_width*GetPixelChannels(image);
3211 case UndefinedMorphology:
3213 break; /* Do nothing */
3215 /* Final mathematics of results (combine with original image?)
3217 ** NOTE: Difference Morphology operators Edge* and *Hat could also
3218 ** be done here but works better with iteration as a image difference
3219 ** in the controling function (below). Thicken and Thinning however
3220 ** should be done here so thay can be iterated correctly.
3223 case HitAndMissMorphology:
3224 case ErodeMorphology:
3225 result = min; /* minimum of neighbourhood */
3227 case DilateMorphology:
3228 result = max; /* maximum of neighbourhood */
3230 case ThinningMorphology:
3231 /* subtract pattern match from original */
3232 result.red -= min.red;
3233 result.green -= min.green;
3234 result.blue -= min.blue;
3235 result.black -= min.black;
3236 result.alpha -= min.alpha;
3238 case ThickenMorphology:
3239 /* Add the pattern matchs to the original */
3240 result.red += min.red;
3241 result.green += min.green;
3242 result.blue += min.blue;
3243 result.black += min.black;
3244 result.alpha += min.alpha;
3247 /* result directly calculated or assigned */
3250 /* Assign the resulting pixel values - Clamping Result */
3252 case UndefinedMorphology:
3253 case ConvolveMorphology:
3254 case DilateIntensityMorphology:
3255 case ErodeIntensityMorphology:
3256 break; /* full pixel was directly assigned - not a channel method */
3258 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3259 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
3260 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3261 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
3262 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3263 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
3264 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3265 (image->colorspace == CMYKColorspace))
3266 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
3267 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3268 (image->matte == MagickTrue))
3269 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3272 /* Count up changed pixels */
3273 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q)) ||
3274 (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q)) ||
3275 (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q)) ||
3276 (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q)) ||
3277 ((image->colorspace == CMYKColorspace) &&
3278 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
3279 changed++; /* The pixel was changed in some way! */
3280 p+=GetPixelChannels(image);
3281 q+=GetPixelChannels(morphology_image);
3283 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3285 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3290 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3291 #pragma omp critical (MagickCore_MorphologyImage)
3293 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3294 if (proceed == MagickFalse)
3298 morphology_view=DestroyCacheView(morphology_view);
3299 image_view=DestroyCacheView(image_view);
3300 return(status ? (ssize_t)changed : -1);
3303 /* This is almost identical to the MorphologyPrimative() function above,
3304 ** but will apply the primitive directly to the actual image using two
3305 ** passes, once in each direction, with the results of the previous (and
3306 ** current) row being re-used.
3308 ** That is after each row is 'Sync'ed' into the image, the next row will
3309 ** make use of those values as part of the calculation of the next row.
3310 ** It then repeats, but going in the oppisite (bottom-up) direction.
3312 ** Because of this 're-use of results' this function can not make use
3313 ** of multi-threaded, parellel processing.
3315 static ssize_t MorphologyPrimitiveDirect(Image *image,
3316 const MorphologyMethod method,const KernelInfo *kernel,
3317 ExceptionInfo *exception)
3340 assert(image != (Image *) NULL);
3341 assert(image->signature == MagickSignature);
3342 assert(kernel != (KernelInfo *) NULL);
3343 assert(kernel->signature == MagickSignature);
3344 assert(exception != (ExceptionInfo *) NULL);
3345 assert(exception->signature == MagickSignature);
3347 /* Some methods (including convolve) needs use a reflected kernel.
3348 * Adjust 'origin' offsets to loop though kernel as a reflection.
3353 case DistanceMorphology:
3354 case VoronoiMorphology:
3355 /* kernel needs to used with reflection about origin */
3356 offx = (ssize_t) kernel->width-offx-1;
3357 offy = (ssize_t) kernel->height-offy-1;
3360 case ?????Morphology:
3361 /* kernel is used as is, without reflection */
3365 assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
3369 /* DO NOT THREAD THIS CODE! */
3370 /* two views into same image (virtual, and actual) */
3371 virt_view=AcquireCacheView(image);
3372 auth_view=AcquireCacheView(image);
3373 virt_width=image->columns+kernel->width-1;
3375 for (y=0; y < (ssize_t) image->rows; y++)
3377 register const Quantum
3389 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3390 ** we read using virtual to get virtual pixel handling, but write back
3391 ** into the same image.
3393 ** Only top half of kernel is processed as we do a single pass downward
3394 ** through the image iterating the distance function as we go.
3396 if (status == MagickFalse)
3398 p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
3400 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3402 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3404 if (status == MagickFalse)
3407 /* offset to origin in 'p'. while 'q' points to it directly */
3408 r = (ssize_t) virt_width*offy + offx;
3410 for (x=0; x < (ssize_t) image->columns; x++)
3418 register const double
3421 register const Quantum
3427 /* Starting Defaults */
3428 GetPixelInfo(image,&result);
3429 GetPixelInfoPixel(image,q,&result);
3430 if ( method != VoronoiMorphology )
3431 result.alpha = QuantumRange - result.alpha;
3434 case DistanceMorphology:
3435 /* Add kernel Value and select the minimum value found. */
3436 k = &kernel->values[ kernel->width*kernel->height-1 ];
3438 for (v=0; v <= (ssize_t) offy; v++) {
3439 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3440 if ( IsNan(*k) ) continue;
3441 Minimize(result.red, (*k)+
3442 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3443 Minimize(result.green, (*k)+
3444 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3445 Minimize(result.blue, (*k)+
3446 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3447 if (image->colorspace == CMYKColorspace)
3448 Minimize(result.black,(*k)+
3449 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3450 Minimize(result.alpha, (*k)+
3451 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3453 k_pixels += virt_width*GetPixelChannels(image);
3455 /* repeat with the just processed pixels of this row */
3456 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3457 k_pixels = q-offx*GetPixelChannels(image);
3458 for (u=0; u < (ssize_t) offx; u++, k--) {
3459 if ( x+u-offx < 0 ) continue; /* off the edge! */
3460 if ( IsNan(*k) ) continue;
3461 Minimize(result.red, (*k)+
3462 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3463 Minimize(result.green, (*k)+
3464 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3465 Minimize(result.blue, (*k)+
3466 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3467 if (image->colorspace == CMYKColorspace)
3468 Minimize(result.black,(*k)+
3469 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3470 Minimize(result.alpha,(*k)+
3471 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3474 case VoronoiMorphology:
3475 /* Apply Distance to 'Matte' channel, while coping the color
3476 ** values of the closest pixel.
3478 ** This is experimental, and realy the 'alpha' component should
3479 ** be completely separate 'masking' channel so that alpha can
3480 ** also be used as part of the results.
3482 k = &kernel->values[ kernel->width*kernel->height-1 ];
3484 for (v=0; v <= (ssize_t) offy; v++) {
3485 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3486 if ( IsNan(*k) ) continue;
3487 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3489 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3494 k_pixels += virt_width*GetPixelChannels(image);
3496 /* repeat with the just processed pixels of this row */
3497 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3498 k_pixels = q-offx*GetPixelChannels(image);
3499 for (u=0; u < (ssize_t) offx; u++, k--) {
3500 if ( x+u-offx < 0 ) continue; /* off the edge! */
3501 if ( IsNan(*k) ) continue;
3502 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3504 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3511 /* result directly calculated or assigned */
3514 /* Assign the resulting pixel values - Clamping Result */
3516 case VoronoiMorphology:
3517 SetPixelInfoPixel(image,&result,q);
3520 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3521 SetPixelRed(image,ClampToQuantum(result.red),q);
3522 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3523 SetPixelGreen(image,ClampToQuantum(result.green),q);
3524 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3525 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3526 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3527 (image->colorspace == CMYKColorspace))
3528 SetPixelBlack(image,ClampToQuantum(result.black),q);
3529 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3530 (image->matte == MagickTrue))
3531 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3534 /* Count up changed pixels */
3535 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q)) ||
3536 (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q)) ||
3537 (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q)) ||
3538 (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q)) ||
3539 ((image->colorspace == CMYKColorspace) &&
3540 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3541 changed++; /* The pixel was changed in some way! */
3543 p+=GetPixelChannels(image); /* increment pixel buffers */
3544 q+=GetPixelChannels(image);
3547 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3549 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3550 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3556 /* Do the reversed pass through the image */
3557 for (y=(ssize_t)image->rows-1; y >= 0; y--)
3559 register const Quantum
3571 if (status == MagickFalse)
3573 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3574 ** we read using virtual to get virtual pixel handling, but write back
3575 ** into the same image.
3577 ** Only the bottom half of the kernel will be processes as we
3580 p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
3581 kernel->y+1,exception);
3582 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3584 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3586 if (status == MagickFalse)
3589 /* adjust positions to end of row */
3590 p += (image->columns-1)*GetPixelChannels(image);
3591 q += (image->columns-1)*GetPixelChannels(image);
3593 /* offset to origin in 'p'. while 'q' points to it directly */
3596 for (x=(ssize_t)image->columns-1; x >= 0; x--)
3604 register const double
3607 register const Quantum
3613 /* Default - previously modified pixel */
3614 GetPixelInfo(image,&result);
3615 GetPixelInfoPixel(image,q,&result);
3616 if ( method != VoronoiMorphology )
3617 result.alpha = QuantumRange - result.alpha;
3620 case DistanceMorphology:
3621 /* Add kernel Value and select the minimum value found. */
3622 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3624 for (v=offy; v < (ssize_t) kernel->height; v++) {
3625 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3626 if ( IsNan(*k) ) continue;
3627 Minimize(result.red, (*k)+
3628 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3629 Minimize(result.green, (*k)+
3630 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3631 Minimize(result.blue, (*k)+
3632 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3633 if ( image->colorspace == CMYKColorspace)
3634 Minimize(result.black,(*k)+
3635 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3636 Minimize(result.alpha, (*k)+
3637 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3639 k_pixels += virt_width*GetPixelChannels(image);
3641 /* repeat with the just processed pixels of this row */
3642 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3643 k_pixels = q-offx*GetPixelChannels(image);
3644 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3645 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3646 if ( IsNan(*k) ) continue;
3647 Minimize(result.red, (*k)+
3648 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3649 Minimize(result.green, (*k)+
3650 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3651 Minimize(result.blue, (*k)+
3652 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3653 if ( image->colorspace == CMYKColorspace)
3654 Minimize(result.black, (*k)+
3655 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3656 Minimize(result.alpha, (*k)+
3657 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3660 case VoronoiMorphology:
3661 /* Apply Distance to 'Matte' channel, coping the closest color.
3663 ** This is experimental, and realy the 'alpha' component should
3664 ** be completely separate 'masking' channel.
3666 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3668 for (v=offy; v < (ssize_t) kernel->height; v++) {
3669 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3670 if ( IsNan(*k) ) continue;
3671 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3673 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3678 k_pixels += virt_width*GetPixelChannels(image);
3680 /* repeat with the just processed pixels of this row */
3681 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3682 k_pixels = q-offx*GetPixelChannels(image);
3683 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3684 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3685 if ( IsNan(*k) ) continue;
3686 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3688 GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
3695 /* result directly calculated or assigned */
3698 /* Assign the resulting pixel values - Clamping Result */
3700 case VoronoiMorphology:
3701 SetPixelInfoPixel(image,&result,q);
3704 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3705 SetPixelRed(image,ClampToQuantum(result.red),q);
3706 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3707 SetPixelGreen(image,ClampToQuantum(result.green),q);
3708 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3709 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3710 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3711 (image->colorspace == CMYKColorspace))
3712 SetPixelBlack(image,ClampToQuantum(result.black),q);
3713 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3714 (image->matte == MagickTrue))
3715 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3718 /* Count up changed pixels */
3719 if ( (GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q))
3720 || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q))
3721 || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q))
3722 || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q))
3723 || ((image->colorspace == CMYKColorspace) &&
3724 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3725 changed++; /* The pixel was changed in some way! */
3727 p-=GetPixelChannels(image); /* go backward through pixel buffers */
3728 q-=GetPixelChannels(image);
3730 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3732 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3733 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3739 auth_view=DestroyCacheView(auth_view);
3740 virt_view=DestroyCacheView(virt_view);
3741 return(status ? (ssize_t) changed : -1);
3744 /* Apply a Morphology by calling one of the above low level primitive
3745 ** application functions. This function handles any iteration loops,
3746 ** composition or re-iteration of results, and compound morphology methods
3747 ** that is based on multiple low-level (staged) morphology methods.
3749 ** Basically this provides the complex glue between the requested morphology
3750 ** method and raw low-level implementation (above).
3752 MagickPrivate Image *MorphologyApply(const Image *image,
3753 const MorphologyMethod method, const ssize_t iterations,
3754 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3755 ExceptionInfo *exception)
3761 *curr_image, /* Image we are working with or iterating */
3762 *work_image, /* secondary image for primitive iteration */
3763 *save_image, /* saved image - for 'edge' method only */
3764 *rslt_image; /* resultant image - after multi-kernel handling */
3767 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3768 *norm_kernel, /* the current normal un-reflected kernel */
3769 *rflt_kernel, /* the current reflected kernel (if needed) */
3770 *this_kernel; /* the kernel being applied */
3773 primitive; /* the current morphology primitive being applied */
3776 rslt_compose; /* multi-kernel compose method for results to use */
3779 special, /* do we use a direct modify function? */
3780 verbose; /* verbose output of results */
3783 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3784 method_limit, /* maximum number of compound method iterations */
3785 kernel_number, /* Loop 2: the kernel number being applied */
3786 stage_loop, /* Loop 3: primitive loop for compound morphology */
3787 stage_limit, /* how many primitives are in this compound */
3788 kernel_loop, /* Loop 4: iterate the kernel over image */
3789 kernel_limit, /* number of times to iterate kernel */
3790 count, /* total count of primitive steps applied */
3791 kernel_changed, /* total count of changed using iterated kernel */
3792 method_changed; /* total count of changed over method iteration */
3795 changed; /* number pixels changed by last primitive operation */
3800 assert(image != (Image *) NULL);
3801 assert(image->signature == MagickSignature);
3802 assert(kernel != (KernelInfo *) NULL);
3803 assert(kernel->signature == MagickSignature);
3804 assert(exception != (ExceptionInfo *) NULL);
3805 assert(exception->signature == MagickSignature);
3807 count = 0; /* number of low-level morphology primitives performed */
3808 if ( iterations == 0 )
3809 return((Image *)NULL); /* null operation - nothing to do! */
3811 kernel_limit = (size_t) iterations;
3812 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3813 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3815 verbose = IsStringTrue(GetImageArtifact(image,"verbose"));
3817 /* initialise for cleanup */
3818 curr_image = (Image *) image;
3819 curr_compose = image->compose;
3820 (void) curr_compose;
3821 work_image = save_image = rslt_image = (Image *) NULL;
3822 reflected_kernel = (KernelInfo *) NULL;
3824 /* Initialize specific methods
3825 * + which loop should use the given iteratations
3826 * + how many primitives make up the compound morphology
3827 * + multi-kernel compose method to use (by default)
3829 method_limit = 1; /* just do method once, unless otherwise set */
3830 stage_limit = 1; /* assume method is not a compound */
3831 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3832 rslt_compose = compose; /* and we are composing multi-kernels as given */
3834 case SmoothMorphology: /* 4 primitive compound morphology */
3837 case OpenMorphology: /* 2 primitive compound morphology */
3838 case OpenIntensityMorphology:
3839 case TopHatMorphology:
3840 case CloseMorphology:
3841 case CloseIntensityMorphology:
3842 case BottomHatMorphology:
3843 case EdgeMorphology:
3846 case HitAndMissMorphology:
3847 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3849 case ThinningMorphology:
3850 case ThickenMorphology:
3851 method_limit = kernel_limit; /* iterate the whole method */
3852 kernel_limit = 1; /* do not do kernel iteration */
3854 case DistanceMorphology:
3855 case VoronoiMorphology:
3856 special = MagickTrue; /* use special direct primative */
3862 /* Apply special methods with special requirments
3863 ** For example, single run only, or post-processing requirements
3865 if ( special == MagickTrue )
3867 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3868 if (rslt_image == (Image *) NULL)
3870 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3873 changed = MorphologyPrimitiveDirect(rslt_image, method,
3876 if ( IfMagickTrue(verbose) )
3877 (void) (void) FormatLocaleFile(stderr,
3878 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3879 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3880 1.0,0.0,1.0, (double) changed);
3885 if ( method == VoronoiMorphology ) {
3886 /* Preserve the alpha channel of input image - but turned off */
3887 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3889 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3890 MagickTrue,0,0,exception);
3891 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3897 /* Handle user (caller) specified multi-kernel composition method */
3898 if ( compose != UndefinedCompositeOp )
3899 rslt_compose = compose; /* override default composition for method */
3900 if ( rslt_compose == UndefinedCompositeOp )
3901 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3903 /* Some methods require a reflected kernel to use with primitives.
3904 * Create the reflected kernel for those methods. */
3906 case CorrelateMorphology:
3907 case CloseMorphology:
3908 case CloseIntensityMorphology:
3909 case BottomHatMorphology:
3910 case SmoothMorphology:
3911 reflected_kernel = CloneKernelInfo(kernel);
3912 if (reflected_kernel == (KernelInfo *) NULL)
3914 RotateKernelInfo(reflected_kernel,180);
3920 /* Loops around more primitive morpholgy methods
3921 ** erose, dilate, open, close, smooth, edge, etc...
3923 /* Loop 1: iterate the compound method */
3926 while ( method_loop < method_limit && method_changed > 0 ) {
3930 /* Loop 2: iterate over each kernel in a multi-kernel list */
3931 norm_kernel = (KernelInfo *) kernel;
3932 this_kernel = (KernelInfo *) kernel;
3933 rflt_kernel = reflected_kernel;
3936 while ( norm_kernel != NULL ) {
3938 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3939 stage_loop = 0; /* the compound morphology stage number */
3940 while ( stage_loop < stage_limit ) {
3941 stage_loop++; /* The stage of the compound morphology */
3943 /* Select primitive morphology for this stage of compound method */
3944 this_kernel = norm_kernel; /* default use unreflected kernel */
3945 primitive = method; /* Assume method is a primitive */
3947 case ErodeMorphology: /* just erode */
3948 case EdgeInMorphology: /* erode and image difference */
3949 primitive = ErodeMorphology;
3951 case DilateMorphology: /* just dilate */
3952 case EdgeOutMorphology: /* dilate and image difference */
3953 primitive = DilateMorphology;
3955 case OpenMorphology: /* erode then dialate */
3956 case TopHatMorphology: /* open and image difference */
3957 primitive = ErodeMorphology;
3958 if ( stage_loop == 2 )
3959 primitive = DilateMorphology;
3961 case OpenIntensityMorphology:
3962 primitive = ErodeIntensityMorphology;
3963 if ( stage_loop == 2 )
3964 primitive = DilateIntensityMorphology;
3966 case CloseMorphology: /* dilate, then erode */
3967 case BottomHatMorphology: /* close and image difference */
3968 this_kernel = rflt_kernel; /* use the reflected kernel */
3969 primitive = DilateMorphology;
3970 if ( stage_loop == 2 )
3971 primitive = ErodeMorphology;
3973 case CloseIntensityMorphology:
3974 this_kernel = rflt_kernel; /* use the reflected kernel */
3975 primitive = DilateIntensityMorphology;
3976 if ( stage_loop == 2 )
3977 primitive = ErodeIntensityMorphology;
3979 case SmoothMorphology: /* open, close */
3980 switch ( stage_loop ) {
3981 case 1: /* start an open method, which starts with Erode */
3982 primitive = ErodeMorphology;
3984 case 2: /* now Dilate the Erode */
3985 primitive = DilateMorphology;
3987 case 3: /* Reflect kernel a close */
3988 this_kernel = rflt_kernel; /* use the reflected kernel */
3989 primitive = DilateMorphology;
3991 case 4: /* Finish the Close */
3992 this_kernel = rflt_kernel; /* use the reflected kernel */
3993 primitive = ErodeMorphology;
3997 case EdgeMorphology: /* dilate and erode difference */
3998 primitive = DilateMorphology;
3999 if ( stage_loop == 2 ) {
4000 save_image = curr_image; /* save the image difference */
4001 curr_image = (Image *) image;
4002 primitive = ErodeMorphology;
4005 case CorrelateMorphology:
4006 /* A Correlation is a Convolution with a reflected kernel.
4007 ** However a Convolution is a weighted sum using a reflected
4008 ** kernel. It may seem stange to convert a Correlation into a
4009 ** Convolution as the Correlation is the simplier method, but
4010 ** Convolution is much more commonly used, and it makes sense to
4011 ** implement it directly so as to avoid the need to duplicate the
4012 ** kernel when it is not required (which is typically the
4015 this_kernel = rflt_kernel; /* use the reflected kernel */
4016 primitive = ConvolveMorphology;
4021 assert( this_kernel != (KernelInfo *) NULL );
4023 /* Extra information for debugging compound operations */
4024 if ( IfMagickTrue(verbose) ) {
4025 if ( stage_limit > 1 )
4026 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
4027 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
4028 method_loop,(double) stage_loop);
4029 else if ( primitive != method )
4030 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
4031 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
4037 /* Loop 4: Iterate the kernel with primitive */
4041 while ( kernel_loop < kernel_limit && changed > 0 ) {
4042 kernel_loop++; /* the iteration of this kernel */
4044 /* Create a clone as the destination image, if not yet defined */
4045 if ( work_image == (Image *) NULL )
4047 work_image=CloneImage(image,0,0,MagickTrue,exception);
4048 if (work_image == (Image *) NULL)
4050 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
4052 /* work_image->type=image->type; ??? */
4055 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4057 changed = MorphologyPrimitive(curr_image, work_image, primitive,
4058 this_kernel, bias, exception);
4060 if ( IfMagickTrue(verbose) ) {
4061 if ( kernel_loop > 1 )
4062 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4063 (void) (void) FormatLocaleFile(stderr,
4064 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4065 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4066 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4067 (double) (method_loop+kernel_loop-1),(double) kernel_number,
4068 (double) count,(double) changed);
4072 kernel_changed += changed;
4073 method_changed += changed;
4075 /* prepare next loop */
4076 { Image *tmp = work_image; /* swap images for iteration */
4077 work_image = curr_image;
4080 if ( work_image == image )
4081 work_image = (Image *) NULL; /* replace input 'image' */
4083 } /* End Loop 4: Iterate the kernel with primitive */
4085 if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed )
4086 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
4087 if ( IfMagickTrue(verbose) && stage_loop < stage_limit )
4088 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4091 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4092 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
4093 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
4094 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
4095 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
4098 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4100 /* Final Post-processing for some Compound Methods
4102 ** The removal of any 'Sync' channel flag in the Image Compositon
4103 ** below ensures the methematical compose method is applied in a
4104 ** purely mathematical way, and only to the selected channels.
4105 ** Turn off SVG composition 'alpha blending'.
4108 case EdgeOutMorphology:
4109 case EdgeInMorphology:
4110 case TopHatMorphology:
4111 case BottomHatMorphology:
4112 if ( IfMagickTrue(verbose) )
4113 (void) FormatLocaleFile(stderr,
4114 "\n%s: Difference with original image",CommandOptionToMnemonic(
4115 MagickMorphologyOptions, method) );
4116 (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
4117 MagickTrue,0,0,exception);
4119 case EdgeMorphology:
4120 if ( IfMagickTrue(verbose) )
4121 (void) FormatLocaleFile(stderr,
4122 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
4123 MagickMorphologyOptions, method) );
4124 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
4125 MagickTrue,0,0,exception);
4126 save_image = DestroyImage(save_image); /* finished with save image */
4132 /* multi-kernel handling: re-iterate, or compose results */
4133 if ( kernel->next == (KernelInfo *) NULL )
4134 rslt_image = curr_image; /* just return the resulting image */
4135 else if ( rslt_compose == NoCompositeOp )
4136 { if ( IfMagickTrue(verbose) ) {
4137 if ( this_kernel->next != (KernelInfo *) NULL )
4138 (void) FormatLocaleFile(stderr, " (re-iterate)");
4140 (void) FormatLocaleFile(stderr, " (done)");
4142 rslt_image = curr_image; /* return result, and re-iterate */
4144 else if ( rslt_image == (Image *) NULL)
4145 { if ( IfMagickTrue(verbose) )
4146 (void) FormatLocaleFile(stderr, " (save for compose)");
4147 rslt_image = curr_image;
4148 curr_image = (Image *) image; /* continue with original image */
4151 { /* Add the new 'current' result to the composition
4153 ** The removal of any 'Sync' channel flag in the Image Compositon
4154 ** below ensures the methematical compose method is applied in a
4155 ** purely mathematical way, and only to the selected channels.
4156 ** IE: Turn off SVG composition 'alpha blending'.
4158 if ( IfMagickTrue(verbose) )
4159 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4160 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4161 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4163 curr_image = DestroyImage(curr_image);
4164 curr_image = (Image *) image; /* continue with original image */
4166 if ( IfMagickTrue(verbose) )
4167 (void) FormatLocaleFile(stderr, "\n");
4169 /* loop to the next kernel in a multi-kernel list */
4170 norm_kernel = norm_kernel->next;
4171 if ( rflt_kernel != (KernelInfo *) NULL )
4172 rflt_kernel = rflt_kernel->next;
4174 } /* End Loop 2: Loop over each kernel */
4176 } /* End Loop 1: compound method interation */
4180 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4182 if ( curr_image == rslt_image )
4183 curr_image = (Image *) NULL;
4184 if ( rslt_image != (Image *) NULL )
4185 rslt_image = DestroyImage(rslt_image);
4187 if ( curr_image == rslt_image || curr_image == image )
4188 curr_image = (Image *) NULL;
4189 if ( curr_image != (Image *) NULL )
4190 curr_image = DestroyImage(curr_image);
4191 if ( work_image != (Image *) NULL )
4192 work_image = DestroyImage(work_image);
4193 if ( save_image != (Image *) NULL )
4194 save_image = DestroyImage(save_image);
4195 if ( reflected_kernel != (KernelInfo *) NULL )
4196 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4206 % M o r p h o l o g y I m a g e %
4210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4212 % MorphologyImage() applies a user supplied kernel to the image
4213 % according to the given mophology method.
4215 % This function applies any and all user defined settings before calling
4216 % the above internal function MorphologyApply().
4218 % User defined settings include...
4219 % * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4220 % * Kernel Scale/normalize settings ("-define convolve:scale=??")
4221 % This can also includes the addition of a scaled unity kernel.
4222 % * Show Kernel being applied ("-define showkernel=1")
4224 % Other operators that do not want user supplied options interfering,
4225 % especially "convolve:bias" and "showkernel" should use MorphologyApply()
4228 % The format of the MorphologyImage method is:
4230 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4231 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4233 % A description of each parameter follows:
4235 % o image: the image.
4237 % o method: the morphology method to be applied.
4239 % o iterations: apply the operation this many times (or no change).
4240 % A value of -1 means loop until no change found.
4241 % How this is applied may depend on the morphology method.
4242 % Typically this is a value of 1.
4244 % o kernel: An array of double representing the morphology kernel.
4245 % Warning: kernel may be normalized for the Convolve method.
4247 % o exception: return any errors or warnings in this structure.
4250 MagickExport Image *MorphologyImage(const Image *image,
4251 const MorphologyMethod method,const ssize_t iterations,
4252 const KernelInfo *kernel,ExceptionInfo *exception)
4266 curr_kernel = (KernelInfo *) kernel;
4268 compose = (ssize_t)UndefinedCompositeOp; /* use default for method */
4270 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4271 * This is done BEFORE the ShowKernelInfo() function is called so that
4272 * users can see the results of the 'option:convolve:scale' option.
4274 if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4278 /* Get the bias value as it will be needed */
4279 artifact = GetImageArtifact(image,"convolve:bias");
4280 if ( artifact != (const char *) NULL) {
4281 if (IfMagickFalse(IsGeometry(artifact)))
4282 (void) ThrowMagickException(exception,GetMagickModule(),
4283 OptionWarning,"InvalidSetting","'%s' '%s'",
4284 "convolve:bias",artifact);
4286 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4289 /* Scale kernel according to user wishes */
4290 artifact = GetImageArtifact(image,"convolve:scale");
4291 if ( artifact != (const char *)NULL ) {
4292 if (IfMagickFalse(IsGeometry(artifact)))
4293 (void) ThrowMagickException(exception,GetMagickModule(),
4294 OptionWarning,"InvalidSetting","'%s' '%s'",
4295 "convolve:scale",artifact);
4297 if ( curr_kernel == kernel )
4298 curr_kernel = CloneKernelInfo(kernel);
4299 if (curr_kernel == (KernelInfo *) NULL)
4300 return((Image *) NULL);
4301 ScaleGeometryKernelInfo(curr_kernel, artifact);
4306 /* display the (normalized) kernel via stderr */
4307 if ( IfMagickTrue(GetImageArtifact(image,"showkernel"))
4308 || IfMagickTrue(GetImageArtifact(image,"convolve:showkernel"))
4309 || IfMagickTrue(GetImageArtifact(image,"morphology:showkernel")) )
4310 ShowKernelInfo(curr_kernel);
4312 /* Override the default handling of multi-kernel morphology results
4313 * If 'Undefined' use the default method
4314 * If 'None' (default for 'Convolve') re-iterate previous result
4315 * Otherwise merge resulting images using compose method given.
4316 * Default for 'HitAndMiss' is 'Lighten'.
4323 artifact = GetImageArtifact(image,"morphology:compose");
4324 if ( artifact != (const char *) NULL) {
4325 parse=ParseCommandOption(MagickComposeOptions,
4326 MagickFalse,artifact);
4328 (void) ThrowMagickException(exception,GetMagickModule(),
4329 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4330 "morphology:compose",artifact);
4332 compose=(CompositeOperator)parse;
4335 /* Apply the Morphology */
4336 morphology_image = MorphologyApply(image,method,iterations,
4337 curr_kernel,compose,bias,exception);
4339 /* Cleanup and Exit */
4340 if ( curr_kernel != kernel )
4341 curr_kernel=DestroyKernelInfo(curr_kernel);
4342 return(morphology_image);
4346 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4350 + R o t a t e K e r n e l I n f o %
4354 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4356 % RotateKernelInfo() rotates the kernel by the angle given.
4358 % Currently it is restricted to 90 degree angles, of either 1D kernels
4359 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4360 % It will ignore usless rotations for specific 'named' built-in kernels.
4362 % The format of the RotateKernelInfo method is:
4364 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4366 % A description of each parameter follows:
4368 % o kernel: the Morphology/Convolution kernel
4370 % o angle: angle to rotate in degrees
4372 % This function is currently internal to this module only, but can be exported
4373 % to other modules if needed.
4375 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4377 /* angle the lower kernels first */
4378 if ( kernel->next != (KernelInfo *) NULL)
4379 RotateKernelInfo(kernel->next, angle);
4381 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4383 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4386 /* Modulus the angle */
4387 angle = fmod(angle, 360.0);
4391 if ( 337.5 < angle || angle <= 22.5 )
4392 return; /* Near zero angle - no change! - At least not at this time */
4394 /* Handle special cases */
4395 switch (kernel->type) {
4396 /* These built-in kernels are cylindrical kernels, rotating is useless */
4397 case GaussianKernel:
4402 case LaplacianKernel:
4403 case ChebyshevKernel:
4404 case ManhattanKernel:
4405 case EuclideanKernel:
4408 /* These may be rotatable at non-90 angles in the future */
4409 /* but simply rotating them in multiples of 90 degrees is useless */
4416 /* These only allows a +/-90 degree rotation (by transpose) */
4417 /* A 180 degree rotation is useless */
4419 if ( 135.0 < angle && angle <= 225.0 )
4421 if ( 225.0 < angle && angle <= 315.0 )
4428 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4429 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4431 if ( kernel->width == 3 && kernel->height == 3 )
4432 { /* Rotate a 3x3 square by 45 degree angle */
4433 MagickRealType t = kernel->values[0];
4434 kernel->values[0] = kernel->values[3];
4435 kernel->values[3] = kernel->values[6];
4436 kernel->values[6] = kernel->values[7];
4437 kernel->values[7] = kernel->values[8];
4438 kernel->values[8] = kernel->values[5];
4439 kernel->values[5] = kernel->values[2];
4440 kernel->values[2] = kernel->values[1];
4441 kernel->values[1] = t;
4442 /* rotate non-centered origin */
4443 if ( kernel->x != 1 || kernel->y != 1 ) {
4445 x = (ssize_t) kernel->x-1;
4446 y = (ssize_t) kernel->y-1;
4447 if ( x == y ) x = 0;
4448 else if ( x == 0 ) x = -y;
4449 else if ( x == -y ) y = 0;
4450 else if ( y == 0 ) y = x;
4451 kernel->x = (ssize_t) x+1;
4452 kernel->y = (ssize_t) y+1;
4454 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4455 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4458 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4460 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4462 if ( kernel->width == 1 || kernel->height == 1 )
4463 { /* Do a transpose of a 1 dimensional kernel,
4464 ** which results in a fast 90 degree rotation of some type.
4468 t = (ssize_t) kernel->width;
4469 kernel->width = kernel->height;
4470 kernel->height = (size_t) t;
4472 kernel->x = kernel->y;
4474 if ( kernel->width == 1 ) {
4475 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4476 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4478 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4479 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4482 else if ( kernel->width == kernel->height )
4483 { /* Rotate a square array of values by 90 degrees */
4489 for( i=0, x=kernel->width-1; i<=x; i++, x--)
4490 for( j=0, y=kernel->height-1; j<y; j++, y--)
4491 { t = k[i+j*kernel->width];
4492 k[i+j*kernel->width] = k[j+x*kernel->width];
4493 k[j+x*kernel->width] = k[x+y*kernel->width];
4494 k[x+y*kernel->width] = k[y+i*kernel->width];
4495 k[y+i*kernel->width] = t;
4498 /* rotate the origin - relative to center of array */
4499 { register ssize_t x,y;
4500 x = (ssize_t) (kernel->x*2-kernel->width+1);
4501 y = (ssize_t) (kernel->y*2-kernel->height+1);
4502 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4503 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4505 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4506 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4509 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4511 if ( 135.0 < angle && angle <= 225.0 )
4513 /* For a 180 degree rotation - also know as a reflection
4514 * This is actually a very very common operation!
4515 * Basically all that is needed is a reversal of the kernel data!
4516 * And a reflection of the origon
4529 j=(ssize_t) (kernel->width*kernel->height-1);
4530 for (i=0; i < j; i++, j--)
4531 t=k[i], k[i]=k[j], k[j]=t;
4533 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4534 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4535 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4536 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4538 /* At this point angle should at least between -45 (315) and +45 degrees
4539 * In the future some form of non-orthogonal angled rotates could be
4540 * performed here, posibily with a linear kernel restriction.
4547 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4551 % S c a l e G e o m e t r y K e r n e l I n f o %
4555 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4557 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4558 % provided as a "-set option:convolve:scale {geometry}" user setting,
4559 % and modifies the kernel according to the parsed arguments of that setting.
4561 % The first argument (and any normalization flags) are passed to
4562 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4563 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4564 % into the scaled/normalized kernel.
4566 % The format of the ScaleGeometryKernelInfo method is:
4568 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4569 % const double scaling_factor,const MagickStatusType normalize_flags)
4571 % A description of each parameter follows:
4573 % o kernel: the Morphology/Convolution kernel to modify
4576 % The geometry string to parse, typically from the user provided
4577 % "-set option:convolve:scale {geometry}" setting.
4580 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4581 const char *geometry)
4589 SetGeometryInfo(&args);
4590 flags = ParseGeometry(geometry, &args);
4593 /* For Debugging Geometry Input */
4594 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4595 flags, args.rho, args.sigma, args.xi, args.psi );
4598 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4599 args.rho *= 0.01, args.sigma *= 0.01;
4601 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4603 if ( (flags & SigmaValue) == 0 )
4606 /* Scale/Normalize the input kernel */
4607 ScaleKernelInfo(kernel, args.rho, flags);
4609 /* Add Unity Kernel, for blending with original */
4610 if ( (flags & SigmaValue) != 0 )
4611 UnityAddKernelInfo(kernel, args.sigma);
4616 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4620 % S c a l e K e r n e l I n f o %
4624 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4626 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4627 % without normalization of the sum of the kernel values (as per given flags).
4629 % By default (no flags given) the values within the kernel is scaled
4630 % directly using given scaling factor without change.
4632 % If either of the two 'normalize_flags' are given the kernel will first be
4633 % normalized and then further scaled by the scaling factor value given.
4635 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4636 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4637 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4638 % non-HDRI versions of IM this may cause images to have any negative results
4639 % clipped, unless some 'bias' is used.
4641 % More specifically. Kernels which only contain positive values (such as a
4642 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4643 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4645 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4646 % the kernel will be scaled by the absolute of the sum of kernel values, so
4647 % that it will generally fall within the +/- 1.0 range.
4649 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4650 % will be scaled by just the sum of the postive values, so that its output
4651 % range will again fall into the +/- 1.0 range.
4653 % For special kernels designed for locating shapes using 'Correlate', (often
4654 % only containing +1 and -1 values, representing foreground/brackground
4655 % matching) a special normalization method is provided to scale the positive
4656 % values separately to those of the negative values, so the kernel will be
4657 % forced to become a zero-sum kernel better suited to such searches.
4659 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4660 % attributes within the kernel structure have been correctly set during the
4663 % NOTE: The values used for 'normalize_flags' have been selected specifically
4664 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4665 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4667 % The format of the ScaleKernelInfo method is:
4669 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4670 % const MagickStatusType normalize_flags )
4672 % A description of each parameter follows:
4674 % o kernel: the Morphology/Convolution kernel
4677 % multiply all values (after normalization) by this factor if not
4678 % zero. If the kernel is normalized regardless of any flags.
4680 % o normalize_flags:
4681 % GeometryFlags defining normalization method to use.
4682 % specifically: NormalizeValue, CorrelateNormalizeValue,
4683 % and/or PercentValue
4686 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4687 const double scaling_factor,const GeometryFlags normalize_flags)
4696 /* do the other kernels in a multi-kernel list first */
4697 if ( kernel->next != (KernelInfo *) NULL)
4698 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4700 /* Normalization of Kernel */
4702 if ( (normalize_flags&NormalizeValue) != 0 ) {
4703 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4704 /* non-zero-summing kernel (generally positive) */
4705 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4707 /* zero-summing kernel */
4708 pos_scale = kernel->positive_range;
4710 /* Force kernel into a normalized zero-summing kernel */
4711 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4712 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4713 ? kernel->positive_range : 1.0;
4714 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
4715 ? -kernel->negative_range : 1.0;
4718 neg_scale = pos_scale;
4720 /* finialize scaling_factor for positive and negative components */
4721 pos_scale = scaling_factor/pos_scale;
4722 neg_scale = scaling_factor/neg_scale;
4724 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4725 if ( ! IsNan(kernel->values[i]) )
4726 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4728 /* convolution output range */
4729 kernel->positive_range *= pos_scale;
4730 kernel->negative_range *= neg_scale;
4731 /* maximum and minimum values in kernel */
4732 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4733 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4735 /* swap kernel settings if user's scaling factor is negative */
4736 if ( scaling_factor < MagickEpsilon ) {
4738 t = kernel->positive_range;
4739 kernel->positive_range = kernel->negative_range;
4740 kernel->negative_range = t;
4741 t = kernel->maximum;
4742 kernel->maximum = kernel->minimum;
4743 kernel->minimum = 1;
4750 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4754 % S h o w K e r n e l I n f o %
4758 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4760 % ShowKernelInfo() outputs the details of the given kernel defination to
4761 % standard error, generally due to a users 'showkernel' option request.
4763 % The format of the ShowKernel method is:
4765 % void ShowKernelInfo(const KernelInfo *kernel)
4767 % A description of each parameter follows:
4769 % o kernel: the Morphology/Convolution kernel
4772 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4780 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4782 (void) FormatLocaleFile(stderr, "Kernel");
4783 if ( kernel->next != (KernelInfo *) NULL )
4784 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4785 (void) FormatLocaleFile(stderr, " \"%s",
4786 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4787 if ( fabs(k->angle) > MagickEpsilon )
4788 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4789 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4790 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4791 (void) FormatLocaleFile(stderr,
4792 " with values from %.*lg to %.*lg\n",
4793 GetMagickPrecision(), k->minimum,
4794 GetMagickPrecision(), k->maximum);
4795 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4796 GetMagickPrecision(), k->negative_range,
4797 GetMagickPrecision(), k->positive_range);
4798 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4799 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4800 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4801 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4803 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4804 GetMagickPrecision(), k->positive_range+k->negative_range);
4805 for (i=v=0; v < k->height; v++) {
4806 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4807 for (u=0; u < k->width; u++, i++)
4808 if ( IsNan(k->values[i]) )
4809 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4811 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4812 GetMagickPrecision(), k->values[i]);
4813 (void) FormatLocaleFile(stderr,"\n");
4819 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4823 % U n i t y A d d K e r n a l I n f o %
4827 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4829 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4830 % to the given pre-scaled and normalized Kernel. This in effect adds that
4831 % amount of the original image into the resulting convolution kernel. This
4832 % value is usually provided by the user as a percentage value in the
4833 % 'convolve:scale' setting.
4835 % The resulting effect is to convert the defined kernels into blended
4836 % soft-blurs, unsharp kernels or into sharpening kernels.
4838 % The format of the UnityAdditionKernelInfo method is:
4840 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4842 % A description of each parameter follows:
4844 % o kernel: the Morphology/Convolution kernel
4847 % scaling factor for the unity kernel to be added to
4851 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4854 /* do the other kernels in a multi-kernel list first */
4855 if ( kernel->next != (KernelInfo *) NULL)
4856 UnityAddKernelInfo(kernel->next, scale);
4858 /* Add the scaled unity kernel to the existing kernel */
4859 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4860 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4866 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4870 % Z e r o K e r n e l N a n s %
4874 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4876 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4877 % the kernel with a zero value. This is typically done when the kernel will
4878 % be used in special hardware (GPU) convolution processors, to simply
4881 % The format of the ZeroKernelNans method is:
4883 % void ZeroKernelNans (KernelInfo *kernel)
4885 % A description of each parameter follows:
4887 % o kernel: the Morphology/Convolution kernel
4890 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4895 /* do the other kernels in a multi-kernel list first */
4896 if ( kernel->next != (KernelInfo *) NULL)
4897 ZeroKernelNans(kernel->next);
4899 for (i=0; i < (kernel->width*kernel->height); i++)
4900 if ( IsNan(kernel->values[i]) )
4901 kernel->values[i] = 0.0;