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

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %              EEEEE  N   N  H   H   AAA   N   N   CCCC  EEEEE                %
   7 %              E      NN  N  H   H  A   A  NN  N  C      E                    %
   8 %              EEE    N N N  HHHHH  AAAAA  N N N  C      EEE                  %
   9 %              E      N  NN  H   H  A   A  N  NN  C      E                    %
  10 %              EEEEE  N   N  H   H  A   A  N   N   CCCC  EEEEE                %
  11 %                                                                             %
  12 %                                                                             %
  13 %                    MagickCore Image Enhancement Methods                     %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                                   Cristy                                    %
  17 %                                 July 1992                                   %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2019 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://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/accelerate-private.h"
  45 #include "MagickCore/artifact.h"
  46 #include "MagickCore/attribute.h"
  47 #include "MagickCore/cache.h"
  48 #include "MagickCore/cache-private.h"
  49 #include "MagickCore/cache-view.h"
  50 #include "MagickCore/channel.h"
  51 #include "MagickCore/color.h"
  52 #include "MagickCore/color-private.h"
  53 #include "MagickCore/colorspace.h"
  54 #include "MagickCore/colorspace-private.h"
  55 #include "MagickCore/composite-private.h"
  56 #include "MagickCore/enhance.h"
  57 #include "MagickCore/exception.h"
  58 #include "MagickCore/exception-private.h"
  59 #include "MagickCore/fx.h"
  60 #include "MagickCore/gem.h"
  61 #include "MagickCore/gem-private.h"
  62 #include "MagickCore/geometry.h"
  63 #include "MagickCore/histogram.h"
  64 #include "MagickCore/image.h"
  65 #include "MagickCore/image-private.h"
  66 #include "MagickCore/memory_.h"
  67 #include "MagickCore/monitor.h"
  68 #include "MagickCore/monitor-private.h"
  69 #include "MagickCore/option.h"
  70 #include "MagickCore/pixel.h"
  71 #include "MagickCore/pixel-accessor.h"
  72 #include "MagickCore/quantum.h"
  73 #include "MagickCore/quantum-private.h"
  74 #include "MagickCore/resample.h"
  75 #include "MagickCore/resample-private.h"
  76 #include "MagickCore/resource_.h"
  77 #include "MagickCore/statistic.h"
  78 #include "MagickCore/string_.h"
  79 #include "MagickCore/string-private.h"
  80 #include "MagickCore/thread-private.h"
  81 #include "MagickCore/threshold.h"
  82 #include "MagickCore/token.h"
  83 #include "MagickCore/xml-tree.h"
  84 #include "MagickCore/xml-tree-private.h"
  85 \f
  86 /*
  87 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  88 %                                                                             %
  89 %                                                                             %
  90 %                                                                             %
  91 %     A u t o G a m m a I m a g e                                             %
  92 %                                                                             %
  93 %                                                                             %
  94 %                                                                             %
  95 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  96 %
  97 %  AutoGammaImage() extract the 'mean' from the image and adjust the image
  98 %  to try make set its gamma appropriatally.
  99 %
 100 %  The format of the AutoGammaImage method is:
 101 %
 102 %      MagickBooleanType AutoGammaImage(Image *image,ExceptionInfo *exception)
 103 %
 104 %  A description of each parameter follows:
 105 %
 106 %    o image: The image to auto-level
 107 %
 108 %    o exception: return any errors or warnings in this structure.
 109 %
 110 */
 111 MagickExport MagickBooleanType AutoGammaImage(Image *image,
 112   ExceptionInfo *exception)
 113 {
 114   double
 115     gamma,
 116     log_mean,
 117     mean,
 118     sans;
 119
 120   MagickStatusType
 121     status;
 122
 123   register ssize_t
 124     i;
 125
 126   log_mean=log(0.5);
 127   if (image->channel_mask == DefaultChannels)
 128     {
 129       /*
 130         Apply gamma correction equally across all given channels.
 131       */
 132       (void) GetImageMean(image,&mean,&sans,exception);
 133       gamma=log(mean*QuantumScale)/log_mean;
 134       return(LevelImage(image,0.0,(double) QuantumRange,gamma,exception));
 135     }
 136   /*
 137     Auto-gamma each channel separately.
 138   */
 139   status=MagickTrue;
 140   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
 141   {
 142     ChannelType
 143       channel_mask;
 144
 145     PixelChannel channel = GetPixelChannelChannel(image,i);
 146     PixelTrait traits = GetPixelChannelTraits(image,channel);
 147     if ((traits & UpdatePixelTrait) == 0)
 148       continue;
 149     channel_mask=SetImageChannelMask(image,(ChannelType) (1UL << i));
 150     status=GetImageMean(image,&mean,&sans,exception);
 151     gamma=log(mean*QuantumScale)/log_mean;
 152     status&=LevelImage(image,0.0,(double) QuantumRange,gamma,exception);
 153     (void) SetImageChannelMask(image,channel_mask);
 154     if (status == MagickFalse)
 155       break;
 156   }
 157   return(status != 0 ? MagickTrue : MagickFalse);
 158 }
 159 \f
 160 /*
 161 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 162 %                                                                             %
 163 %                                                                             %
 164 %                                                                             %
 165 %     A u t o L e v e l I m a g e                                             %
 166 %                                                                             %
 167 %                                                                             %
 168 %                                                                             %
 169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 170 %
 171 %  AutoLevelImage() adjusts the levels of a particular image channel by
 172 %  scaling the minimum and maximum values to the full quantum range.
 173 %
 174 %  The format of the LevelImage method is:
 175 %
 176 %      MagickBooleanType AutoLevelImage(Image *image,ExceptionInfo *exception)
 177 %
 178 %  A description of each parameter follows:
 179 %
 180 %    o image: The image to auto-level
 181 %
 182 %    o exception: return any errors or warnings in this structure.
 183 %
 184 */
 185 MagickExport MagickBooleanType AutoLevelImage(Image *image,
 186   ExceptionInfo *exception)
 187 {
 188   return(MinMaxStretchImage(image,0.0,0.0,1.0,exception));
 189 }
 190 \f
 191 /*
 192 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 193 %                                                                             %
 194 %                                                                             %
 195 %                                                                             %
 196 %     B r i g h t n e s s C o n t r a s t I m a g e                           %
 197 %                                                                             %
 198 %                                                                             %
 199 %                                                                             %
 200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 201 %
 202 %  BrightnessContrastImage() changes the brightness and/or contrast of an
 203 %  image.  It converts the brightness and contrast parameters into slope and
 204 %  intercept and calls a polynomical function to apply to the image.
 205 %
 206 %  The format of the BrightnessContrastImage method is:
 207 %
 208 %      MagickBooleanType BrightnessContrastImage(Image *image,
 209 %        const double brightness,const double contrast,ExceptionInfo *exception)
 210 %
 211 %  A description of each parameter follows:
 212 %
 213 %    o image: the image.
 214 %
 215 %    o brightness: the brightness percent (-100 .. 100).
 216 %
 217 %    o contrast: the contrast percent (-100 .. 100).
 218 %
 219 %    o exception: return any errors or warnings in this structure.
 220 %
 221 */
 222 MagickExport MagickBooleanType BrightnessContrastImage(Image *image,
 223   const double brightness,const double contrast,ExceptionInfo *exception)
 224 {
 225 #define BrightnessContastImageTag  "BrightnessContast/Image"
 226
 227   double
 228     alpha,
 229     coefficients[2],
 230     intercept,
 231     slope;
 232
 233   MagickBooleanType
 234     status;
 235
 236   /*
 237     Compute slope and intercept.
 238   */
 239   assert(image != (Image *) NULL);
 240   assert(image->signature == MagickCoreSignature);
 241   if (image->debug != MagickFalse)
 242     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 243   alpha=contrast;
 244   slope=tan((double) (MagickPI*(alpha/100.0+1.0)/4.0));
 245   if (slope < 0.0)
 246     slope=0.0;
 247   intercept=brightness/100.0+((100-brightness)/200.0)*(1.0-slope);
 248   coefficients[0]=slope;
 249   coefficients[1]=intercept;
 250   status=FunctionImage(image,PolynomialFunction,2,coefficients,exception);
 251   return(status);
 252 }
 253 \f
 254 /*
 255 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 256 %                                                                             %
 257 %                                                                             %
 258 %                                                                             %
 259 %     C L A H E I m a g e                                                     %
 260 %                                                                             %
 261 %                                                                             %
 262 %                                                                             %
 263 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 264 %
 265 %  CLAHEImage() is a variant of adaptive histogram equalization in which the
 266 %  contrast amplification is limited, so as to reduce this problem of noise
 267 %  amplification.
 268 %
 269 %  Adapted from implementation by Karel Zuiderveld, karel@cv.ruu.nl in
 270 %  "Graphics Gems IV", Academic Press, 1994.
 271 %
 272 %  The format of the CLAHEImage method is:
 273 %
 274 %      MagickBooleanType CLAHEImage(Image *image,const size_t width,
 275 %        const size_t height,const size_t number_bins,const double clip_limit,
 276 %        ExceptionInfo *exception)
 277 %
 278 %  A description of each parameter follows:
 279 %
 280 %    o image: the image.
 281 %
 282 %    o width: the width of the tile divisions to use in horizontal direction.
 283 %
 284 %    o height: the height of the tile divisions to use in vertical direction.
 285 %
 286 %    o number_bins: number of bins for histogram ("dynamic range").
 287 %
 288 %    o clip_limit: contrast limit for localised changes in contrast. A limit
 289 %      less than 1 results in standard non-contrast limited AHE.
 290 %
 291 %    o exception: return any errors or warnings in this structure.
 292 %
 293 */
 294
 295 typedef struct _RangeInfo
 296 {
 297   unsigned short
 298     min,
 299     max;
 300 } RangeInfo;
 301
 302 static void ClipCLAHEHistogram(const double clip_limit,const size_t number_bins,
 303   size_t *histogram)
 304 {
 305 #define NumberCLAHEGrays  (65536)
 306
 307   register ssize_t
 308     i;
 309
 310   size_t
 311     cumulative_excess,
 312     previous_excess,
 313     step;
 314
 315   ssize_t
 316     excess;
 317
 318   /*
 319     Compute total number of excess pixels.
 320   */
 321   cumulative_excess=0;
 322   for (i=0; i < (ssize_t) number_bins; i++)
 323   {
 324     excess=(ssize_t) histogram[i]-(ssize_t) clip_limit;
 325     if (excess > 0)
 326       cumulative_excess+=excess;
 327   }
 328   /*
 329     Clip histogram and redistribute excess pixels across all bins.
 330   */
 331   step=cumulative_excess/number_bins;
 332   excess=(ssize_t) (clip_limit-step);
 333   for (i=0; i < (ssize_t) number_bins; i++)
 334   {
 335     if ((double) histogram[i] > clip_limit)
 336       histogram[i]=(size_t) clip_limit;
 337     else
 338       if ((ssize_t) histogram[i] > excess)
 339         {
 340           cumulative_excess-=histogram[i]-excess;
 341           histogram[i]=(size_t) clip_limit;
 342         }
 343       else
 344         {
 345           cumulative_excess-=step;
 346           histogram[i]+=step;
 347         }
 348   }
 349   /*
 350     Redistribute remaining excess.
 351   */
 352   do
 353   {
 354     register size_t
 355       *p;
 356
 357     size_t
 358       *q;
 359
 360     previous_excess=cumulative_excess;
 361     p=histogram;
 362     q=histogram+number_bins;
 363     while ((cumulative_excess != 0) && (p < q))
 364     {
 365       step=number_bins/cumulative_excess;
 366       if (step < 1)
 367         step=1;
 368       for (p=histogram; (p < q) && (cumulative_excess != 0); p+=step)
 369         if ((double) *p < clip_limit)
 370           {
 371             (*p)++;
 372             cumulative_excess--;
 373           }
 374       p++;
 375     }
 376   } while ((cumulative_excess != 0) && (cumulative_excess < previous_excess));
 377 }
 378
 379 static void GenerateCLAHEHistogram(const RectangleInfo *clahe_info,
 380   const RectangleInfo *tile_info,const size_t number_bins,
 381   const unsigned short *lut,const unsigned short *pixels,size_t *histogram)
 382 {
 383   register const unsigned short
 384     *p;
 385
 386   register ssize_t
 387     i;
 388
 389   /*
 390     Classify the pixels into a gray histogram.
 391   */
 392   for (i=0; i < (ssize_t) number_bins; i++)
 393     histogram[i]=0L;
 394   p=pixels;
 395   for (i=0; i < (ssize_t) tile_info->height; i++)
 396   {
 397     const unsigned short
 398       *q;
 399
 400     q=p+tile_info->width;
 401     while (p < q)
 402       histogram[lut[*p++]]++;
 403     q+=clahe_info->width;
 404     p=q-tile_info->width;
 405   }
 406 }
 407
 408 static void InterpolateCLAHE(const RectangleInfo *clahe_info,const size_t *Q12,
 409   const size_t *Q22,const size_t *Q11,const size_t *Q21,
 410   const RectangleInfo *tile,const unsigned short *lut,unsigned short *pixels)
 411 {
 412   ssize_t
 413     y;
 414
 415   unsigned short
 416     intensity;
 417
 418   /*
 419     Bilinear interpolate four tiles to eliminate boundary artifacts.
 420   */
 421   for (y=(ssize_t) tile->height; y > 0; y--)
 422   {
 423     register ssize_t
 424       x;
 425
 426     for (x=(ssize_t) tile->width; x > 0; x--)
 427     {
 428       intensity=lut[*pixels];
 429       *pixels++=(unsigned short ) (PerceptibleReciprocal((double) tile->width*
 430         tile->height)*(y*(x*Q12[intensity]+(tile->width-x)*Q22[intensity])+
 431         (tile->height-y)*(x*Q11[intensity]+(tile->width-x)*Q21[intensity])));
 432     }
 433     pixels+=(clahe_info->width-tile->width);
 434   }
 435 }
 436
 437 static void GenerateCLAHELut(const RangeInfo *range_info,
 438   const size_t number_bins,unsigned short *lut)
 439 {
 440   ssize_t
 441     i;
 442
 443   unsigned short
 444     delta;
 445
 446   /*
 447     Scale input image [intensity min,max] to [0,number_bins-1].
 448   */
 449   delta=(unsigned short) ((range_info->max-range_info->min)/number_bins+1);
 450   for (i=(ssize_t) range_info->min; i <= (ssize_t) range_info->max; i++)
 451     lut[i]=(unsigned short) ((i-range_info->min)/delta);
 452 }
 453
 454 static void MapCLAHEHistogram(const RangeInfo *range_info,
 455   const size_t number_bins,const size_t number_pixels,size_t *histogram)
 456 {
 457   double
 458     scale,
 459     sum;
 460
 461   register ssize_t
 462     i;
 463
 464   /*
 465     Rescale histogram to range [min-intensity .. max-intensity].
 466   */
 467   scale=(double) (range_info->max-range_info->min)/number_pixels;
 468   sum=0.0;
 469   for (i=0; i < (ssize_t) number_bins; i++)
 470   {
 471     sum+=histogram[i];
 472     histogram[i]=(size_t) (range_info->min+scale*sum);
 473     if (histogram[i] > range_info->max)
 474       histogram[i]=range_info->max;
 475   }
 476 }
 477
 478 static MagickBooleanType CLAHE(const RectangleInfo *clahe_info,
 479   const RectangleInfo *tile_info,const RangeInfo *range_info,
 480   const size_t number_bins,const double clip_limit,unsigned short *pixels)
 481 {
 482   MemoryInfo
 483     *tile_cache;
 484
 485   register unsigned short
 486     *p;
 487
 488   size_t
 489     limit,
 490     *tiles;
 491
 492   ssize_t
 493     y;
 494
 495   unsigned short
 496     lut[NumberCLAHEGrays];
 497
 498   /*
 499     Constrast limited adapted histogram equalization.
 500   */
 501   if (clip_limit == 1.0)
 502     return(MagickTrue);
 503   tile_cache=AcquireVirtualMemory((size_t) clahe_info->x*clahe_info->y,
 504     number_bins*sizeof(*tiles));
 505   if (tile_cache == (MemoryInfo *) NULL)
 506     return(MagickFalse);
 507   tiles=(size_t *) GetVirtualMemoryBlob(tile_cache);
 508   limit=(size_t) (clip_limit*(tile_info->width*tile_info->height)/number_bins);
 509   if (limit < 1UL)
 510     limit=1UL;
 511   /*
 512     Generate greylevel mappings for each tile.
 513   */
 514   GenerateCLAHELut(range_info,number_bins,lut);
 515   p=pixels;
 516   for (y=0; y < (ssize_t) clahe_info->y; y++)
 517   {
 518     register ssize_t
 519       x;
 520
 521     for (x=0; x < (ssize_t) clahe_info->x; x++)
 522     {
 523       size_t
 524         *histogram;
 525
 526       histogram=tiles+(number_bins*(y*clahe_info->x+x));
 527       GenerateCLAHEHistogram(clahe_info,tile_info,number_bins,lut,p,histogram);
 528       ClipCLAHEHistogram((double) limit,number_bins,histogram);
 529       MapCLAHEHistogram(range_info,number_bins,tile_info->width*
 530         tile_info->height,histogram);
 531       p+=tile_info->width;
 532     }
 533     p+=clahe_info->width*(tile_info->height-1);
 534   }
 535   /*
 536     Interpolate greylevel mappings to get CLAHE image.
 537   */
 538   p=pixels;
 539   for (y=0; y <= (ssize_t) clahe_info->y; y++)
 540   {
 541     OffsetInfo
 542       offset;
 543
 544     RectangleInfo
 545       tile;
 546
 547     register ssize_t
 548       x;
 549
 550     tile.height=tile_info->height;
 551     tile.y=y-1;
 552     offset.y=tile.y+1;
 553     if (y == 0)
 554       {
 555         /*
 556           Top row.
 557         */
 558         tile.height=tile_info->height >> 1;
 559         tile.y=0;
 560         offset.y=0;
 561       }
 562     else
 563       if (y == (ssize_t) clahe_info->y)
 564         {
 565           /*
 566             Bottom row.
 567           */
 568           tile.height=(tile_info->height+1) >> 1;
 569           tile.y=clahe_info->y-1;
 570           offset.y=tile.y;
 571         }
 572     for (x=0; x <= (ssize_t) clahe_info->x; x++)
 573     {
 574       tile.width=tile_info->width;
 575       tile.x=x-1;
 576       offset.x=tile.x+1;
 577       if (x == 0)
 578         {
 579           /*
 580             Left column.
 581           */
 582           tile.width=tile_info->width >> 1;
 583           tile.x=0;
 584           offset.x=0;
 585         }
 586       else
 587         if (x == (ssize_t) clahe_info->x)
 588           {
 589             /*
 590               Right column.
 591             */
 592             tile.width=(tile_info->width+1) >> 1;
 593             tile.x=clahe_info->x-1;
 594             offset.x=tile.x;
 595           }
 596       InterpolateCLAHE(clahe_info,
 597         tiles+(number_bins*(tile.y*clahe_info->x+tile.x)),     /* Q12 */
 598         tiles+(number_bins*(tile.y*clahe_info->x+offset.x)),   /* Q22 */
 599         tiles+(number_bins*(offset.y*clahe_info->x+tile.x)),   /* Q11 */
 600         tiles+(number_bins*(offset.y*clahe_info->x+offset.x)), /* Q21 */
 601         &tile,lut,p);
 602       p+=tile.width;
 603     }
 604     p+=clahe_info->width*(tile.height-1);
 605   }
 606   tile_cache=RelinquishVirtualMemory(tile_cache);
 607   return(MagickTrue);
 608 }
 609
 610 MagickExport MagickBooleanType CLAHEImage(Image *image,const size_t width,
 611   const size_t height,const size_t number_bins,const double clip_limit,
 612   ExceptionInfo *exception)
 613 {
 614 #define CLAHEImageTag  "CLAHE/Image"
 615
 616   CacheView
 617     *image_view;
 618
 619   ColorspaceType
 620     colorspace;
 621
 622   MagickBooleanType
 623     status;
 624
 625   MagickOffsetType
 626     progress;
 627
 628   MemoryInfo
 629     *pixel_cache;
 630
 631   RangeInfo
 632     range_info;
 633
 634   RectangleInfo
 635     clahe_info,
 636     tile_info;
 637
 638   size_t
 639     n;
 640
 641   ssize_t
 642     y;
 643
 644   unsigned short
 645     *pixels;
 646
 647   /*
 648     Configure CLAHE parameters.
 649   */
 650   assert(image != (Image *) NULL);
 651   assert(image->signature == MagickCoreSignature);
 652   if (image->debug != MagickFalse)
 653     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 654   range_info.min=0;
 655   range_info.max=NumberCLAHEGrays-1;
 656   tile_info.width=width;
 657   if (tile_info.width == 0)
 658     tile_info.width=image->columns >> 3;
 659   tile_info.height=height;
 660   if (tile_info.height == 0)
 661     tile_info.height=image->rows >> 3;
 662   tile_info.x=0;
 663   if ((image->columns % tile_info.width) != 0)
 664     tile_info.x=(ssize_t) tile_info.width-(image->columns % tile_info.width);
 665   tile_info.y=0;
 666   if ((image->rows % tile_info.height) != 0)
 667     tile_info.y=(ssize_t) tile_info.height-(image->rows % tile_info.height);
 668   clahe_info.width=image->columns+tile_info.x;
 669   clahe_info.height=image->rows+tile_info.y;
 670   clahe_info.x=(ssize_t) clahe_info.width/tile_info.width;
 671   clahe_info.y=(ssize_t) clahe_info.height/tile_info.height;
 672   pixel_cache=AcquireVirtualMemory(clahe_info.width,clahe_info.height*
 673     sizeof(*pixels));
 674   if (pixel_cache == (MemoryInfo *) NULL)
 675     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
 676       image->filename);
 677   pixels=(unsigned short *) GetVirtualMemoryBlob(pixel_cache);
 678   colorspace=image->colorspace;
 679   if (TransformImageColorspace(image,LabColorspace,exception) == MagickFalse)
 680     {
 681       pixel_cache=RelinquishVirtualMemory(pixel_cache);
 682       return(MagickFalse);
 683     }
 684   /*
 685     Initialize CLAHE pixels.
 686   */
 687   image_view=AcquireVirtualCacheView(image,exception);
 688   progress=0;
 689   status=MagickTrue;
 690   n=0;
 691   for (y=0; y < (ssize_t) clahe_info.height; y++)
 692   {
 693     register const Quantum
 694       *magick_restrict p;
 695
 696     register ssize_t
 697       x;
 698
 699     if (status == MagickFalse)
 700       continue;
 701     p=GetCacheViewVirtualPixels(image_view,-(tile_info.x >> 1),y-
 702       (tile_info.y >> 1),clahe_info.width,1,exception);
 703     if (p == (const Quantum *) NULL)
 704       {
 705         status=MagickFalse;
 706         continue;
 707       }
 708     for (x=0; x < (ssize_t) clahe_info.width; x++)
 709     {
 710       pixels[n++]=ScaleQuantumToShort(p[0]);
 711       p+=GetPixelChannels(image);
 712     }
 713     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 714       {
 715         MagickBooleanType
 716           proceed;
 717
 718 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 719         #pragma omp atomic
 720 #endif
 721         progress++;
 722         proceed=SetImageProgress(image,CLAHEImageTag,progress,2*
 723           GetPixelChannels(image));
 724         if (proceed == MagickFalse)
 725           status=MagickFalse;
 726       }
 727   }
 728   image_view=DestroyCacheView(image_view);
 729   status=CLAHE(&clahe_info,&tile_info,&range_info,number_bins == 0 ?
 730     (size_t) 128 : MagickMin(number_bins,256),clip_limit,pixels);
 731   if (status == MagickFalse)
 732     (void) ThrowMagickException(exception,GetMagickModule(),
 733       ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 734   /*
 735     Push CLAHE pixels to CLAHE image.
 736   */
 737   image_view=AcquireAuthenticCacheView(image,exception);
 738   n=clahe_info.width*(tile_info.y >> 1);
 739   for (y=0; y < (ssize_t) image->rows; y++)
 740   {
 741     register Quantum
 742       *magick_restrict q;
 743
 744     register ssize_t
 745       x;
 746
 747     if (status == MagickFalse)
 748       continue;
 749     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 750     if (q == (Quantum *) NULL)
 751       {
 752         status=MagickFalse;
 753         continue;
 754       }
 755     n+=tile_info.x >> 1;
 756     for (x=0; x < (ssize_t) image->columns; x++)
 757     {
 758       q[0]=ScaleShortToQuantum(pixels[n++]);
 759       q+=GetPixelChannels(image);
 760     }
 761     n+=(clahe_info.width-image->columns-(tile_info.x >> 1));
 762     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 763       status=MagickFalse;
 764     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 765       {
 766         MagickBooleanType
 767           proceed;
 768
 769 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 770         #pragma omp atomic
 771 #endif
 772         progress++;
 773         proceed=SetImageProgress(image,CLAHEImageTag,progress,2*
 774           GetPixelChannels(image));
 775         if (proceed == MagickFalse)
 776           status=MagickFalse;
 777       }
 778   }
 779   image_view=DestroyCacheView(image_view);
 780   pixel_cache=RelinquishVirtualMemory(pixel_cache);
 781   if (TransformImageColorspace(image,colorspace,exception) == MagickFalse)
 782     status=MagickFalse;
 783   return(status);
 784 }
 785 \f
 786 /*
 787 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 788 %                                                                             %
 789 %                                                                             %
 790 %                                                                             %
 791 %     C l u t I m a g e                                                       %
 792 %                                                                             %
 793 %                                                                             %
 794 %                                                                             %
 795 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 796 %
 797 %  ClutImage() replaces each color value in the given image, by using it as an
 798 %  index to lookup a replacement color value in a Color Look UP Table in the
 799 %  form of an image.  The values are extracted along a diagonal of the CLUT
 800 %  image so either a horizontal or vertial gradient image can be used.
 801 %
 802 %  Typically this is used to either re-color a gray-scale image according to a
 803 %  color gradient in the CLUT image, or to perform a freeform histogram
 804 %  (level) adjustment according to the (typically gray-scale) gradient in the
 805 %  CLUT image.
 806 %
 807 %  When the 'channel' mask includes the matte/alpha transparency channel but
 808 %  one image has no such channel it is assumed that that image is a simple
 809 %  gray-scale image that will effect the alpha channel values, either for
 810 %  gray-scale coloring (with transparent or semi-transparent colors), or
 811 %  a histogram adjustment of existing alpha channel values.   If both images
 812 %  have matte channels, direct and normal indexing is applied, which is rarely
 813 %  used.
 814 %
 815 %  The format of the ClutImage method is:
 816 %
 817 %      MagickBooleanType ClutImage(Image *image,Image *clut_image,
 818 %        const PixelInterpolateMethod method,ExceptionInfo *exception)
 819 %
 820 %  A description of each parameter follows:
 821 %
 822 %    o image: the image, which is replaced by indexed CLUT values
 823 %
 824 %    o clut_image: the color lookup table image for replacement color values.
 825 %
 826 %    o method: the pixel interpolation method.
 827 %
 828 %    o exception: return any errors or warnings in this structure.
 829 %
 830 */
 831 MagickExport MagickBooleanType ClutImage(Image *image,const Image *clut_image,
 832   const PixelInterpolateMethod method,ExceptionInfo *exception)
 833 {
 834 #define ClutImageTag  "Clut/Image"
 835
 836   CacheView
 837     *clut_view,
 838     *image_view;
 839
 840   MagickBooleanType
 841     status;
 842
 843   MagickOffsetType
 844     progress;
 845
 846   PixelInfo
 847     *clut_map;
 848
 849   register ssize_t
 850     i;
 851
 852   ssize_t adjust,
 853     y;
 854
 855   assert(image != (Image *) NULL);
 856   assert(image->signature == MagickCoreSignature);
 857   if (image->debug != MagickFalse)
 858     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 859   assert(clut_image != (Image *) NULL);
 860   assert(clut_image->signature == MagickCoreSignature);
 861   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
 862     return(MagickFalse);
 863   if ((IsGrayColorspace(image->colorspace) != MagickFalse) &&
 864       (IsGrayColorspace(clut_image->colorspace) == MagickFalse))
 865     (void) SetImageColorspace(image,sRGBColorspace,exception);
 866   clut_map=(PixelInfo *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*clut_map));
 867   if (clut_map == (PixelInfo *) NULL)
 868     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
 869       image->filename);
 870   /*
 871     Clut image.
 872   */
 873   status=MagickTrue;
 874   progress=0;
 875   adjust=(ssize_t) (clut_image->interpolate == IntegerInterpolatePixel ? 0 : 1);
 876   clut_view=AcquireVirtualCacheView(clut_image,exception);
 877   for (i=0; i <= (ssize_t) MaxMap; i++)
 878   {
 879     GetPixelInfo(clut_image,clut_map+i);
 880     status=InterpolatePixelInfo(clut_image,clut_view,method,
 881       (double) i*(clut_image->columns-adjust)/MaxMap,(double) i*
 882       (clut_image->rows-adjust)/MaxMap,clut_map+i,exception);
 883     if (status == MagickFalse)
 884       break;
 885   }
 886   clut_view=DestroyCacheView(clut_view);
 887   image_view=AcquireAuthenticCacheView(image,exception);
 888 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 889   #pragma omp parallel for schedule(static) shared(progress,status) \
 890     magick_number_threads(image,image,image->rows,1)
 891 #endif
 892   for (y=0; y < (ssize_t) image->rows; y++)
 893   {
 894     PixelInfo
 895       pixel;
 896
 897     register Quantum
 898       *magick_restrict q;
 899
 900     register ssize_t
 901       x;
 902
 903     if (status == MagickFalse)
 904       continue;
 905     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
 906     if (q == (Quantum *) NULL)
 907       {
 908         status=MagickFalse;
 909         continue;
 910       }
 911     GetPixelInfo(image,&pixel);
 912     for (x=0; x < (ssize_t) image->columns; x++)
 913     {
 914       PixelTrait
 915         traits;
 916
 917       GetPixelInfoPixel(image,q,&pixel);
 918       traits=GetPixelChannelTraits(image,RedPixelChannel);
 919       if ((traits & UpdatePixelTrait) != 0)
 920         pixel.red=clut_map[ScaleQuantumToMap(ClampToQuantum(
 921           pixel.red))].red;
 922       traits=GetPixelChannelTraits(image,GreenPixelChannel);
 923       if ((traits & UpdatePixelTrait) != 0)
 924         pixel.green=clut_map[ScaleQuantumToMap(ClampToQuantum(
 925           pixel.green))].green;
 926       traits=GetPixelChannelTraits(image,BluePixelChannel);
 927       if ((traits & UpdatePixelTrait) != 0)
 928         pixel.blue=clut_map[ScaleQuantumToMap(ClampToQuantum(
 929           pixel.blue))].blue;
 930       traits=GetPixelChannelTraits(image,BlackPixelChannel);
 931       if ((traits & UpdatePixelTrait) != 0)
 932         pixel.black=clut_map[ScaleQuantumToMap(ClampToQuantum(
 933           pixel.black))].black;
 934       traits=GetPixelChannelTraits(image,AlphaPixelChannel);
 935       if ((traits & UpdatePixelTrait) != 0)
 936         pixel.alpha=clut_map[ScaleQuantumToMap(ClampToQuantum(
 937           pixel.alpha))].alpha;
 938       SetPixelViaPixelInfo(image,&pixel,q);
 939       q+=GetPixelChannels(image);
 940     }
 941     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
 942       status=MagickFalse;
 943     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 944       {
 945         MagickBooleanType
 946           proceed;
 947
 948 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 949         #pragma omp atomic
 950 #endif
 951         progress++;
 952         proceed=SetImageProgress(image,ClutImageTag,progress,image->rows);
 953         if (proceed == MagickFalse)
 954           status=MagickFalse;
 955       }
 956   }
 957   image_view=DestroyCacheView(image_view);
 958   clut_map=(PixelInfo *) RelinquishMagickMemory(clut_map);
 959   if ((clut_image->alpha_trait != UndefinedPixelTrait) &&
 960       ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0))
 961     (void) SetImageAlphaChannel(image,ActivateAlphaChannel,exception);
 962   return(status);
 963 }
 964 \f
 965 /*
 966 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 967 %                                                                             %
 968 %                                                                             %
 969 %                                                                             %
 970 %     C o l o r D e c i s i o n L i s t I m a g e                             %
 971 %                                                                             %
 972 %                                                                             %
 973 %                                                                             %
 974 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 975 %
 976 %  ColorDecisionListImage() accepts a lightweight Color Correction Collection
 977 %  (CCC) file which solely contains one or more color corrections and applies
 978 %  the correction to the image.  Here is a sample CCC file:
 979 %
 980 %    <ColorCorrectionCollection xmlns="urn:ASC:CDL:v1.2">
 981 %          <ColorCorrection id="cc03345">
 982 %                <SOPNode>
 983 %                     <Slope> 0.9 1.2 0.5 </Slope>
 984 %                     <Offset> 0.4 -0.5 0.6 </Offset>
 985 %                     <Power> 1.0 0.8 1.5 </Power>
 986 %                </SOPNode>
 987 %                <SATNode>
 988 %                     <Saturation> 0.85 </Saturation>
 989 %                </SATNode>
 990 %          </ColorCorrection>
 991 %    </ColorCorrectionCollection>
 992 %
 993 %  which includes the slop, offset, and power for each of the RGB channels
 994 %  as well as the saturation.
 995 %
 996 %  The format of the ColorDecisionListImage method is:
 997 %
 998 %      MagickBooleanType ColorDecisionListImage(Image *image,
 999 %        const char *color_correction_collection,ExceptionInfo *exception)
1000 %
1001 %  A description of each parameter follows:
1002 %
1003 %    o image: the image.
1004 %
1005 %    o color_correction_collection: the color correction collection in XML.
1006 %
1007 %    o exception: return any errors or warnings in this structure.
1008 %
1009 */
1010 MagickExport MagickBooleanType ColorDecisionListImage(Image *image,
1011   const char *color_correction_collection,ExceptionInfo *exception)
1012 {
1013 #define ColorDecisionListCorrectImageTag  "ColorDecisionList/Image"
1014
1015   typedef struct _Correction
1016   {
1017     double
1018       slope,
1019       offset,
1020       power;
1021   } Correction;
1022
1023   typedef struct _ColorCorrection
1024   {
1025     Correction
1026       red,
1027       green,
1028       blue;
1029
1030     double
1031       saturation;
1032   } ColorCorrection;
1033
1034   CacheView
1035     *image_view;
1036
1037   char
1038     token[MagickPathExtent];
1039
1040   ColorCorrection
1041     color_correction;
1042
1043   const char
1044     *content,
1045     *p;
1046
1047   MagickBooleanType
1048     status;
1049
1050   MagickOffsetType
1051     progress;
1052
1053   PixelInfo
1054     *cdl_map;
1055
1056   register ssize_t
1057     i;
1058
1059   ssize_t
1060     y;
1061
1062   XMLTreeInfo
1063     *cc,
1064     *ccc,
1065     *sat,
1066     *sop;
1067
1068   /*
1069     Allocate and initialize cdl maps.
1070   */
1071   assert(image != (Image *) NULL);
1072   assert(image->signature == MagickCoreSignature);
1073   if (image->debug != MagickFalse)
1074     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1075   if (color_correction_collection == (const char *) NULL)
1076     return(MagickFalse);
1077   ccc=NewXMLTree((const char *) color_correction_collection,exception);
1078   if (ccc == (XMLTreeInfo *) NULL)
1079     return(MagickFalse);
1080   cc=GetXMLTreeChild(ccc,"ColorCorrection");
1081   if (cc == (XMLTreeInfo *) NULL)
1082     {
1083       ccc=DestroyXMLTree(ccc);
1084       return(MagickFalse);
1085     }
1086   color_correction.red.slope=1.0;
1087   color_correction.red.offset=0.0;
1088   color_correction.red.power=1.0;
1089   color_correction.green.slope=1.0;
1090   color_correction.green.offset=0.0;
1091   color_correction.green.power=1.0;
1092   color_correction.blue.slope=1.0;
1093   color_correction.blue.offset=0.0;
1094   color_correction.blue.power=1.0;
1095   color_correction.saturation=0.0;
1096   sop=GetXMLTreeChild(cc,"SOPNode");
1097   if (sop != (XMLTreeInfo *) NULL)
1098     {
1099       XMLTreeInfo
1100         *offset,
1101         *power,
1102         *slope;
1103
1104       slope=GetXMLTreeChild(sop,"Slope");
1105       if (slope != (XMLTreeInfo *) NULL)
1106         {
1107           content=GetXMLTreeContent(slope);
1108           p=(const char *) content;
1109           for (i=0; (*p != '\0') && (i < 3); i++)
1110           {
1111             (void) GetNextToken(p,&p,MagickPathExtent,token);
1112             if (*token == ',')
1113               (void) GetNextToken(p,&p,MagickPathExtent,token);
1114             switch (i)
1115             {
1116               case 0:
1117               {
1118                 color_correction.red.slope=StringToDouble(token,(char **) NULL);
1119                 break;
1120               }
1121               case 1:
1122               {
1123                 color_correction.green.slope=StringToDouble(token,
1124                   (char **) NULL);
1125                 break;
1126               }
1127               case 2:
1128               {
1129                 color_correction.blue.slope=StringToDouble(token,
1130                   (char **) NULL);
1131                 break;
1132               }
1133             }
1134           }
1135         }
1136       offset=GetXMLTreeChild(sop,"Offset");
1137       if (offset != (XMLTreeInfo *) NULL)
1138         {
1139           content=GetXMLTreeContent(offset);
1140           p=(const char *) content;
1141           for (i=0; (*p != '\0') && (i < 3); i++)
1142           {
1143             (void) GetNextToken(p,&p,MagickPathExtent,token);
1144             if (*token == ',')
1145               (void) GetNextToken(p,&p,MagickPathExtent,token);
1146             switch (i)
1147             {
1148               case 0:
1149               {
1150                 color_correction.red.offset=StringToDouble(token,
1151                   (char **) NULL);
1152                 break;
1153               }
1154               case 1:
1155               {
1156                 color_correction.green.offset=StringToDouble(token,
1157                   (char **) NULL);
1158                 break;
1159               }
1160               case 2:
1161               {
1162                 color_correction.blue.offset=StringToDouble(token,
1163                   (char **) NULL);
1164                 break;
1165               }
1166             }
1167           }
1168         }
1169       power=GetXMLTreeChild(sop,"Power");
1170       if (power != (XMLTreeInfo *) NULL)
1171         {
1172           content=GetXMLTreeContent(power);
1173           p=(const char *) content;
1174           for (i=0; (*p != '\0') && (i < 3); i++)
1175           {
1176             (void) GetNextToken(p,&p,MagickPathExtent,token);
1177             if (*token == ',')
1178               (void) GetNextToken(p,&p,MagickPathExtent,token);
1179             switch (i)
1180             {
1181               case 0:
1182               {
1183                 color_correction.red.power=StringToDouble(token,(char **) NULL);
1184                 break;
1185               }
1186               case 1:
1187               {
1188                 color_correction.green.power=StringToDouble(token,
1189                   (char **) NULL);
1190                 break;
1191               }
1192               case 2:
1193               {
1194                 color_correction.blue.power=StringToDouble(token,
1195                   (char **) NULL);
1196                 break;
1197               }
1198             }
1199           }
1200         }
1201     }
1202   sat=GetXMLTreeChild(cc,"SATNode");
1203   if (sat != (XMLTreeInfo *) NULL)
1204     {
1205       XMLTreeInfo
1206         *saturation;
1207
1208       saturation=GetXMLTreeChild(sat,"Saturation");
1209       if (saturation != (XMLTreeInfo *) NULL)
1210         {
1211           content=GetXMLTreeContent(saturation);
1212           p=(const char *) content;
1213           (void) GetNextToken(p,&p,MagickPathExtent,token);
1214           color_correction.saturation=StringToDouble(token,(char **) NULL);
1215         }
1216     }
1217   ccc=DestroyXMLTree(ccc);
1218   if (image->debug != MagickFalse)
1219     {
1220       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1221         "  Color Correction Collection:");
1222       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1223         "  color_correction.red.slope: %g",color_correction.red.slope);
1224       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1225         "  color_correction.red.offset: %g",color_correction.red.offset);
1226       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1227         "  color_correction.red.power: %g",color_correction.red.power);
1228       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1229         "  color_correction.green.slope: %g",color_correction.green.slope);
1230       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1231         "  color_correction.green.offset: %g",color_correction.green.offset);
1232       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1233         "  color_correction.green.power: %g",color_correction.green.power);
1234       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1235         "  color_correction.blue.slope: %g",color_correction.blue.slope);
1236       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1237         "  color_correction.blue.offset: %g",color_correction.blue.offset);
1238       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1239         "  color_correction.blue.power: %g",color_correction.blue.power);
1240       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
1241         "  color_correction.saturation: %g",color_correction.saturation);
1242     }
1243   cdl_map=(PixelInfo *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*cdl_map));
1244   if (cdl_map == (PixelInfo *) NULL)
1245     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1246       image->filename);
1247   for (i=0; i <= (ssize_t) MaxMap; i++)
1248   {
1249     cdl_map[i].red=(double) ScaleMapToQuantum((double)
1250       (MaxMap*(pow(color_correction.red.slope*i/MaxMap+
1251       color_correction.red.offset,color_correction.red.power))));
1252     cdl_map[i].green=(double) ScaleMapToQuantum((double)
1253       (MaxMap*(pow(color_correction.green.slope*i/MaxMap+
1254       color_correction.green.offset,color_correction.green.power))));
1255     cdl_map[i].blue=(double) ScaleMapToQuantum((double)
1256       (MaxMap*(pow(color_correction.blue.slope*i/MaxMap+
1257       color_correction.blue.offset,color_correction.blue.power))));
1258   }
1259   if (image->storage_class == PseudoClass)
1260     for (i=0; i < (ssize_t) image->colors; i++)
1261     {
1262       /*
1263         Apply transfer function to colormap.
1264       */
1265       double
1266         luma;
1267
1268       luma=0.21267f*image->colormap[i].red+0.71526*image->colormap[i].green+
1269         0.07217f*image->colormap[i].blue;
1270       image->colormap[i].red=luma+color_correction.saturation*cdl_map[
1271         ScaleQuantumToMap(ClampToQuantum(image->colormap[i].red))].red-luma;
1272       image->colormap[i].green=luma+color_correction.saturation*cdl_map[
1273         ScaleQuantumToMap(ClampToQuantum(image->colormap[i].green))].green-luma;
1274       image->colormap[i].blue=luma+color_correction.saturation*cdl_map[
1275         ScaleQuantumToMap(ClampToQuantum(image->colormap[i].blue))].blue-luma;
1276     }
1277   /*
1278     Apply transfer function to image.
1279   */
1280   status=MagickTrue;
1281   progress=0;
1282   image_view=AcquireAuthenticCacheView(image,exception);
1283 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1284   #pragma omp parallel for schedule(static) shared(progress,status) \
1285     magick_number_threads(image,image,image->rows,1)
1286 #endif
1287   for (y=0; y < (ssize_t) image->rows; y++)
1288   {
1289     double
1290       luma;
1291
1292     register Quantum
1293       *magick_restrict q;
1294
1295     register ssize_t
1296       x;
1297
1298     if (status == MagickFalse)
1299       continue;
1300     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1301     if (q == (Quantum *) NULL)
1302       {
1303         status=MagickFalse;
1304         continue;
1305       }
1306     for (x=0; x < (ssize_t) image->columns; x++)
1307     {
1308       luma=0.21267f*GetPixelRed(image,q)+0.71526*GetPixelGreen(image,q)+
1309         0.07217f*GetPixelBlue(image,q);
1310       SetPixelRed(image,ClampToQuantum(luma+color_correction.saturation*
1311         (cdl_map[ScaleQuantumToMap(GetPixelRed(image,q))].red-luma)),q);
1312       SetPixelGreen(image,ClampToQuantum(luma+color_correction.saturation*
1313         (cdl_map[ScaleQuantumToMap(GetPixelGreen(image,q))].green-luma)),q);
1314       SetPixelBlue(image,ClampToQuantum(luma+color_correction.saturation*
1315         (cdl_map[ScaleQuantumToMap(GetPixelBlue(image,q))].blue-luma)),q);
1316       q+=GetPixelChannels(image);
1317     }
1318     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1319       status=MagickFalse;
1320     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1321       {
1322         MagickBooleanType
1323           proceed;
1324
1325 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1326         #pragma omp atomic
1327 #endif
1328         progress++;
1329         proceed=SetImageProgress(image,ColorDecisionListCorrectImageTag,
1330           progress,image->rows);
1331         if (proceed == MagickFalse)
1332           status=MagickFalse;
1333       }
1334   }
1335   image_view=DestroyCacheView(image_view);
1336   cdl_map=(PixelInfo *) RelinquishMagickMemory(cdl_map);
1337   return(status);
1338 }
1339 \f
1340 /*
1341 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1342 %                                                                             %
1343 %                                                                             %
1344 %                                                                             %
1345 %     C o n t r a s t I m a g e                                               %
1346 %                                                                             %
1347 %                                                                             %
1348 %                                                                             %
1349 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1350 %
1351 %  ContrastImage() enhances the intensity differences between the lighter and
1352 %  darker elements of the image.  Set sharpen to a MagickTrue to increase the
1353 %  image contrast otherwise the contrast is reduced.
1354 %
1355 %  The format of the ContrastImage method is:
1356 %
1357 %      MagickBooleanType ContrastImage(Image *image,
1358 %        const MagickBooleanType sharpen,ExceptionInfo *exception)
1359 %
1360 %  A description of each parameter follows:
1361 %
1362 %    o image: the image.
1363 %
1364 %    o sharpen: Increase or decrease image contrast.
1365 %
1366 %    o exception: return any errors or warnings in this structure.
1367 %
1368 */
1369
1370 static void Contrast(const int sign,double *red,double *green,double *blue)
1371 {
1372   double
1373     brightness,
1374     hue,
1375     saturation;
1376
1377   /*
1378     Enhance contrast: dark color become darker, light color become lighter.
1379   */
1380   assert(red != (double *) NULL);
1381   assert(green != (double *) NULL);
1382   assert(blue != (double *) NULL);
1383   hue=0.0;
1384   saturation=0.0;
1385   brightness=0.0;
1386   ConvertRGBToHSB(*red,*green,*blue,&hue,&saturation,&brightness);
1387   brightness+=0.5*sign*(0.5*(sin((double) (MagickPI*(brightness-0.5)))+1.0)-
1388     brightness);
1389   if (brightness > 1.0)
1390     brightness=1.0;
1391   else
1392     if (brightness < 0.0)
1393       brightness=0.0;
1394   ConvertHSBToRGB(hue,saturation,brightness,red,green,blue);
1395 }
1396
1397 MagickExport MagickBooleanType ContrastImage(Image *image,
1398   const MagickBooleanType sharpen,ExceptionInfo *exception)
1399 {
1400 #define ContrastImageTag  "Contrast/Image"
1401
1402   CacheView
1403     *image_view;
1404
1405   int
1406     sign;
1407
1408   MagickBooleanType
1409     status;
1410
1411   MagickOffsetType
1412     progress;
1413
1414   register ssize_t
1415     i;
1416
1417   ssize_t
1418     y;
1419
1420   assert(image != (Image *) NULL);
1421   assert(image->signature == MagickCoreSignature);
1422 #if defined(MAGICKCORE_OPENCL_SUPPORT)
1423   if (AccelerateContrastImage(image,sharpen,exception) != MagickFalse)
1424     return(MagickTrue);
1425 #endif
1426   if (image->debug != MagickFalse)
1427     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1428   sign=sharpen != MagickFalse ? 1 : -1;
1429   if (image->storage_class == PseudoClass)
1430     {
1431       /*
1432         Contrast enhance colormap.
1433       */
1434       for (i=0; i < (ssize_t) image->colors; i++)
1435       {
1436         double
1437           blue,
1438           green,
1439           red;
1440
1441         red=(double) image->colormap[i].red;
1442         green=(double) image->colormap[i].green;
1443         blue=(double) image->colormap[i].blue;
1444         Contrast(sign,&red,&green,&blue);
1445         image->colormap[i].red=(MagickRealType) red;
1446         image->colormap[i].green=(MagickRealType) green;
1447         image->colormap[i].blue=(MagickRealType) blue;
1448       }
1449     }
1450   /*
1451     Contrast enhance image.
1452   */
1453   status=MagickTrue;
1454   progress=0;
1455   image_view=AcquireAuthenticCacheView(image,exception);
1456 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1457   #pragma omp parallel for schedule(static) shared(progress,status) \
1458     magick_number_threads(image,image,image->rows,1)
1459 #endif
1460   for (y=0; y < (ssize_t) image->rows; y++)
1461   {
1462     double
1463       blue,
1464       green,
1465       red;
1466
1467     register Quantum
1468       *magick_restrict q;
1469
1470     register ssize_t
1471       x;
1472
1473     if (status == MagickFalse)
1474       continue;
1475     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1476     if (q == (Quantum *) NULL)
1477       {
1478         status=MagickFalse;
1479         continue;
1480       }
1481     for (x=0; x < (ssize_t) image->columns; x++)
1482     {
1483       red=(double) GetPixelRed(image,q);
1484       green=(double) GetPixelGreen(image,q);
1485       blue=(double) GetPixelBlue(image,q);
1486       Contrast(sign,&red,&green,&blue);
1487       SetPixelRed(image,ClampToQuantum(red),q);
1488       SetPixelGreen(image,ClampToQuantum(green),q);
1489       SetPixelBlue(image,ClampToQuantum(blue),q);
1490       q+=GetPixelChannels(image);
1491     }
1492     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1493       status=MagickFalse;
1494     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1495       {
1496         MagickBooleanType
1497           proceed;
1498
1499 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1500         #pragma omp atomic
1501 #endif
1502         progress++;
1503         proceed=SetImageProgress(image,ContrastImageTag,progress,image->rows);
1504         if (proceed == MagickFalse)
1505           status=MagickFalse;
1506       }
1507   }
1508   image_view=DestroyCacheView(image_view);
1509   return(status);
1510 }
1511 \f
1512 /*
1513 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1514 %                                                                             %
1515 %                                                                             %
1516 %                                                                             %
1517 %     C o n t r a s t S t r e t c h I m a g e                                 %
1518 %                                                                             %
1519 %                                                                             %
1520 %                                                                             %
1521 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1522 %
1523 %  ContrastStretchImage() is a simple image enhancement technique that attempts
1524 %  to improve the contrast in an image by 'stretching' the range of intensity
1525 %  values it contains to span a desired range of values. It differs from the
1526 %  more sophisticated histogram equalization in that it can only apply a
1527 %  linear scaling function to the image pixel values.  As a result the
1528 %  'enhancement' is less harsh.
1529 %
1530 %  The format of the ContrastStretchImage method is:
1531 %
1532 %      MagickBooleanType ContrastStretchImage(Image *image,
1533 %        const char *levels,ExceptionInfo *exception)
1534 %
1535 %  A description of each parameter follows:
1536 %
1537 %    o image: the image.
1538 %
1539 %    o black_point: the black point.
1540 %
1541 %    o white_point: the white point.
1542 %
1543 %    o levels: Specify the levels where the black and white points have the
1544 %      range of 0 to number-of-pixels (e.g. 1%, 10x90%, etc.).
1545 %
1546 %    o exception: return any errors or warnings in this structure.
1547 %
1548 */
1549 MagickExport MagickBooleanType ContrastStretchImage(Image *image,
1550   const double black_point,const double white_point,ExceptionInfo *exception)
1551 {
1552 #define MaxRange(color)  ((double) ScaleQuantumToMap((Quantum) (color)))
1553 #define ContrastStretchImageTag  "ContrastStretch/Image"
1554
1555   CacheView
1556     *image_view;
1557
1558   double
1559     *black,
1560     *histogram,
1561     *stretch_map,
1562     *white;
1563
1564   MagickBooleanType
1565     status;
1566
1567   MagickOffsetType
1568     progress;
1569
1570   register ssize_t
1571     i;
1572
1573   ssize_t
1574     y;
1575
1576   /*
1577     Allocate histogram and stretch map.
1578   */
1579   assert(image != (Image *) NULL);
1580   assert(image->signature == MagickCoreSignature);
1581   if (image->debug != MagickFalse)
1582     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1583   if (SetImageGray(image,exception) != MagickFalse)
1584     (void) SetImageColorspace(image,GRAYColorspace,exception);
1585   black=(double *) AcquireQuantumMemory(MaxPixelChannels,sizeof(*black));
1586   white=(double *) AcquireQuantumMemory(MaxPixelChannels,sizeof(*white));
1587   histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,MaxPixelChannels*
1588     sizeof(*histogram));
1589   stretch_map=(double *) AcquireQuantumMemory(MaxMap+1UL,MaxPixelChannels*
1590     sizeof(*stretch_map));
1591   if ((black == (double *) NULL) || (white == (double *) NULL) ||
1592       (histogram == (double *) NULL) || (stretch_map == (double *) NULL))
1593     {
1594       if (stretch_map != (double *) NULL)
1595         stretch_map=(double *) RelinquishMagickMemory(stretch_map);
1596       if (histogram != (double *) NULL)
1597         histogram=(double *) RelinquishMagickMemory(histogram);
1598       if (white != (double *) NULL)
1599         white=(double *) RelinquishMagickMemory(white);
1600       if (black != (double *) NULL)
1601         black=(double *) RelinquishMagickMemory(black);
1602       ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
1603         image->filename);
1604     }
1605   /*
1606     Form histogram.
1607   */
1608   status=MagickTrue;
1609   (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
1610     sizeof(*histogram));
1611   image_view=AcquireVirtualCacheView(image,exception);
1612   for (y=0; y < (ssize_t) image->rows; y++)
1613   {
1614     register const Quantum
1615       *magick_restrict p;
1616
1617     register ssize_t
1618       x;
1619
1620     if (status == MagickFalse)
1621       continue;
1622     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1623     if (p == (const Quantum *) NULL)
1624       {
1625         status=MagickFalse;
1626         continue;
1627       }
1628     for (x=0; x < (ssize_t) image->columns; x++)
1629     {
1630       double
1631         pixel;
1632
1633       pixel=GetPixelIntensity(image,p);
1634       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1635       {
1636         if (image->channel_mask != DefaultChannels)
1637           pixel=(double) p[i];
1638         histogram[GetPixelChannels(image)*ScaleQuantumToMap(
1639           ClampToQuantum(pixel))+i]++;
1640       }
1641       p+=GetPixelChannels(image);
1642     }
1643   }
1644   image_view=DestroyCacheView(image_view);
1645   /*
1646     Find the histogram boundaries by locating the black/white levels.
1647   */
1648   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1649   {
1650     double
1651       intensity;
1652
1653     register ssize_t
1654       j;
1655
1656     black[i]=0.0;
1657     white[i]=MaxRange(QuantumRange);
1658     intensity=0.0;
1659     for (j=0; j <= (ssize_t) MaxMap; j++)
1660     {
1661       intensity+=histogram[GetPixelChannels(image)*j+i];
1662       if (intensity > black_point)
1663         break;
1664     }
1665     black[i]=(double) j;
1666     intensity=0.0;
1667     for (j=(ssize_t) MaxMap; j != 0; j--)
1668     {
1669       intensity+=histogram[GetPixelChannels(image)*j+i];
1670       if (intensity > ((double) image->columns*image->rows-white_point))
1671         break;
1672     }
1673     white[i]=(double) j;
1674   }
1675   histogram=(double *) RelinquishMagickMemory(histogram);
1676   /*
1677     Stretch the histogram to create the stretched image mapping.
1678   */
1679   (void) memset(stretch_map,0,(MaxMap+1)*GetPixelChannels(image)*
1680     sizeof(*stretch_map));
1681   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1682   {
1683     register ssize_t
1684       j;
1685
1686     for (j=0; j <= (ssize_t) MaxMap; j++)
1687     {
1688       double
1689         gamma;
1690
1691       gamma=PerceptibleReciprocal(white[i]-black[i]);
1692       if (j < (ssize_t) black[i])
1693         stretch_map[GetPixelChannels(image)*j+i]=0.0;
1694       else
1695         if (j > (ssize_t) white[i])
1696           stretch_map[GetPixelChannels(image)*j+i]=(double) QuantumRange;
1697         else
1698           if (black[i] != white[i])
1699             stretch_map[GetPixelChannels(image)*j+i]=(double) ScaleMapToQuantum(
1700               (double) (MaxMap*gamma*(j-black[i])));
1701     }
1702   }
1703   if (image->storage_class == PseudoClass)
1704     {
1705       register ssize_t
1706         j;
1707
1708       /*
1709         Stretch-contrast colormap.
1710       */
1711       for (j=0; j < (ssize_t) image->colors; j++)
1712       {
1713         if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
1714           {
1715             i=GetPixelChannelOffset(image,RedPixelChannel);
1716             image->colormap[j].red=stretch_map[GetPixelChannels(image)*
1717               ScaleQuantumToMap(ClampToQuantum(image->colormap[j].red))+i];
1718           }
1719         if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
1720           {
1721             i=GetPixelChannelOffset(image,GreenPixelChannel);
1722             image->colormap[j].green=stretch_map[GetPixelChannels(image)*
1723               ScaleQuantumToMap(ClampToQuantum(image->colormap[j].green))+i];
1724           }
1725         if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
1726           {
1727             i=GetPixelChannelOffset(image,BluePixelChannel);
1728             image->colormap[j].blue=stretch_map[GetPixelChannels(image)*
1729               ScaleQuantumToMap(ClampToQuantum(image->colormap[j].blue))+i];
1730           }
1731         if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
1732           {
1733             i=GetPixelChannelOffset(image,AlphaPixelChannel);
1734             image->colormap[j].alpha=stretch_map[GetPixelChannels(image)*
1735               ScaleQuantumToMap(ClampToQuantum(image->colormap[j].alpha))+i];
1736           }
1737       }
1738     }
1739   /*
1740     Stretch-contrast image.
1741   */
1742   status=MagickTrue;
1743   progress=0;
1744   image_view=AcquireAuthenticCacheView(image,exception);
1745 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1746   #pragma omp parallel for schedule(static) shared(progress,status) \
1747     magick_number_threads(image,image,image->rows,1)
1748 #endif
1749   for (y=0; y < (ssize_t) image->rows; y++)
1750   {
1751     register Quantum
1752       *magick_restrict q;
1753
1754     register ssize_t
1755       x;
1756
1757     if (status == MagickFalse)
1758       continue;
1759     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
1760     if (q == (Quantum *) NULL)
1761       {
1762         status=MagickFalse;
1763         continue;
1764       }
1765     for (x=0; x < (ssize_t) image->columns; x++)
1766     {
1767       register ssize_t
1768         j;
1769
1770       for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
1771       {
1772         PixelChannel channel = GetPixelChannelChannel(image,j);
1773         PixelTrait traits = GetPixelChannelTraits(image,channel);
1774         if ((traits & UpdatePixelTrait) == 0)
1775           continue;
1776         if (black[j] == white[j])
1777           continue;
1778         q[j]=ClampToQuantum(stretch_map[GetPixelChannels(image)*
1779           ScaleQuantumToMap(q[j])+j]);
1780       }
1781       q+=GetPixelChannels(image);
1782     }
1783     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
1784       status=MagickFalse;
1785     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1786       {
1787         MagickBooleanType
1788           proceed;
1789
1790 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1791         #pragma omp atomic
1792 #endif
1793         progress++;
1794         proceed=SetImageProgress(image,ContrastStretchImageTag,progress,
1795           image->rows);
1796         if (proceed == MagickFalse)
1797           status=MagickFalse;
1798       }
1799   }
1800   image_view=DestroyCacheView(image_view);
1801   stretch_map=(double *) RelinquishMagickMemory(stretch_map);
1802   white=(double *) RelinquishMagickMemory(white);
1803   black=(double *) RelinquishMagickMemory(black);
1804   return(status);
1805 }
1806 \f
1807 /*
1808 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1809 %                                                                             %
1810 %                                                                             %
1811 %                                                                             %
1812 %     E n h a n c e I m a g e                                                 %
1813 %                                                                             %
1814 %                                                                             %
1815 %                                                                             %
1816 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1817 %
1818 %  EnhanceImage() applies a digital filter that improves the quality of a
1819 %  noisy image.
1820 %
1821 %  The format of the EnhanceImage method is:
1822 %
1823 %      Image *EnhanceImage(const Image *image,ExceptionInfo *exception)
1824 %
1825 %  A description of each parameter follows:
1826 %
1827 %    o image: the image.
1828 %
1829 %    o exception: return any errors or warnings in this structure.
1830 %
1831 */
1832 MagickExport Image *EnhanceImage(const Image *image,ExceptionInfo *exception)
1833 {
1834 #define EnhanceImageTag  "Enhance/Image"
1835 #define EnhancePixel(weight) \
1836   mean=QuantumScale*((double) GetPixelRed(image,r)+pixel.red)/2.0; \
1837   distance=QuantumScale*((double) GetPixelRed(image,r)-pixel.red); \
1838   distance_squared=(4.0+mean)*distance*distance; \
1839   mean=QuantumScale*((double) GetPixelGreen(image,r)+pixel.green)/2.0; \
1840   distance=QuantumScale*((double) GetPixelGreen(image,r)-pixel.green); \
1841   distance_squared+=(7.0-mean)*distance*distance; \
1842   mean=QuantumScale*((double) GetPixelBlue(image,r)+pixel.blue)/2.0; \
1843   distance=QuantumScale*((double) GetPixelBlue(image,r)-pixel.blue); \
1844   distance_squared+=(5.0-mean)*distance*distance; \
1845   mean=QuantumScale*((double) GetPixelBlack(image,r)+pixel.black)/2.0; \
1846   distance=QuantumScale*((double) GetPixelBlack(image,r)-pixel.black); \
1847   distance_squared+=(5.0-mean)*distance*distance; \
1848   mean=QuantumScale*((double) GetPixelAlpha(image,r)+pixel.alpha)/2.0; \
1849   distance=QuantumScale*((double) GetPixelAlpha(image,r)-pixel.alpha); \
1850   distance_squared+=(5.0-mean)*distance*distance; \
1851   if (distance_squared < 0.069) \
1852     { \
1853       aggregate.red+=(weight)*GetPixelRed(image,r); \
1854       aggregate.green+=(weight)*GetPixelGreen(image,r); \
1855       aggregate.blue+=(weight)*GetPixelBlue(image,r); \
1856       aggregate.black+=(weight)*GetPixelBlack(image,r); \
1857       aggregate.alpha+=(weight)*GetPixelAlpha(image,r); \
1858       total_weight+=(weight); \
1859     } \
1860   r+=GetPixelChannels(image);
1861
1862   CacheView
1863     *enhance_view,
1864     *image_view;
1865
1866   Image
1867     *enhance_image;
1868
1869   MagickBooleanType
1870     status;
1871
1872   MagickOffsetType
1873     progress;
1874
1875   ssize_t
1876     y;
1877
1878   /*
1879     Initialize enhanced image attributes.
1880   */
1881   assert(image != (const Image *) NULL);
1882   assert(image->signature == MagickCoreSignature);
1883   if (image->debug != MagickFalse)
1884     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1885   assert(exception != (ExceptionInfo *) NULL);
1886   assert(exception->signature == MagickCoreSignature);
1887   enhance_image=CloneImage(image,0,0,MagickTrue,
1888     exception);
1889   if (enhance_image == (Image *) NULL)
1890     return((Image *) NULL);
1891   if (SetImageStorageClass(enhance_image,DirectClass,exception) == MagickFalse)
1892     {
1893       enhance_image=DestroyImage(enhance_image);
1894       return((Image *) NULL);
1895     }
1896   /*
1897     Enhance image.
1898   */
1899   status=MagickTrue;
1900   progress=0;
1901   image_view=AcquireVirtualCacheView(image,exception);
1902   enhance_view=AcquireAuthenticCacheView(enhance_image,exception);
1903 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1904   #pragma omp parallel for schedule(static) shared(progress,status) \
1905     magick_number_threads(image,enhance_image,image->rows,1)
1906 #endif
1907   for (y=0; y < (ssize_t) image->rows; y++)
1908   {
1909     PixelInfo
1910       pixel;
1911
1912     register const Quantum
1913       *magick_restrict p;
1914
1915     register Quantum
1916       *magick_restrict q;
1917
1918     register ssize_t
1919       x;
1920
1921     ssize_t
1922       center;
1923
1924     if (status == MagickFalse)
1925       continue;
1926     p=GetCacheViewVirtualPixels(image_view,-2,y-2,image->columns+4,5,exception);
1927     q=QueueCacheViewAuthenticPixels(enhance_view,0,y,enhance_image->columns,1,
1928       exception);
1929     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1930       {
1931         status=MagickFalse;
1932         continue;
1933       }
1934     center=(ssize_t) GetPixelChannels(image)*(2*(image->columns+4)+2);
1935     GetPixelInfo(image,&pixel);
1936     for (x=0; x < (ssize_t) image->columns; x++)
1937     {
1938       double
1939         distance,
1940         distance_squared,
1941         mean,
1942         total_weight;
1943
1944       PixelInfo
1945         aggregate;
1946
1947       register const Quantum
1948         *magick_restrict r;
1949
1950       GetPixelInfo(image,&aggregate);
1951       total_weight=0.0;
1952       GetPixelInfoPixel(image,p+center,&pixel);
1953       r=p;
1954       EnhancePixel(5.0); EnhancePixel(8.0); EnhancePixel(10.0);
1955         EnhancePixel(8.0); EnhancePixel(5.0);
1956       r=p+GetPixelChannels(image)*(image->columns+4);
1957       EnhancePixel(8.0); EnhancePixel(20.0); EnhancePixel(40.0);
1958         EnhancePixel(20.0); EnhancePixel(8.0);
1959       r=p+2*GetPixelChannels(image)*(image->columns+4);
1960       EnhancePixel(10.0); EnhancePixel(40.0); EnhancePixel(80.0);
1961         EnhancePixel(40.0); EnhancePixel(10.0);
1962       r=p+3*GetPixelChannels(image)*(image->columns+4);
1963       EnhancePixel(8.0); EnhancePixel(20.0); EnhancePixel(40.0);
1964         EnhancePixel(20.0); EnhancePixel(8.0);
1965       r=p+4*GetPixelChannels(image)*(image->columns+4);
1966       EnhancePixel(5.0); EnhancePixel(8.0); EnhancePixel(10.0);
1967         EnhancePixel(8.0); EnhancePixel(5.0);
1968       if (total_weight > MagickEpsilon)
1969         {
1970           pixel.red=((aggregate.red+total_weight/2.0)/total_weight);
1971           pixel.green=((aggregate.green+total_weight/2.0)/total_weight);
1972           pixel.blue=((aggregate.blue+total_weight/2.0)/total_weight);
1973           pixel.black=((aggregate.black+total_weight/2.0)/total_weight);
1974           pixel.alpha=((aggregate.alpha+total_weight/2.0)/total_weight);
1975         }
1976       SetPixelViaPixelInfo(enhance_image,&pixel,q);
1977       p+=GetPixelChannels(image);
1978       q+=GetPixelChannels(enhance_image);
1979     }
1980     if (SyncCacheViewAuthenticPixels(enhance_view,exception) == MagickFalse)
1981       status=MagickFalse;
1982     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1983       {
1984         MagickBooleanType
1985           proceed;
1986
1987 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1988         #pragma omp atomic
1989 #endif
1990         progress++;
1991         proceed=SetImageProgress(image,EnhanceImageTag,progress,image->rows);
1992         if (proceed == MagickFalse)
1993           status=MagickFalse;
1994       }
1995   }
1996   enhance_view=DestroyCacheView(enhance_view);
1997   image_view=DestroyCacheView(image_view);
1998   if (status == MagickFalse)
1999     enhance_image=DestroyImage(enhance_image);
2000   return(enhance_image);
2001 }
2002 \f
2003 /*
2004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2005 %                                                                             %
2006 %                                                                             %
2007 %                                                                             %
2008 %     E q u a l i z e I m a g e                                               %
2009 %                                                                             %
2010 %                                                                             %
2011 %                                                                             %
2012 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2013 %
2014 %  EqualizeImage() applies a histogram equalization to the image.
2015 %
2016 %  The format of the EqualizeImage method is:
2017 %
2018 %      MagickBooleanType EqualizeImage(Image *image,ExceptionInfo *exception)
2019 %
2020 %  A description of each parameter follows:
2021 %
2022 %    o image: the image.
2023 %
2024 %    o exception: return any errors or warnings in this structure.
2025 %
2026 */
2027 MagickExport MagickBooleanType EqualizeImage(Image *image,
2028   ExceptionInfo *exception)
2029 {
2030 #define EqualizeImageTag  "Equalize/Image"
2031
2032   CacheView
2033     *image_view;
2034
2035   double
2036     black[CompositePixelChannel+1],
2037     *equalize_map,
2038     *histogram,
2039     *map,
2040     white[CompositePixelChannel+1];
2041
2042   MagickBooleanType
2043     status;
2044
2045   MagickOffsetType
2046     progress;
2047
2048   register ssize_t
2049     i;
2050
2051   ssize_t
2052     y;
2053
2054   /*
2055     Allocate and initialize histogram arrays.
2056   */
2057   assert(image != (Image *) NULL);
2058   assert(image->signature == MagickCoreSignature);
2059 #if defined(MAGICKCORE_OPENCL_SUPPORT)
2060   if (AccelerateEqualizeImage(image,exception) != MagickFalse)
2061     return(MagickTrue);
2062 #endif
2063   if (image->debug != MagickFalse)
2064     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2065   equalize_map=(double *) AcquireQuantumMemory(MaxMap+1UL,MaxPixelChannels*
2066     sizeof(*equalize_map));
2067   histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,MaxPixelChannels*
2068     sizeof(*histogram));
2069   map=(double *) AcquireQuantumMemory(MaxMap+1UL,MaxPixelChannels*sizeof(*map));
2070   if ((equalize_map == (double *) NULL) || (histogram == (double *) NULL) ||
2071       (map == (double *) NULL))
2072     {
2073       if (map != (double *) NULL)
2074         map=(double *) RelinquishMagickMemory(map);
2075       if (histogram != (double *) NULL)
2076         histogram=(double *) RelinquishMagickMemory(histogram);
2077       if (equalize_map != (double *) NULL)
2078         equalize_map=(double *) RelinquishMagickMemory(equalize_map);
2079       ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
2080         image->filename);
2081     }
2082   /*
2083     Form histogram.
2084   */
2085   status=MagickTrue;
2086   (void) memset(histogram,0,(MaxMap+1)*GetPixelChannels(image)*
2087     sizeof(*histogram));
2088   image_view=AcquireVirtualCacheView(image,exception);
2089   for (y=0; y < (ssize_t) image->rows; y++)
2090   {
2091     register const Quantum
2092       *magick_restrict p;
2093
2094     register ssize_t
2095       x;
2096
2097     if (status == MagickFalse)
2098       continue;
2099     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2100     if (p == (const Quantum *) NULL)
2101       {
2102         status=MagickFalse;
2103         continue;
2104       }
2105     for (x=0; x < (ssize_t) image->columns; x++)
2106     {
2107       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2108       {
2109         double
2110           intensity;
2111
2112         intensity=(double) p[i];
2113         if ((image->channel_mask & SyncChannels) != 0)
2114           intensity=GetPixelIntensity(image,p);
2115         histogram[GetPixelChannels(image)*ScaleQuantumToMap(
2116           ClampToQuantum(intensity))+i]++;
2117       }
2118       p+=GetPixelChannels(image);
2119     }
2120   }
2121   image_view=DestroyCacheView(image_view);
2122   /*
2123     Integrate the histogram to get the equalization map.
2124   */
2125   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2126   {
2127     double
2128       intensity;
2129
2130     register ssize_t
2131       j;
2132
2133     intensity=0.0;
2134     for (j=0; j <= (ssize_t) MaxMap; j++)
2135     {
2136       intensity+=histogram[GetPixelChannels(image)*j+i];
2137       map[GetPixelChannels(image)*j+i]=intensity;
2138     }
2139   }
2140   (void) memset(equalize_map,0,(MaxMap+1)*GetPixelChannels(image)*
2141     sizeof(*equalize_map));
2142   (void) memset(black,0,sizeof(*black));
2143   (void) memset(white,0,sizeof(*white));
2144   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2145   {
2146     register ssize_t
2147       j;
2148
2149     black[i]=map[i];
2150     white[i]=map[GetPixelChannels(image)*MaxMap+i];
2151     if (black[i] != white[i])
2152       for (j=0; j <= (ssize_t) MaxMap; j++)
2153         equalize_map[GetPixelChannels(image)*j+i]=(double)
2154           ScaleMapToQuantum((double) ((MaxMap*(map[
2155           GetPixelChannels(image)*j+i]-black[i]))/(white[i]-black[i])));
2156   }
2157   histogram=(double *) RelinquishMagickMemory(histogram);
2158   map=(double *) RelinquishMagickMemory(map);
2159   if (image->storage_class == PseudoClass)
2160     {
2161       register ssize_t
2162         j;
2163
2164       /*
2165         Equalize colormap.
2166       */
2167       for (j=0; j < (ssize_t) image->colors; j++)
2168       {
2169         if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2170           {
2171             PixelChannel channel = GetPixelChannelChannel(image,
2172               RedPixelChannel);
2173             if (black[channel] != white[channel])
2174               image->colormap[j].red=equalize_map[GetPixelChannels(image)*
2175                 ScaleQuantumToMap(ClampToQuantum(image->colormap[j].red))+
2176                 channel];
2177           }
2178         if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2179           {
2180             PixelChannel channel = GetPixelChannelChannel(image,
2181               GreenPixelChannel);
2182             if (black[channel] != white[channel])
2183               image->colormap[j].green=equalize_map[GetPixelChannels(image)*
2184                 ScaleQuantumToMap(ClampToQuantum(image->colormap[j].green))+
2185                 channel];
2186           }
2187         if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2188           {
2189             PixelChannel channel = GetPixelChannelChannel(image,
2190               BluePixelChannel);
2191             if (black[channel] != white[channel])
2192               image->colormap[j].blue=equalize_map[GetPixelChannels(image)*
2193                 ScaleQuantumToMap(ClampToQuantum(image->colormap[j].blue))+
2194                 channel];
2195           }
2196         if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
2197           {
2198             PixelChannel channel = GetPixelChannelChannel(image,
2199               AlphaPixelChannel);
2200             if (black[channel] != white[channel])
2201               image->colormap[j].alpha=equalize_map[GetPixelChannels(image)*
2202                 ScaleQuantumToMap(ClampToQuantum(image->colormap[j].alpha))+
2203                 channel];
2204           }
2205       }
2206     }
2207   /*
2208     Equalize image.
2209   */
2210   progress=0;
2211   image_view=AcquireAuthenticCacheView(image,exception);
2212 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2213   #pragma omp parallel for schedule(static) shared(progress,status) \
2214     magick_number_threads(image,image,image->rows,1)
2215 #endif
2216   for (y=0; y < (ssize_t) image->rows; y++)
2217   {
2218     register Quantum
2219       *magick_restrict q;
2220
2221     register ssize_t
2222       x;
2223
2224     if (status == MagickFalse)
2225       continue;
2226     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2227     if (q == (Quantum *) NULL)
2228       {
2229         status=MagickFalse;
2230         continue;
2231       }
2232     for (x=0; x < (ssize_t) image->columns; x++)
2233     {
2234       register ssize_t
2235         j;
2236
2237       for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
2238       {
2239         PixelChannel channel = GetPixelChannelChannel(image,j);
2240         PixelTrait traits = GetPixelChannelTraits(image,channel);
2241         if (((traits & UpdatePixelTrait) == 0) || (black[j] == white[j]))
2242           continue;
2243         q[j]=ClampToQuantum(equalize_map[GetPixelChannels(image)*
2244           ScaleQuantumToMap(q[j])+j]);
2245       }
2246       q+=GetPixelChannels(image);
2247     }
2248     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2249       status=MagickFalse;
2250     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2251       {
2252         MagickBooleanType
2253           proceed;
2254
2255 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2256         #pragma omp atomic
2257 #endif
2258         progress++;
2259         proceed=SetImageProgress(image,EqualizeImageTag,progress,image->rows);
2260         if (proceed == MagickFalse)
2261           status=MagickFalse;
2262       }
2263   }
2264   image_view=DestroyCacheView(image_view);
2265   equalize_map=(double *) RelinquishMagickMemory(equalize_map);
2266   return(status);
2267 }
2268 \f
2269 /*
2270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2271 %                                                                             %
2272 %                                                                             %
2273 %                                                                             %
2274 %     G a m m a I m a g e                                                     %
2275 %                                                                             %
2276 %                                                                             %
2277 %                                                                             %
2278 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2279 %
2280 %  GammaImage() gamma-corrects a particular image channel.  The same
2281 %  image viewed on different devices will have perceptual differences in the
2282 %  way the image's intensities are represented on the screen.  Specify
2283 %  individual gamma levels for the red, green, and blue channels, or adjust
2284 %  all three with the gamma parameter.  Values typically range from 0.8 to 2.3.
2285 %
2286 %  You can also reduce the influence of a particular channel with a gamma
2287 %  value of 0.
2288 %
2289 %  The format of the GammaImage method is:
2290 %
2291 %      MagickBooleanType GammaImage(Image *image,const double gamma,
2292 %        ExceptionInfo *exception)
2293 %
2294 %  A description of each parameter follows:
2295 %
2296 %    o image: the image.
2297 %
2298 %    o level: the image gamma as a string (e.g. 1.6,1.2,1.0).
2299 %
2300 %    o gamma: the image gamma.
2301 %
2302 */
2303
2304 static inline double gamma_pow(const double value,const double gamma)
2305 {
2306   return(value < 0.0 ? value : pow(value,gamma));
2307 }
2308
2309 MagickExport MagickBooleanType GammaImage(Image *image,const double gamma,
2310   ExceptionInfo *exception)
2311 {
2312 #define GammaImageTag  "Gamma/Image"
2313
2314   CacheView
2315     *image_view;
2316
2317   MagickBooleanType
2318     status;
2319
2320   MagickOffsetType
2321     progress;
2322
2323   Quantum
2324     *gamma_map;
2325
2326   register ssize_t
2327     i;
2328
2329   ssize_t
2330     y;
2331
2332   /*
2333     Allocate and initialize gamma maps.
2334   */
2335   assert(image != (Image *) NULL);
2336   assert(image->signature == MagickCoreSignature);
2337   if (image->debug != MagickFalse)
2338     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2339   if (gamma == 1.0)
2340     return(MagickTrue);
2341   gamma_map=(Quantum *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*gamma_map));
2342   if (gamma_map == (Quantum *) NULL)
2343     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
2344       image->filename);
2345   (void) memset(gamma_map,0,(MaxMap+1)*sizeof(*gamma_map));
2346   if (gamma != 0.0)
2347     for (i=0; i <= (ssize_t) MaxMap; i++)
2348       gamma_map[i]=ScaleMapToQuantum((double) (MaxMap*pow((double) i/
2349         MaxMap,PerceptibleReciprocal(gamma))));
2350   if (image->storage_class == PseudoClass)
2351     for (i=0; i < (ssize_t) image->colors; i++)
2352     {
2353       /*
2354         Gamma-correct colormap.
2355       */
2356       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2357         image->colormap[i].red=(double) gamma_map[ScaleQuantumToMap(
2358           ClampToQuantum(image->colormap[i].red))];
2359       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2360         image->colormap[i].green=(double) gamma_map[ScaleQuantumToMap(
2361           ClampToQuantum(image->colormap[i].green))];
2362       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2363         image->colormap[i].blue=(double) gamma_map[ScaleQuantumToMap(
2364           ClampToQuantum(image->colormap[i].blue))];
2365       if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
2366         image->colormap[i].alpha=(double) gamma_map[ScaleQuantumToMap(
2367           ClampToQuantum(image->colormap[i].alpha))];
2368     }
2369   /*
2370     Gamma-correct image.
2371   */
2372   status=MagickTrue;
2373   progress=0;
2374   image_view=AcquireAuthenticCacheView(image,exception);
2375 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2376   #pragma omp parallel for schedule(static) shared(progress,status) \
2377     magick_number_threads(image,image,image->rows,1)
2378 #endif
2379   for (y=0; y < (ssize_t) image->rows; y++)
2380   {
2381     register Quantum
2382       *magick_restrict q;
2383
2384     register ssize_t
2385       x;
2386
2387     if (status == MagickFalse)
2388       continue;
2389     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2390     if (q == (Quantum *) NULL)
2391       {
2392         status=MagickFalse;
2393         continue;
2394       }
2395     for (x=0; x < (ssize_t) image->columns; x++)
2396     {
2397       register ssize_t
2398         j;
2399
2400       for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
2401       {
2402         PixelChannel channel = GetPixelChannelChannel(image,j);
2403         PixelTrait traits = GetPixelChannelTraits(image,channel);
2404         if ((traits & UpdatePixelTrait) == 0)
2405           continue;
2406         q[j]=gamma_map[ScaleQuantumToMap(ClampToQuantum((MagickRealType)
2407           q[j]))];
2408       }
2409       q+=GetPixelChannels(image);
2410     }
2411     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2412       status=MagickFalse;
2413     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2414       {
2415         MagickBooleanType
2416           proceed;
2417
2418 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2419         #pragma omp atomic
2420 #endif
2421         progress++;
2422         proceed=SetImageProgress(image,GammaImageTag,progress,image->rows);
2423         if (proceed == MagickFalse)
2424           status=MagickFalse;
2425       }
2426   }
2427   image_view=DestroyCacheView(image_view);
2428   gamma_map=(Quantum *) RelinquishMagickMemory(gamma_map);
2429   if (image->gamma != 0.0)
2430     image->gamma*=gamma;
2431   return(status);
2432 }
2433 \f
2434 /*
2435 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2436 %                                                                             %
2437 %                                                                             %
2438 %                                                                             %
2439 %     G r a y s c a l e I m a g e                                             %
2440 %                                                                             %
2441 %                                                                             %
2442 %                                                                             %
2443 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2444 %
2445 %  GrayscaleImage() converts the image to grayscale.
2446 %
2447 %  The format of the GrayscaleImage method is:
2448 %
2449 %      MagickBooleanType GrayscaleImage(Image *image,
2450 %        const PixelIntensityMethod method ,ExceptionInfo *exception)
2451 %
2452 %  A description of each parameter follows:
2453 %
2454 %    o image: the image.
2455 %
2456 %    o method: the pixel intensity method.
2457 %
2458 %    o exception: return any errors or warnings in this structure.
2459 %
2460 */
2461 MagickExport MagickBooleanType GrayscaleImage(Image *image,
2462   const PixelIntensityMethod method,ExceptionInfo *exception)
2463 {
2464 #define GrayscaleImageTag  "Grayscale/Image"
2465
2466   CacheView
2467     *image_view;
2468
2469   MagickBooleanType
2470     status;
2471
2472   MagickOffsetType
2473     progress;
2474
2475   ssize_t
2476     y;
2477
2478   assert(image != (Image *) NULL);
2479   assert(image->signature == MagickCoreSignature);
2480   if (image->debug != MagickFalse)
2481     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2482   if (image->storage_class == PseudoClass)
2483     {
2484       if (SyncImage(image,exception) == MagickFalse)
2485         return(MagickFalse);
2486       if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2487         return(MagickFalse);
2488     }
2489 #if defined(MAGICKCORE_OPENCL_SUPPORT)
2490   if (AccelerateGrayscaleImage(image,method,exception) != MagickFalse)
2491     {
2492       image->intensity=method;
2493       image->type=GrayscaleType;
2494       if ((method == Rec601LuminancePixelIntensityMethod) ||
2495           (method == Rec709LuminancePixelIntensityMethod))
2496         return(SetImageColorspace(image,LinearGRAYColorspace,exception));
2497       return(SetImageColorspace(image,GRAYColorspace,exception));
2498     }
2499 #endif
2500   /*
2501     Grayscale image.
2502   */
2503   status=MagickTrue;
2504   progress=0;
2505   image_view=AcquireAuthenticCacheView(image,exception);
2506 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2507   #pragma omp parallel for schedule(static) shared(progress,status) \
2508     magick_number_threads(image,image,image->rows,1)
2509 #endif
2510   for (y=0; y < (ssize_t) image->rows; y++)
2511   {
2512     register Quantum
2513       *magick_restrict q;
2514
2515     register ssize_t
2516       x;
2517
2518     if (status == MagickFalse)
2519       continue;
2520     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2521     if (q == (Quantum *) NULL)
2522       {
2523         status=MagickFalse;
2524         continue;
2525       }
2526     for (x=0; x < (ssize_t) image->columns; x++)
2527     {
2528       MagickRealType
2529         blue,
2530         green,
2531         red,
2532         intensity;
2533
2534       red=(MagickRealType) GetPixelRed(image,q);
2535       green=(MagickRealType) GetPixelGreen(image,q);
2536       blue=(MagickRealType) GetPixelBlue(image,q);
2537       intensity=0.0;
2538       switch (method)
2539       {
2540         case AveragePixelIntensityMethod:
2541         {
2542           intensity=(red+green+blue)/3.0;
2543           break;
2544         }
2545         case BrightnessPixelIntensityMethod:
2546         {
2547           intensity=MagickMax(MagickMax(red,green),blue);
2548           break;
2549         }
2550         case LightnessPixelIntensityMethod:
2551         {
2552           intensity=(MagickMin(MagickMin(red,green),blue)+
2553             MagickMax(MagickMax(red,green),blue))/2.0;
2554           break;
2555         }
2556         case MSPixelIntensityMethod:
2557         {
2558           intensity=(MagickRealType) (((double) red*red+green*green+
2559             blue*blue)/3.0);
2560           break;
2561         }
2562         case Rec601LumaPixelIntensityMethod:
2563         {
2564           if (image->colorspace == RGBColorspace)
2565             {
2566               red=EncodePixelGamma(red);
2567               green=EncodePixelGamma(green);
2568               blue=EncodePixelGamma(blue);
2569             }
2570           intensity=0.298839*red+0.586811*green+0.114350*blue;
2571           break;
2572         }
2573         case Rec601LuminancePixelIntensityMethod:
2574         {
2575           if (image->colorspace == sRGBColorspace)
2576             {
2577               red=DecodePixelGamma(red);
2578               green=DecodePixelGamma(green);
2579               blue=DecodePixelGamma(blue);
2580             }
2581           intensity=0.298839*red+0.586811*green+0.114350*blue;
2582           break;
2583         }
2584         case Rec709LumaPixelIntensityMethod:
2585         default:
2586         {
2587           if (image->colorspace == RGBColorspace)
2588             {
2589               red=EncodePixelGamma(red);
2590               green=EncodePixelGamma(green);
2591               blue=EncodePixelGamma(blue);
2592             }
2593           intensity=0.212656*red+0.715158*green+0.072186*blue;
2594           break;
2595         }
2596         case Rec709LuminancePixelIntensityMethod:
2597         {
2598           if (image->colorspace == sRGBColorspace)
2599             {
2600               red=DecodePixelGamma(red);
2601               green=DecodePixelGamma(green);
2602               blue=DecodePixelGamma(blue);
2603             }
2604           intensity=0.212656*red+0.715158*green+0.072186*blue;
2605           break;
2606         }
2607         case RMSPixelIntensityMethod:
2608         {
2609           intensity=(MagickRealType) (sqrt((double) red*red+green*green+
2610             blue*blue)/sqrt(3.0));
2611           break;
2612         }
2613       }
2614       SetPixelGray(image,ClampToQuantum(intensity),q);
2615       q+=GetPixelChannels(image);
2616     }
2617     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2618       status=MagickFalse;
2619     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2620       {
2621         MagickBooleanType
2622           proceed;
2623
2624 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2625         #pragma omp atomic
2626 #endif
2627         progress++;
2628         proceed=SetImageProgress(image,GrayscaleImageTag,progress,image->rows);
2629         if (proceed == MagickFalse)
2630           status=MagickFalse;
2631       }
2632   }
2633   image_view=DestroyCacheView(image_view);
2634   image->intensity=method;
2635   image->type=GrayscaleType;
2636   if ((method == Rec601LuminancePixelIntensityMethod) ||
2637       (method == Rec709LuminancePixelIntensityMethod))
2638     return(SetImageColorspace(image,LinearGRAYColorspace,exception));
2639   return(SetImageColorspace(image,GRAYColorspace,exception));
2640 }
2641 \f
2642 /*
2643 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2644 %                                                                             %
2645 %                                                                             %
2646 %                                                                             %
2647 %     H a l d C l u t I m a g e                                               %
2648 %                                                                             %
2649 %                                                                             %
2650 %                                                                             %
2651 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2652 %
2653 %  HaldClutImage() applies a Hald color lookup table to the image.  A Hald
2654 %  color lookup table is a 3-dimensional color cube mapped to 2 dimensions.
2655 %  Create it with the HALD coder.  You can apply any color transformation to
2656 %  the Hald image and then use this method to apply the transform to the
2657 %  image.
2658 %
2659 %  The format of the HaldClutImage method is:
2660 %
2661 %      MagickBooleanType HaldClutImage(Image *image,Image *hald_image,
2662 %        ExceptionInfo *exception)
2663 %
2664 %  A description of each parameter follows:
2665 %
2666 %    o image: the image, which is replaced by indexed CLUT values
2667 %
2668 %    o hald_image: the color lookup table image for replacement color values.
2669 %
2670 %    o exception: return any errors or warnings in this structure.
2671 %
2672 */
2673 MagickExport MagickBooleanType HaldClutImage(Image *image,
2674   const Image *hald_image,ExceptionInfo *exception)
2675 {
2676 #define HaldClutImageTag  "Clut/Image"
2677
2678   typedef struct _HaldInfo
2679   {
2680     double
2681       x,
2682       y,
2683       z;
2684   } HaldInfo;
2685
2686   CacheView
2687     *hald_view,
2688     *image_view;
2689
2690   double
2691     width;
2692
2693   MagickBooleanType
2694     status;
2695
2696   MagickOffsetType
2697     progress;
2698
2699   PixelInfo
2700     zero;
2701
2702   size_t
2703     cube_size,
2704     length,
2705     level;
2706
2707   ssize_t
2708     y;
2709
2710   assert(image != (Image *) NULL);
2711   assert(image->signature == MagickCoreSignature);
2712   if (image->debug != MagickFalse)
2713     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2714   assert(hald_image != (Image *) NULL);
2715   assert(hald_image->signature == MagickCoreSignature);
2716   if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse)
2717     return(MagickFalse);
2718   if (image->alpha_trait == UndefinedPixelTrait)
2719     (void) SetImageAlphaChannel(image,OpaqueAlphaChannel,exception);
2720   /*
2721     Hald clut image.
2722   */
2723   status=MagickTrue;
2724   progress=0;
2725   length=(size_t) MagickMin((MagickRealType) hald_image->columns,
2726     (MagickRealType) hald_image->rows);
2727   for (level=2; (level*level*level) < length; level++) ;
2728   level*=level;
2729   cube_size=level*level;
2730   width=(double) hald_image->columns;
2731   GetPixelInfo(hald_image,&zero);
2732   hald_view=AcquireVirtualCacheView(hald_image,exception);
2733   image_view=AcquireAuthenticCacheView(image,exception);
2734 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2735   #pragma omp parallel for schedule(static) shared(progress,status) \
2736     magick_number_threads(image,image,image->rows,1)
2737 #endif
2738   for (y=0; y < (ssize_t) image->rows; y++)
2739   {
2740     register Quantum
2741       *magick_restrict q;
2742
2743     register ssize_t
2744       x;
2745
2746     if (status == MagickFalse)
2747       continue;
2748     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2749     if (q == (Quantum *) NULL)
2750       {
2751         status=MagickFalse;
2752         continue;
2753       }
2754     for (x=0; x < (ssize_t) image->columns; x++)
2755     {
2756       double
2757         offset,
2758         area;
2759
2760       HaldInfo
2761         point;
2762
2763       PixelInfo
2764         pixel,
2765         pixel1,
2766         pixel2,
2767         pixel3,
2768         pixel4;
2769
2770       point.x=QuantumScale*(level-1.0)*GetPixelRed(image,q);
2771       point.y=QuantumScale*(level-1.0)*GetPixelGreen(image,q);
2772       point.z=QuantumScale*(level-1.0)*GetPixelBlue(image,q);
2773       offset=point.x+level*floor(point.y)+cube_size*floor(point.z);
2774       point.x-=floor(point.x);
2775       point.y-=floor(point.y);
2776       point.z-=floor(point.z);
2777       pixel1=zero;
2778       status=InterpolatePixelInfo(hald_image,hald_view,hald_image->interpolate,
2779         fmod(offset,width),floor(offset/width),&pixel1,exception);
2780       if (status == MagickFalse)
2781         break;
2782       pixel2=zero;
2783       status=InterpolatePixelInfo(hald_image,hald_view,hald_image->interpolate,
2784         fmod(offset+level,width),floor((offset+level)/width),&pixel2,exception);
2785       if (status == MagickFalse)
2786         break;
2787       pixel3=zero;
2788       area=point.y;
2789       if (hald_image->interpolate == NearestInterpolatePixel)
2790         area=(point.y < 0.5) ? 0 : 1;
2791       CompositePixelInfoAreaBlend(&pixel1,pixel1.alpha,&pixel2,pixel2.alpha,
2792         area,&pixel3);
2793       offset+=cube_size;
2794       status=InterpolatePixelInfo(hald_image,hald_view,hald_image->interpolate,
2795         fmod(offset,width),floor(offset/width),&pixel1,exception);
2796       if (status == MagickFalse)
2797         break;
2798       status=InterpolatePixelInfo(hald_image,hald_view,hald_image->interpolate,
2799         fmod(offset+level,width),floor((offset+level)/width),&pixel2,exception);
2800       if (status == MagickFalse)
2801         break;
2802       pixel4=zero;
2803       CompositePixelInfoAreaBlend(&pixel1,pixel1.alpha,&pixel2,pixel2.alpha,
2804         area,&pixel4);
2805       pixel=zero;
2806       area=point.z;
2807       if (hald_image->interpolate==NearestInterpolatePixel)
2808         area=(point.z<0.5)?0:1;
2809       CompositePixelInfoAreaBlend(&pixel3,pixel3.alpha,&pixel4,pixel4.alpha,
2810         area,&pixel);
2811       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2812         SetPixelRed(image,ClampToQuantum(pixel.red),q);
2813       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2814         SetPixelGreen(image,ClampToQuantum(pixel.green),q);
2815       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2816         SetPixelBlue(image,ClampToQuantum(pixel.blue),q);
2817       if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2818           (image->colorspace == CMYKColorspace))
2819         SetPixelBlack(image,ClampToQuantum(pixel.black),q);
2820       if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2821           (image->alpha_trait != UndefinedPixelTrait))
2822         SetPixelAlpha(image,ClampToQuantum(pixel.alpha),q);
2823       q+=GetPixelChannels(image);
2824     }
2825     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2826       status=MagickFalse;
2827     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2828       {
2829         MagickBooleanType
2830           proceed;
2831
2832 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2833         #pragma omp atomic
2834 #endif
2835         progress++;
2836         proceed=SetImageProgress(image,HaldClutImageTag,progress,image->rows);
2837         if (proceed == MagickFalse)
2838           status=MagickFalse;
2839       }
2840   }
2841   hald_view=DestroyCacheView(hald_view);
2842   image_view=DestroyCacheView(image_view);
2843   return(status);
2844 }
2845 \f
2846 /*
2847 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2848 %                                                                             %
2849 %                                                                             %
2850 %                                                                             %
2851 %     L e v e l I m a g e                                                     %
2852 %                                                                             %
2853 %                                                                             %
2854 %                                                                             %
2855 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2856 %
2857 %  LevelImage() adjusts the levels of a particular image channel by
2858 %  scaling the colors falling between specified white and black points to
2859 %  the full available quantum range.
2860 %
2861 %  The parameters provided represent the black, and white points.  The black
2862 %  point specifies the darkest color in the image. Colors darker than the
2863 %  black point are set to zero.  White point specifies the lightest color in
2864 %  the image.  Colors brighter than the white point are set to the maximum
2865 %  quantum value.
2866 %
2867 %  If a '!' flag is given, map black and white colors to the given levels
2868 %  rather than mapping those levels to black and white.  See
2869 %  LevelizeImage() below.
2870 %
2871 %  Gamma specifies a gamma correction to apply to the image.
2872 %
2873 %  The format of the LevelImage method is:
2874 %
2875 %      MagickBooleanType LevelImage(Image *image,const double black_point,
2876 %        const double white_point,const double gamma,ExceptionInfo *exception)
2877 %
2878 %  A description of each parameter follows:
2879 %
2880 %    o image: the image.
2881 %
2882 %    o black_point: The level to map zero (black) to.
2883 %
2884 %    o white_point: The level to map QuantumRange (white) to.
2885 %
2886 %    o exception: return any errors or warnings in this structure.
2887 %
2888 */
2889
2890 static inline double LevelPixel(const double black_point,
2891   const double white_point,const double gamma,const double pixel)
2892 {
2893   double
2894     level_pixel,
2895     scale;
2896
2897   scale=PerceptibleReciprocal(white_point-black_point);
2898   level_pixel=QuantumRange*gamma_pow(scale*((double) pixel-black_point),
2899     PerceptibleReciprocal(gamma));
2900   return(level_pixel);
2901 }
2902
2903 MagickExport MagickBooleanType LevelImage(Image *image,const double black_point,
2904   const double white_point,const double gamma,ExceptionInfo *exception)
2905 {
2906 #define LevelImageTag  "Level/Image"
2907
2908   CacheView
2909     *image_view;
2910
2911   MagickBooleanType
2912     status;
2913
2914   MagickOffsetType
2915     progress;
2916
2917   register ssize_t
2918     i;
2919
2920   ssize_t
2921     y;
2922
2923   /*
2924     Allocate and initialize levels map.
2925   */
2926   assert(image != (Image *) NULL);
2927   assert(image->signature == MagickCoreSignature);
2928   if (image->debug != MagickFalse)
2929     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2930   if (image->storage_class == PseudoClass)
2931     for (i=0; i < (ssize_t) image->colors; i++)
2932     {
2933       /*
2934         Level colormap.
2935       */
2936       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2937         image->colormap[i].red=(double) ClampToQuantum(LevelPixel(black_point,
2938           white_point,gamma,image->colormap[i].red));
2939       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2940         image->colormap[i].green=(double) ClampToQuantum(LevelPixel(black_point,
2941           white_point,gamma,image->colormap[i].green));
2942       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2943         image->colormap[i].blue=(double) ClampToQuantum(LevelPixel(black_point,
2944           white_point,gamma,image->colormap[i].blue));
2945       if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
2946         image->colormap[i].alpha=(double) ClampToQuantum(LevelPixel(black_point,
2947           white_point,gamma,image->colormap[i].alpha));
2948     }
2949   /*
2950     Level image.
2951   */
2952   status=MagickTrue;
2953   progress=0;
2954   image_view=AcquireAuthenticCacheView(image,exception);
2955 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2956   #pragma omp parallel for schedule(static) shared(progress,status) \
2957     magick_number_threads(image,image,image->rows,1)
2958 #endif
2959   for (y=0; y < (ssize_t) image->rows; y++)
2960   {
2961     register Quantum
2962       *magick_restrict q;
2963
2964     register ssize_t
2965       x;
2966
2967     if (status == MagickFalse)
2968       continue;
2969     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
2970     if (q == (Quantum *) NULL)
2971       {
2972         status=MagickFalse;
2973         continue;
2974       }
2975     for (x=0; x < (ssize_t) image->columns; x++)
2976     {
2977       register ssize_t
2978         j;
2979
2980       for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
2981       {
2982         PixelChannel channel = GetPixelChannelChannel(image,j);
2983         PixelTrait traits = GetPixelChannelTraits(image,channel);
2984         if ((traits & UpdatePixelTrait) == 0)
2985           continue;
2986         q[j]=ClampToQuantum(LevelPixel(black_point,white_point,gamma,
2987           (double) q[j]));
2988       }
2989       q+=GetPixelChannels(image);
2990     }
2991     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
2992       status=MagickFalse;
2993     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2994       {
2995         MagickBooleanType
2996           proceed;
2997
2998 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2999         #pragma omp atomic
3000 #endif
3001         progress++;
3002         proceed=SetImageProgress(image,LevelImageTag,progress,image->rows);
3003         if (proceed == MagickFalse)
3004           status=MagickFalse;
3005       }
3006   }
3007   image_view=DestroyCacheView(image_view);
3008   (void) ClampImage(image,exception);
3009   return(status);
3010 }
3011 \f
3012 /*
3013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3014 %                                                                             %
3015 %                                                                             %
3016 %                                                                             %
3017 %     L e v e l i z e I m a g e                                               %
3018 %                                                                             %
3019 %                                                                             %
3020 %                                                                             %
3021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3022 %
3023 %  LevelizeImage() applies the reversed LevelImage() operation to just
3024 %  the specific channels specified.  It compresses the full range of color
3025 %  values, so that they lie between the given black and white points. Gamma is
3026 %  applied before the values are mapped.
3027 %
3028 %  LevelizeImage() can be called with by using a +level command line
3029 %  API option, or using a '!' on a -level or LevelImage() geometry string.
3030 %
3031 %  It can be used to de-contrast a greyscale image to the exact levels
3032 %  specified.  Or by using specific levels for each channel of an image you
3033 %  can convert a gray-scale image to any linear color gradient, according to
3034 %  those levels.
3035 %
3036 %  The format of the LevelizeImage method is:
3037 %
3038 %      MagickBooleanType LevelizeImage(Image *image,const double black_point,
3039 %        const double white_point,const double gamma,ExceptionInfo *exception)
3040 %
3041 %  A description of each parameter follows:
3042 %
3043 %    o image: the image.
3044 %
3045 %    o black_point: The level to map zero (black) to.
3046 %
3047 %    o white_point: The level to map QuantumRange (white) to.
3048 %
3049 %    o gamma: adjust gamma by this factor before mapping values.
3050 %
3051 %    o exception: return any errors or warnings in this structure.
3052 %
3053 */
3054 MagickExport MagickBooleanType LevelizeImage(Image *image,
3055   const double black_point,const double white_point,const double gamma,
3056   ExceptionInfo *exception)
3057 {
3058 #define LevelizeImageTag  "Levelize/Image"
3059 #define LevelizeValue(x) ClampToQuantum(((MagickRealType) gamma_pow((double) \
3060   (QuantumScale*(x)),gamma))*(white_point-black_point)+black_point)
3061
3062   CacheView
3063     *image_view;
3064
3065   MagickBooleanType
3066     status;
3067
3068   MagickOffsetType
3069     progress;
3070
3071   register ssize_t
3072     i;
3073
3074   ssize_t
3075     y;
3076
3077   /*
3078     Allocate and initialize levels map.
3079   */
3080   assert(image != (Image *) NULL);
3081   assert(image->signature == MagickCoreSignature);
3082   if (image->debug != MagickFalse)
3083     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3084   if (image->storage_class == PseudoClass)
3085     for (i=0; i < (ssize_t) image->colors; i++)
3086     {
3087       /*
3088         Level colormap.
3089       */
3090       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3091         image->colormap[i].red=(double) LevelizeValue(image->colormap[i].red);
3092       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3093         image->colormap[i].green=(double) LevelizeValue(
3094           image->colormap[i].green);
3095       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3096         image->colormap[i].blue=(double) LevelizeValue(image->colormap[i].blue);
3097       if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
3098         image->colormap[i].alpha=(double) LevelizeValue(
3099           image->colormap[i].alpha);
3100     }
3101   /*
3102     Level image.
3103   */
3104   status=MagickTrue;
3105   progress=0;
3106   image_view=AcquireAuthenticCacheView(image,exception);
3107 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3108   #pragma omp parallel for schedule(static) shared(progress,status) \
3109     magick_number_threads(image,image,image->rows,1)
3110 #endif
3111   for (y=0; y < (ssize_t) image->rows; y++)
3112   {
3113     register Quantum
3114       *magick_restrict q;
3115
3116     register ssize_t
3117       x;
3118
3119     if (status == MagickFalse)
3120       continue;
3121     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3122     if (q == (Quantum *) NULL)
3123       {
3124         status=MagickFalse;
3125         continue;
3126       }
3127     for (x=0; x < (ssize_t) image->columns; x++)
3128     {
3129       register ssize_t
3130         j;
3131
3132       for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
3133       {
3134         PixelChannel channel = GetPixelChannelChannel(image,j);
3135         PixelTrait traits = GetPixelChannelTraits(image,channel);
3136         if ((traits & UpdatePixelTrait) == 0)
3137           continue;
3138         q[j]=LevelizeValue(q[j]);
3139       }
3140       q+=GetPixelChannels(image);
3141     }
3142     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3143       status=MagickFalse;
3144     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3145       {
3146         MagickBooleanType
3147           proceed;
3148
3149 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3150         #pragma omp atomic
3151 #endif
3152         progress++;
3153         proceed=SetImageProgress(image,LevelizeImageTag,progress,image->rows);
3154         if (proceed == MagickFalse)
3155           status=MagickFalse;
3156       }
3157   }
3158   image_view=DestroyCacheView(image_view);
3159   return(status);
3160 }
3161 \f
3162 /*
3163 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3164 %                                                                             %
3165 %                                                                             %
3166 %                                                                             %
3167 %     L e v e l I m a g e C o l o r s                                         %
3168 %                                                                             %
3169 %                                                                             %
3170 %                                                                             %
3171 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3172 %
3173 %  LevelImageColors() maps the given color to "black" and "white" values,
3174 %  linearly spreading out the colors, and level values on a channel by channel
3175 %  bases, as per LevelImage().  The given colors allows you to specify
3176 %  different level ranges for each of the color channels separately.
3177 %
3178 %  If the boolean 'invert' is set true the image values will modifyed in the
3179 %  reverse direction. That is any existing "black" and "white" colors in the
3180 %  image will become the color values given, with all other values compressed
3181 %  appropriatally.  This effectivally maps a greyscale gradient into the given
3182 %  color gradient.
3183 %
3184 %  The format of the LevelImageColors method is:
3185 %
3186 %    MagickBooleanType LevelImageColors(Image *image,
3187 %      const PixelInfo *black_color,const PixelInfo *white_color,
3188 %      const MagickBooleanType invert,ExceptionInfo *exception)
3189 %
3190 %  A description of each parameter follows:
3191 %
3192 %    o image: the image.
3193 %
3194 %    o black_color: The color to map black to/from
3195 %
3196 %    o white_point: The color to map white to/from
3197 %
3198 %    o invert: if true map the colors (levelize), rather than from (level)
3199 %
3200 %    o exception: return any errors or warnings in this structure.
3201 %
3202 */
3203 MagickExport MagickBooleanType LevelImageColors(Image *image,
3204   const PixelInfo *black_color,const PixelInfo *white_color,
3205   const MagickBooleanType invert,ExceptionInfo *exception)
3206 {
3207   ChannelType
3208     channel_mask;
3209
3210   MagickStatusType
3211     status;
3212
3213   /*
3214     Allocate and initialize levels map.
3215   */
3216   assert(image != (Image *) NULL);
3217   assert(image->signature == MagickCoreSignature);
3218   if (image->debug != MagickFalse)
3219     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3220   if ((IsGrayColorspace(image->colorspace) != MagickFalse) &&
3221       ((IsGrayColorspace(black_color->colorspace) == MagickFalse) ||
3222        (IsGrayColorspace(white_color->colorspace) == MagickFalse)))
3223     (void) SetImageColorspace(image,sRGBColorspace,exception);
3224   status=MagickTrue;
3225   if (invert == MagickFalse)
3226     {
3227       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3228         {
3229           channel_mask=SetImageChannelMask(image,RedChannel);
3230           status&=LevelImage(image,black_color->red,white_color->red,1.0,
3231             exception);
3232           (void) SetImageChannelMask(image,channel_mask);
3233         }
3234       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3235         {
3236           channel_mask=SetImageChannelMask(image,GreenChannel);
3237           status&=LevelImage(image,black_color->green,white_color->green,1.0,
3238             exception);
3239           (void) SetImageChannelMask(image,channel_mask);
3240         }
3241       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3242         {
3243           channel_mask=SetImageChannelMask(image,BlueChannel);
3244           status&=LevelImage(image,black_color->blue,white_color->blue,1.0,
3245             exception);
3246           (void) SetImageChannelMask(image,channel_mask);
3247         }
3248       if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3249           (image->colorspace == CMYKColorspace))
3250         {
3251           channel_mask=SetImageChannelMask(image,BlackChannel);
3252           status&=LevelImage(image,black_color->black,white_color->black,1.0,
3253             exception);
3254           (void) SetImageChannelMask(image,channel_mask);
3255         }
3256       if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3257           (image->alpha_trait != UndefinedPixelTrait))
3258         {
3259           channel_mask=SetImageChannelMask(image,AlphaChannel);
3260           status&=LevelImage(image,black_color->alpha,white_color->alpha,1.0,
3261             exception);
3262           (void) SetImageChannelMask(image,channel_mask);
3263         }
3264     }
3265   else
3266     {
3267       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3268         {
3269           channel_mask=SetImageChannelMask(image,RedChannel);
3270           status&=LevelizeImage(image,black_color->red,white_color->red,1.0,
3271             exception);
3272           (void) SetImageChannelMask(image,channel_mask);
3273         }
3274       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3275         {
3276           channel_mask=SetImageChannelMask(image,GreenChannel);
3277           status&=LevelizeImage(image,black_color->green,white_color->green,1.0,
3278             exception);
3279           (void) SetImageChannelMask(image,channel_mask);
3280         }
3281       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3282         {
3283           channel_mask=SetImageChannelMask(image,BlueChannel);
3284           status&=LevelizeImage(image,black_color->blue,white_color->blue,1.0,
3285             exception);
3286           (void) SetImageChannelMask(image,channel_mask);
3287         }
3288       if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3289           (image->colorspace == CMYKColorspace))
3290         {
3291           channel_mask=SetImageChannelMask(image,BlackChannel);
3292           status&=LevelizeImage(image,black_color->black,white_color->black,1.0,
3293             exception);
3294           (void) SetImageChannelMask(image,channel_mask);
3295         }
3296       if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3297           (image->alpha_trait != UndefinedPixelTrait))
3298         {
3299           channel_mask=SetImageChannelMask(image,AlphaChannel);
3300           status&=LevelizeImage(image,black_color->alpha,white_color->alpha,1.0,
3301             exception);
3302           (void) SetImageChannelMask(image,channel_mask);
3303         }
3304     }
3305   return(status != 0 ? MagickTrue : MagickFalse);
3306 }
3307 \f
3308 /*
3309 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3310 %                                                                             %
3311 %                                                                             %
3312 %                                                                             %
3313 %     L i n e a r S t r e t c h I m a g e                                     %
3314 %                                                                             %
3315 %                                                                             %
3316 %                                                                             %
3317 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3318 %
3319 %  LinearStretchImage() discards any pixels below the black point and above
3320 %  the white point and levels the remaining pixels.
3321 %
3322 %  The format of the LinearStretchImage method is:
3323 %
3324 %      MagickBooleanType LinearStretchImage(Image *image,
3325 %        const double black_point,const double white_point,
3326 %        ExceptionInfo *exception)
3327 %
3328 %  A description of each parameter follows:
3329 %
3330 %    o image: the image.
3331 %
3332 %    o black_point: the black point.
3333 %
3334 %    o white_point: the white point.
3335 %
3336 %    o exception: return any errors or warnings in this structure.
3337 %
3338 */
3339 MagickExport MagickBooleanType LinearStretchImage(Image *image,
3340   const double black_point,const double white_point,ExceptionInfo *exception)
3341 {
3342 #define LinearStretchImageTag  "LinearStretch/Image"
3343
3344   CacheView
3345     *image_view;
3346
3347   double
3348     *histogram,
3349     intensity;
3350
3351   MagickBooleanType
3352     status;
3353
3354   ssize_t
3355     black,
3356     white,
3357     y;
3358
3359   /*
3360     Allocate histogram and linear map.
3361   */
3362   assert(image != (Image *) NULL);
3363   assert(image->signature == MagickCoreSignature);
3364   histogram=(double *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*histogram));
3365   if (histogram == (double *) NULL)
3366     ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
3367       image->filename);
3368   /*
3369     Form histogram.
3370   */
3371   (void) memset(histogram,0,(MaxMap+1)*sizeof(*histogram));
3372   image_view=AcquireVirtualCacheView(image,exception);
3373   for (y=0; y < (ssize_t) image->rows; y++)
3374   {
3375     register const Quantum
3376       *magick_restrict p;
3377
3378     register ssize_t
3379       x;
3380
3381     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3382     if (p == (const Quantum *) NULL)
3383       break;
3384     for (x=0; x < (ssize_t) image->columns; x++)
3385     {
3386       intensity=GetPixelIntensity(image,p);
3387       histogram[ScaleQuantumToMap(ClampToQuantum(intensity))]++;
3388       p+=GetPixelChannels(image);
3389     }
3390   }
3391   image_view=DestroyCacheView(image_view);
3392   /*
3393     Find the histogram boundaries by locating the black and white point levels.
3394   */
3395   intensity=0.0;
3396   for (black=0; black < (ssize_t) MaxMap; black++)
3397   {
3398     intensity+=histogram[black];
3399     if (intensity >= black_point)
3400       break;
3401   }
3402   intensity=0.0;
3403   for (white=(ssize_t) MaxMap; white != 0; white--)
3404   {
3405     intensity+=histogram[white];
3406     if (intensity >= white_point)
3407       break;
3408   }
3409   histogram=(double *) RelinquishMagickMemory(histogram);
3410   status=LevelImage(image,(double) ScaleMapToQuantum((MagickRealType) black),
3411     (double) ScaleMapToQuantum((MagickRealType) white),1.0,exception);
3412   return(status);
3413 }
3414
3415
3416 /*
3417 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3418 %                                                                             %
3419 %                                                                             %
3420 %                                                                             %
3421 %     M o d u l a t e I m a g e                                               %
3422 %                                                                             %
3423 %                                                                             %
3424 %                                                                             %
3425 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3426 %
3427 %  ModulateImage() lets you control the brightness, saturation, and hue
3428 %  of an image.  Modulate represents the brightness, saturation, and hue
3429 %  as one parameter (e.g. 90,150,100).  If the image colorspace is HSL, the
3430 %  modulation is lightness, saturation, and hue.  For HWB, use blackness,
3431 %  whiteness, and hue. And for HCL, use chrome, luma, and hue.
3432 %
3433 %  The format of the ModulateImage method is:
3434 %
3435 %      MagickBooleanType ModulateImage(Image *image,const char *modulate,
3436 %        ExceptionInfo *exception)
3437 %
3438 %  A description of each parameter follows:
3439 %
3440 %    o image: the image.
3441 %
3442 %    o modulate: Define the percent change in brightness, saturation, and hue.
3443 %
3444 %    o exception: return any errors or warnings in this structure.
3445 %
3446 */
3447
3448 static inline void ModulateHCL(const double percent_hue,
3449   const double percent_chroma,const double percent_luma,double *red,
3450   double *green,double *blue)
3451 {
3452   double
3453     hue,
3454     luma,
3455     chroma;
3456
3457   /*
3458     Increase or decrease color luma, chroma, or hue.
3459   */
3460   ConvertRGBToHCL(*red,*green,*blue,&hue,&chroma,&luma);
3461   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3462   chroma*=0.01*percent_chroma;
3463   luma*=0.01*percent_luma;
3464   ConvertHCLToRGB(hue,chroma,luma,red,green,blue);
3465 }
3466
3467 static inline void ModulateHCLp(const double percent_hue,
3468   const double percent_chroma,const double percent_luma,double *red,
3469   double *green,double *blue)
3470 {
3471   double
3472     hue,
3473     luma,
3474     chroma;
3475
3476   /*
3477     Increase or decrease color luma, chroma, or hue.
3478   */
3479   ConvertRGBToHCLp(*red,*green,*blue,&hue,&chroma,&luma);
3480   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3481   chroma*=0.01*percent_chroma;
3482   luma*=0.01*percent_luma;
3483   ConvertHCLpToRGB(hue,chroma,luma,red,green,blue);
3484 }
3485
3486 static inline void ModulateHSB(const double percent_hue,
3487   const double percent_saturation,const double percent_brightness,double *red,
3488   double *green,double *blue)
3489 {
3490   double
3491     brightness,
3492     hue,
3493     saturation;
3494
3495   /*
3496     Increase or decrease color brightness, saturation, or hue.
3497   */
3498   ConvertRGBToHSB(*red,*green,*blue,&hue,&saturation,&brightness);
3499   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3500   saturation*=0.01*percent_saturation;
3501   brightness*=0.01*percent_brightness;
3502   ConvertHSBToRGB(hue,saturation,brightness,red,green,blue);
3503 }
3504
3505 static inline void ModulateHSI(const double percent_hue,
3506   const double percent_saturation,const double percent_intensity,double *red,
3507   double *green,double *blue)
3508 {
3509   double
3510     intensity,
3511     hue,
3512     saturation;
3513
3514   /*
3515     Increase or decrease color intensity, saturation, or hue.
3516   */
3517   ConvertRGBToHSI(*red,*green,*blue,&hue,&saturation,&intensity);
3518   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3519   saturation*=0.01*percent_saturation;
3520   intensity*=0.01*percent_intensity;
3521   ConvertHSIToRGB(hue,saturation,intensity,red,green,blue);
3522 }
3523
3524 static inline void ModulateHSL(const double percent_hue,
3525   const double percent_saturation,const double percent_lightness,double *red,
3526   double *green,double *blue)
3527 {
3528   double
3529     hue,
3530     lightness,
3531     saturation;
3532
3533   /*
3534     Increase or decrease color lightness, saturation, or hue.
3535   */
3536   ConvertRGBToHSL(*red,*green,*blue,&hue,&saturation,&lightness);
3537   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3538   saturation*=0.01*percent_saturation;
3539   lightness*=0.01*percent_lightness;
3540   ConvertHSLToRGB(hue,saturation,lightness,red,green,blue);
3541 }
3542
3543 static inline void ModulateHSV(const double percent_hue,
3544   const double percent_saturation,const double percent_value,double *red,
3545   double *green,double *blue)
3546 {
3547   double
3548     hue,
3549     saturation,
3550     value;
3551
3552   /*
3553     Increase or decrease color value, saturation, or hue.
3554   */
3555   ConvertRGBToHSV(*red,*green,*blue,&hue,&saturation,&value);
3556   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3557   saturation*=0.01*percent_saturation;
3558   value*=0.01*percent_value;
3559   ConvertHSVToRGB(hue,saturation,value,red,green,blue);
3560 }
3561
3562 static inline void ModulateHWB(const double percent_hue,
3563   const double percent_whiteness,const double percent_blackness,double *red,
3564   double *green,double *blue)
3565 {
3566   double
3567     blackness,
3568     hue,
3569     whiteness;
3570
3571   /*
3572     Increase or decrease color blackness, whiteness, or hue.
3573   */
3574   ConvertRGBToHWB(*red,*green,*blue,&hue,&whiteness,&blackness);
3575   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3576   blackness*=0.01*percent_blackness;
3577   whiteness*=0.01*percent_whiteness;
3578   ConvertHWBToRGB(hue,whiteness,blackness,red,green,blue);
3579 }
3580
3581 static inline void ModulateLCHab(const double percent_luma,
3582   const double percent_chroma,const double percent_hue,double *red,
3583   double *green,double *blue)
3584 {
3585   double
3586     hue,
3587     luma,
3588     chroma;
3589
3590   /*
3591     Increase or decrease color luma, chroma, or hue.
3592   */
3593   ConvertRGBToLCHab(*red,*green,*blue,&luma,&chroma,&hue);
3594   luma*=0.01*percent_luma;
3595   chroma*=0.01*percent_chroma;
3596   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3597   ConvertLCHabToRGB(luma,chroma,hue,red,green,blue);
3598 }
3599
3600 static inline void ModulateLCHuv(const double percent_luma,
3601   const double percent_chroma,const double percent_hue,double *red,
3602   double *green,double *blue)
3603 {
3604   double
3605     hue,
3606     luma,
3607     chroma;
3608
3609   /*
3610     Increase or decrease color luma, chroma, or hue.
3611   */
3612   ConvertRGBToLCHuv(*red,*green,*blue,&luma,&chroma,&hue);
3613   luma*=0.01*percent_luma;
3614   chroma*=0.01*percent_chroma;
3615   hue+=fmod((percent_hue-100.0),200.0)/200.0;
3616   ConvertLCHuvToRGB(luma,chroma,hue,red,green,blue);
3617 }
3618
3619 MagickExport MagickBooleanType ModulateImage(Image *image,const char *modulate,
3620   ExceptionInfo *exception)
3621 {
3622 #define ModulateImageTag  "Modulate/Image"
3623
3624   CacheView
3625     *image_view;
3626
3627   ColorspaceType
3628     colorspace;
3629
3630   const char
3631     *artifact;
3632
3633   double
3634     percent_brightness,
3635     percent_hue,
3636     percent_saturation;
3637
3638   GeometryInfo
3639     geometry_info;
3640
3641   MagickBooleanType
3642     status;
3643
3644   MagickOffsetType
3645     progress;
3646
3647   MagickStatusType
3648     flags;
3649
3650   register ssize_t
3651     i;
3652
3653   ssize_t
3654     y;
3655
3656   /*
3657     Initialize modulate table.
3658   */
3659   assert(image != (Image *) NULL);
3660   assert(image->signature == MagickCoreSignature);
3661   if (image->debug != MagickFalse)
3662     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3663   if (modulate == (char *) NULL)
3664     return(MagickFalse);
3665   if (IssRGBCompatibleColorspace(image->colorspace) == MagickFalse)
3666     (void) SetImageColorspace(image,sRGBColorspace,exception);
3667   flags=ParseGeometry(modulate,&geometry_info);
3668   percent_brightness=geometry_info.rho;
3669   percent_saturation=geometry_info.sigma;
3670   if ((flags & SigmaValue) == 0)
3671     percent_saturation=100.0;
3672   percent_hue=geometry_info.xi;
3673   if ((flags & XiValue) == 0)
3674     percent_hue=100.0;
3675   colorspace=UndefinedColorspace;
3676   artifact=GetImageArtifact(image,"modulate:colorspace");
3677   if (artifact != (const char *) NULL)
3678     colorspace=(ColorspaceType) ParseCommandOption(MagickColorspaceOptions,
3679       MagickFalse,artifact);
3680   if (image->storage_class == PseudoClass)
3681     for (i=0; i < (ssize_t) image->colors; i++)
3682     {
3683       double
3684         blue,
3685         green,
3686         red;
3687
3688       /*
3689         Modulate image colormap.
3690       */
3691       red=(double) image->colormap[i].red;
3692       green=(double) image->colormap[i].green;
3693       blue=(double) image->colormap[i].blue;
3694       switch (colorspace)
3695       {
3696         case HCLColorspace:
3697         {
3698           ModulateHCL(percent_hue,percent_saturation,percent_brightness,
3699             &red,&green,&blue);
3700           break;
3701         }
3702         case HCLpColorspace:
3703         {
3704           ModulateHCLp(percent_hue,percent_saturation,percent_brightness,
3705             &red,&green,&blue);
3706           break;
3707         }
3708         case HSBColorspace:
3709         {
3710           ModulateHSB(percent_hue,percent_saturation,percent_brightness,
3711             &red,&green,&blue);
3712           break;
3713         }
3714         case HSIColorspace:
3715         {
3716           ModulateHSI(percent_hue,percent_saturation,percent_brightness,
3717             &red,&green,&blue);
3718           break;
3719         }
3720         case HSLColorspace:
3721         default:
3722         {
3723           ModulateHSL(percent_hue,percent_saturation,percent_brightness,
3724             &red,&green,&blue);
3725           break;
3726         }
3727         case HSVColorspace:
3728         {
3729           ModulateHSV(percent_hue,percent_saturation,percent_brightness,
3730             &red,&green,&blue);
3731           break;
3732         }
3733         case HWBColorspace:
3734         {
3735           ModulateHWB(percent_hue,percent_saturation,percent_brightness,
3736             &red,&green,&blue);
3737           break;
3738         }
3739         case LCHColorspace:
3740         case LCHabColorspace:
3741         {
3742           ModulateLCHab(percent_brightness,percent_saturation,percent_hue,
3743             &red,&green,&blue);
3744           break;
3745         }
3746         case LCHuvColorspace:
3747         {
3748           ModulateLCHuv(percent_brightness,percent_saturation,percent_hue,
3749             &red,&green,&blue);
3750           break;
3751         }
3752       }
3753       image->colormap[i].red=red;
3754       image->colormap[i].green=green;
3755       image->colormap[i].blue=blue;
3756     }
3757   /*
3758     Modulate image.
3759   */
3760 #if defined(MAGICKCORE_OPENCL_SUPPORT)
3761   if (AccelerateModulateImage(image,percent_brightness,percent_hue,
3762         percent_saturation,colorspace,exception) != MagickFalse)
3763     return(MagickTrue);
3764 #endif
3765   status=MagickTrue;
3766   progress=0;
3767   image_view=AcquireAuthenticCacheView(image,exception);
3768 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3769   #pragma omp parallel for schedule(static) shared(progress,status) \
3770     magick_number_threads(image,image,image->rows,1)
3771 #endif
3772   for (y=0; y < (ssize_t) image->rows; y++)
3773   {
3774     register Quantum
3775       *magick_restrict q;
3776
3777     register ssize_t
3778       x;
3779
3780     if (status == MagickFalse)
3781       continue;
3782     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
3783     if (q == (Quantum *) NULL)
3784       {
3785         status=MagickFalse;
3786         continue;
3787       }
3788     for (x=0; x < (ssize_t) image->columns; x++)
3789     {
3790       double
3791         blue,
3792         green,
3793         red;
3794
3795       red=(double) GetPixelRed(image,q);
3796       green=(double) GetPixelGreen(image,q);
3797       blue=(double) GetPixelBlue(image,q);
3798       switch (colorspace)
3799       {
3800         case HCLColorspace:
3801         {
3802           ModulateHCL(percent_hue,percent_saturation,percent_brightness,
3803             &red,&green,&blue);
3804           break;
3805         }
3806         case HCLpColorspace:
3807         {
3808           ModulateHCLp(percent_hue,percent_saturation,percent_brightness,
3809             &red,&green,&blue);
3810           break;
3811         }
3812         case HSBColorspace:
3813         {
3814           ModulateHSB(percent_hue,percent_saturation,percent_brightness,
3815             &red,&green,&blue);
3816           break;
3817         }
3818         case HSLColorspace:
3819         default:
3820         {
3821           ModulateHSL(percent_hue,percent_saturation,percent_brightness,
3822             &red,&green,&blue);
3823           break;
3824         }
3825         case HSVColorspace:
3826         {
3827           ModulateHSV(percent_hue,percent_saturation,percent_brightness,
3828             &red,&green,&blue);
3829           break;
3830         }
3831         case HWBColorspace:
3832         {
3833           ModulateHWB(percent_hue,percent_saturation,percent_brightness,
3834             &red,&green,&blue);
3835           break;
3836         }
3837         case LCHabColorspace:
3838         {
3839           ModulateLCHab(percent_brightness,percent_saturation,percent_hue,
3840             &red,&green,&blue);
3841           break;
3842         }
3843         case LCHColorspace:
3844         case LCHuvColorspace:
3845         {
3846           ModulateLCHuv(percent_brightness,percent_saturation,percent_hue,
3847             &red,&green,&blue);
3848           break;
3849         }
3850       }
3851       SetPixelRed(image,ClampToQuantum(red),q);
3852       SetPixelGreen(image,ClampToQuantum(green),q);
3853       SetPixelBlue(image,ClampToQuantum(blue),q);
3854       q+=GetPixelChannels(image);
3855     }
3856     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
3857       status=MagickFalse;
3858     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3859       {
3860         MagickBooleanType
3861           proceed;
3862
3863 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3864         #pragma omp atomic
3865 #endif
3866         progress++;
3867         proceed=SetImageProgress(image,ModulateImageTag,progress,image->rows);
3868         if (proceed == MagickFalse)
3869           status=MagickFalse;
3870       }
3871   }
3872   image_view=DestroyCacheView(image_view);
3873   return(status);
3874 }
3875 \f
3876 /*
3877 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3878 %                                                                             %
3879 %                                                                             %
3880 %                                                                             %
3881 %     N e g a t e I m a g e                                                   %
3882 %                                                                             %
3883 %                                                                             %
3884 %                                                                             %
3885 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3886 %
3887 %  NegateImage() negates the colors in the reference image.  The grayscale
3888 %  option means that only grayscale values within the image are negated.
3889 %
3890 %  The format of the NegateImage method is:
3891 %
3892 %      MagickBooleanType NegateImage(Image *image,
3893 %        const MagickBooleanType grayscale,ExceptionInfo *exception)
3894 %
3895 %  A description of each parameter follows:
3896 %
3897 %    o image: the image.
3898 %
3899 %    o grayscale: If MagickTrue, only negate grayscale pixels within the image.
3900 %
3901 %    o exception: return any errors or warnings in this structure.
3902 %
3903 */
3904 MagickExport MagickBooleanType NegateImage(Image *image,
3905   const MagickBooleanType grayscale,ExceptionInfo *exception)
3906 {
3907 #define NegateImageTag  "Negate/Image"
3908
3909   CacheView
3910     *image_view;
3911
3912   MagickBooleanType
3913     status;
3914
3915   MagickOffsetType
3916     progress;
3917
3918   register ssize_t
3919     i;
3920
3921   ssize_t
3922     y;
3923
3924   assert(image != (Image *) NULL);
3925   assert(image->signature == MagickCoreSignature);
3926   if (image->debug != MagickFalse)
3927     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3928   if (image->storage_class == PseudoClass)
3929     for (i=0; i < (ssize_t) image->colors; i++)
3930     {
3931       /*
3932         Negate colormap.
3933       */
3934       if( grayscale != MagickFalse )
3935         if ((image->colormap[i].red != image->colormap[i].green) ||
3936             (image->colormap[i].green != image->colormap[i].blue))
3937           continue;
3938       if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3939         image->colormap[i].red=QuantumRange-image->colormap[i].red;
3940       if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3941         image->colormap[i].green=QuantumRange-image->colormap[i].green;
3942       if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3943         image->colormap[i].blue=QuantumRange-image->colormap[i].blue;
3944     }
3945   /*
3946     Negate image.
3947   */
3948   status=MagickTrue;
3949   progress=0;
3950   image_view=AcquireAuthenticCacheView(image,exception);
3951   if( grayscale != MagickFalse )
3952     {
3953       for (y=0; y < (ssize_t) image->rows; y++)
3954       {
3955         MagickBooleanType
3956           sync;
3957
3958         register Quantum
3959           *magick_restrict q;
3960
3961         register ssize_t
3962           x;
3963
3964         if (status == MagickFalse)
3965           continue;
3966         q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,
3967           exception);
3968         if (q == (Quantum *) NULL)
3969           {
3970             status=MagickFalse;
3971             continue;
3972           }
3973         for (x=0; x < (ssize_t) image->columns; x++)
3974         {
3975           register ssize_t
3976             j;
3977
3978           if (IsPixelGray(image,q) != MagickFalse)
3979             {
3980               q+=GetPixelChannels(image);
3981               continue;
3982             }
3983           for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
3984           {
3985             PixelChannel channel = GetPixelChannelChannel(image,j);
3986             PixelTrait traits = GetPixelChannelTraits(image,channel);
3987             if ((traits & UpdatePixelTrait) == 0)
3988               continue;
3989             q[j]=QuantumRange-q[j];
3990           }
3991           q+=GetPixelChannels(image);
3992         }
3993         sync=SyncCacheViewAuthenticPixels(image_view,exception);
3994         if (sync == MagickFalse)
3995           status=MagickFalse;
3996         if (image->progress_monitor != (MagickProgressMonitor) NULL)
3997           {
3998             MagickBooleanType
3999               proceed;
4000
4001 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4002             #pragma omp atomic
4003 #endif
4004             progress++;
4005             proceed=SetImageProgress(image,NegateImageTag,progress,image->rows);
4006             if (proceed == MagickFalse)
4007               status=MagickFalse;
4008           }
4009       }
4010       image_view=DestroyCacheView(image_view);
4011       return(MagickTrue);
4012     }
4013   /*
4014     Negate image.
4015   */
4016 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4017   #pragma omp parallel for schedule(static) shared(progress,status) \
4018     magick_number_threads(image,image,image->rows,1)
4019 #endif
4020   for (y=0; y < (ssize_t) image->rows; y++)
4021   {
4022     register Quantum
4023       *magick_restrict q;
4024
4025     register ssize_t
4026       x;
4027
4028     if (status == MagickFalse)
4029       continue;
4030     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
4031     if (q == (Quantum *) NULL)
4032       {
4033         status=MagickFalse;
4034         continue;
4035       }
4036     for (x=0; x < (ssize_t) image->columns; x++)
4037     {
4038       register ssize_t
4039         j;
4040
4041       for (j=0; j < (ssize_t) GetPixelChannels(image); j++)
4042       {
4043         PixelChannel channel = GetPixelChannelChannel(image,j);
4044         PixelTrait traits = GetPixelChannelTraits(image,channel);
4045         if ((traits & UpdatePixelTrait) == 0)
4046           continue;
4047         q[j]=QuantumRange-q[j];
4048       }
4049       q+=GetPixelChannels(image);
4050     }
4051     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
4052       status=MagickFalse;
4053     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4054       {
4055         MagickBooleanType
4056           proceed;
4057
4058 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4059         #pragma omp atomic
4060 #endif
4061         progress++;
4062         proceed=SetImageProgress(image,NegateImageTag,progress,image->rows);
4063         if (proceed == MagickFalse)
4064           status=MagickFalse;
4065       }
4066   }
4067   image_view=DestroyCacheView(image_view);
4068   return(status);
4069 }
4070 \f
4071 /*
4072 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4073 %                                                                             %
4074 %                                                                             %
4075 %                                                                             %
4076 %     N o r m a l i z e I m a g e                                             %
4077 %                                                                             %
4078 %                                                                             %
4079 %                                                                             %
4080 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4081 %
4082 %  The NormalizeImage() method enhances the contrast of a color image by
4083 %  mapping the darkest 2 percent of all pixel to black and the brightest
4084 %  1 percent to white.
4085 %
4086 %  The format of the NormalizeImage method is:
4087 %
4088 %      MagickBooleanType NormalizeImage(Image *image,ExceptionInfo *exception)
4089 %
4090 %  A description of each parameter follows:
4091 %
4092 %    o image: the image.
4093 %
4094 %    o exception: return any errors or warnings in this structure.
4095 %
4096 */
4097 MagickExport MagickBooleanType NormalizeImage(Image *image,
4098   ExceptionInfo *exception)
4099 {
4100   double
4101     black_point,
4102     white_point;
4103
4104   black_point=(double) image->columns*image->rows*0.0015;
4105   white_point=(double) image->columns*image->rows*0.9995;
4106   return(ContrastStretchImage(image,black_point,white_point,exception));
4107 }
4108 \f
4109 /*
4110 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4111 %                                                                             %
4112 %                                                                             %
4113 %                                                                             %
4114 %     S i g m o i d a l C o n t r a s t I m a g e                             %
4115 %                                                                             %
4116 %                                                                             %
4117 %                                                                             %
4118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4119 %
4120 %  SigmoidalContrastImage() adjusts the contrast of an image with a non-linear
4121 %  sigmoidal contrast algorithm.  Increase the contrast of the image using a
4122 %  sigmoidal transfer function without saturating highlights or shadows.
4123 %  Contrast indicates how much to increase the contrast (0 is none; 3 is
4124 %  typical; 20 is pushing it); mid-point indicates where midtones fall in the
4125 %  resultant image (0 is white; 50% is middle-gray; 100% is black).  Set
4126 %  sharpen to MagickTrue to increase the image contrast otherwise the contrast
4127 %  is reduced.
4128 %
4129 %  The format of the SigmoidalContrastImage method is:
4130 %
4131 %      MagickBooleanType SigmoidalContrastImage(Image *image,
4132 %        const MagickBooleanType sharpen,const char *levels,
4133 %        ExceptionInfo *exception)
4134 %
4135 %  A description of each parameter follows:
4136 %
4137 %    o image: the image.
4138 %
4139 %    o sharpen: Increase or decrease image contrast.
4140 %
4141 %    o contrast: strength of the contrast, the larger the number the more
4142 %      'threshold-like' it becomes.
4143 %
4144 %    o midpoint: midpoint of the function as a color value 0 to QuantumRange.
4145 %
4146 %    o exception: return any errors or warnings in this structure.
4147 %
4148 */
4149
4150 /*
4151   ImageMagick 6 has a version of this function which uses LUTs.
4152 */
4153
4154 /*
4155   Sigmoidal function Sigmoidal with inflexion point moved to b and "slope
4156   constant" set to a.
4157
4158   The first version, based on the hyperbolic tangent tanh, when combined with
4159   the scaling step, is an exact arithmetic clone of the the sigmoid function
4160   based on the logistic curve. The equivalence is based on the identity
4161
4162     1/(1+exp(-t)) = (1+tanh(t/2))/2
4163
4164   (http://de.wikipedia.org/wiki/Sigmoidfunktion) and the fact that the
4165   scaled sigmoidal derivation is invariant under affine transformations of
4166   the ordinate.
4167
4168   The tanh version is almost certainly more accurate and cheaper.  The 0.5
4169   factor in the argument is to clone the legacy ImageMagick behavior. The
4170   reason for making the define depend on atanh even though it only uses tanh
4171   has to do with the construction of the inverse of the scaled sigmoidal.
4172 */
4173 #if defined(MAGICKCORE_HAVE_ATANH)
4174 #define Sigmoidal(a,b,x) ( tanh((0.5*(a))*((x)-(b))) )
4175 #else
4176 #define Sigmoidal(a,b,x) ( 1.0/(1.0+exp((a)*((b)-(x)))) )
4177 #endif
4178 /*
4179   Scaled sigmoidal function:
4180
4181     ( Sigmoidal(a,b,x) - Sigmoidal(a,b,0) ) /
4182     ( Sigmoidal(a,b,1) - Sigmoidal(a,b,0) )
4183
4184   See http://osdir.com/ml/video.image-magick.devel/2005-04/msg00006.html and
4185   http://www.cs.dartmouth.edu/farid/downloads/tutorials/fip.pdf.  The limit
4186   of ScaledSigmoidal as a->0 is the identity, but a=0 gives a division by
4187   zero. This is fixed below by exiting immediately when contrast is small,
4188   leaving the image (or colormap) unmodified. This appears to be safe because
4189   the series expansion of the logistic sigmoidal function around x=b is
4190
4191   1/2-a*(b-x)/4+...
4192
4193   so that the key denominator s(1)-s(0) is about a/4 (a/2 with tanh).
4194 */
4195 #define ScaledSigmoidal(a,b,x) (                    \
4196   (Sigmoidal((a),(b),(x))-Sigmoidal((a),(b),0.0)) / \
4197   (Sigmoidal((a),(b),1.0)-Sigmoidal((a),(b),0.0)) )
4198 /*
4199   Inverse of ScaledSigmoidal, used for +sigmoidal-contrast.  Because b
4200   may be 0 or 1, the argument of the hyperbolic tangent (resp. logistic
4201   sigmoidal) may be outside of the interval (-1,1) (resp. (0,1)), even
4202   when creating a LUT from in gamut values, hence the branching.  In
4203   addition, HDRI may have out of gamut values.
4204   InverseScaledSigmoidal is not a two-sided inverse of ScaledSigmoidal:
4205   It is only a right inverse. This is unavoidable.
4206 */
4207 static inline double InverseScaledSigmoidal(const double a,const double b,
4208   const double x)
4209 {
4210   const double sig0=Sigmoidal(a,b,0.0);
4211   const double sig1=Sigmoidal(a,b,1.0);
4212   const double argument=(sig1-sig0)*x+sig0;
4213   const double clamped=
4214     (
4215 #if defined(MAGICKCORE_HAVE_ATANH)
4216       argument < -1+MagickEpsilon
4217       ?
4218       -1+MagickEpsilon
4219       :
4220       ( argument > 1-MagickEpsilon ? 1-MagickEpsilon : argument )
4221     );
4222   return(b+(2.0/a)*atanh(clamped));
4223 #else
4224       argument < MagickEpsilon
4225       ?
4226       MagickEpsilon
4227       :
4228       ( argument > 1-MagickEpsilon ? 1-MagickEpsilon : argument )
4229     );
4230   return(b-log(1.0/clamped-1.0)/a);
4231 #endif
4232 }
4233
4234 MagickExport MagickBooleanType SigmoidalContrastImage(Image *image,
4235   const MagickBooleanType sharpen,const double contrast,const double midpoint,
4236   ExceptionInfo *exception)
4237 {
4238 #define SigmoidalContrastImageTag  "SigmoidalContrast/Image"
4239 #define ScaledSig(x) ( ClampToQuantum(QuantumRange* \
4240   ScaledSigmoidal(contrast,QuantumScale*midpoint,QuantumScale*(x))) )
4241 #define InverseScaledSig(x) ( ClampToQuantum(QuantumRange* \
4242   InverseScaledSigmoidal(contrast,QuantumScale*midpoint,QuantumScale*(x))) )
4243
4244   CacheView
4245     *image_view;
4246
4247   MagickBooleanType
4248     status;
4249
4250   MagickOffsetType
4251     progress;
4252
4253   ssize_t
4254     y;
4255
4256   /*
4257     Convenience macros.
4258   */
4259   assert(image != (Image *) NULL);
4260   assert(image->signature == MagickCoreSignature);
4261   if (image->debug != MagickFalse)
4262     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4263   /*
4264     Side effect: may clamp values unless contrast<MagickEpsilon, in which
4265     case nothing is done.
4266   */
4267   if (contrast < MagickEpsilon)
4268     return(MagickTrue);
4269   /*
4270     Sigmoidal-contrast enhance colormap.
4271   */
4272   if (image->storage_class == PseudoClass)
4273     {
4274       register ssize_t
4275         i;
4276
4277       if( sharpen != MagickFalse )
4278         for (i=0; i < (ssize_t) image->colors; i++)
4279         {
4280           if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
4281             image->colormap[i].red=(MagickRealType) ScaledSig(
4282               image->colormap[i].red);
4283           if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
4284             image->colormap[i].green=(MagickRealType) ScaledSig(
4285               image->colormap[i].green);
4286           if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
4287             image->colormap[i].blue=(MagickRealType) ScaledSig(
4288               image->colormap[i].blue);
4289           if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
4290             image->colormap[i].alpha=(MagickRealType) ScaledSig(
4291               image->colormap[i].alpha);
4292         }
4293       else
4294         for (i=0; i < (ssize_t) image->colors; i++)
4295         {
4296           if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
4297             image->colormap[i].red=(MagickRealType) InverseScaledSig(
4298               image->colormap[i].red);
4299           if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
4300             image->colormap[i].green=(MagickRealType) InverseScaledSig(
4301               image->colormap[i].green);
4302           if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
4303             image->colormap[i].blue=(MagickRealType) InverseScaledSig(
4304               image->colormap[i].blue);
4305           if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
4306             image->colormap[i].alpha=(MagickRealType) InverseScaledSig(
4307               image->colormap[i].alpha);
4308         }
4309     }
4310   /*
4311     Sigmoidal-contrast enhance image.
4312   */
4313   status=MagickTrue;
4314   progress=0;
4315   image_view=AcquireAuthenticCacheView(image,exception);
4316 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4317   #pragma omp parallel for schedule(static) shared(progress,status) \
4318     magick_number_threads(image,image,image->rows,1)
4319 #endif
4320   for (y=0; y < (ssize_t) image->rows; y++)
4321   {
4322     register Quantum
4323       *magick_restrict q;
4324
4325     register ssize_t
4326       x;
4327
4328     if (status == MagickFalse)
4329       continue;
4330     q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
4331     if (q == (Quantum *) NULL)
4332       {
4333         status=MagickFalse;
4334         continue;
4335       }
4336     for (x=0; x < (ssize_t) image->columns; x++)
4337     {
4338       register ssize_t
4339         i;
4340
4341       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
4342       {
4343         PixelChannel channel = GetPixelChannelChannel(image,i);
4344         PixelTrait traits = GetPixelChannelTraits(image,channel);
4345         if ((traits & UpdatePixelTrait) == 0)
4346           continue;
4347         if( sharpen != MagickFalse )
4348           q[i]=ScaledSig(q[i]);
4349         else
4350           q[i]=InverseScaledSig(q[i]);
4351       }
4352       q+=GetPixelChannels(image);
4353     }
4354     if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
4355       status=MagickFalse;
4356     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4357       {
4358         MagickBooleanType
4359           proceed;
4360
4361 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4362         #pragma omp atomic
4363 #endif
4364         progress++;
4365         proceed=SetImageProgress(image,SigmoidalContrastImageTag,progress,
4366           image->rows);
4367         if (proceed == MagickFalse)
4368           status=MagickFalse;
4369       }
4370   }
4371   image_view=DestroyCacheView(image_view);
4372   return(status);
4373 }