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