2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %
13 % MagickCore Methods to Segment an Image with Thresholding Fuzzy c-Means %
20 % Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % http://www.imagemagick.org/script/license.php %
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. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
46 % The fuzzy c-Means algorithm can be summarized as follows:
48 % o Build a histogram, one for each color component of the image.
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
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
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
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.
72 % Segment is strongly based on software written by Andy Gallo,
73 % University of Delaware.
75 % The following reference was used in creating this program:
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
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"
90 #include "MagickCore/colorspace-private.h"
91 #include "MagickCore/exception.h"
92 #include "MagickCore/exception-private.h"
93 #include "MagickCore/image.h"
94 #include "MagickCore/image-private.h"
95 #include "MagickCore/memory_.h"
96 #include "MagickCore/monitor.h"
97 #include "MagickCore/monitor-private.h"
98 #include "MagickCore/pixel-accessor.h"
99 #include "MagickCore/quantize.h"
100 #include "MagickCore/quantum.h"
101 #include "MagickCore/quantum-private.h"
102 #include "MagickCore/segment.h"
103 #include "MagickCore/string_.h"
108 #define MaxDimension 3
109 #define DeltaTau 0.5f
110 #if defined(FastClassify)
111 #define WeightingExponent 2.0
112 #define SegmentPower(ratio) (ratio)
114 #define WeightingExponent 2.5
115 #define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0)));
120 Typedef declarations.
122 typedef struct _ExtentPacket
133 typedef struct _Cluster
148 typedef struct _IntervalTree
166 typedef struct _ZeroCrossing
177 Constant declarations.
189 static MagickRealType
190 OptimalTau(const ssize_t *,const double,const double,const double,
191 const double,short *);
194 DefineRegion(const short *,ExtentPacket *);
197 InitializeHistogram(const Image *,ssize_t **,ExceptionInfo *),
198 ScaleSpace(const ssize_t *,const MagickRealType,MagickRealType *),
199 ZeroCrossHistogram(MagickRealType *,const MagickRealType,short *);
202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
212 % Classify() defines one or more classes. Each pixel is thresholded to
213 % determine which class it belongs to. If the class is not identified it is
214 % assigned to the closest class based on the fuzzy c-Means technique.
216 % The format of the Classify method is:
218 % MagickBooleanType Classify(Image *image,short **extrema,
219 % const MagickRealType cluster_threshold,
220 % const MagickRealType weighting_exponent,
221 % const MagickBooleanType verbose,ExceptionInfo *exception)
223 % A description of each parameter follows.
225 % o image: the image.
227 % o extrema: Specifies a pointer to an array of integers. They
228 % represent the peaks and valleys of the histogram for each color
231 % o cluster_threshold: This MagickRealType represents the minimum number of
232 % pixels contained in a hexahedra before it can be considered valid
233 % (expressed as a percentage).
235 % o weighting_exponent: Specifies the membership weighting exponent.
237 % o verbose: A value greater than zero prints detailed information about
238 % the identified classes.
240 % o exception: return any errors or warnings in this structure.
243 static MagickBooleanType Classify(Image *image,short **extrema,
244 const MagickRealType cluster_threshold,
245 const MagickRealType weighting_exponent,const MagickBooleanType verbose,
246 ExceptionInfo *exception)
248 #define SegmentImageTag "Segment/Image"
276 register MagickRealType
289 cluster=(Cluster *) NULL;
290 head=(Cluster *) NULL;
291 (void) ResetMagickMemory(&red,0,sizeof(red));
292 (void) ResetMagickMemory(&green,0,sizeof(green));
293 (void) ResetMagickMemory(&blue,0,sizeof(blue));
294 while (DefineRegion(extrema[Red],&red) != 0)
297 while (DefineRegion(extrema[Green],&green) != 0)
300 while (DefineRegion(extrema[Blue],&blue) != 0)
303 Allocate a new class.
305 if (head != (Cluster *) NULL)
307 cluster->next=(Cluster *) AcquireMagickMemory(
308 sizeof(*cluster->next));
309 cluster=cluster->next;
313 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
316 if (cluster == (Cluster *) NULL)
317 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
320 Initialize a new class.
324 cluster->green=green;
326 cluster->next=(Cluster *) NULL;
330 if (head == (Cluster *) NULL)
333 No classes were identified-- create one.
335 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
336 if (cluster == (Cluster *) NULL)
337 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
340 Initialize a new class.
344 cluster->green=green;
346 cluster->next=(Cluster *) NULL;
350 Count the pixels for each cluster.
355 image_view=AcquireCacheView(image);
356 for (y=0; y < (ssize_t) image->rows; y++)
358 register const Quantum
364 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
365 if (p == (const Quantum *) NULL)
367 for (x=0; x < (ssize_t) image->columns; x++)
369 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
370 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >=
371 (cluster->red.left-SafeMargin)) &&
372 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <=
373 (cluster->red.right+SafeMargin)) &&
374 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >=
375 (cluster->green.left-SafeMargin)) &&
376 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <=
377 (cluster->green.right+SafeMargin)) &&
378 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >=
379 (cluster->blue.left-SafeMargin)) &&
380 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <=
381 (cluster->blue.right+SafeMargin)))
387 cluster->red.center+=(MagickRealType) ScaleQuantumToChar(
388 GetPixelRed(image,p));
389 cluster->green.center+=(MagickRealType) ScaleQuantumToChar(
390 GetPixelGreen(image,p));
391 cluster->blue.center+=(MagickRealType) ScaleQuantumToChar(
392 GetPixelBlue(image,p));
396 p+=GetPixelChannels(image);
398 if (image->progress_monitor != (MagickProgressMonitor) NULL)
403 #if defined(MAGICKCORE_OPENMP_SUPPORT)
404 #pragma omp critical (MagickCore_Classify)
406 proceed=SetImageProgress(image,SegmentImageTag,progress++,
408 if (proceed == MagickFalse)
412 image_view=DestroyCacheView(image_view);
414 Remove clusters that do not meet minimum cluster threshold.
419 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
421 next_cluster=cluster->next;
422 if ((cluster->count > 0) &&
423 (cluster->count >= (count*cluster_threshold/100.0)))
429 cluster->red.center/=cluster->count;
430 cluster->green.center/=cluster->count;
431 cluster->blue.center/=cluster->count;
433 last_cluster=cluster;
442 last_cluster->next=next_cluster;
443 cluster=(Cluster *) RelinquishMagickMemory(cluster);
445 number_clusters=(size_t) count;
446 if (verbose != MagickFalse)
449 Print cluster statistics.
451 (void) FormatLocaleFile(stdout,"Fuzzy C-means Statistics\n");
452 (void) FormatLocaleFile(stdout,"===================\n\n");
453 (void) FormatLocaleFile(stdout,"\tCluster Threshold = %g\n",(double)
455 (void) FormatLocaleFile(stdout,"\tWeighting Exponent = %g\n",(double)
457 (void) FormatLocaleFile(stdout,"\tTotal Number of Clusters = %.20g\n\n",
458 (double) number_clusters);
460 Print the total number of points per cluster.
462 (void) FormatLocaleFile(stdout,"\n\nNumber of Vectors Per Cluster\n");
463 (void) FormatLocaleFile(stdout,"=============================\n\n");
464 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
465 (void) FormatLocaleFile(stdout,"Cluster #%.20g = %.20g\n",(double)
466 cluster->id,(double) cluster->count);
468 Print the cluster extents.
470 (void) FormatLocaleFile(stdout,
471 "\n\n\nCluster Extents: (Vector Size: %d)\n",MaxDimension);
472 (void) FormatLocaleFile(stdout,"================");
473 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
475 (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
477 (void) FormatLocaleFile(stdout,
478 "%.20g-%.20g %.20g-%.20g %.20g-%.20g\n",(double)
479 cluster->red.left,(double) cluster->red.right,(double)
480 cluster->green.left,(double) cluster->green.right,(double)
481 cluster->blue.left,(double) cluster->blue.right);
484 Print the cluster center values.
486 (void) FormatLocaleFile(stdout,
487 "\n\n\nCluster Center Values: (Vector Size: %d)\n",MaxDimension);
488 (void) FormatLocaleFile(stdout,"=====================");
489 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
491 (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
493 (void) FormatLocaleFile(stdout,"%g %g %g\n",(double)
494 cluster->red.center,(double) cluster->green.center,(double)
495 cluster->blue.center);
497 (void) FormatLocaleFile(stdout,"\n");
499 if (number_clusters > 256)
500 ThrowBinaryException(ImageError,"TooManyClusters",image->filename);
502 Speed up distance calculations.
504 squares=(MagickRealType *) AcquireQuantumMemory(513UL,sizeof(*squares));
505 if (squares == (MagickRealType *) NULL)
506 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
509 for (i=(-255); i <= 255; i++)
510 squares[i]=(MagickRealType) i*(MagickRealType) i;
512 Allocate image colormap.
514 if (AcquireImageColormap(image,number_clusters,exception) == MagickFalse)
515 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
518 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
520 image->colormap[i].red=ScaleCharToQuantum((unsigned char)
521 (cluster->red.center+0.5));
522 image->colormap[i].green=ScaleCharToQuantum((unsigned char)
523 (cluster->green.center+0.5));
524 image->colormap[i].blue=ScaleCharToQuantum((unsigned char)
525 (cluster->blue.center+0.5));
529 Do course grain classes.
531 exception=(&image->exception);
532 image_view=AcquireCacheView(image);
533 #if defined(MAGICKCORE_OPENMP_SUPPORT)
534 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
536 for (y=0; y < (ssize_t) image->rows; y++)
541 register const PixelPacket
550 if (status == MagickFalse)
552 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
553 if (q == (Quantum *) NULL)
558 for (x=0; x < (ssize_t) image->columns; x++)
560 SetPixelIndex(image,0,q);
561 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
563 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,q)) >=
564 (cluster->red.left-SafeMargin)) &&
565 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,q)) <=
566 (cluster->red.right+SafeMargin)) &&
567 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q)) >=
568 (cluster->green.left-SafeMargin)) &&
569 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q)) <=
570 (cluster->green.right+SafeMargin)) &&
571 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q)) >=
572 (cluster->blue.left-SafeMargin)) &&
573 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q)) <=
574 (cluster->blue.right+SafeMargin)))
579 SetPixelIndex(image,(Quantum) cluster->id,q);
583 if (cluster == (Cluster *) NULL)
597 Compute fuzzy membership.
600 for (j=0; j < (ssize_t) image->colors; j++)
604 distance_squared=squares[(ssize_t) ScaleQuantumToChar(
605 GetPixelRed(image,q))-(ssize_t)
606 ScaleQuantumToChar(p->red)]+squares[(ssize_t)
607 ScaleQuantumToChar(GetPixelGreen(image,q))-(ssize_t)
608 ScaleQuantumToChar(p->green)]+squares[(ssize_t)
609 ScaleQuantumToChar(GetPixelBlue(image,q))-(ssize_t)
610 ScaleQuantumToChar(p->blue)];
611 numerator=distance_squared;
612 for (k=0; k < (ssize_t) image->colors; k++)
615 distance_squared=squares[(ssize_t) ScaleQuantumToChar(
616 GetPixelRed(image,q))-(ssize_t)
617 ScaleQuantumToChar(p->red)]+squares[(ssize_t)
618 ScaleQuantumToChar(GetPixelGreen(image,q))-(ssize_t)
619 ScaleQuantumToChar(p->green)]+squares[(ssize_t)
620 ScaleQuantumToChar(GetPixelBlue(image,q))-(ssize_t)
621 ScaleQuantumToChar(p->blue)];
622 ratio=numerator/distance_squared;
623 sum+=SegmentPower(ratio);
625 if ((sum != 0.0) && ((1.0/sum) > local_minima))
630 local_minima=1.0/sum;
631 SetPixelIndex(image,(Quantum) j,q);
635 q+=GetPixelChannels(image);
637 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
639 if (image->progress_monitor != (MagickProgressMonitor) NULL)
644 #if defined(MAGICKCORE_OPENMP_SUPPORT)
645 #pragma omp critical (MagickCore_Classify)
647 proceed=SetImageProgress(image,SegmentImageTag,progress++,
649 if (proceed == MagickFalse)
653 image_view=DestroyCacheView(image_view);
654 status&=SyncImage(image);
656 Relinquish resources.
658 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
660 next_cluster=cluster->next;
661 cluster=(Cluster *) RelinquishMagickMemory(cluster);
664 free_squares=squares;
665 free_squares=(MagickRealType *) RelinquishMagickMemory(free_squares);
670 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
674 + C o n s o l i d a t e C r o s s i n g s %
678 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
680 % ConsolidateCrossings() guarantees that an even number of zero crossings
681 % always lie between two crossings.
683 % The format of the ConsolidateCrossings method is:
685 % ConsolidateCrossings(ZeroCrossing *zero_crossing,
686 % const size_t number_crossings)
688 % A description of each parameter follows.
690 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
692 % o number_crossings: This size_t specifies the number of elements
693 % in the zero_crossing array.
697 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
704 static inline ssize_t MagickMax(const ssize_t x,const ssize_t y)
711 static inline ssize_t MagickMin(const ssize_t x,const ssize_t y)
718 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
719 const size_t number_crossings)
735 Consolidate zero crossings.
737 for (i=(ssize_t) number_crossings-1; i >= 0; i--)
738 for (j=0; j <= 255; j++)
740 if (zero_crossing[i].crossings[j] == 0)
743 Find the entry that is closest to j and still preserves the
744 property that there are an even number of crossings between
747 for (k=j-1; k > 0; k--)
748 if (zero_crossing[i+1].crossings[k] != 0)
752 for (k=j+1; k < 255; k++)
753 if (zero_crossing[i+1].crossings[k] != 0)
755 right=MagickMin(k,255);
757 K is the zero crossing just left of j.
759 for (k=j-1; k > 0; k--)
760 if (zero_crossing[i].crossings[k] != 0)
765 Check center for an even number of crossings between k and j.
768 if (zero_crossing[i+1].crossings[j] != 0)
771 for (l=k+1; l < center; l++)
772 if (zero_crossing[i+1].crossings[l] != 0)
774 if (((count % 2) == 0) && (center != k))
778 Check left for an even number of crossings between k and j.
783 for (l=k+1; l < left; l++)
784 if (zero_crossing[i+1].crossings[l] != 0)
786 if (((count % 2) == 0) && (left != k))
790 Check right for an even number of crossings between k and j.
795 for (l=k+1; l < right; l++)
796 if (zero_crossing[i+1].crossings[l] != 0)
798 if (((count % 2) == 0) && (right != k))
801 l=(ssize_t) zero_crossing[i].crossings[j];
802 zero_crossing[i].crossings[j]=0;
804 zero_crossing[i].crossings[correct]=(short) l;
809 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
813 + D e f i n e R e g i o n %
817 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
819 % DefineRegion() defines the left and right boundaries of a peak region.
821 % The format of the DefineRegion method is:
823 % ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
825 % A description of each parameter follows.
827 % o extrema: Specifies a pointer to an array of integers. They
828 % represent the peaks and valleys of the histogram for each color
831 % o extents: This pointer to an ExtentPacket represent the extends
832 % of a particular peak or valley of a color component.
835 static ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
838 Initialize to default values.
844 Find the left side (maxima).
846 for ( ; extents->index <= 255; extents->index++)
847 if (extrema[extents->index] > 0)
849 if (extents->index > 255)
850 return(MagickFalse); /* no left side - no region exists */
851 extents->left=extents->index;
853 Find the right side (minima).
855 for ( ; extents->index <= 255; extents->index++)
856 if (extrema[extents->index] < 0)
858 extents->right=extents->index-1;
863 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
867 + D e r i v a t i v e H i s t o g r a m %
871 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
873 % DerivativeHistogram() determines the derivative of the histogram using
874 % central differencing.
876 % The format of the DerivativeHistogram method is:
878 % DerivativeHistogram(const MagickRealType *histogram,
879 % MagickRealType *derivative)
881 % A description of each parameter follows.
883 % o histogram: Specifies an array of MagickRealTypes representing the number
884 % of pixels for each intensity of a particular color component.
886 % o derivative: This array of MagickRealTypes is initialized by
887 % DerivativeHistogram to the derivative of the histogram using central
891 static void DerivativeHistogram(const MagickRealType *histogram,
892 MagickRealType *derivative)
899 Compute endpoints using second order polynomial interpolation.
902 derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);
903 derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);
905 Compute derivative using central differencing.
907 for (i=1; i < n; i++)
908 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
913 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
917 + G e t I m a g e D y n a m i c T h r e s h o l d %
921 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
923 % GetImageDynamicThreshold() returns the dynamic threshold for an image.
925 % The format of the GetImageDynamicThreshold method is:
927 % MagickBooleanType GetImageDynamicThreshold(const Image *image,
928 % const double cluster_threshold,const double smooth_threshold,
929 % PixelInfo *pixel,ExceptionInfo *exception)
931 % A description of each parameter follows.
933 % o image: the image.
935 % o cluster_threshold: This MagickRealType represents the minimum number of
936 % pixels contained in a hexahedra before it can be considered valid
937 % (expressed as a percentage).
939 % o smooth_threshold: the smoothing threshold eliminates noise in the second
940 % derivative of the histogram. As the value is increased, you can expect a
941 % smoother second derivative.
943 % o pixel: return the dynamic threshold here.
945 % o exception: return any errors or warnings in this structure.
948 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,
949 const double cluster_threshold,const double smooth_threshold,
950 PixelInfo *pixel,ExceptionInfo *exception)
971 register const Quantum
979 *extrema[MaxDimension];
983 *histogram[MaxDimension],
987 Allocate histogram and extrema.
989 assert(image != (Image *) NULL);
990 assert(image->signature == MagickSignature);
991 if (image->debug != MagickFalse)
992 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
993 GetPixelInfo(image,pixel);
994 for (i=0; i < MaxDimension; i++)
996 histogram[i]=(ssize_t *) AcquireQuantumMemory(256UL,sizeof(**histogram));
997 extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
998 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1000 for (i-- ; i >= 0; i--)
1002 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1003 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1005 (void) ThrowMagickException(exception,GetMagickModule(),
1006 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1007 return(MagickFalse);
1011 Initialize histogram.
1013 InitializeHistogram(image,histogram,exception);
1014 (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
1015 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
1016 (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
1017 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
1018 (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
1019 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
1023 cluster=(Cluster *) NULL;
1024 head=(Cluster *) NULL;
1025 (void) ResetMagickMemory(&red,0,sizeof(red));
1026 (void) ResetMagickMemory(&green,0,sizeof(green));
1027 (void) ResetMagickMemory(&blue,0,sizeof(blue));
1028 while (DefineRegion(extrema[Red],&red) != 0)
1031 while (DefineRegion(extrema[Green],&green) != 0)
1034 while (DefineRegion(extrema[Blue],&blue) != 0)
1037 Allocate a new class.
1039 if (head != (Cluster *) NULL)
1041 cluster->next=(Cluster *) AcquireMagickMemory(
1042 sizeof(*cluster->next));
1043 cluster=cluster->next;
1047 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1050 if (cluster == (Cluster *) NULL)
1052 (void) ThrowMagickException(exception,GetMagickModule(),
1053 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1055 return(MagickFalse);
1058 Initialize a new class.
1062 cluster->green=green;
1064 cluster->next=(Cluster *) NULL;
1068 if (head == (Cluster *) NULL)
1071 No classes were identified-- create one.
1073 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1074 if (cluster == (Cluster *) NULL)
1076 (void) ThrowMagickException(exception,GetMagickModule(),
1077 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1078 return(MagickFalse);
1081 Initialize a new class.
1085 cluster->green=green;
1087 cluster->next=(Cluster *) NULL;
1091 Count the pixels for each cluster.
1094 for (y=0; y < (ssize_t) image->rows; y++)
1096 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1097 if (p == (const Quantum *) NULL)
1099 for (x=0; x < (ssize_t) image->columns; x++)
1101 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1102 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >=
1103 (cluster->red.left-SafeMargin)) &&
1104 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <=
1105 (cluster->red.right+SafeMargin)) &&
1106 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >=
1107 (cluster->green.left-SafeMargin)) &&
1108 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <=
1109 (cluster->green.right+SafeMargin)) &&
1110 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >=
1111 (cluster->blue.left-SafeMargin)) &&
1112 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <=
1113 (cluster->blue.right+SafeMargin)))
1119 cluster->red.center+=(MagickRealType) ScaleQuantumToChar(
1120 GetPixelRed(image,p));
1121 cluster->green.center+=(MagickRealType) ScaleQuantumToChar(
1122 GetPixelGreen(image,p));
1123 cluster->blue.center+=(MagickRealType) ScaleQuantumToChar(
1124 GetPixelBlue(image,p));
1128 p+=GetPixelChannels(image);
1130 proceed=SetImageProgress(image,SegmentImageTag,(MagickOffsetType) y,
1132 if (proceed == MagickFalse)
1136 Remove clusters that do not meet minimum cluster threshold.
1141 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1143 next_cluster=cluster->next;
1144 if ((cluster->count > 0) &&
1145 (cluster->count >= (count*cluster_threshold/100.0)))
1151 cluster->red.center/=cluster->count;
1152 cluster->green.center/=cluster->count;
1153 cluster->blue.center/=cluster->count;
1155 last_cluster=cluster;
1161 if (cluster == head)
1164 last_cluster->next=next_cluster;
1165 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1172 for (cluster=object; cluster->next != (Cluster *) NULL; )
1174 if (cluster->count < object->count)
1176 cluster=cluster->next;
1178 background=head->next;
1179 for (cluster=background; cluster->next != (Cluster *) NULL; )
1181 if (cluster->count > background->count)
1183 cluster=cluster->next;
1186 threshold=(background->red.center+object->red.center)/2.0;
1187 pixel->red=(MagickRealType) ScaleCharToQuantum((unsigned char)
1189 threshold=(background->green.center+object->green.center)/2.0;
1190 pixel->green=(MagickRealType) ScaleCharToQuantum((unsigned char)
1192 threshold=(background->blue.center+object->blue.center)/2.0;
1193 pixel->blue=(MagickRealType) ScaleCharToQuantum((unsigned char)
1196 Relinquish resources.
1198 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1200 next_cluster=cluster->next;
1201 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1203 for (i=0; i < MaxDimension; i++)
1205 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1206 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1216 + I n i t i a l i z e H i s t o g r a m %
1220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1222 % InitializeHistogram() computes the histogram for an image.
1224 % The format of the InitializeHistogram method is:
1226 % InitializeHistogram(const Image *image,ssize_t **histogram)
1228 % A description of each parameter follows.
1230 % o image: Specifies a pointer to an Image structure; returned from
1233 % o histogram: Specifies an array of integers representing the number
1234 % of pixels for each intensity of a particular color component.
1237 static void InitializeHistogram(const Image *image,ssize_t **histogram,
1238 ExceptionInfo *exception)
1240 register const Quantum
1251 Initialize histogram.
1253 for (i=0; i <= 255; i++)
1255 histogram[Red][i]=0;
1256 histogram[Green][i]=0;
1257 histogram[Blue][i]=0;
1259 for (y=0; y < (ssize_t) image->rows; y++)
1261 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1262 if (p == (const Quantum *) NULL)
1264 for (x=0; x < (ssize_t) image->columns; x++)
1266 histogram[Red][(ssize_t) ScaleQuantumToChar(GetPixelRed(image,p))]++;
1267 histogram[Green][(ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p))]++;
1268 histogram[Blue][(ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p))]++;
1269 p+=GetPixelChannels(image);
1275 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1279 + I n i t i a l i z e I n t e r v a l T r e e %
1283 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1285 % InitializeIntervalTree() initializes an interval tree from the lists of
1288 % The format of the InitializeIntervalTree method is:
1290 % InitializeIntervalTree(IntervalTree **list,ssize_t *number_nodes,
1291 % IntervalTree *node)
1293 % A description of each parameter follows.
1295 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
1297 % o number_crossings: This size_t specifies the number of elements
1298 % in the zero_crossing array.
1302 static void InitializeList(IntervalTree **list,ssize_t *number_nodes,
1305 if (node == (IntervalTree *) NULL)
1307 if (node->child == (IntervalTree *) NULL)
1308 list[(*number_nodes)++]=node;
1309 InitializeList(list,number_nodes,node->sibling);
1310 InitializeList(list,number_nodes,node->child);
1313 static void MeanStability(IntervalTree *node)
1315 register IntervalTree
1318 if (node == (IntervalTree *) NULL)
1320 node->mean_stability=0.0;
1322 if (child != (IntervalTree *) NULL)
1327 register MagickRealType
1332 for ( ; child != (IntervalTree *) NULL; child=child->sibling)
1334 sum+=child->stability;
1337 node->mean_stability=sum/(MagickRealType) count;
1339 MeanStability(node->sibling);
1340 MeanStability(node->child);
1343 static void Stability(IntervalTree *node)
1345 if (node == (IntervalTree *) NULL)
1347 if (node->child == (IntervalTree *) NULL)
1348 node->stability=0.0;
1350 node->stability=node->tau-(node->child)->tau;
1351 Stability(node->sibling);
1352 Stability(node->child);
1355 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing,
1356 const size_t number_crossings)
1374 Allocate interval tree.
1376 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1378 if (list == (IntervalTree **) NULL)
1379 return((IntervalTree *) NULL);
1381 The root is the entire histogram.
1383 root=(IntervalTree *) AcquireMagickMemory(sizeof(*root));
1384 root->child=(IntervalTree *) NULL;
1385 root->sibling=(IntervalTree *) NULL;
1389 for (i=(-1); i < (ssize_t) number_crossings; i++)
1392 Initialize list with all nodes with no children.
1395 InitializeList(list,&number_nodes,root);
1399 for (j=0; j < number_nodes; j++)
1404 for (k=head->left+1; k < head->right; k++)
1406 if (zero_crossing[i+1].crossings[k] != 0)
1410 node->child=(IntervalTree *) AcquireMagickMemory(
1411 sizeof(*node->child));
1416 node->sibling=(IntervalTree *) AcquireMagickMemory(
1417 sizeof(*node->sibling));
1420 node->tau=zero_crossing[i+1].tau;
1421 node->child=(IntervalTree *) NULL;
1422 node->sibling=(IntervalTree *) NULL;
1428 if (left != head->left)
1430 node->sibling=(IntervalTree *) AcquireMagickMemory(
1431 sizeof(*node->sibling));
1433 node->tau=zero_crossing[i+1].tau;
1434 node->child=(IntervalTree *) NULL;
1435 node->sibling=(IntervalTree *) NULL;
1437 node->right=head->right;
1442 Determine the stability: difference between a nodes tau and its child.
1444 Stability(root->child);
1445 MeanStability(root->child);
1446 list=(IntervalTree **) RelinquishMagickMemory(list);
1451 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1455 + O p t i m a l T a u %
1459 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1461 % OptimalTau() finds the optimal tau for each band of the histogram.
1463 % The format of the OptimalTau method is:
1465 % MagickRealType OptimalTau(const ssize_t *histogram,const double max_tau,
1466 % const double min_tau,const double delta_tau,
1467 % const double smooth_threshold,short *extrema)
1469 % A description of each parameter follows.
1471 % o histogram: Specifies an array of integers representing the number
1472 % of pixels for each intensity of a particular color component.
1474 % o extrema: Specifies a pointer to an array of integers. They
1475 % represent the peaks and valleys of the histogram for each color
1480 static void ActiveNodes(IntervalTree **list,ssize_t *number_nodes,
1483 if (node == (IntervalTree *) NULL)
1485 if (node->stability >= node->mean_stability)
1487 list[(*number_nodes)++]=node;
1488 ActiveNodes(list,number_nodes,node->sibling);
1492 ActiveNodes(list,number_nodes,node->sibling);
1493 ActiveNodes(list,number_nodes,node->child);
1497 static void FreeNodes(IntervalTree *node)
1499 if (node == (IntervalTree *) NULL)
1501 FreeNodes(node->sibling);
1502 FreeNodes(node->child);
1503 node=(IntervalTree *) RelinquishMagickMemory(node);
1506 static MagickRealType OptimalTau(const ssize_t *histogram,const double max_tau,
1507 const double min_tau,const double delta_tau,const double smooth_threshold,
1543 Allocate interval tree.
1545 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1547 if (list == (IntervalTree **) NULL)
1550 Allocate zero crossing list.
1552 count=(size_t) ((max_tau-min_tau)/delta_tau)+2;
1553 zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count,
1554 sizeof(*zero_crossing));
1555 if (zero_crossing == (ZeroCrossing *) NULL)
1557 for (i=0; i < (ssize_t) count; i++)
1558 zero_crossing[i].tau=(-1.0);
1560 Initialize zero crossing list.
1562 derivative=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*derivative));
1563 second_derivative=(MagickRealType *) AcquireQuantumMemory(256,
1564 sizeof(*second_derivative));
1565 if ((derivative == (MagickRealType *) NULL) ||
1566 (second_derivative == (MagickRealType *) NULL))
1567 ThrowFatalException(ResourceLimitFatalError,
1568 "UnableToAllocateDerivatives");
1570 for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
1572 zero_crossing[i].tau=tau;
1573 ScaleSpace(histogram,tau,zero_crossing[i].histogram);
1574 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1575 DerivativeHistogram(derivative,second_derivative);
1576 ZeroCrossHistogram(second_derivative,smooth_threshold,
1577 zero_crossing[i].crossings);
1581 Add an entry for the original histogram.
1583 zero_crossing[i].tau=0.0;
1584 for (j=0; j <= 255; j++)
1585 zero_crossing[i].histogram[j]=(MagickRealType) histogram[j];
1586 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1587 DerivativeHistogram(derivative,second_derivative);
1588 ZeroCrossHistogram(second_derivative,smooth_threshold,
1589 zero_crossing[i].crossings);
1590 number_crossings=(size_t) i;
1591 derivative=(MagickRealType *) RelinquishMagickMemory(derivative);
1592 second_derivative=(MagickRealType *)
1593 RelinquishMagickMemory(second_derivative);
1595 Ensure the scale-space fingerprints form lines in scale-space, not loops.
1597 ConsolidateCrossings(zero_crossing,number_crossings);
1599 Force endpoints to be included in the interval.
1601 for (i=0; i <= (ssize_t) number_crossings; i++)
1603 for (j=0; j < 255; j++)
1604 if (zero_crossing[i].crossings[j] != 0)
1606 zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]);
1607 for (j=255; j > 0; j--)
1608 if (zero_crossing[i].crossings[j] != 0)
1610 zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1613 Initialize interval tree.
1615 root=InitializeIntervalTree(zero_crossing,number_crossings);
1616 if (root == (IntervalTree *) NULL)
1619 Find active nodes: stability is greater (or equal) to the mean stability of
1623 ActiveNodes(list,&number_nodes,root->child);
1627 for (i=0; i <= 255; i++)
1629 for (i=0; i < number_nodes; i++)
1632 Find this tau in zero crossings list.
1636 for (j=0; j <= (ssize_t) number_crossings; j++)
1637 if (zero_crossing[j].tau == node->tau)
1640 Find the value of the peak.
1642 peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1645 value=zero_crossing[k].histogram[index];
1646 for (x=node->left; x <= node->right; x++)
1648 if (peak != MagickFalse)
1650 if (zero_crossing[k].histogram[x] > value)
1652 value=zero_crossing[k].histogram[x];
1657 if (zero_crossing[k].histogram[x] < value)
1659 value=zero_crossing[k].histogram[x];
1663 for (x=node->left; x <= node->right; x++)
1667 if (peak != MagickFalse)
1668 extrema[x]=(short) index;
1670 extrema[x]=(short) (-index);
1674 Determine the average tau.
1677 for (i=0; i < number_nodes; i++)
1678 average_tau+=list[i]->tau;
1679 average_tau/=(MagickRealType) number_nodes;
1681 Relinquish resources.
1684 zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1685 list=(IntervalTree **) RelinquishMagickMemory(list);
1686 return(average_tau);
1690 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1694 + S c a l e S p a c e %
1698 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1700 % ScaleSpace() performs a scale-space filter on the 1D histogram.
1702 % The format of the ScaleSpace method is:
1704 % ScaleSpace(const ssize_t *histogram,const MagickRealType tau,
1705 % MagickRealType *scale_histogram)
1707 % A description of each parameter follows.
1709 % o histogram: Specifies an array of MagickRealTypes representing the number
1710 % of pixels for each intensity of a particular color component.
1714 static void ScaleSpace(const ssize_t *histogram,const MagickRealType tau,
1715 MagickRealType *scale_histogram)
1727 gamma=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*gamma));
1728 if (gamma == (MagickRealType *) NULL)
1729 ThrowFatalException(ResourceLimitFatalError,
1730 "UnableToAllocateGammaMap");
1731 alpha=1.0/(tau*sqrt(2.0*MagickPI));
1732 beta=(-1.0/(2.0*tau*tau));
1733 for (x=0; x <= 255; x++)
1735 for (x=0; x <= 255; x++)
1737 gamma[x]=exp((double) beta*x*x);
1738 if (gamma[x] < MagickEpsilon)
1741 for (x=0; x <= 255; x++)
1744 for (u=0; u <= 255; u++)
1745 sum+=(MagickRealType) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1746 scale_histogram[x]=alpha*sum;
1748 gamma=(MagickRealType *) RelinquishMagickMemory(gamma);
1752 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1756 % S e g m e n t I m a g e %
1760 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1762 % SegmentImage() segment an image by analyzing the histograms of the color
1763 % components and identifying units that are homogeneous with the fuzzy
1764 % C-means technique.
1766 % The format of the SegmentImage method is:
1768 % MagickBooleanType SegmentImage(Image *image,
1769 % const ColorspaceType colorspace,const MagickBooleanType verbose,
1770 % const double cluster_threshold,const double smooth_threshold,
1771 % ExceptionInfo *exception)
1773 % A description of each parameter follows.
1775 % o image: the image.
1777 % o colorspace: Indicate the colorspace.
1779 % o verbose: Set to MagickTrue to print detailed information about the
1780 % identified classes.
1782 % o cluster_threshold: This represents the minimum number of pixels
1783 % contained in a hexahedra before it can be considered valid (expressed
1786 % o smooth_threshold: the smoothing threshold eliminates noise in the second
1787 % derivative of the histogram. As the value is increased, you can expect a
1788 % smoother second derivative.
1790 % o exception: return any errors or warnings in this structure.
1793 MagickExport MagickBooleanType SegmentImage(Image *image,
1794 const ColorspaceType colorspace,const MagickBooleanType verbose,
1795 const double cluster_threshold,const double smooth_threshold,
1796 ExceptionInfo *exception)
1805 *extrema[MaxDimension];
1808 *histogram[MaxDimension];
1811 Allocate histogram and extrema.
1813 assert(image != (Image *) NULL);
1814 assert(image->signature == MagickSignature);
1815 if (image->debug != MagickFalse)
1816 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1817 for (i=0; i < MaxDimension; i++)
1819 histogram[i]=(ssize_t *) AcquireQuantumMemory(256,sizeof(**histogram));
1820 extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1821 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1823 for (i-- ; i >= 0; i--)
1825 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1826 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1828 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1832 if (IsRGBColorspace(colorspace) == MagickFalse)
1833 (void) TransformImageColorspace(image,colorspace);
1835 Initialize histogram.
1837 InitializeHistogram(image,histogram,&image->exception);
1838 (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,
1839 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Red]);
1840 (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,
1841 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Green]);
1842 (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,
1843 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Blue]);
1845 Classify using the fuzzy c-Means technique.
1847 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose,
1849 if (IsRGBColorspace(colorspace) == MagickFalse)
1850 (void) TransformImageColorspace(image,colorspace);
1852 Relinquish resources.
1854 for (i=0; i < MaxDimension; i++)
1856 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1857 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1863 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1867 + Z e r o C r o s s H i s t o g r a m %
1871 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1873 % ZeroCrossHistogram() find the zero crossings in a histogram and marks
1874 % directions as: 1 is negative to positive; 0 is zero crossing; and -1
1875 % is positive to negative.
1877 % The format of the ZeroCrossHistogram method is:
1879 % ZeroCrossHistogram(MagickRealType *second_derivative,
1880 % const MagickRealType smooth_threshold,short *crossings)
1882 % A description of each parameter follows.
1884 % o second_derivative: Specifies an array of MagickRealTypes representing the
1885 % second derivative of the histogram of a particular color component.
1887 % o crossings: This array of integers is initialized with
1888 % -1, 0, or 1 representing the slope of the first derivative of the
1889 % of a particular color component.
1892 static void ZeroCrossHistogram(MagickRealType *second_derivative,
1893 const MagickRealType smooth_threshold,short *crossings)
1902 Merge low numbers to zero to help prevent noise.
1904 for (i=0; i <= 255; i++)
1905 if ((second_derivative[i] < smooth_threshold) &&
1906 (second_derivative[i] >= -smooth_threshold))
1907 second_derivative[i]=0.0;
1909 Mark zero crossings.
1912 for (i=0; i <= 255; i++)
1915 if (second_derivative[i] < 0.0)
1922 if (second_derivative[i] > 0.0)