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-2012 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/pixel-private.h"
100 #include "MagickCore/quantize.h"
101 #include "MagickCore/quantum.h"
102 #include "MagickCore/quantum-private.h"
103 #include "MagickCore/resource_.h"
104 #include "MagickCore/segment.h"
105 #include "MagickCore/string_.h"
106 #include "MagickCore/thread-private.h"
111 #define MaxDimension 3
112 #define DeltaTau 0.5f
113 #if defined(FastClassify)
114 #define WeightingExponent 2.0
115 #define SegmentPower(ratio) (ratio)
117 #define WeightingExponent 2.5
118 #define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0)));
123 Typedef declarations.
125 typedef struct _ExtentPacket
136 typedef struct _Cluster
151 typedef struct _IntervalTree
169 typedef struct _ZeroCrossing
180 Constant declarations.
193 OptimalTau(const ssize_t *,const double,const double,const double,
194 const double,short *);
197 DefineRegion(const short *,ExtentPacket *);
200 InitializeHistogram(const Image *,ssize_t **,ExceptionInfo *),
201 ScaleSpace(const ssize_t *,const double,double *),
202 ZeroCrossHistogram(double *,const double,short *);
205 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
219 % The format of the Classify method is:
221 % MagickBooleanType Classify(Image *image,short **extrema,
222 % const double cluster_threshold,
223 % const double weighting_exponent,
224 % const MagickBooleanType verbose,ExceptionInfo *exception)
226 % A description of each parameter follows.
228 % o image: the image.
230 % o extrema: Specifies a pointer to an array of integers. They
231 % represent the peaks and valleys of the histogram for each color
234 % o cluster_threshold: This double represents the minimum number of
235 % pixels contained in a hexahedra before it can be considered valid
236 % (expressed as a percentage).
238 % o weighting_exponent: Specifies the membership weighting exponent.
240 % o verbose: A value greater than zero prints detailed information about
241 % the identified classes.
243 % o exception: return any errors or warnings in this structure.
246 static MagickBooleanType Classify(Image *image,short **extrema,
247 const double cluster_threshold,
248 const double weighting_exponent,const MagickBooleanType verbose,
249 ExceptionInfo *exception)
251 #define SegmentImageTag "Segment/Image"
292 cluster=(Cluster *) NULL;
293 head=(Cluster *) NULL;
294 (void) ResetMagickMemory(&red,0,sizeof(red));
295 (void) ResetMagickMemory(&green,0,sizeof(green));
296 (void) ResetMagickMemory(&blue,0,sizeof(blue));
297 while (DefineRegion(extrema[Red],&red) != 0)
300 while (DefineRegion(extrema[Green],&green) != 0)
303 while (DefineRegion(extrema[Blue],&blue) != 0)
306 Allocate a new class.
308 if (head != (Cluster *) NULL)
310 cluster->next=(Cluster *) AcquireMagickMemory(
311 sizeof(*cluster->next));
312 cluster=cluster->next;
316 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
319 if (cluster == (Cluster *) NULL)
320 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
323 Initialize a new class.
327 cluster->green=green;
329 cluster->next=(Cluster *) NULL;
333 if (head == (Cluster *) NULL)
336 No classes were identified-- create one.
338 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
339 if (cluster == (Cluster *) NULL)
340 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
343 Initialize a new class.
347 cluster->green=green;
349 cluster->next=(Cluster *) NULL;
353 Count the pixels for each cluster.
358 image_view=AcquireVirtualCacheView(image,exception);
359 for (y=0; y < (ssize_t) image->rows; y++)
361 register const Quantum
367 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
368 if (p == (const Quantum *) NULL)
370 for (x=0; x < (ssize_t) image->columns; x++)
372 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
373 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >=
374 (cluster->red.left-SafeMargin)) &&
375 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <=
376 (cluster->red.right+SafeMargin)) &&
377 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >=
378 (cluster->green.left-SafeMargin)) &&
379 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <=
380 (cluster->green.right+SafeMargin)) &&
381 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >=
382 (cluster->blue.left-SafeMargin)) &&
383 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <=
384 (cluster->blue.right+SafeMargin)))
390 cluster->red.center+=(double) ScaleQuantumToChar(
391 GetPixelRed(image,p));
392 cluster->green.center+=(double) ScaleQuantumToChar(
393 GetPixelGreen(image,p));
394 cluster->blue.center+=(double) ScaleQuantumToChar(
395 GetPixelBlue(image,p));
399 p+=GetPixelChannels(image);
401 if (image->progress_monitor != (MagickProgressMonitor) NULL)
406 #if defined(MAGICKCORE_OPENMP_SUPPORT)
407 #pragma omp critical (MagickCore_Classify)
409 proceed=SetImageProgress(image,SegmentImageTag,progress++,2*
411 if (proceed == MagickFalse)
415 image_view=DestroyCacheView(image_view);
417 Remove clusters that do not meet minimum cluster threshold.
422 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
424 next_cluster=cluster->next;
425 if ((cluster->count > 0) &&
426 (cluster->count >= (count*cluster_threshold/100.0)))
432 cluster->red.center/=cluster->count;
433 cluster->green.center/=cluster->count;
434 cluster->blue.center/=cluster->count;
436 last_cluster=cluster;
445 last_cluster->next=next_cluster;
446 cluster=(Cluster *) RelinquishMagickMemory(cluster);
448 number_clusters=(size_t) count;
449 if (verbose != MagickFalse)
452 Print cluster statistics.
454 (void) FormatLocaleFile(stdout,"Fuzzy C-means Statistics\n");
455 (void) FormatLocaleFile(stdout,"===================\n\n");
456 (void) FormatLocaleFile(stdout,"\tCluster Threshold = %g\n",(double)
458 (void) FormatLocaleFile(stdout,"\tWeighting Exponent = %g\n",(double)
460 (void) FormatLocaleFile(stdout,"\tTotal Number of Clusters = %.20g\n\n",
461 (double) number_clusters);
463 Print the total number of points per cluster.
465 (void) FormatLocaleFile(stdout,"\n\nNumber of Vectors Per Cluster\n");
466 (void) FormatLocaleFile(stdout,"=============================\n\n");
467 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
468 (void) FormatLocaleFile(stdout,"Cluster #%.20g = %.20g\n",(double)
469 cluster->id,(double) cluster->count);
471 Print the cluster extents.
473 (void) FormatLocaleFile(stdout,
474 "\n\n\nCluster Extents: (Vector Size: %d)\n",MaxDimension);
475 (void) FormatLocaleFile(stdout,"================");
476 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
478 (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
480 (void) FormatLocaleFile(stdout,
481 "%.20g-%.20g %.20g-%.20g %.20g-%.20g\n",(double)
482 cluster->red.left,(double) cluster->red.right,(double)
483 cluster->green.left,(double) cluster->green.right,(double)
484 cluster->blue.left,(double) cluster->blue.right);
487 Print the cluster center values.
489 (void) FormatLocaleFile(stdout,
490 "\n\n\nCluster Center Values: (Vector Size: %d)\n",MaxDimension);
491 (void) FormatLocaleFile(stdout,"=====================");
492 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
494 (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)
496 (void) FormatLocaleFile(stdout,"%g %g %g\n",(double)
497 cluster->red.center,(double) cluster->green.center,(double)
498 cluster->blue.center);
500 (void) FormatLocaleFile(stdout,"\n");
502 if (number_clusters > 256)
503 ThrowBinaryException(ImageError,"TooManyClusters",image->filename);
505 Speed up distance calculations.
507 squares=(double *) AcquireQuantumMemory(513UL,sizeof(*squares));
508 if (squares == (double *) NULL)
509 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
512 for (i=(-255); i <= 255; i++)
513 squares[i]=(double) i*(double) i;
515 Allocate image colormap.
517 if (AcquireImageColormap(image,number_clusters,exception) == MagickFalse)
518 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
521 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
523 image->colormap[i].red=(double) ScaleCharToQuantum((unsigned char)
524 (cluster->red.center+0.5));
525 image->colormap[i].green=(double) ScaleCharToQuantum((unsigned char)
526 (cluster->green.center+0.5));
527 image->colormap[i].blue=(double) ScaleCharToQuantum((unsigned char)
528 (cluster->blue.center+0.5));
532 Do course grain classes.
534 image_view=AcquireAuthenticCacheView(image,exception);
535 #if defined(MAGICKCORE_OPENMP_SUPPORT)
536 #pragma omp parallel for schedule(static,4) shared(progress,status) \
537 dynamic_number_threads(image,image->columns,image->rows,1)
539 for (y=0; y < (ssize_t) image->rows; y++)
544 register const PixelInfo
553 if (status == MagickFalse)
555 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
556 if (q == (Quantum *) NULL)
561 for (x=0; x < (ssize_t) image->columns; x++)
563 SetPixelIndex(image,0,q);
564 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
566 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,q)) >=
567 (cluster->red.left-SafeMargin)) &&
568 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,q)) <=
569 (cluster->red.right+SafeMargin)) &&
570 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q)) >=
571 (cluster->green.left-SafeMargin)) &&
572 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q)) <=
573 (cluster->green.right+SafeMargin)) &&
574 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q)) >=
575 (cluster->blue.left-SafeMargin)) &&
576 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q)) <=
577 (cluster->blue.right+SafeMargin)))
582 SetPixelIndex(image,(Quantum) cluster->id,q);
586 if (cluster == (Cluster *) NULL)
600 Compute fuzzy membership.
603 for (j=0; j < (ssize_t) image->colors; j++)
607 distance_squared=squares[(ssize_t) ScaleQuantumToChar(
608 GetPixelRed(image,q))-(ssize_t)
609 ScaleQuantumToChar(ClampToQuantum(p->red))]+squares[(ssize_t)
610 ScaleQuantumToChar(GetPixelGreen(image,q))-(ssize_t)
611 ScaleQuantumToChar(ClampToQuantum(p->green))]+squares[(ssize_t)
612 ScaleQuantumToChar(GetPixelBlue(image,q))-(ssize_t)
613 ScaleQuantumToChar(ClampToQuantum(p->blue))];
614 numerator=distance_squared;
615 for (k=0; k < (ssize_t) image->colors; k++)
618 distance_squared=squares[(ssize_t) ScaleQuantumToChar(
619 GetPixelRed(image,q))-(ssize_t)
620 ScaleQuantumToChar(ClampToQuantum(p->red))]+squares[
621 (ssize_t) ScaleQuantumToChar(GetPixelGreen(image,q))-(ssize_t)
622 ScaleQuantumToChar(ClampToQuantum(p->green))]+squares[
623 (ssize_t) ScaleQuantumToChar(GetPixelBlue(image,q))-(ssize_t)
624 ScaleQuantumToChar(ClampToQuantum(p->blue))];
625 ratio=numerator/distance_squared;
626 sum+=SegmentPower(ratio);
628 if ((sum != 0.0) && ((1.0/sum) > local_minima))
633 local_minima=1.0/sum;
634 SetPixelIndex(image,(Quantum) j,q);
638 q+=GetPixelChannels(image);
640 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
642 if (image->progress_monitor != (MagickProgressMonitor) NULL)
647 #if defined(MAGICKCORE_OPENMP_SUPPORT)
648 #pragma omp critical (MagickCore_Classify)
650 proceed=SetImageProgress(image,SegmentImageTag,progress++,
652 if (proceed == MagickFalse)
656 image_view=DestroyCacheView(image_view);
657 status&=SyncImage(image,exception);
659 Relinquish resources.
661 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
663 next_cluster=cluster->next;
664 cluster=(Cluster *) RelinquishMagickMemory(cluster);
667 free_squares=squares;
668 free_squares=(double *) RelinquishMagickMemory(free_squares);
673 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
677 + C o n s o l i d a t e C r o s s i n g s %
681 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
683 % ConsolidateCrossings() guarantees that an even number of zero crossings
684 % always lie between two crossings.
686 % The format of the ConsolidateCrossings method is:
688 % ConsolidateCrossings(ZeroCrossing *zero_crossing,
689 % const size_t number_crossings)
691 % A description of each parameter follows.
693 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
695 % o number_crossings: This size_t specifies the number of elements
696 % in the zero_crossing array.
700 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
707 static inline ssize_t MagickMax(const ssize_t x,const ssize_t y)
714 static inline ssize_t MagickMin(const ssize_t x,const ssize_t y)
721 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
722 const size_t number_crossings)
738 Consolidate zero crossings.
740 for (i=(ssize_t) number_crossings-1; i >= 0; i--)
741 for (j=0; j <= 255; j++)
743 if (zero_crossing[i].crossings[j] == 0)
746 Find the entry that is closest to j and still preserves the
747 property that there are an even number of crossings between
750 for (k=j-1; k > 0; k--)
751 if (zero_crossing[i+1].crossings[k] != 0)
755 for (k=j+1; k < 255; k++)
756 if (zero_crossing[i+1].crossings[k] != 0)
758 right=MagickMin(k,255);
760 K is the zero crossing just left of j.
762 for (k=j-1; k > 0; k--)
763 if (zero_crossing[i].crossings[k] != 0)
768 Check center for an even number of crossings between k and j.
771 if (zero_crossing[i+1].crossings[j] != 0)
774 for (l=k+1; l < center; l++)
775 if (zero_crossing[i+1].crossings[l] != 0)
777 if (((count % 2) == 0) && (center != k))
781 Check left for an even number of crossings between k and j.
786 for (l=k+1; l < left; l++)
787 if (zero_crossing[i+1].crossings[l] != 0)
789 if (((count % 2) == 0) && (left != k))
793 Check right for an even number of crossings between k and j.
798 for (l=k+1; l < right; l++)
799 if (zero_crossing[i+1].crossings[l] != 0)
801 if (((count % 2) == 0) && (right != k))
804 l=(ssize_t) zero_crossing[i].crossings[j];
805 zero_crossing[i].crossings[j]=0;
807 zero_crossing[i].crossings[correct]=(short) l;
812 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
816 + D e f i n e R e g i o n %
820 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
822 % DefineRegion() defines the left and right boundaries of a peak region.
824 % The format of the DefineRegion method is:
826 % ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
828 % A description of each parameter follows.
830 % o extrema: Specifies a pointer to an array of integers. They
831 % represent the peaks and valleys of the histogram for each color
834 % o extents: This pointer to an ExtentPacket represent the extends
835 % of a particular peak or valley of a color component.
838 static ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
841 Initialize to default values.
847 Find the left side (maxima).
849 for ( ; extents->index <= 255; extents->index++)
850 if (extrema[extents->index] > 0)
852 if (extents->index > 255)
853 return(MagickFalse); /* no left side - no region exists */
854 extents->left=extents->index;
856 Find the right side (minima).
858 for ( ; extents->index <= 255; extents->index++)
859 if (extrema[extents->index] < 0)
861 extents->right=extents->index-1;
866 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
870 + D e r i v a t i v e H i s t o g r a m %
874 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
876 % DerivativeHistogram() determines the derivative of the histogram using
877 % central differencing.
879 % The format of the DerivativeHistogram method is:
881 % DerivativeHistogram(const double *histogram,
882 % double *derivative)
884 % A description of each parameter follows.
886 % o histogram: Specifies an array of doubles representing the number
887 % of pixels for each intensity of a particular color component.
889 % o derivative: This array of doubles is initialized by
890 % DerivativeHistogram to the derivative of the histogram using central
894 static void DerivativeHistogram(const double *histogram,
902 Compute endpoints using second order polynomial interpolation.
905 derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);
906 derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);
908 Compute derivative using central differencing.
910 for (i=1; i < n; i++)
911 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
916 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
920 + G e t I m a g e D y n a m i c T h r e s h o l d %
924 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
926 % GetImageDynamicThreshold() returns the dynamic threshold for an image.
928 % The format of the GetImageDynamicThreshold method is:
930 % MagickBooleanType GetImageDynamicThreshold(const Image *image,
931 % const double cluster_threshold,const double smooth_threshold,
932 % PixelInfo *pixel,ExceptionInfo *exception)
934 % A description of each parameter follows.
936 % o image: the image.
938 % o cluster_threshold: This double represents the minimum number of
939 % pixels contained in a hexahedra before it can be considered valid
940 % (expressed as a percentage).
942 % o smooth_threshold: the smoothing threshold eliminates noise in the second
943 % derivative of the histogram. As the value is increased, you can expect a
944 % smoother second derivative.
946 % o pixel: return the dynamic threshold here.
948 % o exception: return any errors or warnings in this structure.
951 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,
952 const double cluster_threshold,const double smooth_threshold,
953 PixelInfo *pixel,ExceptionInfo *exception)
974 register const Quantum
982 *extrema[MaxDimension];
986 *histogram[MaxDimension],
990 Allocate histogram and extrema.
992 assert(image != (Image *) NULL);
993 assert(image->signature == MagickSignature);
994 if (image->debug != MagickFalse)
995 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
996 GetPixelInfo(image,pixel);
997 for (i=0; i < MaxDimension; i++)
999 histogram[i]=(ssize_t *) AcquireQuantumMemory(256UL,sizeof(**histogram));
1000 extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
1001 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1003 for (i-- ; i >= 0; i--)
1005 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1006 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1008 (void) ThrowMagickException(exception,GetMagickModule(),
1009 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
1010 return(MagickFalse);
1014 Initialize histogram.
1016 InitializeHistogram(image,histogram,exception);
1017 (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
1018 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
1019 (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
1020 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
1021 (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
1022 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
1026 cluster=(Cluster *) NULL;
1027 head=(Cluster *) NULL;
1028 (void) ResetMagickMemory(&red,0,sizeof(red));
1029 (void) ResetMagickMemory(&green,0,sizeof(green));
1030 (void) ResetMagickMemory(&blue,0,sizeof(blue));
1031 while (DefineRegion(extrema[Red],&red) != 0)
1034 while (DefineRegion(extrema[Green],&green) != 0)
1037 while (DefineRegion(extrema[Blue],&blue) != 0)
1040 Allocate a new class.
1042 if (head != (Cluster *) NULL)
1044 cluster->next=(Cluster *) AcquireMagickMemory(
1045 sizeof(*cluster->next));
1046 cluster=cluster->next;
1050 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1053 if (cluster == (Cluster *) NULL)
1055 (void) ThrowMagickException(exception,GetMagickModule(),
1056 ResourceLimitError,"MemoryAllocationFailed","'%s'",
1058 return(MagickFalse);
1061 Initialize a new class.
1065 cluster->green=green;
1067 cluster->next=(Cluster *) NULL;
1071 if (head == (Cluster *) NULL)
1074 No classes were identified-- create one.
1076 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1077 if (cluster == (Cluster *) NULL)
1079 (void) ThrowMagickException(exception,GetMagickModule(),
1080 ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
1081 return(MagickFalse);
1084 Initialize a new class.
1088 cluster->green=green;
1090 cluster->next=(Cluster *) NULL;
1094 Count the pixels for each cluster.
1097 for (y=0; y < (ssize_t) image->rows; y++)
1099 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1100 if (p == (const Quantum *) NULL)
1102 for (x=0; x < (ssize_t) image->columns; x++)
1104 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1105 if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >=
1106 (cluster->red.left-SafeMargin)) &&
1107 ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <=
1108 (cluster->red.right+SafeMargin)) &&
1109 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >=
1110 (cluster->green.left-SafeMargin)) &&
1111 ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <=
1112 (cluster->green.right+SafeMargin)) &&
1113 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >=
1114 (cluster->blue.left-SafeMargin)) &&
1115 ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <=
1116 (cluster->blue.right+SafeMargin)))
1122 cluster->red.center+=(double) ScaleQuantumToChar(
1123 GetPixelRed(image,p));
1124 cluster->green.center+=(double) ScaleQuantumToChar(
1125 GetPixelGreen(image,p));
1126 cluster->blue.center+=(double) ScaleQuantumToChar(
1127 GetPixelBlue(image,p));
1131 p+=GetPixelChannels(image);
1133 proceed=SetImageProgress(image,SegmentImageTag,(MagickOffsetType) y,
1135 if (proceed == MagickFalse)
1139 Remove clusters that do not meet minimum cluster threshold.
1144 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1146 next_cluster=cluster->next;
1147 if ((cluster->count > 0) &&
1148 (cluster->count >= (count*cluster_threshold/100.0)))
1154 cluster->red.center/=cluster->count;
1155 cluster->green.center/=cluster->count;
1156 cluster->blue.center/=cluster->count;
1158 last_cluster=cluster;
1164 if (cluster == head)
1167 last_cluster->next=next_cluster;
1168 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1175 for (cluster=object; cluster->next != (Cluster *) NULL; )
1177 if (cluster->count < object->count)
1179 cluster=cluster->next;
1181 background=head->next;
1182 for (cluster=background; cluster->next != (Cluster *) NULL; )
1184 if (cluster->count > background->count)
1186 cluster=cluster->next;
1189 threshold=(background->red.center+object->red.center)/2.0;
1190 pixel->red=(double) ScaleCharToQuantum((unsigned char)
1192 threshold=(background->green.center+object->green.center)/2.0;
1193 pixel->green=(double) ScaleCharToQuantum((unsigned char)
1195 threshold=(background->blue.center+object->blue.center)/2.0;
1196 pixel->blue=(double) ScaleCharToQuantum((unsigned char)
1199 Relinquish resources.
1201 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1203 next_cluster=cluster->next;
1204 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1206 for (i=0; i < MaxDimension; i++)
1208 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1209 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1215 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1219 + I n i t i a l i z e H i s t o g r a m %
1223 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1225 % InitializeHistogram() computes the histogram for an image.
1227 % The format of the InitializeHistogram method is:
1229 % InitializeHistogram(const Image *image,ssize_t **histogram)
1231 % A description of each parameter follows.
1233 % o image: Specifies a pointer to an Image structure; returned from
1236 % o histogram: Specifies an array of integers representing the number
1237 % of pixels for each intensity of a particular color component.
1240 static void InitializeHistogram(const Image *image,ssize_t **histogram,
1241 ExceptionInfo *exception)
1243 register const Quantum
1254 Initialize histogram.
1256 for (i=0; i <= 255; i++)
1258 histogram[Red][i]=0;
1259 histogram[Green][i]=0;
1260 histogram[Blue][i]=0;
1262 for (y=0; y < (ssize_t) image->rows; y++)
1264 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1265 if (p == (const Quantum *) NULL)
1267 for (x=0; x < (ssize_t) image->columns; x++)
1269 histogram[Red][(ssize_t) ScaleQuantumToChar(GetPixelRed(image,p))]++;
1270 histogram[Green][(ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p))]++;
1271 histogram[Blue][(ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p))]++;
1272 p+=GetPixelChannels(image);
1278 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1282 + I n i t i a l i z e I n t e r v a l T r e e %
1286 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1288 % InitializeIntervalTree() initializes an interval tree from the lists of
1291 % The format of the InitializeIntervalTree method is:
1293 % InitializeIntervalTree(IntervalTree **list,ssize_t *number_nodes,
1294 % IntervalTree *node)
1296 % A description of each parameter follows.
1298 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
1300 % o number_crossings: This size_t specifies the number of elements
1301 % in the zero_crossing array.
1305 static void InitializeList(IntervalTree **list,ssize_t *number_nodes,
1308 if (node == (IntervalTree *) NULL)
1310 if (node->child == (IntervalTree *) NULL)
1311 list[(*number_nodes)++]=node;
1312 InitializeList(list,number_nodes,node->sibling);
1313 InitializeList(list,number_nodes,node->child);
1316 static void MeanStability(IntervalTree *node)
1318 register IntervalTree
1321 if (node == (IntervalTree *) NULL)
1323 node->mean_stability=0.0;
1325 if (child != (IntervalTree *) NULL)
1335 for ( ; child != (IntervalTree *) NULL; child=child->sibling)
1337 sum+=child->stability;
1340 node->mean_stability=sum/(double) count;
1342 MeanStability(node->sibling);
1343 MeanStability(node->child);
1346 static void Stability(IntervalTree *node)
1348 if (node == (IntervalTree *) NULL)
1350 if (node->child == (IntervalTree *) NULL)
1351 node->stability=0.0;
1353 node->stability=node->tau-(node->child)->tau;
1354 Stability(node->sibling);
1355 Stability(node->child);
1358 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing,
1359 const size_t number_crossings)
1377 Allocate interval tree.
1379 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1381 if (list == (IntervalTree **) NULL)
1382 return((IntervalTree *) NULL);
1384 The root is the entire histogram.
1386 root=(IntervalTree *) AcquireMagickMemory(sizeof(*root));
1387 root->child=(IntervalTree *) NULL;
1388 root->sibling=(IntervalTree *) NULL;
1392 for (i=(-1); i < (ssize_t) number_crossings; i++)
1395 Initialize list with all nodes with no children.
1398 InitializeList(list,&number_nodes,root);
1402 for (j=0; j < number_nodes; j++)
1407 for (k=head->left+1; k < head->right; k++)
1409 if (zero_crossing[i+1].crossings[k] != 0)
1413 node->child=(IntervalTree *) AcquireMagickMemory(
1414 sizeof(*node->child));
1419 node->sibling=(IntervalTree *) AcquireMagickMemory(
1420 sizeof(*node->sibling));
1423 node->tau=zero_crossing[i+1].tau;
1424 node->child=(IntervalTree *) NULL;
1425 node->sibling=(IntervalTree *) NULL;
1431 if (left != head->left)
1433 node->sibling=(IntervalTree *) AcquireMagickMemory(
1434 sizeof(*node->sibling));
1436 node->tau=zero_crossing[i+1].tau;
1437 node->child=(IntervalTree *) NULL;
1438 node->sibling=(IntervalTree *) NULL;
1440 node->right=head->right;
1445 Determine the stability: difference between a nodes tau and its child.
1447 Stability(root->child);
1448 MeanStability(root->child);
1449 list=(IntervalTree **) RelinquishMagickMemory(list);
1454 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1458 + O p t i m a l T a u %
1462 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1464 % OptimalTau() finds the optimal tau for each band of the histogram.
1466 % The format of the OptimalTau method is:
1468 % double OptimalTau(const ssize_t *histogram,const double max_tau,
1469 % const double min_tau,const double delta_tau,
1470 % const double smooth_threshold,short *extrema)
1472 % A description of each parameter follows.
1474 % o histogram: Specifies an array of integers representing the number
1475 % of pixels for each intensity of a particular color component.
1477 % o extrema: Specifies a pointer to an array of integers. They
1478 % represent the peaks and valleys of the histogram for each color
1483 static void ActiveNodes(IntervalTree **list,ssize_t *number_nodes,
1486 if (node == (IntervalTree *) NULL)
1488 if (node->stability >= node->mean_stability)
1490 list[(*number_nodes)++]=node;
1491 ActiveNodes(list,number_nodes,node->sibling);
1495 ActiveNodes(list,number_nodes,node->sibling);
1496 ActiveNodes(list,number_nodes,node->child);
1500 static void FreeNodes(IntervalTree *node)
1502 if (node == (IntervalTree *) NULL)
1504 FreeNodes(node->sibling);
1505 FreeNodes(node->child);
1506 node=(IntervalTree *) RelinquishMagickMemory(node);
1509 static double OptimalTau(const ssize_t *histogram,const double max_tau,
1510 const double min_tau,const double delta_tau,const double smooth_threshold,
1546 Allocate interval tree.
1548 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1550 if (list == (IntervalTree **) NULL)
1553 Allocate zero crossing list.
1555 count=(size_t) ((max_tau-min_tau)/delta_tau)+2;
1556 zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count,
1557 sizeof(*zero_crossing));
1558 if (zero_crossing == (ZeroCrossing *) NULL)
1560 for (i=0; i < (ssize_t) count; i++)
1561 zero_crossing[i].tau=(-1.0);
1563 Initialize zero crossing list.
1565 derivative=(double *) AcquireQuantumMemory(256,sizeof(*derivative));
1566 second_derivative=(double *) AcquireQuantumMemory(256,
1567 sizeof(*second_derivative));
1568 if ((derivative == (double *) NULL) ||
1569 (second_derivative == (double *) NULL))
1570 ThrowFatalException(ResourceLimitFatalError,
1571 "UnableToAllocateDerivatives");
1573 for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
1575 zero_crossing[i].tau=tau;
1576 ScaleSpace(histogram,tau,zero_crossing[i].histogram);
1577 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1578 DerivativeHistogram(derivative,second_derivative);
1579 ZeroCrossHistogram(second_derivative,smooth_threshold,
1580 zero_crossing[i].crossings);
1584 Add an entry for the original histogram.
1586 zero_crossing[i].tau=0.0;
1587 for (j=0; j <= 255; j++)
1588 zero_crossing[i].histogram[j]=(double) histogram[j];
1589 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1590 DerivativeHistogram(derivative,second_derivative);
1591 ZeroCrossHistogram(second_derivative,smooth_threshold,
1592 zero_crossing[i].crossings);
1593 number_crossings=(size_t) i;
1594 derivative=(double *) RelinquishMagickMemory(derivative);
1595 second_derivative=(double *)
1596 RelinquishMagickMemory(second_derivative);
1598 Ensure the scale-space fingerprints form lines in scale-space, not loops.
1600 ConsolidateCrossings(zero_crossing,number_crossings);
1602 Force endpoints to be included in the interval.
1604 for (i=0; i <= (ssize_t) number_crossings; i++)
1606 for (j=0; j < 255; j++)
1607 if (zero_crossing[i].crossings[j] != 0)
1609 zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]);
1610 for (j=255; j > 0; j--)
1611 if (zero_crossing[i].crossings[j] != 0)
1613 zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1616 Initialize interval tree.
1618 root=InitializeIntervalTree(zero_crossing,number_crossings);
1619 if (root == (IntervalTree *) NULL)
1622 Find active nodes: stability is greater (or equal) to the mean stability of
1626 ActiveNodes(list,&number_nodes,root->child);
1630 for (i=0; i <= 255; i++)
1632 for (i=0; i < number_nodes; i++)
1635 Find this tau in zero crossings list.
1639 for (j=0; j <= (ssize_t) number_crossings; j++)
1640 if (zero_crossing[j].tau == node->tau)
1643 Find the value of the peak.
1645 peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1648 value=zero_crossing[k].histogram[index];
1649 for (x=node->left; x <= node->right; x++)
1651 if (peak != MagickFalse)
1653 if (zero_crossing[k].histogram[x] > value)
1655 value=zero_crossing[k].histogram[x];
1660 if (zero_crossing[k].histogram[x] < value)
1662 value=zero_crossing[k].histogram[x];
1666 for (x=node->left; x <= node->right; x++)
1670 if (peak != MagickFalse)
1671 extrema[x]=(short) index;
1673 extrema[x]=(short) (-index);
1677 Determine the average tau.
1680 for (i=0; i < number_nodes; i++)
1681 average_tau+=list[i]->tau;
1682 average_tau/=(double) number_nodes;
1684 Relinquish resources.
1687 zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1688 list=(IntervalTree **) RelinquishMagickMemory(list);
1689 return(average_tau);
1693 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1697 + S c a l e S p a c e %
1701 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1703 % ScaleSpace() performs a scale-space filter on the 1D histogram.
1705 % The format of the ScaleSpace method is:
1707 % ScaleSpace(const ssize_t *histogram,const double tau,
1708 % double *scale_histogram)
1710 % A description of each parameter follows.
1712 % o histogram: Specifies an array of doubles representing the number
1713 % of pixels for each intensity of a particular color component.
1717 static void ScaleSpace(const ssize_t *histogram,const double tau,
1718 double *scale_histogram)
1730 gamma=(double *) AcquireQuantumMemory(256,sizeof(*gamma));
1731 if (gamma == (double *) NULL)
1732 ThrowFatalException(ResourceLimitFatalError,
1733 "UnableToAllocateGammaMap");
1734 alpha=MagickEpsilonReciprocal(tau*sqrt(2.0*MagickPI));
1735 beta=(-1.0*MagickEpsilonReciprocal(2.0*tau*tau));
1736 for (x=0; x <= 255; x++)
1738 for (x=0; x <= 255; x++)
1740 gamma[x]=exp((double) beta*x*x);
1741 if (gamma[x] < MagickEpsilon)
1744 for (x=0; x <= 255; x++)
1747 for (u=0; u <= 255; u++)
1748 sum+=(double) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1749 scale_histogram[x]=alpha*sum;
1751 gamma=(double *) RelinquishMagickMemory(gamma);
1755 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1759 % S e g m e n t I m a g e %
1763 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1765 % SegmentImage() segment an image by analyzing the histograms of the color
1766 % components and identifying units that are homogeneous with the fuzzy
1767 % C-means technique.
1769 % The format of the SegmentImage method is:
1771 % MagickBooleanType SegmentImage(Image *image,
1772 % const ColorspaceType colorspace,const MagickBooleanType verbose,
1773 % const double cluster_threshold,const double smooth_threshold,
1774 % ExceptionInfo *exception)
1776 % A description of each parameter follows.
1778 % o image: the image.
1780 % o colorspace: Indicate the colorspace.
1782 % o verbose: Set to MagickTrue to print detailed information about the
1783 % identified classes.
1785 % o cluster_threshold: This represents the minimum number of pixels
1786 % contained in a hexahedra before it can be considered valid (expressed
1789 % o smooth_threshold: the smoothing threshold eliminates noise in the second
1790 % derivative of the histogram. As the value is increased, you can expect a
1791 % smoother second derivative.
1793 % o exception: return any errors or warnings in this structure.
1796 MagickExport MagickBooleanType SegmentImage(Image *image,
1797 const ColorspaceType colorspace,const MagickBooleanType verbose,
1798 const double cluster_threshold,const double smooth_threshold,
1799 ExceptionInfo *exception)
1802 previous_colorspace;
1811 *extrema[MaxDimension];
1814 *histogram[MaxDimension];
1817 Allocate histogram and extrema.
1819 assert(image != (Image *) NULL);
1820 assert(image->signature == MagickSignature);
1821 if (image->debug != MagickFalse)
1822 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1823 for (i=0; i < MaxDimension; i++)
1825 histogram[i]=(ssize_t *) AcquireQuantumMemory(256,sizeof(**histogram));
1826 extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1827 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1829 for (i-- ; i >= 0; i--)
1831 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1832 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1834 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1839 Initialize histogram.
1841 previous_colorspace=image->colorspace;
1842 (void) TransformImageColorspace(image,colorspace,exception);
1843 InitializeHistogram(image,histogram,exception);
1844 (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,
1845 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Red]);
1846 (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,
1847 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Green]);
1848 (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,
1849 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Blue]);
1851 Classify using the fuzzy c-Means technique.
1853 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose,
1855 (void) TransformImageColorspace(image,previous_colorspace,exception);
1857 Relinquish resources.
1859 for (i=0; i < MaxDimension; i++)
1861 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1862 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1868 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1872 + Z e r o C r o s s H i s t o g r a m %
1876 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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.
1882 % The format of the ZeroCrossHistogram method is:
1884 % ZeroCrossHistogram(double *second_derivative,
1885 % const double smooth_threshold,short *crossings)
1887 % A description of each parameter follows.
1889 % o second_derivative: Specifies an array of doubles representing the
1890 % second derivative of the histogram of a particular color component.
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.
1897 static void ZeroCrossHistogram(double *second_derivative,
1898 const double smooth_threshold,short *crossings)
1907 Merge low numbers to zero to help prevent noise.
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;
1914 Mark zero crossings.
1917 for (i=0; i <= 255; i++)
1920 if (second_derivative[i] < 0.0)
1927 if (second_derivative[i] > 0.0)