/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

%               SSSSS  EEEEE   GGGG  M   M  EEEEE  N   N  TTTTT               %

%               SS     E      G      MM MM  E      NN  N    T                 %

%                SSS   EEE    G GGG  M M M  EEE    N N N    T                 %

%                  SS  E      G   G  M   M  E      N  NN    T                 %

%               SSSSS  EEEEE   GGGG  M   M  EEEEE  N   N    T                 %

%                                                                             %

%                                                                             %

%    MagickCore Methods to Segment an Image with Thresholding Fuzzy c-Means   %

%                                                                             %

%                              Software Design                                %

%                                   Cristy                                    %

%                                April 1993                                   %

%                                                                             %

%                                                                             %

%  Copyright @ 1999 ImageMagick Studio LLC, a non-profit organization         %

%  dedicated to making software imaging solutions freely available.           %

%                                                                             %

%  You may not use this file except in compliance with the License.  You may  %

%  obtain a copy of the License at                                            %

%                                                                             %

%    https://imagemagick.org/license/                                         %

%                                                                             %

%  Unless required by applicable law or agreed to in writing, software        %

%  distributed under the License is distributed on an "AS IS" BASIS,          %

%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %

%  See the License for the specific language governing permissions and        %

%  limitations under the License.                                             %

%                                                                             %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%  Segment segments an image by analyzing the histograms of the color

%  components and identifying units that are homogeneous with the fuzzy

%  c-means technique.  The scale-space filter analyzes the histograms of

%  the three color components of the image and identifies a set of

%  classes.  The extents of each class is used to coarsely segment the

%  image with thresholding.  The color associated with each class is

%  determined by the mean color of all pixels within the extents of a

%  particular class.  Finally, any unclassified pixels are assigned to

%  the closest class with the fuzzy c-means technique.

%

%  The fuzzy c-Means algorithm can be summarized as follows:

%

%    o Build a histogram, one for each color component of the image.

%

%    o For each histogram, successively apply the scale-space filter and

%      build an interval tree of zero crossings in the second derivative

%      at each scale.  Analyze this scale-space ''fingerprint'' to

%      determine which peaks and valleys in the histogram are most

%      predominant.

%

%    o The fingerprint defines intervals on the axis of the histogram.

%      Each interval contains either a minima or a maxima in the original

%      signal.  If each color component lies within the maxima interval,

%      that pixel is considered ''classified'' and is assigned an unique

%      class number.

%

%    o Any pixel that fails to be classified in the above thresholding

%      pass is classified using the fuzzy c-Means technique.  It is

%      assigned to one of the classes discovered in the histogram analysis

%      phase.

%

%  The fuzzy c-Means technique attempts to cluster a pixel by finding

%  the local minima of the generalized within group sum of squared error

%  objective function.  A pixel is assigned to the closest class of

%  which the fuzzy membership has a maximum value.

%

%  Segment is strongly based on software written by Andy Gallo,

%  University of Delaware.

%

%  The following reference was used in creating this program:

%

%    Young Won Lim, Sang Uk Lee, "On The Color Image Segmentation

%    Algorithm Based on the Thresholding and the Fuzzy c-Means

%    Techniques", Pattern Recognition, Volume 23, Number 9, pages

%    935-952, 1990.

%

%

*/

#include "MagickCore/studio.h"

#include "MagickCore/cache.h"

#include "MagickCore/color.h"

#include "MagickCore/colormap.h"

#include "MagickCore/colorspace.h"

#include "MagickCore/colorspace-private.h"

#include "MagickCore/exception.h"

#include "MagickCore/exception-private.h"

#include "MagickCore/image.h"

#include "MagickCore/image-private.h"

#include "MagickCore/memory_.h"

#include "MagickCore/memory-private.h"

#include "MagickCore/monitor.h"

#include "MagickCore/monitor-private.h"

#include "MagickCore/pixel-accessor.h"

#include "MagickCore/quantize.h"

#include "MagickCore/quantum.h"

#include "MagickCore/quantum-private.h"

#include "MagickCore/resource_.h"

#include "MagickCore/segment.h"

#include "MagickCore/string_.h"

#include "MagickCore/thread-private.h"

/*

  Define declarations.

*/

#define MaxDimension  3

#define DeltaTau  0.5f

#if defined(FastClassify)

#define WeightingExponent  2.0

#define SegmentPower(ratio) (ratio)

#else

#define WeightingExponent  2.5

#define SegmentPower(ratio) pow(ratio,(double) (1.0/(weighting_exponent-1.0)));

