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