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
    proceed;

  double
    threshold;

  register const Quantum
    *p;

  register ssize_t
    i,
    x;

  short
    *extrema[MaxDimension];

  ssize_t
    count,
    *histogram[MaxDimension],
    y;

  /*
    Allocate histogram and extrema.
  */
  assert(image != (Image *) NULL);
  assert(image->signature == MagickCoreSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  GetPixelInfo(image,pixel);
  for (i=0; i < MaxDimension; i++)
  {
    histogram[i]=(ssize_t *) AcquireQuantumMemory(256UL,sizeof(**histogram));
    extrema[i]=(short *) AcquireQuantumMemory(256UL,sizeof(**histogram));
    if ((histogram[i] == (ssize_t *) NULL) || (extrema[i] == (short *) NULL))
      {
        for (i-- ; i >= 0; i--)
        {
          extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
          histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
        }
        (void) ThrowMagickException(exception,GetMagickModule(),
          ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
        return(MagickFalse);
      }
  }
  /*
    Initialize histogram.
  */
  InitializeHistogram(image,histogram,exception);
  (void) OptimalTau(histogram[Red],Tau,0.2f,DeltaTau,
    (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Red]);
  (void) OptimalTau(histogram[Green],Tau,0.2f,DeltaTau,
    (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Green]);
  (void) OptimalTau(histogram[Blue],Tau,0.2f,DeltaTau,
    (smooth_threshold == 0.0f ? 1.0f : smooth_threshold),extrema[Blue]);
  /*
    Form clusters.
  */
  cluster=(Cluster *) NULL;
  head=(Cluster *) 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 *) AcquireMagickMemory(
              sizeof(*cluster->next));
            cluster=cluster->next;
          }
        else
          {
            cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
            head=cluster;
          }
        if (cluster == (Cluster *) NULL)
          {
            (void) ThrowMagickException(exception,GetMagickModule(),
              ResourceLimitError,"MemoryAllocationFailed","`%s'",
              image->filename);
            return(MagickFalse);
          }
        /*
          Initialize a new class.
        */
        cluster->count=0;
        cluster->red=red;
        cluster->green=green;
        cluster->blue=blue;
        cluster->next=(Cluster *) NULL;
      }
    }
  }
  if (head == (Cluster *) NULL)
    {
      /*
        No classes were identified-- create one.
      */
      cluster=(Cluster *) AcquireMagickMemory(sizeof(*cluster));
      if (cluster == (Cluster *) NULL)
        {
          (void) ThrowMagickException(exception,GetMagickModule(),
            ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
          return(MagickFalse);
        }
      /*
        Initialize a new class.
      */
      cluster->count=0;
      cluster->red=red;
      cluster->green=green;
      cluster->blue=blue;
      cluster->next=(Cluster *) NULL;
      head=cluster;
    }
  /*
    Count the pixels for each cluster.
  */
  count=0;
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    p=GetVirtualPixels(image,0,y,image->columns,1,exception);
    if (p == (const Quantum *) NULL)
      break;
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      for (cluster=head; cluster != (Cluster *) NULL; cluster=cluster->next)
        if (((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) >=
             (cluster->red.left-SafeMargin)) &&
            ((ssize_t) ScaleQuantumToChar(GetPixelRed(image,p)) <=
             (cluster->red.right+SafeMargin)) &&
            ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) >=
             (cluster->green.left-SafeMargin)) &&
            ((ssize_t) ScaleQuantumToChar(GetPixelGreen(image,p)) <=
             (cluster->green.right+SafeMargin)) &&
            ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) >=
             (cluster->blue.left-SafeMargin)) &&
            ((ssize_t) ScaleQuantumToChar(GetPixelBlue(image,p)) <=
             (cluster->blue.right+SafeMargin)))
          {
            /*
              Count this pixel.
            */
            count++;
            cluster->red.center+=(double) ScaleQuantumToChar(
              GetPixelRed(image,p));
            cluster->green.center+=(double) ScaleQuantumToChar(
              GetPixelGreen(image,p));
            cluster->blue.center+=(double) ScaleQuantumToChar(
              GetPixelBlue(image,p));
            cluster->count++;
            break;
          }
      p+=GetPixelChannels(image);
    }
    proceed=SetImageProgress(image,SegmentImageTag,(MagickOffsetType) y,
      2*image->rows);
    if (proceed == MagickFalse)
      break;
  }
  /*
    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);
  }
  object=head;
  background=head;
  if (count > 1)
    {
      object=head->next;
      for (cluster=object; cluster->next != (Cluster *) NULL; )
      {
        if (cluster->count < object->count)
          object=cluster;
        cluster=cluster->next;
      }
      background=head->next;
      for (cluster=background; cluster->next != (Cluster *) NULL; )
      {
        if (cluster->count > background->count)
          background=cluster;
        cluster=cluster->next;
      }
    }
  if (background != (Cluster *) NULL)
    {
      threshold=(background->red.center+object->red.center)/2.0;
      pixel->red=(double) ScaleCharToQuantum((unsigned char)
        (threshold+0.5));
      threshold=(background->green.center+object->green.center)/2.0;
      pixel->green=(double) ScaleCharToQuantum((unsigned char)
        (threshold+0.5));
      threshold=(background->blue.center+object->blue.center)/2.0;
      pixel->blue=(double) ScaleCharToQuantum((unsigned char)
        (threshold+0.5));
    }
  /*
    Relinquish resources.
  */
  for (cluster=head; cluster != (Cluster *) NULL; cluster=next_cluster)
  {
    next_cluster=cluster->next;
    cluster=(Cluster *) RelinquishMagickMemory(cluster);
  }
  for (i=0; i < MaxDimension; i++)
  {
    extrema[i]=(short *) RelinquishMagickMemory(extrema[i]);
    histogram[i]=(ssize_t *) RelinquishMagickMemory(histogram[i]);
  }
  return(MagickTrue);
}