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