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"
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 *) 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(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",(double)
449 (void) fprintf(stdout,"\tWeighting Exponent = %g\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,"%g %g %g\n",(double) cluster->red.center,
484 (double) cluster->green.center,(double) cluster->blue.center);
486 (void) fprintf(stdout,"\n");
488 if (number_clusters > 256)
489 ThrowBinaryException(ImageError,"TooManyClusters",image->filename);
491 Speed up distance calculations.
493 squares=(MagickRealType *) AcquireQuantumMemory(513UL,sizeof(*squares));
494 if (squares == (MagickRealType *) NULL)
495 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
498 for (i=(-255); i <= 255; i++)
499 squares[i]=(MagickRealType) i*(MagickRealType) i;
501 Allocate image colormap.
503 if (AcquireImageColormap(image,number_clusters) == MagickFalse)
504 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
507 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
509 image->colormap[i].red=ScaleCharToQuantum((unsigned char)
510 (cluster->red.center+0.5));
511 image->colormap[i].green=ScaleCharToQuantum((unsigned char)
512 (cluster->green.center+0.5));
513 image->colormap[i].blue=ScaleCharToQuantum((unsigned char)
514 (cluster->blue.center+0.5));
518 Do course grain classes.
520 exception=(&image->exception);
521 image_view=AcquireCacheView(image);
522 #if defined(MAGICKCORE_OPENMP_SUPPORT)
523 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
525 for (y=0; y < (long) image->rows; y++)
530 register const PixelPacket
542 if (status == MagickFalse)
544 q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
545 if (q == (PixelPacket *) NULL)
550 indexes=GetCacheViewAuthenticIndexQueue(image_view);
551 for (x=0; x < (long) image->columns; x++)
553 indexes[x]=(IndexPacket) 0;
554 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
556 if (((long) ScaleQuantumToChar(q->red) >=
557 (cluster->red.left-SafeMargin)) &&
558 ((long) ScaleQuantumToChar(q->red) <=
559 (cluster->red.right+SafeMargin)) &&
560 ((long) ScaleQuantumToChar(q->green) >=
561 (cluster->green.left-SafeMargin)) &&
562 ((long) ScaleQuantumToChar(q->green) <=
563 (cluster->green.right+SafeMargin)) &&
564 ((long) ScaleQuantumToChar(q->blue) >=
565 (cluster->blue.left-SafeMargin)) &&
566 ((long) ScaleQuantumToChar(q->blue) <=
567 (cluster->blue.right+SafeMargin)))
572 indexes[x]=(IndexPacket) cluster->id;
576 if (cluster == (Cluster *) NULL)
590 Compute fuzzy membership.
593 for (j=0; j < (long) image->colors; j++)
597 distance_squared=squares[(long) ScaleQuantumToChar(q->red)-
598 (long) ScaleQuantumToChar(p->red)]+
599 squares[(long) ScaleQuantumToChar(q->green)-
600 (long) ScaleQuantumToChar(p->green)]+
601 squares[(long) ScaleQuantumToChar(q->blue)-
602 (long) ScaleQuantumToChar(p->blue)];
603 numerator=distance_squared;
604 for (k=0; k < (long) image->colors; k++)
607 distance_squared=squares[(long) ScaleQuantumToChar(q->red)-
608 (long) ScaleQuantumToChar(p->red)]+
609 squares[(long) ScaleQuantumToChar(q->green)-
610 (long) ScaleQuantumToChar(p->green)]+
611 squares[(long) ScaleQuantumToChar(q->blue)-
612 (long) ScaleQuantumToChar(p->blue)];
613 ratio=numerator/distance_squared;
614 sum+=SegmentPower(ratio);
616 if ((sum != 0.0) && ((1.0/sum) > local_minima))
621 local_minima=1.0/sum;
622 indexes[x]=(IndexPacket) j;
628 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
630 if (image->progress_monitor != (MagickProgressMonitor) NULL)
635 #if defined(MAGICKCORE_OPENMP_SUPPORT)
636 #pragma omp critical (MagickCore_Classify)
638 proceed=SetImageProgress(image,SegmentImageTag,progress++,
640 if (proceed == MagickFalse)
644 image_view=DestroyCacheView(image_view);
645 status&=SyncImage(image);
647 Relinquish resources.
649 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
651 next_cluster=cluster->next;
652 cluster=(Cluster *) RelinquishMagickMemory(cluster);
655 free_squares=squares;
656 free_squares=(MagickRealType *) RelinquishMagickMemory(free_squares);
661 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
665 + C o n s o l i d a t e C r o s s i n g s %
669 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
671 % ConsolidateCrossings() guarantees that an even number of zero crossings
672 % always lie between two crossings.
674 % The format of the ConsolidateCrossings method is:
676 % ConsolidateCrossings(ZeroCrossing *zero_crossing,
677 % const unsigned long number_crossings)
679 % A description of each parameter follows.
681 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
683 % o number_crossings: This unsigned long specifies the number of elements
684 % in the zero_crossing array.
688 static inline long MagickAbsoluteValue(const long x)
695 static inline long MagickMax(const long x,const long y)
702 static inline long MagickMin(const long x,const long y)
709 static void ConsolidateCrossings(ZeroCrossing *zero_crossing,
710 const unsigned long number_crossings)
726 Consolidate zero crossings.
728 for (i=(long) number_crossings-1; i >= 0; i--)
729 for (j=0; j <= 255; j++)
731 if (zero_crossing[i].crossings[j] == 0)
734 Find the entry that is closest to j and still preserves the
735 property that there are an even number of crossings between
738 for (k=j-1; k > 0; k--)
739 if (zero_crossing[i+1].crossings[k] != 0)
743 for (k=j+1; k < 255; k++)
744 if (zero_crossing[i+1].crossings[k] != 0)
746 right=MagickMin(k,255);
748 K is the zero crossing just left of j.
750 for (k=j-1; k > 0; k--)
751 if (zero_crossing[i].crossings[k] != 0)
756 Check center for an even number of crossings between k and j.
759 if (zero_crossing[i+1].crossings[j] != 0)
762 for (l=k+1; l < center; l++)
763 if (zero_crossing[i+1].crossings[l] != 0)
765 if (((count % 2) == 0) && (center != k))
769 Check left for an even number of crossings between k and j.
774 for (l=k+1; l < left; l++)
775 if (zero_crossing[i+1].crossings[l] != 0)
777 if (((count % 2) == 0) && (left != k))
781 Check right for an even number of crossings between k and j.
786 for (l=k+1; l < right; l++)
787 if (zero_crossing[i+1].crossings[l] != 0)
789 if (((count % 2) == 0) && (right != k))
792 l=zero_crossing[i].crossings[j];
793 zero_crossing[i].crossings[j]=0;
795 zero_crossing[i].crossings[correct]=(short) l;
800 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
804 + D e f i n e R e g i o n %
808 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
810 % DefineRegion() defines the left and right boundaries of a peak region.
812 % The format of the DefineRegion method is:
814 % long DefineRegion(const short *extrema,ExtentPacket *extents)
816 % A description of each parameter follows.
818 % o extrema: Specifies a pointer to an array of integers. They
819 % represent the peaks and valleys of the histogram for each color
822 % o extents: This pointer to an ExtentPacket represent the extends
823 % of a particular peak or valley of a color component.
826 static long DefineRegion(const short *extrema,ExtentPacket *extents)
829 Initialize to default values.
835 Find the left side (maxima).
837 for ( ; extents->index <= 255; extents->index++)
838 if (extrema[extents->index] > 0)
840 if (extents->index > 255)
841 return(MagickFalse); /* no left side - no region exists */
842 extents->left=extents->index;
844 Find the right side (minima).
846 for ( ; extents->index <= 255; extents->index++)
847 if (extrema[extents->index] < 0)
849 extents->right=extents->index-1;
854 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
858 + D e r i v a t i v e H i s t o g r a m %
862 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
864 % DerivativeHistogram() determines the derivative of the histogram using
865 % central differencing.
867 % The format of the DerivativeHistogram method is:
869 % DerivativeHistogram(const MagickRealType *histogram,
870 % MagickRealType *derivative)
872 % A description of each parameter follows.
874 % o histogram: Specifies an array of MagickRealTypes representing the number
875 % of pixels for each intensity of a particular color component.
877 % o derivative: This array of MagickRealTypes is initialized by
878 % DerivativeHistogram to the derivative of the histogram using central
882 static void DerivativeHistogram(const MagickRealType *histogram,
883 MagickRealType *derivative)
890 Compute endpoints using second order polynomial interpolation.
893 derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);
894 derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);
896 Compute derivative using central differencing.
898 for (i=1; i < n; i++)
899 derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;
904 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
908 + G e t I m a g e D y n a m i c T h r e s h o l d %
912 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
914 % GetImageDynamicThreshold() returns the dynamic threshold for an image.
916 % The format of the GetImageDynamicThreshold method is:
918 % MagickBooleanType GetImageDynamicThreshold(const Image *image,
919 % const double cluster_threshold,const double smooth_threshold,
920 % MagickPixelPacket *pixel,ExceptionInfo *exception)
922 % A description of each parameter follows.
924 % o image: the image.
926 % o cluster_threshold: This MagickRealType represents the minimum number of
927 % pixels contained in a hexahedra before it can be considered valid
928 % (expressed as a percentage).
930 % o smooth_threshold: the smoothing threshold eliminates noise in the second
931 % derivative of the histogram. As the value is increased, you can expect a
932 % smoother second derivative.
934 % o pixel: return the dynamic threshold here.
936 % o exception: return any errors or warnings in this structure.
939 MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,
940 const double cluster_threshold,const double smooth_threshold,
941 MagickPixelPacket *pixel,ExceptionInfo *exception)
958 *histogram[MaxDimension],
967 register const PixelPacket
975 *extrema[MaxDimension];
978 Allocate histogram and extrema.
980 assert(image != (Image *) NULL);
981 assert(image->signature == MagickSignature);
982 if (image->debug != MagickFalse)
983 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
984 GetMagickPixelPacket(image,pixel);
985 for (i=0; i < MaxDimension; i++)
987 histogram[i]=(long *) AcquireQuantumMemory(256UL,sizeof(**histogram));
988 extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
989 if ((histogram[i] == (long *) NULL) || (extrema[i] == (short *) NULL))
991 for (i-- ; i >= 0; i--)
993 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
994 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
996 (void) ThrowMagickException(exception,GetMagickModule(),
997 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1002 Initialize histogram.
1004 InitializeHistogram(image,histogram,exception);
1005 (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
1006 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
1007 (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
1008 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
1009 (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
1010 (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
1014 cluster=(Cluster *) NULL;
1015 head=(Cluster *) NULL;
1016 (void) ResetMagickMemory(&red,0,sizeof(red));
1017 (void) ResetMagickMemory(&green,0,sizeof(green));
1018 (void) ResetMagickMemory(&blue,0,sizeof(blue));
1019 while (DefineRegion(extrema[Red],&red) != 0)
1022 while (DefineRegion(extrema[Green],&green) != 0)
1025 while (DefineRegion(extrema[Blue],&blue) != 0)
1028 Allocate a new class.
1030 if (head != (Cluster *) NULL)
1032 cluster->next=(Cluster *) AcquireMagickMemory(
1033 sizeof(*cluster->next));
1034 cluster=cluster->next;
1038 cluster=(Cluster *) AcquireAlignedMemory(1,sizeof(*cluster));
1041 if (cluster == (Cluster *) NULL)
1043 (void) ThrowMagickException(exception,GetMagickModule(),
1044 ResourceLimitError,"MemoryAllocationFailed","`%s'",
1046 return(MagickFalse);
1049 Initialize a new class.
1053 cluster->green=green;
1055 cluster->next=(Cluster *) NULL;
1059 if (head == (Cluster *) NULL)
1062 No classes were identified-- create one.
1064 cluster=(Cluster *) AcquireAlignedMemory(1,sizeof(*cluster));
1065 if (cluster == (Cluster *) NULL)
1067 (void) ThrowMagickException(exception,GetMagickModule(),
1068 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1069 return(MagickFalse);
1072 Initialize a new class.
1076 cluster->green=green;
1078 cluster->next=(Cluster *) NULL;
1082 Count the pixels for each cluster.
1085 for (y=0; y < (long) image->rows; y++)
1087 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1088 if (p == (const PixelPacket *) NULL)
1090 for (x=0; x < (long) image->columns; x++)
1092 for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
1093 if (((long) ScaleQuantumToChar(p->red) >=
1094 (cluster->red.left-SafeMargin)) &&
1095 ((long) ScaleQuantumToChar(p->red) <=
1096 (cluster->red.right+SafeMargin)) &&
1097 ((long) ScaleQuantumToChar(p->green) >=
1098 (cluster->green.left-SafeMargin)) &&
1099 ((long) ScaleQuantumToChar(p->green) <=
1100 (cluster->green.right+SafeMargin)) &&
1101 ((long) ScaleQuantumToChar(p->blue) >=
1102 (cluster->blue.left-SafeMargin)) &&
1103 ((long) ScaleQuantumToChar(p->blue) <=
1104 (cluster->blue.right+SafeMargin)))
1110 cluster->red.center+=(MagickRealType)
1111 ScaleQuantumToChar(p->red);
1112 cluster->green.center+=(MagickRealType)
1113 ScaleQuantumToChar(p->green);
1114 cluster->blue.center+=(MagickRealType)
1115 ScaleQuantumToChar(p->blue);
1121 proceed=SetImageProgress(image,SegmentImageTag,y,2*image->rows);
1122 if (proceed == MagickFalse)
1126 Remove clusters that do not meet minimum cluster threshold.
1131 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1133 next_cluster=cluster->next;
1134 if ((cluster->count > 0) &&
1135 (cluster->count >= (count*cluster_threshold/100.0)))
1141 cluster->red.center/=cluster->count;
1142 cluster->green.center/=cluster->count;
1143 cluster->blue.center/=cluster->count;
1145 last_cluster=cluster;
1151 if (cluster == head)
1154 last_cluster->next=next_cluster;
1155 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1162 for (cluster=object; cluster->next != (Cluster *) NULL; )
1164 if (cluster->count < object->count)
1166 cluster=cluster->next;
1168 background=head->next;
1169 for (cluster=background; cluster->next != (Cluster *) NULL; )
1171 if (cluster->count > background->count)
1173 cluster=cluster->next;
1176 threshold=(background->red.center+object->red.center)/2.0;
1177 pixel->red=(MagickRealType) ScaleCharToQuantum((unsigned char)
1179 threshold=(background->green.center+object->green.center)/2.0;
1180 pixel->green=(MagickRealType) ScaleCharToQuantum((unsigned char)
1182 threshold=(background->blue.center+object->blue.center)/2.0;
1183 pixel->blue=(MagickRealType) ScaleCharToQuantum((unsigned char)
1186 Relinquish resources.
1188 for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
1190 next_cluster=cluster->next;
1191 cluster=(Cluster *) RelinquishMagickMemory(cluster);
1193 for (i=0; i < MaxDimension; i++)
1195 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1196 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1202 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1206 + I n i t i a l i z e H i s t o g r a m %
1210 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1212 % InitializeHistogram() computes the histogram for an image.
1214 % The format of the InitializeHistogram method is:
1216 % InitializeHistogram(const Image *image,long **histogram)
1218 % A description of each parameter follows.
1220 % o image: Specifies a pointer to an Image structure; returned from
1223 % o histogram: Specifies an array of integers representing the number
1224 % of pixels for each intensity of a particular color component.
1227 static void InitializeHistogram(const Image *image,long **histogram,
1228 ExceptionInfo *exception)
1233 register const PixelPacket
1241 Initialize histogram.
1243 for (i=0; i <= 255; i++)
1245 histogram[Red][i]=0;
1246 histogram[Green][i]=0;
1247 histogram[Blue][i]=0;
1249 for (y=0; y < (long) image->rows; y++)
1251 p=GetVirtualPixels(image,0,y,image->columns,1,exception);
1252 if (p == (const PixelPacket *) NULL)
1254 for (x=0; x < (long) image->columns; x++)
1256 histogram[Red][(long) ScaleQuantumToChar(p->red)]++;
1257 histogram[Green][(long) ScaleQuantumToChar(p->green)]++;
1258 histogram[Blue][(long) ScaleQuantumToChar(p->blue)]++;
1265 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1269 + I n i t i a l i z e I n t e r v a l T r e e %
1273 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1275 % InitializeIntervalTree() initializes an interval tree from the lists of
1278 % The format of the InitializeIntervalTree method is:
1280 % InitializeIntervalTree(IntervalTree **list,long *number_nodes,
1281 % IntervalTree *node)
1283 % A description of each parameter follows.
1285 % o zero_crossing: Specifies an array of structures of type ZeroCrossing.
1287 % o number_crossings: This unsigned long specifies the number of elements
1288 % in the zero_crossing array.
1292 static void InitializeList(IntervalTree **list,long *number_nodes,
1295 if (node == (IntervalTree *) NULL)
1297 if (node->child == (IntervalTree *) NULL)
1298 list[(*number_nodes)++]=node;
1299 InitializeList(list,number_nodes,node->sibling);
1300 InitializeList(list,number_nodes,node->child);
1303 static void MeanStability(IntervalTree *node)
1305 register IntervalTree
1308 if (node == (IntervalTree *) NULL)
1310 node->mean_stability=0.0;
1312 if (child != (IntervalTree *) NULL)
1317 register MagickRealType
1322 for ( ; child != (IntervalTree *) NULL; child=child->sibling)
1324 sum+=child->stability;
1327 node->mean_stability=sum/(MagickRealType) count;
1329 MeanStability(node->sibling);
1330 MeanStability(node->child);
1333 static void Stability(IntervalTree *node)
1335 if (node == (IntervalTree *) NULL)
1337 if (node->child == (IntervalTree *) NULL)
1338 node->stability=0.0;
1340 node->stability=node->tau-(node->child)->tau;
1341 Stability(node->sibling);
1342 Stability(node->child);
1345 static IntervalTree *InitializeIntervalTree(const ZeroCrossing *zero_crossing,
1346 const unsigned long number_crossings)
1364 Allocate interval tree.
1366 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1368 if (list == (IntervalTree **) NULL)
1369 return((IntervalTree *) NULL);
1371 The root is the entire histogram.
1373 root=(IntervalTree *) AcquireAlignedMemory(1,sizeof(*root));
1374 root->child=(IntervalTree *) NULL;
1375 root->sibling=(IntervalTree *) NULL;
1379 for (i=(-1); i < (long) number_crossings; i++)
1382 Initialize list with all nodes with no children.
1385 InitializeList(list,&number_nodes,root);
1389 for (j=0; j < number_nodes; j++)
1394 for (k=head->left+1; k < head->right; k++)
1396 if (zero_crossing[i+1].crossings[k] != 0)
1400 node->child=(IntervalTree *) AcquireMagickMemory(
1401 sizeof(*node->child));
1406 node->sibling=(IntervalTree *) AcquireMagickMemory(
1407 sizeof(*node->sibling));
1410 node->tau=zero_crossing[i+1].tau;
1411 node->child=(IntervalTree *) NULL;
1412 node->sibling=(IntervalTree *) NULL;
1418 if (left != head->left)
1420 node->sibling=(IntervalTree *) AcquireMagickMemory(
1421 sizeof(*node->sibling));
1423 node->tau=zero_crossing[i+1].tau;
1424 node->child=(IntervalTree *) NULL;
1425 node->sibling=(IntervalTree *) NULL;
1427 node->right=head->right;
1432 Determine the stability: difference between a nodes tau and its child.
1434 Stability(root->child);
1435 MeanStability(root->child);
1436 list=(IntervalTree **) RelinquishMagickMemory(list);
1441 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1445 + O p t i m a l T a u %
1449 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1451 % OptimalTau() finds the optimal tau for each band of the histogram.
1453 % The format of the OptimalTau method is:
1455 % MagickRealType OptimalTau(const long *histogram,const double max_tau,
1456 % const double min_tau,const double delta_tau,
1457 % const double smooth_threshold,short *extrema)
1459 % A description of each parameter follows.
1461 % o histogram: Specifies an array of integers representing the number
1462 % of pixels for each intensity of a particular color component.
1464 % o extrema: Specifies a pointer to an array of integers. They
1465 % represent the peaks and valleys of the histogram for each color
1470 static void ActiveNodes(IntervalTree **list,long *number_nodes,
1473 if (node == (IntervalTree *) NULL)
1475 if (node->stability >= node->mean_stability)
1477 list[(*number_nodes)++]=node;
1478 ActiveNodes(list,number_nodes,node->sibling);
1482 ActiveNodes(list,number_nodes,node->sibling);
1483 ActiveNodes(list,number_nodes,node->child);
1487 static void FreeNodes(IntervalTree *node)
1489 if (node == (IntervalTree *) NULL)
1491 FreeNodes(node->sibling);
1492 FreeNodes(node->child);
1493 node=(IntervalTree *) RelinquishMagickMemory(node);
1496 static MagickRealType OptimalTau(const long *histogram,const double max_tau,
1497 const double min_tau,const double delta_tau,const double smooth_threshold,
1533 Allocate interval tree.
1535 list=(IntervalTree **) AcquireQuantumMemory((size_t) TreeLength,
1537 if (list == (IntervalTree **) NULL)
1540 Allocate zero crossing list.
1542 count=(unsigned long) ((max_tau-min_tau)/delta_tau)+2;
1543 zero_crossing=(ZeroCrossing *) AcquireQuantumMemory((size_t) count,
1544 sizeof(*zero_crossing));
1545 if (zero_crossing == (ZeroCrossing *) NULL)
1547 for (i=0; i < (long) count; i++)
1548 zero_crossing[i].tau=(-1.0);
1550 Initialize zero crossing list.
1552 derivative=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*derivative));
1553 second_derivative=(MagickRealType *) AcquireQuantumMemory(256,
1554 sizeof(*second_derivative));
1555 if ((derivative == (MagickRealType *) NULL) ||
1556 (second_derivative == (MagickRealType *) NULL))
1557 ThrowFatalException(ResourceLimitFatalError,
1558 "UnableToAllocateDerivatives");
1560 for (tau=max_tau; tau >= min_tau; tau-=delta_tau)
1562 zero_crossing[i].tau=tau;
1563 ScaleSpace(histogram,tau,zero_crossing[i].histogram);
1564 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1565 DerivativeHistogram(derivative,second_derivative);
1566 ZeroCrossHistogram(second_derivative,smooth_threshold,
1567 zero_crossing[i].crossings);
1571 Add an entry for the original histogram.
1573 zero_crossing[i].tau=0.0;
1574 for (j=0; j <= 255; j++)
1575 zero_crossing[i].histogram[j]=(MagickRealType) histogram[j];
1576 DerivativeHistogram(zero_crossing[i].histogram,derivative);
1577 DerivativeHistogram(derivative,second_derivative);
1578 ZeroCrossHistogram(second_derivative,smooth_threshold,
1579 zero_crossing[i].crossings);
1580 number_crossings=(unsigned long) i;
1581 derivative=(MagickRealType *) RelinquishMagickMemory(derivative);
1582 second_derivative=(MagickRealType *)
1583 RelinquishMagickMemory(second_derivative);
1585 Ensure the scale-space fingerprints form lines in scale-space, not loops.
1587 ConsolidateCrossings(zero_crossing,number_crossings);
1589 Force endpoints to be included in the interval.
1591 for (i=0; i <= (long) number_crossings; i++)
1593 for (j=0; j < 255; j++)
1594 if (zero_crossing[i].crossings[j] != 0)
1596 zero_crossing[i].crossings[0]=(-zero_crossing[i].crossings[j]);
1597 for (j=255; j > 0; j--)
1598 if (zero_crossing[i].crossings[j] != 0)
1600 zero_crossing[i].crossings[255]=(-zero_crossing[i].crossings[j]);
1603 Initialize interval tree.
1605 root=InitializeIntervalTree(zero_crossing,number_crossings);
1606 if (root == (IntervalTree *) NULL)
1609 Find active nodes: stability is greater (or equal) to the mean stability of
1613 ActiveNodes(list,&number_nodes,root->child);
1617 for (i=0; i <= 255; i++)
1619 for (i=0; i < number_nodes; i++)
1622 Find this tau in zero crossings list.
1626 for (j=0; j <= (long) number_crossings; j++)
1627 if (zero_crossing[j].tau == node->tau)
1630 Find the value of the peak.
1632 peak=zero_crossing[k].crossings[node->right] == -1 ? MagickTrue :
1635 value=zero_crossing[k].histogram[index];
1636 for (x=node->left; x <= node->right; x++)
1638 if (peak != MagickFalse)
1640 if (zero_crossing[k].histogram[x] > value)
1642 value=zero_crossing[k].histogram[x];
1647 if (zero_crossing[k].histogram[x] < value)
1649 value=zero_crossing[k].histogram[x];
1653 for (x=node->left; x <= node->right; x++)
1657 if (peak != MagickFalse)
1658 extrema[x]=(short) index;
1660 extrema[x]=(short) (-index);
1664 Determine the average tau.
1667 for (i=0; i < number_nodes; i++)
1668 average_tau+=list[i]->tau;
1669 average_tau/=(MagickRealType) number_nodes;
1671 Relinquish resources.
1674 zero_crossing=(ZeroCrossing *) RelinquishMagickMemory(zero_crossing);
1675 list=(IntervalTree **) RelinquishMagickMemory(list);
1676 return(average_tau);
1680 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1684 + S c a l e S p a c e %
1688 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1690 % ScaleSpace() performs a scale-space filter on the 1D histogram.
1692 % The format of the ScaleSpace method is:
1694 % ScaleSpace(const long *histogram,const MagickRealType tau,
1695 % MagickRealType *scale_histogram)
1697 % A description of each parameter follows.
1699 % o histogram: Specifies an array of MagickRealTypes representing the number
1700 % of pixels for each intensity of a particular color component.
1704 static void ScaleSpace(const long *histogram,const MagickRealType tau,
1705 MagickRealType *scale_histogram)
1717 gamma=(MagickRealType *) AcquireQuantumMemory(256,sizeof(*gamma));
1718 if (gamma == (MagickRealType *) NULL)
1719 ThrowFatalException(ResourceLimitFatalError,
1720 "UnableToAllocateGammaMap");
1721 alpha=1.0/(tau*sqrt(2.0*MagickPI));
1722 beta=(-1.0/(2.0*tau*tau));
1723 for (x=0; x <= 255; x++)
1725 for (x=0; x <= 255; x++)
1727 gamma[x]=exp((double) beta*x*x);
1728 if (gamma[x] < MagickEpsilon)
1731 for (x=0; x <= 255; x++)
1734 for (u=0; u <= 255; u++)
1735 sum+=(MagickRealType) histogram[u]*gamma[MagickAbsoluteValue(x-u)];
1736 scale_histogram[x]=alpha*sum;
1738 gamma=(MagickRealType *) RelinquishMagickMemory(gamma);
1742 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1746 % S e g m e n t I m a g e %
1750 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1752 % SegmentImage() segment an image by analyzing the histograms of the color
1753 % components and identifying units that are homogeneous with the fuzzy
1754 % C-means technique.
1756 % The format of the SegmentImage method is:
1758 % MagickBooleanType SegmentImage(Image *image,
1759 % const ColorspaceType colorspace,const MagickBooleanType verbose,
1760 % const double cluster_threshold,const double smooth_threshold)
1762 % A description of each parameter follows.
1764 % o image: the image.
1766 % o colorspace: Indicate the colorspace.
1768 % o verbose: Set to MagickTrue to print detailed information about the
1769 % identified classes.
1771 % o cluster_threshold: This represents the minimum number of pixels
1772 % contained in a hexahedra before it can be considered valid (expressed
1775 % o smooth_threshold: the smoothing threshold eliminates noise in the second
1776 % derivative of the histogram. As the value is increased, you can expect a
1777 % smoother second derivative.
1780 MagickExport MagickBooleanType SegmentImage(Image *image,
1781 const ColorspaceType colorspace,const MagickBooleanType verbose,
1782 const double cluster_threshold,const double smooth_threshold)
1785 *histogram[MaxDimension];
1794 *extrema[MaxDimension];
1797 Allocate histogram and extrema.
1799 assert(image != (Image *) NULL);
1800 assert(image->signature == MagickSignature);
1801 if (image->debug != MagickFalse)
1802 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1803 for (i=0; i < MaxDimension; i++)
1805 histogram[i]=(long *) AcquireQuantumMemory(256,sizeof(**histogram));
1806 extrema[i]=(short *) AcquireQuantumMemory(256,sizeof(**extrema));
1807 if ((histogram[i] == (long *) NULL) || (extrema[i] == (short *) NULL))
1809 for (i-- ; i >= 0; i--)
1811 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1812 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1814 ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1818 if (colorspace != RGBColorspace)
1819 (void) TransformImageColorspace(image,colorspace);
1821 Initialize histogram.
1823 InitializeHistogram(image,histogram,&image->exception);
1824 (void) OptimalTau(histogram[Red],Tau,0.2,DeltaTau,
1825 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Red]);
1826 (void) OptimalTau(histogram[Green],Tau,0.2,DeltaTau,
1827 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Green]);
1828 (void) OptimalTau(histogram[Blue],Tau,0.2,DeltaTau,
1829 smooth_threshold == 0.0 ? 1.0 : smooth_threshold,extrema[Blue]);
1831 Classify using the fuzzy c-Means technique.
1833 status=Classify(image,extrema,cluster_threshold,WeightingExponent,verbose);
1834 if (colorspace != RGBColorspace)
1835 (void) TransformImageColorspace(image,colorspace);
1837 Relinquish resources.
1839 for (i=0; i < MaxDimension; i++)
1841 extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
1842 histogram[i]=(long *) RelinquishMagickMemory(histogram[i]);
1848 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1852 + Z e r o C r o s s H i s t o g r a m %
1856 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1858 % ZeroCrossHistogram() find the zero crossings in a histogram and marks
1859 % directions as: 1 is negative to positive; 0 is zero crossing; and -1
1860 % is positive to negative.
1862 % The format of the ZeroCrossHistogram method is:
1864 % ZeroCrossHistogram(MagickRealType *second_derivative,
1865 % const MagickRealType smooth_threshold,short *crossings)
1867 % A description of each parameter follows.
1869 % o second_derivative: Specifies an array of MagickRealTypes representing the
1870 % second derivative of the histogram of a particular color component.
1872 % o crossings: This array of integers is initialized with
1873 % -1, 0, or 1 representing the slope of the first derivative of the
1874 % of a particular color component.
1877 static void ZeroCrossHistogram(MagickRealType *second_derivative,
1878 const MagickRealType smooth_threshold,short *crossings)
1887 Merge low numbers to zero to help prevent noise.
1889 for (i=0; i <= 255; i++)
1890 if ((second_derivative[i] < smooth_threshold) &&
1891 (second_derivative[i] >= -smooth_threshold))
1892 second_derivative[i]=0.0;
1894 Mark zero crossings.
1897 for (i=0; i <= 255; i++)
1900 if (second_derivative[i] < 0.0)
1907 if (second_derivative[i] > 0.0)