#endif

#define Tau  5.2f

/*

  Typedef declarations.

*/

typedef struct _ExtentPacket

{

  double

    center;

  ssize_t

    index,

    left,

    right;

} ExtentPacket;

typedef struct _Cluster

{

  struct _Cluster

    *next;

  ExtentPacket

    red,

    green,

    blue;

  ssize_t

    count,

    id;

} Cluster;

typedef struct _IntervalTree

{

  double

    tau;

  ssize_t

    left,

    right;

  double

    mean_stability,

    stability;

  struct _IntervalTree

    *sibling,

    *child;

} IntervalTree;

typedef struct _ZeroCrossing

{

  double

    tau,

    histogram[256];

  short

    crossings[256];

} ZeroCrossing;

/*

  Constant declarations.

*/

static const int

  Blue = 2,

  Green = 1,

  Red = 0,

  SafeMargin = 3,

  TreeLength = 600;

/*

  Method prototypes.

*/

static double

  OptimalTau(const ssize_t *,const double,const double,const double,

    const double,short *);

static ssize_t

  DefineRegion(const short *,ExtentPacket *);

static void

  FreeNodes(IntervalTree *),

  InitializeHistogram(const Image *,ssize_t **,ExceptionInfo *),

  ScaleSpace(const ssize_t *,const double,double *),

  ZeroCrossHistogram(double *,const double,short *);

/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

+   C l a s s i f y                                                           %

%                                                                             %

%                                                                             %

%                                                                             %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%  Classify() defines one or more classes.  Each pixel is thresholded to

%  determine which class it belongs to.  If the class is not identified it is

%  assigned to the closest class based on the fuzzy c-Means technique.

%

%  The format of the Classify method is:

%

%      MagickBooleanType Classify(Image *image,short **extrema,

%        const double cluster_threshold,const double weighting_exponent,

%        const MagickBooleanType verbose,ExceptionInfo *exception)

%

%  A description of each parameter follows.

%

%    o image: the image.

%

%    o extrema:  Specifies a pointer to an array of integers.  They

%      represent the peaks and valleys of the histogram for each color

%      component.

%

%    o cluster_threshold:  This double represents the minimum number of

%      pixels contained in a hexahedra before it can be considered valid

%      (expressed as a percentage).

%

%    o weighting_exponent: Specifies the membership weighting exponent.

%

%    o verbose:  A value greater than zero prints detailed information about

%      the identified classes.

%

%    o exception: return any errors or warnings in this structure.

%

*/

static MagickBooleanType Classify(Image *image,short **extrema,

  const double cluster_threshold,const double weighting_exponent,

  const MagickBooleanType verbose,ExceptionInfo *exception)

