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