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 RotateKernelInfo(KernelInfo *, double);
113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117 % A c q u i r e K e r n e l I n f o %
121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123 % AcquireKernelInfo() takes the given string (generally supplied by the
124 % user) and converts it into a Morphology/Convolution Kernel. This allows
125 % users to specify a kernel from a number of pre-defined kernels, or to fully
126 % specify their own kernel for a specific Convolution or Morphology
129 % The kernel so generated can be any rectangular array of floating point
130 % values (doubles) with the 'control point' or 'pixel being affected'
131 % anywhere within that array of values.
133 % Previously IM was restricted to a square of odd size using the exact
134 % center as origin, this is no longer the case, and any rectangular kernel
135 % with any value being declared the origin. This in turn allows the use of
136 % highly asymmetrical kernels.
138 % The floating point values in the kernel can also include a special value
139 % known as 'nan' or 'not a number' to indicate that this value is not part
140 % of the kernel array. This allows you to shaped the kernel within its
141 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
142 % shape. However at least one non-nan value must be provided for correct
143 % working of a kernel.
145 % The returned kernel should be freed using the DestroyKernelInfo() when you
146 % are finished with it. Do not free this memory yourself.
148 % Input kernel defintion strings can consist of any of three types.
151 % Select from one of the built in kernels, using the name and
152 % geometry arguments supplied. See AcquireKernelBuiltIn()
154 % "WxH[+X+Y]:num, num, num ..."
155 % a kernel of size W by H, with W*H floating point numbers following.
156 % the 'center' can be optionally be defined at +X+Y (such that +0+0
157 % is top left corner). If not defined the pixel in the center, for
158 % odd sizes, or to the immediate top or left of center for even sizes
159 % is automatically selected.
161 % "num, num, num, num, ..."
162 % list of floating point numbers defining an 'old style' odd sized
163 % square kernel. At least 9 values should be provided for a 3x3
164 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
165 % Values can be space or comma separated. This is not recommended.
167 % Note that 'name' kernels will start with an alphabetic character while the
168 % new kernel specification has a ':' character in its specification string.
169 % If neither is the case, it is assumed an old style of a simple list of
170 % numbers generating a odd-sized square kernel has been given.
172 % You can define a 'list of kernels' which can be used by some morphology
173 % operators A list is defined as a semi-colon seperated list kernels.
175 % " kernel ; kernel ; kernel ; "
177 % Extra ';' characters are simply ignored.
179 % The format of the AcquireKernal method is:
181 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
183 % A description of each parameter follows:
185 % o kernel_string: the Morphology/Convolution kernel wanted.
189 /* This was separated so that it could be used as a separate
190 ** array input handling function, such as for -color-matrix
192 static KernelInfo *ParseKernelArray(const char *kernel_string)
198 token[MaxTextExtent];
208 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
210 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
211 if (kernel == (KernelInfo *)NULL)
213 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
214 kernel->minimum = kernel->maximum = 0.0;
215 kernel->negative_range = kernel->positive_range = 0.0;
216 kernel->type = UserDefinedKernel;
217 kernel->next = (KernelInfo *) NULL;
218 kernel->signature = MagickSignature;
220 /* find end of this specific kernel definition string */
221 end = strchr(kernel_string, ';');
222 if ( end == (char *) NULL )
223 end = strchr(kernel_string, '\0');
225 /* Has a ':' in argument - New user kernel specification */
226 p = strchr(kernel_string, ':');
227 if ( p != (char *) NULL && p < end)
235 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
236 memcpy(token, kernel_string, (size_t) (p-kernel_string));
237 token[p-kernel_string] = '\0';
238 SetGeometryInfo(&args);
239 flags = ParseGeometry(token, &args);
241 /* Size handling and checks of geometry settings */
242 if ( (flags & WidthValue) == 0 ) /* if no width then */
243 args.rho = args.sigma; /* then width = height */
244 if ( args.rho < 1.0 ) /* if width too small */
245 args.rho = 1.0; /* then width = 1 */
246 if ( args.sigma < 1.0 ) /* if height too small */
247 args.sigma = args.rho; /* then height = width */
248 kernel->width = (unsigned long)args.rho;
249 kernel->height = (unsigned long)args.sigma;
251 /* Offset Handling and Checks */
252 if ( args.xi < 0.0 || args.psi < 0.0 )
253 return(DestroyKernelInfo(kernel));
254 kernel->x = ((flags & XValue)!=0) ? (long)args.xi
255 : (long) (kernel->width-1)/2;
256 kernel->y = ((flags & YValue)!=0) ? (long)args.psi
257 : (long) (kernel->height-1)/2;
258 if ( kernel->x >= (long) kernel->width ||
259 kernel->y >= (long) kernel->height )
260 return(DestroyKernelInfo(kernel));
262 p++; /* advance beyond the ':' */
265 { /* ELSE - Old old specification, forming odd-square kernel */
266 /* count up number of values given */
267 p=(const char *) kernel_string;
268 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
269 p++; /* ignore "'" chars for convolve filter usage - Cristy */
270 for (i=0; p < end; i++)
272 GetMagickToken(p,&p,token);
274 GetMagickToken(p,&p,token);
276 /* set the size of the kernel - old sized square */
277 kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
278 kernel->x = kernel->y = (long) (kernel->width-1)/2;
279 p=(const char *) kernel_string;
280 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
281 p++; /* ignore "'" chars for convolve filter usage - Cristy */
284 /* Read in the kernel values from rest of input string argument */
285 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
286 kernel->height*sizeof(double));
287 if (kernel->values == (double *) NULL)
288 return(DestroyKernelInfo(kernel));
290 kernel->minimum = +MagickHuge;
291 kernel->maximum = -MagickHuge;
292 kernel->negative_range = kernel->positive_range = 0.0;
294 for (i=0; (i < (long) (kernel->width*kernel->height)) && (p < end); i++)
296 GetMagickToken(p,&p,token);
298 GetMagickToken(p,&p,token);
299 if ( LocaleCompare("nan",token) == 0
300 || LocaleCompare("-",token) == 0 ) {
301 kernel->values[i] = nan; /* do not include this value in kernel */
304 kernel->values[i] = StringToDouble(token);
305 ( kernel->values[i] < 0)
306 ? ( kernel->negative_range += kernel->values[i] )
307 : ( kernel->positive_range += kernel->values[i] );
308 Minimize(kernel->minimum, kernel->values[i]);
309 Maximize(kernel->maximum, kernel->values[i]);
313 /* sanity check -- no more values in kernel definition */
314 GetMagickToken(p,&p,token);
315 if ( *token != '\0' && *token != ';' && *token != '\'' )
316 return(DestroyKernelInfo(kernel));
319 /* this was the old method of handling a incomplete kernel */
320 if ( i < (long) (kernel->width*kernel->height) ) {
321 Minimize(kernel->minimum, kernel->values[i]);
322 Maximize(kernel->maximum, kernel->values[i]);
323 for ( ; i < (long) (kernel->width*kernel->height); i++)
324 kernel->values[i]=0.0;
327 /* Number of values for kernel was not enough - Report Error */
328 if ( i < (long) (kernel->width*kernel->height) )
329 return(DestroyKernelInfo(kernel));
332 /* check that we recieved at least one real (non-nan) value! */
333 if ( kernel->minimum == MagickHuge )
334 return(DestroyKernelInfo(kernel));
339 static KernelInfo *ParseNamedKernel(const char *kernel_string)
342 token[MaxTextExtent];
357 /* Parse special 'named' kernel */
358 GetMagickToken(kernel_string,&p,token);
359 type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
360 if ( type < 0 || type == UserDefinedKernel )
361 return((KernelInfo *)NULL); /* not a valid named kernel */
363 while (((isspace((int) ((unsigned char) *p)) != 0) ||
364 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
367 end = strchr(p, ';'); /* end of this kernel defintion */
368 if ( end == (char *) NULL )
369 end = strchr(p, '\0');
371 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
372 memcpy(token, p, (size_t) (end-p));
374 SetGeometryInfo(&args);
375 flags = ParseGeometry(token, &args);
377 /* special handling of missing values in input string */
379 case RectangleKernel:
380 if ( (flags & WidthValue) == 0 ) /* if no width then */
381 args.rho = args.sigma; /* then width = height */
382 if ( args.rho < 1.0 ) /* if width too small */
383 args.rho = 3; /* then width = 3 */
384 if ( args.sigma < 1.0 ) /* if height too small */
385 args.sigma = args.rho; /* then height = width */
386 if ( (flags & XValue) == 0 ) /* center offset if not defined */
387 args.xi = (double)(((long)args.rho-1)/2);
388 if ( (flags & YValue) == 0 )
389 args.psi = (double)(((long)args.sigma-1)/2);
395 /* If no scale given (a 0 scale is valid! - set it to 1.0 */
396 if ( (flags & HeightValue) == 0 )
399 case ChebyshevKernel:
400 case ManhattenKernel:
401 case EuclideanKernel:
402 if ( (flags & HeightValue) == 0 )
403 args.sigma = 100.0; /* default distance scaling */
404 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
405 args.sigma = QuantumRange/args.sigma; /* maximum pixel distance */
406 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
407 args.sigma *= QuantumRange/100.0; /* percentage of color range */
413 return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
416 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
425 token[MaxTextExtent];
431 kernel = last_kernel = NULL;
433 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
435 /* ignore multiple ';' following each other */
436 if ( *token != ';' ) {
438 /* tokens starting with alpha is a Named kernel */
439 if (isalpha((int) *token) == 0)
440 new_kernel = ParseKernelArray(p);
441 else /* otherwise a user defined kernel array */
442 new_kernel = ParseNamedKernel(p);
444 if ( last_kernel == (KernelInfo *) NULL ) {
445 /* first kernel: initialize the kernel list */
446 if ( new_kernel == (KernelInfo *) NULL )
447 return((KernelInfo *) NULL);
448 kernel = last_kernel = new_kernel;
451 /* append kernel to end of the list */
452 if ( new_kernel == (KernelInfo *) NULL )
453 return(DestroyKernelInfo(kernel));
454 last_kernel->next = new_kernel;
455 last_kernel = new_kernel;
459 /* look for the next kernel in list */
461 if ( p == (char *) NULL )
471 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
475 % A c q u i r e K e r n e l B u i l t I n %
479 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
481 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
482 % kernels used for special purposes such as gaussian blurring, skeleton
483 % pruning, and edge distance determination.
485 % They take a KernelType, and a set of geometry style arguments, which were
486 % typically decoded from a user supplied string, or from a more complex
487 % Morphology Method that was requested.
489 % The format of the AcquireKernalBuiltIn method is:
491 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
492 % const GeometryInfo args)
494 % A description of each parameter follows:
496 % o type: the pre-defined type of kernel wanted
498 % o args: arguments defining or modifying the kernel
500 % Convolution Kernels
502 % Gaussian "{radius},{sigma}"
503 % Generate a two-dimentional gaussian kernel, as used by -gaussian
504 % A sigma is required, (with the 'x'), due to historical reasons.
506 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
507 % the final size of the resulting kernel to a square 2*radius+1 in size.
508 % The radius should be at least 2 times that of the sigma value, or
509 % sever clipping and aliasing may result. If not given or set to 0 the
510 % radius will be determined so as to produce the best minimal error
511 % result, which is usally much larger than is normally needed.
513 % Blur "{radius},{sigma},{angle}"
514 % As per Gaussian, but generates a 1 dimensional or linear gaussian
515 % blur, at the angle given (current restricted to orthogonal angles).
516 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
518 % NOTE that two such blurs perpendicular to each other is equivelent to
519 % -blur and the previous gaussian, but is often 10 or more times faster.
521 % Comet "{width},{sigma},{angle}"
522 % Blur in one direction only, mush like how a bright object leaves
523 % a comet like trail. The Kernel is actually half a gaussian curve,
524 % Adding two such blurs in oppiste directions produces a Linear Blur.
526 % NOTE: that the first argument is the width of the kernel and not the
527 % radius of the kernel.
529 % # Still to be implemented...
531 % # Sharpen "{radius},{sigma}
532 % # Negated Gaussian (center zeroed and re-normalized),
533 % # with a 2 unit positive peak. -- Check On line documentation
535 % # Laplacian "{radius},{sigma}"
536 % # Laplacian (a mexican hat like) Function
538 % # LOG "{radius},{sigma1},{sigma2}
539 % # Laplacian of Gaussian
541 % # DOG "{radius},{sigma1},{sigma2}
542 % # Difference of two Gaussians
546 % # Set kernel values using a resize filter, and given scale (sigma)
547 % # Cylindrical or Linear. Is this posible with an image?
552 % Rectangle "{geometry}"
553 % Simply generate a rectangle of 1's with the size given. You can also
554 % specify the location of the 'control point', otherwise the closest
555 % pixel to the center of the rectangle is selected.
557 % Properly centered and odd sized rectangles work the best.
559 % Diamond "[{radius}[,{scale}]]"
560 % Generate a diamond shaped kernel with given radius to the points.
561 % Kernel size will again be radius*2+1 square and defaults to radius 1,
562 % generating a 3x3 kernel that is slightly larger than a square.
564 % Square "[{radius}[,{scale}]]"
565 % Generate a square shaped kernel of size radius*2+1, and defaulting
566 % to a 3x3 (radius 1).
568 % Note that using a larger radius for the "Square" or the "Diamond"
569 % is also equivelent to iterating the basic morphological method
570 % that many times. However However iterating with the smaller radius 1
571 % default is actually faster than using a larger kernel radius.
573 % Disk "[{radius}[,{scale}]]
574 % Generate a binary disk of the radius given, radius may be a float.
575 % Kernel size will be ceil(radius)*2+1 square.
576 % NOTE: Here are some disk shapes of specific interest
577 % "disk:1" => "diamond" or "cross:1"
578 % "disk:1.5" => "square"
579 % "disk:2" => "diamond:2"
580 % "disk:2.5" => a general disk shape of radius 2
581 % "disk:2.9" => "square:2"
582 % "disk:3.5" => default - octagonal/disk shape of radius 3
583 % "disk:4.2" => roughly octagonal shape of radius 4
584 % "disk:4.3" => a general disk shape of radius 4
585 % After this all the kernel shape becomes more and more circular.
587 % Because a "disk" is more circular when using a larger radius, using a
588 % larger radius is preferred over iterating the morphological operation.
590 % Plus "[{radius}[,{scale}]]"
591 % Generate a kernel in the shape of a 'plus' sign. The length of each
592 % arm is also the radius, which defaults to 2.
594 % This kernel is not a good general morphological kernel, but is used
595 % more for highlighting and marking any single pixels in an image using,
596 % a "Dilate" or "Erode" method as appropriate.
598 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
600 % Note that unlike other kernels iterating a plus does not produce the
601 % same result as using a larger radius for the cross.
603 % Distance Measuring Kernels
605 % Chebyshev "[{radius}][x{scale}[%!]]"
606 % Manhatten "[{radius}][x{scale}[%!]]"
607 % Euclidean "[{radius}][x{scale}[%!]]"
609 % Different types of distance measuring methods, which are used with the
610 % a 'Distance' morphology method for generating a gradient based on
611 % distance from an edge of a binary shape, though there is a technique
612 % for handling a anti-aliased shape.
614 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
615 % one to any neighbour, orthogonal or diagonal. One why of thinking of
616 % it is the number of squares a 'King' or 'Queen' in chess needs to
617 % traverse reach any other position on a chess board. It results in a
618 % 'square' like distance function, but one where diagonals are closer
621 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
622 % Cab metric), is the distance needed when you can only travel in
623 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
624 % in chess would travel. It results in a diamond like distances, where
625 % diagonals are further than expected.
627 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
628 % However by default the kernel size only has a radius of 1, which
629 % limits the distance to 'Knight' like moves, with only orthogonal and
630 % diagonal measurements being correct. As such for the default kernel
631 % you will get octagonal like distance function, which is reasonally
634 % However if you use a larger radius such as "Euclidean:4" you will
635 % get a much smoother distance gradient from the edge of the shape.
636 % Of course a larger kernel is slower to use, and generally not needed.
638 % To allow the use of fractional distances that you get with diagonals
639 % the actual distance is scaled by a fixed value which the user can
640 % provide. This is not actually nessary for either ""Chebyshev" or
641 % "Manhatten" distance kernels, but is done for all three distance
642 % kernels. If no scale is provided it is set to a value of 100,
643 % allowing for a maximum distance measurement of 655 pixels using a Q16
644 % version of IM, from any edge. However for small images this can
645 % result in quite a dark gradient.
647 % See the 'Distance' Morphological Method, for information of how it is
650 % # Hit-n-Miss Kernel-Lists -- Still to be implemented
652 % # specifically for Pruning, Thinning, Thickening
656 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
657 const GeometryInfo *args)
670 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
672 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
673 if (kernel == (KernelInfo *) NULL)
675 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
676 kernel->minimum = kernel->maximum = 0.0;
677 kernel->negative_range = kernel->positive_range = 0.0;
679 kernel->next = (KernelInfo *) NULL;
680 kernel->signature = MagickSignature;
683 /* Convolution Kernels */
686 sigma = fabs(args->sigma);
688 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
690 kernel->width = kernel->height =
691 GetOptimalKernelWidth2D(args->rho,sigma);
692 kernel->x = kernel->y = (long) (kernel->width-1)/2;
693 kernel->negative_range = kernel->positive_range = 0.0;
694 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
695 kernel->height*sizeof(double));
696 if (kernel->values == (double *) NULL)
697 return(DestroyKernelInfo(kernel));
699 sigma = 2.0*sigma*sigma; /* simplify the expression */
700 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
701 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
702 kernel->positive_range += (
704 exp(-((double)(u*u+v*v))/sigma)
705 /* / (MagickPI*sigma) */ );
707 kernel->maximum = kernel->values[
708 kernel->y*kernel->width+kernel->x ];
710 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
716 sigma = fabs(args->sigma);
718 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
720 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
721 kernel->x = (long) (kernel->width-1)/2;
724 kernel->negative_range = kernel->positive_range = 0.0;
725 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
726 kernel->height*sizeof(double));
727 if (kernel->values == (double *) NULL)
728 return(DestroyKernelInfo(kernel));
732 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
733 ** It generates a gaussian 3 times the width, and compresses it into
734 ** the expected range. This produces a closer normalization of the
735 ** resulting kernel, especially for very low sigma values.
736 ** As such while wierd it is prefered.
738 ** I am told this method originally came from Photoshop.
740 sigma *= KernelRank; /* simplify expanded curve */
741 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
742 (void) ResetMagickMemory(kernel->values,0, (size_t)
743 kernel->width*sizeof(double));
744 for ( u=-v; u <= v; u++) {
745 kernel->values[(u+v)/KernelRank] +=
746 exp(-((double)(u*u))/(2.0*sigma*sigma))
747 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
749 for (i=0; i < (long) kernel->width; i++)
750 kernel->positive_range += kernel->values[i];
752 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
753 kernel->positive_range += (
755 exp(-((double)(u*u))/(2.0*sigma*sigma))
756 /* / (MagickSQ2PI*sigma) */ );
759 kernel->maximum = kernel->values[ kernel->x ];
760 /* Note that neither methods above generate a normalized kernel,
761 ** though it gets close. The kernel may be 'clipped' by a user defined
762 ** radius, producing a smaller (darker) kernel. Also for very small
763 ** sigma's (> 0.1) the central value becomes larger than one, and thus
764 ** producing a very bright kernel.
767 /* Normalize the 1D Gaussian Kernel
769 ** Because of this the divisor in the above kernel generator is
770 ** not needed, so is not done above.
772 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
774 /* rotate the kernel by given angle */
775 RotateKernelInfo(kernel, args->xi);
780 sigma = fabs(args->sigma);
782 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
784 if ( args->rho < 1.0 )
785 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
787 kernel->width = (unsigned long)args->rho;
788 kernel->x = kernel->y = 0;
790 kernel->negative_range = kernel->positive_range = 0.0;
791 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
792 kernel->height*sizeof(double));
793 if (kernel->values == (double *) NULL)
794 return(DestroyKernelInfo(kernel));
796 /* A comet blur is half a gaussian curve, so that the object is
797 ** blurred in one direction only. This may not be quite the right
798 ** curve so may change in the future. The function must be normalised.
802 sigma *= KernelRank; /* simplify expanded curve */
803 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
804 (void) ResetMagickMemory(kernel->values,0, (size_t)
805 kernel->width*sizeof(double));
806 for ( u=0; u < v; u++) {
807 kernel->values[u/KernelRank] +=
808 exp(-((double)(u*u))/(2.0*sigma*sigma))
809 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
811 for (i=0; i < (long) kernel->width; i++)
812 kernel->positive_range += kernel->values[i];
814 for ( i=0; i < (long) kernel->width; i++)
815 kernel->positive_range += (
817 exp(-((double)(i*i))/(2.0*sigma*sigma))
818 /* / (MagickSQ2PI*sigma) */ );
821 kernel->maximum = kernel->values[0];
823 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
824 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
827 /* Boolean Kernels */
828 case RectangleKernel:
832 if ( type == SquareKernel )
835 kernel->width = kernel->height = 3; /* default radius = 1 */
837 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
838 kernel->x = kernel->y = (long) (kernel->width-1)/2;
842 /* NOTE: user defaults set in "AcquireKernelInfo()" */
843 if ( args->rho < 1.0 || args->sigma < 1.0 )
844 return(DestroyKernelInfo(kernel)); /* invalid args given */
845 kernel->width = (unsigned long)args->rho;
846 kernel->height = (unsigned long)args->sigma;
847 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
848 args->psi < 0.0 || args->psi > (double)kernel->height )
849 return(DestroyKernelInfo(kernel)); /* invalid args given */
850 kernel->x = (long) args->xi;
851 kernel->y = (long) args->psi;
854 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
855 kernel->height*sizeof(double));
856 if (kernel->values == (double *) NULL)
857 return(DestroyKernelInfo(kernel));
859 /* set all kernel values to 1.0 */
860 u=(long) kernel->width*kernel->height;
861 for ( i=0; i < u; i++)
862 kernel->values[i] = scale;
863 kernel->minimum = kernel->maximum = scale; /* a flat shape */
864 kernel->positive_range = scale*u;
870 kernel->width = kernel->height = 3; /* default radius = 1 */
872 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
873 kernel->x = kernel->y = (long) (kernel->width-1)/2;
875 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
876 kernel->height*sizeof(double));
877 if (kernel->values == (double *) NULL)
878 return(DestroyKernelInfo(kernel));
880 /* set all kernel values within diamond area to scale given */
881 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
882 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
883 if ((labs(u)+labs(v)) <= (long)kernel->x)
884 kernel->positive_range += kernel->values[i] = args->sigma;
886 kernel->values[i] = nan;
887 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
895 limit = (long)(args->rho*args->rho);
896 if (args->rho < 0.1) /* default radius approx 3.5 */
897 kernel->width = kernel->height = 7L, limit = 10L;
899 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
900 kernel->x = kernel->y = (long) (kernel->width-1)/2;
902 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
903 kernel->height*sizeof(double));
904 if (kernel->values == (double *) NULL)
905 return(DestroyKernelInfo(kernel));
907 /* set all kernel values within disk area to 1.0 */
908 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
909 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
910 if ((u*u+v*v) <= limit)
911 kernel->positive_range += kernel->values[i] = args->sigma;
913 kernel->values[i] = nan;
914 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
920 kernel->width = kernel->height = 5; /* default radius 2 */
922 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
923 kernel->x = kernel->y = (long) (kernel->width-1)/2;
925 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
926 kernel->height*sizeof(double));
927 if (kernel->values == (double *) NULL)
928 return(DestroyKernelInfo(kernel));
930 /* set all kernel values along axises to 1.0 */
931 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
932 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
933 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
934 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
935 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
938 /* Distance Measuring Kernels */
939 case ChebyshevKernel:
942 kernel->width = kernel->height = 3; /* default radius = 1 */
944 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
945 kernel->x = kernel->y = (long) (kernel->width-1)/2;
947 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
948 kernel->height*sizeof(double));
949 if (kernel->values == (double *) NULL)
950 return(DestroyKernelInfo(kernel));
952 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
953 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
954 kernel->positive_range += ( kernel->values[i] =
955 args->sigma*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
956 kernel->maximum = kernel->values[0];
959 case ManhattenKernel:
962 kernel->width = kernel->height = 3; /* default radius = 1 */
964 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
965 kernel->x = kernel->y = (long) (kernel->width-1)/2;
967 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
968 kernel->height*sizeof(double));
969 if (kernel->values == (double *) NULL)
970 return(DestroyKernelInfo(kernel));
972 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
973 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
974 kernel->positive_range += ( kernel->values[i] =
975 args->sigma*(labs(u)+labs(v)) );
976 kernel->maximum = kernel->values[0];
979 case EuclideanKernel:
982 kernel->width = kernel->height = 3; /* default radius = 1 */
984 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
985 kernel->x = kernel->y = (long) (kernel->width-1)/2;
987 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
988 kernel->height*sizeof(double));
989 if (kernel->values == (double *) NULL)
990 return(DestroyKernelInfo(kernel));
992 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
993 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
994 kernel->positive_range += ( kernel->values[i] =
995 args->sigma*sqrt((double)(u*u+v*v)) );
996 kernel->maximum = kernel->values[0];
999 /* Undefined Kernels */
1000 case LaplacianKernel:
1003 perror("Kernel Type has not been defined yet");
1006 /* Generate a No-Op minimal kernel - 1x1 pixel */
1007 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
1008 if (kernel->values == (double *) NULL)
1009 return(DestroyKernelInfo(kernel));
1010 kernel->width = kernel->height = 1;
1011 kernel->x = kernel->x = 0;
1012 kernel->type = UndefinedKernel;
1014 kernel->positive_range =
1015 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
1023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1027 % C l o n e K e r n e l I n f o %
1031 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1033 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
1034 % can be modified without effecting the original. The cloned kernel should
1035 % be destroyed using DestoryKernelInfo() when no longer needed.
1037 % The format of the CloneKernelInfo method is:
1039 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1041 % A description of each parameter follows:
1043 % o kernel: the Morphology/Convolution kernel to be cloned
1046 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1054 assert(kernel != (KernelInfo *) NULL);
1055 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1056 if (new_kernel == (KernelInfo *) NULL)
1058 *new_kernel=(*kernel); /* copy values in structure */
1060 /* replace the values with a copy of the values */
1061 new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1062 kernel->height*sizeof(double));
1063 if (new_kernel->values == (double *) NULL)
1064 return(DestroyKernelInfo(new_kernel));
1065 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1066 new_kernel->values[i]=kernel->values[i];
1068 /* Also clone the next kernel in the kernel list */
1069 if ( kernel->next != (KernelInfo *) NULL ) {
1070 new_kernel->next = CloneKernelInfo(kernel->next);
1071 if ( new_kernel->next == (KernelInfo *) NULL )
1072 return(DestroyKernelInfo(new_kernel));
1079 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1083 % D e s t r o y K e r n e l I n f o %
1087 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1089 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
1092 % The format of the DestroyKernelInfo method is:
1094 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1096 % A description of each parameter follows:
1098 % o kernel: the Morphology/Convolution kernel to be destroyed
1102 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1104 assert(kernel != (KernelInfo *) NULL);
1106 if ( kernel->next != (KernelInfo *) NULL )
1107 kernel->next = DestroyKernelInfo(kernel->next);
1109 kernel->values = (double *)RelinquishMagickMemory(kernel->values);
1110 kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
1115 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1119 % M o r p h o l o g y I m a g e C h a n n e l %
1123 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1125 % MorphologyImageChannel() applies a user supplied kernel to the image
1126 % according to the given mophology method.
1128 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1129 % by the kernel generator.
1131 % The format of the MorphologyImage method is:
1133 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
1134 % const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
1135 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1136 % channel,MorphologyMethod method,const long iterations,
1137 % KernelInfo *kernel,ExceptionInfo *exception)
1139 % A description of each parameter follows:
1141 % o image: the image.
1143 % o method: the morphology method to be applied.
1145 % o iterations: apply the operation this many times (or no change).
1146 % A value of -1 means loop until no change found.
1147 % How this is applied may depend on the morphology method.
1148 % Typically this is a value of 1.
1150 % o channel: the channel type.
1152 % o kernel: An array of double representing the morphology kernel.
1153 % Warning: kernel may be normalized for the Convolve method.
1155 % o exception: return any errors or warnings in this structure.
1158 % TODO: bias and auto-scale handling of the kernel for convolution
1159 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1160 % by the kernel generator.
1165 /* Internal function
1166 * Apply the low-level Morphology Method Primatives using the given Kernel
1167 * Returning the number of pixels that changed.
1168 * Two pre-created images must be provided, no image is created.
1170 static unsigned long MorphologyPrimative(const Image *image, Image
1171 *result_image, const MorphologyMethod method, const ChannelType channel,
1172 const KernelInfo *kernel, ExceptionInfo *exception)
1174 #define MorphologyTag "Morphology/Image"
1191 /* Only the most basic morphology is actually performed by this routine */
1194 Apply Basic Morphology to Image.
1200 GetMagickPixelPacket(image,&bias);
1201 SetMagickPixelPacketBias(image,&bias);
1202 /* Future: handle auto-bias from user, based on kernel input */
1204 p_view=AcquireCacheView(image);
1205 q_view=AcquireCacheView(result_image);
1207 /* Some methods (including convolve) needs use a reflected kernel.
1208 * Adjust 'origin' offsets for this reflected kernel.
1213 case ConvolveMorphology:
1214 case DilateMorphology:
1215 case DilateIntensityMorphology:
1216 case DistanceMorphology:
1217 /* kernel needs to used with reflection about origin */
1218 offx = (long) kernel->width-offx-1;
1219 offy = (long) kernel->height-offy-1;
1221 case ErodeMorphology:
1222 case ErodeIntensityMorphology:
1223 case HitAndMissMorphology:
1224 case ThinningMorphology:
1225 case ThickenMorphology:
1226 /* kernel is user as is, without reflection */
1229 perror("Not a low level Morphology Method");
1233 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1234 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1236 for (y=0; y < (long) image->rows; y++)
1241 register const PixelPacket
1244 register const IndexPacket
1245 *restrict p_indexes;
1247 register PixelPacket
1250 register IndexPacket
1251 *restrict q_indexes;
1259 if (status == MagickFalse)
1261 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1262 image->columns+kernel->width, kernel->height, exception);
1263 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1265 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1270 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1271 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1272 r = (image->columns+kernel->width)*offy+offx; /* constant */
1274 for (x=0; x < (long) image->columns; x++)
1282 register const double
1285 register const PixelPacket
1288 register const IndexPacket
1289 *restrict k_indexes;
1296 /* Copy input to ouput image for unused channels
1297 * This removes need for 'cloning' a new image every iteration
1300 if (image->colorspace == CMYKColorspace)
1301 q_indexes[x] = p_indexes[r];
1308 min.index = (MagickRealType) QuantumRange;
1313 max.index = (MagickRealType) 0;
1314 /* original pixel value */
1315 result.red = (MagickRealType) p[r].red;
1316 result.green = (MagickRealType) p[r].green;
1317 result.blue = (MagickRealType) p[r].blue;
1318 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1319 if ( image->colorspace == CMYKColorspace)
1320 result.index = (MagickRealType) p_indexes[r];
1323 case ConvolveMorphology:
1324 /* Set the user defined bias of the weighted average output
1326 ** FUTURE: provide some way for internal functions to disable
1327 ** user provided bias and scaling effects.
1331 case DilateIntensityMorphology:
1332 case ErodeIntensityMorphology:
1333 result.red = 0.0; /* flag indicating when first match found */
1340 case ConvolveMorphology:
1341 /* Weighted Average of pixels using reflected kernel
1343 ** NOTE for correct working of this operation for asymetrical
1344 ** kernels, the kernel needs to be applied in its reflected form.
1345 ** That is its values needs to be reversed.
1347 ** Correlation is actually the same as this but without reflecting
1348 ** the kernel, and thus 'lower-level' that Convolution. However
1349 ** as Convolution is the more common method used, and it does not
1350 ** really cost us much in terms of processing to use a reflected
1351 ** kernel, so it is Convolution that is implemented.
1353 ** Correlation will have its kernel reflected before calling
1354 ** this function to do a Convolve.
1356 ** For more details of Correlation vs Convolution see
1357 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
1359 if (((channel & SyncChannels) != 0 ) &&
1360 (image->matte == MagickTrue))
1361 { /* Channel has a 'Sync' Flag, and Alpha Channel enabled.
1362 ** Weight the color channels with Alpha Channel so that
1363 ** transparent pixels are not part of the results.
1366 alpha, /* color channel weighting : kernel*alpha */
1367 gamma; /* divisor, sum of weighting values */
1370 k = &kernel->values[ kernel->width*kernel->height-1 ];
1372 k_indexes = p_indexes;
1373 for (v=0; v < (long) kernel->height; v++) {
1374 for (u=0; u < (long) kernel->width; u++, k--) {
1375 if ( IsNan(*k) ) continue;
1376 alpha=(*k)*(QuantumScale*(QuantumRange-
1377 k_pixels[u].opacity));
1379 result.red += alpha*k_pixels[u].red;
1380 result.green += alpha*k_pixels[u].green;
1381 result.blue += alpha*k_pixels[u].blue;
1382 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1383 if ( image->colorspace == CMYKColorspace)
1384 result.index += alpha*k_indexes[u];
1386 k_pixels += image->columns+kernel->width;
1387 k_indexes += image->columns+kernel->width;
1389 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1390 result.red *= gamma;
1391 result.green *= gamma;
1392 result.blue *= gamma;
1393 result.opacity *= gamma;
1394 result.index *= gamma;
1398 /* No 'Sync' flag, or no Alpha involved.
1399 ** Convolution is simple individual channel weigthed sum.
1401 k = &kernel->values[ kernel->width*kernel->height-1 ];
1403 k_indexes = p_indexes;
1404 for (v=0; v < (long) kernel->height; v++) {
1405 for (u=0; u < (long) kernel->width; u++, k--) {
1406 if ( IsNan(*k) ) continue;
1407 result.red += (*k)*k_pixels[u].red;
1408 result.green += (*k)*k_pixels[u].green;
1409 result.blue += (*k)*k_pixels[u].blue;
1410 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1411 if ( image->colorspace == CMYKColorspace)
1412 result.index += (*k)*k_indexes[u];
1414 k_pixels += image->columns+kernel->width;
1415 k_indexes += image->columns+kernel->width;
1420 case ErodeMorphology:
1421 /* Minimum Value within kernel neighbourhood
1423 ** NOTE that the kernel is not reflected for this operation!
1425 ** NOTE: in normal Greyscale Morphology, the kernel value should
1426 ** be added to the real value, this is currently not done, due to
1427 ** the nature of the boolean kernels being used.
1431 k_indexes = p_indexes;
1432 for (v=0; v < (long) kernel->height; v++) {
1433 for (u=0; u < (long) kernel->width; u++, k++) {
1434 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1435 Minimize(min.red, (double) k_pixels[u].red);
1436 Minimize(min.green, (double) k_pixels[u].green);
1437 Minimize(min.blue, (double) k_pixels[u].blue);
1438 Minimize(min.opacity,
1439 QuantumRange-(double) k_pixels[u].opacity);
1440 if ( image->colorspace == CMYKColorspace)
1441 Minimize(min.index, (double) k_indexes[u]);
1443 k_pixels += image->columns+kernel->width;
1444 k_indexes += image->columns+kernel->width;
1449 case DilateMorphology:
1450 /* Maximum Value within kernel neighbourhood
1452 ** NOTE for correct working of this operation for asymetrical
1453 ** kernels, the kernel needs to be applied in its reflected form.
1454 ** That is its values needs to be reversed.
1456 ** NOTE: in normal Greyscale Morphology, the kernel value should
1457 ** be added to the real value, this is currently not done, due to
1458 ** the nature of the boolean kernels being used.
1461 k = &kernel->values[ kernel->width*kernel->height-1 ];
1463 k_indexes = p_indexes;
1464 for (v=0; v < (long) kernel->height; v++) {
1465 for (u=0; u < (long) kernel->width; u++, k--) {
1466 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1467 Maximize(max.red, (double) k_pixels[u].red);
1468 Maximize(max.green, (double) k_pixels[u].green);
1469 Maximize(max.blue, (double) k_pixels[u].blue);
1470 Maximize(max.opacity,
1471 QuantumRange-(double) k_pixels[u].opacity);
1472 if ( image->colorspace == CMYKColorspace)
1473 Maximize(max.index, (double) k_indexes[u]);
1475 k_pixels += image->columns+kernel->width;
1476 k_indexes += image->columns+kernel->width;
1480 case HitAndMissMorphology:
1481 case ThinningMorphology:
1482 case ThickenMorphology:
1483 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
1485 ** NOTE that the kernel is not reflected for this operation,
1486 ** and consists of both foreground and background pixel
1487 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
1488 ** with either Nan or 0.5 values for don't care.
1490 ** Note that this can produce negative results, though really
1491 ** only a positive match has any real value.
1495 k_indexes = p_indexes;
1496 for (v=0; v < (long) kernel->height; v++) {
1497 for (u=0; u < (long) kernel->width; u++, k++) {
1498 if ( IsNan(*k) ) continue;
1500 { /* minimim of foreground pixels */
1501 Minimize(min.red, (double) k_pixels[u].red);
1502 Minimize(min.green, (double) k_pixels[u].green);
1503 Minimize(min.blue, (double) k_pixels[u].blue);
1504 Minimize(min.opacity,
1505 QuantumRange-(double) k_pixels[u].opacity);
1506 if ( image->colorspace == CMYKColorspace)
1507 Minimize(min.index, (double) k_indexes[u]);
1509 else if ( (*k) < 0.3 )
1510 { /* maximum of background pixels */
1511 Maximize(max.red, (double) k_pixels[u].red);
1512 Maximize(max.green, (double) k_pixels[u].green);
1513 Maximize(max.blue, (double) k_pixels[u].blue);
1514 Maximize(max.opacity,
1515 QuantumRange-(double) k_pixels[u].opacity);
1516 if ( image->colorspace == CMYKColorspace)
1517 Maximize(max.index, (double) k_indexes[u]);
1520 k_pixels += image->columns+kernel->width;
1521 k_indexes += image->columns+kernel->width;
1523 /* Pattern Match only if min fg larger than min bg pixels */
1524 min.red -= max.red; Maximize( min.red, 0.0 );
1525 min.green -= max.green; Maximize( min.green, 0.0 );
1526 min.blue -= max.blue; Maximize( min.blue, 0.0 );
1527 min.opacity -= max.opacity; Maximize( min.opacity, 0.0 );
1528 min.index -= max.index; Maximize( min.index, 0.0 );
1531 case ErodeIntensityMorphology:
1532 /* Select Pixel with Minimum Intensity within kernel neighbourhood
1534 ** WARNING: the intensity test fails for CMYK and does not
1535 ** take into account the moderating effect of teh alpha channel
1536 ** on the intensity.
1538 ** NOTE that the kernel is not reflected for this operation!
1542 k_indexes = p_indexes;
1543 for (v=0; v < (long) kernel->height; v++) {
1544 for (u=0; u < (long) kernel->width; u++, k++) {
1545 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1546 if ( result.red == 0.0 ||
1547 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1548 /* copy the whole pixel - no channel selection */
1550 if ( result.red > 0.0 ) changed++;
1554 k_pixels += image->columns+kernel->width;
1555 k_indexes += image->columns+kernel->width;
1559 case DilateIntensityMorphology:
1560 /* Select Pixel with Maximum Intensity within kernel neighbourhood
1562 ** WARNING: the intensity test fails for CMYK and does not
1563 ** take into account the moderating effect of teh alpha channel
1564 ** on the intensity.
1566 ** NOTE for correct working of this operation for asymetrical
1567 ** kernels, the kernel needs to be applied in its reflected form.
1568 ** That is its values needs to be reversed.
1570 k = &kernel->values[ kernel->width*kernel->height-1 ];
1572 k_indexes = p_indexes;
1573 for (v=0; v < (long) kernel->height; v++) {
1574 for (u=0; u < (long) kernel->width; u++, k--) {
1575 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1576 if ( result.red == 0.0 ||
1577 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1578 /* copy the whole pixel - no channel selection */
1580 if ( result.red > 0.0 ) changed++;
1584 k_pixels += image->columns+kernel->width;
1585 k_indexes += image->columns+kernel->width;
1590 case DistanceMorphology:
1591 /* Add kernel Value and select the minimum value found.
1592 ** The result is a iterative distance from edge of image shape.
1594 ** All Distance Kernels are symetrical, but that may not always
1595 ** be the case. For example how about a distance from left edges?
1596 ** To work correctly with asymetrical kernels the reflected kernel
1597 ** needs to be applied.
1599 ** Actually this is really a GreyErode with a negative kernel!
1602 k = &kernel->values[ kernel->width*kernel->height-1 ];
1604 k_indexes = p_indexes;
1605 for (v=0; v < (long) kernel->height; v++) {
1606 for (u=0; u < (long) kernel->width; u++, k--) {
1607 if ( IsNan(*k) ) continue;
1608 Minimize(result.red, (*k)+k_pixels[u].red);
1609 Minimize(result.green, (*k)+k_pixels[u].green);
1610 Minimize(result.blue, (*k)+k_pixels[u].blue);
1611 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1612 if ( image->colorspace == CMYKColorspace)
1613 Minimize(result.index, (*k)+k_indexes[u]);
1615 k_pixels += image->columns+kernel->width;
1616 k_indexes += image->columns+kernel->width;
1620 case UndefinedMorphology:
1622 break; /* Do nothing */
1624 /* Final mathematics of results (combine with original image?)
1626 ** NOTE: Difference Morphology operators Edge* and *Hat could also
1627 ** be done here but works better with iteration as a image difference
1628 ** in the controling function (below). Thicken and Thinning however
1629 ** should be done here so thay can be iterated correctly.
1632 case HitAndMissMorphology:
1633 case ErodeMorphology:
1634 result = min; /* minimum of neighbourhood */
1636 case DilateMorphology:
1637 result = max; /* maximum of neighbourhood */
1639 case ThinningMorphology:
1640 /* subtract pattern match from original */
1641 result.red -= min.red;
1642 result.green -= min.green;
1643 result.blue -= min.blue;
1644 result.opacity -= min.opacity;
1645 result.index -= min.index;
1647 case ThickenMorphology:
1648 /* Union with original image (maximize) - or should this be + */
1649 Maximize( result.red, min.red );
1650 Maximize( result.green, min.green );
1651 Maximize( result.blue, min.blue );
1652 Maximize( result.opacity, min.opacity );
1653 Maximize( result.index, min.index );
1656 /* result directly calculated or assigned */
1659 /* Assign the resulting pixel values - Clamping Result */
1661 case UndefinedMorphology:
1662 case DilateIntensityMorphology:
1663 case ErodeIntensityMorphology:
1664 break; /* full pixel was directly assigned - not a channel method */
1666 if ((channel & RedChannel) != 0)
1667 q->red = ClampToQuantum(result.red);
1668 if ((channel & GreenChannel) != 0)
1669 q->green = ClampToQuantum(result.green);
1670 if ((channel & BlueChannel) != 0)
1671 q->blue = ClampToQuantum(result.blue);
1672 if ((channel & OpacityChannel) != 0
1673 && image->matte == MagickTrue )
1674 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1675 if ((channel & IndexChannel) != 0
1676 && image->colorspace == CMYKColorspace)
1677 q_indexes[x] = ClampToQuantum(result.index);
1680 /* Count up changed pixels */
1681 if ( ( p[r].red != q->red )
1682 || ( p[r].green != q->green )
1683 || ( p[r].blue != q->blue )
1684 || ( p[r].opacity != q->opacity )
1685 || ( image->colorspace == CMYKColorspace &&
1686 p_indexes[r] != q_indexes[x] ) )
1687 changed++; /* The pixel had some value changed! */
1691 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1692 if (sync == MagickFalse)
1694 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1699 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1700 #pragma omp critical (MagickCore_MorphologyImage)
1702 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1703 if (proceed == MagickFalse)
1707 result_image->type=image->type;
1708 q_view=DestroyCacheView(q_view);
1709 p_view=DestroyCacheView(p_view);
1710 return(status ? (unsigned long) changed : 0);
1714 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1715 method, const long iterations,const KernelInfo *kernel, ExceptionInfo
1721 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1722 iterations,kernel,exception);
1723 return(morphology_image);
1727 MagickExport Image *MorphologyImageChannel(const Image *image,
1728 const ChannelType channel,const MorphologyMethod method,
1729 const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
1750 assert(image != (Image *) NULL);
1751 assert(image->signature == MagickSignature);
1752 assert(kernel != (KernelInfo *) NULL);
1753 assert(kernel->signature == MagickSignature);
1754 assert(exception != (ExceptionInfo *) NULL);
1755 assert(exception->signature == MagickSignature);
1757 if ( iterations == 0 )
1758 return((Image *)NULL); /* null operation - nothing to do! */
1760 /* kernel must be valid at this point
1761 * (except maybe for posible future morphology methods like "Prune"
1763 assert(kernel != (KernelInfo *)NULL);
1765 new_image = (Image *) NULL; /* for cleanup */
1766 old_image = (Image *) NULL;
1767 grad_image = (Image *) NULL;
1768 curr_kernel = (KernelInfo *) NULL;
1769 count = 0; /* interation count */
1770 changed = 1; /* was last run succesfull */
1771 curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
1772 curr_method = method; /* to be changed as nessary */
1774 limit = (unsigned long) iterations;
1775 if ( iterations < 0 )
1776 limit = image->columns > image->rows ? image->columns : image->rows;
1778 /* Third-level morphology methods */
1779 switch( curr_method ) {
1780 case EdgeMorphology:
1781 grad_image = MorphologyImageChannel(image, channel,
1782 DilateMorphology, iterations, curr_kernel, exception);
1783 if ( grad_image == (Image *) NULL )
1786 case EdgeInMorphology:
1787 curr_method = ErodeMorphology;
1789 case EdgeOutMorphology:
1790 curr_method = DilateMorphology;
1792 case TopHatMorphology:
1793 curr_method = OpenMorphology;
1795 case BottomHatMorphology:
1796 curr_method = CloseMorphology;
1799 break; /* not a third-level method */
1802 /* Second-level morphology methods */
1803 switch( curr_method ) {
1804 case OpenMorphology:
1805 /* Open is a Erode then a Dilate without reflection */
1806 new_image = MorphologyImageChannel(image, channel,
1807 ErodeMorphology, iterations, curr_kernel, exception);
1808 if (new_image == (Image *) NULL)
1810 curr_method = DilateMorphology;
1812 case OpenIntensityMorphology:
1813 new_image = MorphologyImageChannel(image, channel,
1814 ErodeIntensityMorphology, iterations, curr_kernel, exception);
1815 if (new_image == (Image *) NULL)
1817 curr_method = DilateIntensityMorphology;
1820 case CloseMorphology:
1821 /* Close is a Dilate then Erode using reflected kernel */
1822 /* A reflected kernel is needed for a Close */
1823 if ( curr_kernel == kernel )
1824 curr_kernel = CloneKernelInfo(kernel);
1825 if (curr_kernel == (KernelInfo *) NULL)
1827 RotateKernelInfo(curr_kernel,180);
1828 new_image = MorphologyImageChannel(image, channel,
1829 DilateMorphology, iterations, curr_kernel, exception);
1830 if (new_image == (Image *) NULL)
1832 curr_method = ErodeMorphology;
1834 case CloseIntensityMorphology:
1835 /* A reflected kernel is needed for a Close */
1836 if ( curr_kernel == kernel )
1837 curr_kernel = CloneKernelInfo(kernel);
1838 if (curr_kernel == (KernelInfo *) NULL)
1840 RotateKernelInfo(curr_kernel,180);
1841 new_image = MorphologyImageChannel(image, channel,
1842 DilateIntensityMorphology, iterations, curr_kernel, exception);
1843 if (new_image == (Image *) NULL)
1845 curr_method = ErodeIntensityMorphology;
1848 case CorrelateMorphology:
1849 /* A Correlation is actually a Convolution with a reflected kernel.
1850 ** However a Convolution is a weighted sum with a reflected kernel.
1851 ** It may seem stange to convert a Correlation into a Convolution
1852 ** as the Correleation is the simplier method, but Convolution is
1853 ** much more commonly used, and it makes sense to implement it directly
1854 ** so as to avoid the need to duplicate the kernel when it is not
1855 ** required (which is typically the default).
1857 if ( curr_kernel == kernel )
1858 curr_kernel = CloneKernelInfo(kernel);
1859 if (curr_kernel == (KernelInfo *) NULL)
1861 RotateKernelInfo(curr_kernel,180);
1862 curr_method = ConvolveMorphology;
1863 /* FALL-THRU into Correlate (weigthed sum without reflection) */
1865 case ConvolveMorphology:
1866 { /* Scale or Normalize kernel, according to user wishes
1867 ** before using it for the Convolve/Correlate method.
1869 ** FUTURE: provide some way for internal functions to disable
1870 ** user bias and scaling effects.
1873 *artifact = GetImageArtifact(image,"convolve:scale");
1874 if ( artifact != (char *)NULL ) {
1880 if ( curr_kernel == kernel )
1881 curr_kernel = CloneKernelInfo(kernel);
1882 if (curr_kernel == (KernelInfo *) NULL)
1885 flags = (GeometryFlags) ParseGeometry(artifact, &args);
1886 ScaleKernelInfo(curr_kernel, args.rho, flags);
1889 /* FALL-THRU to do the first, and typically the only iteration */
1892 /* Do a single iteration using the Low-Level Morphology method!
1893 ** This ensures a "new_image" has been generated, but allows us to skip
1894 ** the creation of 'old_image' if no more iterations are needed.
1896 ** The "curr_method" should also be set to a low-level method that is
1897 ** understood by the MorphologyPrimative() internal function.
1899 new_image=CloneImage(image,0,0,MagickTrue,exception);
1900 if (new_image == (Image *) NULL)
1902 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1904 InheritException(exception,&new_image->exception);
1907 count++; /* iteration count */
1908 changed = MorphologyPrimative(image,new_image,curr_method,channel,
1909 curr_kernel, exception);
1910 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1911 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
1912 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1913 count, 0L, changed);
1917 * + "curr_method" should be set to a low-level morphology method
1918 * + "count=1" if the first iteration of the first kernel has been done.
1919 * + "new_image" is defined and contains the current resulting image
1922 /* The Hit-And-Miss Morphology is not 'iterated' against the resulting
1923 ** image from the previous morphology step. It must always be applied
1924 ** to the original image, with the results collected into a union (maximimum
1925 ** or lighten) of all the results. Multiple kernels may be applied but
1926 ** an iteration of the morphology does nothing, so is ignored.
1928 ** The first kernel is guranteed to have been done, so lets continue the
1929 ** process, with the other kernels in the kernel list.
1931 if ( method == HitAndMissMorphology )
1933 if ( curr_kernel->next != (KernelInfo *) NULL ) {
1934 /* create a second working image */
1935 old_image = CloneImage(image,0,0,MagickTrue,exception);
1936 if (old_image == (Image *) NULL)
1938 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1940 InheritException(exception,&old_image->exception);
1944 /* loop through rest of the kernels */
1945 this_kernel=curr_kernel->next;
1947 while( this_kernel != (KernelInfo *) NULL )
1949 changed = MorphologyPrimative(image,old_image,curr_method,channel,
1950 this_kernel,exception);
1951 (void) CompositeImageChannel(new_image,
1952 (ChannelType) (channel & ~SyncChannels), LightenCompositeOp,
1954 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1955 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
1956 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1957 count, kernel_number, changed);
1958 this_kernel = this_kernel->next;
1961 old_image=DestroyImage(old_image);
1966 /* Repeat the low-level morphology over all kernels
1967 until iteration count limit or no change from any kernel is found */
1968 if ( ( count < limit && changed > 0 ) ||
1969 curr_kernel->next != (KernelInfo *) NULL ) {
1971 /* create a second working image */
1972 old_image = CloneImage(image,0,0,MagickTrue,exception);
1973 if (old_image == (Image *) NULL)
1975 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1977 InheritException(exception,&old_image->exception);
1981 /* reset variables for the first/next iteration, or next kernel) */
1983 this_kernel = curr_kernel;
1984 total_changed = count != 0 ? changed : 0;
1985 if ( count != 0 && this_kernel != (KernelInfo *) NULL ) {
1986 count = 0; /* first iteration is not yet finished! */
1987 this_kernel = curr_kernel->next;
1989 total_changed = changed;
1992 while ( count < limit ) {
1994 while ( this_kernel != (KernelInfo *) NULL ) {
1995 Image *tmp = old_image;
1996 old_image = new_image;
1998 changed = MorphologyPrimative(old_image,new_image,curr_method,channel,
1999 this_kernel,exception);
2000 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2001 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
2002 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2003 count, kernel_number, changed);
2004 total_changed += changed;
2005 this_kernel = this_kernel->next;
2008 if ( total_changed == 0 )
2009 break; /* no changes after processing all kernels - ABORT */
2010 /* prepare for next loop */
2013 this_kernel = curr_kernel;
2015 old_image=DestroyImage(old_image);
2018 /* finished with kernel - destary any copy that was made */
2019 if ( curr_kernel != kernel )
2020 curr_kernel=DestroyKernelInfo(curr_kernel);
2022 /* Third-level Subtractive methods post-processing
2024 ** The removal of any 'Sync' channel flag in the Image Compositon below
2025 ** ensures the compose method is applied in a purely mathematical way, only
2026 ** the selected channels, without any normal 'alpha blending' normally
2027 ** associated with the compose method.
2029 ** Note "method" here is the 'original' morphological method, and not the
2030 ** 'current' morphological method used above to generate "new_image".
2033 case EdgeOutMorphology:
2034 case EdgeInMorphology:
2035 case TopHatMorphology:
2036 case BottomHatMorphology:
2037 /* Get Difference relative to the original image */
2038 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2039 DifferenceCompositeOp, image, 0, 0);
2041 case EdgeMorphology:
2042 /* Difference the Eroded image from the saved Dilated image */
2043 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2044 DifferenceCompositeOp, grad_image, 0, 0);
2045 grad_image=DestroyImage(grad_image);
2053 /* Yes goto's are bad, but in this case it makes cleanup lot more efficient */
2055 if ( new_image != (Image *) NULL )
2056 DestroyImage(new_image);
2058 if ( old_image != (Image *) NULL )
2059 DestroyImage(old_image);
2060 if ( curr_kernel != (KernelInfo *) NULL && curr_kernel != kernel )
2061 curr_kernel=DestroyKernelInfo(curr_kernel);
2066 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2070 + R o t a t e K e r n e l I n f o %
2074 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2076 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
2077 % restricted to 90 degree angles, but this may be improved in the future.
2079 % The format of the RotateKernelInfo method is:
2081 % void RotateKernelInfo(KernelInfo *kernel, double angle)
2083 % A description of each parameter follows:
2085 % o kernel: the Morphology/Convolution kernel
2087 % o angle: angle to rotate in degrees
2089 % This function is only internel to this module, as it is not finalized,
2090 % especially with regard to non-orthogonal angles, and rotation of larger
2093 static void RotateKernelInfo(KernelInfo *kernel, double angle)
2095 /* angle the lower kernels first */
2096 if ( kernel->next != (KernelInfo *) NULL)
2097 RotateKernelInfo(kernel->next, angle);
2099 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
2101 ** TODO: expand beyond simple 90 degree rotates, flips and flops
2104 /* Modulus the angle */
2105 angle = fmod(angle, 360.0);
2109 if ( 315.0 < angle || angle <= 45.0 )
2110 return; /* no change! - At least at this time */
2112 switch (kernel->type) {
2113 /* These built-in kernels are cylindrical kernels, rotating is useless */
2114 case GaussianKernel:
2115 case LaplacianKernel:
2119 case ChebyshevKernel:
2120 case ManhattenKernel:
2121 case EuclideanKernel:
2124 /* These may be rotatable at non-90 angles in the future */
2125 /* but simply rotating them in multiples of 90 degrees is useless */
2131 /* These only allows a +/-90 degree rotation (by transpose) */
2132 /* A 180 degree rotation is useless */
2134 case RectangleKernel:
2135 if ( 135.0 < angle && angle <= 225.0 )
2137 if ( 225.0 < angle && angle <= 315.0 )
2141 /* these are freely rotatable in 90 degree units */
2143 case UndefinedKernel:
2144 case UserDefinedKernel:
2147 if ( 135.0 < angle && angle <= 225.0 )
2149 /* For a 180 degree rotation - also know as a reflection */
2150 /* This is actually a very very common operation! */
2151 /* Basically all that is needed is a reversal of the kernel data! */
2158 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
2159 t=k[i], k[i]=k[j], k[j]=t;
2161 kernel->x = (long) kernel->width - kernel->x - 1;
2162 kernel->y = (long) kernel->height - kernel->y - 1;
2163 angle = fmod(angle+180.0, 360.0);
2165 if ( 45.0 < angle && angle <= 135.0 )
2166 { /* Do a transpose and a flop, of the image, which results in a 90
2167 * degree rotation using two mirror operations.
2169 * WARNING: this assumes the original image was a 1 dimentional image
2170 * but currently that is the only built-ins it is applied to.
2174 t = (long) kernel->width;
2175 kernel->width = kernel->height;
2176 kernel->height = (unsigned long) t;
2178 kernel->x = kernel->y;
2180 angle = fmod(450.0 - angle, 360.0);
2182 /* At this point angle should be between -45 (315) and +45 degrees
2183 * In the future some form of non-orthogonal angled rotates could be
2184 * performed here, posibily with a linear kernel restriction.
2188 Not currently in use!
2189 { /* Do a flop, this assumes kernel is horizontally symetrical.
2190 * Each row of the kernel needs to be reversed!
2199 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2200 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2201 t=k[x], k[x]=k[r], k[r]=t;
2203 kernel->x = kernel->width - kernel->x - 1;
2204 angle = fmod(angle+180.0, 360.0);
2211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2215 % S c a l e K e r n e l I n f o %
2219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2221 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
2222 % without normalization of the sum of the kernel values (as per given flags).
2224 % By default (no flags given) the values within the kernel is scaled
2225 % directly using given scaling factor without change.
2227 % If any 'normalize_flags' are given the kernel will first be normalized and
2228 % then further scaled by the scaling factor value given. A 'PercentValue'
2229 % flag will cause the given scaling factor to be divided by one hundred
2232 % Kernel normalization ('normalize_flags' given) is designed to ensure that
2233 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
2234 % morphology methods will fall into -1.0 to +1.0 range. Note that for
2235 % non-HDRI versions of IM this may cause images to have any negative results
2236 % clipped, unless some 'bias' is used.
2238 % More specifically. Kernels which only contain positive values (such as a
2239 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
2240 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
2242 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
2243 % the kernel will be scaled by the absolute of the sum of kernel values, so
2244 % that it will generally fall within the +/- 1.0 range.
2246 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
2247 % will be scaled by just the sum of the postive values, so that its output
2248 % range will again fall into the +/- 1.0 range.
2250 % For special kernels designed for locating shapes using 'Correlate', (often
2251 % only containing +1 and -1 values, representing foreground/brackground
2252 % matching) a special normalization method is provided to scale the positive
2253 % values seperatally to those of the negative values, so the kernel will be
2254 % forced to become a zero-sum kernel better suited to such searches.
2256 % WARNING: Correct normalization of the kernel assumes that the '*_range'
2257 % attributes within the kernel structure have been correctly set during the
2260 % NOTE: The values used for 'normalize_flags' have been selected specifically
2261 % to match the use of geometry options, so that '!' means NormalizeValue,
2262 % '^' means CorrelateNormalizeValue, and '%' means PercentValue. All other
2263 % GeometryFlags values are ignored.
2265 % The format of the ScaleKernelInfo method is:
2267 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
2268 % const MagickStatusType normalize_flags )
2270 % A description of each parameter follows:
2272 % o kernel: the Morphology/Convolution kernel
2275 % multiply all values (after normalization) by this factor if not
2276 % zero. If the kernel is normalized regardless of any flags.
2278 % o normalize_flags:
2279 % GeometryFlags defining normalization method to use.
2280 % specifically: NormalizeValue, CorrelateNormalizeValue,
2281 % and/or PercentValue
2283 % This function is internal to this module only at this time, but can be
2284 % exported to other modules if needed.
2286 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
2287 const double scaling_factor,const GeometryFlags normalize_flags)
2296 /* scale the lower kernels first */
2297 if ( kernel->next != (KernelInfo *) NULL)
2298 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
2301 if ( (normalize_flags&NormalizeValue) != 0 ) {
2302 /* normalize kernel appropriately */
2303 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
2304 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
2306 pos_scale = kernel->positive_range; /* special zero-summing kernel */
2308 /* force kernel into being a normalized zero-summing kernel */
2309 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
2310 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
2311 ? kernel->positive_range : 1.0;
2312 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
2313 ? -kernel->negative_range : 1.0;
2316 neg_scale = pos_scale;
2318 /* finialize scaling_factor for positive and negative components */
2319 pos_scale = scaling_factor/pos_scale;
2320 neg_scale = scaling_factor/neg_scale;
2321 if ( (normalize_flags&PercentValue) != 0 ) {
2326 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2327 if ( ! IsNan(kernel->values[i]) )
2328 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
2330 /* convolution output range */
2331 kernel->positive_range *= pos_scale;
2332 kernel->negative_range *= neg_scale;
2333 /* maximum and minimum values in kernel */
2334 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
2335 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
2337 /* swap kernel settings if user scaling factor is negative */
2338 if ( scaling_factor < MagickEpsilon ) {
2340 t = kernel->positive_range;
2341 kernel->positive_range = kernel->negative_range;
2342 kernel->negative_range = t;
2343 t = kernel->maximum;
2344 kernel->maximum = kernel->minimum;
2345 kernel->minimum = 1;
2352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2356 + S h o w K e r n e l I n f o %
2360 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2362 % ShowKernelInfo() outputs the details of the given kernel defination to
2363 % standard error, generally due to a users 'showkernel' option request.
2365 % The format of the ShowKernel method is:
2367 % void ShowKernelInfo(KernelInfo *kernel)
2369 % A description of each parameter follows:
2371 % o kernel: the Morphology/Convolution kernel
2373 % This function is internal to this module only at this time. That may change
2376 MagickExport void ShowKernelInfo(KernelInfo *kernel)
2384 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
2386 fprintf(stderr, "Kernel ");
2387 if ( kernel->next != (KernelInfo *) NULL )
2388 fprintf(stderr, " #%ld", c );
2389 fprintf(stderr, " \"%s\" of size %lux%lu%+ld%+ld ",
2390 MagickOptionToMnemonic(MagickKernelOptions, k->type),
2391 kernel->width, k->height,
2392 kernel->x, kernel->y );
2394 " with values from %.*lg to %.*lg\n",
2395 GetMagickPrecision(), k->minimum,
2396 GetMagickPrecision(), k->maximum);
2397 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2398 GetMagickPrecision(), k->negative_range,
2399 GetMagickPrecision(), k->positive_range,
2400 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2401 for (i=v=0; v < (long) k->height; v++) {
2402 fprintf(stderr,"%2ld:",v);
2403 for (u=0; u < (long) k->width; u++, i++)
2404 if ( IsNan(k->values[i]) )
2405 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2407 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2408 GetMagickPrecision(), k->values[i]);
2409 fprintf(stderr,"\n");
2415 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2419 + Z e r o K e r n e l N a n s %
2423 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2425 % ZeroKernelNans() replaces any special 'nan' value that may be present in
2426 % the kernel with a zero value. This is typically done when the kernel will
2427 % be used in special hardware (GPU) convolution processors, to simply
2430 % The format of the ZeroKernelNans method is:
2432 % voidZeroKernelNans (KernelInfo *kernel)
2434 % A description of each parameter follows:
2436 % o kernel: the Morphology/Convolution kernel
2438 % FUTURE: return the information in a string for API usage.
2440 MagickExport void ZeroKernelNans(KernelInfo *kernel)
2445 /* scale the lower kernels first */
2446 if ( kernel->next != (KernelInfo *) NULL)
2447 ZeroKernelNans(kernel->next);
2449 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2450 if ( IsNan(kernel->values[i]) )
2451 kernel->values[i] = 0.0;