{

#define SegmentImageTag  "Segment/Image"

#define ThrowClassifyException(severity,tag,label) \

{\

  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster) \

  { \

    next_cluster=cluster->next; \

    cluster=(Cluster *) RelinquishMagickMemory(cluster); \

  } \

  if (squares != (double *) NULL) \

    { \

      squares-=255; \

      free_squares=squares; \

      free_squares=(double *) RelinquishMagickMemory(free_squares); \

    } \

  ThrowBinaryException(severity,tag,label); \

}

  CacheView

    *image_view;

  Cluster

    *cluster,

    *head,

    *last_cluster,

    *next_cluster;

  double

    *free_squares;

  ExtentPacket

    blue,

    green,

    red;

  MagickOffsetType

    progress;

  MagickStatusType

    status;

  ssize_t

    i;

  double

    *squares;

  size_t

    number_clusters;

  ssize_t

    count,

    y;

  /*

    Form clusters.

  */

  cluster=(Cluster *) NULL;

  head=(Cluster *) NULL;

  squares=(double *) NULL;

  (void) memset(&red,0,sizeof(red));

  (void) memset(&green,0,sizeof(green));

  (void) memset(&blue,0,sizeof(blue));

  while (DefineRegion(extrema[Red],&red) != 0)

  {

    green.index=0;

    while (DefineRegion(extrema[Green],&green) != 0)

    {

      blue.index=0;

      while (DefineRegion(extrema[Blue],&blue) != 0)

      {

        /*

          Allocate a new class.

        */

        if (head != (Cluster *) NULL)

          {

            cluster->next=(Cluster *) AcquireQuantumMemory(1,

              sizeof(*cluster->next));

            cluster=cluster->next;

          }

        else

          {

            cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));

            head=cluster;

          }

        if (cluster == (Cluster *) NULL)

          ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",

            image->filename);

        /*

          Initialize a new class.

        */

        (void) memset(cluster,0,sizeof(*cluster));

        cluster->red=red;

        cluster->green=green;

        cluster->blue=blue;

      }

    }

  }

  if (head == (Cluster *) NULL)

    {

      /*

        No classes were identified-- create one.

      */

      cluster=(Cluster *) AcquireQuantumMemory(1,sizeof(*cluster));

      if (cluster == (Cluster *) NULL)

        ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",

          image->filename);

      /*

        Initialize a new class.

      */

      (void) memset(cluster,0,sizeof(*cluster));

      cluster->red=red;

      cluster->green=green;

      cluster->blue=blue;

      head=cluster;

    }

  /*

    Count the pixels for each cluster.

  */

  status=MagickTrue;

  count=0;

  progress=0;

  image_view=AcquireVirtualCacheView(image,exception);

  for (y=0; y < (ssize_t) image->rows; y++)

  {

    const Quantum

      *p;

    ssize_t

      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);

    if (p == (const Quantum *) NULL)

      break;

    for (x=0; x < (ssize_t) image->columns; x++)

    {

      PixelInfo

        pixel;

      pixel.red=(double) ScaleQuantumToChar(GetPixelRed(image,p));

      pixel.green=(double) ScaleQuantumToChar(GetPixelGreen(image,p));

      pixel.blue=(double) ScaleQuantumToChar(GetPixelBlue(image,p));

      for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)

        if ((pixel.red >= (double) (cluster->red.left-SafeMargin)) &&

            (pixel.red <= (double) (cluster->red.right+SafeMargin)) &&

            (pixel.green >= (double) (cluster->green.left-SafeMargin)) &&

            (pixel.green <= (double) (cluster->green.right+SafeMargin)) &&

            (pixel.blue >= (double) (cluster->blue.left-SafeMargin)) &&

            (pixel.blue <= (double) (cluster->blue.right+SafeMargin)))

          {

            /*

              Count this pixel.

            */

            count++;

            cluster->red.center+=pixel.red;

            cluster->green.center+=pixel.green;

            cluster->blue.center+=pixel.blue;

            cluster->count++;

            break;

          }

      p+=(ptrdiff_t) GetPixelChannels(image);

    }

    if (image->progress_monitor != (MagickProgressMonitor) NULL)

      {

        MagickBooleanType

          proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)

        #pragma omp atomic

