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-2011 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 *) AcquireMagickMemory(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 p = strchr(kernel_string, ':');
265 if ( p != (char *) NULL && p < end)
267 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
268 memcpy(token, kernel_string, (size_t) (p-kernel_string));
269 token[p-kernel_string] = '\0';
270 SetGeometryInfo(&args);
271 flags = ParseGeometry(token, &args);
273 /* Size handling and checks of geometry settings */
274 if ( (flags & WidthValue) == 0 ) /* if no width then */
275 args.rho = args.sigma; /* then width = height */
276 if ( args.rho < 1.0 ) /* if width too small */
277 args.rho = 1.0; /* then width = 1 */
278 if ( args.sigma < 1.0 ) /* if height too small */
279 args.sigma = args.rho; /* then height = width */
280 kernel->width = (size_t)args.rho;
281 kernel->height = (size_t)args.sigma;
283 /* Offset Handling and Checks */
284 if ( args.xi < 0.0 || args.psi < 0.0 )
285 return(DestroyKernelInfo(kernel));
286 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
287 : (ssize_t) (kernel->width-1)/2;
288 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
289 : (ssize_t) (kernel->height-1)/2;
290 if ( kernel->x >= (ssize_t) kernel->width ||
291 kernel->y >= (ssize_t) kernel->height )
292 return(DestroyKernelInfo(kernel));
294 p++; /* advance beyond the ':' */
297 { /* ELSE - Old old specification, forming odd-square kernel */
298 /* count up number of values given */
299 p=(const char *) kernel_string;
300 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
301 p++; /* ignore "'" chars for convolve filter usage - Cristy */
302 for (i=0; p < end; i++)
304 GetMagickToken(p,&p,token);
306 GetMagickToken(p,&p,token);
308 /* set the size of the kernel - old sized square */
309 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
310 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
311 p=(const char *) kernel_string;
312 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
313 p++; /* ignore "'" chars for convolve filter usage - Cristy */
316 /* Read in the kernel values from rest of input string argument */
317 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
318 kernel->height*sizeof(double));
319 if (kernel->values == (double *) NULL)
320 return(DestroyKernelInfo(kernel));
322 kernel->minimum = +MagickHuge;
323 kernel->maximum = -MagickHuge;
324 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; /* do not include this value in kernel */
336 kernel->values[i] = InterpretLocaleValue(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=(double *) AcquireAlignedMemory(1,sizeof(double));
1034 if (kernel->values == (double *) NULL)
1035 return(DestroyKernelInfo(kernel));
1036 kernel->maximum = kernel->values[0] = args->rho;
1040 case GaussianKernel:
1044 sigma = fabs(args->sigma),
1045 sigma2 = fabs(args->xi),
1048 if ( args->rho >= 1.0 )
1049 kernel->width = (size_t)args->rho*2+1;
1050 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1051 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1053 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1054 kernel->height = kernel->width;
1055 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1056 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1057 kernel->height*sizeof(double));
1058 if (kernel->values == (double *) NULL)
1059 return(DestroyKernelInfo(kernel));
1061 /* WARNING: The following generates a 'sampled gaussian' kernel.
1062 * What we really want is a 'discrete gaussian' kernel.
1064 * How to do this is I don't know, but appears to be basied on the
1065 * Error Function 'erf()' (intergral of a gaussian)
1068 if ( type == GaussianKernel || type == DoGKernel )
1069 { /* Calculate a Gaussian, OR positive half of a DoG */
1070 if ( sigma > MagickEpsilon )
1071 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1072 B = (double) (1.0/(Magick2PI*sigma*sigma));
1073 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1074 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1075 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1077 else /* limiting case - a unity (normalized Dirac) kernel */
1078 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1079 kernel->width*kernel->height*sizeof(double));
1080 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1084 if ( type == DoGKernel )
1085 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1086 if ( sigma2 > MagickEpsilon )
1087 { sigma = sigma2; /* simplify loop expressions */
1088 A = 1.0/(2.0*sigma*sigma);
1089 B = (double) (1.0/(Magick2PI*sigma*sigma));
1090 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1091 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1092 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1094 else /* limiting case - a unity (normalized Dirac) kernel */
1095 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1098 if ( type == LoGKernel )
1099 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1100 if ( sigma > MagickEpsilon )
1101 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1102 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1103 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1104 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1105 { R = ((double)(u*u+v*v))*A;
1106 kernel->values[i] = (1-R)*exp(-R)*B;
1109 else /* special case - generate a unity kernel */
1110 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1111 kernel->width*kernel->height*sizeof(double));
1112 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1116 /* Note the above kernels may have been 'clipped' by a user defined
1117 ** radius, producing a smaller (darker) kernel. Also for very small
1118 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1119 ** producing a very bright kernel.
1121 ** Normalization will still be needed.
1124 /* Normalize the 2D Gaussian Kernel
1126 ** NB: a CorrelateNormalize performs a normal Normalize if
1127 ** there are no negative values.
1129 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1130 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1136 sigma = fabs(args->sigma),
1139 if ( args->rho >= 1.0 )
1140 kernel->width = (size_t)args->rho*2+1;
1142 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1144 kernel->x = (ssize_t) (kernel->width-1)/2;
1146 kernel->negative_range = kernel->positive_range = 0.0;
1147 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1148 kernel->height*sizeof(double));
1149 if (kernel->values == (double *) NULL)
1150 return(DestroyKernelInfo(kernel));
1153 #define KernelRank 3
1154 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1155 ** It generates a gaussian 3 times the width, and compresses it into
1156 ** the expected range. This produces a closer normalization of the
1157 ** resulting kernel, especially for very low sigma values.
1158 ** As such while wierd it is prefered.
1160 ** I am told this method originally came from Photoshop.
1162 ** A properly normalized curve is generated (apart from edge clipping)
1163 ** even though we later normalize the result (for edge clipping)
1164 ** to allow the correct generation of a "Difference of Blurs".
1168 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1169 (void) ResetMagickMemory(kernel->values,0, (size_t)
1170 kernel->width*kernel->height*sizeof(double));
1171 /* Calculate a Positive 1D Gaussian */
1172 if ( sigma > MagickEpsilon )
1173 { sigma *= KernelRank; /* simplify loop expressions */
1174 alpha = 1.0/(2.0*sigma*sigma);
1175 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1176 for ( u=-v; u <= v; u++) {
1177 kernel->values[(u+v)/KernelRank] +=
1178 exp(-((double)(u*u))*alpha)*beta;
1181 else /* special case - generate a unity kernel */
1182 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1184 /* Direct calculation without curve averaging */
1186 /* Calculate a Positive Gaussian */
1187 if ( sigma > MagickEpsilon )
1188 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1189 beta = 1.0/(MagickSQ2PI*sigma);
1190 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1191 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1193 else /* special case - generate a unity kernel */
1194 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1195 kernel->width*kernel->height*sizeof(double));
1196 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1199 /* Note the above kernel may have been 'clipped' by a user defined
1200 ** radius, producing a smaller (darker) kernel. Also for very small
1201 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1202 ** producing a very bright kernel.
1204 ** Normalization will still be needed.
1207 /* Normalize the 1D Gaussian Kernel
1209 ** NB: a CorrelateNormalize performs a normal Normalize if
1210 ** there are no negative values.
1212 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1213 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1215 /* rotate the 1D kernel by given angle */
1216 RotateKernelInfo(kernel, args->xi );
1221 sigma = fabs(args->sigma),
1224 if ( args->rho < 1.0 )
1225 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1227 kernel->width = (size_t)args->rho;
1228 kernel->x = kernel->y = 0;
1230 kernel->negative_range = kernel->positive_range = 0.0;
1231 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1232 kernel->height*sizeof(double));
1233 if (kernel->values == (double *) NULL)
1234 return(DestroyKernelInfo(kernel));
1236 /* A comet blur is half a 1D gaussian curve, so that the object is
1237 ** blurred in one direction only. This may not be quite the right
1238 ** curve to use so may change in the future. The function must be
1239 ** normalised after generation, which also resolves any clipping.
1241 ** As we are normalizing and not subtracting gaussians,
1242 ** there is no need for a divisor in the gaussian formula
1244 ** It is less comples
1246 if ( sigma > MagickEpsilon )
1249 #define KernelRank 3
1250 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1251 (void) ResetMagickMemory(kernel->values,0, (size_t)
1252 kernel->width*sizeof(double));
1253 sigma *= KernelRank; /* simplify the loop expression */
1254 A = 1.0/(2.0*sigma*sigma);
1255 /* B = 1.0/(MagickSQ2PI*sigma); */
1256 for ( u=0; u < v; u++) {
1257 kernel->values[u/KernelRank] +=
1258 exp(-((double)(u*u))*A);
1259 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1261 for (i=0; i < (ssize_t) kernel->width; i++)
1262 kernel->positive_range += kernel->values[i];
1264 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1265 /* B = 1.0/(MagickSQ2PI*sigma); */
1266 for ( i=0; i < (ssize_t) kernel->width; i++)
1267 kernel->positive_range +=
1269 exp(-((double)(i*i))*A);
1270 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1273 else /* special case - generate a unity kernel */
1274 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1275 kernel->width*kernel->height*sizeof(double));
1276 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1277 kernel->positive_range = 1.0;
1280 kernel->minimum = 0.0;
1281 kernel->maximum = kernel->values[0];
1282 kernel->negative_range = 0.0;
1284 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1285 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1290 Convolution Kernels - Well Known Named Constant Kernels
1292 case LaplacianKernel:
1293 { switch ( (int) args->rho ) {
1295 default: /* laplacian square filter -- default */
1296 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1298 case 1: /* laplacian diamond filter */
1299 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1302 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1305 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1307 case 5: /* a 5x5 laplacian */
1308 kernel=ParseKernelArray(
1309 "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");
1311 case 7: /* a 7x7 laplacian */
1312 kernel=ParseKernelArray(
1313 "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" );
1315 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1316 kernel=ParseKernelArray(
1317 "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");
1319 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1320 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1321 kernel=ParseKernelArray(
1322 "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");
1325 if (kernel == (KernelInfo *) NULL)
1327 kernel->type = type;
1331 { /* Simple Sobel Kernel */
1332 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1333 if (kernel == (KernelInfo *) NULL)
1335 kernel->type = type;
1336 RotateKernelInfo(kernel, args->rho);
1341 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1342 if (kernel == (KernelInfo *) NULL)
1344 kernel->type = type;
1345 RotateKernelInfo(kernel, args->rho);
1350 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1351 if (kernel == (KernelInfo *) NULL)
1353 kernel->type = type;
1354 RotateKernelInfo(kernel, args->rho);
1359 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1360 if (kernel == (KernelInfo *) NULL)
1362 kernel->type = type;
1363 RotateKernelInfo(kernel, args->rho);
1368 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1369 if (kernel == (KernelInfo *) NULL)
1371 kernel->type = type;
1372 RotateKernelInfo(kernel, args->rho);
1375 case FreiChenKernel:
1376 /* Direction is set to be left to right positive */
1377 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1378 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1379 { switch ( (int) args->rho ) {
1382 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1383 if (kernel == (KernelInfo *) NULL)
1385 kernel->type = type;
1386 kernel->values[3] = +MagickSQ2;
1387 kernel->values[5] = -MagickSQ2;
1388 CalcKernelMetaData(kernel); /* recalculate meta-data */
1391 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1392 if (kernel == (KernelInfo *) NULL)
1394 kernel->type = type;
1395 kernel->values[1] = kernel->values[3] = +MagickSQ2;
1396 kernel->values[5] = kernel->values[7] = -MagickSQ2;
1397 CalcKernelMetaData(kernel); /* recalculate meta-data */
1398 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1401 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
1402 if (kernel == (KernelInfo *) NULL)
1407 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1408 if (kernel == (KernelInfo *) NULL)
1410 kernel->type = type;
1411 kernel->values[3] = +MagickSQ2;
1412 kernel->values[5] = -MagickSQ2;
1413 CalcKernelMetaData(kernel); /* recalculate meta-data */
1414 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1417 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1418 if (kernel == (KernelInfo *) NULL)
1420 kernel->type = type;
1421 kernel->values[1] = +MagickSQ2;
1422 kernel->values[7] = +MagickSQ2;
1423 CalcKernelMetaData(kernel);
1424 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1427 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1428 if (kernel == (KernelInfo *) NULL)
1430 kernel->type = type;
1431 kernel->values[0] = +MagickSQ2;
1432 kernel->values[8] = -MagickSQ2;
1433 CalcKernelMetaData(kernel);
1434 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1437 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1438 if (kernel == (KernelInfo *) NULL)
1440 kernel->type = type;
1441 kernel->values[2] = -MagickSQ2;
1442 kernel->values[6] = +MagickSQ2;
1443 CalcKernelMetaData(kernel);
1444 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1447 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1448 if (kernel == (KernelInfo *) NULL)
1450 kernel->type = type;
1451 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1454 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1455 if (kernel == (KernelInfo *) NULL)
1457 kernel->type = type;
1458 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1461 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1462 if (kernel == (KernelInfo *) NULL)
1464 kernel->type = type;
1465 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1468 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1469 if (kernel == (KernelInfo *) NULL)
1471 kernel->type = type;
1472 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1475 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1476 if (kernel == (KernelInfo *) NULL)
1478 kernel->type = type;
1479 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1482 if ( fabs(args->sigma) > MagickEpsilon )
1483 /* Rotate by correctly supplied 'angle' */
1484 RotateKernelInfo(kernel, args->sigma);
1485 else if ( args->rho > 30.0 || args->rho < -30.0 )
1486 /* Rotate by out of bounds 'type' */
1487 RotateKernelInfo(kernel, args->rho);
1492 Boolean or Shaped Kernels
1496 if (args->rho < 1.0)
1497 kernel->width = kernel->height = 3; /* default radius = 1 */
1499 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1500 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1502 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1503 kernel->height*sizeof(double));
1504 if (kernel->values == (double *) NULL)
1505 return(DestroyKernelInfo(kernel));
1507 /* set all kernel values within diamond area to scale given */
1508 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1509 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1510 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1511 kernel->positive_range += kernel->values[i] = args->sigma;
1513 kernel->values[i] = nan;
1514 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1518 case RectangleKernel:
1521 if ( type == SquareKernel )
1523 if (args->rho < 1.0)
1524 kernel->width = kernel->height = 3; /* default radius = 1 */
1526 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1527 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1528 scale = args->sigma;
1531 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1532 if ( args->rho < 1.0 || args->sigma < 1.0 )
1533 return(DestroyKernelInfo(kernel)); /* invalid args given */
1534 kernel->width = (size_t)args->rho;
1535 kernel->height = (size_t)args->sigma;
1536 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1537 args->psi < 0.0 || args->psi > (double)kernel->height )
1538 return(DestroyKernelInfo(kernel)); /* invalid args given */
1539 kernel->x = (ssize_t) args->xi;
1540 kernel->y = (ssize_t) args->psi;
1543 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1544 kernel->height*sizeof(double));
1545 if (kernel->values == (double *) NULL)
1546 return(DestroyKernelInfo(kernel));
1548 /* set all kernel values to scale given */
1549 u=(ssize_t) (kernel->width*kernel->height);
1550 for ( i=0; i < u; i++)
1551 kernel->values[i] = scale;
1552 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1553 kernel->positive_range = scale*u;
1558 if (args->rho < 1.0)
1559 kernel->width = kernel->height = 5; /* default radius = 2 */
1561 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1562 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1564 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1565 kernel->height*sizeof(double));
1566 if (kernel->values == (double *) NULL)
1567 return(DestroyKernelInfo(kernel));
1569 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1570 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1571 if ( (labs((long) u)+labs((long) v)) <=
1572 ((long)kernel->x + (long)(kernel->x/2)) )
1573 kernel->positive_range += kernel->values[i] = args->sigma;
1575 kernel->values[i] = nan;
1576 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1582 limit = (ssize_t)(args->rho*args->rho);
1584 if (args->rho < 0.4) /* default radius approx 4.3 */
1585 kernel->width = kernel->height = 9L, limit = 18L;
1587 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1588 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1590 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1591 kernel->height*sizeof(double));
1592 if (kernel->values == (double *) NULL)
1593 return(DestroyKernelInfo(kernel));
1595 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1596 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1597 if ((u*u+v*v) <= limit)
1598 kernel->positive_range += kernel->values[i] = args->sigma;
1600 kernel->values[i] = nan;
1601 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1606 if (args->rho < 1.0)
1607 kernel->width = kernel->height = 5; /* default radius 2 */
1609 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1610 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1612 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1613 kernel->height*sizeof(double));
1614 if (kernel->values == (double *) NULL)
1615 return(DestroyKernelInfo(kernel));
1617 /* set all kernel values along axises to given scale */
1618 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1619 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1620 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1621 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1622 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1627 if (args->rho < 1.0)
1628 kernel->width = kernel->height = 5; /* default radius 2 */
1630 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1631 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1633 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1634 kernel->height*sizeof(double));
1635 if (kernel->values == (double *) NULL)
1636 return(DestroyKernelInfo(kernel));
1638 /* set all kernel values along axises to given scale */
1639 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1640 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1641 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1642 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1643 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1657 if (args->rho < args->sigma)
1659 kernel->width = ((size_t)args->sigma)*2+1;
1660 limit1 = (ssize_t)(args->rho*args->rho);
1661 limit2 = (ssize_t)(args->sigma*args->sigma);
1665 kernel->width = ((size_t)args->rho)*2+1;
1666 limit1 = (ssize_t)(args->sigma*args->sigma);
1667 limit2 = (ssize_t)(args->rho*args->rho);
1670 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1672 kernel->height = kernel->width;
1673 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1674 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
1675 kernel->height*sizeof(double));
1676 if (kernel->values == (double *) NULL)
1677 return(DestroyKernelInfo(kernel));
1679 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1680 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1681 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1682 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1683 { ssize_t radius=u*u+v*v;
1684 if (limit1 < radius && radius <= limit2)
1685 kernel->positive_range += kernel->values[i] = (double) scale;
1687 kernel->values[i] = nan;
1689 kernel->minimum = kernel->maximum = (double) scale;
1690 if ( type == PeaksKernel ) {
1691 /* set the central point in the middle */
1692 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1693 kernel->positive_range = 1.0;
1694 kernel->maximum = 1.0;
1700 kernel=AcquireKernelInfo("ThinSE:482");
1701 if (kernel == (KernelInfo *) NULL)
1703 kernel->type = type;
1704 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1709 kernel=AcquireKernelInfo("ThinSE:87");
1710 if (kernel == (KernelInfo *) NULL)
1712 kernel->type = type;
1713 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1716 case DiagonalsKernel:
1718 switch ( (int) args->rho ) {
1723 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1724 if (kernel == (KernelInfo *) NULL)
1726 kernel->type = type;
1727 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1728 if (new_kernel == (KernelInfo *) NULL)
1729 return(DestroyKernelInfo(kernel));
1730 new_kernel->type = type;
1731 LastKernelInfo(kernel)->next = new_kernel;
1732 ExpandMirrorKernelInfo(kernel);
1736 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1739 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1742 if (kernel == (KernelInfo *) NULL)
1744 kernel->type = type;
1745 RotateKernelInfo(kernel, args->sigma);
1748 case LineEndsKernel:
1749 { /* Kernels for finding the end of thin lines */
1750 switch ( (int) args->rho ) {
1753 /* set of kernels to find all end of lines */
1754 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
1756 /* kernel for 4-connected line ends - no rotation */
1757 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1760 /* kernel to add for 8-connected lines - no rotation */
1761 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1764 /* kernel to add for orthogonal line ends - does not find corners */
1765 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1768 /* traditional line end - fails on last T end */
1769 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1772 if (kernel == (KernelInfo *) NULL)
1774 kernel->type = type;
1775 RotateKernelInfo(kernel, args->sigma);
1778 case LineJunctionsKernel:
1779 { /* kernels for finding the junctions of multiple lines */
1780 switch ( (int) args->rho ) {
1783 /* set of kernels to find all line junctions */
1784 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
1787 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1790 /* Diagonal T Junctions */
1791 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1794 /* Orthogonal T Junctions */
1795 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1798 /* Diagonal X Junctions */
1799 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1802 /* Orthogonal X Junctions - minimal diamond kernel */
1803 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1806 if (kernel == (KernelInfo *) NULL)
1808 kernel->type = type;
1809 RotateKernelInfo(kernel, args->sigma);
1813 { /* Ridges - Ridge finding kernels */
1816 switch ( (int) args->rho ) {
1819 kernel=ParseKernelArray("3x1:0,1,0");
1820 if (kernel == (KernelInfo *) NULL)
1822 kernel->type = type;
1823 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1826 kernel=ParseKernelArray("4x1:0,1,1,0");
1827 if (kernel == (KernelInfo *) NULL)
1829 kernel->type = type;
1830 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1832 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1833 /* Unfortunatally we can not yet rotate a non-square kernel */
1834 /* But then we can't flip a non-symetrical kernel either */
1835 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1836 if (new_kernel == (KernelInfo *) NULL)
1837 return(DestroyKernelInfo(kernel));
1838 new_kernel->type = type;
1839 LastKernelInfo(kernel)->next = new_kernel;
1840 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1841 if (new_kernel == (KernelInfo *) NULL)
1842 return(DestroyKernelInfo(kernel));
1843 new_kernel->type = type;
1844 LastKernelInfo(kernel)->next = new_kernel;
1845 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1846 if (new_kernel == (KernelInfo *) NULL)
1847 return(DestroyKernelInfo(kernel));
1848 new_kernel->type = type;
1849 LastKernelInfo(kernel)->next = new_kernel;
1850 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1851 if (new_kernel == (KernelInfo *) NULL)
1852 return(DestroyKernelInfo(kernel));
1853 new_kernel->type = type;
1854 LastKernelInfo(kernel)->next = new_kernel;
1855 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1856 if (new_kernel == (KernelInfo *) NULL)
1857 return(DestroyKernelInfo(kernel));
1858 new_kernel->type = type;
1859 LastKernelInfo(kernel)->next = new_kernel;
1860 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1861 if (new_kernel == (KernelInfo *) NULL)
1862 return(DestroyKernelInfo(kernel));
1863 new_kernel->type = type;
1864 LastKernelInfo(kernel)->next = new_kernel;
1865 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1866 if (new_kernel == (KernelInfo *) NULL)
1867 return(DestroyKernelInfo(kernel));
1868 new_kernel->type = type;
1869 LastKernelInfo(kernel)->next = new_kernel;
1870 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1871 if (new_kernel == (KernelInfo *) NULL)
1872 return(DestroyKernelInfo(kernel));
1873 new_kernel->type = type;
1874 LastKernelInfo(kernel)->next = new_kernel;
1879 case ConvexHullKernel:
1883 /* first set of 8 kernels */
1884 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1885 if (kernel == (KernelInfo *) NULL)
1887 kernel->type = type;
1888 ExpandRotateKernelInfo(kernel, 90.0);
1889 /* append the mirror versions too - no flip function yet */
1890 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1891 if (new_kernel == (KernelInfo *) NULL)
1892 return(DestroyKernelInfo(kernel));
1893 new_kernel->type = type;
1894 ExpandRotateKernelInfo(new_kernel, 90.0);
1895 LastKernelInfo(kernel)->next = new_kernel;
1898 case SkeletonKernel:
1900 switch ( (int) args->rho ) {
1903 /* Traditional Skeleton...
1904 ** A cyclically rotated single kernel
1906 kernel=AcquireKernelInfo("ThinSE:482");
1907 if (kernel == (KernelInfo *) NULL)
1909 kernel->type = type;
1910 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1913 /* HIPR Variation of the cyclic skeleton
1914 ** Corners of the traditional method made more forgiving,
1915 ** but the retain the same cyclic order.
1917 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
1918 if (kernel == (KernelInfo *) NULL)
1920 if (kernel->next == (KernelInfo *) NULL)
1921 return(DestroyKernelInfo(kernel));
1922 kernel->type = type;
1923 kernel->next->type = type;
1924 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1927 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1928 ** "Connectivity-Preserving Morphological Image Thransformations"
1929 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1930 ** http://www.leptonica.com/papers/conn.pdf
1932 kernel=AcquireKernelInfo(
1933 "ThinSE:41; ThinSE:42; ThinSE:43");
1934 if (kernel == (KernelInfo *) NULL)
1936 kernel->type = type;
1937 kernel->next->type = type;
1938 kernel->next->next->type = type;
1939 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1945 { /* Special kernels for general thinning, while preserving connections
1946 ** "Connectivity-Preserving Morphological Image Thransformations"
1947 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1948 ** http://www.leptonica.com/papers/conn.pdf
1950 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
1952 ** Note kernels do not specify the origin pixel, allowing them
1953 ** to be used for both thickening and thinning operations.
1955 switch ( (int) args->rho ) {
1956 /* SE for 4-connected thinning */
1957 case 41: /* SE_4_1 */
1958 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
1960 case 42: /* SE_4_2 */
1961 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
1963 case 43: /* SE_4_3 */
1964 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
1966 case 44: /* SE_4_4 */
1967 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
1969 case 45: /* SE_4_5 */
1970 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
1972 case 46: /* SE_4_6 */
1973 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
1975 case 47: /* SE_4_7 */
1976 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
1978 case 48: /* SE_4_8 */
1979 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
1981 case 49: /* SE_4_9 */
1982 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
1984 /* SE for 8-connected thinning - negatives of the above */
1985 case 81: /* SE_8_0 */
1986 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
1988 case 82: /* SE_8_2 */
1989 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
1991 case 83: /* SE_8_3 */
1992 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
1994 case 84: /* SE_8_4 */
1995 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
1997 case 85: /* SE_8_5 */
1998 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2000 case 86: /* SE_8_6 */
2001 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2003 case 87: /* SE_8_7 */
2004 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2006 case 88: /* SE_8_8 */
2007 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2009 case 89: /* SE_8_9 */
2010 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2012 /* Special combined SE kernels */
2013 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2014 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2016 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2017 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2019 case 481: /* SE_48_1 - General Connected Corner Kernel */
2020 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2023 case 482: /* SE_48_2 - General Edge Kernel */
2024 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2027 if (kernel == (KernelInfo *) NULL)
2029 kernel->type = type;
2030 RotateKernelInfo(kernel, args->sigma);
2034 Distance Measuring Kernels
2036 case ChebyshevKernel:
2038 if (args->rho < 1.0)
2039 kernel->width = kernel->height = 3; /* default radius = 1 */
2041 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2042 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2044 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2045 kernel->height*sizeof(double));
2046 if (kernel->values == (double *) NULL)
2047 return(DestroyKernelInfo(kernel));
2049 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2050 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2051 kernel->positive_range += ( kernel->values[i] =
2052 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2053 kernel->maximum = kernel->values[0];
2056 case ManhattanKernel:
2058 if (args->rho < 1.0)
2059 kernel->width = kernel->height = 3; /* default radius = 1 */
2061 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2062 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2064 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2065 kernel->height*sizeof(double));
2066 if (kernel->values == (double *) NULL)
2067 return(DestroyKernelInfo(kernel));
2069 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2070 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2071 kernel->positive_range += ( kernel->values[i] =
2072 args->sigma*(labs((long) u)+labs((long) v)) );
2073 kernel->maximum = kernel->values[0];
2076 case OctagonalKernel:
2078 if (args->rho < 2.0)
2079 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2081 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2082 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2084 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2085 kernel->height*sizeof(double));
2086 if (kernel->values == (double *) NULL)
2087 return(DestroyKernelInfo(kernel));
2089 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2090 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2093 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2094 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2095 kernel->positive_range += kernel->values[i] =
2096 args->sigma*MagickMax(r1,r2);
2098 kernel->maximum = kernel->values[0];
2101 case EuclideanKernel:
2103 if (args->rho < 1.0)
2104 kernel->width = kernel->height = 3; /* default radius = 1 */
2106 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2107 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2109 kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2110 kernel->height*sizeof(double));
2111 if (kernel->values == (double *) NULL)
2112 return(DestroyKernelInfo(kernel));
2114 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2115 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2116 kernel->positive_range += ( kernel->values[i] =
2117 args->sigma*sqrt((double)(u*u+v*v)) );
2118 kernel->maximum = kernel->values[0];
2123 /* No-Op Kernel - Basically just a single pixel on its own */
2124 kernel=ParseKernelArray("1:1");
2125 if (kernel == (KernelInfo *) NULL)
2127 kernel->type = UndefinedKernel;
2136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2140 % C l o n e K e r n e l I n f o %
2144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2146 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2147 % can be modified without effecting the original. The cloned kernel should
2148 % be destroyed using DestoryKernelInfo() when no longer needed.
2150 % The format of the CloneKernelInfo method is:
2152 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2154 % A description of each parameter follows:
2156 % o kernel: the Morphology/Convolution kernel to be cloned
2159 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2167 assert(kernel != (KernelInfo *) NULL);
2168 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2169 if (new_kernel == (KernelInfo *) NULL)
2171 *new_kernel=(*kernel); /* copy values in structure */
2173 /* replace the values with a copy of the values */
2174 new_kernel->values=(double *) AcquireAlignedMemory(kernel->width,
2175 kernel->height*sizeof(double));
2176 if (new_kernel->values == (double *) NULL)
2177 return(DestroyKernelInfo(new_kernel));
2178 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2179 new_kernel->values[i]=kernel->values[i];
2181 /* Also clone the next kernel in the kernel list */
2182 if ( kernel->next != (KernelInfo *) NULL ) {
2183 new_kernel->next = CloneKernelInfo(kernel->next);
2184 if ( new_kernel->next == (KernelInfo *) NULL )
2185 return(DestroyKernelInfo(new_kernel));
2192 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2196 % D e s t r o y K e r n e l I n f o %
2200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2202 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2205 % The format of the DestroyKernelInfo method is:
2207 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2209 % A description of each parameter follows:
2211 % o kernel: the Morphology/Convolution kernel to be destroyed
2214 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2216 assert(kernel != (KernelInfo *) NULL);
2218 if ( kernel->next != (KernelInfo *) NULL )
2219 kernel->next = DestroyKernelInfo(kernel->next);
2221 kernel->values = (double *)RelinquishMagickMemory(kernel->values);
2222 kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
2227 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2231 + E x p a n d M i r r o r K e r n e l I n f o %
2235 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2237 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2238 % sequence of 90-degree rotated kernels but providing a reflected 180
2239 % rotatation, before the -/+ 90-degree rotations.
2241 % This special rotation order produces a better, more symetrical thinning of
2244 % The format of the ExpandMirrorKernelInfo method is:
2246 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2248 % A description of each parameter follows:
2250 % o kernel: the Morphology/Convolution kernel
2252 % This function is only internel to this module, as it is not finalized,
2253 % especially with regard to non-orthogonal angles, and rotation of larger
2258 static void FlopKernelInfo(KernelInfo *kernel)
2259 { /* Do a Flop by reversing each row. */
2267 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2268 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2269 t=k[x], k[x]=k[r], k[r]=t;
2271 kernel->x = kernel->width - kernel->x - 1;
2272 angle = fmod(angle+180.0, 360.0);
2276 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2284 clone = CloneKernelInfo(last);
2285 RotateKernelInfo(clone, 180); /* flip */
2286 LastKernelInfo(last)->next = clone;
2289 clone = CloneKernelInfo(last);
2290 RotateKernelInfo(clone, 90); /* transpose */
2291 LastKernelInfo(last)->next = clone;
2294 clone = CloneKernelInfo(last);
2295 RotateKernelInfo(clone, 180); /* flop */
2296 LastKernelInfo(last)->next = clone;
2302 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2306 + E x p a n d R o t a t e K e r n e l I n f o %
2310 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2312 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2313 % incrementally by the angle given, until the kernel repeats.
2315 % WARNING: 45 degree rotations only works for 3x3 kernels.
2316 % While 90 degree roatations only works for linear and square kernels
2318 % The format of the ExpandRotateKernelInfo method is:
2320 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2322 % A description of each parameter follows:
2324 % o kernel: the Morphology/Convolution kernel
2326 % o angle: angle to rotate in degrees
2328 % This function is only internel to this module, as it is not finalized,
2329 % especially with regard to non-orthogonal angles, and rotation of larger
2333 /* Internal Routine - Return true if two kernels are the same */
2334 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2335 const KernelInfo *kernel2)
2340 /* check size and origin location */
2341 if ( kernel1->width != kernel2->width
2342 || kernel1->height != kernel2->height
2343 || kernel1->x != kernel2->x
2344 || kernel1->y != kernel2->y )
2347 /* check actual kernel values */
2348 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2349 /* Test for Nan equivalence */
2350 if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
2352 if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
2354 /* Test actual values are equivalent */
2355 if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
2362 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2370 clone = CloneKernelInfo(last);
2371 RotateKernelInfo(clone, angle);
2372 if ( SameKernelInfo(kernel, clone) == MagickTrue )
2374 LastKernelInfo(last)->next = clone;
2377 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2382 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2386 + C a l c M e t a K e r n a l I n f o %
2390 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2392 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2393 % using the kernel values. This should only ne used if it is not possible to
2394 % calculate that meta-data in some easier way.
2396 % It is important that the meta-data is correct before ScaleKernelInfo() is
2397 % used to perform kernel normalization.
2399 % The format of the CalcKernelMetaData method is:
2401 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2403 % A description of each parameter follows:
2405 % o kernel: the Morphology/Convolution kernel to modify
2407 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2408 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2409 % however is not true for flat-shaped morphological kernels.
2411 % WARNING: Only the specific kernel pointed to is modified, not a list of
2414 % This is an internal function and not expected to be useful outside this
2415 % module. This could change however.
2417 static void CalcKernelMetaData(KernelInfo *kernel)
2422 kernel->minimum = kernel->maximum = 0.0;
2423 kernel->negative_range = kernel->positive_range = 0.0;
2424 for (i=0; i < (kernel->width*kernel->height); i++)
2426 if ( fabs(kernel->values[i]) < MagickEpsilon )
2427 kernel->values[i] = 0.0;
2428 ( kernel->values[i] < 0)
2429 ? ( kernel->negative_range += kernel->values[i] )
2430 : ( kernel->positive_range += kernel->values[i] );
2431 Minimize(kernel->minimum, kernel->values[i]);
2432 Maximize(kernel->maximum, kernel->values[i]);
2439 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2443 % M o r p h o l o g y A p p l y %
2447 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2449 % MorphologyApply() applies a morphological method, multiple times using
2450 % a list of multiple kernels.
2452 % It is basically equivalent to as MorphologyImage() (see below) but
2453 % without any user controls. This allows internel programs to use this
2454 % function, to actually perform a specific task without possible interference
2455 % by any API user supplied settings.
2457 % It is MorphologyImage() task to extract any such user controls, and
2458 % pass them to this function for processing.
2460 % More specifically kernels are not normalized/scaled/blended by the
2461 % 'convolve:scale' Image Artifact (setting), nor is the convolve bias
2462 % (-bias setting or image->bias) loooked at, but must be supplied from the
2463 % function arguments.
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 DistanceMorphology:*/
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(dynamic,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->sync == MagickFalse) || (image->matte == MagickFalse) )
2663 { /* No 'Sync' involved.
2664 ** Convolution is simple greyscale channel operation
2666 for (v=0; v < (ssize_t) kernel->height; v++) {
2667 if ( IsNan(*k) ) continue;
2668 result.red += (*k)*GetPixelRed(image,k_pixels);
2669 result.green += (*k)*GetPixelGreen(image,k_pixels);
2670 result.blue += (*k)*GetPixelBlue(image,k_pixels);
2671 if (image->colorspace == CMYKColorspace)
2672 result.black+=(*k)*GetPixelBlack(image,k_pixels);
2673 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2675 k_pixels+=GetPixelChannels(image);
2677 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2678 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
2679 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2680 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
2681 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2682 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
2683 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2684 (image->colorspace == CMYKColorspace))
2685 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
2686 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2687 (image->matte == MagickTrue))
2688 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2691 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2692 ** Weight the color channels with Alpha Channel so that
2693 ** transparent pixels are not part of the results.
2696 alpha, /* alpha weighting of colors : kernel*alpha */
2697 gamma; /* divisor, sum of color weighting values */
2700 for (v=0; v < (ssize_t) kernel->height; v++) {
2701 if ( IsNan(*k) ) continue;
2702 alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels));
2704 result.red += alpha*GetPixelRed(image,k_pixels);
2705 result.green += alpha*GetPixelGreen(image,k_pixels);
2706 result.blue += alpha*GetPixelBlue(image,k_pixels);
2707 if (image->colorspace == CMYKColorspace)
2708 result.black += alpha*GetPixelBlack(image,k_pixels);
2709 result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
2711 k_pixels+=GetPixelChannels(image);
2713 /* Sync'ed channels, all channels are modified */
2714 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2715 SetPixelRed(morphology_image,
2716 ClampToQuantum(gamma*result.red),q);
2717 SetPixelGreen(morphology_image,
2718 ClampToQuantum(gamma*result.green),q);
2719 SetPixelBlue(morphology_image,
2720 ClampToQuantum(gamma*result.blue),q);
2721 if (image->colorspace == CMYKColorspace)
2722 SetPixelBlack(morphology_image,
2723 ClampToQuantum(gamma*result.black),q);
2724 SetPixelAlpha(morphology_image,
2725 ClampToQuantum(result.alpha),q);
2728 /* Count up changed pixels */
2729 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q))
2730 || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q))
2731 || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q))
2732 || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q))
2733 || ((image->colorspace == CMYKColorspace) &&
2734 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
2735 changed++; /* The pixel was changed in some way! */
2736 p+=GetPixelChannels(image);
2737 q+=GetPixelChannels(morphology_image);
2739 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2741 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2746 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2747 #pragma omp critical (MagickCore_MorphologyImage)
2749 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
2750 if (proceed == MagickFalse)
2754 morphology_image->type=image->type;
2755 morphology_view=DestroyCacheView(morphology_view);
2756 image_view=DestroyCacheView(image_view);
2757 return(status ? (ssize_t) changed : 0);
2761 ** Normal handling of horizontal or rectangular kernels (row by row)
2763 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2764 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2766 for (y=0; y < (ssize_t) image->rows; y++)
2768 register const Quantum
2780 if (status == MagickFalse)
2782 p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
2783 kernel->height, exception);
2784 q=GetCacheViewAuthenticPixels(morphology_view,0,y,
2785 morphology_image->columns,1,exception);
2786 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2791 /* offset to origin in 'p'. while 'q' points to it directly */
2792 r = virt_width*offy + offx;
2794 for (x=0; x < (ssize_t) image->columns; x++)
2802 register const double
2805 register const Quantum
2813 /* Copy input image to the output image for unused channels
2814 * This removes need for 'cloning' a new image every iteration
2816 SetPixelRed(morphology_image,GetPixelRed(image,p+r*
2817 GetPixelChannels(image)),q);
2818 SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
2819 GetPixelChannels(image)),q);
2820 SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
2821 GetPixelChannels(image)),q);
2822 if (image->colorspace == CMYKColorspace)
2823 SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
2824 GetPixelChannels(image)),q);
2831 min.black = (MagickRealType) QuantumRange;
2836 max.black = (MagickRealType) 0;
2837 /* default result is the original pixel value */
2838 result.red = (MagickRealType) GetPixelRed(image,p+r*GetPixelChannels(image));
2839 result.green = (MagickRealType) GetPixelGreen(image,p+r*GetPixelChannels(image));
2840 result.blue = (MagickRealType) GetPixelBlue(image,p+r*GetPixelChannels(image));
2842 if (image->colorspace == CMYKColorspace)
2843 result.black = (MagickRealType) GetPixelBlack(image,p+r*GetPixelChannels(image));
2844 result.alpha=(MagickRealType) GetPixelAlpha(image,p+r*GetPixelChannels(image));
2847 case ConvolveMorphology:
2848 /* Set the bias of the weighted average output */
2853 result.black = bias;
2855 case DilateIntensityMorphology:
2856 case ErodeIntensityMorphology:
2857 /* use a boolean flag indicating when first match found */
2858 result.red = 0.0; /* result is not used otherwise */
2865 case ConvolveMorphology:
2866 /* Weighted Average of pixels using reflected kernel
2868 ** NOTE for correct working of this operation for asymetrical
2869 ** kernels, the kernel needs to be applied in its reflected form.
2870 ** That is its values needs to be reversed.
2872 ** Correlation is actually the same as this but without reflecting
2873 ** the kernel, and thus 'lower-level' that Convolution. However
2874 ** as Convolution is the more common method used, and it does not
2875 ** really cost us much in terms of processing to use a reflected
2876 ** kernel, so it is Convolution that is implemented.
2878 ** Correlation will have its kernel reflected before calling
2879 ** this function to do a Convolve.
2881 ** For more details of Correlation vs Convolution see
2882 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2884 k = &kernel->values[ kernel->width*kernel->height-1 ];
2886 if ( (image->sync == MagickFalse) ||
2887 (image->matte == MagickFalse) )
2888 { /* No 'Sync' involved.
2889 ** Convolution is simple greyscale channel operation
2891 for (v=0; v < (ssize_t) kernel->height; v++) {
2892 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2893 if ( IsNan(*k) ) continue;
2895 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2896 result.green += (*k)*
2897 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2898 result.blue += (*k)*
2899 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2900 if (image->colorspace == CMYKColorspace)
2901 result.black += (*k)*
2902 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2903 result.alpha += (*k)*
2904 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2906 k_pixels += virt_width*GetPixelChannels(image);
2908 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2909 SetPixelRed(morphology_image,ClampToQuantum(result.red),
2911 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2912 SetPixelGreen(morphology_image,ClampToQuantum(result.green),
2914 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2915 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
2917 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2918 (image->colorspace == CMYKColorspace))
2919 SetPixelBlack(morphology_image,ClampToQuantum(result.black),
2921 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2922 (image->matte == MagickTrue))
2923 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
2927 { /* Channel 'Sync' Flag, and Alpha Channel enabled.
2928 ** Weight the color channels with Alpha Channel so that
2929 ** transparent pixels are not part of the results.
2932 alpha, /* alpha weighting of colors : kernel*alpha */
2933 gamma; /* divisor, sum of color weighting values */
2936 for (v=0; v < (ssize_t) kernel->height; v++) {
2937 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
2938 if ( IsNan(*k) ) continue;
2939 alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels+u*
2940 GetPixelChannels(image)));
2942 result.red += alpha*
2943 GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
2944 result.green += alpha*
2945 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
2946 result.blue += alpha*
2947 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
2948 if (image->colorspace == CMYKColorspace)
2949 result.black+=alpha*
2950 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
2951 result.alpha += (*k)*
2952 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
2954 k_pixels += virt_width*GetPixelChannels(image);
2956 /* Sync'ed channels, all channels are modified */
2957 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
2958 SetPixelRed(morphology_image,
2959 ClampToQuantum(gamma*result.red),q);
2960 SetPixelGreen(morphology_image,
2961 ClampToQuantum(gamma*result.green),q);
2962 SetPixelBlue(morphology_image,
2963 ClampToQuantum(gamma*result.blue),q);
2964 if (image->colorspace == CMYKColorspace)
2965 SetPixelBlack(morphology_image,
2966 ClampToQuantum(gamma*result.black),q);
2967 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
2971 case ErodeMorphology:
2972 /* Minimum Value within kernel neighbourhood
2974 ** NOTE that the kernel is not reflected for this operation!
2976 ** NOTE: in normal Greyscale Morphology, the kernel value should
2977 ** be added to the real value, this is currently not done, due to
2978 ** the nature of the boolean kernels being used.
2982 for (v=0; v < (ssize_t) kernel->height; v++) {
2983 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
2984 if ( IsNan(*k) || (*k) < 0.5 ) continue;
2985 Minimize(min.red, (double)
2986 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
2987 Minimize(min.green, (double)
2988 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
2989 Minimize(min.blue, (double)
2990 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
2991 if (image->colorspace == CMYKColorspace)
2992 Minimize(min.black,(double)
2993 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
2994 Minimize(min.alpha,(double)
2995 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
2997 k_pixels += virt_width*GetPixelChannels(image);
3001 case DilateMorphology:
3002 /* Maximum Value within kernel neighbourhood
3004 ** NOTE for correct working of this operation for asymetrical
3005 ** kernels, the kernel needs to be applied in its reflected form.
3006 ** That is its values needs to be reversed.
3008 ** NOTE: in normal Greyscale Morphology, the kernel value should
3009 ** be added to the real value, this is currently not done, due to
3010 ** the nature of the boolean kernels being used.
3013 k = &kernel->values[ kernel->width*kernel->height-1 ];
3015 for (v=0; v < (ssize_t) kernel->height; v++) {
3016 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3017 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3018 Maximize(max.red, (double)
3019 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3020 Maximize(max.green, (double)
3021 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3022 Maximize(max.blue, (double)
3023 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3024 if (image->colorspace == CMYKColorspace)
3025 Maximize(max.black, (double)
3026 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3027 Maximize(max.alpha,(double)
3028 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3030 k_pixels += virt_width*GetPixelChannels(image);
3034 case HitAndMissMorphology:
3035 case ThinningMorphology:
3036 case ThickenMorphology:
3037 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
3039 ** NOTE that the kernel is not reflected for this operation,
3040 ** and consists of both foreground and background pixel
3041 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
3042 ** with either Nan or 0.5 values for don't care.
3044 ** Note that this will never produce a meaningless negative
3045 ** result. Such results can cause Thinning/Thicken to not work
3046 ** correctly when used against a greyscale image.
3050 for (v=0; v < (ssize_t) kernel->height; v++) {
3051 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3052 if ( IsNan(*k) ) continue;
3054 { /* minimim of foreground pixels */
3055 Minimize(min.red, (double)
3056 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3057 Minimize(min.green, (double)
3058 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3059 Minimize(min.blue, (double)
3060 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3061 if ( image->colorspace == CMYKColorspace)
3062 Minimize(min.black,(double)
3063 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3064 Minimize(min.alpha,(double)
3065 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3067 else if ( (*k) < 0.3 )
3068 { /* maximum of background pixels */
3069 Maximize(max.red, (double)
3070 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3071 Maximize(max.green, (double)
3072 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3073 Maximize(max.blue, (double)
3074 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3075 if (image->colorspace == CMYKColorspace)
3076 Maximize(max.black, (double)
3077 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3078 Maximize(max.alpha,(double)
3079 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3082 k_pixels += virt_width*GetPixelChannels(image);
3084 /* Pattern Match if difference is positive */
3085 min.red -= max.red; Maximize( min.red, 0.0 );
3086 min.green -= max.green; Maximize( min.green, 0.0 );
3087 min.blue -= max.blue; Maximize( min.blue, 0.0 );
3088 min.black -= max.black; Maximize( min.black, 0.0 );
3089 min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
3092 case ErodeIntensityMorphology:
3093 /* Select Pixel with Minimum Intensity within kernel neighbourhood
3095 ** WARNING: the intensity test fails for CMYK and does not
3096 ** take into account the moderating effect of the alpha channel
3097 ** on the intensity.
3099 ** NOTE that the kernel is not reflected for this operation!
3103 for (v=0; v < (ssize_t) kernel->height; v++) {
3104 for (u=0; u < (ssize_t) kernel->width; u++, k++) {
3105 if ( IsNan(*k) || (*k) < 0.5 ) continue;
3106 if ( result.red == 0.0 ||
3107 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) < GetPixelIntensity(morphology_image,q) ) {
3108 /* copy the whole pixel - no channel selection */
3109 SetPixelRed(morphology_image,GetPixelRed(image,
3110 k_pixels+u*GetPixelChannels(image)),q);
3111 SetPixelGreen(morphology_image,GetPixelGreen(image,
3112 k_pixels+u*GetPixelChannels(image)),q);
3113 SetPixelBlue(morphology_image,GetPixelBlue(image,
3114 k_pixels+u*GetPixelChannels(image)),q);
3115 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3116 k_pixels+u*GetPixelChannels(image)),q);
3117 if ( result.red > 0.0 ) changed++;
3121 k_pixels += virt_width*GetPixelChannels(image);
3125 case DilateIntensityMorphology:
3126 /* Select Pixel with Maximum Intensity within kernel neighbourhood
3128 ** WARNING: the intensity test fails for CMYK and does not
3129 ** take into account the moderating effect of the alpha channel
3130 ** on the intensity (yet).
3132 ** NOTE for correct working of this operation for asymetrical
3133 ** kernels, the kernel needs to be applied in its reflected form.
3134 ** That is its values needs to be reversed.
3136 k = &kernel->values[ kernel->width*kernel->height-1 ];
3138 for (v=0; v < (ssize_t) kernel->height; v++) {
3139 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3140 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
3141 if ( result.red == 0.0 ||
3142 GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) > GetPixelIntensity(morphology_image,q) ) {
3143 /* copy the whole pixel - no channel selection */
3144 SetPixelRed(morphology_image,GetPixelRed(image,
3145 k_pixels+u*GetPixelChannels(image)),q);
3146 SetPixelGreen(morphology_image,GetPixelGreen(image,
3147 k_pixels+u*GetPixelChannels(image)),q);
3148 SetPixelBlue(morphology_image,GetPixelBlue(image,
3149 k_pixels+u*GetPixelChannels(image)),q);
3150 SetPixelAlpha(morphology_image,GetPixelAlpha(image,
3151 k_pixels+u*GetPixelChannels(image)),q);
3152 if ( result.red > 0.0 ) changed++;
3156 k_pixels += virt_width*GetPixelChannels(image);
3160 This code has been obsoleted by the MorphologyPrimitiveDirect() function.
3161 However it is still (almost) correct coding for Grayscale Morphology.
3164 GrayErode is equivalent but with kernel values subtracted from pixels
3165 without the kernel rotation
3166 GreyDilate is equivalent but using Maximum() instead of Minimum()
3167 using kernel rotation
3169 It has thus been preserved for future implementation of those methods.
3171 case DistanceMorphology:
3172 /* Add kernel Value and select the minimum value found.
3173 ** The result is a iterative distance from edge of image shape.
3175 ** All Distance Kernels are symetrical, but that may not always
3176 ** be the case. For example how about a distance from left edges?
3177 ** To work correctly with asymetrical kernels the reflected kernel
3178 ** needs to be applied.
3180 k = &kernel->values[ kernel->width*kernel->height-1 ];
3182 for (v=0; v < (ssize_t) kernel->height; v++) {
3183 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3184 if ( IsNan(*k) ) continue;
3185 Minimize(result.red, (*k)+k_pixels[u].red);
3186 Minimize(result.green, (*k)+k_pixels[u].green);
3187 Minimize(result.blue, (*k)+k_pixels[u].blue);
3188 Minimize(result.alpha, (*k)+k_pixels[u].alpha);
3189 if ( image->colorspace == CMYKColorspace)
3190 Minimize(result.black,(*k)+GetPixelBlack(p_image,k_indexes+u));
3192 k_pixels += virt_width*GetPixelChannels(image);
3196 case UndefinedMorphology:
3198 break; /* Do nothing */
3200 /* Final mathematics of results (combine with original image?)
3202 ** NOTE: Difference Morphology operators Edge* and *Hat could also
3203 ** be done here but works better with iteration as a image difference
3204 ** in the controling function (below). Thicken and Thinning however
3205 ** should be done here so thay can be iterated correctly.
3208 case HitAndMissMorphology:
3209 case ErodeMorphology:
3210 result = min; /* minimum of neighbourhood */
3212 case DilateMorphology:
3213 result = max; /* maximum of neighbourhood */
3215 case ThinningMorphology:
3216 /* subtract pattern match from original */
3217 result.red -= min.red;
3218 result.green -= min.green;
3219 result.blue -= min.blue;
3220 result.black -= min.black;
3221 result.alpha -= min.alpha;
3223 case ThickenMorphology:
3224 /* Add the pattern matchs to the original */
3225 result.red += min.red;
3226 result.green += min.green;
3227 result.blue += min.blue;
3228 result.black += min.black;
3229 result.alpha += min.alpha;
3232 /* result directly calculated or assigned */
3235 /* Assign the resulting pixel values - Clamping Result */
3237 case UndefinedMorphology:
3238 case ConvolveMorphology:
3239 case DilateIntensityMorphology:
3240 case ErodeIntensityMorphology:
3241 break; /* full pixel was directly assigned - not a channel method */
3243 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3244 SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
3245 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3246 SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
3247 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3248 SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
3249 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3250 (image->colorspace == CMYKColorspace))
3251 SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
3252 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3253 (image->matte == MagickTrue))
3254 SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
3257 /* Count up changed pixels */
3258 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q)) ||
3259 (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q)) ||
3260 (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q)) ||
3261 (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q)) ||
3262 ((image->colorspace == CMYKColorspace) &&
3263 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
3264 changed++; /* The pixel was changed in some way! */
3265 p+=GetPixelChannels(image);
3266 q+=GetPixelChannels(morphology_image);
3268 if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3270 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3275 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3276 #pragma omp critical (MagickCore_MorphologyImage)
3278 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3279 if (proceed == MagickFalse)
3283 morphology_view=DestroyCacheView(morphology_view);
3284 image_view=DestroyCacheView(image_view);
3285 return(status ? (ssize_t)changed : -1);
3288 /* This is almost identical to the MorphologyPrimative() function above,
3289 ** but will apply the primitive directly to the image in two passes.
3291 ** That is after each row is 'Sync'ed' into the image, the next row will
3292 ** make use of those values as part of the calculation of the next row.
3293 ** It then repeats, but going in the oppisite (bottom-up) direction.
3295 ** Because of this 'iterative' handling this function can not make use
3296 ** of multi-threaded, parellel processing.
3298 static ssize_t MorphologyPrimitiveDirect(Image *image,
3299 const MorphologyMethod method,const KernelInfo *kernel,
3300 ExceptionInfo *exception)
3323 assert(image != (Image *) NULL);
3324 assert(image->signature == MagickSignature);
3325 assert(kernel != (KernelInfo *) NULL);
3326 assert(kernel->signature == MagickSignature);
3327 assert(exception != (ExceptionInfo *) NULL);
3328 assert(exception->signature == MagickSignature);
3330 /* Some methods (including convolve) needs use a reflected kernel.
3331 * Adjust 'origin' offsets to loop though kernel as a reflection.
3336 case DistanceMorphology:
3337 case VoronoiMorphology:
3338 /* kernel needs to used with reflection about origin */
3339 offx = (ssize_t) kernel->width-offx-1;
3340 offy = (ssize_t) kernel->height-offy-1;
3343 case ?????Morphology:
3344 /* kernel is used as is, without reflection */
3348 assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
3352 /* DO NOT THREAD THIS CODE! */
3353 /* two views into same image (virtual, and actual) */
3354 virt_view=AcquireCacheView(image);
3355 auth_view=AcquireCacheView(image);
3356 virt_width=image->columns+kernel->width-1;
3358 for (y=0; y < (ssize_t) image->rows; y++)
3360 register const Quantum
3372 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3373 ** we read using virtual to get virtual pixel handling, but write back
3374 ** into the same image.
3376 ** Only top half of kernel is processed as we do a single pass downward
3377 ** through the image iterating the distance function as we go.
3379 if (status == MagickFalse)
3381 p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
3383 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3385 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3387 if (status == MagickFalse)
3390 /* offset to origin in 'p'. while 'q' points to it directly */
3391 r = (ssize_t) virt_width*offy + offx;
3393 for (x=0; x < (ssize_t) image->columns; x++)
3401 register const double
3404 register const Quantum
3410 /* Starting Defaults */
3411 GetPixelInfo(image,&result);
3412 SetPixelInfo(image,q,&result);
3413 if ( method != VoronoiMorphology )
3414 result.alpha = QuantumRange - result.alpha;
3417 case DistanceMorphology:
3418 /* Add kernel Value and select the minimum value found. */
3419 k = &kernel->values[ kernel->width*kernel->height-1 ];
3421 for (v=0; v <= (ssize_t) offy; v++) {
3422 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3423 if ( IsNan(*k) ) continue;
3424 Minimize(result.red, (*k)+
3425 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3426 Minimize(result.green, (*k)+
3427 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3428 Minimize(result.blue, (*k)+
3429 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3430 if (image->colorspace == CMYKColorspace)
3431 Minimize(result.black,(*k)+
3432 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3433 Minimize(result.alpha, (*k)+
3434 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3436 k_pixels += virt_width*GetPixelChannels(image);
3438 /* repeat with the just processed pixels of this row */
3439 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3440 k_pixels = q-offx*GetPixelChannels(image);
3441 for (u=0; u < (ssize_t) offx; u++, k--) {
3442 if ( x+u-offx < 0 ) continue; /* off the edge! */
3443 if ( IsNan(*k) ) continue;
3444 Minimize(result.red, (*k)+
3445 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3446 Minimize(result.green, (*k)+
3447 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3448 Minimize(result.blue, (*k)+
3449 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3450 if (image->colorspace == CMYKColorspace)
3451 Minimize(result.black,(*k)+
3452 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3453 Minimize(result.alpha,(*k)+
3454 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3457 case VoronoiMorphology:
3458 /* Apply Distance to 'Matte' channel, coping the closest color.
3460 ** This is experimental, and realy the 'alpha' component should
3461 ** be completely separate 'masking' channel.
3463 k = &kernel->values[ kernel->width*kernel->height-1 ];
3465 for (v=0; v <= (ssize_t) offy; v++) {
3466 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3467 if ( IsNan(*k) ) continue;
3468 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3470 SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3475 k_pixels += virt_width*GetPixelChannels(image);
3477 /* repeat with the just processed pixels of this row */
3478 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3479 k_pixels = q-offx*GetPixelChannels(image);
3480 for (u=0; u < (ssize_t) offx; u++, k--) {
3481 if ( x+u-offx < 0 ) continue; /* off the edge! */
3482 if ( IsNan(*k) ) continue;
3483 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3485 SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3492 /* result directly calculated or assigned */
3495 /* Assign the resulting pixel values - Clamping Result */
3497 case VoronoiMorphology:
3498 SetPixelPixelInfo(image,&result,q);
3501 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3502 SetPixelRed(image,ClampToQuantum(result.red),q);
3503 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3504 SetPixelGreen(image,ClampToQuantum(result.green),q);
3505 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3506 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3507 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3508 (image->colorspace == CMYKColorspace))
3509 SetPixelBlack(image,ClampToQuantum(result.black),q);
3510 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3511 (image->matte == MagickTrue))
3512 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3515 /* Count up changed pixels */
3516 if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q)) ||
3517 (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q)) ||
3518 (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q)) ||
3519 (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q)) ||
3520 ((image->colorspace == CMYKColorspace) &&
3521 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3522 changed++; /* The pixel was changed in some way! */
3524 p+=GetPixelChannels(image); /* increment pixel buffers */
3525 q+=GetPixelChannels(image);
3528 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3530 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3531 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3537 /* Do the reversed pass through the image */
3538 for (y=(ssize_t)image->rows-1; y >= 0; y--)
3540 register const Quantum
3552 if (status == MagickFalse)
3554 /* NOTE read virtual pixels, and authentic pixels, from the same image!
3555 ** we read using virtual to get virtual pixel handling, but write back
3556 ** into the same image.
3558 ** Only the bottom half of the kernel will be processes as we
3561 p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
3562 kernel->y+1,exception);
3563 q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
3565 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3567 if (status == MagickFalse)
3570 /* adjust positions to end of row */
3571 p += (image->columns-1)*GetPixelChannels(image);
3572 q += (image->columns-1)*GetPixelChannels(image);
3574 /* offset to origin in 'p'. while 'q' points to it directly */
3577 for (x=(ssize_t)image->columns-1; x >= 0; x--)
3585 register const double
3588 register const Quantum
3594 /* Default - previously modified pixel */
3595 GetPixelInfo(image,&result);
3596 SetPixelInfo(image,q,&result);
3597 if ( method != VoronoiMorphology )
3598 result.alpha = QuantumRange - result.alpha;
3601 case DistanceMorphology:
3602 /* Add kernel Value and select the minimum value found. */
3603 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3605 for (v=offy; v < (ssize_t) kernel->height; v++) {
3606 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3607 if ( IsNan(*k) ) continue;
3608 Minimize(result.red, (*k)+
3609 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3610 Minimize(result.green, (*k)+
3611 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3612 Minimize(result.blue, (*k)+
3613 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3614 if ( image->colorspace == CMYKColorspace)
3615 Minimize(result.black,(*k)+
3616 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3617 Minimize(result.alpha, (*k)+
3618 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3620 k_pixels += virt_width*GetPixelChannels(image);
3622 /* repeat with the just processed pixels of this row */
3623 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3624 k_pixels = q-offx*GetPixelChannels(image);
3625 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3626 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3627 if ( IsNan(*k) ) continue;
3628 Minimize(result.red, (*k)+
3629 GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
3630 Minimize(result.green, (*k)+
3631 GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
3632 Minimize(result.blue, (*k)+
3633 GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
3634 if ( image->colorspace == CMYKColorspace)
3635 Minimize(result.black, (*k)+
3636 GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
3637 Minimize(result.alpha, (*k)+
3638 GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
3641 case VoronoiMorphology:
3642 /* Apply Distance to 'Matte' channel, coping the closest color.
3644 ** This is experimental, and realy the 'alpha' component should
3645 ** be completely separate 'masking' channel.
3647 k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
3649 for (v=offy; v < (ssize_t) kernel->height; v++) {
3650 for (u=0; u < (ssize_t) kernel->width; u++, k--) {
3651 if ( IsNan(*k) ) continue;
3652 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3654 SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3659 k_pixels += virt_width*GetPixelChannels(image);
3661 /* repeat with the just processed pixels of this row */
3662 k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
3663 k_pixels = q-offx*GetPixelChannels(image);
3664 for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
3665 if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
3666 if ( IsNan(*k) ) continue;
3667 if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
3669 SetPixelInfo(image,k_pixels+u*GetPixelChannels(image),
3676 /* result directly calculated or assigned */
3679 /* Assign the resulting pixel values - Clamping Result */
3681 case VoronoiMorphology:
3682 SetPixelPixelInfo(image,&result,q);
3685 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3686 SetPixelRed(image,ClampToQuantum(result.red),q);
3687 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3688 SetPixelGreen(image,ClampToQuantum(result.green),q);
3689 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3690 SetPixelBlue(image,ClampToQuantum(result.blue),q);
3691 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3692 (image->colorspace == CMYKColorspace))
3693 SetPixelBlack(image,ClampToQuantum(result.black),q);
3694 if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
3695 (image->matte == MagickTrue))
3696 SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
3699 /* Count up changed pixels */
3700 if ( (GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q))
3701 || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q))
3702 || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q))
3703 || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q))
3704 || ((image->colorspace == CMYKColorspace) &&
3705 (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
3706 changed++; /* The pixel was changed in some way! */
3708 p-=GetPixelChannels(image); /* go backward through pixel buffers */
3709 q-=GetPixelChannels(image);
3711 if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
3713 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3714 if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
3720 auth_view=DestroyCacheView(auth_view);
3721 virt_view=DestroyCacheView(virt_view);
3722 return(status ? (ssize_t) changed : -1);
3725 /* Apply a Morphology by calling theabove low level primitive application
3726 ** functions. This function handles any iteration loops, composition or
3727 ** re-iteration of results, and compound morphology methods that is based
3728 ** on multiple low-level (staged) morphology methods.
3730 ** Basically this provides the complex grue between the requested morphology
3731 ** method and raw low-level implementation (above).
3733 MagickPrivate Image *MorphologyApply(const Image *image,
3734 const MorphologyMethod method, const ssize_t iterations,
3735 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3736 ExceptionInfo *exception)
3742 *curr_image, /* Image we are working with or iterating */
3743 *work_image, /* secondary image for primitive iteration */
3744 *save_image, /* saved image - for 'edge' method only */
3745 *rslt_image; /* resultant image - after multi-kernel handling */
3748 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3749 *norm_kernel, /* the current normal un-reflected kernel */
3750 *rflt_kernel, /* the current reflected kernel (if needed) */
3751 *this_kernel; /* the kernel being applied */
3754 primitive; /* the current morphology primitive being applied */
3757 rslt_compose; /* multi-kernel compose method for results to use */
3760 special, /* do we use a direct modify function? */
3761 verbose; /* verbose output of results */
3764 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3765 method_limit, /* maximum number of compound method iterations */
3766 kernel_number, /* Loop 2: the kernel number being applied */
3767 stage_loop, /* Loop 3: primitive loop for compound morphology */
3768 stage_limit, /* how many primitives are in this compound */
3769 kernel_loop, /* Loop 4: iterate the kernel over image */
3770 kernel_limit, /* number of times to iterate kernel */
3771 count, /* total count of primitive steps applied */
3772 kernel_changed, /* total count of changed using iterated kernel */
3773 method_changed; /* total count of changed over method iteration */
3776 changed; /* number pixels changed by last primitive operation */
3781 assert(image != (Image *) NULL);
3782 assert(image->signature == MagickSignature);
3783 assert(kernel != (KernelInfo *) NULL);
3784 assert(kernel->signature == MagickSignature);
3785 assert(exception != (ExceptionInfo *) NULL);
3786 assert(exception->signature == MagickSignature);
3788 count = 0; /* number of low-level morphology primitives performed */
3789 if ( iterations == 0 )
3790 return((Image *)NULL); /* null operation - nothing to do! */
3792 kernel_limit = (size_t) iterations;
3793 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3794 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3796 verbose = IsMagickTrue(GetImageArtifact(image,"verbose"));
3798 /* initialise for cleanup */
3799 curr_image = (Image *) image;
3800 curr_compose = image->compose;
3801 (void) curr_compose;
3802 work_image = save_image = rslt_image = (Image *) NULL;
3803 reflected_kernel = (KernelInfo *) NULL;
3805 /* Initialize specific methods
3806 * + which loop should use the given iteratations
3807 * + how many primitives make up the compound morphology
3808 * + multi-kernel compose method to use (by default)
3810 method_limit = 1; /* just do method once, unless otherwise set */
3811 stage_limit = 1; /* assume method is not a compound */
3812 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3813 rslt_compose = compose; /* and we are composing multi-kernels as given */
3815 case SmoothMorphology: /* 4 primitive compound morphology */
3818 case OpenMorphology: /* 2 primitive compound morphology */
3819 case OpenIntensityMorphology:
3820 case TopHatMorphology:
3821 case CloseMorphology:
3822 case CloseIntensityMorphology:
3823 case BottomHatMorphology:
3824 case EdgeMorphology:
3827 case HitAndMissMorphology:
3828 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3830 case ThinningMorphology:
3831 case ThickenMorphology:
3832 method_limit = kernel_limit; /* iterate the whole method */
3833 kernel_limit = 1; /* do not do kernel iteration */
3835 case DistanceMorphology:
3836 case VoronoiMorphology:
3837 special = MagickTrue;
3843 /* Apply special methods with special requirments
3844 ** For example, single run only, or post-processing requirements
3846 if ( special == MagickTrue )
3848 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3849 if (rslt_image == (Image *) NULL)
3851 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3854 changed = MorphologyPrimitiveDirect(rslt_image, method,
3857 if ( verbose == MagickTrue )
3858 (void) (void) FormatLocaleFile(stderr,
3859 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3860 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3861 1.0,0.0,1.0, (double) changed);
3866 if ( method == VoronoiMorphology ) {
3867 /* Preserve the alpha channel of input image - but turned off */
3868 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3870 (void) CompositeImage(rslt_image, CopyOpacityCompositeOp, image, 0, 0);
3871 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3877 /* Handle user (caller) specified multi-kernel composition method */
3878 if ( compose != UndefinedCompositeOp )
3879 rslt_compose = compose; /* override default composition for method */
3880 if ( rslt_compose == UndefinedCompositeOp )
3881 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3883 /* Some methods require a reflected kernel to use with primitives.
3884 * Create the reflected kernel for those methods. */
3886 case CorrelateMorphology:
3887 case CloseMorphology:
3888 case CloseIntensityMorphology:
3889 case BottomHatMorphology:
3890 case SmoothMorphology:
3891 reflected_kernel = CloneKernelInfo(kernel);
3892 if (reflected_kernel == (KernelInfo *) NULL)
3894 RotateKernelInfo(reflected_kernel,180);
3900 /* Loops around more primitive morpholgy methods
3901 ** erose, dilate, open, close, smooth, edge, etc...
3903 /* Loop 1: iterate the compound method */
3906 while ( method_loop < method_limit && method_changed > 0 ) {
3910 /* Loop 2: iterate over each kernel in a multi-kernel list */
3911 norm_kernel = (KernelInfo *) kernel;
3912 this_kernel = (KernelInfo *) kernel;
3913 rflt_kernel = reflected_kernel;
3916 while ( norm_kernel != NULL ) {
3918 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3919 stage_loop = 0; /* the compound morphology stage number */
3920 while ( stage_loop < stage_limit ) {
3921 stage_loop++; /* The stage of the compound morphology */
3923 /* Select primitive morphology for this stage of compound method */
3924 this_kernel = norm_kernel; /* default use unreflected kernel */
3925 primitive = method; /* Assume method is a primitive */
3927 case ErodeMorphology: /* just erode */
3928 case EdgeInMorphology: /* erode and image difference */
3929 primitive = ErodeMorphology;
3931 case DilateMorphology: /* just dilate */
3932 case EdgeOutMorphology: /* dilate and image difference */
3933 primitive = DilateMorphology;
3935 case OpenMorphology: /* erode then dialate */
3936 case TopHatMorphology: /* open and image difference */
3937 primitive = ErodeMorphology;
3938 if ( stage_loop == 2 )
3939 primitive = DilateMorphology;
3941 case OpenIntensityMorphology:
3942 primitive = ErodeIntensityMorphology;
3943 if ( stage_loop == 2 )
3944 primitive = DilateIntensityMorphology;
3946 case CloseMorphology: /* dilate, then erode */
3947 case BottomHatMorphology: /* close and image difference */
3948 this_kernel = rflt_kernel; /* use the reflected kernel */
3949 primitive = DilateMorphology;
3950 if ( stage_loop == 2 )
3951 primitive = ErodeMorphology;
3953 case CloseIntensityMorphology:
3954 this_kernel = rflt_kernel; /* use the reflected kernel */
3955 primitive = DilateIntensityMorphology;
3956 if ( stage_loop == 2 )
3957 primitive = ErodeIntensityMorphology;
3959 case SmoothMorphology: /* open, close */
3960 switch ( stage_loop ) {
3961 case 1: /* start an open method, which starts with Erode */
3962 primitive = ErodeMorphology;
3964 case 2: /* now Dilate the Erode */
3965 primitive = DilateMorphology;
3967 case 3: /* Reflect kernel a close */
3968 this_kernel = rflt_kernel; /* use the reflected kernel */
3969 primitive = DilateMorphology;
3971 case 4: /* Finish the Close */
3972 this_kernel = rflt_kernel; /* use the reflected kernel */
3973 primitive = ErodeMorphology;
3977 case EdgeMorphology: /* dilate and erode difference */
3978 primitive = DilateMorphology;
3979 if ( stage_loop == 2 ) {
3980 save_image = curr_image; /* save the image difference */
3981 curr_image = (Image *) image;
3982 primitive = ErodeMorphology;
3985 case CorrelateMorphology:
3986 /* A Correlation is a Convolution with a reflected kernel.
3987 ** However a Convolution is a weighted sum using a reflected
3988 ** kernel. It may seem stange to convert a Correlation into a
3989 ** Convolution as the Correlation is the simplier method, but
3990 ** Convolution is much more commonly used, and it makes sense to
3991 ** implement it directly so as to avoid the need to duplicate the
3992 ** kernel when it is not required (which is typically the
3995 this_kernel = rflt_kernel; /* use the reflected kernel */
3996 primitive = ConvolveMorphology;
4001 assert( this_kernel != (KernelInfo *) NULL );
4003 /* Extra information for debugging compound operations */
4004 if ( verbose == MagickTrue ) {
4005 if ( stage_limit > 1 )
4006 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
4007 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
4008 method_loop,(double) stage_loop);
4009 else if ( primitive != method )
4010 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
4011 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
4017 /* Loop 4: Iterate the kernel with primitive */
4021 while ( kernel_loop < kernel_limit && changed > 0 ) {
4022 kernel_loop++; /* the iteration of this kernel */
4024 /* Create a clone as the destination image, if not yet defined */
4025 if ( work_image == (Image *) NULL )
4027 work_image=CloneImage(image,0,0,MagickTrue,exception);
4028 if (work_image == (Image *) NULL)
4030 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
4032 /* work_image->type=image->type; ??? */
4035 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
4037 changed = MorphologyPrimitive(curr_image, work_image, primitive,
4038 this_kernel, bias, exception);
4040 if ( verbose == MagickTrue ) {
4041 if ( kernel_loop > 1 )
4042 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
4043 (void) (void) FormatLocaleFile(stderr,
4044 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
4045 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
4046 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
4047 (double) (method_loop+kernel_loop-1),(double) kernel_number,
4048 (double) count,(double) changed);
4052 kernel_changed += changed;
4053 method_changed += changed;
4055 /* prepare next loop */
4056 { Image *tmp = work_image; /* swap images for iteration */
4057 work_image = curr_image;
4060 if ( work_image == image )
4061 work_image = (Image *) NULL; /* replace input 'image' */
4063 } /* End Loop 4: Iterate the kernel with primitive */
4065 if ( verbose == MagickTrue && kernel_changed != (size_t)changed )
4066 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
4067 if ( verbose == MagickTrue && stage_loop < stage_limit )
4068 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
4071 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
4072 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
4073 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
4074 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
4075 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
4078 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
4080 /* Final Post-processing for some Compound Methods
4082 ** The removal of any 'Sync' channel flag in the Image Compositon
4083 ** below ensures the methematical compose method is applied in a
4084 ** purely mathematical way, and only to the selected channels.
4085 ** Turn off SVG composition 'alpha blending'.
4088 case EdgeOutMorphology:
4089 case EdgeInMorphology:
4090 case TopHatMorphology:
4091 case BottomHatMorphology:
4092 if ( verbose == MagickTrue )
4093 (void) FormatLocaleFile(stderr, "\n%s: Difference with original image",
4094 CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4095 curr_image->sync=MagickFalse;
4096 (void) CompositeImage(curr_image,DifferenceCompositeOp,image,0,0);
4098 case EdgeMorphology:
4099 if ( verbose == MagickTrue )
4100 (void) FormatLocaleFile(stderr, "\n%s: Difference of Dilate and Erode",
4101 CommandOptionToMnemonic(MagickMorphologyOptions, method) );
4102 curr_image->sync=MagickFalse;
4103 (void) CompositeImage(curr_image,DifferenceCompositeOp,save_image,0,
4105 save_image = DestroyImage(save_image); /* finished with save image */
4111 /* multi-kernel handling: re-iterate, or compose results */
4112 if ( kernel->next == (KernelInfo *) NULL )
4113 rslt_image = curr_image; /* just return the resulting image */
4114 else if ( rslt_compose == NoCompositeOp )
4115 { if ( verbose == MagickTrue ) {
4116 if ( this_kernel->next != (KernelInfo *) NULL )
4117 (void) FormatLocaleFile(stderr, " (re-iterate)");
4119 (void) FormatLocaleFile(stderr, " (done)");
4121 rslt_image = curr_image; /* return result, and re-iterate */
4123 else if ( rslt_image == (Image *) NULL)
4124 { if ( verbose == MagickTrue )
4125 (void) FormatLocaleFile(stderr, " (save for compose)");
4126 rslt_image = curr_image;
4127 curr_image = (Image *) image; /* continue with original image */
4130 { /* Add the new 'current' result to the composition
4132 ** The removal of any 'Sync' channel flag in the Image Compositon
4133 ** below ensures the methematical compose method is applied in a
4134 ** purely mathematical way, and only to the selected channels.
4135 ** IE: Turn off SVG composition 'alpha blending'.
4137 if ( verbose == MagickTrue )
4138 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4139 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4140 rslt_image->sync=MagickFalse;
4141 (void) CompositeImage(rslt_image, rslt_compose, curr_image, 0, 0);
4142 curr_image = DestroyImage(curr_image);
4143 curr_image = (Image *) image; /* continue with original image */
4145 if ( verbose == MagickTrue )
4146 (void) FormatLocaleFile(stderr, "\n");
4148 /* loop to the next kernel in a multi-kernel list */
4149 norm_kernel = norm_kernel->next;
4150 if ( rflt_kernel != (KernelInfo *) NULL )
4151 rflt_kernel = rflt_kernel->next;
4153 } /* End Loop 2: Loop over each kernel */
4155 } /* End Loop 1: compound method interation */
4159 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4161 if ( curr_image == rslt_image )
4162 curr_image = (Image *) NULL;
4163 if ( rslt_image != (Image *) NULL )
4164 rslt_image = DestroyImage(rslt_image);
4166 if ( curr_image == rslt_image || curr_image == image )
4167 curr_image = (Image *) NULL;
4168 if ( curr_image != (Image *) NULL )
4169 curr_image = DestroyImage(curr_image);
4170 if ( work_image != (Image *) NULL )
4171 work_image = DestroyImage(work_image);
4172 if ( save_image != (Image *) NULL )
4173 save_image = DestroyImage(save_image);
4174 if ( reflected_kernel != (KernelInfo *) NULL )
4175 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4185 % M o r p h o l o g y I m a g e %
4189 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4191 % MorphologyImage() applies a user supplied kernel to the image
4192 % according to the given mophology method.
4194 % This function applies any and all user defined settings before calling
4195 % the above internal function MorphologyApply().
4197 % User defined settings include...
4198 % * Output Bias for Convolution and correlation ("-bias")
4199 % * Kernel Scale/normalize settings ("-set 'option:convolve:scale'")
4200 % This can also includes the addition of a scaled unity kernel.
4201 % * Show Kernel being applied ("-set option:showkernel 1")
4203 % The format of the MorphologyImage method is:
4205 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4206 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4208 % Image *MorphologyImage(const Image *image, const ChannelType
4209 % channel,MorphologyMethod method,const ssize_t iterations,
4210 % KernelInfo *kernel,ExceptionInfo *exception)
4212 % A description of each parameter follows:
4214 % o image: the image.
4216 % o method: the morphology method to be applied.
4218 % o iterations: apply the operation this many times (or no change).
4219 % A value of -1 means loop until no change found.
4220 % How this is applied may depend on the morphology method.
4221 % Typically this is a value of 1.
4223 % o kernel: An array of double representing the morphology kernel.
4224 % Warning: kernel may be normalized for the Convolve method.
4226 % o exception: return any errors or warnings in this structure.
4229 MagickExport Image *MorphologyImage(const Image *image,
4230 const MorphologyMethod method,const ssize_t iterations,
4231 const KernelInfo *kernel,ExceptionInfo *exception)
4243 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4244 * This is done BEFORE the ShowKernelInfo() function is called so that
4245 * users can see the results of the 'option:convolve:scale' option.
4247 curr_kernel = (KernelInfo *) kernel;
4248 if ( method == ConvolveMorphology || method == CorrelateMorphology )
4252 artifact = GetImageArtifact(image,"convolve:scale");
4253 if ( artifact != (const char *)NULL ) {
4254 if ( curr_kernel == kernel )
4255 curr_kernel = CloneKernelInfo(kernel);
4256 if (curr_kernel == (KernelInfo *) NULL) {
4257 curr_kernel=DestroyKernelInfo(curr_kernel);
4258 return((Image *) NULL);
4260 ScaleGeometryKernelInfo(curr_kernel, artifact);
4264 /* display the (normalized) kernel via stderr */
4265 if ( IsMagickTrue(GetImageArtifact(image,"showkernel"))
4266 || IsMagickTrue(GetImageArtifact(image,"convolve:showkernel"))
4267 || IsMagickTrue(GetImageArtifact(image,"morphology:showkernel")) )
4268 ShowKernelInfo(curr_kernel);
4270 /* Override the default handling of multi-kernel morphology results
4271 * If 'Undefined' use the default method
4272 * If 'None' (default for 'Convolve') re-iterate previous result
4273 * Otherwise merge resulting images using compose method given.
4274 * Default for 'HitAndMiss' is 'Lighten'.
4278 artifact = GetImageArtifact(image,"morphology:compose");
4279 compose = UndefinedCompositeOp; /* use default for method */
4280 if ( artifact != (const char *) NULL)
4281 compose=(CompositeOperator) ParseCommandOption(MagickComposeOptions,
4282 MagickFalse,artifact);
4284 /* Apply the Morphology */
4285 morphology_image = MorphologyApply(image, method, iterations,
4286 curr_kernel, compose, image->bias, exception);
4288 /* Cleanup and Exit */
4289 if ( curr_kernel != kernel )
4290 curr_kernel=DestroyKernelInfo(curr_kernel);
4291 return(morphology_image);
4295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4299 + R o t a t e K e r n e l I n f o %
4303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4305 % RotateKernelInfo() rotates the kernel by the angle given.
4307 % Currently it is restricted to 90 degree angles, of either 1D kernels
4308 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4309 % It will ignore usless rotations for specific 'named' built-in kernels.
4311 % The format of the RotateKernelInfo method is:
4313 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4315 % A description of each parameter follows:
4317 % o kernel: the Morphology/Convolution kernel
4319 % o angle: angle to rotate in degrees
4321 % This function is currently internal to this module only, but can be exported
4322 % to other modules if needed.
4324 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4326 /* angle the lower kernels first */
4327 if ( kernel->next != (KernelInfo *) NULL)
4328 RotateKernelInfo(kernel->next, angle);
4330 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4332 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4335 /* Modulus the angle */
4336 angle = fmod(angle, 360.0);
4340 if ( 337.5 < angle || angle <= 22.5 )
4341 return; /* Near zero angle - no change! - At least not at this time */
4343 /* Handle special cases */
4344 switch (kernel->type) {
4345 /* These built-in kernels are cylindrical kernels, rotating is useless */
4346 case GaussianKernel:
4351 case LaplacianKernel:
4352 case ChebyshevKernel:
4353 case ManhattanKernel:
4354 case EuclideanKernel:
4357 /* These may be rotatable at non-90 angles in the future */
4358 /* but simply rotating them in multiples of 90 degrees is useless */
4365 /* These only allows a +/-90 degree rotation (by transpose) */
4366 /* A 180 degree rotation is useless */
4368 if ( 135.0 < angle && angle <= 225.0 )
4370 if ( 225.0 < angle && angle <= 315.0 )
4377 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4378 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4380 if ( kernel->width == 3 && kernel->height == 3 )
4381 { /* Rotate a 3x3 square by 45 degree angle */
4382 MagickRealType t = kernel->values[0];
4383 kernel->values[0] = kernel->values[3];
4384 kernel->values[3] = kernel->values[6];
4385 kernel->values[6] = kernel->values[7];
4386 kernel->values[7] = kernel->values[8];
4387 kernel->values[8] = kernel->values[5];
4388 kernel->values[5] = kernel->values[2];
4389 kernel->values[2] = kernel->values[1];
4390 kernel->values[1] = t;
4391 /* rotate non-centered origin */
4392 if ( kernel->x != 1 || kernel->y != 1 ) {
4394 x = (ssize_t) kernel->x-1;
4395 y = (ssize_t) kernel->y-1;
4396 if ( x == y ) x = 0;
4397 else if ( x == 0 ) x = -y;
4398 else if ( x == -y ) y = 0;
4399 else if ( y == 0 ) y = x;
4400 kernel->x = (ssize_t) x+1;
4401 kernel->y = (ssize_t) y+1;
4403 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4404 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4407 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4409 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4411 if ( kernel->width == 1 || kernel->height == 1 )
4412 { /* Do a transpose of a 1 dimensional kernel,
4413 ** which results in a fast 90 degree rotation of some type.
4417 t = (ssize_t) kernel->width;
4418 kernel->width = kernel->height;
4419 kernel->height = (size_t) t;
4421 kernel->x = kernel->y;
4423 if ( kernel->width == 1 ) {
4424 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4425 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4427 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4428 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4431 else if ( kernel->width == kernel->height )
4432 { /* Rotate a square array of values by 90 degrees */
4435 register MagickRealType
4438 for( i=0, x=kernel->width-1; i<=x; i++, x--)
4439 for( j=0, y=kernel->height-1; j<y; j++, y--)
4440 { t = k[i+j*kernel->width];
4441 k[i+j*kernel->width] = k[j+x*kernel->width];
4442 k[j+x*kernel->width] = k[x+y*kernel->width];
4443 k[x+y*kernel->width] = k[y+i*kernel->width];
4444 k[y+i*kernel->width] = t;
4447 /* rotate the origin - relative to center of array */
4448 { register ssize_t x,y;
4449 x = (ssize_t) (kernel->x*2-kernel->width+1);
4450 y = (ssize_t) (kernel->y*2-kernel->height+1);
4451 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4452 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4454 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4455 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4458 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4460 if ( 135.0 < angle && angle <= 225.0 )
4462 /* For a 180 degree rotation - also know as a reflection
4463 * This is actually a very very common operation!
4464 * Basically all that is needed is a reversal of the kernel data!
4465 * And a reflection of the origon
4473 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
4474 t=k[i], k[i]=k[j], k[j]=t;
4476 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4477 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4478 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4479 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4481 /* At this point angle should at least between -45 (315) and +45 degrees
4482 * In the future some form of non-orthogonal angled rotates could be
4483 * performed here, posibily with a linear kernel restriction.
4490 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4494 % S c a l e G e o m e t r y K e r n e l I n f o %
4498 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4500 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4501 % provided as a "-set option:convolve:scale {geometry}" user setting,
4502 % and modifies the kernel according to the parsed arguments of that setting.
4504 % The first argument (and any normalization flags) are passed to
4505 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4506 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4507 % into the scaled/normalized kernel.
4509 % The format of the ScaleGeometryKernelInfo method is:
4511 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4512 % const double scaling_factor,const MagickStatusType normalize_flags)
4514 % A description of each parameter follows:
4516 % o kernel: the Morphology/Convolution kernel to modify
4519 % The geometry string to parse, typically from the user provided
4520 % "-set option:convolve:scale {geometry}" setting.
4523 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4524 const char *geometry)
4531 SetGeometryInfo(&args);
4532 flags = (GeometryFlags) ParseGeometry(geometry, &args);
4535 /* For Debugging Geometry Input */
4536 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4537 flags, args.rho, args.sigma, args.xi, args.psi );
4540 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4541 args.rho *= 0.01, args.sigma *= 0.01;
4543 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4545 if ( (flags & SigmaValue) == 0 )
4548 /* Scale/Normalize the input kernel */
4549 ScaleKernelInfo(kernel, args.rho, flags);
4551 /* Add Unity Kernel, for blending with original */
4552 if ( (flags & SigmaValue) != 0 )
4553 UnityAddKernelInfo(kernel, args.sigma);
4558 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4562 % S c a l e K e r n e l I n f o %
4566 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4568 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4569 % without normalization of the sum of the kernel values (as per given flags).
4571 % By default (no flags given) the values within the kernel is scaled
4572 % directly using given scaling factor without change.
4574 % If either of the two 'normalize_flags' are given the kernel will first be
4575 % normalized and then further scaled by the scaling factor value given.
4577 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4578 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4579 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4580 % non-HDRI versions of IM this may cause images to have any negative results
4581 % clipped, unless some 'bias' is used.
4583 % More specifically. Kernels which only contain positive values (such as a
4584 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4585 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4587 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4588 % the kernel will be scaled by the absolute of the sum of kernel values, so
4589 % that it will generally fall within the +/- 1.0 range.
4591 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4592 % will be scaled by just the sum of the postive values, so that its output
4593 % range will again fall into the +/- 1.0 range.
4595 % For special kernels designed for locating shapes using 'Correlate', (often
4596 % only containing +1 and -1 values, representing foreground/brackground
4597 % matching) a special normalization method is provided to scale the positive
4598 % values separately to those of the negative values, so the kernel will be
4599 % forced to become a zero-sum kernel better suited to such searches.
4601 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4602 % attributes within the kernel structure have been correctly set during the
4605 % NOTE: The values used for 'normalize_flags' have been selected specifically
4606 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4607 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4609 % The format of the ScaleKernelInfo method is:
4611 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4612 % const MagickStatusType normalize_flags )
4614 % A description of each parameter follows:
4616 % o kernel: the Morphology/Convolution kernel
4619 % multiply all values (after normalization) by this factor if not
4620 % zero. If the kernel is normalized regardless of any flags.
4622 % o normalize_flags:
4623 % GeometryFlags defining normalization method to use.
4624 % specifically: NormalizeValue, CorrelateNormalizeValue,
4625 % and/or PercentValue
4628 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4629 const double scaling_factor,const GeometryFlags normalize_flags)
4638 /* do the other kernels in a multi-kernel list first */
4639 if ( kernel->next != (KernelInfo *) NULL)
4640 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4642 /* Normalization of Kernel */
4644 if ( (normalize_flags&NormalizeValue) != 0 ) {
4645 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
4646 /* non-zero-summing kernel (generally positive) */
4647 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4649 /* zero-summing kernel */
4650 pos_scale = kernel->positive_range;
4652 /* Force kernel into a normalized zero-summing kernel */
4653 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4654 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
4655 ? kernel->positive_range : 1.0;
4656 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
4657 ? -kernel->negative_range : 1.0;
4660 neg_scale = pos_scale;
4662 /* finialize scaling_factor for positive and negative components */
4663 pos_scale = scaling_factor/pos_scale;
4664 neg_scale = scaling_factor/neg_scale;
4666 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4667 if ( ! IsNan(kernel->values[i]) )
4668 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4670 /* convolution output range */
4671 kernel->positive_range *= pos_scale;
4672 kernel->negative_range *= neg_scale;
4673 /* maximum and minimum values in kernel */
4674 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4675 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4677 /* swap kernel settings if user's scaling factor is negative */
4678 if ( scaling_factor < MagickEpsilon ) {
4680 t = kernel->positive_range;
4681 kernel->positive_range = kernel->negative_range;
4682 kernel->negative_range = t;
4683 t = kernel->maximum;
4684 kernel->maximum = kernel->minimum;
4685 kernel->minimum = 1;
4692 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4696 % S h o w K e r n e l I n f o %
4700 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4702 % ShowKernelInfo() outputs the details of the given kernel defination to
4703 % standard error, generally due to a users 'showkernel' option request.
4705 % The format of the ShowKernel method is:
4707 % void ShowKernelInfo(const KernelInfo *kernel)
4709 % A description of each parameter follows:
4711 % o kernel: the Morphology/Convolution kernel
4714 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4722 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4724 (void) FormatLocaleFile(stderr, "Kernel");
4725 if ( kernel->next != (KernelInfo *) NULL )
4726 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4727 (void) FormatLocaleFile(stderr, " \"%s",
4728 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4729 if ( fabs(k->angle) > MagickEpsilon )
4730 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4731 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4732 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4733 (void) FormatLocaleFile(stderr,
4734 " with values from %.*lg to %.*lg\n",
4735 GetMagickPrecision(), k->minimum,
4736 GetMagickPrecision(), k->maximum);
4737 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4738 GetMagickPrecision(), k->negative_range,
4739 GetMagickPrecision(), k->positive_range);
4740 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4741 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4742 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4743 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4745 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4746 GetMagickPrecision(), k->positive_range+k->negative_range);
4747 for (i=v=0; v < k->height; v++) {
4748 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4749 for (u=0; u < k->width; u++, i++)
4750 if ( IsNan(k->values[i]) )
4751 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4753 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4754 GetMagickPrecision(), k->values[i]);
4755 (void) FormatLocaleFile(stderr,"\n");
4761 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4765 % U n i t y A d d K e r n a l I n f o %
4769 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4771 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4772 % to the given pre-scaled and normalized Kernel. This in effect adds that
4773 % amount of the original image into the resulting convolution kernel. This
4774 % value is usually provided by the user as a percentage value in the
4775 % 'convolve:scale' setting.
4777 % The resulting effect is to convert the defined kernels into blended
4778 % soft-blurs, unsharp kernels or into sharpening kernels.
4780 % The format of the UnityAdditionKernelInfo method is:
4782 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4784 % A description of each parameter follows:
4786 % o kernel: the Morphology/Convolution kernel
4789 % scaling factor for the unity kernel to be added to
4793 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4796 /* do the other kernels in a multi-kernel list first */
4797 if ( kernel->next != (KernelInfo *) NULL)
4798 UnityAddKernelInfo(kernel->next, scale);
4800 /* Add the scaled unity kernel to the existing kernel */
4801 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4802 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4808 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4812 % Z e r o K e r n e l N a n s %
4816 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4818 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4819 % the kernel with a zero value. This is typically done when the kernel will
4820 % be used in special hardware (GPU) convolution processors, to simply
4823 % The format of the ZeroKernelNans method is:
4825 % void ZeroKernelNans (KernelInfo *kernel)
4827 % A description of each parameter follows:
4829 % o kernel: the Morphology/Convolution kernel
4832 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4837 /* do the other kernels in a multi-kernel list first */
4838 if ( kernel->next != (KernelInfo *) NULL)
4839 ZeroKernelNans(kernel->next);
4841 for (i=0; i < (kernel->width*kernel->height); i++)
4842 if ( IsNan(kernel->values[i]) )
4843 kernel->values[i] = 0.0;