]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/threshold.c
(no commit message)
[imagemagick] / MagickCore / threshold.c
index af85140139b1b3e50420c68091aa1508e1566021..e7f1103bd77590e4750febc61c3021f2f973555c 100644 (file)
 %                      MagickCore Image Threshold Methods                     %
 %                                                                             %
 %                               Software Design                               %
-%                                 John Cristy                                 %
+%                                    Cristy                                   %
 %                                 October 1996                                %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2014 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  %
@@ -48,6 +48,7 @@
 #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"
@@ -109,6 +110,28 @@ struct _ThresholdMap
     *levels;
 };
 \f
+/*
+  Static declarations.
+*/
+static const char
+  *MinimalThresholdMap =
+    "<?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>";
+\f
 /*
   Forward declarations.
 */
@@ -183,8 +206,6 @@ MagickExport Image *AdaptiveThresholdImage(const Image *image,
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
-  if ((width % 2) == 0)
-    ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
   threshold_image=CloneImage(image,image->columns,image->rows,MagickTrue,
     exception);
   if (threshold_image == (Image *) NULL)
@@ -201,24 +222,33 @@ MagickExport Image *AdaptiveThresholdImage(const Image *image,
   status=MagickTrue;
   progress=0;
   number_pixels=(MagickSizeType) width*height;
-  image_view=AcquireCacheView(image);
-  threshold_view=AcquireCacheView(threshold_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  threshold_view=AcquireAuthenticCacheView(threshold_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    magick_threads(image,threshold_image,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
+    double
+      channel_bias[MaxPixelChannels],
+      channel_sum[MaxPixelChannels];
+
     register const Quantum
-      *restrict p;
+      *restrict p,
+      *restrict pixels;
 
     register Quantum
       *restrict q;
 
     register ssize_t
+      i,
       x;
 
     ssize_t
-      center;
+      center,
+      u,
+      v;
 
     if (status == MagickFalse)
       continue;
@@ -233,57 +263,68 @@ MagickExport Image *AdaptiveThresholdImage(const Image *image,
       }
     center=(ssize_t) GetPixelChannels(image)*(image->columns+width)*(height/2L)+
       GetPixelChannels(image)*(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) ||
+          (GetPixelReadMask(image,p) == 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]+=pixels[i];
+          channel_sum[channel]+=pixels[i];
+          pixels+=GetPixelChannels(image);
+        }
+        pixels+=(image->columns-1)*GetPixelChannels(image);
+      }
+    }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
-      register ssize_t
-        i;
-
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
-          mean,
-          pixel;
-
-        PixelChannel
-          channel;
-
-        PixelTrait
-          threshold_traits,
-          traits;
-
-        register const Quantum
-          *restrict pixels;
-
-        register ssize_t
-          u;
-
-        ssize_t
-          v;
+        double
+          mean;
 
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
-        threshold_traits=GetPixelChannelMapTraits(threshold_image,channel);
+        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)
+        if (((threshold_traits & CopyPixelTrait) != 0) ||
+            (GetPixelReadMask(image,p) == 0))
           {
             SetPixelChannel(threshold_image,channel,p[center+i],q);
             continue;
           }
+        channel_sum[channel]-=channel_bias[channel];
+        channel_bias[channel]=0.0;
         pixels=p;
-        pixel=0.0;
         for (v=0; v < (ssize_t) height; v++)
         {
-          for (u=0; u < (ssize_t) width; u++)
-          {
-            pixel+=pixels[i];
-            pixels+=GetPixelChannels(image);
-          }
-          pixels+=image->columns*GetPixelChannels(image);
+          channel_bias[channel]+=pixels[i];
+          pixels+=(width-1)*GetPixelChannels(image);
+          channel_sum[channel]+=pixels[i];
+          pixels+=(image->columns-1)*GetPixelChannels(image);
         }
-        mean=(MagickRealType) (pixel/number_pixels+bias);
-        SetPixelChannel(threshold_image,channel,(Quantum) ((MagickRealType)
+        mean=(double) (channel_sum[channel]/number_pixels+bias);
+        SetPixelChannel(threshold_image,channel,(Quantum) ((double)
           p[center+i] <= mean ? 0 : QuantumRange),q);
       }
       p+=GetPixelChannels(image);
@@ -297,7 +338,7 @@ MagickExport Image *AdaptiveThresholdImage(const Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_AdaptiveThresholdImage)
+        #pragma omp critical (MagickCore_AdaptiveThresholdImage)
 #endif
         proceed=SetImageProgress(image,AdaptiveThresholdImageTag,progress++,
           image->rows);
@@ -377,14 +418,17 @@ MagickExport MagickBooleanType BilevelImage(Image *image,const double threshold,
     (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=AcquireCacheView(image);
+  image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,8) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    magick_threads(image,image,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -404,22 +448,27 @@ MagickExport MagickBooleanType BilevelImage(Image *image,const double threshold,
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        pixel;
+
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,q) == 0)
+        {
+          q+=GetPixelChannels(image);
+          continue;
+        }
+      pixel=GetPixelIntensity(image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          traits;
-
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
-        if ((traits & UpdatePixelTrait) != 0)
-          q[i]=(Quantum) ((MagickRealType) q[i] <= threshold ? 0 :
-            QuantumRange);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        if ((traits & UpdatePixelTrait) == 0)
+          continue;
+        if (image->channel_mask != DefaultChannels)
+          pixel=(double) q[i];
+        q[i]=(Quantum) (pixel <= threshold ? 0 : QuantumRange);
       }
       q+=GetPixelChannels(image);
     }
@@ -431,7 +480,7 @@ MagickExport MagickBooleanType BilevelImage(Image *image,const double threshold,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_BilevelImage)
+        #pragma omp critical (MagickCore_BilevelImage)
 #endif
         proceed=SetImageProgress(image,ThresholdImageTag,progress++,
           image->rows);
@@ -489,15 +538,12 @@ MagickExport MagickBooleanType BlackThresholdImage(Image *image,
   MagickOffsetType
     progress;
 
-  MagickRealType
-    threshold[5];
+  PixelInfo
+    threshold;
 
   MagickStatusType
     flags;
 
-  register ssize_t
-    i;
-
   ssize_t
     y;
 
@@ -509,28 +555,45 @@ MagickExport MagickBooleanType BlackThresholdImage(Image *image,
     return(MagickTrue);
   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
     return(MagickFalse);
+  if (IsGrayColorspace(image->colorspace) != MagickFalse)
+    (void) SetImageColorspace(image,sRGBColorspace,exception);
+  GetPixelInfo(image,&threshold);
   flags=ParseGeometry(thresholds,&geometry_info);
-  for (i=0; i < 5; i++)
-    threshold[i]=geometry_info.rho;
+  threshold.red=geometry_info.rho;
+  threshold.green=geometry_info.rho;
+  threshold.blue=geometry_info.rho;
+  threshold.black=geometry_info.rho;
+  threshold.alpha=100.0;
   if ((flags & SigmaValue) != 0)
-    threshold[1]=geometry_info.sigma;
+    threshold.green=geometry_info.sigma;
   if ((flags & XiValue) != 0)
-    threshold[2]=geometry_info.xi;
+    threshold.blue=geometry_info.xi;
   if ((flags & PsiValue) != 0)
-    threshold[3]=geometry_info.psi;
-  if ((flags & ChiValue) != 0)
-    threshold[4]=geometry_info.chi;
+    threshold.alpha=geometry_info.psi;
+  if (threshold.colorspace == CMYKColorspace)
+    {
+      if ((flags & PsiValue) != 0)
+        threshold.black=geometry_info.psi;
+      if ((flags & ChiValue) != 0)
+        threshold.alpha=geometry_info.chi;
+    }
   if ((flags & PercentValue) != 0)
-    for (i=0; i < 5; i++)
-      threshold[i]*=(QuantumRange/100.0);
+    {
+      threshold.red*=(MagickRealType) (QuantumRange/100.0);
+      threshold.green*=(MagickRealType) (QuantumRange/100.0);
+      threshold.blue*=(MagickRealType) (QuantumRange/100.0);
+      threshold.black*=(MagickRealType) (QuantumRange/100.0);
+      threshold.alpha*=(MagickRealType) (QuantumRange/100.0);
+    }
   /*
     White threshold image.
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
+  image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,8) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    magick_threads(image,image,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -550,27 +613,28 @@ MagickExport MagickBooleanType BlackThresholdImage(Image *image,
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        pixel;
+
       register ssize_t
         i;
 
-      ssize_t
-        n;
-
-      n=0;
+      if (GetPixelReadMask(image,q) == 0)
+        {
+          q+=GetPixelChannels(image);
+          continue;
+        }
+      pixel=GetPixelIntensity(image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          traits;
-
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
         if ((traits & UpdatePixelTrait) == 0)
           continue;
-        if ((MagickRealType) q[i] < threshold[n++ % 5])
-          q[i]=QuantumRange;
+        if (image->channel_mask != DefaultChannels)
+          pixel=(double) q[i];
+        if (pixel <= GetPixelInfoChannel(&threshold,channel))
+          q[i]=(Quantum) 0;
       }
       q+=GetPixelChannels(image);
     }
@@ -582,7 +646,7 @@ MagickExport MagickBooleanType BlackThresholdImage(Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_WhiteThresholdImage)
+        #pragma omp critical (MagickCore_BlackThresholdImage)
 #endif
         proceed=SetImageProgress(image,ThresholdImageTag,progress++,
           image->rows);
@@ -605,7 +669,9 @@ MagickExport MagickBooleanType BlackThresholdImage(Image *image,
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  ClampImage() restricts the color range from 0 to the quantum depth.
+%  ClampImage() set each pixel whose value is below zero to zero and any the
+%  pixel whose value is above the quantum range to the quantum range (e.g.
+%  65535) otherwise the pixel value remains unchanged.
 %
 %  The format of the ClampImage method is:
 %
@@ -619,16 +685,16 @@ MagickExport MagickBooleanType BlackThresholdImage(Image *image,
 %
 */
 
-static inline Quantum ClampToUnsignedQuantum(const Quantum quantum)
+static inline Quantum ClampPixel(const MagickRealType value)
 {
-#if defined(MAGICKCORE_HDRI_SUPPORT)
-  if (quantum <= 0)
-    return(0);
-  if (quantum >= QuantumRange)
-    return(QuantumRange);
-  return(quantum);
+#if !defined(MAGICKCORE_HDRI_SUPPORT)
+  return((Quantum) value);
 #else
-  return(quantum);
+  if (value < 0.0f)
+    return(0.0);
+  if (value >= (MagickRealType) QuantumRange)
+    return((Quantum) QuantumRange);
+  return(value);
 #endif
 }
 
@@ -663,10 +729,10 @@ MagickExport MagickBooleanType ClampImage(Image *image,ExceptionInfo *exception)
       q=image->colormap;
       for (i=0; i < (ssize_t) image->colors; i++)
       {
-        q->red=ClampToUnsignedQuantum(q->red);
-        q->green=ClampToUnsignedQuantum(q->green);
-        q->blue=ClampToUnsignedQuantum(q->blue);
-        q->alpha=ClampToUnsignedQuantum(q->alpha);
+        q->red=(double) ClampPixel(q->red);
+        q->green=(double) ClampPixel(q->green);
+        q->blue=(double) ClampPixel(q->blue);
+        q->alpha=(double) ClampPixel(q->alpha);
         q++;
       }
       return(SyncImage(image,exception));
@@ -676,9 +742,10 @@ MagickExport MagickBooleanType ClampImage(Image *image,ExceptionInfo *exception)
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
+  image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,8) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    magick_threads(image,image,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -701,19 +768,18 @@ MagickExport MagickBooleanType ClampImage(Image *image,ExceptionInfo *exception)
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,q) == 0)
+        {
+          q+=GetPixelChannels(image);
+          continue;
+        }
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          traits;
-
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
-        if (traits == UndefinedPixelTrait)
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        if ((traits & UpdatePixelTrait) == 0)
           continue;
-        q[i]=ClampToUnsignedQuantum(q[i]);
+        q[i]=ClampPixel(q[i]);
       }
       q+=GetPixelChannels(image);
     }