#endif

        progress++;

        proceed=SetImageProgress(image,SegmentImageTag,progress,2*image->rows);

        if (proceed == MagickFalse)

          status=MagickFalse;

      }

  }

  image_view=DestroyCacheView(image_view);

  /*

    Remove clusters that do not meet minimum cluster threshold.

  */

  count=0;

  last_cluster=head;

  next_cluster=head;

  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)

  {

    next_cluster=cluster->next;

    if ((cluster->count > 0) &&

        (cluster->count >= (count*cluster_threshold/100.0)))

      {

        /*

          Initialize cluster.

        */

        cluster->id=count;

        cluster->red.center/=cluster->count;

        cluster->green.center/=cluster->count;

        cluster->blue.center/=cluster->count;

        count++;

        last_cluster=cluster;

        continue;

      }

    /*

      Delete cluster.

    */

    if (cluster == head)

      head=next_cluster;

    else

      last_cluster->next=next_cluster;

    cluster=(Cluster *) RelinquishMagickMemory(cluster);

  }

  number_clusters=(size_t) count;

  if (verbose != MagickFalse)

    {

      /*

        Print cluster statistics.

      */

      (void) FormatLocaleFile(stdout,"Fuzzy C-means Statistics\n");

      (void) FormatLocaleFile(stdout,"===================\n\n");

      (void) FormatLocaleFile(stdout,"\tCluster Threshold = %g\n",(double)

        cluster_threshold);

      (void) FormatLocaleFile(stdout,"\tWeighting Exponent = %g\n",(double)

        weighting_exponent);

      (void) FormatLocaleFile(stdout,"\tTotal Number of Clusters = %.20g\n\n",

        (double) number_clusters);

      /*

        Print the total number of points per cluster.

      */

      (void) FormatLocaleFile(stdout,"\n\nNumber of Vectors Per Cluster\n");

      (void) FormatLocaleFile(stdout,"=============================\n\n");

      for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)

        (void) FormatLocaleFile(stdout,"Cluster #%.20g = %.20g\n",(double)

          cluster->id,(double) cluster->count);

      /*

        Print the cluster extents.

      */

      (void) FormatLocaleFile(stdout,

        "\n\n\nCluster Extents:        (Vector Size: %d)\n",MaxDimension);

      (void) FormatLocaleFile(stdout,"================");

      for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)

      {

        (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)

          cluster->id);

        (void) FormatLocaleFile(stdout,

          "%.20g-%.20g  %.20g-%.20g  %.20g-%.20g\n",(double)

          cluster->red.left,(double) cluster->red.right,(double)

          cluster->green.left,(double) cluster->green.right,(double)

          cluster->blue.left,(double) cluster->blue.right);

      }

      /*

        Print the cluster center values.

      */

      (void) FormatLocaleFile(stdout,

        "\n\n\nCluster Center Values:        (Vector Size: %d)\n",MaxDimension);

      (void) FormatLocaleFile(stdout,"=====================");

      for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)

      {

        (void) FormatLocaleFile(stdout,"\n\nCluster #%.20g\n\n",(double)

          cluster->id);

        (void) FormatLocaleFile(stdout,"%g  %g  %g\n",(double)

          cluster->red.center,(double) cluster->green.center,(double)

          cluster->blue.center);

      }

      (void) FormatLocaleFile(stdout,"\n");

    }

  if (number_clusters > 256)

    ThrowClassifyException(ImageError,"TooManyClusters",image->filename);

  /*

    Speed up distance calculations.

  */

  squares=(double *) AcquireQuantumMemory(513UL,sizeof(*squares));

  if (squares == (double *) NULL)

    ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",

      image->filename);

  squares+=255;

  for (i=(-255); i <= 255; i++)

    squares[i]=(double) i*(double) i;

  /*

    Allocate image colormap.

  */

  if (AcquireImageColormap(image,number_clusters,exception) == MagickFalse)

    ThrowClassifyException(ResourceLimitError,"MemoryAllocationFailed",

      image->filename);

  i=0;

  for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)

  {

    image->colormap[i].red=(double) ScaleCharToQuantum((unsigned char)

      (cluster->red.center+0.5));

    image->colormap[i].green=(double) ScaleCharToQuantum((unsigned char)

      (cluster->green.center+0.5));

    image->colormap[i].blue=(double) ScaleCharToQuantum((unsigned char)

      (cluster->blue.center+0.5));

    i++;

  }

  /*

    Do course grain classes.

  */

  image_view=AcquireAuthenticCacheView(image,exception);

#if defined(MAGICKCORE_OPENMP_SUPPORT)

  #pragma omp parallel for schedule(static) shared(progress,status) \

    magick_number_threads(image,image,image->rows,1)

#endif

  for (y=0; y < (ssize_t) image->rows; y++)

  {

    Cluster

      *c;

    const PixelInfo

      *magick_restrict p;

    ssize_t

      x;

    Quantum

      *magick_restrict q;

    if (status == MagickFalse)

      continue;

    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);

    if (q == (Quantum *) NULL)

      {

        status=MagickFalse;

        continue;

      }

    for (x=0; x < (ssize_t) image->columns; x++)

    {

      PixelInfo

        pixel;

      SetPixelIndex(image,(Quantum) 0,q);

      pixel.red=(double) ScaleQuantumToChar(GetPixelRed(image,q));

      pixel.green=(double) ScaleQuantumToChar(GetPixelGreen(image,q));

      pixel.blue=(double) ScaleQuantumToChar(GetPixelBlue(image,q));

      for (c=head; c != (Cluster *) NULL; c=c->next)

      {

        if ((pixel.red >= (double) (c->red.left-SafeMargin)) &&

            (pixel.red <= (double) (c->red.right+SafeMargin)) &&

            (pixel.green >= (double) (c->green.left-SafeMargin)) &&

            (pixel.green <= (double) (c->green.right+SafeMargin)) &&

            (pixel.blue >= (double) (c->blue.left-SafeMargin)) &&

            (pixel.blue <= (double) (c->blue.right+SafeMargin)))

          {

            /*

              Classify this pixel.

            */

            SetPixelIndex(image,(Quantum) c->id,q);

            break;

          }

      }

      if (c == (Cluster *) NULL)

        {

          double

            distance_squared,

            local_minima,

            numerator,

            ratio,

            sum;

          ssize_t

            j,

            k;

          /*

            Compute fuzzy membership.

          */

          local_minima=0.0;

          for (j=0; j < (ssize_t) image->colors; j++)

          {

            sum=0.0;

            p=image->colormap+j;

            distance_squared=

              squares[(ssize_t) (pixel.red-ScaleQuantumToChar((const Quantum) p->red))]+

              squares[(ssize_t) (pixel.green-ScaleQuantumToChar((const Quantum) p->green))]+

              squares[(ssize_t) (pixel.blue-ScaleQuantumToChar((const Quantum) p->blue))];

            numerator=distance_squared;

            for (k=0; k < (ssize_t) image->colors; k++)

            {

              p=image->colormap+k;

              distance_squared=

                squares[(ssize_t) (pixel.red-ScaleQuantumToChar((const Quantum) p->red))]+

                squares[(ssize_t) (pixel.green-ScaleQuantumToChar((const Quantum) p->green))]+

                squares[(ssize_t) (pixel.blue-ScaleQuantumToChar((const Quantum) p->blue))];

              ratio=numerator/distance_squared;

              sum+=SegmentPower(ratio);

            }

            if ((sum != 0.0) && ((1.0/sum) > local_minima))

              {

                /*

                  Classify this pixel.

                */

                local_minima=1.0/sum;

                SetPixelIndex(image,(Quantum) j,q);

              }

          }

        }

      q+=(ptrdiff_t) GetPixelChannels(image);

    }

    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)

      status=MagickFalse;

    if (image->progress_monitor != (MagickProgressMonitor) NULL)

      {

        MagickBooleanType

          proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)

        #pragma omp atomic

