MagickCore 7.1.1
Convert, Edit, Or Compose Bitmap Images
Loading...
Searching...
No Matches
morphology.c
1/*
2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3% %
4% %
5% %
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 %
11% %
12% %
13% MagickCore Morphology Methods %
14% %
15% Software Design %
16% Anthony Thyssen %
17% January 2010 %
18% %
19% %
20% Copyright @ 2010 ImageMagick Studio LLC, a non-profit organization %
21% dedicated to making software imaging solutions freely available. %
22% %
23% You may not use this file except in compliance with the License. You may %
24% obtain a copy of the License at %
25% %
26% https://imagemagick.org/script/license.php %
27% %
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. %
33% %
34%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35%
36% Morphology is the application of various kernels, of any size or shape, to an
37% image in various ways (typically binary, but not always).
38%
39% Convolution (weighted sum or average) is just one specific type of
40% morphology. Just one that is very common for image blurring and sharpening
41% effects. Not only 2D Gaussian blurring, but also 2-pass 1D Blurring.
42%
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.
47*/
48
49/*
50 Include declarations.
51*/
52#include "MagickCore/studio.h"
53#include "MagickCore/artifact.h"
54#include "MagickCore/cache-view.h"
55#include "MagickCore/channel.h"
56#include "MagickCore/color-private.h"
57#include "MagickCore/enhance.h"
58#include "MagickCore/exception.h"
59#include "MagickCore/exception-private.h"
60#include "MagickCore/gem.h"
61#include "MagickCore/gem-private.h"
62#include "MagickCore/image.h"
63#include "MagickCore/image-private.h"
64#include "MagickCore/linked-list.h"
65#include "MagickCore/list.h"
66#include "MagickCore/magick.h"
67#include "MagickCore/memory_.h"
68#include "MagickCore/memory-private.h"
69#include "MagickCore/monitor-private.h"
70#include "MagickCore/morphology.h"
71#include "MagickCore/morphology-private.h"
72#include "MagickCore/option.h"
73#include "MagickCore/pixel-accessor.h"
74#include "MagickCore/prepress.h"
75#include "MagickCore/quantize.h"
76#include "MagickCore/resource_.h"
77#include "MagickCore/registry.h"
78#include "MagickCore/semaphore.h"
79#include "MagickCore/splay-tree.h"
80#include "MagickCore/statistic.h"
81#include "MagickCore/string_.h"
82#include "MagickCore/string-private.h"
83#include "MagickCore/thread-private.h"
84#include "MagickCore/token.h"
85#include "MagickCore/utility.h"
86#include "MagickCore/utility-private.h"
87
88/*
89 Other global definitions used by module.
90*/
91#define Minimize(assign,value) assign=MagickMin(assign,value)
92#define Maximize(assign,value) assign=MagickMax(assign,value)
93
94/* Integer Factorial Function - for a Binomial kernel */
95#if 1
96static inline size_t fact(size_t n)
97{
98 size_t f,l;
99 for(f=1, l=2; l <= n; f=f*l, l++);
100 return(f);
101}
102#elif 1 /* glibc floating point alternatives */
103#define fact(n) ((size_t)tgamma((double)n+1))
104#else
105#define fact(n) ((size_t)lgamma((double)n+1))
106#endif
107
108
109/* Currently these are only internal to this module */
110static void
111 CalcKernelMetaData(KernelInfo *),
112 ExpandMirrorKernelInfo(KernelInfo *),
113 ExpandRotateKernelInfo(KernelInfo *, const double),
114 RotateKernelInfo(KernelInfo *, double);
115
116
117/* Quick function to find last kernel in a kernel list */
118static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
119{
120 while (kernel->next != (KernelInfo *) NULL)
121 kernel=kernel->next;
122 return(kernel);
123}
124
125/*
126%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127% %
128% %
129% %
130% A c q u i r e K e r n e l I n f o %
131% %
132% %
133% %
134%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
135%
136% AcquireKernelInfo() takes the given string (generally supplied by the
137% user) and converts it into a Morphology/Convolution Kernel. This allows
138% users to specify a kernel from a number of pre-defined kernels, or to fully
139% specify their own kernel for a specific Convolution or Morphology
140% Operation.
141%
142% The kernel so generated can be any rectangular array of floating point
143% values (doubles) with the 'control point' or 'pixel being affected'
144% anywhere within that array of values.
145%
146% Previously IM was restricted to a square of odd size using the exact
147% center as origin, this is no longer the case, and any rectangular kernel
148% with any value being declared the origin. This in turn allows the use of
149% highly asymmetrical kernels.
150%
151% The floating point values in the kernel can also include a special value
152% known as 'nan' or 'not a number' to indicate that this value is not part
153% of the kernel array. This allows you to shaped the kernel within its
154% rectangular area. That is 'nan' values provide a 'mask' for the kernel
155% shape. However at least one non-nan value must be provided for correct
156% working of a kernel.
157%
158% The returned kernel should be freed using the DestroyKernelInfo() when you
159% are finished with it. Do not free this memory yourself.
160%
161% Input kernel definition strings can consist of any of three types.
162%
163% "name:args[[@><]"
164% Select from one of the built in kernels, using the name and
165% geometry arguments supplied. See AcquireKernelBuiltIn()
166%
167% "WxH[+X+Y][@><]:num, num, num ..."
168% a kernel of size W by H, with W*H floating point numbers following.
169% the 'center' can be optionally be defined at +X+Y (such that +0+0
170% is top left corner). If not defined the pixel in the center, for
171% odd sizes, or to the immediate top or left of center for even sizes
172% is automatically selected.
173%
174% "num, num, num, num, ..."
175% list of floating point numbers defining an 'old style' odd sized
176% square kernel. At least 9 values should be provided for a 3x3
177% square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
178% Values can be space or comma separated. This is not recommended.
179%
180% You can define a 'list of kernels' which can be used by some morphology
181% operators A list is defined as a semi-colon separated list kernels.
182%
183% " kernel ; kernel ; kernel ; "
184%
185% Any extra ';' characters, at start, end or between kernel definitions are
186% simply ignored.
187%
188% The special flags will expand a single kernel, into a list of rotated
189% kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
190% cyclic rotations, while a '>' will generate a list of 90-degree rotations.
191% The '<' also expands using 90-degree rotates, but giving a 180-degree
192% reflected kernel before the +/- 90-degree rotations, which can be important
193% for Thinning operations.
194%
195% Note that 'name' kernels will start with an alphabetic character while the
196% new kernel specification has a ':' character in its specification string.
197% If neither is the case, it is assumed an old style of a simple list of
198% numbers generating a odd-sized square kernel has been given.
199%
200% The format of the AcquireKernel method is:
201%
202% KernelInfo *AcquireKernelInfo(const char *kernel_string)
203%
204% A description of each parameter follows:
205%
206% o kernel_string: the Morphology/Convolution kernel wanted.
207%
208*/
209
210/* This was separated so that it could be used as a separate
211** array input handling function, such as for -color-matrix
212*/
213static KernelInfo *ParseKernelArray(const char *kernel_string)
214{
216 *kernel;
217
218 char
219 token[MagickPathExtent];
220
221 const char
222 *p,
223 *end;
224
225 ssize_t
226 i;
227
228 double
229 nan = sqrt(-1.0); /* Special Value : Not A Number */
230
231 MagickStatusType
232 flags;
233
235 args;
236
237 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
238 if (kernel == (KernelInfo *) NULL)
239 return(kernel);
240 (void) memset(kernel,0,sizeof(*kernel));
241 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
242 kernel->negative_range = kernel->positive_range = 0.0;
243 kernel->type = UserDefinedKernel;
244 kernel->next = (KernelInfo *) NULL;
245 kernel->signature=MagickCoreSignature;
246 if (kernel_string == (const char *) NULL)
247 return(kernel);
248
249 /* find end of this specific kernel definition string */
250 end = strchr(kernel_string, ';');
251 if ( end == (char *) NULL )
252 end = strchr(kernel_string, '\0');
253
254 /* clear flags - for Expanding kernel lists through rotations */
255 flags = NoValue;
256
257 /* Has a ':' in argument - New user kernel specification
258 FUTURE: this split on ':' could be done by StringToken()
259 */
260 p = strchr(kernel_string, ':');
261 if ( p != (char *) NULL && p < end)
262 {
263 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
264 (void) memcpy(token, kernel_string, (size_t) (p-kernel_string));
265 token[p-kernel_string] = '\0';
266 SetGeometryInfo(&args);
267 flags = ParseGeometry(token, &args);
268
269 /* Size handling and checks of geometry settings */
270 if ( (flags & WidthValue) == 0 ) /* if no width then */
271 args.rho = args.sigma; /* then width = height */
272 if ( args.rho < 1.0 ) /* if width too small */
273 args.rho = 1.0; /* then width = 1 */
274 if ( args.sigma < 1.0 ) /* if height too small */
275 args.sigma = args.rho; /* then height = width */
276 kernel->width = (size_t)args.rho;
277 kernel->height = (size_t)args.sigma;
278
279 /* Offset Handling and Checks */
280 if ( args.xi < 0.0 || args.psi < 0.0 )
281 return(DestroyKernelInfo(kernel));
282 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
283 : (ssize_t) (kernel->width-1)/2;
284 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
285 : (ssize_t) (kernel->height-1)/2;
286 if ( kernel->x >= (ssize_t) kernel->width ||
287 kernel->y >= (ssize_t) kernel->height )
288 return(DestroyKernelInfo(kernel));
289
290 p++; /* advance beyond the ':' */
291 }
292 else
293 { /* ELSE - Old old specification, forming odd-square kernel */
294 /* count up number of values given */
295 p=(const char *) kernel_string;
296 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
297 p++; /* ignore "'" chars for convolve filter usage - Cristy */
298 for (i=0; p < end; i++)
299 {
300 (void) GetNextToken(p,&p,MagickPathExtent,token);
301 if (*token == ',')
302 (void) GetNextToken(p,&p,MagickPathExtent,token);
303 }
304 /* set the size of the kernel - old sized square */
305 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
306 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
307 p=(const char *) kernel_string;
308 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
309 p++; /* ignore "'" chars for convolve filter usage - Cristy */
310 }
311
312 /* Read in the kernel values from rest of input string argument */
313 kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
314 kernel->width,kernel->height*sizeof(*kernel->values)));
315 if (kernel->values == (MagickRealType *) NULL)
316 return(DestroyKernelInfo(kernel));
317 kernel->minimum=MagickMaximumValue;
318 kernel->maximum=(-MagickMaximumValue);
319 kernel->negative_range = kernel->positive_range = 0.0;
320 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
321 {
322 (void) GetNextToken(p,&p,MagickPathExtent,token);
323 if (*token == ',')
324 (void) GetNextToken(p,&p,MagickPathExtent,token);
325 if ( LocaleCompare("nan",token) == 0
326 || LocaleCompare("-",token) == 0 ) {
327 kernel->values[i] = nan; /* this value is not part of neighbourhood */
328 }
329 else {
330 kernel->values[i] = StringToDouble(token,(char **) NULL);
331 ( kernel->values[i] < 0)
332 ? ( kernel->negative_range += kernel->values[i] )
333 : ( kernel->positive_range += kernel->values[i] );
334 Minimize(kernel->minimum, kernel->values[i]);
335 Maximize(kernel->maximum, kernel->values[i]);
336 }
337 }
338
339 /* sanity check -- no more values in kernel definition */
340 (void) GetNextToken(p,&p,MagickPathExtent,token);
341 if ( *token != '\0' && *token != ';' && *token != '\'' )
342 return(DestroyKernelInfo(kernel));
343
344#if 0
345 /* this was the old method of handling a incomplete kernel */
346 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
347 Minimize(kernel->minimum, kernel->values[i]);
348 Maximize(kernel->maximum, kernel->values[i]);
349 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
350 kernel->values[i]=0.0;
351 }
352#else
353 /* Number of values for kernel was not enough - Report Error */
354 if ( i < (ssize_t) (kernel->width*kernel->height) )
355 return(DestroyKernelInfo(kernel));
356#endif
357
358 /* check that we received at least one real (non-nan) value! */
359 if (kernel->minimum == MagickMaximumValue)
360 return(DestroyKernelInfo(kernel));
361
362 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
363 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
364 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
365 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
366 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
367 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
368
369 return(kernel);
370}
371
372static KernelInfo *ParseKernelName(const char *kernel_string,
373 ExceptionInfo *exception)
374{
375 char
376 token[MagickPathExtent] = "";
377
378 const char
379 *p,
380 *end;
381
383 args;
384
386 *kernel;
387
388 MagickStatusType
389 flags;
390
391 ssize_t
392 type;
393
394 /* Parse special 'named' kernel */
395 (void) GetNextToken(kernel_string,&p,MagickPathExtent,token);
396 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
397 if ( type < 0 || type == UserDefinedKernel )
398 return((KernelInfo *) NULL); /* not a valid named kernel */
399
400 while (((isspace((int) ((unsigned char) *p)) != 0) ||
401 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
402 p++;
403
404 end = strchr(p, ';'); /* end of this kernel definition */
405 if ( end == (char *) NULL )
406 end = strchr(p, '\0');
407
408 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
409 (void) memcpy(token, p, (size_t) (end-p));
410 token[end-p] = '\0';
411 SetGeometryInfo(&args);
412 flags = ParseGeometry(token, &args);
413
414#if 0
415 /* For Debugging Geometry Input */
416 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
417 flags, args.rho, args.sigma, args.xi, args.psi );
418#endif
419
420 /* special handling of missing values in input string */
421 switch( type ) {
422 /* Shape Kernel Defaults */
423 case UnityKernel:
424 if ( (flags & WidthValue) == 0 )
425 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
426 break;
427 case SquareKernel:
428 case DiamondKernel:
429 case OctagonKernel:
430 case DiskKernel:
431 case PlusKernel:
432 case CrossKernel:
433 if ( (flags & HeightValue) == 0 )
434 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
435 break;
436 case RingKernel:
437 if ( (flags & XValue) == 0 )
438 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
439 break;
440 case RectangleKernel: /* Rectangle - set size defaults */
441 if ( (flags & WidthValue) == 0 ) /* if no width then */
442 args.rho = args.sigma; /* then width = height */
443 if ( args.rho < 1.0 ) /* if width too small */
444 args.rho = 3; /* then width = 3 */
445 if ( args.sigma < 1.0 ) /* if height too small */
446 args.sigma = args.rho; /* then height = width */
447 if ( (flags & XValue) == 0 ) /* center offset if not defined */
448 args.xi = (double)(((ssize_t)args.rho-1)/2);
449 if ( (flags & YValue) == 0 )
450 args.psi = (double)(((ssize_t)args.sigma-1)/2);
451 break;
452 /* Distance Kernel Defaults */
453 case ChebyshevKernel:
454 case ManhattanKernel:
455 case OctagonalKernel:
456 case EuclideanKernel:
457 if ( (flags & HeightValue) == 0 ) /* no distance scale */
458 args.sigma = 100.0; /* default distance scaling */
459 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
460 args.sigma = (double) QuantumRange/(args.sigma+1); /* maximum pixel distance */
461 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
462 args.sigma *= (double) QuantumRange/100.0; /* percentage of color range */
463 break;
464 default:
465 break;
466 }
467
468 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args, exception);
469 if ( kernel == (KernelInfo *) NULL )
470 return(kernel);
471
472 /* global expand to rotated kernel list - only for single kernels */
473 if ( kernel->next == (KernelInfo *) NULL ) {
474 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
475 ExpandRotateKernelInfo(kernel, 45.0);
476 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
477 ExpandRotateKernelInfo(kernel, 90.0);
478 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
479 ExpandMirrorKernelInfo(kernel);
480 }
481
482 return(kernel);
483}
484
485MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string,
486 ExceptionInfo *exception)
487{
489 *kernel,
490 *new_kernel;
491
492 char
493 *kernel_cache,
494 token[MagickPathExtent];
495
496 const char
497 *p;
498
499 if (kernel_string == (const char *) NULL)
500 return(ParseKernelArray(kernel_string));
501 p=kernel_string;
502 kernel_cache=(char *) NULL;
503 if (*kernel_string == '@')
504 {
505 kernel_cache=FileToString(kernel_string,~0UL,exception);
506 if (kernel_cache == (char *) NULL)
507 return((KernelInfo *) NULL);
508 p=(const char *) kernel_cache;
509 }
510 kernel=NULL;
511 while (GetNextToken(p,(const char **) NULL,MagickPathExtent,token), *token != '\0')
512 {
513 /* ignore extra or multiple ';' kernel separators */
514 if (*token != ';')
515 {
516 /* tokens starting with alpha is a Named kernel */
517 if (isalpha((int) ((unsigned char) *token)) != 0)
518 new_kernel=ParseKernelName(p,exception);
519 else /* otherwise a user defined kernel array */
520 new_kernel=ParseKernelArray(p);
521
522 /* Error handling -- this is not proper error handling! */
523 if (new_kernel == (KernelInfo *) NULL)
524 {
525 if (kernel != (KernelInfo *) NULL)
526 kernel=DestroyKernelInfo(kernel);
527 return((KernelInfo *) NULL);
528 }
529
530 /* initialise or append the kernel list */
531 if (kernel == (KernelInfo *) NULL)
532 kernel=new_kernel;
533 else
534 LastKernelInfo(kernel)->next=new_kernel;
535 }
536
537 /* look for the next kernel in list */
538 p=strchr(p,';');
539 if (p == (char *) NULL)
540 break;
541 p++;
542 }
543 if (kernel_cache != (char *) NULL)
544 kernel_cache=DestroyString(kernel_cache);
545 return(kernel);
546}
547
548/*
549%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
550% %
551% %
552% %
553% A c q u i r e K e r n e l B u i l t I n %
554% %
555% %
556% %
557%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
558%
559% AcquireKernelBuiltIn() returned one of the 'named' built-in types of
560% kernels used for special purposes such as gaussian blurring, skeleton
561% pruning, and edge distance determination.
562%
563% They take a KernelType, and a set of geometry style arguments, which were
564% typically decoded from a user supplied string, or from a more complex
565% Morphology Method that was requested.
566%
567% The format of the AcquireKernelBuiltIn method is:
568%
569% KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
570% const GeometryInfo args)
571%
572% A description of each parameter follows:
573%
574% o type: the pre-defined type of kernel wanted
575%
576% o args: arguments defining or modifying the kernel
577%
578% Convolution Kernels
579%
580% Unity
581% The a No-Op or Scaling single element kernel.
582%
583% Gaussian:{radius},{sigma}
584% Generate a two-dimensional gaussian kernel, as used by -gaussian.
585% The sigma for the curve is required. The resulting kernel is
586% normalized,
587%
588% If 'sigma' is zero, you get a single pixel on a field of zeros.
589%
590% NOTE: that the 'radius' is optional, but if provided can limit (clip)
591% the final size of the resulting kernel to a square 2*radius+1 in size.
592% The radius should be at least 2 times that of the sigma value, or
593% sever clipping and aliasing may result. If not given or set to 0 the
594% radius will be determined so as to produce the best minimal error
595% result, which is usually much larger than is normally needed.
596%
597% LoG:{radius},{sigma}
598% "Laplacian of a Gaussian" or "Mexican Hat" Kernel.
599% The supposed ideal edge detection, zero-summing kernel.
600%
601% An alternative to this kernel is to use a "DoG" with a sigma ratio of
602% approx 1.6 (according to wikipedia).
603%
604% DoG:{radius},{sigma1},{sigma2}
605% "Difference of Gaussians" Kernel.
606% As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
607% from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
608% The result is a zero-summing kernel.
609%
610% Blur:{radius},{sigma}[,{angle}]
611% Generates a 1 dimensional or linear gaussian blur, at the angle given
612% (current restricted to orthogonal angles). If a 'radius' is given the
613% kernel is clipped to a width of 2*radius+1. Kernel can be rotated
614% by a 90 degree angle.
615%
616% If 'sigma' is zero, you get a single pixel on a field of zeros.
617%
618% Note that two convolutions with two "Blur" kernels perpendicular to
619% each other, is equivalent to a far larger "Gaussian" kernel with the
620% same sigma value, However it is much faster to apply. This is how the
621% "-blur" operator actually works.
622%
623% Comet:{width},{sigma},{angle}
624% Blur in one direction only, much like how a bright object leaves
625% a comet like trail. The Kernel is actually half a gaussian curve,
626% Adding two such blurs in opposite directions produces a Blur Kernel.
627% Angle can be rotated in multiples of 90 degrees.
628%
629% Note that the first argument is the width of the kernel and not the
630% radius of the kernel.
631%
632% Binomial:[{radius}]
633% Generate a discrete kernel using a 2 dimensional Pascal's Triangle
634% of values. Used for special forma of image filters.
635%
636% # Still to be implemented...
637% #
638% # Filter2D
639% # Filter1D
640% # Set kernel values using a resize filter, and given scale (sigma)
641% # Cylindrical or Linear. Is this possible with an image?
642% #
643%
644% Named Constant Convolution Kernels
645%
646% All these are unscaled, zero-summing kernels by default. As such for
647% non-HDRI version of ImageMagick some form of normalization, user scaling,
648% and biasing the results is recommended, to prevent the resulting image
649% being 'clipped'.
650%
651% The 3x3 kernels (most of these) can be circularly rotated in multiples of
652% 45 degrees to generate the 8 angled variants of each of the kernels.
653%
654% Laplacian:{type}
655% Discrete Laplacian Kernels, (without normalization)
656% Type 0 : 3x3 with center:8 surrounded by -1 (8 neighbourhood)
657% Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
658% Type 2 : 3x3 with center:4 edge:1 corner:-2
659% Type 3 : 3x3 with center:4 edge:-2 corner:1
660% Type 5 : 5x5 laplacian
661% Type 7 : 7x7 laplacian
662% Type 15 : 5x5 LoG (sigma approx 1.4)
663% Type 19 : 9x9 LoG (sigma approx 1.4)
664%
665% Sobel:{angle}
666% Sobel 'Edge' convolution kernel (3x3)
667% | -1, 0, 1 |
668% | -2, 0,-2 |
669% | -1, 0, 1 |
670%
671% Roberts:{angle}
672% Roberts convolution kernel (3x3)
673% | 0, 0, 0 |
674% | -1, 1, 0 |
675% | 0, 0, 0 |
676%
677% Prewitt:{angle}
678% Prewitt Edge convolution kernel (3x3)
679% | -1, 0, 1 |
680% | -1, 0, 1 |
681% | -1, 0, 1 |
682%
683% Compass:{angle}
684% Prewitt's "Compass" convolution kernel (3x3)
685% | -1, 1, 1 |
686% | -1,-2, 1 |
687% | -1, 1, 1 |
688%
689% Kirsch:{angle}
690% Kirsch's "Compass" convolution kernel (3x3)
691% | -3,-3, 5 |
692% | -3, 0, 5 |
693% | -3,-3, 5 |
694%
695% FreiChen:{angle}
696% Frei-Chen Edge Detector is based on a kernel that is similar to
697% the Sobel Kernel, but is designed to be isotropic. That is it takes
698% into account the distance of the diagonal in the kernel.
699%
700% | 1, 0, -1 |
701% | sqrt(2), 0, -sqrt(2) |
702% | 1, 0, -1 |
703%
704% FreiChen:{type},{angle}
705%
706% Frei-Chen Pre-weighted kernels...
707%
708% Type 0: default un-normalized version shown above.
709%
710% Type 1: Orthogonal Kernel (same as type 11 below)
711% | 1, 0, -1 |
712% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
713% | 1, 0, -1 |
714%
715% Type 2: Diagonal form of Kernel...
716% | 1, sqrt(2), 0 |
717% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
718% | 0, -sqrt(2) -1 |
719%
720% However this kernel is als at the heart of the FreiChen Edge Detection
721% Process which uses a set of 9 specially weighted kernel. These 9
722% kernels not be normalized, but directly applied to the image. The
723% results is then added together, to produce the intensity of an edge in
724% a specific direction. The square root of the pixel value can then be
725% taken as the cosine of the edge, and at least 2 such runs at 90 degrees
726% from each other, both the direction and the strength of the edge can be
727% determined.
728%
729% Type 10: All 9 of the following pre-weighted kernels...
730%
731% Type 11: | 1, 0, -1 |
732% | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
733% | 1, 0, -1 |
734%
735% Type 12: | 1, sqrt(2), 1 |
736% | 0, 0, 0 | / 2*sqrt(2)
737% | 1, sqrt(2), 1 |
738%
739% Type 13: | sqrt(2), -1, 0 |
740% | -1, 0, 1 | / 2*sqrt(2)
741% | 0, 1, -sqrt(2) |
742%
743% Type 14: | 0, 1, -sqrt(2) |
744% | -1, 0, 1 | / 2*sqrt(2)
745% | sqrt(2), -1, 0 |
746%
747% Type 15: | 0, -1, 0 |
748% | 1, 0, 1 | / 2
749% | 0, -1, 0 |
750%
751% Type 16: | 1, 0, -1 |
752% | 0, 0, 0 | / 2
753% | -1, 0, 1 |
754%
755% Type 17: | 1, -2, 1 |
756% | -2, 4, -2 | / 6
757% | -1, -2, 1 |
758%
759% Type 18: | -2, 1, -2 |
760% | 1, 4, 1 | / 6
761% | -2, 1, -2 |
762%
763% Type 19: | 1, 1, 1 |
764% | 1, 1, 1 | / 3
765% | 1, 1, 1 |
766%
767% The first 4 are for edge detection, the next 4 are for line detection
768% and the last is to add a average component to the results.
769%
770% Using a special type of '-1' will return all 9 pre-weighted kernels
771% as a multi-kernel list, so that you can use them directly (without
772% normalization) with the special "-set option:morphology:compose Plus"
773% setting to apply the full FreiChen Edge Detection Technique.
774%
775% If 'type' is large it will be taken to be an actual rotation angle for
776% the default FreiChen (type 0) kernel. As such FreiChen:45 will look
777% like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
778%
779% WARNING: The above was layed out as per
780% http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
781% But rotated 90 degrees so direction is from left rather than the top.
782% I have yet to find any secondary confirmation of the above. The only
783% other source found was actual source code at
784% http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
785% Neither paper defines the kernels in a way that looks logical or
786% correct when taken as a whole.
787%
788% Boolean Kernels
789%
790% Diamond:[{radius}[,{scale}]]
791% Generate a diamond shaped kernel with given radius to the points.
792% Kernel size will again be radius*2+1 square and defaults to radius 1,
793% generating a 3x3 kernel that is slightly larger than a square.
794%
795% Square:[{radius}[,{scale}]]
796% Generate a square shaped kernel of size radius*2+1, and defaulting
797% to a 3x3 (radius 1).
798%
799% Octagon:[{radius}[,{scale}]]
800% Generate octagonal shaped kernel of given radius and constant scale.
801% Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
802% in "Diamond" kernel.
803%
804% Disk:[{radius}[,{scale}]]
805% Generate a binary disk, thresholded at the radius given, the radius
806% may be a float-point value. Final Kernel size is floor(radius)*2+1
807% square. A radius of 5.3 is the default.
808%
809% NOTE: That a low radii Disk kernels produce the same results as
810% many of the previously defined kernels, but differ greatly at larger
811% radii. Here is a table of equivalences...
812% "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
813% "Disk:1.5" => "Square"
814% "Disk:2" => "Diamond:2"
815% "Disk:2.5" => "Octagon"
816% "Disk:2.9" => "Square:2"
817% "Disk:3.5" => "Octagon:3"
818% "Disk:4.5" => "Octagon:4"
819% "Disk:5.4" => "Octagon:5"
820% "Disk:6.4" => "Octagon:6"
821% All other Disk shapes are unique to this kernel, but because a "Disk"
822% is more circular when using a larger radius, using a larger radius is
823% preferred over iterating the morphological operation.
824%
825% Rectangle:{geometry}
826% Simply generate a rectangle of 1's with the size given. You can also
827% specify the location of the 'control point', otherwise the closest
828% pixel to the center of the rectangle is selected.
829%
830% Properly centered and odd sized rectangles work the best.
831%
832% Symbol Dilation Kernels
833%
834% These kernel is not a good general morphological kernel, but is used
835% more for highlighting and marking any single pixels in an image using,
836% a "Dilate" method as appropriate.
837%
838% For the same reasons iterating these kernels does not produce the
839% same result as using a larger radius for the symbol.
840%
841% Plus:[{radius}[,{scale}]]
842% Cross:[{radius}[,{scale}]]
843% Generate a kernel in the shape of a 'plus' or a 'cross' with
844% a each arm the length of the given radius (default 2).
845%
846% NOTE: "plus:1" is equivalent to a "Diamond" kernel.
847%
848% Ring:{radius1},{radius2}[,{scale}]
849% A ring of the values given that falls between the two radii.
850% Defaults to a ring of approximately 3 radius in a 7x7 kernel.
851% This is the 'edge' pixels of the default "Disk" kernel,
852% More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
853%
854% Hit and Miss Kernels
855%
856% Peak:radius1,radius2
857% Find any peak larger than the pixels the fall between the two radii.
858% The default ring of pixels is as per "Ring".
859% Edges
860% Find flat orthogonal edges of a binary shape
861% Corners
862% Find 90 degree corners of a binary shape
863% Diagonals:type
864% A special kernel to thin the 'outside' of diagonals
865% LineEnds:type
866% Find end points of lines (for pruning a skeleton)
867% Two types of lines ends (default to both) can be searched for
868% Type 0: All line ends
869% Type 1: single kernel for 4-connected line ends
870% Type 2: single kernel for simple line ends
871% LineJunctions
872% Find three line junctions (within a skeleton)
873% Type 0: all line junctions
874% Type 1: Y Junction kernel
875% Type 2: Diagonal T Junction kernel
876% Type 3: Orthogonal T Junction kernel
877% Type 4: Diagonal X Junction kernel
878% Type 5: Orthogonal + Junction kernel
879% Ridges:type
880% Find single pixel ridges or thin lines
881% Type 1: Fine single pixel thick lines and ridges
882% Type 2: Find two pixel thick lines and ridges
883% ConvexHull
884% Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
885% Skeleton:type
886% Traditional skeleton generating kernels.
887% Type 1: Traditional Skeleton kernel (4 connected skeleton)
888% Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
889% Type 3: Thinning skeleton based on a research paper by
890% Dan S. Bloomberg (Default Type)
891% ThinSE:type
892% A huge variety of Thinning Kernels designed to preserve connectivity.
893% many other kernel sets use these kernels as source definitions.
894% Type numbers are 41-49, 81-89, 481, and 482 which are based on
895% the super and sub notations used in the source research paper.
896%
897% Distance Measuring Kernels
898%
899% Different types of distance measuring methods, which are used with the
900% a 'Distance' morphology method for generating a gradient based on
901% distance from an edge of a binary shape, though there is a technique
902% for handling a anti-aliased shape.
903%
904% See the 'Distance' Morphological Method, for information of how it is
905% applied.
906%
907% Chebyshev:[{radius}][x{scale}[%!]]
908% Chebyshev Distance (also known as Tchebychev or Chessboard distance)
909% is a value of one to any neighbour, orthogonal or diagonal. One why
910% of thinking of it is the number of squares a 'King' or 'Queen' in
911% chess needs to traverse reach any other position on a chess board.
912% It results in a 'square' like distance function, but one where
913% diagonals are given a value that is closer than expected.
914%
915% Manhattan:[{radius}][x{scale}[%!]]
916% Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
917% Cab distance metric), it is the distance needed when you can only
918% travel in horizontal or vertical directions only. It is the
919% distance a 'Rook' in chess would have to travel, and results in a
920% diamond like distances, where diagonals are further than expected.
921%
922% Octagonal:[{radius}][x{scale}[%!]]
923% An interleaving of Manhattan and Chebyshev metrics producing an
924% increasing octagonally shaped distance. Distances matches those of
925% the "Octagon" shaped kernel of the same radius. The minimum radius
926% and default is 2, producing a 5x5 kernel.
927%
928% Euclidean:[{radius}][x{scale}[%!]]
929% Euclidean distance is the 'direct' or 'as the crow flys' distance.
930% However by default the kernel size only has a radius of 1, which
931% limits the distance to 'Knight' like moves, with only orthogonal and
932% diagonal measurements being correct. As such for the default kernel
933% you will get octagonal like distance function.
934%
935% However using a larger radius such as "Euclidean:4" you will get a
936% much smoother distance gradient from the edge of the shape. Especially
937% if the image is pre-processed to include any anti-aliasing pixels.
938% Of course a larger kernel is slower to use, and not always needed.
939%
940% The first three Distance Measuring Kernels will only generate distances
941% of exact multiples of {scale} in binary images. As such you can use a
942% scale of 1 without loosing any information. However you also need some
943% scaling when handling non-binary anti-aliased shapes.
944%
945% The "Euclidean" Distance Kernel however does generate a non-integer
946% fractional results, and as such scaling is vital even for binary shapes.
947%
948*/
949
950MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
951 const GeometryInfo *args,ExceptionInfo *exception)
952{
954 *kernel;
955
956 ssize_t
957 i;
958
959 ssize_t
960 u,
961 v;
962
963 double
964 nan = sqrt(-1.0); /* Special Value : Not A Number */
965
966 /* Generate a new empty kernel if needed */
967 kernel=(KernelInfo *) NULL;
968 switch(type) {
969 case UndefinedKernel: /* These should not call this function */
970 case UserDefinedKernel:
971 ThrowMagickException(exception,GetMagickModule(),OptionWarning,
972 "InvalidOption","`%s'","Should not call this function");
973 return((KernelInfo *) NULL);
974 case LaplacianKernel: /* Named Discrete Convolution Kernels */
975 case SobelKernel: /* these are defined using other kernels */
976 case RobertsKernel:
977 case PrewittKernel:
978 case CompassKernel:
979 case KirschKernel:
980 case FreiChenKernel:
981 case EdgesKernel: /* Hit and Miss kernels */
982 case CornersKernel:
983 case DiagonalsKernel:
984 case LineEndsKernel:
985 case LineJunctionsKernel:
986 case RidgesKernel:
987 case ConvexHullKernel:
988 case SkeletonKernel:
989 case ThinSEKernel:
990 break; /* A pre-generated kernel is not needed */
991#if 0
992 /* set to 1 to do a compile-time check that we haven't missed anything */
993 case UnityKernel:
994 case GaussianKernel:
995 case DoGKernel:
996 case LoGKernel:
997 case BlurKernel:
998 case CometKernel:
999 case BinomialKernel:
1000 case DiamondKernel:
1001 case SquareKernel:
1002 case RectangleKernel:
1003 case OctagonKernel:
1004 case DiskKernel:
1005 case PlusKernel:
1006 case CrossKernel:
1007 case RingKernel:
1008 case PeaksKernel:
1009 case ChebyshevKernel:
1010 case ManhattanKernel:
1011 case OctagonalKernel:
1012 case EuclideanKernel:
1013#else
1014 default:
1015#endif
1016 /* Generate the base Kernel Structure */
1017 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1018 if (kernel == (KernelInfo *) NULL)
1019 return(kernel);
1020 (void) memset(kernel,0,sizeof(*kernel));
1021 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1022 kernel->negative_range = kernel->positive_range = 0.0;
1023 kernel->type = type;
1024 kernel->next = (KernelInfo *) NULL;
1025 kernel->signature=MagickCoreSignature;
1026 break;
1027 }
1028
1029 switch(type) {
1030 /*
1031 Convolution Kernels
1032 */
1033 case UnityKernel:
1034 {
1035 kernel->height = kernel->width = (size_t) 1;
1036 kernel->x = kernel->y = (ssize_t) 0;
1037 kernel->values=(MagickRealType *) MagickAssumeAligned(
1038 AcquireAlignedMemory(1,sizeof(*kernel->values)));
1039 if (kernel->values == (MagickRealType *) NULL)
1040 return(DestroyKernelInfo(kernel));
1041 kernel->maximum = kernel->values[0] = args->rho;
1042 break;
1043 }
1044 break;
1045 case GaussianKernel:
1046 case DoGKernel:
1047 case LoGKernel:
1048 { double
1049 sigma = fabs(args->sigma),
1050 sigma2 = fabs(args->xi),
1051 A, B, R;
1052
1053 if ( args->rho >= 1.0 )
1054 kernel->width = (size_t)args->rho*2+1;
1055 else if ( (type != DoGKernel) || (sigma >= sigma2) )
1056 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1057 else
1058 kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1059 kernel->height = kernel->width;
1060 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1061 kernel->values=(MagickRealType *) MagickAssumeAligned(
1062 AcquireAlignedMemory(kernel->width,kernel->height*
1063 sizeof(*kernel->values)));
1064 if (kernel->values == (MagickRealType *) NULL)
1065 return(DestroyKernelInfo(kernel));
1066
1067 /* WARNING: The following generates a 'sampled gaussian' kernel.
1068 * What we really want is a 'discrete gaussian' kernel.
1069 *
1070 * How to do this is I don't know, but appears to be basied on the
1071 * Error Function 'erf()' (integral of a gaussian)
1072 */
1073
1074 if ( type == GaussianKernel || type == DoGKernel )
1075 { /* Calculate a Gaussian, OR positive half of a DoG */
1076 if ( sigma > MagickEpsilon )
1077 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1078 B = (double) (1.0/(Magick2PI*sigma*sigma));
1079 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1080 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1081 kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1082 }
1083 else /* limiting case - a unity (normalized Dirac) kernel */
1084 { (void) memset(kernel->values,0, (size_t)
1085 kernel->width*kernel->height*sizeof(*kernel->values));
1086 kernel->values[kernel->x+kernel->y*(ssize_t) kernel->width] = 1.0;
1087 }
1088 }
1089
1090 if ( type == DoGKernel )
1091 { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1092 if ( sigma2 > MagickEpsilon )
1093 { sigma = sigma2; /* simplify loop expressions */
1094 A = 1.0/(2.0*sigma*sigma);
1095 B = (double) (1.0/(Magick2PI*sigma*sigma));
1096 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1097 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1098 kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1099 }
1100 else /* limiting case - a unity (normalized Dirac) kernel */
1101 kernel->values[kernel->x+kernel->y*(ssize_t) kernel->width] -= 1.0;
1102 }
1103
1104 if ( type == LoGKernel )
1105 { /* Calculate a Laplacian of a Gaussian - Or Mexican Hat */
1106 if ( sigma > MagickEpsilon )
1107 { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1108 B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1109 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1110 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1111 { R = ((double)(u*u+v*v))*A;
1112 kernel->values[i] = (1-R)*exp(-R)*B;
1113 }
1114 }
1115 else /* special case - generate a unity kernel */
1116 { (void) memset(kernel->values,0, (size_t)
1117 kernel->width*kernel->height*sizeof(*kernel->values));
1118 kernel->values[kernel->x+kernel->y*(ssize_t) kernel->width] = 1.0;
1119 }
1120 }
1121
1122 /* Note the above kernels may have been 'clipped' by a user defined
1123 ** radius, producing a smaller (darker) kernel. Also for very small
1124 ** sigma's (> 0.1) the central value becomes larger than one, and thus
1125 ** producing a very bright kernel.
1126 **
1127 ** Normalization will still be needed.
1128 */
1129
1130 /* Normalize the 2D Gaussian Kernel
1131 **
1132 ** NB: a CorrelateNormalize performs a normal Normalize if
1133 ** there are no negative values.
1134 */
1135 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1136 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1137
1138 break;
1139 }
1140 case BlurKernel:
1141 { double
1142 sigma = fabs(args->sigma),
1143 alpha, beta;
1144
1145 if ( args->rho >= 1.0 )
1146 kernel->width = (size_t)args->rho*2+1;
1147 else
1148 kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1149 kernel->height = 1;
1150 kernel->x = (ssize_t) (kernel->width-1)/2;
1151 kernel->y = 0;
1152 kernel->negative_range = kernel->positive_range = 0.0;
1153 kernel->values=(MagickRealType *) MagickAssumeAligned(
1154 AcquireAlignedMemory(kernel->width,kernel->height*
1155 sizeof(*kernel->values)));
1156 if (kernel->values == (MagickRealType *) NULL)
1157 return(DestroyKernelInfo(kernel));
1158
1159#if 1
1160#define KernelRank 3
1161 /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1162 ** It generates a gaussian 3 times the width, and compresses it into
1163 ** the expected range. This produces a closer normalization of the
1164 ** resulting kernel, especially for very low sigma values.
1165 ** As such while wierd it is prefered.
1166 **
1167 ** I am told this method originally came from Photoshop.
1168 **
1169 ** A properly normalized curve is generated (apart from edge clipping)
1170 ** even though we later normalize the result (for edge clipping)
1171 ** to allow the correct generation of a "Difference of Blurs".
1172 */
1173
1174 /* initialize */
1175 v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1176 (void) memset(kernel->values,0, (size_t)
1177 kernel->width*kernel->height*sizeof(*kernel->values));
1178 /* Calculate a Positive 1D Gaussian */
1179 if ( sigma > MagickEpsilon )
1180 { sigma *= KernelRank; /* simplify loop expressions */
1181 alpha = 1.0/(2.0*sigma*sigma);
1182 beta= (double) (1.0/(MagickSQ2PI*sigma ));
1183 for ( u=-v; u <= v; u++) {
1184 kernel->values[(u+v)/KernelRank] +=
1185 exp(-((double)(u*u))*alpha)*beta;
1186 }
1187 }
1188 else /* special case - generate a unity kernel */
1189 kernel->values[kernel->x+kernel->y*(ssize_t) kernel->width] = 1.0;
1190#else
1191 /* Direct calculation without curve averaging
1192 This is equivalent to a KernelRank of 1 */
1193
1194 /* Calculate a Positive Gaussian */
1195 if ( sigma > MagickEpsilon )
1196 { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1197 beta = 1.0/(MagickSQ2PI*sigma);
1198 for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1199 kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1200 }
1201 else /* special case - generate a unity kernel */
1202 { (void) memset(kernel->values,0, (size_t)
1203 kernel->width*kernel->height*sizeof(*kernel->values));
1204 kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1205 }
1206#endif
1207 /* Note the above kernel may have been 'clipped' by a user defined
1208 ** radius, producing a smaller (darker) kernel. Also for very small
1209 ** sigma's (> 0.1) the central value becomes larger than one, as a
1210 ** result of not generating a actual 'discrete' kernel, and thus
1211 ** producing a very bright 'impulse'.
1212 **
1213 ** Because of these two factors Normalization is required!
1214 */
1215
1216 /* Normalize the 1D Gaussian Kernel
1217 **
1218 ** NB: a CorrelateNormalize performs a normal Normalize if
1219 ** there are no negative values.
1220 */
1221 CalcKernelMetaData(kernel); /* the other kernel meta-data */
1222 ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
1223
1224 /* rotate the 1D kernel by given angle */
1225 RotateKernelInfo(kernel, args->xi );
1226 break;
1227 }
1228 case CometKernel:
1229 { double
1230 sigma = fabs(args->sigma),
1231 A;
1232
1233 if ( args->rho < 1.0 )
1234 kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1235 else
1236 kernel->width = (size_t)args->rho;
1237 kernel->x = kernel->y = 0;
1238 kernel->height = 1;
1239 kernel->negative_range = kernel->positive_range = 0.0;
1240 kernel->values=(MagickRealType *) MagickAssumeAligned(
1241 AcquireAlignedMemory(kernel->width,kernel->height*
1242 sizeof(*kernel->values)));
1243 if (kernel->values == (MagickRealType *) NULL)
1244 return(DestroyKernelInfo(kernel));
1245
1246 /* A comet blur is half a 1D gaussian curve, so that the object is
1247 ** blurred in one direction only. This may not be quite the right
1248 ** curve to use so may change in the future. The function must be
1249 ** normalised after generation, which also resolves any clipping.
1250 **
1251 ** As we are normalizing and not subtracting gaussians,
1252 ** there is no need for a divisor in the gaussian formula
1253 **
1254 ** It is less complex
1255 */
1256 if ( sigma > MagickEpsilon )
1257 {
1258#if 1
1259#define KernelRank 3
1260 v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1261 (void) memset(kernel->values,0, (size_t)
1262 kernel->width*sizeof(*kernel->values));
1263 sigma *= KernelRank; /* simplify the loop expression */
1264 A = 1.0/(2.0*sigma*sigma);
1265 /* B = 1.0/(MagickSQ2PI*sigma); */
1266 for ( u=0; u < v; u++) {
1267 kernel->values[u/KernelRank] +=
1268 exp(-((double)(u*u))*A);
1269 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1270 }
1271 for (i=0; i < (ssize_t) kernel->width; i++)
1272 kernel->positive_range += kernel->values[i];
1273#else
1274 A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1275 /* B = 1.0/(MagickSQ2PI*sigma); */
1276 for ( i=0; i < (ssize_t) kernel->width; i++)
1277 kernel->positive_range +=
1278 kernel->values[i] = exp(-((double)(i*i))*A);
1279 /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1280#endif
1281 }
1282 else /* special case - generate a unity kernel */
1283 { (void) memset(kernel->values,0, (size_t)
1284 kernel->width*kernel->height*sizeof(*kernel->values));
1285 kernel->values[kernel->x+kernel->y*(ssize_t) kernel->width] = 1.0;
1286 kernel->positive_range = 1.0;
1287 }
1288
1289 kernel->minimum = 0.0;
1290 kernel->maximum = kernel->values[0];
1291 kernel->negative_range = 0.0;
1292
1293 ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1294 RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1295 break;
1296 }
1297 case BinomialKernel:
1298 {
1299 size_t
1300 order_f;
1301
1302 if (args->rho < 1.0)
1303 kernel->width = kernel->height = 3; /* default radius = 1 */
1304 else
1305 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1306 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1307
1308 order_f = fact(kernel->width-1);
1309
1310 kernel->values=(MagickRealType *) MagickAssumeAligned(
1311 AcquireAlignedMemory(kernel->width,kernel->height*
1312 sizeof(*kernel->values)));
1313 if (kernel->values == (MagickRealType *) NULL)
1314 return(DestroyKernelInfo(kernel));
1315
1316 /* set all kernel values within diamond area to scale given */
1317 for ( i=0, v=0; v < (ssize_t)kernel->height; v++)
1318 { size_t
1319 alpha = order_f / ( fact((size_t) v) * fact(kernel->height-(size_t) v-1) );
1320 for ( u=0; u < (ssize_t)kernel->width; u++, i++)
1321 kernel->positive_range += kernel->values[i] = (double)
1322 (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-(size_t) u-1) ));
1323 }
1324 kernel->minimum = 1.0;
1325 kernel->maximum = kernel->values[kernel->x+kernel->y*(ssize_t) kernel->width];
1326 kernel->negative_range = 0.0;
1327 break;
1328 }
1329
1330 /*
1331 Convolution Kernels - Well Known Named Constant Kernels
1332 */
1333 case LaplacianKernel:
1334 { switch ( (int) args->rho ) {
1335 case 0:
1336 default: /* laplacian square filter -- default */
1337 kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1338 break;
1339 case 1: /* laplacian diamond filter */
1340 kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1341 break;
1342 case 2:
1343 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1344 break;
1345 case 3:
1346 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1347 break;
1348 case 5: /* a 5x5 laplacian */
1349 kernel=ParseKernelArray(
1350 "5: -4,-1,0,-1,-4 -1,2,3,2,-1 0,3,4,3,0 -1,2,3,2,-1 -4,-1,0,-1,-4");
1351 break;
1352 case 7: /* a 7x7 laplacian */
1353 kernel=ParseKernelArray(
1354 "7:-10,-5,-2,-1,-2,-5,-10 -5,0,3,4,3,0,-5 -2,3,6,7,6,3,-2 -1,4,7,8,7,4,-1 -2,3,6,7,6,3,-2 -5,0,3,4,3,0,-5 -10,-5,-2,-1,-2,-5,-10" );
1355 break;
1356 case 15: /* a 5x5 LoG (sigma approx 1.4) */
1357 kernel=ParseKernelArray(
1358 "5: 0,0,-1,0,0 0,-1,-2,-1,0 -1,-2,16,-2,-1 0,-1,-2,-1,0 0,0,-1,0,0");
1359 break;
1360 case 19: /* a 9x9 LoG (sigma approx 1.4) */
1361 /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1362 kernel=ParseKernelArray(
1363 "9: 0,-1,-1,-2,-2,-2,-1,-1,0 -1,-2,-4,-5,-5,-5,-4,-2,-1 -1,-4,-5,-3,-0,-3,-5,-4,-1 -2,-5,-3,12,24,12,-3,-5,-2 -2,-5,-0,24,40,24,-0,-5,-2 -2,-5,-3,12,24,12,-3,-5,-2 -1,-4,-5,-3,-0,-3,-5,-4,-1 -1,-2,-4,-5,-5,-5,-4,-2,-1 0,-1,-1,-2,-2,-2,-1,-1,0");
1364 break;
1365 }
1366 if (kernel == (KernelInfo *) NULL)
1367 return(kernel);
1368 kernel->type = type;
1369 break;
1370 }
1371 case SobelKernel:
1372 { /* Simple Sobel Kernel */
1373 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1374 if (kernel == (KernelInfo *) NULL)
1375 return(kernel);
1376 kernel->type = type;
1377 RotateKernelInfo(kernel, args->rho);
1378 break;
1379 }
1380 case RobertsKernel:
1381 {
1382 kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1383 if (kernel == (KernelInfo *) NULL)
1384 return(kernel);
1385 kernel->type = type;
1386 RotateKernelInfo(kernel, args->rho);
1387 break;
1388 }
1389 case PrewittKernel:
1390 {
1391 kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1392 if (kernel == (KernelInfo *) NULL)
1393 return(kernel);
1394 kernel->type = type;
1395 RotateKernelInfo(kernel, args->rho);
1396 break;
1397 }
1398 case CompassKernel:
1399 {
1400 kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1401 if (kernel == (KernelInfo *) NULL)
1402 return(kernel);
1403 kernel->type = type;
1404 RotateKernelInfo(kernel, args->rho);
1405 break;
1406 }
1407 case KirschKernel:
1408 {
1409 kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1410 if (kernel == (KernelInfo *) NULL)
1411 return(kernel);
1412 kernel->type = type;
1413 RotateKernelInfo(kernel, args->rho);
1414 break;
1415 }
1416 case FreiChenKernel:
1417 /* Direction is set to be left to right positive */
1418 /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1419 /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1420 { switch ( (int) args->rho ) {
1421 default:
1422 case 0:
1423 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1424 if (kernel == (KernelInfo *) NULL)
1425 return(kernel);
1426 kernel->type = type;
1427 kernel->values[3] = +(MagickRealType) MagickSQ2;
1428 kernel->values[5] = -(MagickRealType) MagickSQ2;
1429 CalcKernelMetaData(kernel); /* recalculate meta-data */
1430 break;
1431 case 2:
1432 kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1433 if (kernel == (KernelInfo *) NULL)
1434 return(kernel);
1435 kernel->type = type;
1436 kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2;
1437 kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2;
1438 CalcKernelMetaData(kernel); /* recalculate meta-data */
1439 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1440 break;
1441 case 10:
1442 {
1443 kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19",exception);
1444 if (kernel == (KernelInfo *) NULL)
1445 return(kernel);
1446 break;
1447 }
1448 case 1:
1449 case 11:
1450 kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1451 if (kernel == (KernelInfo *) NULL)
1452 return(kernel);
1453 kernel->type = type;
1454 kernel->values[3] = +(MagickRealType) MagickSQ2;
1455 kernel->values[5] = -(MagickRealType) MagickSQ2;
1456 CalcKernelMetaData(kernel); /* recalculate meta-data */
1457 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1458 break;
1459 case 12:
1460 kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1461 if (kernel == (KernelInfo *) NULL)
1462 return(kernel);
1463 kernel->type = type;
1464 kernel->values[1] = +(MagickRealType) MagickSQ2;
1465 kernel->values[7] = +(MagickRealType) MagickSQ2;
1466 CalcKernelMetaData(kernel);
1467 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1468 break;
1469 case 13:
1470 kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1471 if (kernel == (KernelInfo *) NULL)
1472 return(kernel);
1473 kernel->type = type;
1474 kernel->values[0] = +(MagickRealType) MagickSQ2;
1475 kernel->values[8] = -(MagickRealType) MagickSQ2;
1476 CalcKernelMetaData(kernel);
1477 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1478 break;
1479 case 14:
1480 kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1481 if (kernel == (KernelInfo *) NULL)
1482 return(kernel);
1483 kernel->type = type;
1484 kernel->values[2] = -(MagickRealType) MagickSQ2;
1485 kernel->values[6] = +(MagickRealType) MagickSQ2;
1486 CalcKernelMetaData(kernel);
1487 ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1488 break;
1489 case 15:
1490 kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1491 if (kernel == (KernelInfo *) NULL)
1492 return(kernel);
1493 kernel->type = type;
1494 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1495 break;
1496 case 16:
1497 kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1498 if (kernel == (KernelInfo *) NULL)
1499 return(kernel);
1500 kernel->type = type;
1501 ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1502 break;
1503 case 17:
1504 kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1505 if (kernel == (KernelInfo *) NULL)
1506 return(kernel);
1507 kernel->type = type;
1508 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1509 break;
1510 case 18:
1511 kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1512 if (kernel == (KernelInfo *) NULL)
1513 return(kernel);
1514 kernel->type = type;
1515 ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1516 break;
1517 case 19:
1518 kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1519 if (kernel == (KernelInfo *) NULL)
1520 return(kernel);
1521 kernel->type = type;
1522 ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1523 break;
1524 }
1525 if ( fabs(args->sigma) >= MagickEpsilon )
1526 /* Rotate by correctly supplied 'angle' */
1527 RotateKernelInfo(kernel, args->sigma);
1528 else if ( args->rho > 30.0 || args->rho < -30.0 )
1529 /* Rotate by out of bounds 'type' */
1530 RotateKernelInfo(kernel, args->rho);
1531 break;
1532 }
1533
1534 /*
1535 Boolean or Shaped Kernels
1536 */
1537 case DiamondKernel:
1538 {
1539 if (args->rho < 1.0)
1540 kernel->width = kernel->height = 3; /* default radius = 1 */
1541 else
1542 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1543 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1544
1545 kernel->values=(MagickRealType *) MagickAssumeAligned(
1546 AcquireAlignedMemory(kernel->width,kernel->height*
1547 sizeof(*kernel->values)));
1548 if (kernel->values == (MagickRealType *) NULL)
1549 return(DestroyKernelInfo(kernel));
1550
1551 /* set all kernel values within diamond area to scale given */
1552 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1553 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1554 if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1555 kernel->positive_range += kernel->values[i] = args->sigma;
1556 else
1557 kernel->values[i] = nan;
1558 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1559 break;
1560 }
1561 case SquareKernel:
1562 case RectangleKernel:
1563 { double
1564 scale;
1565 if ( type == SquareKernel )
1566 {
1567 if (args->rho < 1.0)
1568 kernel->width = kernel->height = 3; /* default radius = 1 */
1569 else
1570 kernel->width = kernel->height = (size_t) (2*args->rho+1);
1571 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1572 scale = args->sigma;
1573 }
1574 else {
1575 /* NOTE: user defaults set in "AcquireKernelInfo()" */
1576 if ( args->rho < 1.0 || args->sigma < 1.0 )
1577 return(DestroyKernelInfo(kernel)); /* invalid args given */
1578 kernel->width = (size_t)args->rho;
1579 kernel->height = (size_t)args->sigma;
1580 if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1581 args->psi < 0.0 || args->psi > (double)kernel->height )
1582 return(DestroyKernelInfo(kernel)); /* invalid args given */
1583 kernel->x = (ssize_t) args->xi;
1584 kernel->y = (ssize_t) args->psi;
1585 scale = 1.0;
1586 }
1587 kernel->values=(MagickRealType *) MagickAssumeAligned(
1588 AcquireAlignedMemory(kernel->width,kernel->height*
1589 sizeof(*kernel->values)));
1590 if (kernel->values == (MagickRealType *) NULL)
1591 return(DestroyKernelInfo(kernel));
1592
1593 /* set all kernel values to scale given */
1594 u=(ssize_t) (kernel->width*kernel->height);
1595 for ( i=0; i < u; i++)
1596 kernel->values[i] = scale;
1597 kernel->minimum = kernel->maximum = scale; /* a flat shape */
1598 kernel->positive_range = scale*u;
1599 break;
1600 }
1601 case OctagonKernel:
1602 {
1603 if (args->rho < 1.0)
1604 kernel->width = kernel->height = 5; /* default radius = 2 */
1605 else
1606 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1607 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1608
1609 kernel->values=(MagickRealType *) MagickAssumeAligned(
1610 AcquireAlignedMemory(kernel->width,kernel->height*
1611 sizeof(*kernel->values)));
1612 if (kernel->values == (MagickRealType *) NULL)
1613 return(DestroyKernelInfo(kernel));
1614
1615 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1616 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1617 if ( (labs((long) u)+labs((long) v)) <=
1618 ((long)kernel->x + (long)(kernel->x/2)) )
1619 kernel->positive_range += kernel->values[i] = args->sigma;
1620 else
1621 kernel->values[i] = nan;
1622 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1623 break;
1624 }
1625 case DiskKernel:
1626 {
1627 ssize_t
1628 limit = (ssize_t)(args->rho*args->rho);
1629
1630 if (args->rho < 0.4) /* default radius approx 4.3 */
1631 kernel->width = kernel->height = 9L, limit = 18L;
1632 else
1633 kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1634 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1635
1636 kernel->values=(MagickRealType *) MagickAssumeAligned(
1637 AcquireAlignedMemory(kernel->width,kernel->height*
1638 sizeof(*kernel->values)));
1639 if (kernel->values == (MagickRealType *) NULL)
1640 return(DestroyKernelInfo(kernel));
1641
1642 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1643 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1644 if ((u*u+v*v) <= limit)
1645 kernel->positive_range += kernel->values[i] = args->sigma;
1646 else
1647 kernel->values[i] = nan;
1648 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1649 break;
1650 }
1651 case PlusKernel:
1652 {
1653 if (args->rho < 1.0)
1654 kernel->width = kernel->height = 5; /* default radius 2 */
1655 else
1656 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1657 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1658
1659 kernel->values=(MagickRealType *) MagickAssumeAligned(
1660 AcquireAlignedMemory(kernel->width,kernel->height*
1661 sizeof(*kernel->values)));
1662 if (kernel->values == (MagickRealType *) NULL)
1663 return(DestroyKernelInfo(kernel));
1664
1665 /* set all kernel values along axises to given scale */
1666 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1667 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1668 kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1669 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1670 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1671 break;
1672 }
1673 case CrossKernel:
1674 {
1675 if (args->rho < 1.0)
1676 kernel->width = kernel->height = 5; /* default radius 2 */
1677 else
1678 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1679 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1680
1681 kernel->values=(MagickRealType *) MagickAssumeAligned(
1682 AcquireAlignedMemory(kernel->width,kernel->height*
1683 sizeof(*kernel->values)));
1684 if (kernel->values == (MagickRealType *) NULL)
1685 return(DestroyKernelInfo(kernel));
1686
1687 /* set all kernel values along axises to given scale */
1688 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1689 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1690 kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1691 kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1692 kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1693 break;
1694 }
1695 /*
1696 HitAndMiss Kernels
1697 */
1698 case RingKernel:
1699 case PeaksKernel:
1700 {
1701 ssize_t
1702 limit1,
1703 limit2,
1704 scale;
1705
1706 if (args->rho < args->sigma)
1707 {
1708 kernel->width = ((size_t)args->sigma)*2+1;
1709 limit1 = (ssize_t)(args->rho*args->rho);
1710 limit2 = (ssize_t)(args->sigma*args->sigma);
1711 }
1712 else
1713 {
1714 kernel->width = ((size_t)args->rho)*2+1;
1715 limit1 = (ssize_t)(args->sigma*args->sigma);
1716 limit2 = (ssize_t)(args->rho*args->rho);
1717 }
1718 if ( limit2 <= 0 )
1719 kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1720
1721 kernel->height = kernel->width;
1722 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1723 kernel->values=(MagickRealType *) MagickAssumeAligned(
1724 AcquireAlignedMemory(kernel->width,kernel->height*
1725 sizeof(*kernel->values)));
1726 if (kernel->values == (MagickRealType *) NULL)
1727 return(DestroyKernelInfo(kernel));
1728
1729 /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1730 scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1731 for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1732 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1733 { ssize_t radius=u*u+v*v;
1734 if (limit1 < radius && radius <= limit2)
1735 kernel->positive_range += kernel->values[i] = (double) scale;
1736 else
1737 kernel->values[i] = nan;
1738 }
1739 kernel->minimum = kernel->maximum = (double) scale;
1740 if ( type == PeaksKernel ) {
1741 /* set the central point in the middle */
1742 kernel->values[kernel->x+kernel->y*(ssize_t) kernel->width] = 1.0;
1743 kernel->positive_range = 1.0;
1744 kernel->maximum = 1.0;
1745 }
1746 break;
1747 }
1748 case EdgesKernel:
1749 {
1750 kernel=AcquireKernelInfo("ThinSE:482",exception);
1751 if (kernel == (KernelInfo *) NULL)
1752 return(kernel);
1753 kernel->type = type;
1754 ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1755 break;
1756 }
1757 case CornersKernel:
1758 {
1759 kernel=AcquireKernelInfo("ThinSE:87",exception);
1760 if (kernel == (KernelInfo *) NULL)
1761 return(kernel);
1762 kernel->type = type;
1763 ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1764 break;
1765 }
1766 case DiagonalsKernel:
1767 {
1768 switch ( (int) args->rho ) {
1769 case 0:
1770 default:
1771 { KernelInfo
1772 *new_kernel;
1773 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1774 if (kernel == (KernelInfo *) NULL)
1775 return(kernel);
1776 kernel->type = type;
1777 new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1778 if (new_kernel == (KernelInfo *) NULL)
1779 return(DestroyKernelInfo(kernel));
1780 new_kernel->type = type;
1781 LastKernelInfo(kernel)->next = new_kernel;
1782 ExpandMirrorKernelInfo(kernel);
1783 return(kernel);
1784 }
1785 case 1:
1786 kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1787 break;
1788 case 2:
1789 kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1790 break;
1791 }
1792 if (kernel == (KernelInfo *) NULL)
1793 return(kernel);
1794 kernel->type = type;
1795 RotateKernelInfo(kernel, args->sigma);
1796 break;
1797 }
1798 case LineEndsKernel:
1799 { /* Kernels for finding the end of thin lines */
1800 switch ( (int) args->rho ) {
1801 case 0:
1802 default:
1803 /* set of kernels to find all end of lines */
1804 return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>",exception));
1805 case 1:
1806 /* kernel for 4-connected line ends - no rotation */
1807 kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1808 break;
1809 case 2:
1810 /* kernel to add for 8-connected lines - no rotation */
1811 kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1812 break;
1813 case 3:
1814 /* kernel to add for orthogonal line ends - does not find corners */
1815 kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1816 break;
1817 case 4:
1818 /* traditional line end - fails on last T end */
1819 kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1820 break;
1821 }
1822 if (kernel == (KernelInfo *) NULL)
1823 return(kernel);
1824 kernel->type = type;
1825 RotateKernelInfo(kernel, args->sigma);
1826 break;
1827 }
1828 case LineJunctionsKernel:
1829 { /* kernels for finding the junctions of multiple lines */
1830 switch ( (int) args->rho ) {
1831 case 0:
1832 default:
1833 /* set of kernels to find all line junctions */
1834 return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>",exception));
1835 case 1:
1836 /* Y Junction */
1837 kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1838 break;
1839 case 2:
1840 /* Diagonal T Junctions */
1841 kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1842 break;
1843 case 3:
1844 /* Orthogonal T Junctions */
1845 kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1846 break;
1847 case 4:
1848 /* Diagonal X Junctions */
1849 kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1850 break;
1851 case 5:
1852 /* Orthogonal X Junctions - minimal diamond kernel */
1853 kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1854 break;
1855 }
1856 if (kernel == (KernelInfo *) NULL)
1857 return(kernel);
1858 kernel->type = type;
1859 RotateKernelInfo(kernel, args->sigma);
1860 break;
1861 }
1862 case RidgesKernel:
1863 { /* Ridges - Ridge finding kernels */
1865 *new_kernel;
1866 switch ( (int) args->rho ) {
1867 case 1:
1868 default:
1869 kernel=ParseKernelArray("3x1:0,1,0");
1870 if (kernel == (KernelInfo *) NULL)
1871 return(kernel);
1872 kernel->type = type;
1873 ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1874 break;
1875 case 2:
1876 kernel=ParseKernelArray("4x1:0,1,1,0");
1877 if (kernel == (KernelInfo *) NULL)
1878 return(kernel);
1879 kernel->type = type;
1880 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1881
1882 /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1883 /* Unfortunately we can not yet rotate a non-square kernel */
1884 /* But then we can't flip a non-symmetrical kernel either */
1885 new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1886 if (new_kernel == (KernelInfo *) NULL)
1887 return(DestroyKernelInfo(kernel));
1888 new_kernel->type = type;
1889 LastKernelInfo(kernel)->next = new_kernel;
1890 new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1891 if (new_kernel == (KernelInfo *) NULL)
1892 return(DestroyKernelInfo(kernel));
1893 new_kernel->type = type;
1894 LastKernelInfo(kernel)->next = new_kernel;
1895 new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1896 if (new_kernel == (KernelInfo *) NULL)
1897 return(DestroyKernelInfo(kernel));
1898 new_kernel->type = type;
1899 LastKernelInfo(kernel)->next = new_kernel;
1900 new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1901 if (new_kernel == (KernelInfo *) NULL)
1902 return(DestroyKernelInfo(kernel));
1903 new_kernel->type = type;
1904 LastKernelInfo(kernel)->next = new_kernel;
1905 new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1906 if (new_kernel == (KernelInfo *) NULL)
1907 return(DestroyKernelInfo(kernel));
1908 new_kernel->type = type;
1909 LastKernelInfo(kernel)->next = new_kernel;
1910 new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1911 if (new_kernel == (KernelInfo *) NULL)
1912 return(DestroyKernelInfo(kernel));
1913 new_kernel->type = type;
1914 LastKernelInfo(kernel)->next = new_kernel;
1915 new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1916 if (new_kernel == (KernelInfo *) NULL)
1917 return(DestroyKernelInfo(kernel));
1918 new_kernel->type = type;
1919 LastKernelInfo(kernel)->next = new_kernel;
1920 new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1921 if (new_kernel == (KernelInfo *) NULL)
1922 return(DestroyKernelInfo(kernel));
1923 new_kernel->type = type;
1924 LastKernelInfo(kernel)->next = new_kernel;
1925 break;
1926 }
1927 break;
1928 }
1929 case ConvexHullKernel:
1930 {
1932 *new_kernel;
1933 /* first set of 8 kernels */
1934 kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1935 if (kernel == (KernelInfo *) NULL)
1936 return(kernel);
1937 kernel->type = type;
1938 ExpandRotateKernelInfo(kernel, 90.0);
1939 /* append the mirror versions too - no flip function yet */
1940 new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1941 if (new_kernel == (KernelInfo *) NULL)
1942 return(DestroyKernelInfo(kernel));
1943 new_kernel->type = type;
1944 ExpandRotateKernelInfo(new_kernel, 90.0);
1945 LastKernelInfo(kernel)->next = new_kernel;
1946 break;
1947 }
1948 case SkeletonKernel:
1949 {
1950 switch ( (int) args->rho ) {
1951 case 1:
1952 default:
1953 /* Traditional Skeleton...
1954 ** A cyclically rotated single kernel
1955 */
1956 kernel=AcquireKernelInfo("ThinSE:482",exception);
1957 if (kernel == (KernelInfo *) NULL)
1958 return(kernel);
1959 kernel->type = type;
1960 ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1961 break;
1962 case 2:
1963 /* HIPR Variation of the cyclic skeleton
1964 ** Corners of the traditional method made more forgiving,
1965 ** but the retain the same cyclic order.
1966 */
1967 kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;",exception);
1968 if (kernel == (KernelInfo *) NULL)
1969 return(kernel);
1970 if (kernel->next == (KernelInfo *) NULL)
1971 return(DestroyKernelInfo(kernel));
1972 kernel->type = type;
1973 kernel->next->type = type;
1974 ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1975 break;
1976 case 3:
1977 /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1978 ** "Connectivity-Preserving Morphological Image Transformations"
1979 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1980 ** http://www.leptonica.com/papers/conn.pdf
1981 */
1982 kernel=AcquireKernelInfo("ThinSE:41; ThinSE:42; ThinSE:43",
1983 exception);
1984 if (kernel == (KernelInfo *) NULL)
1985 return(kernel);
1986 if (kernel->next == (KernelInfo *) NULL)
1987 return(DestroyKernelInfo(kernel));
1988 if (kernel->next->next == (KernelInfo *) NULL)
1989 return(DestroyKernelInfo(kernel));
1990 kernel->type = type;
1991 kernel->next->type = type;
1992 kernel->next->next->type = type;
1993 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1994 break;
1995 }
1996 break;
1997 }
1998 case ThinSEKernel:
1999 { /* Special kernels for general thinning, while preserving connections
2000 ** "Connectivity-Preserving Morphological Image Transformations"
2001 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
2002 ** http://www.leptonica.com/papers/conn.pdf
2003 ** And
2004 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2005 **
2006 ** Note kernels do not specify the origin pixel, allowing them
2007 ** to be used for both thickening and thinning operations.
2008 */
2009 switch ( (int) args->rho ) {
2010 /* SE for 4-connected thinning */
2011 case 41: /* SE_4_1 */
2012 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
2013 break;
2014 case 42: /* SE_4_2 */
2015 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
2016 break;
2017 case 43: /* SE_4_3 */
2018 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
2019 break;
2020 case 44: /* SE_4_4 */
2021 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
2022 break;
2023 case 45: /* SE_4_5 */
2024 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
2025 break;
2026 case 46: /* SE_4_6 */
2027 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
2028 break;
2029 case 47: /* SE_4_7 */
2030 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
2031 break;
2032 case 48: /* SE_4_8 */
2033 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
2034 break;
2035 case 49: /* SE_4_9 */
2036 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
2037 break;
2038 /* SE for 8-connected thinning - negatives of the above */
2039 case 81: /* SE_8_0 */
2040 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
2041 break;
2042 case 82: /* SE_8_2 */
2043 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
2044 break;
2045 case 83: /* SE_8_3 */
2046 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
2047 break;
2048 case 84: /* SE_8_4 */
2049 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
2050 break;
2051 case 85: /* SE_8_5 */
2052 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2053 break;
2054 case 86: /* SE_8_6 */
2055 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2056 break;
2057 case 87: /* SE_8_7 */
2058 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2059 break;
2060 case 88: /* SE_8_8 */
2061 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2062 break;
2063 case 89: /* SE_8_9 */
2064 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2065 break;
2066 /* Special combined SE kernels */
2067 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2068 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2069 break;
2070 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2071 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2072 break;
2073 case 481: /* SE_48_1 - General Connected Corner Kernel */
2074 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2075 break;
2076 default:
2077 case 482: /* SE_48_2 - General Edge Kernel */
2078 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2079 break;
2080 }
2081 if (kernel == (KernelInfo *) NULL)
2082 return(kernel);
2083 kernel->type = type;
2084 RotateKernelInfo(kernel, args->sigma);
2085 break;
2086 }
2087 /*
2088 Distance Measuring Kernels
2089 */
2090 case ChebyshevKernel:
2091 {
2092 if (args->rho < 1.0)
2093 kernel->width = kernel->height = 3; /* default radius = 1 */
2094 else
2095 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2096 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2097
2098 kernel->values=(MagickRealType *) MagickAssumeAligned(
2099 AcquireAlignedMemory(kernel->width,kernel->height*
2100 sizeof(*kernel->values)));
2101 if (kernel->values == (MagickRealType *) NULL)
2102 return(DestroyKernelInfo(kernel));
2103
2104 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2105 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2106 kernel->positive_range += ( kernel->values[i] =
2107 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2108 kernel->maximum = kernel->values[0];
2109 break;
2110 }
2111 case ManhattanKernel:
2112 {
2113 if (args->rho < 1.0)
2114 kernel->width = kernel->height = 3; /* default radius = 1 */
2115 else
2116 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2117 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2118
2119 kernel->values=(MagickRealType *) MagickAssumeAligned(
2120 AcquireAlignedMemory(kernel->width,kernel->height*
2121 sizeof(*kernel->values)));
2122 if (kernel->values == (MagickRealType *) NULL)
2123 return(DestroyKernelInfo(kernel));
2124
2125 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2126 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2127 kernel->positive_range += ( kernel->values[i] =
2128 args->sigma*(labs((long) u)+labs((long) v)) );
2129 kernel->maximum = kernel->values[0];
2130 break;
2131 }
2132 case OctagonalKernel:
2133 {
2134 if (args->rho < 2.0)
2135 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2136 else
2137 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2138 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2139
2140 kernel->values=(MagickRealType *) MagickAssumeAligned(
2141 AcquireAlignedMemory(kernel->width,kernel->height*
2142 sizeof(*kernel->values)));
2143 if (kernel->values == (MagickRealType *) NULL)
2144 return(DestroyKernelInfo(kernel));
2145
2146 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2147 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2148 {
2149 double
2150 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2151 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2152 kernel->positive_range += kernel->values[i] =
2153 args->sigma*MagickMax(r1,r2);
2154 }
2155 kernel->maximum = kernel->values[0];
2156 break;
2157 }
2158 case EuclideanKernel:
2159 {
2160 if (args->rho < 1.0)
2161 kernel->width = kernel->height = 3; /* default radius = 1 */
2162 else
2163 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2164 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2165
2166 kernel->values=(MagickRealType *) MagickAssumeAligned(
2167 AcquireAlignedMemory(kernel->width,kernel->height*
2168 sizeof(*kernel->values)));
2169 if (kernel->values == (MagickRealType *) NULL)
2170 return(DestroyKernelInfo(kernel));
2171
2172 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2173 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2174 kernel->positive_range += ( kernel->values[i] =
2175 args->sigma*sqrt((double) (u*u+v*v)) );
2176 kernel->maximum = kernel->values[0];
2177 break;
2178 }
2179 default:
2180 {
2181 /* No-Op Kernel - Basically just a single pixel on its own */
2182 kernel=ParseKernelArray("1:1");
2183 if (kernel == (KernelInfo *) NULL)
2184 return(kernel);
2185 kernel->type = UndefinedKernel;
2186 break;
2187 }
2188 break;
2189 }
2190 return(kernel);
2191}
2192
2193/*
2194%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2195% %
2196% %
2197% %
2198% C l o n e K e r n e l I n f o %
2199% %
2200% %
2201% %
2202%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2203%
2204% CloneKernelInfo() creates a new clone of the given Kernel List so that its
2205% can be modified without effecting the original. The cloned kernel should
2206% be destroyed using DestroyKernelInfo() when no longer needed.
2207%
2208% The format of the CloneKernelInfo method is:
2209%
2210% KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2211%
2212% A description of each parameter follows:
2213%
2214% o kernel: the Morphology/Convolution kernel to be cloned
2215%
2216*/
2217MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2218{
2219 ssize_t
2220 i;
2221
2223 *new_kernel;
2224
2225 assert(kernel != (KernelInfo *) NULL);
2226 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2227 if (new_kernel == (KernelInfo *) NULL)
2228 return(new_kernel);
2229 *new_kernel=(*kernel); /* copy values in structure */
2230
2231 /* replace the values with a copy of the values */
2232 new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2233 AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2234 if (new_kernel->values == (MagickRealType *) NULL)
2235 return(DestroyKernelInfo(new_kernel));
2236 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2237 new_kernel->values[i]=kernel->values[i];
2238
2239 /* Also clone the next kernel in the kernel list */
2240 if ( kernel->next != (KernelInfo *) NULL ) {
2241 new_kernel->next = CloneKernelInfo(kernel->next);
2242 if ( new_kernel->next == (KernelInfo *) NULL )
2243 return(DestroyKernelInfo(new_kernel));
2244 }
2245
2246 return(new_kernel);
2247}
2248
2249/*
2250%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2251% %
2252% %
2253% %
2254% D e s t r o y K e r n e l I n f o %
2255% %
2256% %
2257% %
2258%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2259%
2260% DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2261% kernel.
2262%
2263% The format of the DestroyKernelInfo method is:
2264%
2265% KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2266%
2267% A description of each parameter follows:
2268%
2269% o kernel: the Morphology/Convolution kernel to be destroyed
2270%
2271*/
2272MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2273{
2274 assert(kernel != (KernelInfo *) NULL);
2275 if (kernel->next != (KernelInfo *) NULL)
2276 kernel->next=DestroyKernelInfo(kernel->next);
2277 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2278 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2279 return(kernel);
2280}
2281
2282/*
2283%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2284% %
2285% %
2286% %
2287+ E x p a n d M i r r o r K e r n e l I n f o %
2288% %
2289% %
2290% %
2291%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2292%
2293% ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2294% sequence of 90-degree rotated kernels but providing a reflected 180
2295% rotation, before the -/+ 90-degree rotations.
2296%
2297% This special rotation order produces a better, more symmetrical thinning of
2298% objects.
2299%
2300% The format of the ExpandMirrorKernelInfo method is:
2301%
2302% void ExpandMirrorKernelInfo(KernelInfo *kernel)
2303%
2304% A description of each parameter follows:
2305%
2306% o kernel: the Morphology/Convolution kernel
2307%
2308% This function is only internal to this module, as it is not finalized,
2309% especially with regard to non-orthogonal angles, and rotation of larger
2310% 2D kernels.
2311*/
2312
2313#if 0
2314static void FlopKernelInfo(KernelInfo *kernel)
2315 { /* Do a Flop by reversing each row. */
2316 size_t
2317 y;
2318 ssize_t
2319 x,r;
2320 double
2321 *k,t;
2322
2323 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2324 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2325 t=k[x], k[x]=k[r], k[r]=t;
2326
2327 kernel->x = kernel->width - kernel->x - 1;
2328 angle = fmod(angle+180.0, 360.0);
2329 }
2330#endif
2331
2332static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2333{
2335 *clone,
2336 *last;
2337
2338 last = kernel;
2339
2340 clone = CloneKernelInfo(last);
2341 if (clone == (KernelInfo *) NULL)
2342 return;
2343 RotateKernelInfo(clone, 180); /* flip */
2344 LastKernelInfo(last)->next = clone;
2345 last = clone;
2346
2347 clone = CloneKernelInfo(last);
2348 if (clone == (KernelInfo *) NULL)
2349 return;
2350 RotateKernelInfo(clone, 90); /* transpose */
2351 LastKernelInfo(last)->next = clone;
2352 last = clone;
2353
2354 clone = CloneKernelInfo(last);
2355 if (clone == (KernelInfo *) NULL)
2356 return;
2357 RotateKernelInfo(clone, 180); /* flop */
2358 LastKernelInfo(last)->next = clone;
2359
2360 return;
2361}
2362
2363/*
2364%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2365% %
2366% %
2367% %
2368+ E x p a n d R o t a t e K e r n e l I n f o %
2369% %
2370% %
2371% %
2372%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2373%
2374% ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2375% incrementally by the angle given, until the kernel repeats.
2376%
2377% WARNING: 45 degree rotations only works for 3x3 kernels.
2378% While 90 degree rotations only works for linear and square kernels
2379%
2380% The format of the ExpandRotateKernelInfo method is:
2381%
2382% void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2383%
2384% A description of each parameter follows:
2385%
2386% o kernel: the Morphology/Convolution kernel
2387%
2388% o angle: angle to rotate in degrees
2389%
2390% This function is only internal to this module, as it is not finalized,
2391% especially with regard to non-orthogonal angles, and rotation of larger
2392% 2D kernels.
2393*/
2394
2395/* Internal Routine - Return true if two kernels are the same */
2396static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2397 const KernelInfo *kernel2)
2398{
2399 size_t
2400 i;
2401
2402 /* check size and origin location */
2403 if ( kernel1->width != kernel2->width
2404 || kernel1->height != kernel2->height
2405 || kernel1->x != kernel2->x
2406 || kernel1->y != kernel2->y )
2407 return MagickFalse;
2408
2409 /* check actual kernel values */
2410 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2411 /* Test for Nan equivalence */
2412 if ( IsNaN(kernel1->values[i]) && !IsNaN(kernel2->values[i]) )
2413 return MagickFalse;
2414 if ( IsNaN(kernel2->values[i]) && !IsNaN(kernel1->values[i]) )
2415 return MagickFalse;
2416 /* Test actual values are equivalent */
2417 if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2418 return MagickFalse;
2419 }
2420
2421 return MagickTrue;
2422}
2423
2424static void ExpandRotateKernelInfo(KernelInfo *kernel,const double angle)
2425{
2427 *clone_info,
2428 *last;
2429
2430 clone_info=(KernelInfo *) NULL;
2431 last=kernel;
2432DisableMSCWarning(4127)
2433 while (1) {
2434RestoreMSCWarning
2435 clone_info=CloneKernelInfo(last);
2436 if (clone_info == (KernelInfo *) NULL)
2437 break;
2438 RotateKernelInfo(clone_info,angle);
2439 if (SameKernelInfo(kernel,clone_info) != MagickFalse)
2440 break;
2441 LastKernelInfo(last)->next=clone_info;
2442 last=clone_info;
2443 }
2444 if (clone_info != (KernelInfo *) NULL)
2445 clone_info=DestroyKernelInfo(clone_info); /* kernel repeated - junk */
2446 return;
2447}
2448
2449/*
2450%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2451% %
2452% %
2453% %
2454+ C a l c M e t a K e r n a l I n f o %
2455% %
2456% %
2457% %
2458%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2459%
2460% CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2461% using the kernel values. This should only ne used if it is not possible to
2462% calculate that meta-data in some easier way.
2463%
2464% It is important that the meta-data is correct before ScaleKernelInfo() is
2465% used to perform kernel normalization.
2466%
2467% The format of the CalcKernelMetaData method is:
2468%
2469% void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2470%
2471% A description of each parameter follows:
2472%
2473% o kernel: the Morphology/Convolution kernel to modify
2474%
2475% WARNING: Minimum and Maximum values are assumed to include zero, even if
2476% zero is not part of the kernel (as in Gaussian Derived kernels). This
2477% however is not true for flat-shaped morphological kernels.
2478%
2479% WARNING: Only the specific kernel pointed to is modified, not a list of
2480% multiple kernels.
2481%
2482% This is an internal function and not expected to be useful outside this
2483% module. This could change however.
2484*/
2485static void CalcKernelMetaData(KernelInfo *kernel)
2486{
2487 size_t
2488 i;
2489
2490 kernel->minimum = kernel->maximum = 0.0;
2491 kernel->negative_range = kernel->positive_range = 0.0;
2492 for (i=0; i < (kernel->width*kernel->height); i++)
2493 {
2494 if ( fabs(kernel->values[i]) < MagickEpsilon )
2495 kernel->values[i] = 0.0;
2496 ( kernel->values[i] < 0)
2497 ? ( kernel->negative_range += kernel->values[i] )
2498 : ( kernel->positive_range += kernel->values[i] );
2499 Minimize(kernel->minimum, kernel->values[i]);
2500 Maximize(kernel->maximum, kernel->values[i]);
2501 }
2502
2503 return;
2504}
2505
2506/*
2507%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2508% %
2509% %
2510% %
2511% M o r p h o l o g y A p p l y %
2512% %
2513% %
2514% %
2515%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2516%
2517% MorphologyApply() applies a morphological method, multiple times using
2518% a list of multiple kernels. This is the method that should be called by
2519% other 'operators' that internally use morphology operations as part of
2520% their processing.
2521%
2522% It is basically equivalent to as MorphologyImage() (see below) but without
2523% any user controls. This allows internal programs to use this method to
2524% perform a specific task without possible interference by any API user
2525% supplied settings.
2526%
2527% It is MorphologyImage() task to extract any such user controls, and
2528% pass them to this function for processing.
2529%
2530% More specifically all given kernels should already be scaled, normalised,
2531% and blended appropriately before being parred to this routine. The
2532% appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2533%
2534% The format of the MorphologyApply method is:
2535%
2536% Image *MorphologyApply(const Image *image,MorphologyMethod method,
2537% const ssize_t iterations,const KernelInfo *kernel,
2538% const CompositeMethod compose,const double bias,
2539% ExceptionInfo *exception)
2540%
2541% A description of each parameter follows:
2542%
2543% o image: the source image
2544%
2545% o method: the morphology method to be applied.
2546%
2547% o iterations: apply the operation this many times (or no change).
2548% A value of -1 means loop until no change found.
2549% How this is applied may depend on the morphology method.
2550% Typically this is a value of 1.
2551%
2552% o channel: the channel type.
2553%
2554% o kernel: An array of double representing the morphology kernel.
2555%
2556% o compose: How to handle or merge multi-kernel results.
2557% If 'UndefinedCompositeOp' use default for the Morphology method.
2558% If 'NoCompositeOp' force image to be re-iterated by each kernel.
2559% Otherwise merge the results using the compose method given.
2560%
2561% o bias: Convolution Output Bias.
2562%
2563% o exception: return any errors or warnings in this structure.
2564%
2565*/
2566static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2567 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2568 ExceptionInfo *exception)
2569{
2570#define MorphologyTag "Morphology/Image"
2571
2572 CacheView
2573 *image_view,
2574 *morphology_view;
2575
2576 MagickBooleanType
2577 status;
2578
2579 MagickOffsetType
2580 progress;
2581
2583 offset;
2584
2585 ssize_t
2586 j,
2587 y;
2588
2589 size_t
2590 changed,
2591 *changes,
2592 width;
2593
2594 /*
2595 Some methods (including convolve) needs to use a reflected kernel.
2596 Adjust 'origin' offsets to loop though kernel as a reflection.
2597 */
2598 assert(image != (Image *) NULL);
2599 assert(image->signature == MagickCoreSignature);
2600 assert(morphology_image != (Image *) NULL);
2601 assert(morphology_image->signature == MagickCoreSignature);
2602 assert(kernel != (KernelInfo *) NULL);
2603 assert(kernel->signature == MagickCoreSignature);
2604 assert(exception != (ExceptionInfo *) NULL);
2605 assert(exception->signature == MagickCoreSignature);
2606 status=MagickTrue;
2607 progress=0;
2608 image_view=AcquireVirtualCacheView(image,exception);
2609 morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2610 width=image->columns+kernel->width-1;
2611 offset.x=0;
2612 offset.y=0;
2613 switch (method)
2614 {
2615 case ConvolveMorphology:
2616 case DilateMorphology:
2617 case DilateIntensityMorphology:
2618 case IterativeDistanceMorphology:
2619 {
2620 /*
2621 Kernel needs to use a reflection about origin.
2622 */
2623 offset.x=(ssize_t) kernel->width-kernel->x-1;
2624 offset.y=(ssize_t) kernel->height-kernel->y-1;
2625 break;
2626 }
2627 case ErodeMorphology:
2628 case ErodeIntensityMorphology:
2629 case HitAndMissMorphology:
2630 case ThinningMorphology:
2631 case ThickenMorphology:
2632 {
2633 /*
2634 Use kernel as is, not reflection required.
2635 */
2636 offset.x=kernel->x;
2637 offset.y=kernel->y;
2638 break;
2639 }
2640 default:
2641 {
2642 ThrowMagickException(exception,GetMagickModule(),OptionWarning,
2643 "InvalidOption","`%s'","not a primitive morphology method");
2644 break;
2645 }
2646 }
2647 changed=0;
2648 changes=(size_t *) AcquireQuantumMemory(GetOpenMPMaximumThreads(),
2649 sizeof(*changes));
2650 if (changes == (size_t *) NULL)
2651 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2652 for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
2653 changes[j]=0;
2654 if ((method == ConvolveMorphology) && (kernel->width == 1))
2655 {
2656 ssize_t
2657 x;
2658
2659 /*
2660 Special handling (for speed) of vertical (blur) kernels. This performs
2661 its handling in columns rather than in rows. This is only done
2662 for convolve as it is the only method that generates very large 1-D
2663 vertical kernels (such as a 'BlurKernel')
2664 */
2665#if defined(MAGICKCORE_OPENMP_SUPPORT)
2666 #pragma omp parallel for schedule(static) shared(progress,status) \
2667 magick_number_threads(image,morphology_image,image->columns,1)
2668#endif
2669 for (x=0; x < (ssize_t) image->columns; x++)
2670 {
2671 const int
2672 id = GetOpenMPThreadId();
2673
2674 const Quantum
2675 *magick_restrict p;
2676
2677 Quantum
2678 *magick_restrict q;
2679
2680 ssize_t
2681 center,
2682 r;
2683
2684 if (status == MagickFalse)
2685 continue;
2686 p=GetCacheViewVirtualPixels(image_view,x,-offset.y,1,image->rows+
2687 kernel->height-1,exception);
2688 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2689 morphology_image->rows,exception);
2690 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2691 {
2692 status=MagickFalse;
2693 continue;
2694 }
2695 center=(ssize_t) GetPixelChannels(image)*offset.y;
2696 for (r=0; r < (ssize_t) image->rows; r++)
2697 {
2698 ssize_t
2699 i;
2700
2701 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2702 {
2703 double
2704 alpha,
2705 gamma,
2706 pixel;
2707
2708 PixelChannel
2709 channel;
2710
2711 PixelTrait
2712 morphology_traits,
2713 traits;
2714
2715 const MagickRealType
2716 *magick_restrict k;
2717
2718 const Quantum
2719 *magick_restrict pixels;
2720
2721 ssize_t
2722 v;
2723
2724 size_t
2725 count;
2726
2727 channel=GetPixelChannelChannel(image,i);
2728 traits=GetPixelChannelTraits(image,channel);
2729 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2730 if ((traits == UndefinedPixelTrait) ||
2731 (morphology_traits == UndefinedPixelTrait))
2732 continue;
2733 if ((traits & CopyPixelTrait) != 0)
2734 {
2735 SetPixelChannel(morphology_image,channel,p[center+i],q);
2736 continue;
2737 }
2738 k=(&kernel->values[kernel->height-1]);
2739 pixels=p;
2740 pixel=bias;
2741 gamma=1.0;
2742 count=0;
2743 if (((image->alpha_trait & BlendPixelTrait) == 0) ||
2744 ((morphology_traits & BlendPixelTrait) == 0))
2745 for (v=0; v < (ssize_t) kernel->height; v++)
2746 {
2747 if (!IsNaN(*k))
2748 {
2749 pixel+=(*k)*(double) pixels[i];
2750 count++;
2751 }
2752 k--;
2753 pixels+=GetPixelChannels(image);
2754 }
2755 else
2756 {
2757 gamma=0.0;
2758 for (v=0; v < (ssize_t) kernel->height; v++)
2759 {
2760 if (!IsNaN(*k))
2761 {
2762 alpha=(double) (QuantumScale*(double)
2763 GetPixelAlpha(image,pixels));
2764 pixel+=alpha*(*k)*(double) pixels[i];
2765 gamma+=alpha*(*k);
2766 count++;
2767 }
2768 k--;
2769 pixels+=GetPixelChannels(image);
2770 }
2771 }
2772 if (fabs(pixel-(double) p[center+i]) >= MagickEpsilon)
2773 changes[id]++;
2774 gamma=PerceptibleReciprocal(gamma);
2775 if (count != 0)
2776 gamma*=(double) kernel->height/count;
2777 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2778 pixel),q);
2779 }
2780 p+=GetPixelChannels(image);
2781 q+=GetPixelChannels(morphology_image);
2782 }
2783 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2784 status=MagickFalse;
2785 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2786 {
2787 MagickBooleanType
2788 proceed;
2789
2790#if defined(MAGICKCORE_OPENMP_SUPPORT)
2791 #pragma omp atomic
2792#endif
2793 progress++;
2794 proceed=SetImageProgress(image,MorphologyTag,progress,
2795 image->columns);
2796 if (proceed == MagickFalse)
2797 status=MagickFalse;
2798 }
2799 }
2800 morphology_image->type=image->type;
2801 morphology_view=DestroyCacheView(morphology_view);
2802 image_view=DestroyCacheView(image_view);
2803 for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
2804 changed+=changes[j];
2805 changes=(size_t *) RelinquishMagickMemory(changes);
2806 return(status ? (ssize_t) (changed/GetImageChannels(image)) : 0);
2807 }
2808 /*
2809 Normal handling of horizontal or rectangular kernels (row by row).
2810 */
2811#if defined(MAGICKCORE_OPENMP_SUPPORT)
2812 #pragma omp parallel for schedule(static) shared(progress,status) \
2813 magick_number_threads(image,morphology_image,image->rows,1)
2814#endif
2815 for (y=0; y < (ssize_t) image->rows; y++)
2816 {
2817 const int
2818 id = GetOpenMPThreadId();
2819
2820 const Quantum
2821 *magick_restrict p;
2822
2823 Quantum
2824 *magick_restrict q;
2825
2826 ssize_t
2827 x;
2828
2829 ssize_t
2830 center;
2831
2832 if (status == MagickFalse)
2833 continue;
2834 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,
2835 kernel->height,exception);
2836 q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2837 1,exception);
2838 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2839 {
2840 status=MagickFalse;
2841 continue;
2842 }
2843 center=(ssize_t) ((ssize_t) GetPixelChannels(image)*(ssize_t) width*
2844 offset.y+(ssize_t) GetPixelChannels(image)*offset.x);
2845 for (x=0; x < (ssize_t) image->columns; x++)
2846 {
2847 ssize_t
2848 i;
2849
2850 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2851 {
2852 double
2853 alpha,
2854 gamma,
2855 intensity,
2856 maximum,
2857 minimum,
2858 pixel;
2859
2860 PixelChannel
2861 channel;
2862
2863 PixelTrait
2864 morphology_traits,
2865 traits;
2866
2867 const MagickRealType
2868 *magick_restrict k;
2869
2870 const Quantum
2871 *magick_restrict pixels,
2872 *magick_restrict quantum_pixels;
2873
2874 ssize_t
2875 u;
2876
2877 ssize_t
2878 v;
2879
2880 channel=GetPixelChannelChannel(image,i);
2881 traits=GetPixelChannelTraits(image,channel);
2882 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2883 if ((traits == UndefinedPixelTrait) ||
2884 (morphology_traits == UndefinedPixelTrait))
2885 continue;
2886 if ((traits & CopyPixelTrait) != 0)
2887 {
2888 SetPixelChannel(morphology_image,channel,p[center+i],q);
2889 continue;
2890 }
2891 pixels=p;
2892 quantum_pixels=(const Quantum *) NULL;
2893 maximum=0.0;
2894 minimum=(double) QuantumRange;
2895 switch (method)
2896 {
2897 case ConvolveMorphology:
2898 {
2899 pixel=bias;
2900 break;
2901 }
2902 case DilateMorphology:
2903 case ErodeIntensityMorphology:
2904 {
2905 pixel=0.0;
2906 break;
2907 }
2908 default:
2909 {
2910 pixel=(double) p[center+i];
2911 break;
2912 }
2913 }
2914 gamma=1.0;
2915 switch (method)
2916 {
2917 case ConvolveMorphology:
2918 {
2919 /*
2920 Weighted Average of pixels using reflected kernel
2921
2922 For correct working of this operation for asymmetrical kernels,
2923 the kernel needs to be applied in its reflected form. That is
2924 its values needs to be reversed.
2925
2926 Correlation is actually the same as this but without reflecting
2927 the kernel, and thus 'lower-level' that Convolution. However as
2928 Convolution is the more common method used, and it does not
2929 really cost us much in terms of processing to use a reflected
2930 kernel, so it is Convolution that is implemented.
2931
2932 Correlation will have its kernel reflected before calling this
2933 function to do a Convolve.
2934
2935 For more details of Correlation vs Convolution see
2936 http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2937 */
2938 k=(&kernel->values[kernel->width*kernel->height-1]);
2939 if (((image->alpha_trait & BlendPixelTrait) == 0) ||
2940 ((morphology_traits & BlendPixelTrait) == 0))
2941 {
2942 /*
2943 No alpha blending.
2944 */
2945 for (v=0; v < (ssize_t) kernel->height; v++)
2946 {
2947 for (u=0; u < (ssize_t) kernel->width; u++)
2948 {
2949 if (!IsNaN(*k))
2950 pixel+=(*k)*(double) pixels[i];
2951 k--;
2952 pixels+=GetPixelChannels(image);
2953 }
2954 pixels+=(image->columns-1)*GetPixelChannels(image);
2955 }
2956 break;
2957 }
2958 /*
2959 Alpha blending.
2960 */
2961 gamma=0.0;
2962 for (v=0; v < (ssize_t) kernel->height; v++)
2963 {
2964 for (u=0; u < (ssize_t) kernel->width; u++)
2965 {
2966 if (!IsNaN(*k))
2967 {
2968 alpha=(double) (QuantumScale*(double)
2969 GetPixelAlpha(image,pixels));
2970 pixel+=alpha*(*k)*(double) pixels[i];
2971 gamma+=alpha*(*k);
2972 }
2973 k--;
2974 pixels+=GetPixelChannels(image);
2975 }
2976 pixels+=(image->columns-1)*GetPixelChannels(image);
2977 }
2978 break;
2979 }
2980 case ErodeMorphology:
2981 {
2982 /*
2983 Minimum value within kernel neighbourhood.
2984
2985 The kernel is not reflected for this operation. In normal
2986 Greyscale Morphology, the kernel value should be added
2987 to the real value, this is currently not done, due to the
2988 nature of the boolean kernels being used.
2989 */
2990 k=kernel->values;
2991 for (v=0; v < (ssize_t) kernel->height; v++)
2992 {
2993 for (u=0; u < (ssize_t) kernel->width; u++)
2994 {
2995 if (!IsNaN(*k) && (*k >= 0.5))
2996 {
2997 if ((double) pixels[i] < pixel)
2998 pixel=(double) pixels[i];
2999 }
3000 k++;
3001 pixels+=GetPixelChannels(image);
3002 }
3003 pixels+=(image->columns-1)*GetPixelChannels(image);
3004 }
3005 break;
3006 }
3007 case DilateMorphology:
3008 {
3009 /*
3010 Maximum value within kernel neighbourhood.
3011
3012 For correct working of this operation for asymmetrical kernels,
3013 the kernel needs to be applied in its reflected form. That is
3014 its values needs to be reversed.
3015
3016 In normal Greyscale Morphology, the kernel value should be
3017 added to the real value, this is currently not done, due to the
3018 nature of the boolean kernels being used.
3019 */
3020 k=(&kernel->values[kernel->width*kernel->height-1]);
3021 for (v=0; v < (ssize_t) kernel->height; v++)
3022 {
3023 for (u=0; u < (ssize_t) kernel->width; u++)
3024 {
3025 if (!IsNaN(*k) && (*k > 0.5))
3026 {
3027 if ((double) pixels[i] > pixel)
3028 pixel=(double) pixels[i];
3029 }
3030 k--;
3031 pixels+=GetPixelChannels(image);
3032 }
3033 pixels+=(image->columns-1)*GetPixelChannels(image);
3034 }
3035 break;
3036 }
3037 case HitAndMissMorphology:
3038 case ThinningMorphology:
3039 case ThickenMorphology:
3040 {
3041 /*
3042 Minimum of foreground pixel minus maximum of background pixels.
3043
3044 The kernel is not reflected for this operation, and consists
3045 of both foreground and background pixel neighbourhoods, 0.0 for
3046 background, and 1.0 for foreground with either Nan or 0.5 values
3047 for don't care.
3048
3049 This never produces a meaningless negative result. Such results
3050 cause Thinning/Thicken to not work correctly when used against a
3051 greyscale image.
3052 */
3053 k=kernel->values;
3054 for (v=0; v < (ssize_t) kernel->height; v++)
3055 {
3056 for (u=0; u < (ssize_t) kernel->width; u++)
3057 {
3058 if (!IsNaN(*k))
3059 {
3060 if (*k > 0.7)
3061 {
3062 if ((double) pixels[i] < minimum)
3063 minimum=(double) pixels[i];
3064 }
3065 else
3066 if (*k < 0.3)
3067 {
3068 if ((double) pixels[i] > maximum)
3069 maximum=(double) pixels[i];
3070 }
3071 }
3072 k++;
3073 pixels+=GetPixelChannels(image);
3074 }
3075 pixels+=(image->columns-1)*GetPixelChannels(image);
3076 }
3077 minimum-=maximum;
3078 if (minimum < 0.0)
3079 minimum=0.0;
3080 pixel=minimum;
3081 if (method == ThinningMorphology)
3082 pixel=(double) p[center+i]-minimum;
3083 else
3084 if (method == ThickenMorphology)
3085 pixel=(double) p[center+i]+minimum;
3086 break;
3087 }
3088 case ErodeIntensityMorphology:
3089 {
3090 /*
3091 Select pixel with minimum intensity within kernel neighbourhood.
3092
3093 The kernel is not reflected for this operation.
3094 */
3095 k=kernel->values;
3096 for (v=0; v < (ssize_t) kernel->height; v++)
3097 {
3098 for (u=0; u < (ssize_t) kernel->width; u++)
3099 {
3100 if (!IsNaN(*k) && (*k >= 0.5))
3101 {
3102 intensity=(double) GetPixelIntensity(image,pixels);
3103 if (intensity < minimum)
3104 {
3105 quantum_pixels=pixels;
3106 pixel=(double) pixels[i];
3107 minimum=intensity;
3108 }
3109 }
3110 k++;
3111 pixels+=GetPixelChannels(image);
3112 }
3113 pixels+=(image->columns-1)*GetPixelChannels(image);
3114 }
3115 break;
3116 }
3117 case DilateIntensityMorphology:
3118 {
3119 /*
3120 Select pixel with maximum intensity within kernel neighbourhood.
3121
3122 The kernel is not reflected for this operation.
3123 */
3124 k=(&kernel->values[kernel->width*kernel->height-1]);
3125 for (v=0; v < (ssize_t) kernel->height; v++)
3126 {
3127 for (u=0; u < (ssize_t) kernel->width; u++)
3128 {
3129 if (!IsNaN(*k) && (*k >= 0.5))
3130 {
3131 intensity=(double) GetPixelIntensity(image,pixels);
3132 if (intensity > maximum)
3133 {
3134 pixel=(double) pixels[i];
3135 quantum_pixels=pixels;
3136 maximum=intensity;
3137 }
3138 }
3139 k--;
3140 pixels+=GetPixelChannels(image);
3141 }
3142 pixels+=(image->columns-1)*GetPixelChannels(image);
3143 }
3144 break;
3145 }
3146 case IterativeDistanceMorphology:
3147 {
3148 /*
3149 Compute th iterative distance from black edge of a white image
3150 shape. Essentially white values are decreased to the smallest
3151 'distance from edge' it can find.
3152
3153 It works by adding kernel values to the neighbourhood, and
3154 select the minimum value found. The kernel is rotated before
3155 use, so kernel distances match resulting distances, when a user
3156 provided asymmetric kernel is applied.
3157
3158 This code is nearly identical to True GrayScale Morphology but
3159 not quite.
3160
3161 GreyDilate Kernel values added, maximum value found Kernel is
3162 rotated before use.
3163
3164 GrayErode: Kernel values subtracted and minimum value found No
3165 kernel rotation used.
3166
3167 Note the Iterative Distance method is essentially a
3168 GrayErode, but with negative kernel values, and kernel rotation
3169 applied.
3170 */
3171 k=(&kernel->values[kernel->width*kernel->height-1]);
3172 for (v=0; v < (ssize_t) kernel->height; v++)
3173 {
3174 for (u=0; u < (ssize_t) kernel->width; u++)
3175 {
3176 if (!IsNaN(*k))
3177 {
3178 if (((double) pixels[i]+(*k)) < pixel)
3179 pixel=(double) pixels[i]+(*k);
3180 }
3181 k--;
3182 pixels+=GetPixelChannels(image);
3183 }
3184 pixels+=(image->columns-1)*GetPixelChannels(image);
3185 }
3186 break;
3187 }
3188 case UndefinedMorphology:
3189 default:
3190 break;
3191 }
3192 if (quantum_pixels != (const Quantum *) NULL)
3193 {
3194 SetPixelChannel(morphology_image,channel,quantum_pixels[i],q);
3195 continue;
3196 }
3197 gamma=PerceptibleReciprocal(gamma);
3198 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*pixel),q);
3199 if (fabs(pixel-(double) p[center+i]) >= MagickEpsilon)
3200 changes[id]++;
3201 }
3202 p+=GetPixelChannels(image);
3203 q+=GetPixelChannels(morphology_image);
3204 }
3205 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3206 status=MagickFalse;
3207 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3208 {
3209 MagickBooleanType
3210 proceed;
3211
3212#if defined(MAGICKCORE_OPENMP_SUPPORT)
3213 #pragma omp atomic
3214#endif
3215 progress++;
3216 proceed=SetImageProgress(image,MorphologyTag,progress,image->rows);
3217 if (proceed == MagickFalse)
3218 status=MagickFalse;
3219 }
3220 }
3221 morphology_view=DestroyCacheView(morphology_view);
3222 image_view=DestroyCacheView(image_view);
3223 for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
3224 changed+=changes[j];
3225 changes=(size_t *) RelinquishMagickMemory(changes);
3226 return(status ? (ssize_t) (changed/GetImageChannels(image)) : -1);
3227}
3228
3229/*
3230 This is almost identical to the MorphologyPrimitive() function above, but
3231 applies the primitive directly to the actual image using two passes, once in
3232 each direction, with the results of the previous (and current) row being
3233 re-used.
3234
3235 That is after each row is 'Sync'ed' into the image, the next row makes use of
3236 those values as part of the calculation of the next row. It repeats, but
3237 going in the opposite (bottom-up) direction.
3238
3239 Because of this 're-use of results' this function can not make use of multi-
3240 threaded, parallel processing.
3241*/
3242static ssize_t MorphologyPrimitiveDirect(Image *image,
3243 const MorphologyMethod method,const KernelInfo *kernel,
3244 ExceptionInfo *exception)
3245{
3246 CacheView
3247 *morphology_view,
3248 *image_view;
3249
3250 MagickBooleanType
3251 status;
3252
3253 MagickOffsetType
3254 progress;
3255
3257 offset;
3258
3259 size_t
3260 width,
3261 changed;
3262
3263 ssize_t
3264 y;
3265
3266 assert(image != (Image *) NULL);
3267 assert(image->signature == MagickCoreSignature);
3268 assert(kernel != (KernelInfo *) NULL);
3269 assert(kernel->signature == MagickCoreSignature);
3270 assert(exception != (ExceptionInfo *) NULL);
3271 assert(exception->signature == MagickCoreSignature);
3272 status=MagickTrue;
3273 changed=0;
3274 progress=0;
3275 switch(method)
3276 {
3277 case DistanceMorphology:
3278 case VoronoiMorphology:
3279 {
3280 /*
3281 Kernel reflected about origin.
3282 */
3283 offset.x=(ssize_t) kernel->width-kernel->x-1;
3284 offset.y=(ssize_t) kernel->height-kernel->y-1;
3285 break;
3286 }
3287 default:
3288 {
3289 offset.x=kernel->x;
3290 offset.y=kernel->y;
3291 break;
3292 }
3293 }
3294 /*
3295 Two views into same image, do not thread.
3296 */
3297 image_view=AcquireVirtualCacheView(image,exception);
3298 morphology_view=AcquireAuthenticCacheView(image,exception);
3299 width=image->columns+kernel->width-1;
3300 for (y=0; y < (ssize_t) image->rows; y++)
3301 {
3302 const Quantum
3303 *magick_restrict p;
3304
3305 Quantum
3306 *magick_restrict q;
3307
3308 ssize_t
3309 x;
3310
3311 /*
3312 Read virtual pixels, and authentic pixels, from the same image! We read
3313 using virtual to get virtual pixel handling, but write back into the same
3314 image.
3315
3316 Only top half of kernel is processed as we do a single pass downward
3317 through the image iterating the distance function as we go.
3318 */
3319 if (status == MagickFalse)
3320 continue;
3321 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t)
3322 offset.y+1,exception);
3323 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3324 exception);
3325 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3326 {
3327 status=MagickFalse;
3328 continue;
3329 }
3330 for (x=0; x < (ssize_t) image->columns; x++)
3331 {
3332 ssize_t
3333 i;
3334
3335 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3336 {
3337 double
3338 pixel;
3339
3340 PixelChannel
3341 channel;
3342
3343 PixelTrait
3344 traits;
3345
3346 const MagickRealType
3347 *magick_restrict k;
3348
3349 const Quantum
3350 *magick_restrict pixels;
3351
3352 ssize_t
3353 u;
3354
3355 ssize_t
3356 v;
3357
3358 channel=GetPixelChannelChannel(image,i);
3359 traits=GetPixelChannelTraits(image,channel);
3360 if (traits == UndefinedPixelTrait)
3361 continue;
3362 if ((traits & CopyPixelTrait) != 0)
3363 continue;
3364 pixels=p;
3365 pixel=(double) QuantumRange;
3366 switch (method)
3367 {
3368 case DistanceMorphology:
3369 {
3370 k=(&kernel->values[kernel->width*kernel->height-1]);
3371 for (v=0; v <= offset.y; v++)
3372 {
3373 for (u=0; u < (ssize_t) kernel->width; u++)
3374 {
3375 if (!IsNaN(*k))
3376 {
3377 if (((double) pixels[i]+(*k)) < pixel)
3378 pixel=(double) pixels[i]+(*k);
3379 }
3380 k--;
3381 pixels+=GetPixelChannels(image);
3382 }
3383 pixels+=(image->columns-1)*GetPixelChannels(image);
3384 }
3385 k=(&kernel->values[(ssize_t) kernel->width*(kernel->y+1)-1]);
3386 pixels=q-offset.x*(ssize_t) GetPixelChannels(image);
3387 for (u=0; u < offset.x; u++)
3388 {
3389 if (!IsNaN(*k) && ((x+u-offset.x) >= 0))
3390 {
3391 if (((double) pixels[i]+(*k)) < pixel)
3392 pixel=(double) pixels[i]+(*k);
3393 }
3394 k--;
3395 pixels+=GetPixelChannels(image);
3396 }
3397 break;
3398 }
3399 case VoronoiMorphology:
3400 {
3401 k=(&kernel->values[kernel->width*kernel->height-1]);
3402 for (v=0; v < offset.y; v++)
3403 {
3404 for (u=0; u < (ssize_t) kernel->width; u++)
3405 {
3406 if (!IsNaN(*k))
3407 {
3408 if (((double) pixels[i]+(*k)) < pixel)
3409 pixel=(double) pixels[i]+(*k);
3410 }
3411 k--;
3412 pixels+=GetPixelChannels(image);
3413 }
3414 pixels+=(image->columns-1)*GetPixelChannels(image);
3415 }
3416 k=(&kernel->values[(ssize_t) kernel->width*(kernel->y+1)-1]);
3417 pixels=q-offset.x*(ssize_t) GetPixelChannels(image);
3418 for (u=0; u < offset.x; u++)
3419 {
3420 if (!IsNaN(*k) && ((x+u-offset.x) >= 0))
3421 {
3422 if (((double) pixels[i]+(*k)) < pixel)
3423 pixel=(double) pixels[i]+(*k);
3424 }
3425 k--;
3426 pixels+=GetPixelChannels(image);
3427 }
3428 break;
3429 }
3430 default:
3431 break;
3432 }
3433 if (fabs(pixel-(double) q[i]) > MagickEpsilon)
3434 changed++;
3435 q[i]=ClampToQuantum(pixel);
3436 }
3437 p+=GetPixelChannels(image);
3438 q+=GetPixelChannels(image);
3439 }
3440 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3441 status=MagickFalse;
3442 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3443 {
3444 MagickBooleanType
3445 proceed;
3446
3447#if defined(MAGICKCORE_OPENMP_SUPPORT)
3448 #pragma omp atomic
3449#endif
3450 progress++;
3451 proceed=SetImageProgress(image,MorphologyTag,progress,2*image->rows);
3452 if (proceed == MagickFalse)
3453 status=MagickFalse;
3454 }
3455 }
3456 morphology_view=DestroyCacheView(morphology_view);
3457 image_view=DestroyCacheView(image_view);
3458 /*
3459 Do the reverse pass through the image.
3460 */
3461 image_view=AcquireVirtualCacheView(image,exception);
3462 morphology_view=AcquireAuthenticCacheView(image,exception);
3463 for (y=(ssize_t) image->rows-1; y >= 0; y--)
3464 {
3465 const Quantum
3466 *magick_restrict p;
3467
3468 Quantum
3469 *magick_restrict q;
3470
3471 ssize_t
3472 x;
3473
3474 /*
3475 Read virtual pixels, and authentic pixels, from the same image. We
3476 read using virtual to get virtual pixel handling, but write back
3477 into the same image.
3478
3479 Only the bottom half of the kernel is processed as we up the image.
3480 */
3481 if (status == MagickFalse)
3482 continue;
3483 p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t)
3484 kernel->y+1,exception);
3485 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3486 exception);
3487 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3488 {
3489 status=MagickFalse;
3490 continue;
3491 }
3492 p+=(image->columns-1)*GetPixelChannels(image);
3493 q+=(image->columns-1)*GetPixelChannels(image);
3494 for (x=(ssize_t) image->columns-1; x >= 0; x--)
3495 {
3496 ssize_t
3497 i;
3498
3499 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3500 {
3501 double
3502 pixel;
3503
3504 PixelChannel
3505 channel;
3506
3507 PixelTrait
3508 traits;
3509
3510 const MagickRealType
3511 *magick_restrict k;
3512
3513 const Quantum
3514 *magick_restrict pixels;
3515
3516 ssize_t
3517 u;
3518
3519 ssize_t
3520 v;
3521
3522 channel=GetPixelChannelChannel(image,i);
3523 traits=GetPixelChannelTraits(image,channel);
3524 if (traits == UndefinedPixelTrait)
3525 continue;
3526 if ((traits & CopyPixelTrait) != 0)
3527 continue;
3528 pixels=p;
3529 pixel=(double) QuantumRange;
3530 switch (method)
3531 {
3532 case DistanceMorphology:
3533 {
3534 k=(&kernel->values[(ssize_t) kernel->width*(kernel->y+1)-1]);
3535 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3536 {
3537 for (u=0; u < (ssize_t) kernel->width; u++)
3538 {
3539 if (!IsNaN(*k))
3540 {
3541 if (((double) pixels[i]+(*k)) < pixel)
3542 pixel=(double) pixels[i]+(*k);
3543 }
3544 k--;
3545 pixels+=GetPixelChannels(image);
3546 }
3547 pixels+=(image->columns-1)*GetPixelChannels(image);
3548 }
3549 k=(&kernel->values[(ssize_t) kernel->width*kernel->y+kernel->x-1]);
3550 pixels=q;
3551 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3552 {
3553 pixels+=GetPixelChannels(image);
3554 if (!IsNaN(*k) && ((x+u-offset.x) < (ssize_t) image->columns))
3555 {
3556 if (((double) pixels[i]+(*k)) < pixel)
3557 pixel=(double) pixels[i]+(*k);
3558 }
3559 k--;
3560 }
3561 break;
3562 }
3563 case VoronoiMorphology:
3564 {
3565 k=(&kernel->values[(ssize_t) kernel->width*(kernel->y+1)-1]);
3566 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3567 {
3568 for (u=0; u < (ssize_t) kernel->width; u++)
3569 {
3570 if (!IsNaN(*k))
3571 {
3572 if (((double) pixels[i]+(*k)) < pixel)
3573 pixel=(double) pixels[i]+(*k);
3574 }
3575 k--;
3576 pixels+=GetPixelChannels(image);
3577 }
3578 pixels+=(image->columns-1)*GetPixelChannels(image);
3579 }
3580 k=(&kernel->values[(ssize_t) kernel->width*(kernel->y+1)-1]);
3581 pixels=q;
3582 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3583 {
3584 pixels+=GetPixelChannels(image);
3585 if (!IsNaN(*k) && ((x+u-offset.x) < (ssize_t) image->columns))
3586 {
3587 if (((double) pixels[i]+(*k)) < pixel)
3588 pixel=(double) pixels[i]+(*k);
3589 }
3590 k--;
3591 }
3592 break;
3593 }
3594 default:
3595 break;
3596 }
3597 if (fabs(pixel-(double) q[i]) > MagickEpsilon)
3598 changed++;
3599 q[i]=ClampToQuantum(pixel);
3600 }
3601 p-=GetPixelChannels(image);
3602 q-=GetPixelChannels(image);
3603 }
3604 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3605 status=MagickFalse;
3606 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3607 {
3608 MagickBooleanType
3609 proceed;
3610
3611#if defined(MAGICKCORE_OPENMP_SUPPORT)
3612 #pragma omp atomic
3613#endif
3614 progress++;
3615 proceed=SetImageProgress(image,MorphologyTag,progress,2*image->rows);
3616 if (proceed == MagickFalse)
3617 status=MagickFalse;
3618 }
3619 }
3620 morphology_view=DestroyCacheView(morphology_view);
3621 image_view=DestroyCacheView(image_view);
3622 return(status ? (ssize_t) (changed/GetImageChannels(image)) : -1);
3623}
3624
3625/*
3626 Apply a Morphology by calling one of the above low level primitive
3627 application functions. This function handles any iteration loops,
3628 composition or re-iteration of results, and compound morphology methods that
3629 is based on multiple low-level (staged) morphology methods.
3630
3631 Basically this provides the complex glue between the requested morphology
3632 method and raw low-level implementation (above).
3633*/
3634MagickPrivate Image *MorphologyApply(const Image *image,
3635 const MorphologyMethod method, const ssize_t iterations,
3636 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3637 ExceptionInfo *exception)
3638{
3639 CompositeOperator
3640 curr_compose;
3641
3642 Image
3643 *curr_image, /* Image we are working with or iterating */
3644 *work_image, /* secondary image for primitive iteration */
3645 *save_image, /* saved image - for 'edge' method only */
3646 *rslt_image; /* resultant image - after multi-kernel handling */
3647
3649 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3650 *norm_kernel, /* the current normal un-reflected kernel */
3651 *rflt_kernel, /* the current reflected kernel (if needed) */
3652 *this_kernel; /* the kernel being applied */
3653
3654 MorphologyMethod
3655 primitive; /* the current morphology primitive being applied */
3656
3657 CompositeOperator
3658 rslt_compose; /* multi-kernel compose method for results to use */
3659
3660 MagickBooleanType
3661 special, /* do we use a direct modify function? */
3662 verbose; /* verbose output of results */
3663
3664 size_t
3665 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3666 method_limit, /* maximum number of compound method iterations */
3667 kernel_number, /* Loop 2: the kernel number being applied */
3668 stage_loop, /* Loop 3: primitive loop for compound morphology */
3669 stage_limit, /* how many primitives are in this compound */
3670 kernel_loop, /* Loop 4: iterate the kernel over image */
3671 kernel_limit, /* number of times to iterate kernel */
3672 count, /* total count of primitive steps applied */
3673 kernel_changed, /* total count of changed using iterated kernel */
3674 method_changed; /* total count of changed over method iteration */
3675
3676 ssize_t
3677 changed; /* number pixels changed by last primitive operation */
3678
3679 char
3680 v_info[MagickPathExtent];
3681
3682 assert(image != (Image *) NULL);
3683 assert(image->signature == MagickCoreSignature);
3684 assert(kernel != (KernelInfo *) NULL);
3685 assert(kernel->signature == MagickCoreSignature);
3686 assert(exception != (ExceptionInfo *) NULL);
3687 assert(exception->signature == MagickCoreSignature);
3688
3689 count = 0; /* number of low-level morphology primitives performed */
3690 if ( iterations == 0 )
3691 return((Image *) NULL); /* null operation - nothing to do! */
3692
3693 kernel_limit = (size_t) iterations;
3694 if ( iterations < 0 ) /* negative interactions = infinite (well almost) */
3695 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3696
3697 verbose = IsStringTrue(GetImageArtifact(image,"debug"));
3698
3699 /* initialise for cleanup */
3700 curr_image = (Image *) image;
3701 curr_compose = image->compose;
3702 (void) curr_compose;
3703 work_image = save_image = rslt_image = (Image *) NULL;
3704 reflected_kernel = (KernelInfo *) NULL;
3705
3706 /* Initialize specific methods
3707 * + which loop should use the given iterations
3708 * + how many primitives make up the compound morphology
3709 * + multi-kernel compose method to use (by default)
3710 */
3711 method_limit = 1; /* just do method once, unless otherwise set */
3712 stage_limit = 1; /* assume method is not a compound */
3713 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3714 rslt_compose = compose; /* and we are composing multi-kernels as given */
3715 switch( method ) {
3716 case SmoothMorphology: /* 4 primitive compound morphology */
3717 stage_limit = 4;
3718 break;
3719 case OpenMorphology: /* 2 primitive compound morphology */
3720 case OpenIntensityMorphology:
3721 case TopHatMorphology:
3722 case CloseMorphology:
3723 case CloseIntensityMorphology:
3724 case BottomHatMorphology:
3725 case EdgeMorphology:
3726 stage_limit = 2;
3727 break;
3728 case HitAndMissMorphology:
3729 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3730 magick_fallthrough;
3731 case ThinningMorphology:
3732 case ThickenMorphology:
3733 method_limit = kernel_limit; /* iterate the whole method */
3734 kernel_limit = 1; /* do not do kernel iteration */
3735 break;
3736 case DistanceMorphology:
3737 case VoronoiMorphology:
3738 special = MagickTrue; /* use special direct primitive */
3739 break;
3740 default:
3741 break;
3742 }
3743
3744 /* Apply special methods with special requirements
3745 ** For example, single run only, or post-processing requirements
3746 */
3747 if ( special != MagickFalse )
3748 {
3749 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3750 if (rslt_image == (Image *) NULL)
3751 goto error_cleanup;
3752 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3753 goto error_cleanup;
3754
3755 changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception);
3756
3757 if (verbose != MagickFalse)
3758 (void) (void) FormatLocaleFile(stderr,
3759 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3760 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3761 1.0,0.0,1.0, (double) changed);
3762
3763 if ( changed < 0 )
3764 goto error_cleanup;
3765
3766 if ( method == VoronoiMorphology ) {
3767 /* Preserve the alpha channel of input image - but turned it off */
3768 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3769 exception);
3770 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3771 MagickTrue,0,0,exception);
3772 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3773 exception);
3774 }
3775 goto exit_cleanup;
3776 }
3777
3778 /* Handle user (caller) specified multi-kernel composition method */
3779 if ( compose != UndefinedCompositeOp )
3780 rslt_compose = compose; /* override default composition for method */
3781 if ( rslt_compose == UndefinedCompositeOp )
3782 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3783
3784 /* Some methods require a reflected kernel to use with primitives.
3785 * Create the reflected kernel for those methods. */
3786 switch ( method ) {
3787 case CorrelateMorphology:
3788 case CloseMorphology:
3789 case CloseIntensityMorphology:
3790 case BottomHatMorphology:
3791 case SmoothMorphology:
3792 reflected_kernel = CloneKernelInfo(kernel);
3793 if (reflected_kernel == (KernelInfo *) NULL)
3794 goto error_cleanup;
3795 RotateKernelInfo(reflected_kernel,180);
3796 break;
3797 default:
3798 break;
3799 }
3800
3801 /* Loops around more primitive morphology methods
3802 ** erose, dilate, open, close, smooth, edge, etc...
3803 */
3804 /* Loop 1: iterate the compound method */
3805 method_loop = 0;
3806 method_changed = 1;
3807 while ( method_loop < method_limit && method_changed > 0 ) {
3808 method_loop++;
3809 method_changed = 0;
3810
3811 /* Loop 2: iterate over each kernel in a multi-kernel list */
3812 norm_kernel = (KernelInfo *) kernel;
3813 this_kernel = (KernelInfo *) kernel;
3814 rflt_kernel = reflected_kernel;
3815
3816 kernel_number = 0;
3817 while ( norm_kernel != NULL ) {
3818
3819 /* Loop 3: Compound Morphology Staging - Select Primitive to apply */
3820 stage_loop = 0; /* the compound morphology stage number */
3821 while ( stage_loop < stage_limit ) {
3822 stage_loop++; /* The stage of the compound morphology */
3823
3824 /* Select primitive morphology for this stage of compound method */
3825 this_kernel = norm_kernel; /* default use unreflected kernel */
3826 primitive = method; /* Assume method is a primitive */
3827 switch( method ) {
3828 case ErodeMorphology: /* just erode */
3829 case EdgeInMorphology: /* erode and image difference */
3830 primitive = ErodeMorphology;
3831 break;
3832 case DilateMorphology: /* just dilate */
3833 case EdgeOutMorphology: /* dilate and image difference */
3834 primitive = DilateMorphology;
3835 break;
3836 case OpenMorphology: /* erode then dilate */
3837 case TopHatMorphology: /* open and image difference */
3838 primitive = ErodeMorphology;
3839 if ( stage_loop == 2 )
3840 primitive = DilateMorphology;
3841 break;
3842 case OpenIntensityMorphology:
3843 primitive = ErodeIntensityMorphology;
3844 if ( stage_loop == 2 )
3845 primitive = DilateIntensityMorphology;
3846 break;
3847 case CloseMorphology: /* dilate, then erode */
3848 case BottomHatMorphology: /* close and image difference */
3849 this_kernel = rflt_kernel; /* use the reflected kernel */
3850 primitive = DilateMorphology;
3851 if ( stage_loop == 2 )
3852 primitive = ErodeMorphology;
3853 break;
3854 case CloseIntensityMorphology:
3855 this_kernel = rflt_kernel; /* use the reflected kernel */
3856 primitive = DilateIntensityMorphology;
3857 if ( stage_loop == 2 )
3858 primitive = ErodeIntensityMorphology;
3859 break;
3860 case SmoothMorphology: /* open, close */
3861 switch ( stage_loop ) {
3862 case 1: /* start an open method, which starts with Erode */
3863 primitive = ErodeMorphology;
3864 break;
3865 case 2: /* now Dilate the Erode */
3866 primitive = DilateMorphology;
3867 break;
3868 case 3: /* Reflect kernel a close */
3869 this_kernel = rflt_kernel; /* use the reflected kernel */
3870 primitive = DilateMorphology;
3871 break;
3872 case 4: /* Finish the Close */
3873 this_kernel = rflt_kernel; /* use the reflected kernel */
3874 primitive = ErodeMorphology;
3875 break;
3876 }
3877 break;
3878 case EdgeMorphology: /* dilate and erode difference */
3879 primitive = DilateMorphology;
3880 if ( stage_loop == 2 ) {
3881 save_image = curr_image; /* save the image difference */
3882 curr_image = (Image *) image;
3883 primitive = ErodeMorphology;
3884 }
3885 break;
3886 case CorrelateMorphology:
3887 /* A Correlation is a Convolution with a reflected kernel.
3888 ** However a Convolution is a weighted sum using a reflected
3889 ** kernel. It may seem strange to convert a Correlation into a
3890 ** Convolution as the Correlation is the simpler method, but
3891 ** Convolution is much more commonly used, and it makes sense to
3892 ** implement it directly so as to avoid the need to duplicate the
3893 ** kernel when it is not required (which is typically the
3894 ** default).
3895 */
3896 this_kernel = rflt_kernel; /* use the reflected kernel */
3897 primitive = ConvolveMorphology;
3898 break;
3899 default:
3900 break;
3901 }
3902 assert( this_kernel != (KernelInfo *) NULL );
3903
3904 /* Extra information for debugging compound operations */
3905 if (verbose != MagickFalse) {
3906 if ( stage_limit > 1 )
3907 (void) FormatLocaleString(v_info,MagickPathExtent,"%s:%.20g.%.20g -> ",
3908 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
3909 method_loop,(double) stage_loop);
3910 else if ( primitive != method )
3911 (void) FormatLocaleString(v_info, MagickPathExtent, "%s:%.20g -> ",
3912 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
3913 method_loop);
3914 else
3915 v_info[0] = '\0';
3916 }
3917
3918 /* Loop 4: Iterate the kernel with primitive */
3919 kernel_loop = 0;
3920 kernel_changed = 0;
3921 changed = 1;
3922 while ( kernel_loop < kernel_limit && changed > 0 ) {
3923 kernel_loop++; /* the iteration of this kernel */
3924
3925 /* Create a clone as the destination image, if not yet defined */
3926 if ( work_image == (Image *) NULL )
3927 {
3928 work_image=CloneImage(image,0,0,MagickTrue,exception);
3929 if (work_image == (Image *) NULL)
3930 goto error_cleanup;
3931 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
3932 goto error_cleanup;
3933 }
3934
3935 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3936 count++;
3937 changed = MorphologyPrimitive(curr_image, work_image, primitive,
3938 this_kernel, bias, exception);
3939 if (verbose != MagickFalse) {
3940 if ( kernel_loop > 1 )
3941 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
3942 (void) (void) FormatLocaleFile(stderr,
3943 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3944 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
3945 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
3946 (double) (method_loop+kernel_loop-1),(double) kernel_number,
3947 (double) count,(double) changed);
3948 }
3949 if ( changed < 0 )
3950 goto error_cleanup;
3951 kernel_changed = (size_t) ((ssize_t) kernel_changed+changed);
3952 method_changed = (size_t) ((ssize_t) method_changed+changed);
3953
3954 /* prepare next loop */
3955 { Image *tmp = work_image; /* swap images for iteration */
3956 work_image = curr_image;
3957 curr_image = tmp;
3958 }
3959 if ( work_image == image )
3960 work_image = (Image *) NULL; /* replace input 'image' */
3961
3962 } /* End Loop 4: Iterate the kernel with primitive */
3963
3964 if (verbose != MagickFalse && kernel_changed != (size_t)changed)
3965 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
3966 if (verbose != MagickFalse && stage_loop < stage_limit)
3967 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
3968
3969#if 0
3970 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3971 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
3972 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
3973 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
3974 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
3975#endif
3976
3977 } /* End Loop 3: Primitive (staging) Loop for Compound Methods */
3978
3979 /* Final Post-processing for some Compound Methods
3980 **
3981 ** The removal of any 'Sync' channel flag in the Image Composition
3982 ** below ensures the mathematical compose method is applied in a
3983 ** purely mathematical way, and only to the selected channels.
3984 ** Turn off SVG composition 'alpha blending'.
3985 */
3986 switch( method ) {
3987 case EdgeOutMorphology:
3988 case EdgeInMorphology:
3989 case TopHatMorphology:
3990 case BottomHatMorphology:
3991 if (verbose != MagickFalse)
3992 (void) FormatLocaleFile(stderr,
3993 "\n%s: Difference with original image",CommandOptionToMnemonic(
3994 MagickMorphologyOptions, method) );
3995 (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
3996 MagickTrue,0,0,exception);
3997 break;
3998 case EdgeMorphology:
3999 if (verbose != MagickFalse)
4000 (void) FormatLocaleFile(stderr,
4001 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
4002 MagickMorphologyOptions, method) );
4003 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
4004 MagickTrue,0,0,exception);
4005 save_image = DestroyImage(save_image); /* finished with save image */
4006 break;
4007 default:
4008 break;
4009 }
4010
4011 /* multi-kernel handling: re-iterate, or compose results */
4012 if ( kernel->next == (KernelInfo *) NULL )
4013 rslt_image = curr_image; /* just return the resulting image */
4014 else if ( rslt_compose == NoCompositeOp )
4015 { if (verbose != MagickFalse) {
4016 if ( this_kernel->next != (KernelInfo *) NULL )
4017 (void) FormatLocaleFile(stderr, " (re-iterate)");
4018 else
4019 (void) FormatLocaleFile(stderr, " (done)");
4020 }
4021 rslt_image = curr_image; /* return result, and re-iterate */
4022 }
4023 else if ( rslt_image == (Image *) NULL)
4024 { if (verbose != MagickFalse)
4025 (void) FormatLocaleFile(stderr, " (save for compose)");
4026 rslt_image = curr_image;
4027 curr_image = (Image *) image; /* continue with original image */
4028 }
4029 else
4030 { /* Add the new 'current' result to the composition
4031 **
4032 ** The removal of any 'Sync' channel flag in the Image Composition
4033 ** below ensures the mathematical compose method is applied in a
4034 ** purely mathematical way, and only to the selected channels.
4035 ** IE: Turn off SVG composition 'alpha blending'.
4036 */
4037 if (verbose != MagickFalse)
4038 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4039 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4040 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4041 0,0,exception);
4042 curr_image = DestroyImage(curr_image);
4043 curr_image = (Image *) image; /* continue with original image */
4044 }
4045 if (verbose != MagickFalse)
4046 (void) FormatLocaleFile(stderr, "\n");
4047
4048 /* loop to the next kernel in a multi-kernel list */
4049 norm_kernel = norm_kernel->next;
4050 if ( rflt_kernel != (KernelInfo *) NULL )
4051 rflt_kernel = rflt_kernel->next;
4052 kernel_number++;
4053 } /* End Loop 2: Loop over each kernel */
4054
4055 } /* End Loop 1: compound method interaction */
4056
4057 goto exit_cleanup;
4058
4059 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4060error_cleanup:
4061 if ( curr_image == rslt_image )
4062 curr_image = (Image *) NULL;
4063 if ( rslt_image != (Image *) NULL )
4064 rslt_image = DestroyImage(rslt_image);
4065exit_cleanup:
4066 if ( curr_image == rslt_image || curr_image == image )
4067 curr_image = (Image *) NULL;
4068 if ( curr_image != (Image *) NULL )
4069 curr_image = DestroyImage(curr_image);
4070 if ( work_image != (Image *) NULL )
4071 work_image = DestroyImage(work_image);
4072 if ( save_image != (Image *) NULL )
4073 save_image = DestroyImage(save_image);
4074 if ( reflected_kernel != (KernelInfo *) NULL )
4075 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4076 return(rslt_image);
4077}
4078
4079
4080/*
4081%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4082% %
4083% %
4084% %
4085% M o r p h o l o g y I m a g e %
4086% %
4087% %
4088% %
4089%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4090%
4091% MorphologyImage() applies a user supplied kernel to the image according to
4092% the given morphology method.
4093%
4094% This function applies any and all user defined settings before calling
4095% the above internal function MorphologyApply().
4096%
4097% User defined settings include...
4098% * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4099% * Kernel Scale/normalize settings ("-define convolve:scale=??")
4100% This can also includes the addition of a scaled unity kernel.
4101% * Show Kernel being applied ("-define morphology:showKernel=1")
4102%
4103% Other operators that do not want user supplied options interfering,
4104% especially "convolve:bias" and "morphology:showKernel" should use
4105% MorphologyApply() directly.
4106%
4107% The format of the MorphologyImage method is:
4108%
4109% Image *MorphologyImage(const Image *image,MorphologyMethod method,
4110% const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4111%
4112% A description of each parameter follows:
4113%
4114% o image: the image.
4115%
4116% o method: the morphology method to be applied.
4117%
4118% o iterations: apply the operation this many times (or no change).
4119% A value of -1 means loop until no change found.
4120% How this is applied may depend on the morphology method.
4121% Typically this is a value of 1.
4122%
4123% o kernel: An array of double representing the morphology kernel.
4124% Warning: kernel may be normalized for the Convolve method.
4125%
4126% o exception: return any errors or warnings in this structure.
4127%
4128*/
4129MagickExport Image *MorphologyImage(const Image *image,
4130 const MorphologyMethod method,const ssize_t iterations,
4131 const KernelInfo *kernel,ExceptionInfo *exception)
4132{
4133 const char
4134 *artifact;
4135
4136 CompositeOperator
4137 compose;
4138
4139 double
4140 bias;
4141
4142 Image
4143 *morphology_image;
4144
4146 *curr_kernel;
4147
4148 assert(image != (const Image *) NULL);
4149 assert(image->signature == MagickCoreSignature);
4150 assert(exception != (ExceptionInfo *) NULL);
4151 assert(exception->signature == MagickCoreSignature);
4152 if (IsEventLogging() != MagickFalse)
4153 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4154 curr_kernel = (KernelInfo *) kernel;
4155 bias=0.0;
4156 compose = UndefinedCompositeOp; /* use default for method */
4157
4158 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4159 * This is done BEFORE the ShowKernelInfo() function is called so that
4160 * users can see the results of the 'option:convolve:scale' option.
4161 */
4162 if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4163 /* Get the bias value as it will be needed */
4164 artifact = GetImageArtifact(image,"convolve:bias");
4165 if ( artifact != (const char *) NULL) {
4166 if (IsGeometry(artifact) == MagickFalse)
4167 (void) ThrowMagickException(exception,GetMagickModule(),
4168 OptionWarning,"InvalidSetting","'%s' '%s'",
4169 "convolve:bias",artifact);
4170 else
4171 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4172 }
4173
4174 /* Scale kernel according to user wishes */
4175 artifact = GetImageArtifact(image,"convolve:scale");
4176 if ( artifact != (const char *) NULL ) {
4177 if (IsGeometry(artifact) == MagickFalse)
4178 (void) ThrowMagickException(exception,GetMagickModule(),
4179 OptionWarning,"InvalidSetting","'%s' '%s'",
4180 "convolve:scale",artifact);
4181 else {
4182 if ( curr_kernel == kernel )
4183 curr_kernel = CloneKernelInfo(kernel);
4184 if (curr_kernel == (KernelInfo *) NULL)
4185 return((Image *) NULL);
4186 ScaleGeometryKernelInfo(curr_kernel, artifact);
4187 }
4188 }
4189 }
4190
4191 /* display the (normalized) kernel via stderr */
4192 artifact=GetImageArtifact(image,"morphology:showKernel");
4193 if (IsStringTrue(artifact) != MagickFalse)
4194 ShowKernelInfo(curr_kernel);
4195
4196 /* Override the default handling of multi-kernel morphology results
4197 * If 'Undefined' use the default method
4198 * If 'None' (default for 'Convolve') re-iterate previous result
4199 * Otherwise merge resulting images using compose method given.
4200 * Default for 'HitAndMiss' is 'Lighten'.
4201 */
4202 {
4203 ssize_t
4204 parse;
4205
4206 artifact = GetImageArtifact(image,"morphology:compose");
4207 if ( artifact != (const char *) NULL) {
4208 parse=ParseCommandOption(MagickComposeOptions,
4209 MagickFalse,artifact);
4210 if ( parse < 0 )
4211 (void) ThrowMagickException(exception,GetMagickModule(),
4212 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4213 "morphology:compose",artifact);
4214 else
4215 compose=(CompositeOperator)parse;
4216 }
4217 }
4218 /* Apply the Morphology */
4219 morphology_image = MorphologyApply(image,method,iterations,
4220 curr_kernel,compose,bias,exception);
4221
4222 /* Cleanup and Exit */
4223 if ( curr_kernel != kernel )
4224 curr_kernel=DestroyKernelInfo(curr_kernel);
4225 return(morphology_image);
4226}
4227
4228/*
4229%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4230% %
4231% %
4232% %
4233+ R o t a t e K e r n e l I n f o %
4234% %
4235% %
4236% %
4237%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4238%
4239% RotateKernelInfo() rotates the kernel by the angle given.
4240%
4241% Currently it is restricted to 90 degree angles, of either 1D kernels
4242% or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4243% It will ignore useless rotations for specific 'named' built-in kernels.
4244%
4245% The format of the RotateKernelInfo method is:
4246%
4247% void RotateKernelInfo(KernelInfo *kernel, double angle)
4248%
4249% A description of each parameter follows:
4250%
4251% o kernel: the Morphology/Convolution kernel
4252%
4253% o angle: angle to rotate in degrees
4254%
4255% This function is currently internal to this module only, but can be exported
4256% to other modules if needed.
4257*/
4258static void RotateKernelInfo(KernelInfo *kernel, double angle)
4259{
4260 /* angle the lower kernels first */
4261 if ( kernel->next != (KernelInfo *) NULL)
4262 RotateKernelInfo(kernel->next, angle);
4263
4264 /* WARNING: Currently assumes the kernel (rightly) is horizontally symmetrical
4265 **
4266 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4267 */
4268
4269 /* Modulus the angle */
4270 angle = fmod(angle, 360.0);
4271 if ( angle < 0 )
4272 angle += 360.0;
4273
4274 if ( 337.5 < angle || angle <= 22.5 )
4275 return; /* Near zero angle - no change! - At least not at this time */
4276
4277 /* Handle special cases */
4278 switch (kernel->type) {
4279 /* These built-in kernels are cylindrical kernels, rotating is useless */
4280 case GaussianKernel:
4281 case DoGKernel:
4282 case LoGKernel:
4283 case DiskKernel:
4284 case PeaksKernel:
4285 case LaplacianKernel:
4286 case ChebyshevKernel:
4287 case ManhattanKernel:
4288 case EuclideanKernel:
4289 return;
4290
4291 /* These may be rotatable at non-90 angles in the future */
4292 /* but simply rotating them in multiples of 90 degrees is useless */
4293 case SquareKernel:
4294 case DiamondKernel:
4295 case PlusKernel:
4296 case CrossKernel:
4297 return;
4298
4299 /* These only allows a +/-90 degree rotation (by transpose) */
4300 /* A 180 degree rotation is useless */
4301 case BlurKernel:
4302 if ( 135.0 < angle && angle <= 225.0 )
4303 return;
4304 if ( 225.0 < angle && angle <= 315.0 )
4305 angle -= 180;
4306 break;
4307
4308 default:
4309 break;
4310 }
4311 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4312 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4313 {
4314 if ( kernel->width == 3 && kernel->height == 3 )
4315 { /* Rotate a 3x3 square by 45 degree angle */
4316 double t = kernel->values[0];
4317 kernel->values[0] = kernel->values[3];
4318 kernel->values[3] = kernel->values[6];
4319 kernel->values[6] = kernel->values[7];
4320 kernel->values[7] = kernel->values[8];
4321 kernel->values[8] = kernel->values[5];
4322 kernel->values[5] = kernel->values[2];
4323 kernel->values[2] = kernel->values[1];
4324 kernel->values[1] = t;
4325 /* rotate non-centered origin */
4326 if ( kernel->x != 1 || kernel->y != 1 ) {
4327 ssize_t x,y;
4328 x = (ssize_t) kernel->x-1;
4329 y = (ssize_t) kernel->y-1;
4330 if ( x == y ) x = 0;
4331 else if ( x == 0 ) x = -y;
4332 else if ( x == -y ) y = 0;
4333 else if ( y == 0 ) y = x;
4334 kernel->x = (ssize_t) x+1;
4335 kernel->y = (ssize_t) y+1;
4336 }
4337 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4338 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4339 }
4340 else
4341 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4342 }
4343 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4344 {
4345 if ( kernel->width == 1 || kernel->height == 1 )
4346 { /* Do a transpose of a 1 dimensional kernel,
4347 ** which results in a fast 90 degree rotation of some type.
4348 */
4349 ssize_t
4350 t;
4351 t = (ssize_t) kernel->width;
4352 kernel->width = kernel->height;
4353 kernel->height = (size_t) t;
4354 t = kernel->x;
4355 kernel->x = kernel->y;
4356 kernel->y = t;
4357 if ( kernel->width == 1 ) {
4358 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4359 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4360 } else {
4361 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4362 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4363 }
4364 }
4365 else if ( kernel->width == kernel->height )
4366 { /* Rotate a square array of values by 90 degrees */
4367 { ssize_t
4368 i,j,x,y;
4369
4370 MagickRealType
4371 *k,t;
4372
4373 k=kernel->values;
4374 for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--)
4375 for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--)
4376 { t = k[i+j*(ssize_t) kernel->width];
4377 k[i+j*(ssize_t) kernel->width] = k[j+x*(ssize_t) kernel->width];
4378 k[j+x*(ssize_t) kernel->width] = k[x+y*(ssize_t) kernel->width];
4379 k[x+y*(ssize_t) kernel->width] = k[y+i*(ssize_t) kernel->width];
4380 k[y+i*(ssize_t) kernel->width] = t;
4381 }
4382 }
4383 /* rotate the origin - relative to center of array */
4384 { ssize_t x,y;
4385 x = (ssize_t) (kernel->x*2-(ssize_t) kernel->width+1);
4386 y = (ssize_t) (kernel->y*2-(ssize_t) kernel->height+1);
4387 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4388 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4389 }
4390 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4391 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4392 }
4393 else
4394 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4395 }
4396 if ( 135.0 < angle && angle <= 225.0 )
4397 {
4398 /* For a 180 degree rotation - also know as a reflection
4399 * This is actually a very very common operation!
4400 * Basically all that is needed is a reversal of the kernel data!
4401 * And a reflection of the origin
4402 */
4403 MagickRealType
4404 t;
4405
4406 MagickRealType
4407 *k;
4408
4409 ssize_t
4410 i,
4411 j;
4412
4413 k=kernel->values;
4414 j=(ssize_t) (kernel->width*kernel->height-1);
4415 for (i=0; i < j; i++, j--)
4416 t=k[i], k[i]=k[j], k[j]=t;
4417
4418 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4419 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4420 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4421 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4422 }
4423 /* At this point angle should at least between -45 (315) and +45 degrees
4424 * In the future some form of non-orthogonal angled rotates could be
4425 * performed here, possibly with a linear kernel restriction.
4426 */
4427
4428 return;
4429}
4430
4431/*
4432%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4433% %
4434% %
4435% %
4436% S c a l e G e o m e t r y K e r n e l I n f o %
4437% %
4438% %
4439% %
4440%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4441%
4442% ScaleGeometryKernelInfo() takes a geometry argument string, typically
4443% provided as a "-set option:convolve:scale {geometry}" user setting,
4444% and modifies the kernel according to the parsed arguments of that setting.
4445%
4446% The first argument (and any normalization flags) are passed to
4447% ScaleKernelInfo() to scale/normalize the kernel. The second argument
4448% is then passed to UnityAddKernelInfo() to add a scaled unity kernel
4449% into the scaled/normalized kernel.
4450%
4451% The format of the ScaleGeometryKernelInfo method is:
4452%
4453% void ScaleGeometryKernelInfo(KernelInfo *kernel,
4454% const double scaling_factor,const MagickStatusType normalize_flags)
4455%
4456% A description of each parameter follows:
4457%
4458% o kernel: the Morphology/Convolution kernel to modify
4459%
4460% o geometry:
4461% The geometry string to parse, typically from the user provided
4462% "-set option:convolve:scale {geometry}" setting.
4463%
4464*/
4465MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4466 const char *geometry)
4467{
4468 MagickStatusType
4469 flags;
4470
4472 args;
4473
4474 SetGeometryInfo(&args);
4475 flags = ParseGeometry(geometry, &args);
4476
4477#if 0
4478 /* For Debugging Geometry Input */
4479 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4480 flags, args.rho, args.sigma, args.xi, args.psi );
4481#endif
4482
4483 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4484 args.rho *= 0.01, args.sigma *= 0.01;
4485
4486 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4487 args.rho = 1.0;
4488 if ( (flags & SigmaValue) == 0 )
4489 args.sigma = 0.0;
4490
4491 /* Scale/Normalize the input kernel */
4492 ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4493
4494 /* Add Unity Kernel, for blending with original */
4495 if ( (flags & SigmaValue) != 0 )
4496 UnityAddKernelInfo(kernel, args.sigma);
4497
4498 return;
4499}
4500/*
4501%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4502% %
4503% %
4504% %
4505% S c a l e K e r n e l I n f o %
4506% %
4507% %
4508% %
4509%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4510%
4511% ScaleKernelInfo() scales the given kernel list by the given amount, with or
4512% without normalization of the sum of the kernel values (as per given flags).
4513%
4514% By default (no flags given) the values within the kernel is scaled
4515% directly using given scaling factor without change.
4516%
4517% If either of the two 'normalize_flags' are given the kernel will first be
4518% normalized and then further scaled by the scaling factor value given.
4519%
4520% Kernel normalization ('normalize_flags' given) is designed to ensure that
4521% any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4522% morphology methods will fall into -1.0 to +1.0 range. Note that for
4523% non-HDRI versions of IM this may cause images to have any negative results
4524% clipped, unless some 'bias' is used.
4525%
4526% More specifically. Kernels which only contain positive values (such as a
4527% 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4528% ensuring a 0.0 to +1.0 output range for non-HDRI images.
4529%
4530% For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4531% the kernel will be scaled by the absolute of the sum of kernel values, so
4532% that it will generally fall within the +/- 1.0 range.
4533%
4534% For kernels whose values sum to zero, (such as 'Laplacian' kernels) kernel
4535% will be scaled by just the sum of the positive values, so that its output
4536% range will again fall into the +/- 1.0 range.
4537%
4538% For special kernels designed for locating shapes using 'Correlate', (often
4539% only containing +1 and -1 values, representing foreground/background
4540% matching) a special normalization method is provided to scale the positive
4541% values separately to those of the negative values, so the kernel will be
4542% forced to become a zero-sum kernel better suited to such searches.
4543%
4544% WARNING: Correct normalization of the kernel assumes that the '*_range'
4545% attributes within the kernel structure have been correctly set during the
4546% kernels creation.
4547%
4548% NOTE: The values used for 'normalize_flags' have been selected specifically
4549% to match the use of geometry options, so that '!' means NormalizeValue, '^'
4550% means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4551%
4552% The format of the ScaleKernelInfo method is:
4553%
4554% void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4555% const MagickStatusType normalize_flags )
4556%
4557% A description of each parameter follows:
4558%
4559% o kernel: the Morphology/Convolution kernel
4560%
4561% o scaling_factor:
4562% multiply all values (after normalization) by this factor if not
4563% zero. If the kernel is normalized regardless of any flags.
4564%
4565% o normalize_flags:
4566% GeometryFlags defining normalization method to use.
4567% specifically: NormalizeValue, CorrelateNormalizeValue,
4568% and/or PercentValue
4569%
4570*/
4571MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4572 const double scaling_factor,const GeometryFlags normalize_flags)
4573{
4574 double
4575 pos_scale,
4576 neg_scale;
4577
4578 ssize_t
4579 i;
4580
4581 /* do the other kernels in a multi-kernel list first */
4582 if ( kernel->next != (KernelInfo *) NULL)
4583 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4584
4585 /* Normalization of Kernel */
4586 pos_scale = 1.0;
4587 if ( (normalize_flags&NormalizeValue) != 0 ) {
4588 if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4589 /* non-zero-summing kernel (generally positive) */
4590 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4591 else
4592 /* zero-summing kernel */
4593 pos_scale = kernel->positive_range;
4594 }
4595 /* Force kernel into a normalized zero-summing kernel */
4596 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4597 pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4598 ? kernel->positive_range : 1.0;
4599 neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4600 ? -kernel->negative_range : 1.0;
4601 }
4602 else
4603 neg_scale = pos_scale;
4604
4605 /* finalize scaling_factor for positive and negative components */
4606 pos_scale = scaling_factor/pos_scale;
4607 neg_scale = scaling_factor/neg_scale;
4608
4609 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4610 if (!IsNaN(kernel->values[i]))
4611 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4612
4613 /* convolution output range */
4614 kernel->positive_range *= pos_scale;
4615 kernel->negative_range *= neg_scale;
4616 /* maximum and minimum values in kernel */
4617 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4618 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4619
4620 /* swap kernel settings if user's scaling factor is negative */
4621 if ( scaling_factor < MagickEpsilon ) {
4622 double t;
4623 t = kernel->positive_range;
4624 kernel->positive_range = kernel->negative_range;
4625 kernel->negative_range = t;
4626 t = kernel->maximum;
4627 kernel->maximum = kernel->minimum;
4628 kernel->minimum = 1;
4629 }
4630
4631 return;
4632}
4633
4634/*
4635%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4636% %
4637% %
4638% %
4639% S h o w K e r n e l I n f o %
4640% %
4641% %
4642% %
4643%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4644%
4645% ShowKernelInfo() outputs the details of the given kernel definition to
4646% standard error, generally due to a users 'morphology:showKernel' option
4647% request.
4648%
4649% The format of the ShowKernel method is:
4650%
4651% void ShowKernelInfo(const KernelInfo *kernel)
4652%
4653% A description of each parameter follows:
4654%
4655% o kernel: the Morphology/Convolution kernel
4656%
4657*/
4658MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4659{
4660 const KernelInfo
4661 *k;
4662
4663 size_t
4664 c, i, u, v;
4665
4666 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4667
4668 (void) FormatLocaleFile(stderr, "Kernel");
4669 if ( kernel->next != (KernelInfo *) NULL )
4670 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4671 (void) FormatLocaleFile(stderr, " \"%s",
4672 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4673 if ( fabs(k->angle) >= MagickEpsilon )
4674 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4675 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4676 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4677 (void) FormatLocaleFile(stderr,
4678 " with values from %.*lg to %.*lg\n",
4679 GetMagickPrecision(), k->minimum,
4680 GetMagickPrecision(), k->maximum);
4681 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4682 GetMagickPrecision(), k->negative_range,
4683 GetMagickPrecision(), k->positive_range);
4684 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4685 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4686 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4687 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4688 else
4689 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4690 GetMagickPrecision(), k->positive_range+k->negative_range);
4691 for (i=v=0; v < k->height; v++) {
4692 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4693 for (u=0; u < k->width; u++, i++)
4694 if (IsNaN(k->values[i]))
4695 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4696 else
4697 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4698 GetMagickPrecision(), (double) k->values[i]);
4699 (void) FormatLocaleFile(stderr,"\n");
4700 }
4701 }
4702}
4703
4704/*
4705%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4706% %
4707% %
4708% %
4709% U n i t y A d d K e r n a l I n f o %
4710% %
4711% %
4712% %
4713%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4714%
4715% UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4716% to the given pre-scaled and normalized Kernel. This in effect adds that
4717% amount of the original image into the resulting convolution kernel. This
4718% value is usually provided by the user as a percentage value in the
4719% 'convolve:scale' setting.
4720%
4721% The resulting effect is to convert the defined kernels into blended
4722% soft-blurs, unsharp kernels or into sharpening kernels.
4723%
4724% The format of the UnityAdditionKernelInfo method is:
4725%
4726% void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4727%
4728% A description of each parameter follows:
4729%
4730% o kernel: the Morphology/Convolution kernel
4731%
4732% o scale:
4733% scaling factor for the unity kernel to be added to
4734% the given kernel.
4735%
4736*/
4737MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4738 const double scale)
4739{
4740 /* do the other kernels in a multi-kernel list first */
4741 if ( kernel->next != (KernelInfo *) NULL)
4742 UnityAddKernelInfo(kernel->next, scale);
4743
4744 /* Add the scaled unity kernel to the existing kernel */
4745 kernel->values[kernel->x+kernel->y*(ssize_t) kernel->width] += scale;
4746 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4747
4748 return;
4749}
4750
4751/*
4752%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4753% %
4754% %
4755% %
4756% Z e r o K e r n e l N a n s %
4757% %
4758% %
4759% %
4760%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4761%
4762% ZeroKernelNans() replaces any special 'nan' value that may be present in
4763% the kernel with a zero value. This is typically done when the kernel will
4764% be used in special hardware (GPU) convolution processors, to simply
4765% matters.
4766%
4767% The format of the ZeroKernelNans method is:
4768%
4769% void ZeroKernelNans (KernelInfo *kernel)
4770%
4771% A description of each parameter follows:
4772%
4773% o kernel: the Morphology/Convolution kernel
4774%
4775*/
4776MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4777{
4778 size_t
4779 i;
4780
4781 /* do the other kernels in a multi-kernel list first */
4782 if (kernel->next != (KernelInfo *) NULL)
4783 ZeroKernelNans(kernel->next);
4784
4785 for (i=0; i < (kernel->width*kernel->height); i++)
4786 if (IsNaN(kernel->values[i]))
4787 kernel->values[i]=0.0;
4788
4789 return;
4790}