@@ -725,10 +791,9 @@ MagickExport MagickBooleanType ClampImage(Image *image,ExceptionInfo *exception)
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_ClampImage)
+        #pragma omp critical (MagickCore_ClampImage)
 #endif
-        proceed=SetImageProgress(image,ClampImageTag,progress++,
-          image->rows);
+        proceed=SetImageProgress(image,ClampImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
           status=MagickFalse;
       }
@@ -810,12 +875,19 @@ MagickExport ThresholdMap *GetThresholdMap(const char *map_id,
   ThresholdMap
     *map;
 
-  map=(ThresholdMap *)NULL;
+  map=GetThresholdMapFile(MinimalThresholdMap,"built-in",map_id,exception);
+  if (map != (ThresholdMap *) NULL)
+    return(map);
   options=GetConfigureOptions(ThresholdsFilename,exception);
-  while ((option=(const StringInfo *) GetNextValueInLinkedList(options)) !=
-         (const StringInfo *) NULL && (map == (ThresholdMap *) NULL))
+  option=(const StringInfo *) GetNextValueInLinkedList(options);
+  while (option != (const StringInfo *) NULL)
+  {
     map=GetThresholdMapFile((const char *) GetStringInfoDatum(option),
       GetStringInfoPath(option),map_id,exception);
+    if (map != (ThresholdMap *) NULL)
+      break;
+    option=(const StringInfo *) GetNextValueInLinkedList(options);
+  }
   options=DestroyConfigureOptions(options);
   return(map);
 }
