/*

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

%                                                                             %

%                                                                             %

%                                                                             %

%       TTTTT  H   H  RRRR   EEEEE  SSSSS  H   H   OOO   L      DDDD          %

%         T    H   H  R   R  E      SS     H   H  O   O  L      D   D         %

%         T    HHHHH  RRRR   EEE     SSS   HHHHH  O   O  L      D   D         %

%         T    H   H  R R    E         SS  H   H  O   O  L      D   D         %

%         T    H   H  R  R   EEEEE  SSSSS  H   H   OOO   LLLLL  DDDD          %

%                                                                             %

%                                                                             %

%                      MagickCore Image Threshold Methods                     %

%                                                                             %

%                               Software Design                               %

%                                    Cristy                                   %

%                                 October 1996                                %

%                                                                             %

%                                                                             %

%  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/script/license.php                               %

%                                                                             %

%  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.                                             %

%                                                                             %

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

%

%

%

*/

/*

  Include declarations.

*/

#include "MagickCore/studio.h"

#include "MagickCore/artifact.h"

#include "MagickCore/blob.h"

#include "MagickCore/cache-view.h"

#include "MagickCore/color.h"

#include "MagickCore/color-private.h"

#include "MagickCore/colormap.h"

#include "MagickCore/colorspace.h"

#include "MagickCore/colorspace-private.h"

#include "MagickCore/configure.h"

#include "MagickCore/constitute.h"

#include "MagickCore/decorate.h"

#include "MagickCore/draw.h"

#include "MagickCore/enhance.h"

#include "MagickCore/exception.h"

#include "MagickCore/exception-private.h"

#include "MagickCore/effect.h"

#include "MagickCore/fx.h"

#include "MagickCore/gem.h"

#include "MagickCore/gem-private.h"

#include "MagickCore/geometry.h"

#include "MagickCore/image-private.h"

#include "MagickCore/list.h"

#include "MagickCore/log.h"

#include "MagickCore/memory_.h"

#include "MagickCore/monitor.h"

#include "MagickCore/monitor-private.h"

#include "MagickCore/montage.h"

#include "MagickCore/option.h"

#include "MagickCore/pixel-accessor.h"

#include "MagickCore/property.h"

#include "MagickCore/quantize.h"

#include "MagickCore/quantum.h"

#include "MagickCore/quantum-private.h"

#include "MagickCore/random_.h"

#include "MagickCore/random-private.h"

#include "MagickCore/resize.h"

#include "MagickCore/resource_.h"

#include "MagickCore/segment.h"

#include "MagickCore/shear.h"

#include "MagickCore/signature-private.h"

#include "MagickCore/string_.h"

#include "MagickCore/string-private.h"

#include "MagickCore/thread-private.h"

#include "MagickCore/threshold.h"

#include "MagickCore/token.h"

#include "MagickCore/transform.h"

#include "MagickCore/xml-tree.h"

#include "MagickCore/xml-tree-private.h"

/*

  Define declarations.

*/

#define ThresholdsFilename  "thresholds.xml"

/*

  Typedef declarations.

*/

struct _ThresholdMap

{

  char

    *map_id,

    *description;

  size_t

    width,

    height;

  ssize_t

    divisor,

    *levels;

};

/*

  Static declarations.

*/

#if MAGICKCORE_ZERO_CONFIGURATION_SUPPORT

  #include "MagickCore/threshold-map.h"

#else

static const char *const

  BuiltinMap=

    "<?xml version=\"1.0\"?>"

    "<thresholds>"

    "  <threshold map=\"threshold\" alias=\"1x1\">"

    "    <description>Threshold 1x1 (non-dither)</description>"

    "    <levels width=\"1\" height=\"1\" divisor=\"2\">"

    "        1"

    "    </levels>"

    "  </threshold>"

    "  <threshold map=\"checks\" alias=\"2x1\">"

    "    <description>Checkerboard 2x1 (dither)</description>"

    "    <levels width=\"2\" height=\"2\" divisor=\"3\">"

    "       1 2"

    "       2 1"

    "    </levels>"

    "  </threshold>"

    "</thresholds>";

#endif

/*

  Forward declarations.

*/

static ThresholdMap

  *GetThresholdMapFile(const char *,const char *,const char *,ExceptionInfo *);

