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 "magick/studio.h"
86 #include "magick/cache.h"
87 #include "magick/color.h"
88 #include "magick/colormap.h"
89 #include "magick/colorspace.h"
90 #include "magick/exception.h"
91 #include "magick/exception-private.h"
92 #include "magick/image.h"
93 #include "magick/image-private.h"
94 #include "magick/memory_.h"
95 #include "magick/monitor.h"
96 #include "magick/monitor-private.h"
97 #include "magick/quantize.h"
98 #include "magick/quantum.h"
99 #include "magick/quantum-private.h"
100 #include "magick/segment.h"
101 #include "magick/string_.h"
106 #define MaxDimension 3
107 #define DeltaTau 0.5f
108 #if defined(FastClassify)
109 #define WeightingExponent 2.0
110 #define SegmentPower(ratio) (ratio)
112 #define WeightingExponent 2.5
113 #define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0)));
118 Typedef declarations.
120 typedef struct _ExtentPacket
131 typedef struct _Cluster
146 typedef struct _IntervalTree
164 typedef struct _ZeroCrossing
175 Constant declarations.
187 static MagickRealType
188 OptimalTau(const ssize_t *,const double,const double,const double,
189 const double,short *);
192 DefineRegion(const short *,ExtentPacket *);
195 InitializeHistogram(const Image *,ssize_t **,ExceptionInfo *),
196 ScaleSpace(const ssize_t *,const MagickRealType,MagickRealType *),
197 ZeroCrossHistogram(MagickRealType *,const MagickRealType,short *);
200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
208 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
210 % Classify() defines one or more classes. Each pixel is thresholded to
211 % determine which class it belongs to. If the class is not identified it is
212 % assigned to the closest class based on the fuzzy c-Means technique.
214 % The format of the Classify method is:
216 % MagickBooleanType Classify(Image *image,short **extrema,
217 % const MagickRealType cluster_threshold,
218 % const MagickRealType weighting_exponent,
219 % const MagickBooleanType verbose)
221 % A description of each parameter follows.
223 % o image: the image.
225 % o extrema: Specifies a pointer to an array of integers. They
226 % represent the peaks and valleys of the histogram for each color
229 % o cluster_threshold: This MagickRealType represents the minimum number of
230 % pixels contained in a hexahedra before it can be considered valid
231 % (expressed as a percentage).
233 % o weighting_exponent: Specifies the membership weighting exponent.
235 % o verbose: A value greater than zero prints detailed information about
236 % the identified classes.
239 static MagickBooleanType Classify(Image *image,short **extrema,
240 const MagickRealType cluster_threshold,
241 const MagickRealType weighting_exponent,const MagickBooleanType verbose)
243 #define SegmentImageTag "Segment/Image"
274 register MagickRealType
287 cluster=(Cluster *) NULL;
288 head=(Cluster *) NULL;
289 (void) ResetMagickMemory(&red,0,sizeof(red));
290 (void) ResetMagickMemory(&green,0,sizeof(green));
291 (void) ResetMagickMemory(&blue,0,sizeof(blue));
292 while (DefineRegion(extrema[Red],&red) != 0)
295 while (DefineRegion(extrema[Green],&green) != 0)
298 while (DefineRegion(extrema[Blue],&blue) != 0)
301 Allocate a new class.
303 if (head != (Cluster *) NULL)
305 cluster->next=(Cluster *) AcquireMagickMemory(
306 sizeof(*cluster->next));
307 cluster=cluster->next;
311 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
314 if (cluster == (Cluster *) NULL)
315 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
318 Initialize a new class.
322 cluster->green=green;
324 cluster->next=(Cluster *) NULL;
328 if (head == (Cluster *) NULL)
331 No classes were identified-- create one.
333 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
334 if (cluster == (Cluster *) NULL)
335 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
338 Initialize a new class.
342 cluster->green=green;
344 cluster->next=(Cluster *) NULL;
348 Count the pixels for each cluster.
353 exception=(&image->exception);
354 image_view=AcquireCacheView(image);
355 for (y=0; y < (ssize_t) image->rows; y++)
357 register const PixelPacket
363 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
364 if (p == (const PixelPacket *) NULL)
366 for (x=0; x < (ssize_t) image->columns; x++)
368 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
369 if (((ssize_t) ScaleQuantumToChar(GetRedPixelComponent(p)) >=
370 (cluster->red.left-SafeMargin)) &&
371 ((ssize_t) ScaleQuantumToChar(GetRedPixelComponent(p)) <=
372 (cluster->red.right+SafeMargin)) &&
373 ((ssize_t) ScaleQuantumToChar(GetGreenPixelComponent(p)) >=
374 (cluster->green.left-SafeMargin)) &&
375 ((ssize_t) ScaleQuantumToChar(GetGreenPixelComponent(p)) <=
376 (cluster->green.right+SafeMargin)) &&
377 ((ssize_t) ScaleQuantumToChar(GetBluePixelComponent(p)) >=
378 (cluster->blue.left-SafeMargin)) &&
379 ((ssize_t) ScaleQuantumToChar(GetBluePixelComponent(p)) <=
380 (cluster->blue.right+SafeMargin)))
386 cluster->red.center+=(MagickRealType) ScaleQuantumToChar(GetRedPixelComponent(p));
387 cluster->green.center+=(MagickRealType)
388 ScaleQuantumToChar(GetGreenPixelComponent(p));
389 cluster->blue.center+=(MagickRealType) ScaleQuantumToChar(GetBluePixelComponent(p));
395 if (image->progress_monitor != (MagickProgressMonitor) NULL)
400 #if defined(MAGICKCORE_OPENMP_SUPPORT)
401 #pragma omp critical (MagickCore_Classify)
403 proceed=SetImageProgress(image,SegmentImageTag,progress++,
405 if (proceed == MagickFalse)
409 image_view=DestroyCacheView(image_view);
411 Remove clusters that do not meet minimum cluster threshold.
416 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
418 next_cluster=cluster->next;
419 if ((cluster->count > 0) &&
420 (cluster->count >= (count*cluster_threshold/100.0)))
426 cluster->red.center/=cluster->count;
427 cluster->green.center/=cluster->count;
428 cluster->blue.center/=cluster->count;
430 last_cluster=cluster;
439 last_cluster->next=next_cluster;
440 cluster=(Cluster *) RelinquishMagickMemory(cluster);
442 number_clusters=(size_t) count;
443 if (verbose != MagickFalse)
446 Print cluster statistics.
448 (void) fprintf(stdout,"Fuzzy C-means Statistics\n");
449 (void) fprintf(stdout,"===================\n\n");
450 (void) fprintf(stdout,"\tCluster Threshold = %g\n",(double)
452 (void) fprintf(stdout,"\tWeighting Exponent = %g\n",(double)
454 (void) fprintf(stdout,"\tTotal Number of Clusters = %.20g\n\n",(double)
457 Print the total number of points per cluster.
459 (void) fprintf(stdout,"\n\nNumber of Vectors Per Cluster\n");
460 (void) fprintf(stdout,"=============================\n\n");
461 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
462 (void) fprintf(stdout,"Cluster #%.20g = %.20g\n",(double) cluster->id,
463 (double) cluster->count);
465 Print the cluster extents.
467 (void) fprintf(stdout,
468 "\n\n\nCluster Extents: (Vector Size: %d)\n",MaxDimension);
469 (void) fprintf(stdout,"================");
470 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
472 (void) fprintf(stdout,"\n\nCluster #%.20g\n\n",(double) cluster->id);
473 (void) fprintf(stdout,"%.20g-%.20g %.20g-%.20g %.20g-%.20g\n",(double)
474 cluster->red.left,(double) cluster->red.right,(double)
475 cluster->green.left,(double) cluster->green.right,(double)
476 cluster->blue.left,(double) cluster->blue.right);
479 Print the cluster center values.
481 (void) fprintf(stdout,
482 "\n\n\nCluster Center Values: (Vector Size: %d)\n",MaxDimension);
483 (void) fprintf(stdout,"=====================");
484 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
486 (void) fprintf(stdout,"\n\nCluster #%.20g\n\n",(double) cluster->id);
487 (void) fprintf(stdout,"%g %g %g\n",(double)
488 cluster->red.center,(double) cluster->green.center,(double)
489 cluster->blue.center);
491 (void) fprintf(stdout,"\n");
493 if (number_clusters > 256)
494 ThrowBinaryException(ImageError,"TooManyClusters",image->filename);
496 Speed up distance calculations.
498 squares=(MagickRealType *) AcquireQuantumMemory(513UL,sizeof(*squares));
499 if (squares == (MagickRealType *) NULL)
500 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
503 for (i=(-255); i <= 255; i++)
504 squares[i]=(MagickRealType) i*(MagickRealType) i;
506 Allocate image colormap.
508 if (AcquireImageColormap(image,number_clusters) == MagickFalse)
509 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
512 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
514 image->colormap[i].red=ScaleCharToQuantum((unsigned char)
515 (cluster->red.center+0.5));
516 image->colormap[i].green=ScaleCharToQuantum((unsigned char)
517 (cluster->green.center+0.5));
518 image->colormap[i].blue=ScaleCharToQuantum((unsigned char)
519 (cluster->blue.center+0.5));
523 Do course grain classes.
525 exception=(&image->exception);
526 image_view=AcquireCacheView(image);
527 #if defined(MAGICKCORE_OPENMP_SUPPORT)
528 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
530 for (y=0; y < (ssize_t) image->rows; y++)
535 register const PixelPacket
547 if (status == MagickFalse)
549 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
550 if (q == (PixelPacket *) NULL)
555 indexes=GetCacheViewAuthenticIndexQueue(image_view);
556 for (x=0; x < (ssize_t) image->columns; x++)
558 SetIndexPixelComponent(indexes+x,0);
559 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
561 if (((ssize_t) ScaleQuantumToChar(q->red) >=
562 (cluster->red.left-SafeMargin)) &&
563 ((ssize_t) ScaleQuantumToChar(q->red) <=
564 (cluster->red.right+SafeMargin)) &&
565 ((ssize_t) ScaleQuantumToChar(q->green) >=
566 (cluster->green.left-SafeMargin)) &&
567 ((ssize_t) ScaleQuantumToChar(q->green) <=
568 (cluster->green.right+SafeMargin)) &&
569 ((ssize_t) ScaleQuantumToChar(q->blue) >=
570 (cluster->blue.left-SafeMargin)) &&
571 ((ssize_t) ScaleQuantumToChar(q->blue) <=
572 (cluster->blue.right+SafeMargin)))
577 SetIndexPixelComponent(indexes+x,cluster->id);
581 if (cluster == (Cluster *) NULL)
595 Compute fuzzy membership.
598 for (j=0; j < (ssize_t) image->colors; j++)
602 distance_squared=squares[(ssize_t) ScaleQuantumToChar(q->red)-
603 (ssize_t) ScaleQuantumToChar(GetRedPixelComponent(p))]+
604 squares[(ssize_t) ScaleQuantumToChar(q->green)-
605 (ssize_t) ScaleQuantumToChar(GetGreenPixelComponent(p))]+
606 squares[(ssize_t) ScaleQuantumToChar(q->blue)-
607 (ssize_t) ScaleQuantumToChar(GetBluePixelComponent(p))];
608 numerator=distance_squared;
609 for (k=0; k < (ssize_t) image->colors; k++)
612 distance_squared=squares[(ssize_t) ScaleQuantumToChar(q->red)-
613 (ssize_t) ScaleQuantumToChar(GetRedPixelComponent(p))]+
614 squares[(ssize_t) ScaleQuantumToChar(q->green)-
615 (ssize_t) ScaleQuantumToChar(GetGreenPixelComponent(p))]+
616 squares[(ssize_t) ScaleQuantumToChar(q->blue)-
617 (ssize_t) ScaleQuantumToChar(GetBluePixelComponent(p))];
618 ratio=numerator/distance_squared;
619 sum+=SegmentPower(ratio);
621 if ((sum != 0.0) && ((1.0/sum) > local_minima))
626 local_minima=1.0/sum;
627 SetIndexPixelComponent(indexes+x,j);
633 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
635 if (image->progress_monitor != (MagickProgressMonitor) NULL)
640 #if defined(MAGICKCORE_OPENMP_SUPPORT)
641 #pragma omp critical (MagickCore_Classify)
643 proceed=SetImageProgress(image,SegmentImageTag,progress++,
645 if (proceed == MagickFalse)
649 image_view=DestroyCacheView(image_view);
650 status&=SyncImage(image);
652 Relinquish resources.
654 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
656 next_cluster=cluster->next;
657 cluster=(Cluster *) RelinquishMagickMemory(cluster);
660 free_squares=squares;
661 free_squares=(MagickRealType *) RelinquishMagickMemory(free_squares);
666 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
670 + C o n s o l i d a t e C r o s s i n g s %
674 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
676 % ConsolidateCrossings() guarantees that an even number of zero crossings
677 % always lie between two crossings.
679 % The format of the ConsolidateCrossings method is:
681 % ConsolidateCrossings(ZeroCrossing *zero_crossing,
682 % const size_t number_crossings)
684 % A description of each parameter follows.
686 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
688 % o number_crossings: This size_t specifies the number of elements
689 % in the zero_crossing array.
693 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
700 static inline ssize_t MagickMax(const ssize_t x,const ssize_t y)
707 static inline ssize_t MagickMin(const ssize_t x,const ssize_t y)
714 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
715 const size_t number_crossings)
731 Consolidate zero crossings.
733 for (i=(ssize_t) number_crossings-1; i >= 0; i--)
734 for (j=0; j <= 255; j++)
736 if (zero_crossing[i].crossings[j] == 0)
739 Find the entry that is closest to j and still preserves the
740 property that there are an even number of crossings between
743 for (k=j-1; k > 0; k--)
744 if (zero_crossing[i+1].crossings[k] != 0)
748 for (k=j+1; k < 255; k++)
749 if (zero_crossing[i+1].crossings[k] != 0)
751 right=MagickMin(k,255);
753 K is the zero crossing just left of j.
755 for (k=j-1; k > 0; k--)
756 if (zero_crossing[i].crossings[k] != 0)
761 Check center for an even number of crossings between k and j.
764 if (zero_crossing[i+1].crossings[j] != 0)
767 for (l=k+1; l < center; l++)
768 if (zero_crossing[i+1].crossings[l] != 0)
770 if (((count % 2) == 0) && (center != k))
774 Check left for an even number of crossings between k and j.
779 for (l=k+1; l < left; l++)
780 if (zero_crossing[i+1].crossings[l] != 0)
782 if (((count % 2) == 0) && (left != k))
786 Check right for an even number of crossings between k and j.
791 for (l=k+1; l < right; l++)
792 if (zero_crossing[i+1].crossings[l] != 0)
794 if (((count % 2) == 0) && (right != k))
797 l=(ssize_t) zero_crossing[i].crossings[j];
798 zero_crossing[i].crossings[j]=0;
800 zero_crossing[i].crossings[correct]=(short) l;
805 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
809 + D e f i n e R e g i o n %
813 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
815 % DefineRegion() defines the left and right boundaries of a peak region.
817 % The format of the DefineRegion method is:
819 % ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
821 % A description of each parameter follows.
823 % o extrema: Specifies a pointer to an array of integers. They
824 % represent the peaks and valleys of the histogram for each color
827 % o extents: This pointer to an ExtentPacket represent the extends
828 % of a particular peak or valley of a color component.
831 static ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)
834 Initialize to default values.
840 Find the left side (maxima).
842 for ( ; extents->index <= 255; extents->index++)
843 if (extrema[extents->index] > 0)
845 if (extents->index > 255)
846 return(MagickFalse); /* no left side - no region exists */
847 extents->left=extents->index;
849 Find the right side (minima).
851 for ( ; extents->index <= 255; extents->index++)
852 if (extrema[extents->index] < 0)
854 extents->right=extents->index-1;
859 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
863 + D e r i v a t i v e H i s t o g r a m %
867 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
869 % DerivativeHistogram() determines the derivative of the histogram using
870 % central differencing.
872 % The format of the DerivativeHistogram method is:
874 % DerivativeHistogram(const MagickRealType *histogram,
875 % MagickRealType *derivative)
877 % A description of each parameter follows.
879 % o histogram: Specifies an array of MagickRealTypes representing the number
880 % of pixels for each intensity of a particular color component.
882 % o derivative: This array of MagickRealTypes is initialized by
883 % DerivativeHistogram to the derivative of the histogram using central
887 static void DerivativeHistogram(const MagickRealType *histogram,
888 MagickRealType *derivative)
895 Compute endpoints using second order polynomial interpolation.
898 derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);
899 derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);
901 Compute derivative using central differencing.
903 for (i=1; i < n; i++)
904 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
909 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
913 + G e t I m a g e D y n a m i c T h r e s h o l d %
917 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
919 % GetImageDynamicThreshold() returns the dynamic threshold for an image.
921 % The format of the GetImageDynamicThreshold method is:
923 % MagickBooleanType GetImageDynamicThreshold(const Image *image,
924 % const double cluster_threshold,const double smooth_threshold,
925 % MagickPixelPacket *pixel,ExceptionInfo *exception)
927 % A description of each parameter follows.
929 % o image: the image.
931 % o cluster_threshold: This MagickRealType represents the minimum number of
932 % pixels contained in a hexahedra before it can be considered valid
933 % (expressed as a percentage).
935 % o smooth_threshold: the smoothing threshold eliminates noise in the second
936 % derivative of the histogram. As the value is increased, you can expect a
937 % smoother second derivative.
939 % o pixel: return the dynamic threshold here.
941 % o exception: return any errors or warnings in this structure.
944 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,
945 const double cluster_threshold,const double smooth_threshold,
946 MagickPixelPacket *pixel,ExceptionInfo *exception)
967 register const PixelPacket
975 *extrema[MaxDimension];
979 *histogram[MaxDimension],
983 Allocate histogram and extrema.
985 assert(image != (Image *) NULL);
986 assert(image->signature == MagickSignature);
987 if (image->debug != MagickFalse)
988 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
989 GetMagickPixelPacket(image,pixel);
990 for (i=0; i < MaxDimension; i++)
992 histogram[i]=(ssize_t *) AcquireQuantumMemory(256UL,sizeof(**histogram));
993 extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
994 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
996 for (i-- ; i >= 0; i--)
998 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
999 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1001 (void) ThrowMagickException(exception,GetMagickModule(),
1002 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1003 return(MagickFalse);
1007 Initialize histogram.
1009 InitializeHistogram(image,histogram,exception);
1010 (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
1011 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
1012 (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
1013 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
1014 (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
1015 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
1019 cluster=(Cluster *) NULL;
1020 head=(Cluster *) NULL;
1021 (void) ResetMagickMemory(&red,0,sizeof(red));
1022 (void) ResetMagickMemory(&green,0,sizeof(green));
1023 (void) ResetMagickMemory(&blue,0,sizeof(blue));
1024 while (DefineRegion(extrema[Red],&red) != 0)
1027 while (DefineRegion(extrema[Green],&green) != 0)
1030 while (DefineRegion(extrema[Blue],&blue) != 0)
1033 Allocate a new class.
1035 if (head != (Cluster *) NULL)
1037 cluster->next=(Cluster *) AcquireMagickMemory(
1038 sizeof(*cluster->next));
1039 cluster=cluster->next;
1043 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1046 if (cluster == (Cluster *) NULL)
1048 (void) ThrowMagickException(exception,GetMagickModule(),
1049 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1051 return(MagickFalse);
1054 Initialize a new class.
1058 cluster->green=green;
1060 cluster->next=(Cluster *) NULL;
1064 if (head == (Cluster *) NULL)
1067 No classes were identified-- create one.
1069 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1070 if (cluster == (Cluster *) NULL)
1072 (void) ThrowMagickException(exception,GetMagickModule(),
1073 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1074 return(MagickFalse);
1077 Initialize a new class.
1081 cluster->green=green;
1083 cluster->next=(Cluster *) NULL;
1087 Count the pixels for each cluster.
1090 for (y=0; y < (ssize_t) image->rows; y++)
1092 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1093 if (p == (const PixelPacket *) NULL)
1095 for (x=0; x < (ssize_t) image->columns; x++)
1097 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1098 if (((ssize_t) ScaleQuantumToChar(GetRedPixelComponent(p)) >=
1099 (cluster->red.left-SafeMargin)) &&
1100 ((ssize_t) ScaleQuantumToChar(GetRedPixelComponent(p)) <=
1101 (cluster->red.right+SafeMargin)) &&
1102 ((ssize_t) ScaleQuantumToChar(GetGreenPixelComponent(p)) >=
1103 (cluster->green.left-SafeMargin)) &&
1104 ((ssize_t) ScaleQuantumToChar(GetGreenPixelComponent(p)) <=
1105 (cluster->green.right+SafeMargin)) &&
1106 ((ssize_t) ScaleQuantumToChar(GetBluePixelComponent(p)) >=
1107 (cluster->blue.left-SafeMargin)) &&
1108 ((ssize_t) ScaleQuantumToChar(GetBluePixelComponent(p)) <=
1109 (cluster->blue.right+SafeMargin)))
1115 cluster->red.center+=(MagickRealType)
1116 ScaleQuantumToChar(GetRedPixelComponent(p));
1117 cluster->green.center+=(MagickRealType)
1118 ScaleQuantumToChar(GetGreenPixelComponent(p));
1119 cluster->blue.center+=(MagickRealType)
1120 ScaleQuantumToChar(GetBluePixelComponent(p));
1126 proceed=SetImageProgress(image,SegmentImageTag,(MagickOffsetType) y,
1128 if (proceed == MagickFalse)
1132 Remove clusters that do not meet minimum cluster threshold.
1137 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1139 next_cluster=cluster->next;
1140 if ((cluster->count > 0) &&
1141 (cluster->count >= (count*cluster_threshold/100.0)))
1147 cluster->red.center/=cluster->count;
1148 cluster->green.center/=cluster->count;
1149 cluster->blue.center/=cluster->count;
1151 last_cluster=cluster;
1157 if (cluster == head)
1160 last_cluster->next=next_cluster;
1161 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1168 for (cluster=object; cluster->next != (Cluster *) NULL; )
1170 if (cluster->count < object->count)
1172 cluster=cluster->next;
1174 background=head->next;
1175 for (cluster=background; cluster->next != (Cluster *) NULL; )
1177 if (cluster->count > background->count)
1179 cluster=cluster->next;
1182 threshold=(background->red.center+object->red.center)/2.0;
1183 pixel->red=(MagickRealType) ScaleCharToQuantum((unsigned char)
1185 threshold=(background->green.center+object->green.center)/2.0;
1186 pixel->green=(MagickRealType) ScaleCharToQuantum((unsigned char)
1188 threshold=(background->blue.center+object->blue.center)/2.0;
1189 pixel->blue=(MagickRealType) ScaleCharToQuantum((unsigned char)
1192 Relinquish resources.
1194 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1196 next_cluster=cluster->next;
1197 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1199 for (i=0; i < MaxDimension; i++)
1201 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1202 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1208 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1212 + I n i t i a l i z e H i s t o g r a m %
1216 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1218 % InitializeHistogram() computes the histogram for an image.
1220 % The format of the InitializeHistogram method is:
1222 % InitializeHistogram(const Image *image,ssize_t **histogram)
1224 % A description of each parameter follows.
1226 % o image: Specifies a pointer to an Image structure; returned from
1229 % o histogram: Specifies an array of integers representing the number
1230 % of pixels for each intensity of a particular color component.
1233 static void InitializeHistogram(const Image *image,ssize_t **histogram,
1234 ExceptionInfo *exception)
1236 register const PixelPacket
1247 Initialize histogram.
1249 for (i=0; i <= 255; i++)
1251 histogram[Red][i]=0;
1252 histogram[Green][i]=0;
1253 histogram[Blue][i]=0;
1255 for (y=0; y < (ssize_t) image->rows; y++)
1257 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1258 if (p == (const PixelPacket *) NULL)
1260 for (x=0; x < (ssize_t) image->columns; x++)
1262 histogram[Red][(ssize_t) ScaleQuantumToChar(GetRedPixelComponent(p))]++;
1263 histogram[Green][(ssize_t) ScaleQuantumToChar(GetGreenPixelComponent(p))]++;
1264 histogram[Blue][(ssize_t) ScaleQuantumToChar(GetBluePixelComponent(p))]++;
1271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1275 + I n i t i a l i z e I n t e r v a l T r e e %
1279 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1281 % InitializeIntervalTree() initializes an interval tree from the lists of
1284 % The format of the InitializeIntervalTree method is:
1286 % InitializeIntervalTree(IntervalTree **list,ssize_t *number_nodes,
1287 % IntervalTree *node)
1289 % A description of each parameter follows.
1291 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
1293 % o number_crossings: This size_t specifies the number of elements
1294 % in the zero_crossing array.
1298 static void InitializeList(IntervalTree **list,ssize_t *number_nodes,
1301 if (node == (IntervalTree *) NULL)
1303 if (node->child == (IntervalTree *) NULL)
1304 list[(*number_nodes)++]=node;
1305 InitializeList(list,number_nodes,node->sibling);
1306 InitializeList(list,number_nodes,node->child);
1309 static void MeanStability(IntervalTree *node)
1311 register IntervalTree
1314 if (node == (IntervalTree *) NULL)
1316 node->mean_stability=0.0;
1318 if (child != (IntervalTree *) NULL)
1323 register MagickRealType
1328 for ( ; child != (IntervalTree *) NULL; child=child->sibling)
1330 sum+=child->stability;
1333 node->mean_stability=sum/(MagickRealType) count;
1335 MeanStability(node->sibling);
1336 MeanStability(node->child);
1339 static void Stability(IntervalTree *node)
1341 if (node == (IntervalTree *) NULL)
1343 if (node->child == (IntervalTree *) NULL)
1344 node->stability=0.0;
1346 node->stability=node->tau-(node->child)->tau;
1347 Stability(node->sibling);
1348 Stability(node->child);
1351 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing,
1352 const size_t number_crossings)
1370 Allocate interval tree.
1372 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1374 if (list == (IntervalTree **) NULL)
1375 return((IntervalTree *) NULL);
1377 The root is the entire histogram.
1379 root=(IntervalTree *) AcquireMagickMemory(sizeof(*root));
1380 root->child=(IntervalTree *) NULL;
1381 root->sibling=(IntervalTree *) NULL;
1385 for (i=(-1); i < (ssize_t) number_crossings; i++)
1388 Initialize list with all nodes with no children.
1391 InitializeList(list,&number_nodes,root);
1395 for (j=0; j < number_nodes; j++)
1400 for (k=head->left+1; k < head->right; k++)
1402 if (zero_crossing[i+1].crossings[k] != 0)
1406 node->child=(IntervalTree *) AcquireMagickMemory(
1407 sizeof(*node->child));
1412 node->sibling=(IntervalTree *) AcquireMagickMemory(
1413 sizeof(*node->sibling));
1416 node->tau=zero_crossing[i+1].tau;
1417 node->child=(IntervalTree *) NULL;
1418 node->sibling=(IntervalTree *) NULL;
1424 if (left != head->left)
1426 node->sibling=(IntervalTree *) AcquireMagickMemory(
1427 sizeof(*node->sibling));
1429 node->tau=zero_crossing[i+1].tau;
1430 node->child=(IntervalTree *) NULL;
1431 node->sibling=(IntervalTree *) NULL;
1433 node->right=head->right;
1438 Determine the stability: difference between a nodes tau and its child.
1440 Stability(root->child);
1441 MeanStability(root->child);
1442 list=(IntervalTree **) RelinquishMagickMemory(list);
1447 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1451 + O p t i m a l T a u %
1455 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1457 % OptimalTau() finds the optimal tau for each band of the histogram.
1459 % The format of the OptimalTau method is:
1461 % MagickRealType OptimalTau(const ssize_t *histogram,const double max_tau,
1462 % const double min_tau,const double delta_tau,
1463 % const double smooth_threshold,short *extrema)
1465 % A description of each parameter follows.
1467 % o histogram: Specifies an array of integers representing the number
1468 % of pixels for each intensity of a particular color component.
1470 % o extrema: Specifies a pointer to an array of integers. They
1471 % represent the peaks and valleys of the histogram for each color
1476 static void ActiveNodes(IntervalTree **list,ssize_t *number_nodes,
1479 if (node == (IntervalTree *) NULL)
1481 if (node->stability >= node->mean_stability)
1483 list[(*number_nodes)++]=node;
1484 ActiveNodes(list,number_nodes,node->sibling);
1488 ActiveNodes(list,number_nodes,node->sibling);
1489 ActiveNodes(list,number_nodes,node->child);
1493 static void FreeNodes(IntervalTree *node)
1495 if (node == (IntervalTree *) NULL)
1497 FreeNodes(node->sibling);
1498 FreeNodes(node->child);
1499 node=(IntervalTree *) RelinquishMagickMemory(node);
1502 static MagickRealType OptimalTau(const ssize_t *histogram,const double max_tau,
1503 const double min_tau,const double delta_tau,const double smooth_threshold,
1539 Allocate interval tree.
1541 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1543 if (list == (IntervalTree **) NULL)
1546 Allocate zero crossing list.
1548 count=(size_t) ((max_tau-min_tau)/delta_tau)+2;
1549 zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count,
1550 sizeof(*zero_crossing));
1551 if (zero_crossing == (ZeroCrossing *) NULL)
1553 for (i=0; i < (ssize_t) count; i++)
1554 zero_crossing[i].tau=(-1.0);
1556 Initialize zero crossing list.
1558 derivative=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*derivative));
1559 second_derivative=(MagickRealType *) AcquireQuantumMemory(256,
1560 sizeof(*second_derivative));
1561 if ((derivative == (MagickRealType *) NULL) ||
1562 (second_derivative == (MagickRealType *) NULL))
1563 ThrowFatalException(ResourceLimitFatalError,
1564 "UnableToAllocateDerivatives");
1566 for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
1568 zero_crossing[i].tau=tau;
1569 ScaleSpace(histogram,tau,zero_crossing[i].histogram);
1570 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1571 DerivativeHistogram(derivative,second_derivative);
1572 ZeroCrossHistogram(second_derivative,smooth_threshold,
1573 zero_crossing[i].crossings);
1577 Add an entry for the original histogram.
1579 zero_crossing[i].tau=0.0;
1580 for (j=0; j <= 255; j++)
1581 zero_crossing[i].histogram[j]=(MagickRealType) histogram[j];
1582 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1583 DerivativeHistogram(derivative,second_derivative);
1584 ZeroCrossHistogram(second_derivative,smooth_threshold,
1585 zero_crossing[i].crossings);
1586 number_crossings=(size_t) i;
1587 derivative=(MagickRealType *) RelinquishMagickMemory(derivative);
1588 second_derivative=(MagickRealType *)
1589 RelinquishMagickMemory(second_derivative);
1591 Ensure the scale-space fingerprints form lines in scale-space, not loops.
1593 ConsolidateCrossings(zero_crossing,number_crossings);
1595 Force endpoints to be included in the interval.
1597 for (i=0; i <= (ssize_t) number_crossings; i++)
1599 for (j=0; j < 255; j++)
1600 if (zero_crossing[i].crossings[j] != 0)
1602 zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]);
1603 for (j=255; j > 0; j--)
1604 if (zero_crossing[i].crossings[j] != 0)
1606 zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1609 Initialize interval tree.
1611 root=InitializeIntervalTree(zero_crossing,number_crossings);
1612 if (root == (IntervalTree *) NULL)
1615 Find active nodes: stability is greater (or equal) to the mean stability of
1619 ActiveNodes(list,&number_nodes,root->child);
1623 for (i=0; i <= 255; i++)
1625 for (i=0; i < number_nodes; i++)
1628 Find this tau in zero crossings list.
1632 for (j=0; j <= (ssize_t) number_crossings; j++)
1633 if (zero_crossing[j].tau == node->tau)
1636 Find the value of the peak.
1638 peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1641 value=zero_crossing[k].histogram[index];
1642 for (x=node->left; x <= node->right; x++)
1644 if (peak != MagickFalse)
1646 if (zero_crossing[k].histogram[x] > value)
1648 value=zero_crossing[k].histogram[x];
1653 if (zero_crossing[k].histogram[x] < value)
1655 value=zero_crossing[k].histogram[x];
1659 for (x=node->left; x <= node->right; x++)
1663 if (peak != MagickFalse)
1664 extrema[x]=(short) index;
1666 extrema[x]=(short) (-index);
1670 Determine the average tau.
1673 for (i=0; i < number_nodes; i++)
1674 average_tau+=list[i]->tau;
1675 average_tau/=(MagickRealType) number_nodes;
1677 Relinquish resources.
1680 zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1681 list=(IntervalTree **) RelinquishMagickMemory(list);
1682 return(average_tau);
1686 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1690 + S c a l e S p a c e %
1694 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1696 % ScaleSpace() performs a scale-space filter on the 1D histogram.
1698 % The format of the ScaleSpace method is:
1700 % ScaleSpace(const ssize_t *histogram,const MagickRealType tau,
1701 % MagickRealType *scale_histogram)
1703 % A description of each parameter follows.
1705 % o histogram: Specifies an array of MagickRealTypes representing the number
1706 % of pixels for each intensity of a particular color component.
1710 static void ScaleSpace(const ssize_t *histogram,const MagickRealType tau,
1711 MagickRealType *scale_histogram)
1723 gamma=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*gamma));
1724 if (gamma == (MagickRealType *) NULL)
1725 ThrowFatalException(ResourceLimitFatalError,
1726 "UnableToAllocateGammaMap");
1727 alpha=1.0/(tau*sqrt(2.0*MagickPI));
1728 beta=(-1.0/(2.0*tau*tau));
1729 for (x=0; x <= 255; x++)
1731 for (x=0; x <= 255; x++)
1733 gamma[x]=exp((double) beta*x*x);
1734 if (gamma[x] < MagickEpsilon)
1737 for (x=0; x <= 255; x++)
1740 for (u=0; u <= 255; u++)
1741 sum+=(MagickRealType) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1742 scale_histogram[x]=alpha*sum;
1744 gamma=(MagickRealType *) RelinquishMagickMemory(gamma);
1748 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1752 % S e g m e n t I m a g e %
1756 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1758 % SegmentImage() segment an image by analyzing the histograms of the color
1759 % components and identifying units that are homogeneous with the fuzzy
1760 % C-means technique.
1762 % The format of the SegmentImage method is:
1764 % MagickBooleanType SegmentImage(Image *image,
1765 % const ColorspaceType colorspace,const MagickBooleanType verbose,
1766 % const double cluster_threshold,const double smooth_threshold)
1768 % A description of each parameter follows.
1770 % o image: the image.
1772 % o colorspace: Indicate the colorspace.
1774 % o verbose: Set to MagickTrue to print detailed information about the
1775 % identified classes.
1777 % o cluster_threshold: This represents the minimum number of pixels
1778 % contained in a hexahedra before it can be considered valid (expressed
1781 % o smooth_threshold: the smoothing threshold eliminates noise in the second
1782 % derivative of the histogram. As the value is increased, you can expect a
1783 % smoother second derivative.
1786 MagickExport MagickBooleanType SegmentImage(Image *image,
1787 const ColorspaceType colorspace,const MagickBooleanType verbose,
1788 const double cluster_threshold,const double smooth_threshold)
1797 *extrema[MaxDimension];
1800 *histogram[MaxDimension];
1803 Allocate histogram and extrema.
1805 assert(image != (Image *) NULL);
1806 assert(image->signature == MagickSignature);
1807 if (image->debug != MagickFalse)
1808 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1809 for (i=0; i < MaxDimension; i++)
1811 histogram[i]=(ssize_t *) AcquireQuantumMemory(256,sizeof(**histogram));
1812 extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1813 if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
1815 for (i-- ; i >= 0; i--)
1817 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1818 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1820 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1824 if (colorspace != RGBColorspace)
1825 (void) TransformImageColorspace(image,colorspace);
1827 Initialize histogram.
1829 InitializeHistogram(image,histogram,&image->exception);
1830 (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,
1831 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Red]);
1832 (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,
1833 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Green]);
1834 (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,
1835 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Blue]);
1837 Classify using the fuzzy c-Means technique.
1839 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose);
1840 if (colorspace != RGBColorspace)
1841 (void) TransformImageColorspace(image,colorspace);
1843 Relinquish resources.
1845 for (i=0; i < MaxDimension; i++)
1847 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1848 histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
1854 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1858 + Z e r o C r o s s H i s t o g r a m %
1862 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1864 % ZeroCrossHistogram() find the zero crossings in a histogram and marks
1865 % directions as: 1 is negative to positive; 0 is zero crossing; and -1
1866 % is positive to negative.
1868 % The format of the ZeroCrossHistogram method is:
1870 % ZeroCrossHistogram(MagickRealType *second_derivative,
1871 % const MagickRealType smooth_threshold,short *crossings)
1873 % A description of each parameter follows.
1875 % o second_derivative: Specifies an array of MagickRealTypes representing the
1876 % second derivative of the histogram of a particular color component.
1878 % o crossings: This array of integers is initialized with
1879 % -1, 0, or 1 representing the slope of the first derivative of the
1880 % of a particular color component.
1883 static void ZeroCrossHistogram(MagickRealType *second_derivative,
1884 const MagickRealType smooth_threshold,short *crossings)
1893 Merge low numbers to zero to help prevent noise.
1895 for (i=0; i <= 255; i++)
1896 if ((second_derivative[i] < smooth_threshold) &&
1897 (second_derivative[i] >= -smooth_threshold))
1898 second_derivative[i]=0.0;
1900 Mark zero crossings.
1903 for (i=0; i <= 255; i++)
1906 if (second_derivative[i] < 0.0)
1913 if (second_derivative[i] > 0.0)