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-2010 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 "magick/studio.h"
86 #include "magick/cache.h"
87 #include "magick/color.h"
88 #include "magick/colorspace.h"
89 #include "magick/exception.h"
90 #include "magick/exception-private.h"
91 #include "magick/image.h"
92 #include "magick/image-private.h"
93 #include "magick/memory_.h"
94 #include "magick/monitor.h"
95 #include "magick/monitor-private.h"
96 #include "magick/quantize.h"
97 #include "magick/quantum.h"
98 #include "magick/quantum-private.h"
99 #include "magick/segment.h"
100 #include "magick/string_.h"
105 #define MaxDimension 3
106 #define DeltaTau 0.5f
107 #if defined(FastClassify)
108 #define WeightingExponent 2.0
109 #define SegmentPower(ratio) (ratio)
111 #define WeightingExponent 2.5
112 #define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0)));
117 Typedef declarations.
119 typedef struct _ExtentPacket
130 typedef struct _Cluster
145 typedef struct _IntervalTree
163 typedef struct _ZeroCrossing
174 Constant declarations.
186 static MagickRealType
187 OptimalTau(const long *,const double,const double,const double,
188 const double,short *);
191 DefineRegion(const short *,ExtentPacket *);
194 InitializeHistogram(const Image *,long **,ExceptionInfo *),
195 ScaleSpace(const long *,const MagickRealType,MagickRealType *),
196 ZeroCrossHistogram(MagickRealType *,const MagickRealType,short *);
199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
207 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
209 % Classify() defines one or more classes. Each pixel is thresholded to
210 % determine which class it belongs to. If the class is not identified it is
211 % assigned to the closest class based on the fuzzy c-Means technique.
213 % The format of the Classify method is:
215 % MagickBooleanType Classify(Image *image,short **extrema,
216 % const MagickRealType cluster_threshold,
217 % const MagickRealType weighting_exponent,
218 % const MagickBooleanType verbose)
220 % A description of each parameter follows.
222 % o image: the image.
224 % o extrema: Specifies a pointer to an array of integers. They
225 % represent the peaks and valleys of the histogram for each color
228 % o cluster_threshold: This MagickRealType represents the minimum number of
229 % pixels contained in a hexahedra before it can be considered valid
230 % (expressed as a percentage).
232 % o weighting_exponent: Specifies the membership weighting exponent.
234 % o verbose: A value greater than zero prints detailed information about
235 % the identified classes.
238 static MagickBooleanType Classify(Image *image,short **extrema,
239 const MagickRealType cluster_threshold,
240 const MagickRealType weighting_exponent,const MagickBooleanType verbose)
242 #define SegmentImageTag "Segment/Image"
275 register MagickRealType
284 cluster=(Cluster *) NULL;
285 head=(Cluster *) NULL;
286 (void) ResetMagickMemory(&red,0,sizeof(red));
287 (void) ResetMagickMemory(&green,0,sizeof(green));
288 (void) ResetMagickMemory(&blue,0,sizeof(blue));
289 while (DefineRegion(extrema[Red],&red) != 0)
292 while (DefineRegion(extrema[Green],&green) != 0)
295 while (DefineRegion(extrema[Blue],&blue) != 0)
298 Allocate a new class.
300 if (head != (Cluster *) NULL)
302 cluster->next=(Cluster *) AcquireMagickMemory(
303 sizeof(*cluster->next));
304 cluster=cluster->next;
308 cluster=(Cluster *) AcquireAlignedMemory(1,sizeof(*cluster));
311 if (cluster == (Cluster *) NULL)
312 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
315 Initialize a new class.
319 cluster->green=green;
321 cluster->next=(Cluster *) NULL;
325 if (head == (Cluster *) NULL)
328 No classes were identified-- create one.
330 cluster=(Cluster *) AcquireAlignedMemory(1,sizeof(*cluster));
331 if (cluster == (Cluster *) NULL)
332 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
335 Initialize a new class.
339 cluster->green=green;
341 cluster->next=(Cluster *) NULL;
345 Count the pixels for each cluster.
350 exception=(&image->exception);
351 image_view=AcquireCacheView(image);
352 for (y=0; y < (long) image->rows; y++)
354 register const PixelPacket
360 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
361 if (p == (const PixelPacket *) NULL)
363 for (x=0; x < (long) image->columns; x++)
365 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
366 if (((long) ScaleQuantumToChar(GetRedSample(p)) >=
367 (cluster->red.left-SafeMargin)) &&
368 ((long) ScaleQuantumToChar(GetRedSample(p)) <=
369 (cluster->red.right+SafeMargin)) &&
370 ((long) ScaleQuantumToChar(GetGreenSample(p)) >=
371 (cluster->green.left-SafeMargin)) &&
372 ((long) ScaleQuantumToChar(GetGreenSample(p)) <=
373 (cluster->green.right+SafeMargin)) &&
374 ((long) ScaleQuantumToChar(GetBlueSample(p)) >=
375 (cluster->blue.left-SafeMargin)) &&
376 ((long) ScaleQuantumToChar(GetBlueSample(p)) <=
377 (cluster->blue.right+SafeMargin)))
383 cluster->red.center+=(MagickRealType) ScaleQuantumToChar(GetRedSample(p));
384 cluster->green.center+=(MagickRealType)
385 ScaleQuantumToChar(GetGreenSample(p));
386 cluster->blue.center+=(MagickRealType) ScaleQuantumToChar(GetBlueSample(p));
392 if (image->progress_monitor != (MagickProgressMonitor) NULL)
397 #if defined(MAGICKCORE_OPENMP_SUPPORT)
398 #pragma omp critical (MagickCore_Classify)
400 proceed=SetImageProgress(image,SegmentImageTag,progress++,
402 if (proceed == MagickFalse)
406 image_view=DestroyCacheView(image_view);
408 Remove clusters that do not meet minimum cluster threshold.
413 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
415 next_cluster=cluster->next;
416 if ((cluster->count > 0) &&
417 (cluster->count >= (count*cluster_threshold/100.0)))
423 cluster->red.center/=cluster->count;
424 cluster->green.center/=cluster->count;
425 cluster->blue.center/=cluster->count;
427 last_cluster=cluster;
436 last_cluster->next=next_cluster;
437 cluster=(Cluster *) RelinquishMagickMemory(cluster);
439 number_clusters=(unsigned long) count;
440 if (verbose != MagickFalse)
443 Print cluster statistics.
445 (void) fprintf(stdout,"Fuzzy C-means Statistics\n");
446 (void) fprintf(stdout,"===================\n\n");
447 (void) fprintf(stdout,"\tCluster Threshold = %.15g\n",(double)
449 (void) fprintf(stdout,"\tWeighting Exponent = %.15g\n",(double)
451 (void) fprintf(stdout,"\tTotal Number of Clusters = %lu\n\n",
454 Print the total number of points per cluster.
456 (void) fprintf(stdout,"\n\nNumber of Vectors Per Cluster\n");
457 (void) fprintf(stdout,"=============================\n\n");
458 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
459 (void) fprintf(stdout,"Cluster #%ld = %ld\n",cluster->id,
462 Print the cluster extents.
464 (void) fprintf(stdout,
465 "\n\n\nCluster Extents: (Vector Size: %d)\n",MaxDimension);
466 (void) fprintf(stdout,"================");
467 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
469 (void) fprintf(stdout,"\n\nCluster #%ld\n\n",cluster->id);
470 (void) fprintf(stdout,"%ld-%ld %ld-%ld %ld-%ld\n",cluster->red.left,
471 cluster->red.right,cluster->green.left,cluster->green.right,
472 cluster->blue.left,cluster->blue.right);
475 Print the cluster center values.
477 (void) fprintf(stdout,
478 "\n\n\nCluster Center Values: (Vector Size: %d)\n",MaxDimension);
479 (void) fprintf(stdout,"=====================");
480 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
482 (void) fprintf(stdout,"\n\nCluster #%ld\n\n",cluster->id);
483 (void) fprintf(stdout,"%.15g %.15g %.15g\n",(double)
484 cluster->red.center,(double) cluster->green.center,(double)
485 cluster->blue.center);
487 (void) fprintf(stdout,"\n");
489 if (number_clusters > 256)
490 ThrowBinaryException(ImageError,"TooManyClusters",image->filename);
492 Speed up distance calculations.
494 squares=(MagickRealType *) AcquireQuantumMemory(513UL,sizeof(*squares));
495 if (squares == (MagickRealType *) NULL)
496 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
499 for (i=(-255); i <= 255; i++)
500 squares[i]=(MagickRealType) i*(MagickRealType) i;
502 Allocate image colormap.
504 if (AcquireImageColormap(image,number_clusters) == MagickFalse)
505 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
508 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
510 image->colormap[i].red=ScaleCharToQuantum((unsigned char)
511 (cluster->red.center+0.5));
512 image->colormap[i].green=ScaleCharToQuantum((unsigned char)
513 (cluster->green.center+0.5));
514 image->colormap[i].blue=ScaleCharToQuantum((unsigned char)
515 (cluster->blue.center+0.5));
519 Do course grain classes.
521 exception=(&image->exception);
522 image_view=AcquireCacheView(image);
523 #if defined(MAGICKCORE_OPENMP_SUPPORT)
524 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
526 for (y=0; y < (long) image->rows; y++)
531 register const PixelPacket
543 if (status == MagickFalse)
545 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
546 if (q == (PixelPacket *) NULL)
551 indexes=GetCacheViewAuthenticIndexQueue(image_view);
552 for (x=0; x < (long) image->columns; x++)
554 indexes[x]=(IndexPacket) 0;
555 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
557 if (((long) ScaleQuantumToChar(q->red) >=
558 (cluster->red.left-SafeMargin)) &&
559 ((long) ScaleQuantumToChar(q->red) <=
560 (cluster->red.right+SafeMargin)) &&
561 ((long) ScaleQuantumToChar(q->green) >=
562 (cluster->green.left-SafeMargin)) &&
563 ((long) ScaleQuantumToChar(q->green) <=
564 (cluster->green.right+SafeMargin)) &&
565 ((long) ScaleQuantumToChar(q->blue) >=
566 (cluster->blue.left-SafeMargin)) &&
567 ((long) ScaleQuantumToChar(q->blue) <=
568 (cluster->blue.right+SafeMargin)))
573 indexes[x]=(IndexPacket) cluster->id;
577 if (cluster == (Cluster *) NULL)
591 Compute fuzzy membership.
594 for (j=0; j < (long) image->colors; j++)
598 distance_squared=squares[(long) ScaleQuantumToChar(q->red)-
599 (long) ScaleQuantumToChar(GetRedSample(p))]+
600 squares[(long) ScaleQuantumToChar(q->green)-
601 (long) ScaleQuantumToChar(GetGreenSample(p))]+
602 squares[(long) ScaleQuantumToChar(q->blue)-
603 (long) ScaleQuantumToChar(GetBlueSample(p))];
604 numerator=distance_squared;
605 for (k=0; k < (long) image->colors; k++)
608 distance_squared=squares[(long) ScaleQuantumToChar(q->red)-
609 (long) ScaleQuantumToChar(GetRedSample(p))]+
610 squares[(long) ScaleQuantumToChar(q->green)-
611 (long) ScaleQuantumToChar(GetGreenSample(p))]+
612 squares[(long) ScaleQuantumToChar(q->blue)-
613 (long) ScaleQuantumToChar(GetBlueSample(p))];
614 ratio=numerator/distance_squared;
615 sum+=SegmentPower(ratio);
617 if ((sum != 0.0) && ((1.0/sum) > local_minima))
622 local_minima=1.0/sum;
623 indexes[x]=(IndexPacket) j;
629 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
631 if (image->progress_monitor != (MagickProgressMonitor) NULL)
636 #if defined(MAGICKCORE_OPENMP_SUPPORT)
637 #pragma omp critical (MagickCore_Classify)
639 proceed=SetImageProgress(image,SegmentImageTag,progress++,
641 if (proceed == MagickFalse)
645 image_view=DestroyCacheView(image_view);
646 status&=SyncImage(image);
648 Relinquish resources.
650 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
652 next_cluster=cluster->next;
653 cluster=(Cluster *) RelinquishMagickMemory(cluster);
656 free_squares=squares;
657 free_squares=(MagickRealType *) RelinquishMagickMemory(free_squares);
662 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
666 + C o n s o l i d a t e C r o s s i n g s %
670 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
672 % ConsolidateCrossings() guarantees that an even number of zero crossings
673 % always lie between two crossings.
675 % The format of the ConsolidateCrossings method is:
677 % ConsolidateCrossings(ZeroCrossing *zero_crossing,
678 % const unsigned long number_crossings)
680 % A description of each parameter follows.
682 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
684 % o number_crossings: This unsigned long specifies the number of elements
685 % in the zero_crossing array.
689 static inline long MagickAbsoluteValue(const long x)
696 static inline long MagickMax(const long x,const long y)
703 static inline long MagickMin(const long x,const long y)
710 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
711 const unsigned long number_crossings)
727 Consolidate zero crossings.
729 for (i=(long) number_crossings-1; i >= 0; i--)
730 for (j=0; j <= 255; j++)
732 if (zero_crossing[i].crossings[j] == 0)
735 Find the entry that is closest to j and still preserves the
736 property that there are an even number of crossings between
739 for (k=j-1; k > 0; k--)
740 if (zero_crossing[i+1].crossings[k] != 0)
744 for (k=j+1; k < 255; k++)
745 if (zero_crossing[i+1].crossings[k] != 0)
747 right=MagickMin(k,255);
749 K is the zero crossing just left of j.
751 for (k=j-1; k > 0; k--)
752 if (zero_crossing[i].crossings[k] != 0)
757 Check center for an even number of crossings between k and j.
760 if (zero_crossing[i+1].crossings[j] != 0)
763 for (l=k+1; l < center; l++)
764 if (zero_crossing[i+1].crossings[l] != 0)
766 if (((count % 2) == 0) && (center != k))
770 Check left for an even number of crossings between k and j.
775 for (l=k+1; l < left; l++)
776 if (zero_crossing[i+1].crossings[l] != 0)
778 if (((count % 2) == 0) && (left != k))
782 Check right for an even number of crossings between k and j.
787 for (l=k+1; l < right; l++)
788 if (zero_crossing[i+1].crossings[l] != 0)
790 if (((count % 2) == 0) && (right != k))
793 l=zero_crossing[i].crossings[j];
794 zero_crossing[i].crossings[j]=0;
796 zero_crossing[i].crossings[correct]=(short) l;
801 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
805 + D e f i n e R e g i o n %
809 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
811 % DefineRegion() defines the left and right boundaries of a peak region.
813 % The format of the DefineRegion method is:
815 % long DefineRegion(const short *extrema,ExtentPacket *extents)
817 % A description of each parameter follows.
819 % o extrema: Specifies a pointer to an array of integers. They
820 % represent the peaks and valleys of the histogram for each color
823 % o extents: This pointer to an ExtentPacket represent the extends
824 % of a particular peak or valley of a color component.
827 static long DefineRegion(const short *extrema,ExtentPacket *extents)
830 Initialize to default values.
836 Find the left side (maxima).
838 for ( ; extents->index <= 255; extents->index++)
839 if (extrema[extents->index] > 0)
841 if (extents->index > 255)
842 return(MagickFalse); /* no left side - no region exists */
843 extents->left=extents->index;
845 Find the right side (minima).
847 for ( ; extents->index <= 255; extents->index++)
848 if (extrema[extents->index] < 0)
850 extents->right=extents->index-1;
855 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
859 + D e r i v a t i v e H i s t o g r a m %
863 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
865 % DerivativeHistogram() determines the derivative of the histogram using
866 % central differencing.
868 % The format of the DerivativeHistogram method is:
870 % DerivativeHistogram(const MagickRealType *histogram,
871 % MagickRealType *derivative)
873 % A description of each parameter follows.
875 % o histogram: Specifies an array of MagickRealTypes representing the number
876 % of pixels for each intensity of a particular color component.
878 % o derivative: This array of MagickRealTypes is initialized by
879 % DerivativeHistogram to the derivative of the histogram using central
883 static void DerivativeHistogram(const MagickRealType *histogram,
884 MagickRealType *derivative)
891 Compute endpoints using second order polynomial interpolation.
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]);
897 Compute derivative using central differencing.
899 for (i=1; i < n; i++)
900 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
905 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %
913 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
915 % GetImageDynamicThreshold() returns the dynamic threshold for an image.
917 % The format of the GetImageDynamicThreshold method is:
919 % MagickBooleanType GetImageDynamicThreshold(const Image *image,
920 % const double cluster_threshold,const double smooth_threshold,
921 % MagickPixelPacket *pixel,ExceptionInfo *exception)
923 % A description of each parameter follows.
925 % o image: the image.
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).
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.
935 % o pixel: return the dynamic threshold here.
937 % o exception: return any errors or warnings in this structure.
940 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,
941 const double cluster_threshold,const double smooth_threshold,
942 MagickPixelPacket *pixel,ExceptionInfo *exception)
959 *histogram[MaxDimension],
968 register const PixelPacket
976 *extrema[MaxDimension];
979 Allocate histogram and extrema.
981 assert(image != (Image *) NULL);
982 assert(image->signature == MagickSignature);
983 if (image->debug != MagickFalse)
984 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
985 GetMagickPixelPacket(image,pixel);
986 for (i=0; i < MaxDimension; i++)
988 histogram[i]=(long *) AcquireQuantumMemory(256UL,sizeof(**histogram));
989 extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
990 if ((histogram[i] == (long *) NULL) || (extrema[i] == (short *) NULL))
992 for (i-- ; i >= 0; i--)
994 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
995 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
997 (void) ThrowMagickException(exception,GetMagickModule(),
998 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1003 Initialize histogram.
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]);
1015 cluster=(Cluster *) NULL;
1016 head=(Cluster *) NULL;
1017 (void) ResetMagickMemory(&red,0,sizeof(red));
1018 (void) ResetMagickMemory(&green,0,sizeof(green));
1019 (void) ResetMagickMemory(&blue,0,sizeof(blue));
1020 while (DefineRegion(extrema[Red],&red) != 0)
1023 while (DefineRegion(extrema[Green],&green) != 0)
1026 while (DefineRegion(extrema[Blue],&blue) != 0)
1029 Allocate a new class.
1031 if (head != (Cluster *) NULL)
1033 cluster->next=(Cluster *) AcquireMagickMemory(
1034 sizeof(*cluster->next));
1035 cluster=cluster->next;
1039 cluster=(Cluster *) AcquireAlignedMemory(1,sizeof(*cluster));
1042 if (cluster == (Cluster *) NULL)
1044 (void) ThrowMagickException(exception,GetMagickModule(),
1045 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1047 return(MagickFalse);
1050 Initialize a new class.
1054 cluster->green=green;
1056 cluster->next=(Cluster *) NULL;
1060 if (head == (Cluster *) NULL)
1063 No classes were identified-- create one.
1065 cluster=(Cluster *) AcquireAlignedMemory(1,sizeof(*cluster));
1066 if (cluster == (Cluster *) NULL)
1068 (void) ThrowMagickException(exception,GetMagickModule(),
1069 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1070 return(MagickFalse);
1073 Initialize a new class.
1077 cluster->green=green;
1079 cluster->next=(Cluster *) NULL;
1083 Count the pixels for each cluster.
1086 for (y=0; y < (long) image->rows; y++)
1088 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1089 if (p == (const PixelPacket *) NULL)
1091 for (x=0; x < (long) image->columns; x++)
1093 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1094 if (((long) ScaleQuantumToChar(GetRedSample(p)) >=
1095 (cluster->red.left-SafeMargin)) &&
1096 ((long) ScaleQuantumToChar(GetRedSample(p)) <=
1097 (cluster->red.right+SafeMargin)) &&
1098 ((long) ScaleQuantumToChar(GetGreenSample(p)) >=
1099 (cluster->green.left-SafeMargin)) &&
1100 ((long) ScaleQuantumToChar(GetGreenSample(p)) <=
1101 (cluster->green.right+SafeMargin)) &&
1102 ((long) ScaleQuantumToChar(GetBlueSample(p)) >=
1103 (cluster->blue.left-SafeMargin)) &&
1104 ((long) ScaleQuantumToChar(GetBlueSample(p)) <=
1105 (cluster->blue.right+SafeMargin)))
1111 cluster->red.center+=(MagickRealType)
1112 ScaleQuantumToChar(GetRedSample(p));
1113 cluster->green.center+=(MagickRealType)
1114 ScaleQuantumToChar(GetGreenSample(p));
1115 cluster->blue.center+=(MagickRealType)
1116 ScaleQuantumToChar(GetBlueSample(p));
1122 proceed=SetImageProgress(image,SegmentImageTag,y,2*image->rows);
1123 if (proceed == MagickFalse)
1127 Remove clusters that do not meet minimum cluster threshold.
1132 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1134 next_cluster=cluster->next;
1135 if ((cluster->count > 0) &&
1136 (cluster->count >= (count*cluster_threshold/100.0)))
1142 cluster->red.center/=cluster->count;
1143 cluster->green.center/=cluster->count;
1144 cluster->blue.center/=cluster->count;
1146 last_cluster=cluster;
1152 if (cluster == head)
1155 last_cluster->next=next_cluster;
1156 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1163 for (cluster=object; cluster->next != (Cluster *) NULL; )
1165 if (cluster->count < object->count)
1167 cluster=cluster->next;
1169 background=head->next;
1170 for (cluster=background; cluster->next != (Cluster *) NULL; )
1172 if (cluster->count > background->count)
1174 cluster=cluster->next;
1177 threshold=(background->red.center+object->red.center)/2.0;
1178 pixel->red=(MagickRealType) ScaleCharToQuantum((unsigned char)
1180 threshold=(background->green.center+object->green.center)/2.0;
1181 pixel->green=(MagickRealType) ScaleCharToQuantum((unsigned char)
1183 threshold=(background->blue.center+object->blue.center)/2.0;
1184 pixel->blue=(MagickRealType) ScaleCharToQuantum((unsigned char)
1187 Relinquish resources.
1189 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1191 next_cluster=cluster->next;
1192 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1194 for (i=0; i < MaxDimension; i++)
1196 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1197 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1203 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1207 + I n i t i a l i z e H i s t o g r a m %
1211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1213 % InitializeHistogram() computes the histogram for an image.
1215 % The format of the InitializeHistogram method is:
1217 % InitializeHistogram(const Image *image,long **histogram)
1219 % A description of each parameter follows.
1221 % o image: Specifies a pointer to an Image structure; returned from
1224 % o histogram: Specifies an array of integers representing the number
1225 % of pixels for each intensity of a particular color component.
1228 static void InitializeHistogram(const Image *image,long **histogram,
1229 ExceptionInfo *exception)
1234 register const PixelPacket
1242 Initialize histogram.
1244 for (i=0; i <= 255; i++)
1246 histogram[Red][i]=0;
1247 histogram[Green][i]=0;
1248 histogram[Blue][i]=0;
1250 for (y=0; y < (long) image->rows; y++)
1252 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1253 if (p == (const PixelPacket *) NULL)
1255 for (x=0; x < (long) image->columns; x++)
1257 histogram[Red][(long) ScaleQuantumToChar(GetRedSample(p))]++;
1258 histogram[Green][(long) ScaleQuantumToChar(GetGreenSample(p))]++;
1259 histogram[Blue][(long) ScaleQuantumToChar(GetBlueSample(p))]++;
1266 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1270 + I n i t i a l i z e I n t e r v a l T r e e %
1274 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1276 % InitializeIntervalTree() initializes an interval tree from the lists of
1279 % The format of the InitializeIntervalTree method is:
1281 % InitializeIntervalTree(IntervalTree **list,long *number_nodes,
1282 % IntervalTree *node)
1284 % A description of each parameter follows.
1286 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
1288 % o number_crossings: This unsigned long specifies the number of elements
1289 % in the zero_crossing array.
1293 static void InitializeList(IntervalTree **list,long *number_nodes,
1296 if (node == (IntervalTree *) NULL)
1298 if (node->child == (IntervalTree *) NULL)
1299 list[(*number_nodes)++]=node;
1300 InitializeList(list,number_nodes,node->sibling);
1301 InitializeList(list,number_nodes,node->child);
1304 static void MeanStability(IntervalTree *node)
1306 register IntervalTree
1309 if (node == (IntervalTree *) NULL)
1311 node->mean_stability=0.0;
1313 if (child != (IntervalTree *) NULL)
1318 register MagickRealType
1323 for ( ; child != (IntervalTree *) NULL; child=child->sibling)
1325 sum+=child->stability;
1328 node->mean_stability=sum/(MagickRealType) count;
1330 MeanStability(node->sibling);
1331 MeanStability(node->child);
1334 static void Stability(IntervalTree *node)
1336 if (node == (IntervalTree *) NULL)
1338 if (node->child == (IntervalTree *) NULL)
1339 node->stability=0.0;
1341 node->stability=node->tau-(node->child)->tau;
1342 Stability(node->sibling);
1343 Stability(node->child);
1346 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing,
1347 const unsigned long number_crossings)
1365 Allocate interval tree.
1367 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1369 if (list == (IntervalTree **) NULL)
1370 return((IntervalTree *) NULL);
1372 The root is the entire histogram.
1374 root=(IntervalTree *) AcquireAlignedMemory(1,sizeof(*root));
1375 root->child=(IntervalTree *) NULL;
1376 root->sibling=(IntervalTree *) NULL;
1380 for (i=(-1); i < (long) number_crossings; i++)
1383 Initialize list with all nodes with no children.
1386 InitializeList(list,&number_nodes,root);
1390 for (j=0; j < number_nodes; j++)
1395 for (k=head->left+1; k < head->right; k++)
1397 if (zero_crossing[i+1].crossings[k] != 0)
1401 node->child=(IntervalTree *) AcquireMagickMemory(
1402 sizeof(*node->child));
1407 node->sibling=(IntervalTree *) AcquireMagickMemory(
1408 sizeof(*node->sibling));
1411 node->tau=zero_crossing[i+1].tau;
1412 node->child=(IntervalTree *) NULL;
1413 node->sibling=(IntervalTree *) NULL;
1419 if (left != head->left)
1421 node->sibling=(IntervalTree *) AcquireMagickMemory(
1422 sizeof(*node->sibling));
1424 node->tau=zero_crossing[i+1].tau;
1425 node->child=(IntervalTree *) NULL;
1426 node->sibling=(IntervalTree *) NULL;
1428 node->right=head->right;
1433 Determine the stability: difference between a nodes tau and its child.
1435 Stability(root->child);
1436 MeanStability(root->child);
1437 list=(IntervalTree **) RelinquishMagickMemory(list);
1442 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1446 + O p t i m a l T a u %
1450 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1452 % OptimalTau() finds the optimal tau for each band of the histogram.
1454 % The format of the OptimalTau method is:
1456 % MagickRealType OptimalTau(const long *histogram,const double max_tau,
1457 % const double min_tau,const double delta_tau,
1458 % const double smooth_threshold,short *extrema)
1460 % A description of each parameter follows.
1462 % o histogram: Specifies an array of integers representing the number
1463 % of pixels for each intensity of a particular color component.
1465 % o extrema: Specifies a pointer to an array of integers. They
1466 % represent the peaks and valleys of the histogram for each color
1471 static void ActiveNodes(IntervalTree **list,long *number_nodes,
1474 if (node == (IntervalTree *) NULL)
1476 if (node->stability >= node->mean_stability)
1478 list[(*number_nodes)++]=node;
1479 ActiveNodes(list,number_nodes,node->sibling);
1483 ActiveNodes(list,number_nodes,node->sibling);
1484 ActiveNodes(list,number_nodes,node->child);
1488 static void FreeNodes(IntervalTree *node)
1490 if (node == (IntervalTree *) NULL)
1492 FreeNodes(node->sibling);
1493 FreeNodes(node->child);
1494 node=(IntervalTree *) RelinquishMagickMemory(node);
1497 static MagickRealType OptimalTau(const long *histogram,const double max_tau,
1498 const double min_tau,const double delta_tau,const double smooth_threshold,
1534 Allocate interval tree.
1536 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1538 if (list == (IntervalTree **) NULL)
1541 Allocate zero crossing list.
1543 count=(unsigned long) ((max_tau-min_tau)/delta_tau)+2;
1544 zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count,
1545 sizeof(*zero_crossing));
1546 if (zero_crossing == (ZeroCrossing *) NULL)
1548 for (i=0; i < (long) count; i++)
1549 zero_crossing[i].tau=(-1.0);
1551 Initialize zero crossing list.
1553 derivative=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*derivative));
1554 second_derivative=(MagickRealType *) AcquireQuantumMemory(256,
1555 sizeof(*second_derivative));
1556 if ((derivative == (MagickRealType *) NULL) ||
1557 (second_derivative == (MagickRealType *) NULL))
1558 ThrowFatalException(ResourceLimitFatalError,
1559 "UnableToAllocateDerivatives");
1561 for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
1563 zero_crossing[i].tau=tau;
1564 ScaleSpace(histogram,tau,zero_crossing[i].histogram);
1565 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1566 DerivativeHistogram(derivative,second_derivative);
1567 ZeroCrossHistogram(second_derivative,smooth_threshold,
1568 zero_crossing[i].crossings);
1572 Add an entry for the original histogram.
1574 zero_crossing[i].tau=0.0;
1575 for (j=0; j <= 255; j++)
1576 zero_crossing[i].histogram[j]=(MagickRealType) histogram[j];
1577 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1578 DerivativeHistogram(derivative,second_derivative);
1579 ZeroCrossHistogram(second_derivative,smooth_threshold,
1580 zero_crossing[i].crossings);
1581 number_crossings=(unsigned long) i;
1582 derivative=(MagickRealType *) RelinquishMagickMemory(derivative);
1583 second_derivative=(MagickRealType *)
1584 RelinquishMagickMemory(second_derivative);
1586 Ensure the scale-space fingerprints form lines in scale-space, not loops.
1588 ConsolidateCrossings(zero_crossing,number_crossings);
1590 Force endpoints to be included in the interval.
1592 for (i=0; i <= (long) number_crossings; i++)
1594 for (j=0; j < 255; j++)
1595 if (zero_crossing[i].crossings[j] != 0)
1597 zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]);
1598 for (j=255; j > 0; j--)
1599 if (zero_crossing[i].crossings[j] != 0)
1601 zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1604 Initialize interval tree.
1606 root=InitializeIntervalTree(zero_crossing,number_crossings);
1607 if (root == (IntervalTree *) NULL)
1610 Find active nodes: stability is greater (or equal) to the mean stability of
1614 ActiveNodes(list,&number_nodes,root->child);
1618 for (i=0; i <= 255; i++)
1620 for (i=0; i < number_nodes; i++)
1623 Find this tau in zero crossings list.
1627 for (j=0; j <= (long) number_crossings; j++)
1628 if (zero_crossing[j].tau == node->tau)
1631 Find the value of the peak.
1633 peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1636 value=zero_crossing[k].histogram[index];
1637 for (x=node->left; x <= node->right; x++)
1639 if (peak != MagickFalse)
1641 if (zero_crossing[k].histogram[x] > value)
1643 value=zero_crossing[k].histogram[x];
1648 if (zero_crossing[k].histogram[x] < value)
1650 value=zero_crossing[k].histogram[x];
1654 for (x=node->left; x <= node->right; x++)
1658 if (peak != MagickFalse)
1659 extrema[x]=(short) index;
1661 extrema[x]=(short) (-index);
1665 Determine the average tau.
1668 for (i=0; i < number_nodes; i++)
1669 average_tau+=list[i]->tau;
1670 average_tau/=(MagickRealType) number_nodes;
1672 Relinquish resources.
1675 zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1676 list=(IntervalTree **) RelinquishMagickMemory(list);
1677 return(average_tau);
1681 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1685 + S c a l e S p a c e %
1689 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1691 % ScaleSpace() performs a scale-space filter on the 1D histogram.
1693 % The format of the ScaleSpace method is:
1695 % ScaleSpace(const long *histogram,const MagickRealType tau,
1696 % MagickRealType *scale_histogram)
1698 % A description of each parameter follows.
1700 % o histogram: Specifies an array of MagickRealTypes representing the number
1701 % of pixels for each intensity of a particular color component.
1705 static void ScaleSpace(const long *histogram,const MagickRealType tau,
1706 MagickRealType *scale_histogram)
1718 gamma=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*gamma));
1719 if (gamma == (MagickRealType *) NULL)
1720 ThrowFatalException(ResourceLimitFatalError,
1721 "UnableToAllocateGammaMap");
1722 alpha=1.0/(tau*sqrt(2.0*MagickPI));
1723 beta=(-1.0/(2.0*tau*tau));
1724 for (x=0; x <= 255; x++)
1726 for (x=0; x <= 255; x++)
1728 gamma[x]=exp((double) beta*x*x);
1729 if (gamma[x] < MagickEpsilon)
1732 for (x=0; x <= 255; x++)
1735 for (u=0; u <= 255; u++)
1736 sum+=(MagickRealType) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1737 scale_histogram[x]=alpha*sum;
1739 gamma=(MagickRealType *) RelinquishMagickMemory(gamma);
1743 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1747 % S e g m e n t I m a g e %
1751 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1753 % SegmentImage() segment an image by analyzing the histograms of the color
1754 % components and identifying units that are homogeneous with the fuzzy
1755 % C-means technique.
1757 % The format of the SegmentImage method is:
1759 % MagickBooleanType SegmentImage(Image *image,
1760 % const ColorspaceType colorspace,const MagickBooleanType verbose,
1761 % const double cluster_threshold,const double smooth_threshold)
1763 % A description of each parameter follows.
1765 % o image: the image.
1767 % o colorspace: Indicate the colorspace.
1769 % o verbose: Set to MagickTrue to print detailed information about the
1770 % identified classes.
1772 % o cluster_threshold: This represents the minimum number of pixels
1773 % contained in a hexahedra before it can be considered valid (expressed
1776 % o smooth_threshold: the smoothing threshold eliminates noise in the second
1777 % derivative of the histogram. As the value is increased, you can expect a
1778 % smoother second derivative.
1781 MagickExport MagickBooleanType SegmentImage(Image *image,
1782 const ColorspaceType colorspace,const MagickBooleanType verbose,
1783 const double cluster_threshold,const double smooth_threshold)
1786 *histogram[MaxDimension];
1795 *extrema[MaxDimension];
1798 Allocate histogram and extrema.
1800 assert(image != (Image *) NULL);
1801 assert(image->signature == MagickSignature);
1802 if (image->debug != MagickFalse)
1803 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1804 for (i=0; i < MaxDimension; i++)
1806 histogram[i]=(long *) AcquireQuantumMemory(256,sizeof(**histogram));
1807 extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1808 if ((histogram[i] == (long *) NULL) || (extrema[i] == (short *) NULL))
1810 for (i-- ; i >= 0; i--)
1812 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1813 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1815 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1819 if (colorspace != RGBColorspace)
1820 (void) TransformImageColorspace(image,colorspace);
1822 Initialize histogram.
1824 InitializeHistogram(image,histogram,&image->exception);
1825 (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,
1826 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Red]);
1827 (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,
1828 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Green]);
1829 (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,
1830 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Blue]);
1832 Classify using the fuzzy c-Means technique.
1834 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose);
1835 if (colorspace != RGBColorspace)
1836 (void) TransformImageColorspace(image,colorspace);
1838 Relinquish resources.
1840 for (i=0; i < MaxDimension; i++)
1842 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1843 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1849 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1853 + Z e r o C r o s s H i s t o g r a m %
1857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1859 % ZeroCrossHistogram() find the zero crossings in a histogram and marks
1860 % directions as: 1 is negative to positive; 0 is zero crossing; and -1
1861 % is positive to negative.
1863 % The format of the ZeroCrossHistogram method is:
1865 % ZeroCrossHistogram(MagickRealType *second_derivative,
1866 % const MagickRealType smooth_threshold,short *crossings)
1868 % A description of each parameter follows.
1870 % o second_derivative: Specifies an array of MagickRealTypes representing the
1871 % second derivative of the histogram of a particular color component.
1873 % o crossings: This array of integers is initialized with
1874 % -1, 0, or 1 representing the slope of the first derivative of the
1875 % of a particular color component.
1878 static void ZeroCrossHistogram(MagickRealType *second_derivative,
1879 const MagickRealType smooth_threshold,short *crossings)
1888 Merge low numbers to zero to help prevent noise.
1890 for (i=0; i <= 255; i++)
1891 if ((second_derivative[i] < smooth_threshold) &&
1892 (second_derivative[i] >= -smooth_threshold))
1893 second_derivative[i]=0.0;
1895 Mark zero crossings.
1898 for (i=0; i <= 255; i++)
1901 if (second_derivative[i] < 0.0)
1908 if (second_derivative[i] > 0.0)