#endif

        progress++;

        proceed=SetImageProgress(image,SegmentImageTag,progress,2*image->rows);

        if (proceed == MagickFalse)

          status=MagickFalse;

      }

  }

  image_view=DestroyCacheView(image_view);

  status&=(MagickStatusType) SyncImage(image,exception);

  /*

    Relinquish resources.

  */

  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)

  {

    next_cluster=cluster->next;

    cluster=(Cluster *) RelinquishMagickMemory(cluster);

  }

  squares-=255;

  free_squares=squares;

  free_squares=(double *) RelinquishMagickMemory(free_squares);

  return(MagickTrue);

}

/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

+   C o n s o l i d a t e C r o s s i n g s                                   %

%                                                                             %

%                                                                             %

%                                                                             %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%  ConsolidateCrossings() guarantees that an even number of zero crossings

%  always lie between two crossings.

%

%  The format of the ConsolidateCrossings method is:

%

%      ConsolidateCrossings(ZeroCrossing *zero_crossing,

%        const size_t number_crossings)

%

%  A description of each parameter follows.

%

%    o zero_crossing: Specifies an array of structures of type ZeroCrossing.

%

%    o number_crossings: This size_t specifies the number of elements

%      in the zero_crossing array.

%

*/

static void ConsolidateCrossings(ZeroCrossing *zero_crossing,

  const size_t number_crossings)

{

  ssize_t

    i,

    j,

    k,

    l;

  ssize_t

    center,

    correct,

    count,

    left,

    right;

  /*

    Consolidate zero crossings.

  */

  for (i=(ssize_t) number_crossings-1; i >= 0; i--)

    for (j=0; j <= 255; j++)

    {

      if (zero_crossing[i].crossings[j] == 0)

        continue;

      /*

        Find the entry that is closest to j and still preserves the

        property that there are an even number of crossings between

        intervals.

      */

      for (k=j-1; k > 0; k--)

        if (zero_crossing[i+1].crossings[k] != 0)

          break;

      left=MagickMax(k,0);

      center=j;

      for (k=j+1; k < 255; k++)

        if (zero_crossing[i+1].crossings[k] != 0)

          break;

      right=MagickMin(k,255);

      /*

        K is the zero crossing just left of j.

      */

      for (k=j-1; k > 0; k--)

        if (zero_crossing[i].crossings[k] != 0)

          break;

      if (k < 0)

        k=0;

      /*

        Check center for an even number of crossings between k and j.

      */

      correct=(-1);

      if (zero_crossing[i+1].crossings[j] != 0)

        {

          count=0;

          for (l=k+1; l < center; l++)

            if (zero_crossing[i+1].crossings[l] != 0)

              count++;

          if (((count % 2) == 0) && (center != k))

            correct=center;

        }

      /*

        Check left for an even number of crossings between k and j.

      */

      if (correct == -1)

        {

          count=0;

          for (l=k+1; l < left; l++)

            if (zero_crossing[i+1].crossings[l] != 0)

              count++;

          if (((count % 2) == 0) && (left != k))

            correct=left;

        }

      /*

        Check right for an even number of crossings between k and j.

      */

      if (correct == -1)

        {

          count=0;

          for (l=k+1; l < right; l++)

            if (zero_crossing[i+1].crossings[l] != 0)

              count++;

          if (((count % 2) == 0) && (right != k))

            correct=right;

        }

      l=(ssize_t) zero_crossing[i].crossings[j];

      zero_crossing[i].crossings[j]=0;

      if (correct != -1)

        zero_crossing[i].crossings[correct]=(short) l;

    }

}

