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-2010 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 "magick/studio.h"
53 #include "magick/artifact.h"
54 #include "magick/cache-view.h"
55 #include "magick/color-private.h"
56 #include "magick/enhance.h"
57 #include "magick/exception.h"
58 #include "magick/exception-private.h"
59 #include "magick/gem.h"
60 #include "magick/hashmap.h"
61 #include "magick/image.h"
62 #include "magick/image-private.h"
63 #include "magick/list.h"
64 #include "magick/magick.h"
65 #include "magick/memory_.h"
66 #include "magick/monitor-private.h"
67 #include "magick/morphology.h"
68 #include "magick/option.h"
69 #include "magick/pixel-private.h"
70 #include "magick/prepress.h"
71 #include "magick/quantize.h"
72 #include "magick/registry.h"
73 #include "magick/semaphore.h"
74 #include "magick/splay-tree.h"
75 #include "magick/statistic.h"
76 #include "magick/string_.h"
77 #include "magick/string-private.h"
78 #include "magick/token.h"
81 The following test is for special floating point numbers of value NaN (not
82 a number), that may be used within a Kernel Definition. NaN's are defined
83 as part of the IEEE standard for floating point number representation.
85 These are used a Kernel value of NaN means that that kernel position is not
86 part of the normal convolution or morphology process, and thus allowing the
87 use of 'shaped' kernels.
89 Special properities two NaN's are never equal, even if they are from the
90 same variable That is the IsNaN() macro is only true if the value is NaN.
92 #define IsNan(a) ((a)!=(a))
95 Other global definitions used by module.
97 static inline double MagickMin(const double x,const double y)
99 return( x < y ? x : y);
101 static inline double MagickMax(const double x,const double y)
103 return( x > y ? x : y);
105 #define Minimize(assign,value) assign=MagickMin(assign,value)
106 #define Maximize(assign,value) assign=MagickMax(assign,value)
108 /* Currently these are only internal to this module */
110 ExpandKernelInfo(KernelInfo *, double),
111 RotateKernelInfo(KernelInfo *, double);
114 /* Quick function to find last kernel in a kernel list */
115 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
117 while (kernel->next != (KernelInfo *) NULL)
118 kernel = kernel->next;
124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 % A c q u i r e K e r n e l I n f o %
132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134 % AcquireKernelInfo() takes the given string (generally supplied by the
135 % user) and converts it into a Morphology/Convolution Kernel. This allows
136 % users to specify a kernel from a number of pre-defined kernels, or to fully
137 % specify their own kernel for a specific Convolution or Morphology
140 % The kernel so generated can be any rectangular array of floating point
141 % values (doubles) with the 'control point' or 'pixel being affected'
142 % anywhere within that array of values.
144 % Previously IM was restricted to a square of odd size using the exact
145 % center as origin, this is no longer the case, and any rectangular kernel
146 % with any value being declared the origin. This in turn allows the use of
147 % highly asymmetrical kernels.
149 % The floating point values in the kernel can also include a special value
150 % known as 'nan' or 'not a number' to indicate that this value is not part
151 % of the kernel array. This allows you to shaped the kernel within its
152 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
153 % shape. However at least one non-nan value must be provided for correct
154 % working of a kernel.
156 % The returned kernel should be freed using the DestroyKernelInfo() when you
157 % are finished with it. Do not free this memory yourself.
159 % Input kernel defintion strings can consist of any of three types.
162 % Select from one of the built in kernels, using the name and
163 % geometry arguments supplied. See AcquireKernelBuiltIn()
165 % "WxH[+X+Y]:num, num, num ..."
166 % a kernel of size W by H, with W*H floating point numbers following.
167 % the 'center' can be optionally be defined at +X+Y (such that +0+0
168 % is top left corner). If not defined the pixel in the center, for
169 % odd sizes, or to the immediate top or left of center for even sizes
170 % is automatically selected.
172 % "num, num, num, num, ..."
173 % list of floating point numbers defining an 'old style' odd sized
174 % square kernel. At least 9 values should be provided for a 3x3
175 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
176 % Values can be space or comma separated. This is not recommended.
178 % Note that 'name' kernels will start with an alphabetic character while the
179 % new kernel specification has a ':' character in its specification string.
180 % If neither is the case, it is assumed an old style of a simple list of
181 % numbers generating a odd-sized square kernel has been given.
183 % You can define a 'list of kernels' which can be used by some morphology
184 % operators A list is defined as a semi-colon seperated list kernels.
186 % " kernel ; kernel ; kernel ; "
188 % Extra ';' characters are simply ignored.
190 % The format of the AcquireKernal method is:
192 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
194 % A description of each parameter follows:
196 % o kernel_string: the Morphology/Convolution kernel wanted.
200 /* This was separated so that it could be used as a separate
201 ** array input handling function, such as for -color-matrix
203 static KernelInfo *ParseKernelArray(const char *kernel_string)
209 token[MaxTextExtent];
219 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
221 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
222 if (kernel == (KernelInfo *)NULL)
224 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
225 kernel->minimum = kernel->maximum = 0.0;
226 kernel->negative_range = kernel->positive_range = 0.0;
227 kernel->type = UserDefinedKernel;
228 kernel->next = (KernelInfo *) NULL;
229 kernel->signature = MagickSignature;
231 /* find end of this specific kernel definition string */
232 end = strchr(kernel_string, ';');
233 if ( end == (char *) NULL )
234 end = strchr(kernel_string, '\0');
236 /* Has a ':' in argument - New user kernel specification */
237 p = strchr(kernel_string, ':');
238 if ( p != (char *) NULL && p < end)
246 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
247 memcpy(token, kernel_string, (size_t) (p-kernel_string));
248 token[p-kernel_string] = '\0';
249 SetGeometryInfo(&args);
250 flags = ParseGeometry(token, &args);
252 /* Size handling and checks of geometry settings */
253 if ( (flags & WidthValue) == 0 ) /* if no width then */
254 args.rho = args.sigma; /* then width = height */
255 if ( args.rho < 1.0 ) /* if width too small */
256 args.rho = 1.0; /* then width = 1 */
257 if ( args.sigma < 1.0 ) /* if height too small */
258 args.sigma = args.rho; /* then height = width */
259 kernel->width = (unsigned long)args.rho;
260 kernel->height = (unsigned long)args.sigma;
262 /* Offset Handling and Checks */
263 if ( args.xi < 0.0 || args.psi < 0.0 )
264 return(DestroyKernelInfo(kernel));
265 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
266 : (long) (kernel->width-1)/2;
267 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
268 : (long) (kernel->height-1)/2;
269 if ( kernel->x >= (long) kernel->width ||
270 kernel->y >= (long) kernel->height )
271 return(DestroyKernelInfo(kernel));
273 p++; /* advance beyond the ':' */
276 { /* ELSE - Old old specification, forming odd-square kernel */
277 /* count up number of values given */
278 p=(const char *) kernel_string;
279 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
280 p++; /* ignore "'" chars for convolve filter usage - Cristy */
281 for (i=0; p < end; i++)
283 GetMagickToken(p,&p,token);
285 GetMagickToken(p,&p,token);
287 /* set the size of the kernel - old sized square */
288 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
289 kernel->x = kernel->y = (long) (kernel->width-1)/2;
290 p=(const char *) kernel_string;
291 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
292 p++; /* ignore "'" chars for convolve filter usage - Cristy */
295 /* Read in the kernel values from rest of input string argument */
296 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
297 kernel->height*sizeof(double));
298 if (kernel->values == (double *) NULL)
299 return(DestroyKernelInfo(kernel));
301 kernel->minimum = +MagickHuge;
302 kernel->maximum = -MagickHuge;
303 kernel->negative_range = kernel->positive_range = 0.0;
305 for (i=0; (i < (long) (kernel->width*kernel->height)) && (p < end); i++)
307 GetMagickToken(p,&p,token);
309 GetMagickToken(p,&p,token);
310 if ( LocaleCompare("nan",token) == 0
311 || LocaleCompare("-",token) == 0 ) {
312 kernel->values[i] = nan; /* do not include this value in kernel */
315 kernel->values[i] = StringToDouble(token);
316 ( kernel->values[i] < 0)
317 ? ( kernel->negative_range += kernel->values[i] )
318 : ( kernel->positive_range += kernel->values[i] );
319 Minimize(kernel->minimum, kernel->values[i]);
320 Maximize(kernel->maximum, kernel->values[i]);
324 /* sanity check -- no more values in kernel definition */
325 GetMagickToken(p,&p,token);
326 if ( *token != '\0' && *token != ';' && *token != '\'' )
327 return(DestroyKernelInfo(kernel));
330 /* this was the old method of handling a incomplete kernel */
331 if ( i < (long) (kernel->width*kernel->height) ) {
332 Minimize(kernel->minimum, kernel->values[i]);
333 Maximize(kernel->maximum, kernel->values[i]);
334 for ( ; i < (long) (kernel->width*kernel->height); i++)
335 kernel->values[i]=0.0;
338 /* Number of values for kernel was not enough - Report Error */
339 if ( i < (long) (kernel->width*kernel->height) )
340 return(DestroyKernelInfo(kernel));
343 /* check that we recieved at least one real (non-nan) value! */
344 if ( kernel->minimum == MagickHuge )
345 return(DestroyKernelInfo(kernel));
350 static KernelInfo *ParseNamedKernel(const char *kernel_string)
353 token[MaxTextExtent];
368 /* Parse special 'named' kernel */
369 GetMagickToken(kernel_string,&p,token);
370 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
371 if ( type < 0 || type == UserDefinedKernel )
372 return((KernelInfo *)NULL); /* not a valid named kernel */
374 while (((isspace((int) ((unsigned char) *p)) != 0) ||
375 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
378 end = strchr(p, ';'); /* end of this kernel defintion */
379 if ( end == (char *) NULL )
380 end = strchr(p, '\0');
382 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
383 memcpy(token, p, (size_t) (end-p));
385 SetGeometryInfo(&args);
386 flags = ParseGeometry(token, &args);
389 /* For Debugging Geometry Input */
390 fprintf(stderr, "Geometry = %04x : %lf x %lf %+lf %+lf\n",
391 flags, args.rho, args.sigma, args.xi, args.psi );
394 /* special handling of missing values in input string */
396 case RectangleKernel:
397 if ( (flags & WidthValue) == 0 ) /* if no width then */
398 args.rho = args.sigma; /* then width = height */
399 if ( args.rho < 1.0 ) /* if width too small */
400 args.rho = 3; /* then width = 3 */
401 if ( args.sigma < 1.0 ) /* if height too small */
402 args.sigma = args.rho; /* then height = width */
403 if ( (flags & XValue) == 0 ) /* center offset if not defined */
404 args.xi = (double)(((long)args.rho-1)/2);
405 if ( (flags & YValue) == 0 )
406 args.psi = (double)(((long)args.sigma-1)/2);
413 /* If no scale given (a 0 scale is valid! - set it to 1.0 */
414 if ( (flags & HeightValue) == 0 )
417 case ChebyshevKernel:
418 case ManhattenKernel:
419 case EuclideanKernel:
420 if ( (flags & HeightValue) == 0 )
421 args.sigma = 100.0; /* default distance scaling */
422 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
423 args.sigma = QuantumRange/args.sigma; /* maximum pixel distance */
424 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
425 args.sigma *= QuantumRange/100.0; /* percentage of color range */
431 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
434 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
443 token[MaxTextExtent];
452 kernel = last_kernel = NULL;
455 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
457 /* ignore multiple ';' following each other */
458 if ( *token != ';' ) {
460 /* tokens starting with alpha is a Named kernel */
461 if (isalpha((int) *token) == 0)
462 new_kernel = ParseKernelArray(p);
463 else /* otherwise a user defined kernel array */
464 new_kernel = ParseNamedKernel(p);
466 /* Error handling -- this is not proper error handling! */
467 if ( new_kernel == (KernelInfo *) NULL ) {
468 fprintf(stderr, "Failed to parse kernel number #%lu\n", kernel_number);
469 if ( kernel != (KernelInfo *) NULL )
470 kernel=DestroyKernelInfo(kernel);
471 return((KernelInfo *) NULL);
474 /* initialise or append the kernel list */
475 if ( kernel == (KernelInfo *) NULL )
478 last_kernel->next = new_kernel;
479 last_kernel = LastKernelInfo(new_kernel);
482 /* look for the next kernel in list */
484 if ( p == (char *) NULL )
494 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
498 % A c q u i r e K e r n e l B u i l t I n %
502 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
504 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
505 % kernels used for special purposes such as gaussian blurring, skeleton
506 % pruning, and edge distance determination.
508 % They take a KernelType, and a set of geometry style arguments, which were
509 % typically decoded from a user supplied string, or from a more complex
510 % Morphology Method that was requested.
512 % The format of the AcquireKernalBuiltIn method is:
514 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
515 % const GeometryInfo args)
517 % A description of each parameter follows:
519 % o type: the pre-defined type of kernel wanted
521 % o args: arguments defining or modifying the kernel
523 % Convolution Kernels
525 % Gaussian:{radius},{sigma}
526 % Generate a two-dimentional gaussian kernel, as used by -gaussian.
527 % A sigma is required.
529 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
530 % the final size of the resulting kernel to a square 2*radius+1 in size.
531 % The radius should be at least 2 times that of the sigma value, or
532 % sever clipping and aliasing may result. If not given or set to 0 the
533 % radius will be determined so as to produce the best minimal error
534 % result, which is usally much larger than is normally needed.
536 % Blur:{radius},{sigma},{angle}
537 % As per Gaussian, but generates a 1 dimensional or linear gaussian
538 % blur, at the angle given (current restricted to orthogonal angles).
539 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
540 % Angle can be rotated in multiples of 90 degrees.
542 % Note that two such blurs perpendicular to each other is equivelent to
543 % the far large "Gaussian" kernel, but much faster to apply. This is how
544 % the -blur operator works.
546 % Comet:{width},{sigma},{angle}
547 % Blur in one direction only, much like how a bright object leaves
548 % a comet like trail. The Kernel is actually half a gaussian curve,
549 % Adding two such blurs in opposite directions produces a Blur Kernel.
550 % Angle can be rotated in multiples of 90 degrees.
552 % Note that the first argument is the width of the kernel and not the
553 % radius of the kernel.
555 % # Still to be implemented...
557 % # DOG:{radius},{sigma1},{sigma2}
558 % # Difference of two Gaussians
562 % # Set kernel values using a resize filter, and given scale (sigma)
563 % # Cylindrical or Linear. Is this posible with an image?
566 % Named Constant Convolution Kernels
569 % The 3x3 sobel convolution kernel. Angle may be given in multiples
570 % of 45 degrees. Kernel is unscaled by default so some normalization
571 % may be required to ensure results are not clipped.
572 % Default kernel is -1,0,1
577 % Generate Lapacian kernel of the type specified. (1 is the default)
578 % Type 0 default square laplacian 3x3: all -1's with central 8
579 % Type 1 3x3: central 4 edge -1 corner 0
580 % Type 2 3x3: central 4 edge 1 corner -2
581 % Type 3 a 5x5 laplacian
582 % Type 4 a 7x7 laplacian
586 % Rectangle:{geometry}
587 % Simply generate a rectangle of 1's with the size given. You can also
588 % specify the location of the 'control point', otherwise the closest
589 % pixel to the center of the rectangle is selected.
591 % Properly centered and odd sized rectangles work the best.
593 % Diamond:[{radius}[,{scale}]]
594 % Generate a diamond shaped kernel with given radius to the points.
595 % Kernel size will again be radius*2+1 square and defaults to radius 1,
596 % generating a 3x3 kernel that is slightly larger than a square.
598 % Square:[{radius}[,{scale}]]
599 % Generate a square shaped kernel of size radius*2+1, and defaulting
600 % to a 3x3 (radius 1).
602 % Note that using a larger radius for the "Square" or the "Diamond"
603 % is also equivelent to iterating the basic morphological method
604 % that many times. However However iterating with the smaller radius 1
605 % default is actually faster than using a larger kernel radius.
607 % Disk:[{radius}[,{scale}]]
608 % Generate a binary disk of the radius given, radius may be a float.
609 % Kernel size will be ceil(radius)*2+1 square.
610 % NOTE: Here are some disk shapes of specific interest
611 % "disk:1" => "diamond" or "cross:1"
612 % "disk:1.5" => "square"
613 % "disk:2" => "diamond:2"
614 % "disk:2.5" => a general disk shape of radius 2
615 % "disk:2.9" => "square:2"
616 % "disk:3.5" => default - octagonal/disk shape of radius 3
617 % "disk:4.2" => roughly octagonal shape of radius 4
618 % "disk:4.3" => a general disk shape of radius 4
619 % After this all the kernel shape becomes more and more circular.
621 % Because a "disk" is more circular when using a larger radius, using a
622 % larger radius is preferred over iterating the morphological operation.
624 % Plus:[{radius}[,{scale}]]
625 % Cross:[{radius}[,{scale}]]
626 % Generate a kernel in the shape of a 'plus' or a cross. The length of
627 % each arm is also the radius, which defaults to 2.
629 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
631 % This kernel is not a good general morphological kernel, but is used
632 % more for highlighting and marking any single pixels in an image using,
633 % a "Dilate" or "Erode" method as appropriate.
635 % For the same reasons iterating these kernels does not produce the
636 % same result as using a larger radius for the symbol.
638 % Hit and Miss Kernels
640 % Peak:radius1,radius2
641 % Find a foreground inside a background ring of the given radii.
643 % Find corners of a binary shape
645 % Find end points of lines (for pruning a skeletion)
647 % Find three line junctions (in a skeletion)
649 % Octagonal thicken kernel, to generate convex hulls of 45 degrees
651 % Thinning kernel, which leaves behind a skeletion of a shape
653 % Distance Measuring Kernels
655 % Chebyshev:[{radius}][x{scale}[%!]]
656 % Manhatten:[{radius}][x{scale}[%!]]
657 % Euclidean:[{radius}][x{scale}[%!]]
659 % Different types of distance measuring methods, which are used with the
660 % a 'Distance' morphology method for generating a gradient based on
661 % distance from an edge of a binary shape, though there is a technique
662 % for handling a anti-aliased shape.
664 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
665 % one to any neighbour, orthogonal or diagonal. One why of thinking of
666 % it is the number of squares a 'King' or 'Queen' in chess needs to
667 % traverse reach any other position on a chess board. It results in a
668 % 'square' like distance function, but one where diagonals are closer
671 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
672 % Cab metric), is the distance needed when you can only travel in
673 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
674 % in chess would travel. It results in a diamond like distances, where
675 % diagonals are further than expected.
677 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
678 % However by default the kernel size only has a radius of 1, which
679 % limits the distance to 'Knight' like moves, with only orthogonal and
680 % diagonal measurements being correct. As such for the default kernel
681 % you will get octagonal like distance function, which is reasonally
684 % However if you use a larger radius such as "Euclidean:4" you will
685 % get a much smoother distance gradient from the edge of the shape.
686 % Of course a larger kernel is slower to use, and generally not needed.
688 % To allow the use of fractional distances that you get with diagonals
689 % the actual distance is scaled by a fixed value which the user can
690 % provide. This is not actually nessary for either ""Chebyshev" or
691 % "Manhatten" distance kernels, but is done for all three distance
692 % kernels. If no scale is provided it is set to a value of 100,
693 % allowing for a maximum distance measurement of 655 pixels using a Q16
694 % version of IM, from any edge. However for small images this can
695 % result in quite a dark gradient.
697 % See the 'Distance' Morphological Method, for information of how it is
702 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
703 const GeometryInfo *args)
716 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
718 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
719 if (kernel == (KernelInfo *) NULL)
721 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
722 kernel->minimum = kernel->maximum = 0.0;
723 kernel->negative_range = kernel->positive_range = 0.0;
725 kernel->next = (KernelInfo *) NULL;
726 kernel->signature = MagickSignature;
729 /* Convolution Kernels */
732 sigma = fabs(args->sigma);
734 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
736 kernel->width = kernel->height =
737 GetOptimalKernelWidth2D(args->rho,sigma);
738 kernel->x = kernel->y = (long) (kernel->width-1)/2;
739 kernel->negative_range = kernel->positive_range = 0.0;
740 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
741 kernel->height*sizeof(double));
742 if (kernel->values == (double *) NULL)
743 return(DestroyKernelInfo(kernel));
745 sigma = 2.0*sigma*sigma; /* simplify the expression */
746 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
747 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
748 kernel->positive_range += (
750 exp(-((double)(u*u+v*v))/sigma)
751 /* / (MagickPI*sigma) */ );
753 kernel->maximum = kernel->values[
754 kernel->y*kernel->width+kernel->x ];
756 /* Normalize the 2D Gaussian Kernel
758 ** As it is normalized the divisor in the above kernel generator is
759 ** not needed, so is not done above.
761 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
767 sigma = fabs(args->sigma);
769 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
771 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
772 kernel->x = (long) (kernel->width-1)/2;
775 kernel->negative_range = kernel->positive_range = 0.0;
776 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
777 kernel->height*sizeof(double));
778 if (kernel->values == (double *) NULL)
779 return(DestroyKernelInfo(kernel));
783 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
784 ** It generates a gaussian 3 times the width, and compresses it into
785 ** the expected range. This produces a closer normalization of the
786 ** resulting kernel, especially for very low sigma values.
787 ** As such while wierd it is prefered.
789 ** I am told this method originally came from Photoshop.
791 sigma *= KernelRank; /* simplify expanded curve */
792 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
793 (void) ResetMagickMemory(kernel->values,0, (size_t)
794 kernel->width*sizeof(double));
795 for ( u=-v; u <= v; u++) {
796 kernel->values[(u+v)/KernelRank] +=
797 exp(-((double)(u*u))/(2.0*sigma*sigma))
798 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
800 for (i=0; i < (long) kernel->width; i++)
801 kernel->positive_range += kernel->values[i];
803 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
804 kernel->positive_range += (
806 exp(-((double)(u*u))/(2.0*sigma*sigma))
807 /* / (MagickSQ2PI*sigma) */ );
810 kernel->maximum = kernel->values[ kernel->x ];
811 /* Note that neither methods above generate a normalized kernel,
812 ** though it gets close. The kernel may be 'clipped' by a user defined
813 ** radius, producing a smaller (darker) kernel. Also for very small
814 ** sigma's (> 0.1) the central value becomes larger than one, and thus
815 ** producing a very bright kernel.
818 /* Normalize the 1D Gaussian Kernel
820 ** As it is normalized the divisor in the above kernel generator is
821 ** not needed, so is not done above.
823 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
825 /* rotate the kernel by given angle */
826 RotateKernelInfo(kernel, args->xi);
831 sigma = fabs(args->sigma);
833 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
835 if ( args->rho < 1.0 )
836 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
838 kernel->width = (unsigned long)args->rho;
839 kernel->x = kernel->y = 0;
841 kernel->negative_range = kernel->positive_range = 0.0;
842 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
843 kernel->height*sizeof(double));
844 if (kernel->values == (double *) NULL)
845 return(DestroyKernelInfo(kernel));
847 /* A comet blur is half a gaussian curve, so that the object is
848 ** blurred in one direction only. This may not be quite the right
849 ** curve to use so may change in the future. The function must be
850 ** normalised after generation, which also resolves any clipping.
854 sigma *= KernelRank; /* simplify expanded curve */
855 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
856 (void) ResetMagickMemory(kernel->values,0, (size_t)
857 kernel->width*sizeof(double));
858 for ( u=0; u < v; u++) {
859 kernel->values[u/KernelRank] +=
860 exp(-((double)(u*u))/(2.0*sigma*sigma))
861 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
863 for (i=0; i < (long) kernel->width; i++)
864 kernel->positive_range += kernel->values[i];
866 for ( i=0; i < (long) kernel->width; i++)
867 kernel->positive_range += (
869 exp(-((double)(i*i))/(2.0*sigma*sigma))
870 /* / (MagickSQ2PI*sigma) */ );
873 kernel->maximum = kernel->values[0];
875 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
876 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
879 /* Convolution Kernels - Well Known Constants */
881 { kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
882 kernel=ParseKernelArray("3x3:-1,0,1 -2,0,2 -1,0,1");
883 if (kernel == (KernelInfo *) NULL)
886 RotateKernelInfo(kernel, args->rho); /* Rotate by angle */
889 case LaplacianKernel:
890 { kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
891 switch ( (int) args->rho ) {
893 default: /* default laplacian 'edge' filter */
894 kernel=ParseKernelArray("3x3: -1,-1,-1 -1,8,-1 -1,-1,-1");
897 kernel=ParseKernelArray("3x3: 0,-1,0 -1,4,-1 0,-1,0");
900 kernel=ParseKernelArray("3x3: -2,1,-2 1,4,1 -2,1,-2");
902 case 3: /* a 5x5 laplacian */
903 kernel=ParseKernelArray(
904 "5x5: -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");
906 case 4: /* a 7x7 laplacian */
907 kernel=ParseKernelArray(
908 "7x7:-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" );
911 if (kernel == (KernelInfo *) NULL)
916 /* Boolean Kernels */
917 case RectangleKernel:
921 if ( type == SquareKernel )
924 kernel->width = kernel->height = 3; /* default radius = 1 */
926 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
927 kernel->x = kernel->y = (long) (kernel->width-1)/2;
931 /* NOTE: user defaults set in "AcquireKernelInfo()" */
932 if ( args->rho < 1.0 || args->sigma < 1.0 )
933 return(DestroyKernelInfo(kernel)); /* invalid args given */
934 kernel->width = (unsigned long)args->rho;
935 kernel->height = (unsigned long)args->sigma;
936 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
937 args->psi < 0.0 || args->psi > (double)kernel->height )
938 return(DestroyKernelInfo(kernel)); /* invalid args given */
939 kernel->x = (long) args->xi;
940 kernel->y = (long) args->psi;
943 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
944 kernel->height*sizeof(double));
945 if (kernel->values == (double *) NULL)
946 return(DestroyKernelInfo(kernel));
948 /* set all kernel values to scale given */
949 u=(long) kernel->width*kernel->height;
950 for ( i=0; i < u; i++)
951 kernel->values[i] = scale;
952 kernel->minimum = kernel->maximum = scale; /* a flat shape */
953 kernel->positive_range = scale*u;
959 kernel->width = kernel->height = 3; /* default radius = 1 */
961 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
962 kernel->x = kernel->y = (long) (kernel->width-1)/2;
964 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
965 kernel->height*sizeof(double));
966 if (kernel->values == (double *) NULL)
967 return(DestroyKernelInfo(kernel));
969 /* set all kernel values within diamond area to scale given */
970 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
971 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
972 if ((labs(u)+labs(v)) <= (long)kernel->x)
973 kernel->positive_range += kernel->values[i] = args->sigma;
975 kernel->values[i] = nan;
976 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
984 limit = (long)(args->rho*args->rho);
985 if (args->rho < 0.1) /* default radius approx 3.5 */
986 kernel->width = kernel->height = 7L, limit = 10L;
988 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
989 kernel->x = kernel->y = (long) (kernel->width-1)/2;
991 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
992 kernel->height*sizeof(double));
993 if (kernel->values == (double *) NULL)
994 return(DestroyKernelInfo(kernel));
996 /* set all kernel values within disk area to scale given */
997 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
998 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
999 if ((u*u+v*v) <= limit)
1000 kernel->positive_range += kernel->values[i] = args->sigma;
1002 kernel->values[i] = nan;
1003 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1008 if (args->rho < 1.0)
1009 kernel->width = kernel->height = 5; /* default radius 2 */
1011 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
1012 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1014 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1015 kernel->height*sizeof(double));
1016 if (kernel->values == (double *) NULL)
1017 return(DestroyKernelInfo(kernel));
1019 /* set all kernel values along axises to given scale */
1020 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1021 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1022 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1023 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1024 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1029 if (args->rho < 1.0)
1030 kernel->width = kernel->height = 5; /* default radius 2 */
1032 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
1033 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1035 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1036 kernel->height*sizeof(double));
1037 if (kernel->values == (double *) NULL)
1038 return(DestroyKernelInfo(kernel));
1040 /* set all kernel values along axises to given scale */
1041 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1042 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1043 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1044 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1045 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1048 /* HitAndMiss Kernels */
1055 if (args->rho < args->sigma)
1057 kernel->width = ((unsigned long)args->sigma)*2+1;
1058 limit1 = (long)args->rho*args->rho;
1059 limit2 = (long)args->sigma*args->sigma;
1063 kernel->width = ((unsigned long)args->rho)*2+1;
1064 limit1 = (long)args->sigma*args->sigma;
1065 limit2 = (long)args->rho*args->rho;
1067 if ( limit2 <= 0 ) /* default outer radius approx 3.5 */
1068 kernel->width = 7L, limit2 = 11L;
1069 kernel->height = kernel->width;
1070 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1071 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1072 kernel->height*sizeof(double));
1073 if (kernel->values == (double *) NULL)
1074 return(DestroyKernelInfo(kernel));
1076 /* set a ring of background points */
1077 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
1078 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1079 { long radius=u*u+v*v;
1080 if ( limit1 <= radius && radius <= limit2)
1081 kernel->values[i] = 0.0;
1083 kernel->values[i] = nan;
1085 /* central point is always foreground */
1086 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1087 kernel->positive_range = 1.0;
1088 kernel->maximum = 1.0;
1093 kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
1094 kernel=ParseKernelArray("3x3: 0,0,- 0,1,1 -,1,-");
1095 if (kernel == (KernelInfo *) NULL)
1097 kernel->type = type;
1098 ExpandKernelInfo(kernel, 90.0); /* Create a list of 4 rotated kernels */
1101 case LineEndsKernel:
1103 kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
1104 kernel=ParseKernelArray("3x3: 0,-,- 0,1,0 0,0,0");
1105 if (kernel == (KernelInfo *) NULL)
1107 kernel->type = type;
1108 ExpandKernelInfo(kernel, 45.0); /* Create a list of 8 rotated kernels */
1111 case LineJunctionsKernel:
1115 kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
1116 /* first set of 4 kernels */
1117 kernel=ParseKernelArray("3x3: -,1,- -,1,- 1,-,1");
1118 if (kernel == (KernelInfo *) NULL)
1120 kernel->type = type;
1121 ExpandKernelInfo(kernel, 90.0);
1122 /* append second set of 4 kernels */
1123 new_kernel=ParseKernelArray("3x3: 1,-,- -,1,- 1,-,1");
1124 if (new_kernel == (KernelInfo *) NULL)
1125 return(DestroyKernelInfo(kernel));
1126 kernel->type = type;
1127 ExpandKernelInfo(new_kernel, 90.0);
1128 LastKernelInfo(kernel)->next = new_kernel;
1129 /* append Thrid set of 4 kernels */
1130 new_kernel=ParseKernelArray("3x3: -,1,- -,1,1 1,-,-");
1131 if (new_kernel == (KernelInfo *) NULL)
1132 return(DestroyKernelInfo(kernel));
1133 kernel->type = type;
1134 ExpandKernelInfo(new_kernel, 90.0);
1135 LastKernelInfo(kernel)->next = new_kernel;
1138 case ConvexHullKernel:
1142 kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
1143 /* first set of 4 kernels */
1144 kernel=ParseKernelArray("3x3: 1,1,- 1,0,- 1,-,0");
1145 if (kernel == (KernelInfo *) NULL)
1147 kernel->type = type;
1148 ExpandKernelInfo(kernel, 90.0);
1149 /* append second set of 4 kernels */
1150 new_kernel=ParseKernelArray("3x3: -,1,1 -,0,1 0,-,1");
1151 if (new_kernel == (KernelInfo *) NULL)
1152 return(DestroyKernelInfo(kernel));
1153 kernel->type = type;
1154 ExpandKernelInfo(new_kernel, 90.0);
1155 LastKernelInfo(kernel)->next = new_kernel;
1158 case SkeletonKernel:
1162 kernel=DestroyKernelInfo(kernel); /* default kernel is not needed */
1163 /* first set of 4 kernels - corners */
1164 kernel=ParseKernelArray("3x3: 0,0,- 0,1,1 -,1,-");
1165 if (kernel == (KernelInfo *) NULL)
1167 kernel->type = type;
1168 ExpandKernelInfo(kernel, 90);
1169 /* append second set of 4 kernels - edge middles */
1170 new_kernel=ParseKernelArray("3x3: 0,0,0 -,1,- 1,1,1");
1171 if (new_kernel == (KernelInfo *) NULL)
1172 return(DestroyKernelInfo(kernel));
1173 kernel->type = type;
1174 ExpandKernelInfo(new_kernel, 90);
1175 LastKernelInfo(kernel)->next = new_kernel;
1178 /* Distance Measuring Kernels */
1179 case ChebyshevKernel:
1181 if (args->rho < 1.0)
1182 kernel->width = kernel->height = 3; /* default radius = 1 */
1184 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
1185 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1187 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1188 kernel->height*sizeof(double));
1189 if (kernel->values == (double *) NULL)
1190 return(DestroyKernelInfo(kernel));
1192 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1193 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1194 kernel->positive_range += ( kernel->values[i] =
1195 args->sigma*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
1196 kernel->maximum = kernel->values[0];
1199 case ManhattenKernel:
1201 if (args->rho < 1.0)
1202 kernel->width = kernel->height = 3; /* default radius = 1 */
1204 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
1205 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1207 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1208 kernel->height*sizeof(double));
1209 if (kernel->values == (double *) NULL)
1210 return(DestroyKernelInfo(kernel));
1212 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1213 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1214 kernel->positive_range += ( kernel->values[i] =
1215 args->sigma*(labs(u)+labs(v)) );
1216 kernel->maximum = kernel->values[0];
1219 case EuclideanKernel:
1221 if (args->rho < 1.0)
1222 kernel->width = kernel->height = 3; /* default radius = 1 */
1224 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
1225 kernel->x = kernel->y = (long) (kernel->width-1)/2;
1227 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1228 kernel->height*sizeof(double));
1229 if (kernel->values == (double *) NULL)
1230 return(DestroyKernelInfo(kernel));
1232 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
1233 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1234 kernel->positive_range += ( kernel->values[i] =
1235 args->sigma*sqrt((double)(u*u+v*v)) );
1236 kernel->maximum = kernel->values[0];
1239 /* Undefined Kernels */
1241 perror("Kernel Type has not been defined yet\n");
1244 /* Generate a No-Op minimal kernel - 1x1 pixel */
1245 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
1246 if (kernel->values == (double *) NULL)
1247 return(DestroyKernelInfo(kernel));
1248 kernel->width = kernel->height = 1;
1249 kernel->x = kernel->x = 0;
1250 kernel->type = UndefinedKernel;
1252 kernel->positive_range =
1253 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
1261 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1265 % C l o n e K e r n e l I n f o %
1269 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1271 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
1272 % can be modified without effecting the original. The cloned kernel should
1273 % be destroyed using DestoryKernelInfo() when no longer needed.
1275 % The format of the CloneKernelInfo method is:
1277 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1279 % A description of each parameter follows:
1281 % o kernel: the Morphology/Convolution kernel to be cloned
1284 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1292 assert(kernel != (KernelInfo *) NULL);
1293 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1294 if (new_kernel == (KernelInfo *) NULL)
1296 *new_kernel=(*kernel); /* copy values in structure */
1298 /* replace the values with a copy of the values */
1299 new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1300 kernel->height*sizeof(double));
1301 if (new_kernel->values == (double *) NULL)
1302 return(DestroyKernelInfo(new_kernel));
1303 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1304 new_kernel->values[i]=kernel->values[i];
1306 /* Also clone the next kernel in the kernel list */
1307 if ( kernel->next != (KernelInfo *) NULL ) {
1308 new_kernel->next = CloneKernelInfo(kernel->next);
1309 if ( new_kernel->next == (KernelInfo *) NULL )
1310 return(DestroyKernelInfo(new_kernel));
1317 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1321 % D e s t r o y K e r n e l I n f o %
1325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1327 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
1330 % The format of the DestroyKernelInfo method is:
1332 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1334 % A description of each parameter follows:
1336 % o kernel: the Morphology/Convolution kernel to be destroyed
1340 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1342 assert(kernel != (KernelInfo *) NULL);
1344 if ( kernel->next != (KernelInfo *) NULL )
1345 kernel->next = DestroyKernelInfo(kernel->next);
1347 kernel->values = (double *)RelinquishMagickMemory(kernel->values);
1348 kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
1353 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1357 % E x p a n d K e r n e l I n f o %
1361 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1363 % ExpandKernelInfo() takes a single kernel, and expands it into a list
1364 % of kernels each incrementally rotated the angle given.
1366 % WARNING: 45 degree rotations only works for 3x3 kernels.
1367 % While 90 degree roatations only works for linear and square kernels
1369 % The format of the RotateKernelInfo method is:
1371 % void ExpandKernelInfo(KernelInfo *kernel, double angle)
1373 % A description of each parameter follows:
1375 % o kernel: the Morphology/Convolution kernel
1377 % o angle: angle to rotate in degrees
1379 % This function is only internel to this module, as it is not finalized,
1380 % especially with regard to non-orthogonal angles, and rotation of larger
1383 static void ExpandKernelInfo(KernelInfo *kernel, double angle)
1393 for (a=angle; a<355.0; a+=angle) {
1394 clone = CloneKernelInfo(last);
1395 RotateKernelInfo(clone, angle);
1402 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1406 % M o r p h o l o g y I m a g e C h a n n e l %
1410 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1412 % MorphologyImageChannel() applies a user supplied kernel to the image
1413 % according to the given mophology method.
1415 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1416 % by the kernel generator.
1418 % The format of the MorphologyImage method is:
1420 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
1421 % const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
1422 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1423 % channel,MorphologyMethod method,const long iterations,
1424 % KernelInfo *kernel,ExceptionInfo *exception)
1426 % A description of each parameter follows:
1428 % o image: the image.
1430 % o method: the morphology method to be applied.
1432 % o iterations: apply the operation this many times (or no change).
1433 % A value of -1 means loop until no change found.
1434 % How this is applied may depend on the morphology method.
1435 % Typically this is a value of 1.
1437 % o channel: the channel type.
1439 % o kernel: An array of double representing the morphology kernel.
1440 % Warning: kernel may be normalized for the Convolve method.
1442 % o exception: return any errors or warnings in this structure.
1445 % TODO: bias and auto-scale handling of the kernel for convolution
1446 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1447 % by the kernel generator.
1452 /* Internal function
1453 * Apply the low-level Morphology Method Primatives using the given Kernel
1454 * Returning the number of pixels that changed.
1455 * Two pre-created images must be provided, no image is created.
1457 static unsigned long MorphologyPrimative(const Image *image, Image
1458 *result_image, const MorphologyMethod method, const ChannelType channel,
1459 const KernelInfo *kernel, ExceptionInfo *exception)
1461 #define MorphologyTag "Morphology/Image"
1478 /* Only the most basic morphology is actually performed by this routine */
1481 Apply Basic Morphology to Image.
1487 GetMagickPixelPacket(image,&bias);
1488 SetMagickPixelPacketBias(image,&bias);
1489 /* Future: handle auto-bias from user, based on kernel input */
1491 p_view=AcquireCacheView(image);
1492 q_view=AcquireCacheView(result_image);
1494 /* Some methods (including convolve) needs use a reflected kernel.
1495 * Adjust 'origin' offsets for this reflected kernel.
1500 case ConvolveMorphology:
1501 case DilateMorphology:
1502 case DilateIntensityMorphology:
1503 case DistanceMorphology:
1504 /* kernel needs to used with reflection about origin */
1505 offx = (long) kernel->width-offx-1;
1506 offy = (long) kernel->height-offy-1;
1508 case ErodeMorphology:
1509 case ErodeIntensityMorphology:
1510 case HitAndMissMorphology:
1511 case ThinningMorphology:
1512 case ThickenMorphology:
1513 /* kernel is user as is, without reflection */
1516 perror("Not a low level Morphology Method");
1520 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1521 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1523 for (y=0; y < (long) image->rows; y++)
1528 register const PixelPacket
1531 register const IndexPacket
1532 *restrict p_indexes;
1534 register PixelPacket
1537 register IndexPacket
1538 *restrict q_indexes;
1546 if (status == MagickFalse)
1548 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1549 image->columns+kernel->width, kernel->height, exception);
1550 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1552 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1557 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1558 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1559 r = (image->columns+kernel->width)*offy+offx; /* constant */
1561 for (x=0; x < (long) image->columns; x++)
1569 register const double
1572 register const PixelPacket
1575 register const IndexPacket
1576 *restrict k_indexes;
1583 /* Copy input to ouput image for unused channels
1584 * This removes need for 'cloning' a new image every iteration
1587 if (image->colorspace == CMYKColorspace)
1588 q_indexes[x] = p_indexes[r];
1595 min.index = (MagickRealType) QuantumRange;
1600 max.index = (MagickRealType) 0;
1601 /* original pixel value */
1602 result.red = (MagickRealType) p[r].red;
1603 result.green = (MagickRealType) p[r].green;
1604 result.blue = (MagickRealType) p[r].blue;
1605 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1607 if ( image->colorspace == CMYKColorspace)
1608 result.index = (MagickRealType) p_indexes[r];
1611 case ConvolveMorphology:
1612 /* Set the user defined bias of the weighted average output
1614 ** FUTURE: provide some way for internal functions to disable
1615 ** user provided bias and scaling effects.
1619 case DilateIntensityMorphology:
1620 case ErodeIntensityMorphology:
1621 result.red = 0.0; /* flag indicating when first match found */
1628 case ConvolveMorphology:
1629 /* Weighted Average of pixels using reflected kernel
1631 ** NOTE for correct working of this operation for asymetrical
1632 ** kernels, the kernel needs to be applied in its reflected form.
1633 ** That is its values needs to be reversed.
1635 ** Correlation is actually the same as this but without reflecting
1636 ** the kernel, and thus 'lower-level' that Convolution. However
1637 ** as Convolution is the more common method used, and it does not
1638 ** really cost us much in terms of processing to use a reflected
1639 ** kernel, so it is Convolution that is implemented.
1641 ** Correlation will have its kernel reflected before calling
1642 ** this function to do a Convolve.
1644 ** For more details of Correlation vs Convolution see
1645 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
1647 if (((channel & SyncChannels) != 0 ) &&
1648 (image->matte == MagickTrue))
1649 { /* Channel has a 'Sync' Flag, and Alpha Channel enabled.
1650 ** Weight the color channels with Alpha Channel so that
1651 ** transparent pixels are not part of the results.
1654 alpha, /* color channel weighting : kernel*alpha */
1655 gamma; /* divisor, sum of weighting values */
1658 k = &kernel->values[ kernel->width*kernel->height-1 ];
1660 k_indexes = p_indexes;
1661 for (v=0; v < (long) kernel->height; v++) {
1662 for (u=0; u < (long) kernel->width; u++, k--) {
1663 if ( IsNan(*k) ) continue;
1664 alpha=(*k)*(QuantumScale*(QuantumRange-
1665 k_pixels[u].opacity));
1667 result.red += alpha*k_pixels[u].red;
1668 result.green += alpha*k_pixels[u].green;
1669 result.blue += alpha*k_pixels[u].blue;
1670 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1671 if ( image->colorspace == CMYKColorspace)
1672 result.index += alpha*k_indexes[u];
1674 k_pixels += image->columns+kernel->width;
1675 k_indexes += image->columns+kernel->width;
1677 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1678 result.red *= gamma;
1679 result.green *= gamma;
1680 result.blue *= gamma;
1681 result.opacity *= gamma;
1682 result.index *= gamma;
1686 /* No 'Sync' flag, or no Alpha involved.
1687 ** Convolution is simple individual channel weigthed sum.
1689 k = &kernel->values[ kernel->width*kernel->height-1 ];
1691 k_indexes = p_indexes;
1692 for (v=0; v < (long) kernel->height; v++) {
1693 for (u=0; u < (long) kernel->width; u++, k--) {
1694 if ( IsNan(*k) ) continue;
1695 result.red += (*k)*k_pixels[u].red;
1696 result.green += (*k)*k_pixels[u].green;
1697 result.blue += (*k)*k_pixels[u].blue;
1698 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1699 if ( image->colorspace == CMYKColorspace)
1700 result.index += (*k)*k_indexes[u];
1702 k_pixels += image->columns+kernel->width;
1703 k_indexes += image->columns+kernel->width;
1708 case ErodeMorphology:
1709 /* Minimum Value within kernel neighbourhood
1711 ** NOTE that the kernel is not reflected for this operation!
1713 ** NOTE: in normal Greyscale Morphology, the kernel value should
1714 ** be added to the real value, this is currently not done, due to
1715 ** the nature of the boolean kernels being used.
1719 k_indexes = p_indexes;
1720 for (v=0; v < (long) kernel->height; v++) {
1721 for (u=0; u < (long) kernel->width; u++, k++) {
1722 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1723 Minimize(min.red, (double) k_pixels[u].red);
1724 Minimize(min.green, (double) k_pixels[u].green);
1725 Minimize(min.blue, (double) k_pixels[u].blue);
1726 Minimize(min.opacity,
1727 QuantumRange-(double) k_pixels[u].opacity);
1728 if ( image->colorspace == CMYKColorspace)
1729 Minimize(min.index, (double) k_indexes[u]);
1731 k_pixels += image->columns+kernel->width;
1732 k_indexes += image->columns+kernel->width;
1737 case DilateMorphology:
1738 /* Maximum Value within kernel neighbourhood
1740 ** NOTE for correct working of this operation for asymetrical
1741 ** kernels, the kernel needs to be applied in its reflected form.
1742 ** That is its values needs to be reversed.
1744 ** NOTE: in normal Greyscale Morphology, the kernel value should
1745 ** be added to the real value, this is currently not done, due to
1746 ** the nature of the boolean kernels being used.
1749 k = &kernel->values[ kernel->width*kernel->height-1 ];
1751 k_indexes = p_indexes;
1752 for (v=0; v < (long) kernel->height; v++) {
1753 for (u=0; u < (long) kernel->width; u++, k--) {
1754 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1755 Maximize(max.red, (double) k_pixels[u].red);
1756 Maximize(max.green, (double) k_pixels[u].green);
1757 Maximize(max.blue, (double) k_pixels[u].blue);
1758 Maximize(max.opacity,
1759 QuantumRange-(double) k_pixels[u].opacity);
1760 if ( image->colorspace == CMYKColorspace)
1761 Maximize(max.index, (double) k_indexes[u]);
1763 k_pixels += image->columns+kernel->width;
1764 k_indexes += image->columns+kernel->width;
1768 case HitAndMissMorphology:
1769 case ThinningMorphology:
1770 case ThickenMorphology:
1771 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
1773 ** NOTE that the kernel is not reflected for this operation,
1774 ** and consists of both foreground and background pixel
1775 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
1776 ** with either Nan or 0.5 values for don't care.
1778 ** Note that this can produce negative results, though really
1779 ** only a positive match has any real value.
1783 k_indexes = p_indexes;
1784 for (v=0; v < (long) kernel->height; v++) {
1785 for (u=0; u < (long) kernel->width; u++, k++) {
1786 if ( IsNan(*k) ) continue;
1788 { /* minimim of foreground pixels */
1789 Minimize(min.red, (double) k_pixels[u].red);
1790 Minimize(min.green, (double) k_pixels[u].green);
1791 Minimize(min.blue, (double) k_pixels[u].blue);
1792 Minimize(min.opacity,
1793 QuantumRange-(double) k_pixels[u].opacity);
1794 if ( image->colorspace == CMYKColorspace)
1795 Minimize(min.index, (double) k_indexes[u]);
1797 else if ( (*k) < 0.3 )
1798 { /* maximum of background pixels */
1799 Maximize(max.red, (double) k_pixels[u].red);
1800 Maximize(max.green, (double) k_pixels[u].green);
1801 Maximize(max.blue, (double) k_pixels[u].blue);
1802 Maximize(max.opacity,
1803 QuantumRange-(double) k_pixels[u].opacity);
1804 if ( image->colorspace == CMYKColorspace)
1805 Maximize(max.index, (double) k_indexes[u]);
1808 k_pixels += image->columns+kernel->width;
1809 k_indexes += image->columns+kernel->width;
1811 /* Pattern Match only if min fg larger than min bg pixels */
1812 min.red -= max.red; Maximize( min.red, 0.0 );
1813 min.green -= max.green; Maximize( min.green, 0.0 );
1814 min.blue -= max.blue; Maximize( min.blue, 0.0 );
1815 min.opacity -= max.opacity; Maximize( min.opacity, 0.0 );
1816 min.index -= max.index; Maximize( min.index, 0.0 );
1819 case ErodeIntensityMorphology:
1820 /* Select Pixel with Minimum Intensity within kernel neighbourhood
1822 ** WARNING: the intensity test fails for CMYK and does not
1823 ** take into account the moderating effect of teh alpha channel
1824 ** on the intensity.
1826 ** NOTE that the kernel is not reflected for this operation!
1830 k_indexes = p_indexes;
1831 for (v=0; v < (long) kernel->height; v++) {
1832 for (u=0; u < (long) kernel->width; u++, k++) {
1833 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1834 if ( result.red == 0.0 ||
1835 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1836 /* copy the whole pixel - no channel selection */
1838 if ( result.red > 0.0 ) changed++;
1842 k_pixels += image->columns+kernel->width;
1843 k_indexes += image->columns+kernel->width;
1847 case DilateIntensityMorphology:
1848 /* Select Pixel with Maximum Intensity within kernel neighbourhood
1850 ** WARNING: the intensity test fails for CMYK and does not
1851 ** take into account the moderating effect of teh alpha channel
1852 ** on the intensity.
1854 ** NOTE for correct working of this operation for asymetrical
1855 ** kernels, the kernel needs to be applied in its reflected form.
1856 ** That is its values needs to be reversed.
1858 k = &kernel->values[ kernel->width*kernel->height-1 ];
1860 k_indexes = p_indexes;
1861 for (v=0; v < (long) kernel->height; v++) {
1862 for (u=0; u < (long) kernel->width; u++, k--) {
1863 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1864 if ( result.red == 0.0 ||
1865 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1866 /* copy the whole pixel - no channel selection */
1868 if ( result.red > 0.0 ) changed++;
1872 k_pixels += image->columns+kernel->width;
1873 k_indexes += image->columns+kernel->width;
1878 case DistanceMorphology:
1879 /* Add kernel Value and select the minimum value found.
1880 ** The result is a iterative distance from edge of image shape.
1882 ** All Distance Kernels are symetrical, but that may not always
1883 ** be the case. For example how about a distance from left edges?
1884 ** To work correctly with asymetrical kernels the reflected kernel
1885 ** needs to be applied.
1887 ** Actually this is really a GreyErode with a negative kernel!
1890 k = &kernel->values[ kernel->width*kernel->height-1 ];
1892 k_indexes = p_indexes;
1893 for (v=0; v < (long) kernel->height; v++) {
1894 for (u=0; u < (long) kernel->width; u++, k--) {
1895 if ( IsNan(*k) ) continue;
1896 Minimize(result.red, (*k)+k_pixels[u].red);
1897 Minimize(result.green, (*k)+k_pixels[u].green);
1898 Minimize(result.blue, (*k)+k_pixels[u].blue);
1899 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1900 if ( image->colorspace == CMYKColorspace)
1901 Minimize(result.index, (*k)+k_indexes[u]);
1903 k_pixels += image->columns+kernel->width;
1904 k_indexes += image->columns+kernel->width;
1908 case UndefinedMorphology:
1910 break; /* Do nothing */
1912 /* Final mathematics of results (combine with original image?)
1914 ** NOTE: Difference Morphology operators Edge* and *Hat could also
1915 ** be done here but works better with iteration as a image difference
1916 ** in the controling function (below). Thicken and Thinning however
1917 ** should be done here so thay can be iterated correctly.
1920 case HitAndMissMorphology:
1921 case ErodeMorphology:
1922 result = min; /* minimum of neighbourhood */
1924 case DilateMorphology:
1925 result = max; /* maximum of neighbourhood */
1927 case ThinningMorphology:
1928 /* subtract pattern match from original */
1929 result.red -= min.red;
1930 result.green -= min.green;
1931 result.blue -= min.blue;
1932 result.opacity -= min.opacity;
1933 result.index -= min.index;
1935 case ThickenMorphology:
1936 /* Union with original image (maximize) - or should this be + */
1937 Maximize( result.red, min.red );
1938 Maximize( result.green, min.green );
1939 Maximize( result.blue, min.blue );
1940 Maximize( result.opacity, min.opacity );
1941 Maximize( result.index, min.index );
1944 /* result directly calculated or assigned */
1947 /* Assign the resulting pixel values - Clamping Result */
1949 case UndefinedMorphology:
1950 case DilateIntensityMorphology:
1951 case ErodeIntensityMorphology:
1952 break; /* full pixel was directly assigned - not a channel method */
1954 if ((channel & RedChannel) != 0)
1955 q->red = ClampToQuantum(result.red);
1956 if ((channel & GreenChannel) != 0)
1957 q->green = ClampToQuantum(result.green);
1958 if ((channel & BlueChannel) != 0)
1959 q->blue = ClampToQuantum(result.blue);
1960 if ((channel & OpacityChannel) != 0
1961 && image->matte == MagickTrue )
1962 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1963 if ((channel & IndexChannel) != 0
1964 && image->colorspace == CMYKColorspace)
1965 q_indexes[x] = ClampToQuantum(result.index);
1968 /* Count up changed pixels */
1969 if ( ( p[r].red != q->red )
1970 || ( p[r].green != q->green )
1971 || ( p[r].blue != q->blue )
1972 || ( p[r].opacity != q->opacity )
1973 || ( image->colorspace == CMYKColorspace &&
1974 p_indexes[r] != q_indexes[x] ) )
1975 changed++; /* The pixel had some value changed! */
1979 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1980 if (sync == MagickFalse)
1982 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1987 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1988 #pragma omp critical (MagickCore_MorphologyImage)
1990 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1991 if (proceed == MagickFalse)
1995 result_image->type=image->type;
1996 q_view=DestroyCacheView(q_view);
1997 p_view=DestroyCacheView(p_view);
1998 return(status ? (unsigned long) changed : 0);
2002 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
2003 method, const long iterations,const KernelInfo *kernel, ExceptionInfo
2009 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
2010 iterations,kernel,exception);
2011 return(morphology_image);
2015 MagickExport Image *MorphologyImageChannel(const Image *image,
2016 const ChannelType channel,const MorphologyMethod method,
2017 const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
2038 assert(image != (Image *) NULL);
2039 assert(image->signature == MagickSignature);
2040 assert(kernel != (KernelInfo *) NULL);
2041 assert(kernel->signature == MagickSignature);
2042 assert(exception != (ExceptionInfo *) NULL);
2043 assert(exception->signature == MagickSignature);
2045 if ( iterations == 0 )
2046 return((Image *)NULL); /* null operation - nothing to do! */
2048 /* kernel must be valid at this point
2049 * (except maybe for posible future morphology methods like "Prune"
2051 assert(kernel != (KernelInfo *)NULL);
2053 new_image = (Image *) NULL; /* for cleanup */
2054 old_image = (Image *) NULL;
2055 grad_image = (Image *) NULL;
2056 curr_kernel = (KernelInfo *) NULL;
2057 count = 0; /* interation count */
2058 changed = 1; /* was last run succesfull */
2059 curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
2060 curr_method = method; /* to be changed as nessary */
2062 limit = (unsigned long) iterations;
2063 if ( iterations < 0 )
2064 limit = image->columns > image->rows ? image->columns : image->rows;
2066 /* Third-level morphology methods */
2067 switch( curr_method ) {
2068 case EdgeMorphology:
2069 grad_image = MorphologyImageChannel(image, channel,
2070 DilateMorphology, iterations, curr_kernel, exception);
2071 if ( grad_image == (Image *) NULL )
2074 case EdgeInMorphology:
2075 curr_method = ErodeMorphology;
2077 case EdgeOutMorphology:
2078 curr_method = DilateMorphology;
2080 case TopHatMorphology:
2081 curr_method = OpenMorphology;
2083 case BottomHatMorphology:
2084 curr_method = CloseMorphology;
2087 break; /* not a third-level method */
2090 /* Second-level morphology methods */
2091 switch( curr_method ) {
2092 case OpenMorphology:
2093 /* Open is a Erode then a Dilate without reflection */
2094 new_image = MorphologyImageChannel(image, channel,
2095 ErodeMorphology, iterations, curr_kernel, exception);
2096 if (new_image == (Image *) NULL)
2098 curr_method = DilateMorphology;
2100 case OpenIntensityMorphology:
2101 new_image = MorphologyImageChannel(image, channel,
2102 ErodeIntensityMorphology, iterations, curr_kernel, exception);
2103 if (new_image == (Image *) NULL)
2105 curr_method = DilateIntensityMorphology;
2108 case CloseMorphology:
2109 /* Close is a Dilate then Erode using reflected kernel */
2110 /* A reflected kernel is needed for a Close */
2111 if ( curr_kernel == kernel )
2112 curr_kernel = CloneKernelInfo(kernel);
2113 if (curr_kernel == (KernelInfo *) NULL)
2115 RotateKernelInfo(curr_kernel,180);
2116 new_image = MorphologyImageChannel(image, channel,
2117 DilateMorphology, iterations, curr_kernel, exception);
2118 if (new_image == (Image *) NULL)
2120 curr_method = ErodeMorphology;
2122 case CloseIntensityMorphology:
2123 /* A reflected kernel is needed for a Close */
2124 if ( curr_kernel == kernel )
2125 curr_kernel = CloneKernelInfo(kernel);
2126 if (curr_kernel == (KernelInfo *) NULL)
2128 RotateKernelInfo(curr_kernel,180);
2129 new_image = MorphologyImageChannel(image, channel,
2130 DilateIntensityMorphology, iterations, curr_kernel, exception);
2131 if (new_image == (Image *) NULL)
2133 curr_method = ErodeIntensityMorphology;
2136 case CorrelateMorphology:
2137 /* A Correlation is actually a Convolution with a reflected kernel.
2138 ** However a Convolution is a weighted sum with a reflected kernel.
2139 ** It may seem stange to convert a Correlation into a Convolution
2140 ** as the Correleation is the simplier method, but Convolution is
2141 ** much more commonly used, and it makes sense to implement it directly
2142 ** so as to avoid the need to duplicate the kernel when it is not
2143 ** required (which is typically the default).
2145 if ( curr_kernel == kernel )
2146 curr_kernel = CloneKernelInfo(kernel);
2147 if (curr_kernel == (KernelInfo *) NULL)
2149 RotateKernelInfo(curr_kernel,180);
2150 curr_method = ConvolveMorphology;
2151 /* FALL-THRU into Correlate (weigthed sum without reflection) */
2153 case ConvolveMorphology:
2154 { /* Scale or Normalize kernel, according to user wishes
2155 ** before using it for the Convolve/Correlate method.
2157 ** FUTURE: provide some way for internal functions to disable
2158 ** user bias and scaling effects.
2161 *artifact = GetImageArtifact(image,"convolve:scale");
2162 if ( artifact != (char *)NULL ) {
2168 if ( curr_kernel == kernel )
2169 curr_kernel = CloneKernelInfo(kernel);
2170 if (curr_kernel == (KernelInfo *) NULL)
2173 flags = (GeometryFlags) ParseGeometry(artifact, &args);
2174 ScaleKernelInfo(curr_kernel, args.rho, flags);
2177 /* FALL-THRU to do the first, and typically the only iteration */
2180 /* Do a single iteration using the Low-Level Morphology method!
2181 ** This ensures a "new_image" has been generated, but allows us to skip
2182 ** the creation of 'old_image' if no more iterations are needed.
2184 ** The "curr_method" should also be set to a low-level method that is
2185 ** understood by the MorphologyPrimative() internal function.
2187 new_image=CloneImage(image,0,0,MagickTrue,exception);
2188 if (new_image == (Image *) NULL)
2190 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
2192 InheritException(exception,&new_image->exception);
2195 count++; /* iteration count */
2196 changed = MorphologyPrimative(image,new_image,curr_method,channel,
2197 curr_kernel, exception);
2198 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2199 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
2200 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2201 count, 0L, changed);
2205 * + "curr_method" should be set to a low-level morphology method
2206 * + "count=1" if the first iteration of the first kernel has been done.
2207 * + "new_image" is defined and contains the current resulting image
2210 /* The Hit-And-Miss Morphology is not 'iterated' against the resulting
2211 ** image from the previous morphology step. It must always be applied
2212 ** to the original image, with the results collected into a union (maximimum
2213 ** or lighten) of all the results. Multiple kernels may be applied but
2214 ** an iteration of the morphology does nothing, so is ignored.
2216 ** The first kernel is guranteed to have been done, so lets continue the
2217 ** process, with the other kernels in the kernel list.
2219 if ( method == HitAndMissMorphology )
2221 if ( curr_kernel->next != (KernelInfo *) NULL ) {
2222 /* create a second working image */
2223 old_image = CloneImage(image,0,0,MagickTrue,exception);
2224 if (old_image == (Image *) NULL)
2226 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
2228 InheritException(exception,&old_image->exception);
2232 /* loop through rest of the kernels */
2233 this_kernel=curr_kernel->next;
2235 while( this_kernel != (KernelInfo *) NULL )
2237 changed = MorphologyPrimative(image,old_image,curr_method,channel,
2238 this_kernel,exception);
2239 (void) CompositeImageChannel(new_image,
2240 (ChannelType) (channel & ~SyncChannels), LightenCompositeOp,
2242 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2243 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
2244 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2245 count, kernel_number, changed);
2246 this_kernel = this_kernel->next;
2249 old_image=DestroyImage(old_image);
2254 /* Repeat the low-level morphology over all kernels
2255 until iteration count limit or no change from any kernel is found */
2256 if ( ( count < limit && changed > 0 ) ||
2257 curr_kernel->next != (KernelInfo *) NULL ) {
2259 /* create a second working image */
2260 old_image = CloneImage(image,0,0,MagickTrue,exception);
2261 if (old_image == (Image *) NULL)
2263 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
2265 InheritException(exception,&old_image->exception);
2269 /* reset variables for the first/next iteration, or next kernel) */
2271 this_kernel = curr_kernel;
2272 total_changed = count != 0 ? changed : 0;
2273 if ( count != 0 && this_kernel != (KernelInfo *) NULL ) {
2274 count = 0; /* first iteration is not yet finished! */
2275 this_kernel = curr_kernel->next;
2277 total_changed = changed;
2280 while ( count < limit ) {
2282 while ( this_kernel != (KernelInfo *) NULL ) {
2283 Image *tmp = old_image;
2284 old_image = new_image;
2286 changed = MorphologyPrimative(old_image,new_image,curr_method,channel,
2287 this_kernel,exception);
2288 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2289 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
2290 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2291 count, kernel_number, changed);
2292 total_changed += changed;
2293 this_kernel = this_kernel->next;
2296 if ( kernel_number > 1 )
2297 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2298 fprintf(stderr, "Morphology %s:%lu ===> Total Changed %lu\n",
2299 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2300 count, total_changed);
2301 if ( total_changed == 0 )
2302 break; /* no changes after processing all kernels - ABORT */
2303 /* prepare for next loop */
2306 this_kernel = curr_kernel;
2308 old_image=DestroyImage(old_image);
2311 /* finished with kernel - destroy any copy that was made */
2312 if ( curr_kernel != kernel )
2313 curr_kernel=DestroyKernelInfo(curr_kernel);
2315 /* Third-level Subtractive methods post-processing
2317 ** The removal of any 'Sync' channel flag in the Image Compositon below
2318 ** ensures the compose method is applied in a purely mathematical way, only
2319 ** the selected channels, without any normal 'alpha blending' normally
2320 ** associated with the compose method.
2322 ** Note "method" here is the 'original' morphological method, and not the
2323 ** 'current' morphological method used above to generate "new_image".
2326 case EdgeOutMorphology:
2327 case EdgeInMorphology:
2328 case TopHatMorphology:
2329 case BottomHatMorphology:
2330 /* Get Difference relative to the original image */
2331 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2332 DifferenceCompositeOp, image, 0, 0);
2334 case EdgeMorphology:
2335 /* Difference the Eroded image from the saved Dilated image */
2336 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2337 DifferenceCompositeOp, grad_image, 0, 0);
2338 grad_image=DestroyImage(grad_image);
2346 /* Yes goto's are bad, but in this case it makes cleanup lot more efficient */
2348 if ( new_image != (Image *) NULL )
2349 DestroyImage(new_image);
2351 if ( old_image != (Image *) NULL )
2352 DestroyImage(old_image);
2353 if ( curr_kernel != (KernelInfo *) NULL && curr_kernel != kernel )
2354 curr_kernel=DestroyKernelInfo(curr_kernel);
2359 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2363 + R o t a t e K e r n e l I n f o %
2367 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2369 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
2370 % restricted to 90 degree angles, but this may be improved in the future.
2372 % The format of the RotateKernelInfo method is:
2374 % void RotateKernelInfo(KernelInfo *kernel, double angle)
2376 % A description of each parameter follows:
2378 % o kernel: the Morphology/Convolution kernel
2380 % o angle: angle to rotate in degrees
2382 % This function is only internel to this module, as it is not finalized,
2383 % especially with regard to non-orthogonal angles, and rotation of larger
2386 static void RotateKernelInfo(KernelInfo *kernel, double angle)
2388 /* angle the lower kernels first */
2389 if ( kernel->next != (KernelInfo *) NULL)
2390 RotateKernelInfo(kernel->next, angle);
2392 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
2394 ** TODO: expand beyond simple 90 degree rotates, flips and flops
2397 /* Modulus the angle */
2398 angle = fmod(angle, 360.0);
2402 if ( 337.5 < angle || angle <= 22.5 )
2403 return; /* no change! - At least at this time */
2405 /* Handle special cases */
2406 switch (kernel->type) {
2407 /* These built-in kernels are cylindrical kernels, rotating is useless */
2408 case GaussianKernel:
2412 case LaplacianKernel:
2413 case ChebyshevKernel:
2414 case ManhattenKernel:
2415 case EuclideanKernel:
2418 /* These may be rotatable at non-90 angles in the future */
2419 /* but simply rotating them in multiples of 90 degrees is useless */
2426 /* These only allows a +/-90 degree rotation (by transpose) */
2427 /* A 180 degree rotation is useless */
2429 case RectangleKernel:
2430 if ( 135.0 < angle && angle <= 225.0 )
2432 if ( 225.0 < angle && angle <= 315.0 )
2439 /* Attempt rotations by 45 degrees */
2440 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
2442 if ( kernel->width == 3 && kernel->height == 3 )
2443 { /* Rotate a 3x3 square by 45 degree angle */
2444 MagickRealType t = kernel->values[0];
2445 kernel->values[0] = kernel->values[1];
2446 kernel->values[1] = kernel->values[2];
2447 kernel->values[2] = kernel->values[5];
2448 kernel->values[5] = kernel->values[8];
2449 kernel->values[8] = kernel->values[7];
2450 kernel->values[7] = kernel->values[6];
2451 kernel->values[6] = kernel->values[3];
2452 kernel->values[3] = t;
2453 angle = fmod(angle+45.0, 360.0);
2456 perror("Unable to rotate non-3x3 kernel by 45 degrees");
2458 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
2460 if ( kernel->width == 1 || kernel->height == 1 )
2461 { /* Do a transpose of the image, which results in a 90
2462 ** degree rotation of a 1 dimentional kernel
2466 t = (long) kernel->width;
2467 kernel->width = kernel->height;
2468 kernel->height = (unsigned long) t;
2470 kernel->x = kernel->y;
2472 if ( kernel->width == 1 )
2473 angle = fmod(angle+270.0, 360.0); /* angle -90 degrees */
2475 angle = fmod(angle+90.0, 360.0); /* angle +90 degrees */
2477 else if ( kernel->width == kernel->height )
2478 { /* Rotate a square array of values by 90 degrees */
2479 register unsigned long
2481 register MagickRealType
2484 for( i=0, x=kernel->width-1; i<=x; i++, x--)
2485 for( j=0, y=kernel->height-1; j<y; j++, y--)
2486 { t = k[i+j*kernel->width];
2487 k[i+j*kernel->width] = k[y+i*kernel->width];
2488 k[y+i*kernel->width] = k[x+y*kernel->width];
2489 k[x+y*kernel->width] = k[j+x*kernel->width];
2490 k[j+x*kernel->width] = t;
2492 angle = fmod(angle+90.0, 360.0); /* angle +90 degrees */
2495 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
2497 if ( 135.0 < angle && angle <= 225.0 )
2499 /* For a 180 degree rotation - also know as a reflection */
2500 /* This is actually a very very common operation! */
2501 /* Basically all that is needed is a reversal of the kernel data! */
2508 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
2509 t=k[i], k[i]=k[j], k[j]=t;
2511 kernel->x = (long) kernel->width - kernel->x - 1;
2512 kernel->y = (long) kernel->height - kernel->y - 1;
2513 angle = fmod(angle+180.0, 360.0); /* angle+180 degrees */
2515 /* At this point angle should at least between -45 (315) and +45 degrees
2516 * In the future some form of non-orthogonal angled rotates could be
2517 * performed here, posibily with a linear kernel restriction.
2521 { /* Do a Flop by reversing each row.
2530 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2531 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2532 t=k[x], k[x]=k[r], k[r]=t;
2534 kernel->x = kernel->width - kernel->x - 1;
2535 angle = fmod(angle+180.0, 360.0);
2542 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2546 % S c a l e K e r n e l I n f o %
2550 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2552 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
2553 % without normalization of the sum of the kernel values (as per given flags).
2555 % By default (no flags given) the values within the kernel is scaled
2556 % directly using given scaling factor without change.
2558 % If any 'normalize_flags' are given the kernel will first be normalized and
2559 % then further scaled by the scaling factor value given. A 'PercentValue'
2560 % flag will cause the given scaling factor to be divided by one hundred
2563 % Kernel normalization ('normalize_flags' given) is designed to ensure that
2564 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
2565 % morphology methods will fall into -1.0 to +1.0 range. Note that for
2566 % non-HDRI versions of IM this may cause images to have any negative results
2567 % clipped, unless some 'bias' is used.
2569 % More specifically. Kernels which only contain positive values (such as a
2570 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
2571 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
2573 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
2574 % the kernel will be scaled by the absolute of the sum of kernel values, so
2575 % that it will generally fall within the +/- 1.0 range.
2577 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
2578 % will be scaled by just the sum of the postive values, so that its output
2579 % range will again fall into the +/- 1.0 range.
2581 % For special kernels designed for locating shapes using 'Correlate', (often
2582 % only containing +1 and -1 values, representing foreground/brackground
2583 % matching) a special normalization method is provided to scale the positive
2584 % values seperatally to those of the negative values, so the kernel will be
2585 % forced to become a zero-sum kernel better suited to such searches.
2587 % WARNING: Correct normalization of the kernel assumes that the '*_range'
2588 % attributes within the kernel structure have been correctly set during the
2591 % NOTE: The values used for 'normalize_flags' have been selected specifically
2592 % to match the use of geometry options, so that '!' means NormalizeValue,
2593 % '^' means CorrelateNormalizeValue, and '%' means PercentValue. All other
2594 % GeometryFlags values are ignored.
2596 % The format of the ScaleKernelInfo method is:
2598 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
2599 % const MagickStatusType normalize_flags )
2601 % A description of each parameter follows:
2603 % o kernel: the Morphology/Convolution kernel
2606 % multiply all values (after normalization) by this factor if not
2607 % zero. If the kernel is normalized regardless of any flags.
2609 % o normalize_flags:
2610 % GeometryFlags defining normalization method to use.
2611 % specifically: NormalizeValue, CorrelateNormalizeValue,
2612 % and/or PercentValue
2614 % This function is internal to this module only at this time, but can be
2615 % exported to other modules if needed.
2617 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
2618 const double scaling_factor,const GeometryFlags normalize_flags)
2627 /* scale the lower kernels first */
2628 if ( kernel->next != (KernelInfo *) NULL)
2629 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
2632 if ( (normalize_flags&NormalizeValue) != 0 ) {
2633 /* normalize kernel appropriately */
2634 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
2635 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
2637 pos_scale = kernel->positive_range; /* special zero-summing kernel */
2639 /* force kernel into being a normalized zero-summing kernel */
2640 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
2641 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
2642 ? kernel->positive_range : 1.0;
2643 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
2644 ? -kernel->negative_range : 1.0;
2647 neg_scale = pos_scale;
2649 /* finialize scaling_factor for positive and negative components */
2650 pos_scale = scaling_factor/pos_scale;
2651 neg_scale = scaling_factor/neg_scale;
2652 if ( (normalize_flags&PercentValue) != 0 ) {
2657 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2658 if ( ! IsNan(kernel->values[i]) )
2659 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
2661 /* convolution output range */
2662 kernel->positive_range *= pos_scale;
2663 kernel->negative_range *= neg_scale;
2664 /* maximum and minimum values in kernel */
2665 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
2666 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
2668 /* swap kernel settings if user scaling factor is negative */
2669 if ( scaling_factor < MagickEpsilon ) {
2671 t = kernel->positive_range;
2672 kernel->positive_range = kernel->negative_range;
2673 kernel->negative_range = t;
2674 t = kernel->maximum;
2675 kernel->maximum = kernel->minimum;
2676 kernel->minimum = 1;
2683 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2687 + S h o w K e r n e l I n f o %
2691 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2693 % ShowKernelInfo() outputs the details of the given kernel defination to
2694 % standard error, generally due to a users 'showkernel' option request.
2696 % The format of the ShowKernel method is:
2698 % void ShowKernelInfo(KernelInfo *kernel)
2700 % A description of each parameter follows:
2702 % o kernel: the Morphology/Convolution kernel
2704 % This function is internal to this module only at this time. That may change
2707 MagickExport void ShowKernelInfo(KernelInfo *kernel)
2715 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
2717 fprintf(stderr, "Kernel ");
2718 if ( kernel->next != (KernelInfo *) NULL )
2719 fprintf(stderr, " #%ld", c );
2720 fprintf(stderr, " \"%s\" of size %lux%lu%+ld%+ld ",
2721 MagickOptionToMnemonic(MagickKernelOptions, k->type),
2722 kernel->width, k->height,
2723 kernel->x, kernel->y );
2725 " with values from %.*lg to %.*lg\n",
2726 GetMagickPrecision(), k->minimum,
2727 GetMagickPrecision(), k->maximum);
2728 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2729 GetMagickPrecision(), k->negative_range,
2730 GetMagickPrecision(), k->positive_range,
2731 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2732 for (i=v=0; v < (long) k->height; v++) {
2733 fprintf(stderr,"%2ld:",v);
2734 for (u=0; u < (long) k->width; u++, i++)
2735 if ( IsNan(k->values[i]) )
2736 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2738 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2739 GetMagickPrecision(), k->values[i]);
2740 fprintf(stderr,"\n");
2746 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2750 + Z e r o K e r n e l N a n s %
2754 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2756 % ZeroKernelNans() replaces any special 'nan' value that may be present in
2757 % the kernel with a zero value. This is typically done when the kernel will
2758 % be used in special hardware (GPU) convolution processors, to simply
2761 % The format of the ZeroKernelNans method is:
2763 % voidZeroKernelNans (KernelInfo *kernel)
2765 % A description of each parameter follows:
2767 % o kernel: the Morphology/Convolution kernel
2769 % FUTURE: return the information in a string for API usage.
2771 MagickExport void ZeroKernelNans(KernelInfo *kernel)
2776 /* scale the lower kernels first */
2777 if ( kernel->next != (KernelInfo *) NULL)
2778 ZeroKernelNans(kernel->next);
2780 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2781 if ( IsNan(kernel->values[i]) )
2782 kernel->values[i] = 0.0;