/*

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

%                                                                             %

%                                                                             %

%                                                                             %

%     A d a p t i v e T h r e s h o l d I m a g e                             %

%                                                                             %

%                                                                             %

%                                                                             %

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

%

%  AdaptiveThresholdImage() selects an individual threshold for each pixel

%  based on the range of intensity values in its local neighborhood.  This

%  allows for thresholding of an image whose global intensity histogram

%  doesn't contain distinctive peaks.

%

%  The format of the AdaptiveThresholdImage method is:

%

%      Image *AdaptiveThresholdImage(const Image *image,const size_t width,

%        const size_t height,const double bias,ExceptionInfo *exception)

%

%  A description of each parameter follows:

%

%    o image: the image.

%

%    o width: the width of the local neighborhood.

%

%    o height: the height of the local neighborhood.

%

%    o bias: the mean bias.

%

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

%

*/

MagickExport Image *AdaptiveThresholdImage(const Image *image,

  const size_t width,const size_t height,const double bias,

  ExceptionInfo *exception)

{

#define AdaptiveThresholdImageTag  "AdaptiveThreshold/Image"

  CacheView

    *image_view,

    *threshold_view;

  Image

    *threshold_image;

  MagickBooleanType

    status;

  MagickOffsetType

    progress;

  MagickSizeType

    number_pixels;

  ssize_t

    y;

  /*

    Initialize threshold image attributes.

  */

  assert(image != (Image *) NULL);

  assert(image->signature == MagickCoreSignature);

  assert(exception != (ExceptionInfo *) NULL);

  assert(exception->signature == MagickCoreSignature);

  if (IsEventLogging() != MagickFalse)

    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);

  threshold_image=CloneImage(image,0,0,MagickTrue,exception);

  if (threshold_image == (Image *) NULL)

    return((Image *) NULL);

  if ((width == 0) || (height == 0))

    return(threshold_image);

  status=SetImageStorageClass(threshold_image,DirectClass,exception);

  if (status == MagickFalse)

    {

      threshold_image=DestroyImage(threshold_image);

      return((Image *) NULL);

    }

  /*

    Threshold image.

  */

  status=MagickTrue;

  progress=0;

  number_pixels=(MagickSizeType) width*height;

  image_view=AcquireVirtualCacheView(image,exception);

  threshold_view=AcquireAuthenticCacheView(threshold_image,exception);

#if defined(MAGICKCORE_OPENMP_SUPPORT)

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

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

