2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 % M M OOO RRRR PPPP H H OOO L OOO GGGG Y Y %
7 % MM MM O O R R P P H H O O L O O G Y Y %
8 % M M M O O RRRR PPPP HHHHH O O L O O G GGG Y %
9 % M M O O R R P H H O O L O O G G Y %
10 % M M OOO R R P H H OOO LLLLL OOO GGG Y %
13 % MagickCore Morphology Methods %
20 % Copyright 1999-2015 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % http://www.imagemagick.org/script/license.php %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 % Morpology is the the application of various kernels, of any size and even
37 % shape, to a image in various ways (typically binary, but not always).
39 % Convolution (weighted sum or average) is just one specific type of
40 % morphology. Just one that is very common for image bluring and sharpening
41 % effects. Not only 2D Gaussian blurring, but also 2-pass 1D Blurring.
43 % This module provides not only a general morphology function, and the ability
44 % to apply more advanced or iterative morphologies, but also functions for the
45 % generation of many different types of kernel arrays from user supplied
46 % arguments. Prehaps even the generation of a kernel from a small image.
52 #include "MagickCore/studio.h"
53 #include "MagickCore/artifact.h"
54 #include "MagickCore/cache-view.h"
55 #include "MagickCore/channel.h"
56 #include "MagickCore/color-private.h"
57 #include "MagickCore/enhance.h"
58 #include "MagickCore/exception.h"
59 #include "MagickCore/exception-private.h"
60 #include "MagickCore/gem.h"
61 #include "MagickCore/gem-private.h"
62 #include "MagickCore/hashmap.h"
63 #include "MagickCore/image.h"
64 #include "MagickCore/image-private.h"
65 #include "MagickCore/list.h"
66 #include "MagickCore/magick.h"
67 #include "MagickCore/memory_.h"
68 #include "MagickCore/memory-private.h"
69 #include "MagickCore/monitor-private.h"
70 #include "MagickCore/morphology.h"
71 #include "MagickCore/morphology-private.h"
72 #include "MagickCore/option.h"
73 #include "MagickCore/pixel-accessor.h"
74 #include "MagickCore/pixel-private.h"
75 #include "MagickCore/prepress.h"
76 #include "MagickCore/quantize.h"
77 #include "MagickCore/resource_.h"
78 #include "MagickCore/registry.h"
79 #include "MagickCore/semaphore.h"
80 #include "MagickCore/splay-tree.h"
81 #include "MagickCore/statistic.h"
82 #include "MagickCore/string_.h"
83 #include "MagickCore/string-private.h"
84 #include "MagickCore/thread-private.h"
85 #include "MagickCore/token.h"
86 #include "MagickCore/utility.h"
87 #include "MagickCore/utility-private.h"
90 Other global definitions used by module.
92 static inline double MagickMin(const double x,const double y)
94 return( x < y ? x : y);
96 static inline double MagickMax(const double x,const double y)
98 return( x > y ? x : y);
100 #define Minimize(assign,value) assign=MagickMin(assign,value)
101 #define Maximize(assign,value) assign=MagickMax(assign,value)
103 /* Integer Factorial Function - for a Binomial kernel */
105 static inline size_t fact(size_t n)
108 for(f=1, l=2; l <= n; f=f*l, l++);
111 #elif 1 /* glibc floating point alternatives */
112 #define fact(n) ((size_t)tgamma((double)n+1))
114 #define fact(n) ((size_t)lgamma((double)n+1))
118 /* Currently these are only internal to this module */
120 CalcKernelMetaData(KernelInfo *),
121 ExpandMirrorKernelInfo(KernelInfo *),
122 ExpandRotateKernelInfo(KernelInfo *, const double),
123 RotateKernelInfo(KernelInfo *, double);
126 /* Quick function to find last kernel in a kernel list */
127 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
129 while (kernel->next != (KernelInfo *) NULL)
135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
139 % A c q u i r e K e r n e l I n f o %
143 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145 % AcquireKernelInfo() takes the given string (generally supplied by the
146 % user) and converts it into a Morphology/Convolution Kernel. This allows
147 % users to specify a kernel from a number of pre-defined kernels, or to fully
148 % specify their own kernel for a specific Convolution or Morphology
151 % The kernel so generated can be any rectangular array of floating point
152 % values (doubles) with the 'control point' or 'pixel being affected'
153 % anywhere within that array of values.
155 % Previously IM was restricted to a square of odd size using the exact
156 % center as origin, this is no longer the case, and any rectangular kernel
157 % with any value being declared the origin. This in turn allows the use of
158 % highly asymmetrical kernels.
160 % The floating point values in the kernel can also include a special value
161 % known as 'nan' or 'not a number' to indicate that this value is not part
162 % of the kernel array. This allows you to shaped the kernel within its
163 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
164 % shape. However at least one non-nan value must be provided for correct
165 % working of a kernel.
167 % The returned kernel should be freed using the DestroyKernelInfo() when you
168 % are finished with it. Do not free this memory yourself.
170 % Input kernel defintion strings can consist of any of three types.
173 % Select from one of the built in kernels, using the name and
174 % geometry arguments supplied. See AcquireKernelBuiltIn()
176 % "WxH[+X+Y][@><]:num, num, num ..."
177 % a kernel of size W by H, with W*H floating point numbers following.
178 % the 'center' can be optionally be defined at +X+Y (such that +0+0
179 % is top left corner). If not defined the pixel in the center, for
180 % odd sizes, or to the immediate top or left of center for even sizes
181 % is automatically selected.
183 % "num, num, num, num, ..."
184 % list of floating point numbers defining an 'old style' odd sized
185 % square kernel. At least 9 values should be provided for a 3x3
186 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
187 % Values can be space or comma separated. This is not recommended.
189 % You can define a 'list of kernels' which can be used by some morphology
190 % operators A list is defined as a semi-colon separated list kernels.
192 % " kernel ; kernel ; kernel ; "
194 % Any extra ';' characters, at start, end or between kernel defintions are
197 % The special flags will expand a single kernel, into a list of rotated
198 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
199 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
200 % The '<' also exands using 90-degree rotates, but giving a 180-degree
201 % reflected kernel before the +/- 90-degree rotations, which can be important
202 % for Thinning operations.
204 % Note that 'name' kernels will start with an alphabetic character while the
205 % new kernel specification has a ':' character in its specification string.
206 % If neither is the case, it is assumed an old style of a simple list of
207 % numbers generating a odd-sized square kernel has been given.
209 % The format of the AcquireKernal method is:
211 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
213 % A description of each parameter follows:
215 % o kernel_string: the Morphology/Convolution kernel wanted.
219 /* This was separated so that it could be used as a separate
220 ** array input handling function, such as for -color-matrix
222 static KernelInfo *ParseKernelArray(const char *kernel_string)
228 token[MaxTextExtent];
238 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
246 kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
247 if (kernel == (KernelInfo *)NULL)
249 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
250 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
251 kernel->negative_range = kernel->positive_range = 0.0;
252 kernel->type = UserDefinedKernel;
253 kernel->next = (KernelInfo *) NULL;
254 kernel->signature = MagickSignature;
255 if (kernel_string == (const char *) NULL)
258 /* find end of this specific kernel definition string */
259 end = strchr(kernel_string, ';');
260 if ( end == (char *) NULL )
261 end = strchr(kernel_string, '\0');
263 /* clear flags - for Expanding kernel lists thorugh rotations */
266 /* Has a ':' in argument - New user kernel specification
267 FUTURE: this split on ':' could be done by StringToken()
269 p = strchr(kernel_string, ':');
270 if ( p != (char *) NULL && p < end)
272 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
273 memcpy(token, kernel_string, (size_t) (p-kernel_string));
274 token[p-kernel_string] = '\0';
275 SetGeometryInfo(&args);
276 flags = ParseGeometry(token, &args);
278 /* Size handling and checks of geometry settings */
279 if ( (flags & WidthValue) == 0 ) /* if no width then */
280 args.rho = args.sigma; /* then width = height */
281 if ( args.rho < 1.0 ) /* if width too small */
282 args.rho = 1.0; /* then width = 1 */
283 if ( args.sigma < 1.0 ) /* if height too small */
284 args.sigma = args.rho; /* then height = width */
285 kernel->width = (size_t)args.rho;
286 kernel->height = (size_t)args.sigma;
288 /* Offset Handling and Checks */
289 if ( args.xi < 0.0 || args.psi < 0.0 )
290 return(DestroyKernelInfo(kernel));
291 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
292 : (ssize_t) (kernel->width-1)/2;
293 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
294 : (ssize_t) (kernel->height-1)/2;
295 if ( kernel->x >= (ssize_t) kernel->width ||
296 kernel->y >= (ssize_t) kernel->height )
297 return(DestroyKernelInfo(kernel));
299 p++; /* advance beyond the ':' */
302 { /* ELSE - Old old specification, forming odd-square kernel */
303 /* count up number of values given */
304 p=(const char *) kernel_string;
305 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
306 p++; /* ignore "'" chars for convolve filter usage - Cristy */
307 for (i=0; p < end; i++)
309 GetMagickToken(p,&p,token);
311 GetMagickToken(p,&p,token);
313 /* set the size of the kernel - old sized square */
314 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
315 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
316 p=(const char *) kernel_string;
317 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
318 p++; /* ignore "'" chars for convolve filter usage - Cristy */
321 /* Read in the kernel values from rest of input string argument */
322 kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
323 kernel->width,kernel->height*sizeof(*kernel->values)));
324 if (kernel->values == (MagickRealType *) NULL)
325 return(DestroyKernelInfo(kernel));
326 kernel->minimum=MagickMaximumValue;
327 kernel->maximum=(-MagickMaximumValue);
328 kernel->negative_range = kernel->positive_range = 0.0;
329 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
331 GetMagickToken(p,&p,token);
333 GetMagickToken(p,&p,token);
334 if ( LocaleCompare("nan",token) == 0
335 || LocaleCompare("-",token) == 0 ) {
336 kernel->values[i] = nan; /* this value is not part of neighbourhood */
339 kernel->values[i] = StringToDouble(token,(char **) NULL);
340 ( kernel->values[i] < 0)
341 ? ( kernel->negative_range += kernel->values[i] )
342 : ( kernel->positive_range += kernel->values[i] );
343 Minimize(kernel->minimum, kernel->values[i]);
344 Maximize(kernel->maximum, kernel->values[i]);
348 /* sanity check -- no more values in kernel definition */
349 GetMagickToken(p,&p,token);
350 if ( *token != '\0' && *token != ';' && *token != '\'' )
351 return(DestroyKernelInfo(kernel));
354 /* this was the old method of handling a incomplete kernel */
355 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
356 Minimize(kernel->minimum, kernel->values[i]);
357 Maximize(kernel->maximum, kernel->values[i]);
358 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
359 kernel->values[i]=0.0;
362 /* Number of values for kernel was not enough - Report Error */
363 if ( i < (ssize_t) (kernel->width*kernel->height) )
364 return(DestroyKernelInfo(kernel));
367 /* check that we recieved at least one real (non-nan) value! */
368 if (kernel->minimum == MagickMaximumValue)
369 return(DestroyKernelInfo(kernel));
371 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
372 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
373 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
374 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
375 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
376 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
381 static KernelInfo *ParseKernelName(const char *kernel_string,
382 ExceptionInfo *exception)
385 token[MaxTextExtent];
403 /* Parse special 'named' kernel */
404 GetMagickToken(kernel_string,&p,token);
405 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
406 if ( type < 0 || type == UserDefinedKernel )
407 return((KernelInfo *)NULL); /* not a valid named kernel */
409 while (((isspace((int) ((unsigned char) *p)) != 0) ||
410 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
413 end = strchr(p, ';'); /* end of this kernel defintion */
414 if ( end == (char *) NULL )
415 end = strchr(p, '\0');
417 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
418 memcpy(token, p, (size_t) (end-p));
420 SetGeometryInfo(&args);
421 flags = ParseGeometry(token, &args);
424 /* For Debugging Geometry Input */
425 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
426 flags, args.rho, args.sigma, args.xi, args.psi );
429 /* special handling of missing values in input string */
431 /* Shape Kernel Defaults */
433 if ( (flags & WidthValue) == 0 )
434 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
442 if ( (flags & HeightValue) == 0 )
443 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
446 if ( (flags & XValue) == 0 )
447 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
449 case RectangleKernel: /* Rectangle - set size defaults */
450 if ( (flags & WidthValue) == 0 ) /* if no width then */
451 args.rho = args.sigma; /* then width = height */
452 if ( args.rho < 1.0 ) /* if width too small */
453 args.rho = 3; /* then width = 3 */
454 if ( args.sigma < 1.0 ) /* if height too small */
455 args.sigma = args.rho; /* then height = width */
456 if ( (flags & XValue) == 0 ) /* center offset if not defined */
457 args.xi = (double)(((ssize_t)args.rho-1)/2);
458 if ( (flags & YValue) == 0 )
459 args.psi = (double)(((ssize_t)args.sigma-1)/2);
461 /* Distance Kernel Defaults */
462 case ChebyshevKernel:
463 case ManhattanKernel:
464 case OctagonalKernel:
465 case EuclideanKernel:
466 if ( (flags & HeightValue) == 0 ) /* no distance scale */
467 args.sigma = 100.0; /* default distance scaling */
468 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
469 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
470 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
471 args.sigma *= QuantumRange/100.0; /* percentage of color range */
477 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args, exception);
478 if ( kernel == (KernelInfo *) NULL )
481 /* global expand to rotated kernel list - only for single kernels */
482 if ( kernel->next == (KernelInfo *) NULL ) {
483 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
484 ExpandRotateKernelInfo(kernel, 45.0);
485 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
486 ExpandRotateKernelInfo(kernel, 90.0);
487 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
488 ExpandMirrorKernelInfo(kernel);
494 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string,
495 ExceptionInfo *exception)
503 token[MaxTextExtent];
508 if (kernel_string == (const char *) NULL)
509 return(ParseKernelArray(kernel_string));
511 kernel_cache=(char *) NULL;
512 if (*kernel_string == '@')
514 kernel_cache=FileToString(kernel_string+1,~0UL,exception);
515 if (kernel_cache == (char *) NULL)
516 return((KernelInfo *) NULL);
517 p=(const char *) kernel_cache;
520 while (GetMagickToken(p,NULL,token), *token != '\0')
522 /* ignore extra or multiple ';' kernel separators */
525 /* tokens starting with alpha is a Named kernel */
526 if (isalpha((int) ((unsigned char) *token)) != 0)
527 new_kernel=ParseKernelName(p,exception);
528 else /* otherwise a user defined kernel array */
529 new_kernel=ParseKernelArray(p);
531 /* Error handling -- this is not proper error handling! */
532 if (new_kernel == (KernelInfo *) NULL)
534 if (kernel != (KernelInfo *) NULL)
535 kernel=DestroyKernelInfo(kernel);
536 return((KernelInfo *) NULL);
539 /* initialise or append the kernel list */
540 if (kernel == (KernelInfo *) NULL)
543 LastKernelInfo(kernel)->next=new_kernel;
546 /* look for the next kernel in list */
548 if (p == (char *) NULL)
552 if (kernel_cache != (char *) NULL)
553 kernel_cache=DestroyString(kernel_cache);
558 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
562 % A c q u i r e K e r n e l B u i l t I n %
566 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
568 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
569 % kernels used for special purposes such as gaussian blurring, skeleton
570 % pruning, and edge distance determination.
572 % They take a KernelType, and a set of geometry style arguments, which were
573 % typically decoded from a user supplied string, or from a more complex
574 % Morphology Method that was requested.
576 % The format of the AcquireKernalBuiltIn method is:
578 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
579 % const GeometryInfo args)
581 % A description of each parameter follows:
583 % o type: the pre-defined type of kernel wanted
585 % o args: arguments defining or modifying the kernel
587 % Convolution Kernels
590 % The a No-Op or Scaling single element kernel.
592 % Gaussian:{radius},{sigma}
593 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
594 % The sigma for the curve is required. The resulting kernel is
597 % If 'sigma' is zero, you get a single pixel on a field of zeros.
599 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
600 % the final size of the resulting kernel to a square 2*radius+1 in size.
601 % The radius should be at least 2 times that of the sigma value, or
602 % sever clipping and aliasing may result. If not given or set to 0 the
603 % radius will be determined so as to produce the best minimal error
604 % result, which is usally much larger than is normally needed.
606 % LoG:{radius},{sigma}
607 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
608 % The supposed ideal edge detection, zero-summing kernel.
610 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
611 % approx 1.6 (according to wikipedia).
613 % DoG:{radius},{sigma1},{sigma2}
614 % "Difference of Gaussians" Kernel.
615 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
616 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
617 % The result is a zero-summing kernel.
619 % Blur:{radius},{sigma}[,{angle}]
620 % Generates a 1 dimensional or linear gaussian blur, at the angle given
621 % (current restricted to orthogonal angles). If a 'radius' is given the
622 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
623 % by a 90 degree angle.
625 % If 'sigma' is zero, you get a single pixel on a field of zeros.
627 % Note that two convolutions with two "Blur" kernels perpendicular to
628 % each other, is equivalent to a far larger "Gaussian" kernel with the
629 % same sigma value, However it is much faster to apply. This is how the
630 % "-blur" operator actually works.
632 % Comet:{width},{sigma},{angle}
633 % Blur in one direction only, much like how a bright object leaves
634 % a comet like trail. The Kernel is actually half a gaussian curve,
635 % Adding two such blurs in opposite directions produces a Blur Kernel.
636 % Angle can be rotated in multiples of 90 degrees.
638 % Note that the first argument is the width of the kernel and not the
639 % radius of the kernel.
641 % Binomial:[{radius}]
642 % Generate a discrete kernel using a 2 dimentional Pascel's Triangle
643 % of values. Used for special forma of image filters.
645 % # Still to be implemented...
649 % # Set kernel values using a resize filter, and given scale (sigma)
650 % # Cylindrical or Linear. Is this possible with an image?
653 % Named Constant Convolution Kernels
655 % All these are unscaled, zero-summing kernels by default. As such for
656 % non-HDRI version of ImageMagick some form of normalization, user scaling,
657 % and biasing the results is recommended, to prevent the resulting image
660 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
661 % 45 degrees to generate the 8 angled varients of each of the kernels.
664 % Discrete Lapacian Kernels, (without normalization)
665 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
666 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
667 % Type 2 : 3x3 with center:4 edge:1 corner:-2
668 % Type 3 : 3x3 with center:4 edge:-2 corner:1
669 % Type 5 : 5x5 laplacian
670 % Type 7 : 7x7 laplacian
671 % Type 15 : 5x5 LoG (sigma approx 1.4)
672 % Type 19 : 9x9 LoG (sigma approx 1.4)
675 % Sobel 'Edge' convolution kernel (3x3)
681 % Roberts convolution kernel (3x3)
687 % Prewitt Edge convolution kernel (3x3)
693 % Prewitt's "Compass" convolution kernel (3x3)
699 % Kirsch's "Compass" convolution kernel (3x3)
705 % Frei-Chen Edge Detector is based on a kernel that is similar to
706 % the Sobel Kernel, but is designed to be isotropic. That is it takes
707 % into account the distance of the diagonal in the kernel.
710 % | sqrt(2), 0, -sqrt(2) |
713 % FreiChen:{type},{angle}
715 % Frei-Chen Pre-weighted kernels...
717 % Type 0: default un-nomalized version shown above.
719 % Type 1: Orthogonal Kernel (same as type 11 below)
721 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
724 % Type 2: Diagonal form of Kernel...
726 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
729 % However this kernel is als at the heart of the FreiChen Edge Detection
730 % Process which uses a set of 9 specially weighted kernel. These 9
731 % kernels not be normalized, but directly applied to the image. The
732 % results is then added together, to produce the intensity of an edge in
733 % a specific direction. The square root of the pixel value can then be
734 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
735 % from each other, both the direction and the strength of the edge can be
738 % Type 10: All 9 of the following pre-weighted kernels...
740 % Type 11: | 1, 0, -1 |
741 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
744 % Type 12: | 1, sqrt(2), 1 |
745 % | 0, 0, 0 | / 2*sqrt(2)
748 % Type 13: | sqrt(2), -1, 0 |
749 % | -1, 0, 1 | / 2*sqrt(2)
752 % Type 14: | 0, 1, -sqrt(2) |
753 % | -1, 0, 1 | / 2*sqrt(2)
756 % Type 15: | 0, -1, 0 |
760 % Type 16: | 1, 0, -1 |
764 % Type 17: | 1, -2, 1 |
768 % Type 18: | -2, 1, -2 |
772 % Type 19: | 1, 1, 1 |
776 % The first 4 are for edge detection, the next 4 are for line detection
777 % and the last is to add a average component to the results.
779 % Using a special type of '-1' will return all 9 pre-weighted kernels
780 % as a multi-kernel list, so that you can use them directly (without
781 % normalization) with the special "-set option:morphology:compose Plus"
782 % setting to apply the full FreiChen Edge Detection Technique.
784 % If 'type' is large it will be taken to be an actual rotation angle for
785 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
786 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
788 % WARNING: The above was layed out as per
789 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
790 % But rotated 90 degrees so direction is from left rather than the top.
791 % I have yet to find any secondary confirmation of the above. The only
792 % other source found was actual source code at
793 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
794 % Neigher paper defineds the kernels in a way that looks locical or
795 % correct when taken as a whole.
799 % Diamond:[{radius}[,{scale}]]
800 % Generate a diamond shaped kernel with given radius to the points.
801 % Kernel size will again be radius*2+1 square and defaults to radius 1,
802 % generating a 3x3 kernel that is slightly larger than a square.
804 % Square:[{radius}[,{scale}]]
805 % Generate a square shaped kernel of size radius*2+1, and defaulting
806 % to a 3x3 (radius 1).
808 % Octagon:[{radius}[,{scale}]]
809 % Generate octagonal shaped kernel of given radius and constant scale.
810 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
811 % in "Diamond" kernel.
813 % Disk:[{radius}[,{scale}]]
814 % Generate a binary disk, thresholded at the radius given, the radius
815 % may be a float-point value. Final Kernel size is floor(radius)*2+1
816 % square. A radius of 5.3 is the default.
818 % NOTE: That a low radii Disk kernels produce the same results as
819 % many of the previously defined kernels, but differ greatly at larger
820 % radii. Here is a table of equivalences...
821 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
822 % "Disk:1.5" => "Square"
823 % "Disk:2" => "Diamond:2"
824 % "Disk:2.5" => "Octagon"
825 % "Disk:2.9" => "Square:2"
826 % "Disk:3.5" => "Octagon:3"
827 % "Disk:4.5" => "Octagon:4"
828 % "Disk:5.4" => "Octagon:5"
829 % "Disk:6.4" => "Octagon:6"
830 % All other Disk shapes are unique to this kernel, but because a "Disk"
831 % is more circular when using a larger radius, using a larger radius is
832 % preferred over iterating the morphological operation.
834 % Rectangle:{geometry}
835 % Simply generate a rectangle of 1's with the size given. You can also
836 % specify the location of the 'control point', otherwise the closest
837 % pixel to the center of the rectangle is selected.
839 % Properly centered and odd sized rectangles work the best.
841 % Symbol Dilation Kernels
843 % These kernel is not a good general morphological kernel, but is used
844 % more for highlighting and marking any single pixels in an image using,
845 % a "Dilate" method as appropriate.
847 % For the same reasons iterating these kernels does not produce the
848 % same result as using a larger radius for the symbol.
850 % Plus:[{radius}[,{scale}]]
851 % Cross:[{radius}[,{scale}]]
852 % Generate a kernel in the shape of a 'plus' or a 'cross' with
853 % a each arm the length of the given radius (default 2).
855 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
857 % Ring:{radius1},{radius2}[,{scale}]
858 % A ring of the values given that falls between the two radii.
859 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
860 % This is the 'edge' pixels of the default "Disk" kernel,
861 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
863 % Hit and Miss Kernels
865 % Peak:radius1,radius2
866 % Find any peak larger than the pixels the fall between the two radii.
867 % The default ring of pixels is as per "Ring".
869 % Find flat orthogonal edges of a binary shape
871 % Find 90 degree corners of a binary shape
873 % A special kernel to thin the 'outside' of diagonals
875 % Find end points of lines (for pruning a skeletion)
876 % Two types of lines ends (default to both) can be searched for
877 % Type 0: All line ends
878 % Type 1: single kernel for 4-conneected line ends
879 % Type 2: single kernel for simple line ends
881 % Find three line junctions (within a skeletion)
882 % Type 0: all line junctions
883 % Type 1: Y Junction kernel
884 % Type 2: Diagonal T Junction kernel
885 % Type 3: Orthogonal T Junction kernel
886 % Type 4: Diagonal X Junction kernel
887 % Type 5: Orthogonal + Junction kernel
889 % Find single pixel ridges or thin lines
890 % Type 1: Fine single pixel thick lines and ridges
891 % Type 2: Find two pixel thick lines and ridges
893 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
895 % Traditional skeleton generating kernels.
896 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
897 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
898 % Type 3: Thinning skeleton based on a ressearch paper by
899 % Dan S. Bloomberg (Default Type)
901 % A huge variety of Thinning Kernels designed to preserve conectivity.
902 % many other kernel sets use these kernels as source definitions.
903 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
904 % the super and sub notations used in the source research paper.
906 % Distance Measuring Kernels
908 % Different types of distance measuring methods, which are used with the
909 % a 'Distance' morphology method for generating a gradient based on
910 % distance from an edge of a binary shape, though there is a technique
911 % for handling a anti-aliased shape.
913 % See the 'Distance' Morphological Method, for information of how it is
916 % Chebyshev:[{radius}][x{scale}[%!]]
917 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
918 % is a value of one to any neighbour, orthogonal or diagonal. One why
919 % of thinking of it is the number of squares a 'King' or 'Queen' in
920 % chess needs to traverse reach any other position on a chess board.
921 % It results in a 'square' like distance function, but one where
922 % diagonals are given a value that is closer than expected.
924 % Manhattan:[{radius}][x{scale}[%!]]
925 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
926 % Cab distance metric), it is the distance needed when you can only
927 % travel in horizontal or vertical directions only. It is the
928 % distance a 'Rook' in chess would have to travel, and results in a
929 % diamond like distances, where diagonals are further than expected.
931 % Octagonal:[{radius}][x{scale}[%!]]
932 % An interleving of Manhatten and Chebyshev metrics producing an
933 % increasing octagonally shaped distance. Distances matches those of
934 % the "Octagon" shaped kernel of the same radius. The minimum radius
935 % and default is 2, producing a 5x5 kernel.
937 % Euclidean:[{radius}][x{scale}[%!]]
938 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
939 % However by default the kernel size only has a radius of 1, which
940 % limits the distance to 'Knight' like moves, with only orthogonal and
941 % diagonal measurements being correct. As such for the default kernel
942 % you will get octagonal like distance function.
944 % However using a larger radius such as "Euclidean:4" you will get a
945 % much smoother distance gradient from the edge of the shape. Especially
946 % if the image is pre-processed to include any anti-aliasing pixels.
947 % Of course a larger kernel is slower to use, and not always needed.
949 % The first three Distance Measuring Kernels will only generate distances
950 % of exact multiples of {scale} in binary images. As such you can use a
951 % scale of 1 without loosing any information. However you also need some
952 % scaling when handling non-binary anti-aliased shapes.
954 % The "Euclidean" Distance Kernel however does generate a non-integer
955 % fractional results, and as such scaling is vital even for binary shapes.
959 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
960 const GeometryInfo *args,ExceptionInfo *exception)
973 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
975 /* Generate a new empty kernel if needed */
976 kernel=(KernelInfo *) NULL;
978 case UndefinedKernel: /* These should not call this function */
979 case UserDefinedKernel:
980 assert("Should not call this function" != (char *)NULL);
982 case LaplacianKernel: /* Named Descrete Convolution Kernels */
983 case SobelKernel: /* these are defined using other kernels */
989 case EdgesKernel: /* Hit and Miss kernels */
991 case DiagonalsKernel:
993 case LineJunctionsKernel:
995 case ConvexHullKernel:
998 break; /* A pre-generated kernel is not needed */
1000 /* set to 1 to do a compile-time check that we haven't missed anything */
1002 case GaussianKernel:
1007 case BinomialKernel:
1010 case RectangleKernel:
1017 case ChebyshevKernel:
1018 case ManhattanKernel:
1019 case OctangonalKernel:
1020 case EuclideanKernel:
1024 /* Generate the base Kernel Structure */
1025 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1026 if (kernel == (KernelInfo *) NULL)
1028 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1029 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1030 kernel->negative_range = kernel->positive_range = 0.0;
1031 kernel->type = type;
1032 kernel->next = (KernelInfo *) NULL;
1033 kernel->signature = MagickSignature;
1043 kernel->height = kernel->width = (size_t) 1;
1044 kernel->x = kernel->y = (ssize_t) 0;
1045 kernel->values=(MagickRealType *) MagickAssumeAligned(
1046 AcquireAlignedMemory(1,sizeof(*kernel->values)));
1047 if (kernel->values == (MagickRealType *) NULL)
1048 return(DestroyKernelInfo(kernel));
1049 kernel->maximum = kernel->values[0] = args->rho;
1053 case GaussianKernel:
1057 sigma = fabs(args->sigma),
1058 sigma2 = fabs(args->xi),
1061 if ( args->rho >= 1.0 )
1062 kernel->width = (size_t)args->rho*2+1;
1063 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1064 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1066 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1067 kernel->height = kernel->width;
1068 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1069 kernel->values=(MagickRealType *) MagickAssumeAligned(
1070 AcquireAlignedMemory(kernel->width,kernel->height*
1071 sizeof(*kernel->values)));
1072 if (kernel->values == (MagickRealType *) NULL)
1073 return(DestroyKernelInfo(kernel));
1075 /* WARNING: The following generates a 'sampled gaussian' kernel.
1076 * What we really want is a 'discrete gaussian' kernel.
1078 * How to do this is I don't know, but appears to be basied on the
1079 * Error Function 'erf()' (intergral of a gaussian)
1082 if ( type == GaussianKernel || type == DoGKernel )
1083 { /* Calculate a Gaussian, OR positive half of a DoG */
1084 if ( sigma > MagickEpsilon )
1085 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1086 B = (double) (1.0/(Magick2PI*sigma*sigma));
1087 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1088 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1089 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1091 else /* limiting case - a unity (normalized Dirac) kernel */
1092 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1093 kernel->width*kernel->height*sizeof(*kernel->values));
1094 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1098 if ( type == DoGKernel )
1099 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1100 if ( sigma2 > MagickEpsilon )
1101 { sigma = sigma2; /* simplify loop expressions */
1102 A = 1.0/(2.0*sigma*sigma);
1103 B = (double) (1.0/(Magick2PI*sigma*sigma));
1104 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1105 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1106 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1108 else /* limiting case - a unity (normalized Dirac) kernel */
1109 kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1112 if ( type == LoGKernel )
1113 { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1114 if ( sigma > MagickEpsilon )
1115 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1116 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1117 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1118 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1119 { R = ((double)(u*u+v*v))*A;
1120 kernel->values[i] = (1-R)*exp(-R)*B;
1123 else /* special case - generate a unity kernel */
1124 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1125 kernel->width*kernel->height*sizeof(*kernel->values));
1126 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1130 /* Note the above kernels may have been 'clipped' by a user defined
1131 ** radius, producing a smaller (darker) kernel. Also for very small
1132 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1133 ** producing a very bright kernel.
1135 ** Normalization will still be needed.
1138 /* Normalize the 2D Gaussian Kernel
1140 ** NB: a CorrelateNormalize performs a normal Normalize if
1141 ** there are no negative values.
1143 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1144 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1150 sigma = fabs(args->sigma),
1153 if ( args->rho >= 1.0 )
1154 kernel->width = (size_t)args->rho*2+1;
1156 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1158 kernel->x = (ssize_t) (kernel->width-1)/2;
1160 kernel->negative_range = kernel->positive_range = 0.0;
1161 kernel->values=(MagickRealType *) MagickAssumeAligned(
1162 AcquireAlignedMemory(kernel->width,kernel->height*
1163 sizeof(*kernel->values)));
1164 if (kernel->values == (MagickRealType *) NULL)
1165 return(DestroyKernelInfo(kernel));
1168 #define KernelRank 3
1169 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1170 ** It generates a gaussian 3 times the width, and compresses it into
1171 ** the expected range. This produces a closer normalization of the
1172 ** resulting kernel, especially for very low sigma values.
1173 ** As such while wierd it is prefered.
1175 ** I am told this method originally came from Photoshop.
1177 ** A properly normalized curve is generated (apart from edge clipping)
1178 ** even though we later normalize the result (for edge clipping)
1179 ** to allow the correct generation of a "Difference of Blurs".
1183 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1184 (void) ResetMagickMemory(kernel->values,0, (size_t)
1185 kernel->width*kernel->height*sizeof(*kernel->values));
1186 /* Calculate a Positive 1D Gaussian */
1187 if ( sigma > MagickEpsilon )
1188 { sigma *= KernelRank; /* simplify loop expressions */
1189 alpha = 1.0/(2.0*sigma*sigma);
1190 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1191 for ( u=-v; u <= v; u++) {
1192 kernel->values[(u+v)/KernelRank] +=
1193 exp(-((double)(u*u))*alpha)*beta;
1196 else /* special case - generate a unity kernel */
1197 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1199 /* Direct calculation without curve averaging
1200 This is equivelent to a KernelRank of 1 */
1202 /* Calculate a Positive Gaussian */
1203 if ( sigma > MagickEpsilon )
1204 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1205 beta = 1.0/(MagickSQ2PI*sigma);
1206 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1207 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1209 else /* special case - generate a unity kernel */
1210 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1211 kernel->width*kernel->height*sizeof(*kernel->values));
1212 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1215 /* Note the above kernel may have been 'clipped' by a user defined
1216 ** radius, producing a smaller (darker) kernel. Also for very small
1217 ** sigma's (> 0.1) the central value becomes larger than one, as a
1218 ** result of not generating a actual 'discrete' kernel, and thus
1219 ** producing a very bright 'impulse'.
1221 ** Becuase of these two factors Normalization is required!
1224 /* Normalize the 1D Gaussian Kernel
1226 ** NB: a CorrelateNormalize performs a normal Normalize if
1227 ** there are no negative values.
1229 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1230 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1232 /* rotate the 1D kernel by given angle */
1233 RotateKernelInfo(kernel, args->xi );
1238 sigma = fabs(args->sigma),
1241 if ( args->rho < 1.0 )
1242 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1244 kernel->width = (size_t)args->rho;
1245 kernel->x = kernel->y = 0;
1247 kernel->negative_range = kernel->positive_range = 0.0;
1248 kernel->values=(MagickRealType *) MagickAssumeAligned(
1249 AcquireAlignedMemory(kernel->width,kernel->height*
1250 sizeof(*kernel->values)));
1251 if (kernel->values == (MagickRealType *) NULL)
1252 return(DestroyKernelInfo(kernel));
1254 /* A comet blur is half a 1D gaussian curve, so that the object is
1255 ** blurred in one direction only. This may not be quite the right
1256 ** curve to use so may change in the future. The function must be
1257 ** normalised after generation, which also resolves any clipping.
1259 ** As we are normalizing and not subtracting gaussians,
1260 ** there is no need for a divisor in the gaussian formula
1262 ** It is less comples
1264 if ( sigma > MagickEpsilon )
1267 #define KernelRank 3
1268 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1269 (void) ResetMagickMemory(kernel->values,0, (size_t)
1270 kernel->width*sizeof(*kernel->values));
1271 sigma *= KernelRank; /* simplify the loop expression */
1272 A = 1.0/(2.0*sigma*sigma);
1273 /* B = 1.0/(MagickSQ2PI*sigma); */
1274 for ( u=0; u < v; u++) {
1275 kernel->values[u/KernelRank] +=
1276 exp(-((double)(u*u))*A);
1277 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1279 for (i=0; i < (ssize_t) kernel->width; i++)
1280 kernel->positive_range += kernel->values[i];
1282 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1283 /* B = 1.0/(MagickSQ2PI*sigma); */
1284 for ( i=0; i < (ssize_t) kernel->width; i++)
1285 kernel->positive_range +=
1286 kernel->values[i] = exp(-((double)(i*i))*A);
1287 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1290 else /* special case - generate a unity kernel */
1291 { (void) ResetMagickMemory(kernel->values,0, (size_t)
1292 kernel->width*kernel->height*sizeof(*kernel->values));
1293 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1294 kernel->positive_range = 1.0;
1297 kernel->minimum = 0.0;
1298 kernel->maximum = kernel->values[0];
1299 kernel->negative_range = 0.0;
1301 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1302 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1305 case BinomialKernel:
1310 if (args->rho < 1.0)
1311 kernel->width = kernel->height = 3; /* default radius = 1 */
1313 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1314 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1316 order_f = fact(kernel->width-1);
1318 kernel->values=(MagickRealType *) MagickAssumeAligned(
1319 AcquireAlignedMemory(kernel->width,kernel->height*
1320 sizeof(*kernel->values)));
1321 if (kernel->values == (MagickRealType *) NULL)
1322 return(DestroyKernelInfo(kernel));
1324 /* set all kernel values within diamond area to scale given */
1325 for ( i=0, v=0; v < (ssize_t)kernel->height; v++)
1327 alpha = order_f / ( fact((size_t) v) * fact(kernel->height-v-1) );
1328 for ( u=0; u < (ssize_t)kernel->width; u++, i++)
1329 kernel->positive_range += kernel->values[i] = (double)
1330 (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-u-1) ));
1332 kernel->minimum = 1.0;
1333 kernel->maximum = kernel->values[kernel->x+kernel->y*kernel->width];
1334 kernel->negative_range = 0.0;
1339 Convolution Kernels - Well Known Named Constant Kernels
1341 case LaplacianKernel:
1342 { switch ( (int) args->rho ) {
1344 default: /* laplacian square filter -- default */
1345 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1347 case 1: /* laplacian diamond filter */
1348 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1351 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1354 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1356 case 5: /* a 5x5 laplacian */
1357 kernel=ParseKernelArray(
1358 "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");
1360 case 7: /* a 7x7 laplacian */
1361 kernel=ParseKernelArray(
1362 "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" );
1364 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1365 kernel=ParseKernelArray(
1366 "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");
1368 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1369 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1370 kernel=ParseKernelArray(
1371 "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");
1374 if (kernel == (KernelInfo *) NULL)
1376 kernel->type = type;
1380 { /* Simple Sobel Kernel */
1381 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1382 if (kernel == (KernelInfo *) NULL)
1384 kernel->type = type;
1385 RotateKernelInfo(kernel, args->rho);
1390 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1391 if (kernel == (KernelInfo *) NULL)
1393 kernel->type = type;
1394 RotateKernelInfo(kernel, args->rho);
1399 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1400 if (kernel == (KernelInfo *) NULL)
1402 kernel->type = type;
1403 RotateKernelInfo(kernel, args->rho);
1408 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1409 if (kernel == (KernelInfo *) NULL)
1411 kernel->type = type;
1412 RotateKernelInfo(kernel, args->rho);
1417 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1418 if (kernel == (KernelInfo *) NULL)
1420 kernel->type = type;
1421 RotateKernelInfo(kernel, args->rho);
1424 case FreiChenKernel:
1425 /* Direction is set to be left to right positive */
1426 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1427 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1428 { switch ( (int) args->rho ) {
1431 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1432 if (kernel == (KernelInfo *) NULL)
1434 kernel->type = type;
1435 kernel->values[3] = +(MagickRealType) MagickSQ2;
1436 kernel->values[5] = -(MagickRealType) MagickSQ2;
1437 CalcKernelMetaData(kernel); /* recalculate meta-data */
1440 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1441 if (kernel == (KernelInfo *) NULL)
1443 kernel->type = type;
1444 kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2;
1445 kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2;
1446 CalcKernelMetaData(kernel); /* recalculate meta-data */
1447 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1451 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19",exception);
1452 if (kernel == (KernelInfo *) NULL)
1458 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1459 if (kernel == (KernelInfo *) NULL)
1461 kernel->type = type;
1462 kernel->values[3] = +(MagickRealType) MagickSQ2;
1463 kernel->values[5] = -(MagickRealType) MagickSQ2;
1464 CalcKernelMetaData(kernel); /* recalculate meta-data */
1465 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1468 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1469 if (kernel == (KernelInfo *) NULL)
1471 kernel->type = type;
1472 kernel->values[1] = +(MagickRealType) MagickSQ2;
1473 kernel->values[7] = +(MagickRealType) MagickSQ2;
1474 CalcKernelMetaData(kernel);
1475 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1478 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1479 if (kernel == (KernelInfo *) NULL)
1481 kernel->type = type;
1482 kernel->values[0] = +(MagickRealType) MagickSQ2;
1483 kernel->values[8] = -(MagickRealType) MagickSQ2;
1484 CalcKernelMetaData(kernel);
1485 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1488 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1489 if (kernel == (KernelInfo *) NULL)
1491 kernel->type = type;
1492 kernel->values[2] = -(MagickRealType) MagickSQ2;
1493 kernel->values[6] = +(MagickRealType) MagickSQ2;
1494 CalcKernelMetaData(kernel);
1495 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1498 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1499 if (kernel == (KernelInfo *) NULL)
1501 kernel->type = type;
1502 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1505 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1506 if (kernel == (KernelInfo *) NULL)
1508 kernel->type = type;
1509 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1512 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1513 if (kernel == (KernelInfo *) NULL)
1515 kernel->type = type;
1516 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1519 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1520 if (kernel == (KernelInfo *) NULL)
1522 kernel->type = type;
1523 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1526 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1527 if (kernel == (KernelInfo *) NULL)
1529 kernel->type = type;
1530 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1533 if ( fabs(args->sigma) >= MagickEpsilon )
1534 /* Rotate by correctly supplied 'angle' */
1535 RotateKernelInfo(kernel, args->sigma);
1536 else if ( args->rho > 30.0 || args->rho < -30.0 )
1537 /* Rotate by out of bounds 'type' */
1538 RotateKernelInfo(kernel, args->rho);
1543 Boolean or Shaped Kernels
1547 if (args->rho < 1.0)
1548 kernel->width = kernel->height = 3; /* default radius = 1 */
1550 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1551 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1553 kernel->values=(MagickRealType *) MagickAssumeAligned(
1554 AcquireAlignedMemory(kernel->width,kernel->height*
1555 sizeof(*kernel->values)));
1556 if (kernel->values == (MagickRealType *) NULL)
1557 return(DestroyKernelInfo(kernel));
1559 /* set all kernel values within diamond area to scale given */
1560 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1561 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1562 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1563 kernel->positive_range += kernel->values[i] = args->sigma;
1565 kernel->values[i] = nan;
1566 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1570 case RectangleKernel:
1573 if ( type == SquareKernel )
1575 if (args->rho < 1.0)
1576 kernel->width = kernel->height = 3; /* default radius = 1 */
1578 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1579 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1580 scale = args->sigma;
1583 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1584 if ( args->rho < 1.0 || args->sigma < 1.0 )
1585 return(DestroyKernelInfo(kernel)); /* invalid args given */
1586 kernel->width = (size_t)args->rho;
1587 kernel->height = (size_t)args->sigma;
1588 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1589 args->psi < 0.0 || args->psi > (double)kernel->height )
1590 return(DestroyKernelInfo(kernel)); /* invalid args given */
1591 kernel->x = (ssize_t) args->xi;
1592 kernel->y = (ssize_t) args->psi;
1595 kernel->values=(MagickRealType *) MagickAssumeAligned(
1596 AcquireAlignedMemory(kernel->width,kernel->height*
1597 sizeof(*kernel->values)));
1598 if (kernel->values == (MagickRealType *) NULL)
1599 return(DestroyKernelInfo(kernel));
1601 /* set all kernel values to scale given */
1602 u=(ssize_t) (kernel->width*kernel->height);
1603 for ( i=0; i < u; i++)
1604 kernel->values[i] = scale;
1605 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1606 kernel->positive_range = scale*u;
1611 if (args->rho < 1.0)
1612 kernel->width = kernel->height = 5; /* default radius = 2 */
1614 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1615 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1617 kernel->values=(MagickRealType *) MagickAssumeAligned(
1618 AcquireAlignedMemory(kernel->width,kernel->height*
1619 sizeof(*kernel->values)));
1620 if (kernel->values == (MagickRealType *) NULL)
1621 return(DestroyKernelInfo(kernel));
1623 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1624 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1625 if ( (labs((long) u)+labs((long) v)) <=
1626 ((long)kernel->x + (long)(kernel->x/2)) )
1627 kernel->positive_range += kernel->values[i] = args->sigma;
1629 kernel->values[i] = nan;
1630 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1636 limit = (ssize_t)(args->rho*args->rho);
1638 if (args->rho < 0.4) /* default radius approx 4.3 */
1639 kernel->width = kernel->height = 9L, limit = 18L;
1641 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1642 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1644 kernel->values=(MagickRealType *) MagickAssumeAligned(
1645 AcquireAlignedMemory(kernel->width,kernel->height*
1646 sizeof(*kernel->values)));
1647 if (kernel->values == (MagickRealType *) NULL)
1648 return(DestroyKernelInfo(kernel));
1650 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1651 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1652 if ((u*u+v*v) <= limit)
1653 kernel->positive_range += kernel->values[i] = args->sigma;
1655 kernel->values[i] = nan;
1656 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1661 if (args->rho < 1.0)
1662 kernel->width = kernel->height = 5; /* default radius 2 */
1664 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1665 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1667 kernel->values=(MagickRealType *) MagickAssumeAligned(
1668 AcquireAlignedMemory(kernel->width,kernel->height*
1669 sizeof(*kernel->values)));
1670 if (kernel->values == (MagickRealType *) NULL)
1671 return(DestroyKernelInfo(kernel));
1673 /* set all kernel values along axises to given scale */
1674 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1675 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1676 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1677 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1678 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1683 if (args->rho < 1.0)
1684 kernel->width = kernel->height = 5; /* default radius 2 */
1686 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1687 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1689 kernel->values=(MagickRealType *) MagickAssumeAligned(
1690 AcquireAlignedMemory(kernel->width,kernel->height*
1691 sizeof(*kernel->values)));
1692 if (kernel->values == (MagickRealType *) NULL)
1693 return(DestroyKernelInfo(kernel));
1695 /* set all kernel values along axises to given scale */
1696 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1697 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1698 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1699 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1700 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1714 if (args->rho < args->sigma)
1716 kernel->width = ((size_t)args->sigma)*2+1;
1717 limit1 = (ssize_t)(args->rho*args->rho);
1718 limit2 = (ssize_t)(args->sigma*args->sigma);
1722 kernel->width = ((size_t)args->rho)*2+1;
1723 limit1 = (ssize_t)(args->sigma*args->sigma);
1724 limit2 = (ssize_t)(args->rho*args->rho);
1727 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1729 kernel->height = kernel->width;
1730 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1731 kernel->values=(MagickRealType *) MagickAssumeAligned(
1732 AcquireAlignedMemory(kernel->width,kernel->height*
1733 sizeof(*kernel->values)));
1734 if (kernel->values == (MagickRealType *) NULL)
1735 return(DestroyKernelInfo(kernel));
1737 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1738 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1739 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1740 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1741 { ssize_t radius=u*u+v*v;
1742 if (limit1 < radius && radius <= limit2)
1743 kernel->positive_range += kernel->values[i] = (double) scale;
1745 kernel->values[i] = nan;
1747 kernel->minimum = kernel->maximum = (double) scale;
1748 if ( type == PeaksKernel ) {
1749 /* set the central point in the middle */
1750 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1751 kernel->positive_range = 1.0;
1752 kernel->maximum = 1.0;
1758 kernel=AcquireKernelInfo("ThinSE:482",exception);
1759 if (kernel == (KernelInfo *) NULL)
1761 kernel->type = type;
1762 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1767 kernel=AcquireKernelInfo("ThinSE:87",exception);
1768 if (kernel == (KernelInfo *) NULL)
1770 kernel->type = type;
1771 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1774 case DiagonalsKernel:
1776 switch ( (int) args->rho ) {
1781 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1782 if (kernel == (KernelInfo *) NULL)
1784 kernel->type = type;
1785 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1786 if (new_kernel == (KernelInfo *) NULL)
1787 return(DestroyKernelInfo(kernel));
1788 new_kernel->type = type;
1789 LastKernelInfo(kernel)->next = new_kernel;
1790 ExpandMirrorKernelInfo(kernel);
1794 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1797 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1800 if (kernel == (KernelInfo *) NULL)
1802 kernel->type = type;
1803 RotateKernelInfo(kernel, args->sigma);
1806 case LineEndsKernel:
1807 { /* Kernels for finding the end of thin lines */
1808 switch ( (int) args->rho ) {
1811 /* set of kernels to find all end of lines */
1812 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>",exception));
1814 /* kernel for 4-connected line ends - no rotation */
1815 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1818 /* kernel to add for 8-connected lines - no rotation */
1819 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1822 /* kernel to add for orthogonal line ends - does not find corners */
1823 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1826 /* traditional line end - fails on last T end */
1827 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1830 if (kernel == (KernelInfo *) NULL)
1832 kernel->type = type;
1833 RotateKernelInfo(kernel, args->sigma);
1836 case LineJunctionsKernel:
1837 { /* kernels for finding the junctions of multiple lines */
1838 switch ( (int) args->rho ) {
1841 /* set of kernels to find all line junctions */
1842 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>",exception));
1845 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1848 /* Diagonal T Junctions */
1849 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1852 /* Orthogonal T Junctions */
1853 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1856 /* Diagonal X Junctions */
1857 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1860 /* Orthogonal X Junctions - minimal diamond kernel */
1861 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1864 if (kernel == (KernelInfo *) NULL)
1866 kernel->type = type;
1867 RotateKernelInfo(kernel, args->sigma);
1871 { /* Ridges - Ridge finding kernels */
1874 switch ( (int) args->rho ) {
1877 kernel=ParseKernelArray("3x1:0,1,0");
1878 if (kernel == (KernelInfo *) NULL)
1880 kernel->type = type;
1881 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1884 kernel=ParseKernelArray("4x1:0,1,1,0");
1885 if (kernel == (KernelInfo *) NULL)
1887 kernel->type = type;
1888 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1890 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1891 /* Unfortunatally we can not yet rotate a non-square kernel */
1892 /* But then we can't flip a non-symetrical kernel either */
1893 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1894 if (new_kernel == (KernelInfo *) NULL)
1895 return(DestroyKernelInfo(kernel));
1896 new_kernel->type = type;
1897 LastKernelInfo(kernel)->next = new_kernel;
1898 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1899 if (new_kernel == (KernelInfo *) NULL)
1900 return(DestroyKernelInfo(kernel));
1901 new_kernel->type = type;
1902 LastKernelInfo(kernel)->next = new_kernel;
1903 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1904 if (new_kernel == (KernelInfo *) NULL)
1905 return(DestroyKernelInfo(kernel));
1906 new_kernel->type = type;
1907 LastKernelInfo(kernel)->next = new_kernel;
1908 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1909 if (new_kernel == (KernelInfo *) NULL)
1910 return(DestroyKernelInfo(kernel));
1911 new_kernel->type = type;
1912 LastKernelInfo(kernel)->next = new_kernel;
1913 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1914 if (new_kernel == (KernelInfo *) NULL)
1915 return(DestroyKernelInfo(kernel));
1916 new_kernel->type = type;
1917 LastKernelInfo(kernel)->next = new_kernel;
1918 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1919 if (new_kernel == (KernelInfo *) NULL)
1920 return(DestroyKernelInfo(kernel));
1921 new_kernel->type = type;
1922 LastKernelInfo(kernel)->next = new_kernel;
1923 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1924 if (new_kernel == (KernelInfo *) NULL)
1925 return(DestroyKernelInfo(kernel));
1926 new_kernel->type = type;
1927 LastKernelInfo(kernel)->next = new_kernel;
1928 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1929 if (new_kernel == (KernelInfo *) NULL)
1930 return(DestroyKernelInfo(kernel));
1931 new_kernel->type = type;
1932 LastKernelInfo(kernel)->next = new_kernel;
1937 case ConvexHullKernel:
1941 /* first set of 8 kernels */
1942 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1943 if (kernel == (KernelInfo *) NULL)
1945 kernel->type = type;
1946 ExpandRotateKernelInfo(kernel, 90.0);
1947 /* append the mirror versions too - no flip function yet */
1948 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1949 if (new_kernel == (KernelInfo *) NULL)
1950 return(DestroyKernelInfo(kernel));
1951 new_kernel->type = type;
1952 ExpandRotateKernelInfo(new_kernel, 90.0);
1953 LastKernelInfo(kernel)->next = new_kernel;
1956 case SkeletonKernel:
1958 switch ( (int) args->rho ) {
1961 /* Traditional Skeleton...
1962 ** A cyclically rotated single kernel
1964 kernel=AcquireKernelInfo("ThinSE:482",exception);
1965 if (kernel == (KernelInfo *) NULL)
1967 kernel->type = type;
1968 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1971 /* HIPR Variation of the cyclic skeleton
1972 ** Corners of the traditional method made more forgiving,
1973 ** but the retain the same cyclic order.
1975 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;",exception);
1976 if (kernel == (KernelInfo *) NULL)
1978 if (kernel->next == (KernelInfo *) NULL)
1979 return(DestroyKernelInfo(kernel));
1980 kernel->type = type;
1981 kernel->next->type = type;
1982 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1985 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1986 ** "Connectivity-Preserving Morphological Image Thransformations"
1987 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1988 ** http://www.leptonica.com/papers/conn.pdf
1990 kernel=AcquireKernelInfo("ThinSE:41; ThinSE:42; ThinSE:43",
1992 if (kernel == (KernelInfo *) NULL)
1994 kernel->type = type;
1995 kernel->next->type = type;
1996 kernel->next->next->type = type;
1997 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
2003 { /* Special kernels for general thinning, while preserving connections
2004 ** "Connectivity-Preserving Morphological Image Thransformations"
2005 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
2006 ** http://www.leptonica.com/papers/conn.pdf
2008 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2010 ** Note kernels do not specify the origin pixel, allowing them
2011 ** to be used for both thickening and thinning operations.
2013 switch ( (int) args->rho ) {
2014 /* SE for 4-connected thinning */
2015 case 41: /* SE_4_1 */
2016 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
2018 case 42: /* SE_4_2 */
2019 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
2021 case 43: /* SE_4_3 */
2022 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
2024 case 44: /* SE_4_4 */
2025 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
2027 case 45: /* SE_4_5 */
2028 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
2030 case 46: /* SE_4_6 */
2031 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
2033 case 47: /* SE_4_7 */
2034 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
2036 case 48: /* SE_4_8 */
2037 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
2039 case 49: /* SE_4_9 */
2040 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
2042 /* SE for 8-connected thinning - negatives of the above */
2043 case 81: /* SE_8_0 */
2044 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
2046 case 82: /* SE_8_2 */
2047 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
2049 case 83: /* SE_8_3 */
2050 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
2052 case 84: /* SE_8_4 */
2053 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
2055 case 85: /* SE_8_5 */
2056 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2058 case 86: /* SE_8_6 */
2059 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2061 case 87: /* SE_8_7 */
2062 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2064 case 88: /* SE_8_8 */
2065 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2067 case 89: /* SE_8_9 */
2068 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2070 /* Special combined SE kernels */
2071 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2072 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2074 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2075 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2077 case 481: /* SE_48_1 - General Connected Corner Kernel */
2078 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2081 case 482: /* SE_48_2 - General Edge Kernel */
2082 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2085 if (kernel == (KernelInfo *) NULL)
2087 kernel->type = type;
2088 RotateKernelInfo(kernel, args->sigma);
2092 Distance Measuring Kernels
2094 case ChebyshevKernel:
2096 if (args->rho < 1.0)
2097 kernel->width = kernel->height = 3; /* default radius = 1 */
2099 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2100 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2102 kernel->values=(MagickRealType *) MagickAssumeAligned(
2103 AcquireAlignedMemory(kernel->width,kernel->height*
2104 sizeof(*kernel->values)));
2105 if (kernel->values == (MagickRealType *) NULL)
2106 return(DestroyKernelInfo(kernel));
2108 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2109 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2110 kernel->positive_range += ( kernel->values[i] =
2111 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2112 kernel->maximum = kernel->values[0];
2115 case ManhattanKernel:
2117 if (args->rho < 1.0)
2118 kernel->width = kernel->height = 3; /* default radius = 1 */
2120 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2121 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2123 kernel->values=(MagickRealType *) MagickAssumeAligned(
2124 AcquireAlignedMemory(kernel->width,kernel->height*
2125 sizeof(*kernel->values)));
2126 if (kernel->values == (MagickRealType *) NULL)
2127 return(DestroyKernelInfo(kernel));
2129 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2130 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2131 kernel->positive_range += ( kernel->values[i] =
2132 args->sigma*(labs((long) u)+labs((long) v)) );
2133 kernel->maximum = kernel->values[0];
2136 case OctagonalKernel:
2138 if (args->rho < 2.0)
2139 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2141 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2142 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2144 kernel->values=(MagickRealType *) MagickAssumeAligned(
2145 AcquireAlignedMemory(kernel->width,kernel->height*
2146 sizeof(*kernel->values)));
2147 if (kernel->values == (MagickRealType *) NULL)
2148 return(DestroyKernelInfo(kernel));
2150 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2151 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2154 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2155 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2156 kernel->positive_range += kernel->values[i] =
2157 args->sigma*MagickMax(r1,r2);
2159 kernel->maximum = kernel->values[0];
2162 case EuclideanKernel:
2164 if (args->rho < 1.0)
2165 kernel->width = kernel->height = 3; /* default radius = 1 */
2167 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2168 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2170 kernel->values=(MagickRealType *) MagickAssumeAligned(
2171 AcquireAlignedMemory(kernel->width,kernel->height*
2172 sizeof(*kernel->values)));
2173 if (kernel->values == (MagickRealType *) NULL)
2174 return(DestroyKernelInfo(kernel));
2176 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2177 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2178 kernel->positive_range += ( kernel->values[i] =
2179 args->sigma*sqrt((double)(u*u+v*v)) );
2180 kernel->maximum = kernel->values[0];
2185 /* No-Op Kernel - Basically just a single pixel on its own */
2186 kernel=ParseKernelArray("1:1");
2187 if (kernel == (KernelInfo *) NULL)
2189 kernel->type = UndefinedKernel;
2198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2202 % C l o n e K e r n e l I n f o %
2206 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2208 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2209 % can be modified without effecting the original. The cloned kernel should
2210 % be destroyed using DestoryKernelInfo() when no longer needed.
2212 % The format of the CloneKernelInfo method is:
2214 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2216 % A description of each parameter follows:
2218 % o kernel: the Morphology/Convolution kernel to be cloned
2221 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2229 assert(kernel != (KernelInfo *) NULL);
2230 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2231 if (new_kernel == (KernelInfo *) NULL)
2233 *new_kernel=(*kernel); /* copy values in structure */
2235 /* replace the values with a copy of the values */
2236 new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2237 AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2238 if (new_kernel->values == (MagickRealType *) NULL)
2239 return(DestroyKernelInfo(new_kernel));
2240 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2241 new_kernel->values[i]=kernel->values[i];
2243 /* Also clone the next kernel in the kernel list */
2244 if ( kernel->next != (KernelInfo *) NULL ) {
2245 new_kernel->next = CloneKernelInfo(kernel->next);
2246 if ( new_kernel->next == (KernelInfo *) NULL )
2247 return(DestroyKernelInfo(new_kernel));
2254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2258 % D e s t r o y K e r n e l I n f o %
2262 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2264 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2267 % The format of the DestroyKernelInfo method is:
2269 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2271 % A description of each parameter follows:
2273 % o kernel: the Morphology/Convolution kernel to be destroyed
2276 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2278 assert(kernel != (KernelInfo *) NULL);
2279 if (kernel->next != (KernelInfo *) NULL)
2280 kernel->next=DestroyKernelInfo(kernel->next);
2281 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2282 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2291 + E x p a n d M i r r o r K e r n e l I n f o %
2295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2297 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2298 % sequence of 90-degree rotated kernels but providing a reflected 180
2299 % rotatation, before the -/+ 90-degree rotations.
2301 % This special rotation order produces a better, more symetrical thinning of
2304 % The format of the ExpandMirrorKernelInfo method is:
2306 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2308 % A description of each parameter follows:
2310 % o kernel: the Morphology/Convolution kernel
2312 % This function is only internel to this module, as it is not finalized,
2313 % especially with regard to non-orthogonal angles, and rotation of larger
2318 static void FlopKernelInfo(KernelInfo *kernel)
2319 { /* Do a Flop by reversing each row. */
2327 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2328 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2329 t=k[x], k[x]=k[r], k[r]=t;
2331 kernel->x = kernel->width - kernel->x - 1;
2332 angle = fmod(angle+180.0, 360.0);
2336 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2344 clone = CloneKernelInfo(last);
2345 RotateKernelInfo(clone, 180); /* flip */
2346 LastKernelInfo(last)->next = clone;
2349 clone = CloneKernelInfo(last);
2350 RotateKernelInfo(clone, 90); /* transpose */
2351 LastKernelInfo(last)->next = clone;
2354 clone = CloneKernelInfo(last);
2355 RotateKernelInfo(clone, 180); /* flop */
2356 LastKernelInfo(last)->next = clone;
2362 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2366 + E x p a n d R o t a t e K e r n e l I n f o %
2370 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2372 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2373 % incrementally by the angle given, until the kernel repeats.
2375 % WARNING: 45 degree rotations only works for 3x3 kernels.
2376 % While 90 degree roatations only works for linear and square kernels
2378 % The format of the ExpandRotateKernelInfo method is:
2380 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2382 % A description of each parameter follows:
2384 % o kernel: the Morphology/Convolution kernel
2386 % o angle: angle to rotate in degrees
2388 % This function is only internel to this module, as it is not finalized,
2389 % especially with regard to non-orthogonal angles, and rotation of larger
2393 /* Internal Routine - Return true if two kernels are the same */
2394 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2395 const KernelInfo *kernel2)
2400 /* check size and origin location */
2401 if ( kernel1->width != kernel2->width
2402 || kernel1->height != kernel2->height
2403 || kernel1->x != kernel2->x
2404 || kernel1->y != kernel2->y )
2407 /* check actual kernel values */
2408 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2409 /* Test for Nan equivalence */
2410 if ( IfNaN(kernel1->values[i]) && !IfNaN(kernel2->values[i]) )
2412 if ( IfNaN(kernel2->values[i]) && !IfNaN(kernel1->values[i]) )
2414 /* Test actual values are equivalent */
2415 if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2422 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2429 DisableMSCWarning(4127)
2432 clone = CloneKernelInfo(last);
2433 RotateKernelInfo(clone, angle);
2434 if ( SameKernelInfo(kernel, clone) != MagickFalse )
2436 LastKernelInfo(last)->next = clone;
2439 clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
2444 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2448 + C a l c M e t a K e r n a l I n f o %
2452 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2454 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2455 % using the kernel values. This should only ne used if it is not possible to
2456 % calculate that meta-data in some easier way.
2458 % It is important that the meta-data is correct before ScaleKernelInfo() is
2459 % used to perform kernel normalization.
2461 % The format of the CalcKernelMetaData method is:
2463 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2465 % A description of each parameter follows:
2467 % o kernel: the Morphology/Convolution kernel to modify
2469 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2470 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2471 % however is not true for flat-shaped morphological kernels.
2473 % WARNING: Only the specific kernel pointed to is modified, not a list of
2476 % This is an internal function and not expected to be useful outside this
2477 % module. This could change however.
2479 static void CalcKernelMetaData(KernelInfo *kernel)
2484 kernel->minimum = kernel->maximum = 0.0;
2485 kernel->negative_range = kernel->positive_range = 0.0;
2486 for (i=0; i < (kernel->width*kernel->height); i++)
2488 if ( fabs(kernel->values[i]) < MagickEpsilon )
2489 kernel->values[i] = 0.0;
2490 ( kernel->values[i] < 0)
2491 ? ( kernel->negative_range += kernel->values[i] )
2492 : ( kernel->positive_range += kernel->values[i] );
2493 Minimize(kernel->minimum, kernel->values[i]);
2494 Maximize(kernel->maximum, kernel->values[i]);
2501 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2505 % M o r p h o l o g y A p p l y %
2509 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2511 % MorphologyApply() applies a morphological method, multiple times using
2512 % a list of multiple kernels. This is the method that should be called by
2513 % other 'operators' that internally use morphology operations as part of
2516 % It is basically equivalent to as MorphologyImage() (see below) but without
2517 % any user controls. This allows internel programs to use this method to
2518 % perform a specific task without possible interference by any API user
2519 % supplied settings.
2521 % It is MorphologyImage() task to extract any such user controls, and
2522 % pass them to this function for processing.
2524 % More specifically all given kernels should already be scaled, normalised,
2525 % and blended appropriatally before being parred to this routine. The
2526 % appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2528 % The format of the MorphologyApply method is:
2530 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2531 % const ssize_t iterations,const KernelInfo *kernel,
2532 % const CompositeMethod compose,const double bias,
2533 % ExceptionInfo *exception)
2535 % A description of each parameter follows:
2537 % o image: the source image
2539 % o method: the morphology method to be applied.
2541 % o iterations: apply the operation this many times (or no change).
2542 % A value of -1 means loop until no change found.
2543 % How this is applied may depend on the morphology method.
2544 % Typically this is a value of 1.
2546 % o channel: the channel type.
2548 % o kernel: An array of double representing the morphology kernel.
2550 % o compose: How to handle or merge multi-kernel results.
2551 % If 'UndefinedCompositeOp' use default for the Morphology method.
2552 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2553 % Otherwise merge the results using the compose method given.
2555 % o bias: Convolution Output Bias.
2557 % o exception: return any errors or warnings in this structure.
2560 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2561 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2562 ExceptionInfo *exception)
2564 #define MorphologyTag "Morphology/Image"
2590 assert(image != (Image *) NULL);
2591 assert(image->signature == MagickSignature);
2592 assert(morphology_image != (Image *) NULL);
2593 assert(morphology_image->signature == MagickSignature);
2594 assert(kernel != (KernelInfo *) NULL);
2595 assert(kernel->signature == MagickSignature);
2596 assert(exception != (ExceptionInfo *) NULL);
2597 assert(exception->signature == MagickSignature);
2600 image_view=AcquireVirtualCacheView(image,exception);
2601 morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2602 width=image->columns+kernel->width-1;
2607 case ConvolveMorphology:
2608 case DilateMorphology:
2609 case DilateIntensityMorphology:
2610 case IterativeDistanceMorphology:
2613 Kernel needs to used with reflection about origin.
2615 offset.x=(ssize_t) kernel->width-kernel->x-1;
2616 offset.y=(ssize_t) kernel->height-kernel->y-1;
2619 case ErodeMorphology:
2620 case ErodeIntensityMorphology:
2621 case HitAndMissMorphology:
2622 case ThinningMorphology:
2623 case ThickenMorphology:
2631 assert("Not a Primitive Morphology Method" != (char *) NULL);
2636 changes=(size_t *) AcquireQuantumMemory(GetOpenMPMaximumThreads(),
2638 if (changes == (size_t *) NULL)
2639 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2640 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2642 if ((method == ConvolveMorphology) && (kernel->width == 1))
2648 Special handling (for speed) of vertical (blur) kernels. This performs
2649 its handling in columns rather than in rows. This is only done
2650 for convolve as it is the only method that generates very large 1-D
2651 vertical kernels (such as a 'BlurKernel')
2653 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2654 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2655 magick_threads(image,morphology_image,image->columns,1)
2657 for (x=0; x < (ssize_t) image->columns; x++)
2660 id = GetOpenMPThreadId();
2662 register const Quantum
2674 if (status == MagickFalse)
2676 p=GetCacheViewVirtualPixels(image_view,x,-offset.y,1,image->rows+
2677 kernel->height-1,exception);
2678 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2679 morphology_image->rows,exception);
2680 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2685 center=(ssize_t) GetPixelChannels(image)*offset.y;
2686 for (y=0; y < (ssize_t) image->rows; y++)
2691 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2705 register const MagickRealType
2708 register const Quantum
2720 channel=GetPixelChannelChannel(image,i);
2721 traits=GetPixelChannelTraits(image,channel);
2722 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2723 if ((traits == UndefinedPixelTrait) ||
2724 (morphology_traits == UndefinedPixelTrait))
2726 if (((traits & CopyPixelTrait) != 0) ||
2727 (GetPixelReadMask(image,p+center) == 0))
2729 SetPixelChannel(morphology_image,channel,p[center+i],q);
2732 k=(&kernel->values[kernel->width*kernel->height-1]);
2737 if ((morphology_traits & BlendPixelTrait) == 0)
2738 for (v=0; v < (ssize_t) kernel->height; v++)
2740 for (u=0; u < (ssize_t) kernel->width; u++)
2742 if (IfNaN(*k) == MagickFalse)
2744 pixel+=(*k)*pixels[i];
2749 pixels+=GetPixelChannels(image);
2753 for (v=0; v < (ssize_t) kernel->height; v++)
2755 for (u=0; u < (ssize_t) kernel->width; u++)
2757 if (IfNaN(*k) == MagickFalse)
2759 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2760 pixel+=alpha*(*k)*pixels[i];
2765 pixels+=GetPixelChannels(image);
2768 if (fabs(pixel-p[center+i]) > MagickEpsilon)
2770 gamma=PerceptibleReciprocal(gamma);
2772 gamma*=(double) kernel->height*kernel->width/count;
2773 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2776 p+=GetPixelChannels(image);
2777 q+=GetPixelChannels(morphology_image);
2779 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2781 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2786 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2787 #pragma omp critical (MagickCore_MorphologyPrimitive)
2789 proceed=SetImageProgress(image,MorphologyTag,progress++,
2791 if (proceed == MagickFalse)
2795 morphology_image->type=image->type;
2796 morphology_view=DestroyCacheView(morphology_view);
2797 image_view=DestroyCacheView(image_view);
2798 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
2799 changed+=changes[i];
2800 changes=(size_t *) RelinquishMagickMemory(changes);
2801 return(status ? (ssize_t) changed : 0);
2804 Normal handling of horizontal or rectangular kernels (row by row).
2806 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2807 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2808 magick_threads(image,morphology_image,image->rows,1)
2810 for (y=0; y < (ssize_t) image->rows; y++)
2813 id = GetOpenMPThreadId();
2815 register const Quantum
2827 if (status == MagickFalse)
2829 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,
2830 kernel->height,exception);
2831 q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2833 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2838 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
2839 GetPixelChannels(image)*offset.x);
2840 for (x=0; x < (ssize_t) image->columns; x++)
2845 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2861 register const MagickRealType
2864 register const Quantum
2876 channel=GetPixelChannelChannel(image,i);
2877 traits=GetPixelChannelTraits(image,channel);
2878 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2879 if ((traits == UndefinedPixelTrait) ||
2880 (morphology_traits == UndefinedPixelTrait))
2882 if (((traits & CopyPixelTrait) != 0) ||
2883 (GetPixelReadMask(image,p+center) == 0))
2885 SetPixelChannel(morphology_image,channel,p[center+i],q);
2890 minimum=(double) QuantumRange;
2891 count=kernel->width*kernel->height;
2894 case ConvolveMorphology: pixel=bias; break;
2895 case HitAndMissMorphology: pixel=(double) QuantumRange; break;
2896 case ThinningMorphology: pixel=(double) QuantumRange; break;
2897 case ThickenMorphology: pixel=(double) QuantumRange; break;
2898 case ErodeMorphology: pixel=(double) QuantumRange; break;
2899 case DilateMorphology: pixel=0.0; break;
2900 case ErodeIntensityMorphology:
2901 case DilateIntensityMorphology:
2902 case IterativeDistanceMorphology:
2904 pixel=(double) p[center+i];
2907 default: pixel=0; break;
2912 case ConvolveMorphology:
2915 Weighted Average of pixels using reflected kernel
2917 For correct working of this operation for asymetrical kernels,
2918 the kernel needs to be applied in its reflected form. That is
2919 its values needs to be reversed.
2921 Correlation is actually the same as this but without reflecting
2922 the kernel, and thus 'lower-level' that Convolution. However as
2923 Convolution is the more common method used, and it does not
2924 really cost us much in terms of processing to use a reflected
2925 kernel, so it is Convolution that is implemented.
2927 Correlation will have its kernel reflected before calling this
2928 function to do a Convolve.
2930 For more details of Correlation vs Convolution see
2931 http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2933 k=(&kernel->values[kernel->width*kernel->height-1]);
2935 if ((morphology_traits & BlendPixelTrait) == 0)
2940 for (v=0; v < (ssize_t) kernel->height; v++)
2942 for (u=0; u < (ssize_t) kernel->width; u++)
2944 if (IfNaN(*k) == MagickFalse)
2946 pixel+=(*k)*pixels[i];
2950 pixels+=GetPixelChannels(image);
2952 pixels+=(image->columns-1)*GetPixelChannels(image);
2960 for (v=0; v < (ssize_t) kernel->height; v++)
2962 for (u=0; u < (ssize_t) kernel->width; u++)
2964 if (IfNaN(*k) == MagickFalse)
2966 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2967 pixel+=alpha*(*k)*pixels[i];
2972 pixels+=GetPixelChannels(image);
2974 pixels+=(image->columns-1)*GetPixelChannels(image);
2978 case ErodeMorphology:
2981 Minimum value within kernel neighbourhood.
2983 The kernel is not reflected for this operation. In normal
2984 Greyscale Morphology, the kernel value should be added
2985 to the real value, this is currently not done, due to the
2986 nature of the boolean kernels being used.
2989 for (v=0; v < (ssize_t) kernel->height; v++)
2991 for (u=0; u < (ssize_t) kernel->width; u++)
2993 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
2995 if ((double) pixels[i] < pixel)
2996 pixel=(double) pixels[i];
2999 pixels+=GetPixelChannels(image);
3001 pixels+=(image->columns-1)*GetPixelChannels(image);
3005 case DilateMorphology:
3008 Maximum value within kernel neighbourhood.
3010 For correct working of this operation for asymetrical kernels,
3011 the kernel needs to be applied in its reflected form. That is
3012 its values needs to be reversed.
3014 In normal Greyscale Morphology, the kernel value should be
3015 added to the real value, this is currently not done, due to the
3016 nature of the boolean kernels being used.
3019 k=(&kernel->values[kernel->width*kernel->height-1]);
3020 for (v=0; v < (ssize_t) kernel->height; v++)
3022 for (u=0; u < (ssize_t) kernel->width; u++)
3024 if ((IfNaN(*k) == MagickFalse) && (*k > 0.5))
3026 if ((double) pixels[i] > pixel)
3027 pixel=(double) pixels[i];
3030 pixels+=GetPixelChannels(image);
3032 pixels+=(image->columns-1)*GetPixelChannels(image);
3036 case HitAndMissMorphology:
3037 case ThinningMorphology:
3038 case ThickenMorphology:
3041 Minimum of foreground pixel minus maxumum of background pixels.
3043 The kernel is not reflected for this operation, and consists
3044 of both foreground and background pixel neighbourhoods, 0.0 for
3045 background, and 1.0 for foreground with either Nan or 0.5 values
3048 This never produces a meaningless negative result. Such results
3049 cause Thinning/Thicken to not work correctly when used against a
3054 for (v=0; v < (ssize_t) kernel->height; v++)
3056 for (u=0; u < (ssize_t) kernel->width; u++)
3058 if (IfNaN(*k) == MagickFalse)
3062 if ((double) pixels[i] < pixel)
3063 pixel=(double) pixels[i];
3068 if ((double) pixels[i] > maximum)
3069 maximum=(double) pixels[i];
3074 pixels+=GetPixelChannels(image);
3076 pixels+=(image->columns-1)*GetPixelChannels(image);
3081 if (method == ThinningMorphology)
3082 pixel=(double) p[center+i]-pixel;
3084 if (method == ThickenMorphology)
3085 pixel+=(double) p[center+i]+pixel;
3088 case ErodeIntensityMorphology:
3091 Select pixel with minimum intensity within kernel neighbourhood.
3093 The kernel is not reflected for this operation.
3097 for (v=0; v < (ssize_t) kernel->height; v++)
3099 for (u=0; u < (ssize_t) kernel->width; u++)
3101 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3103 if (GetPixelIntensity(image,pixels) < minimum)
3105 pixel=(double) pixels[i];
3106 minimum=GetPixelIntensity(image,pixels);
3111 pixels+=GetPixelChannels(image);
3113 pixels+=(image->columns-1)*GetPixelChannels(image);
3117 case DilateIntensityMorphology:
3120 Select pixel with maximum intensity within kernel neighbourhood.
3122 The kernel is not reflected for this operation.
3125 k=(&kernel->values[kernel->width*kernel->height-1]);
3126 for (v=0; v < (ssize_t) kernel->height; v++)
3128 for (u=0; u < (ssize_t) kernel->width; u++)
3130 if ((IfNaN(*k) == MagickFalse) && (*k >= 0.5))
3132 if (GetPixelIntensity(image,pixels) > maximum)
3134 pixel=(double) pixels[i];
3135 maximum=GetPixelIntensity(image,pixels);
3140 pixels+=GetPixelChannels(image);
3142 pixels+=(image->columns-1)*GetPixelChannels(image);
3146 case IterativeDistanceMorphology:
3149 Compute th iterative distance from black edge of a white image
3150 shape. Essentually white values are decreased to the smallest
3151 'distance from edge' it can find.
3153 It works by adding kernel values to the neighbourhood, and and
3154 select the minimum value found. The kernel is rotated before
3155 use, so kernel distances match resulting distances, when a user
3156 provided asymmetric kernel is applied.
3158 This code is nearly identical to True GrayScale Morphology but
3161 GreyDilate Kernel values added, maximum value found Kernel is
3164 GrayErode: Kernel values subtracted and minimum value found No
3165 kernel rotation used.
3167 Note the the Iterative Distance method is essentially a
3168 GrayErode, but with negative kernel values, and kernel rotation
3172 k=(&kernel->values[kernel->width*kernel->height-1]);
3173 for (v=0; v < (ssize_t) kernel->height; v++)
3175 for (u=0; u < (ssize_t) kernel->width; u++)
3177 if (IfNaN(*k) == MagickFalse)
3179 if ((pixels[i]+(*k)) < pixel)
3180 pixel=(double) pixels[i]+(*k);
3184 pixels+=GetPixelChannels(image);
3186 pixels+=(image->columns-1)*GetPixelChannels(image);
3190 case UndefinedMorphology:
3194 if (fabs(pixel-p[center+i]) > MagickEpsilon)
3196 gamma=PerceptibleReciprocal(gamma);
3198 gamma*=(double) kernel->height*kernel->width/count;
3199 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*pixel),q);
3201 p+=GetPixelChannels(image);
3202 q+=GetPixelChannels(morphology_image);
3204 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3206 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3211 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3212 #pragma omp critical (MagickCore_MorphologyPrimitive)
3214 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3215 if (proceed == MagickFalse)
3219 morphology_view=DestroyCacheView(morphology_view);
3220 image_view=DestroyCacheView(image_view);
3221 for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
3222 changed+=changes[i];
3223 changes=(size_t *) RelinquishMagickMemory(changes);
3224 return(status ? (ssize_t) changed : -1);
3228 This is almost identical to the MorphologyPrimative() function above, but
3229 applies the primitive directly to the actual image using two passes, once in
3230 each direction, with the results of the previous (and current) row being
3233 That is after each row is 'Sync'ed' into the image, the next row makes use of
3234 those values as part of the calculation of the next row. It repeats, but
3235 going in the oppisite (bottom-up) direction.
3237 Because of this 're-use of results' this function can not make use of multi-
3238 threaded, parellel processing.
3240 static ssize_t MorphologyPrimitiveDirect(Image *image,
3241 const MorphologyMethod method,const KernelInfo *kernel,
3242 ExceptionInfo *exception)
3264 assert(image != (Image *) NULL);
3265 assert(image->signature == MagickSignature);
3266 assert(kernel != (KernelInfo *) NULL);
3267 assert(kernel->signature == MagickSignature);
3268 assert(exception != (ExceptionInfo *) NULL);
3269 assert(exception->signature == MagickSignature);
3275 case DistanceMorphology:
3276 case VoronoiMorphology:
3279 Kernel reflected about origin.
3281 offset.x=(ssize_t) kernel->width-kernel->x-1;
3282 offset.y=(ssize_t) kernel->height-kernel->y-1;
3293 Two views into same image, do not thread.
3295 image_view=AcquireVirtualCacheView(image,exception);
3296 morphology_view=AcquireAuthenticCacheView(image,exception);
3297 width=image->columns+kernel->width-1;
3298 for (y=0; y < (ssize_t) image->rows; y++)
3300 register const Quantum
3313 Read virtual pixels, and authentic pixels, from the same image! We read
3314 using virtual to get virtual pixel handling, but write back into the same
3317 Only top half of kernel is processed as we do a single pass downward
3318 through the image iterating the distance function as we go.
3320 if (status == MagickFalse)
3322 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t)
3323 offset.y+1,exception);
3324 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3326 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3331 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
3332 GetPixelChannels(image)*offset.x);
3333 for (x=0; x < (ssize_t) image->columns; x++)
3338 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3346 register const MagickRealType
3349 register const Quantum
3358 traits=GetPixelChannelTraits(image,(PixelChannel) i);
3359 if (traits == UndefinedPixelTrait)
3361 if (((traits & CopyPixelTrait) != 0) ||
3362 (GetPixelReadMask(image,p+center) == 0))
3365 pixel=(double) QuantumRange;
3368 case DistanceMorphology:
3370 k=(&kernel->values[kernel->width*kernel->height-1]);
3371 for (v=0; v <= offset.y; v++)
3373 for (u=0; u < (ssize_t) kernel->width; u++)
3375 if (IfNaN(*k) == MagickFalse)
3377 if ((pixels[i]+(*k)) < pixel)
3378 pixel=(double) pixels[i]+(*k);
3381 pixels+=GetPixelChannels(image);
3383 pixels+=(image->columns-1)*GetPixelChannels(image);
3385 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3386 pixels=q-offset.x*GetPixelChannels(image);
3387 for (u=0; u < offset.x; u++)
3389 if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3391 if ((pixels[i]+(*k)) < pixel)
3392 pixel=(double) pixels[i]+(*k);
3395 pixels+=GetPixelChannels(image);
3399 case VoronoiMorphology:
3401 k=(&kernel->values[kernel->width*kernel->height-1]);
3402 for (v=0; v < offset.y; v++)
3404 for (u=0; u < (ssize_t) kernel->width; u++)
3406 if (IfNaN(*k) == MagickFalse)
3408 if ((pixels[i]+(*k)) < pixel)
3409 pixel=(double) pixels[i]+(*k);
3412 pixels+=GetPixelChannels(image);
3414 pixels+=(image->columns-1)*GetPixelChannels(image);
3416 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3417 pixels=q-offset.x*GetPixelChannels(image);
3418 for (u=0; u < offset.x; u++)
3420 if ((IfNaN(*k) == MagickFalse) && ((x+u-offset.x) >= 0))
3422 if ((pixels[i]+(*k)) < pixel)
3423 pixel=(double) pixels[i]+(*k);
3426 pixels+=GetPixelChannels(image);
3433 if (fabs(pixel-q[i]) > MagickEpsilon)
3435 q[i]=ClampToQuantum(pixel);
3437 p+=GetPixelChannels(image);
3438 q+=GetPixelChannels(image);
3440 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3442 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3447 proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3448 if (proceed == MagickFalse)
3452 morphology_view=DestroyCacheView(morphology_view);
3453 image_view=DestroyCacheView(image_view);
3455 Do the reverse pass through the image.
3457 image_view=AcquireVirtualCacheView(image,exception);
3458 morphology_view=AcquireAuthenticCacheView(image,exception);
3459 for (y=(ssize_t) image->rows-1; y >= 0; y--)
3461 register const Quantum
3474 Read virtual pixels, and authentic pixels, from the same image. We
3475 read using virtual to get virtual pixel handling, but write back
3476 into the same image.
3478 Only the bottom half of the kernel is processed as we up the image.
3480 if (status == MagickFalse)
3482 p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t)
3483 kernel->y+1,exception);
3484 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3486 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3491 p+=(image->columns-1)*GetPixelChannels(image);
3492 q+=(image->columns-1)*GetPixelChannels(image);
3493 center=(ssize_t) (offset.x*GetPixelChannels(image));
3494 for (x=(ssize_t) image->columns-1; x >= 0; x--)
3499 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3507 register const MagickRealType
3510 register const Quantum
3519 traits=GetPixelChannelTraits(image,(PixelChannel) i);
3520 if (traits == UndefinedPixelTrait)
3522 if (((traits & CopyPixelTrait) != 0) ||
3523 (GetPixelReadMask(image,p+center) == 0))
3526 pixel=(double) QuantumRange;
3529 case DistanceMorphology:
3531 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3532 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3534 for (u=0; u < (ssize_t) kernel->width; u++)
3536 if (IfNaN(*k) == MagickFalse)
3538 if ((pixels[i]+(*k)) < pixel)
3539 pixel=(double) pixels[i]+(*k);
3542 pixels+=GetPixelChannels(image);
3544 pixels+=(image->columns-1)*GetPixelChannels(image);
3546 k=(&kernel->values[kernel->width*kernel->y+kernel->x-1]);
3548 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3550 pixels+=GetPixelChannels(image);
3551 if ((IfNaN(*k) == MagickFalse) &&
3552 ((x+u-offset.x) < (ssize_t) image->columns))
3554 if ((pixels[i]+(*k)) < pixel)
3555 pixel=(double) pixels[i]+(*k);
3561 case VoronoiMorphology:
3563 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3564 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3566 for (u=0; u < (ssize_t) kernel->width; u++)
3568 if (IfNaN(*k) == MagickFalse)
3570 if ((pixels[i]+(*k)) < pixel)
3571 pixel=(double) pixels[i]+(*k);
3574 pixels+=GetPixelChannels(image);
3576 pixels+=(image->columns-1)*GetPixelChannels(image);
3578 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3580 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3582 pixels+=GetPixelChannels(image);
3583 if ((IfNaN(*k) == MagickFalse) &&
3584 ((x+u-offset.x) < (ssize_t) image->columns))
3586 if ((pixels[i]+(*k)) < pixel)
3587 pixel=(double) pixels[i]+(*k);
3596 if (fabs(pixel-q[i]) > MagickEpsilon)
3598 q[i]=ClampToQuantum(pixel);
3600 p-=GetPixelChannels(image);
3601 q-=GetPixelChannels(image);
3603 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3605 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3610 proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3611 if (proceed == MagickFalse)
3615 morphology_view=DestroyCacheView(morphology_view);
3616 image_view=DestroyCacheView(image_view);
3617 return(status ? (ssize_t) changed : -1);
3621 Apply a Morphology by calling one of the above low level primitive
3622 application functions. This function handles any iteration loops,
3623 composition or re-iteration of results, and compound morphology methods that
3624 is based on multiple low-level (staged) morphology methods.
3626 Basically this provides the complex glue between the requested morphology
3627 method and raw low-level implementation (above).
3629 MagickPrivate Image *MorphologyApply(const Image *image,
3630 const MorphologyMethod method, const ssize_t iterations,
3631 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3632 ExceptionInfo *exception)
3638 *curr_image, /* Image we are working with or iterating */
3639 *work_image, /* secondary image for primitive iteration */
3640 *save_image, /* saved image - for 'edge' method only */
3641 *rslt_image; /* resultant image - after multi-kernel handling */
3644 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3645 *norm_kernel, /* the current normal un-reflected kernel */
3646 *rflt_kernel, /* the current reflected kernel (if needed) */
3647 *this_kernel; /* the kernel being applied */
3650 primitive; /* the current morphology primitive being applied */
3653 rslt_compose; /* multi-kernel compose method for results to use */
3656 special, /* do we use a direct modify function? */
3657 verbose; /* verbose output of results */
3660 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3661 method_limit, /* maximum number of compound method iterations */
3662 kernel_number, /* Loop 2: the kernel number being applied */
3663 stage_loop, /* Loop 3: primitive loop for compound morphology */
3664 stage_limit, /* how many primitives are in this compound */
3665 kernel_loop, /* Loop 4: iterate the kernel over image */
3666 kernel_limit, /* number of times to iterate kernel */
3667 count, /* total count of primitive steps applied */
3668 kernel_changed, /* total count of changed using iterated kernel */
3669 method_changed; /* total count of changed over method iteration */
3672 changed; /* number pixels changed by last primitive operation */
3675 v_info[MaxTextExtent];
3677 assert(image != (Image *) NULL);
3678 assert(image->signature == MagickSignature);
3679 assert(kernel != (KernelInfo *) NULL);
3680 assert(kernel->signature == MagickSignature);
3681 assert(exception != (ExceptionInfo *) NULL);
3682 assert(exception->signature == MagickSignature);
3684 count = 0; /* number of low-level morphology primitives performed */
3685 if ( iterations == 0 )
3686 return((Image *)NULL); /* null operation - nothing to do! */
3688 kernel_limit = (size_t) iterations;
3689 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3690 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3692 verbose = IsStringTrue(GetImageArtifact(image,"debug"));
3694 /* initialise for cleanup */
3695 curr_image = (Image *) image;
3696 curr_compose = image->compose;
3697 (void) curr_compose;
3698 work_image = save_image = rslt_image = (Image *) NULL;
3699 reflected_kernel = (KernelInfo *) NULL;
3701 /* Initialize specific methods
3702 * + which loop should use the given iteratations
3703 * + how many primitives make up the compound morphology
3704 * + multi-kernel compose method to use (by default)
3706 method_limit = 1; /* just do method once, unless otherwise set */
3707 stage_limit = 1; /* assume method is not a compound */
3708 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3709 rslt_compose = compose; /* and we are composing multi-kernels as given */
3711 case SmoothMorphology: /* 4 primitive compound morphology */
3714 case OpenMorphology: /* 2 primitive compound morphology */
3715 case OpenIntensityMorphology:
3716 case TopHatMorphology:
3717 case CloseMorphology:
3718 case CloseIntensityMorphology:
3719 case BottomHatMorphology:
3720 case EdgeMorphology:
3723 case HitAndMissMorphology:
3724 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3726 case ThinningMorphology:
3727 case ThickenMorphology:
3728 method_limit = kernel_limit; /* iterate the whole method */
3729 kernel_limit = 1; /* do not do kernel iteration */
3731 case DistanceMorphology:
3732 case VoronoiMorphology:
3733 special = MagickTrue; /* use special direct primative */
3739 /* Apply special methods with special requirments
3740 ** For example, single run only, or post-processing requirements
3742 if ( special != MagickFalse )
3744 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3745 if (rslt_image == (Image *) NULL)
3747 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3750 changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception);
3752 if ( IfMagickTrue(verbose) )
3753 (void) (void) FormatLocaleFile(stderr,
3754 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3755 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3756 1.0,0.0,1.0, (double) changed);
3761 if ( method == VoronoiMorphology ) {
3762 /* Preserve the alpha channel of input image - but turned it off */
3763 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3765 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3766 MagickTrue,0,0,exception);
3767 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3773 /* Handle user (caller) specified multi-kernel composition method */
3774 if ( compose != UndefinedCompositeOp )
3775 rslt_compose = compose; /* override default composition for method */
3776 if ( rslt_compose == UndefinedCompositeOp )
3777 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3779 /* Some methods require a reflected kernel to use with primitives.
3780 * Create the reflected kernel for those methods. */
3782 case CorrelateMorphology:
3783 case CloseMorphology:
3784 case CloseIntensityMorphology:
3785 case BottomHatMorphology:
3786 case SmoothMorphology:
3787 reflected_kernel = CloneKernelInfo(kernel);
3788 if (reflected_kernel == (KernelInfo *) NULL)
3790 RotateKernelInfo(reflected_kernel,180);
3796 /* Loops around more primitive morpholgy methods
3797 ** erose, dilate, open, close, smooth, edge, etc...
3799 /* Loop 1: iterate the compound method */
3802 while ( method_loop < method_limit && method_changed > 0 ) {
3806 /* Loop 2: iterate over each kernel in a multi-kernel list */
3807 norm_kernel = (KernelInfo *) kernel;
3808 this_kernel = (KernelInfo *) kernel;
3809 rflt_kernel = reflected_kernel;
3812 while ( norm_kernel != NULL ) {
3814 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3815 stage_loop = 0; /* the compound morphology stage number */
3816 while ( stage_loop < stage_limit ) {
3817 stage_loop++; /* The stage of the compound morphology */
3819 /* Select primitive morphology for this stage of compound method */
3820 this_kernel = norm_kernel; /* default use unreflected kernel */
3821 primitive = method; /* Assume method is a primitive */
3823 case ErodeMorphology: /* just erode */
3824 case EdgeInMorphology: /* erode and image difference */
3825 primitive = ErodeMorphology;
3827 case DilateMorphology: /* just dilate */
3828 case EdgeOutMorphology: /* dilate and image difference */
3829 primitive = DilateMorphology;
3831 case OpenMorphology: /* erode then dialate */
3832 case TopHatMorphology: /* open and image difference */
3833 primitive = ErodeMorphology;
3834 if ( stage_loop == 2 )
3835 primitive = DilateMorphology;
3837 case OpenIntensityMorphology:
3838 primitive = ErodeIntensityMorphology;
3839 if ( stage_loop == 2 )
3840 primitive = DilateIntensityMorphology;
3842 case CloseMorphology: /* dilate, then erode */
3843 case BottomHatMorphology: /* close and image difference */
3844 this_kernel = rflt_kernel; /* use the reflected kernel */
3845 primitive = DilateMorphology;
3846 if ( stage_loop == 2 )
3847 primitive = ErodeMorphology;
3849 case CloseIntensityMorphology:
3850 this_kernel = rflt_kernel; /* use the reflected kernel */
3851 primitive = DilateIntensityMorphology;
3852 if ( stage_loop == 2 )
3853 primitive = ErodeIntensityMorphology;
3855 case SmoothMorphology: /* open, close */
3856 switch ( stage_loop ) {
3857 case 1: /* start an open method, which starts with Erode */
3858 primitive = ErodeMorphology;
3860 case 2: /* now Dilate the Erode */
3861 primitive = DilateMorphology;
3863 case 3: /* Reflect kernel a close */
3864 this_kernel = rflt_kernel; /* use the reflected kernel */
3865 primitive = DilateMorphology;
3867 case 4: /* Finish the Close */
3868 this_kernel = rflt_kernel; /* use the reflected kernel */
3869 primitive = ErodeMorphology;
3873 case EdgeMorphology: /* dilate and erode difference */
3874 primitive = DilateMorphology;
3875 if ( stage_loop == 2 ) {
3876 save_image = curr_image; /* save the image difference */
3877 curr_image = (Image *) image;
3878 primitive = ErodeMorphology;
3881 case CorrelateMorphology:
3882 /* A Correlation is a Convolution with a reflected kernel.
3883 ** However a Convolution is a weighted sum using a reflected
3884 ** kernel. It may seem stange to convert a Correlation into a
3885 ** Convolution as the Correlation is the simplier method, but
3886 ** Convolution is much more commonly used, and it makes sense to
3887 ** implement it directly so as to avoid the need to duplicate the
3888 ** kernel when it is not required (which is typically the
3891 this_kernel = rflt_kernel; /* use the reflected kernel */
3892 primitive = ConvolveMorphology;
3897 assert( this_kernel != (KernelInfo *) NULL );
3899 /* Extra information for debugging compound operations */
3900 if ( IfMagickTrue(verbose) ) {
3901 if ( stage_limit > 1 )
3902 (void) FormatLocaleString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
3903 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
3904 method_loop,(double) stage_loop);
3905 else if ( primitive != method )
3906 (void) FormatLocaleString(v_info, MaxTextExtent, "%s:%.20g -> ",
3907 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
3913 /* Loop 4: Iterate the kernel with primitive */
3917 while ( kernel_loop < kernel_limit && changed > 0 ) {
3918 kernel_loop++; /* the iteration of this kernel */
3920 /* Create a clone as the destination image, if not yet defined */
3921 if ( work_image == (Image *) NULL )
3923 work_image=CloneImage(image,0,0,MagickTrue,exception);
3924 if (work_image == (Image *) NULL)
3926 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
3930 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3932 changed = MorphologyPrimitive(curr_image, work_image, primitive,
3933 this_kernel, bias, exception);
3934 if ( IfMagickTrue(verbose) ) {
3935 if ( kernel_loop > 1 )
3936 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
3937 (void) (void) FormatLocaleFile(stderr,
3938 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3939 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
3940 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
3941 (double) (method_loop+kernel_loop-1),(double) kernel_number,
3942 (double) count,(double) changed);
3946 kernel_changed += changed;
3947 method_changed += changed;
3949 /* prepare next loop */
3950 { Image *tmp = work_image; /* swap images for iteration */
3951 work_image = curr_image;
3954 if ( work_image == image )
3955 work_image = (Image *) NULL; /* replace input 'image' */
3957 } /* End Loop 4: Iterate the kernel with primitive */
3959 if ( IfMagickTrue(verbose) && kernel_changed != (size_t)changed )
3960 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
3961 if ( IfMagickTrue(verbose) && stage_loop < stage_limit )
3962 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
3965 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3966 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
3967 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
3968 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
3969 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
3972 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3974 /* Final Post-processing for some Compound Methods
3976 ** The removal of any 'Sync' channel flag in the Image Compositon
3977 ** below ensures the methematical compose method is applied in a
3978 ** purely mathematical way, and only to the selected channels.
3979 ** Turn off SVG composition 'alpha blending'.
3982 case EdgeOutMorphology:
3983 case EdgeInMorphology:
3984 case TopHatMorphology:
3985 case BottomHatMorphology:
3986 if ( IfMagickTrue(verbose) )
3987 (void) FormatLocaleFile(stderr,
3988 "\n%s: Difference with original image",CommandOptionToMnemonic(
3989 MagickMorphologyOptions, method) );
3990 (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
3991 MagickTrue,0,0,exception);
3993 case EdgeMorphology:
3994 if ( IfMagickTrue(verbose) )
3995 (void) FormatLocaleFile(stderr,
3996 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
3997 MagickMorphologyOptions, method) );
3998 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
3999 MagickTrue,0,0,exception);
4000 save_image = DestroyImage(save_image); /* finished with save image */
4006 /* multi-kernel handling: re-iterate, or compose results */
4007 if ( kernel->next == (KernelInfo *) NULL )
4008 rslt_image = curr_image; /* just return the resulting image */
4009 else if ( rslt_compose == NoCompositeOp )
4010 { if ( IfMagickTrue(verbose) ) {
4011 if ( this_kernel->next != (KernelInfo *) NULL )
4012 (void) FormatLocaleFile(stderr, " (re-iterate)");
4014 (void) FormatLocaleFile(stderr, " (done)");
4016 rslt_image = curr_image; /* return result, and re-iterate */
4018 else if ( rslt_image == (Image *) NULL)
4019 { if ( IfMagickTrue(verbose) )
4020 (void) FormatLocaleFile(stderr, " (save for compose)");
4021 rslt_image = curr_image;
4022 curr_image = (Image *) image; /* continue with original image */
4025 { /* Add the new 'current' result to the composition
4027 ** The removal of any 'Sync' channel flag in the Image Compositon
4028 ** below ensures the methematical compose method is applied in a
4029 ** purely mathematical way, and only to the selected channels.
4030 ** IE: Turn off SVG composition 'alpha blending'.
4032 if ( IfMagickTrue(verbose) )
4033 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4034 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4035 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4037 curr_image = DestroyImage(curr_image);
4038 curr_image = (Image *) image; /* continue with original image */
4040 if ( IfMagickTrue(verbose) )
4041 (void) FormatLocaleFile(stderr, "\n");
4043 /* loop to the next kernel in a multi-kernel list */
4044 norm_kernel = norm_kernel->next;
4045 if ( rflt_kernel != (KernelInfo *) NULL )
4046 rflt_kernel = rflt_kernel->next;
4048 } /* End Loop 2: Loop over each kernel */
4050 } /* End Loop 1: compound method interation */
4054 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4056 if ( curr_image == rslt_image )
4057 curr_image = (Image *) NULL;
4058 if ( rslt_image != (Image *) NULL )
4059 rslt_image = DestroyImage(rslt_image);
4061 if ( curr_image == rslt_image || curr_image == image )
4062 curr_image = (Image *) NULL;
4063 if ( curr_image != (Image *) NULL )
4064 curr_image = DestroyImage(curr_image);
4065 if ( work_image != (Image *) NULL )
4066 work_image = DestroyImage(work_image);
4067 if ( save_image != (Image *) NULL )
4068 save_image = DestroyImage(save_image);
4069 if ( reflected_kernel != (KernelInfo *) NULL )
4070 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4076 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4080 % M o r p h o l o g y I m a g e %
4084 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4086 % MorphologyImage() applies a user supplied kernel to the image according to
4087 % the given mophology method.
4089 % This function applies any and all user defined settings before calling
4090 % the above internal function MorphologyApply().
4092 % User defined settings include...
4093 % * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4094 % * Kernel Scale/normalize settings ("-define convolve:scale=??")
4095 % This can also includes the addition of a scaled unity kernel.
4096 % * Show Kernel being applied ("-define showkernel=1")
4098 % Other operators that do not want user supplied options interfering,
4099 % especially "convolve:bias" and "showkernel" should use MorphologyApply()
4102 % The format of the MorphologyImage method is:
4104 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4105 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4107 % A description of each parameter follows:
4109 % o image: the image.
4111 % o method: the morphology method to be applied.
4113 % o iterations: apply the operation this many times (or no change).
4114 % A value of -1 means loop until no change found.
4115 % How this is applied may depend on the morphology method.
4116 % Typically this is a value of 1.
4118 % o kernel: An array of double representing the morphology kernel.
4119 % Warning: kernel may be normalized for the Convolve method.
4121 % o exception: return any errors or warnings in this structure.
4124 MagickExport Image *MorphologyImage(const Image *image,
4125 const MorphologyMethod method,const ssize_t iterations,
4126 const KernelInfo *kernel,ExceptionInfo *exception)
4140 curr_kernel = (KernelInfo *) kernel;
4142 compose = UndefinedCompositeOp; /* use default for method */
4144 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4145 * This is done BEFORE the ShowKernelInfo() function is called so that
4146 * users can see the results of the 'option:convolve:scale' option.
4148 if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4152 /* Get the bias value as it will be needed */
4153 artifact = GetImageArtifact(image,"convolve:bias");
4154 if ( artifact != (const char *) NULL) {
4155 if (IfMagickFalse(IsGeometry(artifact)))
4156 (void) ThrowMagickException(exception,GetMagickModule(),
4157 OptionWarning,"InvalidSetting","'%s' '%s'",
4158 "convolve:bias",artifact);
4160 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4163 /* Scale kernel according to user wishes */
4164 artifact = GetImageArtifact(image,"convolve:scale");
4165 if ( artifact != (const char *)NULL ) {
4166 if (IfMagickFalse(IsGeometry(artifact)))
4167 (void) ThrowMagickException(exception,GetMagickModule(),
4168 OptionWarning,"InvalidSetting","'%s' '%s'",
4169 "convolve:scale",artifact);
4171 if ( curr_kernel == kernel )
4172 curr_kernel = CloneKernelInfo(kernel);
4173 if (curr_kernel == (KernelInfo *) NULL)
4174 return((Image *) NULL);
4175 ScaleGeometryKernelInfo(curr_kernel, artifact);
4180 /* display the (normalized) kernel via stderr */
4181 if ( IfStringTrue(GetImageArtifact(image,"showkernel"))
4182 || IfStringTrue(GetImageArtifact(image,"convolve:showkernel"))
4183 || IfStringTrue(GetImageArtifact(image,"morphology:showkernel")) )
4184 ShowKernelInfo(curr_kernel);
4186 /* Override the default handling of multi-kernel morphology results
4187 * If 'Undefined' use the default method
4188 * If 'None' (default for 'Convolve') re-iterate previous result
4189 * Otherwise merge resulting images using compose method given.
4190 * Default for 'HitAndMiss' is 'Lighten'.
4197 artifact = GetImageArtifact(image,"morphology:compose");
4198 if ( artifact != (const char *) NULL) {
4199 parse=ParseCommandOption(MagickComposeOptions,
4200 MagickFalse,artifact);
4202 (void) ThrowMagickException(exception,GetMagickModule(),
4203 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4204 "morphology:compose",artifact);
4206 compose=(CompositeOperator)parse;
4209 /* Apply the Morphology */
4210 morphology_image = MorphologyApply(image,method,iterations,
4211 curr_kernel,compose,bias,exception);
4213 /* Cleanup and Exit */
4214 if ( curr_kernel != kernel )
4215 curr_kernel=DestroyKernelInfo(curr_kernel);
4216 return(morphology_image);
4220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4224 + R o t a t e K e r n e l I n f o %
4228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4230 % RotateKernelInfo() rotates the kernel by the angle given.
4232 % Currently it is restricted to 90 degree angles, of either 1D kernels
4233 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4234 % It will ignore usless rotations for specific 'named' built-in kernels.
4236 % The format of the RotateKernelInfo method is:
4238 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4240 % A description of each parameter follows:
4242 % o kernel: the Morphology/Convolution kernel
4244 % o angle: angle to rotate in degrees
4246 % This function is currently internal to this module only, but can be exported
4247 % to other modules if needed.
4249 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4251 /* angle the lower kernels first */
4252 if ( kernel->next != (KernelInfo *) NULL)
4253 RotateKernelInfo(kernel->next, angle);
4255 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4257 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4260 /* Modulus the angle */
4261 angle = fmod(angle, 360.0);
4265 if ( 337.5 < angle || angle <= 22.5 )
4266 return; /* Near zero angle - no change! - At least not at this time */
4268 /* Handle special cases */
4269 switch (kernel->type) {
4270 /* These built-in kernels are cylindrical kernels, rotating is useless */
4271 case GaussianKernel:
4276 case LaplacianKernel:
4277 case ChebyshevKernel:
4278 case ManhattanKernel:
4279 case EuclideanKernel:
4282 /* These may be rotatable at non-90 angles in the future */
4283 /* but simply rotating them in multiples of 90 degrees is useless */
4290 /* These only allows a +/-90 degree rotation (by transpose) */
4291 /* A 180 degree rotation is useless */
4293 if ( 135.0 < angle && angle <= 225.0 )
4295 if ( 225.0 < angle && angle <= 315.0 )
4302 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4303 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4305 if ( kernel->width == 3 && kernel->height == 3 )
4306 { /* Rotate a 3x3 square by 45 degree angle */
4307 double t = kernel->values[0];
4308 kernel->values[0] = kernel->values[3];
4309 kernel->values[3] = kernel->values[6];
4310 kernel->values[6] = kernel->values[7];
4311 kernel->values[7] = kernel->values[8];
4312 kernel->values[8] = kernel->values[5];
4313 kernel->values[5] = kernel->values[2];
4314 kernel->values[2] = kernel->values[1];
4315 kernel->values[1] = t;
4316 /* rotate non-centered origin */
4317 if ( kernel->x != 1 || kernel->y != 1 ) {
4319 x = (ssize_t) kernel->x-1;
4320 y = (ssize_t) kernel->y-1;
4321 if ( x == y ) x = 0;
4322 else if ( x == 0 ) x = -y;
4323 else if ( x == -y ) y = 0;
4324 else if ( y == 0 ) y = x;
4325 kernel->x = (ssize_t) x+1;
4326 kernel->y = (ssize_t) y+1;
4328 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4329 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4332 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4334 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4336 if ( kernel->width == 1 || kernel->height == 1 )
4337 { /* Do a transpose of a 1 dimensional kernel,
4338 ** which results in a fast 90 degree rotation of some type.
4342 t = (ssize_t) kernel->width;
4343 kernel->width = kernel->height;
4344 kernel->height = (size_t) t;
4346 kernel->x = kernel->y;
4348 if ( kernel->width == 1 ) {
4349 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4350 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4352 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4353 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4356 else if ( kernel->width == kernel->height )
4357 { /* Rotate a square array of values by 90 degrees */
4361 register MagickRealType
4365 for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--)
4366 for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--)
4367 { t = k[i+j*kernel->width];
4368 k[i+j*kernel->width] = k[j+x*kernel->width];
4369 k[j+x*kernel->width] = k[x+y*kernel->width];
4370 k[x+y*kernel->width] = k[y+i*kernel->width];
4371 k[y+i*kernel->width] = t;
4374 /* rotate the origin - relative to center of array */
4375 { register ssize_t x,y;
4376 x = (ssize_t) (kernel->x*2-kernel->width+1);
4377 y = (ssize_t) (kernel->y*2-kernel->height+1);
4378 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4379 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4381 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4382 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4385 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4387 if ( 135.0 < angle && angle <= 225.0 )
4389 /* For a 180 degree rotation - also know as a reflection
4390 * This is actually a very very common operation!
4391 * Basically all that is needed is a reversal of the kernel data!
4392 * And a reflection of the origon
4397 register MagickRealType
4405 j=(ssize_t) (kernel->width*kernel->height-1);
4406 for (i=0; i < j; i++, j--)
4407 t=k[i], k[i]=k[j], k[j]=t;
4409 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4410 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4411 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4412 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4414 /* At this point angle should at least between -45 (315) and +45 degrees
4415 * In the future some form of non-orthogonal angled rotates could be
4416 * performed here, posibily with a linear kernel restriction.
4423 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4427 % S c a l e G e o m e t r y K e r n e l I n f o %
4431 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4433 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4434 % provided as a "-set option:convolve:scale {geometry}" user setting,
4435 % and modifies the kernel according to the parsed arguments of that setting.
4437 % The first argument (and any normalization flags) are passed to
4438 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4439 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4440 % into the scaled/normalized kernel.
4442 % The format of the ScaleGeometryKernelInfo method is:
4444 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4445 % const double scaling_factor,const MagickStatusType normalize_flags)
4447 % A description of each parameter follows:
4449 % o kernel: the Morphology/Convolution kernel to modify
4452 % The geometry string to parse, typically from the user provided
4453 % "-set option:convolve:scale {geometry}" setting.
4456 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4457 const char *geometry)
4465 SetGeometryInfo(&args);
4466 flags = ParseGeometry(geometry, &args);
4469 /* For Debugging Geometry Input */
4470 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4471 flags, args.rho, args.sigma, args.xi, args.psi );
4474 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4475 args.rho *= 0.01, args.sigma *= 0.01;
4477 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4479 if ( (flags & SigmaValue) == 0 )
4482 /* Scale/Normalize the input kernel */
4483 ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4485 /* Add Unity Kernel, for blending with original */
4486 if ( (flags & SigmaValue) != 0 )
4487 UnityAddKernelInfo(kernel, args.sigma);
4492 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4496 % S c a l e K e r n e l I n f o %
4500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4502 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4503 % without normalization of the sum of the kernel values (as per given flags).
4505 % By default (no flags given) the values within the kernel is scaled
4506 % directly using given scaling factor without change.
4508 % If either of the two 'normalize_flags' are given the kernel will first be
4509 % normalized and then further scaled by the scaling factor value given.
4511 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4512 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4513 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4514 % non-HDRI versions of IM this may cause images to have any negative results
4515 % clipped, unless some 'bias' is used.
4517 % More specifically. Kernels which only contain positive values (such as a
4518 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4519 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4521 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4522 % the kernel will be scaled by the absolute of the sum of kernel values, so
4523 % that it will generally fall within the +/- 1.0 range.
4525 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4526 % will be scaled by just the sum of the postive values, so that its output
4527 % range will again fall into the +/- 1.0 range.
4529 % For special kernels designed for locating shapes using 'Correlate', (often
4530 % only containing +1 and -1 values, representing foreground/brackground
4531 % matching) a special normalization method is provided to scale the positive
4532 % values separately to those of the negative values, so the kernel will be
4533 % forced to become a zero-sum kernel better suited to such searches.
4535 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4536 % attributes within the kernel structure have been correctly set during the
4539 % NOTE: The values used for 'normalize_flags' have been selected specifically
4540 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4541 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4543 % The format of the ScaleKernelInfo method is:
4545 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4546 % const MagickStatusType normalize_flags )
4548 % A description of each parameter follows:
4550 % o kernel: the Morphology/Convolution kernel
4553 % multiply all values (after normalization) by this factor if not
4554 % zero. If the kernel is normalized regardless of any flags.
4556 % o normalize_flags:
4557 % GeometryFlags defining normalization method to use.
4558 % specifically: NormalizeValue, CorrelateNormalizeValue,
4559 % and/or PercentValue
4562 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4563 const double scaling_factor,const GeometryFlags normalize_flags)
4572 /* do the other kernels in a multi-kernel list first */
4573 if ( kernel->next != (KernelInfo *) NULL)
4574 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4576 /* Normalization of Kernel */
4578 if ( (normalize_flags&NormalizeValue) != 0 ) {
4579 if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4580 /* non-zero-summing kernel (generally positive) */
4581 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4583 /* zero-summing kernel */
4584 pos_scale = kernel->positive_range;
4586 /* Force kernel into a normalized zero-summing kernel */
4587 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4588 pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4589 ? kernel->positive_range : 1.0;
4590 neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4591 ? -kernel->negative_range : 1.0;
4594 neg_scale = pos_scale;
4596 /* finialize scaling_factor for positive and negative components */
4597 pos_scale = scaling_factor/pos_scale;
4598 neg_scale = scaling_factor/neg_scale;
4600 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4601 if ( ! IfNaN(kernel->values[i]) )
4602 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4604 /* convolution output range */
4605 kernel->positive_range *= pos_scale;
4606 kernel->negative_range *= neg_scale;
4607 /* maximum and minimum values in kernel */
4608 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4609 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4611 /* swap kernel settings if user's scaling factor is negative */
4612 if ( scaling_factor < MagickEpsilon ) {
4614 t = kernel->positive_range;
4615 kernel->positive_range = kernel->negative_range;
4616 kernel->negative_range = t;
4617 t = kernel->maximum;
4618 kernel->maximum = kernel->minimum;
4619 kernel->minimum = 1;
4626 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4630 % S h o w K e r n e l I n f o %
4634 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4636 % ShowKernelInfo() outputs the details of the given kernel defination to
4637 % standard error, generally due to a users 'showkernel' option request.
4639 % The format of the ShowKernel method is:
4641 % void ShowKernelInfo(const KernelInfo *kernel)
4643 % A description of each parameter follows:
4645 % o kernel: the Morphology/Convolution kernel
4648 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4656 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4658 (void) FormatLocaleFile(stderr, "Kernel");
4659 if ( kernel->next != (KernelInfo *) NULL )
4660 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4661 (void) FormatLocaleFile(stderr, " \"%s",
4662 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4663 if ( fabs(k->angle) >= MagickEpsilon )
4664 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4665 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4666 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4667 (void) FormatLocaleFile(stderr,
4668 " with values from %.*lg to %.*lg\n",
4669 GetMagickPrecision(), k->minimum,
4670 GetMagickPrecision(), k->maximum);
4671 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4672 GetMagickPrecision(), k->negative_range,
4673 GetMagickPrecision(), k->positive_range);
4674 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4675 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4676 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4677 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4679 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4680 GetMagickPrecision(), k->positive_range+k->negative_range);
4681 for (i=v=0; v < k->height; v++) {
4682 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4683 for (u=0; u < k->width; u++, i++)
4684 if ( IfNaN(k->values[i]) )
4685 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4687 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4688 GetMagickPrecision(), (double) k->values[i]);
4689 (void) FormatLocaleFile(stderr,"\n");
4695 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4699 % U n i t y A d d K e r n a l I n f o %
4703 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4705 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4706 % to the given pre-scaled and normalized Kernel. This in effect adds that
4707 % amount of the original image into the resulting convolution kernel. This
4708 % value is usually provided by the user as a percentage value in the
4709 % 'convolve:scale' setting.
4711 % The resulting effect is to convert the defined kernels into blended
4712 % soft-blurs, unsharp kernels or into sharpening kernels.
4714 % The format of the UnityAdditionKernelInfo method is:
4716 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4718 % A description of each parameter follows:
4720 % o kernel: the Morphology/Convolution kernel
4723 % scaling factor for the unity kernel to be added to
4727 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4730 /* do the other kernels in a multi-kernel list first */
4731 if ( kernel->next != (KernelInfo *) NULL)
4732 UnityAddKernelInfo(kernel->next, scale);
4734 /* Add the scaled unity kernel to the existing kernel */
4735 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4736 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4742 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4746 % Z e r o K e r n e l N a n s %
4750 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4752 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4753 % the kernel with a zero value. This is typically done when the kernel will
4754 % be used in special hardware (GPU) convolution processors, to simply
4757 % The format of the ZeroKernelNans method is:
4759 % void ZeroKernelNans (KernelInfo *kernel)
4761 % A description of each parameter follows:
4763 % o kernel: the Morphology/Convolution kernel
4766 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4771 /* do the other kernels in a multi-kernel list first */
4772 if ( kernel->next != (KernelInfo *) NULL)
4773 ZeroKernelNans(kernel->next);
4775 for (i=0; i < (kernel->width*kernel->height); i++)
4776 if ( IfNaN(kernel->values[i]) )
4777 kernel->values[i] = 0.0;