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