#endif

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

  {

    double

      channel_bias[MaxPixelChannels],

      channel_sum[MaxPixelChannels];

    const Quantum

      *magick_restrict p,

      *magick_restrict pixels;

    Quantum

      *magick_restrict q;

    ssize_t

      center,

      i,

      u,

      v,

      x;

    if (status == MagickFalse)

      continue;

    p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)

      (height/2L),image->columns+width,height,exception);

    q=QueueCacheViewAuthenticPixels(threshold_view,0,y,threshold_image->columns,

      1,exception);

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

      {

        status=MagickFalse;

        continue;

      }

    center=(ssize_t) GetPixelChannels(image)*((ssize_t) image->columns+

      (ssize_t) width)*((ssize_t) height/2L)+(ssize_t) GetPixelChannels(image)*

      ((ssize_t) width/2);

    for (i=0; i < (ssize_t) GetPixelChannels(image); i++)

    {

      PixelChannel channel = GetPixelChannelChannel(image,i);

      PixelTrait traits = GetPixelChannelTraits(image,channel);

      PixelTrait threshold_traits=GetPixelChannelTraits(threshold_image,

        channel);

      if ((traits == UndefinedPixelTrait) ||

          (threshold_traits == UndefinedPixelTrait))

        continue;

      if ((threshold_traits & CopyPixelTrait) != 0)

        {

          SetPixelChannel(threshold_image,channel,p[center+i],q);

          continue;

        }

      pixels=p;

      channel_bias[channel]=0.0;

      channel_sum[channel]=0.0;

      for (v=0; v < (ssize_t) height; v++)

      {

        for (u=0; u < (ssize_t) width; u++)

        {

          if (u == (ssize_t) (width-1))

            channel_bias[channel]+=(double) pixels[i];

          channel_sum[channel]+=(double) pixels[i];

          pixels+=(ptrdiff_t) GetPixelChannels(image);

        }

        pixels+=(ptrdiff_t) GetPixelChannels(image)*image->columns;

      }

    }

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

    {

      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)

      {

        double

          mean;

        PixelChannel channel = GetPixelChannelChannel(image,i);

        PixelTrait traits = GetPixelChannelTraits(image,channel);

        PixelTrait threshold_traits=GetPixelChannelTraits(threshold_image,

          channel);

        if ((traits == UndefinedPixelTrait) ||

            (threshold_traits == UndefinedPixelTrait))

          continue;

        if ((threshold_traits & CopyPixelTrait) != 0)

          {

            SetPixelChannel(threshold_image,channel,p[center+i],q);

            continue;

          }

        channel_sum[channel]-=channel_bias[channel];

        channel_bias[channel]=0.0;

        pixels=p;

        for (v=0; v < (ssize_t) height; v++)

        {

          channel_bias[channel]+=(double) pixels[i];

          pixels+=(width-1)*GetPixelChannels(image);

          channel_sum[channel]+=(double) pixels[i];

          pixels+=(ptrdiff_t) GetPixelChannels(image)*(image->columns+1);

        }

        mean=(double) (channel_sum[channel]/number_pixels+bias);

        SetPixelChannel(threshold_image,channel,(Quantum) ((double)

          p[center+i] <= mean ? 0 : QuantumRange),q);

      }

      p+=(ptrdiff_t) GetPixelChannels(image);

      q+=(ptrdiff_t) GetPixelChannels(threshold_image);

    }

    if (SyncCacheViewAuthenticPixels(threshold_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,AdaptiveThresholdImageTag,progress,

          image->rows);

        if (proceed == MagickFalse)

          status=MagickFalse;

      }

  }

  threshold_image->type=image->type;

  threshold_view=DestroyCacheView(threshold_view);

  image_view=DestroyCacheView(image_view);

  if (status == MagickFalse)

    threshold_image=DestroyImage(threshold_image);

  return(threshold_image);

}

/*

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

%                                                                             %

%                                                                             %

%                                                                             %

%     A u t o T h r e s h o l d I m a g e                                     %

%                                                                             %

%                                                                             %

%                                                                             %

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

%

%  AutoThresholdImage()  automatically performs image thresholding

%  dependent on which method you specify.

%

%  The format of the AutoThresholdImage method is:

%

%      MagickBooleanType AutoThresholdImage(Image *image,

%        const AutoThresholdMethod method,ExceptionInfo *exception)

%

%  A description of each parameter follows:

%

%    o image: The image to auto-threshold.

%

%    o method: choose from Kapur, OTSU, or Triangle.

%

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

%

*/

static double KapurThreshold(const Image *image,const double *histogram,

  ExceptionInfo *exception)

