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 1999-2021 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 bluring 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/pixel-private.h"
75 #include "MagickCore/prepress.h"
76 #include "MagickCore/quantize.h"
77 #include "MagickCore/resource_.h"
78 #include "MagickCore/registry.h"
79 #include "MagickCore/semaphore.h"
80 #include "MagickCore/splay-tree.h"
81 #include "MagickCore/statistic.h"
82 #include "MagickCore/string_.h"
83 #include "MagickCore/string-private.h"
84 #include "MagickCore/thread-private.h"
85 #include "MagickCore/token.h"
86 #include "MagickCore/utility.h"
87 #include "MagickCore/utility-private.h"
88
89 /*
90 Other global definitions used by module.
91 */
92 #define Minimize(assign,value) assign=MagickMin(assign,value)
93 #define Maximize(assign,value) assign=MagickMax(assign,value)
94
95 /* Integer Factorial Function - for a Binomial kernel */
96 #if 1
fact(size_t n)97 static inline size_t fact(size_t n)
98 {
99 size_t f,l;
100 for(f=1, l=2; l <= n; f=f*l, l++);
101 return(f);
102 }
103 #elif 1 /* glibc floating point alternatives */
104 #define fact(n) ((size_t)tgamma((double)n+1))
105 #else
106 #define fact(n) ((size_t)lgamma((double)n+1))
107 #endif
108
109
110 /* Currently these are only internal to this module */
111 static void
112 CalcKernelMetaData(KernelInfo *),
113 ExpandMirrorKernelInfo(KernelInfo *),
114 ExpandRotateKernelInfo(KernelInfo *, const double),
115 RotateKernelInfo(KernelInfo *, double);
116
117
118 /* Quick function to find last kernel in a kernel list */
LastKernelInfo(KernelInfo * kernel)119 static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
120 {
121 while (kernel->next != (KernelInfo *) NULL)
122 kernel=kernel->next;
123 return(kernel);
124 }
125
126 /*
127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
128 % %
129 % %
130 % %
131 % A c q u i r e K e r n e l I n f o %
132 % %
133 % %
134 % %
135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
136 %
137 % AcquireKernelInfo() takes the given string (generally supplied by the
138 % user) and converts it into a Morphology/Convolution Kernel. This allows
139 % users to specify a kernel from a number of pre-defined kernels, or to fully
140 % specify their own kernel for a specific Convolution or Morphology
141 % Operation.
142 %
143 % The kernel so generated can be any rectangular array of floating point
144 % values (doubles) with the 'control point' or 'pixel being affected'
145 % anywhere within that array of values.
146 %
147 % Previously IM was restricted to a square of odd size using the exact
148 % center as origin, this is no longer the case, and any rectangular kernel
149 % with any value being declared the origin. This in turn allows the use of
150 % highly asymmetrical kernels.
151 %
152 % The floating point values in the kernel can also include a special value
153 % known as 'nan' or 'not a number' to indicate that this value is not part
154 % of the kernel array. This allows you to shaped the kernel within its
155 % rectangular area. That is 'nan' values provide a 'mask' for the kernel
156 % shape. However at least one non-nan value must be provided for correct
157 % working of a kernel.
158 %
159 % The returned kernel should be freed using the DestroyKernelInfo() when you
160 % are finished with it. Do not free this memory yourself.
161 %
162 % Input kernel defintion strings can consist of any of three types.
163 %
164 % "name:args[[@><]"
165 % Select from one of the built in kernels, using the name and
166 % geometry arguments supplied. See AcquireKernelBuiltIn()
167 %
168 % "WxH[+X+Y][@><]:num, num, num ..."
169 % a kernel of size W by H, with W*H floating point numbers following.
170 % the 'center' can be optionally be defined at +X+Y (such that +0+0
171 % is top left corner). If not defined the pixel in the center, for
172 % odd sizes, or to the immediate top or left of center for even sizes
173 % is automatically selected.
174 %
175 % "num, num, num, num, ..."
176 % list of floating point numbers defining an 'old style' odd sized
177 % square kernel. At least 9 values should be provided for a 3x3
178 % square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
179 % Values can be space or comma separated. This is not recommended.
180 %
181 % You can define a 'list of kernels' which can be used by some morphology
182 % operators A list is defined as a semi-colon separated list kernels.
183 %
184 % " kernel ; kernel ; kernel ; "
185 %
186 % Any extra ';' characters, at start, end or between kernel defintions are
187 % simply ignored.
188 %
189 % The special flags will expand a single kernel, into a list of rotated
190 % kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
191 % cyclic rotations, while a '>' will generate a list of 90-degree rotations.
192 % The '<' also exands using 90-degree rotates, but giving a 180-degree
193 % reflected kernel before the +/- 90-degree rotations, which can be important
194 % for Thinning operations.
195 %
196 % Note that 'name' kernels will start with an alphabetic character while the
197 % new kernel specification has a ':' character in its specification string.
198 % If neither is the case, it is assumed an old style of a simple list of
199 % numbers generating a odd-sized square kernel has been given.
200 %
201 % The format of the AcquireKernal method is:
202 %
203 % KernelInfo *AcquireKernelInfo(const char *kernel_string)
204 %
205 % A description of each parameter follows:
206 %
207 % o kernel_string: the Morphology/Convolution kernel wanted.
208 %
209 */
210
211 /* This was separated so that it could be used as a separate
212 ** array input handling function, such as for -color-matrix
213 */
ParseKernelArray(const char * kernel_string)214 static KernelInfo *ParseKernelArray(const char *kernel_string)
215 {
216 KernelInfo
217 *kernel;
218
219 char
220 token[MagickPathExtent];
221
222 const char
223 *p,
224 *end;
225
226 ssize_t
227 i;
228
229 double
230 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
231
232 MagickStatusType
233 flags;
234
235 GeometryInfo
236 args;
237
238 kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
239 if (kernel == (KernelInfo *) NULL)
240 return(kernel);
241 (void) memset(kernel,0,sizeof(*kernel));
242 kernel->minimum = kernel->maximum = kernel->angle = 0.0;
243 kernel->negative_range = kernel->positive_range = 0.0;
244 kernel->type = UserDefinedKernel;
245 kernel->next = (KernelInfo *) NULL;
246 kernel->signature=MagickCoreSignature;
247 if (kernel_string == (const char *) NULL)
248 return(kernel);
249
250 /* find end of this specific kernel definition string */
251 end = strchr(kernel_string, ';');
252 if ( end == (char *) NULL )
253 end = strchr(kernel_string, '\0');
254
255 /* clear flags - for Expanding kernel lists thorugh rotations */
256 flags = NoValue;
257
258 /* Has a ':' in argument - New user kernel specification
259 FUTURE: this split on ':' could be done by StringToken()
260 */
261 p = strchr(kernel_string, ':');
262 if ( p != (char *) NULL && p < end)
263 {
264 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
265 memcpy(token, kernel_string, (size_t) (p-kernel_string));
266 token[p-kernel_string] = '\0';
267 SetGeometryInfo(&args);
268 flags = ParseGeometry(token, &args);
269
270 /* Size handling and checks of geometry settings */
271 if ( (flags & WidthValue) == 0 ) /* if no width then */
272 args.rho = args.sigma; /* then width = height */
273 if ( args.rho < 1.0 ) /* if width too small */
274 args.rho = 1.0; /* then width = 1 */
275 if ( args.sigma < 1.0 ) /* if height too small */
276 args.sigma = args.rho; /* then height = width */
277 kernel->width = (size_t)args.rho;
278 kernel->height = (size_t)args.sigma;
279
280 /* Offset Handling and Checks */
281 if ( args.xi < 0.0 || args.psi < 0.0 )
282 return(DestroyKernelInfo(kernel));
283 kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
284 : (ssize_t) (kernel->width-1)/2;
285 kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
286 : (ssize_t) (kernel->height-1)/2;
287 if ( kernel->x >= (ssize_t) kernel->width ||
288 kernel->y >= (ssize_t) kernel->height )
289 return(DestroyKernelInfo(kernel));
290
291 p++; /* advance beyond the ':' */
292 }
293 else
294 { /* ELSE - Old old specification, forming odd-square kernel */
295 /* count up number of values given */
296 p=(const char *) kernel_string;
297 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
298 p++; /* ignore "'" chars for convolve filter usage - Cristy */
299 for (i=0; p < end; i++)
300 {
301 (void) GetNextToken(p,&p,MagickPathExtent,token);
302 if (*token == ',')
303 (void) GetNextToken(p,&p,MagickPathExtent,token);
304 }
305 /* set the size of the kernel - old sized square */
306 kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
307 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
308 p=(const char *) kernel_string;
309 while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
310 p++; /* ignore "'" chars for convolve filter usage - Cristy */
311 }
312
313 /* Read in the kernel values from rest of input string argument */
314 kernel->values=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
315 kernel->width,kernel->height*sizeof(*kernel->values)));
316 if (kernel->values == (MagickRealType *) NULL)
317 return(DestroyKernelInfo(kernel));
318 kernel->minimum=MagickMaximumValue;
319 kernel->maximum=(-MagickMaximumValue);
320 kernel->negative_range = kernel->positive_range = 0.0;
321 for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
322 {
323 (void) GetNextToken(p,&p,MagickPathExtent,token);
324 if (*token == ',')
325 (void) GetNextToken(p,&p,MagickPathExtent,token);
326 if ( LocaleCompare("nan",token) == 0
327 || LocaleCompare("-",token) == 0 ) {
328 kernel->values[i] = nan; /* this value is not part of neighbourhood */
329 }
330 else {
331 kernel->values[i] = StringToDouble(token,(char **) NULL);
332 ( kernel->values[i] < 0)
333 ? ( kernel->negative_range += kernel->values[i] )
334 : ( kernel->positive_range += kernel->values[i] );
335 Minimize(kernel->minimum, kernel->values[i]);
336 Maximize(kernel->maximum, kernel->values[i]);
337 }
338 }
339
340 /* sanity check -- no more values in kernel definition */
341 (void) GetNextToken(p,&p,MagickPathExtent,token);
342 if ( *token != '\0' && *token != ';' && *token != '\'' )
343 return(DestroyKernelInfo(kernel));
344
345 #if 0
346 /* this was the old method of handling a incomplete kernel */
347 if ( i < (ssize_t) (kernel->width*kernel->height) ) {
348 Minimize(kernel->minimum, kernel->values[i]);
349 Maximize(kernel->maximum, kernel->values[i]);
350 for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
351 kernel->values[i]=0.0;
352 }
353 #else
354 /* Number of values for kernel was not enough - Report Error */
355 if ( i < (ssize_t) (kernel->width*kernel->height) )
356 return(DestroyKernelInfo(kernel));
357 #endif
358
359 /* check that we recieved at least one real (non-nan) value! */
360 if (kernel->minimum == MagickMaximumValue)
361 return(DestroyKernelInfo(kernel));
362
363 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel size */
364 ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
365 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
366 ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
367 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
368 ExpandMirrorKernelInfo(kernel); /* 90 degree mirror rotate */
369
370 return(kernel);
371 }
372
ParseKernelName(const char * kernel_string,ExceptionInfo * exception)373 static KernelInfo *ParseKernelName(const char *kernel_string,
374 ExceptionInfo *exception)
375 {
376 char
377 token[MagickPathExtent];
378
379 const char
380 *p,
381 *end;
382
383 GeometryInfo
384 args;
385
386 KernelInfo
387 *kernel;
388
389 MagickStatusType
390 flags;
391
392 ssize_t
393 type;
394
395 /* Parse special 'named' kernel */
396 (void) GetNextToken(kernel_string,&p,MagickPathExtent,token);
397 type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
398 if ( type < 0 || type == UserDefinedKernel )
399 return((KernelInfo *) NULL); /* not a valid named kernel */
400
401 while (((isspace((int) ((unsigned char) *p)) != 0) ||
402 (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
403 p++;
404
405 end = strchr(p, ';'); /* end of this kernel defintion */
406 if ( end == (char *) NULL )
407 end = strchr(p, '\0');
408
409 /* ParseGeometry() needs the geometry separated! -- Arrgghh */
410 memcpy(token, p, (size_t) (end-p));
411 token[end-p] = '\0';
412 SetGeometryInfo(&args);
413 flags = ParseGeometry(token, &args);
414
415 #if 0
416 /* For Debugging Geometry Input */
417 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
418 flags, args.rho, args.sigma, args.xi, args.psi );
419 #endif
420
421 /* special handling of missing values in input string */
422 switch( type ) {
423 /* Shape Kernel Defaults */
424 case UnityKernel:
425 if ( (flags & WidthValue) == 0 )
426 args.rho = 1.0; /* Default scale = 1.0, zero is valid */
427 break;
428 case SquareKernel:
429 case DiamondKernel:
430 case OctagonKernel:
431 case DiskKernel:
432 case PlusKernel:
433 case CrossKernel:
434 if ( (flags & HeightValue) == 0 )
435 args.sigma = 1.0; /* Default scale = 1.0, zero is valid */
436 break;
437 case RingKernel:
438 if ( (flags & XValue) == 0 )
439 args.xi = 1.0; /* Default scale = 1.0, zero is valid */
440 break;
441 case RectangleKernel: /* Rectangle - set size defaults */
442 if ( (flags & WidthValue) == 0 ) /* if no width then */
443 args.rho = args.sigma; /* then width = height */
444 if ( args.rho < 1.0 ) /* if width too small */
445 args.rho = 3; /* then width = 3 */
446 if ( args.sigma < 1.0 ) /* if height too small */
447 args.sigma = args.rho; /* then height = width */
448 if ( (flags & XValue) == 0 ) /* center offset if not defined */
449 args.xi = (double)(((ssize_t)args.rho-1)/2);
450 if ( (flags & YValue) == 0 )
451 args.psi = (double)(((ssize_t)args.sigma-1)/2);
452 break;
453 /* Distance Kernel Defaults */
454 case ChebyshevKernel:
455 case ManhattanKernel:
456 case OctagonalKernel:
457 case EuclideanKernel:
458 if ( (flags & HeightValue) == 0 ) /* no distance scale */
459 args.sigma = 100.0; /* default distance scaling */
460 else if ( (flags & AspectValue ) != 0 ) /* '!' flag */
461 args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
462 else if ( (flags & PercentValue ) != 0 ) /* '%' flag */
463 args.sigma *= QuantumRange/100.0; /* percentage of color range */
464 break;
465 default:
466 break;
467 }
468
469 kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args, exception);
470 if ( kernel == (KernelInfo *) NULL )
471 return(kernel);
472
473 /* global expand to rotated kernel list - only for single kernels */
474 if ( kernel->next == (KernelInfo *) NULL ) {
475 if ( (flags & AreaValue) != 0 ) /* '@' symbol in kernel args */
476 ExpandRotateKernelInfo(kernel, 45.0);
477 else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
478 ExpandRotateKernelInfo(kernel, 90.0);
479 else if ( (flags & LessValue) != 0 ) /* '<' symbol in kernel args */
480 ExpandMirrorKernelInfo(kernel);
481 }
482
483 return(kernel);
484 }
485
AcquireKernelInfo(const char * kernel_string,ExceptionInfo * exception)486 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string,
487 ExceptionInfo *exception)
488 {
489 KernelInfo
490 *kernel,
491 *new_kernel;
492
493 char
494 *kernel_cache,
495 token[MagickPathExtent];
496
497 const char
498 *p;
499
500 if (kernel_string == (const char *) NULL)
501 return(ParseKernelArray(kernel_string));
502 p=kernel_string;
503 kernel_cache=(char *) NULL;
504 if (*kernel_string == '@')
505 {
506 kernel_cache=FileToString(kernel_string+1,~0UL,exception);
507 if (kernel_cache == (char *) NULL)
508 return((KernelInfo *) NULL);
509 p=(const char *) kernel_cache;
510 }
511 kernel=NULL;
512 while (GetNextToken(p,(const char **) NULL,MagickPathExtent,token), *token != '\0')
513 {
514 /* ignore extra or multiple ';' kernel separators */
515 if (*token != ';')
516 {
517 /* tokens starting with alpha is a Named kernel */
518 if (isalpha((int) ((unsigned char) *token)) != 0)
519 new_kernel=ParseKernelName(p,exception);
520 else /* otherwise a user defined kernel array */
521 new_kernel=ParseKernelArray(p);
522
523 /* Error handling -- this is not proper error handling! */
524 if (new_kernel == (KernelInfo *) NULL)
525 {
526 if (kernel != (KernelInfo *) NULL)
527 kernel=DestroyKernelInfo(kernel);
528 return((KernelInfo *) NULL);
529 }
530
531 /* initialise or append the kernel list */
532 if (kernel == (KernelInfo *) NULL)
533 kernel=new_kernel;
534 else
535 LastKernelInfo(kernel)->next=new_kernel;
536 }
537
538 /* look for the next kernel in list */
539 p=strchr(p,';');
540 if (p == (char *) NULL)
541 break;
542 p++;
543 }
544 if (kernel_cache != (char *) NULL)
545 kernel_cache=DestroyString(kernel_cache);
546 return(kernel);
547 }
548
549 /*
550 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
551 % %
552 % %
553 % %
554 % A c q u i r e K e r n e l B u i l t I n %
555 % %
556 % %
557 % %
558 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
559 %
560 % AcquireKernelBuiltIn() returned one of the 'named' built-in types of
561 % kernels used for special purposes such as gaussian blurring, skeleton
562 % pruning, and edge distance determination.
563 %
564 % They take a KernelType, and a set of geometry style arguments, which were
565 % typically decoded from a user supplied string, or from a more complex
566 % Morphology Method that was requested.
567 %
568 % The format of the AcquireKernalBuiltIn method is:
569 %
570 % KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
571 % const GeometryInfo args)
572 %
573 % A description of each parameter follows:
574 %
575 % o type: the pre-defined type of kernel wanted
576 %
577 % o args: arguments defining or modifying the kernel
578 %
579 % Convolution Kernels
580 %
581 % Unity
582 % The a No-Op or Scaling single element kernel.
583 %
584 % Gaussian:{radius},{sigma}
585 % Generate a two-dimensional gaussian kernel, as used by -gaussian.
586 % The sigma for the curve is required. The resulting kernel is
587 % normalized,
588 %
589 % If 'sigma' is zero, you get a single pixel on a field of zeros.
590 %
591 % NOTE: that the 'radius' is optional, but if provided can limit (clip)
592 % the final size of the resulting kernel to a square 2*radius+1 in size.
593 % The radius should be at least 2 times that of the sigma value, or
594 % sever clipping and aliasing may result. If not given or set to 0 the
595 % radius will be determined so as to produce the best minimal error
596 % result, which is usally much larger than is normally needed.
597 %
598 % LoG:{radius},{sigma}
599 % "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
600 % The supposed ideal edge detection, zero-summing kernel.
601 %
602 % An alturnative to this kernel is to use a "DoG" with a sigma ratio of
603 % approx 1.6 (according to wikipedia).
604 %
605 % DoG:{radius},{sigma1},{sigma2}
606 % "Difference of Gaussians" Kernel.
607 % As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
608 % from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
609 % The result is a zero-summing kernel.
610 %
611 % Blur:{radius},{sigma}[,{angle}]
612 % Generates a 1 dimensional or linear gaussian blur, at the angle given
613 % (current restricted to orthogonal angles). If a 'radius' is given the
614 % kernel is clipped to a width of 2*radius+1. Kernel can be rotated
615 % by a 90 degree angle.
616 %
617 % If 'sigma' is zero, you get a single pixel on a field of zeros.
618 %
619 % Note that two convolutions with two "Blur" kernels perpendicular to
620 % each other, is equivalent to a far larger "Gaussian" kernel with the
621 % same sigma value, However it is much faster to apply. This is how the
622 % "-blur" operator actually works.
623 %
624 % Comet:{width},{sigma},{angle}
625 % Blur in one direction only, much like how a bright object leaves
626 % a comet like trail. The Kernel is actually half a gaussian curve,
627 % Adding two such blurs in opposite directions produces a Blur Kernel.
628 % Angle can be rotated in multiples of 90 degrees.
629 %
630 % Note that the first argument is the width of the kernel and not the
631 % radius of the kernel.
632 %
633 % Binomial:[{radius}]
634 % Generate a discrete kernel using a 2 dimentional Pascel's Triangle
635 % of values. Used for special forma of image filters.
636 %
637 % # Still to be implemented...
638 % #
639 % # Filter2D
640 % # Filter1D
641 % # Set kernel values using a resize filter, and given scale (sigma)
642 % # Cylindrical or Linear. Is this possible with an image?
643 % #
644 %
645 % Named Constant Convolution Kernels
646 %
647 % All these are unscaled, zero-summing kernels by default. As such for
648 % non-HDRI version of ImageMagick some form of normalization, user scaling,
649 % and biasing the results is recommended, to prevent the resulting image
650 % being 'clipped'.
651 %
652 % The 3x3 kernels (most of these) can be circularly rotated in multiples of
653 % 45 degrees to generate the 8 angled varients of each of the kernels.
654 %
655 % Laplacian:{type}
656 % Discrete Lapacian Kernels, (without normalization)
657 % Type 0 : 3x3 with center:8 surounded by -1 (8 neighbourhood)
658 % Type 1 : 3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
659 % Type 2 : 3x3 with center:4 edge:1 corner:-2
660 % Type 3 : 3x3 with center:4 edge:-2 corner:1
661 % Type 5 : 5x5 laplacian
662 % Type 7 : 7x7 laplacian
663 % Type 15 : 5x5 LoG (sigma approx 1.4)
664 % Type 19 : 9x9 LoG (sigma approx 1.4)
665 %
666 % Sobel:{angle}
667 % Sobel 'Edge' convolution kernel (3x3)
668 % | -1, 0, 1 |
669 % | -2, 0,-2 |
670 % | -1, 0, 1 |
671 %
672 % Roberts:{angle}
673 % Roberts convolution kernel (3x3)
674 % | 0, 0, 0 |
675 % | -1, 1, 0 |
676 % | 0, 0, 0 |
677 %
678 % Prewitt:{angle}
679 % Prewitt Edge convolution kernel (3x3)
680 % | -1, 0, 1 |
681 % | -1, 0, 1 |
682 % | -1, 0, 1 |
683 %
684 % Compass:{angle}
685 % Prewitt's "Compass" convolution kernel (3x3)
686 % | -1, 1, 1 |
687 % | -1,-2, 1 |
688 % | -1, 1, 1 |
689 %
690 % Kirsch:{angle}
691 % Kirsch's "Compass" convolution kernel (3x3)
692 % | -3,-3, 5 |
693 % | -3, 0, 5 |
694 % | -3,-3, 5 |
695 %
696 % FreiChen:{angle}
697 % Frei-Chen Edge Detector is based on a kernel that is similar to
698 % the Sobel Kernel, but is designed to be isotropic. That is it takes
699 % into account the distance of the diagonal in the kernel.
700 %
701 % | 1, 0, -1 |
702 % | sqrt(2), 0, -sqrt(2) |
703 % | 1, 0, -1 |
704 %
705 % FreiChen:{type},{angle}
706 %
707 % Frei-Chen Pre-weighted kernels...
708 %
709 % Type 0: default un-nomalized version shown above.
710 %
711 % Type 1: Orthogonal Kernel (same as type 11 below)
712 % | 1, 0, -1 |
713 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
714 % | 1, 0, -1 |
715 %
716 % Type 2: Diagonal form of Kernel...
717 % | 1, sqrt(2), 0 |
718 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
719 % | 0, -sqrt(2) -1 |
720 %
721 % However this kernel is als at the heart of the FreiChen Edge Detection
722 % Process which uses a set of 9 specially weighted kernel. These 9
723 % kernels not be normalized, but directly applied to the image. The
724 % results is then added together, to produce the intensity of an edge in
725 % a specific direction. The square root of the pixel value can then be
726 % taken as the cosine of the edge, and at least 2 such runs at 90 degrees
727 % from each other, both the direction and the strength of the edge can be
728 % determined.
729 %
730 % Type 10: All 9 of the following pre-weighted kernels...
731 %
732 % Type 11: | 1, 0, -1 |
733 % | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
734 % | 1, 0, -1 |
735 %
736 % Type 12: | 1, sqrt(2), 1 |
737 % | 0, 0, 0 | / 2*sqrt(2)
738 % | 1, sqrt(2), 1 |
739 %
740 % Type 13: | sqrt(2), -1, 0 |
741 % | -1, 0, 1 | / 2*sqrt(2)
742 % | 0, 1, -sqrt(2) |
743 %
744 % Type 14: | 0, 1, -sqrt(2) |
745 % | -1, 0, 1 | / 2*sqrt(2)
746 % | sqrt(2), -1, 0 |
747 %
748 % Type 15: | 0, -1, 0 |
749 % | 1, 0, 1 | / 2
750 % | 0, -1, 0 |
751 %
752 % Type 16: | 1, 0, -1 |
753 % | 0, 0, 0 | / 2
754 % | -1, 0, 1 |
755 %
756 % Type 17: | 1, -2, 1 |
757 % | -2, 4, -2 | / 6
758 % | -1, -2, 1 |
759 %
760 % Type 18: | -2, 1, -2 |
761 % | 1, 4, 1 | / 6
762 % | -2, 1, -2 |
763 %
764 % Type 19: | 1, 1, 1 |
765 % | 1, 1, 1 | / 3
766 % | 1, 1, 1 |
767 %
768 % The first 4 are for edge detection, the next 4 are for line detection
769 % and the last is to add a average component to the results.
770 %
771 % Using a special type of '-1' will return all 9 pre-weighted kernels
772 % as a multi-kernel list, so that you can use them directly (without
773 % normalization) with the special "-set option:morphology:compose Plus"
774 % setting to apply the full FreiChen Edge Detection Technique.
775 %
776 % If 'type' is large it will be taken to be an actual rotation angle for
777 % the default FreiChen (type 0) kernel. As such FreiChen:45 will look
778 % like a Sobel:45 but with 'sqrt(2)' instead of '2' values.
779 %
780 % WARNING: The above was layed out as per
781 % http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
782 % But rotated 90 degrees so direction is from left rather than the top.
783 % I have yet to find any secondary confirmation of the above. The only
784 % other source found was actual source code at
785 % http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
786 % Neigher paper defineds the kernels in a way that looks locical or
787 % correct when taken as a whole.
788 %
789 % Boolean Kernels
790 %
791 % Diamond:[{radius}[,{scale}]]
792 % Generate a diamond shaped kernel with given radius to the points.
793 % Kernel size will again be radius*2+1 square and defaults to radius 1,
794 % generating a 3x3 kernel that is slightly larger than a square.
795 %
796 % Square:[{radius}[,{scale}]]
797 % Generate a square shaped kernel of size radius*2+1, and defaulting
798 % to a 3x3 (radius 1).
799 %
800 % Octagon:[{radius}[,{scale}]]
801 % Generate octagonal shaped kernel of given radius and constant scale.
802 % Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
803 % in "Diamond" kernel.
804 %
805 % Disk:[{radius}[,{scale}]]
806 % Generate a binary disk, thresholded at the radius given, the radius
807 % may be a float-point value. Final Kernel size is floor(radius)*2+1
808 % square. A radius of 5.3 is the default.
809 %
810 % NOTE: That a low radii Disk kernels produce the same results as
811 % many of the previously defined kernels, but differ greatly at larger
812 % radii. Here is a table of equivalences...
813 % "Disk:1" => "Diamond", "Octagon:1", or "Cross:1"
814 % "Disk:1.5" => "Square"
815 % "Disk:2" => "Diamond:2"
816 % "Disk:2.5" => "Octagon"
817 % "Disk:2.9" => "Square:2"
818 % "Disk:3.5" => "Octagon:3"
819 % "Disk:4.5" => "Octagon:4"
820 % "Disk:5.4" => "Octagon:5"
821 % "Disk:6.4" => "Octagon:6"
822 % All other Disk shapes are unique to this kernel, but because a "Disk"
823 % is more circular when using a larger radius, using a larger radius is
824 % preferred over iterating the morphological operation.
825 %
826 % Rectangle:{geometry}
827 % Simply generate a rectangle of 1's with the size given. You can also
828 % specify the location of the 'control point', otherwise the closest
829 % pixel to the center of the rectangle is selected.
830 %
831 % Properly centered and odd sized rectangles work the best.
832 %
833 % Symbol Dilation Kernels
834 %
835 % These kernel is not a good general morphological kernel, but is used
836 % more for highlighting and marking any single pixels in an image using,
837 % a "Dilate" method as appropriate.
838 %
839 % For the same reasons iterating these kernels does not produce the
840 % same result as using a larger radius for the symbol.
841 %
842 % Plus:[{radius}[,{scale}]]
843 % Cross:[{radius}[,{scale}]]
844 % Generate a kernel in the shape of a 'plus' or a 'cross' with
845 % a each arm the length of the given radius (default 2).
846 %
847 % NOTE: "plus:1" is equivalent to a "Diamond" kernel.
848 %
849 % Ring:{radius1},{radius2}[,{scale}]
850 % A ring of the values given that falls between the two radii.
851 % Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
852 % This is the 'edge' pixels of the default "Disk" kernel,
853 % More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
854 %
855 % Hit and Miss Kernels
856 %
857 % Peak:radius1,radius2
858 % Find any peak larger than the pixels the fall between the two radii.
859 % The default ring of pixels is as per "Ring".
860 % Edges
861 % Find flat orthogonal edges of a binary shape
862 % Corners
863 % Find 90 degree corners of a binary shape
864 % Diagonals:type
865 % A special kernel to thin the 'outside' of diagonals
866 % LineEnds:type
867 % Find end points of lines (for pruning a skeletion)
868 % Two types of lines ends (default to both) can be searched for
869 % Type 0: All line ends
870 % Type 1: single kernel for 4-conneected line ends
871 % Type 2: single kernel for simple line ends
872 % LineJunctions
873 % Find three line junctions (within a skeletion)
874 % Type 0: all line junctions
875 % Type 1: Y Junction kernel
876 % Type 2: Diagonal T Junction kernel
877 % Type 3: Orthogonal T Junction kernel
878 % Type 4: Diagonal X Junction kernel
879 % Type 5: Orthogonal + Junction kernel
880 % Ridges:type
881 % Find single pixel ridges or thin lines
882 % Type 1: Fine single pixel thick lines and ridges
883 % Type 2: Find two pixel thick lines and ridges
884 % ConvexHull
885 % Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
886 % Skeleton:type
887 % Traditional skeleton generating kernels.
888 % Type 1: Tradional Skeleton kernel (4 connected skeleton)
889 % Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
890 % Type 3: Thinning skeleton based on a ressearch paper by
891 % Dan S. Bloomberg (Default Type)
892 % ThinSE:type
893 % A huge variety of Thinning Kernels designed to preserve conectivity.
894 % many other kernel sets use these kernels as source definitions.
895 % Type numbers are 41-49, 81-89, 481, and 482 which are based on
896 % the super and sub notations used in the source research paper.
897 %
898 % Distance Measuring Kernels
899 %
900 % Different types of distance measuring methods, which are used with the
901 % a 'Distance' morphology method for generating a gradient based on
902 % distance from an edge of a binary shape, though there is a technique
903 % for handling a anti-aliased shape.
904 %
905 % See the 'Distance' Morphological Method, for information of how it is
906 % applied.
907 %
908 % Chebyshev:[{radius}][x{scale}[%!]]
909 % Chebyshev Distance (also known as Tchebychev or Chessboard distance)
910 % is a value of one to any neighbour, orthogonal or diagonal. One why
911 % of thinking of it is the number of squares a 'King' or 'Queen' in
912 % chess needs to traverse reach any other position on a chess board.
913 % It results in a 'square' like distance function, but one where
914 % diagonals are given a value that is closer than expected.
915 %
916 % Manhattan:[{radius}][x{scale}[%!]]
917 % Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
918 % Cab distance metric), it is the distance needed when you can only
919 % travel in horizontal or vertical directions only. It is the
920 % distance a 'Rook' in chess would have to travel, and results in a
921 % diamond like distances, where diagonals are further than expected.
922 %
923 % Octagonal:[{radius}][x{scale}[%!]]
924 % An interleving of Manhatten and Chebyshev metrics producing an
925 % increasing octagonally shaped distance. Distances matches those of
926 % the "Octagon" shaped kernel of the same radius. The minimum radius
927 % and default is 2, producing a 5x5 kernel.
928 %
929 % Euclidean:[{radius}][x{scale}[%!]]
930 % Euclidean distance is the 'direct' or 'as the crow flys' distance.
931 % However by default the kernel size only has a radius of 1, which
932 % limits the distance to 'Knight' like moves, with only orthogonal and
933 % diagonal measurements being correct. As such for the default kernel
934 % you will get octagonal like distance function.
935 %
936 % However using a larger radius such as "Euclidean:4" you will get a
937 % much smoother distance gradient from the edge of the shape. Especially
938 % if the image is pre-processed to include any anti-aliasing pixels.
939 % Of course a larger kernel is slower to use, and not always needed.
940 %
941 % The first three Distance Measuring Kernels will only generate distances
942 % of exact multiples of {scale} in binary images. As such you can use a
943 % scale of 1 without loosing any information. However you also need some
944 % scaling when handling non-binary anti-aliased shapes.
945 %
946 % The "Euclidean" Distance Kernel however does generate a non-integer
947 % fractional results, and as such scaling is vital even for binary shapes.
948 %
949 */
950
AcquireKernelBuiltIn(const KernelInfoType type,const GeometryInfo * args,ExceptionInfo * exception)951 MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
952 const GeometryInfo *args,ExceptionInfo *exception)
953 {
954 KernelInfo
955 *kernel;
956
957 ssize_t
958 i;
959
960 ssize_t
961 u,
962 v;
963
964 double
965 nan = sqrt((double)-1.0); /* Special Value : Not A Number */
966
967 /* Generate a new empty kernel if needed */
968 kernel=(KernelInfo *) NULL;
969 switch(type) {
970 case UndefinedKernel: /* These should not call this function */
971 case UserDefinedKernel:
972 assert("Should not call this function" != (char *) NULL);
973 break;
974 case LaplacianKernel: /* Named Descrete 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 OctangonalKernel:
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()' (intergral 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*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*kernel->width] -= 1.0;
1102 }
1103
1104 if ( type == LoGKernel )
1105 { /* Calculate a Laplacian of a Gaussian - Or Mexician 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*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*kernel->width] = 1.0;
1190 #else
1191 /* Direct calculation without curve averaging
1192 This is equivelent 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 ** Becuase 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 comples
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*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-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-u-1) ));
1323 }
1324 kernel->minimum = 1.0;
1325 kernel->maximum = kernel->values[kernel->x+kernel->y*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*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 */
1864 KernelInfo
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 /* Unfortunatally we can not yet rotate a non-square kernel */
1884 /* But then we can't flip a non-symetrical 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 {
1931 KernelInfo
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 Thransformations"
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 kernel->type = type;
1987 kernel->next->type = type;
1988 kernel->next->next->type = type;
1989 ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1990 break;
1991 }
1992 break;
1993 }
1994 case ThinSEKernel:
1995 { /* Special kernels for general thinning, while preserving connections
1996 ** "Connectivity-Preserving Morphological Image Thransformations"
1997 ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1998 ** http://www.leptonica.com/papers/conn.pdf
1999 ** And
2000 ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2001 **
2002 ** Note kernels do not specify the origin pixel, allowing them
2003 ** to be used for both thickening and thinning operations.
2004 */
2005 switch ( (int) args->rho ) {
2006 /* SE for 4-connected thinning */
2007 case 41: /* SE_4_1 */
2008 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
2009 break;
2010 case 42: /* SE_4_2 */
2011 kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
2012 break;
2013 case 43: /* SE_4_3 */
2014 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
2015 break;
2016 case 44: /* SE_4_4 */
2017 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
2018 break;
2019 case 45: /* SE_4_5 */
2020 kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
2021 break;
2022 case 46: /* SE_4_6 */
2023 kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
2024 break;
2025 case 47: /* SE_4_7 */
2026 kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
2027 break;
2028 case 48: /* SE_4_8 */
2029 kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
2030 break;
2031 case 49: /* SE_4_9 */
2032 kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
2033 break;
2034 /* SE for 8-connected thinning - negatives of the above */
2035 case 81: /* SE_8_0 */
2036 kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
2037 break;
2038 case 82: /* SE_8_2 */
2039 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
2040 break;
2041 case 83: /* SE_8_3 */
2042 kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
2043 break;
2044 case 84: /* SE_8_4 */
2045 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
2046 break;
2047 case 85: /* SE_8_5 */
2048 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2049 break;
2050 case 86: /* SE_8_6 */
2051 kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2052 break;
2053 case 87: /* SE_8_7 */
2054 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2055 break;
2056 case 88: /* SE_8_8 */
2057 kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2058 break;
2059 case 89: /* SE_8_9 */
2060 kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2061 break;
2062 /* Special combined SE kernels */
2063 case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2064 kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2065 break;
2066 case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2067 kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2068 break;
2069 case 481: /* SE_48_1 - General Connected Corner Kernel */
2070 kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2071 break;
2072 default:
2073 case 482: /* SE_48_2 - General Edge Kernel */
2074 kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2075 break;
2076 }
2077 if (kernel == (KernelInfo *) NULL)
2078 return(kernel);
2079 kernel->type = type;
2080 RotateKernelInfo(kernel, args->sigma);
2081 break;
2082 }
2083 /*
2084 Distance Measuring Kernels
2085 */
2086 case ChebyshevKernel:
2087 {
2088 if (args->rho < 1.0)
2089 kernel->width = kernel->height = 3; /* default radius = 1 */
2090 else
2091 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2092 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2093
2094 kernel->values=(MagickRealType *) MagickAssumeAligned(
2095 AcquireAlignedMemory(kernel->width,kernel->height*
2096 sizeof(*kernel->values)));
2097 if (kernel->values == (MagickRealType *) NULL)
2098 return(DestroyKernelInfo(kernel));
2099
2100 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2101 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2102 kernel->positive_range += ( kernel->values[i] =
2103 args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2104 kernel->maximum = kernel->values[0];
2105 break;
2106 }
2107 case ManhattanKernel:
2108 {
2109 if (args->rho < 1.0)
2110 kernel->width = kernel->height = 3; /* default radius = 1 */
2111 else
2112 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2113 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2114
2115 kernel->values=(MagickRealType *) MagickAssumeAligned(
2116 AcquireAlignedMemory(kernel->width,kernel->height*
2117 sizeof(*kernel->values)));
2118 if (kernel->values == (MagickRealType *) NULL)
2119 return(DestroyKernelInfo(kernel));
2120
2121 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2122 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2123 kernel->positive_range += ( kernel->values[i] =
2124 args->sigma*(labs((long) u)+labs((long) v)) );
2125 kernel->maximum = kernel->values[0];
2126 break;
2127 }
2128 case OctagonalKernel:
2129 {
2130 if (args->rho < 2.0)
2131 kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2132 else
2133 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2134 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2135
2136 kernel->values=(MagickRealType *) MagickAssumeAligned(
2137 AcquireAlignedMemory(kernel->width,kernel->height*
2138 sizeof(*kernel->values)));
2139 if (kernel->values == (MagickRealType *) NULL)
2140 return(DestroyKernelInfo(kernel));
2141
2142 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2143 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2144 {
2145 double
2146 r1 = MagickMax(fabs((double)u),fabs((double)v)),
2147 r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2148 kernel->positive_range += kernel->values[i] =
2149 args->sigma*MagickMax(r1,r2);
2150 }
2151 kernel->maximum = kernel->values[0];
2152 break;
2153 }
2154 case EuclideanKernel:
2155 {
2156 if (args->rho < 1.0)
2157 kernel->width = kernel->height = 3; /* default radius = 1 */
2158 else
2159 kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2160 kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2161
2162 kernel->values=(MagickRealType *) MagickAssumeAligned(
2163 AcquireAlignedMemory(kernel->width,kernel->height*
2164 sizeof(*kernel->values)));
2165 if (kernel->values == (MagickRealType *) NULL)
2166 return(DestroyKernelInfo(kernel));
2167
2168 for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2169 for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2170 kernel->positive_range += ( kernel->values[i] =
2171 args->sigma*sqrt((double)(u*u+v*v)) );
2172 kernel->maximum = kernel->values[0];
2173 break;
2174 }
2175 default:
2176 {
2177 /* No-Op Kernel - Basically just a single pixel on its own */
2178 kernel=ParseKernelArray("1:1");
2179 if (kernel == (KernelInfo *) NULL)
2180 return(kernel);
2181 kernel->type = UndefinedKernel;
2182 break;
2183 }
2184 break;
2185 }
2186 return(kernel);
2187 }
2188
2189 /*
2190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2191 % %
2192 % %
2193 % %
2194 % C l o n e K e r n e l I n f o %
2195 % %
2196 % %
2197 % %
2198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2199 %
2200 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2201 % can be modified without effecting the original. The cloned kernel should
2202 % be destroyed using DestoryKernelInfo() when no longer needed.
2203 %
2204 % The format of the CloneKernelInfo method is:
2205 %
2206 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2207 %
2208 % A description of each parameter follows:
2209 %
2210 % o kernel: the Morphology/Convolution kernel to be cloned
2211 %
2212 */
CloneKernelInfo(const KernelInfo * kernel)2213 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2214 {
2215 ssize_t
2216 i;
2217
2218 KernelInfo
2219 *new_kernel;
2220
2221 assert(kernel != (KernelInfo *) NULL);
2222 new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2223 if (new_kernel == (KernelInfo *) NULL)
2224 return(new_kernel);
2225 *new_kernel=(*kernel); /* copy values in structure */
2226
2227 /* replace the values with a copy of the values */
2228 new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2229 AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2230 if (new_kernel->values == (MagickRealType *) NULL)
2231 return(DestroyKernelInfo(new_kernel));
2232 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2233 new_kernel->values[i]=kernel->values[i];
2234
2235 /* Also clone the next kernel in the kernel list */
2236 if ( kernel->next != (KernelInfo *) NULL ) {
2237 new_kernel->next = CloneKernelInfo(kernel->next);
2238 if ( new_kernel->next == (KernelInfo *) NULL )
2239 return(DestroyKernelInfo(new_kernel));
2240 }
2241
2242 return(new_kernel);
2243 }
2244
2245 /*
2246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2247 % %
2248 % %
2249 % %
2250 % D e s t r o y K e r n e l I n f o %
2251 % %
2252 % %
2253 % %
2254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2255 %
2256 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2257 % kernel.
2258 %
2259 % The format of the DestroyKernelInfo method is:
2260 %
2261 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2262 %
2263 % A description of each parameter follows:
2264 %
2265 % o kernel: the Morphology/Convolution kernel to be destroyed
2266 %
2267 */
DestroyKernelInfo(KernelInfo * kernel)2268 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2269 {
2270 assert(kernel != (KernelInfo *) NULL);
2271 if (kernel->next != (KernelInfo *) NULL)
2272 kernel->next=DestroyKernelInfo(kernel->next);
2273 kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2274 kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2275 return(kernel);
2276 }
2277
2278 /*
2279 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2280 % %
2281 % %
2282 % %
2283 + E x p a n d M i r r o r K e r n e l I n f o %
2284 % %
2285 % %
2286 % %
2287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2288 %
2289 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2290 % sequence of 90-degree rotated kernels but providing a reflected 180
2291 % rotatation, before the -/+ 90-degree rotations.
2292 %
2293 % This special rotation order produces a better, more symetrical thinning of
2294 % objects.
2295 %
2296 % The format of the ExpandMirrorKernelInfo method is:
2297 %
2298 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2299 %
2300 % A description of each parameter follows:
2301 %
2302 % o kernel: the Morphology/Convolution kernel
2303 %
2304 % This function is only internel to this module, as it is not finalized,
2305 % especially with regard to non-orthogonal angles, and rotation of larger
2306 % 2D kernels.
2307 */
2308
2309 #if 0
2310 static void FlopKernelInfo(KernelInfo *kernel)
2311 { /* Do a Flop by reversing each row. */
2312 size_t
2313 y;
2314 ssize_t
2315 x,r;
2316 double
2317 *k,t;
2318
2319 for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2320 for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2321 t=k[x], k[x]=k[r], k[r]=t;
2322
2323 kernel->x = kernel->width - kernel->x - 1;
2324 angle = fmod(angle+180.0, 360.0);
2325 }
2326 #endif
2327
ExpandMirrorKernelInfo(KernelInfo * kernel)2328 static void ExpandMirrorKernelInfo(KernelInfo *kernel)
2329 {
2330 KernelInfo
2331 *clone,
2332 *last;
2333
2334 last = kernel;
2335
2336 clone = CloneKernelInfo(last);
2337 if (clone == (KernelInfo *) NULL)
2338 return;
2339 RotateKernelInfo(clone, 180); /* flip */
2340 LastKernelInfo(last)->next = clone;
2341 last = clone;
2342
2343 clone = CloneKernelInfo(last);
2344 if (clone == (KernelInfo *) NULL)
2345 return;
2346 RotateKernelInfo(clone, 90); /* transpose */
2347 LastKernelInfo(last)->next = clone;
2348 last = clone;
2349
2350 clone = CloneKernelInfo(last);
2351 if (clone == (KernelInfo *) NULL)
2352 return;
2353 RotateKernelInfo(clone, 180); /* flop */
2354 LastKernelInfo(last)->next = clone;
2355
2356 return;
2357 }
2358
2359 /*
2360 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2361 % %
2362 % %
2363 % %
2364 + E x p a n d R o t a t e K e r n e l I n f o %
2365 % %
2366 % %
2367 % %
2368 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2369 %
2370 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2371 % incrementally by the angle given, until the kernel repeats.
2372 %
2373 % WARNING: 45 degree rotations only works for 3x3 kernels.
2374 % While 90 degree roatations only works for linear and square kernels
2375 %
2376 % The format of the ExpandRotateKernelInfo method is:
2377 %
2378 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2379 %
2380 % A description of each parameter follows:
2381 %
2382 % o kernel: the Morphology/Convolution kernel
2383 %
2384 % o angle: angle to rotate in degrees
2385 %
2386 % This function is only internel to this module, as it is not finalized,
2387 % especially with regard to non-orthogonal angles, and rotation of larger
2388 % 2D kernels.
2389 */
2390
2391 /* Internal Routine - Return true if two kernels are the same */
SameKernelInfo(const KernelInfo * kernel1,const KernelInfo * kernel2)2392 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
2393 const KernelInfo *kernel2)
2394 {
2395 size_t
2396 i;
2397
2398 /* check size and origin location */
2399 if ( kernel1->width != kernel2->width
2400 || kernel1->height != kernel2->height
2401 || kernel1->x != kernel2->x
2402 || kernel1->y != kernel2->y )
2403 return MagickFalse;
2404
2405 /* check actual kernel values */
2406 for (i=0; i < (kernel1->width*kernel1->height); i++) {
2407 /* Test for Nan equivalence */
2408 if ( IsNaN(kernel1->values[i]) && !IsNaN(kernel2->values[i]) )
2409 return MagickFalse;
2410 if ( IsNaN(kernel2->values[i]) && !IsNaN(kernel1->values[i]) )
2411 return MagickFalse;
2412 /* Test actual values are equivalent */
2413 if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2414 return MagickFalse;
2415 }
2416
2417 return MagickTrue;
2418 }
2419
ExpandRotateKernelInfo(KernelInfo * kernel,const double angle)2420 static void ExpandRotateKernelInfo(KernelInfo *kernel,const double angle)
2421 {
2422 KernelInfo
2423 *clone_info,
2424 *last;
2425
2426 clone_info=(KernelInfo *) NULL;
2427 last=kernel;
2428 DisableMSCWarning(4127)
2429 while (1) {
2430 RestoreMSCWarning
2431 clone_info=CloneKernelInfo(last);
2432 if (clone_info == (KernelInfo *) NULL)
2433 break;
2434 RotateKernelInfo(clone_info,angle);
2435 if (SameKernelInfo(kernel,clone_info) != MagickFalse)
2436 break;
2437 LastKernelInfo(last)->next=clone_info;
2438 last=clone_info;
2439 }
2440 if (clone_info != (KernelInfo *) NULL)
2441 clone_info=DestroyKernelInfo(clone_info); /* kernel repeated - junk */
2442 return;
2443 }
2444
2445 /*
2446 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2447 % %
2448 % %
2449 % %
2450 + C a l c M e t a K e r n a l I n f o %
2451 % %
2452 % %
2453 % %
2454 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2455 %
2456 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2457 % using the kernel values. This should only ne used if it is not possible to
2458 % calculate that meta-data in some easier way.
2459 %
2460 % It is important that the meta-data is correct before ScaleKernelInfo() is
2461 % used to perform kernel normalization.
2462 %
2463 % The format of the CalcKernelMetaData method is:
2464 %
2465 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2466 %
2467 % A description of each parameter follows:
2468 %
2469 % o kernel: the Morphology/Convolution kernel to modify
2470 %
2471 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2472 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2473 % however is not true for flat-shaped morphological kernels.
2474 %
2475 % WARNING: Only the specific kernel pointed to is modified, not a list of
2476 % multiple kernels.
2477 %
2478 % This is an internal function and not expected to be useful outside this
2479 % module. This could change however.
2480 */
CalcKernelMetaData(KernelInfo * kernel)2481 static void CalcKernelMetaData(KernelInfo *kernel)
2482 {
2483 size_t
2484 i;
2485
2486 kernel->minimum = kernel->maximum = 0.0;
2487 kernel->negative_range = kernel->positive_range = 0.0;
2488 for (i=0; i < (kernel->width*kernel->height); i++)
2489 {
2490 if ( fabs(kernel->values[i]) < MagickEpsilon )
2491 kernel->values[i] = 0.0;
2492 ( kernel->values[i] < 0)
2493 ? ( kernel->negative_range += kernel->values[i] )
2494 : ( kernel->positive_range += kernel->values[i] );
2495 Minimize(kernel->minimum, kernel->values[i]);
2496 Maximize(kernel->maximum, kernel->values[i]);
2497 }
2498
2499 return;
2500 }
2501
2502 /*
2503 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2504 % %
2505 % %
2506 % %
2507 % M o r p h o l o g y A p p l y %
2508 % %
2509 % %
2510 % %
2511 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2512 %
2513 % MorphologyApply() applies a morphological method, multiple times using
2514 % a list of multiple kernels. This is the method that should be called by
2515 % other 'operators' that internally use morphology operations as part of
2516 % their processing.
2517 %
2518 % It is basically equivalent to as MorphologyImage() (see below) but without
2519 % any user controls. This allows internel programs to use this method to
2520 % perform a specific task without possible interference by any API user
2521 % supplied settings.
2522 %
2523 % It is MorphologyImage() task to extract any such user controls, and
2524 % pass them to this function for processing.
2525 %
2526 % More specifically all given kernels should already be scaled, normalised,
2527 % and blended appropriatally before being parred to this routine. The
2528 % appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2529 %
2530 % The format of the MorphologyApply method is:
2531 %
2532 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2533 % const ssize_t iterations,const KernelInfo *kernel,
2534 % const CompositeMethod compose,const double bias,
2535 % ExceptionInfo *exception)
2536 %
2537 % A description of each parameter follows:
2538 %
2539 % o image: the source image
2540 %
2541 % o method: the morphology method to be applied.
2542 %
2543 % o iterations: apply the operation this many times (or no change).
2544 % A value of -1 means loop until no change found.
2545 % How this is applied may depend on the morphology method.
2546 % Typically this is a value of 1.
2547 %
2548 % o channel: the channel type.
2549 %
2550 % o kernel: An array of double representing the morphology kernel.
2551 %
2552 % o compose: How to handle or merge multi-kernel results.
2553 % If 'UndefinedCompositeOp' use default for the Morphology method.
2554 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2555 % Otherwise merge the results using the compose method given.
2556 %
2557 % o bias: Convolution Output Bias.
2558 %
2559 % o exception: return any errors or warnings in this structure.
2560 %
2561 */
MorphologyPrimitive(const Image * image,Image * morphology_image,const MorphologyMethod method,const KernelInfo * kernel,const double bias,ExceptionInfo * exception)2562 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2563 const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2564 ExceptionInfo *exception)
2565 {
2566 #define MorphologyTag "Morphology/Image"
2567
2568 CacheView
2569 *image_view,
2570 *morphology_view;
2571
2572 OffsetInfo
2573 offset;
2574
2575 ssize_t
2576 j,
2577 y;
2578
2579 size_t
2580 *changes,
2581 changed,
2582 width;
2583
2584 MagickBooleanType
2585 status;
2586
2587 MagickOffsetType
2588 progress;
2589
2590 assert(image != (Image *) NULL);
2591 assert(image->signature == MagickCoreSignature);
2592 assert(morphology_image != (Image *) NULL);
2593 assert(morphology_image->signature == MagickCoreSignature);
2594 assert(kernel != (KernelInfo *) NULL);
2595 assert(kernel->signature == MagickCoreSignature);
2596 assert(exception != (ExceptionInfo *) NULL);
2597 assert(exception->signature == MagickCoreSignature);
2598 status=MagickTrue;
2599 progress=0;
2600 image_view=AcquireVirtualCacheView(image,exception);
2601 morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2602 width=image->columns+kernel->width-1;
2603 offset.x=0;
2604 offset.y=0;
2605 switch (method)
2606 {
2607 case ConvolveMorphology:
2608 case DilateMorphology:
2609 case DilateIntensityMorphology:
2610 case IterativeDistanceMorphology:
2611 {
2612 /*
2613 Kernel needs to used with reflection about origin.
2614 */
2615 offset.x=(ssize_t) kernel->width-kernel->x-1;
2616 offset.y=(ssize_t) kernel->height-kernel->y-1;
2617 break;
2618 }
2619 case ErodeMorphology:
2620 case ErodeIntensityMorphology:
2621 case HitAndMissMorphology:
2622 case ThinningMorphology:
2623 case ThickenMorphology:
2624 {
2625 offset.x=kernel->x;
2626 offset.y=kernel->y;
2627 break;
2628 }
2629 default:
2630 {
2631 assert("Not a Primitive Morphology Method" != (char *) NULL);
2632 break;
2633 }
2634 }
2635 changed=0;
2636 changes=(size_t *) AcquireQuantumMemory(GetOpenMPMaximumThreads(),
2637 sizeof(*changes));
2638 if (changes == (size_t *) NULL)
2639 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2640 for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
2641 changes[j]=0;
2642
2643 if ((method == ConvolveMorphology) && (kernel->width == 1))
2644 {
2645 ssize_t
2646 x;
2647
2648 /*
2649 Special handling (for speed) of vertical (blur) kernels. This performs
2650 its handling in columns rather than in rows. This is only done
2651 for convolve as it is the only method that generates very large 1-D
2652 vertical kernels (such as a 'BlurKernel')
2653 */
2654 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2655 #pragma omp parallel for schedule(static) shared(progress,status) \
2656 magick_number_threads(image,morphology_image,image->columns,1)
2657 #endif
2658 for (x=0; x < (ssize_t) image->columns; x++)
2659 {
2660 const int
2661 id = GetOpenMPThreadId();
2662
2663 const Quantum
2664 *magick_restrict p;
2665
2666 Quantum
2667 *magick_restrict q;
2668
2669 ssize_t
2670 r;
2671
2672 ssize_t
2673 center;
2674
2675 if (status == MagickFalse)
2676 continue;
2677 p=GetCacheViewVirtualPixels(image_view,x,-offset.y,1,image->rows+
2678 kernel->height-1,exception);
2679 q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2680 morphology_image->rows,exception);
2681 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2682 {
2683 status=MagickFalse;
2684 continue;
2685 }
2686 center=(ssize_t) GetPixelChannels(image)*offset.y;
2687 for (r=0; r < (ssize_t) image->rows; r++)
2688 {
2689 ssize_t
2690 i;
2691
2692 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2693 {
2694 double
2695 alpha,
2696 gamma,
2697 pixel;
2698
2699 PixelChannel
2700 channel;
2701
2702 PixelTrait
2703 morphology_traits,
2704 traits;
2705
2706 const MagickRealType
2707 *magick_restrict k;
2708
2709 const Quantum
2710 *magick_restrict pixels;
2711
2712 ssize_t
2713 v;
2714
2715 size_t
2716 count;
2717
2718 channel=GetPixelChannelChannel(image,i);
2719 traits=GetPixelChannelTraits(image,channel);
2720 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2721 if ((traits == UndefinedPixelTrait) ||
2722 (morphology_traits == UndefinedPixelTrait))
2723 continue;
2724 if ((traits & CopyPixelTrait) != 0)
2725 {
2726 SetPixelChannel(morphology_image,channel,p[center+i],q);
2727 continue;
2728 }
2729 k=(&kernel->values[kernel->height-1]);
2730 pixels=p;
2731 pixel=bias;
2732 gamma=1.0;
2733 count=0;
2734 if (((image->alpha_trait & BlendPixelTrait) == 0) ||
2735 ((morphology_traits & BlendPixelTrait) == 0))
2736 for (v=0; v < (ssize_t) kernel->height; v++)
2737 {
2738 if (!IsNaN(*k))
2739 {
2740 pixel+=(*k)*pixels[i];
2741 count++;
2742 }
2743 k--;
2744 pixels+=GetPixelChannels(image);
2745 }
2746 else
2747 {
2748 gamma=0.0;
2749 for (v=0; v < (ssize_t) kernel->height; v++)
2750 {
2751 if (!IsNaN(*k))
2752 {
2753 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2754 pixel+=alpha*(*k)*pixels[i];
2755 gamma+=alpha*(*k);
2756 count++;
2757 }
2758 k--;
2759 pixels+=GetPixelChannels(image);
2760 }
2761 }
2762 if (fabs(pixel-p[center+i]) > MagickEpsilon)
2763 changes[id]++;
2764 gamma=PerceptibleReciprocal(gamma);
2765 if (count != 0)
2766 gamma*=(double) kernel->height/count;
2767 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2768 pixel),q);
2769 }
2770 p+=GetPixelChannels(image);
2771 q+=GetPixelChannels(morphology_image);
2772 }
2773 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2774 status=MagickFalse;
2775 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2776 {
2777 MagickBooleanType
2778 proceed;
2779
2780 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2781 #pragma omp atomic
2782 #endif
2783 progress++;
2784 proceed=SetImageProgress(image,MorphologyTag,progress,image->rows);
2785 if (proceed == MagickFalse)
2786 status=MagickFalse;
2787 }
2788 }
2789 morphology_image->type=image->type;
2790 morphology_view=DestroyCacheView(morphology_view);
2791 image_view=DestroyCacheView(image_view);
2792 for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
2793 changed+=changes[j];
2794 changes=(size_t *) RelinquishMagickMemory(changes);
2795 return(status ? (ssize_t) changed : 0);
2796 }
2797 /*
2798 Normal handling of horizontal or rectangular kernels (row by row).
2799 */
2800 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2801 #pragma omp parallel for schedule(static) shared(progress,status) \
2802 magick_number_threads(image,morphology_image,image->rows,1)
2803 #endif
2804 for (y=0; y < (ssize_t) image->rows; y++)
2805 {
2806 const int
2807 id = GetOpenMPThreadId();
2808
2809 const Quantum
2810 *magick_restrict p;
2811
2812 Quantum
2813 *magick_restrict q;
2814
2815 ssize_t
2816 x;
2817
2818 ssize_t
2819 center;
2820
2821 if (status == MagickFalse)
2822 continue;
2823 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,
2824 kernel->height,exception);
2825 q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2826 1,exception);
2827 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2828 {
2829 status=MagickFalse;
2830 continue;
2831 }
2832 center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
2833 GetPixelChannels(image)*offset.x);
2834 for (x=0; x < (ssize_t) image->columns; x++)
2835 {
2836 ssize_t
2837 i;
2838
2839 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2840 {
2841 double
2842 alpha,
2843 gamma,
2844 intensity,
2845 maximum,
2846 minimum,
2847 pixel;
2848
2849 PixelChannel
2850 channel;
2851
2852 PixelTrait
2853 morphology_traits,
2854 traits;
2855
2856 const MagickRealType
2857 *magick_restrict k;
2858
2859 const Quantum
2860 *magick_restrict pixels,
2861 *magick_restrict quantum_pixels;
2862
2863 ssize_t
2864 u;
2865
2866 size_t
2867 count;
2868
2869 ssize_t
2870 v;
2871
2872 channel=GetPixelChannelChannel(image,i);
2873 traits=GetPixelChannelTraits(image,channel);
2874 morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2875 if ((traits == UndefinedPixelTrait) ||
2876 (morphology_traits == UndefinedPixelTrait))
2877 continue;
2878 if ((traits & CopyPixelTrait) != 0)
2879 {
2880 SetPixelChannel(morphology_image,channel,p[center+i],q);
2881 continue;
2882 }
2883 pixels=p;
2884 quantum_pixels=(const Quantum *) NULL;
2885 maximum=0.0;
2886 minimum=(double) QuantumRange;
2887 switch (method)
2888 {
2889 case ConvolveMorphology:
2890 {
2891 pixel=bias;
2892 break;
2893 }
2894 case DilateMorphology:
2895 case ErodeIntensityMorphology:
2896 {
2897 pixel=0.0;
2898 break;
2899 }
2900 case HitAndMissMorphology:
2901 case ErodeMorphology:
2902 {
2903 pixel=QuantumRange;
2904 break;
2905 }
2906 default:
2907 {
2908 pixel=(double) p[center+i];
2909 break;
2910 }
2911 }
2912 count=0;
2913 gamma=1.0;
2914 switch (method)
2915 {
2916 case ConvolveMorphology:
2917 {
2918 /*
2919 Weighted Average of pixels using reflected kernel
2920
2921 For correct working of this operation for asymetrical kernels,
2922 the kernel needs to be applied in its reflected form. That is
2923 its values needs to be reversed.
2924
2925 Correlation is actually the same as this but without reflecting
2926 the kernel, and thus 'lower-level' that Convolution. However as
2927 Convolution is the more common method used, and it does not
2928 really cost us much in terms of processing to use a reflected
2929 kernel, so it is Convolution that is implemented.
2930
2931 Correlation will have its kernel reflected before calling this
2932 function to do a Convolve.
2933
2934 For more details of Correlation vs Convolution see
2935 http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2936 */
2937 k=(&kernel->values[kernel->width*kernel->height-1]);
2938 if (((image->alpha_trait & BlendPixelTrait) == 0) ||
2939 ((morphology_traits & BlendPixelTrait) == 0))
2940 {
2941 /*
2942 No alpha blending.
2943 */
2944 for (v=0; v < (ssize_t) kernel->height; v++)
2945 {
2946 for (u=0; u < (ssize_t) kernel->width; u++)
2947 {
2948 if (!IsNaN(*k))
2949 {
2950 pixel+=(*k)*pixels[i];
2951 count++;
2952 }
2953 k--;
2954 pixels+=GetPixelChannels(image);
2955 }
2956 pixels+=(image->columns-1)*GetPixelChannels(image);
2957 }
2958 break;
2959 }
2960 /*
2961 Alpha blending.
2962 */
2963 gamma=0.0;
2964 for (v=0; v < (ssize_t) kernel->height; v++)
2965 {
2966 for (u=0; u < (ssize_t) kernel->width; u++)
2967 {
2968 if (!IsNaN(*k))
2969 {
2970 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2971 pixel+=alpha*(*k)*pixels[i];
2972 gamma+=alpha*(*k);
2973 count++;
2974 }
2975 k--;
2976 pixels+=GetPixelChannels(image);
2977 }
2978 pixels+=(image->columns-1)*GetPixelChannels(image);
2979 }
2980 break;
2981 }
2982 case ErodeMorphology:
2983 {
2984 /*
2985 Minimum value within kernel neighbourhood.
2986
2987 The kernel is not reflected for this operation. In normal
2988 Greyscale Morphology, the kernel value should be added
2989 to the real value, this is currently not done, due to the
2990 nature of the boolean kernels being used.
2991 */
2992 k=kernel->values;
2993 for (v=0; v < (ssize_t) kernel->height; v++)
2994 {
2995 for (u=0; u < (ssize_t) kernel->width; u++)
2996 {
2997 if (!IsNaN(*k) && (*k >= 0.5))
2998 {
2999 if ((double) pixels[i] < pixel)
3000 pixel=(double) pixels[i];
3001 }
3002 k++;
3003 pixels+=GetPixelChannels(image);
3004 }
3005 pixels+=(image->columns-1)*GetPixelChannels(image);
3006 }
3007 break;
3008 }
3009 case DilateMorphology:
3010 {
3011 /*
3012 Maximum value within kernel neighbourhood.
3013
3014 For correct working of this operation for asymetrical kernels,
3015 the kernel needs to be applied in its reflected form. That is
3016 its values needs to be reversed.
3017
3018 In normal Greyscale Morphology, the kernel value should be
3019 added to the real value, this is currently not done, due to the
3020 nature of the boolean kernels being used.
3021 */
3022 k=(&kernel->values[kernel->width*kernel->height-1]);
3023 for (v=0; v < (ssize_t) kernel->height; v++)
3024 {
3025 for (u=0; u < (ssize_t) kernel->width; u++)
3026 {
3027 if (!IsNaN(*k) && (*k > 0.5))
3028 {
3029 if ((double) pixels[i] > pixel)
3030 pixel=(double) pixels[i];
3031 }
3032 k--;
3033 pixels+=GetPixelChannels(image);
3034 }
3035 pixels+=(image->columns-1)*GetPixelChannels(image);
3036 }
3037 break;
3038 }
3039 case HitAndMissMorphology:
3040 case ThinningMorphology:
3041 case ThickenMorphology:
3042 {
3043 /*
3044 Minimum of foreground pixel minus maxumum of background pixels.
3045
3046 The kernel is not reflected for this operation, and consists
3047 of both foreground and background pixel neighbourhoods, 0.0 for
3048 background, and 1.0 for foreground with either Nan or 0.5 values
3049 for don't care.
3050
3051 This never produces a meaningless negative result. Such results
3052 cause Thinning/Thicken to not work correctly when used against a
3053 greyscale image.
3054 */
3055 k=kernel->values;
3056 for (v=0; v < (ssize_t) kernel->height; v++)
3057 {
3058 for (u=0; u < (ssize_t) kernel->width; u++)
3059 {
3060 if (!IsNaN(*k))
3061 {
3062 if (*k > 0.7)
3063 {
3064 if ((double) pixels[i] < pixel)
3065 pixel=(double) pixels[i];
3066 }
3067 else
3068 if (*k < 0.3)
3069 {
3070 if ((double) pixels[i] > maximum)
3071 maximum=(double) pixels[i];
3072 }
3073 count++;
3074 }
3075 k++;
3076 pixels+=GetPixelChannels(image);
3077 }
3078 pixels+=(image->columns-1)*GetPixelChannels(image);
3079 }
3080 pixel-=maximum;
3081 if (pixel < 0.0)
3082 pixel=0.0;
3083 if (method == ThinningMorphology)
3084 pixel=(double) p[center+i]-pixel;
3085 else
3086 if (method == ThickenMorphology)
3087 pixel+=(double) p[center+i]+pixel;
3088 break;
3089 }
3090 case ErodeIntensityMorphology:
3091 {
3092 /*
3093 Select pixel with minimum intensity within kernel neighbourhood.
3094
3095 The kernel is not reflected for this operation.
3096 */
3097 k=kernel->values;
3098 for (v=0; v < (ssize_t) kernel->height; v++)
3099 {
3100 for (u=0; u < (ssize_t) kernel->width; u++)
3101 {
3102 if (!IsNaN(*k) && (*k >= 0.5))
3103 {
3104 intensity=(double) GetPixelIntensity(image,pixels);
3105 if (intensity < minimum)
3106 {
3107 quantum_pixels=pixels;
3108 pixel=(double) pixels[i];
3109 minimum=intensity;
3110 }
3111 count++;
3112 }
3113 k++;
3114 pixels+=GetPixelChannels(image);
3115 }
3116 pixels+=(image->columns-1)*GetPixelChannels(image);
3117 }
3118 break;
3119 }
3120 case DilateIntensityMorphology:
3121 {
3122 /*
3123 Select pixel with maximum intensity within kernel neighbourhood.
3124
3125 The kernel is not reflected for this operation.
3126 */
3127 k=(&kernel->values[kernel->width*kernel->height-1]);
3128 for (v=0; v < (ssize_t) kernel->height; v++)
3129 {
3130 for (u=0; u < (ssize_t) kernel->width; u++)
3131 {
3132 if (!IsNaN(*k) && (*k >= 0.5))
3133 {
3134 intensity=(double) GetPixelIntensity(image,pixels);
3135 if (intensity > maximum)
3136 {
3137 pixel=(double) pixels[i];
3138 quantum_pixels=pixels;
3139 maximum=intensity;
3140 }
3141 count++;
3142 }
3143 k--;
3144 pixels+=GetPixelChannels(image);
3145 }
3146 pixels+=(image->columns-1)*GetPixelChannels(image);
3147 }
3148 break;
3149 }
3150 case IterativeDistanceMorphology:
3151 {
3152 /*
3153 Compute th iterative distance from black edge of a white image
3154 shape. Essentially white values are decreased to the smallest
3155 'distance from edge' it can find.
3156
3157 It works by adding kernel values to the neighbourhood, and
3158 select the minimum value found. The kernel is rotated before
3159 use, so kernel distances match resulting distances, when a user
3160 provided asymmetric kernel is applied.
3161
3162 This code is nearly identical to True GrayScale Morphology but
3163 not quite.
3164
3165 GreyDilate Kernel values added, maximum value found Kernel is
3166 rotated before use.
3167
3168 GrayErode: Kernel values subtracted and minimum value found No
3169 kernel rotation used.
3170
3171 Note the Iterative Distance method is essentially a
3172 GrayErode, but with negative kernel values, and kernel rotation
3173 applied.
3174 */
3175 k=(&kernel->values[kernel->width*kernel->height-1]);
3176 for (v=0; v < (ssize_t) kernel->height; v++)
3177 {
3178 for (u=0; u < (ssize_t) kernel->width; u++)
3179 {
3180 if (!IsNaN(*k))
3181 {
3182 if ((pixels[i]+(*k)) < pixel)
3183 pixel=(double) pixels[i]+(*k);
3184 count++;
3185 }
3186 k--;
3187 pixels+=GetPixelChannels(image);
3188 }
3189 pixels+=(image->columns-1)*GetPixelChannels(image);
3190 }
3191 break;
3192 }
3193 case UndefinedMorphology:
3194 default:
3195 break;
3196 }
3197 if (fabs(pixel-p[center+i]) > MagickEpsilon)
3198 changes[id]++;
3199 if (quantum_pixels != (const Quantum *) NULL)
3200 {
3201 SetPixelChannel(morphology_image,channel,quantum_pixels[i],q);
3202 continue;
3203 }
3204 gamma=PerceptibleReciprocal(gamma);
3205 if (count != 0)
3206 gamma*=(double) kernel->height*kernel->width/count;
3207 SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*pixel),q);
3208 }
3209 p+=GetPixelChannels(image);
3210 q+=GetPixelChannels(morphology_image);
3211 }
3212 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3213 status=MagickFalse;
3214 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3215 {
3216 MagickBooleanType
3217 proceed;
3218
3219 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3220 #pragma omp atomic
3221 #endif
3222 progress++;
3223 proceed=SetImageProgress(image,MorphologyTag,progress,image->rows);
3224 if (proceed == MagickFalse)
3225 status=MagickFalse;
3226 }
3227 }
3228 morphology_view=DestroyCacheView(morphology_view);
3229 image_view=DestroyCacheView(image_view);
3230 for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
3231 changed+=changes[j];
3232 changes=(size_t *) RelinquishMagickMemory(changes);
3233 return(status ? (ssize_t) changed : -1);
3234 }
3235
3236 /*
3237 This is almost identical to the MorphologyPrimative() function above, but
3238 applies the primitive directly to the actual image using two passes, once in
3239 each direction, with the results of the previous (and current) row being
3240 re-used.
3241
3242 That is after each row is 'Sync'ed' into the image, the next row makes use of
3243 those values as part of the calculation of the next row. It repeats, but
3244 going in the oppisite (bottom-up) direction.
3245
3246 Because of this 're-use of results' this function can not make use of multi-
3247 threaded, parellel processing.
3248 */
MorphologyPrimitiveDirect(Image * image,const MorphologyMethod method,const KernelInfo * kernel,ExceptionInfo * exception)3249 static ssize_t MorphologyPrimitiveDirect(Image *image,
3250 const MorphologyMethod method,const KernelInfo *kernel,
3251 ExceptionInfo *exception)
3252 {
3253 CacheView
3254 *morphology_view,
3255 *image_view;
3256
3257 MagickBooleanType
3258 status;
3259
3260 MagickOffsetType
3261 progress;
3262
3263 OffsetInfo
3264 offset;
3265
3266 size_t
3267 width,
3268 changed;
3269
3270 ssize_t
3271 y;
3272
3273 assert(image != (Image *) NULL);
3274 assert(image->signature == MagickCoreSignature);
3275 assert(kernel != (KernelInfo *) NULL);
3276 assert(kernel->signature == MagickCoreSignature);
3277 assert(exception != (ExceptionInfo *) NULL);
3278 assert(exception->signature == MagickCoreSignature);
3279 status=MagickTrue;
3280 changed=0;
3281 progress=0;
3282 switch(method)
3283 {
3284 case DistanceMorphology:
3285 case VoronoiMorphology:
3286 {
3287 /*
3288 Kernel reflected about origin.
3289 */
3290 offset.x=(ssize_t) kernel->width-kernel->x-1;
3291 offset.y=(ssize_t) kernel->height-kernel->y-1;
3292 break;
3293 }
3294 default:
3295 {
3296 offset.x=kernel->x;
3297 offset.y=kernel->y;
3298 break;
3299 }
3300 }
3301 /*
3302 Two views into same image, do not thread.
3303 */
3304 image_view=AcquireVirtualCacheView(image,exception);
3305 morphology_view=AcquireAuthenticCacheView(image,exception);
3306 width=image->columns+kernel->width-1;
3307 for (y=0; y < (ssize_t) image->rows; y++)
3308 {
3309 const Quantum
3310 *magick_restrict p;
3311
3312 Quantum
3313 *magick_restrict q;
3314
3315 ssize_t
3316 x;
3317
3318 /*
3319 Read virtual pixels, and authentic pixels, from the same image! We read
3320 using virtual to get virtual pixel handling, but write back into the same
3321 image.
3322
3323 Only top half of kernel is processed as we do a single pass downward
3324 through the image iterating the distance function as we go.
3325 */
3326 if (status == MagickFalse)
3327 continue;
3328 p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t)
3329 offset.y+1,exception);
3330 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3331 exception);
3332 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3333 {
3334 status=MagickFalse;
3335 continue;
3336 }
3337 for (x=0; x < (ssize_t) image->columns; x++)
3338 {
3339 ssize_t
3340 i;
3341
3342 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3343 {
3344 double
3345 pixel;
3346
3347 PixelChannel
3348 channel;
3349
3350 PixelTrait
3351 traits;
3352
3353 const MagickRealType
3354 *magick_restrict k;
3355
3356 const Quantum
3357 *magick_restrict pixels;
3358
3359 ssize_t
3360 u;
3361
3362 ssize_t
3363 v;
3364
3365 channel=GetPixelChannelChannel(image,i);
3366 traits=GetPixelChannelTraits(image,channel);
3367 if (traits == UndefinedPixelTrait)
3368 continue;
3369 if ((traits & CopyPixelTrait) != 0)
3370 continue;
3371 pixels=p;
3372 pixel=(double) QuantumRange;
3373 switch (method)
3374 {
3375 case DistanceMorphology:
3376 {
3377 k=(&kernel->values[kernel->width*kernel->height-1]);
3378 for (v=0; v <= offset.y; v++)
3379 {
3380 for (u=0; u < (ssize_t) kernel->width; u++)
3381 {
3382 if (!IsNaN(*k))
3383 {
3384 if ((pixels[i]+(*k)) < pixel)
3385 pixel=(double) pixels[i]+(*k);
3386 }
3387 k--;
3388 pixels+=GetPixelChannels(image);
3389 }
3390 pixels+=(image->columns-1)*GetPixelChannels(image);
3391 }
3392 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3393 pixels=q-offset.x*GetPixelChannels(image);
3394 for (u=0; u < offset.x; u++)
3395 {
3396 if (!IsNaN(*k) && ((x+u-offset.x) >= 0))
3397 {
3398 if ((pixels[i]+(*k)) < pixel)
3399 pixel=(double) pixels[i]+(*k);
3400 }
3401 k--;
3402 pixels+=GetPixelChannels(image);
3403 }
3404 break;
3405 }
3406 case VoronoiMorphology:
3407 {
3408 k=(&kernel->values[kernel->width*kernel->height-1]);
3409 for (v=0; v < offset.y; v++)
3410 {
3411 for (u=0; u < (ssize_t) kernel->width; u++)
3412 {
3413 if (!IsNaN(*k))
3414 {
3415 if ((pixels[i]+(*k)) < pixel)
3416 pixel=(double) pixels[i]+(*k);
3417 }
3418 k--;
3419 pixels+=GetPixelChannels(image);
3420 }
3421 pixels+=(image->columns-1)*GetPixelChannels(image);
3422 }
3423 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3424 pixels=q-offset.x*GetPixelChannels(image);
3425 for (u=0; u < offset.x; u++)
3426 {
3427 if (!IsNaN(*k) && ((x+u-offset.x) >= 0))
3428 {
3429 if ((pixels[i]+(*k)) < pixel)
3430 pixel=(double) pixels[i]+(*k);
3431 }
3432 k--;
3433 pixels+=GetPixelChannels(image);
3434 }
3435 break;
3436 }
3437 default:
3438 break;
3439 }
3440 if (fabs(pixel-q[i]) > MagickEpsilon)
3441 changed++;
3442 q[i]=ClampToQuantum(pixel);
3443 }
3444 p+=GetPixelChannels(image);
3445 q+=GetPixelChannels(image);
3446 }
3447 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3448 status=MagickFalse;
3449 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3450 {
3451 MagickBooleanType
3452 proceed;
3453
3454 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3455 #pragma omp atomic
3456 #endif
3457 progress++;
3458 proceed=SetImageProgress(image,MorphologyTag,progress,2*image->rows);
3459 if (proceed == MagickFalse)
3460 status=MagickFalse;
3461 }
3462 }
3463 morphology_view=DestroyCacheView(morphology_view);
3464 image_view=DestroyCacheView(image_view);
3465 /*
3466 Do the reverse pass through the image.
3467 */
3468 image_view=AcquireVirtualCacheView(image,exception);
3469 morphology_view=AcquireAuthenticCacheView(image,exception);
3470 for (y=(ssize_t) image->rows-1; y >= 0; y--)
3471 {
3472 const Quantum
3473 *magick_restrict p;
3474
3475 Quantum
3476 *magick_restrict q;
3477
3478 ssize_t
3479 x;
3480
3481 /*
3482 Read virtual pixels, and authentic pixels, from the same image. We
3483 read using virtual to get virtual pixel handling, but write back
3484 into the same image.
3485
3486 Only the bottom half of the kernel is processed as we up the image.
3487 */
3488 if (status == MagickFalse)
3489 continue;
3490 p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t)
3491 kernel->y+1,exception);
3492 q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3493 exception);
3494 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3495 {
3496 status=MagickFalse;
3497 continue;
3498 }
3499 p+=(image->columns-1)*GetPixelChannels(image);
3500 q+=(image->columns-1)*GetPixelChannels(image);
3501 for (x=(ssize_t) image->columns-1; x >= 0; x--)
3502 {
3503 ssize_t
3504 i;
3505
3506 for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3507 {
3508 double
3509 pixel;
3510
3511 PixelChannel
3512 channel;
3513
3514 PixelTrait
3515 traits;
3516
3517 const MagickRealType
3518 *magick_restrict k;
3519
3520 const Quantum
3521 *magick_restrict pixels;
3522
3523 ssize_t
3524 u;
3525
3526 ssize_t
3527 v;
3528
3529 channel=GetPixelChannelChannel(image,i);
3530 traits=GetPixelChannelTraits(image,channel);
3531 if (traits == UndefinedPixelTrait)
3532 continue;
3533 if ((traits & CopyPixelTrait) != 0)
3534 continue;
3535 pixels=p;
3536 pixel=(double) QuantumRange;
3537 switch (method)
3538 {
3539 case DistanceMorphology:
3540 {
3541 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3542 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3543 {
3544 for (u=0; u < (ssize_t) kernel->width; u++)
3545 {
3546 if (!IsNaN(*k))
3547 {
3548 if ((pixels[i]+(*k)) < pixel)
3549 pixel=(double) pixels[i]+(*k);
3550 }
3551 k--;
3552 pixels+=GetPixelChannels(image);
3553 }
3554 pixels+=(image->columns-1)*GetPixelChannels(image);
3555 }
3556 k=(&kernel->values[kernel->width*kernel->y+kernel->x-1]);
3557 pixels=q;
3558 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3559 {
3560 pixels+=GetPixelChannels(image);
3561 if (!IsNaN(*k) && ((x+u-offset.x) < (ssize_t) image->columns))
3562 {
3563 if ((pixels[i]+(*k)) < pixel)
3564 pixel=(double) pixels[i]+(*k);
3565 }
3566 k--;
3567 }
3568 break;
3569 }
3570 case VoronoiMorphology:
3571 {
3572 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3573 for (v=offset.y; v < (ssize_t) kernel->height; v++)
3574 {
3575 for (u=0; u < (ssize_t) kernel->width; u++)
3576 {
3577 if (!IsNaN(*k))
3578 {
3579 if ((pixels[i]+(*k)) < pixel)
3580 pixel=(double) pixels[i]+(*k);
3581 }
3582 k--;
3583 pixels+=GetPixelChannels(image);
3584 }
3585 pixels+=(image->columns-1)*GetPixelChannels(image);
3586 }
3587 k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3588 pixels=q;
3589 for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3590 {
3591 pixels+=GetPixelChannels(image);
3592 if (!IsNaN(*k) && ((x+u-offset.x) < (ssize_t) image->columns))
3593 {
3594 if ((pixels[i]+(*k)) < pixel)
3595 pixel=(double) pixels[i]+(*k);
3596 }
3597 k--;
3598 }
3599 break;
3600 }
3601 default:
3602 break;
3603 }
3604 if (fabs(pixel-q[i]) > MagickEpsilon)
3605 changed++;
3606 q[i]=ClampToQuantum(pixel);
3607 }
3608 p-=GetPixelChannels(image);
3609 q-=GetPixelChannels(image);
3610 }
3611 if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3612 status=MagickFalse;
3613 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3614 {
3615 MagickBooleanType
3616 proceed;
3617
3618 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3619 #pragma omp atomic
3620 #endif
3621 progress++;
3622 proceed=SetImageProgress(image,MorphologyTag,progress,2*image->rows);
3623 if (proceed == MagickFalse)
3624 status=MagickFalse;
3625 }
3626 }
3627 morphology_view=DestroyCacheView(morphology_view);
3628 image_view=DestroyCacheView(image_view);
3629 return(status ? (ssize_t) changed : -1);
3630 }
3631
3632 /*
3633 Apply a Morphology by calling one of the above low level primitive
3634 application functions. This function handles any iteration loops,
3635 composition or re-iteration of results, and compound morphology methods that
3636 is based on multiple low-level (staged) morphology methods.
3637
3638 Basically this provides the complex glue between the requested morphology
3639 method and raw low-level implementation (above).
3640 */
MorphologyApply(const Image * image,const MorphologyMethod method,const ssize_t iterations,const KernelInfo * kernel,const CompositeOperator compose,const double bias,ExceptionInfo * exception)3641 MagickPrivate Image *MorphologyApply(const Image *image,
3642 const MorphologyMethod method, const ssize_t iterations,
3643 const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3644 ExceptionInfo *exception)
3645 {
3646 CompositeOperator
3647 curr_compose;
3648
3649 Image
3650 *curr_image, /* Image we are working with or iterating */
3651 *work_image, /* secondary image for primitive iteration */
3652 *save_image, /* saved image - for 'edge' method only */
3653 *rslt_image; /* resultant image - after multi-kernel handling */
3654
3655 KernelInfo
3656 *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3657 *norm_kernel, /* the current normal un-reflected kernel */
3658 *rflt_kernel, /* the current reflected kernel (if needed) */
3659 *this_kernel; /* the kernel being applied */
3660
3661 MorphologyMethod
3662 primitive; /* the current morphology primitive being applied */
3663
3664 CompositeOperator
3665 rslt_compose; /* multi-kernel compose method for results to use */
3666
3667 MagickBooleanType
3668 special, /* do we use a direct modify function? */
3669 verbose; /* verbose output of results */
3670
3671 size_t
3672 method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3673 method_limit, /* maximum number of compound method iterations */
3674 kernel_number, /* Loop 2: the kernel number being applied */
3675 stage_loop, /* Loop 3: primitive loop for compound morphology */
3676 stage_limit, /* how many primitives are in this compound */
3677 kernel_loop, /* Loop 4: iterate the kernel over image */
3678 kernel_limit, /* number of times to iterate kernel */
3679 count, /* total count of primitive steps applied */
3680 kernel_changed, /* total count of changed using iterated kernel */
3681 method_changed; /* total count of changed over method iteration */
3682
3683 ssize_t
3684 changed; /* number pixels changed by last primitive operation */
3685
3686 char
3687 v_info[MagickPathExtent];
3688
3689 assert(image != (Image *) NULL);
3690 assert(image->signature == MagickCoreSignature);
3691 assert(kernel != (KernelInfo *) NULL);
3692 assert(kernel->signature == MagickCoreSignature);
3693 assert(exception != (ExceptionInfo *) NULL);
3694 assert(exception->signature == MagickCoreSignature);
3695
3696 count = 0; /* number of low-level morphology primitives performed */
3697 if ( iterations == 0 )
3698 return((Image *) NULL); /* null operation - nothing to do! */
3699
3700 kernel_limit = (size_t) iterations;
3701 if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3702 kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3703
3704 verbose = IsStringTrue(GetImageArtifact(image,"debug"));
3705
3706 /* initialise for cleanup */
3707 curr_image = (Image *) image;
3708 curr_compose = image->compose;
3709 (void) curr_compose;
3710 work_image = save_image = rslt_image = (Image *) NULL;
3711 reflected_kernel = (KernelInfo *) NULL;
3712
3713 /* Initialize specific methods
3714 * + which loop should use the given iteratations
3715 * + how many primitives make up the compound morphology
3716 * + multi-kernel compose method to use (by default)
3717 */
3718 method_limit = 1; /* just do method once, unless otherwise set */
3719 stage_limit = 1; /* assume method is not a compound */
3720 special = MagickFalse; /* assume it is NOT a direct modify primitive */
3721 rslt_compose = compose; /* and we are composing multi-kernels as given */
3722 switch( method ) {
3723 case SmoothMorphology: /* 4 primitive compound morphology */
3724 stage_limit = 4;
3725 break;
3726 case OpenMorphology: /* 2 primitive compound morphology */
3727 case OpenIntensityMorphology:
3728 case TopHatMorphology:
3729 case CloseMorphology:
3730 case CloseIntensityMorphology:
3731 case BottomHatMorphology:
3732 case EdgeMorphology:
3733 stage_limit = 2;
3734 break;
3735 case HitAndMissMorphology:
3736 rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3737 /* FALL THUR */
3738 case ThinningMorphology:
3739 case ThickenMorphology:
3740 method_limit = kernel_limit; /* iterate the whole method */
3741 kernel_limit = 1; /* do not do kernel iteration */
3742 break;
3743 case DistanceMorphology:
3744 case VoronoiMorphology:
3745 special = MagickTrue; /* use special direct primative */
3746 break;
3747 default:
3748 break;
3749 }
3750
3751 /* Apply special methods with special requirments
3752 ** For example, single run only, or post-processing requirements
3753 */
3754 if ( special != MagickFalse )
3755 {
3756 rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3757 if (rslt_image == (Image *) NULL)
3758 goto error_cleanup;
3759 if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3760 goto error_cleanup;
3761
3762 changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception);
3763
3764 if (verbose != MagickFalse)
3765 (void) (void) FormatLocaleFile(stderr,
3766 "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3767 CommandOptionToMnemonic(MagickMorphologyOptions, method),
3768 1.0,0.0,1.0, (double) changed);
3769
3770 if ( changed < 0 )
3771 goto error_cleanup;
3772
3773 if ( method == VoronoiMorphology ) {
3774 /* Preserve the alpha channel of input image - but turned it off */
3775 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3776 exception);
3777 (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3778 MagickTrue,0,0,exception);
3779 (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3780 exception);
3781 }
3782 goto exit_cleanup;
3783 }
3784
3785 /* Handle user (caller) specified multi-kernel composition method */
3786 if ( compose != UndefinedCompositeOp )
3787 rslt_compose = compose; /* override default composition for method */
3788 if ( rslt_compose == UndefinedCompositeOp )
3789 rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3790
3791 /* Some methods require a reflected kernel to use with primitives.
3792 * Create the reflected kernel for those methods. */
3793 switch ( method ) {
3794 case CorrelateMorphology:
3795 case CloseMorphology:
3796 case CloseIntensityMorphology:
3797 case BottomHatMorphology:
3798 case SmoothMorphology:
3799 reflected_kernel = CloneKernelInfo(kernel);
3800 if (reflected_kernel == (KernelInfo *) NULL)
3801 goto error_cleanup;
3802 RotateKernelInfo(reflected_kernel,180);
3803 break;
3804 default:
3805 break;
3806 }
3807
3808 /* Loops around more primitive morpholgy methods
3809 ** erose, dilate, open, close, smooth, edge, etc...
3810 */
3811 /* Loop 1: iterate the compound method */
3812 method_loop = 0;
3813 method_changed = 1;
3814 while ( method_loop < method_limit && method_changed > 0 ) {
3815 method_loop++;
3816 method_changed = 0;
3817
3818 /* Loop 2: iterate over each kernel in a multi-kernel list */
3819 norm_kernel = (KernelInfo *) kernel;
3820 this_kernel = (KernelInfo *) kernel;
3821 rflt_kernel = reflected_kernel;
3822
3823 kernel_number = 0;
3824 while ( norm_kernel != NULL ) {
3825
3826 /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3827 stage_loop = 0; /* the compound morphology stage number */
3828 while ( stage_loop < stage_limit ) {
3829 stage_loop++; /* The stage of the compound morphology */
3830
3831 /* Select primitive morphology for this stage of compound method */
3832 this_kernel = norm_kernel; /* default use unreflected kernel */
3833 primitive = method; /* Assume method is a primitive */
3834 switch( method ) {
3835 case ErodeMorphology: /* just erode */
3836 case EdgeInMorphology: /* erode and image difference */
3837 primitive = ErodeMorphology;
3838 break;
3839 case DilateMorphology: /* just dilate */
3840 case EdgeOutMorphology: /* dilate and image difference */
3841 primitive = DilateMorphology;
3842 break;
3843 case OpenMorphology: /* erode then dialate */
3844 case TopHatMorphology: /* open and image difference */
3845 primitive = ErodeMorphology;
3846 if ( stage_loop == 2 )
3847 primitive = DilateMorphology;
3848 break;
3849 case OpenIntensityMorphology:
3850 primitive = ErodeIntensityMorphology;
3851 if ( stage_loop == 2 )
3852 primitive = DilateIntensityMorphology;
3853 break;
3854 case CloseMorphology: /* dilate, then erode */
3855 case BottomHatMorphology: /* close and image difference */
3856 this_kernel = rflt_kernel; /* use the reflected kernel */
3857 primitive = DilateMorphology;
3858 if ( stage_loop == 2 )
3859 primitive = ErodeMorphology;
3860 break;
3861 case CloseIntensityMorphology:
3862 this_kernel = rflt_kernel; /* use the reflected kernel */
3863 primitive = DilateIntensityMorphology;
3864 if ( stage_loop == 2 )
3865 primitive = ErodeIntensityMorphology;
3866 break;
3867 case SmoothMorphology: /* open, close */
3868 switch ( stage_loop ) {
3869 case 1: /* start an open method, which starts with Erode */
3870 primitive = ErodeMorphology;
3871 break;
3872 case 2: /* now Dilate the Erode */
3873 primitive = DilateMorphology;
3874 break;
3875 case 3: /* Reflect kernel a close */
3876 this_kernel = rflt_kernel; /* use the reflected kernel */
3877 primitive = DilateMorphology;
3878 break;
3879 case 4: /* Finish the Close */
3880 this_kernel = rflt_kernel; /* use the reflected kernel */
3881 primitive = ErodeMorphology;
3882 break;
3883 }
3884 break;
3885 case EdgeMorphology: /* dilate and erode difference */
3886 primitive = DilateMorphology;
3887 if ( stage_loop == 2 ) {
3888 save_image = curr_image; /* save the image difference */
3889 curr_image = (Image *) image;
3890 primitive = ErodeMorphology;
3891 }
3892 break;
3893 case CorrelateMorphology:
3894 /* A Correlation is a Convolution with a reflected kernel.
3895 ** However a Convolution is a weighted sum using a reflected
3896 ** kernel. It may seem stange to convert a Correlation into a
3897 ** Convolution as the Correlation is the simplier method, but
3898 ** Convolution is much more commonly used, and it makes sense to
3899 ** implement it directly so as to avoid the need to duplicate the
3900 ** kernel when it is not required (which is typically the
3901 ** default).
3902 */
3903 this_kernel = rflt_kernel; /* use the reflected kernel */
3904 primitive = ConvolveMorphology;
3905 break;
3906 default:
3907 break;
3908 }
3909 assert( this_kernel != (KernelInfo *) NULL );
3910
3911 /* Extra information for debugging compound operations */
3912 if (verbose != MagickFalse) {
3913 if ( stage_limit > 1 )
3914 (void) FormatLocaleString(v_info,MagickPathExtent,"%s:%.20g.%.20g -> ",
3915 CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
3916 method_loop,(double) stage_loop);
3917 else if ( primitive != method )
3918 (void) FormatLocaleString(v_info, MagickPathExtent, "%s:%.20g -> ",
3919 CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
3920 method_loop);
3921 else
3922 v_info[0] = '\0';
3923 }
3924
3925 /* Loop 4: Iterate the kernel with primitive */
3926 kernel_loop = 0;
3927 kernel_changed = 0;
3928 changed = 1;
3929 while ( kernel_loop < kernel_limit && changed > 0 ) {
3930 kernel_loop++; /* the iteration of this kernel */
3931
3932 /* Create a clone as the destination image, if not yet defined */
3933 if ( work_image == (Image *) NULL )
3934 {
3935 work_image=CloneImage(image,0,0,MagickTrue,exception);
3936 if (work_image == (Image *) NULL)
3937 goto error_cleanup;
3938 if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
3939 goto error_cleanup;
3940 }
3941
3942 /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3943 count++;
3944 changed = MorphologyPrimitive(curr_image, work_image, primitive,
3945 this_kernel, bias, exception);
3946 if (verbose != MagickFalse) {
3947 if ( kernel_loop > 1 )
3948 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
3949 (void) (void) FormatLocaleFile(stderr,
3950 "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3951 v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
3952 primitive),(this_kernel == rflt_kernel ) ? "*" : "",
3953 (double) (method_loop+kernel_loop-1),(double) kernel_number,
3954 (double) count,(double) changed);
3955 }
3956 if ( changed < 0 )
3957 goto error_cleanup;
3958 kernel_changed += changed;
3959 method_changed += changed;
3960
3961 /* prepare next loop */
3962 { Image *tmp = work_image; /* swap images for iteration */
3963 work_image = curr_image;
3964 curr_image = tmp;
3965 }
3966 if ( work_image == image )
3967 work_image = (Image *) NULL; /* replace input 'image' */
3968
3969 } /* End Loop 4: Iterate the kernel with primitive */
3970
3971 if (verbose != MagickFalse && kernel_changed != (size_t)changed)
3972 (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
3973 if (verbose != MagickFalse && stage_loop < stage_limit)
3974 (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
3975
3976 #if 0
3977 (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3978 (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
3979 (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
3980 (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
3981 (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
3982 #endif
3983
3984 } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3985
3986 /* Final Post-processing for some Compound Methods
3987 **
3988 ** The removal of any 'Sync' channel flag in the Image Compositon
3989 ** below ensures the methematical compose method is applied in a
3990 ** purely mathematical way, and only to the selected channels.
3991 ** Turn off SVG composition 'alpha blending'.
3992 */
3993 switch( method ) {
3994 case EdgeOutMorphology:
3995 case EdgeInMorphology:
3996 case TopHatMorphology:
3997 case BottomHatMorphology:
3998 if (verbose != MagickFalse)
3999 (void) FormatLocaleFile(stderr,
4000 "\n%s: Difference with original image",CommandOptionToMnemonic(
4001 MagickMorphologyOptions, method) );
4002 (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
4003 MagickTrue,0,0,exception);
4004 break;
4005 case EdgeMorphology:
4006 if (verbose != MagickFalse)
4007 (void) FormatLocaleFile(stderr,
4008 "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
4009 MagickMorphologyOptions, method) );
4010 (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
4011 MagickTrue,0,0,exception);
4012 save_image = DestroyImage(save_image); /* finished with save image */
4013 break;
4014 default:
4015 break;
4016 }
4017
4018 /* multi-kernel handling: re-iterate, or compose results */
4019 if ( kernel->next == (KernelInfo *) NULL )
4020 rslt_image = curr_image; /* just return the resulting image */
4021 else if ( rslt_compose == NoCompositeOp )
4022 { if (verbose != MagickFalse) {
4023 if ( this_kernel->next != (KernelInfo *) NULL )
4024 (void) FormatLocaleFile(stderr, " (re-iterate)");
4025 else
4026 (void) FormatLocaleFile(stderr, " (done)");
4027 }
4028 rslt_image = curr_image; /* return result, and re-iterate */
4029 }
4030 else if ( rslt_image == (Image *) NULL)
4031 { if (verbose != MagickFalse)
4032 (void) FormatLocaleFile(stderr, " (save for compose)");
4033 rslt_image = curr_image;
4034 curr_image = (Image *) image; /* continue with original image */
4035 }
4036 else
4037 { /* Add the new 'current' result to the composition
4038 **
4039 ** The removal of any 'Sync' channel flag in the Image Compositon
4040 ** below ensures the methematical compose method is applied in a
4041 ** purely mathematical way, and only to the selected channels.
4042 ** IE: Turn off SVG composition 'alpha blending'.
4043 */
4044 if (verbose != MagickFalse)
4045 (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4046 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
4047 (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4048 0,0,exception);
4049 curr_image = DestroyImage(curr_image);
4050 curr_image = (Image *) image; /* continue with original image */
4051 }
4052 if (verbose != MagickFalse)
4053 (void) FormatLocaleFile(stderr, "\n");
4054
4055 /* loop to the next kernel in a multi-kernel list */
4056 norm_kernel = norm_kernel->next;
4057 if ( rflt_kernel != (KernelInfo *) NULL )
4058 rflt_kernel = rflt_kernel->next;
4059 kernel_number++;
4060 } /* End Loop 2: Loop over each kernel */
4061
4062 } /* End Loop 1: compound method interation */
4063
4064 goto exit_cleanup;
4065
4066 /* Yes goto's are bad, but it makes cleanup lot more efficient */
4067 error_cleanup:
4068 if ( curr_image == rslt_image )
4069 curr_image = (Image *) NULL;
4070 if ( rslt_image != (Image *) NULL )
4071 rslt_image = DestroyImage(rslt_image);
4072 exit_cleanup:
4073 if ( curr_image == rslt_image || curr_image == image )
4074 curr_image = (Image *) NULL;
4075 if ( curr_image != (Image *) NULL )
4076 curr_image = DestroyImage(curr_image);
4077 if ( work_image != (Image *) NULL )
4078 work_image = DestroyImage(work_image);
4079 if ( save_image != (Image *) NULL )
4080 save_image = DestroyImage(save_image);
4081 if ( reflected_kernel != (KernelInfo *) NULL )
4082 reflected_kernel = DestroyKernelInfo(reflected_kernel);
4083 return(rslt_image);
4084 }
4085
4086
4087 /*
4088 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4089 % %
4090 % %
4091 % %
4092 % M o r p h o l o g y I m a g e %
4093 % %
4094 % %
4095 % %
4096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4097 %
4098 % MorphologyImage() applies a user supplied kernel to the image according to
4099 % the given mophology method.
4100 %
4101 % This function applies any and all user defined settings before calling
4102 % the above internal function MorphologyApply().
4103 %
4104 % User defined settings include...
4105 % * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4106 % * Kernel Scale/normalize settings ("-define convolve:scale=??")
4107 % This can also includes the addition of a scaled unity kernel.
4108 % * Show Kernel being applied ("-define morphology:showKernel=1")
4109 %
4110 % Other operators that do not want user supplied options interfering,
4111 % especially "convolve:bias" and "morphology:showKernel" should use
4112 % MorphologyApply() directly.
4113 %
4114 % The format of the MorphologyImage method is:
4115 %
4116 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4117 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4118 %
4119 % A description of each parameter follows:
4120 %
4121 % o image: the image.
4122 %
4123 % o method: the morphology method to be applied.
4124 %
4125 % o iterations: apply the operation this many times (or no change).
4126 % A value of -1 means loop until no change found.
4127 % How this is applied may depend on the morphology method.
4128 % Typically this is a value of 1.
4129 %
4130 % o kernel: An array of double representing the morphology kernel.
4131 % Warning: kernel may be normalized for the Convolve method.
4132 %
4133 % o exception: return any errors or warnings in this structure.
4134 %
4135 */
MorphologyImage(const Image * image,const MorphologyMethod method,const ssize_t iterations,const KernelInfo * kernel,ExceptionInfo * exception)4136 MagickExport Image *MorphologyImage(const Image *image,
4137 const MorphologyMethod method,const ssize_t iterations,
4138 const KernelInfo *kernel,ExceptionInfo *exception)
4139 {
4140 const char
4141 *artifact;
4142
4143 CompositeOperator
4144 compose;
4145
4146 double
4147 bias;
4148
4149 Image
4150 *morphology_image;
4151
4152 KernelInfo
4153 *curr_kernel;
4154
4155 curr_kernel = (KernelInfo *) kernel;
4156 bias=0.0;
4157 compose = UndefinedCompositeOp; /* use default for method */
4158
4159 /* Apply Convolve/Correlate Normalization and Scaling Factors.
4160 * This is done BEFORE the ShowKernelInfo() function is called so that
4161 * users can see the results of the 'option:convolve:scale' option.
4162 */
4163 if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4164 /* Get the bias value as it will be needed */
4165 artifact = GetImageArtifact(image,"convolve:bias");
4166 if ( artifact != (const char *) NULL) {
4167 if (IsGeometry(artifact) == MagickFalse)
4168 (void) ThrowMagickException(exception,GetMagickModule(),
4169 OptionWarning,"InvalidSetting","'%s' '%s'",
4170 "convolve:bias",artifact);
4171 else
4172 bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4173 }
4174
4175 /* Scale kernel according to user wishes */
4176 artifact = GetImageArtifact(image,"convolve:scale");
4177 if ( artifact != (const char *) NULL ) {
4178 if (IsGeometry(artifact) == MagickFalse)
4179 (void) ThrowMagickException(exception,GetMagickModule(),
4180 OptionWarning,"InvalidSetting","'%s' '%s'",
4181 "convolve:scale",artifact);
4182 else {
4183 if ( curr_kernel == kernel )
4184 curr_kernel = CloneKernelInfo(kernel);
4185 if (curr_kernel == (KernelInfo *) NULL)
4186 return((Image *) NULL);
4187 ScaleGeometryKernelInfo(curr_kernel, artifact);
4188 }
4189 }
4190 }
4191
4192 /* display the (normalized) kernel via stderr */
4193 artifact=GetImageArtifact(image,"morphology:showKernel");
4194 if (IsStringTrue(artifact) != MagickFalse)
4195 ShowKernelInfo(curr_kernel);
4196
4197 /* Override the default handling of multi-kernel morphology results
4198 * If 'Undefined' use the default method
4199 * If 'None' (default for 'Convolve') re-iterate previous result
4200 * Otherwise merge resulting images using compose method given.
4201 * Default for 'HitAndMiss' is 'Lighten'.
4202 */
4203 {
4204 ssize_t
4205 parse;
4206
4207 artifact = GetImageArtifact(image,"morphology:compose");
4208 if ( artifact != (const char *) NULL) {
4209 parse=ParseCommandOption(MagickComposeOptions,
4210 MagickFalse,artifact);
4211 if ( parse < 0 )
4212 (void) ThrowMagickException(exception,GetMagickModule(),
4213 OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4214 "morphology:compose",artifact);
4215 else
4216 compose=(CompositeOperator)parse;
4217 }
4218 }
4219 /* Apply the Morphology */
4220 morphology_image = MorphologyApply(image,method,iterations,
4221 curr_kernel,compose,bias,exception);
4222
4223 /* Cleanup and Exit */
4224 if ( curr_kernel != kernel )
4225 curr_kernel=DestroyKernelInfo(curr_kernel);
4226 return(morphology_image);
4227 }
4228
4229 /*
4230 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4231 % %
4232 % %
4233 % %
4234 + R o t a t e K e r n e l I n f o %
4235 % %
4236 % %
4237 % %
4238 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4239 %
4240 % RotateKernelInfo() rotates the kernel by the angle given.
4241 %
4242 % Currently it is restricted to 90 degree angles, of either 1D kernels
4243 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4244 % It will ignore usless rotations for specific 'named' built-in kernels.
4245 %
4246 % The format of the RotateKernelInfo method is:
4247 %
4248 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4249 %
4250 % A description of each parameter follows:
4251 %
4252 % o kernel: the Morphology/Convolution kernel
4253 %
4254 % o angle: angle to rotate in degrees
4255 %
4256 % This function is currently internal to this module only, but can be exported
4257 % to other modules if needed.
4258 */
RotateKernelInfo(KernelInfo * kernel,double angle)4259 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4260 {
4261 /* angle the lower kernels first */
4262 if ( kernel->next != (KernelInfo *) NULL)
4263 RotateKernelInfo(kernel->next, angle);
4264
4265 /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4266 **
4267 ** TODO: expand beyond simple 90 degree rotates, flips and flops
4268 */
4269
4270 /* Modulus the angle */
4271 angle = fmod(angle, 360.0);
4272 if ( angle < 0 )
4273 angle += 360.0;
4274
4275 if ( 337.5 < angle || angle <= 22.5 )
4276 return; /* Near zero angle - no change! - At least not at this time */
4277
4278 /* Handle special cases */
4279 switch (kernel->type) {
4280 /* These built-in kernels are cylindrical kernels, rotating is useless */
4281 case GaussianKernel:
4282 case DoGKernel:
4283 case LoGKernel:
4284 case DiskKernel:
4285 case PeaksKernel:
4286 case LaplacianKernel:
4287 case ChebyshevKernel:
4288 case ManhattanKernel:
4289 case EuclideanKernel:
4290 return;
4291
4292 /* These may be rotatable at non-90 angles in the future */
4293 /* but simply rotating them in multiples of 90 degrees is useless */
4294 case SquareKernel:
4295 case DiamondKernel:
4296 case PlusKernel:
4297 case CrossKernel:
4298 return;
4299
4300 /* These only allows a +/-90 degree rotation (by transpose) */
4301 /* A 180 degree rotation is useless */
4302 case BlurKernel:
4303 if ( 135.0 < angle && angle <= 225.0 )
4304 return;
4305 if ( 225.0 < angle && angle <= 315.0 )
4306 angle -= 180;
4307 break;
4308
4309 default:
4310 break;
4311 }
4312 /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4313 if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4314 {
4315 if ( kernel->width == 3 && kernel->height == 3 )
4316 { /* Rotate a 3x3 square by 45 degree angle */
4317 double t = kernel->values[0];
4318 kernel->values[0] = kernel->values[3];
4319 kernel->values[3] = kernel->values[6];
4320 kernel->values[6] = kernel->values[7];
4321 kernel->values[7] = kernel->values[8];
4322 kernel->values[8] = kernel->values[5];
4323 kernel->values[5] = kernel->values[2];
4324 kernel->values[2] = kernel->values[1];
4325 kernel->values[1] = t;
4326 /* rotate non-centered origin */
4327 if ( kernel->x != 1 || kernel->y != 1 ) {
4328 ssize_t x,y;
4329 x = (ssize_t) kernel->x-1;
4330 y = (ssize_t) kernel->y-1;
4331 if ( x == y ) x = 0;
4332 else if ( x == 0 ) x = -y;
4333 else if ( x == -y ) y = 0;
4334 else if ( y == 0 ) y = x;
4335 kernel->x = (ssize_t) x+1;
4336 kernel->y = (ssize_t) y+1;
4337 }
4338 angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4339 kernel->angle = fmod(kernel->angle+45.0, 360.0);
4340 }
4341 else
4342 perror("Unable to rotate non-3x3 kernel by 45 degrees");
4343 }
4344 if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4345 {
4346 if ( kernel->width == 1 || kernel->height == 1 )
4347 { /* Do a transpose of a 1 dimensional kernel,
4348 ** which results in a fast 90 degree rotation of some type.
4349 */
4350 ssize_t
4351 t;
4352 t = (ssize_t) kernel->width;
4353 kernel->width = kernel->height;
4354 kernel->height = (size_t) t;
4355 t = kernel->x;
4356 kernel->x = kernel->y;
4357 kernel->y = t;
4358 if ( kernel->width == 1 ) {
4359 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4360 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4361 } else {
4362 angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4363 kernel->angle = fmod(kernel->angle+270.0, 360.0);
4364 }
4365 }
4366 else if ( kernel->width == kernel->height )
4367 { /* Rotate a square array of values by 90 degrees */
4368 { ssize_t
4369 i,j,x,y;
4370
4371 MagickRealType
4372 *k,t;
4373
4374 k=kernel->values;
4375 for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--)
4376 for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--)
4377 { t = k[i+j*kernel->width];
4378 k[i+j*kernel->width] = k[j+x*kernel->width];
4379 k[j+x*kernel->width] = k[x+y*kernel->width];
4380 k[x+y*kernel->width] = k[y+i*kernel->width];
4381 k[y+i*kernel->width] = t;
4382 }
4383 }
4384 /* rotate the origin - relative to center of array */
4385 { ssize_t x,y;
4386 x = (ssize_t) (kernel->x*2-kernel->width+1);
4387 y = (ssize_t) (kernel->y*2-kernel->height+1);
4388 kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4389 kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4390 }
4391 angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4392 kernel->angle = fmod(kernel->angle+90.0, 360.0);
4393 }
4394 else
4395 perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4396 }
4397 if ( 135.0 < angle && angle <= 225.0 )
4398 {
4399 /* For a 180 degree rotation - also know as a reflection
4400 * This is actually a very very common operation!
4401 * Basically all that is needed is a reversal of the kernel data!
4402 * And a reflection of the origon
4403 */
4404 MagickRealType
4405 t;
4406
4407 MagickRealType
4408 *k;
4409
4410 ssize_t
4411 i,
4412 j;
4413
4414 k=kernel->values;
4415 j=(ssize_t) (kernel->width*kernel->height-1);
4416 for (i=0; i < j; i++, j--)
4417 t=k[i], k[i]=k[j], k[j]=t;
4418
4419 kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4420 kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4421 angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4422 kernel->angle = fmod(kernel->angle+180.0, 360.0);
4423 }
4424 /* At this point angle should at least between -45 (315) and +45 degrees
4425 * In the future some form of non-orthogonal angled rotates could be
4426 * performed here, posibily with a linear kernel restriction.
4427 */
4428
4429 return;
4430 }
4431
4432 /*
4433 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4434 % %
4435 % %
4436 % %
4437 % S c a l e G e o m e t r y K e r n e l I n f o %
4438 % %
4439 % %
4440 % %
4441 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4442 %
4443 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4444 % provided as a "-set option:convolve:scale {geometry}" user setting,
4445 % and modifies the kernel according to the parsed arguments of that setting.
4446 %
4447 % The first argument (and any normalization flags) are passed to
4448 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4449 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4450 % into the scaled/normalized kernel.
4451 %
4452 % The format of the ScaleGeometryKernelInfo method is:
4453 %
4454 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4455 % const double scaling_factor,const MagickStatusType normalize_flags)
4456 %
4457 % A description of each parameter follows:
4458 %
4459 % o kernel: the Morphology/Convolution kernel to modify
4460 %
4461 % o geometry:
4462 % The geometry string to parse, typically from the user provided
4463 % "-set option:convolve:scale {geometry}" setting.
4464 %
4465 */
ScaleGeometryKernelInfo(KernelInfo * kernel,const char * geometry)4466 MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
4467 const char *geometry)
4468 {
4469 MagickStatusType
4470 flags;
4471
4472 GeometryInfo
4473 args;
4474
4475 SetGeometryInfo(&args);
4476 flags = ParseGeometry(geometry, &args);
4477
4478 #if 0
4479 /* For Debugging Geometry Input */
4480 (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4481 flags, args.rho, args.sigma, args.xi, args.psi );
4482 #endif
4483
4484 if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4485 args.rho *= 0.01, args.sigma *= 0.01;
4486
4487 if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4488 args.rho = 1.0;
4489 if ( (flags & SigmaValue) == 0 )
4490 args.sigma = 0.0;
4491
4492 /* Scale/Normalize the input kernel */
4493 ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4494
4495 /* Add Unity Kernel, for blending with original */
4496 if ( (flags & SigmaValue) != 0 )
4497 UnityAddKernelInfo(kernel, args.sigma);
4498
4499 return;
4500 }
4501 /*
4502 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4503 % %
4504 % %
4505 % %
4506 % S c a l e K e r n e l I n f o %
4507 % %
4508 % %
4509 % %
4510 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4511 %
4512 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4513 % without normalization of the sum of the kernel values (as per given flags).
4514 %
4515 % By default (no flags given) the values within the kernel is scaled
4516 % directly using given scaling factor without change.
4517 %
4518 % If either of the two 'normalize_flags' are given the kernel will first be
4519 % normalized and then further scaled by the scaling factor value given.
4520 %
4521 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4522 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4523 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4524 % non-HDRI versions of IM this may cause images to have any negative results
4525 % clipped, unless some 'bias' is used.
4526 %
4527 % More specifically. Kernels which only contain positive values (such as a
4528 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4529 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4530 %
4531 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4532 % the kernel will be scaled by the absolute of the sum of kernel values, so
4533 % that it will generally fall within the +/- 1.0 range.
4534 %
4535 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4536 % will be scaled by just the sum of the postive values, so that its output
4537 % range will again fall into the +/- 1.0 range.
4538 %
4539 % For special kernels designed for locating shapes using 'Correlate', (often
4540 % only containing +1 and -1 values, representing foreground/brackground
4541 % matching) a special normalization method is provided to scale the positive
4542 % values separately to those of the negative values, so the kernel will be
4543 % forced to become a zero-sum kernel better suited to such searches.
4544 %
4545 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4546 % attributes within the kernel structure have been correctly set during the
4547 % kernels creation.
4548 %
4549 % NOTE: The values used for 'normalize_flags' have been selected specifically
4550 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4551 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4552 %
4553 % The format of the ScaleKernelInfo method is:
4554 %
4555 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4556 % const MagickStatusType normalize_flags )
4557 %
4558 % A description of each parameter follows:
4559 %
4560 % o kernel: the Morphology/Convolution kernel
4561 %
4562 % o scaling_factor:
4563 % multiply all values (after normalization) by this factor if not
4564 % zero. If the kernel is normalized regardless of any flags.
4565 %
4566 % o normalize_flags:
4567 % GeometryFlags defining normalization method to use.
4568 % specifically: NormalizeValue, CorrelateNormalizeValue,
4569 % and/or PercentValue
4570 %
4571 */
ScaleKernelInfo(KernelInfo * kernel,const double scaling_factor,const GeometryFlags normalize_flags)4572 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
4573 const double scaling_factor,const GeometryFlags normalize_flags)
4574 {
4575 double
4576 pos_scale,
4577 neg_scale;
4578
4579 ssize_t
4580 i;
4581
4582 /* do the other kernels in a multi-kernel list first */
4583 if ( kernel->next != (KernelInfo *) NULL)
4584 ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4585
4586 /* Normalization of Kernel */
4587 pos_scale = 1.0;
4588 if ( (normalize_flags&NormalizeValue) != 0 ) {
4589 if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4590 /* non-zero-summing kernel (generally positive) */
4591 pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4592 else
4593 /* zero-summing kernel */
4594 pos_scale = kernel->positive_range;
4595 }
4596 /* Force kernel into a normalized zero-summing kernel */
4597 if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4598 pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4599 ? kernel->positive_range : 1.0;
4600 neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4601 ? -kernel->negative_range : 1.0;
4602 }
4603 else
4604 neg_scale = pos_scale;
4605
4606 /* finialize scaling_factor for positive and negative components */
4607 pos_scale = scaling_factor/pos_scale;
4608 neg_scale = scaling_factor/neg_scale;
4609
4610 for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4611 if (!IsNaN(kernel->values[i]))
4612 kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4613
4614 /* convolution output range */
4615 kernel->positive_range *= pos_scale;
4616 kernel->negative_range *= neg_scale;
4617 /* maximum and minimum values in kernel */
4618 kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4619 kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4620
4621 /* swap kernel settings if user's scaling factor is negative */
4622 if ( scaling_factor < MagickEpsilon ) {
4623 double t;
4624 t = kernel->positive_range;
4625 kernel->positive_range = kernel->negative_range;
4626 kernel->negative_range = t;
4627 t = kernel->maximum;
4628 kernel->maximum = kernel->minimum;
4629 kernel->minimum = 1;
4630 }
4631
4632 return;
4633 }
4634
4635 /*
4636 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4637 % %
4638 % %
4639 % %
4640 % S h o w K e r n e l I n f o %
4641 % %
4642 % %
4643 % %
4644 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4645 %
4646 % ShowKernelInfo() outputs the details of the given kernel defination to
4647 % standard error, generally due to a users 'morphology:showKernel' option
4648 % request.
4649 %
4650 % The format of the ShowKernel method is:
4651 %
4652 % void ShowKernelInfo(const KernelInfo *kernel)
4653 %
4654 % A description of each parameter follows:
4655 %
4656 % o kernel: the Morphology/Convolution kernel
4657 %
4658 */
ShowKernelInfo(const KernelInfo * kernel)4659 MagickPrivate void ShowKernelInfo(const KernelInfo *kernel)
4660 {
4661 const KernelInfo
4662 *k;
4663
4664 size_t
4665 c, i, u, v;
4666
4667 for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4668
4669 (void) FormatLocaleFile(stderr, "Kernel");
4670 if ( kernel->next != (KernelInfo *) NULL )
4671 (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4672 (void) FormatLocaleFile(stderr, " \"%s",
4673 CommandOptionToMnemonic(MagickKernelOptions, k->type) );
4674 if ( fabs(k->angle) >= MagickEpsilon )
4675 (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4676 (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4677 k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4678 (void) FormatLocaleFile(stderr,
4679 " with values from %.*lg to %.*lg\n",
4680 GetMagickPrecision(), k->minimum,
4681 GetMagickPrecision(), k->maximum);
4682 (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4683 GetMagickPrecision(), k->negative_range,
4684 GetMagickPrecision(), k->positive_range);
4685 if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4686 (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4687 else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4688 (void) FormatLocaleFile(stderr, " (Normalized)\n");
4689 else
4690 (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4691 GetMagickPrecision(), k->positive_range+k->negative_range);
4692 for (i=v=0; v < k->height; v++) {
4693 (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4694 for (u=0; u < k->width; u++, i++)
4695 if (IsNaN(k->values[i]))
4696 (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4697 else
4698 (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4699 GetMagickPrecision(), (double) k->values[i]);
4700 (void) FormatLocaleFile(stderr,"\n");
4701 }
4702 }
4703 }
4704
4705 /*
4706 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4707 % %
4708 % %
4709 % %
4710 % U n i t y A d d K e r n a l I n f o %
4711 % %
4712 % %
4713 % %
4714 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4715 %
4716 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4717 % to the given pre-scaled and normalized Kernel. This in effect adds that
4718 % amount of the original image into the resulting convolution kernel. This
4719 % value is usually provided by the user as a percentage value in the
4720 % 'convolve:scale' setting.
4721 %
4722 % The resulting effect is to convert the defined kernels into blended
4723 % soft-blurs, unsharp kernels or into sharpening kernels.
4724 %
4725 % The format of the UnityAdditionKernelInfo method is:
4726 %
4727 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4728 %
4729 % A description of each parameter follows:
4730 %
4731 % o kernel: the Morphology/Convolution kernel
4732 %
4733 % o scale:
4734 % scaling factor for the unity kernel to be added to
4735 % the given kernel.
4736 %
4737 */
UnityAddKernelInfo(KernelInfo * kernel,const double scale)4738 MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
4739 const double scale)
4740 {
4741 /* do the other kernels in a multi-kernel list first */
4742 if ( kernel->next != (KernelInfo *) NULL)
4743 UnityAddKernelInfo(kernel->next, scale);
4744
4745 /* Add the scaled unity kernel to the existing kernel */
4746 kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4747 CalcKernelMetaData(kernel); /* recalculate the meta-data */
4748
4749 return;
4750 }
4751
4752 /*
4753 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4754 % %
4755 % %
4756 % %
4757 % Z e r o K e r n e l N a n s %
4758 % %
4759 % %
4760 % %
4761 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4762 %
4763 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4764 % the kernel with a zero value. This is typically done when the kernel will
4765 % be used in special hardware (GPU) convolution processors, to simply
4766 % matters.
4767 %
4768 % The format of the ZeroKernelNans method is:
4769 %
4770 % void ZeroKernelNans (KernelInfo *kernel)
4771 %
4772 % A description of each parameter follows:
4773 %
4774 % o kernel: the Morphology/Convolution kernel
4775 %
4776 */
ZeroKernelNans(KernelInfo * kernel)4777 MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
4778 {
4779 size_t
4780 i;
4781
4782 /* do the other kernels in a multi-kernel list first */
4783 if (kernel->next != (KernelInfo *) NULL)
4784 ZeroKernelNans(kernel->next);
4785
4786 for (i=0; i < (kernel->width*kernel->height); i++)
4787 if (IsNaN(kernel->values[i]))
4788 kernel->values[i]=0.0;
4789
4790 return;
4791 }
4792