granicus.if.org Git - imagemagick/blob - MagickCore/threshold.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %       TTTTT  H   H  RRRR   EEEEE  SSSSS  H   H   OOO   L      DDDD          %
   7 %         T    H   H  R   R  E      SS     H   H  O   O  L      D   D         %
   8 %         T    HHHHH  RRRR   EEE     SSS   HHHHH  O   O  L      D   D         %
   9 %         T    H   H  R R    E         SS  H   H  O   O  L      D   D         %
  10 %         T    H   H  R  R   EEEEE  SSSSS  H   H   OOO   LLLLL  DDDD          %
  11 %                                                                             %
  12 %                                                                             %
  13 %                      MagickCore Image Threshold Methods                     %
  14 %                                                                             %
  15 %                               Software Design                               %
  16 %                                    Cristy                                   %
  17 %                                 October 1996                                %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2017 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    https://www.imagemagick.org/script/license.php                           %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "MagickCore/studio.h"
  44 #include "MagickCore/property.h"
  45 #include "MagickCore/blob.h"
  46 #include "MagickCore/cache-view.h"
  47 #include "MagickCore/color.h"
  48 #include "MagickCore/color-private.h"
  49 #include "MagickCore/colormap.h"
  50 #include "MagickCore/colorspace.h"
  51 #include "MagickCore/colorspace-private.h"
  52 #include "MagickCore/configure.h"
  53 #include "MagickCore/constitute.h"
  54 #include "MagickCore/decorate.h"
  55 #include "MagickCore/draw.h"
  56 #include "MagickCore/enhance.h"
  57 #include "MagickCore/exception.h"
  58 #include "MagickCore/exception-private.h"
  59 #include "MagickCore/effect.h"
  60 #include "MagickCore/fx.h"
  61 #include "MagickCore/gem.h"
  62 #include "MagickCore/geometry.h"
  63 #include "MagickCore/image-private.h"
  64 #include "MagickCore/list.h"
  65 #include "MagickCore/log.h"
  66 #include "MagickCore/memory_.h"
  67 #include "MagickCore/monitor.h"
  68 #include "MagickCore/monitor-private.h"
  69 #include "MagickCore/montage.h"
  70 #include "MagickCore/option.h"
  71 #include "MagickCore/pixel-accessor.h"
  72 #include "MagickCore/pixel-private.h"
  73 #include "MagickCore/quantize.h"
  74 #include "MagickCore/quantum.h"
  75 #include "MagickCore/quantum-private.h"
  76 #include "MagickCore/random_.h"
  77 #include "MagickCore/random-private.h"
  78 #include "MagickCore/resize.h"
  79 #include "MagickCore/resource_.h"
  80 #include "MagickCore/segment.h"
  81 #include "MagickCore/shear.h"
  82 #include "MagickCore/signature-private.h"
  83 #include "MagickCore/string_.h"
  84 #include "MagickCore/string-private.h"
  85 #include "MagickCore/thread-private.h"
  86 #include "MagickCore/threshold.h"
  87 #include "MagickCore/token.h"
  88 #include "MagickCore/transform.h"
  89 #include "MagickCore/xml-tree.h"
  90 #include "MagickCore/xml-tree-private.h"
  91 \f
  92 /*
  93   Define declarations.
  94 */
  95 #define ThresholdsFilename  "thresholds.xml"
  96 \f
  97 /*
  98   Typedef declarations.
  99 */
 100 struct _ThresholdMap
 101 {
 102   char
 103     *map_id,
 104     *description;
 105
 106   size_t
 107     width,
 108     height;
 109
 110   ssize_t
 111     divisor,
 112     *levels;
 113 };
 114 \f
 115 /*
 116   Static declarations.
 117 */
 118 static const char
 119   *MinimalThresholdMap =
 120     "<?xml version=\"1.0\"?>"
 121     "<thresholds>"
 122     "  <threshold map=\"threshold\" alias=\"1x1\">"
 123     "    <description>Threshold 1x1 (non-dither)</description>"
 124     "    <levels width=\"1\" height=\"1\" divisor=\"2\">"
 125     "        1"
 126     "    </levels>"
 127     "  </threshold>"
 128     "  <threshold map=\"checks\" alias=\"2x1\">"
 129     "    <description>Checkerboard 2x1 (dither)</description>"
 130     "    <levels width=\"2\" height=\"2\" divisor=\"3\">"
 131     "       1 2"
 132     "       2 1"
 133     "    </levels>"
 134     "  </threshold>"
 135     "</thresholds>";
 136 \f
 137 /*
 138   Forward declarations.
 139 */
 140 static ThresholdMap
 141   *GetThresholdMapFile(const char *,const char *,const char *,ExceptionInfo *);
 142 \f
 143 /*
 144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 145 %                                                                             %
 146 %                                                                             %
 147 %                                                                             %
 148 %     A d a p t i v e T h r e s h o l d I m a g e                             %
 149 %                                                                             %
 150 %                                                                             %
 151 %                                                                             %
 152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 153 %
 154 %  AdaptiveThresholdImage() selects an individual threshold for each pixel
 155 %  based on the range of intensity values in its local neighborhood.  This
 156 %  allows for thresholding of an image whose global intensity histogram
 157 %  doesn't contain distinctive peaks.
 158 %
 159 %  The format of the AdaptiveThresholdImage method is:
 160 %
 161 %      Image *AdaptiveThresholdImage(const Image *image,const size_t width,
 162 %        const size_t height,const double bias,ExceptionInfo *exception)
 163 %
 164 %  A description of each parameter follows:
 165 %
 166 %    o image: the image.
 167 %
 168 %    o width: the width of the local neighborhood.
 169 %
 170 %    o height: the height of the local neighborhood.
 171 %
 172 %    o bias: the mean bias.
 173 %
 174 %    o exception: return any errors or warnings in this structure.
 175 %
 176 */
 177 MagickExport Image *AdaptiveThresholdImage(const Image *image,
 178   const size_t width,const size_t height,const double bias,
 179   ExceptionInfo *exception)
 180 {
 181 #define AdaptiveThresholdImageTag  "AdaptiveThreshold/Image"
 182
 183   CacheView
 184     *image_view,
 185     *threshold_view;
 186
 187   Image
 188     *threshold_image;
 189
 190   MagickBooleanType
 191     status;
 192
 193   MagickOffsetType
 194     progress;
 195
 196   MagickSizeType
 197     number_pixels;
 198
 199   ssize_t
 200     y;
 201
 202   /*
 203     Initialize threshold image attributes.
 204   */
 205   assert(image != (Image *) NULL);
 206   assert(image->signature == MagickCoreSignature);
 207   if (image->debug != MagickFalse)
 208     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 209   assert(exception != (ExceptionInfo *) NULL);
 210   assert(exception->signature == MagickCoreSignature);
 211   threshold_image=CloneImage(image,image->columns,image->rows,MagickTrue,
 212     exception);
 213   if (threshold_image == (Image *) NULL)
 214     return((Image *) NULL);
 215   status=SetImageStorageClass(threshold_image,DirectClass,exception);
 216   if (status == MagickFalse)
 217     {
 218       threshold_image=DestroyImage(threshold_image);
 219       return((Image *) NULL);
 220     }
 221   /*
 222     Threshold image.
 223   */
 224   status=MagickTrue;
 225   progress=0;
 226   number_pixels=(MagickSizeType) width*height;
 227   image_view=AcquireVirtualCacheView(image,exception);
 228   threshold_view=AcquireAuthenticCacheView(threshold_image,exception);
 229 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 230   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 231     magick_threads(image,threshold_image,image->rows,1)
 232 #endif
 233   for (y=0; y < (ssize_t) image->rows; y++)
 234   {
 235     double
 236       channel_bias[MaxPixelChannels],
 237       channel_sum[MaxPixelChannels];
 238
 239     register const Quantum
 240       *magick_restrict p,
 241       *magick_restrict pixels;
 242
 243     register Quantum
 244       *magick_restrict q;
 245
 246     register ssize_t
 247       i,
 248       x;
 249
 250     ssize_t
 251       center,
 252       u,
 253       v;
 254
 255     if (status == MagickFalse)
 256       continue;
 257     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
 258       (height/2L),image->columns+width,height,exception);
 259     q=QueueCacheViewAuthenticPixels(threshold_view,0,y,threshold_image->columns,
 260       1,exception);
 261     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
 262       {
 263         status=MagickFalse;
 264         continue;
 265       }
 266     center=(ssize_t) GetPixelChannels(image)*(image->columns+width)*(height/2L)+
 267       GetPixelChannels(image)*(width/2);
 268     for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 269     {
 270       PixelChannel channel=GetPixelChannelChannel(image,i);
 271       PixelTrait traits=GetPixelChannelTraits(image,channel);
 272       PixelTrait threshold_traits=GetPixelChannelTraits(threshold_image,
 273         channel);
 274       if ((traits == UndefinedPixelTrait) ||
 275           (threshold_traits == UndefinedPixelTrait))
 276         continue;
 277       if (((threshold_traits & CopyPixelTrait) != 0) ||
 278           (GetPixelWriteMask(image,p) == 0))
 279         {
 280           SetPixelChannel(threshold_image,channel,p[center+i],q);
 281           continue;
 282         }
 283       pixels=p;
 284       channel_bias[channel]=0.0;
 285       channel_sum[channel]=0.0;
 286       for (v=0; v < (ssize_t) height; v++)
 287       {
 288         for (u=0; u < (ssize_t) width; u++)
 289         {
 290           if (u == (ssize_t) (width-1))
 291             channel_bias[channel]+=pixels[i];
 292           channel_sum[channel]+=pixels[i];
 293           pixels+=GetPixelChannels(image);
 294         }
 295         pixels+=GetPixelChannels(image)*image->columns;
 296       }
 297     }
 298     for (x=0; x < (ssize_t) image->columns; x++)
 299     {
 300       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 301       {
 302         double
 303           mean;
 304
 305         PixelChannel channel=GetPixelChannelChannel(image,i);
 306         PixelTrait traits=GetPixelChannelTraits(image,channel);
 307         PixelTrait threshold_traits=GetPixelChannelTraits(threshold_image,
 308           channel);
 309         if ((traits == UndefinedPixelTrait) ||
 310             (threshold_traits == UndefinedPixelTrait))
 311           continue;
 312         if (((threshold_traits & CopyPixelTrait) != 0) ||
 313             (GetPixelWriteMask(image,p) == 0))
 314           {
 315             SetPixelChannel(threshold_image,channel,p[center+i],q);
 316             continue;
 317           }
 318         channel_sum[channel]-=channel_bias[channel];
 319         channel_bias[channel]=0.0;
 320         pixels=p;
 321         for (v=0; v < (ssize_t) height; v++)
 322         {
 323           channel_bias[channel]+=pixels[i];
 324           pixels+=(width-1)*GetPixelChannels(image);
 325           channel_sum[channel]+=pixels[i];
 326           pixels+=GetPixelChannels(image)*(image->columns+1);
 327         }
 328         mean=(double) (channel_sum[channel]/number_pixels+bias);
 329         SetPixelChannel(threshold_image,channel,(Quantum) ((double)
 330           p[center+i] <= mean ? 0 : QuantumRange),q);
 331       }
 332       p+=GetPixelChannels(image);
 333       q+=GetPixelChannels(threshold_image);
 334     }
 335     if (SyncCacheViewAuthenticPixels(threshold_view,exception) == MagickFalse)
 336       status=MagickFalse;
 337     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 338       {
 339         MagickBooleanType
 340           proceed;
 341
 342 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 343         #pragma omp critical (MagickCore_AdaptiveThresholdImage)
 344 #endif
 345         proceed=SetImageProgress(image,AdaptiveThresholdImageTag,progress++,
 346           image->rows);
 347         if (proceed == MagickFalse)
 348           status=MagickFalse;
 349       }
 350   }
 351   threshold_image->type=image->type;
 352   threshold_view=DestroyCacheView(threshold_view);
 353   image_view=DestroyCacheView(image_view);
 354   if (status == MagickFalse)
 355     threshold_image=DestroyImage(threshold_image);
 356   return(threshold_image);
 357 }
 358 \f
 359 /*
 360 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 361 %                                                                             %
 362 %                                                                             %
 363 %                                                                             %
 364 %     A u t o T h r e s o l d I m a g e                                       %
 365 %                                                                             %
 366 %                                                                             %
 367 %                                                                             %
 368 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 369 %
 370 %  AutoThresholdImage() automatically selects a threshold and replaces each
 371 %  pixel in the image with a black pixel if the image intentsity is less than
 372 %  the selected threshold otherwise white.
 373 %
 374 %  The format of the AutoThresholdImage method is:
 375 %
 376 %      MagickBooleanType AutoThresholdImage(Image *image,
 377 %        const AutoThresholdMethod method,ExceptionInfo *exception)
 378 %
 379 %  A description of each parameter follows:
 380 %
 381 %    o image: The image to auto-threshold.
 382 %
 383 %    o method: choose from Kapur, OTSU, or Triangle.
 384 %
 385 %    o exception: return any errors or warnings in this structure.
 386 %
 387 */
 388
 389 static double KapurThreshold(const Image *image,const double *histogram,
 390   ExceptionInfo *exception)
 391 {
 392 #define MaxIntensity  255
 393
 394   double
 395     *black_entropy,
 396     *cumulative_histogram,
 397     entropy,
 398     epsilon,
 399     maximum_entropy,
 400     *white_entropy;
 401
 402   register ssize_t
 403     i,
 404     j;
 405
 406   size_t
 407     threshold;
 408
 409   /*
 410     Compute optimal threshold from the entopy of the histogram.
 411   */
 412   cumulative_histogram=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
 413     sizeof(*cumulative_histogram));
 414   black_entropy=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
 415     sizeof(*black_entropy));
 416   white_entropy=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
 417     sizeof(*white_entropy));
 418   if ((cumulative_histogram == (double *) NULL) ||
 419       (black_entropy == (double *) NULL) || (white_entropy == (double *) NULL))
 420     {
 421       if (white_entropy != (double *) NULL)
 422         white_entropy=(double *) RelinquishMagickMemory(white_entropy);
 423       if (black_entropy != (double *) NULL)
 424         black_entropy=(double *) RelinquishMagickMemory(black_entropy);
 425       if (cumulative_histogram != (double *) NULL)
 426         cumulative_histogram=(double *)
 427           RelinquishMagickMemory(cumulative_histogram);
 428       (void) ThrowMagickException(exception,GetMagickModule(),
 429         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 430       return(-1.0);
 431     }
 432    /*
 433      Entropy for black and white parts of the histogram.
 434    */
 435    cumulative_histogram[0]=histogram[0];
 436    for (i=1; i <= MaxIntensity; i++)
 437      cumulative_histogram[i]=cumulative_histogram[i-1]+histogram[i];
 438    epsilon=MagickMinimumValue;
 439    for (j=0; j <= MaxIntensity; j++)
 440    {
 441      /*
 442        Black entropy.
 443      */
 444      black_entropy[j]=0.0;
 445      if (cumulative_histogram[j] > epsilon)
 446        {
 447          entropy=0.0;
 448          for (i=0; i <= j; i++)
 449            if (histogram[i] > epsilon)
 450              entropy-=histogram[i]/cumulative_histogram[j]*
 451                log(histogram[i]/cumulative_histogram[j]);
 452          black_entropy[j]=entropy;
 453        }
 454      /*
 455        White entropy.
 456      */
 457      white_entropy[j]=0.0;
 458      if ((1.0-cumulative_histogram[j]) > epsilon)
 459        {
 460          entropy=0.0;
 461          for (i=j+1; i <= MaxIntensity; i++)
 462            if (histogram[i] > epsilon)
 463              entropy-=histogram[i]/(1.0-cumulative_histogram[j])*
 464                log(histogram[i]/(1.0-cumulative_histogram[j]));
 465          white_entropy[j]=entropy;
 466        }
 467    }
 468   /*
 469     Find histogram bin with maximum entropy.
 470   */
 471   maximum_entropy=black_entropy[0]+white_entropy[0];
 472   threshold=0;
 473   for (j=1; j <= MaxIntensity; j++)
 474     if ((black_entropy[j]+white_entropy[j]) > maximum_entropy)
 475       {
 476         maximum_entropy=black_entropy[j]+white_entropy[j];
 477         threshold=(ssize_t) j;
 478       }
 479   /*
 480     Free resources.
 481   */
 482   white_entropy=(double *) RelinquishMagickMemory(white_entropy);
 483   black_entropy=(double *) RelinquishMagickMemory(black_entropy);
 484   cumulative_histogram=(double *) RelinquishMagickMemory(cumulative_histogram);
 485   return(100.0*threshold/MaxIntensity);
 486 }
 487
 488 static double OTSUThreshold(const Image *image,const double *histogram,
 489   ExceptionInfo *exception)
 490 {
 491   double
 492     max_sigma,
 493     *myu,
 494     *omega,
 495     *probability,
 496     *sigma,
 497     threshold;
 498
 499   register ssize_t
 500     i;
 501
 502   /*
 503     Compute optimal threshold from maximization of inter-class variance.
 504   */
 505   myu=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*myu));
 506   omega=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*omega));
 507   probability=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
 508     sizeof(*probability));
 509   sigma=(double *) AcquireQuantumMemory(MaxIntensity+1UL,sizeof(*sigma));
 510   if ((myu == (double *) NULL) || (omega == (double *) NULL) ||
 511       (probability == (double *) NULL) || (sigma == (double *) NULL))
 512     {
 513       if (sigma != (double *) NULL)
 514         sigma=(double *) RelinquishMagickMemory(sigma);
 515       if (probability != (double *) NULL)
 516         probability=(double *) RelinquishMagickMemory(probability);
 517       if (omega != (double *) NULL)
 518         omega=(double *) RelinquishMagickMemory(omega);
 519       if (myu != (double *) NULL)
 520         myu=(double *) RelinquishMagickMemory(myu);
 521       (void) ThrowMagickException(exception,GetMagickModule(),
 522         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 523       return(-1.0);
 524     }
 525   /*
 526     Calculate probability density.
 527   */
 528   for (i=0; i <= (ssize_t) MaxIntensity; i++)
 529     probability[i]=histogram[i];
 530   /*
 531     Generate probability of graylevels and mean value for separation.
 532   */
 533   omega[0]=probability[0];
 534   myu[0]=0.0;
 535   for (i=0; i <= (ssize_t) MaxIntensity; i++)
 536   {
 537     omega[i]=omega[i-1]+probability[i];
 538     myu[i]=myu[i-1]+i*probability[i];
 539   }
 540   /*
 541     Sigma maximization: inter-class variance and compute optimal threshold.
 542   */
 543   threshold=0;
 544   max_sigma=0.0;
 545   for (i=0; i < (ssize_t) MaxIntensity; i++)
 546   {
 547     sigma[i]=0.0;
 548     if ((omega[i] != 0.0) && (omega[i] != 1.0))
 549       sigma[i]=pow(myu[MaxIntensity]*omega[i]-myu[i],2.0)/(omega[i]*(1.0-
 550         omega[i]));
 551     if (sigma[i] > max_sigma)
 552       {
 553         max_sigma=sigma[i];
 554         threshold=(double) i;
 555       }
 556   }
 557   /*
 558     Free resources.
 559   */
 560   myu=(double *) RelinquishMagickMemory(myu);
 561   omega=(double *) RelinquishMagickMemory(omega);
 562   probability=(double *) RelinquishMagickMemory(probability);
 563   sigma=(double *) RelinquishMagickMemory(sigma);
 564   return(100.0*threshold/MaxIntensity);
 565 }
 566
 567 static double TriangleThreshold(const Image *image,const double *histogram,
 568   ExceptionInfo *exception)
 569 {
 570   double
 571     a,
 572     b,
 573     c,
 574     count,
 575     distance,
 576     inverse_ratio,
 577     max_distance,
 578     segment,
 579     x1,
 580     x2,
 581     y1,
 582     y2;
 583
 584   register ssize_t
 585     i;
 586
 587   ssize_t
 588     end,
 589     max,
 590     start,
 591     threshold;
 592
 593   /*
 594     Compute optimal threshold with triangle algorithm.
 595   */
 596   start=0;  /* find start bin, first bin not zero count */
 597   for (i=0; i <= (ssize_t) MaxIntensity; i++)
 598     if (histogram[i] > 0.0)
 599       {
 600         start=i;
 601         break;
 602       }
 603   end=0;  /* find end bin, last bin not zero count */
 604   for (i=(ssize_t) MaxIntensity; i >= 0; i--)
 605     if (histogram[i] > 0.0)
 606       {
 607         end=i;
 608         break;
 609       }
 610   max=0;  /* find max bin, bin with largest count */
 611   count=0.0;
 612   for (i=0; i <= (ssize_t) MaxIntensity; i++)
 613     if (histogram[i] > count)
 614       {
 615         max=i;
 616         count=histogram[i];
 617       }
 618   /*
 619     Compute threshold at split point.
 620   */
 621   x1=(double) max;
 622   y1=histogram[max];
 623   x2=(double) end;
 624   if ((max-start) >= (end-max))
 625     x2=(double) start;
 626   y2=0.0;
 627   a=y1-y2;
 628   b=x2-x1;
 629   c=(-1.0)*(a*x1+b*y1);
 630   inverse_ratio=1.0/sqrt(a*a+b*b+c*c);
 631   threshold=0;
 632   max_distance=0.0;
 633   if (x2 == (double) start)
 634     for (i=start; i < max; i++)
 635     {
 636       segment=inverse_ratio*(a*i+b*histogram[i]+c);
 637       distance=sqrt(segment*segment);
 638       if ((distance > max_distance) && (segment > 0.0))
 639         {
 640           threshold=i;
 641           max_distance=distance;
 642         }
 643     }
 644   else
 645     for (i=end; i > max; i--)
 646     {
 647       segment=inverse_ratio*(a*i+b*histogram[i]+c);
 648       distance=sqrt(segment*segment);
 649       if ((distance > max_distance) && (segment < 0.0))
 650         {
 651           threshold=i;
 652           max_distance=distance;
 653         }
 654     }
 655   return(100.0*threshold/MaxIntensity);
 656 }
 657
 658 MagickExport MagickBooleanType AutoThresholdImage(Image *image,
 659   const AutoThresholdMethod method,ExceptionInfo *exception)
 660 {
 661   CacheView
 662     *image_view;
 663
 664   char
 665     property[MagickPathExtent];
 666
 667   double
 668     gamma,
 669     *histogram,
 670     sum,
 671     threshold;
 672
 673   MagickBooleanType
 674     status;
 675
 676   register ssize_t
 677     i;
 678
 679   ssize_t
 680     y;
 681
 682   /*
 683     Form histogram.
 684   */
 685   assert(image != (Image *) NULL);
 686   assert(image->signature == MagickCoreSignature);
 687   if (image->debug != MagickFalse)
 688     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 689   histogram=(double *) AcquireQuantumMemory(MaxIntensity+1UL,
 690     sizeof(*histogram));
 691   if (histogram == (double *) NULL)
 692     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
 693       image->filename);
 694   status=MagickTrue;
 695   (void) ResetMagickMemory(histogram,0,(MaxIntensity+1UL)*sizeof(*histogram));
 696   image_view=AcquireVirtualCacheView(image,exception);
 697   for (y=0; y < (ssize_t) image->rows; y++)
 698   {
 699     register const Quantum
 700       *magick_restrict p;
 701
 702     register ssize_t
 703       x;
 704
 705     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 706     if (p == (const Quantum *) NULL)
 707       break;
 708     for (x=0; x < (ssize_t) image->columns; x++)
 709     {
 710       double intensity = GetPixelIntensity(image,p);
 711       histogram[ScaleQuantumToChar(ClampToQuantum(intensity))]++;
 712       p+=GetPixelChannels(image);
 713     }
 714   }
 715   image_view=DestroyCacheView(image_view);
 716   /*
 717     Normalize histogram.
 718   */
 719   sum=0.0;
 720   for (i=0; i <= (ssize_t) MaxIntensity; i++)
 721     sum+=histogram[i];
 722   gamma=PerceptibleReciprocal(sum);
 723   for (i=0; i <= (ssize_t) MaxIntensity; i++)
 724     histogram[i]=gamma*histogram[i];
 725   /*
 726     Discover threshold from histogram.
 727   */
 728   switch (method)
 729   {
 730     case KapurThresholdMethod:
 731     {
 732       threshold=KapurThreshold(image,histogram,exception);
 733       break;
 734     }
 735     case OTSUThresholdMethod:
 736     default:
 737     {
 738       threshold=OTSUThreshold(image,histogram,exception);
 739       break;
 740     }
 741     case TriangleThresholdMethod:
 742     {
 743       threshold=TriangleThreshold(image,histogram,exception);
 744       break;
 745     }
 746   }
 747   histogram=(double *) RelinquishMagickMemory(histogram);
 748   if (threshold < 0.0)
 749     status=MagickFalse;
 750   if (status == MagickFalse)
 751     return(MagickFalse);
 752   /*
 753     Threshold image.
 754   */
 755   (void) FormatLocaleString(property,MagickPathExtent,"%g%%",threshold);
 756   (void) SetImageProperty(image,"auto-threshold:threshold",property,exception);
 757   return(BilevelImage(image,QuantumRange*threshold/100.0,exception));
 758 }
 759 \f
 760 /*
 761 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 762 %                                                                             %
 763 %                                                                             %
 764 %                                                                             %
 765 %     B i l e v e l I m a g e                                                 %
 766 %                                                                             %
 767 %                                                                             %
 768 %                                                                             %
 769 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 770 %
 771 %  BilevelImage() changes the value of individual pixels based on the
 772 %  intensity of each pixel channel.  The result is a high-contrast image.
 773 %
 774 %  More precisely each channel value of the image is 'thresholded' so that if
 775 %  it is equal to or less than the given value it is set to zero, while any
 776 %  value greater than that give is set to it maximum or QuantumRange.
 777 %
 778 %  This function is what is used to implement the "-threshold" operator for
 779 %  the command line API.
 780 %
 781 %  If the default channel setting is given the image is thresholded using just
 782 %  the gray 'intensity' of the image, rather than the individual channels.
 783 %
 784 %  The format of the BilevelImage method is:
 785 %
 786 %      MagickBooleanType BilevelImage(Image *image,const double threshold,
 787 %        ExceptionInfo *exception)
 788 %
 789 %  A description of each parameter follows:
 790 %
 791 %    o image: the image.
 792 %
 793 %    o threshold: define the threshold values.
 794 %
 795 %    o exception: return any errors or warnings in this structure.
 796 %
 797 %  Aside: You can get the same results as operator using LevelImages()
 798 %  with the 'threshold' value for both the black_point and the white_point.
 799 %
 800 */
 801 MagickExport MagickBooleanType BilevelImage(Image *image,const double threshold,
 802   ExceptionInfo *exception)
 803 {
 804 #define ThresholdImageTag  "Threshold/Image"
 805
 806   CacheView
 807     *image_view;
 808
 809   MagickBooleanType
 810     status;
 811
 812   MagickOffsetType
 813     progress;
 814
 815   ssize_t
 816     y;
 817
 818   assert(image != (Image *) NULL);
 819   assert(image->signature == MagickCoreSignature);
 820   if (image->debug != MagickFalse)
 821     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 822   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 823     return(MagickFalse);
 824   if (IsGrayColorspace(image->colorspace) != MagickFalse)
 825     (void) SetImageColorspace(image,sRGBColorspace,exception);
 826   /*
 827     Bilevel threshold image.
 828   */
 829   status=MagickTrue;
 830   progress=0;
 831   image_view=AcquireAuthenticCacheView(image,exception);
 832 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 833   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 834     magick_threads(image,image,image->rows,1)
 835 #endif
 836   for (y=0; y < (ssize_t) image->rows; y++)
 837   {
 838     register ssize_t
 839       x;
 840
 841     register Quantum
 842       *magick_restrict q;
 843
 844     if (status == MagickFalse)
 845       continue;
 846     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 847     if (q == (Quantum *) NULL)
 848       {
 849         status=MagickFalse;
 850         continue;
 851       }
 852     for (x=0; x < (ssize_t) image->columns; x++)
 853     {
 854       double
 855         pixel;
 856
 857       register ssize_t
 858         i;
 859
 860       if (GetPixelWriteMask(image,q) == 0)
 861         {
 862           q+=GetPixelChannels(image);
 863           continue;
 864         }
 865       pixel=GetPixelIntensity(image,q);
 866       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 867       {
 868         PixelChannel channel=GetPixelChannelChannel(image,i);
 869         PixelTrait traits=GetPixelChannelTraits(image,channel);
 870         if ((traits & UpdatePixelTrait) == 0)
 871           continue;
 872         if (image->channel_mask != DefaultChannels)
 873           pixel=(double) q[i];
 874         q[i]=(Quantum) (pixel <= threshold ? 0 : QuantumRange);
 875       }
 876       q+=GetPixelChannels(image);
 877     }
 878     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 879       status=MagickFalse;
 880     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 881       {
 882         MagickBooleanType
 883           proceed;
 884
 885 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 886         #pragma omp critical (MagickCore_BilevelImage)
 887 #endif
 888         proceed=SetImageProgress(image,ThresholdImageTag,progress++,
 889           image->rows);
 890         if (proceed == MagickFalse)
 891           status=MagickFalse;
 892       }
 893   }
 894   image_view=DestroyCacheView(image_view);
 895   return(status);
 896 }
 897 \f
 898 /*
 899 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 900 %                                                                             %
 901 %                                                                             %
 902 %                                                                             %
 903 %     B l a c k T h r e s h o l d I m a g e                                   %
 904 %                                                                             %
 905 %                                                                             %
 906 %                                                                             %
 907 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 908 %
 909 %  BlackThresholdImage() is like ThresholdImage() but forces all pixels below
 910 %  the threshold into black while leaving all pixels at or above the threshold
 911 %  unchanged.
 912 %
 913 %  The format of the BlackThresholdImage method is:
 914 %
 915 %      MagickBooleanType BlackThresholdImage(Image *image,
 916 %        const char *threshold,ExceptionInfo *exception)
 917 %
 918 %  A description of each parameter follows:
 919 %
 920 %    o image: the image.
 921 %
 922 %    o threshold: define the threshold value.
 923 %
 924 %    o exception: return any errors or warnings in this structure.
 925 %
 926 */
 927 MagickExport MagickBooleanType BlackThresholdImage(Image *image,
 928   const char *thresholds,ExceptionInfo *exception)
 929 {
 930 #define ThresholdImageTag  "Threshold/Image"
 931
 932   CacheView
 933     *image_view;
 934
 935   GeometryInfo
 936     geometry_info;
 937
 938   MagickBooleanType
 939     status;
 940
 941   MagickOffsetType
 942     progress;
 943
 944   PixelInfo
 945     threshold;
 946
 947   MagickStatusType
 948     flags;
 949
 950   ssize_t
 951     y;
 952
 953   assert(image != (Image *) NULL);
 954   assert(image->signature == MagickCoreSignature);
 955   if (image->debug != MagickFalse)
 956     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 957   if (thresholds == (const char *) NULL)
 958     return(MagickTrue);
 959   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 960     return(MagickFalse);
 961   if (IsGrayColorspace(image->colorspace) != MagickFalse)
 962     (void) SetImageColorspace(image,sRGBColorspace,exception);
 963   GetPixelInfo(image,&threshold);
 964   flags=ParseGeometry(thresholds,&geometry_info);
 965   threshold.red=geometry_info.rho;
 966   threshold.green=geometry_info.rho;
 967   threshold.blue=geometry_info.rho;
 968   threshold.black=geometry_info.rho;
 969   threshold.alpha=100.0;
 970   if ((flags & SigmaValue) != 0)
 971     threshold.green=geometry_info.sigma;
 972   if ((flags & XiValue) != 0)
 973     threshold.blue=geometry_info.xi;
 974   if ((flags & PsiValue) != 0)
 975     threshold.alpha=geometry_info.psi;
 976   if (threshold.colorspace == CMYKColorspace)
 977     {
 978       if ((flags & PsiValue) != 0)
 979         threshold.black=geometry_info.psi;
 980       if ((flags & ChiValue) != 0)
 981         threshold.alpha=geometry_info.chi;
 982     }
 983   if ((flags & PercentValue) != 0)
 984     {
 985       threshold.red*=(MagickRealType) (QuantumRange/100.0);
 986       threshold.green*=(MagickRealType) (QuantumRange/100.0);
 987       threshold.blue*=(MagickRealType) (QuantumRange/100.0);
 988       threshold.black*=(MagickRealType) (QuantumRange/100.0);
 989       threshold.alpha*=(MagickRealType) (QuantumRange/100.0);
 990     }
 991   /*
 992     White threshold image.
 993   */
 994   status=MagickTrue;
 995   progress=0;
 996   image_view=AcquireAuthenticCacheView(image,exception);
 997 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 998   #pragma omp parallel for schedule(static,4) shared(progress,status) \
 999     magick_threads(image,image,image->rows,1)
1000 #endif
1001   for (y=0; y < (ssize_t) image->rows; y++)
1002   {
1003     register ssize_t
1004       x;
1005
1006     register Quantum
1007       *magick_restrict q;
1008
1009     if (status == MagickFalse)
1010       continue;
1011     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1012     if (q == (Quantum *) NULL)
1013       {
1014         status=MagickFalse;
1015         continue;
1016       }
1017     for (x=0; x < (ssize_t) image->columns; x++)
1018     {
1019       double
1020         pixel;
1021
1022       register ssize_t
1023         i;
1024
1025       if (GetPixelWriteMask(image,q) == 0)
1026         {
1027           q+=GetPixelChannels(image);
1028           continue;
1029         }
1030       pixel=GetPixelIntensity(image,q);
1031       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1032       {
1033         PixelChannel channel=GetPixelChannelChannel(image,i);
1034         PixelTrait traits=GetPixelChannelTraits(image,channel);
1035         if ((traits & UpdatePixelTrait) == 0)
1036           continue;
1037         if (image->channel_mask != DefaultChannels)
1038           pixel=(double) q[i];
1039         if (pixel < GetPixelInfoChannel(&threshold,channel))
1040           q[i]=(Quantum) 0;
1041       }
1042       q+=GetPixelChannels(image);
1043     }
1044     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1045       status=MagickFalse;
1046     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1047       {
1048         MagickBooleanType
1049           proceed;
1050
1051 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1052         #pragma omp critical (MagickCore_BlackThresholdImage)
1053 #endif
1054         proceed=SetImageProgress(image,ThresholdImageTag,progress++,
1055           image->rows);
1056         if (proceed == MagickFalse)
1057           status=MagickFalse;
1058       }
1059   }
1060   image_view=DestroyCacheView(image_view);
1061   return(status);
1062 }
1063 \f
1064 /*
1065 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1066 %                                                                             %
1067 %                                                                             %
1068 %                                                                             %
1069 %     C l a m p I m a g e                                                     %
1070 %                                                                             %
1071 %                                                                             %
1072 %                                                                             %
1073 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1074 %
1075 %  ClampImage() set each pixel whose value is below zero to zero and any the
1076 %  pixel whose value is above the quantum range to the quantum range (e.g.
1077 %  65535) otherwise the pixel value remains unchanged.
1078 %
1079 %  The format of the ClampImage method is:
1080 %
1081 %      MagickBooleanType ClampImage(Image *image,ExceptionInfo *exception)
1082 %
1083 %  A description of each parameter follows:
1084 %
1085 %    o image: the image.
1086 %
1087 %    o exception: return any errors or warnings in this structure.
1088 %
1089 */
1090
1091 MagickExport MagickBooleanType ClampImage(Image *image,ExceptionInfo *exception)
1092 {
1093 #define ClampImageTag  "Clamp/Image"
1094
1095   CacheView
1096     *image_view;
1097
1098   MagickBooleanType
1099     status;
1100
1101   MagickOffsetType
1102     progress;
1103
1104   ssize_t
1105     y;
1106
1107   assert(image != (Image *) NULL);
1108   assert(image->signature == MagickCoreSignature);
1109   if (image->debug != MagickFalse)
1110     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1111   if (image->storage_class == PseudoClass)
1112     {
1113       register ssize_t
1114         i;
1115
1116       register PixelInfo
1117         *magick_restrict q;
1118
1119       q=image->colormap;
1120       for (i=0; i < (ssize_t) image->colors; i++)
1121       {
1122         q->red=(double) ClampPixel(q->red);
1123         q->green=(double) ClampPixel(q->green);
1124         q->blue=(double) ClampPixel(q->blue);
1125         q->alpha=(double) ClampPixel(q->alpha);
1126         q++;
1127       }
1128       return(SyncImage(image,exception));
1129     }
1130   /*
1131     Clamp image.
1132   */
1133   status=MagickTrue;
1134   progress=0;
1135   image_view=AcquireAuthenticCacheView(image,exception);
1136 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1137   #pragma omp parallel for schedule(static,4) shared(progress,status) \
1138     magick_threads(image,image,image->rows,1)
1139 #endif
1140   for (y=0; y < (ssize_t) image->rows; y++)
1141   {
1142     register ssize_t
1143       x;
1144
1145     register Quantum
1146       *magick_restrict q;
1147
1148     if (status == MagickFalse)
1149       continue;
1150     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1151     if (q == (Quantum *) NULL)
1152       {
1153         status=MagickFalse;
1154         continue;
1155       }
1156     for (x=0; x < (ssize_t) image->columns; x++)
1157     {
1158       register ssize_t
1159         i;
1160
1161       if (GetPixelWriteMask(image,q) == 0)
1162         {
1163           q+=GetPixelChannels(image);
1164           continue;
1165         }
1166       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1167       {
1168         PixelChannel channel=GetPixelChannelChannel(image,i);
1169         PixelTrait traits=GetPixelChannelTraits(image,channel);
1170         if ((traits & UpdatePixelTrait) == 0)
1171           continue;
1172         q[i]=ClampPixel((MagickRealType) q[i]);
1173       }
1174       q+=GetPixelChannels(image);
1175     }
1176     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1177       status=MagickFalse;
1178     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1179       {
1180         MagickBooleanType
1181           proceed;
1182
1183 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1184         #pragma omp critical (MagickCore_ClampImage)
1185 #endif
1186         proceed=SetImageProgress(image,ClampImageTag,progress++,image->rows);
1187         if (proceed == MagickFalse)
1188           status=MagickFalse;
1189       }
1190   }
1191   image_view=DestroyCacheView(image_view);
1192   return(status);
1193 }
1194 \f
1195 /*
1196 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1197 %                                                                             %
1198 %                                                                             %
1199 %                                                                             %
1200 %  D e s t r o y T h r e s h o l d M a p                                      %
1201 %                                                                             %
1202 %                                                                             %
1203 %                                                                             %
1204 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1205 %
1206 %  DestroyThresholdMap() de-allocate the given ThresholdMap
1207 %
1208 %  The format of the ListThresholdMaps method is:
1209 %
1210 %      ThresholdMap *DestroyThresholdMap(Threshold *map)
1211 %
1212 %  A description of each parameter follows.
1213 %
1214 %    o map:    Pointer to the Threshold map to destroy
1215 %
1216 */
1217 MagickExport ThresholdMap *DestroyThresholdMap(ThresholdMap *map)
1218 {
1219   assert(map != (ThresholdMap *) NULL);
1220   if (map->map_id != (char *) NULL)
1221     map->map_id=DestroyString(map->map_id);
1222   if (map->description != (char *) NULL)
1223     map->description=DestroyString(map->description);
1224   if (map->levels != (ssize_t *) NULL)
1225     map->levels=(ssize_t *) RelinquishMagickMemory(map->levels);
1226   map=(ThresholdMap *) RelinquishMagickMemory(map);
1227   return(map);
1228 }
1229 \f
1230 /*
1231 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1232 %                                                                             %
1233 %                                                                             %
1234 %                                                                             %
1235 %  G e t T h r e s h o l d M a p                                              %
1236 %                                                                             %
1237 %                                                                             %
1238 %                                                                             %
1239 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1240 %
1241 %  GetThresholdMap() loads and searches one or more threshold map files for the
1242 %  map matching the given name or alias.
1243 %
1244 %  The format of the GetThresholdMap method is:
1245 %
1246 %      ThresholdMap *GetThresholdMap(const char *map_id,
1247 %        ExceptionInfo *exception)
1248 %
1249 %  A description of each parameter follows.
1250 %
1251 %    o map_id:  ID of the map to look for.
1252 %
1253 %    o exception: return any errors or warnings in this structure.
1254 %
1255 */
1256 MagickExport ThresholdMap *GetThresholdMap(const char *map_id,
1257   ExceptionInfo *exception)
1258 {
1259   ThresholdMap
1260     *map;
1261
1262   map=GetThresholdMapFile(MinimalThresholdMap,"built-in",map_id,exception);
1263   if (map != (ThresholdMap *) NULL)
1264     return(map);
1265 #if !defined(MAGICKCORE_ZERO_CONFIGURATION_SUPPORT)
1266   {
1267     const StringInfo
1268       *option;
1269
1270     LinkedListInfo
1271       *options;
1272
1273     options=GetConfigureOptions(ThresholdsFilename,exception);
1274     option=(const StringInfo *) GetNextValueInLinkedList(options);
1275     while (option != (const StringInfo *) NULL)
1276     {
1277       map=GetThresholdMapFile((const char *) GetStringInfoDatum(option),
1278         GetStringInfoPath(option),map_id,exception);
1279       if (map != (ThresholdMap *) NULL)
1280         break;
1281       option=(const StringInfo *) GetNextValueInLinkedList(options);
1282     }
1283     options=DestroyConfigureOptions(options);
1284   }
1285 #endif
1286   return(map);
1287 }
1288 \f
1289 /*
1290 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1291 %                                                                             %
1292 %                                                                             %
1293 %                                                                             %
1294 +  G e t T h r e s h o l d M a p F i l e                                      %
1295 %                                                                             %
1296 %                                                                             %
1297 %                                                                             %
1298 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1299 %
1300 %  GetThresholdMapFile() look for a given threshold map name or alias in the
1301 %  given XML file data, and return the allocated the map when found.
1302 %
1303 %  The format of the ListThresholdMaps method is:
1304 %
1305 %      ThresholdMap *GetThresholdMap(const char *xml,const char *filename,
1306 %         const char *map_id,ExceptionInfo *exception)
1307 %
1308 %  A description of each parameter follows.
1309 %
1310 %    o xml:  The threshold map list in XML format.
1311 %
1312 %    o filename:  The threshold map XML filename.
1313 %
1314 %    o map_id:  ID of the map to look for in XML list.
1315 %
1316 %    o exception: return any errors or warnings in this structure.
1317 %
1318 */
1319 static ThresholdMap *GetThresholdMapFile(const char *xml,const char *filename,
1320   const char *map_id,ExceptionInfo *exception)
1321 {
1322   char
1323     *p;
1324
1325   const char
1326     *attribute,
1327     *content;
1328
1329   double
1330     value;
1331
1332   register ssize_t
1333     i;
1334
1335   ThresholdMap
1336     *map;
1337
1338   XMLTreeInfo
1339     *description,
1340     *levels,
1341     *threshold,
1342     *thresholds;
1343
1344   (void) LogMagickEvent(ConfigureEvent,GetMagickModule(),
1345     "Loading threshold map file \"%s\" ...",filename);
1346   map=(ThresholdMap *) NULL;
1347   thresholds=NewXMLTree(xml,exception);
1348   if (thresholds == (XMLTreeInfo *) NULL)
1349     return(map);
1350   for (threshold=GetXMLTreeChild(thresholds,"threshold");
1351        threshold != (XMLTreeInfo *) NULL;
1352        threshold=GetNextXMLTreeTag(threshold))
1353   {
1354     attribute=GetXMLTreeAttribute(threshold,"map");
1355     if ((attribute != (char *) NULL) && (LocaleCompare(map_id,attribute) == 0))
1356       break;
1357     attribute=GetXMLTreeAttribute(threshold,"alias");
1358     if ((attribute != (char *) NULL) && (LocaleCompare(map_id,attribute) == 0))
1359       break;
1360   }
1361   if (threshold == (XMLTreeInfo *) NULL)
1362     {
1363       thresholds=DestroyXMLTree(thresholds);
1364       return(map);
1365     }
1366   description=GetXMLTreeChild(threshold,"description");
1367   if (description == (XMLTreeInfo *) NULL)
1368     {
1369       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1370         "XmlMissingElement", "<description>, map \"%s\"",map_id);
1371       thresholds=DestroyXMLTree(thresholds);
1372       return(map);
1373     }
1374   levels=GetXMLTreeChild(threshold,"levels");
1375   if (levels == (XMLTreeInfo *) NULL)
1376     {
1377       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1378         "XmlMissingElement", "<levels>, map \"%s\"", map_id);
1379       thresholds=DestroyXMLTree(thresholds);
1380       return(map);
1381     }
1382   map=(ThresholdMap *) AcquireMagickMemory(sizeof(ThresholdMap));
1383   if (map == (ThresholdMap *) NULL)
1384     ThrowFatalException(ResourceLimitFatalError,"UnableToAcquireThresholdMap");
1385   map->map_id=(char *) NULL;
1386   map->description=(char *) NULL;
1387   map->levels=(ssize_t *) NULL;
1388   attribute=GetXMLTreeAttribute(threshold,"map");
1389   if (attribute != (char *) NULL)
1390     map->map_id=ConstantString(attribute);
1391   content=GetXMLTreeContent(description);
1392   if (content != (char *) NULL)
1393     map->description=ConstantString(content);
1394   attribute=GetXMLTreeAttribute(levels,"width");
1395   if (attribute == (char *) NULL)
1396     {
1397       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1398         "XmlMissingAttribute", "<levels width>, map \"%s\"",map_id);
1399       thresholds=DestroyXMLTree(thresholds);
1400       map=DestroyThresholdMap(map);
1401       return(map);
1402     }
1403   map->width=StringToUnsignedLong(attribute);
1404   if (map->width == 0)
1405     {
1406       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1407        "XmlInvalidAttribute", "<levels width>, map \"%s\"",map_id);
1408       thresholds=DestroyXMLTree(thresholds);
1409       map=DestroyThresholdMap(map);
1410       return(map);
1411     }
1412   attribute=GetXMLTreeAttribute(levels,"height");
1413   if (attribute == (char *) NULL)
1414     {
1415       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1416         "XmlMissingAttribute", "<levels height>, map \"%s\"",map_id);
1417       thresholds=DestroyXMLTree(thresholds);
1418       map=DestroyThresholdMap(map);
1419       return(map);
1420     }
1421   map->height=StringToUnsignedLong(attribute);
1422   if (map->height == 0)
1423     {
1424       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1425         "XmlInvalidAttribute", "<levels height>, map \"%s\"",map_id);
1426       thresholds=DestroyXMLTree(thresholds);
1427       map=DestroyThresholdMap(map);
1428       return(map);
1429     }
1430   attribute=GetXMLTreeAttribute(levels,"divisor");
1431   if (attribute == (char *) NULL)
1432     {
1433       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1434         "XmlMissingAttribute", "<levels divisor>, map \"%s\"",map_id);
1435       thresholds=DestroyXMLTree(thresholds);
1436       map=DestroyThresholdMap(map);
1437       return(map);
1438     }
1439   map->divisor=(ssize_t) StringToLong(attribute);
1440   if (map->divisor < 2)
1441     {
1442       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1443         "XmlInvalidAttribute", "<levels divisor>, map \"%s\"",map_id);
1444       thresholds=DestroyXMLTree(thresholds);
1445       map=DestroyThresholdMap(map);
1446       return(map);
1447     }
1448   content=GetXMLTreeContent(levels);
1449   if (content == (char *) NULL)
1450     {
1451       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1452         "XmlMissingContent", "<levels>, map \"%s\"",map_id);
1453       thresholds=DestroyXMLTree(thresholds);
1454       map=DestroyThresholdMap(map);
1455       return(map);
1456     }
1457   map->levels=(ssize_t *) AcquireQuantumMemory((size_t) map->width,map->height*
1458     sizeof(*map->levels));
1459   if (map->levels == (ssize_t *) NULL)
1460     ThrowFatalException(ResourceLimitFatalError,"UnableToAcquireThresholdMap");
1461   for (i=0; i < (ssize_t) (map->width*map->height); i++)
1462   {
1463     map->levels[i]=(ssize_t) strtol(content,&p,10);
1464     if (p == content)
1465       {
1466         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1467           "XmlInvalidContent", "<level> too few values, map \"%s\"",map_id);
1468         thresholds=DestroyXMLTree(thresholds);
1469         map=DestroyThresholdMap(map);
1470         return(map);
1471       }
1472     if ((map->levels[i] < 0) || (map->levels[i] > map->divisor))
1473       {
1474         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1475           "XmlInvalidContent", "<level> %.20g out of range, map \"%s\"",
1476           (double) map->levels[i],map_id);
1477         thresholds=DestroyXMLTree(thresholds);
1478         map=DestroyThresholdMap(map);
1479         return(map);
1480       }
1481     content=p;
1482   }
1483   value=(double) strtol(content,&p,10);
1484   (void) value;
1485   if (p != content)
1486     {
1487       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1488         "XmlInvalidContent", "<level> too many values, map \"%s\"",map_id);
1489      thresholds=DestroyXMLTree(thresholds);
1490      map=DestroyThresholdMap(map);
1491      return(map);
1492    }
1493   thresholds=DestroyXMLTree(thresholds);
1494   return(map);
1495 }
1496 \f
1497 /*
1498 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1499 %                                                                             %
1500 %                                                                             %
1501 %                                                                             %
1502 +  L i s t T h r e s h o l d M a p F i l e                                    %
1503 %                                                                             %
1504 %                                                                             %
1505 %                                                                             %
1506 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1507 %
1508 %  ListThresholdMapFile() lists the threshold maps and their descriptions
1509 %  in the given XML file data.
1510 %
1511 %  The format of the ListThresholdMaps method is:
1512 %
1513 %      MagickBooleanType ListThresholdMaps(FILE *file,const char*xml,
1514 %         const char *filename,ExceptionInfo *exception)
1515 %
1516 %  A description of each parameter follows.
1517 %
1518 %    o file:  An pointer to the output FILE.
1519 %
1520 %    o xml:  The threshold map list in XML format.
1521 %
1522 %    o filename:  The threshold map XML filename.
1523 %
1524 %    o exception: return any errors or warnings in this structure.
1525 %
1526 */
1527 MagickBooleanType ListThresholdMapFile(FILE *file,const char *xml,
1528   const char *filename,ExceptionInfo *exception)
1529 {
1530   const char
1531     *alias,
1532     *content,
1533     *map;
1534
1535   XMLTreeInfo
1536     *description,
1537     *threshold,
1538     *thresholds;
1539
1540   assert( xml != (char *) NULL );
1541   assert( file != (FILE *) NULL );
1542   (void) LogMagickEvent(ConfigureEvent,GetMagickModule(),
1543     "Loading threshold map file \"%s\" ...",filename);
1544   thresholds=NewXMLTree(xml,exception);
1545   if ( thresholds == (XMLTreeInfo *) NULL )
1546     return(MagickFalse);
1547   (void) FormatLocaleFile(file,"%-16s %-12s %s\n","Map","Alias","Description");
1548   (void) FormatLocaleFile(file,
1549     "----------------------------------------------------\n");
1550   threshold=GetXMLTreeChild(thresholds,"threshold");
1551   for ( ; threshold != (XMLTreeInfo *) NULL;
1552           threshold=GetNextXMLTreeTag(threshold))
1553   {
1554     map=GetXMLTreeAttribute(threshold,"map");
1555     if (map == (char *) NULL)
1556       {
1557         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1558           "XmlMissingAttribute", "<map>");
1559         thresholds=DestroyXMLTree(thresholds);
1560         return(MagickFalse);
1561       }
1562     alias=GetXMLTreeAttribute(threshold,"alias");
1563     description=GetXMLTreeChild(threshold,"description");
1564     if (description == (XMLTreeInfo *) NULL)
1565       {
1566         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1567           "XmlMissingElement", "<description>, map \"%s\"",map);
1568         thresholds=DestroyXMLTree(thresholds);
1569         return(MagickFalse);
1570       }
1571     content=GetXMLTreeContent(description);
1572     if (content == (char *) NULL)
1573       {
1574         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1575           "XmlMissingContent", "<description>, map \"%s\"", map);
1576         thresholds=DestroyXMLTree(thresholds);
1577         return(MagickFalse);
1578       }
1579     (void) FormatLocaleFile(file,"%-16s %-12s %s\n",map,alias ? alias : "",
1580       content);
1581   }
1582   thresholds=DestroyXMLTree(thresholds);
1583   return(MagickTrue);
1584 }
1585 \f
1586 /*
1587 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1588 %                                                                             %
1589 %                                                                             %
1590 %                                                                             %
1591 %  L i s t T h r e s h o l d M a p s                                          %
1592 %                                                                             %
1593 %                                                                             %
1594 %                                                                             %
1595 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1596 %
1597 %  ListThresholdMaps() lists the threshold maps and their descriptions
1598 %  as defined by "threshold.xml" to a file.
1599 %
1600 %  The format of the ListThresholdMaps method is:
1601 %
1602 %      MagickBooleanType ListThresholdMaps(FILE *file,ExceptionInfo *exception)
1603 %
1604 %  A description of each parameter follows.
1605 %
1606 %    o file:  An pointer to the output FILE.
1607 %
1608 %    o exception: return any errors or warnings in this structure.
1609 %
1610 */
1611 MagickExport MagickBooleanType ListThresholdMaps(FILE *file,
1612   ExceptionInfo *exception)
1613 {
1614   const StringInfo
1615     *option;
1616
1617   LinkedListInfo
1618     *options;
1619
1620   MagickStatusType
1621     status;
1622
1623   status=MagickTrue;
1624   if (file == (FILE *) NULL)
1625     file=stdout;
1626   options=GetConfigureOptions(ThresholdsFilename,exception);
1627   (void) FormatLocaleFile(file,
1628     "\n   Threshold Maps for Ordered Dither Operations\n");
1629   option=(const StringInfo *) GetNextValueInLinkedList(options);
1630   while (option != (const StringInfo *) NULL)
1631   {
1632     (void) FormatLocaleFile(file,"\nPath: %s\n\n",GetStringInfoPath(option));
1633     status&=ListThresholdMapFile(file,(const char *) GetStringInfoDatum(option),
1634       GetStringInfoPath(option),exception);
1635     option=(const StringInfo *) GetNextValueInLinkedList(options);
1636   }
1637   options=DestroyConfigureOptions(options);
1638   return(status != 0 ? MagickTrue : MagickFalse);
1639 }
1640 \f
1641 /*
1642 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1643 %                                                                             %
1644 %                                                                             %
1645 %                                                                             %
1646 %     O r d e r e d D i t h e r I m a g e                                     %
1647 %                                                                             %
1648 %                                                                             %
1649 %                                                                             %
1650 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1651 %
1652 %  OrderedDitherImage() will perform a ordered dither based on a number
1653 %  of pre-defined dithering threshold maps, but over multiple intensity
1654 %  levels, which can be different for different channels, according to the
1655 %  input argument.
1656 %
1657 %  The format of the OrderedDitherImage method is:
1658 %
1659 %      MagickBooleanType OrderedDitherImage(Image *image,
1660 %        const char *threshold_map,ExceptionInfo *exception)
1661 %
1662 %  A description of each parameter follows:
1663 %
1664 %    o image: the image.
1665 %
1666 %    o threshold_map: A string containing the name of the threshold dither
1667 %      map to use, followed by zero or more numbers representing the number
1668 %      of color levels tho dither between.
1669 %
1670 %      Any level number less than 2 will be equivalent to 2, and means only
1671 %      binary dithering will be applied to each color channel.
1672 %
1673 %      No numbers also means a 2 level (bitmap) dither will be applied to all
1674 %      channels, while a single number is the number of levels applied to each
1675 %      channel in sequence.  More numbers will be applied in turn to each of
1676 %      the color channels.
1677 %
1678 %      For example: "o3x3,6" will generate a 6 level posterization of the
1679 %      image with a ordered 3x3 diffused pixel dither being applied between
1680 %      each level. While checker,8,8,4 will produce a 332 colormaped image
1681 %      with only a single checkerboard hash pattern (50% grey) between each
1682 %      color level, to basically double the number of color levels with
1683 %      a bare minimim of dithering.
1684 %
1685 %    o exception: return any errors or warnings in this structure.
1686 %
1687 */
1688 MagickExport MagickBooleanType OrderedDitherImage(Image *image,
1689   const char *threshold_map,ExceptionInfo *exception)
1690 {
1691 #define DitherImageTag  "Dither/Image"
1692
1693   CacheView
1694     *image_view;
1695
1696   char
1697     token[MagickPathExtent];
1698
1699   const char
1700     *p;
1701
1702   double
1703     levels[CompositePixelChannel];
1704
1705   MagickBooleanType
1706     status;
1707
1708   MagickOffsetType
1709     progress;
1710
1711   register ssize_t
1712     i;
1713
1714   ssize_t
1715     y;
1716
1717   ThresholdMap
1718     *map;
1719
1720   assert(image != (Image *) NULL);
1721   assert(image->signature == MagickCoreSignature);
1722   if (image->debug != MagickFalse)
1723     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1724   assert(exception != (ExceptionInfo *) NULL);
1725   assert(exception->signature == MagickCoreSignature);
1726   if (threshold_map == (const char *) NULL)
1727     return(MagickTrue);
1728   p=(char *) threshold_map;
1729   while (((isspace((int) ((unsigned char) *p)) != 0) || (*p == ',')) &&
1730          (*p != '\0'))
1731     p++;
1732   threshold_map=p;
1733   while (((isspace((int) ((unsigned char) *p)) == 0) && (*p != ',')) &&
1734          (*p != '\0'))
1735   {
1736     if ((p-threshold_map) >= (MagickPathExtent-1))
1737       break;
1738     token[p-threshold_map]=(*p);
1739     p++;
1740   }
1741   token[p-threshold_map]='\0';
1742   map=GetThresholdMap(token,exception);
1743   if (map == (ThresholdMap *) NULL)
1744     {
1745       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1746         "InvalidArgument","%s : '%s'","ordered-dither",threshold_map);
1747       return(MagickFalse);
1748     }
1749   for (i=0; i < MaxPixelChannels; i++)
1750     levels[i]=2.0;
1751   p=strchr((char *) threshold_map,',');
1752   if ((p != (char *) NULL) && (isdigit((int) ((unsigned char) *(++p))) != 0))
1753     {
1754       GetNextToken(p,&p,MagickPathExtent,token);
1755       for (i=0; (i < MaxPixelChannels); i++)
1756         levels[i]=StringToDouble(token,(char **) NULL);
1757       for (i=0; (*p != '\0') && (i < MaxPixelChannels); i++)
1758       {
1759         GetNextToken(p,&p,MagickPathExtent,token);
1760         if (*token == ',')
1761           GetNextToken(p,&p,MagickPathExtent,token);
1762         levels[i]=StringToDouble(token,(char **) NULL);
1763       }
1764     }
1765   for (i=0; i < MaxPixelChannels; i++)
1766     if (fabs(levels[i]) >= 1)
1767       levels[i]-=1.0;
1768   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
1769     return(MagickFalse);
1770   status=MagickTrue;
1771   progress=0;
1772   image_view=AcquireAuthenticCacheView(image,exception);
1773 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1774   #pragma omp parallel for schedule(static,4) shared(progress,status) \
1775     magick_threads(image,image,image->rows,1)
1776 #endif
1777   for (y=0; y < (ssize_t) image->rows; y++)
1778   {
1779     register ssize_t
1780       x;
1781
1782     register Quantum
1783       *magick_restrict q;
1784
1785     if (status == MagickFalse)
1786       continue;
1787     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1788     if (q == (Quantum *) NULL)
1789       {
1790         status=MagickFalse;
1791         continue;
1792       }
1793     for (x=0; x < (ssize_t) image->columns; x++)
1794     {
1795       register ssize_t
1796         i;
1797
1798       ssize_t
1799         n;
1800
1801       n=0;
1802       if (GetPixelWriteMask(image,q) == 0)
1803         {
1804           q+=GetPixelChannels(image);
1805           continue;
1806         }
1807       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1808       {
1809         ssize_t
1810           level,
1811           threshold;
1812
1813         PixelChannel channel=GetPixelChannelChannel(image,i);
1814         PixelTrait traits=GetPixelChannelTraits(image,channel);
1815         if ((traits & UpdatePixelTrait) == 0)
1816           continue;
1817         if (fabs(levels[n]) < MagickEpsilon)
1818           {
1819             n++;
1820             continue;
1821           }
1822         threshold=(ssize_t) (QuantumScale*q[i]*(levels[n]*(map->divisor-1)+1));
1823         level=threshold/(map->divisor-1);
1824         threshold-=level*(map->divisor-1);
1825         q[i]=ClampToQuantum((double) (level+(threshold >=
1826           map->levels[(x % map->width)+map->width*(y % map->height)]))*
1827           QuantumRange/levels[n]);
1828         n++;
1829       }
1830       q+=GetPixelChannels(image);
1831     }
1832     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1833       status=MagickFalse;
1834     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1835       {
1836         MagickBooleanType
1837           proceed;
1838
1839 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1840         #pragma omp critical (MagickCore_OrderedDitherImage)
1841 #endif
1842         proceed=SetImageProgress(image,DitherImageTag,progress++,image->rows);
1843         if (proceed == MagickFalse)
1844           status=MagickFalse;
1845       }
1846   }
1847   image_view=DestroyCacheView(image_view);
1848   map=DestroyThresholdMap(map);
1849   return(MagickTrue);
1850 }
1851 \f
1852 /*
1853 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1854 %                                                                             %
1855 %                                                                             %
1856 %                                                                             %
1857 %     P e r c e p t i b l e I m a g e                                         %
1858 %                                                                             %
1859 %                                                                             %
1860 %                                                                             %
1861 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1862 %
1863 %  PerceptibleImage() set each pixel whose value is less than |epsilon| to
1864 %  epsilon or -epsilon (whichever is closer) otherwise the pixel value remains
1865 %  unchanged.
1866 %
1867 %  The format of the PerceptibleImage method is:
1868 %
1869 %      MagickBooleanType PerceptibleImage(Image *image,const double epsilon,
1870 %        ExceptionInfo *exception)
1871 %
1872 %  A description of each parameter follows:
1873 %
1874 %    o image: the image.
1875 %
1876 %    o epsilon: the epsilon threshold (e.g. 1.0e-9).
1877 %
1878 %    o exception: return any errors or warnings in this structure.
1879 %
1880 */
1881
1882 static inline Quantum PerceptibleThreshold(const Quantum quantum,
1883   const double epsilon)
1884 {
1885   double
1886     sign;
1887
1888   sign=(double) quantum < 0.0 ? -1.0 : 1.0;
1889   if ((sign*quantum) >= epsilon)
1890     return(quantum);
1891   return((Quantum) (sign*epsilon));
1892 }
1893
1894 MagickExport MagickBooleanType PerceptibleImage(Image *image,
1895   const double epsilon,ExceptionInfo *exception)
1896 {
1897 #define PerceptibleImageTag  "Perceptible/Image"
1898
1899   CacheView
1900     *image_view;
1901
1902   MagickBooleanType
1903     status;
1904
1905   MagickOffsetType
1906     progress;
1907
1908   ssize_t
1909     y;
1910
1911   assert(image != (Image *) NULL);
1912   assert(image->signature == MagickCoreSignature);
1913   if (image->debug != MagickFalse)
1914     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1915   if (image->storage_class == PseudoClass)
1916     {
1917       register ssize_t
1918         i;
1919
1920       register PixelInfo
1921         *magick_restrict q;
1922
1923       q=image->colormap;
1924       for (i=0; i < (ssize_t) image->colors; i++)
1925       {
1926         q->red=(double) PerceptibleThreshold(ClampToQuantum(q->red),
1927           epsilon);
1928         q->green=(double) PerceptibleThreshold(ClampToQuantum(q->green),
1929           epsilon);
1930         q->blue=(double) PerceptibleThreshold(ClampToQuantum(q->blue),
1931           epsilon);
1932         q->alpha=(double) PerceptibleThreshold(ClampToQuantum(q->alpha),
1933           epsilon);
1934         q++;
1935       }
1936       return(SyncImage(image,exception));
1937     }
1938   /*
1939     Perceptible image.
1940   */
1941   status=MagickTrue;
1942   progress=0;
1943   image_view=AcquireAuthenticCacheView(image,exception);
1944 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1945   #pragma omp parallel for schedule(static,4) shared(progress,status) \
1946     magick_threads(image,image,image->rows,1)
1947 #endif
1948   for (y=0; y < (ssize_t) image->rows; y++)
1949   {
1950     register ssize_t
1951       x;
1952
1953     register Quantum
1954       *magick_restrict q;
1955
1956     if (status == MagickFalse)
1957       continue;
1958     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1959     if (q == (Quantum *) NULL)
1960       {
1961         status=MagickFalse;
1962         continue;
1963       }
1964     for (x=0; x < (ssize_t) image->columns; x++)
1965     {
1966       register ssize_t
1967         i;
1968
1969       if (GetPixelWriteMask(image,q) == 0)
1970         {
1971           q+=GetPixelChannels(image);
1972           continue;
1973         }
1974       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1975       {
1976         PixelChannel channel=GetPixelChannelChannel(image,i);
1977         PixelTrait traits=GetPixelChannelTraits(image,channel);
1978         if (traits == UndefinedPixelTrait)
1979           continue;
1980         q[i]=PerceptibleThreshold(q[i],epsilon);
1981       }
1982       q+=GetPixelChannels(image);
1983     }
1984     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1985       status=MagickFalse;
1986     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1987       {
1988         MagickBooleanType
1989           proceed;
1990
1991 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1992         #pragma omp critical (MagickCore_PerceptibleImage)
1993 #endif
1994         proceed=SetImageProgress(image,PerceptibleImageTag,progress++,image->rows);
1995         if (proceed == MagickFalse)
1996           status=MagickFalse;
1997       }
1998   }
1999   image_view=DestroyCacheView(image_view);
2000   return(status);
2001 }
2002 \f
2003 /*
2004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2005 %                                                                             %
2006 %                                                                             %
2007 %                                                                             %
2008 %     R a n d o m T h r e s h o l d I m a g e                                 %
2009 %                                                                             %
2010 %                                                                             %
2011 %                                                                             %
2012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2013 %
2014 %  RandomThresholdImage() changes the value of individual pixels based on the
2015 %  intensity of each pixel compared to a random threshold.  The result is a
2016 %  low-contrast, two color image.
2017 %
2018 %  The format of the RandomThresholdImage method is:
2019 %
2020 %      MagickBooleanType RandomThresholdImage(Image *image,
2021 %        const char *thresholds,ExceptionInfo *exception)
2022 %
2023 %  A description of each parameter follows:
2024 %
2025 %    o image: the image.
2026 %
2027 %    o low,high: Specify the high and low thresholds. These values range from
2028 %      0 to QuantumRange.
2029 %
2030 %    o exception: return any errors or warnings in this structure.
2031 %
2032 */
2033 MagickExport MagickBooleanType RandomThresholdImage(Image *image,
2034   const double min_threshold, const double max_threshold,ExceptionInfo *exception)
2035 {
2036 #define ThresholdImageTag  "Threshold/Image"
2037
2038   CacheView
2039     *image_view;
2040
2041   MagickBooleanType
2042     status;
2043
2044   MagickOffsetType
2045     progress;
2046
2047   PixelInfo
2048     threshold;
2049
2050   RandomInfo
2051     **magick_restrict random_info;
2052
2053   ssize_t
2054     y;
2055
2056 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2057   unsigned long
2058     key;
2059 #endif
2060
2061   assert(image != (Image *) NULL);
2062   assert(image->signature == MagickCoreSignature);
2063   if (image->debug != MagickFalse)
2064     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2065   assert(exception != (ExceptionInfo *) NULL);
2066   assert(exception->signature == MagickCoreSignature);
2067   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2068     return(MagickFalse);
2069   GetPixelInfo(image,&threshold);
2070   /*
2071     Random threshold image.
2072   */
2073   status=MagickTrue;
2074   progress=0;
2075   random_info=AcquireRandomInfoThreadSet();
2076   image_view=AcquireAuthenticCacheView(image,exception);
2077 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2078   key=GetRandomSecretKey(random_info[0]);
2079   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2080     magick_threads(image,image,image->rows,key == ~0UL)
2081 #endif
2082   for (y=0; y < (ssize_t) image->rows; y++)
2083   {
2084     const int
2085       id = GetOpenMPThreadId();
2086
2087     register Quantum
2088       *magick_restrict q;
2089
2090     register ssize_t
2091       x;
2092
2093     if (status == MagickFalse)
2094       continue;
2095     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2096     if (q == (Quantum *) NULL)
2097       {
2098         status=MagickFalse;
2099         continue;
2100       }
2101     for (x=0; x < (ssize_t) image->columns; x++)
2102     {
2103       register ssize_t
2104         i;
2105
2106       if (GetPixelWriteMask(image,q) == 0)
2107         {
2108           q+=GetPixelChannels(image);
2109           continue;
2110         }
2111       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2112       {
2113         double
2114           threshold;
2115
2116         PixelChannel channel=GetPixelChannelChannel(image,i);
2117         PixelTrait traits=GetPixelChannelTraits(image,channel);
2118         if ((traits & UpdatePixelTrait) == 0)
2119           continue;
2120         if ((double) q[i] < min_threshold)
2121           threshold=min_threshold;
2122         else
2123           if ((double) q[i] > max_threshold)
2124             threshold=max_threshold;
2125           else
2126             threshold=(double) (QuantumRange*
2127               GetPseudoRandomValue(random_info[id]));
2128         q[i]=(double) q[i] <= threshold ? 0 : QuantumRange;
2129       }
2130       q+=GetPixelChannels(image);
2131     }
2132     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2133       status=MagickFalse;
2134     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2135       {
2136         MagickBooleanType
2137           proceed;
2138
2139 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2140         #pragma omp critical (MagickCore_RandomThresholdImage)
2141 #endif
2142         proceed=SetImageProgress(image,ThresholdImageTag,progress++,
2143           image->rows);
2144         if (proceed == MagickFalse)
2145           status=MagickFalse;
2146       }
2147   }
2148   image_view=DestroyCacheView(image_view);
2149   random_info=DestroyRandomInfoThreadSet(random_info);
2150   return(status);
2151 }
2152 \f
2153 /*
2154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2155 %                                                                             %
2156 %                                                                             %
2157 %                                                                             %
2158 %     W h i t e T h r e s h o l d I m a g e                                   %
2159 %                                                                             %
2160 %                                                                             %
2161 %                                                                             %
2162 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2163 %
2164 %  WhiteThresholdImage() is like ThresholdImage() but forces all pixels above
2165 %  the threshold into white while leaving all pixels at or below the threshold
2166 %  unchanged.
2167 %
2168 %  The format of the WhiteThresholdImage method is:
2169 %
2170 %      MagickBooleanType WhiteThresholdImage(Image *image,
2171 %        const char *threshold,ExceptionInfo *exception)
2172 %
2173 %  A description of each parameter follows:
2174 %
2175 %    o image: the image.
2176 %
2177 %    o threshold: Define the threshold value.
2178 %
2179 %    o exception: return any errors or warnings in this structure.
2180 %
2181 */
2182 MagickExport MagickBooleanType WhiteThresholdImage(Image *image,
2183   const char *thresholds,ExceptionInfo *exception)
2184 {
2185 #define ThresholdImageTag  "Threshold/Image"
2186
2187   CacheView
2188     *image_view;
2189
2190   GeometryInfo
2191     geometry_info;
2192
2193   MagickBooleanType
2194     status;
2195
2196   MagickOffsetType
2197     progress;
2198
2199   PixelInfo
2200     threshold;
2201
2202   MagickStatusType
2203     flags;
2204
2205   ssize_t
2206     y;
2207
2208   assert(image != (Image *) NULL);
2209   assert(image->signature == MagickCoreSignature);
2210   if (image->debug != MagickFalse)
2211     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2212   if (thresholds == (const char *) NULL)
2213     return(MagickTrue);
2214   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2215     return(MagickFalse);
2216   if (IsGrayColorspace(image->colorspace) != MagickFalse)
2217     (void) TransformImageColorspace(image,sRGBColorspace,exception);
2218   GetPixelInfo(image,&threshold);
2219   flags=ParseGeometry(thresholds,&geometry_info);
2220   threshold.red=geometry_info.rho;
2221   threshold.green=geometry_info.rho;
2222   threshold.blue=geometry_info.rho;
2223   threshold.black=geometry_info.rho;
2224   threshold.alpha=100.0;
2225   if ((flags & SigmaValue) != 0)
2226     threshold.green=geometry_info.sigma;
2227   if ((flags & XiValue) != 0)
2228     threshold.blue=geometry_info.xi;
2229   if ((flags & PsiValue) != 0)
2230     threshold.alpha=geometry_info.psi;
2231   if (threshold.colorspace == CMYKColorspace)
2232     {
2233       if ((flags & PsiValue) != 0)
2234         threshold.black=geometry_info.psi;
2235       if ((flags & ChiValue) != 0)
2236         threshold.alpha=geometry_info.chi;
2237     }
2238   if ((flags & PercentValue) != 0)
2239     {
2240       threshold.red*=(MagickRealType) (QuantumRange/100.0);
2241       threshold.green*=(MagickRealType) (QuantumRange/100.0);
2242       threshold.blue*=(MagickRealType) (QuantumRange/100.0);
2243       threshold.black*=(MagickRealType) (QuantumRange/100.0);
2244       threshold.alpha*=(MagickRealType) (QuantumRange/100.0);
2245     }
2246   /*
2247     White threshold image.
2248   */
2249   status=MagickTrue;
2250   progress=0;
2251   image_view=AcquireAuthenticCacheView(image,exception);
2252 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2253   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2254     magick_threads(image,image,image->rows,1)
2255 #endif
2256   for (y=0; y < (ssize_t) image->rows; y++)
2257   {
2258     register ssize_t
2259       x;
2260
2261     register Quantum
2262       *magick_restrict q;
2263
2264     if (status == MagickFalse)
2265       continue;
2266     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2267     if (q == (Quantum *) NULL)
2268       {
2269         status=MagickFalse;
2270         continue;
2271       }
2272     for (x=0; x < (ssize_t) image->columns; x++)
2273     {
2274       double
2275         pixel;
2276
2277       register ssize_t
2278         i;
2279
2280       if (GetPixelWriteMask(image,q) == 0)
2281         {
2282           q+=GetPixelChannels(image);
2283           continue;
2284         }
2285       pixel=GetPixelIntensity(image,q);
2286       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2287       {
2288         PixelChannel channel=GetPixelChannelChannel(image,i);
2289         PixelTrait traits=GetPixelChannelTraits(image,channel);
2290         if ((traits & UpdatePixelTrait) == 0)
2291           continue;
2292         if (image->channel_mask != DefaultChannels)
2293           pixel=(double) q[i];
2294         if (pixel > GetPixelInfoChannel(&threshold,channel))
2295           q[i]=QuantumRange;
2296       }
2297       q+=GetPixelChannels(image);
2298     }
2299     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2300       status=MagickFalse;
2301     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2302       {
2303         MagickBooleanType
2304           proceed;
2305
2306 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2307         #pragma omp critical (MagickCore_WhiteThresholdImage)
2308 #endif
2309         proceed=SetImageProgress(image,ThresholdImageTag,progress++,
2310           image->rows);
2311         if (proceed == MagickFalse)
2312           status=MagickFalse;
2313       }
2314   }
2315   image_view=DestroyCacheView(image_view);
2316   return(status);
2317 }