{

#define MaxIntensity  255

  double

    *black_entropy,

    *cumulative_histogram,

    entropy,

    epsilon,

    maximum_entropy,

    *white_entropy;

  ssize_t

    i,

    j;

  size_t

    threshold;

  /*

    Compute optimal threshold from the entropy of the histogram.

  */

  cumulative_histogram=(double *) AcquireQuantumMemory(MaxIntensity+1UL,

    sizeof(*cumulative_histogram));

  black_entropy=(double *) AcquireQuantumMemory(MaxIntensity+1UL,

    sizeof(*black_entropy));

  white_entropy=(double *) AcquireQuantumMemory(MaxIntensity+1UL,

    sizeof(*white_entropy));

  if ((cumulative_histogram == (double *) NULL) ||

      (black_entropy == (double *) NULL) || (white_entropy == (double *) NULL))

    {

      if (white_entropy != (double *) NULL)

        white_entropy=(double *) RelinquishMagickMemory(white_entropy);

      if (black_entropy != (double *) NULL)

        black_entropy=(double *) RelinquishMagickMemory(black_entropy);

      if (cumulative_histogram != (double *) NULL)

        cumulative_histogram=(double *)

          RelinquishMagickMemory(cumulative_histogram);

      (void) ThrowMagickException(exception,GetMagickModule(),

        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);

      return(-1.0);

    }

   /*

     Entropy for black and white parts of the histogram.

   */

   cumulative_histogram[0]=histogram[0];

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

     cumulative_histogram[i]=cumulative_histogram[i-1]+histogram[i];

   epsilon=MagickMinimumValue;

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

   {

     /*

       Black entropy.

     */

     black_entropy[j]=0.0;

     if (cumulative_histogram[j] > epsilon)

       {

         entropy=0.0;

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

           if (histogram[i] > epsilon)

             entropy-=histogram[i]/cumulative_histogram[j]*

               log(histogram[i]/cumulative_histogram[j]);

         black_entropy[j]=entropy;

       }

     /*

       White entropy.

     */

     white_entropy[j]=0.0;

     if ((1.0-cumulative_histogram[j]) > epsilon)

       {

         entropy=0.0;

         for (i=j+1; i <= MaxIntensity; i++)

           if (histogram[i] > epsilon)

             entropy-=histogram[i]/(1.0-cumulative_histogram[j])*

               log(histogram[i]/(1.0-cumulative_histogram[j]));

         white_entropy[j]=entropy;

       }

   }

  /*

    Find histogram bin with maximum entropy.

  */

  maximum_entropy=black_entropy[0]+white_entropy[0];

  threshold=0;

  for (j=1; j <= MaxIntensity; j++)

    if ((black_entropy[j]+white_entropy[j]) > maximum_entropy)

      {

        maximum_entropy=black_entropy[j]+white_entropy[j];

        threshold=(size_t) j;

      }

  /*

    Free resources.

  */

  white_entropy=(double *) RelinquishMagickMemory(white_entropy);

  black_entropy=(double *) RelinquishMagickMemory(black_entropy);

  cumulative_histogram=(double *) RelinquishMagickMemory(cumulative_histogram);

  return(100.0*threshold/MaxIntensity);

}

static double OTSUThreshold(const Image *image,const double *histogram,

  ExceptionInfo *exception)

{

  double

    max_sigma,

    *myu,

    *omega,

    *probability,

    *sigma,

    threshold;

  ssize_t

    i;

  /*

    Compute optimal threshold from maximization of inter-class variance.

  */

  myu=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*myu));

  omega=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*omega));

  probability=(double *) AcquireQuantumMemory(MaxIntensity+1UL,

    sizeof(*probability));

  sigma=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*sigma));

  if ((myu == (double *) NULL) || (omega == (double *) NULL) ||

      (probability == (double *) NULL) || (sigma == (double *) NULL))

    {

      if (sigma != (double *) NULL)

        sigma=(double *) RelinquishMagickMemory(sigma);

      if (probability != (double *) NULL)

        probability=(double *) RelinquishMagickMemory(probability);

      if (omega != (double *) NULL)

        omega=(double *) RelinquishMagickMemory(omega);

      if (myu != (double *) NULL)

        myu=(double *) RelinquishMagickMemory(myu);

      (void) ThrowMagickException(exception,GetMagickModule(),

        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);

      return(-1.0);

    }

  /*

    Calculate probability density.

  */

  for (i=0; i <= (ssize_t) MaxIntensity; i++)

    probability[i]=histogram[i];

  /*

    Generate probability of graylevels and mean value for separation.

  */

  omega[0]=probability[0];

  myu[0]=0.0;

  for (i=1; i <= (ssize_t) MaxIntensity; i++)

  {

    omega[i]=omega[i-1]+probability[i];

    myu[i]=myu[i-1]+i*probability[i];

  }

  /*

    Sigma maximization: inter-class variance and compute optimal threshold.

  */

  threshold=0;

  max_sigma=0.0;

  for (i=0; i < (ssize_t) MaxIntensity; i++)

  {

    sigma[i]=0.0;

    if ((omega[i] != 0.0) && (omega[i] != 1.0))

      sigma[i]=pow(myu[MaxIntensity]*omega[i]-myu[i],2.0)/(omega[i]*(1.0-

        omega[i]));

    if (sigma[i] > max_sigma)

      {

        max_sigma=sigma[i];

        threshold=(double) i;

      }

  }

  /*

    Free resources.

  */

  myu=(double *) RelinquishMagickMemory(myu);

  omega=(double *) RelinquishMagickMemory(omega);

  probability=(double *) RelinquishMagickMemory(probability);

  sigma=(double *) RelinquishMagickMemory(sigma);

  return(100.0*threshold/MaxIntensity);

}

