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];
434 kernel = last_kernel = NULL;
437 while ( GetMagickToken(p,NULL,token), *token != '\0' ) {
439 /* ignore multiple ';' following each other */
440 if ( *token != ';' ) {
442 /* tokens starting with alpha is a Named kernel */
443 if (isalpha((int) *token) == 0)
444 new_kernel = ParseKernelArray(p);
445 else /* otherwise a user defined kernel array */
446 new_kernel = ParseNamedKernel(p);
448 /* Error handling -- this is not proper error handling! */
449 if ( new_kernel == (KernelInfo *) NULL ) {
450 fprintf(stderr, "Failed to parse kernel number #%lu\n", kernel_number);
451 if ( kernel != (KernelInfo *) NULL )
452 kernel=DestroyKernelInfo(kernel);
453 return((KernelInfo *) NULL);
456 /* initialise or append the kernel list */
457 if ( last_kernel == (KernelInfo *) NULL )
458 kernel = last_kernel = new_kernel;
460 last_kernel->next = new_kernel;
461 last_kernel = new_kernel;
465 /* look for the next kernel in list */
467 if ( p == (char *) NULL )
477 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
481 % A c q u i r e K e r n e l B u i l t I n %
485 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
487 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
488 % kernels used for special purposes such as gaussian blurring, skeleton
489 % pruning, and edge distance determination.
491 % They take a KernelType, and a set of geometry style arguments, which were
492 % typically decoded from a user supplied string, or from a more complex
493 % Morphology Method that was requested.
495 % The format of the AcquireKernalBuiltIn method is:
497 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
498 % const GeometryInfo args)
500 % A description of each parameter follows:
502 % o type: the pre-defined type of kernel wanted
504 % o args: arguments defining or modifying the kernel
506 % Convolution Kernels
508 % Gaussian "{radius},{sigma}"
509 % Generate a two-dimentional gaussian kernel, as used by -gaussian
510 % A sigma is required, (with the 'x'), due to historical reasons.
512 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
513 % the final size of the resulting kernel to a square 2*radius+1 in size.
514 % The radius should be at least 2 times that of the sigma value, or
515 % sever clipping and aliasing may result. If not given or set to 0 the
516 % radius will be determined so as to produce the best minimal error
517 % result, which is usally much larger than is normally needed.
519 % Blur "{radius},{sigma},{angle}"
520 % As per Gaussian, but generates a 1 dimensional or linear gaussian
521 % blur, at the angle given (current restricted to orthogonal angles).
522 % If a 'radius' is given the kernel is clipped to a width of 2*radius+1.
524 % NOTE that two such blurs perpendicular to each other is equivelent to
525 % -blur and the previous gaussian, but is often 10 or more times faster.
527 % Comet "{width},{sigma},{angle}"
528 % Blur in one direction only, mush like how a bright object leaves
529 % a comet like trail. The Kernel is actually half a gaussian curve,
530 % Adding two such blurs in oppiste directions produces a Linear Blur.
532 % NOTE: that the first argument is the width of the kernel and not the
533 % radius of the kernel.
535 % # Still to be implemented...
537 % # Sharpen "{radius},{sigma}
538 % # Negated Gaussian (center zeroed and re-normalized),
539 % # with a 2 unit positive peak. -- Check On line documentation
541 % # Laplacian "{radius},{sigma}"
542 % # Laplacian (a mexican hat like) Function
544 % # LOG "{radius},{sigma1},{sigma2}
545 % # Laplacian of Gaussian
547 % # DOG "{radius},{sigma1},{sigma2}
548 % # Difference of two Gaussians
552 % # Set kernel values using a resize filter, and given scale (sigma)
553 % # Cylindrical or Linear. Is this posible with an image?
558 % Rectangle "{geometry}"
559 % Simply generate a rectangle of 1's with the size given. You can also
560 % specify the location of the 'control point', otherwise the closest
561 % pixel to the center of the rectangle is selected.
563 % Properly centered and odd sized rectangles work the best.
565 % Diamond "[{radius}[,{scale}]]"
566 % Generate a diamond shaped kernel with given radius to the points.
567 % Kernel size will again be radius*2+1 square and defaults to radius 1,
568 % generating a 3x3 kernel that is slightly larger than a square.
570 % Square "[{radius}[,{scale}]]"
571 % Generate a square shaped kernel of size radius*2+1, and defaulting
572 % to a 3x3 (radius 1).
574 % Note that using a larger radius for the "Square" or the "Diamond"
575 % is also equivelent to iterating the basic morphological method
576 % that many times. However However iterating with the smaller radius 1
577 % default is actually faster than using a larger kernel radius.
579 % Disk "[{radius}[,{scale}]]
580 % Generate a binary disk of the radius given, radius may be a float.
581 % Kernel size will be ceil(radius)*2+1 square.
582 % NOTE: Here are some disk shapes of specific interest
583 % "disk:1" => "diamond" or "cross:1"
584 % "disk:1.5" => "square"
585 % "disk:2" => "diamond:2"
586 % "disk:2.5" => a general disk shape of radius 2
587 % "disk:2.9" => "square:2"
588 % "disk:3.5" => default - octagonal/disk shape of radius 3
589 % "disk:4.2" => roughly octagonal shape of radius 4
590 % "disk:4.3" => a general disk shape of radius 4
591 % After this all the kernel shape becomes more and more circular.
593 % Because a "disk" is more circular when using a larger radius, using a
594 % larger radius is preferred over iterating the morphological operation.
596 % Plus "[{radius}[,{scale}]]"
597 % Generate a kernel in the shape of a 'plus' sign. The length of each
598 % arm is also the radius, which defaults to 2.
600 % This kernel is not a good general morphological kernel, but is used
601 % more for highlighting and marking any single pixels in an image using,
602 % a "Dilate" or "Erode" method as appropriate.
604 % NOTE: "plus:1" is equivelent to a "Diamond" kernel.
606 % Note that unlike other kernels iterating a plus does not produce the
607 % same result as using a larger radius for the cross.
609 % Distance Measuring Kernels
611 % Chebyshev "[{radius}][x{scale}[%!]]"
612 % Manhatten "[{radius}][x{scale}[%!]]"
613 % Euclidean "[{radius}][x{scale}[%!]]"
615 % Different types of distance measuring methods, which are used with the
616 % a 'Distance' morphology method for generating a gradient based on
617 % distance from an edge of a binary shape, though there is a technique
618 % for handling a anti-aliased shape.
620 % Chebyshev Distance (also known as Tchebychev Distance) is a value of
621 % one to any neighbour, orthogonal or diagonal. One why of thinking of
622 % it is the number of squares a 'King' or 'Queen' in chess needs to
623 % traverse reach any other position on a chess board. It results in a
624 % 'square' like distance function, but one where diagonals are closer
627 % Manhatten Distance (also known as Rectilinear Distance, or the Taxi
628 % Cab metric), is the distance needed when you can only travel in
629 % orthogonal (horizontal or vertical) only. It is the distance a 'Rook'
630 % in chess would travel. It results in a diamond like distances, where
631 % diagonals are further than expected.
633 % Euclidean Distance is the 'direct' or 'as the crow flys distance.
634 % However by default the kernel size only has a radius of 1, which
635 % limits the distance to 'Knight' like moves, with only orthogonal and
636 % diagonal measurements being correct. As such for the default kernel
637 % you will get octagonal like distance function, which is reasonally
640 % However if you use a larger radius such as "Euclidean:4" you will
641 % get a much smoother distance gradient from the edge of the shape.
642 % Of course a larger kernel is slower to use, and generally not needed.
644 % To allow the use of fractional distances that you get with diagonals
645 % the actual distance is scaled by a fixed value which the user can
646 % provide. This is not actually nessary for either ""Chebyshev" or
647 % "Manhatten" distance kernels, but is done for all three distance
648 % kernels. If no scale is provided it is set to a value of 100,
649 % allowing for a maximum distance measurement of 655 pixels using a Q16
650 % version of IM, from any edge. However for small images this can
651 % result in quite a dark gradient.
653 % See the 'Distance' Morphological Method, for information of how it is
656 % # Hit-n-Miss Kernel-Lists -- Still to be implemented
658 % # specifically for Pruning, Thinning, Thickening
662 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
663 const GeometryInfo *args)
676 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
678 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
679 if (kernel == (KernelInfo *) NULL)
681 (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
682 kernel->minimum = kernel->maximum = 0.0;
683 kernel->negative_range = kernel->positive_range = 0.0;
685 kernel->next = (KernelInfo *) NULL;
686 kernel->signature = MagickSignature;
689 /* Convolution Kernels */
692 sigma = fabs(args->sigma);
694 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
696 kernel->width = kernel->height =
697 GetOptimalKernelWidth2D(args->rho,sigma);
698 kernel->x = kernel->y = (long) (kernel->width-1)/2;
699 kernel->negative_range = kernel->positive_range = 0.0;
700 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
701 kernel->height*sizeof(double));
702 if (kernel->values == (double *) NULL)
703 return(DestroyKernelInfo(kernel));
705 sigma = 2.0*sigma*sigma; /* simplify the expression */
706 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
707 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
708 kernel->positive_range += (
710 exp(-((double)(u*u+v*v))/sigma)
711 /* / (MagickPI*sigma) */ );
713 kernel->maximum = kernel->values[
714 kernel->y*kernel->width+kernel->x ];
716 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
722 sigma = fabs(args->sigma);
724 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
726 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
727 kernel->x = (long) (kernel->width-1)/2;
730 kernel->negative_range = kernel->positive_range = 0.0;
731 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
732 kernel->height*sizeof(double));
733 if (kernel->values == (double *) NULL)
734 return(DestroyKernelInfo(kernel));
738 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
739 ** It generates a gaussian 3 times the width, and compresses it into
740 ** the expected range. This produces a closer normalization of the
741 ** resulting kernel, especially for very low sigma values.
742 ** As such while wierd it is prefered.
744 ** I am told this method originally came from Photoshop.
746 sigma *= KernelRank; /* simplify expanded curve */
747 v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
748 (void) ResetMagickMemory(kernel->values,0, (size_t)
749 kernel->width*sizeof(double));
750 for ( u=-v; u <= v; u++) {
751 kernel->values[(u+v)/KernelRank] +=
752 exp(-((double)(u*u))/(2.0*sigma*sigma))
753 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
755 for (i=0; i < (long) kernel->width; i++)
756 kernel->positive_range += kernel->values[i];
758 for ( i=0, u=-kernel->x; i < kernel->width; i++, u++)
759 kernel->positive_range += (
761 exp(-((double)(u*u))/(2.0*sigma*sigma))
762 /* / (MagickSQ2PI*sigma) */ );
765 kernel->maximum = kernel->values[ kernel->x ];
766 /* Note that neither methods above generate a normalized kernel,
767 ** though it gets close. The kernel may be 'clipped' by a user defined
768 ** radius, producing a smaller (darker) kernel. Also for very small
769 ** sigma's (> 0.1) the central value becomes larger than one, and thus
770 ** producing a very bright kernel.
773 /* Normalize the 1D Gaussian Kernel
775 ** Because of this the divisor in the above kernel generator is
776 ** not needed, so is not done above.
778 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
780 /* rotate the kernel by given angle */
781 RotateKernelInfo(kernel, args->xi);
786 sigma = fabs(args->sigma);
788 sigma = (sigma <= MagickEpsilon) ? 1.0 : sigma;
790 if ( args->rho < 1.0 )
791 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
793 kernel->width = (unsigned long)args->rho;
794 kernel->x = kernel->y = 0;
796 kernel->negative_range = kernel->positive_range = 0.0;
797 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
798 kernel->height*sizeof(double));
799 if (kernel->values == (double *) NULL)
800 return(DestroyKernelInfo(kernel));
802 /* A comet blur is half a gaussian curve, so that the object is
803 ** blurred in one direction only. This may not be quite the right
804 ** curve so may change in the future. The function must be normalised.
808 sigma *= KernelRank; /* simplify expanded curve */
809 v = (long) kernel->width*KernelRank; /* start/end points to fit range */
810 (void) ResetMagickMemory(kernel->values,0, (size_t)
811 kernel->width*sizeof(double));
812 for ( u=0; u < v; u++) {
813 kernel->values[u/KernelRank] +=
814 exp(-((double)(u*u))/(2.0*sigma*sigma))
815 /* / (MagickSQ2PI*sigma/KernelRank) */ ;
817 for (i=0; i < (long) kernel->width; i++)
818 kernel->positive_range += kernel->values[i];
820 for ( i=0; i < (long) kernel->width; i++)
821 kernel->positive_range += (
823 exp(-((double)(i*i))/(2.0*sigma*sigma))
824 /* / (MagickSQ2PI*sigma) */ );
827 kernel->maximum = kernel->values[0];
829 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
830 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
833 /* Boolean Kernels */
834 case RectangleKernel:
838 if ( type == SquareKernel )
841 kernel->width = kernel->height = 3; /* default radius = 1 */
843 kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
844 kernel->x = kernel->y = (long) (kernel->width-1)/2;
848 /* NOTE: user defaults set in "AcquireKernelInfo()" */
849 if ( args->rho < 1.0 || args->sigma < 1.0 )
850 return(DestroyKernelInfo(kernel)); /* invalid args given */
851 kernel->width = (unsigned long)args->rho;
852 kernel->height = (unsigned long)args->sigma;
853 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
854 args->psi < 0.0 || args->psi > (double)kernel->height )
855 return(DestroyKernelInfo(kernel)); /* invalid args given */
856 kernel->x = (long) args->xi;
857 kernel->y = (long) args->psi;
860 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
861 kernel->height*sizeof(double));
862 if (kernel->values == (double *) NULL)
863 return(DestroyKernelInfo(kernel));
865 /* set all kernel values to 1.0 */
866 u=(long) kernel->width*kernel->height;
867 for ( i=0; i < u; i++)
868 kernel->values[i] = scale;
869 kernel->minimum = kernel->maximum = scale; /* a flat shape */
870 kernel->positive_range = scale*u;
876 kernel->width = kernel->height = 3; /* default radius = 1 */
878 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
879 kernel->x = kernel->y = (long) (kernel->width-1)/2;
881 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
882 kernel->height*sizeof(double));
883 if (kernel->values == (double *) NULL)
884 return(DestroyKernelInfo(kernel));
886 /* set all kernel values within diamond area to scale given */
887 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
888 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
889 if ((labs(u)+labs(v)) <= (long)kernel->x)
890 kernel->positive_range += kernel->values[i] = args->sigma;
892 kernel->values[i] = nan;
893 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
901 limit = (long)(args->rho*args->rho);
902 if (args->rho < 0.1) /* default radius approx 3.5 */
903 kernel->width = kernel->height = 7L, limit = 10L;
905 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
906 kernel->x = kernel->y = (long) (kernel->width-1)/2;
908 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
909 kernel->height*sizeof(double));
910 if (kernel->values == (double *) NULL)
911 return(DestroyKernelInfo(kernel));
913 /* set all kernel values within disk area to 1.0 */
914 for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
915 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
916 if ((u*u+v*v) <= limit)
917 kernel->positive_range += kernel->values[i] = args->sigma;
919 kernel->values[i] = nan;
920 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
926 kernel->width = kernel->height = 5; /* default radius 2 */
928 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
929 kernel->x = kernel->y = (long) (kernel->width-1)/2;
931 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
932 kernel->height*sizeof(double));
933 if (kernel->values == (double *) NULL)
934 return(DestroyKernelInfo(kernel));
936 /* set all kernel values along axises to 1.0 */
937 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
938 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
939 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
940 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
941 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
944 /* Distance Measuring Kernels */
945 case ChebyshevKernel:
948 kernel->width = kernel->height = 3; /* default radius = 1 */
950 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
951 kernel->x = kernel->y = (long) (kernel->width-1)/2;
953 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
954 kernel->height*sizeof(double));
955 if (kernel->values == (double *) NULL)
956 return(DestroyKernelInfo(kernel));
958 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
959 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
960 kernel->positive_range += ( kernel->values[i] =
961 args->sigma*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
962 kernel->maximum = kernel->values[0];
965 case ManhattenKernel:
968 kernel->width = kernel->height = 3; /* default radius = 1 */
970 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
971 kernel->x = kernel->y = (long) (kernel->width-1)/2;
973 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
974 kernel->height*sizeof(double));
975 if (kernel->values == (double *) NULL)
976 return(DestroyKernelInfo(kernel));
978 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
979 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
980 kernel->positive_range += ( kernel->values[i] =
981 args->sigma*(labs(u)+labs(v)) );
982 kernel->maximum = kernel->values[0];
985 case EuclideanKernel:
988 kernel->width = kernel->height = 3; /* default radius = 1 */
990 kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
991 kernel->x = kernel->y = (long) (kernel->width-1)/2;
993 kernel->values=(double *) AcquireQuantumMemory(kernel->width,
994 kernel->height*sizeof(double));
995 if (kernel->values == (double *) NULL)
996 return(DestroyKernelInfo(kernel));
998 for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
999 for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
1000 kernel->positive_range += ( kernel->values[i] =
1001 args->sigma*sqrt((double)(u*u+v*v)) );
1002 kernel->maximum = kernel->values[0];
1005 /* Undefined Kernels */
1006 case LaplacianKernel:
1009 perror("Kernel Type has not been defined yet");
1012 /* Generate a No-Op minimal kernel - 1x1 pixel */
1013 kernel->values=(double *)AcquireQuantumMemory((size_t)1,sizeof(double));
1014 if (kernel->values == (double *) NULL)
1015 return(DestroyKernelInfo(kernel));
1016 kernel->width = kernel->height = 1;
1017 kernel->x = kernel->x = 0;
1018 kernel->type = UndefinedKernel;
1020 kernel->positive_range =
1021 kernel->values[0] = 1.0; /* a flat single-point no-op kernel! */
1029 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1033 % C l o n e K e r n e l I n f o %
1037 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1039 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
1040 % can be modified without effecting the original. The cloned kernel should
1041 % be destroyed using DestoryKernelInfo() when no longer needed.
1043 % The format of the CloneKernelInfo method is:
1045 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1047 % A description of each parameter follows:
1049 % o kernel: the Morphology/Convolution kernel to be cloned
1052 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
1060 assert(kernel != (KernelInfo *) NULL);
1061 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1062 if (new_kernel == (KernelInfo *) NULL)
1064 *new_kernel=(*kernel); /* copy values in structure */
1066 /* replace the values with a copy of the values */
1067 new_kernel->values=(double *) AcquireQuantumMemory(kernel->width,
1068 kernel->height*sizeof(double));
1069 if (new_kernel->values == (double *) NULL)
1070 return(DestroyKernelInfo(new_kernel));
1071 for (i=0; i < (long) (kernel->width*kernel->height); i++)
1072 new_kernel->values[i]=kernel->values[i];
1074 /* Also clone the next kernel in the kernel list */
1075 if ( kernel->next != (KernelInfo *) NULL ) {
1076 new_kernel->next = CloneKernelInfo(kernel->next);
1077 if ( new_kernel->next == (KernelInfo *) NULL )
1078 return(DestroyKernelInfo(new_kernel));
1085 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1089 % D e s t r o y K e r n e l I n f o %
1093 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1095 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
1098 % The format of the DestroyKernelInfo method is:
1100 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1102 % A description of each parameter follows:
1104 % o kernel: the Morphology/Convolution kernel to be destroyed
1108 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
1110 assert(kernel != (KernelInfo *) NULL);
1112 if ( kernel->next != (KernelInfo *) NULL )
1113 kernel->next = DestroyKernelInfo(kernel->next);
1115 kernel->values = (double *)RelinquishMagickMemory(kernel->values);
1116 kernel = (KernelInfo *) RelinquishMagickMemory(kernel);
1121 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1125 % M o r p h o l o g y I m a g e C h a n n e l %
1129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1131 % MorphologyImageChannel() applies a user supplied kernel to the image
1132 % according to the given mophology method.
1134 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1135 % by the kernel generator.
1137 % The format of the MorphologyImage method is:
1139 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
1140 % const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
1141 % Image *MorphologyImageChannel(const Image *image, const ChannelType
1142 % channel,MorphologyMethod method,const long iterations,
1143 % KernelInfo *kernel,ExceptionInfo *exception)
1145 % A description of each parameter follows:
1147 % o image: the image.
1149 % o method: the morphology method to be applied.
1151 % o iterations: apply the operation this many times (or no change).
1152 % A value of -1 means loop until no change found.
1153 % How this is applied may depend on the morphology method.
1154 % Typically this is a value of 1.
1156 % o channel: the channel type.
1158 % o kernel: An array of double representing the morphology kernel.
1159 % Warning: kernel may be normalized for the Convolve method.
1161 % o exception: return any errors or warnings in this structure.
1164 % TODO: bias and auto-scale handling of the kernel for convolution
1165 % The given kernel is assumed to have been pre-scaled appropriatally, usally
1166 % by the kernel generator.
1171 /* Internal function
1172 * Apply the low-level Morphology Method Primatives using the given Kernel
1173 * Returning the number of pixels that changed.
1174 * Two pre-created images must be provided, no image is created.
1176 static unsigned long MorphologyPrimative(const Image *image, Image
1177 *result_image, const MorphologyMethod method, const ChannelType channel,
1178 const KernelInfo *kernel, ExceptionInfo *exception)
1180 #define MorphologyTag "Morphology/Image"
1197 /* Only the most basic morphology is actually performed by this routine */
1200 Apply Basic Morphology to Image.
1206 GetMagickPixelPacket(image,&bias);
1207 SetMagickPixelPacketBias(image,&bias);
1208 /* Future: handle auto-bias from user, based on kernel input */
1210 p_view=AcquireCacheView(image);
1211 q_view=AcquireCacheView(result_image);
1213 /* Some methods (including convolve) needs use a reflected kernel.
1214 * Adjust 'origin' offsets for this reflected kernel.
1219 case ConvolveMorphology:
1220 case DilateMorphology:
1221 case DilateIntensityMorphology:
1222 case DistanceMorphology:
1223 /* kernel needs to used with reflection about origin */
1224 offx = (long) kernel->width-offx-1;
1225 offy = (long) kernel->height-offy-1;
1227 case ErodeMorphology:
1228 case ErodeIntensityMorphology:
1229 case HitAndMissMorphology:
1230 case ThinningMorphology:
1231 case ThickenMorphology:
1232 /* kernel is user as is, without reflection */
1235 perror("Not a low level Morphology Method");
1239 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1240 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1242 for (y=0; y < (long) image->rows; y++)
1247 register const PixelPacket
1250 register const IndexPacket
1251 *restrict p_indexes;
1253 register PixelPacket
1256 register IndexPacket
1257 *restrict q_indexes;
1265 if (status == MagickFalse)
1267 p=GetCacheViewVirtualPixels(p_view, -offx, y-offy,
1268 image->columns+kernel->width, kernel->height, exception);
1269 q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
1271 if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
1276 p_indexes=GetCacheViewVirtualIndexQueue(p_view);
1277 q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
1278 r = (image->columns+kernel->width)*offy+offx; /* constant */
1280 for (x=0; x < (long) image->columns; x++)
1288 register const double
1291 register const PixelPacket
1294 register const IndexPacket
1295 *restrict k_indexes;
1302 /* Copy input to ouput image for unused channels
1303 * This removes need for 'cloning' a new image every iteration
1306 if (image->colorspace == CMYKColorspace)
1307 q_indexes[x] = p_indexes[r];
1314 min.index = (MagickRealType) QuantumRange;
1319 max.index = (MagickRealType) 0;
1320 /* original pixel value */
1321 result.red = (MagickRealType) p[r].red;
1322 result.green = (MagickRealType) p[r].green;
1323 result.blue = (MagickRealType) p[r].blue;
1324 result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
1326 if ( image->colorspace == CMYKColorspace)
1327 result.index = (MagickRealType) p_indexes[r];
1330 case ConvolveMorphology:
1331 /* Set the user defined bias of the weighted average output
1333 ** FUTURE: provide some way for internal functions to disable
1334 ** user provided bias and scaling effects.
1338 case DilateIntensityMorphology:
1339 case ErodeIntensityMorphology:
1340 result.red = 0.0; /* flag indicating when first match found */
1347 case ConvolveMorphology:
1348 /* Weighted Average of pixels using reflected kernel
1350 ** NOTE for correct working of this operation for asymetrical
1351 ** kernels, the kernel needs to be applied in its reflected form.
1352 ** That is its values needs to be reversed.
1354 ** Correlation is actually the same as this but without reflecting
1355 ** the kernel, and thus 'lower-level' that Convolution. However
1356 ** as Convolution is the more common method used, and it does not
1357 ** really cost us much in terms of processing to use a reflected
1358 ** kernel, so it is Convolution that is implemented.
1360 ** Correlation will have its kernel reflected before calling
1361 ** this function to do a Convolve.
1363 ** For more details of Correlation vs Convolution see
1364 ** http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
1366 if (((channel & SyncChannels) != 0 ) &&
1367 (image->matte == MagickTrue))
1368 { /* Channel has a 'Sync' Flag, and Alpha Channel enabled.
1369 ** Weight the color channels with Alpha Channel so that
1370 ** transparent pixels are not part of the results.
1373 alpha, /* color channel weighting : kernel*alpha */
1374 gamma; /* divisor, sum of weighting values */
1377 k = &kernel->values[ kernel->width*kernel->height-1 ];
1379 k_indexes = p_indexes;
1380 for (v=0; v < (long) kernel->height; v++) {
1381 for (u=0; u < (long) kernel->width; u++, k--) {
1382 if ( IsNan(*k) ) continue;
1383 alpha=(*k)*(QuantumScale*(QuantumRange-
1384 k_pixels[u].opacity));
1386 result.red += alpha*k_pixels[u].red;
1387 result.green += alpha*k_pixels[u].green;
1388 result.blue += alpha*k_pixels[u].blue;
1389 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1390 if ( image->colorspace == CMYKColorspace)
1391 result.index += alpha*k_indexes[u];
1393 k_pixels += image->columns+kernel->width;
1394 k_indexes += image->columns+kernel->width;
1396 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
1397 result.red *= gamma;
1398 result.green *= gamma;
1399 result.blue *= gamma;
1400 result.opacity *= gamma;
1401 result.index *= gamma;
1405 /* No 'Sync' flag, or no Alpha involved.
1406 ** Convolution is simple individual channel weigthed sum.
1408 k = &kernel->values[ kernel->width*kernel->height-1 ];
1410 k_indexes = p_indexes;
1411 for (v=0; v < (long) kernel->height; v++) {
1412 for (u=0; u < (long) kernel->width; u++, k--) {
1413 if ( IsNan(*k) ) continue;
1414 result.red += (*k)*k_pixels[u].red;
1415 result.green += (*k)*k_pixels[u].green;
1416 result.blue += (*k)*k_pixels[u].blue;
1417 result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
1418 if ( image->colorspace == CMYKColorspace)
1419 result.index += (*k)*k_indexes[u];
1421 k_pixels += image->columns+kernel->width;
1422 k_indexes += image->columns+kernel->width;
1427 case ErodeMorphology:
1428 /* Minimum Value within kernel neighbourhood
1430 ** NOTE that the kernel is not reflected for this operation!
1432 ** NOTE: in normal Greyscale Morphology, the kernel value should
1433 ** be added to the real value, this is currently not done, due to
1434 ** the nature of the boolean kernels being used.
1438 k_indexes = p_indexes;
1439 for (v=0; v < (long) kernel->height; v++) {
1440 for (u=0; u < (long) kernel->width; u++, k++) {
1441 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1442 Minimize(min.red, (double) k_pixels[u].red);
1443 Minimize(min.green, (double) k_pixels[u].green);
1444 Minimize(min.blue, (double) k_pixels[u].blue);
1445 Minimize(min.opacity,
1446 QuantumRange-(double) k_pixels[u].opacity);
1447 if ( image->colorspace == CMYKColorspace)
1448 Minimize(min.index, (double) k_indexes[u]);
1450 k_pixels += image->columns+kernel->width;
1451 k_indexes += image->columns+kernel->width;
1456 case DilateMorphology:
1457 /* Maximum Value within kernel neighbourhood
1459 ** NOTE for correct working of this operation for asymetrical
1460 ** kernels, the kernel needs to be applied in its reflected form.
1461 ** That is its values needs to be reversed.
1463 ** NOTE: in normal Greyscale Morphology, the kernel value should
1464 ** be added to the real value, this is currently not done, due to
1465 ** the nature of the boolean kernels being used.
1468 k = &kernel->values[ kernel->width*kernel->height-1 ];
1470 k_indexes = p_indexes;
1471 for (v=0; v < (long) kernel->height; v++) {
1472 for (u=0; u < (long) kernel->width; u++, k--) {
1473 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1474 Maximize(max.red, (double) k_pixels[u].red);
1475 Maximize(max.green, (double) k_pixels[u].green);
1476 Maximize(max.blue, (double) k_pixels[u].blue);
1477 Maximize(max.opacity,
1478 QuantumRange-(double) k_pixels[u].opacity);
1479 if ( image->colorspace == CMYKColorspace)
1480 Maximize(max.index, (double) k_indexes[u]);
1482 k_pixels += image->columns+kernel->width;
1483 k_indexes += image->columns+kernel->width;
1487 case HitAndMissMorphology:
1488 case ThinningMorphology:
1489 case ThickenMorphology:
1490 /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
1492 ** NOTE that the kernel is not reflected for this operation,
1493 ** and consists of both foreground and background pixel
1494 ** neighbourhoods, 0.0 for background, and 1.0 for foreground
1495 ** with either Nan or 0.5 values for don't care.
1497 ** Note that this can produce negative results, though really
1498 ** only a positive match has any real value.
1502 k_indexes = p_indexes;
1503 for (v=0; v < (long) kernel->height; v++) {
1504 for (u=0; u < (long) kernel->width; u++, k++) {
1505 if ( IsNan(*k) ) continue;
1507 { /* minimim of foreground pixels */
1508 Minimize(min.red, (double) k_pixels[u].red);
1509 Minimize(min.green, (double) k_pixels[u].green);
1510 Minimize(min.blue, (double) k_pixels[u].blue);
1511 Minimize(min.opacity,
1512 QuantumRange-(double) k_pixels[u].opacity);
1513 if ( image->colorspace == CMYKColorspace)
1514 Minimize(min.index, (double) k_indexes[u]);
1516 else if ( (*k) < 0.3 )
1517 { /* maximum of background pixels */
1518 Maximize(max.red, (double) k_pixels[u].red);
1519 Maximize(max.green, (double) k_pixels[u].green);
1520 Maximize(max.blue, (double) k_pixels[u].blue);
1521 Maximize(max.opacity,
1522 QuantumRange-(double) k_pixels[u].opacity);
1523 if ( image->colorspace == CMYKColorspace)
1524 Maximize(max.index, (double) k_indexes[u]);
1527 k_pixels += image->columns+kernel->width;
1528 k_indexes += image->columns+kernel->width;
1530 /* Pattern Match only if min fg larger than min bg pixels */
1531 min.red -= max.red; Maximize( min.red, 0.0 );
1532 min.green -= max.green; Maximize( min.green, 0.0 );
1533 min.blue -= max.blue; Maximize( min.blue, 0.0 );
1534 min.opacity -= max.opacity; Maximize( min.opacity, 0.0 );
1535 min.index -= max.index; Maximize( min.index, 0.0 );
1538 case ErodeIntensityMorphology:
1539 /* Select Pixel with Minimum Intensity within kernel neighbourhood
1541 ** WARNING: the intensity test fails for CMYK and does not
1542 ** take into account the moderating effect of teh alpha channel
1543 ** on the intensity.
1545 ** NOTE that the kernel is not reflected for this operation!
1549 k_indexes = p_indexes;
1550 for (v=0; v < (long) kernel->height; v++) {
1551 for (u=0; u < (long) kernel->width; u++, k++) {
1552 if ( IsNan(*k) || (*k) < 0.5 ) continue;
1553 if ( result.red == 0.0 ||
1554 PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
1555 /* copy the whole pixel - no channel selection */
1557 if ( result.red > 0.0 ) changed++;
1561 k_pixels += image->columns+kernel->width;
1562 k_indexes += image->columns+kernel->width;
1566 case DilateIntensityMorphology:
1567 /* Select Pixel with Maximum Intensity within kernel neighbourhood
1569 ** WARNING: the intensity test fails for CMYK and does not
1570 ** take into account the moderating effect of teh alpha channel
1571 ** on the intensity.
1573 ** NOTE for correct working of this operation for asymetrical
1574 ** kernels, the kernel needs to be applied in its reflected form.
1575 ** That is its values needs to be reversed.
1577 k = &kernel->values[ kernel->width*kernel->height-1 ];
1579 k_indexes = p_indexes;
1580 for (v=0; v < (long) kernel->height; v++) {
1581 for (u=0; u < (long) kernel->width; u++, k--) {
1582 if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
1583 if ( result.red == 0.0 ||
1584 PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
1585 /* copy the whole pixel - no channel selection */
1587 if ( result.red > 0.0 ) changed++;
1591 k_pixels += image->columns+kernel->width;
1592 k_indexes += image->columns+kernel->width;
1597 case DistanceMorphology:
1598 /* Add kernel Value and select the minimum value found.
1599 ** The result is a iterative distance from edge of image shape.
1601 ** All Distance Kernels are symetrical, but that may not always
1602 ** be the case. For example how about a distance from left edges?
1603 ** To work correctly with asymetrical kernels the reflected kernel
1604 ** needs to be applied.
1606 ** Actually this is really a GreyErode with a negative kernel!
1609 k = &kernel->values[ kernel->width*kernel->height-1 ];
1611 k_indexes = p_indexes;
1612 for (v=0; v < (long) kernel->height; v++) {
1613 for (u=0; u < (long) kernel->width; u++, k--) {
1614 if ( IsNan(*k) ) continue;
1615 Minimize(result.red, (*k)+k_pixels[u].red);
1616 Minimize(result.green, (*k)+k_pixels[u].green);
1617 Minimize(result.blue, (*k)+k_pixels[u].blue);
1618 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
1619 if ( image->colorspace == CMYKColorspace)
1620 Minimize(result.index, (*k)+k_indexes[u]);
1622 k_pixels += image->columns+kernel->width;
1623 k_indexes += image->columns+kernel->width;
1627 case UndefinedMorphology:
1629 break; /* Do nothing */
1631 /* Final mathematics of results (combine with original image?)
1633 ** NOTE: Difference Morphology operators Edge* and *Hat could also
1634 ** be done here but works better with iteration as a image difference
1635 ** in the controling function (below). Thicken and Thinning however
1636 ** should be done here so thay can be iterated correctly.
1639 case HitAndMissMorphology:
1640 case ErodeMorphology:
1641 result = min; /* minimum of neighbourhood */
1643 case DilateMorphology:
1644 result = max; /* maximum of neighbourhood */
1646 case ThinningMorphology:
1647 /* subtract pattern match from original */
1648 result.red -= min.red;
1649 result.green -= min.green;
1650 result.blue -= min.blue;
1651 result.opacity -= min.opacity;
1652 result.index -= min.index;
1654 case ThickenMorphology:
1655 /* Union with original image (maximize) - or should this be + */
1656 Maximize( result.red, min.red );
1657 Maximize( result.green, min.green );
1658 Maximize( result.blue, min.blue );
1659 Maximize( result.opacity, min.opacity );
1660 Maximize( result.index, min.index );
1663 /* result directly calculated or assigned */
1666 /* Assign the resulting pixel values - Clamping Result */
1668 case UndefinedMorphology:
1669 case DilateIntensityMorphology:
1670 case ErodeIntensityMorphology:
1671 break; /* full pixel was directly assigned - not a channel method */
1673 if ((channel & RedChannel) != 0)
1674 q->red = ClampToQuantum(result.red);
1675 if ((channel & GreenChannel) != 0)
1676 q->green = ClampToQuantum(result.green);
1677 if ((channel & BlueChannel) != 0)
1678 q->blue = ClampToQuantum(result.blue);
1679 if ((channel & OpacityChannel) != 0
1680 && image->matte == MagickTrue )
1681 q->opacity = ClampToQuantum(QuantumRange-result.opacity);
1682 if ((channel & IndexChannel) != 0
1683 && image->colorspace == CMYKColorspace)
1684 q_indexes[x] = ClampToQuantum(result.index);
1687 /* Count up changed pixels */
1688 if ( ( p[r].red != q->red )
1689 || ( p[r].green != q->green )
1690 || ( p[r].blue != q->blue )
1691 || ( p[r].opacity != q->opacity )
1692 || ( image->colorspace == CMYKColorspace &&
1693 p_indexes[r] != q_indexes[x] ) )
1694 changed++; /* The pixel had some value changed! */
1698 sync=SyncCacheViewAuthenticPixels(q_view,exception);
1699 if (sync == MagickFalse)
1701 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1706 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1707 #pragma omp critical (MagickCore_MorphologyImage)
1709 proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
1710 if (proceed == MagickFalse)
1714 result_image->type=image->type;
1715 q_view=DestroyCacheView(q_view);
1716 p_view=DestroyCacheView(p_view);
1717 return(status ? (unsigned long) changed : 0);
1721 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
1722 method, const long iterations,const KernelInfo *kernel, ExceptionInfo
1728 morphology_image=MorphologyImageChannel(image,DefaultChannels,method,
1729 iterations,kernel,exception);
1730 return(morphology_image);
1734 MagickExport Image *MorphologyImageChannel(const Image *image,
1735 const ChannelType channel,const MorphologyMethod method,
1736 const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
1757 assert(image != (Image *) NULL);
1758 assert(image->signature == MagickSignature);
1759 assert(kernel != (KernelInfo *) NULL);
1760 assert(kernel->signature == MagickSignature);
1761 assert(exception != (ExceptionInfo *) NULL);
1762 assert(exception->signature == MagickSignature);
1764 if ( iterations == 0 )
1765 return((Image *)NULL); /* null operation - nothing to do! */
1767 /* kernel must be valid at this point
1768 * (except maybe for posible future morphology methods like "Prune"
1770 assert(kernel != (KernelInfo *)NULL);
1772 new_image = (Image *) NULL; /* for cleanup */
1773 old_image = (Image *) NULL;
1774 grad_image = (Image *) NULL;
1775 curr_kernel = (KernelInfo *) NULL;
1776 count = 0; /* interation count */
1777 changed = 1; /* was last run succesfull */
1778 curr_kernel = (KernelInfo *)kernel; /* allow kernel and method */
1779 curr_method = method; /* to be changed as nessary */
1781 limit = (unsigned long) iterations;
1782 if ( iterations < 0 )
1783 limit = image->columns > image->rows ? image->columns : image->rows;
1785 /* Third-level morphology methods */
1786 switch( curr_method ) {
1787 case EdgeMorphology:
1788 grad_image = MorphologyImageChannel(image, channel,
1789 DilateMorphology, iterations, curr_kernel, exception);
1790 if ( grad_image == (Image *) NULL )
1793 case EdgeInMorphology:
1794 curr_method = ErodeMorphology;
1796 case EdgeOutMorphology:
1797 curr_method = DilateMorphology;
1799 case TopHatMorphology:
1800 curr_method = OpenMorphology;
1802 case BottomHatMorphology:
1803 curr_method = CloseMorphology;
1806 break; /* not a third-level method */
1809 /* Second-level morphology methods */
1810 switch( curr_method ) {
1811 case OpenMorphology:
1812 /* Open is a Erode then a Dilate without reflection */
1813 new_image = MorphologyImageChannel(image, channel,
1814 ErodeMorphology, iterations, curr_kernel, exception);
1815 if (new_image == (Image *) NULL)
1817 curr_method = DilateMorphology;
1819 case OpenIntensityMorphology:
1820 new_image = MorphologyImageChannel(image, channel,
1821 ErodeIntensityMorphology, iterations, curr_kernel, exception);
1822 if (new_image == (Image *) NULL)
1824 curr_method = DilateIntensityMorphology;
1827 case CloseMorphology:
1828 /* Close is a Dilate then Erode using reflected kernel */
1829 /* A reflected kernel is needed for a Close */
1830 if ( curr_kernel == kernel )
1831 curr_kernel = CloneKernelInfo(kernel);
1832 if (curr_kernel == (KernelInfo *) NULL)
1834 RotateKernelInfo(curr_kernel,180);
1835 new_image = MorphologyImageChannel(image, channel,
1836 DilateMorphology, iterations, curr_kernel, exception);
1837 if (new_image == (Image *) NULL)
1839 curr_method = ErodeMorphology;
1841 case CloseIntensityMorphology:
1842 /* A reflected kernel is needed for a Close */
1843 if ( curr_kernel == kernel )
1844 curr_kernel = CloneKernelInfo(kernel);
1845 if (curr_kernel == (KernelInfo *) NULL)
1847 RotateKernelInfo(curr_kernel,180);
1848 new_image = MorphologyImageChannel(image, channel,
1849 DilateIntensityMorphology, iterations, curr_kernel, exception);
1850 if (new_image == (Image *) NULL)
1852 curr_method = ErodeIntensityMorphology;
1855 case CorrelateMorphology:
1856 /* A Correlation is actually a Convolution with a reflected kernel.
1857 ** However a Convolution is a weighted sum with a reflected kernel.
1858 ** It may seem stange to convert a Correlation into a Convolution
1859 ** as the Correleation is the simplier method, but Convolution is
1860 ** much more commonly used, and it makes sense to implement it directly
1861 ** so as to avoid the need to duplicate the kernel when it is not
1862 ** required (which is typically the default).
1864 if ( curr_kernel == kernel )
1865 curr_kernel = CloneKernelInfo(kernel);
1866 if (curr_kernel == (KernelInfo *) NULL)
1868 RotateKernelInfo(curr_kernel,180);
1869 curr_method = ConvolveMorphology;
1870 /* FALL-THRU into Correlate (weigthed sum without reflection) */
1872 case ConvolveMorphology:
1873 { /* Scale or Normalize kernel, according to user wishes
1874 ** before using it for the Convolve/Correlate method.
1876 ** FUTURE: provide some way for internal functions to disable
1877 ** user bias and scaling effects.
1880 *artifact = GetImageArtifact(image,"convolve:scale");
1881 if ( artifact != (char *)NULL ) {
1887 if ( curr_kernel == kernel )
1888 curr_kernel = CloneKernelInfo(kernel);
1889 if (curr_kernel == (KernelInfo *) NULL)
1892 flags = (GeometryFlags) ParseGeometry(artifact, &args);
1893 ScaleKernelInfo(curr_kernel, args.rho, flags);
1896 /* FALL-THRU to do the first, and typically the only iteration */
1899 /* Do a single iteration using the Low-Level Morphology method!
1900 ** This ensures a "new_image" has been generated, but allows us to skip
1901 ** the creation of 'old_image' if no more iterations are needed.
1903 ** The "curr_method" should also be set to a low-level method that is
1904 ** understood by the MorphologyPrimative() internal function.
1906 new_image=CloneImage(image,0,0,MagickTrue,exception);
1907 if (new_image == (Image *) NULL)
1909 if (SetImageStorageClass(new_image,DirectClass) == MagickFalse)
1911 InheritException(exception,&new_image->exception);
1914 count++; /* iteration count */
1915 changed = MorphologyPrimative(image,new_image,curr_method,channel,
1916 curr_kernel, exception);
1917 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1918 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
1919 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1920 count, 0L, changed);
1924 * + "curr_method" should be set to a low-level morphology method
1925 * + "count=1" if the first iteration of the first kernel has been done.
1926 * + "new_image" is defined and contains the current resulting image
1929 /* The Hit-And-Miss Morphology is not 'iterated' against the resulting
1930 ** image from the previous morphology step. It must always be applied
1931 ** to the original image, with the results collected into a union (maximimum
1932 ** or lighten) of all the results. Multiple kernels may be applied but
1933 ** an iteration of the morphology does nothing, so is ignored.
1935 ** The first kernel is guranteed to have been done, so lets continue the
1936 ** process, with the other kernels in the kernel list.
1938 if ( method == HitAndMissMorphology )
1940 if ( curr_kernel->next != (KernelInfo *) NULL ) {
1941 /* create a second working image */
1942 old_image = CloneImage(image,0,0,MagickTrue,exception);
1943 if (old_image == (Image *) NULL)
1945 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1947 InheritException(exception,&old_image->exception);
1951 /* loop through rest of the kernels */
1952 this_kernel=curr_kernel->next;
1954 while( this_kernel != (KernelInfo *) NULL )
1956 changed = MorphologyPrimative(image,old_image,curr_method,channel,
1957 this_kernel,exception);
1958 (void) CompositeImageChannel(new_image,
1959 (ChannelType) (channel & ~SyncChannels), LightenCompositeOp,
1961 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
1962 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
1963 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
1964 count, kernel_number, changed);
1965 this_kernel = this_kernel->next;
1968 old_image=DestroyImage(old_image);
1973 /* Repeat the low-level morphology over all kernels
1974 until iteration count limit or no change from any kernel is found */
1975 if ( ( count < limit && changed > 0 ) ||
1976 curr_kernel->next != (KernelInfo *) NULL ) {
1978 /* create a second working image */
1979 old_image = CloneImage(image,0,0,MagickTrue,exception);
1980 if (old_image == (Image *) NULL)
1982 if (SetImageStorageClass(old_image,DirectClass) == MagickFalse)
1984 InheritException(exception,&old_image->exception);
1988 /* reset variables for the first/next iteration, or next kernel) */
1990 this_kernel = curr_kernel;
1991 total_changed = count != 0 ? changed : 0;
1992 if ( count != 0 && this_kernel != (KernelInfo *) NULL ) {
1993 count = 0; /* first iteration is not yet finished! */
1994 this_kernel = curr_kernel->next;
1996 total_changed = changed;
1999 while ( count < limit ) {
2001 while ( this_kernel != (KernelInfo *) NULL ) {
2002 Image *tmp = old_image;
2003 old_image = new_image;
2005 changed = MorphologyPrimative(old_image,new_image,curr_method,channel,
2006 this_kernel,exception);
2007 if ( GetImageArtifact(image,"verbose") != (const char *) NULL )
2008 fprintf(stderr, "Morphology %s:%lu.%lu => Changed %lu\n",
2009 MagickOptionToMnemonic(MagickMorphologyOptions, curr_method),
2010 count, kernel_number, changed);
2011 total_changed += changed;
2012 this_kernel = this_kernel->next;
2015 if ( total_changed == 0 )
2016 break; /* no changes after processing all kernels - ABORT */
2017 /* prepare for next loop */
2020 this_kernel = curr_kernel;
2022 old_image=DestroyImage(old_image);
2025 /* finished with kernel - destary any copy that was made */
2026 if ( curr_kernel != kernel )
2027 curr_kernel=DestroyKernelInfo(curr_kernel);
2029 /* Third-level Subtractive methods post-processing
2031 ** The removal of any 'Sync' channel flag in the Image Compositon below
2032 ** ensures the compose method is applied in a purely mathematical way, only
2033 ** the selected channels, without any normal 'alpha blending' normally
2034 ** associated with the compose method.
2036 ** Note "method" here is the 'original' morphological method, and not the
2037 ** 'current' morphological method used above to generate "new_image".
2040 case EdgeOutMorphology:
2041 case EdgeInMorphology:
2042 case TopHatMorphology:
2043 case BottomHatMorphology:
2044 /* Get Difference relative to the original image */
2045 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2046 DifferenceCompositeOp, image, 0, 0);
2048 case EdgeMorphology:
2049 /* Difference the Eroded image from the saved Dilated image */
2050 (void) CompositeImageChannel(new_image, (ChannelType) (channel & ~SyncChannels),
2051 DifferenceCompositeOp, grad_image, 0, 0);
2052 grad_image=DestroyImage(grad_image);
2060 /* Yes goto's are bad, but in this case it makes cleanup lot more efficient */
2062 if ( new_image != (Image *) NULL )
2063 DestroyImage(new_image);
2065 if ( old_image != (Image *) NULL )
2066 DestroyImage(old_image);
2067 if ( curr_kernel != (KernelInfo *) NULL && curr_kernel != kernel )
2068 curr_kernel=DestroyKernelInfo(curr_kernel);
2073 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2077 + R o t a t e K e r n e l I n f o %
2081 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2083 % RotateKernelInfo() rotates the kernel by the angle given. Currently it is
2084 % restricted to 90 degree angles, but this may be improved in the future.
2086 % The format of the RotateKernelInfo method is:
2088 % void RotateKernelInfo(KernelInfo *kernel, double angle)
2090 % A description of each parameter follows:
2092 % o kernel: the Morphology/Convolution kernel
2094 % o angle: angle to rotate in degrees
2096 % This function is only internel to this module, as it is not finalized,
2097 % especially with regard to non-orthogonal angles, and rotation of larger
2100 static void RotateKernelInfo(KernelInfo *kernel, double angle)
2102 /* angle the lower kernels first */
2103 if ( kernel->next != (KernelInfo *) NULL)
2104 RotateKernelInfo(kernel->next, angle);
2106 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
2108 ** TODO: expand beyond simple 90 degree rotates, flips and flops
2111 /* Modulus the angle */
2112 angle = fmod(angle, 360.0);
2116 if ( 315.0 < angle || angle <= 45.0 )
2117 return; /* no change! - At least at this time */
2119 switch (kernel->type) {
2120 /* These built-in kernels are cylindrical kernels, rotating is useless */
2121 case GaussianKernel:
2122 case LaplacianKernel:
2126 case ChebyshevKernel:
2127 case ManhattenKernel:
2128 case EuclideanKernel:
2131 /* These may be rotatable at non-90 angles in the future */
2132 /* but simply rotating them in multiples of 90 degrees is useless */
2138 /* These only allows a +/-90 degree rotation (by transpose) */
2139 /* A 180 degree rotation is useless */
2141 case RectangleKernel:
2142 if ( 135.0 < angle && angle <= 225.0 )
2144 if ( 225.0 < angle && angle <= 315.0 )
2148 /* these are freely rotatable in 90 degree units */
2150 case UndefinedKernel:
2151 case UserDefinedKernel:
2154 if ( 135.0 < angle && angle <= 225.0 )
2156 /* For a 180 degree rotation - also know as a reflection */
2157 /* This is actually a very very common operation! */
2158 /* Basically all that is needed is a reversal of the kernel data! */
2165 for ( i=0, j=kernel->width*kernel->height-1; i<j; i++, j--)
2166 t=k[i], k[i]=k[j], k[j]=t;
2168 kernel->x = (long) kernel->width - kernel->x - 1;
2169 kernel->y = (long) kernel->height - kernel->y - 1;
2170 angle = fmod(angle+180.0, 360.0);
2172 if ( 45.0 < angle && angle <= 135.0 )
2173 { /* Do a transpose and a flop, of the image, which results in a 90
2174 * degree rotation using two mirror operations.
2176 * WARNING: this assumes the original image was a 1 dimentional image
2177 * but currently that is the only built-ins it is applied to.
2181 t = (long) kernel->width;
2182 kernel->width = kernel->height;
2183 kernel->height = (unsigned long) t;
2185 kernel->x = kernel->y;
2187 angle = fmod(450.0 - angle, 360.0);
2189 /* At this point angle should be between -45 (315) and +45 degrees
2190 * In the future some form of non-orthogonal angled rotates could be
2191 * performed here, posibily with a linear kernel restriction.
2195 Not currently in use!
2196 { /* Do a flop, this assumes kernel is horizontally symetrical.
2197 * Each row of the kernel needs to be reversed!
2206 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2207 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2208 t=k[x], k[x]=k[r], k[r]=t;
2210 kernel->x = kernel->width - kernel->x - 1;
2211 angle = fmod(angle+180.0, 360.0);
2218 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2222 % S c a l e K e r n e l I n f o %
2226 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2228 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
2229 % without normalization of the sum of the kernel values (as per given flags).
2231 % By default (no flags given) the values within the kernel is scaled
2232 % directly using given scaling factor without change.
2234 % If any 'normalize_flags' are given the kernel will first be normalized and
2235 % then further scaled by the scaling factor value given. A 'PercentValue'
2236 % flag will cause the given scaling factor to be divided by one hundred
2239 % Kernel normalization ('normalize_flags' given) is designed to ensure that
2240 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
2241 % morphology methods will fall into -1.0 to +1.0 range. Note that for
2242 % non-HDRI versions of IM this may cause images to have any negative results
2243 % clipped, unless some 'bias' is used.
2245 % More specifically. Kernels which only contain positive values (such as a
2246 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
2247 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
2249 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
2250 % the kernel will be scaled by the absolute of the sum of kernel values, so
2251 % that it will generally fall within the +/- 1.0 range.
2253 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
2254 % will be scaled by just the sum of the postive values, so that its output
2255 % range will again fall into the +/- 1.0 range.
2257 % For special kernels designed for locating shapes using 'Correlate', (often
2258 % only containing +1 and -1 values, representing foreground/brackground
2259 % matching) a special normalization method is provided to scale the positive
2260 % values seperatally to those of the negative values, so the kernel will be
2261 % forced to become a zero-sum kernel better suited to such searches.
2263 % WARNING: Correct normalization of the kernel assumes that the '*_range'
2264 % attributes within the kernel structure have been correctly set during the
2267 % NOTE: The values used for 'normalize_flags' have been selected specifically
2268 % to match the use of geometry options, so that '!' means NormalizeValue,
2269 % '^' means CorrelateNormalizeValue, and '%' means PercentValue. All other
2270 % GeometryFlags values are ignored.
2272 % The format of the ScaleKernelInfo method is:
2274 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
2275 % const MagickStatusType normalize_flags )
2277 % A description of each parameter follows:
2279 % o kernel: the Morphology/Convolution kernel
2282 % multiply all values (after normalization) by this factor if not
2283 % zero. If the kernel is normalized regardless of any flags.
2285 % o normalize_flags:
2286 % GeometryFlags defining normalization method to use.
2287 % specifically: NormalizeValue, CorrelateNormalizeValue,
2288 % and/or PercentValue
2290 % This function is internal to this module only at this time, but can be
2291 % exported to other modules if needed.
2293 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
2294 const double scaling_factor,const GeometryFlags normalize_flags)
2303 /* scale the lower kernels first */
2304 if ( kernel->next != (KernelInfo *) NULL)
2305 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
2308 if ( (normalize_flags&NormalizeValue) != 0 ) {
2309 /* normalize kernel appropriately */
2310 if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
2311 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
2313 pos_scale = kernel->positive_range; /* special zero-summing kernel */
2315 /* force kernel into being a normalized zero-summing kernel */
2316 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
2317 pos_scale = ( fabs(kernel->positive_range) > MagickEpsilon )
2318 ? kernel->positive_range : 1.0;
2319 neg_scale = ( fabs(kernel->negative_range) > MagickEpsilon )
2320 ? -kernel->negative_range : 1.0;
2323 neg_scale = pos_scale;
2325 /* finialize scaling_factor for positive and negative components */
2326 pos_scale = scaling_factor/pos_scale;
2327 neg_scale = scaling_factor/neg_scale;
2328 if ( (normalize_flags&PercentValue) != 0 ) {
2333 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2334 if ( ! IsNan(kernel->values[i]) )
2335 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
2337 /* convolution output range */
2338 kernel->positive_range *= pos_scale;
2339 kernel->negative_range *= neg_scale;
2340 /* maximum and minimum values in kernel */
2341 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
2342 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
2344 /* swap kernel settings if user scaling factor is negative */
2345 if ( scaling_factor < MagickEpsilon ) {
2347 t = kernel->positive_range;
2348 kernel->positive_range = kernel->negative_range;
2349 kernel->negative_range = t;
2350 t = kernel->maximum;
2351 kernel->maximum = kernel->minimum;
2352 kernel->minimum = 1;
2359 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2363 + S h o w K e r n e l I n f o %
2367 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2369 % ShowKernelInfo() outputs the details of the given kernel defination to
2370 % standard error, generally due to a users 'showkernel' option request.
2372 % The format of the ShowKernel method is:
2374 % void ShowKernelInfo(KernelInfo *kernel)
2376 % A description of each parameter follows:
2378 % o kernel: the Morphology/Convolution kernel
2380 % This function is internal to this module only at this time. That may change
2383 MagickExport void ShowKernelInfo(KernelInfo *kernel)
2391 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
2393 fprintf(stderr, "Kernel ");
2394 if ( kernel->next != (KernelInfo *) NULL )
2395 fprintf(stderr, " #%ld", c );
2396 fprintf(stderr, " \"%s\" of size %lux%lu%+ld%+ld ",
2397 MagickOptionToMnemonic(MagickKernelOptions, k->type),
2398 kernel->width, k->height,
2399 kernel->x, kernel->y );
2401 " with values from %.*lg to %.*lg\n",
2402 GetMagickPrecision(), k->minimum,
2403 GetMagickPrecision(), k->maximum);
2404 fprintf(stderr, "Forming convolution output range from %.*lg to %.*lg%s\n",
2405 GetMagickPrecision(), k->negative_range,
2406 GetMagickPrecision(), k->positive_range,
2407 /*kernel->normalized == MagickTrue ? " (normalized)" : */ "" );
2408 for (i=v=0; v < (long) k->height; v++) {
2409 fprintf(stderr,"%2ld:",v);
2410 for (u=0; u < (long) k->width; u++, i++)
2411 if ( IsNan(k->values[i]) )
2412 fprintf(stderr," %*s", GetMagickPrecision()+2, "nan");
2414 fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
2415 GetMagickPrecision(), k->values[i]);
2416 fprintf(stderr,"\n");
2422 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2426 + Z e r o K e r n e l N a n s %
2430 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2432 % ZeroKernelNans() replaces any special 'nan' value that may be present in
2433 % the kernel with a zero value. This is typically done when the kernel will
2434 % be used in special hardware (GPU) convolution processors, to simply
2437 % The format of the ZeroKernelNans method is:
2439 % voidZeroKernelNans (KernelInfo *kernel)
2441 % A description of each parameter follows:
2443 % o kernel: the Morphology/Convolution kernel
2445 % FUTURE: return the information in a string for API usage.
2447 MagickExport void ZeroKernelNans(KernelInfo *kernel)
2452 /* scale the lower kernels first */
2453 if ( kernel->next != (KernelInfo *) NULL)
2454 ZeroKernelNans(kernel->next);
2456 for (i=0; i < (long) (kernel->width*kernel->height); i++)
2457 if ( IsNan(kernel->values[i]) )
2458 kernel->values[i] = 0.0;