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