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