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