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