@@ -850,8 +922,8 @@ MagickExport ThresholdMap *GetThresholdMap(const char *map_id,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-static ThresholdMap *GetThresholdMapFile(const char *xml,
-  const char *filename,const char *map_id,ExceptionInfo *exception)
+static ThresholdMap *GetThresholdMapFile(const char *xml,const char *filename,
+  const char *map_id,ExceptionInfo *exception)
 {
   char
     *p;
@@ -893,7 +965,10 @@ static ThresholdMap *GetThresholdMapFile(const char *xml,
       break;
   }
   if (threshold == (XMLTreeInfo *) NULL)
-    return(map);
+    {
+      thresholds=DestroyXMLTree(thresholds);
+      return(map);
+    }
   description=GetXMLTreeChild(threshold,"description");
   if (description == (XMLTreeInfo *) NULL)
     {
@@ -1157,12 +1232,13 @@ MagickExport MagickBooleanType ListThresholdMaps(FILE *file,
   options=GetConfigureOptions(ThresholdsFilename,exception);
   (void) FormatLocaleFile(file,
     "\n   Threshold Maps for Ordered Dither Operations\n");
-  while ((option=(const StringInfo *) GetNextValueInLinkedList(options)) !=
-         (const StringInfo *) NULL)
+  option=(const StringInfo *) GetNextValueInLinkedList(options);
+  while (option != (const StringInfo *) NULL)
   {
-    (void) FormatLocaleFile(file,"\nPATH: %s\n\n",GetStringInfoPath(option));
-    status|=ListThresholdMapFile(file,(const char *) GetStringInfoDatum(option),
+    (void) FormatLocaleFile(file,"\nPath: %s\n\n",GetStringInfoPath(option));
+    status&=ListThresholdMapFile(file,(const char *) GetStringInfoDatum(option),
       GetStringInfoPath(option),exception);
+    option=(const StringInfo *) GetNextValueInLinkedList(options);
   }
   options=DestroyConfigureOptions(options);
   return(status != 0 ? MagickTrue : MagickFalse);
@@ -1229,15 +1305,15 @@ MagickExport MagickBooleanType OrderedPosterizeImage(Image *image,
   const char
     *p;
 
+  double
+    levels[CompositePixelChannel];
+
   MagickBooleanType
     status;
 
   MagickOffsetType
     progress;
 
-  MagickRealType
-    levels[CompositePixelChannel];
-
   register ssize_t
     i;
 
@@ -1294,9 +1370,10 @@ MagickExport MagickBooleanType OrderedPosterizeImage(Image *image,
     return(MagickFalse);
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
+  image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,8) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    magick_threads(image,image,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1323,20 +1400,19 @@ MagickExport MagickBooleanType OrderedPosterizeImage(Image *image,
         n;
 
       n=0;
+      if (GetPixelReadMask(image,q) == 0)
+        {
+          q+=GetPixelChannels(image);
+          continue;
+        }
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          traits;
-
         ssize_t
           level,
           threshold;
 
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
         if ((traits & UpdatePixelTrait) == 0)
           continue;
         if (fabs(levels[n++]) < MagickEpsilon)
@@ -1344,8 +1420,9 @@ MagickExport MagickBooleanType OrderedPosterizeImage(Image *image,
         threshold=(ssize_t) (QuantumScale*q[i]*(levels[n]*(map->divisor-1)+1));
         level=threshold/(map->divisor-1);
         threshold-=level*(map->divisor-1);
-        q[i]=RoundToQuantum((level+(threshold >= map->levels[(x % map->width)+
-          map->width*(y % map->height)]))*QuantumRange/levels[n]);
+        q[i]=ClampToQuantum((double) (level+(threshold >=
+          map->levels[(x % map->width)+map->width*(y % map->height)]))*
+          QuantumRange/levels[n]);
         n++;
       }
       q+=GetPixelChannels(image);
@@ -1358,7 +1435,7 @@ MagickExport MagickBooleanType OrderedPosterizeImage(Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-#pragma omp critical (MagickCore_OrderedPosterizeImage)
+        #pragma omp critical (MagickCore_OrderedPosterizeImage)
 #endif
         proceed=SetImageProgress(image,DitherImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
@@ -1375,6 +1452,157 @@ MagickExport MagickBooleanType OrderedPosterizeImage(Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
+%     P e r c e p t i b l e I m a g e                                         %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  PerceptibleImage() set each pixel whose value is less than |epsilon| to
+%  epsilon or -epsilon (whichever is closer) otherwise the pixel value remains
+%  unchanged.
+%
+%  The format of the PerceptibleImage method is:
+%
+%      MagickBooleanType PerceptibleImage(Image *image,const double epsilon,
+%        ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o epsilon: the epsilon threshold (e.g. 1.0e-9).
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+
+static inline Quantum PerceptibleThreshold(const Quantum quantum,
+  const double epsilon)
+{
+  double
+    sign;
+
+  sign=(double) quantum < 0.0 ? -1.0 : 1.0;
+  if ((sign*quantum) >= epsilon)
+    return(quantum);
+  return((Quantum) (sign*epsilon));
+}
+
+MagickExport MagickBooleanType PerceptibleImage(Image *image,
+  const double epsilon,ExceptionInfo *exception)
+{
+#define PerceptibleImageTag  "Perceptible/Image"
+
+  CacheView
+    *image_view;
+
+  MagickBooleanType
+    status;
+
+  MagickOffsetType
+    progress;
+
+  ssize_t
+    y;
+
+  assert(image != (Image *) NULL);
+  assert(image->signature == MagickSignature);
+  if (image->debug != MagickFalse)
+    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  if (image->storage_class == PseudoClass)
+    {
+      register ssize_t
+        i;
+
+      register PixelInfo
+        *restrict q;
+
+      q=image->colormap;
+      for (i=0; i < (ssize_t) image->colors; i++)
+      {
+        q->red=(double) PerceptibleThreshold(ClampToQuantum(q->red),
+          epsilon);
+        q->green=(double) PerceptibleThreshold(ClampToQuantum(q->green),
+          epsilon);
+        q->blue=(double) PerceptibleThreshold(ClampToQuantum(q->blue),
+          epsilon);
+        q->alpha=(double) PerceptibleThreshold(ClampToQuantum(q->alpha),
+          epsilon);
+        q++;
+      }
+      return(SyncImage(image,exception));
+    }
+  /*
+    Perceptible image.
+  */
+  status=MagickTrue;
+  progress=0;
+  image_view=AcquireAuthenticCacheView(image,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    magick_threads(image,image,image->rows,1)
+#endif
+  for (y=0; y < (ssize_t) image->rows; y++)
+  {
+    register ssize_t
+      x;
+
+    register Quantum
+      *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++)
+    {
+      register ssize_t
+        i;
+
+      if (GetPixelReadMask(image,q) == 0)
+        {
+          q+=GetPixelChannels(image);
+          continue;
+        }
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
+        if (traits == UndefinedPixelTrait)
+          continue;
+        q[i]=PerceptibleThreshold(q[i],epsilon);
+      }
+      q+=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 critical (MagickCore_PerceptibleImage)
+#endif
+        proceed=SetImageProgress(image,PerceptibleImageTag,progress++,image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
+  }
+  image_view=DestroyCacheView(image_view);
+  return(status);
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
 %     R a n d o m T h r e s h o l d I m a g e                                 %
 %                                                                             %
 %                                                                             %
@@ -1409,6 +1637,10 @@ MagickExport MagickBooleanType RandomThresholdImage(Image *image,
   CacheView
     *image_view;
 
+  double
+    min_threshold,
+    max_threshold;
+
   GeometryInfo
     geometry_info;
 
@@ -1424,16 +1656,17 @@ MagickExport MagickBooleanType RandomThresholdImage(Image *image,
   PixelInfo
     threshold;
 
-  MagickRealType
-    min_threshold,
-    max_threshold;
-
   RandomInfo
     **restrict random_info;
 
   ssize_t
     y;
 
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  unsigned long
+    key;
+#endif
+
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
@@ -1442,9 +1675,11 @@ MagickExport MagickBooleanType RandomThresholdImage(Image *image,
   assert(exception->signature == MagickSignature);
   if (thresholds == (const char *) NULL)
     return(MagickTrue);
+  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
+    return(MagickFalse);
   GetPixelInfo(image,&threshold);
   min_threshold=0.0;
-  max_threshold=(MagickRealType) QuantumRange;
+  max_threshold=(double) QuantumRange;
   flags=ParseGeometry(thresholds,&geometry_info);
   min_threshold=geometry_info.rho;
   max_threshold=geometry_info.sigma;
@@ -1452,20 +1687,22 @@ MagickExport MagickBooleanType RandomThresholdImage(Image *image,
     max_threshold=min_threshold;
   if (strchr(thresholds,'%') != (char *) NULL)
     {
-      max_threshold*=(MagickRealType) (0.01*QuantumRange);
-      min_threshold*=(MagickRealType) (0.01*QuantumRange);
+      max_threshold*=(double) (0.01*QuantumRange);
+      min_threshold*=(double) (0.01*QuantumRange);
     }
   /*
     Random threshold image.
   */
   status=MagickTrue;
   progress=0;
-  if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
-    return(MagickFalse);
   random_info=AcquireRandomInfoThreadSet();
-  image_view=AcquireCacheView(image);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,8) shared(progress,status)
+  key=GetRandomSecretKey(random_info[0]);
+#endif
+  image_view=AcquireAuthenticCacheView(image,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    magick_threads(image,image,image->rows,key == ~0UL)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1491,31 +1728,29 @@ MagickExport MagickBooleanType RandomThresholdImage(Image *image,
       register ssize_t
         i;
 
+      if (GetPixelReadMask(image,q) == 0)
+        {
+          q+=GetPixelChannels(image);
+          continue;
+        }
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           threshold;
 
-        PixelChannel
-          channel;
-
-        PixelTrait
-          traits;
-
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
         if ((traits & UpdatePixelTrait) == 0)
           continue;
-        if ((MagickRealType) q[i] < min_threshold)
+        if ((double) q[i] < min_threshold)
           threshold=min_threshold;
         else
-          if ((MagickRealType) q[i] > max_threshold)
+          if ((double) q[i] > max_threshold)
             threshold=max_threshold;
           else
-            threshold=(MagickRealType) (QuantumRange*
+            threshold=(double) (QuantumRange*
               GetPseudoRandomValue(random_info[id]));
-          q[i]=(Quantum) ((MagickRealType) q[i] <= threshold ? 0 :
-            QuantumRange);
+        q[i]=(double) q[i] <= threshold ? 0 : QuantumRange;
       }
       q+=GetPixelChannels(image);
     }
@@ -1527,7 +1762,7 @@ MagickExport MagickBooleanType RandomThresholdImage(Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_RandomThresholdImage)
+        #pragma omp critical (MagickCore_RandomThresholdImage)
 #endif
         proceed=SetImageProgress(image,ThresholdImageTag,progress++,
           image->rows);
@@ -1586,15 +1821,12 @@ MagickExport MagickBooleanType WhiteThresholdImage(Image *image,
   MagickOffsetType
     progress;
 
-  MagickRealType
-    threshold[5];
+  PixelInfo
+    threshold;
 
   MagickStatusType
     flags;
 
-  register ssize_t
-    i;
-
   ssize_t
     y;
 
@@ -1606,28 +1838,45 @@ MagickExport MagickBooleanType WhiteThresholdImage(Image *image,
     return(MagickTrue);
   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
     return(MagickFalse);
+  if (IsGrayColorspace(image->colorspace) != MagickFalse)
+    (void) TransformImageColorspace(image,sRGBColorspace,exception);
+  GetPixelInfo(image,&threshold);
   flags=ParseGeometry(thresholds,&geometry_info);
-  for (i=0; i < 5; i++)
-    threshold[i]=geometry_info.rho;
+  threshold.red=geometry_info.rho;
+  threshold.green=geometry_info.rho;
+  threshold.blue=geometry_info.rho;
+  threshold.black=geometry_info.rho;
+  threshold.alpha=100.0;
   if ((flags & SigmaValue) != 0)
-    threshold[1]=geometry_info.sigma;
+    threshold.green=geometry_info.sigma;
   if ((flags & XiValue) != 0)
-    threshold[2]=geometry_info.xi;
+    threshold.blue=geometry_info.xi;
   if ((flags & PsiValue) != 0)
-    threshold[3]=geometry_info.psi;
-  if ((flags & ChiValue) != 0)
-    threshold[4]=geometry_info.chi;
+    threshold.alpha=geometry_info.psi;
+  if (threshold.colorspace == CMYKColorspace)
+    {
+      if ((flags & PsiValue) != 0)
+        threshold.black=geometry_info.psi;
+      if ((flags & ChiValue) != 0)
+        threshold.alpha=geometry_info.chi;
+    }
   if ((flags & PercentValue) != 0)
-    for (i=0; i < 5; i++)
-      threshold[i]*=(QuantumRange/100.0);
+    {
+      threshold.red*=(MagickRealType) (QuantumRange/100.0);
+      threshold.green*=(MagickRealType) (QuantumRange/100.0);
+      threshold.blue*=(MagickRealType) (QuantumRange/100.0);
+      threshold.black*=(MagickRealType) (QuantumRange/100.0);
+      threshold.alpha*=(MagickRealType) (QuantumRange/100.0);
+    }
   /*
     White threshold image.
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
+  image_view=AcquireAuthenticCacheView(image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,8) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    magick_threads(image,image,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -1647,26 +1896,27 @@ MagickExport MagickBooleanType WhiteThresholdImage(Image *image,
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      double
+        pixel;
+
       register ssize_t
         i;
 
-      ssize_t
-        n;
-
-      n=0;
+      if (GetPixelReadMask(image,q) == 0)
+        {
+          q+=GetPixelChannels(image);
+          continue;
+        }
+      pixel=GetPixelIntensity(image,q);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        PixelChannel
-          channel;
-
-        PixelTrait
-          traits;
-
-        channel=GetPixelChannelMapChannel(image,i);
-        traits=GetPixelChannelMapTraits(image,channel);
+        PixelChannel channel=GetPixelChannelChannel(image,i);
+        PixelTrait traits=GetPixelChannelTraits(image,channel);
         if ((traits & UpdatePixelTrait) == 0)
           continue;
-        if ((MagickRealType) q[i] > threshold[n++ % 5])
+        if (image->channel_mask != DefaultChannels)
+          pixel=(double) q[i];
+        if (pixel > GetPixelInfoChannel(&threshold,channel))
           q[i]=QuantumRange;
       }
       q+=GetPixelChannels(image);
@@ -1679,7 +1929,7 @@ MagickExport MagickBooleanType WhiteThresholdImage(Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_WhiteThresholdImage)
+        #pragma omp critical (MagickCore_WhiteThresholdImage)
 #endif
         proceed=SetImageProgress(image,ThresholdImageTag,progress++,
           image->rows);