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