static double TriangleThreshold(const double *histogram)

{

  double

    a,

    b,

    c,

    count,

    distance,

    inverse_ratio,

    max_distance,

    segment,

    x1,

    x2,

    y1,

    y2;

  ssize_t

    i;

  ssize_t

    end,

    max,

    start,

    threshold;

  /*

    Compute optimal threshold with triangle algorithm.

  */

  start=0;  /* find start bin, first bin not zero count */

  for (i=0; i <= (ssize_t) MaxIntensity; i++)

    if (histogram[i] > 0.0)

      {

        start=i;

        break;

      }

  end=0;  /* find end bin, last bin not zero count */

  for (i=(ssize_t) MaxIntensity; i >= 0; i--)

    if (histogram[i] > 0.0)

      {

        end=i;

        break;

      }

  max=0;  /* find max bin, bin with largest count */

  count=0.0;

  for (i=0; i <= (ssize_t) MaxIntensity; i++)

    if (histogram[i] > count)

      {

        max=i;

        count=histogram[i];

      }

  /*

    Compute threshold at split point.

  */

  x1=(double) max;

  y1=histogram[max];

  x2=(double) end;

  if ((max-start) >= (end-max))

    x2=(double) start;

  y2=0.0;

  a=y1-y2;

  b=x2-x1;

  c=(-1.0)*(a*x1+b*y1);

  inverse_ratio=1.0/sqrt(a*a+b*b+c*c);

  threshold=0;

  max_distance=0.0;

  if (x2 == (double) start)

    for (i=start; i < max; i++)

    {

      segment=inverse_ratio*(a*i+b*histogram[i]+c);

      distance=sqrt(segment*segment);

      if ((distance > max_distance) && (segment > 0.0))

        {

          threshold=i;

          max_distance=distance;

        }

    }

  else

    for (i=end; i > max; i--)

    {

      segment=inverse_ratio*(a*i+b*histogram[i]+c);

      distance=sqrt(segment*segment);

      if ((distance > max_distance) && (segment < 0.0))

        {

          threshold=i;

          max_distance=distance;

        }

    }

  return(100.0*threshold/MaxIntensity);

}

MagickExport MagickBooleanType AutoThresholdImage(Image *image,

  const AutoThresholdMethod method,ExceptionInfo *exception)

{

  CacheView

    *image_view;

  char

    property[MagickPathExtent];

  double

    gamma,

    *histogram,

    sum,

    threshold;

  MagickBooleanType

    status;

  ssize_t

    i;

  ssize_t

    y;

  /*

    Form histogram.

  */

  assert(image != (Image *) NULL);

  assert(image->signature == MagickCoreSignature);

  if (IsEventLogging() != MagickFalse)

    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);

  histogram=(double *) AcquireQuantumMemory(MaxIntensity+1UL,

    sizeof(*histogram));

  if (histogram == (double *) NULL)

    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",

      image->filename);

  status=MagickTrue;

  (void) memset(histogram,0,(MaxIntensity+1UL)*sizeof(*histogram));

  image_view=AcquireVirtualCacheView(image,exception);

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

  {

    const Quantum

      *magick_restrict 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++)

    {

      double intensity = GetPixelIntensity(image,p);

      histogram[ScaleQuantumToChar(ClampToQuantum(intensity))]++;

      p+=(ptrdiff_t) GetPixelChannels(image);

    }

  }

  image_view=DestroyCacheView(image_view);

  /*

    Normalize histogram.

  */

  sum=0.0;

  for (i=0; i <= (ssize_t) MaxIntensity; i++)

    sum+=histogram[i];

  gamma=MagickSafeReciprocal(sum);

  for (i=0; i <= (ssize_t) MaxIntensity; i++)

    histogram[i]=gamma*histogram[i];

  /*

    Discover threshold from histogram.

  */

  switch (method)

  {

    case KapurThresholdMethod:

    {

      threshold=KapurThreshold(image,histogram,exception);

      break;

    }

    case OTSUThresholdMethod:

    default:

    {

      threshold=OTSUThreshold(image,histogram,exception);

      break;

    }

    case TriangleThresholdMethod:

    {

      threshold=TriangleThreshold(histogram);

      break;

    }

  }

  histogram=(double *) RelinquishMagickMemory(histogram);

  if (threshold < 0.0)

    status=MagickFalse;

  if (status == MagickFalse)

    return(MagickFalse);

  /*

    Threshold image.

  */

  (void) FormatLocaleString(property,MagickPathExtent,"%g%%",threshold);

  (void) SetImageProperty(image,"auto-threshold:threshold",property,exception);

  if (IsStringTrue(GetImageArtifact(image,"auto-threshold:verbose")) != MagickFalse)

    (void) FormatLocaleFile(stdout,"%.*g%%\n",GetMagickPrecision(),threshold);

  return(BilevelImage(image,(double) QuantumRange*threshold/100.0,exception));

}

