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-2009 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"
272 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 *) AcquireMagickMemory(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 *) AcquireMagickMemory(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(p->red) >=
367 (cluster->red.left-SafeMargin)) &&
368 ((long) ScaleQuantumToChar(p->red) <=
369 (cluster->red.right+SafeMargin)) &&
370 ((long) ScaleQuantumToChar(p->green) >=
371 (cluster->green.left-SafeMargin)) &&
372 ((long) ScaleQuantumToChar(p->green) <=
373 (cluster->green.right+SafeMargin)) &&
374 ((long) ScaleQuantumToChar(p->blue) >=
375 (cluster->blue.left-SafeMargin)) &&
376 ((long) ScaleQuantumToChar(p->blue) <=
377 (cluster->blue.right+SafeMargin)))
383 cluster->red.center+=(MagickRealType) ScaleQuantumToChar(p->red);
384 cluster->green.center+=(MagickRealType)
385 ScaleQuantumToChar(p->green);
386 cluster->blue.center+=(MagickRealType) ScaleQuantumToChar(p->blue);
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 = %g\n",cluster_threshold);
448 (void) fprintf(stdout,"\tWeighting Exponent = %g\n",weighting_exponent);
449 (void) fprintf(stdout,"\tTotal Number of Clusters = %lu\n\n",
452 Print the total number of points per cluster.
454 (void) fprintf(stdout,"\n\nNumber of Vectors Per Cluster\n");
455 (void) fprintf(stdout,"=============================\n\n");
456 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
457 (void) fprintf(stdout,"Cluster #%ld = %ld\n",cluster->id,
460 Print the cluster extents.
462 (void) fprintf(stdout,
463 "\n\n\nCluster Extents: (Vector Size: %d)\n",MaxDimension);
464 (void) fprintf(stdout,"================");
465 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
467 (void) fprintf(stdout,"\n\nCluster #%ld\n\n",cluster->id);
468 (void) fprintf(stdout,"%ld-%ld %ld-%ld %ld-%ld\n",cluster->red.left,
469 cluster->red.right,cluster->green.left,cluster->green.right,
470 cluster->blue.left,cluster->blue.right);
473 Print the cluster center values.
475 (void) fprintf(stdout,
476 "\n\n\nCluster Center Values: (Vector Size: %d)\n",MaxDimension);
477 (void) fprintf(stdout,"=====================");
478 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
480 (void) fprintf(stdout,"\n\nCluster #%ld\n\n",cluster->id);
481 (void) fprintf(stdout,"%g %g %g\n",(double) cluster->red.center,
482 (double) cluster->green.center,(double) cluster->blue.center);
484 (void) fprintf(stdout,"\n");
486 if (number_clusters > 256)
487 ThrowBinaryException(ImageError,"TooManyClusters",image->filename);
489 Speed up distance calculations.
491 squares=(MagickRealType *) AcquireQuantumMemory(513UL,sizeof(*squares));
492 if (squares == (MagickRealType *) NULL)
493 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
496 for (i=(-255); i <= 255; i++)
497 squares[i]=(MagickRealType) i*(MagickRealType) i;
499 Allocate image colormap.
501 if (AcquireImageColormap(image,number_clusters) == MagickFalse)
502 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
505 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
507 image->colormap[i].red=ScaleCharToQuantum((unsigned char)
508 (cluster->red.center+0.5));
509 image->colormap[i].green=ScaleCharToQuantum((unsigned char)
510 (cluster->green.center+0.5));
511 image->colormap[i].blue=ScaleCharToQuantum((unsigned char)
512 (cluster->blue.center+0.5));
516 Do course grain classes.
518 exception=(&image->exception);
519 image_view=AcquireCacheView(image);
520 #if defined(MAGICKCORE_OPENMP_SUPPORT)
521 #pragma omp parallel for schedule(dynamic,8) shared(progress,status)
523 for (y=0; y < (long) image->rows; y++)
528 register const PixelPacket
540 if (status == MagickFalse)
542 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
543 if (q == (PixelPacket *) NULL)
548 indexes=GetCacheViewAuthenticIndexQueue(image_view);
549 for (x=0; x < (long) image->columns; x++)
551 indexes[x]=(IndexPacket) 0;
552 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
554 if (((long) ScaleQuantumToChar(q->red) >=
555 (cluster->red.left-SafeMargin)) &&
556 ((long) ScaleQuantumToChar(q->red) <=
557 (cluster->red.right+SafeMargin)) &&
558 ((long) ScaleQuantumToChar(q->green) >=
559 (cluster->green.left-SafeMargin)) &&
560 ((long) ScaleQuantumToChar(q->green) <=
561 (cluster->green.right+SafeMargin)) &&
562 ((long) ScaleQuantumToChar(q->blue) >=
563 (cluster->blue.left-SafeMargin)) &&
564 ((long) ScaleQuantumToChar(q->blue) <=
565 (cluster->blue.right+SafeMargin)))
570 indexes[x]=(IndexPacket) cluster->id;
574 if (cluster == (Cluster *) NULL)
588 Compute fuzzy membership.
591 for (j=0; j < (long) image->colors; j++)
595 distance_squared=squares[(long) ScaleQuantumToChar(q->red)-
596 (long) ScaleQuantumToChar(p->red)]+
597 squares[(long) ScaleQuantumToChar(q->green)-
598 (long) ScaleQuantumToChar(p->green)]+
599 squares[(long) ScaleQuantumToChar(q->blue)-
600 (long) ScaleQuantumToChar(p->blue)];
601 numerator=distance_squared;
602 for (k=0; k < (long) image->colors; k++)
605 distance_squared=squares[(long) ScaleQuantumToChar(q->red)-
606 (long) ScaleQuantumToChar(p->red)]+
607 squares[(long) ScaleQuantumToChar(q->green)-
608 (long) ScaleQuantumToChar(p->green)]+
609 squares[(long) ScaleQuantumToChar(q->blue)-
610 (long) ScaleQuantumToChar(p->blue)];
611 ratio=numerator/distance_squared;
612 sum+=SegmentPower(ratio);
614 if ((sum != 0.0) && ((1.0/sum) > local_minima))
619 local_minima=1.0/sum;
620 indexes[x]=(IndexPacket) j;
626 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
628 if (image->progress_monitor != (MagickProgressMonitor) NULL)
633 #if defined(MAGICKCORE_OPENMP_SUPPORT)
634 #pragma omp critical (MagickCore_Classify)
636 proceed=SetImageProgress(image,SegmentImageTag,progress++,
638 if (proceed == MagickFalse)
642 image_view=DestroyCacheView(image_view);
643 status&=SyncImage(image);
645 Relinquish resources.
647 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
649 next_cluster=cluster->next;
650 cluster=(Cluster *) RelinquishMagickMemory(cluster);
653 free_squares=squares;
654 free_squares=(MagickRealType *) RelinquishMagickMemory(free_squares);
659 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
663 + C o n s o l i d a t e C r o s s i n g s %
667 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
669 % ConsolidateCrossings() guarantees that an even number of zero crossings
670 % always lie between two crossings.
672 % The format of the ConsolidateCrossings method is:
674 % ConsolidateCrossings(ZeroCrossing *zero_crossing,
675 % const unsigned long number_crossings)
677 % A description of each parameter follows.
679 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
681 % o number_crossings: This unsigned long specifies the number of elements
682 % in the zero_crossing array.
686 static inline long MagickAbsoluteValue(const long x)
693 static inline long MagickMax(const long x,const long y)
700 static inline long MagickMin(const long x,const long y)
707 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
708 const unsigned long number_crossings)
724 Consolidate zero crossings.
726 for (i=(long) number_crossings-1; i >= 0; i--)
727 for (j=0; j <= 255; j++)
729 if (zero_crossing[i].crossings[j] == 0)
732 Find the entry that is closest to j and still preserves the
733 property that there are an even number of crossings between
736 for (k=j-1; k > 0; k--)
737 if (zero_crossing[i+1].crossings[k] != 0)
741 for (k=j+1; k < 255; k++)
742 if (zero_crossing[i+1].crossings[k] != 0)
744 right=MagickMin(k,255);
746 K is the zero crossing just left of j.
748 for (k=j-1; k > 0; k--)
749 if (zero_crossing[i].crossings[k] != 0)
754 Check center for an even number of crossings between k and j.
757 if (zero_crossing[i+1].crossings[j] != 0)
760 for (l=k+1; l < center; l++)
761 if (zero_crossing[i+1].crossings[l] != 0)
763 if (((count % 2) == 0) && (center != k))
767 Check left for an even number of crossings between k and j.
772 for (l=k+1; l < left; l++)
773 if (zero_crossing[i+1].crossings[l] != 0)
775 if (((count % 2) == 0) && (left != k))
779 Check right for an even number of crossings between k and j.
784 for (l=k+1; l < right; l++)
785 if (zero_crossing[i+1].crossings[l] != 0)
787 if (((count % 2) == 0) && (right != k))
790 l=zero_crossing[i].crossings[j];
791 zero_crossing[i].crossings[j]=0;
793 zero_crossing[i].crossings[correct]=(short) l;
798 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
802 + D e f i n e R e g i o n %
806 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
808 % DefineRegion() defines the left and right boundaries of a peak region.
810 % The format of the DefineRegion method is:
812 % long DefineRegion(const short *extrema,ExtentPacket *extents)
814 % A description of each parameter follows.
816 % o extrema: Specifies a pointer to an array of integers. They
817 % represent the peaks and valleys of the histogram for each color
820 % o extents: This pointer to an ExtentPacket represent the extends
821 % of a particular peak or valley of a color component.
824 static long DefineRegion(const short *extrema,ExtentPacket *extents)
827 Initialize to default values.
833 Find the left side (maxima).
835 for ( ; extents->index <= 255; extents->index++)
836 if (extrema[extents->index] > 0)
838 if (extents->index > 255)
839 return(MagickFalse); /* no left side - no region exists */
840 extents->left=extents->index;
842 Find the right side (minima).
844 for ( ; extents->index <= 255; extents->index++)
845 if (extrema[extents->index] < 0)
847 extents->right=extents->index-1;
852 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
856 + D e r i v a t i v e H i s t o g r a m %
860 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
862 % DerivativeHistogram() determines the derivative of the histogram using
863 % central differencing.
865 % The format of the DerivativeHistogram method is:
867 % DerivativeHistogram(const MagickRealType *histogram,
868 % MagickRealType *derivative)
870 % A description of each parameter follows.
872 % o histogram: Specifies an array of MagickRealTypes representing the number
873 % of pixels for each intensity of a particular color component.
875 % o derivative: This array of MagickRealTypes is initialized by
876 % DerivativeHistogram to the derivative of the histogram using central
880 static void DerivativeHistogram(const MagickRealType *histogram,
881 MagickRealType *derivative)
888 Compute endpoints using second order polynomial interpolation.
891 derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);
892 derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);
894 Compute derivative using central differencing.
896 for (i=1; i < n; i++)
897 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
902 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
906 + G e t I m a g e D y n a m i c T h r e s h o l d %
910 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
912 % GetImageDynamicThreshold() returns the dynamic threshold for an image.
914 % The format of the GetImageDynamicThreshold method is:
916 % MagickBooleanType GetImageDynamicThreshold(const Image *image,
917 % const double cluster_threshold,const double smooth_threshold,
918 % MagickPixelPacket *pixel,ExceptionInfo *exception)
920 % A description of each parameter follows.
922 % o image: the image.
924 % o cluster_threshold: This MagickRealType represents the minimum number of
925 % pixels contained in a hexahedra before it can be considered valid
926 % (expressed as a percentage).
928 % o smooth_threshold: the smoothing threshold eliminates noise in the second
929 % derivative of the histogram. As the value is increased, you can expect a
930 % smoother second derivative.
932 % o pixel: return the dynamic threshold here.
934 % o exception: return any errors or warnings in this structure.
937 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,
938 const double cluster_threshold,const double smooth_threshold,
939 MagickPixelPacket *pixel,ExceptionInfo *exception)
956 *histogram[MaxDimension],
965 register const PixelPacket
973 *extrema[MaxDimension];
976 Allocate histogram and extrema.
978 assert(image != (Image *) NULL);
979 assert(image->signature == MagickSignature);
980 if (image->debug != MagickFalse)
981 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
982 GetMagickPixelPacket(image,pixel);
983 for (i=0; i < MaxDimension; i++)
985 histogram[i]=(long *) AcquireQuantumMemory(256UL,sizeof(**histogram));
986 extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
987 if ((histogram[i] == (long *) NULL) || (extrema[i] == (short *) NULL))
989 for (i-- ; i >= 0; i--)
991 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
992 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
994 (void) ThrowMagickException(exception,GetMagickModule(),
995 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1000 Initialize histogram.
1002 InitializeHistogram(image,histogram,exception);
1003 (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
1004 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
1005 (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
1006 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
1007 (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
1008 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
1012 cluster=(Cluster *) NULL;
1013 head=(Cluster *) NULL;
1014 (void) ResetMagickMemory(&red,0,sizeof(red));
1015 (void) ResetMagickMemory(&green,0,sizeof(green));
1016 (void) ResetMagickMemory(&blue,0,sizeof(blue));
1017 while (DefineRegion(extrema[Red],&red) != 0)
1020 while (DefineRegion(extrema[Green],&green) != 0)
1023 while (DefineRegion(extrema[Blue],&blue) != 0)
1026 Allocate a new class.
1028 if (head != (Cluster *) NULL)
1030 cluster->next=(Cluster *) AcquireMagickMemory(
1031 sizeof(*cluster->next));
1032 cluster=cluster->next;
1036 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1039 if (cluster == (Cluster *) NULL)
1041 (void) ThrowMagickException(exception,GetMagickModule(),
1042 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1044 return(MagickFalse);
1047 Initialize a new class.
1051 cluster->green=green;
1053 cluster->next=(Cluster *) NULL;
1057 if (head == (Cluster *) NULL)
1060 No classes were identified-- create one.
1062 cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
1063 if (cluster == (Cluster *) NULL)
1065 (void) ThrowMagickException(exception,GetMagickModule(),
1066 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1067 return(MagickFalse);
1070 Initialize a new class.
1074 cluster->green=green;
1076 cluster->next=(Cluster *) NULL;
1080 Count the pixels for each cluster.
1083 for (y=0; y < (long) image->rows; y++)
1085 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1086 if (p == (const PixelPacket *) NULL)
1088 for (x=0; x < (long) image->columns; x++)
1090 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1091 if (((long) ScaleQuantumToChar(p->red) >=
1092 (cluster->red.left-SafeMargin)) &&
1093 ((long) ScaleQuantumToChar(p->red) <=
1094 (cluster->red.right+SafeMargin)) &&
1095 ((long) ScaleQuantumToChar(p->green) >=
1096 (cluster->green.left-SafeMargin)) &&
1097 ((long) ScaleQuantumToChar(p->green) <=
1098 (cluster->green.right+SafeMargin)) &&
1099 ((long) ScaleQuantumToChar(p->blue) >=
1100 (cluster->blue.left-SafeMargin)) &&
1101 ((long) ScaleQuantumToChar(p->blue) <=
1102 (cluster->blue.right+SafeMargin)))
1108 cluster->red.center+=(MagickRealType)
1109 ScaleQuantumToChar(p->red);
1110 cluster->green.center+=(MagickRealType)
1111 ScaleQuantumToChar(p->green);
1112 cluster->blue.center+=(MagickRealType)
1113 ScaleQuantumToChar(p->blue);
1119 proceed=SetImageProgress(image,SegmentImageTag,y,2*image->rows);
1120 if (proceed == MagickFalse)
1124 Remove clusters that do not meet minimum cluster threshold.
1129 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1131 next_cluster=cluster->next;
1132 if ((cluster->count > 0) &&
1133 (cluster->count >= (count*cluster_threshold/100.0)))
1139 cluster->red.center/=cluster->count;
1140 cluster->green.center/=cluster->count;
1141 cluster->blue.center/=cluster->count;
1143 last_cluster=cluster;
1149 if (cluster == head)
1152 last_cluster->next=next_cluster;
1153 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1160 for (cluster=object; cluster->next != (Cluster *) NULL; )
1162 if (cluster->count < object->count)
1164 cluster=cluster->next;
1166 background=head->next;
1167 for (cluster=background; cluster->next != (Cluster *) NULL; )
1169 if (cluster->count > background->count)
1171 cluster=cluster->next;
1174 threshold=(background->red.center+object->red.center)/2.0;
1175 pixel->red=(MagickRealType) ScaleCharToQuantum((unsigned char)
1177 threshold=(background->green.center+object->green.center)/2.0;
1178 pixel->green=(MagickRealType) ScaleCharToQuantum((unsigned char)
1180 threshold=(background->blue.center+object->blue.center)/2.0;
1181 pixel->blue=(MagickRealType) ScaleCharToQuantum((unsigned char)
1184 Relinquish resources.
1186 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1188 next_cluster=cluster->next;
1189 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1191 for (i=0; i < MaxDimension; i++)
1193 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1194 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1204 + I n i t i a l i z e H i s t o g r a m %
1208 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1210 % InitializeHistogram() computes the histogram for an image.
1212 % The format of the InitializeHistogram method is:
1214 % InitializeHistogram(const Image *image,long **histogram)
1216 % A description of each parameter follows.
1218 % o image: Specifies a pointer to an Image structure; returned from
1221 % o histogram: Specifies an array of integers representing the number
1222 % of pixels for each intensity of a particular color component.
1225 static void InitializeHistogram(const Image *image,long **histogram,
1226 ExceptionInfo *exception)
1231 register const PixelPacket
1239 Initialize histogram.
1241 for (i=0; i <= 255; i++)
1243 histogram[Red][i]=0;
1244 histogram[Green][i]=0;
1245 histogram[Blue][i]=0;
1247 for (y=0; y < (long) image->rows; y++)
1249 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1250 if (p == (const PixelPacket *) NULL)
1252 for (x=0; x < (long) image->columns; x++)
1254 histogram[Red][(long) ScaleQuantumToChar(p->red)]++;
1255 histogram[Green][(long) ScaleQuantumToChar(p->green)]++;
1256 histogram[Blue][(long) ScaleQuantumToChar(p->blue)]++;
1263 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1267 + I n i t i a l i z e I n t e r v a l T r e e %
1271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1273 % InitializeIntervalTree() initializes an interval tree from the lists of
1276 % The format of the InitializeIntervalTree method is:
1278 % InitializeIntervalTree(IntervalTree **list,long *number_nodes,
1279 % IntervalTree *node)
1281 % A description of each parameter follows.
1283 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
1285 % o number_crossings: This unsigned long specifies the number of elements
1286 % in the zero_crossing array.
1290 static void InitializeList(IntervalTree **list,long *number_nodes,
1293 if (node == (IntervalTree *) NULL)
1295 if (node->child == (IntervalTree *) NULL)
1296 list[(*number_nodes)++]=node;
1297 InitializeList(list,number_nodes,node->sibling);
1298 InitializeList(list,number_nodes,node->child);
1301 static void MeanStability(IntervalTree *node)
1303 register IntervalTree
1306 if (node == (IntervalTree *) NULL)
1308 node->mean_stability=0.0;
1310 if (child != (IntervalTree *) NULL)
1315 register MagickRealType
1320 for ( ; child != (IntervalTree *) NULL; child=child->sibling)
1322 sum+=child->stability;
1325 node->mean_stability=sum/(MagickRealType) count;
1327 MeanStability(node->sibling);
1328 MeanStability(node->child);
1331 static void Stability(IntervalTree *node)
1333 if (node == (IntervalTree *) NULL)
1335 if (node->child == (IntervalTree *) NULL)
1336 node->stability=0.0;
1338 node->stability=node->tau-(node->child)->tau;
1339 Stability(node->sibling);
1340 Stability(node->child);
1343 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing,
1344 const unsigned long number_crossings)
1362 Allocate interval tree.
1364 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1366 if (list == (IntervalTree **) NULL)
1367 return((IntervalTree *) NULL);
1369 The root is the entire histogram.
1371 root=(IntervalTree *) AcquireMagickMemory(sizeof(*root));
1372 root->child=(IntervalTree *) NULL;
1373 root->sibling=(IntervalTree *) NULL;
1377 for (i=(-1); i < (long) number_crossings; i++)
1380 Initialize list with all nodes with no children.
1383 InitializeList(list,&number_nodes,root);
1387 for (j=0; j < number_nodes; j++)
1392 for (k=head->left+1; k < head->right; k++)
1394 if (zero_crossing[i+1].crossings[k] != 0)
1398 node->child=(IntervalTree *) AcquireMagickMemory(
1399 sizeof(*node->child));
1404 node->sibling=(IntervalTree *) AcquireMagickMemory(
1405 sizeof(*node->sibling));
1408 node->tau=zero_crossing[i+1].tau;
1409 node->child=(IntervalTree *) NULL;
1410 node->sibling=(IntervalTree *) NULL;
1416 if (left != head->left)
1418 node->sibling=(IntervalTree *) AcquireMagickMemory(
1419 sizeof(*node->sibling));
1421 node->tau=zero_crossing[i+1].tau;
1422 node->child=(IntervalTree *) NULL;
1423 node->sibling=(IntervalTree *) NULL;
1425 node->right=head->right;
1430 Determine the stability: difference between a nodes tau and its child.
1432 Stability(root->child);
1433 MeanStability(root->child);
1434 list=(IntervalTree **) RelinquishMagickMemory(list);
1439 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1443 + O p t i m a l T a u %
1447 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1449 % OptimalTau() finds the optimal tau for each band of the histogram.
1451 % The format of the OptimalTau method is:
1453 % MagickRealType OptimalTau(const long *histogram,const double max_tau,
1454 % const double min_tau,const double delta_tau,
1455 % const double smooth_threshold,short *extrema)
1457 % A description of each parameter follows.
1459 % o histogram: Specifies an array of integers representing the number
1460 % of pixels for each intensity of a particular color component.
1462 % o extrema: Specifies a pointer to an array of integers. They
1463 % represent the peaks and valleys of the histogram for each color
1468 static void ActiveNodes(IntervalTree **list,long *number_nodes,
1471 if (node == (IntervalTree *) NULL)
1473 if (node->stability >= node->mean_stability)
1475 list[(*number_nodes)++]=node;
1476 ActiveNodes(list,number_nodes,node->sibling);
1480 ActiveNodes(list,number_nodes,node->sibling);
1481 ActiveNodes(list,number_nodes,node->child);
1485 static void FreeNodes(IntervalTree *node)
1487 if (node == (IntervalTree *) NULL)
1489 FreeNodes(node->sibling);
1490 FreeNodes(node->child);
1491 node=(IntervalTree *) RelinquishMagickMemory(node);
1494 static MagickRealType OptimalTau(const long *histogram,const double max_tau,
1495 const double min_tau,const double delta_tau,const double smooth_threshold,
1531 Allocate interval tree.
1533 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1535 if (list == (IntervalTree **) NULL)
1538 Allocate zero crossing list.
1540 count=(unsigned long) ((max_tau-min_tau)/delta_tau)+2;
1541 zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count,
1542 sizeof(*zero_crossing));
1543 if (zero_crossing == (ZeroCrossing *) NULL)
1545 for (i=0; i < (long) count; i++)
1546 zero_crossing[i].tau=(-1.0);
1548 Initialize zero crossing list.
1550 derivative=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*derivative));
1551 second_derivative=(MagickRealType *) AcquireQuantumMemory(256,
1552 sizeof(*second_derivative));
1553 if ((derivative == (MagickRealType *) NULL) ||
1554 (second_derivative == (MagickRealType *) NULL))
1555 ThrowFatalException(ResourceLimitFatalError,
1556 "UnableToAllocateDerivatives");
1558 for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
1560 zero_crossing[i].tau=tau;
1561 ScaleSpace(histogram,tau,zero_crossing[i].histogram);
1562 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1563 DerivativeHistogram(derivative,second_derivative);
1564 ZeroCrossHistogram(second_derivative,smooth_threshold,
1565 zero_crossing[i].crossings);
1569 Add an entry for the original histogram.
1571 zero_crossing[i].tau=0.0;
1572 for (j=0; j <= 255; j++)
1573 zero_crossing[i].histogram[j]=(MagickRealType) histogram[j];
1574 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1575 DerivativeHistogram(derivative,second_derivative);
1576 ZeroCrossHistogram(second_derivative,smooth_threshold,
1577 zero_crossing[i].crossings);
1578 number_crossings=(unsigned long) i;
1579 derivative=(MagickRealType *) RelinquishMagickMemory(derivative);
1580 second_derivative=(MagickRealType *)
1581 RelinquishMagickMemory(second_derivative);
1583 Ensure the scale-space fingerprints form lines in scale-space, not loops.
1585 ConsolidateCrossings(zero_crossing,number_crossings);
1587 Force endpoints to be included in the interval.
1589 for (i=0; i <= (long) number_crossings; i++)
1591 for (j=0; j < 255; j++)
1592 if (zero_crossing[i].crossings[j] != 0)
1594 zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]);
1595 for (j=255; j > 0; j--)
1596 if (zero_crossing[i].crossings[j] != 0)
1598 zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1601 Initialize interval tree.
1603 root=InitializeIntervalTree(zero_crossing,number_crossings);
1604 if (root == (IntervalTree *) NULL)
1607 Find active nodes: stability is greater (or equal) to the mean stability of
1611 ActiveNodes(list,&number_nodes,root->child);
1615 for (i=0; i <= 255; i++)
1617 for (i=0; i < number_nodes; i++)
1620 Find this tau in zero crossings list.
1624 for (j=0; j <= (long) number_crossings; j++)
1625 if (zero_crossing[j].tau == node->tau)
1628 Find the value of the peak.
1630 peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1633 value=zero_crossing[k].histogram[index];
1634 for (x=node->left; x <= node->right; x++)
1636 if (peak != MagickFalse)
1638 if (zero_crossing[k].histogram[x] > value)
1640 value=zero_crossing[k].histogram[x];
1645 if (zero_crossing[k].histogram[x] < value)
1647 value=zero_crossing[k].histogram[x];
1651 for (x=node->left; x <= node->right; x++)
1655 if (peak != MagickFalse)
1656 extrema[x]=(short) index;
1658 extrema[x]=(short) (-index);
1662 Determine the average tau.
1665 for (i=0; i < number_nodes; i++)
1666 average_tau+=list[i]->tau;
1667 average_tau/=(MagickRealType) number_nodes;
1669 Relinquish resources.
1672 zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1673 list=(IntervalTree **) RelinquishMagickMemory(list);
1674 return(average_tau);
1678 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1682 + S c a l e S p a c e %
1686 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1688 % ScaleSpace() performs a scale-space filter on the 1D histogram.
1690 % The format of the ScaleSpace method is:
1692 % ScaleSpace(const long *histogram,const MagickRealType tau,
1693 % MagickRealType *scale_histogram)
1695 % A description of each parameter follows.
1697 % o histogram: Specifies an array of MagickRealTypes representing the number
1698 % of pixels for each intensity of a particular color component.
1702 static void ScaleSpace(const long *histogram,const MagickRealType tau,
1703 MagickRealType *scale_histogram)
1715 gamma=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*gamma));
1716 if (gamma == (MagickRealType *) NULL)
1717 ThrowFatalException(ResourceLimitFatalError,
1718 "UnableToAllocateGammaMap");
1719 alpha=1.0/(tau*sqrt(2.0*MagickPI));
1720 beta=(-1.0/(2.0*tau*tau));
1721 for (x=0; x <= 255; x++)
1723 for (x=0; x <= 255; x++)
1725 gamma[x]=exp((double) beta*x*x);
1726 if (gamma[x] < MagickEpsilon)
1729 for (x=0; x <= 255; x++)
1732 for (u=0; u <= 255; u++)
1733 sum+=(MagickRealType) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1734 scale_histogram[x]=alpha*sum;
1736 gamma=(MagickRealType *) RelinquishMagickMemory(gamma);
1740 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1744 % S e g m e n t I m a g e %
1748 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1750 % SegmentImage() segment an image by analyzing the histograms of the color
1751 % components and identifying units that are homogeneous with the fuzzy
1752 % C-means technique.
1754 % The format of the SegmentImage method is:
1756 % MagickBooleanType SegmentImage(Image *image,
1757 % const ColorspaceType colorspace,const MagickBooleanType verbose,
1758 % const double cluster_threshold,const double smooth_threshold)
1760 % A description of each parameter follows.
1762 % o image: the image.
1764 % o colorspace: Indicate the colorspace.
1766 % o verbose: Set to MagickTrue to print detailed information about the
1767 % identified classes.
1769 % o cluster_threshold: This represents the minimum number of pixels
1770 % contained in a hexahedra before it can be considered valid (expressed
1773 % o smooth_threshold: the smoothing threshold eliminates noise in the second
1774 % derivative of the histogram. As the value is increased, you can expect a
1775 % smoother second derivative.
1778 MagickExport MagickBooleanType SegmentImage(Image *image,
1779 const ColorspaceType colorspace,const MagickBooleanType verbose,
1780 const double cluster_threshold,const double smooth_threshold)
1783 *histogram[MaxDimension];
1792 *extrema[MaxDimension];
1795 Allocate histogram and extrema.
1797 assert(image != (Image *) NULL);
1798 assert(image->signature == MagickSignature);
1799 if (image->debug != MagickFalse)
1800 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1801 for (i=0; i < MaxDimension; i++)
1803 histogram[i]=(long *) AcquireQuantumMemory(256,sizeof(**histogram));
1804 extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1805 if ((histogram[i] == (long *) NULL) || (extrema[i] == (short *) NULL))
1807 for (i-- ; i >= 0; i--)
1809 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1810 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1812 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1816 if (colorspace != RGBColorspace)
1817 (void) TransformImageColorspace(image,colorspace);
1819 Initialize histogram.
1821 InitializeHistogram(image,histogram,&image->exception);
1822 (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,
1823 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Red]);
1824 (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,
1825 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Green]);
1826 (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,
1827 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Blue]);
1829 Classify using the fuzzy c-Means technique.
1831 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose);
1832 if (colorspace != RGBColorspace)
1833 (void) TransformImageColorspace(image,colorspace);
1835 Relinquish resources.
1837 for (i=0; i < MaxDimension; i++)
1839 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1840 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1846 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1850 + Z e r o C r o s s H i s t o g r a m %
1854 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1856 % ZeroCrossHistogram() find the zero crossings in a histogram and marks
1857 % directions as: 1 is negative to positive; 0 is zero crossing; and -1
1858 % is positive to negative.
1860 % The format of the ZeroCrossHistogram method is:
1862 % ZeroCrossHistogram(MagickRealType *second_derivative,
1863 % const MagickRealType smooth_threshold,short *crossings)
1865 % A description of each parameter follows.
1867 % o second_derivative: Specifies an array of MagickRealTypes representing the
1868 % second derivative of the histogram of a particular color component.
1870 % o crossings: This array of integers is initialized with
1871 % -1, 0, or 1 representing the slope of the first derivative of the
1872 % of a particular color component.
1875 static void ZeroCrossHistogram(MagickRealType *second_derivative,
1876 const MagickRealType smooth_threshold,short *crossings)
1885 Merge low numbers to zero to help prevent noise.
1887 for (i=0; i <= 255; i++)
1888 if ((second_derivative[i] < smooth_threshold) &&
1889 (second_derivative[i] >= -smooth_threshold))
1890 second_derivative[i]=0.0;
1892 Mark zero crossings.
1895 for (i=0; i <= 255; i++)
1898 if (second_derivative[i] < 0.0)
1905 if (second_derivative[i] > 0.0)