MagickCore  7.0.7
Convert, Edit, Or Compose Bitmap Images
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-2018 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://www.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  register 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 *) AcquireQuantumMemory(1,sizeof(*kernel));
239  if (kernel == (KernelInfo *) NULL)
240  return(kernel);
241  (void) ResetMagickMemory(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  GetNextToken(p,&p,MagickPathExtent,token);
302  if (*token == ',')
303  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  GetNextToken(p,&p,MagickPathExtent,token);
324  if (*token == ',')
325  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  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  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  register ssize_t
958  i;
959 
960  register ssize_t
961  u,
962  v;
963 
964  double
965  nan = sqrt((double)-1.0); /* Special Value : Not A Number */
966 
967  /* Generate a new empty kernel if needed */
968  kernel=(KernelInfo *) NULL;
969  switch(type) {
970  case UndefinedKernel: /* These should not call this function */
971  case UserDefinedKernel:
972  assert("Should not call this function" != (char *) NULL);
973  break;
974  case LaplacianKernel: /* Named Descrete Convolution Kernels */
975  case SobelKernel: /* these are defined using other kernels */
976  case RobertsKernel:
977  case PrewittKernel:
978  case CompassKernel:
979  case KirschKernel:
980  case FreiChenKernel:
981  case EdgesKernel: /* Hit and Miss kernels */
982  case CornersKernel:
983  case DiagonalsKernel:
984  case LineEndsKernel:
985  case LineJunctionsKernel:
986  case RidgesKernel:
987  case ConvexHullKernel:
988  case SkeletonKernel:
989  case ThinSEKernel:
990  break; /* A pre-generated kernel is not needed */
991 #if 0
992  /* set to 1 to do a compile-time check that we haven't missed anything */
993  case UnityKernel:
994  case GaussianKernel:
995  case DoGKernel:
996  case LoGKernel:
997  case BlurKernel:
998  case CometKernel:
999  case BinomialKernel:
1000  case DiamondKernel:
1001  case SquareKernel:
1002  case RectangleKernel:
1003  case OctagonKernel:
1004  case DiskKernel:
1005  case PlusKernel:
1006  case CrossKernel:
1007  case RingKernel:
1008  case PeaksKernel:
1009  case ChebyshevKernel:
1010  case ManhattanKernel:
1011  case OctangonalKernel:
1012  case EuclideanKernel:
1013 #else
1014  default:
1015 #endif
1016  /* Generate the base Kernel Structure */
1017  kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
1018  if (kernel == (KernelInfo *) NULL)
1019  return(kernel);
1020  (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
1021  kernel->minimum = kernel->maximum = kernel->angle = 0.0;
1022  kernel->negative_range = kernel->positive_range = 0.0;
1023  kernel->type = type;
1024  kernel->next = (KernelInfo *) NULL;
1026  break;
1027  }
1028 
1029  switch(type) {
1030  /*
1031  Convolution Kernels
1032  */
1033  case UnityKernel:
1034  {
1035  kernel->height = kernel->width = (size_t) 1;
1036  kernel->x = kernel->y = (ssize_t) 0;
1038  AcquireAlignedMemory(1,sizeof(*kernel->values)));
1039  if (kernel->values == (MagickRealType *) NULL)
1040  return(DestroyKernelInfo(kernel));
1041  kernel->maximum = kernel->values[0] = args->rho;
1042  break;
1043  }
1044  break;
1045  case GaussianKernel:
1046  case DoGKernel:
1047  case LoGKernel:
1048  { double
1049  sigma = fabs(args->sigma),
1050  sigma2 = fabs(args->xi),
1051  A, B, R;
1052 
1053  if ( args->rho >= 1.0 )
1054  kernel->width = (size_t)args->rho*2+1;
1055  else if ( (type != DoGKernel) || (sigma >= sigma2) )
1056  kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
1057  else
1058  kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
1059  kernel->height = kernel->width;
1060  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1062  AcquireAlignedMemory(kernel->width,kernel->height*
1063  sizeof(*kernel->values)));
1064  if (kernel->values == (MagickRealType *) NULL)
1065  return(DestroyKernelInfo(kernel));
1066 
1067  /* WARNING: The following generates a 'sampled gaussian' kernel.
1068  * What we really want is a 'discrete gaussian' kernel.
1069  *
1070  * How to do this is I don't know, but appears to be basied on the
1071  * Error Function 'erf()' (intergral of a gaussian)
1072  */
1073 
1074  if ( type == GaussianKernel || type == DoGKernel )
1075  { /* Calculate a Gaussian, OR positive half of a DoG */
1076  if ( sigma > MagickEpsilon )
1077  { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1078  B = (double) (1.0/(Magick2PI*sigma*sigma));
1079  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1080  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1081  kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
1082  }
1083  else /* limiting case - a unity (normalized Dirac) kernel */
1084  { (void) ResetMagickMemory(kernel->values,0, (size_t)
1085  kernel->width*kernel->height*sizeof(*kernel->values));
1086  kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1087  }
1088  }
1089 
1090  if ( type == DoGKernel )
1091  { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
1092  if ( sigma2 > MagickEpsilon )
1093  { sigma = sigma2; /* simplify loop expressions */
1094  A = 1.0/(2.0*sigma*sigma);
1095  B = (double) (1.0/(Magick2PI*sigma*sigma));
1096  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1097  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1098  kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
1099  }
1100  else /* limiting case - a unity (normalized Dirac) kernel */
1101  kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
1102  }
1103 
1104  if ( type == LoGKernel )
1105  { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
1106  if ( sigma > MagickEpsilon )
1107  { A = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1108  B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
1109  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1110  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1111  { R = ((double)(u*u+v*v))*A;
1112  kernel->values[i] = (1-R)*exp(-R)*B;
1113  }
1114  }
1115  else /* special case - generate a unity kernel */
1116  { (void) ResetMagickMemory(kernel->values,0, (size_t)
1117  kernel->width*kernel->height*sizeof(*kernel->values));
1118  kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1119  }
1120  }
1121 
1122  /* Note the above kernels may have been 'clipped' by a user defined
1123  ** radius, producing a smaller (darker) kernel. Also for very small
1124  ** sigma's (> 0.1) the central value becomes larger than one, and thus
1125  ** producing a very bright kernel.
1126  **
1127  ** Normalization will still be needed.
1128  */
1129 
1130  /* Normalize the 2D Gaussian Kernel
1131  **
1132  ** NB: a CorrelateNormalize performs a normal Normalize if
1133  ** there are no negative values.
1134  */
1135  CalcKernelMetaData(kernel); /* the other kernel meta-data */
1137 
1138  break;
1139  }
1140  case BlurKernel:
1141  { double
1142  sigma = fabs(args->sigma),
1143  alpha, beta;
1144 
1145  if ( args->rho >= 1.0 )
1146  kernel->width = (size_t)args->rho*2+1;
1147  else
1148  kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
1149  kernel->height = 1;
1150  kernel->x = (ssize_t) (kernel->width-1)/2;
1151  kernel->y = 0;
1152  kernel->negative_range = kernel->positive_range = 0.0;
1154  AcquireAlignedMemory(kernel->width,kernel->height*
1155  sizeof(*kernel->values)));
1156  if (kernel->values == (MagickRealType *) NULL)
1157  return(DestroyKernelInfo(kernel));
1158 
1159 #if 1
1160 #define KernelRank 3
1161  /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
1162  ** It generates a gaussian 3 times the width, and compresses it into
1163  ** the expected range. This produces a closer normalization of the
1164  ** resulting kernel, especially for very low sigma values.
1165  ** As such while wierd it is prefered.
1166  **
1167  ** I am told this method originally came from Photoshop.
1168  **
1169  ** A properly normalized curve is generated (apart from edge clipping)
1170  ** even though we later normalize the result (for edge clipping)
1171  ** to allow the correct generation of a "Difference of Blurs".
1172  */
1173 
1174  /* initialize */
1175  v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
1176  (void) ResetMagickMemory(kernel->values,0, (size_t)
1177  kernel->width*kernel->height*sizeof(*kernel->values));
1178  /* Calculate a Positive 1D Gaussian */
1179  if ( sigma > MagickEpsilon )
1180  { sigma *= KernelRank; /* simplify loop expressions */
1181  alpha = 1.0/(2.0*sigma*sigma);
1182  beta= (double) (1.0/(MagickSQ2PI*sigma ));
1183  for ( u=-v; u <= v; u++) {
1184  kernel->values[(u+v)/KernelRank] +=
1185  exp(-((double)(u*u))*alpha)*beta;
1186  }
1187  }
1188  else /* special case - generate a unity kernel */
1189  kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1190 #else
1191  /* Direct calculation without curve averaging
1192  This is equivelent to a KernelRank of 1 */
1193 
1194  /* Calculate a Positive Gaussian */
1195  if ( sigma > MagickEpsilon )
1196  { alpha = 1.0/(2.0*sigma*sigma); /* simplify loop expressions */
1197  beta = 1.0/(MagickSQ2PI*sigma);
1198  for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1199  kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
1200  }
1201  else /* special case - generate a unity kernel */
1202  { (void) ResetMagickMemory(kernel->values,0, (size_t)
1203  kernel->width*kernel->height*sizeof(*kernel->values));
1204  kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1205  }
1206 #endif
1207  /* Note the above kernel may have been 'clipped' by a user defined
1208  ** radius, producing a smaller (darker) kernel. Also for very small
1209  ** sigma's (> 0.1) the central value becomes larger than one, as a
1210  ** result of not generating a actual 'discrete' kernel, and thus
1211  ** producing a very bright 'impulse'.
1212  **
1213  ** Becuase of these two factors Normalization is required!
1214  */
1215 
1216  /* Normalize the 1D Gaussian Kernel
1217  **
1218  ** NB: a CorrelateNormalize performs a normal Normalize if
1219  ** there are no negative values.
1220  */
1221  CalcKernelMetaData(kernel); /* the other kernel meta-data */
1223 
1224  /* rotate the 1D kernel by given angle */
1225  RotateKernelInfo(kernel, args->xi );
1226  break;
1227  }
1228  case CometKernel:
1229  { double
1230  sigma = fabs(args->sigma),
1231  A;
1232 
1233  if ( args->rho < 1.0 )
1234  kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
1235  else
1236  kernel->width = (size_t)args->rho;
1237  kernel->x = kernel->y = 0;
1238  kernel->height = 1;
1239  kernel->negative_range = kernel->positive_range = 0.0;
1241  AcquireAlignedMemory(kernel->width,kernel->height*
1242  sizeof(*kernel->values)));
1243  if (kernel->values == (MagickRealType *) NULL)
1244  return(DestroyKernelInfo(kernel));
1245 
1246  /* A comet blur is half a 1D gaussian curve, so that the object is
1247  ** blurred in one direction only. This may not be quite the right
1248  ** curve to use so may change in the future. The function must be
1249  ** normalised after generation, which also resolves any clipping.
1250  **
1251  ** As we are normalizing and not subtracting gaussians,
1252  ** there is no need for a divisor in the gaussian formula
1253  **
1254  ** It is less comples
1255  */
1256  if ( sigma > MagickEpsilon )
1257  {
1258 #if 1
1259 #define KernelRank 3
1260  v = (ssize_t) kernel->width*KernelRank; /* start/end points */
1261  (void) ResetMagickMemory(kernel->values,0, (size_t)
1262  kernel->width*sizeof(*kernel->values));
1263  sigma *= KernelRank; /* simplify the loop expression */
1264  A = 1.0/(2.0*sigma*sigma);
1265  /* B = 1.0/(MagickSQ2PI*sigma); */
1266  for ( u=0; u < v; u++) {
1267  kernel->values[u/KernelRank] +=
1268  exp(-((double)(u*u))*A);
1269  /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1270  }
1271  for (i=0; i < (ssize_t) kernel->width; i++)
1272  kernel->positive_range += kernel->values[i];
1273 #else
1274  A = 1.0/(2.0*sigma*sigma); /* simplify the loop expression */
1275  /* B = 1.0/(MagickSQ2PI*sigma); */
1276  for ( i=0; i < (ssize_t) kernel->width; i++)
1277  kernel->positive_range +=
1278  kernel->values[i] = exp(-((double)(i*i))*A);
1279  /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
1280 #endif
1281  }
1282  else /* special case - generate a unity kernel */
1283  { (void) ResetMagickMemory(kernel->values,0, (size_t)
1284  kernel->width*kernel->height*sizeof(*kernel->values));
1285  kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1286  kernel->positive_range = 1.0;
1287  }
1288 
1289  kernel->minimum = 0.0;
1290  kernel->maximum = kernel->values[0];
1291  kernel->negative_range = 0.0;
1292 
1293  ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
1294  RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
1295  break;
1296  }
1297  case BinomialKernel:
1298  {
1299  size_t
1300  order_f;
1301 
1302  if (args->rho < 1.0)
1303  kernel->width = kernel->height = 3; /* default radius = 1 */
1304  else
1305  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1306  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1307 
1308  order_f = fact(kernel->width-1);
1309 
1311  AcquireAlignedMemory(kernel->width,kernel->height*
1312  sizeof(*kernel->values)));
1313  if (kernel->values == (MagickRealType *) NULL)
1314  return(DestroyKernelInfo(kernel));
1315 
1316  /* set all kernel values within diamond area to scale given */
1317  for ( i=0, v=0; v < (ssize_t)kernel->height; v++)
1318  { size_t
1319  alpha = order_f / ( fact((size_t) v) * fact(kernel->height-v-1) );
1320  for ( u=0; u < (ssize_t)kernel->width; u++, i++)
1321  kernel->positive_range += kernel->values[i] = (double)
1322  (alpha * order_f / ( fact((size_t) u) * fact(kernel->height-u-1) ));
1323  }
1324  kernel->minimum = 1.0;
1325  kernel->maximum = kernel->values[kernel->x+kernel->y*kernel->width];
1326  kernel->negative_range = 0.0;
1327  break;
1328  }
1329 
1330  /*
1331  Convolution Kernels - Well Known Named Constant Kernels
1332  */
1333  case LaplacianKernel:
1334  { switch ( (int) args->rho ) {
1335  case 0:
1336  default: /* laplacian square filter -- default */
1337  kernel=ParseKernelArray("3: -1,-1,-1 -1,8,-1 -1,-1,-1");
1338  break;
1339  case 1: /* laplacian diamond filter */
1340  kernel=ParseKernelArray("3: 0,-1,0 -1,4,-1 0,-1,0");
1341  break;
1342  case 2:
1343  kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1344  break;
1345  case 3:
1346  kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 1,-2,1");
1347  break;
1348  case 5: /* a 5x5 laplacian */
1349  kernel=ParseKernelArray(
1350  "5: -4,-1,0,-1,-4 -1,2,3,2,-1 0,3,4,3,0 -1,2,3,2,-1 -4,-1,0,-1,-4");
1351  break;
1352  case 7: /* a 7x7 laplacian */
1353  kernel=ParseKernelArray(
1354  "7:-10,-5,-2,-1,-2,-5,-10 -5,0,3,4,3,0,-5 -2,3,6,7,6,3,-2 -1,4,7,8,7,4,-1 -2,3,6,7,6,3,-2 -5,0,3,4,3,0,-5 -10,-5,-2,-1,-2,-5,-10" );
1355  break;
1356  case 15: /* a 5x5 LoG (sigma approx 1.4) */
1357  kernel=ParseKernelArray(
1358  "5: 0,0,-1,0,0 0,-1,-2,-1,0 -1,-2,16,-2,-1 0,-1,-2,-1,0 0,0,-1,0,0");
1359  break;
1360  case 19: /* a 9x9 LoG (sigma approx 1.4) */
1361  /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
1362  kernel=ParseKernelArray(
1363  "9: 0,-1,-1,-2,-2,-2,-1,-1,0 -1,-2,-4,-5,-5,-5,-4,-2,-1 -1,-4,-5,-3,-0,-3,-5,-4,-1 -2,-5,-3,12,24,12,-3,-5,-2 -2,-5,-0,24,40,24,-0,-5,-2 -2,-5,-3,12,24,12,-3,-5,-2 -1,-4,-5,-3,-0,-3,-5,-4,-1 -1,-2,-4,-5,-5,-5,-4,-2,-1 0,-1,-1,-2,-2,-2,-1,-1,0");
1364  break;
1365  }
1366  if (kernel == (KernelInfo *) NULL)
1367  return(kernel);
1368  kernel->type = type;
1369  break;
1370  }
1371  case SobelKernel:
1372  { /* Simple Sobel Kernel */
1373  kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1374  if (kernel == (KernelInfo *) NULL)
1375  return(kernel);
1376  kernel->type = type;
1377  RotateKernelInfo(kernel, args->rho);
1378  break;
1379  }
1380  case RobertsKernel:
1381  {
1382  kernel=ParseKernelArray("3: 0,0,0 1,-1,0 0,0,0");
1383  if (kernel == (KernelInfo *) NULL)
1384  return(kernel);
1385  kernel->type = type;
1386  RotateKernelInfo(kernel, args->rho);
1387  break;
1388  }
1389  case PrewittKernel:
1390  {
1391  kernel=ParseKernelArray("3: 1,0,-1 1,0,-1 1,0,-1");
1392  if (kernel == (KernelInfo *) NULL)
1393  return(kernel);
1394  kernel->type = type;
1395  RotateKernelInfo(kernel, args->rho);
1396  break;
1397  }
1398  case CompassKernel:
1399  {
1400  kernel=ParseKernelArray("3: 1,1,-1 1,-2,-1 1,1,-1");
1401  if (kernel == (KernelInfo *) NULL)
1402  return(kernel);
1403  kernel->type = type;
1404  RotateKernelInfo(kernel, args->rho);
1405  break;
1406  }
1407  case KirschKernel:
1408  {
1409  kernel=ParseKernelArray("3: 5,-3,-3 5,0,-3 5,-3,-3");
1410  if (kernel == (KernelInfo *) NULL)
1411  return(kernel);
1412  kernel->type = type;
1413  RotateKernelInfo(kernel, args->rho);
1414  break;
1415  }
1416  case FreiChenKernel:
1417  /* Direction is set to be left to right positive */
1418  /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
1419  /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
1420  { switch ( (int) args->rho ) {
1421  default:
1422  case 0:
1423  kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1424  if (kernel == (KernelInfo *) NULL)
1425  return(kernel);
1426  kernel->type = type;
1427  kernel->values[3] = +(MagickRealType) MagickSQ2;
1428  kernel->values[5] = -(MagickRealType) MagickSQ2;
1429  CalcKernelMetaData(kernel); /* recalculate meta-data */
1430  break;
1431  case 2:
1432  kernel=ParseKernelArray("3: 1,2,0 2,0,-2 0,-2,-1");
1433  if (kernel == (KernelInfo *) NULL)
1434  return(kernel);
1435  kernel->type = type;
1436  kernel->values[1] = kernel->values[3]= +(MagickRealType) MagickSQ2;
1437  kernel->values[5] = kernel->values[7]= -(MagickRealType) MagickSQ2;
1438  CalcKernelMetaData(kernel); /* recalculate meta-data */
1439  ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1440  break;
1441  case 10:
1442  {
1443  kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19",exception);
1444  if (kernel == (KernelInfo *) NULL)
1445  return(kernel);
1446  break;
1447  }
1448  case 1:
1449  case 11:
1450  kernel=ParseKernelArray("3: 1,0,-1 2,0,-2 1,0,-1");
1451  if (kernel == (KernelInfo *) NULL)
1452  return(kernel);
1453  kernel->type = type;
1454  kernel->values[3] = +(MagickRealType) MagickSQ2;
1455  kernel->values[5] = -(MagickRealType) MagickSQ2;
1456  CalcKernelMetaData(kernel); /* recalculate meta-data */
1457  ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1458  break;
1459  case 12:
1460  kernel=ParseKernelArray("3: 1,2,1 0,0,0 1,2,1");
1461  if (kernel == (KernelInfo *) NULL)
1462  return(kernel);
1463  kernel->type = type;
1464  kernel->values[1] = +(MagickRealType) MagickSQ2;
1465  kernel->values[7] = +(MagickRealType) MagickSQ2;
1466  CalcKernelMetaData(kernel);
1467  ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1468  break;
1469  case 13:
1470  kernel=ParseKernelArray("3: 2,-1,0 -1,0,1 0,1,-2");
1471  if (kernel == (KernelInfo *) NULL)
1472  return(kernel);
1473  kernel->type = type;
1474  kernel->values[0] = +(MagickRealType) MagickSQ2;
1475  kernel->values[8] = -(MagickRealType) MagickSQ2;
1476  CalcKernelMetaData(kernel);
1477  ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1478  break;
1479  case 14:
1480  kernel=ParseKernelArray("3: 0,1,-2 -1,0,1 2,-1,0");
1481  if (kernel == (KernelInfo *) NULL)
1482  return(kernel);
1483  kernel->type = type;
1484  kernel->values[2] = -(MagickRealType) MagickSQ2;
1485  kernel->values[6] = +(MagickRealType) MagickSQ2;
1486  CalcKernelMetaData(kernel);
1487  ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
1488  break;
1489  case 15:
1490  kernel=ParseKernelArray("3: 0,-1,0 1,0,1 0,-1,0");
1491  if (kernel == (KernelInfo *) NULL)
1492  return(kernel);
1493  kernel->type = type;
1494  ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1495  break;
1496  case 16:
1497  kernel=ParseKernelArray("3: 1,0,-1 0,0,0 -1,0,1");
1498  if (kernel == (KernelInfo *) NULL)
1499  return(kernel);
1500  kernel->type = type;
1501  ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
1502  break;
1503  case 17:
1504  kernel=ParseKernelArray("3: 1,-2,1 -2,4,-2 -1,-2,1");
1505  if (kernel == (KernelInfo *) NULL)
1506  return(kernel);
1507  kernel->type = type;
1508  ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1509  break;
1510  case 18:
1511  kernel=ParseKernelArray("3: -2,1,-2 1,4,1 -2,1,-2");
1512  if (kernel == (KernelInfo *) NULL)
1513  return(kernel);
1514  kernel->type = type;
1515  ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
1516  break;
1517  case 19:
1518  kernel=ParseKernelArray("3: 1,1,1 1,1,1 1,1,1");
1519  if (kernel == (KernelInfo *) NULL)
1520  return(kernel);
1521  kernel->type = type;
1522  ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
1523  break;
1524  }
1525  if ( fabs(args->sigma) >= MagickEpsilon )
1526  /* Rotate by correctly supplied 'angle' */
1527  RotateKernelInfo(kernel, args->sigma);
1528  else if ( args->rho > 30.0 || args->rho < -30.0 )
1529  /* Rotate by out of bounds 'type' */
1530  RotateKernelInfo(kernel, args->rho);
1531  break;
1532  }
1533 
1534  /*
1535  Boolean or Shaped Kernels
1536  */
1537  case DiamondKernel:
1538  {
1539  if (args->rho < 1.0)
1540  kernel->width = kernel->height = 3; /* default radius = 1 */
1541  else
1542  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1543  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1544 
1546  AcquireAlignedMemory(kernel->width,kernel->height*
1547  sizeof(*kernel->values)));
1548  if (kernel->values == (MagickRealType *) NULL)
1549  return(DestroyKernelInfo(kernel));
1550 
1551  /* set all kernel values within diamond area to scale given */
1552  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1553  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1554  if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
1555  kernel->positive_range += kernel->values[i] = args->sigma;
1556  else
1557  kernel->values[i] = nan;
1558  kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1559  break;
1560  }
1561  case SquareKernel:
1562  case RectangleKernel:
1563  { double
1564  scale;
1565  if ( type == SquareKernel )
1566  {
1567  if (args->rho < 1.0)
1568  kernel->width = kernel->height = 3; /* default radius = 1 */
1569  else
1570  kernel->width = kernel->height = (size_t) (2*args->rho+1);
1571  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1572  scale = args->sigma;
1573  }
1574  else {
1575  /* NOTE: user defaults set in "AcquireKernelInfo()" */
1576  if ( args->rho < 1.0 || args->sigma < 1.0 )
1577  return(DestroyKernelInfo(kernel)); /* invalid args given */
1578  kernel->width = (size_t)args->rho;
1579  kernel->height = (size_t)args->sigma;
1580  if ( args->xi < 0.0 || args->xi > (double)kernel->width ||
1581  args->psi < 0.0 || args->psi > (double)kernel->height )
1582  return(DestroyKernelInfo(kernel)); /* invalid args given */
1583  kernel->x = (ssize_t) args->xi;
1584  kernel->y = (ssize_t) args->psi;
1585  scale = 1.0;
1586  }
1588  AcquireAlignedMemory(kernel->width,kernel->height*
1589  sizeof(*kernel->values)));
1590  if (kernel->values == (MagickRealType *) NULL)
1591  return(DestroyKernelInfo(kernel));
1592 
1593  /* set all kernel values to scale given */
1594  u=(ssize_t) (kernel->width*kernel->height);
1595  for ( i=0; i < u; i++)
1596  kernel->values[i] = scale;
1597  kernel->minimum = kernel->maximum = scale; /* a flat shape */
1598  kernel->positive_range = scale*u;
1599  break;
1600  }
1601  case OctagonKernel:
1602  {
1603  if (args->rho < 1.0)
1604  kernel->width = kernel->height = 5; /* default radius = 2 */
1605  else
1606  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1607  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1608 
1610  AcquireAlignedMemory(kernel->width,kernel->height*
1611  sizeof(*kernel->values)));
1612  if (kernel->values == (MagickRealType *) NULL)
1613  return(DestroyKernelInfo(kernel));
1614 
1615  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1616  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1617  if ( (labs((long) u)+labs((long) v)) <=
1618  ((long)kernel->x + (long)(kernel->x/2)) )
1619  kernel->positive_range += kernel->values[i] = args->sigma;
1620  else
1621  kernel->values[i] = nan;
1622  kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1623  break;
1624  }
1625  case DiskKernel:
1626  {
1627  ssize_t
1628  limit = (ssize_t)(args->rho*args->rho);
1629 
1630  if (args->rho < 0.4) /* default radius approx 4.3 */
1631  kernel->width = kernel->height = 9L, limit = 18L;
1632  else
1633  kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
1634  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1635 
1637  AcquireAlignedMemory(kernel->width,kernel->height*
1638  sizeof(*kernel->values)));
1639  if (kernel->values == (MagickRealType *) NULL)
1640  return(DestroyKernelInfo(kernel));
1641 
1642  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1643  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1644  if ((u*u+v*v) <= limit)
1645  kernel->positive_range += kernel->values[i] = args->sigma;
1646  else
1647  kernel->values[i] = nan;
1648  kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1649  break;
1650  }
1651  case PlusKernel:
1652  {
1653  if (args->rho < 1.0)
1654  kernel->width = kernel->height = 5; /* default radius 2 */
1655  else
1656  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1657  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1658 
1660  AcquireAlignedMemory(kernel->width,kernel->height*
1661  sizeof(*kernel->values)));
1662  if (kernel->values == (MagickRealType *) NULL)
1663  return(DestroyKernelInfo(kernel));
1664 
1665  /* set all kernel values along axises to given scale */
1666  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1667  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1668  kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
1669  kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1670  kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1671  break;
1672  }
1673  case CrossKernel:
1674  {
1675  if (args->rho < 1.0)
1676  kernel->width = kernel->height = 5; /* default radius 2 */
1677  else
1678  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
1679  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1680 
1682  AcquireAlignedMemory(kernel->width,kernel->height*
1683  sizeof(*kernel->values)));
1684  if (kernel->values == (MagickRealType *) NULL)
1685  return(DestroyKernelInfo(kernel));
1686 
1687  /* set all kernel values along axises to given scale */
1688  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
1689  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1690  kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
1691  kernel->minimum = kernel->maximum = args->sigma; /* a flat shape */
1692  kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
1693  break;
1694  }
1695  /*
1696  HitAndMiss Kernels
1697  */
1698  case RingKernel:
1699  case PeaksKernel:
1700  {
1701  ssize_t
1702  limit1,
1703  limit2,
1704  scale;
1705 
1706  if (args->rho < args->sigma)
1707  {
1708  kernel->width = ((size_t)args->sigma)*2+1;
1709  limit1 = (ssize_t)(args->rho*args->rho);
1710  limit2 = (ssize_t)(args->sigma*args->sigma);
1711  }
1712  else
1713  {
1714  kernel->width = ((size_t)args->rho)*2+1;
1715  limit1 = (ssize_t)(args->sigma*args->sigma);
1716  limit2 = (ssize_t)(args->rho*args->rho);
1717  }
1718  if ( limit2 <= 0 )
1719  kernel->width = 7L, limit1 = 7L, limit2 = 11L;
1720 
1721  kernel->height = kernel->width;
1722  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
1724  AcquireAlignedMemory(kernel->width,kernel->height*
1725  sizeof(*kernel->values)));
1726  if (kernel->values == (MagickRealType *) NULL)
1727  return(DestroyKernelInfo(kernel));
1728 
1729  /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
1730  scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
1731  for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
1732  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
1733  { ssize_t radius=u*u+v*v;
1734  if (limit1 < radius && radius <= limit2)
1735  kernel->positive_range += kernel->values[i] = (double) scale;
1736  else
1737  kernel->values[i] = nan;
1738  }
1739  kernel->minimum = kernel->maximum = (double) scale;
1740  if ( type == PeaksKernel ) {
1741  /* set the central point in the middle */
1742  kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
1743  kernel->positive_range = 1.0;
1744  kernel->maximum = 1.0;
1745  }
1746  break;
1747  }
1748  case EdgesKernel:
1749  {
1750  kernel=AcquireKernelInfo("ThinSE:482",exception);
1751  if (kernel == (KernelInfo *) NULL)
1752  return(kernel);
1753  kernel->type = type;
1754  ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
1755  break;
1756  }
1757  case CornersKernel:
1758  {
1759  kernel=AcquireKernelInfo("ThinSE:87",exception);
1760  if (kernel == (KernelInfo *) NULL)
1761  return(kernel);
1762  kernel->type = type;
1763  ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
1764  break;
1765  }
1766  case DiagonalsKernel:
1767  {
1768  switch ( (int) args->rho ) {
1769  case 0:
1770  default:
1771  { KernelInfo
1772  *new_kernel;
1773  kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1774  if (kernel == (KernelInfo *) NULL)
1775  return(kernel);
1776  kernel->type = type;
1777  new_kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1778  if (new_kernel == (KernelInfo *) NULL)
1779  return(DestroyKernelInfo(kernel));
1780  new_kernel->type = type;
1781  LastKernelInfo(kernel)->next = new_kernel;
1782  ExpandMirrorKernelInfo(kernel);
1783  return(kernel);
1784  }
1785  case 1:
1786  kernel=ParseKernelArray("3: 0,0,0 0,-,1 1,1,-");
1787  break;
1788  case 2:
1789  kernel=ParseKernelArray("3: 0,0,1 0,-,1 0,1,-");
1790  break;
1791  }
1792  if (kernel == (KernelInfo *) NULL)
1793  return(kernel);
1794  kernel->type = type;
1795  RotateKernelInfo(kernel, args->sigma);
1796  break;
1797  }
1798  case LineEndsKernel:
1799  { /* Kernels for finding the end of thin lines */
1800  switch ( (int) args->rho ) {
1801  case 0:
1802  default:
1803  /* set of kernels to find all end of lines */
1804  return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>",exception));
1805  case 1:
1806  /* kernel for 4-connected line ends - no rotation */
1807  kernel=ParseKernelArray("3: 0,0,- 0,1,1 0,0,-");
1808  break;
1809  case 2:
1810  /* kernel to add for 8-connected lines - no rotation */
1811  kernel=ParseKernelArray("3: 0,0,0 0,1,0 0,0,1");
1812  break;
1813  case 3:
1814  /* kernel to add for orthogonal line ends - does not find corners */
1815  kernel=ParseKernelArray("3: 0,0,0 0,1,1 0,0,0");
1816  break;
1817  case 4:
1818  /* traditional line end - fails on last T end */
1819  kernel=ParseKernelArray("3: 0,0,0 0,1,- 0,0,-");
1820  break;
1821  }
1822  if (kernel == (KernelInfo *) NULL)
1823  return(kernel);
1824  kernel->type = type;
1825  RotateKernelInfo(kernel, args->sigma);
1826  break;
1827  }
1828  case LineJunctionsKernel:
1829  { /* kernels for finding the junctions of multiple lines */
1830  switch ( (int) args->rho ) {
1831  case 0:
1832  default:
1833  /* set of kernels to find all line junctions */
1834  return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>",exception));
1835  case 1:
1836  /* Y Junction */
1837  kernel=ParseKernelArray("3: 1,-,1 -,1,- -,1,-");
1838  break;
1839  case 2:
1840  /* Diagonal T Junctions */
1841  kernel=ParseKernelArray("3: 1,-,- -,1,- 1,-,1");
1842  break;
1843  case 3:
1844  /* Orthogonal T Junctions */
1845  kernel=ParseKernelArray("3: -,-,- 1,1,1 -,1,-");
1846  break;
1847  case 4:
1848  /* Diagonal X Junctions */
1849  kernel=ParseKernelArray("3: 1,-,1 -,1,- 1,-,1");
1850  break;
1851  case 5:
1852  /* Orthogonal X Junctions - minimal diamond kernel */
1853  kernel=ParseKernelArray("3: -,1,- 1,1,1 -,1,-");
1854  break;
1855  }
1856  if (kernel == (KernelInfo *) NULL)
1857  return(kernel);
1858  kernel->type = type;
1859  RotateKernelInfo(kernel, args->sigma);
1860  break;
1861  }
1862  case RidgesKernel:
1863  { /* Ridges - Ridge finding kernels */
1864  KernelInfo
1865  *new_kernel;
1866  switch ( (int) args->rho ) {
1867  case 1:
1868  default:
1869  kernel=ParseKernelArray("3x1:0,1,0");
1870  if (kernel == (KernelInfo *) NULL)
1871  return(kernel);
1872  kernel->type = type;
1873  ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
1874  break;
1875  case 2:
1876  kernel=ParseKernelArray("4x1:0,1,1,0");
1877  if (kernel == (KernelInfo *) NULL)
1878  return(kernel);
1879  kernel->type = type;
1880  ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
1881 
1882  /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
1883  /* Unfortunatally we can not yet rotate a non-square kernel */
1884  /* But then we can't flip a non-symetrical kernel either */
1885  new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
1886  if (new_kernel == (KernelInfo *) NULL)
1887  return(DestroyKernelInfo(kernel));
1888  new_kernel->type = type;
1889  LastKernelInfo(kernel)->next = new_kernel;
1890  new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
1891  if (new_kernel == (KernelInfo *) NULL)
1892  return(DestroyKernelInfo(kernel));
1893  new_kernel->type = type;
1894  LastKernelInfo(kernel)->next = new_kernel;
1895  new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
1896  if (new_kernel == (KernelInfo *) NULL)
1897  return(DestroyKernelInfo(kernel));
1898  new_kernel->type = type;
1899  LastKernelInfo(kernel)->next = new_kernel;
1900  new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
1901  if (new_kernel == (KernelInfo *) NULL)
1902  return(DestroyKernelInfo(kernel));
1903  new_kernel->type = type;
1904  LastKernelInfo(kernel)->next = new_kernel;
1905  new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
1906  if (new_kernel == (KernelInfo *) NULL)
1907  return(DestroyKernelInfo(kernel));
1908  new_kernel->type = type;
1909  LastKernelInfo(kernel)->next = new_kernel;
1910  new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
1911  if (new_kernel == (KernelInfo *) NULL)
1912  return(DestroyKernelInfo(kernel));
1913  new_kernel->type = type;
1914  LastKernelInfo(kernel)->next = new_kernel;
1915  new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
1916  if (new_kernel == (KernelInfo *) NULL)
1917  return(DestroyKernelInfo(kernel));
1918  new_kernel->type = type;
1919  LastKernelInfo(kernel)->next = new_kernel;
1920  new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
1921  if (new_kernel == (KernelInfo *) NULL)
1922  return(DestroyKernelInfo(kernel));
1923  new_kernel->type = type;
1924  LastKernelInfo(kernel)->next = new_kernel;
1925  break;
1926  }
1927  break;
1928  }
1929  case ConvexHullKernel:
1930  {
1931  KernelInfo
1932  *new_kernel;
1933  /* first set of 8 kernels */
1934  kernel=ParseKernelArray("3: 1,1,- 1,0,- 1,-,0");
1935  if (kernel == (KernelInfo *) NULL)
1936  return(kernel);
1937  kernel->type = type;
1938  ExpandRotateKernelInfo(kernel, 90.0);
1939  /* append the mirror versions too - no flip function yet */
1940  new_kernel=ParseKernelArray("3: 1,1,1 1,0,- -,-,0");
1941  if (new_kernel == (KernelInfo *) NULL)
1942  return(DestroyKernelInfo(kernel));
1943  new_kernel->type = type;
1944  ExpandRotateKernelInfo(new_kernel, 90.0);
1945  LastKernelInfo(kernel)->next = new_kernel;
1946  break;
1947  }
1948  case SkeletonKernel:
1949  {
1950  switch ( (int) args->rho ) {
1951  case 1:
1952  default:
1953  /* Traditional Skeleton...
1954  ** A cyclically rotated single kernel
1955  */
1956  kernel=AcquireKernelInfo("ThinSE:482",exception);
1957  if (kernel == (KernelInfo *) NULL)
1958  return(kernel);
1959  kernel->type = type;
1960  ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
1961  break;
1962  case 2:
1963  /* HIPR Variation of the cyclic skeleton
1964  ** Corners of the traditional method made more forgiving,
1965  ** but the retain the same cyclic order.
1966  */
1967  kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;",exception);
1968  if (kernel == (KernelInfo *) NULL)
1969  return(kernel);
1970  if (kernel->next == (KernelInfo *) NULL)
1971  return(DestroyKernelInfo(kernel));
1972  kernel->type = type;
1973  kernel->next->type = type;
1974  ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
1975  break;
1976  case 3:
1977  /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
1978  ** "Connectivity-Preserving Morphological Image Thransformations"
1979  ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1980  ** http://www.leptonica.com/papers/conn.pdf
1981  */
1982  kernel=AcquireKernelInfo("ThinSE:41; ThinSE:42; ThinSE:43",
1983  exception);
1984  if (kernel == (KernelInfo *) NULL)
1985  return(kernel);
1986  kernel->type = type;
1987  kernel->next->type = type;
1988  kernel->next->next->type = type;
1989  ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
1990  break;
1991  }
1992  break;
1993  }
1994  case ThinSEKernel:
1995  { /* Special kernels for general thinning, while preserving connections
1996  ** "Connectivity-Preserving Morphological Image Thransformations"
1997  ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
1998  ** http://www.leptonica.com/papers/conn.pdf
1999  ** And
2000  ** http://tpgit.github.com/Leptonica/ccthin_8c_source.html
2001  **
2002  ** Note kernels do not specify the origin pixel, allowing them
2003  ** to be used for both thickening and thinning operations.
2004  */
2005  switch ( (int) args->rho ) {
2006  /* SE for 4-connected thinning */
2007  case 41: /* SE_4_1 */
2008  kernel=ParseKernelArray("3: -,-,1 0,-,1 -,-,1");
2009  break;
2010  case 42: /* SE_4_2 */
2011  kernel=ParseKernelArray("3: -,-,1 0,-,1 -,0,-");
2012  break;
2013  case 43: /* SE_4_3 */
2014  kernel=ParseKernelArray("3: -,0,- 0,-,1 -,-,1");
2015  break;
2016  case 44: /* SE_4_4 */
2017  kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,-");
2018  break;
2019  case 45: /* SE_4_5 */
2020  kernel=ParseKernelArray("3: -,0,1 0,-,1 -,0,-");
2021  break;
2022  case 46: /* SE_4_6 */
2023  kernel=ParseKernelArray("3: -,0,- 0,-,1 -,0,1");
2024  break;
2025  case 47: /* SE_4_7 */
2026  kernel=ParseKernelArray("3: -,1,1 0,-,1 -,0,-");
2027  break;
2028  case 48: /* SE_4_8 */
2029  kernel=ParseKernelArray("3: -,-,1 0,-,1 0,-,1");
2030  break;
2031  case 49: /* SE_4_9 */
2032  kernel=ParseKernelArray("3: 0,-,1 0,-,1 -,-,1");
2033  break;
2034  /* SE for 8-connected thinning - negatives of the above */
2035  case 81: /* SE_8_0 */
2036  kernel=ParseKernelArray("3: -,1,- 0,-,1 -,1,-");
2037  break;
2038  case 82: /* SE_8_2 */
2039  kernel=ParseKernelArray("3: -,1,- 0,-,1 0,-,-");
2040  break;
2041  case 83: /* SE_8_3 */
2042  kernel=ParseKernelArray("3: 0,-,- 0,-,1 -,1,-");
2043  break;
2044  case 84: /* SE_8_4 */
2045  kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,-");
2046  break;
2047  case 85: /* SE_8_5 */
2048  kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,-");
2049  break;
2050  case 86: /* SE_8_6 */
2051  kernel=ParseKernelArray("3: 0,-,- 0,-,1 0,-,1");
2052  break;
2053  case 87: /* SE_8_7 */
2054  kernel=ParseKernelArray("3: -,1,- 0,-,1 0,0,-");
2055  break;
2056  case 88: /* SE_8_8 */
2057  kernel=ParseKernelArray("3: -,1,- 0,-,1 0,1,-");
2058  break;
2059  case 89: /* SE_8_9 */
2060  kernel=ParseKernelArray("3: 0,1,- 0,-,1 -,1,-");
2061  break;
2062  /* Special combined SE kernels */
2063  case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
2064  kernel=ParseKernelArray("3: -,-,1 0,-,- -,0,-");
2065  break;
2066  case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
2067  kernel=ParseKernelArray("3: -,1,- -,-,1 0,-,-");
2068  break;
2069  case 481: /* SE_48_1 - General Connected Corner Kernel */
2070  kernel=ParseKernelArray("3: -,1,1 0,-,1 0,0,-");
2071  break;
2072  default:
2073  case 482: /* SE_48_2 - General Edge Kernel */
2074  kernel=ParseKernelArray("3: 0,-,1 0,-,1 0,-,1");
2075  break;
2076  }
2077  if (kernel == (KernelInfo *) NULL)
2078  return(kernel);
2079  kernel->type = type;
2080  RotateKernelInfo(kernel, args->sigma);
2081  break;
2082  }
2083  /*
2084  Distance Measuring Kernels
2085  */
2086  case ChebyshevKernel:
2087  {
2088  if (args->rho < 1.0)
2089  kernel->width = kernel->height = 3; /* default radius = 1 */
2090  else
2091  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2092  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2093 
2095  AcquireAlignedMemory(kernel->width,kernel->height*
2096  sizeof(*kernel->values)));
2097  if (kernel->values == (MagickRealType *) NULL)
2098  return(DestroyKernelInfo(kernel));
2099 
2100  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2101  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2102  kernel->positive_range += ( kernel->values[i] =
2103  args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
2104  kernel->maximum = kernel->values[0];
2105  break;
2106  }
2107  case ManhattanKernel:
2108  {
2109  if (args->rho < 1.0)
2110  kernel->width = kernel->height = 3; /* default radius = 1 */
2111  else
2112  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2113  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2114 
2116  AcquireAlignedMemory(kernel->width,kernel->height*
2117  sizeof(*kernel->values)));
2118  if (kernel->values == (MagickRealType *) NULL)
2119  return(DestroyKernelInfo(kernel));
2120 
2121  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2122  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2123  kernel->positive_range += ( kernel->values[i] =
2124  args->sigma*(labs((long) u)+labs((long) v)) );
2125  kernel->maximum = kernel->values[0];
2126  break;
2127  }
2128  case OctagonalKernel:
2129  {
2130  if (args->rho < 2.0)
2131  kernel->width = kernel->height = 5; /* default/minimum radius = 2 */
2132  else
2133  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2134  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2135 
2137  AcquireAlignedMemory(kernel->width,kernel->height*
2138  sizeof(*kernel->values)));
2139  if (kernel->values == (MagickRealType *) NULL)
2140  return(DestroyKernelInfo(kernel));
2141 
2142  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2143  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2144  {
2145  double
2146  r1 = MagickMax(fabs((double)u),fabs((double)v)),
2147  r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
2148  kernel->positive_range += kernel->values[i] =
2149  args->sigma*MagickMax(r1,r2);
2150  }
2151  kernel->maximum = kernel->values[0];
2152  break;
2153  }
2154  case EuclideanKernel:
2155  {
2156  if (args->rho < 1.0)
2157  kernel->width = kernel->height = 3; /* default radius = 1 */
2158  else
2159  kernel->width = kernel->height = ((size_t)args->rho)*2+1;
2160  kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
2161 
2163  AcquireAlignedMemory(kernel->width,kernel->height*
2164  sizeof(*kernel->values)));
2165  if (kernel->values == (MagickRealType *) NULL)
2166  return(DestroyKernelInfo(kernel));
2167 
2168  for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
2169  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
2170  kernel->positive_range += ( kernel->values[i] =
2171  args->sigma*sqrt((double)(u*u+v*v)) );
2172  kernel->maximum = kernel->values[0];
2173  break;
2174  }
2175  default:
2176  {
2177  /* No-Op Kernel - Basically just a single pixel on its own */
2178  kernel=ParseKernelArray("1:1");
2179  if (kernel == (KernelInfo *) NULL)
2180  return(kernel);
2181  kernel->type = UndefinedKernel;
2182  break;
2183  }
2184  break;
2185  }
2186  return(kernel);
2187 }
2188 
2189 /*
2190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2191 % %
2192 % %
2193 % %
2194 % C l o n e K e r n e l I n f o %
2195 % %
2196 % %
2197 % %
2198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2199 %
2200 % CloneKernelInfo() creates a new clone of the given Kernel List so that its
2201 % can be modified without effecting the original. The cloned kernel should
2202 % be destroyed using DestoryKernelInfo() when no longer needed.
2203 %
2204 % The format of the CloneKernelInfo method is:
2205 %
2206 % KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
2207 %
2208 % A description of each parameter follows:
2209 %
2210 % o kernel: the Morphology/Convolution kernel to be cloned
2211 %
2212 */
2214 {
2215  register ssize_t
2216  i;
2217 
2218  KernelInfo
2219  *new_kernel;
2220 
2221  assert(kernel != (KernelInfo *) NULL);
2222  new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
2223  if (new_kernel == (KernelInfo *) NULL)
2224  return(new_kernel);
2225  *new_kernel=(*kernel); /* copy values in structure */
2226 
2227  /* replace the values with a copy of the values */
2228  new_kernel->values=(MagickRealType *) MagickAssumeAligned(
2229  AcquireAlignedMemory(kernel->width,kernel->height*sizeof(*kernel->values)));
2230  if (new_kernel->values == (MagickRealType *) NULL)
2231  return(DestroyKernelInfo(new_kernel));
2232  for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
2233  new_kernel->values[i]=kernel->values[i];
2234 
2235  /* Also clone the next kernel in the kernel list */
2236  if ( kernel->next != (KernelInfo *) NULL ) {
2237  new_kernel->next = CloneKernelInfo(kernel->next);
2238  if ( new_kernel->next == (KernelInfo *) NULL )
2239  return(DestroyKernelInfo(new_kernel));
2240  }
2241 
2242  return(new_kernel);
2243 }
2244 
2245 /*
2246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2247 % %
2248 % %
2249 % %
2250 % D e s t r o y K e r n e l I n f o %
2251 % %
2252 % %
2253 % %
2254 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2255 %
2256 % DestroyKernelInfo() frees the memory used by a Convolution/Morphology
2257 % kernel.
2258 %
2259 % The format of the DestroyKernelInfo method is:
2260 %
2261 % KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
2262 %
2263 % A description of each parameter follows:
2264 %
2265 % o kernel: the Morphology/Convolution kernel to be destroyed
2266 %
2267 */
2269 {
2270  assert(kernel != (KernelInfo *) NULL);
2271  if (kernel->next != (KernelInfo *) NULL)
2272  kernel->next=DestroyKernelInfo(kernel->next);
2273  kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
2274  kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
2275  return(kernel);
2276 }
2277 
2278 /*
2279 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2280 % %
2281 % %
2282 % %
2283 + E x p a n d M i r r o r K e r n e l I n f o %
2284 % %
2285 % %
2286 % %
2287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2288 %
2289 % ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
2290 % sequence of 90-degree rotated kernels but providing a reflected 180
2291 % rotatation, before the -/+ 90-degree rotations.
2292 %
2293 % This special rotation order produces a better, more symetrical thinning of
2294 % objects.
2295 %
2296 % The format of the ExpandMirrorKernelInfo method is:
2297 %
2298 % void ExpandMirrorKernelInfo(KernelInfo *kernel)
2299 %
2300 % A description of each parameter follows:
2301 %
2302 % o kernel: the Morphology/Convolution kernel
2303 %
2304 % This function is only internel to this module, as it is not finalized,
2305 % especially with regard to non-orthogonal angles, and rotation of larger
2306 % 2D kernels.
2307 */
2308 
2309 #if 0
2310 static void FlopKernelInfo(KernelInfo *kernel)
2311  { /* Do a Flop by reversing each row. */
2312  size_t
2313  y;
2314  register ssize_t
2315  x,r;
2316  register double
2317  *k,t;
2318 
2319  for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
2320  for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
2321  t=k[x], k[x]=k[r], k[r]=t;
2322 
2323  kernel->x = kernel->width - kernel->x - 1;
2324  angle = fmod(angle+180.0, 360.0);
2325  }
2326 #endif
2327 
2329 {
2330  KernelInfo
2331  *clone,
2332  *last;
2333 
2334  last = kernel;
2335 
2336  clone = CloneKernelInfo(last);
2337  if (clone == (KernelInfo *) NULL)
2338  return;
2339  RotateKernelInfo(clone, 180); /* flip */
2340  LastKernelInfo(last)->next = clone;
2341  last = clone;
2342 
2343  clone = CloneKernelInfo(last);
2344  if (clone == (KernelInfo *) NULL)
2345  return;
2346  RotateKernelInfo(clone, 90); /* transpose */
2347  LastKernelInfo(last)->next = clone;
2348  last = clone;
2349 
2350  clone = CloneKernelInfo(last);
2351  if (clone == (KernelInfo *) NULL)
2352  return;
2353  RotateKernelInfo(clone, 180); /* flop */
2354  LastKernelInfo(last)->next = clone;
2355 
2356  return;
2357 }
2358 
2359 /*
2360 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2361 % %
2362 % %
2363 % %
2364 + E x p a n d R o t a t e K e r n e l I n f o %
2365 % %
2366 % %
2367 % %
2368 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2369 %
2370 % ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
2371 % incrementally by the angle given, until the kernel repeats.
2372 %
2373 % WARNING: 45 degree rotations only works for 3x3 kernels.
2374 % While 90 degree roatations only works for linear and square kernels
2375 %
2376 % The format of the ExpandRotateKernelInfo method is:
2377 %
2378 % void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
2379 %
2380 % A description of each parameter follows:
2381 %
2382 % o kernel: the Morphology/Convolution kernel
2383 %
2384 % o angle: angle to rotate in degrees
2385 %
2386 % This function is only internel to this module, as it is not finalized,
2387 % especially with regard to non-orthogonal angles, and rotation of larger
2388 % 2D kernels.
2389 */
2390 
2391 /* Internal Routine - Return true if two kernels are the same */
2393  const KernelInfo *kernel2)
2394 {
2395  register size_t
2396  i;
2397 
2398  /* check size and origin location */
2399  if ( kernel1->width != kernel2->width
2400  || kernel1->height != kernel2->height
2401  || kernel1->x != kernel2->x
2402  || kernel1->y != kernel2->y )
2403  return MagickFalse;
2404 
2405  /* check actual kernel values */
2406  for (i=0; i < (kernel1->width*kernel1->height); i++) {
2407  /* Test for Nan equivalence */
2408  if ( IsNaN(kernel1->values[i]) && !IsNaN(kernel2->values[i]) )
2409  return MagickFalse;
2410  if ( IsNaN(kernel2->values[i]) && !IsNaN(kernel1->values[i]) )
2411  return MagickFalse;
2412  /* Test actual values are equivalent */
2413  if ( fabs(kernel1->values[i] - kernel2->values[i]) >= MagickEpsilon )
2414  return MagickFalse;
2415  }
2416 
2417  return MagickTrue;
2418 }
2419 
2420 static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
2421 {
2422  KernelInfo
2423  *clone_info,
2424  *last;
2425 
2426  last=kernel;
2427 DisableMSCWarning(4127)
2428  while (1) {
2430  clone_info=CloneKernelInfo(last);
2431  if (clone_info == (KernelInfo *) NULL)
2432  break;
2433  RotateKernelInfo(clone_info,angle);
2434  if (SameKernelInfo(kernel,clone_info) != MagickFalse)
2435  break;
2436  LastKernelInfo(last)->next=clone_info;
2437  last=clone_info;
2438  }
2439  if (clone_info != (KernelInfo *) NULL)
2440  clone_info=DestroyKernelInfo(clone_info); /* kernel repeated - junk */
2441  return;
2442 }
2443 
2444 /*
2445 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2446 % %
2447 % %
2448 % %
2449 + C a l c M e t a K e r n a l I n f o %
2450 % %
2451 % %
2452 % %
2453 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2454 %
2455 % CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
2456 % using the kernel values. This should only ne used if it is not possible to
2457 % calculate that meta-data in some easier way.
2458 %
2459 % It is important that the meta-data is correct before ScaleKernelInfo() is
2460 % used to perform kernel normalization.
2461 %
2462 % The format of the CalcKernelMetaData method is:
2463 %
2464 % void CalcKernelMetaData(KernelInfo *kernel, const double scale )
2465 %
2466 % A description of each parameter follows:
2467 %
2468 % o kernel: the Morphology/Convolution kernel to modify
2469 %
2470 % WARNING: Minimum and Maximum values are assumed to include zero, even if
2471 % zero is not part of the kernel (as in Gaussian Derived kernels). This
2472 % however is not true for flat-shaped morphological kernels.
2473 %
2474 % WARNING: Only the specific kernel pointed to is modified, not a list of
2475 % multiple kernels.
2476 %
2477 % This is an internal function and not expected to be useful outside this
2478 % module. This could change however.
2479 */
2480 static void CalcKernelMetaData(KernelInfo *kernel)
2481 {
2482  register size_t
2483  i;
2484 
2485  kernel->minimum = kernel->maximum = 0.0;
2486  kernel->negative_range = kernel->positive_range = 0.0;
2487  for (i=0; i < (kernel->width*kernel->height); i++)
2488  {
2489  if ( fabs(kernel->values[i]) < MagickEpsilon )
2490  kernel->values[i] = 0.0;
2491  ( kernel->values[i] < 0)
2492  ? ( kernel->negative_range += kernel->values[i] )
2493  : ( kernel->positive_range += kernel->values[i] );
2494  Minimize(kernel->minimum, kernel->values[i]);
2495  Maximize(kernel->maximum, kernel->values[i]);
2496  }
2497 
2498  return;
2499 }
2500 
2501 /*
2502 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2503 % %
2504 % %
2505 % %
2506 % M o r p h o l o g y A p p l y %
2507 % %
2508 % %
2509 % %
2510 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2511 %
2512 % MorphologyApply() applies a morphological method, multiple times using
2513 % a list of multiple kernels. This is the method that should be called by
2514 % other 'operators' that internally use morphology operations as part of
2515 % their processing.
2516 %
2517 % It is basically equivalent to as MorphologyImage() (see below) but without
2518 % any user controls. This allows internel programs to use this method to
2519 % perform a specific task without possible interference by any API user
2520 % supplied settings.
2521 %
2522 % It is MorphologyImage() task to extract any such user controls, and
2523 % pass them to this function for processing.
2524 %
2525 % More specifically all given kernels should already be scaled, normalised,
2526 % and blended appropriatally before being parred to this routine. The
2527 % appropriate bias, and compose (typically 'UndefinedComposeOp') given.
2528 %
2529 % The format of the MorphologyApply method is:
2530 %
2531 % Image *MorphologyApply(const Image *image,MorphologyMethod method,
2532 % const ssize_t iterations,const KernelInfo *kernel,
2533 % const CompositeMethod compose,const double bias,
2534 % ExceptionInfo *exception)
2535 %
2536 % A description of each parameter follows:
2537 %
2538 % o image: the source image
2539 %
2540 % o method: the morphology method to be applied.
2541 %
2542 % o iterations: apply the operation this many times (or no change).
2543 % A value of -1 means loop until no change found.
2544 % How this is applied may depend on the morphology method.
2545 % Typically this is a value of 1.
2546 %
2547 % o channel: the channel type.
2548 %
2549 % o kernel: An array of double representing the morphology kernel.
2550 %
2551 % o compose: How to handle or merge multi-kernel results.
2552 % If 'UndefinedCompositeOp' use default for the Morphology method.
2553 % If 'NoCompositeOp' force image to be re-iterated by each kernel.
2554 % Otherwise merge the results using the compose method given.
2555 %
2556 % o bias: Convolution Output Bias.
2557 %
2558 % o exception: return any errors or warnings in this structure.
2559 %
2560 */
2561 static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
2562  const MorphologyMethod method,const KernelInfo *kernel,const double bias,
2563  ExceptionInfo *exception)
2564 {
2565 #define MorphologyTag "Morphology/Image"
2566 
2567  CacheView
2568  *image_view,
2569  *morphology_view;
2570 
2571  OffsetInfo
2572  offset;
2573 
2574  register ssize_t
2575  j,
2576  y;
2577 
2578  size_t
2579  *changes,
2580  changed,
2581  width;
2582 
2584  status;
2585 
2587  progress;
2588 
2589  assert(image != (Image *) NULL);
2590  assert(image->signature == MagickCoreSignature);
2591  assert(morphology_image != (Image *) NULL);
2592  assert(morphology_image->signature == MagickCoreSignature);
2593  assert(kernel != (KernelInfo *) NULL);
2594  assert(kernel->signature == MagickCoreSignature);
2595  assert(exception != (ExceptionInfo *) NULL);
2596  assert(exception->signature == MagickCoreSignature);
2597  status=MagickTrue;
2598  progress=0;
2599  image_view=AcquireVirtualCacheView(image,exception);
2600  morphology_view=AcquireAuthenticCacheView(morphology_image,exception);
2601  width=image->columns+kernel->width-1;
2602  offset.x=0;
2603  offset.y=0;
2604  switch (method)
2605  {
2606  case ConvolveMorphology:
2607  case DilateMorphology:
2610  {
2611  /*
2612  Kernel needs to used with reflection about origin.
2613  */
2614  offset.x=(ssize_t) kernel->width-kernel->x-1;
2615  offset.y=(ssize_t) kernel->height-kernel->y-1;
2616  break;
2617  }
2618  case ErodeMorphology:
2620  case HitAndMissMorphology:
2621  case ThinningMorphology:
2622  case ThickenMorphology:
2623  {
2624  offset.x=kernel->x;
2625  offset.y=kernel->y;
2626  break;
2627  }
2628  default:
2629  {
2630  assert("Not a Primitive Morphology Method" != (char *) NULL);
2631  break;
2632  }
2633  }
2634  changed=0;
2635  changes=(size_t *) AcquireQuantumMemory(GetOpenMPMaximumThreads(),
2636  sizeof(*changes));
2637  if (changes == (size_t *) NULL)
2638  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
2639  for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
2640  changes[j]=0;
2641 
2642  if ((method == ConvolveMorphology) && (kernel->width == 1))
2643  {
2644  register ssize_t
2645  x;
2646 
2647  /*
2648  Special handling (for speed) of vertical (blur) kernels. This performs
2649  its handling in columns rather than in rows. This is only done
2650  for convolve as it is the only method that generates very large 1-D
2651  vertical kernels (such as a 'BlurKernel')
2652  */
2653 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2654  #pragma omp parallel for schedule(static,4) shared(progress,status) \
2655  magick_number_threads(image,morphology_image,image->columns,1)
2656 #endif
2657  for (x=0; x < (ssize_t) image->columns; x++)
2658  {
2659  const int
2660  id = GetOpenMPThreadId();
2661 
2662  register const Quantum
2663  *magick_restrict p;
2664 
2665  register Quantum
2666  *magick_restrict q;
2667 
2668  register ssize_t
2669  r;
2670 
2671  ssize_t
2672  center;
2673 
2674  if (status == MagickFalse)
2675  continue;
2676  p=GetCacheViewVirtualPixels(image_view,x,-offset.y,1,image->rows+
2677  kernel->height-1,exception);
2678  q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
2679  morphology_image->rows,exception);
2680  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2681  {
2682  status=MagickFalse;
2683  continue;
2684  }
2685  center=(ssize_t) GetPixelChannels(image)*offset.y;
2686  for (r=0; r < (ssize_t) image->rows; r++)
2687  {
2688  register ssize_t
2689  i;
2690 
2691  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2692  {
2693  double
2694  alpha,
2695  gamma,
2696  pixel;
2697 
2698  PixelChannel
2699  channel;
2700 
2701  PixelTrait
2702  morphology_traits,
2703  traits;
2704 
2705  register const MagickRealType
2706  *magick_restrict k;
2707 
2708  register const Quantum
2709  *magick_restrict pixels;
2710 
2711  register ssize_t
2712  v;
2713 
2714  size_t
2715  count;
2716 
2717  channel=GetPixelChannelChannel(image,i);
2718  traits=GetPixelChannelTraits(image,channel);
2719  morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2720  if ((traits == UndefinedPixelTrait) ||
2721  (morphology_traits == UndefinedPixelTrait))
2722  continue;
2723  if (((traits & CopyPixelTrait) != 0) ||
2724  (GetPixelWriteMask(image,p+center) <= (QuantumRange/2)))
2725  {
2726  SetPixelChannel(morphology_image,channel,p[center+i],q);
2727  continue;
2728  }
2729  k=(&kernel->values[kernel->height-1]);
2730  pixels=p;
2731  pixel=bias;
2732  gamma=0.0;
2733  count=0;
2734  if ((morphology_traits & BlendPixelTrait) == 0)
2735  for (v=0; v < (ssize_t) kernel->height; v++)
2736  {
2737  if (!IsNaN(*k))
2738  {
2739  pixel+=(*k)*pixels[i];
2740  gamma+=(*k);
2741  count++;
2742  }
2743  k--;
2744  pixels+=GetPixelChannels(image);
2745  }
2746  else
2747  for (v=0; v < (ssize_t) kernel->height; v++)
2748  {
2749  if (!IsNaN(*k))
2750  {
2751  alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2752  pixel+=alpha*(*k)*pixels[i];
2753  gamma+=alpha*(*k);
2754  count++;
2755  }
2756  k--;
2757  pixels+=GetPixelChannels(image);
2758  }
2759  if (fabs(pixel-p[center+i]) > MagickEpsilon)
2760  changes[id]++;
2761  gamma=PerceptibleReciprocal(gamma);
2762  if (count != 0)
2763  gamma*=(double) kernel->height/count;
2764  SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*
2765  pixel),q);
2766  }
2767  p+=GetPixelChannels(image);
2768  q+=GetPixelChannels(morphology_image);
2769  }
2770  if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
2771  status=MagickFalse;
2772  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2773  {
2775  proceed;
2776 
2777 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2778  #pragma omp critical (MagickCore_MorphologyPrimitive)
2779 #endif
2780  proceed=SetImageProgress(image,MorphologyTag,progress++,
2781  image->rows);
2782  if (proceed == MagickFalse)
2783  status=MagickFalse;
2784  }
2785  }
2786  morphology_image->type=image->type;
2787  morphology_view=DestroyCacheView(morphology_view);
2788  image_view=DestroyCacheView(image_view);
2789  for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
2790  changed+=changes[j];
2791  changes=(size_t *) RelinquishMagickMemory(changes);
2792  return(status ? (ssize_t) changed : 0);
2793  }
2794  /*
2795  Normal handling of horizontal or rectangular kernels (row by row).
2796  */
2797 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2798  #pragma omp parallel for schedule(static,4) shared(progress,status) \
2799  magick_number_threads(image,morphology_image,image->rows,1)
2800 #endif
2801  for (y=0; y < (ssize_t) image->rows; y++)
2802  {
2803  const int
2804  id = GetOpenMPThreadId();
2805 
2806  register const Quantum
2807  *magick_restrict p;
2808 
2809  register Quantum
2810  *magick_restrict q;
2811 
2812  register ssize_t
2813  x;
2814 
2815  ssize_t
2816  center;
2817 
2818  if (status == MagickFalse)
2819  continue;
2820  p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,
2821  kernel->height,exception);
2822  q=GetCacheViewAuthenticPixels(morphology_view,0,y,morphology_image->columns,
2823  1,exception);
2824  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2825  {
2826  status=MagickFalse;
2827  continue;
2828  }
2829  center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
2830  GetPixelChannels(image)*offset.x);
2831  for (x=0; x < (ssize_t) image->columns; x++)
2832  {
2833  register ssize_t
2834  i;
2835 
2836  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2837  {
2838  double
2839  alpha,
2840  gamma,
2841  intensity,
2842  maximum,
2843  minimum,
2844  pixel;
2845 
2846  PixelChannel
2847  channel;
2848 
2849  PixelTrait
2850  morphology_traits,
2851  traits;
2852 
2853  register const MagickRealType
2854  *magick_restrict k;
2855 
2856  register const Quantum
2857  *magick_restrict pixels;
2858 
2859  register ssize_t
2860  u;
2861 
2862  size_t
2863  count;
2864 
2865  ssize_t
2866  v;
2867 
2868  channel=GetPixelChannelChannel(image,i);
2869  traits=GetPixelChannelTraits(image,channel);
2870  morphology_traits=GetPixelChannelTraits(morphology_image,channel);
2871  if ((traits == UndefinedPixelTrait) ||
2872  (morphology_traits == UndefinedPixelTrait))
2873  continue;
2874  if (((traits & CopyPixelTrait) != 0) ||
2875  (GetPixelWriteMask(image,p+center) <= (QuantumRange/2)))
2876  {
2877  SetPixelChannel(morphology_image,channel,p[center+i],q);
2878  continue;
2879  }
2880  pixels=p;
2881  maximum=0.0;
2882  minimum=(double) QuantumRange;
2883  switch (method)
2884  {
2885  case ConvolveMorphology: pixel=bias; break;
2886  case DilateMorphology:
2888  {
2889  pixel=0.0;
2890  break;
2891  }
2892  default:
2893  {
2894  pixel=(double) p[center+i];
2895  break;
2896  }
2897  }
2898  count=0;
2899  gamma=1.0;
2900  switch (method)
2901  {
2902  case ConvolveMorphology:
2903  {
2904  /*
2905  Weighted Average of pixels using reflected kernel
2906 
2907  For correct working of this operation for asymetrical kernels,
2908  the kernel needs to be applied in its reflected form. That is
2909  its values needs to be reversed.
2910 
2911  Correlation is actually the same as this but without reflecting
2912  the kernel, and thus 'lower-level' that Convolution. However as
2913  Convolution is the more common method used, and it does not
2914  really cost us much in terms of processing to use a reflected
2915  kernel, so it is Convolution that is implemented.
2916 
2917  Correlation will have its kernel reflected before calling this
2918  function to do a Convolve.
2919 
2920  For more details of Correlation vs Convolution see
2921  http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
2922  */
2923  k=(&kernel->values[kernel->width*kernel->height-1]);
2924  if ((morphology_traits & BlendPixelTrait) == 0)
2925  {
2926  /*
2927  No alpha blending.
2928  */
2929  for (v=0; v < (ssize_t) kernel->height; v++)
2930  {
2931  for (u=0; u < (ssize_t) kernel->width; u++)
2932  {
2933  if (!IsNaN(*k))
2934  {
2935  pixel+=(*k)*pixels[i];
2936  count++;
2937  }
2938  k--;
2939  pixels+=GetPixelChannels(image);
2940  }
2941  pixels+=(image->columns-1)*GetPixelChannels(image);
2942  }
2943  break;
2944  }
2945  /*
2946  Alpha blending.
2947  */
2948  gamma=0.0;
2949  for (v=0; v < (ssize_t) kernel->height; v++)
2950  {
2951  for (u=0; u < (ssize_t) kernel->width; u++)
2952  {
2953  if (!IsNaN(*k))
2954  {
2955  alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2956  pixel+=alpha*(*k)*pixels[i];
2957  gamma+=alpha*(*k);
2958  count++;
2959  }
2960  k--;
2961  pixels+=GetPixelChannels(image);
2962  }
2963  pixels+=(image->columns-1)*GetPixelChannels(image);
2964  }
2965  break;
2966  }
2967  case ErodeMorphology:
2968  {
2969  /*
2970  Minimum value within kernel neighbourhood.
2971 
2972  The kernel is not reflected for this operation. In normal
2973  Greyscale Morphology, the kernel value should be added
2974  to the real value, this is currently not done, due to the
2975  nature of the boolean kernels being used.
2976  */
2977  k=kernel->values;
2978  for (v=0; v < (ssize_t) kernel->height; v++)
2979  {
2980  for (u=0; u < (ssize_t) kernel->width; u++)
2981  {
2982  if (!IsNaN(*k) && (*k >= 0.5))
2983  {
2984  if ((double) pixels[i] < pixel)
2985  pixel=(double) pixels[i];
2986  }
2987  k++;
2988  pixels+=GetPixelChannels(image);
2989  }
2990  pixels+=(image->columns-1)*GetPixelChannels(image);
2991  }
2992  break;
2993  }
2994  case DilateMorphology:
2995  {
2996  /*
2997  Maximum value within kernel neighbourhood.
2998 
2999  For correct working of this operation for asymetrical kernels,
3000  the kernel needs to be applied in its reflected form. That is
3001  its values needs to be reversed.
3002 
3003  In normal Greyscale Morphology, the kernel value should be
3004  added to the real value, this is currently not done, due to the
3005  nature of the boolean kernels being used.
3006  */
3007  k=(&kernel->values[kernel->width*kernel->height-1]);
3008  for (v=0; v < (ssize_t) kernel->height; v++)
3009  {
3010  for (u=0; u < (ssize_t) kernel->width; u++)
3011  {
3012  if (!IsNaN(*k) && (*k > 0.5))
3013  {
3014  if ((double) pixels[i] > pixel)
3015  pixel=(double) pixels[i];
3016  }
3017  k--;
3018  pixels+=GetPixelChannels(image);
3019  }
3020  pixels+=(image->columns-1)*GetPixelChannels(image);
3021  }
3022  break;
3023  }
3024  case HitAndMissMorphology:
3025  case ThinningMorphology:
3026  case ThickenMorphology:
3027  {
3028  /*
3029  Minimum of foreground pixel minus maxumum of background pixels.
3030 
3031  The kernel is not reflected for this operation, and consists
3032  of both foreground and background pixel neighbourhoods, 0.0 for
3033  background, and 1.0 for foreground with either Nan or 0.5 values
3034  for don't care.
3035 
3036  This never produces a meaningless negative result. Such results
3037  cause Thinning/Thicken to not work correctly when used against a
3038  greyscale image.
3039  */
3040  k=kernel->values;
3041  for (v=0; v < (ssize_t) kernel->height; v++)
3042  {
3043  for (u=0; u < (ssize_t) kernel->width; u++)
3044  {
3045  if (!IsNaN(*k))
3046  {
3047  if (*k > 0.7)
3048  {
3049  if ((double) pixels[i] < pixel)
3050  pixel=(double) pixels[i];
3051  }
3052  else
3053  if (*k < 0.3)
3054  {
3055  if ((double) pixels[i] > maximum)
3056  maximum=(double) pixels[i];
3057  }
3058  count++;
3059  }
3060  k++;
3061  pixels+=GetPixelChannels(image);
3062  }
3063  pixels+=(image->columns-1)*GetPixelChannels(image);
3064  }
3065  pixel-=maximum;
3066  if (pixel < 0.0)
3067  pixel=0.0;
3068  if (method == ThinningMorphology)
3069  pixel=(double) p[center+i]-pixel;
3070  else
3071  if (method == ThickenMorphology)
3072  pixel+=(double) p[center+i]+pixel;
3073  break;
3074  }
3076  {
3077  /*
3078  Select pixel with minimum intensity within kernel neighbourhood.
3079 
3080  The kernel is not reflected for this operation.
3081  */
3082  k=kernel->values;
3083  for (v=0; v < (ssize_t) kernel->height; v++)
3084  {
3085  for (u=0; u < (ssize_t) kernel->width; u++)
3086  {
3087  if (!IsNaN(*k) && (*k >= 0.5))
3088  {
3089  intensity=(double) GetPixelIntensity(image,pixels);
3090  if (intensity < minimum)
3091  {
3092  pixel=(double) pixels[i];
3093  minimum=intensity;
3094  }
3095  count++;
3096  }
3097  k++;
3098  pixels+=GetPixelChannels(image);
3099  }
3100  pixels+=(image->columns-1)*GetPixelChannels(image);
3101  }
3102  break;
3103  }
3105  {
3106  /*
3107  Select pixel with maximum intensity within kernel neighbourhood.
3108 
3109  The kernel is not reflected for this operation.
3110  */
3111  k=(&kernel->values[kernel->width*kernel->height-1]);
3112  for (v=0; v < (ssize_t) kernel->height; v++)
3113  {
3114  for (u=0; u < (ssize_t) kernel->width; u++)
3115  {
3116  if (!IsNaN(*k) && (*k >= 0.5))
3117  {
3118  intensity=(double) GetPixelIntensity(image,pixels);
3119  if (intensity > maximum)
3120  {
3121  pixel=(double) pixels[i];
3122  maximum=intensity;
3123  }
3124  count++;
3125  }
3126  k--;
3127  pixels+=GetPixelChannels(image);
3128  }
3129  pixels+=(image->columns-1)*GetPixelChannels(image);
3130  }
3131  break;
3132  }
3134  {
3135  /*
3136  Compute th iterative distance from black edge of a white image
3137  shape. Essentually white values are decreased to the smallest
3138  'distance from edge' it can find.
3139 
3140  It works by adding kernel values to the neighbourhood, and and
3141  select the minimum value found. The kernel is rotated before
3142  use, so kernel distances match resulting distances, when a user
3143  provided asymmetric kernel is applied.
3144 
3145  This code is nearly identical to True GrayScale Morphology but
3146  not quite.
3147 
3148  GreyDilate Kernel values added, maximum value found Kernel is
3149  rotated before use.
3150 
3151  GrayErode: Kernel values subtracted and minimum value found No
3152  kernel rotation used.
3153 
3154  Note the the Iterative Distance method is essentially a
3155  GrayErode, but with negative kernel values, and kernel rotation
3156  applied.
3157  */
3158  k=(&kernel->values[kernel->width*kernel->height-1]);
3159  for (v=0; v < (ssize_t) kernel->height; v++)
3160  {
3161  for (u=0; u < (ssize_t) kernel->width; u++)
3162  {
3163  if (!IsNaN(*k))
3164  {
3165  if ((pixels[i]+(*k)) < pixel)
3166  pixel=(double) pixels[i]+(*k);
3167  count++;
3168  }
3169  k--;
3170  pixels+=GetPixelChannels(image);
3171  }
3172  pixels+=(image->columns-1)*GetPixelChannels(image);
3173  }
3174  break;
3175  }
3176  case UndefinedMorphology:
3177  default:
3178  break;
3179  }
3180  if (fabs(pixel-p[center+i]) > MagickEpsilon)
3181  changes[id]++;
3182  gamma=PerceptibleReciprocal(gamma);
3183  if (count != 0)
3184  gamma*=(double) kernel->height*kernel->width/count;
3185  SetPixelChannel(morphology_image,channel,ClampToQuantum(gamma*pixel),q);
3186  }
3187  p+=GetPixelChannels(image);
3188  q+=GetPixelChannels(morphology_image);
3189  }
3190  if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3191  status=MagickFalse;
3192  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3193  {
3195  proceed;
3196 
3197 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3198  #pragma omp critical (MagickCore_MorphologyPrimitive)
3199 #endif
3200  proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
3201  if (proceed == MagickFalse)
3202  status=MagickFalse;
3203  }
3204  }
3205  morphology_view=DestroyCacheView(morphology_view);
3206  image_view=DestroyCacheView(image_view);
3207  for (j=0; j < (ssize_t) GetOpenMPMaximumThreads(); j++)
3208  changed+=changes[j];
3209  changes=(size_t *) RelinquishMagickMemory(changes);
3210  return(status ? (ssize_t) changed : -1);
3211 }
3212 
3213 /*
3214  This is almost identical to the MorphologyPrimative() function above, but
3215  applies the primitive directly to the actual image using two passes, once in
3216  each direction, with the results of the previous (and current) row being
3217  re-used.
3218 
3219  That is after each row is 'Sync'ed' into the image, the next row makes use of
3220  those values as part of the calculation of the next row. It repeats, but
3221  going in the oppisite (bottom-up) direction.
3222 
3223  Because of this 're-use of results' this function can not make use of multi-
3224  threaded, parellel processing.
3225 */
3226 static ssize_t MorphologyPrimitiveDirect(Image *image,
3227  const MorphologyMethod method,const KernelInfo *kernel,
3228  ExceptionInfo *exception)
3229 {
3230  CacheView
3231  *morphology_view,
3232  *image_view;
3233 
3235  status;
3236 
3238  progress;
3239 
3240  OffsetInfo
3241  offset;
3242 
3243  size_t
3244  width,
3245  changed;
3246 
3247  ssize_t
3248  y;
3249 
3250  assert(image != (Image *) NULL);
3251  assert(image->signature == MagickCoreSignature);
3252  assert(kernel != (KernelInfo *) NULL);
3253  assert(kernel->signature == MagickCoreSignature);
3254  assert(exception != (ExceptionInfo *) NULL);
3255  assert(exception->signature == MagickCoreSignature);
3256  status=MagickTrue;
3257  changed=0;
3258  progress=0;
3259  switch(method)
3260  {
3261  case DistanceMorphology:
3262  case VoronoiMorphology:
3263  {
3264  /*
3265  Kernel reflected about origin.
3266  */
3267  offset.x=(ssize_t) kernel->width-kernel->x-1;
3268  offset.y=(ssize_t) kernel->height-kernel->y-1;
3269  break;
3270  }
3271  default:
3272  {
3273  offset.x=kernel->x;
3274  offset.y=kernel->y;
3275  break;
3276  }
3277  }
3278  /*
3279  Two views into same image, do not thread.
3280  */
3281  image_view=AcquireVirtualCacheView(image,exception);
3282  morphology_view=AcquireAuthenticCacheView(image,exception);
3283  width=image->columns+kernel->width-1;
3284  for (y=0; y < (ssize_t) image->rows; y++)
3285  {
3286  register const Quantum
3287  *magick_restrict p;
3288 
3289  register Quantum
3290  *magick_restrict q;
3291 
3292  register ssize_t
3293  x;
3294 
3295  ssize_t
3296  center;
3297 
3298  /*
3299  Read virtual pixels, and authentic pixels, from the same image! We read
3300  using virtual to get virtual pixel handling, but write back into the same
3301  image.
3302 
3303  Only top half of kernel is processed as we do a single pass downward
3304  through the image iterating the distance function as we go.
3305  */
3306  if (status == MagickFalse)
3307  continue;
3308  p=GetCacheViewVirtualPixels(image_view,-offset.x,y-offset.y,width,(size_t)
3309  offset.y+1,exception);
3310  q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3311  exception);
3312  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3313  {
3314  status=MagickFalse;
3315  continue;
3316  }
3317  center=(ssize_t) (GetPixelChannels(image)*width*offset.y+
3318  GetPixelChannels(image)*offset.x);
3319  for (x=0; x < (ssize_t) image->columns; x++)
3320  {
3321  register ssize_t
3322  i;
3323 
3324  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3325  {
3326  double
3327  pixel;
3328 
3329  PixelTrait
3330  traits;
3331 
3332  register const MagickRealType
3333  *magick_restrict k;
3334 
3335  register const Quantum
3336  *magick_restrict pixels;
3337 
3338  register ssize_t
3339  u;
3340 
3341  ssize_t
3342  v;
3343 
3344  traits=GetPixelChannelTraits(image,(PixelChannel) i);
3345  if (traits == UndefinedPixelTrait)
3346  continue;
3347  if (((traits & CopyPixelTrait) != 0) ||
3348  (GetPixelWriteMask(image,p+center) <= (QuantumRange/2)))
3349  continue;
3350  pixels=p;
3351  pixel=(double) QuantumRange;
3352  switch (method)
3353  {
3354  case DistanceMorphology:
3355  {
3356  k=(&kernel->values[kernel->width*kernel->height-1]);
3357  for (v=0; v <= offset.y; v++)
3358  {
3359  for (u=0; u < (ssize_t) kernel->width; u++)
3360  {
3361  if (!IsNaN(*k))
3362  {
3363  if ((pixels[i]+(*k)) < pixel)
3364  pixel=(double) pixels[i]+(*k);
3365  }
3366  k--;
3367  pixels+=GetPixelChannels(image);
3368  }
3369  pixels+=(image->columns-1)*GetPixelChannels(image);
3370  }
3371  k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3372  pixels=q-offset.x*GetPixelChannels(image);
3373  for (u=0; u < offset.x; u++)
3374  {
3375  if (!IsNaN(*k) && ((x+u-offset.x) >= 0))
3376  {
3377  if ((pixels[i]+(*k)) < pixel)
3378  pixel=(double) pixels[i]+(*k);
3379  }
3380  k--;
3381  pixels+=GetPixelChannels(image);
3382  }
3383  break;
3384  }
3385  case VoronoiMorphology:
3386  {
3387  k=(&kernel->values[kernel->width*kernel->height-1]);
3388  for (v=0; v < offset.y; v++)
3389  {
3390  for (u=0; u < (ssize_t) kernel->width; u++)
3391  {
3392  if (!IsNaN(*k))
3393  {
3394  if ((pixels[i]+(*k)) < pixel)
3395  pixel=(double) pixels[i]+(*k);
3396  }
3397  k--;
3398  pixels+=GetPixelChannels(image);
3399  }
3400  pixels+=(image->columns-1)*GetPixelChannels(image);
3401  }
3402  k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3403  pixels=q-offset.x*GetPixelChannels(image);
3404  for (u=0; u < offset.x; u++)
3405  {
3406  if (!IsNaN(*k) && ((x+u-offset.x) >= 0))
3407  {
3408  if ((pixels[i]+(*k)) < pixel)
3409  pixel=(double) pixels[i]+(*k);
3410  }
3411  k--;
3412  pixels+=GetPixelChannels(image);
3413  }
3414  break;
3415  }
3416  default:
3417  break;
3418  }
3419  if (fabs(pixel-q[i]) > MagickEpsilon)
3420  changed++;
3421  q[i]=ClampToQuantum(pixel);
3422  }
3423  p+=GetPixelChannels(image);
3424  q+=GetPixelChannels(image);
3425  }
3426  if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3427  status=MagickFalse;
3428  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3429  {
3431  proceed;
3432 
3433  proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3434  if (proceed == MagickFalse)
3435  status=MagickFalse;
3436  }
3437  }
3438  morphology_view=DestroyCacheView(morphology_view);
3439  image_view=DestroyCacheView(image_view);
3440  /*
3441  Do the reverse pass through the image.
3442  */
3443  image_view=AcquireVirtualCacheView(image,exception);
3444  morphology_view=AcquireAuthenticCacheView(image,exception);
3445  for (y=(ssize_t) image->rows-1; y >= 0; y--)
3446  {
3447  register const Quantum
3448  *magick_restrict p;
3449 
3450  register Quantum
3451  *magick_restrict q;
3452 
3453  register ssize_t
3454  x;
3455 
3456  ssize_t
3457  center;
3458 
3459  /*
3460  Read virtual pixels, and authentic pixels, from the same image. We
3461  read using virtual to get virtual pixel handling, but write back
3462  into the same image.
3463 
3464  Only the bottom half of the kernel is processed as we up the image.
3465  */
3466  if (status == MagickFalse)
3467  continue;
3468  p=GetCacheViewVirtualPixels(image_view,-offset.x,y,width,(size_t)
3469  kernel->y+1,exception);
3470  q=GetCacheViewAuthenticPixels(morphology_view,0,y,image->columns,1,
3471  exception);
3472  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3473  {
3474  status=MagickFalse;
3475  continue;
3476  }
3477  p+=(image->columns-1)*GetPixelChannels(image);
3478  q+=(image->columns-1)*GetPixelChannels(image);
3479  center=(ssize_t) (offset.x*GetPixelChannels(image));
3480  for (x=(ssize_t) image->columns-1; x >= 0; x--)
3481  {
3482  register ssize_t
3483  i;
3484 
3485  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3486  {
3487  double
3488  pixel;
3489 
3490  PixelTrait
3491  traits;
3492 
3493  register const MagickRealType
3494  *magick_restrict k;
3495 
3496  register const Quantum
3497  *magick_restrict pixels;
3498 
3499  register ssize_t
3500  u;
3501 
3502  ssize_t
3503  v;
3504 
3505  traits=GetPixelChannelTraits(image,(PixelChannel) i);
3506  if (traits == UndefinedPixelTrait)
3507  continue;
3508  if (((traits & CopyPixelTrait) != 0) ||
3509  (GetPixelWriteMask(image,p+center) <= (QuantumRange/2)))
3510  continue;
3511  pixels=p;
3512  pixel=(double) QuantumRange;
3513  switch (method)
3514  {
3515  case DistanceMorphology:
3516  {
3517  k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3518  for (v=offset.y; v < (ssize_t) kernel->height; v++)
3519  {
3520  for (u=0; u < (ssize_t) kernel->width; u++)
3521  {
3522  if (!IsNaN(*k))
3523  {
3524  if ((pixels[i]+(*k)) < pixel)
3525  pixel=(double) pixels[i]+(*k);
3526  }
3527  k--;
3528  pixels+=GetPixelChannels(image);
3529  }
3530  pixels+=(image->columns-1)*GetPixelChannels(image);
3531  }
3532  k=(&kernel->values[kernel->width*kernel->y+kernel->x-1]);
3533  pixels=q;
3534  for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3535  {
3536  pixels+=GetPixelChannels(image);
3537  if (!IsNaN(*k) && ((x+u-offset.x) < (ssize_t) image->columns))
3538  {
3539  if ((pixels[i]+(*k)) < pixel)
3540  pixel=(double) pixels[i]+(*k);
3541  }
3542  k--;
3543  }
3544  break;
3545  }
3546  case VoronoiMorphology:
3547  {
3548  k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3549  for (v=offset.y; v < (ssize_t) kernel->height; v++)
3550  {
3551  for (u=0; u < (ssize_t) kernel->width; u++)
3552  {
3553  if (!IsNaN(*k))
3554  {
3555  if ((pixels[i]+(*k)) < pixel)
3556  pixel=(double) pixels[i]+(*k);
3557  }
3558  k--;
3559  pixels+=GetPixelChannels(image);
3560  }
3561  pixels+=(image->columns-1)*GetPixelChannels(image);
3562  }
3563  k=(&kernel->values[kernel->width*(kernel->y+1)-1]);
3564  pixels=q;
3565  for (u=offset.x+1; u < (ssize_t) kernel->width; u++)
3566  {
3567  pixels+=GetPixelChannels(image);
3568  if (!IsNaN(*k) && ((x+u-offset.x) < (ssize_t) image->columns))
3569  {
3570  if ((pixels[i]+(*k)) < pixel)
3571  pixel=(double) pixels[i]+(*k);
3572  }
3573  k--;
3574  }
3575  break;
3576  }
3577  default:
3578  break;
3579  }
3580  if (fabs(pixel-q[i]) > MagickEpsilon)
3581  changed++;
3582  q[i]=ClampToQuantum(pixel);
3583  }
3584  p-=GetPixelChannels(image);
3585  q-=GetPixelChannels(image);
3586  }
3587  if (SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
3588  status=MagickFalse;
3589  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3590  {
3592  proceed;
3593 
3594  proceed=SetImageProgress(image,MorphologyTag,progress++,2*image->rows);
3595  if (proceed == MagickFalse)
3596  status=MagickFalse;
3597  }
3598  }
3599  morphology_view=DestroyCacheView(morphology_view);
3600  image_view=DestroyCacheView(image_view);
3601  return(status ? (ssize_t) changed : -1);
3602 }
3603 
3604 /*
3605  Apply a Morphology by calling one of the above low level primitive
3606  application functions. This function handles any iteration loops,
3607  composition or re-iteration of results, and compound morphology methods that
3608  is based on multiple low-level (staged) morphology methods.
3609 
3610  Basically this provides the complex glue between the requested morphology
3611  method and raw low-level implementation (above).
3612 */
3614  const MorphologyMethod method, const ssize_t iterations,
3615  const KernelInfo *kernel, const CompositeOperator compose,const double bias,
3616  ExceptionInfo *exception)
3617 {
3619  curr_compose;
3620 
3621  Image
3622  *curr_image, /* Image we are working with or iterating */
3623  *work_image, /* secondary image for primitive iteration */
3624  *save_image, /* saved image - for 'edge' method only */
3625  *rslt_image; /* resultant image - after multi-kernel handling */
3626 
3627  KernelInfo
3628  *reflected_kernel, /* A reflected copy of the kernel (if needed) */
3629  *norm_kernel, /* the current normal un-reflected kernel */
3630  *rflt_kernel, /* the current reflected kernel (if needed) */
3631  *this_kernel; /* the kernel being applied */
3632 
3634  primitive; /* the current morphology primitive being applied */
3635 
3637  rslt_compose; /* multi-kernel compose method for results to use */
3638 
3640  special, /* do we use a direct modify function? */
3641  verbose; /* verbose output of results */
3642 
3643  size_t
3644  method_loop, /* Loop 1: number of compound method iterations (norm 1) */
3645  method_limit, /* maximum number of compound method iterations */
3646  kernel_number, /* Loop 2: the kernel number being applied */
3647  stage_loop, /* Loop 3: primitive loop for compound morphology */
3648  stage_limit, /* how many primitives are in this compound */
3649  kernel_loop, /* Loop 4: iterate the kernel over image */
3650  kernel_limit, /* number of times to iterate kernel */
3651  count, /* total count of primitive steps applied */
3652  kernel_changed, /* total count of changed using iterated kernel */
3653  method_changed; /* total count of changed over method iteration */
3654 
3655  ssize_t
3656  changed; /* number pixels changed by last primitive operation */
3657 
3658  char
3659  v_info[MagickPathExtent];
3660 
3661  assert(image != (Image *) NULL);
3662  assert(image->signature == MagickCoreSignature);
3663  assert(kernel != (KernelInfo *) NULL);
3664  assert(kernel->signature == MagickCoreSignature);
3665  assert(exception != (ExceptionInfo *) NULL);
3666  assert(exception->signature == MagickCoreSignature);
3667 
3668  count = 0; /* number of low-level morphology primitives performed */
3669  if ( iterations == 0 )
3670  return((Image *) NULL); /* null operation - nothing to do! */
3671 
3672  kernel_limit = (size_t) iterations;
3673  if ( iterations < 0 ) /* negative interations = infinite (well alomst) */
3674  kernel_limit = image->columns>image->rows ? image->columns : image->rows;
3675 
3676  verbose = IsStringTrue(GetImageArtifact(image,"debug"));
3677 
3678  /* initialise for cleanup */
3679  curr_image = (Image *) image;
3680  curr_compose = image->compose;
3681  (void) curr_compose;
3682  work_image = save_image = rslt_image = (Image *) NULL;
3683  reflected_kernel = (KernelInfo *) NULL;
3684 
3685  /* Initialize specific methods
3686  * + which loop should use the given iteratations
3687  * + how many primitives make up the compound morphology
3688  * + multi-kernel compose method to use (by default)
3689  */
3690  method_limit = 1; /* just do method once, unless otherwise set */
3691  stage_limit = 1; /* assume method is not a compound */
3692  special = MagickFalse; /* assume it is NOT a direct modify primitive */
3693  rslt_compose = compose; /* and we are composing multi-kernels as given */
3694  switch( method ) {
3695  case SmoothMorphology: /* 4 primitive compound morphology */
3696  stage_limit = 4;
3697  break;
3698  case OpenMorphology: /* 2 primitive compound morphology */
3700  case TopHatMorphology:
3701  case CloseMorphology:
3703  case BottomHatMorphology:
3704  case EdgeMorphology:
3705  stage_limit = 2;
3706  break;
3707  case HitAndMissMorphology:
3708  rslt_compose = LightenCompositeOp; /* Union of multi-kernel results */
3709  /* FALL THUR */
3710  case ThinningMorphology:
3711  case ThickenMorphology:
3712  method_limit = kernel_limit; /* iterate the whole method */
3713  kernel_limit = 1; /* do not do kernel iteration */
3714  break;
3715  case DistanceMorphology:
3716  case VoronoiMorphology:
3717  special = MagickTrue; /* use special direct primative */
3718  break;
3719  default:
3720  break;
3721  }
3722 
3723  /* Apply special methods with special requirments
3724  ** For example, single run only, or post-processing requirements
3725  */
3726  if ( special != MagickFalse )
3727  {
3728  rslt_image=CloneImage(image,0,0,MagickTrue,exception);
3729  if (rslt_image == (Image *) NULL)
3730  goto error_cleanup;
3731  if (SetImageStorageClass(rslt_image,DirectClass,exception) == MagickFalse)
3732  goto error_cleanup;
3733 
3734  changed=MorphologyPrimitiveDirect(rslt_image,method,kernel,exception);
3735 
3736  if (verbose != MagickFalse)
3737  (void) (void) FormatLocaleFile(stderr,
3738  "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
3740  1.0,0.0,1.0, (double) changed);
3741 
3742  if ( changed < 0 )
3743  goto error_cleanup;
3744 
3745  if ( method == VoronoiMorphology ) {
3746  /* Preserve the alpha channel of input image - but turned it off */
3747  (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3748  exception);
3749  (void) CompositeImage(rslt_image,image,CopyAlphaCompositeOp,
3750  MagickTrue,0,0,exception);
3751  (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel,
3752  exception);
3753  }
3754  goto exit_cleanup;
3755  }
3756 
3757  /* Handle user (caller) specified multi-kernel composition method */
3758  if ( compose != UndefinedCompositeOp )
3759  rslt_compose = compose; /* override default composition for method */
3760  if ( rslt_compose == UndefinedCompositeOp )
3761  rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
3762 
3763  /* Some methods require a reflected kernel to use with primitives.
3764  * Create the reflected kernel for those methods. */
3765  switch ( method ) {
3766  case CorrelateMorphology:
3767  case CloseMorphology:
3769  case BottomHatMorphology:
3770  case SmoothMorphology:
3771  reflected_kernel = CloneKernelInfo(kernel);
3772  if (reflected_kernel == (KernelInfo *) NULL)
3773  goto error_cleanup;
3774  RotateKernelInfo(reflected_kernel,180);
3775  break;
3776  default:
3777  break;
3778  }
3779 
3780  /* Loops around more primitive morpholgy methods
3781  ** erose, dilate, open, close, smooth, edge, etc...
3782  */
3783  /* Loop 1: iterate the compound method */
3784  method_loop = 0;
3785  method_changed = 1;
3786  while ( method_loop < method_limit && method_changed > 0 ) {
3787  method_loop++;
3788  method_changed = 0;
3789 
3790  /* Loop 2: iterate over each kernel in a multi-kernel list */
3791  norm_kernel = (KernelInfo *) kernel;
3792  this_kernel = (KernelInfo *) kernel;
3793  rflt_kernel = reflected_kernel;
3794 
3795  kernel_number = 0;
3796  while ( norm_kernel != NULL ) {
3797 
3798  /* Loop 3: Compound Morphology Staging - Select Primative to apply */
3799  stage_loop = 0; /* the compound morphology stage number */
3800  while ( stage_loop < stage_limit ) {
3801  stage_loop++; /* The stage of the compound morphology */
3802 
3803  /* Select primitive morphology for this stage of compound method */
3804  this_kernel = norm_kernel; /* default use unreflected kernel */
3805  primitive = method; /* Assume method is a primitive */
3806  switch( method ) {
3807  case ErodeMorphology: /* just erode */
3808  case EdgeInMorphology: /* erode and image difference */
3809  primitive = ErodeMorphology;
3810  break;
3811  case DilateMorphology: /* just dilate */
3812  case EdgeOutMorphology: /* dilate and image difference */
3813  primitive = DilateMorphology;
3814  break;
3815  case OpenMorphology: /* erode then dialate */
3816  case TopHatMorphology: /* open and image difference */
3817  primitive = ErodeMorphology;
3818  if ( stage_loop == 2 )
3819  primitive = DilateMorphology;
3820  break;
3822  primitive = ErodeIntensityMorphology;
3823  if ( stage_loop == 2 )
3824  primitive = DilateIntensityMorphology;
3825  break;
3826  case CloseMorphology: /* dilate, then erode */
3827  case BottomHatMorphology: /* close and image difference */
3828  this_kernel = rflt_kernel; /* use the reflected kernel */
3829  primitive = DilateMorphology;
3830  if ( stage_loop == 2 )
3831  primitive = ErodeMorphology;
3832  break;
3834  this_kernel = rflt_kernel; /* use the reflected kernel */
3835  primitive = DilateIntensityMorphology;
3836  if ( stage_loop == 2 )
3837  primitive = ErodeIntensityMorphology;
3838  break;
3839  case SmoothMorphology: /* open, close */
3840  switch ( stage_loop ) {
3841  case 1: /* start an open method, which starts with Erode */
3842  primitive = ErodeMorphology;
3843  break;
3844  case 2: /* now Dilate the Erode */
3845  primitive = DilateMorphology;
3846  break;
3847  case 3: /* Reflect kernel a close */
3848  this_kernel = rflt_kernel; /* use the reflected kernel */
3849  primitive = DilateMorphology;
3850  break;
3851  case 4: /* Finish the Close */
3852  this_kernel = rflt_kernel; /* use the reflected kernel */
3853  primitive = ErodeMorphology;
3854  break;
3855  }
3856  break;
3857  case EdgeMorphology: /* dilate and erode difference */
3858  primitive = DilateMorphology;
3859  if ( stage_loop == 2 ) {
3860  save_image = curr_image; /* save the image difference */
3861  curr_image = (Image *) image;
3862  primitive = ErodeMorphology;
3863  }
3864  break;
3865  case CorrelateMorphology:
3866  /* A Correlation is a Convolution with a reflected kernel.
3867  ** However a Convolution is a weighted sum using a reflected
3868  ** kernel. It may seem stange to convert a Correlation into a
3869  ** Convolution as the Correlation is the simplier method, but
3870  ** Convolution is much more commonly used, and it makes sense to
3871  ** implement it directly so as to avoid the need to duplicate the
3872  ** kernel when it is not required (which is typically the
3873  ** default).
3874  */
3875  this_kernel = rflt_kernel; /* use the reflected kernel */
3876  primitive = ConvolveMorphology;
3877  break;
3878  default:
3879  break;
3880  }
3881  assert( this_kernel != (KernelInfo *) NULL );
3882 
3883  /* Extra information for debugging compound operations */
3884  if (verbose != MagickFalse) {
3885  if ( stage_limit > 1 )
3886  (void) FormatLocaleString(v_info,MagickPathExtent,"%s:%.20g.%.20g -> ",
3888  method_loop,(double) stage_loop);
3889  else if ( primitive != method )
3890  (void) FormatLocaleString(v_info, MagickPathExtent, "%s:%.20g -> ",
3892  method_loop);
3893  else
3894  v_info[0] = '\0';
3895  }
3896 
3897  /* Loop 4: Iterate the kernel with primitive */
3898  kernel_loop = 0;
3899  kernel_changed = 0;
3900  changed = 1;
3901  while ( kernel_loop < kernel_limit && changed > 0 ) {
3902  kernel_loop++; /* the iteration of this kernel */
3903 
3904  /* Create a clone as the destination image, if not yet defined */
3905  if ( work_image == (Image *) NULL )
3906  {
3907  work_image=CloneImage(image,0,0,MagickTrue,exception);
3908  if (work_image == (Image *) NULL)
3909  goto error_cleanup;
3910  if (SetImageStorageClass(work_image,DirectClass,exception) == MagickFalse)
3911  goto error_cleanup;
3912  }
3913 
3914  /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
3915  count++;
3916  changed = MorphologyPrimitive(curr_image, work_image, primitive,
3917  this_kernel, bias, exception);
3918  if (verbose != MagickFalse) {
3919  if ( kernel_loop > 1 )
3920  (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line from previous */
3921  (void) (void) FormatLocaleFile(stderr,
3922  "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
3924  primitive),(this_kernel == rflt_kernel ) ? "*" : "",
3925  (double) (method_loop+kernel_loop-1),(double) kernel_number,
3926  (double) count,(double) changed);
3927  }
3928  if ( changed < 0 )
3929  goto error_cleanup;
3930  kernel_changed += changed;
3931  method_changed += changed;
3932 
3933  /* prepare next loop */
3934  { Image *tmp = work_image; /* swap images for iteration */
3935  work_image = curr_image;
3936  curr_image = tmp;
3937  }
3938  if ( work_image == image )
3939  work_image = (Image *) NULL; /* replace input 'image' */
3940 
3941  } /* End Loop 4: Iterate the kernel with primitive */
3942 
3943  if (verbose != MagickFalse && kernel_changed != (size_t)changed)
3944  (void) FormatLocaleFile(stderr, " Total %.20g",(double) kernel_changed);
3945  if (verbose != MagickFalse && stage_loop < stage_limit)
3946  (void) FormatLocaleFile(stderr, "\n"); /* add end-of-line before looping */
3947 
3948 #if 0
3949  (void) FormatLocaleFile(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
3950  (void) FormatLocaleFile(stderr, " curr =0x%lx\n", (unsigned long)curr_image);
3951  (void) FormatLocaleFile(stderr, " work =0x%lx\n", (unsigned long)work_image);
3952  (void) FormatLocaleFile(stderr, " save =0x%lx\n", (unsigned long)save_image);
3953  (void) FormatLocaleFile(stderr, " union=0x%lx\n", (unsigned long)rslt_image);
3954 #endif
3955 
3956  } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
3957 
3958  /* Final Post-processing for some Compound Methods
3959  **
3960  ** The removal of any 'Sync' channel flag in the Image Compositon
3961  ** below ensures the methematical compose method is applied in a
3962  ** purely mathematical way, and only to the selected channels.
3963  ** Turn off SVG composition 'alpha blending'.
3964  */
3965  switch( method ) {
3966  case EdgeOutMorphology:
3967  case EdgeInMorphology:
3968  case TopHatMorphology:
3969  case BottomHatMorphology:
3970  if (verbose != MagickFalse)
3971  (void) FormatLocaleFile(stderr,
3972  "\n%s: Difference with original image",CommandOptionToMnemonic(
3973  MagickMorphologyOptions, method) );
3974  (void) CompositeImage(curr_image,image,DifferenceCompositeOp,
3975  MagickTrue,0,0,exception);
3976  break;
3977  case EdgeMorphology:
3978  if (verbose != MagickFalse)
3979  (void) FormatLocaleFile(stderr,
3980  "\n%s: Difference of Dilate and Erode",CommandOptionToMnemonic(
3981  MagickMorphologyOptions, method) );
3982  (void) CompositeImage(curr_image,save_image,DifferenceCompositeOp,
3983  MagickTrue,0,0,exception);
3984  save_image = DestroyImage(save_image); /* finished with save image */
3985  break;
3986  default:
3987  break;
3988  }
3989 
3990  /* multi-kernel handling: re-iterate, or compose results */
3991  if ( kernel->next == (KernelInfo *) NULL )
3992  rslt_image = curr_image; /* just return the resulting image */
3993  else if ( rslt_compose == NoCompositeOp )
3994  { if (verbose != MagickFalse) {
3995  if ( this_kernel->next != (KernelInfo *) NULL )
3996  (void) FormatLocaleFile(stderr, " (re-iterate)");
3997  else
3998  (void) FormatLocaleFile(stderr, " (done)");
3999  }
4000  rslt_image = curr_image; /* return result, and re-iterate */
4001  }
4002  else if ( rslt_image == (Image *) NULL)
4003  { if (verbose != MagickFalse)
4004  (void) FormatLocaleFile(stderr, " (save for compose)");
4005  rslt_image = curr_image;
4006  curr_image = (Image *) image; /* continue with original image */
4007  }
4008  else
4009  { /* Add the new 'current' result to the composition
4010  **
4011  ** The removal of any 'Sync' channel flag in the Image Compositon
4012  ** below ensures the methematical compose method is applied in a
4013  ** purely mathematical way, and only to the selected channels.
4014  ** IE: Turn off SVG composition 'alpha blending'.
4015  */
4016  if (verbose != MagickFalse)
4017  (void) FormatLocaleFile(stderr, " (compose \"%s\")",
4019  (void) CompositeImage(rslt_image,curr_image,rslt_compose,MagickTrue,
4020  0,0,exception);
4021  curr_image = DestroyImage(curr_image);
4022  curr_image = (Image *) image; /* continue with original image */
4023  }
4024  if (verbose != MagickFalse)
4025  (void) FormatLocaleFile(stderr, "\n");
4026 
4027  /* loop to the next kernel in a multi-kernel list */
4028  norm_kernel = norm_kernel->next;
4029  if ( rflt_kernel != (KernelInfo *) NULL )
4030  rflt_kernel = rflt_kernel->next;
4031  kernel_number++;
4032  } /* End Loop 2: Loop over each kernel */
4033 
4034  } /* End Loop 1: compound method interation */
4035 
4036  goto exit_cleanup;
4037 
4038  /* Yes goto's are bad, but it makes cleanup lot more efficient */
4039 error_cleanup:
4040  if ( curr_image == rslt_image )
4041  curr_image = (Image *) NULL;
4042  if ( rslt_image != (Image *) NULL )
4043  rslt_image = DestroyImage(rslt_image);
4044 exit_cleanup:
4045  if ( curr_image == rslt_image || curr_image == image )
4046  curr_image = (Image *) NULL;
4047  if ( curr_image != (Image *) NULL )
4048  curr_image = DestroyImage(curr_image);
4049  if ( work_image != (Image *) NULL )
4050  work_image = DestroyImage(work_image);
4051  if ( save_image != (Image *) NULL )
4052  save_image = DestroyImage(save_image);
4053  if ( reflected_kernel != (KernelInfo *) NULL )
4054  reflected_kernel = DestroyKernelInfo(reflected_kernel);
4055  return(rslt_image);
4056 }
4057 
4058 
4059 /*
4060 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4061 % %
4062 % %
4063 % %
4064 % M o r p h o l o g y I m a g e %
4065 % %
4066 % %
4067 % %
4068 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4069 %
4070 % MorphologyImage() applies a user supplied kernel to the image according to
4071 % the given mophology method.
4072 %
4073 % This function applies any and all user defined settings before calling
4074 % the above internal function MorphologyApply().
4075 %
4076 % User defined settings include...
4077 % * Output Bias for Convolution and correlation ("-define convolve:bias=??")
4078 % * Kernel Scale/normalize settings ("-define convolve:scale=??")
4079 % This can also includes the addition of a scaled unity kernel.
4080 % * Show Kernel being applied ("-define morphology:showkernel=1")
4081 %
4082 % Other operators that do not want user supplied options interfering,
4083 % especially "convolve:bias" and "morphology:showkernel" should use
4084 % MorphologyApply() directly.
4085 %
4086 % The format of the MorphologyImage method is:
4087 %
4088 % Image *MorphologyImage(const Image *image,MorphologyMethod method,
4089 % const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
4090 %
4091 % A description of each parameter follows:
4092 %
4093 % o image: the image.
4094 %
4095 % o method: the morphology method to be applied.
4096 %
4097 % o iterations: apply the operation this many times (or no change).
4098 % A value of -1 means loop until no change found.
4099 % How this is applied may depend on the morphology method.
4100 % Typically this is a value of 1.
4101 %
4102 % o kernel: An array of double representing the morphology kernel.
4103 % Warning: kernel may be normalized for the Convolve method.
4104 %
4105 % o exception: return any errors or warnings in this structure.
4106 %
4107 */
4109  const MorphologyMethod method,const ssize_t iterations,
4110  const KernelInfo *kernel,ExceptionInfo *exception)
4111 {
4112  const char
4113  *artifact;
4114 
4116  compose;
4117 
4118  double
4119  bias;
4120 
4121  Image
4122  *morphology_image;
4123 
4124  KernelInfo
4125  *curr_kernel;
4126 
4127  curr_kernel = (KernelInfo *) kernel;
4128  bias=0.0;
4129  compose = UndefinedCompositeOp; /* use default for method */
4130 
4131  /* Apply Convolve/Correlate Normalization and Scaling Factors.
4132  * This is done BEFORE the ShowKernelInfo() function is called so that
4133  * users can see the results of the 'option:convolve:scale' option.
4134  */
4135  if ( method == ConvolveMorphology || method == CorrelateMorphology ) {
4136  /* Get the bias value as it will be needed */
4137  artifact = GetImageArtifact(image,"convolve:bias");
4138  if ( artifact != (const char *) NULL) {
4139  if (IsGeometry(artifact) == MagickFalse)
4140  (void) ThrowMagickException(exception,GetMagickModule(),
4141  OptionWarning,"InvalidSetting","'%s' '%s'",
4142  "convolve:bias",artifact);
4143  else
4144  bias=StringToDoubleInterval(artifact,(double) QuantumRange+1.0);
4145  }
4146 
4147  /* Scale kernel according to user wishes */
4148  artifact = GetImageArtifact(image,"convolve:scale");
4149  if ( artifact != (const char *) NULL ) {
4150  if (IsGeometry(artifact) == MagickFalse)
4151  (void) ThrowMagickException(exception,GetMagickModule(),
4152  OptionWarning,"InvalidSetting","'%s' '%s'",
4153  "convolve:scale",artifact);
4154  else {
4155  if ( curr_kernel == kernel )
4156  curr_kernel = CloneKernelInfo(kernel);
4157  if (curr_kernel == (KernelInfo *) NULL)
4158  return((Image *) NULL);
4159  ScaleGeometryKernelInfo(curr_kernel, artifact);
4160  }
4161  }
4162  }
4163 
4164  /* display the (normalized) kernel via stderr */
4165  artifact=GetImageArtifact(image,"morphology:showkernel");
4166  if (IsStringTrue(artifact) != MagickFalse)
4167  ShowKernelInfo(curr_kernel);
4168 
4169  /* Override the default handling of multi-kernel morphology results
4170  * If 'Undefined' use the default method
4171  * If 'None' (default for 'Convolve') re-iterate previous result
4172  * Otherwise merge resulting images using compose method given.
4173  * Default for 'HitAndMiss' is 'Lighten'.
4174  */
4175  {
4176  ssize_t
4177  parse;
4178 
4179  artifact = GetImageArtifact(image,"morphology:compose");
4180  if ( artifact != (const char *) NULL) {
4182  MagickFalse,artifact);
4183  if ( parse < 0 )
4184  (void) ThrowMagickException(exception,GetMagickModule(),
4185  OptionWarning,"UnrecognizedComposeOperator","'%s' '%s'",
4186  "morphology:compose",artifact);
4187  else
4188  compose=(CompositeOperator)parse;
4189  }
4190  }
4191  /* Apply the Morphology */
4192  morphology_image = MorphologyApply(image,method,iterations,
4193  curr_kernel,compose,bias,exception);
4194 
4195  /* Cleanup and Exit */
4196  if ( curr_kernel != kernel )
4197  curr_kernel=DestroyKernelInfo(curr_kernel);
4198  return(morphology_image);
4199 }
4200 
4201 /*
4202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4203 % %
4204 % %
4205 % %
4206 + R o t a t e K e r n e l I n f o %
4207 % %
4208 % %
4209 % %
4210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4211 %
4212 % RotateKernelInfo() rotates the kernel by the angle given.
4213 %
4214 % Currently it is restricted to 90 degree angles, of either 1D kernels
4215 % or square kernels. And 'circular' rotations of 45 degrees for 3x3 kernels.
4216 % It will ignore usless rotations for specific 'named' built-in kernels.
4217 %
4218 % The format of the RotateKernelInfo method is:
4219 %
4220 % void RotateKernelInfo(KernelInfo *kernel, double angle)
4221 %
4222 % A description of each parameter follows:
4223 %
4224 % o kernel: the Morphology/Convolution kernel
4225 %
4226 % o angle: angle to rotate in degrees
4227 %
4228 % This function is currently internal to this module only, but can be exported
4229 % to other modules if needed.
4230 */
4231 static void RotateKernelInfo(KernelInfo *kernel, double angle)
4232 {
4233  /* angle the lower kernels first */
4234  if ( kernel->next != (KernelInfo *) NULL)
4235  RotateKernelInfo(kernel->next, angle);
4236 
4237  /* WARNING: Currently assumes the kernel (rightly) is horizontally symetrical
4238  **
4239  ** TODO: expand beyond simple 90 degree rotates, flips and flops
4240  */
4241 
4242  /* Modulus the angle */
4243  angle = fmod(angle, 360.0);
4244  if ( angle < 0 )
4245  angle += 360.0;
4246 
4247  if ( 337.5 < angle || angle <= 22.5 )
4248  return; /* Near zero angle - no change! - At least not at this time */
4249 
4250  /* Handle special cases */
4251  switch (kernel->type) {
4252  /* These built-in kernels are cylindrical kernels, rotating is useless */
4253  case GaussianKernel:
4254  case DoGKernel:
4255  case LoGKernel:
4256  case DiskKernel:
4257  case PeaksKernel:
4258  case LaplacianKernel:
4259  case ChebyshevKernel:
4260  case ManhattanKernel:
4261  case EuclideanKernel:
4262  return;
4263 
4264  /* These may be rotatable at non-90 angles in the future */
4265  /* but simply rotating them in multiples of 90 degrees is useless */
4266  case SquareKernel:
4267  case DiamondKernel:
4268  case PlusKernel:
4269  case CrossKernel:
4270  return;
4271 
4272  /* These only allows a +/-90 degree rotation (by transpose) */
4273  /* A 180 degree rotation is useless */
4274  case BlurKernel:
4275  if ( 135.0 < angle && angle <= 225.0 )
4276  return;
4277  if ( 225.0 < angle && angle <= 315.0 )
4278  angle -= 180;
4279  break;
4280 
4281  default:
4282  break;
4283  }
4284  /* Attempt rotations by 45 degrees -- 3x3 kernels only */
4285  if ( 22.5 < fmod(angle,90.0) && fmod(angle,90.0) <= 67.5 )
4286  {
4287  if ( kernel->width == 3 && kernel->height == 3 )
4288  { /* Rotate a 3x3 square by 45 degree angle */
4289  double t = kernel->values[0];
4290  kernel->values[0] = kernel->values[3];
4291  kernel->values[3] = kernel->values[6];
4292  kernel->values[6] = kernel->values[7];
4293  kernel->values[7] = kernel->values[8];
4294  kernel->values[8] = kernel->values[5];
4295  kernel->values[5] = kernel->values[2];
4296  kernel->values[2] = kernel->values[1];
4297  kernel->values[1] = t;
4298  /* rotate non-centered origin */
4299  if ( kernel->x != 1 || kernel->y != 1 ) {
4300  ssize_t x,y;
4301  x = (ssize_t) kernel->x-1;
4302  y = (ssize_t) kernel->y-1;
4303  if ( x == y ) x = 0;
4304  else if ( x == 0 ) x = -y;
4305  else if ( x == -y ) y = 0;
4306  else if ( y == 0 ) y = x;
4307  kernel->x = (ssize_t) x+1;
4308  kernel->y = (ssize_t) y+1;
4309  }
4310  angle = fmod(angle+315.0, 360.0); /* angle reduced 45 degrees */
4311  kernel->angle = fmod(kernel->angle+45.0, 360.0);
4312  }
4313  else
4314  perror("Unable to rotate non-3x3 kernel by 45 degrees");
4315  }
4316  if ( 45.0 < fmod(angle, 180.0) && fmod(angle,180.0) <= 135.0 )
4317  {
4318  if ( kernel->width == 1 || kernel->height == 1 )
4319  { /* Do a transpose of a 1 dimensional kernel,
4320  ** which results in a fast 90 degree rotation of some type.
4321  */
4322  ssize_t
4323  t;
4324  t = (ssize_t) kernel->width;
4325  kernel->width = kernel->height;
4326  kernel->height = (size_t) t;
4327  t = kernel->x;
4328  kernel->x = kernel->y;
4329  kernel->y = t;
4330  if ( kernel->width == 1 ) {
4331  angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4332  kernel->angle = fmod(kernel->angle+90.0, 360.0);
4333  } else {
4334  angle = fmod(angle+90.0, 360.0); /* angle increased 90 degrees */
4335  kernel->angle = fmod(kernel->angle+270.0, 360.0);
4336  }
4337  }
4338  else if ( kernel->width == kernel->height )
4339  { /* Rotate a square array of values by 90 degrees */
4340  { register ssize_t
4341  i,j,x,y;
4342 
4343  register MagickRealType
4344  *k,t;
4345 
4346  k=kernel->values;
4347  for( i=0, x=(ssize_t) kernel->width-1; i<=x; i++, x--)
4348  for( j=0, y=(ssize_t) kernel->height-1; j<y; j++, y--)
4349  { t = k[i+j*kernel->width];
4350  k[i+j*kernel->width] = k[j+x*kernel->width];
4351  k[j+x*kernel->width] = k[x+y*kernel->width];
4352  k[x+y*kernel->width] = k[y+i*kernel->width];
4353  k[y+i*kernel->width] = t;
4354  }
4355  }
4356  /* rotate the origin - relative to center of array */
4357  { register ssize_t x,y;
4358  x = (ssize_t) (kernel->x*2-kernel->width+1);
4359  y = (ssize_t) (kernel->y*2-kernel->height+1);
4360  kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
4361  kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
4362  }
4363  angle = fmod(angle+270.0, 360.0); /* angle reduced 90 degrees */
4364  kernel->angle = fmod(kernel->angle+90.0, 360.0);
4365  }
4366  else
4367  perror("Unable to rotate a non-square, non-linear kernel 90 degrees");
4368  }
4369  if ( 135.0 < angle && angle <= 225.0 )
4370  {
4371  /* For a 180 degree rotation - also know as a reflection
4372  * This is actually a very very common operation!
4373  * Basically all that is needed is a reversal of the kernel data!
4374  * And a reflection of the origon
4375  */
4377  t;
4378 
4379  register MagickRealType
4380  *k;
4381 
4382  ssize_t
4383  i,
4384  j;
4385 
4386  k=kernel->values;
4387  j=(ssize_t) (kernel->width*kernel->height-1);
4388  for (i=0; i < j; i++, j--)
4389  t=k[i], k[i]=k[j], k[j]=t;
4390 
4391  kernel->x = (ssize_t) kernel->width - kernel->x - 1;
4392  kernel->y = (ssize_t) kernel->height - kernel->y - 1;
4393  angle = fmod(angle-180.0, 360.0); /* angle+180 degrees */
4394  kernel->angle = fmod(kernel->angle+180.0, 360.0);
4395  }
4396  /* At this point angle should at least between -45 (315) and +45 degrees
4397  * In the future some form of non-orthogonal angled rotates could be
4398  * performed here, posibily with a linear kernel restriction.
4399  */
4400 
4401  return;
4402 }
4403 
4404 /*
4405 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4406 % %
4407 % %
4408 % %
4409 % S c a l e G e o m e t r y K e r n e l I n f o %
4410 % %
4411 % %
4412 % %
4413 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4414 %
4415 % ScaleGeometryKernelInfo() takes a geometry argument string, typically
4416 % provided as a "-set option:convolve:scale {geometry}" user setting,
4417 % and modifies the kernel according to the parsed arguments of that setting.
4418 %
4419 % The first argument (and any normalization flags) are passed to
4420 % ScaleKernelInfo() to scale/normalize the kernel. The second argument
4421 % is then passed to UnityAddKernelInfo() to add a scled unity kernel
4422 % into the scaled/normalized kernel.
4423 %
4424 % The format of the ScaleGeometryKernelInfo method is:
4425 %
4426 % void ScaleGeometryKernelInfo(KernelInfo *kernel,
4427 % const double scaling_factor,const MagickStatusType normalize_flags)
4428 %
4429 % A description of each parameter follows:
4430 %
4431 % o kernel: the Morphology/Convolution kernel to modify
4432 %
4433 % o geometry:
4434 % The geometry string to parse, typically from the user provided
4435 % "-set option:convolve:scale {geometry}" setting.
4436 %
4437 */
4439  const char *geometry)
4440 {
4442  flags;
4443 
4444  GeometryInfo
4445  args;
4446 
4447  SetGeometryInfo(&args);
4448  flags = ParseGeometry(geometry, &args);
4449 
4450 #if 0
4451  /* For Debugging Geometry Input */
4452  (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
4453  flags, args.rho, args.sigma, args.xi, args.psi );
4454 #endif
4455 
4456  if ( (flags & PercentValue) != 0 ) /* Handle Percentage flag*/
4457  args.rho *= 0.01, args.sigma *= 0.01;
4458 
4459  if ( (flags & RhoValue) == 0 ) /* Set Defaults for missing args */
4460  args.rho = 1.0;
4461  if ( (flags & SigmaValue) == 0 )
4462  args.sigma = 0.0;
4463 
4464  /* Scale/Normalize the input kernel */
4465  ScaleKernelInfo(kernel, args.rho, (GeometryFlags) flags);
4466 
4467  /* Add Unity Kernel, for blending with original */
4468  if ( (flags & SigmaValue) != 0 )
4469  UnityAddKernelInfo(kernel, args.sigma);
4470 
4471  return;
4472 }
4473 /*
4474 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4475 % %
4476 % %
4477 % %
4478 % S c a l e K e r n e l I n f o %
4479 % %
4480 % %
4481 % %
4482 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4483 %
4484 % ScaleKernelInfo() scales the given kernel list by the given amount, with or
4485 % without normalization of the sum of the kernel values (as per given flags).
4486 %
4487 % By default (no flags given) the values within the kernel is scaled
4488 % directly using given scaling factor without change.
4489 %
4490 % If either of the two 'normalize_flags' are given the kernel will first be
4491 % normalized and then further scaled by the scaling factor value given.
4492 %
4493 % Kernel normalization ('normalize_flags' given) is designed to ensure that
4494 % any use of the kernel scaling factor with 'Convolve' or 'Correlate'
4495 % morphology methods will fall into -1.0 to +1.0 range. Note that for
4496 % non-HDRI versions of IM this may cause images to have any negative results
4497 % clipped, unless some 'bias' is used.
4498 %
4499 % More specifically. Kernels which only contain positive values (such as a
4500 % 'Gaussian' kernel) will be scaled so that those values sum to +1.0,
4501 % ensuring a 0.0 to +1.0 output range for non-HDRI images.
4502 %
4503 % For Kernels that contain some negative values, (such as 'Sharpen' kernels)
4504 % the kernel will be scaled by the absolute of the sum of kernel values, so
4505 % that it will generally fall within the +/- 1.0 range.
4506 %
4507 % For kernels whose values sum to zero, (such as 'Laplician' kernels) kernel
4508 % will be scaled by just the sum of the postive values, so that its output
4509 % range will again fall into the +/- 1.0 range.
4510 %
4511 % For special kernels designed for locating shapes using 'Correlate', (often
4512 % only containing +1 and -1 values, representing foreground/brackground
4513 % matching) a special normalization method is provided to scale the positive
4514 % values separately to those of the negative values, so the kernel will be
4515 % forced to become a zero-sum kernel better suited to such searches.
4516 %
4517 % WARNING: Correct normalization of the kernel assumes that the '*_range'
4518 % attributes within the kernel structure have been correctly set during the
4519 % kernels creation.
4520 %
4521 % NOTE: The values used for 'normalize_flags' have been selected specifically
4522 % to match the use of geometry options, so that '!' means NormalizeValue, '^'
4523 % means CorrelateNormalizeValue. All other GeometryFlags values are ignored.
4524 %
4525 % The format of the ScaleKernelInfo method is:
4526 %
4527 % void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
4528 % const MagickStatusType normalize_flags )
4529 %
4530 % A description of each parameter follows:
4531 %
4532 % o kernel: the Morphology/Convolution kernel
4533 %
4534 % o scaling_factor:
4535 % multiply all values (after normalization) by this factor if not
4536 % zero. If the kernel is normalized regardless of any flags.
4537 %
4538 % o normalize_flags:
4539 % GeometryFlags defining normalization method to use.
4540 % specifically: NormalizeValue, CorrelateNormalizeValue,
4541 % and/or PercentValue
4542 %
4543 */
4545  const double scaling_factor,const GeometryFlags normalize_flags)
4546 {
4547  register double
4548  pos_scale,
4549  neg_scale;
4550 
4551  register ssize_t
4552  i;
4553 
4554  /* do the other kernels in a multi-kernel list first */
4555  if ( kernel->next != (KernelInfo *) NULL)
4556  ScaleKernelInfo(kernel->next, scaling_factor, normalize_flags);
4557 
4558  /* Normalization of Kernel */
4559  pos_scale = 1.0;
4560  if ( (normalize_flags&NormalizeValue) != 0 ) {
4561  if ( fabs(kernel->positive_range + kernel->negative_range) >= MagickEpsilon )
4562  /* non-zero-summing kernel (generally positive) */
4563  pos_scale = fabs(kernel->positive_range + kernel->negative_range);
4564  else
4565  /* zero-summing kernel */
4566  pos_scale = kernel->positive_range;
4567  }
4568  /* Force kernel into a normalized zero-summing kernel */
4569  if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
4570  pos_scale = ( fabs(kernel->positive_range) >= MagickEpsilon )
4571  ? kernel->positive_range : 1.0;
4572  neg_scale = ( fabs(kernel->negative_range) >= MagickEpsilon )
4573  ? -kernel->negative_range : 1.0;
4574  }
4575  else
4576  neg_scale = pos_scale;
4577 
4578  /* finialize scaling_factor for positive and negative components */
4579  pos_scale = scaling_factor/pos_scale;
4580  neg_scale = scaling_factor/neg_scale;
4581 
4582  for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
4583  if (!IsNaN(kernel->values[i]))
4584  kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
4585 
4586  /* convolution output range */
4587  kernel->positive_range *= pos_scale;
4588  kernel->negative_range *= neg_scale;
4589  /* maximum and minimum values in kernel */
4590  kernel->maximum *= (kernel->maximum >= 0.0) ? pos_scale : neg_scale;
4591  kernel->minimum *= (kernel->minimum >= 0.0) ? pos_scale : neg_scale;
4592 
4593  /* swap kernel settings if user's scaling factor is negative */
4594  if ( scaling_factor < MagickEpsilon ) {
4595  double t;
4596  t = kernel->positive_range;
4597  kernel->positive_range = kernel->negative_range;
4598  kernel->negative_range = t;
4599  t = kernel->maximum;
4600  kernel->maximum = kernel->minimum;
4601  kernel->minimum = 1;
4602  }
4603 
4604  return;
4605 }
4606 
4607 /*
4608 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4609 % %
4610 % %
4611 % %
4612 % S h o w K e r n e l I n f o %
4613 % %
4614 % %
4615 % %
4616 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4617 %
4618 % ShowKernelInfo() outputs the details of the given kernel defination to
4619 % standard error, generally due to a users 'morphology:showkernel' option
4620 % request.
4621 %
4622 % The format of the ShowKernel method is:
4623 %
4624 % void ShowKernelInfo(const KernelInfo *kernel)
4625 %
4626 % A description of each parameter follows:
4627 %
4628 % o kernel: the Morphology/Convolution kernel
4629 %
4630 */
4632 {
4633  const KernelInfo
4634  *k;
4635 
4636  size_t
4637  c, i, u, v;
4638 
4639  for (c=0, k=kernel; k != (KernelInfo *) NULL; c++, k=k->next ) {
4640 
4641  (void) FormatLocaleFile(stderr, "Kernel");
4642  if ( kernel->next != (KernelInfo *) NULL )
4643  (void) FormatLocaleFile(stderr, " #%lu", (unsigned long) c );
4644  (void) FormatLocaleFile(stderr, " \"%s",
4646  if ( fabs(k->angle) >= MagickEpsilon )
4647  (void) FormatLocaleFile(stderr, "@%lg", k->angle);
4648  (void) FormatLocaleFile(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long)
4649  k->width,(unsigned long) k->height,(long) k->x,(long) k->y);
4650  (void) FormatLocaleFile(stderr,
4651  " with values from %.*lg to %.*lg\n",
4652  GetMagickPrecision(), k->minimum,
4653  GetMagickPrecision(), k->maximum);
4654  (void) FormatLocaleFile(stderr, "Forming a output range from %.*lg to %.*lg",
4655  GetMagickPrecision(), k->negative_range,
4656  GetMagickPrecision(), k->positive_range);
4657  if ( fabs(k->positive_range+k->negative_range) < MagickEpsilon )
4658  (void) FormatLocaleFile(stderr, " (Zero-Summing)\n");
4659  else if ( fabs(k->positive_range+k->negative_range-1.0) < MagickEpsilon )
4660  (void) FormatLocaleFile(stderr, " (Normalized)\n");
4661  else
4662  (void) FormatLocaleFile(stderr, " (Sum %.*lg)\n",
4663  GetMagickPrecision(), k->positive_range+k->negative_range);
4664  for (i=v=0; v < k->height; v++) {
4665  (void) FormatLocaleFile(stderr, "%2lu:", (unsigned long) v );
4666  for (u=0; u < k->width; u++, i++)
4667  if (IsNaN(k->values[i]))
4668  (void) FormatLocaleFile(stderr," %*s", GetMagickPrecision()+3, "nan");
4669  else
4670  (void) FormatLocaleFile(stderr," %*.*lg", GetMagickPrecision()+3,
4671  GetMagickPrecision(), (double) k->values[i]);
4672  (void) FormatLocaleFile(stderr,"\n");
4673  }
4674  }
4675 }
4676 
4677 /*
4678 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4679 % %
4680 % %
4681 % %
4682 % U n i t y A d d K e r n a l I n f o %
4683 % %
4684 % %
4685 % %
4686 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4687 %
4688 % UnityAddKernelInfo() Adds a given amount of the 'Unity' Convolution Kernel
4689 % to the given pre-scaled and normalized Kernel. This in effect adds that
4690 % amount of the original image into the resulting convolution kernel. This
4691 % value is usually provided by the user as a percentage value in the
4692 % 'convolve:scale' setting.
4693 %
4694 % The resulting effect is to convert the defined kernels into blended
4695 % soft-blurs, unsharp kernels or into sharpening kernels.
4696 %
4697 % The format of the UnityAdditionKernelInfo method is:
4698 %
4699 % void UnityAdditionKernelInfo(KernelInfo *kernel, const double scale )
4700 %
4701 % A description of each parameter follows:
4702 %
4703 % o kernel: the Morphology/Convolution kernel
4704 %
4705 % o scale:
4706 % scaling factor for the unity kernel to be added to
4707 % the given kernel.
4708 %
4709 */
4711  const double scale)
4712 {
4713  /* do the other kernels in a multi-kernel list first */
4714  if ( kernel->next != (KernelInfo *) NULL)
4715  UnityAddKernelInfo(kernel->next, scale);
4716 
4717  /* Add the scaled unity kernel to the existing kernel */
4718  kernel->values[kernel->x+kernel->y*kernel->width] += scale;
4719  CalcKernelMetaData(kernel); /* recalculate the meta-data */
4720 
4721  return;
4722 }
4723 
4724 /*
4725 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4726 % %
4727 % %
4728 % %
4729 % Z e r o K e r n e l N a n s %
4730 % %
4731 % %
4732 % %
4733 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4734 %
4735 % ZeroKernelNans() replaces any special 'nan' value that may be present in
4736 % the kernel with a zero value. This is typically done when the kernel will
4737 % be used in special hardware (GPU) convolution processors, to simply
4738 % matters.
4739 %
4740 % The format of the ZeroKernelNans method is:
4741 %
4742 % void ZeroKernelNans (KernelInfo *kernel)
4743 %
4744 % A description of each parameter follows:
4745 %
4746 % o kernel: the Morphology/Convolution kernel
4747 %
4748 */
4750 {
4751  register size_t
4752  i;
4753 
4754  /* do the other kernels in a multi-kernel list first */
4755  if (kernel->next != (KernelInfo *) NULL)
4756  ZeroKernelNans(kernel->next);
4757 
4758  for (i=0; i < (kernel->width*kernel->height); i++)
4759  if (IsNaN(kernel->values[i]))
4760  kernel->values[i]=0.0;
4761 
4762  return;
4763 }
double psi
Definition: geometry.h:105
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:3613
#define magick_restrict
Definition: MagickCore.h:41
MagickDoubleType MagickRealType
Definition: magick-type.h:118
MagickExport CacheView * DestroyCacheView(CacheView *cache_view)
Definition: cache-view.c:252
static MagickBooleanType SetImageProgress(const Image *image, const char *tag, const MagickOffsetType offset, const MagickSizeType extent)
static void ExpandMirrorKernelInfo(KernelInfo *)
Definition: morphology.c:2328
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:345
#define MagickAssumeAligned(address)
ssize_t y
Definition: geometry.h:115
double positive_range
Definition: morphology.h:119
MagickExport ssize_t ParseCommandOption(const CommandOption option, const MagickBooleanType list, const char *options)
Definition: option.c:2953
size_t height
Definition: morphology.h:108
MagickExport KernelInfo * DestroyKernelInfo(KernelInfo *kernel)
Definition: morphology.c:2268
#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:4108
double rho
Definition: geometry.h:105
MagickExport void SetGeometryInfo(GeometryInfo *geometry_info)
Definition: geometry.c:1605
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:2561
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:3226
#define Minimize(assign, value)
Definition: morphology.c:92
static double StringToDouble(const char *magick_restrict string, char **magick_restrict sentinal)
ssize_t x
Definition: morphology.h:112
static PixelTrait GetPixelChannelTraits(const Image *magick_restrict image, const PixelChannel channel)
#define MagickPI
Definition: image-private.h:30
MagickExport ssize_t FormatLocaleString(char *magick_restrict string, const size_t length, const char *magick_restrict format,...)
Definition: locale.c:473
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:4438
size_t signature
Definition: morphology.h:129
MagickExport void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor, const GeometryFlags normalize_flags)
Definition: morphology.c:4544
#define MagickSQ2
Definition: image-private.h:33
#define MagickEpsilon
Definition: magick-type.h:110
double sigma
Definition: geometry.h:105
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:539
#define Maximize(assign, value)
Definition: morphology.c:93
MagickExport char * FileToString(const char *filename, const size_t extent, ExceptionInfo *exception)
Definition: string.c:985
ssize_t MagickOffsetType
Definition: magick-type.h:127
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:4631
MagickExport MagickBooleanType IsGeometry(const char *geometry)
Definition: geometry.c:593
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:378
MagickExport MagickBooleanType SetImageAlphaChannel(Image *image, const AlphaChannelOption alpha_type, ExceptionInfo *exception)
Definition: channel.c:966
MagickBooleanType
Definition: magick-type.h:156
unsigned int MagickStatusType
Definition: magick-type.h:119
static double PerceptibleReciprocal(const double x)
MagickExport void * ResetMagickMemory(void *memory, int byte, const size_t size)
Definition: memory.c:1164
MagickExport const char * CommandOptionToMnemonic(const CommandOption option, const ssize_t type)
Definition: option.c:2666
#define Magick2PI
Definition: image-private.h:28
static Quantum GetPixelWriteMask(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
MagickExport void * AcquireQuantumMemory(const size_t count, const size_t quantum)
Definition: memory.c:529
MagickPrivate size_t GetOptimalKernelWidth2D(const double, const double)
Definition: gem.c:1635
static int GetOpenMPThreadId(void)
#define MagickSQ2PI
Definition: image-private.h:34
#define RestoreMSCWarning
Definition: studio.h:346
static size_t fact(size_t n)
Definition: morphology.c:97
#define MagickPathExtent
MagickExport void * RelinquishAlignedMemory(void *memory)
Definition: memory.c:1001
MagickExport MagickBooleanType IsStringTrue(const char *value)
Definition: string.c:1464
static void CalcKernelMetaData(KernelInfo *)
Definition: morphology.c:2480
MagickExport int GetMagickPrecision(void)
Definition: magick.c:865
#define MagickMaximumValue
Definition: magick-type.h:111
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:1058
size_t width
Definition: morphology.h:108
size_t signature
Definition: image.h:354
#define QuantumScale
Definition: magick-type.h:113
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:2508
PixelChannel
Definition: pixel.h:66
MagickExport void * AcquireAlignedMemory(const size_t count, const size_t quantum)
Definition: memory.c:238
#define MagickMax(x, y)
Definition: image-private.h:26
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:1409
#define IsNaN(a)
Definition: magick-type.h:179
MagickPrivate void ZeroKernelNans(KernelInfo *kernel)
Definition: morphology.c:4749
#define GetMagickModule()
Definition: log.h:28
static Quantum ClampToQuantum(const MagickRealType value)
Definition: quantum.h:84
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:4231
MagickExport void GetNextToken(const char *start, const char **end, const size_t extent, char *token)
Definition: token.c:171
static double StringToDoubleInterval(const char *string, const double interval)
static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1, const KernelInfo *kernel2)
Definition: morphology.c:2392
unsigned short Quantum
Definition: magick-type.h:82
double xi
Definition: geometry.h:105
MagickExport KernelInfo * CloneKernelInfo(const KernelInfo *kernel)
Definition: morphology.c:2213
MagickExport char * DestroyString(char *string)
Definition: string.c:810
MagickExport void * AcquireMagickMemory(const size_t size)
Definition: memory.c:458
MagickExport MagickStatusType ParseGeometry(const char *geometry, GeometryInfo *geometry_info)
Definition: geometry.c:836
static void SetPixelChannel(const Image *magick_restrict image, const PixelChannel channel, const Quantum quantum, Quantum *magick_restrict pixel)
ssize_t x
Definition: geometry.h:115
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1038
CompositeOperator compose
Definition: image.h:234
MagickExport void UnityAddKernelInfo(KernelInfo *kernel, const double scale)
Definition: morphology.c:4710
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:132
MagickExport MagickRealType GetPixelIntensity(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
Definition: pixel.c:2356
static void ExpandRotateKernelInfo(KernelInfo *, const double)
Definition: morphology.c:2420
#define MorphologyTag
#define KernelRank
MagickExport Image * DestroyImage(Image *image)
Definition: image.c:1182
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:799
#define QuantumRange
Definition: magick-type.h:83
MagickRealType * values
Definition: morphology.h:116