MagickCore  7.1.0
segment.c
Go to the documentation of this file.
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % SSSSS EEEEE GGGG M M EEEEE N N TTTTT %
7 % SS E G MM MM E NN N T %
8 % SSS EEE G GGG M M M EEE N N N T %
9 % SS E G G M M E N NN T %
10 % SSSSS EEEEE GGGG M M EEEEE N N T %
11 % %
12 % %
13 % MagickCore Methods to Segment an Image with Thresholding Fuzzy c-Means %
14 % %
15 % Software Design %
16 % Cristy %
17 % April 1993 %
18 % %
19 % %
20 % Copyright 1999-2021 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 % Segment segments an image by analyzing the histograms of the color
37 % components and identifying units that are homogeneous with the fuzzy
38 % c-means technique. The scale-space filter analyzes the histograms of
39 % the three color components of the image and identifies a set of
40 % classes. The extents of each class is used to coarsely segment the
41 % image with thresholding. The color associated with each class is
42 % determined by the mean color of all pixels within the extents of a
43 % particular class. Finally, any unclassified pixels are assigned to
44 % the closest class with the fuzzy c-means technique.
45 %
46 % The fuzzy c-Means algorithm can be summarized as follows:
47 %
48 % o Build a histogram, one for each color component of the image.
49 %
50 % o For each histogram, successively apply the scale-space filter and
51 % build an interval tree of zero crossings in the second derivative
52 % at each scale. Analyze this scale-space ''fingerprint'' to
53 % determine which peaks and valleys in the histogram are most
54 % predominant.
55 %
56 % o The fingerprint defines intervals on the axis of the histogram.
57 % Each interval contains either a minima or a maxima in the original
58 % signal. If each color component lies within the maxima interval,
59 % that pixel is considered ''classified'' and is assigned an unique
60 % class number.
61 %
62 % o Any pixel that fails to be classified in the above thresholding
63 % pass is classified using the fuzzy c-Means technique. It is
64 % assigned to one of the classes discovered in the histogram analysis
65 % phase.
66 %
67 % The fuzzy c-Means technique attempts to cluster a pixel by finding
68 % the local minima of the generalized within group sum of squared error
69 % objective function. A pixel is assigned to the closest class of
70 % which the fuzzy membership has a maximum value.
71 %
72 % Segment is strongly based on software written by Andy Gallo,
73 % University of Delaware.
74 %
75 % The following reference was used in creating this program:
76 %
77 % Young Won Lim, Sang Uk Lee, "On The Color Image Segmentation
78 % Algorithm Based on the Thresholding and the Fuzzy c-Means
79 % Techniques", Pattern Recognition, Volume 23, Number 9, pages
80 % 935-952, 1990.
81 %
82 %
83 */
84 
85 #include "MagickCore/studio.h"
86 #include "MagickCore/cache.h"
87 #include "MagickCore/color.h"
88 #include "MagickCore/colormap.h"
89 #include "MagickCore/colorspace.h"
91 #include "MagickCore/exception.h"
93 #include "MagickCore/image.h"
95 #include "MagickCore/memory_.h"
97 #include "MagickCore/monitor.h"
100 #include "MagickCore/quantize.h"
101 #include "MagickCore/quantum.h"
103 #include "MagickCore/resource_.h"
104 #include "MagickCore/segment.h"
105 #include "MagickCore/string_.h"
107 
108 /*
109  Define declarations.
110 */
111 #define MaxDimension 3
112 #define DeltaTau 0.5f
113 #if defined(FastClassify)
114 #define WeightingExponent 2.0
115 #define SegmentPower(ratio) (ratio)
116 #else
117 #define WeightingExponent 2.5
118 #define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0)));
119 #endif
120 #define Tau 5.2f
121 
122 /*
123  Typedef declarations.
124 */
125 typedef struct _ExtentPacket
126 {
127  double
129 
130  ssize_t
132  left,
133  right;
134 } ExtentPacket;
135 
136 typedef struct _Cluster
137 {
138  struct _Cluster
139  *next;
140 
143  green,
144  blue;
145 
146  ssize_t
148  id;
149 } Cluster;
150 
151 typedef struct _IntervalTree
152 {
153  double
155 
156  ssize_t
158  right;
159 
160  double
162  stability;
163 
165  *sibling,
166  *child;
167 } IntervalTree;
168 
169 typedef struct _ZeroCrossing
170 {
171  double
173  histogram[256];
174 
175  short
176  crossings[256];
177 } ZeroCrossing;
178 
179 /*
180  Constant declarations.
181 */
182 static const int
183  Blue = 2,
184  Green = 1,
185  Red = 0,
187  TreeLength = 600;
188 
189 /*
190  Method prototypes.
191 */
192 static double
193  OptimalTau(const ssize_t *,const double,const double,const double,
194  const double,short *);
195 
196 static ssize_t
197  DefineRegion(const short *,ExtentPacket *);
198 
199 static void
201  InitializeHistogram(const Image *,ssize_t **,ExceptionInfo *),
202  ScaleSpace(const ssize_t *,const double,double *),
203  ZeroCrossHistogram(double *,const double,short *);
204 
205 /*
206 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
207 % %
208 % %
209 % %
210 + C l a s s i f y %
211 % %
212 % %
213 % %
214 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
215 %
216 % Classify() defines one or more classes. Each pixel is thresholded to
217 % determine which class it belongs to. If the class is not identified it is
218 % assigned to the closest class based on the fuzzy c-Means technique.
219 %
220 % The format of the Classify method is:
221 %
222 % MagickBooleanType Classify(Image *image,short **extrema,
223 % const double cluster_threshold,const double weighting_exponent,
224 % const MagickBooleanType verbose,ExceptionInfo *exception)
225 %
226 % A description of each parameter follows.
227 %
228 % o image: the image.
229 %
230 % o extrema: Specifies a pointer to an array of integers. They
231 % represent the peaks and valleys of the histogram for each color
232 % component.
233 %
234 % o cluster_threshold: This double represents the minimum number of
235 % pixels contained in a hexahedra before it can be considered valid
236 % (expressed as a percentage).
237 %
238 % o weighting_exponent: Specifies the membership weighting exponent.
239 %
240 % o verbose: A value greater than zero prints detailed information about
241 % the identified classes.
242 %
243 % o exception: return any errors or warnings in this structure.
244 %
245 */
246 static MagickBooleanType Classify(Image *image,short **extrema,
247  const double cluster_threshold,const double weighting_exponent,
248  const MagickBooleanType verbose,ExceptionInfo *exception)
249 {
250 #define SegmentImageTag "Segment/Image"
251 #define ThrowClassifyException(severity,tag,label) \
252 {\
253  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster) \
254  { \
255  next_cluster=cluster->next; \
256  cluster=(Cluster *) RelinquishMagickMemory(cluster); \
257  } \
258  if (squares != (double *) NULL) \
259  { \
260  squares-=255; \
261  free_squares=squares; \
262  free_squares=(double *) RelinquishMagickMemory(free_squares); \
263  } \
264  ThrowBinaryException(severity,tag,label); \
265 }
266 
267  CacheView
268  *image_view;
269 
270  Cluster
271  *cluster,
272  *head,
273  *last_cluster,
274  *next_cluster;
275 
276  double
277  *free_squares;
278 
280  blue,
281  green,
282  red;
283 
285  progress;
286 
288  status;
289 
290  ssize_t
291  i;
292 
293  double
294  *squares;
295 
296  size_t
297  number_clusters;
298 
299  ssize_t
300  count,
301  y;
302 
303  /*
304  Form clusters.
305  */
306  cluster=(Cluster *) NULL;
307  head=(Cluster *) NULL;
308  squares=(double *) NULL;
309  (void) memset(&red,0,sizeof(red));
310  (void) memset(&green,0,sizeof(green));
311  (void) memset(&blue,0,sizeof(blue));
312  while (DefineRegion(extrema[Red],&red) != 0)
313  {
314  green.index=0;
315  while (DefineRegion(extrema[Green],&green) != 0)
316  {
317  blue.index=0;
318  while (DefineRegion(extrema[Blue],&blue) != 0)
319  {
320  /*
321  Allocate a new class.
322  */
323  if (head != (Cluster *) NULL)
324  {
325  cluster->next=(Cluster *) AcquireQuantumMemory(1,
326  sizeof(*cluster->next));
327  cluster=cluster->next;
328  }
329  else
330  {
331  cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));
332  head=cluster;
333  }
334  if (cluster == (Cluster *) NULL)
335  ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",
336  image->filename);
337  /*
338  Initialize a new class.
339  */
340  (void) memset(cluster,0,sizeof(*cluster));
341  cluster->red=red;
342  cluster->green=green;
343  cluster->blue=blue;
344  }
345  }
346  }
347  if (head == (Cluster *) NULL)
348  {
349  /*
350  No classes were identified-- create one.
351  */
352  cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));
353  if (cluster == (Cluster *) NULL)
354  ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",
355  image->filename);
356  /*
357  Initialize a new class.
358  */
359  (void) memset(cluster,0,sizeof(*cluster));
360  cluster->red=red;
361  cluster->green=green;
362  cluster->blue=blue;
363  head=cluster;
364  }
365  /*
366  Count the pixels for each cluster.
367  */
368  status=MagickTrue;
369  count=0;
370  progress=0;
371  image_view=AcquireVirtualCacheView(image,exception);
372  for (y=0; y < (ssize_t) image->rows; y++)
373  {
374  const Quantum
375  *p;
376 
377  ssize_t
378  x;
379 
380  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
381  if (p == (const Quantum *) NULL)
382  break;
383  for (x=0; x < (ssize_t) image->columns; x++)
384  {
385  PixelInfo
386  pixel;
387 
388  pixel.red=(double) ScaleQuantumToChar(GetPixelRed(image,p));
389  pixel.green=(double) ScaleQuantumToChar(GetPixelGreen(image,p));
390  pixel.blue=(double) ScaleQuantumToChar(GetPixelBlue(image,p));
391  for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
392  if ((pixel.red >= (double) (cluster->red.left-SafeMargin)) &&
393  (pixel.red <= (double) (cluster->red.right+SafeMargin)) &&
394  (pixel.green >= (double) (cluster->green.left-SafeMargin)) &&
395  (pixel.green <= (double) (cluster->green.right+SafeMargin)) &&
396  (pixel.blue >= (double) (cluster->blue.left-SafeMargin)) &&
397  (pixel.blue <= (double) (cluster->blue.right+SafeMargin)))
398  {
399  /*
400  Count this pixel.
401  */
402  count++;
403  cluster->red.center+=pixel.red;
404  cluster->green.center+=pixel.green;
405  cluster->blue.center+=pixel.blue;
406  cluster->count++;
407  break;
408  }
409  p+=GetPixelChannels(image);
410  }
411  if (image->progress_monitor != (MagickProgressMonitor) NULL)
412  {
414  proceed;
415 
416 #if defined(MAGICKCORE_OPENMP_SUPPORT)
417  #pragma omp atomic
418 #endif
419  progress++;
420  proceed=SetImageProgress(image,SegmentImageTag,progress,2*image->rows);
421  if (proceed == MagickFalse)
422  status=MagickFalse;
423  }
424  }
425  image_view=DestroyCacheView(image_view);
426  /*
427  Remove clusters that do not meet minimum cluster threshold.
428  */
429  count=0;
430  last_cluster=head;
431  next_cluster=head;
432  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
433  {
434  next_cluster=cluster->next;
435  if ((cluster->count > 0) &&
436  (cluster->count >= (count*cluster_threshold/100.0)))
437  {
438  /*
439  Initialize cluster.
440  */
441  cluster->id=count;
442  cluster->red.center/=cluster->count;
443  cluster->green.center/=cluster->count;
444  cluster->blue.center/=cluster->count;
445  count++;
446  last_cluster=cluster;
447  continue;
448  }
449  /*
450  Delete cluster.
451  */
452  if (cluster == head)
453  head=next_cluster;
454  else
455  last_cluster->next=next_cluster;
456  cluster=(Cluster *) RelinquishMagickMemory(cluster);
457  }
458  number_clusters=(size_t) count;
459  if (verbose != MagickFalse)
460  {
461  /*
462  Print cluster statistics.
463  */
464  (void) FormatLocaleFile(stdout,"Fuzzy C-means Statistics\n");
465  (void) FormatLocaleFile(stdout,"===================\n\n");
466  (void) FormatLocaleFile(stdout,"\tCluster Threshold = %g\n",(double)
467  cluster_threshold);
468  (void) FormatLocaleFile(stdout,"\tWeighting Exponent = %g\n",(double)
469  weighting_exponent);
470  (void) FormatLocaleFile(stdout,"\tTotal Number of Clusters = %.20g\n\n",
471  (double) number_clusters);
472  /*
473  Print the total number of points per cluster.
474  */
475  (void) FormatLocaleFile(stdout,"\n\nNumber of Vectors Per Cluster\n");
476  (void) FormatLocaleFile(stdout,"=============================\n\n");
477  for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
478  (void) FormatLocaleFile(stdout,"Cluster #%.20g = %.20g\n",(double)
479  cluster->id,(double) cluster->count);
480  /*
481  Print the cluster extents.
482  */
483  (void) FormatLocaleFile(stdout,
484  "\n\n\nCluster Extents: (Vector Size: %d)\n",MaxDimension);
485  (void) FormatLocaleFile(stdout,"================");
486  for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
487  {
488  (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
489  cluster->id);
490  (void) FormatLocaleFile(stdout,
491  "%.20g-%.20g %.20g-%.20g %.20g-%.20g\n",(double)
492  cluster->red.left,(double) cluster->red.right,(double)
493  cluster->green.left,(double) cluster->green.right,(double)
494  cluster->blue.left,(double) cluster->blue.right);
495  }
496  /*
497  Print the cluster center values.
498  */
499  (void) FormatLocaleFile(stdout,
500  "\n\n\nCluster Center Values: (Vector Size: %d)\n",MaxDimension);
501  (void) FormatLocaleFile(stdout,"=====================");
502  for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
503  {
504  (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
505  cluster->id);
506  (void) FormatLocaleFile(stdout,"%g %g %g\n",(double)
507  cluster->red.center,(double) cluster->green.center,(double)
508  cluster->blue.center);
509  }
510  (void) FormatLocaleFile(stdout,"\n");
511  }
512  if (number_clusters > 256)
513  ThrowClassifyException(ImageError,"TooManyClusters",image->filename);
514  /*
515  Speed up distance calculations.
516  */
517  squares=(double *) AcquireQuantumMemory(513UL,sizeof(*squares));
518  if (squares == (double *) NULL)
519  ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",
520  image->filename);
521  squares+=255;
522  for (i=(-255); i <= 255; i++)
523  squares[i]=(double) i*(double) i;
524  /*
525  Allocate image colormap.
526  */
527  if (AcquireImageColormap(image,number_clusters,exception) == MagickFalse)
528  ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",
529  image->filename);
530  i=0;
531  for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
532  {
533  image->colormap[i].red=(double) ScaleCharToQuantum((unsigned char)
534  (cluster->red.center+0.5));
535  image->colormap[i].green=(double) ScaleCharToQuantum((unsigned char)
536  (cluster->green.center+0.5));
537  image->colormap[i].blue=(double) ScaleCharToQuantum((unsigned char)
538  (cluster->blue.center+0.5));
539  i++;
540  }
541  /*
542  Do course grain classes.
543  */
544  image_view=AcquireAuthenticCacheView(image,exception);
545 #if defined(MAGICKCORE_OPENMP_SUPPORT)
546  #pragma omp parallel for schedule(static) shared(progress,status) \
547  magick_number_threads(image,image,image->rows,1)
548 #endif
549  for (y=0; y < (ssize_t) image->rows; y++)
550  {
551  Cluster
552  *c;
553 
554  const PixelInfo
555  *magick_restrict p;
556 
557  ssize_t
558  x;
559 
560  Quantum
561  *magick_restrict q;
562 
563  if (status == MagickFalse)
564  continue;
565  q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
566  if (q == (Quantum *) NULL)
567  {
568  status=MagickFalse;
569  continue;
570  }
571  for (x=0; x < (ssize_t) image->columns; x++)
572  {
573  PixelInfo
574  pixel;
575 
576  SetPixelIndex(image,(Quantum) 0,q);
577  pixel.red=(double) ScaleQuantumToChar(GetPixelRed(image,q));
578  pixel.green=(double) ScaleQuantumToChar(GetPixelGreen(image,q));
579  pixel.blue=(double) ScaleQuantumToChar(GetPixelBlue(image,q));
580  for (c=head; c != (Cluster *) NULL; c=c->next)
581  {
582  if ((pixel.red >= (double) (c->red.left-SafeMargin)) &&
583  (pixel.red <= (double) (c->red.right+SafeMargin)) &&
584  (pixel.green >= (double) (c->green.left-SafeMargin)) &&
585  (pixel.green <= (double) (c->green.right+SafeMargin)) &&
586  (pixel.blue >= (double) (c->blue.left-SafeMargin)) &&
587  (pixel.blue <= (double) (c->blue.right+SafeMargin)))
588  {
589  /*
590  Classify this pixel.
591  */
592  SetPixelIndex(image,(Quantum) c->id,q);
593  break;
594  }
595  }
596  if (c == (Cluster *) NULL)
597  {
598  double
599  distance_squared,
600  local_minima,
601  numerator,
602  ratio,
603  sum;
604 
605  ssize_t
606  j,
607  k;
608 
609  /*
610  Compute fuzzy membership.
611  */
612  local_minima=0.0;
613  for (j=0; j < (ssize_t) image->colors; j++)
614  {
615  sum=0.0;
616  p=image->colormap+j;
617  distance_squared=
618  squares[(ssize_t) (pixel.red-ScaleQuantumToChar(p->red))]+
619  squares[(ssize_t) (pixel.green-ScaleQuantumToChar(p->green))]+
620  squares[(ssize_t) (pixel.blue-ScaleQuantumToChar(p->blue))];
621  numerator=distance_squared;
622  for (k=0; k < (ssize_t) image->colors; k++)
623  {
624  p=image->colormap+k;
625  distance_squared=
626  squares[(ssize_t) (pixel.red-ScaleQuantumToChar(p->red))]+
627  squares[(ssize_t) (pixel.green-ScaleQuantumToChar(p->green))]+
628  squares[(ssize_t) (pixel.blue-ScaleQuantumToChar(p->blue))];
629  ratio=numerator/distance_squared;
630  sum+=SegmentPower(ratio);
631  }
632  if ((sum != 0.0) && ((1.0/sum) > local_minima))
633  {
634  /*
635  Classify this pixel.
636  */
637  local_minima=1.0/sum;
638  SetPixelIndex(image,(Quantum) j,q);
639  }
640  }
641  }
642  q+=GetPixelChannels(image);
643  }
644  if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
645  status=MagickFalse;
646  if (image->progress_monitor != (MagickProgressMonitor) NULL)
647  {
649  proceed;
650 
651 #if defined(MAGICKCORE_OPENMP_SUPPORT)
652  #pragma omp atomic
653 #endif
654  progress++;
655  proceed=SetImageProgress(image,SegmentImageTag,progress,2*image->rows);
656  if (proceed == MagickFalse)
657  status=MagickFalse;
658  }
659  }
660  image_view=DestroyCacheView(image_view);
661  status&=SyncImage(image,exception);
662  /*
663  Relinquish resources.
664  */
665  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
666  {
667  next_cluster=cluster->next;
668  cluster=(Cluster *) RelinquishMagickMemory(cluster);
669  }
670  squares-=255;
671  free_squares=squares;
672  free_squares=(double *) RelinquishMagickMemory(free_squares);
673  return(MagickTrue);
674 }
675 
676 /*
677 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
678 % %
679 % %
680 % %
681 + C o n s o l i d a t e C r o s s i n g s %
682 % %
683 % %
684 % %
685 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
686 %
687 % ConsolidateCrossings() guarantees that an even number of zero crossings
688 % always lie between two crossings.
689 %
690 % The format of the ConsolidateCrossings method is:
691 %
692 % ConsolidateCrossings(ZeroCrossing *zero_crossing,
693 % const size_t number_crossings)
694 %
695 % A description of each parameter follows.
696 %
697 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
698 %
699 % o number_crossings: This size_t specifies the number of elements
700 % in the zero_crossing array.
701 %
702 */
703 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
704  const size_t number_crossings)
705 {
706  ssize_t
707  i,
708  j,
709  k,
710  l;
711 
712  ssize_t
713  center,
714  correct,
715  count,
716  left,
717  right;
718 
719  /*
720  Consolidate zero crossings.
721  */
722  for (i=(ssize_t) number_crossings-1; i >= 0; i--)
723  for (j=0; j <= 255; j++)
724  {
725  if (zero_crossing[i].crossings[j] == 0)
726  continue;
727  /*
728  Find the entry that is closest to j and still preserves the
729  property that there are an even number of crossings between
730  intervals.
731  */
732  for (k=j-1; k > 0; k--)
733  if (zero_crossing[i+1].crossings[k] != 0)
734  break;
735  left=MagickMax(k,0);
736  center=j;
737  for (k=j+1; k < 255; k++)
738  if (zero_crossing[i+1].crossings[k] != 0)
739  break;
740  right=MagickMin(k,255);
741  /*
742  K is the zero crossing just left of j.
743  */
744  for (k=j-1; k > 0; k--)
745  if (zero_crossing[i].crossings[k] != 0)
746  break;
747  if (k < 0)
748  k=0;
749  /*
750  Check center for an even number of crossings between k and j.
751  */
752  correct=(-1);
753  if (zero_crossing[i+1].crossings[j] != 0)
754  {
755  count=0;
756  for (l=k+1; l < center; l++)
757  if (zero_crossing[i+1].crossings[l] != 0)
758  count++;
759  if (((count % 2) == 0) && (center != k))
760  correct=center;
761  }
762  /*
763  Check left for an even number of crossings between k and j.
764  */
765  if (correct == -1)
766  {
767  count=0;
768  for (l=k+1; l < left; l++)
769  if (zero_crossing[i+1].crossings[l] != 0)
770  count++;
771  if (((count % 2) == 0) && (left != k))
772  correct=left;
773  }
774  /*
775  Check right for an even number of crossings between k and j.
776  */
777  if (correct == -1)
778  {
779  count=0;
780  for (l=k+1; l < right; l++)
781  if (zero_crossing[i+1].crossings[l] != 0)
782  count++;
783  if (((count % 2) == 0) && (right != k))
784  correct=right;
785  }
786  l=(ssize_t) zero_crossing[i].crossings[j];
787  zero_crossing[i].crossings[j]=0;
788  if (correct != -1)
789  zero_crossing[i].crossings[correct]=(short) l;
790  }
791 }
792 
793 /*
794 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
795 % %
796 % %
797 % %
798 + D e f i n e R e g i o n %
799 % %
800 % %
801 % %
802 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
803 %
804 % DefineRegion() defines the left and right boundaries of a peak region.
805 %
806 % The format of the DefineRegion method is:
807 %
808 % ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
809 %
810 % A description of each parameter follows.
811 %
812 % o extrema: Specifies a pointer to an array of integers. They
813 % represent the peaks and valleys of the histogram for each color
814 % component.
815 %
816 % o extents: This pointer to an ExtentPacket represent the extends
817 % of a particular peak or valley of a color component.
818 %
819 */
820 static ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
821 {
822  /*
823  Initialize to default values.
824  */
825  extents->left=0;
826  extents->center=0.0;
827  extents->right=255;
828  /*
829  Find the left side (maxima).
830  */
831  for ( ; extents->index <= 255; extents->index++)
832  if (extrema[extents->index] > 0)
833  break;
834  if (extents->index > 255)
835  return(MagickFalse); /* no left side - no region exists */
836  extents->left=extents->index;
837  /*
838  Find the right side (minima).
839  */
840  for ( ; extents->index <= 255; extents->index++)
841  if (extrema[extents->index] < 0)
842  break;
843  extents->right=extents->index-1;
844  return(MagickTrue);
845 }
846 
847 /*
848 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
849 % %
850 % %
851 % %
852 + D e r i v a t i v e H i s t o g r a m %
853 % %
854 % %
855 % %
856 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
857 %
858 % DerivativeHistogram() determines the derivative of the histogram using
859 % central differencing.
860 %
861 % The format of the DerivativeHistogram method is:
862 %
863 % DerivativeHistogram(const double *histogram,
864 % double *derivative)
865 %
866 % A description of each parameter follows.
867 %
868 % o histogram: Specifies an array of doubles representing the number
869 % of pixels for each intensity of a particular color component.
870 %
871 % o derivative: This array of doubles is initialized by
872 % DerivativeHistogram to the derivative of the histogram using central
873 % differencing.
874 %
875 */
876 static void DerivativeHistogram(const double *histogram,
877  double *derivative)
878 {
879  ssize_t
880  i,
881  n;
882 
883  /*
884  Compute endpoints using second order polynomial interpolation.
885  */
886  n=255;
887  derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);
888  derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);
889  /*
890  Compute derivative using central differencing.
891  */
892  for (i=1; i < n; i++)
893  derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
894  return;
895 }
896 
897 /*
898 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
899 % %
900 % %
901 % %
902 + G e t I m a g e D y n a m i c T h r e s h o l d %
903 % %
904 % %
905 % %
906 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
907 %
908 % GetImageDynamicThreshold() returns the dynamic threshold for an image.
909 %
910 % The format of the GetImageDynamicThreshold method is:
911 %
912 % MagickBooleanType GetImageDynamicThreshold(const Image *image,
913 % const double cluster_threshold,const double smooth_threshold,
914 % PixelInfo *pixel,ExceptionInfo *exception)
915 %
916 % A description of each parameter follows.
917 %
918 % o image: the image.
919 %
920 % o cluster_threshold: This double represents the minimum number of
921 % pixels contained in a hexahedra before it can be considered valid
922 % (expressed as a percentage).
923 %
924 % o smooth_threshold: the smoothing threshold eliminates noise in the second
925 % derivative of the histogram. As the value is increased, you can expect a
926 % smoother second derivative.
927 %
928 % o pixel: return the dynamic threshold here.
929 %
930 % o exception: return any errors or warnings in this structure.
931 %
932 */
934  const double cluster_threshold,const double smooth_threshold,
935  PixelInfo *pixel,ExceptionInfo *exception)
936 {
937  Cluster
938  *background,
939  *cluster,
940  *object,
941  *head,
942  *last_cluster,
943  *next_cluster;
944 
946  blue,
947  green,
948  red;
949 
951  proceed;
952 
953  double
954  threshold;
955 
956  const Quantum
957  *p;
958 
959  ssize_t
960  i,
961  x;
962 
963  short
964  *extrema[MaxDimension];
965 
966  ssize_t
967  count,
968  *histogram[MaxDimension],
969  y;
970 
971  /*
972  Allocate histogram and extrema.
973  */
974  assert(image != (Image *) NULL);
975  assert(image->signature == MagickCoreSignature);
976  if (image->debug != MagickFalse)
977  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
978  GetPixelInfo(image,pixel);
979  for (i=0; i < MaxDimension; i++)
980  {
981  histogram[i]=(ssize_t *) AcquireQuantumMemory(256UL,sizeof(**histogram));
982  extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
983  if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
984  {
985  for (i-- ; i >= 0; i--)
986  {
987  extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
988  histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
989  }
990  (void) ThrowMagickException(exception,GetMagickModule(),
991  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
992  return(MagickFalse);
993  }
994  }
995  /*
996  Initialize histogram.
997  */
998  InitializeHistogram(image,histogram,exception);
999  (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
1000  (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
1001  (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
1002  (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
1003  (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
1004  (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
1005  /*
1006  Form clusters.
1007  */
1008  cluster=(Cluster *) NULL;
1009  head=(Cluster *) NULL;
1010  (void) memset(&red,0,sizeof(red));
1011  (void) memset(&green,0,sizeof(green));
1012  (void) memset(&blue,0,sizeof(blue));
1013  while (DefineRegion(extrema[Red],&red) != 0)
1014  {
1015  green.index=0;
1016  while (DefineRegion(extrema[Green],&green) != 0)
1017  {
1018  blue.index=0;
1019  while (DefineRegion(extrema[Blue],&blue) != 0)
1020  {
1021  /*
1022  Allocate a new class.
1023  */
1024  if (head != (Cluster *) NULL)
1025  {
1026  cluster->next=(Cluster *) AcquireQuantumMemory(1,
1027  sizeof(*cluster->next));
1028  cluster=cluster->next;
1029  }
1030  else
1031  {
1032  cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));
1033  head=cluster;
1034  }
1035  if (cluster == (Cluster *) NULL)
1036  {
1037  (void) ThrowMagickException(exception,GetMagickModule(),
1038  ResourceLimitError,"MemoryAllocationFailed","`%s'",
1039  image->filename);
1040  return(MagickFalse);
1041  }
1042  /*
1043  Initialize a new class.
1044  */
1045  cluster->count=0;
1046  cluster->red=red;
1047  cluster->green=green;
1048  cluster->blue=blue;
1049  cluster->next=(Cluster *) NULL;
1050  }
1051  }
1052  }
1053  if (head == (Cluster *) NULL)
1054  {
1055  /*
1056  No classes were identified-- create one.
1057  */
1058  cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));
1059  if (cluster == (Cluster *) NULL)
1060  {
1061  (void) ThrowMagickException(exception,GetMagickModule(),
1062  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1063  return(MagickFalse);
1064  }
1065  /*
1066  Initialize a new class.
1067  */
1068  cluster->count=0;
1069  cluster->red=red;
1070  cluster->green=green;
1071  cluster->blue=blue;
1072  cluster->next=(Cluster *) NULL;
1073  head=cluster;
1074  }
1075  /*
1076  Count the pixels for each cluster.
1077  */
1078  count=0;
1079  for (y=0; y < (ssize_t) image->rows; y++)
1080  {
1081  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1082  if (p == (const Quantum *) NULL)
1083  break;
1084  for (x=0; x < (ssize_t) image->columns; x++)
1085  {
1086  double
1087  b,
1088  g,
1089  r;
1090 
1091  r=(double) ScaleQuantumToChar(GetPixelRed(image,p));
1092  g=(double) ScaleQuantumToChar(GetPixelGreen(image,p));
1093  b=(double) ScaleQuantumToChar(GetPixelBlue(image,p));
1094  for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1095  if ((r >= (double) (cluster->red.left-SafeMargin)) &&
1096  (r <= (double) (cluster->red.right+SafeMargin)) &&
1097  (g >= (double) (cluster->green.left-SafeMargin)) &&
1098  (g <= (double) (cluster->green.right+SafeMargin)) &&
1099  (b >= (double) (cluster->blue.left-SafeMargin)) &&
1100  (b <= (double) (cluster->blue.right+SafeMargin)))
1101  {
1102  /*
1103  Count this pixel.
1104  */
1105  count++;
1106  cluster->red.center+=r;
1107  cluster->green.center+=g;
1108  cluster->blue.center+=b;
1109  cluster->count++;
1110  break;
1111  }
1112  p+=GetPixelChannels(image);
1113  }
1115  2*image->rows);
1116  if (proceed == MagickFalse)
1117  break;
1118  }
1119  /*
1120  Remove clusters that do not meet minimum cluster threshold.
1121  */
1122  count=0;
1123  last_cluster=head;
1124  next_cluster=head;
1125  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1126  {
1127  next_cluster=cluster->next;
1128  if ((cluster->count > 0) &&
1129  (cluster->count >= (count*cluster_threshold/100.0)))
1130  {
1131  /*
1132  Initialize cluster.
1133  */
1134  cluster->id=count;
1135  cluster->red.center/=cluster->count;
1136  cluster->green.center/=cluster->count;
1137  cluster->blue.center/=cluster->count;
1138  count++;
1139  last_cluster=cluster;
1140  continue;
1141  }
1142  /*
1143  Delete cluster.
1144  */
1145  if (cluster == head)
1146  head=next_cluster;
1147  else
1148  last_cluster->next=next_cluster;
1149  cluster=(Cluster *) RelinquishMagickMemory(cluster);
1150  }
1151  object=head;
1152  background=head;
1153  if (count > 1)
1154  {
1155  object=head->next;
1156  for (cluster=object; cluster->next != (Cluster *) NULL; )
1157  {
1158  if (cluster->count < object->count)
1159  object=cluster;
1160  cluster=cluster->next;
1161  }
1162  background=head->next;
1163  for (cluster=background; cluster->next != (Cluster *) NULL; )
1164  {
1165  if (cluster->count > background->count)
1166  background=cluster;
1167  cluster=cluster->next;
1168  }
1169  }
1170  if (background != (Cluster *) NULL)
1171  {
1172  threshold=(background->red.center+object->red.center)/2.0;
1173  pixel->red=(double) ScaleCharToQuantum((unsigned char)
1174  (threshold+0.5));
1175  threshold=(background->green.center+object->green.center)/2.0;
1176  pixel->green=(double) ScaleCharToQuantum((unsigned char)
1177  (threshold+0.5));
1178  threshold=(background->blue.center+object->blue.center)/2.0;
1179  pixel->blue=(double) ScaleCharToQuantum((unsigned char)
1180  (threshold+0.5));
1181  }
1182  /*
1183  Relinquish resources.
1184  */
1185  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1186  {
1187  next_cluster=cluster->next;
1188  cluster=(Cluster *) RelinquishMagickMemory(cluster);
1189  }
1190  for (i=0; i < MaxDimension; i++)
1191  {
1192  extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1193  histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1194  }
1195  return(MagickTrue);
1196 }
1197 
1198 /*
1199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1200 % %
1201 % %
1202 % %
1203 + I n i t i a l i z e H i s t o g r a m %
1204 % %
1205 % %
1206 % %
1207 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1208 %
1209 % InitializeHistogram() computes the histogram for an image.
1210 %
1211 % The format of the InitializeHistogram method is:
1212 %
1213 % InitializeHistogram(const Image *image,ssize_t **histogram)
1214 %
1215 % A description of each parameter follows.
1216 %
1217 % o image: Specifies a pointer to an Image structure; returned from
1218 % ReadImage.
1219 %
1220 % o histogram: Specifies an array of integers representing the number
1221 % of pixels for each intensity of a particular color component.
1222 %
1223 */
1224 static void InitializeHistogram(const Image *image,ssize_t **histogram,
1225  ExceptionInfo *exception)
1226 {
1227  const Quantum
1228  *p;
1229 
1230  ssize_t
1231  i,
1232  x;
1233 
1234  ssize_t
1235  y;
1236 
1237  /*
1238  Initialize histogram.
1239  */
1240  for (i=0; i <= 255; i++)
1241  {
1242  histogram[Red][i]=0;
1243  histogram[Green][i]=0;
1244  histogram[Blue][i]=0;
1245  }
1246  for (y=0; y < (ssize_t) image->rows; y++)
1247  {
1248  p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1249  if (p == (const Quantum *) NULL)
1250  break;
1251  for (x=0; x < (ssize_t) image->columns; x++)
1252  {
1253  histogram[Red][(ssize_t) ScaleQuantumToChar(GetPixelRed(image,p))]++;
1254  histogram[Green][(ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p))]++;
1255  histogram[Blue][(ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p))]++;
1256  p+=GetPixelChannels(image);
1257  }
1258  }
1259 }
1260 
1261 /*
1262 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1263 % %
1264 % %
1265 % %
1266 + I n i t i a l i z e I n t e r v a l T r e e %
1267 % %
1268 % %
1269 % %
1270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1271 %
1272 % InitializeIntervalTree() initializes an interval tree from the lists of
1273 % zero crossings.
1274 %
1275 % The format of the InitializeIntervalTree method is:
1276 %
1277 % InitializeIntervalTree(IntervalTree **list,ssize_t *number_nodes,
1278 % IntervalTree *node)
1279 %
1280 % A description of each parameter follows.
1281 %
1282 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
1283 %
1284 % o number_crossings: This size_t specifies the number of elements
1285 % in the zero_crossing array.
1286 %
1287 */
1288 
1289 static void InitializeList(IntervalTree **list,ssize_t *number_nodes,
1290  IntervalTree *node)
1291 {
1292  if (node == (IntervalTree *) NULL)
1293  return;
1294  if (node->child == (IntervalTree *) NULL)
1295  list[(*number_nodes)++]=node;
1296  InitializeList(list,number_nodes,node->sibling);
1297  InitializeList(list,number_nodes,node->child);
1298 }
1299 
1300 static void MeanStability(IntervalTree *node)
1301 {
1302  IntervalTree
1303  *child;
1304 
1305  if (node == (IntervalTree *) NULL)
1306  return;
1307  node->mean_stability=0.0;
1308  child=node->child;
1309  if (child != (IntervalTree *) NULL)
1310  {
1311  ssize_t
1312  count;
1313 
1314  double
1315  sum;
1316 
1317  sum=0.0;
1318  count=0;
1319  for ( ; child != (IntervalTree *) NULL; child=child->sibling)
1320  {
1321  sum+=child->stability;
1322  count++;
1323  }
1324  node->mean_stability=sum/(double) count;
1325  }
1326  MeanStability(node->sibling);
1327  MeanStability(node->child);
1328 }
1329 
1330 static void Stability(IntervalTree *node)
1331 {
1332  if (node == (IntervalTree *) NULL)
1333  return;
1334  if (node->child == (IntervalTree *) NULL)
1335  node->stability=0.0;
1336  else
1337  node->stability=node->tau-(node->child)->tau;
1338  Stability(node->sibling);
1339  Stability(node->child);
1340 }
1341 
1342 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing,
1343  const size_t number_crossings)
1344 {
1345  IntervalTree
1346  *head,
1347  **list,
1348  *node,
1349  *root;
1350 
1351  ssize_t
1352  i;
1353 
1354  ssize_t
1355  j,
1356  k,
1357  left,
1358  number_nodes;
1359 
1360  /*
1361  Allocate interval tree.
1362  */
1363  list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1364  sizeof(*list));
1365  if (list == (IntervalTree **) NULL)
1366  return((IntervalTree *) NULL);
1367  /*
1368  The root is the entire histogram.
1369  */
1370  root=(IntervalTree *) AcquireCriticalMemory(sizeof(*root));
1371  root->child=(IntervalTree *) NULL;
1372  root->sibling=(IntervalTree *) NULL;
1373  root->tau=0.0;
1374  root->left=0;
1375  root->right=255;
1376  root->mean_stability=0.0;
1377  root->stability=0.0;
1378  (void) memset(list,0,TreeLength*sizeof(*list));
1379  for (i=(-1); i < (ssize_t) number_crossings; i++)
1380  {
1381  /*
1382  Initialize list with all nodes with no children.
1383  */
1384  number_nodes=0;
1385  InitializeList(list,&number_nodes,root);
1386  /*
1387  Split list.
1388  */
1389  for (j=0; j < number_nodes; j++)
1390  {
1391  head=list[j];
1392  left=head->left;
1393  node=head;
1394  for (k=head->left+1; k < head->right; k++)
1395  {
1396  if (zero_crossing[i+1].crossings[k] != 0)
1397  {
1398  if (node == head)
1399  {
1401  sizeof(*node->child));
1402  node=node->child;
1403  }
1404  else
1405  {
1407  sizeof(*node->sibling));
1408  node=node->sibling;
1409  }
1410  if (node == (IntervalTree *) NULL)
1411  {
1412  list=(IntervalTree **) RelinquishMagickMemory(list);
1413  FreeNodes(root);
1414  return((IntervalTree *) NULL);
1415  }
1416  node->tau=zero_crossing[i+1].tau;
1417  node->child=(IntervalTree *) NULL;
1418  node->sibling=(IntervalTree *) NULL;
1419  node->left=left;
1420  node->right=k;
1421  left=k;
1422  }
1423  }
1424  if (left != head->left)
1425  {
1427  sizeof(*node->sibling));
1428  node=node->sibling;
1429  if (node == (IntervalTree *) NULL)
1430  {
1431  list=(IntervalTree **) RelinquishMagickMemory(list);
1432  FreeNodes(root);
1433  return((IntervalTree *) NULL);
1434  }
1435  node->tau=zero_crossing[i+1].tau;
1436  node->child=(IntervalTree *) NULL;
1437  node->sibling=(IntervalTree *) NULL;
1438  node->left=left;
1439  node->right=head->right;
1440  }
1441  }
1442  }
1443  /*
1444  Determine the stability: difference between a nodes tau and its child.
1445  */
1446  Stability(root->child);
1447  MeanStability(root->child);
1448  list=(IntervalTree **) RelinquishMagickMemory(list);
1449  return(root);
1450 }
1451 
1452 /*
1453 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1454 % %
1455 % %
1456 % %
1457 + O p t i m a l T a u %
1458 % %
1459 % %
1460 % %
1461 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1462 %
1463 % OptimalTau() finds the optimal tau for each band of the histogram.
1464 %
1465 % The format of the OptimalTau method is:
1466 %
1467 % double OptimalTau(const ssize_t *histogram,const double max_tau,
1468 % const double min_tau,const double delta_tau,
1469 % const double smooth_threshold,short *extrema)
1470 %
1471 % A description of each parameter follows.
1472 %
1473 % o histogram: Specifies an array of integers representing the number
1474 % of pixels for each intensity of a particular color component.
1475 %
1476 % o extrema: Specifies a pointer to an array of integers. They
1477 % represent the peaks and valleys of the histogram for each color
1478 % component.
1479 %
1480 */
1481 
1482 static void ActiveNodes(IntervalTree **list,ssize_t *number_nodes,
1483  IntervalTree *node)
1484 {
1485  if (node == (IntervalTree *) NULL)
1486  return;
1487  if (node->stability >= node->mean_stability)
1488  {
1489  list[(*number_nodes)++]=node;
1490  ActiveNodes(list,number_nodes,node->sibling);
1491  }
1492  else
1493  {
1494  ActiveNodes(list,number_nodes,node->sibling);
1495  ActiveNodes(list,number_nodes,node->child);
1496  }
1497 }
1498 
1499 static void FreeNodes(IntervalTree *node)
1500 {
1501  if (node == (IntervalTree *) NULL)
1502  return;
1503  FreeNodes(node->sibling);
1504  FreeNodes(node->child);
1505  node=(IntervalTree *) RelinquishMagickMemory(node);
1506 }
1507 
1508 static double OptimalTau(const ssize_t *histogram,const double max_tau,
1509  const double min_tau,const double delta_tau,const double smooth_threshold,
1510  short *extrema)
1511 {
1512  double
1513  average_tau,
1514  *derivative,
1515  *second_derivative,
1516  tau,
1517  value;
1518 
1519  IntervalTree
1520  **list,
1521  *node,
1522  *root;
1523 
1525  peak;
1526 
1527  ssize_t
1528  i,
1529  x;
1530 
1531  size_t
1532  count,
1533  number_crossings;
1534 
1535  ssize_t
1536  index,
1537  j,
1538  k,
1539  number_nodes;
1540 
1541  ZeroCrossing
1542  *zero_crossing;
1543 
1544  /*
1545  Allocate interval tree.
1546  */
1547  list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1548  sizeof(*list));
1549  if (list == (IntervalTree **) NULL)
1550  return(0.0);
1551  /*
1552  Allocate zero crossing list.
1553  */
1554  count=(size_t) ((max_tau-min_tau)/delta_tau)+2;
1555  zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count,
1556  sizeof(*zero_crossing));
1557  if (zero_crossing == (ZeroCrossing *) NULL)
1558  {
1559  list=(IntervalTree **) RelinquishMagickMemory(list);
1560  return(0.0);
1561  }
1562  for (i=0; i < (ssize_t) count; i++)
1563  zero_crossing[i].tau=(-1.0);
1564  /*
1565  Initialize zero crossing list.
1566  */
1567  derivative=(double *) AcquireCriticalMemory(256*sizeof(*derivative));
1568  second_derivative=(double *) AcquireCriticalMemory(256*
1569  sizeof(*second_derivative));
1570  i=0;
1571  for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
1572  {
1573  zero_crossing[i].tau=tau;
1574  ScaleSpace(histogram,tau,zero_crossing[i].histogram);
1575  DerivativeHistogram(zero_crossing[i].histogram,derivative);
1576  DerivativeHistogram(derivative,second_derivative);
1577  ZeroCrossHistogram(second_derivative,smooth_threshold,
1578  zero_crossing[i].crossings);
1579  i++;
1580  }
1581  /*
1582  Add an entry for the original histogram.
1583  */
1584  zero_crossing[i].tau=0.0;
1585  for (j=0; j <= 255; j++)
1586  zero_crossing[i].histogram[j]=(double) histogram[j];
1587  DerivativeHistogram(zero_crossing[i].histogram,derivative);
1588  DerivativeHistogram(derivative,second_derivative);
1589  ZeroCrossHistogram(second_derivative,smooth_threshold,
1590  zero_crossing[i].crossings);
1591  number_crossings=(size_t) i;
1592  derivative=(double *) RelinquishMagickMemory(derivative);
1593  second_derivative=(double *) RelinquishMagickMemory(second_derivative);
1594  /*
1595  Ensure the scale-space fingerprints form lines in scale-space, not loops.
1596  */
1597  ConsolidateCrossings(zero_crossing,number_crossings);
1598  /*
1599  Force endpoints to be included in the interval.
1600  */
1601  for (i=0; i <= (ssize_t) number_crossings; i++)
1602  {
1603  for (j=0; j < 255; j++)
1604  if (zero_crossing[i].crossings[j] != 0)
1605  break;
1606  zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]);
1607  for (j=255; j > 0; j--)
1608  if (zero_crossing[i].crossings[j] != 0)
1609  break;
1610  zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1611  }
1612  /*
1613  Initialize interval tree.
1614  */
1615  root=InitializeIntervalTree(zero_crossing,number_crossings);
1616  if (root == (IntervalTree *) NULL)
1617  {
1618  zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1619  list=(IntervalTree **) RelinquishMagickMemory(list);
1620  return(0.0);
1621  }
1622  /*
1623  Find active nodes: Stability is greater (or equal) to the mean stability of
1624  its children.
1625  */
1626  number_nodes=0;
1627  ActiveNodes(list,&number_nodes,root->child);
1628  /*
1629  Initialize extrema.
1630  */
1631  for (i=0; i <= 255; i++)
1632  extrema[i]=0;
1633  for (i=0; i < number_nodes; i++)
1634  {
1635  /*
1636  Find this tau in zero crossings list.
1637  */
1638  k=0;
1639  node=list[i];
1640  for (j=0; j <= (ssize_t) number_crossings; j++)
1641  if (zero_crossing[j].tau == node->tau)
1642  k=j;
1643  /*
1644  Find the value of the peak.
1645  */
1646  peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1647  MagickFalse;
1648  index=node->left;
1649  value=zero_crossing[k].histogram[index];
1650  for (x=node->left; x <= node->right; x++)
1651  {
1652  if (peak != MagickFalse)
1653  {
1654  if (zero_crossing[k].histogram[x] > value)
1655  {
1656  value=zero_crossing[k].histogram[x];
1657  index=x;
1658  }
1659  }
1660  else
1661  if (zero_crossing[k].histogram[x] < value)
1662  {
1663  value=zero_crossing[k].histogram[x];
1664  index=x;
1665  }
1666  }
1667  for (x=node->left; x <= node->right; x++)
1668  {
1669  if (index == 0)
1670  index=256;
1671  if (peak != MagickFalse)
1672  extrema[x]=(short) index;
1673  else
1674  extrema[x]=(short) (-index);
1675  }
1676  }
1677  /*
1678  Determine the average tau.
1679  */
1680  average_tau=0.0;
1681  for (i=0; i < number_nodes; i++)
1682  average_tau+=list[i]->tau;
1683  average_tau*=PerceptibleReciprocal((double) number_nodes);
1684  /*
1685  Relinquish resources.
1686  */
1687  FreeNodes(root);
1688  zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1689  list=(IntervalTree **) RelinquishMagickMemory(list);
1690  return(average_tau);
1691 }
1692 
1693 /*
1694 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1695 % %
1696 % %
1697 % %
1698 + S c a l e S p a c e %
1699 % %
1700 % %
1701 % %
1702 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1703 %
1704 % ScaleSpace() performs a scale-space filter on the 1D histogram.
1705 %
1706 % The format of the ScaleSpace method is:
1707 %
1708 % ScaleSpace(const ssize_t *histogram,const double tau,
1709 % double *scale_histogram)
1710 %
1711 % A description of each parameter follows.
1712 %
1713 % o histogram: Specifies an array of doubles representing the number
1714 % of pixels for each intensity of a particular color component.
1715 %
1716 */
1717 static void ScaleSpace(const ssize_t *histogram,const double tau,
1718  double *scale_histogram)
1719 {
1720  double
1721  alpha,
1722  beta,
1723  *gamma,
1724  sum;
1725 
1726  ssize_t
1727  u,
1728  x;
1729 
1730  gamma=(double *) AcquireQuantumMemory(256,sizeof(*gamma));
1731  if (gamma == (double *) NULL)
1732  ThrowFatalException(ResourceLimitFatalError,"UnableToAllocateGammaMap");
1733  alpha=PerceptibleReciprocal(tau*sqrt(2.0*MagickPI));
1734  beta=(-1.0*PerceptibleReciprocal(2.0*tau*tau));
1735  for (x=0; x <= 255; x++)
1736  gamma[x]=0.0;
1737  for (x=0; x <= 255; x++)
1738  {
1739  gamma[x]=exp((double) beta*x*x);
1740  if (gamma[x] < MagickEpsilon)
1741  break;
1742  }
1743  for (x=0; x <= 255; x++)
1744  {
1745  sum=0.0;
1746  for (u=0; u <= 255; u++)
1747  sum+=(double) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1748  scale_histogram[x]=alpha*sum;
1749  }
1750  gamma=(double *) RelinquishMagickMemory(gamma);
1751 }
1752 
1753 /*
1754 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1755 % %
1756 % %
1757 % %
1758 % S e g m e n t I m a g e %
1759 % %
1760 % %
1761 % %
1762 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1763 %
1764 % SegmentImage() segment an image by analyzing the histograms of the color
1765 % components and identifying units that are homogeneous with the fuzzy
1766 % C-means technique.
1767 %
1768 % The format of the SegmentImage method is:
1769 %
1770 % MagickBooleanType SegmentImage(Image *image,
1771 % const ColorspaceType colorspace,const MagickBooleanType verbose,
1772 % const double cluster_threshold,const double smooth_threshold,
1773 % ExceptionInfo *exception)
1774 %
1775 % A description of each parameter follows.
1776 %
1777 % o image: the image.
1778 %
1779 % o colorspace: Indicate the colorspace.
1780 %
1781 % o verbose: Set to MagickTrue to print detailed information about the
1782 % identified classes.
1783 %
1784 % o cluster_threshold: This represents the minimum number of pixels
1785 % contained in a hexahedra before it can be considered valid (expressed
1786 % as a percentage).
1787 %
1788 % o smooth_threshold: the smoothing threshold eliminates noise in the second
1789 % derivative of the histogram. As the value is increased, you can expect a
1790 % smoother second derivative.
1791 %
1792 % o exception: return any errors or warnings in this structure.
1793 %
1794 */
1796  const ColorspaceType colorspace,const MagickBooleanType verbose,
1797  const double cluster_threshold,const double smooth_threshold,
1798  ExceptionInfo *exception)
1799 {
1801  previous_colorspace;
1802 
1804  status;
1805 
1806  ssize_t
1807  i;
1808 
1809  short
1810  *extrema[MaxDimension];
1811 
1812  ssize_t
1813  *histogram[MaxDimension];
1814 
1815  /*
1816  Allocate histogram and extrema.
1817  */
1818  assert(image != (Image *) NULL);
1819  assert(image->signature == MagickCoreSignature);
1820  if (image->debug != MagickFalse)
1821  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1822  for (i=0; i < MaxDimension; i++)
1823  {
1824  histogram[i]=(ssize_t *) AcquireQuantumMemory(256,sizeof(**histogram));
1825  extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1826  if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1827  {
1828  for (i-- ; i >= 0; i--)
1829  {
1830  extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1831  histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1832  }
1833  ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1834  image->filename)
1835  }
1836  }
1837  /*
1838  Initialize histogram.
1839  */
1840  previous_colorspace=image->colorspace;
1841  (void) TransformImageColorspace(image,colorspace,exception);
1842  InitializeHistogram(image,histogram,exception);
1843  (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,smooth_threshold == 0.0 ?
1844  1.0 : smooth_threshold,extrema[Red]);
1845  (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,smooth_threshold == 0.0 ?
1846  1.0 : smooth_threshold,extrema[Green]);
1847  (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,smooth_threshold == 0.0 ?
1848  1.0 : smooth_threshold,extrema[Blue]);
1849  /*
1850  Classify using the fuzzy c-Means technique.
1851  */
1852  status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose,
1853  exception);
1854  (void) TransformImageColorspace(image,previous_colorspace,exception);
1855  /*
1856  Relinquish resources.
1857  */
1858  for (i=0; i < MaxDimension; i++)
1859  {
1860  extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1861  histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1862  }
1863  return(status);
1864 }
1865 
1866 /*
1867 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1868 % %
1869 % %
1870 % %
1871 + Z e r o C r o s s H i s t o g r a m %
1872 % %
1873 % %
1874 % %
1875 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1876 %
1877 % ZeroCrossHistogram() find the zero crossings in a histogram and marks
1878 % directions as: 1 is negative to positive; 0 is zero crossing; and -1
1879 % is positive to negative.
1880 %
1881 % The format of the ZeroCrossHistogram method is:
1882 %
1883 % ZeroCrossHistogram(double *second_derivative,
1884 % const double smooth_threshold,short *crossings)
1885 %
1886 % A description of each parameter follows.
1887 %
1888 % o second_derivative: Specifies an array of doubles representing the
1889 % second derivative of the histogram of a particular color component.
1890 %
1891 % o crossings: This array of integers is initialized with
1892 % -1, 0, or 1 representing the slope of the first derivative of the
1893 % of a particular color component.
1894 %
1895 */
1896 static void ZeroCrossHistogram(double *second_derivative,
1897  const double smooth_threshold,short *crossings)
1898 {
1899  ssize_t
1900  i;
1901 
1902  ssize_t
1903  parity;
1904 
1905  /*
1906  Merge low numbers to zero to help prevent noise.
1907  */
1908  for (i=0; i <= 255; i++)
1909  if ((second_derivative[i] < smooth_threshold) &&
1910  (second_derivative[i] >= -smooth_threshold))
1911  second_derivative[i]=0.0;
1912  /*
1913  Mark zero crossings.
1914  */
1915  parity=0;
1916  for (i=0; i <= 255; i++)
1917  {
1918  crossings[i]=0;
1919  if (second_derivative[i] < 0.0)
1920  {
1921  if (parity > 0)
1922  crossings[i]=(-1);
1923  parity=1;
1924  }
1925  else
1926  if (second_derivative[i] > 0.0)
1927  {
1928  if (parity < 0)
1929  crossings[i]=1;
1930  parity=(-1);
1931  }
1932  }
1933 }
size_t rows
Definition: image.h:172
#define magick_restrict
Definition: MagickCore.h:41
struct _IntervalTree IntervalTree
ExtentPacket blue
Definition: segment.c:142
ExtentPacket green
Definition: segment.c:142
MagickExport CacheView * DestroyCacheView(CacheView *cache_view)
Definition: cache-view.c:252
static void ConsolidateCrossings(ZeroCrossing *zero_crossing, const size_t number_crossings)
Definition: segment.c:703
struct _ExtentPacket ExtentPacket
ssize_t right
Definition: segment.c:157
PixelInfo * colormap
Definition: image.h:179
ExtentPacket red
Definition: segment.c:142
MagickProgressMonitor progress_monitor
Definition: image.h:303
struct _IntervalTree * sibling
Definition: segment.c:164
MagickExport MagickBooleanType SyncImage(Image *image, ExceptionInfo *exception)
Definition: image.c:3899
MagickExport MagickBooleanType TransformImageColorspace(Image *image, const ColorspaceType colorspace, ExceptionInfo *exception)
Definition: colorspace.c:1607
ssize_t right
Definition: segment.c:131
static Quantum GetPixelRed(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
#define ThrowFatalException(severity, tag)
#define WeightingExponent
Definition: segment.c:117
double mean_stability
Definition: segment.c:161
ssize_t index
Definition: segment.c:131
MagickRealType red
Definition: pixel.h:193
#define MagickPI
Definition: image-private.h:42
#define MagickAbsoluteValue(x)
Definition: image-private.h:37
struct _Cluster * next
Definition: segment.c:138
#define Tau
Definition: segment.c:120
MagickExport const Quantum * GetCacheViewVirtualPixels(const CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:651
ssize_t count
Definition: segment.c:147
#define MagickEpsilon
Definition: magick-type.h:114
static void Stability(IntervalTree *node)
Definition: segment.c:1330
#define ThrowBinaryException(severity, tag, context)
Definition: log.h:52
ssize_t MagickOffsetType
Definition: magick-type.h:133
MagickExport void GetPixelInfo(const Image *image, PixelInfo *pixel)
Definition: pixel.c:2170
static void ScaleSpace(const ssize_t *, const double, double *)
Definition: segment.c:1717
Definition: image.h:151
static IntervalTree * InitializeIntervalTree(const ZeroCrossing *zero_crossing, const size_t number_crossings)
Definition: segment.c:1342
static const int Green
Definition: segment.c:184
MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image, const double cluster_threshold, const double smooth_threshold, PixelInfo *pixel, ExceptionInfo *exception)
Definition: segment.c:933
#define MagickCoreSignature
MagickExport Quantum * GetCacheViewAuthenticPixels(CacheView *cache_view, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache-view.c:299
static void MeanStability(IntervalTree *node)
Definition: segment.c:1300
double histogram[256]
Definition: segment.c:172
static void InitializeList(IntervalTree **list, ssize_t *number_nodes, IntervalTree *node)
Definition: segment.c:1289
MagickExport ssize_t FormatLocaleFile(FILE *file, const char *magick_restrict format,...)
Definition: locale.c:370
MagickBooleanType
Definition: magick-type.h:161
unsigned int MagickStatusType
Definition: magick-type.h:125
static double PerceptibleReciprocal(const double x)
MagickExport void * AcquireCriticalMemory(const size_t size)
Definition: memory.c:626
struct _ZeroCrossing ZeroCrossing
double center
Definition: segment.c:128
MagickExport MagickBooleanType SegmentImage(Image *image, const ColorspaceType colorspace, const MagickBooleanType verbose, const double cluster_threshold, const double smooth_threshold, ExceptionInfo *exception)
Definition: segment.c:1795
MagickExport const Quantum * GetVirtualPixels(const Image *image, const ssize_t x, const ssize_t y, const size_t columns, const size_t rows, ExceptionInfo *exception)
Definition: cache.c:3252
MagickExport void * AcquireQuantumMemory(const size_t count, const size_t quantum)
Definition: memory.c:665
struct _Cluster Cluster
#define SegmentImageTag
static Quantum GetPixelGreen(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
static void ActiveNodes(IntervalTree **list, ssize_t *number_nodes, IntervalTree *node)
Definition: segment.c:1482
MagickRealType blue
Definition: pixel.h:193
short crossings[256]
Definition: segment.c:176
MagickExport MagickBooleanType ThrowMagickException(ExceptionInfo *exception, const char *module, const char *function, const size_t line, const ExceptionType severity, const char *tag, const char *format,...)
Definition: exception.c:1145
MagickExport MagickBooleanType LogMagickEvent(const LogEventType type, const char *module, const char *function, const size_t line, const char *format,...)
Definition: log.c:1662
size_t signature
Definition: image.h:354
size_t columns
Definition: image.h:172
static double OptimalTau(const ssize_t *, const double, const double, const double, const double, short *)
Definition: segment.c:1508
MagickBooleanType(* MagickProgressMonitor)(const char *, const MagickOffsetType, const MagickSizeType, void *)
Definition: monitor.h:26
static void ZeroCrossHistogram(double *, const double, short *)
Definition: segment.c:1896
static ssize_t DefineRegion(const short *, ExtentPacket *)
Definition: segment.c:820
#define MagickMax(x, y)
Definition: image-private.h:38
static MagickBooleanType Classify(Image *image, short **extrema, const double cluster_threshold, const double weighting_exponent, const MagickBooleanType verbose, ExceptionInfo *exception)
Definition: segment.c:246
size_t colors
Definition: image.h:172
static const int SafeMargin
Definition: segment.c:186
static size_t GetPixelChannels(const Image *magick_restrict image)
MagickExport MagickBooleanType AcquireImageColormap(Image *image, const size_t colors, ExceptionInfo *exception)
Definition: colormap.c:105
char filename[MagickPathExtent]
Definition: image.h:319
#define SegmentPower(ratio)
Definition: segment.c:118
#define GetMagickModule()
Definition: log.h:28
ssize_t id
Definition: segment.c:147
MagickExport CacheView * AcquireVirtualCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:149
static void DerivativeHistogram(const double *histogram, double *derivative)
Definition: segment.c:876
unsigned short Quantum
Definition: magick-type.h:86
static void FreeNodes(IntervalTree *)
Definition: segment.c:1499
static const int TreeLength
Definition: segment.c:187
#define DeltaTau
Definition: segment.c:112
static void SetPixelIndex(const Image *magick_restrict image, const Quantum index, Quantum *magick_restrict pixel)
static const int Red
Definition: segment.c:185
#define MagickMin(x, y)
Definition: image-private.h:39
#define ThrowClassifyException(severity, tag, label)
ColorspaceType
Definition: colorspace.h:25
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1162
ssize_t left
Definition: segment.c:131
MagickRealType green
Definition: pixel.h:193
static const int Blue
Definition: segment.c:183
ssize_t left
Definition: segment.c:157
#define MagickExport
MagickExport MagickBooleanType SyncCacheViewAuthenticPixels(CacheView *magick_restrict cache_view, ExceptionInfo *exception)
Definition: cache-view.c:1100
MagickExport CacheView * AcquireAuthenticCacheView(const Image *image, ExceptionInfo *exception)
Definition: cache-view.c:112
static void InitializeHistogram(const Image *, ssize_t **, ExceptionInfo *)
Definition: segment.c:1224
double tau
Definition: segment.c:154
double tau
Definition: segment.c:172
static Quantum GetPixelBlue(const Image *magick_restrict image, const Quantum *magick_restrict pixel)
struct _IntervalTree * child
Definition: segment.c:164
#define MaxDimension
Definition: segment.c:111
ColorspaceType colorspace
Definition: image.h:157
MagickExport MagickBooleanType SetImageProgress(const Image *image, const char *tag, const MagickOffsetType offset, const MagickSizeType extent)
Definition: monitor.c:136
double stability
Definition: segment.c:161
MagickBooleanType debug
Definition: image.h:334