1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %                                                                             %
4 %                                                                             %
5 %                                                                             %
6 %                   EEEEE  FFFFF  FFFFF  EEEEE  CCCC  TTTTT                   %
7 %                   E      F      F      E     C        T                     %
8 %                   EEE    FFF    FFF    EEE   C        T                     %
9 %                   E      F      F      E     C        T                     %
10 %                   EEEEE  F      F      EEEEE  CCCC    T                     %
11 %                                                                             %
12 %                                                                             %
13 %                       MagickCore Image Effects Methods                      %
14 %                                                                             %
15 %                               Software Design                               %
16 %                                    Cristy                                   %
17 %                                 October 1996                                %
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/blob.h"
46 #include "MagickCore/cache-view.h"
47 #include "MagickCore/color.h"
48 #include "MagickCore/color-private.h"
49 #include "MagickCore/colorspace.h"
50 #include "MagickCore/constitute.h"
51 #include "MagickCore/decorate.h"
52 #include "MagickCore/distort.h"
53 #include "MagickCore/draw.h"
54 #include "MagickCore/enhance.h"
55 #include "MagickCore/exception.h"
56 #include "MagickCore/exception-private.h"
57 #include "MagickCore/effect.h"
58 #include "MagickCore/fx.h"
59 #include "MagickCore/gem.h"
60 #include "MagickCore/gem-private.h"
61 #include "MagickCore/geometry.h"
62 #include "MagickCore/image-private.h"
63 #include "MagickCore/list.h"
64 #include "MagickCore/log.h"
65 #include "MagickCore/matrix.h"
66 #include "MagickCore/memory_.h"
67 #include "MagickCore/memory-private.h"
68 #include "MagickCore/monitor.h"
69 #include "MagickCore/monitor-private.h"
70 #include "MagickCore/montage.h"
71 #include "MagickCore/morphology.h"
72 #include "MagickCore/morphology-private.h"
73 #include "MagickCore/paint.h"
74 #include "MagickCore/pixel-accessor.h"
75 #include "MagickCore/property.h"
76 #include "MagickCore/quantize.h"
77 #include "MagickCore/quantum.h"
78 #include "MagickCore/quantum-private.h"
79 #include "MagickCore/random_.h"
80 #include "MagickCore/random-private.h"
81 #include "MagickCore/resample.h"
82 #include "MagickCore/resample-private.h"
83 #include "MagickCore/resize.h"
84 #include "MagickCore/resource_.h"
85 #include "MagickCore/segment.h"
86 #include "MagickCore/shear.h"
87 #include "MagickCore/signature-private.h"
88 #include "MagickCore/statistic.h"
89 #include "MagickCore/string_.h"
90 #include "MagickCore/thread-private.h"
91 #include "MagickCore/transform.h"
92 #include "MagickCore/threshold.h"
93 
94 /*
95 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96 %                                                                             %
97 %                                                                             %
98 %                                                                             %
99 %     A d a p t i v e B l u r I m a g e                                       %
100 %                                                                             %
101 %                                                                             %
102 %                                                                             %
103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104 %
105 %  AdaptiveBlurImage() adaptively blurs the image by blurring less
106 %  intensely near image edges and more intensely far from edges.  We blur the
107 %  image with a Gaussian operator of the given radius and standard deviation
108 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
109 %  radius of 0 and AdaptiveBlurImage() selects a suitable radius for you.
110 %
111 %  The format of the AdaptiveBlurImage method is:
112 %
113 %      Image *AdaptiveBlurImage(const Image *image,const double radius,
114 %        const double sigma,ExceptionInfo *exception)
115 %
116 %  A description of each parameter follows:
117 %
118 %    o image: the image.
119 %
120 %    o radius: the radius of the Gaussian, in pixels, not counting the center
121 %      pixel.
122 %
123 %    o sigma: the standard deviation of the Laplacian, in pixels.
124 %
125 %    o exception: return any errors or warnings in this structure.
126 %
127 */
AdaptiveBlurImage(const Image * image,const double radius,const double sigma,ExceptionInfo * exception)128 MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
129   const double sigma,ExceptionInfo *exception)
130 {
131 #define AdaptiveBlurImageTag  "Convolve/Image"
132 #define MagickSigma  (fabs(sigma) < MagickEpsilon ? MagickEpsilon : sigma)
133 
134   CacheView
135     *blur_view,
136     *edge_view,
137     *image_view;
138 
139   double
140     normalize,
141     **kernel;
142 
143   Image
144     *blur_image,
145     *edge_image,
146     *gaussian_image;
147 
148   MagickBooleanType
149     status;
150 
151   MagickOffsetType
152     progress;
153 
154   size_t
155     width;
156 
157   ssize_t
158     w,
159     y;
160 
161   assert(image != (const Image *) NULL);
162   assert(image->signature == MagickCoreSignature);
163   if (image->debug != MagickFalse)
164     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
165   assert(exception != (ExceptionInfo *) NULL);
166   assert(exception->signature == MagickCoreSignature);
167   blur_image=CloneImage(image,0,0,MagickTrue,exception);
168   if (blur_image == (Image *) NULL)
169     return((Image *) NULL);
170   if (fabs(sigma) < MagickEpsilon)
171     return(blur_image);
172   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
173     {
174       blur_image=DestroyImage(blur_image);
175       return((Image *) NULL);
176     }
177   /*
178     Edge detect the image brightness channel, level, blur, and level again.
179   */
180   edge_image=EdgeImage(image,radius,exception);
181   if (edge_image == (Image *) NULL)
182     {
183       blur_image=DestroyImage(blur_image);
184       return((Image *) NULL);
185     }
186   (void) AutoLevelImage(edge_image,exception);
187   gaussian_image=BlurImage(edge_image,radius,sigma,exception);
188   if (gaussian_image != (Image *) NULL)
189     {
190       edge_image=DestroyImage(edge_image);
191       edge_image=gaussian_image;
192     }
193   (void) AutoLevelImage(edge_image,exception);
194   /*
195     Create a set of kernels from maximum (radius,sigma) to minimum.
196   */
197   width=GetOptimalKernelWidth2D(radius,sigma);
198   kernel=(double **) MagickAssumeAligned(AcquireAlignedMemory((size_t) width,
199     sizeof(*kernel)));
200   if (kernel == (double **) NULL)
201     {
202       edge_image=DestroyImage(edge_image);
203       blur_image=DestroyImage(blur_image);
204       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
205     }
206   (void) memset(kernel,0,(size_t) width*sizeof(*kernel));
207   for (w=0; w < (ssize_t) width; w+=2)
208   {
209     ssize_t
210       j,
211       k,
212       u,
213       v;
214 
215     kernel[w]=(double *) MagickAssumeAligned(AcquireAlignedMemory(
216       (size_t) (width-w),(width-w)*sizeof(**kernel)));
217     if (kernel[w] == (double *) NULL)
218       break;
219     normalize=0.0;
220     j=(ssize_t) (width-w-1)/2;
221     k=0;
222     for (v=(-j); v <= j; v++)
223     {
224       for (u=(-j); u <= j; u++)
225       {
226         kernel[w][k]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
227           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
228         normalize+=kernel[w][k];
229         k++;
230       }
231     }
232     kernel[w][(k-1)/2]+=(double) (1.0-normalize);
233     if (sigma < MagickEpsilon)
234       kernel[w][(k-1)/2]=1.0;
235   }
236   if (w < (ssize_t) width)
237     {
238       for (w-=2; w >= 0; w-=2)
239         kernel[w]=(double *) RelinquishAlignedMemory(kernel[w]);
240       kernel=(double **) RelinquishAlignedMemory(kernel);
241       edge_image=DestroyImage(edge_image);
242       blur_image=DestroyImage(blur_image);
243       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
244     }
245   /*
246     Adaptively blur image.
247   */
248   status=MagickTrue;
249   progress=0;
250   image_view=AcquireVirtualCacheView(image,exception);
251   edge_view=AcquireVirtualCacheView(edge_image,exception);
252   blur_view=AcquireAuthenticCacheView(blur_image,exception);
253 #if defined(MAGICKCORE_OPENMP_SUPPORT)
254   #pragma omp parallel for schedule(static) shared(progress,status) \
255     magick_number_threads(image,blur_image,blur_image->rows,1)
256 #endif
257   for (y=0; y < (ssize_t) blur_image->rows; y++)
258   {
259     const Quantum
260       *magick_restrict r;
261 
262     Quantum
263       *magick_restrict q;
264 
265     ssize_t
266       x;
267 
268     if (status == MagickFalse)
269       continue;
270     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
271     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
272       exception);
273     if ((r == (const Quantum *) NULL) || (q == (Quantum *) NULL))
274       {
275         status=MagickFalse;
276         continue;
277       }
278     for (x=0; x < (ssize_t) blur_image->columns; x++)
279     {
280       const Quantum
281         *magick_restrict p;
282 
283       ssize_t
284         i;
285 
286       ssize_t
287         center,
288         j;
289 
290       j=CastDoubleToLong(ceil((double) width*(1.0-QuantumScale*
291         GetPixelIntensity(edge_image,r))-0.5));
292       if (j < 0)
293         j=0;
294       else
295         if (j > (ssize_t) width)
296           j=(ssize_t) width;
297       if ((j & 0x01) != 0)
298         j--;
299       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-j)/2L),y-
300         (ssize_t) ((width-j)/2L),width-j,width-j,exception);
301       if (p == (const Quantum *) NULL)
302         break;
303       center=(ssize_t) GetPixelChannels(image)*(width-j)*((width-j)/2L)+
304         GetPixelChannels(image)*((width-j)/2);
305       for (i=0; i < (ssize_t) GetPixelChannels(blur_image); i++)
306       {
307         double
308           alpha,
309           gamma,
310           pixel;
311 
312         PixelChannel
313           channel;
314 
315         PixelTrait
316           blur_traits,
317           traits;
318 
319         const double
320           *magick_restrict k;
321 
322         const Quantum
323           *magick_restrict pixels;
324 
325         ssize_t
326           u;
327 
328         ssize_t
329           v;
330 
331         channel=GetPixelChannelChannel(image,i);
332         traits=GetPixelChannelTraits(image,channel);
333         blur_traits=GetPixelChannelTraits(blur_image,channel);
334         if ((traits == UndefinedPixelTrait) ||
335             (blur_traits == UndefinedPixelTrait))
336           continue;
337         if ((blur_traits & CopyPixelTrait) != 0)
338           {
339             SetPixelChannel(blur_image,channel,p[center+i],q);
340             continue;
341           }
342         k=kernel[j];
343         pixels=p;
344         pixel=0.0;
345         gamma=0.0;
346         if ((blur_traits & BlendPixelTrait) == 0)
347           {
348             /*
349               No alpha blending.
350             */
351             for (v=0; v < (ssize_t) (width-j); v++)
352             {
353               for (u=0; u < (ssize_t) (width-j); u++)
354               {
355                 pixel+=(*k)*pixels[i];
356                 gamma+=(*k);
357                 k++;
358                 pixels+=GetPixelChannels(image);
359               }
360             }
361             gamma=PerceptibleReciprocal(gamma);
362             SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
363             continue;
364           }
365         /*
366           Alpha blending.
367         */
368         for (v=0; v < (ssize_t) (width-j); v++)
369         {
370           for (u=0; u < (ssize_t) (width-j); u++)
371           {
372             alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
373             pixel+=(*k)*alpha*pixels[i];
374             gamma+=(*k)*alpha;
375             k++;
376             pixels+=GetPixelChannels(image);
377           }
378         }
379         gamma=PerceptibleReciprocal(gamma);
380         SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
381       }
382       q+=GetPixelChannels(blur_image);
383       r+=GetPixelChannels(edge_image);
384     }
385     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
386       status=MagickFalse;
387     if (image->progress_monitor != (MagickProgressMonitor) NULL)
388       {
389         MagickBooleanType
390           proceed;
391 
392 #if defined(MAGICKCORE_OPENMP_SUPPORT)
393         #pragma omp atomic
394 #endif
395         progress++;
396         proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress,
397           image->rows);
398         if (proceed == MagickFalse)
399           status=MagickFalse;
400       }
401   }
402   blur_image->type=image->type;
403   blur_view=DestroyCacheView(blur_view);
404   edge_view=DestroyCacheView(edge_view);
405   image_view=DestroyCacheView(image_view);
406   edge_image=DestroyImage(edge_image);
407   for (w=0; w < (ssize_t) width; w+=2)
408     kernel[w]=(double *) RelinquishAlignedMemory(kernel[w]);
409   kernel=(double **) RelinquishAlignedMemory(kernel);
410   if (status == MagickFalse)
411     blur_image=DestroyImage(blur_image);
412   return(blur_image);
413 }
414 
415 /*
416 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
417 %                                                                             %
418 %                                                                             %
419 %                                                                             %
420 %     A d a p t i v e S h a r p e n I m a g e                                 %
421 %                                                                             %
422 %                                                                             %
423 %                                                                             %
424 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
425 %
426 %  AdaptiveSharpenImage() adaptively sharpens the image by sharpening more
427 %  intensely near image edges and less intensely far from edges. We sharpen the
428 %  image with a Gaussian operator of the given radius and standard deviation
429 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
430 %  radius of 0 and AdaptiveSharpenImage() selects a suitable radius for you.
431 %
432 %  The format of the AdaptiveSharpenImage method is:
433 %
434 %      Image *AdaptiveSharpenImage(const Image *image,const double radius,
435 %        const double sigma,ExceptionInfo *exception)
436 %
437 %  A description of each parameter follows:
438 %
439 %    o image: the image.
440 %
441 %    o radius: the radius of the Gaussian, in pixels, not counting the center
442 %      pixel.
443 %
444 %    o sigma: the standard deviation of the Laplacian, in pixels.
445 %
446 %    o exception: return any errors or warnings in this structure.
447 %
448 */
AdaptiveSharpenImage(const Image * image,const double radius,const double sigma,ExceptionInfo * exception)449 MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
450   const double sigma,ExceptionInfo *exception)
451 {
452 #define AdaptiveSharpenImageTag  "Convolve/Image"
453 #define MagickSigma  (fabs(sigma) < MagickEpsilon ? MagickEpsilon : sigma)
454 
455   CacheView
456     *sharp_view,
457     *edge_view,
458     *image_view;
459 
460   double
461     normalize,
462     **kernel;
463 
464   Image
465     *sharp_image,
466     *edge_image,
467     *gaussian_image;
468 
469   MagickBooleanType
470     status;
471 
472   MagickOffsetType
473     progress;
474 
475   size_t
476     width;
477 
478   ssize_t
479     w,
480     y;
481 
482   assert(image != (const Image *) NULL);
483   assert(image->signature == MagickCoreSignature);
484   if (image->debug != MagickFalse)
485     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
486   assert(exception != (ExceptionInfo *) NULL);
487   assert(exception->signature == MagickCoreSignature);
488   sharp_image=CloneImage(image,0,0,MagickTrue,exception);
489   if (sharp_image == (Image *) NULL)
490     return((Image *) NULL);
491   if (fabs(sigma) < MagickEpsilon)
492     return(sharp_image);
493   if (SetImageStorageClass(sharp_image,DirectClass,exception) == MagickFalse)
494     {
495       sharp_image=DestroyImage(sharp_image);
496       return((Image *) NULL);
497     }
498   /*
499     Edge detect the image brightness channel, level, sharp, and level again.
500   */
501   edge_image=EdgeImage(image,radius,exception);
502   if (edge_image == (Image *) NULL)
503     {
504       sharp_image=DestroyImage(sharp_image);
505       return((Image *) NULL);
506     }
507   (void) AutoLevelImage(edge_image,exception);
508   gaussian_image=BlurImage(edge_image,radius,sigma,exception);
509   if (gaussian_image != (Image *) NULL)
510     {
511       edge_image=DestroyImage(edge_image);
512       edge_image=gaussian_image;
513     }
514   (void) AutoLevelImage(edge_image,exception);
515   /*
516     Create a set of kernels from maximum (radius,sigma) to minimum.
517   */
518   width=GetOptimalKernelWidth2D(radius,sigma);
519   kernel=(double **) MagickAssumeAligned(AcquireAlignedMemory((size_t)
520     width,sizeof(*kernel)));
521   if (kernel == (double **) NULL)
522     {
523       edge_image=DestroyImage(edge_image);
524       sharp_image=DestroyImage(sharp_image);
525       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
526     }
527   (void) memset(kernel,0,(size_t) width*sizeof(*kernel));
528   for (w=0; w < (ssize_t) width; w+=2)
529   {
530     ssize_t
531       j,
532       k,
533       u,
534       v;
535 
536     kernel[w]=(double *) MagickAssumeAligned(AcquireAlignedMemory((size_t)
537       (width-w),(width-w)*sizeof(**kernel)));
538     if (kernel[w] == (double *) NULL)
539       break;
540     normalize=0.0;
541     j=(ssize_t) (width-w-1)/2;
542     k=0;
543     for (v=(-j); v <= j; v++)
544     {
545       for (u=(-j); u <= j; u++)
546       {
547         kernel[w][k]=(double) (-exp(-((double) u*u+v*v)/(2.0*MagickSigma*
548           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
549         normalize+=kernel[w][k];
550         k++;
551       }
552     }
553     kernel[w][(k-1)/2]=(double) ((-2.0)*normalize);
554     if (sigma < MagickEpsilon)
555       kernel[w][(k-1)/2]=1.0;
556   }
557   if (w < (ssize_t) width)
558     {
559       for (w-=2; w >= 0; w-=2)
560         kernel[w]=(double *) RelinquishAlignedMemory(kernel[w]);
561       kernel=(double **) RelinquishAlignedMemory(kernel);
562       edge_image=DestroyImage(edge_image);
563       sharp_image=DestroyImage(sharp_image);
564       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
565     }
566   /*
567     Adaptively sharpen image.
568   */
569   status=MagickTrue;
570   progress=0;
571   image_view=AcquireVirtualCacheView(image,exception);
572   edge_view=AcquireVirtualCacheView(edge_image,exception);
573   sharp_view=AcquireAuthenticCacheView(sharp_image,exception);
574 #if defined(MAGICKCORE_OPENMP_SUPPORT)
575   #pragma omp parallel for schedule(static) shared(progress,status) \
576     magick_number_threads(image,sharp_image,sharp_image->rows,1)
577 #endif
578   for (y=0; y < (ssize_t) sharp_image->rows; y++)
579   {
580     const Quantum
581       *magick_restrict r;
582 
583     Quantum
584       *magick_restrict q;
585 
586     ssize_t
587       x;
588 
589     if (status == MagickFalse)
590       continue;
591     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
592     q=QueueCacheViewAuthenticPixels(sharp_view,0,y,sharp_image->columns,1,
593       exception);
594     if ((r == (const Quantum *) NULL) || (q == (Quantum *) NULL))
595       {
596         status=MagickFalse;
597         continue;
598       }
599     for (x=0; x < (ssize_t) sharp_image->columns; x++)
600     {
601       const Quantum
602         *magick_restrict p;
603 
604       ssize_t
605         i;
606 
607       ssize_t
608         center,
609         j;
610 
611       j=CastDoubleToLong(ceil((double) width*(1.0-QuantumScale*
612         GetPixelIntensity(edge_image,r))-0.5));
613       if (j < 0)
614         j=0;
615       else
616         if (j > (ssize_t) width)
617           j=(ssize_t) width;
618       if ((j & 0x01) != 0)
619         j--;
620       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-j)/2L),y-
621         (ssize_t) ((width-j)/2L),width-j,width-j,exception);
622       if (p == (const Quantum *) NULL)
623         break;
624       center=(ssize_t) GetPixelChannels(image)*(width-j)*((width-j)/2L)+
625         GetPixelChannels(image)*((width-j)/2);
626       for (i=0; i < (ssize_t) GetPixelChannels(sharp_image); i++)
627       {
628         double
629           alpha,
630           gamma,
631           pixel;
632 
633         PixelChannel
634           channel;
635 
636         PixelTrait
637           sharp_traits,
638           traits;
639 
640         const double
641           *magick_restrict k;
642 
643         const Quantum
644           *magick_restrict pixels;
645 
646         ssize_t
647           u;
648 
649         ssize_t
650           v;
651 
652         channel=GetPixelChannelChannel(image,i);
653         traits=GetPixelChannelTraits(image,channel);
654         sharp_traits=GetPixelChannelTraits(sharp_image,channel);
655         if ((traits == UndefinedPixelTrait) ||
656             (sharp_traits == UndefinedPixelTrait))
657           continue;
658         if ((sharp_traits & CopyPixelTrait) != 0)
659           {
660             SetPixelChannel(sharp_image,channel,p[center+i],q);
661             continue;
662           }
663         k=kernel[j];
664         pixels=p;
665         pixel=0.0;
666         gamma=0.0;
667         if ((sharp_traits & BlendPixelTrait) == 0)
668           {
669             /*
670               No alpha blending.
671             */
672             for (v=0; v < (ssize_t) (width-j); v++)
673             {
674               for (u=0; u < (ssize_t) (width-j); u++)
675               {
676                 pixel+=(*k)*pixels[i];
677                 gamma+=(*k);
678                 k++;
679                 pixels+=GetPixelChannels(image);
680               }
681             }
682             gamma=PerceptibleReciprocal(gamma);
683             SetPixelChannel(sharp_image,channel,ClampToQuantum(gamma*pixel),q);
684             continue;
685           }
686         /*
687           Alpha blending.
688         */
689         for (v=0; v < (ssize_t) (width-j); v++)
690         {
691           for (u=0; u < (ssize_t) (width-j); u++)
692           {
693             alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
694             pixel+=(*k)*alpha*pixels[i];
695             gamma+=(*k)*alpha;
696             k++;
697             pixels+=GetPixelChannels(image);
698           }
699         }
700         gamma=PerceptibleReciprocal(gamma);
701         SetPixelChannel(sharp_image,channel,ClampToQuantum(gamma*pixel),q);
702       }
703       q+=GetPixelChannels(sharp_image);
704       r+=GetPixelChannels(edge_image);
705     }
706     if (SyncCacheViewAuthenticPixels(sharp_view,exception) == MagickFalse)
707       status=MagickFalse;
708     if (image->progress_monitor != (MagickProgressMonitor) NULL)
709       {
710         MagickBooleanType
711           proceed;
712 
713 #if defined(MAGICKCORE_OPENMP_SUPPORT)
714         #pragma omp atomic
715 #endif
716         progress++;
717         proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress,
718           image->rows);
719         if (proceed == MagickFalse)
720           status=MagickFalse;
721       }
722   }
723   sharp_image->type=image->type;
724   sharp_view=DestroyCacheView(sharp_view);
725   edge_view=DestroyCacheView(edge_view);
726   image_view=DestroyCacheView(image_view);
727   edge_image=DestroyImage(edge_image);
728   for (w=0; w < (ssize_t) width; w+=2)
729     kernel[w]=(double *) RelinquishAlignedMemory(kernel[w]);
730   kernel=(double **) RelinquishAlignedMemory(kernel);
731   if (status == MagickFalse)
732     sharp_image=DestroyImage(sharp_image);
733   return(sharp_image);
734 }
735 
736 /*
737 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
738 %                                                                             %
739 %                                                                             %
740 %                                                                             %
741 %     B l u r I m a g e                                                       %
742 %                                                                             %
743 %                                                                             %
744 %                                                                             %
745 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
746 %
747 %  BlurImage() blurs an image.  We convolve the image with a Gaussian operator
748 %  of the given radius and standard deviation (sigma).  For reasonable results,
749 %  the radius should be larger than sigma.  Use a radius of 0 and BlurImage()
750 %  selects a suitable radius for you.
751 %
752 %  The format of the BlurImage method is:
753 %
754 %      Image *BlurImage(const Image *image,const double radius,
755 %        const double sigma,ExceptionInfo *exception)
756 %
757 %  A description of each parameter follows:
758 %
759 %    o image: the image.
760 %
761 %    o radius: the radius of the Gaussian, in pixels, not counting the center
762 %      pixel.
763 %
764 %    o sigma: the standard deviation of the Gaussian, in pixels.
765 %
766 %    o exception: return any errors or warnings in this structure.
767 %
768 */
BlurImage(const Image * image,const double radius,const double sigma,ExceptionInfo * exception)769 MagickExport Image *BlurImage(const Image *image,const double radius,
770   const double sigma,ExceptionInfo *exception)
771 {
772   char
773     geometry[MagickPathExtent];
774 
775   KernelInfo
776     *kernel_info;
777 
778   Image
779     *blur_image;
780 
781   assert(image != (const Image *) NULL);
782   assert(image->signature == MagickCoreSignature);
783   if (image->debug != MagickFalse)
784     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
785   assert(exception != (ExceptionInfo *) NULL);
786   assert(exception->signature == MagickCoreSignature);
787 #if defined(MAGICKCORE_OPENCL_SUPPORT)
788   blur_image=AccelerateBlurImage(image,radius,sigma,exception);
789   if (blur_image != (Image *) NULL)
790     return(blur_image);
791 #endif
792   (void) FormatLocaleString(geometry,MagickPathExtent,
793     "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
794   kernel_info=AcquireKernelInfo(geometry,exception);
795   if (kernel_info == (KernelInfo *) NULL)
796     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
797   blur_image=ConvolveImage(image,kernel_info,exception);
798   kernel_info=DestroyKernelInfo(kernel_info);
799   return(blur_image);
800 }
801 
802 /*
803 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
804 %                                                                             %
805 %                                                                             %
806 %                                                                             %
807 %     B i l a t e r a l B l u r I m a g e                                     %
808 %                                                                             %
809 %                                                                             %
810 %                                                                             %
811 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
812 %
813 %  BilateralBlurImage() is a non-linear, edge-preserving, and noise-reducing
814 %  smoothing filter for images.  It replaces the intensity of each pixel with
815 %  a weighted average of intensity values from nearby pixels.  This weight is
816 %  based on a Gaussian distribution.  The weights depend not only on Euclidean
817 %  distance of pixels, but also on the radiometric differences (e.g., range
818 %  differences, such as color intensity, depth distance, etc.). This preserves
819 %  sharp edges.
820 %
821 %  The format of the BilateralBlurImage method is:
822 %
823 %      Image *BilateralBlurImage(const Image *image,const size_t width,
824 %        const size_t height,const double intensity_sigma,
825 %        const double spatial_sigma,ExceptionInfo *exception)
826 %
827 %  A description of each parameter follows:
828 %
829 %    o image: the image.
830 %
831 %    o width: the width of the neighborhood in pixels.
832 %
833 %    o height: the height of the neighborhood in pixels.
834 %
835 %    o intensity_sigma: sigma in the intensity space. A larger value means
836 %      that farther colors within the pixel neighborhood (see spatial_sigma)
837 %      will be mixed together, resulting in larger areas of semi-equal color.
838 %
839 %    o spatial_sigma: sigma in the coordinate space. A larger value means that
840 %      farther pixels influence each other as long as their colors are close
841 %      enough (see intensity_sigma ). When the neigborhood diameter is greater
842 %      than zero, it specifies the neighborhood size regardless of
843 %      spatial_sigma. Otherwise, the neigborhood diameter is proportional to
844 %      spatial_sigma.
845 %
846 %    o exception: return any errors or warnings in this structure.
847 %
848 */
849 
BlurDistance(const ssize_t x,const ssize_t y,const ssize_t u,const ssize_t v)850 static inline double BlurDistance(const ssize_t x,const ssize_t y,
851   const ssize_t u,const ssize_t v)
852 {
853   return(sqrt(((double) x-u)*((double) x-u)+((double) y-v)*((double) y-v)));
854 }
855 
BlurGaussian(const double x,const double sigma)856 static inline double BlurGaussian(const double x,const double sigma)
857 {
858   return(exp(-((double) x*x)*PerceptibleReciprocal(2.0*sigma*sigma))*
859     PerceptibleReciprocal(Magick2PI*sigma*sigma));
860 }
861 
DestroyBilateralThreadSet(const ssize_t number_threads,double ** weights)862 static double **DestroyBilateralThreadSet(const ssize_t number_threads,
863   double **weights)
864 {
865   ssize_t
866     i;
867 
868   assert(weights != (double **) NULL);
869   for (i=0; i <= (ssize_t) number_threads; i++)
870     if (weights[i] != (double *) NULL)
871       weights[i]=(double *) RelinquishMagickMemory(weights[i]);
872   weights=(double **) RelinquishMagickMemory(weights);
873   return(weights);
874 }
875 
AcquireBilateralThreadSet(const size_t number_threads,const size_t width,const size_t height)876 static double **AcquireBilateralThreadSet(const size_t number_threads,
877   const size_t width,const size_t height)
878 {
879   double
880     **weights;
881 
882   ssize_t
883     i;
884 
885   weights=(double **) AcquireQuantumMemory(number_threads+1,sizeof(*weights));
886   if (weights == (double **) NULL)
887     return((double **) NULL);
888   (void) memset(weights,0,number_threads*sizeof(*weights));
889   for (i=0; i <= (ssize_t) number_threads; i++)
890   {
891     weights[i]=(double *) AcquireQuantumMemory(width,height*sizeof(**weights));
892     if (weights[i] == (double *) NULL)
893       return(DestroyBilateralThreadSet(number_threads,weights));
894   }
895   return(weights);
896 }
897 
BilateralBlurImage(const Image * image,const size_t width,const size_t height,const double intensity_sigma,const double spatial_sigma,ExceptionInfo * exception)898 MagickExport Image *BilateralBlurImage(const Image *image,const size_t width,
899   const size_t height,const double intensity_sigma,const double spatial_sigma,
900   ExceptionInfo *exception)
901 {
902 #define MaxIntensity  (255)
903 #define BilateralBlurImageTag  "Blur/Image"
904 
905   CacheView
906     *blur_view,
907     *image_view;
908 
909   double
910     intensity_gaussian[2*(MaxIntensity+1)],
911     *spatial_gaussian,
912     **weights;
913 
914   Image
915     *blur_image;
916 
917   MagickBooleanType
918     status;
919 
920   MagickOffsetType
921     progress;
922 
923   OffsetInfo
924     mid;
925 
926   ssize_t
927     number_threads,
928     w,
929     y;
930 
931   assert(image != (const Image *) NULL);
932   assert(image->signature == MagickCoreSignature);
933   if (image->debug != MagickFalse)
934     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
935   assert(exception != (ExceptionInfo *) NULL);
936   assert(exception->signature == MagickCoreSignature);
937   blur_image=CloneImage(image,0,0,MagickTrue,exception);
938   if (blur_image == (Image *) NULL)
939     return((Image *) NULL);
940   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
941     {
942       blur_image=DestroyImage(blur_image);
943       return((Image *) NULL);
944     }
945   number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
946   weights=AcquireBilateralThreadSet(number_threads,width,height);
947   if (weights == (double **) NULL)
948     {
949       blur_image=DestroyImage(blur_image);
950       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
951     }
952   for (w=(-MaxIntensity); w < MaxIntensity; w++)
953     intensity_gaussian[w+MaxIntensity]=BlurGaussian((double) w,intensity_sigma);
954   spatial_gaussian=weights[number_threads];
955   {
956     ssize_t
957       n,
958       v;
959 
960     n=0;
961     mid.x=(ssize_t) (width/2L);
962     mid.y=(ssize_t) (height/2L);
963     for (v=0; v < (ssize_t) height; v++)
964     {
965       ssize_t
966         u;
967 
968       for (u=0; u < (ssize_t) width; u++)
969         spatial_gaussian[n++]=BlurGaussian(BlurDistance(0,0,u-mid.x,v-mid.y),
970           spatial_sigma);
971     }
972   }
973   /*
974     Bilateral blur image.
975   */
976   status=MagickTrue;
977   progress=0;
978   image_view=AcquireVirtualCacheView(image,exception);
979   blur_view=AcquireAuthenticCacheView(blur_image,exception);
980 #if defined(MAGICKCORE_OPENMP_SUPPORT)
981   #pragma omp parallel for schedule(static) shared(progress,status) \
982     magick_number_threads(image,blur_image,blur_image->rows,1)
983 #endif
984   for (y=0; y < (ssize_t) blur_image->rows; y++)
985   {
986     const int
987       id = GetOpenMPThreadId();
988 
989     Quantum
990       *magick_restrict q;
991 
992     ssize_t
993       x;
994 
995     if (status == MagickFalse)
996       continue;
997     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
998       exception);
999     if (q == (Quantum *) NULL)
1000       {
1001         status=MagickFalse;
1002         continue;
1003       }
1004     for (x=0; x < (ssize_t) blur_image->columns; x++)
1005     {
1006       double
1007         gamma,
1008         pixel;
1009 
1010       const Quantum
1011         *magick_restrict p,
1012         *magick_restrict r;
1013 
1014       ssize_t
1015         i,
1016         u;
1017 
1018       ssize_t
1019         n,
1020         v;
1021 
1022       /*
1023         Tonal weighting preserves edges while smoothing in the flat regions.
1024       */
1025       p=GetCacheViewVirtualPixels(image_view,x-mid.x,y-mid.y,width,height,
1026         exception);
1027       if (p == (const Quantum *) NULL)
1028         break;
1029       p+=(ssize_t) GetPixelChannels(image)*width*mid.y+GetPixelChannels(image)*
1030         mid.x;
1031       n=0;
1032       for (v=0; v < (ssize_t) height; v++)
1033       {
1034         for (u=0; u < (ssize_t) width; u++)
1035         {
1036           double
1037             intensity;
1038 
1039           r=p+(ssize_t) GetPixelChannels(image)*(ssize_t) width*(mid.y-v)+
1040             GetPixelChannels(image)*(mid.x-u);
1041           intensity=ScaleQuantumToChar(GetPixelIntensity(image,r))-
1042             (double) ScaleQuantumToChar(GetPixelIntensity(image,p));
1043           if ((intensity >= -MaxIntensity) && (intensity <= MaxIntensity))
1044             weights[id][n]=intensity_gaussian[(ssize_t) intensity+MaxIntensity]*
1045               spatial_gaussian[n];
1046           else
1047             weights[id][n]=BlurGaussian(intensity,intensity_sigma)*
1048               BlurGaussian(BlurDistance(x,y,x+u-mid.x,y+v-mid.y),spatial_sigma);
1049           n++;
1050         }
1051       }
1052       for (i=0; i < (ssize_t) GetPixelChannels(blur_image); i++)
1053       {
1054         PixelChannel
1055           channel;
1056 
1057         PixelTrait
1058           blur_traits,
1059           traits;
1060 
1061         channel=GetPixelChannelChannel(image,i);
1062         traits=GetPixelChannelTraits(image,channel);
1063         blur_traits=GetPixelChannelTraits(blur_image,channel);
1064         if ((traits == UndefinedPixelTrait) ||
1065             (blur_traits == UndefinedPixelTrait))
1066           continue;
1067         if ((blur_traits & CopyPixelTrait) != 0)
1068           {
1069             SetPixelChannel(blur_image,channel,p[i],q);
1070             continue;
1071           }
1072         pixel=0.0;
1073         gamma=0.0;
1074         n=0;
1075         if ((blur_traits & BlendPixelTrait) == 0)
1076           {
1077             /*
1078               No alpha blending.
1079             */
1080             for (v=0; v < (ssize_t) height; v++)
1081             {
1082               for (u=0; u < (ssize_t) width; u++)
1083               {
1084                 r=p+(ssize_t) GetPixelChannels(image)*width*(mid.y-v)+
1085                   GetPixelChannels(image)*(mid.x-u);
1086                 pixel+=weights[id][n]*r[i];
1087                 gamma+=weights[id][n];
1088                 n++;
1089               }
1090             }
1091             SetPixelChannel(blur_image,channel,ClampToQuantum(
1092               PerceptibleReciprocal(gamma)*pixel),q);
1093             continue;
1094           }
1095         /*
1096           Alpha blending.
1097         */
1098         for (v=0; v < (ssize_t) height; v++)
1099         {
1100           for (u=0; u < (ssize_t) width; u++)
1101           {
1102             double
1103               alpha,
1104               beta;
1105 
1106             r=p+(ssize_t) GetPixelChannels(image)*width*(mid.y-v)+
1107               GetPixelChannels(image)*(mid.x-u);
1108             alpha=(double) (QuantumScale*GetPixelAlpha(image,p));
1109             beta=(double) (QuantumScale*GetPixelAlpha(image,r));
1110             pixel+=weights[id][n]*r[i];
1111             gamma+=weights[id][n]*alpha*beta;
1112             n++;
1113           }
1114         }
1115         SetPixelChannel(blur_image,channel,ClampToQuantum(
1116           PerceptibleReciprocal(gamma)*pixel),q);
1117       }
1118       q+=GetPixelChannels(blur_image);
1119     }
1120     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
1121       status=MagickFalse;
1122     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1123       {
1124         MagickBooleanType
1125           proceed;
1126 
1127 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1128         #pragma omp atomic
1129 #endif
1130         progress++;
1131         proceed=SetImageProgress(image,BilateralBlurImageTag,progress,
1132           image->rows);
1133         if (proceed == MagickFalse)
1134           status=MagickFalse;
1135       }
1136   }
1137   blur_image->type=image->type;
1138   blur_view=DestroyCacheView(blur_view);
1139   image_view=DestroyCacheView(image_view);
1140   weights=DestroyBilateralThreadSet(number_threads,weights);
1141   if (status == MagickFalse)
1142     blur_image=DestroyImage(blur_image);
1143   return(blur_image);
1144 }
1145 
1146 /*
1147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1148 %                                                                             %
1149 %                                                                             %
1150 %                                                                             %
1151 %     C o n v o l v e I m a g e                                               %
1152 %                                                                             %
1153 %                                                                             %
1154 %                                                                             %
1155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1156 %
1157 %  ConvolveImage() applies a custom convolution kernel to the image.
1158 %
1159 %  The format of the ConvolveImage method is:
1160 %
1161 %      Image *ConvolveImage(const Image *image,const KernelInfo *kernel,
1162 %        ExceptionInfo *exception)
1163 %
1164 %  A description of each parameter follows:
1165 %
1166 %    o image: the image.
1167 %
1168 %    o kernel: the filtering kernel.
1169 %
1170 %    o exception: return any errors or warnings in this structure.
1171 %
1172 */
ConvolveImage(const Image * image,const KernelInfo * kernel_info,ExceptionInfo * exception)1173 MagickExport Image *ConvolveImage(const Image *image,
1174   const KernelInfo *kernel_info,ExceptionInfo *exception)
1175 {
1176   Image
1177     *convolve_image;
1178 
1179 #if defined(MAGICKCORE_OPENCL_SUPPORT)
1180   convolve_image=AccelerateConvolveImage(image,kernel_info,exception);
1181   if (convolve_image != (Image *) NULL)
1182     return(convolve_image);
1183 #endif
1184 
1185   convolve_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,
1186     exception);
1187   return(convolve_image);
1188 }
1189 
1190 /*
1191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1192 %                                                                             %
1193 %                                                                             %
1194 %                                                                             %
1195 %     D e s p e c k l e I m a g e                                             %
1196 %                                                                             %
1197 %                                                                             %
1198 %                                                                             %
1199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1200 %
1201 %  DespeckleImage() reduces the speckle noise in an image while perserving the
1202 %  edges of the original image.  A speckle removing filter uses a complementary
1203 %  hulling technique (raising pixels that are darker than their surrounding
1204 %  neighbors, then complementarily lowering pixels that are brighter than their
1205 %  surrounding neighbors) to reduce the speckle index of that image (reference
1206 %  Crimmins speckle removal).
1207 %
1208 %  The format of the DespeckleImage method is:
1209 %
1210 %      Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1211 %
1212 %  A description of each parameter follows:
1213 %
1214 %    o image: the image.
1215 %
1216 %    o exception: return any errors or warnings in this structure.
1217 %
1218 */
1219 
Hull(const Image * image,const ssize_t x_offset,const ssize_t y_offset,const size_t columns,const size_t rows,const int polarity,Quantum * magick_restrict f,Quantum * magick_restrict g)1220 static void Hull(const Image *image,const ssize_t x_offset,
1221   const ssize_t y_offset,const size_t columns,const size_t rows,
1222   const int polarity,Quantum *magick_restrict f,Quantum *magick_restrict g)
1223 {
1224   Quantum
1225     *p,
1226     *q,
1227     *r,
1228     *s;
1229 
1230   ssize_t
1231     y;
1232 
1233   assert(image != (const Image *) NULL);
1234   assert(image->signature == MagickCoreSignature);
1235   if (image->debug != MagickFalse)
1236     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1237   assert(f != (Quantum *) NULL);
1238   assert(g != (Quantum *) NULL);
1239   p=f+(columns+2);
1240   q=g+(columns+2);
1241   r=p+(y_offset*((ssize_t) columns+2)+x_offset);
1242 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1243   #pragma omp parallel for schedule(static) \
1244     magick_number_threads(image,image,rows,1)
1245 #endif
1246   for (y=0; y < (ssize_t) rows; y++)
1247   {
1248     MagickRealType
1249       v;
1250 
1251     ssize_t
1252       i,
1253       x;
1254 
1255     i=(2*y+1)+y*columns;
1256     if (polarity > 0)
1257       for (x=0; x < (ssize_t) columns; x++)
1258       {
1259         v=(MagickRealType) p[i];
1260         if ((MagickRealType) r[i] >= (v+ScaleCharToQuantum(2)))
1261           v+=ScaleCharToQuantum(1);
1262         q[i]=(Quantum) v;
1263         i++;
1264       }
1265     else
1266       for (x=0; x < (ssize_t) columns; x++)
1267       {
1268         v=(MagickRealType) p[i];
1269         if ((MagickRealType) r[i] <= (v-ScaleCharToQuantum(2)))
1270           v-=ScaleCharToQuantum(1);
1271         q[i]=(Quantum) v;
1272         i++;
1273       }
1274   }
1275   p=f+(columns+2);
1276   q=g+(columns+2);
1277   r=q+(y_offset*((ssize_t) columns+2)+x_offset);
1278   s=q-(y_offset*((ssize_t) columns+2)+x_offset);
1279 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1280   #pragma omp parallel for schedule(static) \
1281     magick_number_threads(image,image,rows,1)
1282 #endif
1283   for (y=0; y < (ssize_t) rows; y++)
1284   {
1285     ssize_t
1286       i,
1287       x;
1288 
1289     MagickRealType
1290       v;
1291 
1292     i=(2*y+1)+y*columns;
1293     if (polarity > 0)
1294       for (x=0; x < (ssize_t) columns; x++)
1295       {
1296         v=(MagickRealType) q[i];
1297         if (((MagickRealType) s[i] >= (v+ScaleCharToQuantum(2))) &&
1298             ((MagickRealType) r[i] > v))
1299           v+=ScaleCharToQuantum(1);
1300         p[i]=(Quantum) v;
1301         i++;
1302       }
1303     else
1304       for (x=0; x < (ssize_t) columns; x++)
1305       {
1306         v=(MagickRealType) q[i];
1307         if (((MagickRealType) s[i] <= (v-ScaleCharToQuantum(2))) &&
1308             ((MagickRealType) r[i] < v))
1309           v-=ScaleCharToQuantum(1);
1310         p[i]=(Quantum) v;
1311         i++;
1312       }
1313   }
1314 }
1315 
DespeckleImage(const Image * image,ExceptionInfo * exception)1316 MagickExport Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1317 {
1318 #define DespeckleImageTag  "Despeckle/Image"
1319 
1320   CacheView
1321     *despeckle_view,
1322     *image_view;
1323 
1324   Image
1325     *despeckle_image;
1326 
1327   MagickBooleanType
1328     status;
1329 
1330   MemoryInfo
1331     *buffer_info,
1332     *pixel_info;
1333 
1334   Quantum
1335     *magick_restrict buffer,
1336     *magick_restrict pixels;
1337 
1338   ssize_t
1339     i;
1340 
1341   size_t
1342     length;
1343 
1344   static const ssize_t
1345     X[4] = {0, 1, 1,-1},
1346     Y[4] = {1, 0, 1, 1};
1347 
1348   /*
1349     Allocate despeckled image.
1350   */
1351   assert(image != (const Image *) NULL);
1352   assert(image->signature == MagickCoreSignature);
1353   if (image->debug != MagickFalse)
1354     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1355   assert(exception != (ExceptionInfo *) NULL);
1356   assert(exception->signature == MagickCoreSignature);
1357 #if defined(MAGICKCORE_OPENCL_SUPPORT)
1358   despeckle_image=AccelerateDespeckleImage(image,exception);
1359   if (despeckle_image != (Image *) NULL)
1360     return(despeckle_image);
1361 #endif
1362   despeckle_image=CloneImage(image,0,0,MagickTrue,exception);
1363   if (despeckle_image == (Image *) NULL)
1364     return((Image *) NULL);
1365   status=SetImageStorageClass(despeckle_image,DirectClass,exception);
1366   if (status == MagickFalse)
1367     {
1368       despeckle_image=DestroyImage(despeckle_image);
1369       return((Image *) NULL);
1370     }
1371   /*
1372     Allocate image buffer.
1373   */
1374   length=(size_t) ((image->columns+2)*(image->rows+2));
1375   pixel_info=AcquireVirtualMemory(length,sizeof(*pixels));
1376   buffer_info=AcquireVirtualMemory(length,sizeof(*buffer));
1377   if ((pixel_info == (MemoryInfo *) NULL) ||
1378       (buffer_info == (MemoryInfo *) NULL))
1379     {
1380       if (buffer_info != (MemoryInfo *) NULL)
1381         buffer_info=RelinquishVirtualMemory(buffer_info);
1382       if (pixel_info != (MemoryInfo *) NULL)
1383         pixel_info=RelinquishVirtualMemory(pixel_info);
1384       despeckle_image=DestroyImage(despeckle_image);
1385       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1386     }
1387   pixels=(Quantum *) GetVirtualMemoryBlob(pixel_info);
1388   buffer=(Quantum *) GetVirtualMemoryBlob(buffer_info);
1389   /*
1390     Reduce speckle in the image.
1391   */
1392   status=MagickTrue;
1393   image_view=AcquireVirtualCacheView(image,exception);
1394   despeckle_view=AcquireAuthenticCacheView(despeckle_image,exception);
1395   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1396   {
1397     PixelChannel
1398        channel;
1399 
1400     PixelTrait
1401       despeckle_traits,
1402       traits;
1403 
1404     ssize_t
1405       k,
1406       x;
1407 
1408     ssize_t
1409       j,
1410       y;
1411 
1412     if (status == MagickFalse)
1413       continue;
1414     channel=GetPixelChannelChannel(image,i);
1415     traits=GetPixelChannelTraits(image,channel);
1416     despeckle_traits=GetPixelChannelTraits(despeckle_image,channel);
1417     if ((traits == UndefinedPixelTrait) ||
1418         (despeckle_traits == UndefinedPixelTrait))
1419       continue;
1420     if ((despeckle_traits & CopyPixelTrait) != 0)
1421       continue;
1422     (void) memset(pixels,0,length*sizeof(*pixels));
1423     j=(ssize_t) image->columns+2;
1424     for (y=0; y < (ssize_t) image->rows; y++)
1425     {
1426       const Quantum
1427         *magick_restrict p;
1428 
1429       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1430       if (p == (const Quantum *) NULL)
1431         {
1432           status=MagickFalse;
1433           continue;
1434         }
1435       j++;
1436       for (x=0; x < (ssize_t) image->columns; x++)
1437       {
1438         pixels[j++]=p[i];
1439         p+=GetPixelChannels(image);
1440       }
1441       j++;
1442     }
1443     (void) memset(buffer,0,length*sizeof(*buffer));
1444     for (k=0; k < 4; k++)
1445     {
1446       Hull(image,X[k],Y[k],image->columns,image->rows,1,pixels,buffer);
1447       Hull(image,-X[k],-Y[k],image->columns,image->rows,1,pixels,buffer);
1448       Hull(image,-X[k],-Y[k],image->columns,image->rows,-1,pixels,buffer);
1449       Hull(image,X[k],Y[k],image->columns,image->rows,-1,pixels,buffer);
1450     }
1451     j=(ssize_t) image->columns+2;
1452     for (y=0; y < (ssize_t) image->rows; y++)
1453     {
1454       MagickBooleanType
1455         sync;
1456 
1457       Quantum
1458         *magick_restrict q;
1459 
1460       q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
1461         1,exception);
1462       if (q == (Quantum *) NULL)
1463         {
1464           status=MagickFalse;
1465           continue;
1466         }
1467       j++;
1468       for (x=0; x < (ssize_t) image->columns; x++)
1469       {
1470         SetPixelChannel(despeckle_image,channel,pixels[j++],q);
1471         q+=GetPixelChannels(despeckle_image);
1472       }
1473       sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
1474       if (sync == MagickFalse)
1475         status=MagickFalse;
1476       j++;
1477     }
1478     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1479       {
1480         MagickBooleanType
1481           proceed;
1482 
1483         proceed=SetImageProgress(image,DespeckleImageTag,(MagickOffsetType) i,
1484           GetPixelChannels(image));
1485         if (proceed == MagickFalse)
1486           status=MagickFalse;
1487       }
1488   }
1489   despeckle_view=DestroyCacheView(despeckle_view);
1490   image_view=DestroyCacheView(image_view);
1491   buffer_info=RelinquishVirtualMemory(buffer_info);
1492   pixel_info=RelinquishVirtualMemory(pixel_info);
1493   despeckle_image->type=image->type;
1494   if (status == MagickFalse)
1495     despeckle_image=DestroyImage(despeckle_image);
1496   return(despeckle_image);
1497 }
1498 
1499 /*
1500 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1501 %                                                                             %
1502 %                                                                             %
1503 %                                                                             %
1504 %     E d g e I m a g e                                                       %
1505 %                                                                             %
1506 %                                                                             %
1507 %                                                                             %
1508 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1509 %
1510 %  EdgeImage() finds edges in an image.  Radius defines the radius of the
1511 %  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
1512 %  radius for you.
1513 %
1514 %  The format of the EdgeImage method is:
1515 %
1516 %      Image *EdgeImage(const Image *image,const double radius,
1517 %        ExceptionInfo *exception)
1518 %
1519 %  A description of each parameter follows:
1520 %
1521 %    o image: the image.
1522 %
1523 %    o radius: the radius of the pixel neighborhood.
1524 %
1525 %    o exception: return any errors or warnings in this structure.
1526 %
1527 */
EdgeImage(const Image * image,const double radius,ExceptionInfo * exception)1528 MagickExport Image *EdgeImage(const Image *image,const double radius,
1529   ExceptionInfo *exception)
1530 {
1531   Image
1532     *edge_image;
1533 
1534   KernelInfo
1535     *kernel_info;
1536 
1537   ssize_t
1538     i;
1539 
1540   size_t
1541     width;
1542 
1543   assert(image != (const Image *) NULL);
1544   assert(image->signature == MagickCoreSignature);
1545   if (image->debug != MagickFalse)
1546     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1547   assert(exception != (ExceptionInfo *) NULL);
1548   assert(exception->signature == MagickCoreSignature);
1549   width=GetOptimalKernelWidth1D(radius,0.5);
1550   kernel_info=AcquireKernelInfo((const char *) NULL,exception);
1551   if (kernel_info == (KernelInfo *) NULL)
1552     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1553   (void) memset(kernel_info,0,sizeof(*kernel_info));
1554   kernel_info->width=width;
1555   kernel_info->height=width;
1556   kernel_info->x=(ssize_t) (kernel_info->width-1)/2;
1557   kernel_info->y=(ssize_t) (kernel_info->height-1)/2;
1558   kernel_info->signature=MagickCoreSignature;
1559   kernel_info->values=(MagickRealType *) MagickAssumeAligned(
1560     AcquireAlignedMemory(kernel_info->width,kernel_info->height*
1561     sizeof(*kernel_info->values)));
1562   if (kernel_info->values == (MagickRealType *) NULL)
1563     {
1564       kernel_info=DestroyKernelInfo(kernel_info);
1565       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1566     }
1567   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
1568     kernel_info->values[i]=(-1.0);
1569   kernel_info->values[i/2]=(double) kernel_info->width*kernel_info->height-1.0;
1570   edge_image=ConvolveImage(image,kernel_info,exception);
1571   kernel_info=DestroyKernelInfo(kernel_info);
1572   return(edge_image);
1573 }
1574 
1575 /*
1576 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1577 %                                                                             %
1578 %                                                                             %
1579 %                                                                             %
1580 %     E m b o s s I m a g e                                                   %
1581 %                                                                             %
1582 %                                                                             %
1583 %                                                                             %
1584 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1585 %
1586 %  EmbossImage() returns a grayscale image with a three-dimensional effect.
1587 %  We convolve the image with a Gaussian operator of the given radius and
1588 %  standard deviation (sigma).  For reasonable results, radius should be
1589 %  larger than sigma.  Use a radius of 0 and Emboss() selects a suitable
1590 %  radius for you.
1591 %
1592 %  The format of the EmbossImage method is:
1593 %
1594 %      Image *EmbossImage(const Image *image,const double radius,
1595 %        const double sigma,ExceptionInfo *exception)
1596 %
1597 %  A description of each parameter follows:
1598 %
1599 %    o image: the image.
1600 %
1601 %    o radius: the radius of the pixel neighborhood.
1602 %
1603 %    o sigma: the standard deviation of the Gaussian, in pixels.
1604 %
1605 %    o exception: return any errors or warnings in this structure.
1606 %
1607 */
EmbossImage(const Image * image,const double radius,const double sigma,ExceptionInfo * exception)1608 MagickExport Image *EmbossImage(const Image *image,const double radius,
1609   const double sigma,ExceptionInfo *exception)
1610 {
1611   double
1612     gamma,
1613     normalize;
1614 
1615   Image
1616     *emboss_image;
1617 
1618   KernelInfo
1619     *kernel_info;
1620 
1621   ssize_t
1622     i;
1623 
1624   size_t
1625     width;
1626 
1627   ssize_t
1628     j,
1629     k,
1630     u,
1631     v;
1632 
1633   assert(image != (const Image *) NULL);
1634   assert(image->signature == MagickCoreSignature);
1635   if (image->debug != MagickFalse)
1636     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1637   assert(exception != (ExceptionInfo *) NULL);
1638   assert(exception->signature == MagickCoreSignature);
1639   width=GetOptimalKernelWidth1D(radius,sigma);
1640   kernel_info=AcquireKernelInfo((const char *) NULL,exception);
1641   if (kernel_info == (KernelInfo *) NULL)
1642     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1643   kernel_info->width=width;
1644   kernel_info->height=width;
1645   kernel_info->x=(ssize_t) (width-1)/2;
1646   kernel_info->y=(ssize_t) (width-1)/2;
1647   kernel_info->values=(MagickRealType *) MagickAssumeAligned(
1648     AcquireAlignedMemory(kernel_info->width,kernel_info->width*
1649     sizeof(*kernel_info->values)));
1650   if (kernel_info->values == (MagickRealType *) NULL)
1651     {
1652       kernel_info=DestroyKernelInfo(kernel_info);
1653       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1654     }
1655   j=(ssize_t) (kernel_info->width-1)/2;
1656   k=j;
1657   i=0;
1658   for (v=(-j); v <= j; v++)
1659   {
1660     for (u=(-j); u <= j; u++)
1661     {
1662       kernel_info->values[i]=(MagickRealType) (((u < 0) || (v < 0) ? -8.0 :
1663         8.0)*exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
1664         (2.0*MagickPI*MagickSigma*MagickSigma));
1665       if (u != k)
1666         kernel_info->values[i]=0.0;
1667       i++;
1668     }
1669     k--;
1670   }
1671   normalize=0.0;
1672   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
1673     normalize+=kernel_info->values[i];
1674   gamma=PerceptibleReciprocal(normalize);
1675   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
1676     kernel_info->values[i]*=gamma;
1677   emboss_image=ConvolveImage(image,kernel_info,exception);
1678   kernel_info=DestroyKernelInfo(kernel_info);
1679   if (emboss_image != (Image *) NULL)
1680     (void) EqualizeImage(emboss_image,exception);
1681   return(emboss_image);
1682 }
1683 
1684 /*
1685 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1686 %                                                                             %
1687 %                                                                             %
1688 %                                                                             %
1689 %     G a u s s i a n B l u r I m a g e                                       %
1690 %                                                                             %
1691 %                                                                             %
1692 %                                                                             %
1693 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1694 %
1695 %  GaussianBlurImage() blurs an image.  We convolve the image with a
1696 %  Gaussian operator of the given radius and standard deviation (sigma).
1697 %  For reasonable results, the radius should be larger than sigma.  Use a
1698 %  radius of 0 and GaussianBlurImage() selects a suitable radius for you.
1699 %
1700 %  The format of the GaussianBlurImage method is:
1701 %
1702 %      Image *GaussianBlurImage(const Image *image,onst double radius,
1703 %        const double sigma,ExceptionInfo *exception)
1704 %
1705 %  A description of each parameter follows:
1706 %
1707 %    o image: the image.
1708 %
1709 %    o radius: the radius of the Gaussian, in pixels, not counting the center
1710 %      pixel.
1711 %
1712 %    o sigma: the standard deviation of the Gaussian, in pixels.
1713 %
1714 %    o exception: return any errors or warnings in this structure.
1715 %
1716 */
GaussianBlurImage(const Image * image,const double radius,const double sigma,ExceptionInfo * exception)1717 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
1718   const double sigma,ExceptionInfo *exception)
1719 {
1720   char
1721     geometry[MagickPathExtent];
1722 
1723   KernelInfo
1724     *kernel_info;
1725 
1726   Image
1727     *blur_image;
1728 
1729   assert(image != (const Image *) NULL);
1730   assert(image->signature == MagickCoreSignature);
1731   if (image->debug != MagickFalse)
1732     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1733   assert(exception != (ExceptionInfo *) NULL);
1734   assert(exception->signature == MagickCoreSignature);
1735   (void) FormatLocaleString(geometry,MagickPathExtent,"gaussian:%.20gx%.20g",
1736     radius,sigma);
1737   kernel_info=AcquireKernelInfo(geometry,exception);
1738   if (kernel_info == (KernelInfo *) NULL)
1739     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1740   blur_image=ConvolveImage(image,kernel_info,exception);
1741   kernel_info=DestroyKernelInfo(kernel_info);
1742   return(blur_image);
1743 }
1744 
1745 /*
1746 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1747 %                                                                             %
1748 %                                                                             %
1749 %                                                                             %
1750 %     K u w a h a r a I m a g e                                               %
1751 %                                                                             %
1752 %                                                                             %
1753 %                                                                             %
1754 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1755 %
1756 %  KuwaharaImage() is an edge preserving noise reduction filter.
1757 %
1758 %  The format of the KuwaharaImage method is:
1759 %
1760 %      Image *KuwaharaImage(const Image *image,const double radius,
1761 %        const double sigma,ExceptionInfo *exception)
1762 %
1763 %  A description of each parameter follows:
1764 %
1765 %    o image: the image.
1766 %
1767 %    o radius: the square window radius.
1768 %
1769 %    o sigma: the standard deviation of the Gaussian, in pixels.
1770 %
1771 %    o exception: return any errors or warnings in this structure.
1772 %
1773 */
1774 
GetMeanLuma(const Image * magick_restrict image,const double * magick_restrict pixel)1775 static inline MagickRealType GetMeanLuma(const Image *magick_restrict image,
1776   const double *magick_restrict pixel)
1777 {
1778   return(0.212656f*pixel[image->channel_map[RedPixelChannel].offset]+
1779     0.715158f*pixel[image->channel_map[GreenPixelChannel].offset]+
1780     0.072186f*pixel[image->channel_map[BluePixelChannel].offset]);  /* Rec709 */
1781 }
1782 
KuwaharaImage(const Image * image,const double radius,const double sigma,ExceptionInfo * exception)1783 MagickExport Image *KuwaharaImage(const Image *image,const double radius,
1784   const double sigma,ExceptionInfo *exception)
1785 {
1786 #define KuwaharaImageTag  "Kuwahara/Image"
1787 
1788   CacheView
1789     *image_view,
1790     *kuwahara_view;
1791 
1792   Image
1793     *gaussian_image,
1794     *kuwahara_image;
1795 
1796   MagickBooleanType
1797     status;
1798 
1799   MagickOffsetType
1800     progress;
1801 
1802   size_t
1803     width;
1804 
1805   ssize_t
1806     y;
1807 
1808   /*
1809     Initialize Kuwahara image attributes.
1810   */
1811   assert(image != (Image *) NULL);
1812   assert(image->signature == MagickCoreSignature);
1813   if (image->debug != MagickFalse)
1814     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1815   assert(exception != (ExceptionInfo *) NULL);
1816   assert(exception->signature == MagickCoreSignature);
1817   width=(size_t) radius+1;
1818   gaussian_image=BlurImage(image,radius,sigma,exception);
1819   if (gaussian_image == (Image *) NULL)
1820     return((Image *) NULL);
1821   kuwahara_image=CloneImage(image,0,0,MagickTrue,exception);
1822   if (kuwahara_image == (Image *) NULL)
1823     {
1824       gaussian_image=DestroyImage(gaussian_image);
1825       return((Image *) NULL);
1826     }
1827   if (SetImageStorageClass(kuwahara_image,DirectClass,exception) == MagickFalse)
1828     {
1829       gaussian_image=DestroyImage(gaussian_image);
1830       kuwahara_image=DestroyImage(kuwahara_image);
1831       return((Image *) NULL);
1832     }
1833   /*
1834     Edge preserving noise reduction filter.
1835   */
1836   status=MagickTrue;
1837   progress=0;
1838   image_view=AcquireVirtualCacheView(gaussian_image,exception);
1839   kuwahara_view=AcquireAuthenticCacheView(kuwahara_image,exception);
1840 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1841   #pragma omp parallel for schedule(static) shared(progress,status) \
1842     magick_number_threads(image,kuwahara_image,gaussian_image->rows,1)
1843 #endif
1844   for (y=0; y < (ssize_t) gaussian_image->rows; y++)
1845   {
1846     Quantum
1847       *magick_restrict q;
1848 
1849     ssize_t
1850       x;
1851 
1852     if (status == MagickFalse)
1853       continue;
1854     q=QueueCacheViewAuthenticPixels(kuwahara_view,0,y,kuwahara_image->columns,1,
1855       exception);
1856     if (q == (Quantum *) NULL)
1857       {
1858         status=MagickFalse;
1859         continue;
1860       }
1861     for (x=0; x < (ssize_t) gaussian_image->columns; x++)
1862     {
1863       const Quantum
1864         *magick_restrict p;
1865 
1866       double
1867         min_variance;
1868 
1869       RectangleInfo
1870         quadrant,
1871         target;
1872 
1873       size_t
1874         i;
1875 
1876       min_variance=MagickMaximumValue;
1877       SetGeometry(gaussian_image,&target);
1878       quadrant.width=width;
1879       quadrant.height=width;
1880       for (i=0; i < 4; i++)
1881       {
1882         const Quantum
1883           *magick_restrict k;
1884 
1885         double
1886           mean[MaxPixelChannels],
1887           variance;
1888 
1889         ssize_t
1890           n;
1891 
1892         ssize_t
1893           j;
1894 
1895         quadrant.x=x;
1896         quadrant.y=y;
1897         switch (i)
1898         {
1899           case 0:
1900           {
1901             quadrant.x=x-(ssize_t) (width-1);
1902             quadrant.y=y-(ssize_t) (width-1);
1903             break;
1904           }
1905           case 1:
1906           {
1907             quadrant.y=y-(ssize_t) (width-1);
1908             break;
1909           }
1910           case 2:
1911           {
1912             quadrant.x=x-(ssize_t) (width-1);
1913             break;
1914           }
1915           case 3:
1916           default:
1917             break;
1918         }
1919         p=GetCacheViewVirtualPixels(image_view,quadrant.x,quadrant.y,
1920           quadrant.width,quadrant.height,exception);
1921         if (p == (const Quantum *) NULL)
1922           break;
1923         for (j=0; j < (ssize_t) GetPixelChannels(gaussian_image); j++)
1924           mean[j]=0.0;
1925         k=p;
1926         for (n=0; n < (ssize_t) (width*width); n++)
1927         {
1928           for (j=0; j < (ssize_t) GetPixelChannels(gaussian_image); j++)
1929             mean[j]+=(double) k[j];
1930           k+=GetPixelChannels(gaussian_image);
1931         }
1932         for (j=0; j < (ssize_t) GetPixelChannels(gaussian_image); j++)
1933           mean[j]/=(double) (width*width);
1934         k=p;
1935         variance=0.0;
1936         for (n=0; n < (ssize_t) (width*width); n++)
1937         {
1938           double
1939             luma;
1940 
1941           luma=GetPixelLuma(gaussian_image,k);
1942           variance+=(luma-GetMeanLuma(gaussian_image,mean))*
1943             (luma-GetMeanLuma(gaussian_image,mean));
1944           k+=GetPixelChannels(gaussian_image);
1945         }
1946         if (variance < min_variance)
1947           {
1948             min_variance=variance;
1949             target=quadrant;
1950           }
1951       }
1952       if (i < 4)
1953         {
1954           status=MagickFalse;
1955           break;
1956         }
1957       status=InterpolatePixelChannels(gaussian_image,image_view,kuwahara_image,
1958         UndefinedInterpolatePixel,(double) target.x+target.width/2.0,(double)
1959         target.y+target.height/2.0,q,exception);
1960       if (status == MagickFalse)
1961         break;
1962       q+=GetPixelChannels(kuwahara_image);
1963     }
1964     if (SyncCacheViewAuthenticPixels(kuwahara_view,exception) == MagickFalse)
1965       status=MagickFalse;
1966     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1967       {
1968         MagickBooleanType
1969           proceed;
1970 
1971 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1972         #pragma omp atomic
1973 #endif
1974         progress++;
1975         proceed=SetImageProgress(image,KuwaharaImageTag,progress,image->rows);
1976         if (proceed == MagickFalse)
1977           status=MagickFalse;
1978       }
1979   }
1980   kuwahara_view=DestroyCacheView(kuwahara_view);
1981   image_view=DestroyCacheView(image_view);
1982   gaussian_image=DestroyImage(gaussian_image);
1983   if (status == MagickFalse)
1984     kuwahara_image=DestroyImage(kuwahara_image);
1985   return(kuwahara_image);
1986 }
1987 
1988 /*
1989 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1990 %                                                                             %
1991 %                                                                             %
1992 %                                                                             %
1993 %     L o c a l C o n t r a s t I m a g e                                     %
1994 %                                                                             %
1995 %                                                                             %
1996 %                                                                             %
1997 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1998 %
1999 %  LocalContrastImage() attempts to increase the appearance of large-scale
2000 %  light-dark transitions. Local contrast enhancement works similarly to
2001 %  sharpening with an unsharp mask, however the mask is instead created using
2002 %  an image with a greater blur distance.
2003 %
2004 %  The format of the LocalContrastImage method is:
2005 %
2006 %      Image *LocalContrastImage(const Image *image, const double radius,
2007 %        const double strength,ExceptionInfo *exception)
2008 %
2009 %  A description of each parameter follows:
2010 %
2011 %    o image: the image.
2012 %
2013 %    o radius: the radius of the Gaussian blur, in percentage with 100%
2014 %      resulting in a blur radius of 20% of largest dimension.
2015 %
2016 %    o strength: the strength of the blur mask in percentage.
2017 %
2018 %    o exception: return any errors or warnings in this structure.
2019 %
2020 */
LocalContrastImage(const Image * image,const double radius,const double strength,ExceptionInfo * exception)2021 MagickExport Image *LocalContrastImage(const Image *image,const double radius,
2022   const double strength,ExceptionInfo *exception)
2023 {
2024 #define LocalContrastImageTag  "LocalContrast/Image"
2025 
2026   CacheView
2027     *image_view,
2028     *contrast_view;
2029 
2030   float
2031     *interImage,
2032     *scanline,
2033     totalWeight;
2034 
2035   Image
2036     *contrast_image;
2037 
2038   MagickBooleanType
2039     status;
2040 
2041   MemoryInfo
2042     *scanline_info,
2043     *interImage_info;
2044 
2045   ssize_t
2046     scanLineSize,
2047     width;
2048 
2049   /*
2050     Initialize contrast image attributes.
2051   */
2052   assert(image != (const Image *) NULL);
2053   assert(image->signature == MagickCoreSignature);
2054   if (image->debug != MagickFalse)
2055     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2056   assert(exception != (ExceptionInfo *) NULL);
2057   assert(exception->signature == MagickCoreSignature);
2058 #if defined(MAGICKCORE_OPENCL_SUPPORT)
2059   contrast_image=AccelerateLocalContrastImage(image,radius,strength,exception);
2060   if (contrast_image != (Image *) NULL)
2061     return(contrast_image);
2062 #endif
2063   contrast_image=CloneImage(image,0,0,MagickTrue,exception);
2064   if (contrast_image == (Image *) NULL)
2065     return((Image *) NULL);
2066   if (SetImageStorageClass(contrast_image,DirectClass,exception) == MagickFalse)
2067     {
2068       contrast_image=DestroyImage(contrast_image);
2069       return((Image *) NULL);
2070     }
2071   image_view=AcquireVirtualCacheView(image,exception);
2072   contrast_view=AcquireAuthenticCacheView(contrast_image,exception);
2073   scanLineSize=(ssize_t) MagickMax(image->columns,image->rows);
2074   width=(ssize_t) scanLineSize*0.002f*fabs(radius);
2075   scanLineSize+=(2*width);
2076   scanline_info=AcquireVirtualMemory((size_t) GetOpenMPMaximumThreads()*
2077     scanLineSize,sizeof(*scanline));
2078   if (scanline_info == (MemoryInfo *) NULL)
2079     {
2080       contrast_view=DestroyCacheView(contrast_view);
2081       image_view=DestroyCacheView(image_view);
2082       contrast_image=DestroyImage(contrast_image);
2083       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2084     }
2085   scanline=(float *) GetVirtualMemoryBlob(scanline_info);
2086   /*
2087     Create intermediate buffer.
2088   */
2089   interImage_info=AcquireVirtualMemory(image->rows*(image->columns+(2*width)),
2090     sizeof(*interImage));
2091   if (interImage_info == (MemoryInfo *) NULL)
2092     {
2093       scanline_info=RelinquishVirtualMemory(scanline_info);
2094       contrast_view=DestroyCacheView(contrast_view);
2095       image_view=DestroyCacheView(image_view);
2096       contrast_image=DestroyImage(contrast_image);
2097       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2098     }
2099   interImage=(float *) GetVirtualMemoryBlob(interImage_info);
2100   totalWeight=(float) ((width+1)*(width+1));
2101   /*
2102     Vertical pass.
2103   */
2104   status=MagickTrue;
2105   {
2106     ssize_t
2107       x;
2108 
2109 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2110 #pragma omp parallel for schedule(static) \
2111     magick_number_threads(image,image,image->columns,1)
2112 #endif
2113     for (x=0; x < (ssize_t) image->columns; x++)
2114     {
2115       const int
2116         id = GetOpenMPThreadId();
2117 
2118       const Quantum
2119         *magick_restrict p;
2120 
2121       float
2122         *out,
2123         *pix,
2124         *pixels;
2125 
2126       ssize_t
2127         y;
2128 
2129       ssize_t
2130         i;
2131 
2132       if (status == MagickFalse)
2133         continue;
2134       pixels=scanline;
2135       pixels+=id*scanLineSize;
2136       pix=pixels;
2137       p=GetCacheViewVirtualPixels(image_view,x,-width,1,image->rows+(2*width),
2138         exception);
2139       if (p == (const Quantum *) NULL)
2140         {
2141           status=MagickFalse;
2142           continue;
2143         }
2144       for (y=0; y < (ssize_t) image->rows+(2*width); y++)
2145       {
2146         *pix++=(float)GetPixelLuma(image,p);
2147         p+=image->number_channels;
2148       }
2149       out=interImage+x+width;
2150       for (y=0; y < (ssize_t) image->rows; y++)
2151       {
2152         float
2153           sum,
2154           weight;
2155 
2156         weight=1.0f;
2157         sum=0;
2158         pix=pixels+y;
2159         for (i=0; i < width; i++)
2160         {
2161           sum+=weight*(*pix++);
2162           weight+=1.0f;
2163         }
2164         for (i=width+1; i < (2*width); i++)
2165         {
2166           sum+=weight*(*pix++);
2167           weight-=1.0f;
2168         }
2169         /* write to output */
2170         *out=sum/totalWeight;
2171         /* mirror into padding */
2172         if (x <= width && x != 0)
2173           *(out-(x*2))=*out;
2174         if ((x > (ssize_t) image->columns-width-2) &&
2175             (x != (ssize_t) image->columns-1))
2176           *(out+((image->columns-x-1)*2))=*out;
2177         out+=image->columns+(width*2);
2178       }
2179     }
2180   }
2181   /*
2182     Horizontal pass.
2183   */
2184   {
2185     ssize_t
2186       y;
2187 
2188 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2189 #pragma omp parallel for schedule(static) \
2190     magick_number_threads(image,image,image->rows,1)
2191 #endif
2192     for (y=0; y < (ssize_t) image->rows; y++)
2193     {
2194       const int
2195         id = GetOpenMPThreadId();
2196 
2197       const Quantum
2198         *magick_restrict p;
2199 
2200       float
2201         *pix,
2202         *pixels;
2203 
2204       Quantum
2205         *magick_restrict q;
2206 
2207       ssize_t
2208         x;
2209 
2210       ssize_t
2211         i;
2212 
2213       if (status == MagickFalse)
2214         continue;
2215       pixels=scanline;
2216       pixels+=id*scanLineSize;
2217       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2218       q=GetCacheViewAuthenticPixels(contrast_view,0,y,image->columns,1,
2219         exception);
2220       if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2221         {
2222           status=MagickFalse;
2223           continue;
2224         }
2225       memcpy(pixels,interImage+(y*(image->columns+(2*width))),(image->columns+
2226         (2*width))*sizeof(float));
2227       for (x=0; x < (ssize_t) image->columns; x++)
2228       {
2229         float
2230           mult,
2231           srcVal,
2232           sum,
2233           weight;
2234 
2235         PixelTrait
2236           traits;
2237 
2238         weight=1.0f;
2239         sum=0;
2240         pix=pixels+x;
2241         for (i=0; i < width; i++)
2242         {
2243           sum+=weight*(*pix++);
2244           weight+=1.0f;
2245         }
2246         for (i=width+1; i < (2*width); i++)
2247         {
2248           sum+=weight*(*pix++);
2249           weight-=1.0f;
2250         }
2251         /* Apply and write */
2252         srcVal=(float) GetPixelLuma(image,p);
2253         mult=(srcVal-(sum/totalWeight))*(strength/100.0f);
2254         mult=(srcVal+mult)/srcVal;
2255         traits=GetPixelChannelTraits(image,RedPixelChannel);
2256         if ((traits & UpdatePixelTrait) != 0)
2257           SetPixelRed(contrast_image,ClampToQuantum((MagickRealType)
2258             GetPixelRed(image,p)*mult),q);
2259         traits=GetPixelChannelTraits(image,GreenPixelChannel);
2260         if ((traits & UpdatePixelTrait) != 0)
2261           SetPixelGreen(contrast_image,ClampToQuantum((MagickRealType)
2262             GetPixelGreen(image,p)*mult),q);
2263         traits=GetPixelChannelTraits(image,BluePixelChannel);
2264         if ((traits & UpdatePixelTrait) != 0)
2265           SetPixelBlue(contrast_image,ClampToQuantum((MagickRealType)
2266             GetPixelBlue(image,p)*mult),q);
2267         p+=image->number_channels;
2268         q+=contrast_image->number_channels;
2269       }
2270       if (SyncCacheViewAuthenticPixels(contrast_view,exception) == MagickFalse)
2271         status=MagickFalse;
2272     }
2273   }
2274   scanline_info=RelinquishVirtualMemory(scanline_info);
2275   interImage_info=RelinquishVirtualMemory(interImage_info);
2276   contrast_view=DestroyCacheView(contrast_view);
2277   image_view=DestroyCacheView(image_view);
2278   if (status == MagickFalse)
2279     contrast_image=DestroyImage(contrast_image);
2280   return(contrast_image);
2281 }
2282 
2283 /*
2284 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2285 %                                                                             %
2286 %                                                                             %
2287 %                                                                             %
2288 %     M o t i o n B l u r I m a g e                                           %
2289 %                                                                             %
2290 %                                                                             %
2291 %                                                                             %
2292 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2293 %
2294 %  MotionBlurImage() simulates motion blur.  We convolve the image with a
2295 %  Gaussian operator of the given radius and standard deviation (sigma).
2296 %  For reasonable results, radius should be larger than sigma.  Use a
2297 %  radius of 0 and MotionBlurImage() selects a suitable radius for you.
2298 %  Angle gives the angle of the blurring motion.
2299 %
2300 %  Andrew Protano contributed this effect.
2301 %
2302 %  The format of the MotionBlurImage method is:
2303 %
2304 %    Image *MotionBlurImage(const Image *image,const double radius,
2305 %      const double sigma,const double angle,ExceptionInfo *exception)
2306 %
2307 %  A description of each parameter follows:
2308 %
2309 %    o image: the image.
2310 %
2311 %    o radius: the radius of the Gaussian, in pixels, not counting
2312 %      the center pixel.
2313 %
2314 %    o sigma: the standard deviation of the Gaussian, in pixels.
2315 %
2316 %    o angle: Apply the effect along this angle.
2317 %
2318 %    o exception: return any errors or warnings in this structure.
2319 %
2320 */
2321 
GetMotionBlurKernel(const size_t width,const double sigma)2322 static MagickRealType *GetMotionBlurKernel(const size_t width,
2323   const double sigma)
2324 {
2325   MagickRealType
2326     *kernel,
2327     normalize;
2328 
2329   ssize_t
2330     i;
2331 
2332   /*
2333    Generate a 1-D convolution kernel.
2334   */
2335   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
2336   kernel=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory((size_t)
2337     width,sizeof(*kernel)));
2338   if (kernel == (MagickRealType *) NULL)
2339     return(kernel);
2340   normalize=0.0;
2341   for (i=0; i < (ssize_t) width; i++)
2342   {
2343     kernel[i]=(MagickRealType) (exp((-((double) i*i)/(double) (2.0*MagickSigma*
2344       MagickSigma)))/(MagickSQ2PI*MagickSigma));
2345     normalize+=kernel[i];
2346   }
2347   for (i=0; i < (ssize_t) width; i++)
2348     kernel[i]/=normalize;
2349   return(kernel);
2350 }
2351 
MotionBlurImage(const Image * image,const double radius,const double sigma,const double angle,ExceptionInfo * exception)2352 MagickExport Image *MotionBlurImage(const Image *image,const double radius,
2353   const double sigma,const double angle,ExceptionInfo *exception)
2354 {
2355 #define BlurImageTag  "Blur/Image"
2356 
2357   CacheView
2358     *blur_view,
2359     *image_view,
2360     *motion_view;
2361 
2362   Image
2363     *blur_image;
2364 
2365   MagickBooleanType
2366     status;
2367 
2368   MagickOffsetType
2369     progress;
2370 
2371   MagickRealType
2372     *kernel;
2373 
2374   OffsetInfo
2375     *offset;
2376 
2377   PointInfo
2378     point;
2379 
2380   size_t
2381     width;
2382 
2383   ssize_t
2384     w,
2385     y;
2386 
2387   assert(image != (Image *) NULL);
2388   assert(image->signature == MagickCoreSignature);
2389   if (image->debug != MagickFalse)
2390     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2391   assert(exception != (ExceptionInfo *) NULL);
2392   width=GetOptimalKernelWidth1D(radius,sigma);
2393   kernel=GetMotionBlurKernel(width,sigma);
2394   if (kernel == (MagickRealType *) NULL)
2395     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2396   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
2397   if (offset == (OffsetInfo *) NULL)
2398     {
2399       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2400       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2401     }
2402   point.x=(double) width*sin(DegreesToRadians(angle));
2403   point.y=(double) width*cos(DegreesToRadians(angle));
2404   for (w=0; w < (ssize_t) width; w++)
2405   {
2406     offset[w].x=CastDoubleToLong(ceil((double) (w*point.y)/
2407       hypot(point.x,point.y)-0.5));
2408     offset[w].y=CastDoubleToLong(ceil((double) (w*point.x)/
2409       hypot(point.x,point.y)-0.5));
2410   }
2411   /*
2412     Motion blur image.
2413   */
2414 #if defined(MAGICKCORE_OPENCL_SUPPORT)
2415   blur_image=AccelerateMotionBlurImage(image,kernel,width,offset,exception);
2416   if (blur_image != (Image *) NULL)
2417     {
2418       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2419       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2420       return(blur_image);
2421     }
2422 #endif
2423   blur_image=CloneImage(image,0,0,MagickTrue,exception);
2424   if (blur_image == (Image *) NULL)
2425     {
2426       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2427       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2428       return((Image *) NULL);
2429     }
2430   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
2431     {
2432       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2433       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2434       blur_image=DestroyImage(blur_image);
2435       return((Image *) NULL);
2436     }
2437   status=MagickTrue;
2438   progress=0;
2439   image_view=AcquireVirtualCacheView(image,exception);
2440   motion_view=AcquireVirtualCacheView(image,exception);
2441   blur_view=AcquireAuthenticCacheView(blur_image,exception);
2442 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2443   #pragma omp parallel for schedule(static) shared(progress,status) \
2444     magick_number_threads(image,blur_image,image->rows,1)
2445 #endif
2446   for (y=0; y < (ssize_t) image->rows; y++)
2447   {
2448     const Quantum
2449       *magick_restrict p;
2450 
2451     Quantum
2452       *magick_restrict q;
2453 
2454     ssize_t
2455       x;
2456 
2457     if (status == MagickFalse)
2458       continue;
2459     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2460     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
2461       exception);
2462     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2463       {
2464         status=MagickFalse;
2465         continue;
2466       }
2467     for (x=0; x < (ssize_t) image->columns; x++)
2468     {
2469       ssize_t
2470         i;
2471 
2472       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2473       {
2474         double
2475           alpha,
2476           gamma,
2477           pixel;
2478 
2479         PixelChannel
2480           channel;
2481 
2482         PixelTrait
2483           blur_traits,
2484           traits;
2485 
2486         const Quantum
2487           *magick_restrict r;
2488 
2489         MagickRealType
2490           *magick_restrict k;
2491 
2492         ssize_t
2493           j;
2494 
2495         channel=GetPixelChannelChannel(image,i);
2496         traits=GetPixelChannelTraits(image,channel);
2497         blur_traits=GetPixelChannelTraits(blur_image,channel);
2498         if ((traits == UndefinedPixelTrait) ||
2499             (blur_traits == UndefinedPixelTrait))
2500           continue;
2501         if ((blur_traits & CopyPixelTrait) != 0)
2502           {
2503             SetPixelChannel(blur_image,channel,p[i],q);
2504             continue;
2505           }
2506         k=kernel;
2507         pixel=0.0;
2508         if ((blur_traits & BlendPixelTrait) == 0)
2509           {
2510             for (j=0; j < (ssize_t) width; j++)
2511             {
2512               r=GetCacheViewVirtualPixels(motion_view,x+offset[j].x,y+
2513                 offset[j].y,1,1,exception);
2514               if (r == (const Quantum *) NULL)
2515                 {
2516                   status=MagickFalse;
2517                   continue;
2518                 }
2519               pixel+=(*k)*r[i];
2520               k++;
2521             }
2522             SetPixelChannel(blur_image,channel,ClampToQuantum(pixel),q);
2523             continue;
2524           }
2525         alpha=0.0;
2526         gamma=0.0;
2527         for (j=0; j < (ssize_t) width; j++)
2528         {
2529           r=GetCacheViewVirtualPixels(motion_view,x+offset[j].x,y+offset[j].y,1,
2530             1,exception);
2531           if (r == (const Quantum *) NULL)
2532             {
2533               status=MagickFalse;
2534               continue;
2535             }
2536           alpha=(double) (QuantumScale*GetPixelAlpha(image,r));
2537           pixel+=(*k)*alpha*r[i];
2538           gamma+=(*k)*alpha;
2539           k++;
2540         }
2541         gamma=PerceptibleReciprocal(gamma);
2542         SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
2543       }
2544       p+=GetPixelChannels(image);
2545       q+=GetPixelChannels(blur_image);
2546     }
2547     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
2548       status=MagickFalse;
2549     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2550       {
2551         MagickBooleanType
2552           proceed;
2553 
2554 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2555         #pragma omp atomic
2556 #endif
2557         progress++;
2558         proceed=SetImageProgress(image,BlurImageTag,progress,image->rows);
2559         if (proceed == MagickFalse)
2560           status=MagickFalse;
2561       }
2562   }
2563   blur_view=DestroyCacheView(blur_view);
2564   motion_view=DestroyCacheView(motion_view);
2565   image_view=DestroyCacheView(image_view);
2566   kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2567   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
2568   if (status == MagickFalse)
2569     blur_image=DestroyImage(blur_image);
2570   return(blur_image);
2571 }
2572 
2573 /*
2574 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2575 %                                                                             %
2576 %                                                                             %
2577 %                                                                             %
2578 %     P r e v i e w I m a g e                                                 %
2579 %                                                                             %
2580 %                                                                             %
2581 %                                                                             %
2582 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2583 %
2584 %  PreviewImage() tiles 9 thumbnails of the specified image with an image
2585 %  processing operation applied with varying parameters.  This may be helpful
2586 %  pin-pointing an appropriate parameter for a particular image processing
2587 %  operation.
2588 %
2589 %  The format of the PreviewImages method is:
2590 %
2591 %      Image *PreviewImages(const Image *image,const PreviewType preview,
2592 %        ExceptionInfo *exception)
2593 %
2594 %  A description of each parameter follows:
2595 %
2596 %    o image: the image.
2597 %
2598 %    o preview: the image processing operation.
2599 %
2600 %    o exception: return any errors or warnings in this structure.
2601 %
2602 */
PreviewImage(const Image * image,const PreviewType preview,ExceptionInfo * exception)2603 MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
2604   ExceptionInfo *exception)
2605 {
2606 #define NumberTiles  9
2607 #define PreviewImageTag  "Preview/Image"
2608 #define DefaultPreviewGeometry  "204x204+10+10"
2609 
2610   char
2611     factor[MagickPathExtent],
2612     label[MagickPathExtent];
2613 
2614   double
2615     degrees,
2616     gamma,
2617     percentage,
2618     radius,
2619     sigma,
2620     threshold;
2621 
2622   Image
2623     *images,
2624     *montage_image,
2625     *preview_image,
2626     *thumbnail;
2627 
2628   ImageInfo
2629     *preview_info;
2630 
2631   MagickBooleanType
2632     proceed;
2633 
2634   MontageInfo
2635     *montage_info;
2636 
2637   QuantizeInfo
2638     quantize_info;
2639 
2640   RectangleInfo
2641     geometry;
2642 
2643   ssize_t
2644     i,
2645     x;
2646 
2647   size_t
2648     colors;
2649 
2650   ssize_t
2651     y;
2652 
2653   /*
2654     Open output image file.
2655   */
2656   assert(image != (Image *) NULL);
2657   assert(image->signature == MagickCoreSignature);
2658   if (image->debug != MagickFalse)
2659     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2660   colors=2;
2661   degrees=0.0;
2662   gamma=(-0.2f);
2663   preview_info=AcquireImageInfo();
2664   SetGeometry(image,&geometry);
2665   (void) ParseMetaGeometry(DefaultPreviewGeometry,&geometry.x,&geometry.y,
2666     &geometry.width,&geometry.height);
2667   images=NewImageList();
2668   percentage=12.5;
2669   GetQuantizeInfo(&quantize_info);
2670   radius=0.0;
2671   sigma=1.0;
2672   threshold=0.0;
2673   x=0;
2674   y=0;
2675   for (i=0; i < NumberTiles; i++)
2676   {
2677     thumbnail=ThumbnailImage(image,geometry.width,geometry.height,exception);
2678     if (thumbnail == (Image *) NULL)
2679       break;
2680     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
2681       (void *) NULL);
2682     (void) SetImageProperty(thumbnail,"label",DefaultTileLabel,exception);
2683     if (i == (NumberTiles/2))
2684       {
2685         (void) QueryColorCompliance("#dfdfdf",AllCompliance,
2686           &thumbnail->matte_color,exception);
2687         AppendImageToList(&images,thumbnail);
2688         continue;
2689       }
2690     switch (preview)
2691     {
2692       case RotatePreview:
2693       {
2694         degrees+=45.0;
2695         preview_image=RotateImage(thumbnail,degrees,exception);
2696         (void) FormatLocaleString(label,MagickPathExtent,"rotate %g",degrees);
2697         break;
2698       }
2699       case ShearPreview:
2700       {
2701         degrees+=5.0;
2702         preview_image=ShearImage(thumbnail,degrees,degrees,exception);
2703         (void) FormatLocaleString(label,MagickPathExtent,"shear %gx%g",degrees,
2704           2.0*degrees);
2705         break;
2706       }
2707       case RollPreview:
2708       {
2709         x=(ssize_t) ((i+1)*thumbnail->columns)/NumberTiles;
2710         y=(ssize_t) ((i+1)*thumbnail->rows)/NumberTiles;
2711         preview_image=RollImage(thumbnail,x,y,exception);
2712         (void) FormatLocaleString(label,MagickPathExtent,"roll %+.20gx%+.20g",
2713           (double) x,(double) y);
2714         break;
2715       }
2716       case HuePreview:
2717       {
2718         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2719         if (preview_image == (Image *) NULL)
2720           break;
2721         (void) FormatLocaleString(factor,MagickPathExtent,"100,100,%g",2.0*
2722           percentage);
2723         (void) ModulateImage(preview_image,factor,exception);
2724         (void) FormatLocaleString(label,MagickPathExtent,"modulate %s",factor);
2725         break;
2726       }
2727       case SaturationPreview:
2728       {
2729         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2730         if (preview_image == (Image *) NULL)
2731           break;
2732         (void) FormatLocaleString(factor,MagickPathExtent,"100,%g",2.0*
2733           percentage);
2734         (void) ModulateImage(preview_image,factor,exception);
2735         (void) FormatLocaleString(label,MagickPathExtent,"modulate %s",factor);
2736         break;
2737       }
2738       case BrightnessPreview:
2739       {
2740         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2741         if (preview_image == (Image *) NULL)
2742           break;
2743         (void) FormatLocaleString(factor,MagickPathExtent,"%g",2.0*percentage);
2744         (void) ModulateImage(preview_image,factor,exception);
2745         (void) FormatLocaleString(label,MagickPathExtent,"modulate %s",factor);
2746         break;
2747       }
2748       case GammaPreview:
2749       default:
2750       {
2751         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2752         if (preview_image == (Image *) NULL)
2753           break;
2754         gamma+=0.4f;
2755         (void) GammaImage(preview_image,gamma,exception);
2756         (void) FormatLocaleString(label,MagickPathExtent,"gamma %g",gamma);
2757         break;
2758       }
2759       case SpiffPreview:
2760       {
2761         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2762         if (preview_image != (Image *) NULL)
2763           for (x=0; x < i; x++)
2764             (void) ContrastImage(preview_image,MagickTrue,exception);
2765         (void) FormatLocaleString(label,MagickPathExtent,"contrast (%.20g)",
2766           (double) i+1);
2767         break;
2768       }
2769       case DullPreview:
2770       {
2771         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2772         if (preview_image == (Image *) NULL)
2773           break;
2774         for (x=0; x < i; x++)
2775           (void) ContrastImage(preview_image,MagickFalse,exception);
2776         (void) FormatLocaleString(label,MagickPathExtent,"+contrast (%.20g)",
2777           (double) i+1);
2778         break;
2779       }
2780       case GrayscalePreview:
2781       {
2782         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2783         if (preview_image == (Image *) NULL)
2784           break;
2785         colors<<=1;
2786         quantize_info.number_colors=colors;
2787         quantize_info.colorspace=GRAYColorspace;
2788         (void) QuantizeImage(&quantize_info,preview_image,exception);
2789         (void) FormatLocaleString(label,MagickPathExtent,
2790           "-colorspace gray -colors %.20g",(double) colors);
2791         break;
2792       }
2793       case QuantizePreview:
2794       {
2795         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2796         if (preview_image == (Image *) NULL)
2797           break;
2798         colors<<=1;
2799         quantize_info.number_colors=colors;
2800         (void) QuantizeImage(&quantize_info,preview_image,exception);
2801         (void) FormatLocaleString(label,MagickPathExtent,"colors %.20g",
2802           (double) colors);
2803         break;
2804       }
2805       case DespecklePreview:
2806       {
2807         for (x=0; x < (i-1); x++)
2808         {
2809           preview_image=DespeckleImage(thumbnail,exception);
2810           if (preview_image == (Image *) NULL)
2811             break;
2812           thumbnail=DestroyImage(thumbnail);
2813           thumbnail=preview_image;
2814         }
2815         preview_image=DespeckleImage(thumbnail,exception);
2816         if (preview_image == (Image *) NULL)
2817           break;
2818         (void) FormatLocaleString(label,MagickPathExtent,"despeckle (%.20g)",
2819           (double) i+1);
2820         break;
2821       }
2822       case ReduceNoisePreview:
2823       {
2824         preview_image=StatisticImage(thumbnail,NonpeakStatistic,(size_t)
2825           radius,(size_t) radius,exception);
2826         (void) FormatLocaleString(label,MagickPathExtent,"noise %g",radius);
2827         break;
2828       }
2829       case AddNoisePreview:
2830       {
2831         switch ((int) i)
2832         {
2833           case 0:
2834           {
2835             (void) CopyMagickString(factor,"uniform",MagickPathExtent);
2836             break;
2837           }
2838           case 1:
2839           {
2840             (void) CopyMagickString(factor,"gaussian",MagickPathExtent);
2841             break;
2842           }
2843           case 2:
2844           {
2845             (void) CopyMagickString(factor,"multiplicative",MagickPathExtent);
2846             break;
2847           }
2848           case 3:
2849           {
2850             (void) CopyMagickString(factor,"impulse",MagickPathExtent);
2851             break;
2852           }
2853           case 5:
2854           {
2855             (void) CopyMagickString(factor,"laplacian",MagickPathExtent);
2856             break;
2857           }
2858           case 6:
2859           {
2860             (void) CopyMagickString(factor,"Poisson",MagickPathExtent);
2861             break;
2862           }
2863           default:
2864           {
2865             (void) CopyMagickString(thumbnail->magick,"NULL",MagickPathExtent);
2866             break;
2867           }
2868         }
2869         preview_image=StatisticImage(thumbnail,NonpeakStatistic,(size_t) i,
2870           (size_t) i,exception);
2871         (void) FormatLocaleString(label,MagickPathExtent,"+noise %s",factor);
2872         break;
2873       }
2874       case SharpenPreview:
2875       {
2876         preview_image=SharpenImage(thumbnail,radius,sigma,exception);
2877         (void) FormatLocaleString(label,MagickPathExtent,"sharpen %gx%g",
2878           radius,sigma);
2879         break;
2880       }
2881       case BlurPreview:
2882       {
2883         preview_image=BlurImage(thumbnail,radius,sigma,exception);
2884         (void) FormatLocaleString(label,MagickPathExtent,"blur %gx%g",radius,
2885           sigma);
2886         break;
2887       }
2888       case ThresholdPreview:
2889       {
2890         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2891         if (preview_image == (Image *) NULL)
2892           break;
2893         (void) BilevelImage(thumbnail,(double) (percentage*((double)
2894           QuantumRange+1.0))/100.0,exception);
2895         (void) FormatLocaleString(label,MagickPathExtent,"threshold %g",
2896           (double) (percentage*((double) QuantumRange+1.0))/100.0);
2897         break;
2898       }
2899       case EdgeDetectPreview:
2900       {
2901         preview_image=EdgeImage(thumbnail,radius,exception);
2902         (void) FormatLocaleString(label,MagickPathExtent,"edge %g",radius);
2903         break;
2904       }
2905       case SpreadPreview:
2906       {
2907         preview_image=SpreadImage(thumbnail,image->interpolate,radius,
2908           exception);
2909         (void) FormatLocaleString(label,MagickPathExtent,"spread %g",
2910           radius+0.5);
2911         break;
2912       }
2913       case SolarizePreview:
2914       {
2915         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2916         if (preview_image == (Image *) NULL)
2917           break;
2918         (void) SolarizeImage(preview_image,(double) QuantumRange*percentage/
2919           100.0,exception);
2920         (void) FormatLocaleString(label,MagickPathExtent,"solarize %g",
2921           (QuantumRange*percentage)/100.0);
2922         break;
2923       }
2924       case ShadePreview:
2925       {
2926         degrees+=10.0;
2927         preview_image=ShadeImage(thumbnail,MagickTrue,degrees,degrees,
2928           exception);
2929         (void) FormatLocaleString(label,MagickPathExtent,"shade %gx%g",degrees,
2930           degrees);
2931         break;
2932       }
2933       case RaisePreview:
2934       {
2935         RectangleInfo
2936           raise;
2937 
2938         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2939         if (preview_image == (Image *) NULL)
2940           break;
2941         raise.width=(size_t) (2*i+2);
2942         raise.height=(size_t) (2*i+2);
2943         raise.x=(i-1)/2;
2944         raise.y=(i-1)/2;
2945         (void) RaiseImage(preview_image,&raise,MagickTrue,exception);
2946         (void) FormatLocaleString(label,MagickPathExtent,
2947           "raise %.20gx%.20g%+.20g%+.20g",(double) raise.width,(double)
2948           raise.height,(double) raise.x,(double) raise.y);
2949         break;
2950       }
2951       case SegmentPreview:
2952       {
2953         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2954         if (preview_image == (Image *) NULL)
2955           break;
2956         threshold+=0.4f;
2957         (void) SegmentImage(preview_image,sRGBColorspace,MagickFalse,threshold,
2958           threshold,exception);
2959         (void) FormatLocaleString(label,MagickPathExtent,"segment %gx%g",
2960           threshold,threshold);
2961         break;
2962       }
2963       case SwirlPreview:
2964       {
2965         preview_image=SwirlImage(thumbnail,degrees,image->interpolate,
2966           exception);
2967         (void) FormatLocaleString(label,MagickPathExtent,"swirl %g",degrees);
2968         degrees+=45.0;
2969         break;
2970       }
2971       case ImplodePreview:
2972       {
2973         degrees+=0.1f;
2974         preview_image=ImplodeImage(thumbnail,degrees,image->interpolate,
2975           exception);
2976         (void) FormatLocaleString(label,MagickPathExtent,"implode %g",degrees);
2977         break;
2978       }
2979       case WavePreview:
2980       {
2981         degrees+=5.0f;
2982         preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,
2983           image->interpolate,exception);
2984         (void) FormatLocaleString(label,MagickPathExtent,"wave %gx%g",0.5*
2985           degrees,2.0*degrees);
2986         break;
2987       }
2988       case OilPaintPreview:
2989       {
2990         preview_image=OilPaintImage(thumbnail,(double) radius,(double) sigma,
2991           exception);
2992         (void) FormatLocaleString(label,MagickPathExtent,"charcoal %gx%g",
2993           radius,sigma);
2994         break;
2995       }
2996       case CharcoalDrawingPreview:
2997       {
2998         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
2999           exception);
3000         (void) FormatLocaleString(label,MagickPathExtent,"charcoal %gx%g",
3001           radius,sigma);
3002         break;
3003       }
3004       case JPEGPreview:
3005       {
3006         char
3007           filename[MagickPathExtent];
3008 
3009         int
3010           file;
3011 
3012         MagickBooleanType
3013           status;
3014 
3015         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
3016         if (preview_image == (Image *) NULL)
3017           break;
3018         preview_info->quality=(size_t) percentage;
3019         (void) FormatLocaleString(factor,MagickPathExtent,"%.20g",(double)
3020           preview_info->quality);
3021         file=AcquireUniqueFileResource(filename);
3022         if (file != -1)
3023           file=close(file)-1;
3024         (void) FormatLocaleString(preview_image->filename,MagickPathExtent,
3025           "jpeg:%s",filename);
3026         status=WriteImage(preview_info,preview_image,exception);
3027         if (status != MagickFalse)
3028           {
3029             Image
3030               *quality_image;
3031 
3032             (void) CopyMagickString(preview_info->filename,
3033               preview_image->filename,MagickPathExtent);
3034             quality_image=ReadImage(preview_info,exception);
3035             if (quality_image != (Image *) NULL)
3036               {
3037                 preview_image=DestroyImage(preview_image);
3038                 preview_image=quality_image;
3039               }
3040           }
3041         (void) RelinquishUniqueFileResource(preview_image->filename);
3042         if ((GetBlobSize(preview_image)/1024) >= 1024)
3043           (void) FormatLocaleString(label,MagickPathExtent,"quality %s\n%gmb ",
3044             factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
3045             1024.0/1024.0);
3046         else
3047           if (GetBlobSize(preview_image) >= 1024)
3048             (void) FormatLocaleString(label,MagickPathExtent,
3049               "quality %s\n%gkb ",factor,(double) ((MagickOffsetType)
3050               GetBlobSize(preview_image))/1024.0);
3051           else
3052             (void) FormatLocaleString(label,MagickPathExtent,
3053               "quality %s\n%.20gb ",factor,(double) ((MagickOffsetType)
3054               GetBlobSize(thumbnail)));
3055         break;
3056       }
3057     }
3058     thumbnail=DestroyImage(thumbnail);
3059     percentage+=12.5;
3060     radius+=0.5;
3061     sigma+=0.25;
3062     if (preview_image == (Image *) NULL)
3063       break;
3064     preview_image->alpha_trait=UndefinedPixelTrait;
3065     (void) DeleteImageProperty(preview_image,"label");
3066     (void) SetImageProperty(preview_image,"label",label,exception);
3067     AppendImageToList(&images,preview_image);
3068     proceed=SetImageProgress(image,PreviewImageTag,(MagickOffsetType) i,
3069       NumberTiles);
3070     if (proceed == MagickFalse)
3071       break;
3072   }
3073   if (images == (Image *) NULL)
3074     {
3075       preview_info=DestroyImageInfo(preview_info);
3076       return((Image *) NULL);
3077     }
3078   /*
3079     Create the montage.
3080   */
3081   montage_info=CloneMontageInfo(preview_info,(MontageInfo *) NULL);
3082   (void) CopyMagickString(montage_info->filename,image->filename,
3083     MagickPathExtent);
3084   montage_info->shadow=MagickTrue;
3085   (void) CloneString(&montage_info->tile,"3x3");
3086   (void) CloneString(&montage_info->geometry,DefaultPreviewGeometry);
3087   (void) CloneString(&montage_info->frame,DefaultTileFrame);
3088   montage_image=MontageImages(images,montage_info,exception);
3089   montage_info=DestroyMontageInfo(montage_info);
3090   images=DestroyImageList(images);
3091   if (montage_image == (Image *) NULL)
3092     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3093   if (montage_image->montage != (char *) NULL)
3094     {
3095       /*
3096         Free image directory.
3097       */
3098       montage_image->montage=(char *) RelinquishMagickMemory(
3099         montage_image->montage);
3100       if (image->directory != (char *) NULL)
3101         montage_image->directory=(char *) RelinquishMagickMemory(
3102           montage_image->directory);
3103     }
3104   preview_info=DestroyImageInfo(preview_info);
3105   return(montage_image);
3106 }
3107 
3108 /*
3109 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3110 %                                                                             %
3111 %                                                                             %
3112 %                                                                             %
3113 %     R o t a t i o n a l B l u r I m a g e                                   %
3114 %                                                                             %
3115 %                                                                             %
3116 %                                                                             %
3117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3118 %
3119 %  RotationalBlurImage() applies a radial blur to the image.
3120 %
3121 %  Andrew Protano contributed this effect.
3122 %
3123 %  The format of the RotationalBlurImage method is:
3124 %
3125 %    Image *RotationalBlurImage(const Image *image,const double angle,
3126 %      ExceptionInfo *exception)
3127 %
3128 %  A description of each parameter follows:
3129 %
3130 %    o image: the image.
3131 %
3132 %    o angle: the angle of the radial blur.
3133 %
3134 %    o blur: the blur.
3135 %
3136 %    o exception: return any errors or warnings in this structure.
3137 %
3138 */
RotationalBlurImage(const Image * image,const double angle,ExceptionInfo * exception)3139 MagickExport Image *RotationalBlurImage(const Image *image,const double angle,
3140   ExceptionInfo *exception)
3141 {
3142   CacheView
3143     *blur_view,
3144     *image_view,
3145     *radial_view;
3146 
3147   double
3148     blur_radius,
3149     *cos_theta,
3150     offset,
3151     *sin_theta,
3152     theta;
3153 
3154   Image
3155     *blur_image;
3156 
3157   MagickBooleanType
3158     status;
3159 
3160   MagickOffsetType
3161     progress;
3162 
3163   PointInfo
3164     blur_center;
3165 
3166   size_t
3167     n;
3168 
3169   ssize_t
3170     w,
3171     y;
3172 
3173   /*
3174     Allocate blur image.
3175   */
3176   assert(image != (Image *) NULL);
3177   assert(image->signature == MagickCoreSignature);
3178   if (image->debug != MagickFalse)
3179     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3180   assert(exception != (ExceptionInfo *) NULL);
3181   assert(exception->signature == MagickCoreSignature);
3182 #if defined(MAGICKCORE_OPENCL_SUPPORT)
3183   blur_image=AccelerateRotationalBlurImage(image,angle,exception);
3184   if (blur_image != (Image *) NULL)
3185     return(blur_image);
3186 #endif
3187   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3188   if (blur_image == (Image *) NULL)
3189     return((Image *) NULL);
3190   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
3191     {
3192       blur_image=DestroyImage(blur_image);
3193       return((Image *) NULL);
3194     }
3195   blur_center.x=(double) (image->columns-1)/2.0;
3196   blur_center.y=(double) (image->rows-1)/2.0;
3197   blur_radius=hypot(blur_center.x,blur_center.y);
3198   n=(size_t) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+2UL);
3199   theta=DegreesToRadians(angle)/(double) (n-1);
3200   cos_theta=(double *) AcquireQuantumMemory((size_t) n,sizeof(*cos_theta));
3201   sin_theta=(double *) AcquireQuantumMemory((size_t) n,sizeof(*sin_theta));
3202   if ((cos_theta == (double *) NULL) || (sin_theta == (double *) NULL))
3203     {
3204       if (cos_theta != (double *) NULL)
3205         cos_theta=(double *) RelinquishMagickMemory(cos_theta);
3206       if (sin_theta != (double *) NULL)
3207         sin_theta=(double *) RelinquishMagickMemory(sin_theta);
3208       blur_image=DestroyImage(blur_image);
3209       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3210     }
3211   offset=theta*(double) (n-1)/2.0;
3212   for (w=0; w < (ssize_t) n; w++)
3213   {
3214     cos_theta[w]=cos((double) (theta*w-offset));
3215     sin_theta[w]=sin((double) (theta*w-offset));
3216   }
3217   /*
3218     Radial blur image.
3219   */
3220   status=MagickTrue;
3221   progress=0;
3222   image_view=AcquireVirtualCacheView(image,exception);
3223   radial_view=AcquireVirtualCacheView(image,exception);
3224   blur_view=AcquireAuthenticCacheView(blur_image,exception);
3225 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3226   #pragma omp parallel for schedule(static) shared(progress,status) \
3227     magick_number_threads(image,blur_image,image->rows,1)
3228 #endif
3229   for (y=0; y < (ssize_t) image->rows; y++)
3230   {
3231     const Quantum
3232       *magick_restrict p;
3233 
3234     Quantum
3235       *magick_restrict q;
3236 
3237     ssize_t
3238       x;
3239 
3240     if (status == MagickFalse)
3241       continue;
3242     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3243     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
3244       exception);
3245     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3246       {
3247         status=MagickFalse;
3248         continue;
3249       }
3250     for (x=0; x < (ssize_t) image->columns; x++)
3251     {
3252       double
3253         radius;
3254 
3255       PointInfo
3256         center;
3257 
3258       ssize_t
3259         i;
3260 
3261       size_t
3262         step;
3263 
3264       center.x=(double) x-blur_center.x;
3265       center.y=(double) y-blur_center.y;
3266       radius=hypot((double) center.x,center.y);
3267       if (radius == 0)
3268         step=1;
3269       else
3270         {
3271           step=(size_t) (blur_radius/radius);
3272           if (step == 0)
3273             step=1;
3274           else
3275             if (step >= n)
3276               step=n-1;
3277         }
3278       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3279       {
3280         double
3281           gamma,
3282           pixel;
3283 
3284         PixelChannel
3285           channel;
3286 
3287         PixelTrait
3288           blur_traits,
3289           traits;
3290 
3291         const Quantum
3292           *magick_restrict r;
3293 
3294         ssize_t
3295           j;
3296 
3297         channel=GetPixelChannelChannel(image,i);
3298         traits=GetPixelChannelTraits(image,channel);
3299         blur_traits=GetPixelChannelTraits(blur_image,channel);
3300         if ((traits == UndefinedPixelTrait) ||
3301             (blur_traits == UndefinedPixelTrait))
3302           continue;
3303         if ((blur_traits & CopyPixelTrait) != 0)
3304           {
3305             SetPixelChannel(blur_image,channel,p[i],q);
3306             continue;
3307           }
3308         gamma=0.0;
3309         pixel=0.0;
3310         if ((GetPixelChannelTraits(image,AlphaPixelChannel) == UndefinedPixelTrait) ||
3311             (channel == AlphaPixelChannel))
3312           {
3313             for (j=0; j < (ssize_t) n; j+=(ssize_t) step)
3314             {
3315               r=GetCacheViewVirtualPixels(radial_view, (ssize_t) (blur_center.x+
3316                 center.x*cos_theta[j]-center.y*sin_theta[j]+0.5),(ssize_t)
3317                 (blur_center.y+center.x*sin_theta[j]+center.y*cos_theta[j]+0.5),
3318                 1,1,exception);
3319               if (r == (const Quantum *) NULL)
3320                 {
3321                   status=MagickFalse;
3322                   continue;
3323                 }
3324               pixel+=r[i];
3325               gamma++;
3326             }
3327             gamma=PerceptibleReciprocal(gamma);
3328             SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
3329             continue;
3330           }
3331         for (j=0; j < (ssize_t) n; j+=(ssize_t) step)
3332         {
3333           double
3334             alpha;
3335 
3336           r=GetCacheViewVirtualPixels(radial_view, (ssize_t) (blur_center.x+
3337             center.x*cos_theta[j]-center.y*sin_theta[j]+0.5),(ssize_t)
3338             (blur_center.y+center.x*sin_theta[j]+center.y*cos_theta[j]+0.5),
3339             1,1,exception);
3340           if (r == (const Quantum *) NULL)
3341             {
3342               status=MagickFalse;
3343               continue;
3344             }
3345           alpha=(double) QuantumScale*GetPixelAlpha(image,r);
3346           pixel+=alpha*r[i];
3347           gamma+=alpha;
3348         }
3349         gamma=PerceptibleReciprocal(gamma);
3350         SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
3351       }
3352       p+=GetPixelChannels(image);
3353       q+=GetPixelChannels(blur_image);
3354     }
3355     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
3356       status=MagickFalse;
3357     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3358       {
3359         MagickBooleanType
3360           proceed;
3361 
3362 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3363         #pragma omp atomic
3364 #endif
3365         progress++;
3366         proceed=SetImageProgress(image,BlurImageTag,progress,image->rows);
3367         if (proceed == MagickFalse)
3368           status=MagickFalse;
3369       }
3370   }
3371   blur_view=DestroyCacheView(blur_view);
3372   radial_view=DestroyCacheView(radial_view);
3373   image_view=DestroyCacheView(image_view);
3374   cos_theta=(double *) RelinquishMagickMemory(cos_theta);
3375   sin_theta=(double *) RelinquishMagickMemory(sin_theta);
3376   if (status == MagickFalse)
3377     blur_image=DestroyImage(blur_image);
3378   return(blur_image);
3379 }
3380 
3381 /*
3382 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3383 %                                                                             %
3384 %                                                                             %
3385 %                                                                             %
3386 %     S e l e c t i v e B l u r I m a g e                                     %
3387 %                                                                             %
3388 %                                                                             %
3389 %                                                                             %
3390 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3391 %
3392 %  SelectiveBlurImage() selectively blur pixels within a contrast threshold.
3393 %  It is similar to the unsharpen mask that sharpens everything with contrast
3394 %  above a certain threshold.
3395 %
3396 %  The format of the SelectiveBlurImage method is:
3397 %
3398 %      Image *SelectiveBlurImage(const Image *image,const double radius,
3399 %        const double sigma,const double threshold,ExceptionInfo *exception)
3400 %
3401 %  A description of each parameter follows:
3402 %
3403 %    o image: the image.
3404 %
3405 %    o radius: the radius of the Gaussian, in pixels, not counting the center
3406 %      pixel.
3407 %
3408 %    o sigma: the standard deviation of the Gaussian, in pixels.
3409 %
3410 %    o threshold: only pixels within this contrast threshold are included
3411 %      in the blur operation.
3412 %
3413 %    o exception: return any errors or warnings in this structure.
3414 %
3415 */
SelectiveBlurImage(const Image * image,const double radius,const double sigma,const double threshold,ExceptionInfo * exception)3416 MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
3417   const double sigma,const double threshold,ExceptionInfo *exception)
3418 {
3419 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
3420 
3421   CacheView
3422     *blur_view,
3423     *image_view,
3424     *luminance_view;
3425 
3426   Image
3427     *blur_image,
3428     *luminance_image;
3429 
3430   MagickBooleanType
3431     status;
3432 
3433   MagickOffsetType
3434     progress;
3435 
3436   MagickRealType
3437     *kernel;
3438 
3439   size_t
3440     width;
3441 
3442   ssize_t
3443     center,
3444     y;
3445 
3446   /*
3447     Initialize blur image attributes.
3448   */
3449   assert(image != (Image *) NULL);
3450   assert(image->signature == MagickCoreSignature);
3451   if (image->debug != MagickFalse)
3452     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3453   assert(exception != (ExceptionInfo *) NULL);
3454   assert(exception->signature == MagickCoreSignature);
3455   width=GetOptimalKernelWidth1D(radius,sigma);
3456   kernel=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory((size_t)
3457     width,width*sizeof(*kernel)));
3458   if (kernel == (MagickRealType *) NULL)
3459     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3460   {
3461     ssize_t
3462       i,
3463       j,
3464       v;
3465 
3466     j=(ssize_t) (width-1)/2;
3467     i=0;
3468     for (v=(-j); v <= j; v++)
3469     {
3470       ssize_t
3471         u;
3472 
3473       for (u=(-j); u <= j; u++)
3474         kernel[i++]=(MagickRealType) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
3475           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
3476     }
3477   }
3478   if (image->debug != MagickFalse)
3479     {
3480       char
3481         format[MagickPathExtent],
3482         *message;
3483 
3484       const MagickRealType
3485         *k;
3486 
3487       ssize_t
3488         u,
3489         v;
3490 
3491       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
3492         "  SelectiveBlurImage with %.20gx%.20g kernel:",(double) width,(double)
3493         width);
3494       message=AcquireString("");
3495       k=kernel;
3496       for (v=0; v < (ssize_t) width; v++)
3497       {
3498         *message='\0';
3499         (void) FormatLocaleString(format,MagickPathExtent,"%.20g: ",(double) v);
3500         (void) ConcatenateString(&message,format);
3501         for (u=0; u < (ssize_t) width; u++)
3502         {
3503           (void) FormatLocaleString(format,MagickPathExtent,"%+f ",(double)
3504             *k++);
3505           (void) ConcatenateString(&message,format);
3506         }
3507         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
3508       }
3509       message=DestroyString(message);
3510     }
3511   blur_image=CloneImage(image,0,0,MagickTrue,exception);
3512   if (blur_image == (Image *) NULL)
3513     return((Image *) NULL);
3514   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
3515     {
3516       blur_image=DestroyImage(blur_image);
3517       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
3518       return((Image *) NULL);
3519     }
3520   luminance_image=CloneImage(image,0,0,MagickTrue,exception);
3521   if (luminance_image == (Image *) NULL)
3522     {
3523       blur_image=DestroyImage(blur_image);
3524       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
3525       return((Image *) NULL);
3526     }
3527   status=TransformImageColorspace(luminance_image,GRAYColorspace,exception);
3528   if (status == MagickFalse)
3529     {
3530       luminance_image=DestroyImage(luminance_image);
3531       blur_image=DestroyImage(blur_image);
3532       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
3533       return((Image *) NULL);
3534     }
3535   /*
3536     Threshold blur image.
3537   */
3538   status=MagickTrue;
3539   progress=0;
3540   center=(ssize_t) (GetPixelChannels(image)*(image->columns+width)*
3541     ((width-1)/2L)+GetPixelChannels(image)*((width-1)/2L));
3542   image_view=AcquireVirtualCacheView(image,exception);
3543   luminance_view=AcquireVirtualCacheView(luminance_image,exception);
3544   blur_view=AcquireAuthenticCacheView(blur_image,exception);
3545 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3546   #pragma omp parallel for schedule(static) shared(progress,status) \
3547     magick_number_threads(image,blur_image,image->rows,1)
3548 #endif
3549   for (y=0; y < (ssize_t) image->rows; y++)
3550   {
3551     double
3552       contrast;
3553 
3554     MagickBooleanType
3555       sync;
3556 
3557     const Quantum
3558       *magick_restrict l,
3559       *magick_restrict p;
3560 
3561     Quantum
3562       *magick_restrict q;
3563 
3564     ssize_t
3565       x;
3566 
3567     if (status == MagickFalse)
3568       continue;
3569     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) (width-1)/2L),y-(ssize_t)
3570       ((width-1)/2L),image->columns+width,width,exception);
3571     l=GetCacheViewVirtualPixels(luminance_view,-((ssize_t) (width-1)/2L),y-
3572       (ssize_t) ((width-1)/2L),luminance_image->columns+width,width,exception);
3573     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
3574       exception);
3575     if ((p == (const Quantum *) NULL) || (l == (const Quantum *) NULL) ||
3576         (q == (Quantum *) NULL))
3577       {
3578         status=MagickFalse;
3579         continue;
3580       }
3581     for (x=0; x < (ssize_t) image->columns; x++)
3582     {
3583       double
3584         intensity;
3585 
3586       ssize_t
3587         i;
3588 
3589       intensity=GetPixelIntensity(image,p+center);
3590       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3591       {
3592         double
3593           alpha,
3594           gamma,
3595           pixel;
3596 
3597         PixelChannel
3598           channel;
3599 
3600         PixelTrait
3601           blur_traits,
3602           traits;
3603 
3604         const MagickRealType
3605           *magick_restrict k;
3606 
3607         const Quantum
3608           *magick_restrict luminance_pixels,
3609           *magick_restrict pixels;
3610 
3611         ssize_t
3612           u;
3613 
3614         ssize_t
3615           v;
3616 
3617         channel=GetPixelChannelChannel(image,i);
3618         traits=GetPixelChannelTraits(image,channel);
3619         blur_traits=GetPixelChannelTraits(blur_image,channel);
3620         if ((traits == UndefinedPixelTrait) ||
3621             (blur_traits == UndefinedPixelTrait))
3622           continue;
3623         if ((blur_traits & CopyPixelTrait) != 0)
3624           {
3625             SetPixelChannel(blur_image,channel,p[center+i],q);
3626             continue;
3627           }
3628         k=kernel;
3629         pixel=0.0;
3630         pixels=p;
3631         luminance_pixels=l;
3632         gamma=0.0;
3633         if ((blur_traits & BlendPixelTrait) == 0)
3634           {
3635             for (v=0; v < (ssize_t) width; v++)
3636             {
3637               for (u=0; u < (ssize_t) width; u++)
3638               {
3639                 contrast=GetPixelIntensity(luminance_image,luminance_pixels)-
3640                   intensity;
3641                 if (fabs(contrast) < threshold)
3642                   {
3643                     pixel+=(*k)*pixels[i];
3644                     gamma+=(*k);
3645                   }
3646                 k++;
3647                 pixels+=GetPixelChannels(image);
3648                 luminance_pixels+=GetPixelChannels(luminance_image);
3649               }
3650               pixels+=GetPixelChannels(image)*image->columns;
3651               luminance_pixels+=GetPixelChannels(luminance_image)*
3652                 luminance_image->columns;
3653             }
3654             if (fabs((double) gamma) < MagickEpsilon)
3655               {
3656                 SetPixelChannel(blur_image,channel,p[center+i],q);
3657                 continue;
3658               }
3659             gamma=PerceptibleReciprocal(gamma);
3660             SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
3661             continue;
3662           }
3663         for (v=0; v < (ssize_t) width; v++)
3664         {
3665           for (u=0; u < (ssize_t) width; u++)
3666           {
3667             contrast=GetPixelIntensity(image,pixels)-intensity;
3668             if (fabs(contrast) < threshold)
3669               {
3670                 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
3671                 pixel+=(*k)*alpha*pixels[i];
3672                 gamma+=(*k)*alpha;
3673               }
3674             k++;
3675             pixels+=GetPixelChannels(image);
3676             luminance_pixels+=GetPixelChannels(luminance_image);
3677           }
3678           pixels+=GetPixelChannels(image)*image->columns;
3679           luminance_pixels+=GetPixelChannels(luminance_image)*
3680             luminance_image->columns;
3681         }
3682         if (fabs((double) gamma) < MagickEpsilon)
3683           {
3684             SetPixelChannel(blur_image,channel,p[center+i],q);
3685             continue;
3686           }
3687         gamma=PerceptibleReciprocal(gamma);
3688         SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
3689       }
3690       p+=GetPixelChannels(image);
3691       l+=GetPixelChannels(luminance_image);
3692       q+=GetPixelChannels(blur_image);
3693     }
3694     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
3695     if (sync == MagickFalse)
3696       status=MagickFalse;
3697     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3698       {
3699         MagickBooleanType
3700           proceed;
3701 
3702 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3703         #pragma omp atomic
3704 #endif
3705         progress++;
3706         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress,
3707           image->rows);
3708         if (proceed == MagickFalse)
3709           status=MagickFalse;
3710       }
3711   }
3712   blur_image->type=image->type;
3713   blur_view=DestroyCacheView(blur_view);
3714   luminance_view=DestroyCacheView(luminance_view);
3715   image_view=DestroyCacheView(image_view);
3716   luminance_image=DestroyImage(luminance_image);
3717   kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
3718   if (status == MagickFalse)
3719     blur_image=DestroyImage(blur_image);
3720   return(blur_image);
3721 }
3722 
3723 /*
3724 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3725 %                                                                             %
3726 %                                                                             %
3727 %                                                                             %
3728 %     S h a d e I m a g e                                                     %
3729 %                                                                             %
3730 %                                                                             %
3731 %                                                                             %
3732 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3733 %
3734 %  ShadeImage() shines a distant light on an image to create a
3735 %  three-dimensional effect. You control the positioning of the light with
3736 %  azimuth and elevation; azimuth is measured in degrees off the x axis
3737 %  and elevation is measured in pixels above the Z axis.
3738 %
3739 %  The format of the ShadeImage method is:
3740 %
3741 %      Image *ShadeImage(const Image *image,const MagickBooleanType gray,
3742 %        const double azimuth,const double elevation,ExceptionInfo *exception)
3743 %
3744 %  A description of each parameter follows:
3745 %
3746 %    o image: the image.
3747 %
3748 %    o gray: A value other than zero shades the intensity of each pixel.
3749 %
3750 %    o azimuth, elevation:  Define the light source direction.
3751 %
3752 %    o exception: return any errors or warnings in this structure.
3753 %
3754 */
ShadeImage(const Image * image,const MagickBooleanType gray,const double azimuth,const double elevation,ExceptionInfo * exception)3755 MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
3756   const double azimuth,const double elevation,ExceptionInfo *exception)
3757 {
3758 #define GetShadeIntensity(image,pixel) \
3759   ClampPixel(GetPixelIntensity((image),(pixel)))
3760 #define ShadeImageTag  "Shade/Image"
3761 
3762   CacheView
3763     *image_view,
3764     *shade_view;
3765 
3766   Image
3767     *linear_image,
3768     *shade_image;
3769 
3770   MagickBooleanType
3771     status;
3772 
3773   MagickOffsetType
3774     progress;
3775 
3776   PrimaryInfo
3777     light;
3778 
3779   ssize_t
3780     y;
3781 
3782   /*
3783     Initialize shaded image attributes.
3784   */
3785   assert(image != (const Image *) NULL);
3786   assert(image->signature == MagickCoreSignature);
3787   if (image->debug != MagickFalse)
3788     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3789   assert(exception != (ExceptionInfo *) NULL);
3790   assert(exception->signature == MagickCoreSignature);
3791   linear_image=CloneImage(image,0,0,MagickTrue,exception);
3792   shade_image=CloneImage(image,0,0,MagickTrue,exception);
3793   if ((linear_image == (Image *) NULL) || (shade_image == (Image *) NULL))
3794     {
3795       if (linear_image != (Image *) NULL)
3796         linear_image=DestroyImage(linear_image);
3797       if (shade_image != (Image *) NULL)
3798         shade_image=DestroyImage(shade_image);
3799       return((Image *) NULL);
3800     }
3801   if (SetImageStorageClass(shade_image,DirectClass,exception) == MagickFalse)
3802     {
3803       linear_image=DestroyImage(linear_image);
3804       shade_image=DestroyImage(shade_image);
3805       return((Image *) NULL);
3806     }
3807   /*
3808     Compute the light vector.
3809   */
3810   light.x=(double) QuantumRange*cos(DegreesToRadians(azimuth))*
3811     cos(DegreesToRadians(elevation));
3812   light.y=(double) QuantumRange*sin(DegreesToRadians(azimuth))*
3813     cos(DegreesToRadians(elevation));
3814   light.z=(double) QuantumRange*sin(DegreesToRadians(elevation));
3815   /*
3816     Shade image.
3817   */
3818   status=MagickTrue;
3819   progress=0;
3820   image_view=AcquireVirtualCacheView(linear_image,exception);
3821   shade_view=AcquireAuthenticCacheView(shade_image,exception);
3822 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3823   #pragma omp parallel for schedule(static) shared(progress,status) \
3824     magick_number_threads(linear_image,shade_image,linear_image->rows,1)
3825 #endif
3826   for (y=0; y < (ssize_t) linear_image->rows; y++)
3827   {
3828     double
3829       distance,
3830       normal_distance,
3831       shade;
3832 
3833     PrimaryInfo
3834       normal;
3835 
3836     const Quantum
3837       *magick_restrict center,
3838       *magick_restrict p,
3839       *magick_restrict post,
3840       *magick_restrict pre;
3841 
3842     Quantum
3843       *magick_restrict q;
3844 
3845     ssize_t
3846       x;
3847 
3848     if (status == MagickFalse)
3849       continue;
3850     p=GetCacheViewVirtualPixels(image_view,-1,y-1,linear_image->columns+2,3,
3851       exception);
3852     q=QueueCacheViewAuthenticPixels(shade_view,0,y,shade_image->columns,1,
3853       exception);
3854     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3855       {
3856         status=MagickFalse;
3857         continue;
3858       }
3859     /*
3860       Shade this row of pixels.
3861     */
3862     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
3863     for (x=0; x < (ssize_t) linear_image->columns; x++)
3864     {
3865       ssize_t
3866         i;
3867 
3868       /*
3869         Determine the surface normal and compute shading.
3870       */
3871       pre=p+GetPixelChannels(linear_image);
3872       center=pre+(linear_image->columns+2)*GetPixelChannels(linear_image);
3873       post=center+(linear_image->columns+2)*GetPixelChannels(linear_image);
3874       normal.x=(double) (
3875         GetShadeIntensity(linear_image,pre-GetPixelChannels(linear_image))+
3876         GetShadeIntensity(linear_image,center-GetPixelChannels(linear_image))+
3877         GetShadeIntensity(linear_image,post-GetPixelChannels(linear_image))-
3878         GetShadeIntensity(linear_image,pre+GetPixelChannels(linear_image))-
3879         GetShadeIntensity(linear_image,center+GetPixelChannels(linear_image))-
3880         GetShadeIntensity(linear_image,post+GetPixelChannels(linear_image)));
3881       normal.y=(double) (
3882         GetShadeIntensity(linear_image,post-GetPixelChannels(linear_image))+
3883         GetShadeIntensity(linear_image,post)+
3884         GetShadeIntensity(linear_image,post+GetPixelChannels(linear_image))-
3885         GetShadeIntensity(linear_image,pre-GetPixelChannels(linear_image))-
3886         GetShadeIntensity(linear_image,pre)-
3887         GetShadeIntensity(linear_image,pre+GetPixelChannels(linear_image)));
3888       if ((fabs(normal.x) <= MagickEpsilon) &&
3889           (fabs(normal.y) <= MagickEpsilon))
3890         shade=light.z;
3891       else
3892         {
3893           shade=0.0;
3894           distance=normal.x*light.x+normal.y*light.y+normal.z*light.z;
3895           if (distance > MagickEpsilon)
3896             {
3897               normal_distance=normal.x*normal.x+normal.y*normal.y+
3898                 normal.z*normal.z;
3899               if (normal_distance > (MagickEpsilon*MagickEpsilon))
3900                 shade=distance/sqrt((double) normal_distance);
3901             }
3902         }
3903       for (i=0; i < (ssize_t) GetPixelChannels(linear_image); i++)
3904       {
3905         PixelChannel
3906           channel;
3907 
3908         PixelTrait
3909           shade_traits,
3910           traits;
3911 
3912         channel=GetPixelChannelChannel(linear_image,i);
3913         traits=GetPixelChannelTraits(linear_image,channel);
3914         shade_traits=GetPixelChannelTraits(shade_image,channel);
3915         if ((traits == UndefinedPixelTrait) ||
3916             (shade_traits == UndefinedPixelTrait))
3917           continue;
3918         if ((shade_traits & CopyPixelTrait) != 0)
3919           {
3920             SetPixelChannel(shade_image,channel,center[i],q);
3921             continue;
3922           }
3923         if ((traits & UpdatePixelTrait) == 0)
3924           {
3925             SetPixelChannel(shade_image,channel,center[i],q);
3926             continue;
3927           }
3928         if (gray != MagickFalse)
3929           {
3930             SetPixelChannel(shade_image,channel,ClampToQuantum(shade),q);
3931             continue;
3932           }
3933         SetPixelChannel(shade_image,channel,ClampToQuantum(QuantumScale*shade*
3934           center[i]),q);
3935       }
3936       p+=GetPixelChannels(linear_image);
3937       q+=GetPixelChannels(shade_image);
3938     }
3939     if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
3940       status=MagickFalse;
3941     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3942       {
3943         MagickBooleanType
3944           proceed;
3945 
3946 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3947         #pragma omp atomic
3948 #endif
3949         progress++;
3950         proceed=SetImageProgress(image,ShadeImageTag,progress,image->rows);
3951         if (proceed == MagickFalse)
3952           status=MagickFalse;
3953       }
3954   }
3955   shade_view=DestroyCacheView(shade_view);
3956   image_view=DestroyCacheView(image_view);
3957   linear_image=DestroyImage(linear_image);
3958   if (status == MagickFalse)
3959     shade_image=DestroyImage(shade_image);
3960   return(shade_image);
3961 }
3962 
3963 /*
3964 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3965 %                                                                             %
3966 %                                                                             %
3967 %                                                                             %
3968 %     S h a r p e n I m a g e                                                 %
3969 %                                                                             %
3970 %                                                                             %
3971 %                                                                             %
3972 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3973 %
3974 %  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
3975 %  operator of the given radius and standard deviation (sigma).  For
3976 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
3977 %  and SharpenImage() selects a suitable radius for you.
3978 %
3979 %  Using a separable kernel would be faster, but the negative weights cancel
3980 %  out on the corners of the kernel producing often undesirable ringing in the
3981 %  filtered result; this can be avoided by using a 2D gaussian shaped image
3982 %  sharpening kernel instead.
3983 %
3984 %  The format of the SharpenImage method is:
3985 %
3986 %    Image *SharpenImage(const Image *image,const double radius,
3987 %      const double sigma,ExceptionInfo *exception)
3988 %
3989 %  A description of each parameter follows:
3990 %
3991 %    o image: the image.
3992 %
3993 %    o radius: the radius of the Gaussian, in pixels, not counting the center
3994 %      pixel.
3995 %
3996 %    o sigma: the standard deviation of the Laplacian, in pixels.
3997 %
3998 %    o exception: return any errors or warnings in this structure.
3999 %
4000 */
SharpenImage(const Image * image,const double radius,const double sigma,ExceptionInfo * exception)4001 MagickExport Image *SharpenImage(const Image *image,const double radius,
4002   const double sigma,ExceptionInfo *exception)
4003 {
4004   double
4005     gamma,
4006     normalize;
4007 
4008   Image
4009     *sharp_image;
4010 
4011   KernelInfo
4012     *kernel_info;
4013 
4014   ssize_t
4015     i;
4016 
4017   size_t
4018     width;
4019 
4020   ssize_t
4021     j,
4022     u,
4023     v;
4024 
4025   assert(image != (const Image *) NULL);
4026   assert(image->signature == MagickCoreSignature);
4027   if (image->debug != MagickFalse)
4028     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4029   assert(exception != (ExceptionInfo *) NULL);
4030   assert(exception->signature == MagickCoreSignature);
4031   width=GetOptimalKernelWidth2D(radius,sigma);
4032   kernel_info=AcquireKernelInfo((const char *) NULL,exception);
4033   if (kernel_info == (KernelInfo *) NULL)
4034     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4035   (void) memset(kernel_info,0,sizeof(*kernel_info));
4036   kernel_info->width=width;
4037   kernel_info->height=width;
4038   kernel_info->x=(ssize_t) (width-1)/2;
4039   kernel_info->y=(ssize_t) (width-1)/2;
4040   kernel_info->signature=MagickCoreSignature;
4041   kernel_info->values=(MagickRealType *) MagickAssumeAligned(
4042     AcquireAlignedMemory(kernel_info->width,kernel_info->height*
4043     sizeof(*kernel_info->values)));
4044   if (kernel_info->values == (MagickRealType *) NULL)
4045     {
4046       kernel_info=DestroyKernelInfo(kernel_info);
4047       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
4048     }
4049   normalize=0.0;
4050   j=(ssize_t) (kernel_info->width-1)/2;
4051   i=0;
4052   for (v=(-j); v <= j; v++)
4053   {
4054     for (u=(-j); u <= j; u++)
4055     {
4056       kernel_info->values[i]=(MagickRealType) (-exp(-((double) u*u+v*v)/(2.0*
4057         MagickSigma*MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
4058       normalize+=kernel_info->values[i];
4059       i++;
4060     }
4061   }
4062   kernel_info->values[i/2]=(double) ((-2.0)*normalize);
4063   normalize=0.0;
4064   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
4065     normalize+=kernel_info->values[i];
4066   gamma=PerceptibleReciprocal(normalize);
4067   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
4068     kernel_info->values[i]*=gamma;
4069   sharp_image=ConvolveImage(image,kernel_info,exception);
4070   kernel_info=DestroyKernelInfo(kernel_info);
4071   return(sharp_image);
4072 }
4073 
4074 /*
4075 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4076 %                                                                             %
4077 %                                                                             %
4078 %                                                                             %
4079 %     S p r e a d I m a g e                                                   %
4080 %                                                                             %
4081 %                                                                             %
4082 %                                                                             %
4083 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4084 %
4085 %  SpreadImage() is a special effects method that randomly displaces each
4086 %  pixel in a square area defined by the radius parameter.
4087 %
4088 %  The format of the SpreadImage method is:
4089 %
4090 %      Image *SpreadImage(const Image *image,
4091 %        const PixelInterpolateMethod method,const double radius,
4092 %        ExceptionInfo *exception)
4093 %
4094 %  A description of each parameter follows:
4095 %
4096 %    o image: the image.
4097 %
4098 %    o method:  intepolation method.
4099 %
4100 %    o radius:  choose a random pixel in a neighborhood of this extent.
4101 %
4102 %    o exception: return any errors or warnings in this structure.
4103 %
4104 */
SpreadImage(const Image * image,const PixelInterpolateMethod method,const double radius,ExceptionInfo * exception)4105 MagickExport Image *SpreadImage(const Image *image,
4106   const PixelInterpolateMethod method,const double radius,
4107   ExceptionInfo *exception)
4108 {
4109 #define SpreadImageTag  "Spread/Image"
4110 
4111   CacheView
4112     *image_view,
4113     *spread_view;
4114 
4115   Image
4116     *spread_image;
4117 
4118   MagickBooleanType
4119     status;
4120 
4121   MagickOffsetType
4122     progress;
4123 
4124   RandomInfo
4125     **magick_restrict random_info;
4126 
4127   size_t
4128     width;
4129 
4130   ssize_t
4131     y;
4132 
4133 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4134   unsigned long
4135     key;
4136 #endif
4137 
4138   /*
4139     Initialize spread image attributes.
4140   */
4141   assert(image != (Image *) NULL);
4142   assert(image->signature == MagickCoreSignature);
4143   if (image->debug != MagickFalse)
4144     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4145   assert(exception != (ExceptionInfo *) NULL);
4146   assert(exception->signature == MagickCoreSignature);
4147   spread_image=CloneImage(image,0,0,MagickTrue,exception);
4148   if (spread_image == (Image *) NULL)
4149     return((Image *) NULL);
4150   if (SetImageStorageClass(spread_image,DirectClass,exception) == MagickFalse)
4151     {
4152       spread_image=DestroyImage(spread_image);
4153       return((Image *) NULL);
4154     }
4155   /*
4156     Spread image.
4157   */
4158   status=MagickTrue;
4159   progress=0;
4160   width=GetOptimalKernelWidth1D(radius,0.5);
4161   random_info=AcquireRandomInfoThreadSet();
4162   image_view=AcquireVirtualCacheView(image,exception);
4163   spread_view=AcquireAuthenticCacheView(spread_image,exception);
4164 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4165   key=GetRandomSecretKey(random_info[0]);
4166   #pragma omp parallel for schedule(static) shared(progress,status) \
4167     magick_number_threads(image,spread_image,image->rows,key == ~0UL)
4168 #endif
4169   for (y=0; y < (ssize_t) image->rows; y++)
4170   {
4171     const int
4172       id = GetOpenMPThreadId();
4173 
4174     Quantum
4175       *magick_restrict q;
4176 
4177     ssize_t
4178       x;
4179 
4180     if (status == MagickFalse)
4181       continue;
4182     q=QueueCacheViewAuthenticPixels(spread_view,0,y,spread_image->columns,1,
4183       exception);
4184     if (q == (Quantum *) NULL)
4185       {
4186         status=MagickFalse;
4187         continue;
4188       }
4189     for (x=0; x < (ssize_t) image->columns; x++)
4190     {
4191       PointInfo
4192         point;
4193 
4194       point.x=GetPseudoRandomValue(random_info[id]);
4195       point.y=GetPseudoRandomValue(random_info[id]);
4196       status=InterpolatePixelChannels(image,image_view,spread_image,method,
4197         (double) x+width*(point.x-0.5),(double) y+width*(point.y-0.5),q,
4198         exception);
4199       if (status == MagickFalse)
4200         break;
4201       q+=GetPixelChannels(spread_image);
4202     }
4203     if (SyncCacheViewAuthenticPixels(spread_view,exception) == MagickFalse)
4204       status=MagickFalse;
4205     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4206       {
4207         MagickBooleanType
4208           proceed;
4209 
4210 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4211         #pragma omp atomic
4212 #endif
4213         progress++;
4214         proceed=SetImageProgress(image,SpreadImageTag,progress,image->rows);
4215         if (proceed == MagickFalse)
4216           status=MagickFalse;
4217       }
4218   }
4219   spread_view=DestroyCacheView(spread_view);
4220   image_view=DestroyCacheView(image_view);
4221   random_info=DestroyRandomInfoThreadSet(random_info);
4222   if (status == MagickFalse)
4223     spread_image=DestroyImage(spread_image);
4224   return(spread_image);
4225 }
4226 
4227 /*
4228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4229 %                                                                             %
4230 %                                                                             %
4231 %                                                                             %
4232 %     U n s h a r p M a s k I m a g e                                         %
4233 %                                                                             %
4234 %                                                                             %
4235 %                                                                             %
4236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4237 %
4238 %  UnsharpMaskImage() sharpens one or more image channels.  We convolve the
4239 %  image with a Gaussian operator of the given radius and standard deviation
4240 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
4241 %  radius of 0 and UnsharpMaskImage() selects a suitable radius for you.
4242 %
4243 %  The format of the UnsharpMaskImage method is:
4244 %
4245 %    Image *UnsharpMaskImage(const Image *image,const double radius,
4246 %      const double sigma,const double amount,const double threshold,
4247 %      ExceptionInfo *exception)
4248 %
4249 %  A description of each parameter follows:
4250 %
4251 %    o image: the image.
4252 %
4253 %    o radius: the radius of the Gaussian, in pixels, not counting the center
4254 %      pixel.
4255 %
4256 %    o sigma: the standard deviation of the Gaussian, in pixels.
4257 %
4258 %    o gain: the percentage of the difference between the original and the
4259 %      blur image that is added back into the original.
4260 %
4261 %    o threshold: the threshold in pixels needed to apply the diffence gain.
4262 %
4263 %    o exception: return any errors or warnings in this structure.
4264 %
4265 */
UnsharpMaskImage(const Image * image,const double radius,const double sigma,const double gain,const double threshold,ExceptionInfo * exception)4266 MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
4267   const double sigma,const double gain,const double threshold,
4268   ExceptionInfo *exception)
4269 {
4270 #define SharpenImageTag  "Sharpen/Image"
4271 
4272   CacheView
4273     *image_view,
4274     *unsharp_view;
4275 
4276   Image
4277     *unsharp_image;
4278 
4279   MagickBooleanType
4280     status;
4281 
4282   MagickOffsetType
4283     progress;
4284 
4285   double
4286     quantum_threshold;
4287 
4288   ssize_t
4289     y;
4290 
4291   assert(image != (const Image *) NULL);
4292   assert(image->signature == MagickCoreSignature);
4293   if (image->debug != MagickFalse)
4294     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
4295   assert(exception != (ExceptionInfo *) NULL);
4296 /* This kernel appears to be broken.
4297 #if defined(MAGICKCORE_OPENCL_SUPPORT)
4298   unsharp_image=AccelerateUnsharpMaskImage(image,radius,sigma,gain,threshold,
4299     exception);
4300   if (unsharp_image != (Image *) NULL)
4301     return(unsharp_image);
4302 #endif
4303 */
4304   unsharp_image=BlurImage(image,radius,sigma,exception);
4305   if (unsharp_image == (Image *) NULL)
4306     return((Image *) NULL);
4307   quantum_threshold=(double) QuantumRange*threshold;
4308   /*
4309     Unsharp-mask image.
4310   */
4311   status=MagickTrue;
4312   progress=0;
4313   image_view=AcquireVirtualCacheView(image,exception);
4314   unsharp_view=AcquireAuthenticCacheView(unsharp_image,exception);
4315 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4316   #pragma omp parallel for schedule(static) shared(progress,status) \
4317     magick_number_threads(image,unsharp_image,image->rows,1)
4318 #endif
4319   for (y=0; y < (ssize_t) image->rows; y++)
4320   {
4321     const Quantum
4322       *magick_restrict p;
4323 
4324     Quantum
4325       *magick_restrict q;
4326 
4327     ssize_t
4328       x;
4329 
4330     if (status == MagickFalse)
4331       continue;
4332     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
4333     q=GetCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
4334       exception);
4335     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
4336       {
4337         status=MagickFalse;
4338         continue;
4339       }
4340     for (x=0; x < (ssize_t) image->columns; x++)
4341     {
4342       ssize_t
4343         i;
4344 
4345       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
4346       {
4347         double
4348           pixel;
4349 
4350         PixelChannel
4351           channel;
4352 
4353         PixelTrait
4354           traits,
4355           unsharp_traits;
4356 
4357         channel=GetPixelChannelChannel(image,i);
4358         traits=GetPixelChannelTraits(image,channel);
4359         unsharp_traits=GetPixelChannelTraits(unsharp_image,channel);
4360         if ((traits == UndefinedPixelTrait) ||
4361             (unsharp_traits == UndefinedPixelTrait))
4362           continue;
4363         if ((unsharp_traits & CopyPixelTrait) != 0)
4364           {
4365             SetPixelChannel(unsharp_image,channel,p[i],q);
4366             continue;
4367           }
4368         pixel=p[i]-(double) GetPixelChannel(unsharp_image,channel,q);
4369         if (fabs(2.0*pixel) < quantum_threshold)
4370           pixel=(double) p[i];
4371         else
4372           pixel=(double) p[i]+gain*pixel;
4373         SetPixelChannel(unsharp_image,channel,ClampToQuantum(pixel),q);
4374       }
4375       p+=GetPixelChannels(image);
4376       q+=GetPixelChannels(unsharp_image);
4377     }
4378     if (SyncCacheViewAuthenticPixels(unsharp_view,exception) == MagickFalse)
4379       status=MagickFalse;
4380     if (image->progress_monitor != (MagickProgressMonitor) NULL)
4381       {
4382         MagickBooleanType
4383           proceed;
4384 
4385 #if defined(MAGICKCORE_OPENMP_SUPPORT)
4386         #pragma omp atomic
4387 #endif
4388         progress++;
4389         proceed=SetImageProgress(image,SharpenImageTag,progress,image->rows);
4390         if (proceed == MagickFalse)
4391           status=MagickFalse;
4392       }
4393   }
4394   unsharp_image->type=image->type;
4395   unsharp_view=DestroyCacheView(unsharp_view);
4396   image_view=DestroyCacheView(image_view);
4397   if (status == MagickFalse)
4398     unsharp_image=DestroyImage(unsharp_image);
4399   return(unsharp_image);
4400 }
4401