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