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