/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

+   D e f i n e R e g i o n                                                   %

%                                                                             %

%                                                                             %

%                                                                             %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%  DefineRegion() defines the left and right boundaries of a peak region.

%

%  The format of the DefineRegion method is:

%

%      ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)

%

%  A description of each parameter follows.

%

%    o extrema:  Specifies a pointer to an array of integers.  They

%      represent the peaks and valleys of the histogram for each color

%      component.

%

%    o extents:  This pointer to an ExtentPacket represent the extends

%      of a particular peak or valley of a color component.

%

*/

static ssize_t DefineRegion(const short *extrema,ExtentPacket *extents)

{

  /*

    Initialize to default values.

  */

  extents->left=0;

  extents->center=0.0;

  extents->right=255;

  /*

    Find the left side (maxima).

  */

  for ( ; extents->index <= 255; extents->index++)

    if (extrema[extents->index] > 0)

      break;

  if (extents->index > 255)

    return(MagickFalse);  /* no left side - no region exists */

  extents->left=extents->index;

  /*

    Find the right side (minima).

  */

  for ( ; extents->index <= 255; extents->index++)

    if (extrema[extents->index] < 0)

      break;

  extents->right=extents->index-1;

  return(MagickTrue);

}

/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

+   D e r i v a t i v e H i s t o g r a m                                     %

%                                                                             %

%                                                                             %

%                                                                             %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%  DerivativeHistogram() determines the derivative of the histogram using

%  central differencing.

%

%  The format of the DerivativeHistogram method is:

%

%      DerivativeHistogram(const double *histogram,

%        double *derivative)

%

%  A description of each parameter follows.

%

%    o histogram: Specifies an array of doubles representing the number

%      of pixels for each intensity of a particular color component.

%

%    o derivative: This array of doubles is initialized by

%      DerivativeHistogram to the derivative of the histogram using central

%      differencing.

%

*/

static void DerivativeHistogram(const double *histogram,

  double *derivative)

{

  ssize_t

    i,

    n;

  /*

    Compute endpoints using second order polynomial interpolation.

  */

  n=255;

  derivative[0]=(-1.5*histogram[0]+2.0*histogram[1]-0.5*histogram[2]);

  derivative[n]=(0.5*histogram[n-2]-2.0*histogram[n-1]+1.5*histogram[n]);

  /*

    Compute derivative using central differencing.

  */

  for (i=1; i < n; i++)

    derivative[i]=(histogram[i+1]-histogram[i-1])/2.0;

  return;

}

/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

+  G e t I m a g e D y n a m i c T h r e s h o l d                            %

%                                                                             %

%                                                                             %

%                                                                             %

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

%  GetImageDynamicThreshold() returns the dynamic threshold for an image.

%

%  The format of the GetImageDynamicThreshold method is:

%

%      MagickBooleanType GetImageDynamicThreshold(const Image *image,

%        const double cluster_threshold,const double smooth_threshold,

%        PixelInfo *pixel,ExceptionInfo *exception)

%

%  A description of each parameter follows.

%

%    o image: the image.

%

%    o cluster_threshold:  This double represents the minimum number of

%      pixels contained in a hexahedra before it can be considered valid

%      (expressed as a percentage).

%

%    o smooth_threshold: the smoothing threshold eliminates noise in the second

%      derivative of the histogram.  As the value is increased, you can expect a

%      smoother second derivative.

%

%    o pixel: return the dynamic threshold here.

%

%    o exception: return any errors or warnings in this structure.

%

*/

MagickExport MagickBooleanType GetImageDynamicThreshold(const Image *image,

  const double cluster_threshold,const double smooth_threshold,

  PixelInfo *pixel,ExceptionInfo *exception)

{

  Cluster

    *background,

    *cluster,

    *object,

    *head,

    *last_cluster,

    *next_cluster;

  ExtentPacket

    blue,

    green,

    red;

  MagickBooleanType