/*

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

%                                                                             %

%                                                                             %

%                                                                             %

%     B i l e v e l I m a g e                                                 %

%                                                                             %

%                                                                             %

%                                                                             %

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

%

%  BilevelImage() changes the value of individual pixels based on the

%  intensity of each pixel channel.  The result is a high-contrast image.

%

%  More precisely each channel value of the image is 'thresholded' so that if

%  it is equal to or less than the given value it is set to zero, while any

%  value greater than that give is set to it maximum or QuantumRange.

%

%  This function is what is used to implement the "-threshold" operator for

%  the command line API.

%

%  If the default channel setting is given the image is thresholded using just

%  the gray 'intensity' of the image, rather than the individual channels.

%

%  The format of the BilevelImage method is:

%

%      MagickBooleanType BilevelImage(Image *image,const double threshold,

%        ExceptionInfo *exception)

%

%  A description of each parameter follows:

%

%    o image: the image.

%

%    o threshold: define the threshold values.

%

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

%

%  Aside: You can get the same results as operator using LevelImages()

%  with the 'threshold' value for both the black_point and the white_point.

%

*/

MagickExport MagickBooleanType BilevelImage(Image *image,const double threshold,

  ExceptionInfo *exception)

{

#define ThresholdImageTag  "Threshold/Image"

  CacheView

    *image_view;

  MagickBooleanType

    status;

  MagickOffsetType

    progress;

  ssize_t

    y;

  assert(image != (Image *) NULL);

  assert(image->signature == MagickCoreSignature);

  if (IsEventLogging() != MagickFalse)

    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);

  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)

    return(MagickFalse);

  if (IsGrayColorspace(image->colorspace) == MagickFalse)

    (void) SetImageColorspace(image,sRGBColorspace,exception);

  /*

    Bilevel threshold image.

  */

  status=MagickTrue;

  progress=0;

  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,2)

#endif

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

  {

    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++)

    {

      double

        pixel;

      ssize_t

        i;

      pixel=GetPixelIntensity(image,q);

      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)

      {

        PixelChannel channel = GetPixelChannelChannel(image,i);

        PixelTrait traits = GetPixelChannelTraits(image,channel);

        if ((traits & UpdatePixelTrait) == 0)

          continue;

        if (image->channel_mask != AllChannels)

          pixel=(double) q[i];

        q[i]=(Quantum) (pixel <= threshold ? 0 : QuantumRange);

      }

      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,ThresholdImageTag,progress++,

          image->rows);

        if (proceed == MagickFalse)

          status=MagickFalse;

      }

  }

  image_view=DestroyCacheView(image_view);

  return(status);

}

/*

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

%                                                                             %

%                                                                             %

%                                                                             %

%     B l a c k T h r e s h o l d I m a g e                                   %

%                                                                             %

%                                                                             %

%                                                                             %

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

%

%  BlackThresholdImage() is like ThresholdImage() but forces all pixels below

%  the threshold into black while leaving all pixels at or above the threshold

%  unchanged.

%

%  The format of the BlackThresholdImage method is:

%

%      MagickBooleanType BlackThresholdImage(Image *image,

%        const char *threshold,ExceptionInfo *exception)

%

%  A description of each parameter follows:

%

%    o image: the image.

%

%    o threshold: define the threshold value.

%

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

%

*/

MagickExport MagickBooleanType BlackThresholdImage(Image *image,

  const char *thresholds,ExceptionInfo *exception)

{

#define ThresholdImageTag  "Threshold/Image"

  CacheView

    *image_view;