]> granicus.if.org Git - imagemagick/blob - MagickCore/effect.c
(no commit message)
[imagemagick] / MagickCore / effect.c
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 %                                 John Cristy                                 %
17 %                                 October 1996                                %
18 %                                                                             %
19 %                                                                             %
20 %  Copyright 1999-2013 ImageMagick Studio LLC, a non-profit organization      %
21 %  dedicated to making software imaging solutions freely available.           %
22 %                                                                             %
23 %  You may not use this file except in compliance with the License.  You may  %
24 %  obtain a copy of the License at                                            %
25 %                                                                             %
26 %    http://www.imagemagick.org/script/license.php                            %
27 %                                                                             %
28 %  Unless required by applicable law or agreed to in writing, software        %
29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
31 %  See the License for the specific language governing permissions and        %
32 %  limitations under the License.                                             %
33 %                                                                             %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39 \f
40 /*
41   Include declarations.
42 */
43 #include "MagickCore/studio.h"
44 #include "MagickCore/accelerate.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/memory_.h"
66 #include "MagickCore/memory-private.h"
67 #include "MagickCore/monitor.h"
68 #include "MagickCore/monitor-private.h"
69 #include "MagickCore/montage.h"
70 #include "MagickCore/morphology.h"
71 #include "MagickCore/morphology-private.h"
72 #include "MagickCore/paint.h"
73 #include "MagickCore/pixel-accessor.h"
74 #include "MagickCore/pixel-private.h"
75 #include "MagickCore/property.h"
76 #include "MagickCore/quantize.h"
77 #include "MagickCore/quantum.h"
78 #include "MagickCore/quantum-private.h"
79 #include "MagickCore/random_.h"
80 #include "MagickCore/random-private.h"
81 #include "MagickCore/resample.h"
82 #include "MagickCore/resample-private.h"
83 #include "MagickCore/resize.h"
84 #include "MagickCore/resource_.h"
85 #include "MagickCore/segment.h"
86 #include "MagickCore/shear.h"
87 #include "MagickCore/signature-private.h"
88 #include "MagickCore/statistic.h"
89 #include "MagickCore/string_.h"
90 #include "MagickCore/thread-private.h"
91 #include "MagickCore/transform.h"
92 #include "MagickCore/threshold.h"
93 \f
94 /*
95 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96 %                                                                             %
97 %                                                                             %
98 %                                                                             %
99 %     A d a p t i v e B l u r I m a g e                                       %
100 %                                                                             %
101 %                                                                             %
102 %                                                                             %
103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
104 %
105 %  AdaptiveBlurImage() adaptively blurs the image by blurring less
106 %  intensely near image edges and more intensely far from edges.  We blur the
107 %  image with a Gaussian operator of the given radius and standard deviation
108 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
109 %  radius of 0 and AdaptiveBlurImage() selects a suitable radius for you.
110 %
111 %  The format of the AdaptiveBlurImage method is:
112 %
113 %      Image *AdaptiveBlurImage(const Image *image,const double radius,
114 %        const double sigma,ExceptionInfo *exception)
115 %
116 %  A description of each parameter follows:
117 %
118 %    o image: the image.
119 %
120 %    o radius: the radius of the Gaussian, in pixels, not counting the center
121 %      pixel.
122 %
123 %    o sigma: the standard deviation of the Laplacian, in pixels.
124 %
125 %    o exception: return any errors or warnings in this structure.
126 %
127 */
128
129 MagickExport MagickBooleanType AdaptiveLevelImage(Image *image,
130   const char *levels,ExceptionInfo *exception)
131 {
132   double
133     black_point,
134     gamma,
135     white_point;
136
137   GeometryInfo
138     geometry_info;
139
140   MagickBooleanType
141     status;
142
143   MagickStatusType
144     flags;
145
146   /*
147     Parse levels.
148   */
149   if (levels == (char *) NULL)
150     return(MagickFalse);
151   flags=ParseGeometry(levels,&geometry_info);
152   black_point=geometry_info.rho;
153   white_point=(double) QuantumRange;
154   if ((flags & SigmaValue) != 0)
155     white_point=geometry_info.sigma;
156   gamma=1.0;
157   if ((flags & XiValue) != 0)
158     gamma=geometry_info.xi;
159   if ((flags & PercentValue) != 0)
160     {
161       black_point*=(double) image->columns*image->rows/100.0;
162       white_point*=(double) image->columns*image->rows/100.0;
163     }
164   if ((flags & SigmaValue) == 0)
165     white_point=(double) QuantumRange-black_point;
166   if ((flags & AspectValue ) == 0)
167     status=LevelImage(image,black_point,white_point,gamma,exception);
168   else
169     status=LevelizeImage(image,black_point,white_point,gamma,exception);
170   return(status);
171 }
172
173 MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
174   const double sigma,ExceptionInfo *exception)
175 {
176 #define AdaptiveBlurImageTag  "Convolve/Image"
177 #define MagickSigma  (fabs(sigma) < MagickEpsilon ? MagickEpsilon : sigma)
178
179   CacheView
180     *blur_view,
181     *edge_view,
182     *image_view;
183
184   double
185     normalize;
186
187   Image
188     *blur_image,
189     *edge_image,
190     *gaussian_image;
191
192   MagickBooleanType
193     status;
194
195   MagickOffsetType
196     progress;
197
198   MagickRealType
199     **kernel;
200
201   register ssize_t
202     i;
203
204   size_t
205     width;
206
207   ssize_t
208     j,
209     k,
210     u,
211     v,
212     y;
213
214   assert(image != (const Image *) NULL);
215   assert(image->signature == MagickSignature);
216   if (image->debug != MagickFalse)
217     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
218   assert(exception != (ExceptionInfo *) NULL);
219   assert(exception->signature == MagickSignature);
220   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
221   if (blur_image == (Image *) NULL)
222     return((Image *) NULL);
223   if (fabs(sigma) < MagickEpsilon)
224     return(blur_image);
225   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
226     {
227       blur_image=DestroyImage(blur_image);
228       return((Image *) NULL);
229     }
230   /*
231     Edge detect the image brighness channel, level, blur, and level again.
232   */
233   edge_image=EdgeImage(image,radius,exception);
234   if (edge_image == (Image *) NULL)
235     {
236       blur_image=DestroyImage(blur_image);
237       return((Image *) NULL);
238     }
239   (void) AdaptiveLevelImage(edge_image,"20%,95%",exception);
240   gaussian_image=BlurImage(edge_image,radius,sigma,exception);
241   if (gaussian_image != (Image *) NULL)
242     {
243       edge_image=DestroyImage(edge_image);
244       edge_image=gaussian_image;
245     }
246   (void) AdaptiveLevelImage(edge_image,"10%,95%",exception);
247   /*
248     Create a set of kernels from maximum (radius,sigma) to minimum.
249   */
250   width=GetOptimalKernelWidth2D(radius,sigma);
251   kernel=(MagickRealType **) MagickAssumeAligned(AcquireAlignedMemory((size_t)
252     width,sizeof(*kernel)));
253   if (kernel == (MagickRealType  **) NULL)
254     {
255       edge_image=DestroyImage(edge_image);
256       blur_image=DestroyImage(blur_image);
257       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
258     }
259   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
260   for (i=0; i < (ssize_t) width; i+=2)
261   {
262     kernel[i]=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
263       (size_t) (width-i),(width-i)*sizeof(**kernel)));
264     if (kernel[i] == (MagickRealType *) NULL)
265       break;
266     normalize=0.0;
267     j=(ssize_t) (width-i-1)/2;
268     k=0;
269     for (v=(-j); v <= j; v++)
270     {
271       for (u=(-j); u <= j; u++)
272       {
273         kernel[i][k]=(MagickRealType) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
274           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
275         normalize+=kernel[i][k];
276         k++;
277       }
278     }
279     kernel[i][(j-1)/2]+=(1.0-normalize);
280     if (sigma < MagickEpsilon)
281       kernel[i][(j-1)/2]=1.0;
282   }
283   if (i < (ssize_t) width)
284     {
285       for (i-=2; i >= 0; i-=2)
286         kernel[i]=(MagickRealType *) RelinquishAlignedMemory(kernel[i]);
287       kernel=(MagickRealType **) RelinquishAlignedMemory(kernel);
288       edge_image=DestroyImage(edge_image);
289       blur_image=DestroyImage(blur_image);
290       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
291     }
292   /*
293     Adaptively blur image.
294   */
295   status=MagickTrue;
296   progress=0;
297   image_view=AcquireVirtualCacheView(image,exception);
298   edge_view=AcquireVirtualCacheView(edge_image,exception);
299   blur_view=AcquireAuthenticCacheView(blur_image,exception);
300 #if defined(MAGICKCORE_OPENMP_SUPPORT)
301   #pragma omp parallel for schedule(static,4) shared(progress,status) \
302     magick_threads(image,blur_image,blur_image->rows,1)
303 #endif
304   for (y=0; y < (ssize_t) blur_image->rows; y++)
305   {
306     register const Quantum
307       *restrict r;
308
309     register Quantum
310       *restrict q;
311
312     register ssize_t
313       x;
314
315     if (status == MagickFalse)
316       continue;
317     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
318     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
319       exception);
320     if ((r == (const Quantum *) NULL) || (q == (Quantum *) NULL))
321       {
322         status=MagickFalse;
323         continue;
324       }
325     for (x=0; x < (ssize_t) blur_image->columns; x++)
326     {
327       register const Quantum
328         *restrict p;
329
330       register ssize_t
331         i;
332
333       ssize_t
334         center,
335         j;
336
337       j=(ssize_t) ceil((double) width*QuantumScale*
338         GetPixelIntensity(edge_image,r)-0.5);
339       if (j < 0)
340         j=0;
341       else
342         if (j > (ssize_t) width)
343           j=(ssize_t) width;
344       if ((j & 0x01) != 0)
345         j--;
346       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-j)/2L),y-
347         (ssize_t) ((width-j)/2L),width-j,width-j,exception);
348       if (p == (const Quantum *) NULL)
349         break;
350       center=(ssize_t) GetPixelChannels(image)*(width-j)*((width-j)/2L)+
351         GetPixelChannels(image)*((width-j)/2L);
352       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
353       {
354         double
355           alpha,
356           gamma,
357           pixel;
358
359         PixelChannel
360           channel;
361
362         PixelTrait
363           blur_traits,
364           traits;
365
366         register const MagickRealType
367           *restrict k;
368
369         register const Quantum
370           *restrict pixels;
371
372         register ssize_t
373           u;
374
375         ssize_t
376           v;
377
378         channel=GetPixelChannelChannel(image,i);
379         traits=GetPixelChannelTraits(image,channel);
380         blur_traits=GetPixelChannelTraits(blur_image,channel);
381         if ((traits == UndefinedPixelTrait) ||
382             (blur_traits == UndefinedPixelTrait))
383           continue;
384         if (((blur_traits & CopyPixelTrait) != 0) ||
385             (GetPixelReadMask(image,p+center) == 0))
386           {
387             SetPixelChannel(blur_image,channel,p[center+i],q);
388             continue;
389           }
390         k=kernel[j];
391         pixels=p;
392         pixel=0.0;
393         gamma=0.0;
394         if ((blur_traits & BlendPixelTrait) == 0)
395           {
396             /*
397               No alpha blending.
398             */
399             for (v=0; v < (ssize_t) (width-j); v++)
400             {
401               for (u=0; u < (ssize_t) (width-j); u++)
402               {
403                 pixel+=(*k)*pixels[i];
404                 gamma+=(*k);
405                 k++;
406                 pixels+=GetPixelChannels(image);
407               }
408             }
409             gamma=PerceptibleReciprocal(gamma);
410             SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
411             continue;
412           }
413         /*
414           Alpha blending.
415         */
416         for (v=0; v < (ssize_t) (width-j); v++)
417         {
418           for (u=0; u < (ssize_t) (width-j); u++)
419           {
420             alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
421             pixel+=(*k)*alpha*pixels[i];
422             gamma+=(*k)*alpha;
423             k++;
424             pixels+=GetPixelChannels(image);
425           }
426         }
427         gamma=PerceptibleReciprocal(gamma);
428         SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
429       }
430       q+=GetPixelChannels(blur_image);
431       r+=GetPixelChannels(edge_image);
432     }
433     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
434       status=MagickFalse;
435     if (image->progress_monitor != (MagickProgressMonitor) NULL)
436       {
437         MagickBooleanType
438           proceed;
439
440 #if defined(MAGICKCORE_OPENMP_SUPPORT)
441         #pragma omp critical (MagickCore_AdaptiveBlurImage)
442 #endif
443         proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress++,
444           image->rows);
445         if (proceed == MagickFalse)
446           status=MagickFalse;
447       }
448   }
449   blur_image->type=image->type;
450   blur_view=DestroyCacheView(blur_view);
451   edge_view=DestroyCacheView(edge_view);
452   image_view=DestroyCacheView(image_view);
453   edge_image=DestroyImage(edge_image);
454   for (i=0; i < (ssize_t) width;  i+=2)
455     kernel[i]=(MagickRealType *) RelinquishAlignedMemory(kernel[i]);
456   kernel=(MagickRealType **) RelinquishAlignedMemory(kernel);
457   if (status == MagickFalse)
458     blur_image=DestroyImage(blur_image);
459   return(blur_image);
460 }
461 \f
462 /*
463 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
464 %                                                                             %
465 %                                                                             %
466 %                                                                             %
467 %     A d a p t i v e S h a r p e n I m a g e                                 %
468 %                                                                             %
469 %                                                                             %
470 %                                                                             %
471 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
472 %
473 %  AdaptiveSharpenImage() adaptively sharpens the image by sharpening more
474 %  intensely near image edges and less intensely far from edges. We sharpen the
475 %  image with a Gaussian operator of the given radius and standard deviation
476 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
477 %  radius of 0 and AdaptiveSharpenImage() selects a suitable radius for you.
478 %
479 %  The format of the AdaptiveSharpenImage method is:
480 %
481 %      Image *AdaptiveSharpenImage(const Image *image,const double radius,
482 %        const double sigma,ExceptionInfo *exception)
483 %
484 %  A description of each parameter follows:
485 %
486 %    o image: the image.
487 %
488 %    o radius: the radius of the Gaussian, in pixels, not counting the center
489 %      pixel.
490 %
491 %    o sigma: the standard deviation of the Laplacian, in pixels.
492 %
493 %    o exception: return any errors or warnings in this structure.
494 %
495 */
496 MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
497   const double sigma,ExceptionInfo *exception)
498 {
499 #define AdaptiveSharpenImageTag  "Convolve/Image"
500 #define MagickSigma  (fabs(sigma) < MagickEpsilon ? MagickEpsilon : sigma)
501
502   CacheView
503     *sharp_view,
504     *edge_view,
505     *image_view;
506
507   double
508     normalize;
509
510   Image
511     *sharp_image,
512     *edge_image,
513     *gaussian_image;
514
515   MagickBooleanType
516     status;
517
518   MagickOffsetType
519     progress;
520
521   MagickRealType
522     **kernel;
523
524   register ssize_t
525     i;
526
527   size_t
528     width;
529
530   ssize_t
531     j,
532     k,
533     u,
534     v,
535     y;
536
537   assert(image != (const Image *) NULL);
538   assert(image->signature == MagickSignature);
539   if (image->debug != MagickFalse)
540     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
541   assert(exception != (ExceptionInfo *) NULL);
542   assert(exception->signature == MagickSignature);
543   sharp_image=CloneImage(image,0,0,MagickTrue,exception);
544   if (sharp_image == (Image *) NULL)
545     return((Image *) NULL);
546   if (fabs(sigma) < MagickEpsilon)
547     return(sharp_image);
548   if (SetImageStorageClass(sharp_image,DirectClass,exception) == MagickFalse)
549     {
550       sharp_image=DestroyImage(sharp_image);
551       return((Image *) NULL);
552     }
553   /*
554     Edge detect the image brighness channel, level, sharp, and level again.
555   */
556   edge_image=EdgeImage(image,radius,exception);
557   if (edge_image == (Image *) NULL)
558     {
559       sharp_image=DestroyImage(sharp_image);
560       return((Image *) NULL);
561     }
562   (void) AdaptiveLevelImage(edge_image,"20%,95%",exception);
563   gaussian_image=BlurImage(edge_image,radius,sigma,exception);
564   if (gaussian_image != (Image *) NULL)
565     {
566       edge_image=DestroyImage(edge_image);
567       edge_image=gaussian_image;
568     }
569   (void) AdaptiveLevelImage(edge_image,"10%,95%",exception);
570   /*
571     Create a set of kernels from maximum (radius,sigma) to minimum.
572   */
573   width=GetOptimalKernelWidth2D(radius,sigma);
574   kernel=(MagickRealType **) MagickAssumeAligned(AcquireAlignedMemory((size_t)
575     width,sizeof(*kernel)));
576   if (kernel == (MagickRealType **) NULL)
577     {
578       edge_image=DestroyImage(edge_image);
579       sharp_image=DestroyImage(sharp_image);
580       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
581     }
582   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
583   for (i=0; i < (ssize_t) width; i+=2)
584   {
585     kernel[i]=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
586       (size_t) (width-i),(width-i)*sizeof(**kernel)));
587     if (kernel[i] == (MagickRealType *) NULL)
588       break;
589     normalize=0.0;
590     j=(ssize_t) (width-i)/2;
591     k=0;
592     for (v=(-j); v <= j; v++)
593     {
594       for (u=(-j); u <= j; u++)
595       {
596         kernel[i][k]=(MagickRealType) (-exp(-((double) u*u+v*v)/(2.0*
597           MagickSigma*MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
598         normalize+=kernel[i][k];
599         k++;
600       }
601     }
602     kernel[i][(k-1)/2]=(double) ((-2.0)*normalize);
603     if (sigma < MagickEpsilon)
604       kernel[i][(k-1)/2]=1.0;
605   }
606   if (i < (ssize_t) width)
607     {
608       for (i-=2; i >= 0; i-=2)
609         kernel[i]=(MagickRealType *) RelinquishAlignedMemory(kernel[i]);
610       kernel=(MagickRealType **) RelinquishAlignedMemory(kernel);
611       edge_image=DestroyImage(edge_image);
612       sharp_image=DestroyImage(sharp_image);
613       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
614     }
615   /*
616     Adaptively sharpen image.
617   */
618   status=MagickTrue;
619   progress=0;
620   image_view=AcquireVirtualCacheView(image,exception);
621   edge_view=AcquireVirtualCacheView(edge_image,exception);
622   sharp_view=AcquireAuthenticCacheView(sharp_image,exception);
623 #if defined(MAGICKCORE_OPENMP_SUPPORT)
624   #pragma omp parallel for schedule(static,4) shared(progress,status) \
625     magick_threads(image,sharp_image,sharp_image->rows,1)
626 #endif
627   for (y=0; y < (ssize_t) sharp_image->rows; y++)
628   {
629     register const Quantum
630       *restrict r;
631
632     register Quantum
633       *restrict q;
634
635     register ssize_t
636       x;
637
638     if (status == MagickFalse)
639       continue;
640     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
641     q=QueueCacheViewAuthenticPixels(sharp_view,0,y,sharp_image->columns,1,
642       exception);
643     if ((r == (const Quantum *) NULL) || (q == (Quantum *) NULL))
644       {
645         status=MagickFalse;
646         continue;
647       }
648     for (x=0; x < (ssize_t) sharp_image->columns; x++)
649     {
650       register const Quantum
651         *restrict p;
652
653       register ssize_t
654         i;
655
656       ssize_t
657         center,
658         j;
659
660       j=(ssize_t) ceil((double) width*(1.0-QuantumScale*
661         GetPixelIntensity(edge_image,r))-0.5);
662       if (j < 0)
663         j=0;
664       else
665         if (j > (ssize_t) width)
666           j=(ssize_t) width;
667       if ((j & 0x01) != 0)
668         j--;
669       p=GetCacheViewVirtualPixels(image_view,x-((ssize_t) (width-j)/2L),y-
670         (ssize_t) ((width-j)/2L),width-j,width-j,exception);
671       if (p == (const Quantum *) NULL)
672         break;
673       center=(ssize_t) GetPixelChannels(image)*(width-j)*((width-j)/2L)+
674         GetPixelChannels(image)*((width-j)/2);
675       for (i=0; i < (ssize_t) GetPixelChannels(sharp_image); i++)
676       {
677         double
678           alpha,
679           gamma,
680           pixel;
681
682         PixelChannel
683           channel;
684
685         PixelTrait
686           sharp_traits,
687           traits;
688
689         register const MagickRealType
690           *restrict k;
691
692         register const Quantum
693           *restrict pixels;
694
695         register ssize_t
696           u;
697
698         ssize_t
699           v;
700
701         channel=GetPixelChannelChannel(image,i);
702         traits=GetPixelChannelTraits(image,channel);
703         sharp_traits=GetPixelChannelTraits(sharp_image,channel);
704         if ((traits == UndefinedPixelTrait) ||
705             (sharp_traits == UndefinedPixelTrait))
706           continue;
707         if (((sharp_traits & CopyPixelTrait) != 0) ||
708             (GetPixelReadMask(image,p+center) == 0))
709           {
710             SetPixelChannel(sharp_image,channel,p[center+i],q);
711             continue;
712           }
713         k=kernel[j];
714         pixels=p;
715         pixel=0.0;
716         gamma=0.0;
717         if ((sharp_traits & BlendPixelTrait) == 0)
718           {
719             /*
720               No alpha blending.
721             */
722             for (v=0; v < (ssize_t) (width-j); v++)
723             {
724               for (u=0; u < (ssize_t) (width-j); u++)
725               {
726                 pixel+=(*k)*pixels[i];
727                 gamma+=(*k);
728                 k++;
729                 pixels+=GetPixelChannels(image);
730               }
731             }
732             gamma=PerceptibleReciprocal(gamma);
733             SetPixelChannel(sharp_image,channel,ClampToQuantum(gamma*pixel),q);
734             continue;
735           }
736         /*
737           Alpha blending.
738         */
739         for (v=0; v < (ssize_t) (width-j); v++)
740         {
741           for (u=0; u < (ssize_t) (width-j); u++)
742           {
743             alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
744             pixel+=(*k)*alpha*pixels[i];
745             gamma+=(*k)*alpha;
746             k++;
747             pixels+=GetPixelChannels(image);
748           }
749         }
750         gamma=PerceptibleReciprocal(gamma);
751         SetPixelChannel(sharp_image,channel,ClampToQuantum(gamma*pixel),q);
752       }
753       q+=GetPixelChannels(sharp_image);
754       r+=GetPixelChannels(edge_image);
755     }
756     if (SyncCacheViewAuthenticPixels(sharp_view,exception) == MagickFalse)
757       status=MagickFalse;
758     if (image->progress_monitor != (MagickProgressMonitor) NULL)
759       {
760         MagickBooleanType
761           proceed;
762
763 #if defined(MAGICKCORE_OPENMP_SUPPORT)
764         #pragma omp critical (MagickCore_AdaptiveSharpenImage)
765 #endif
766         proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress++,
767           image->rows);
768         if (proceed == MagickFalse)
769           status=MagickFalse;
770       }
771   }
772   sharp_image->type=image->type;
773   sharp_view=DestroyCacheView(sharp_view);
774   edge_view=DestroyCacheView(edge_view);
775   image_view=DestroyCacheView(image_view);
776   edge_image=DestroyImage(edge_image);
777   for (i=0; i < (ssize_t) width;  i+=2)
778     kernel[i]=(MagickRealType *) RelinquishAlignedMemory(kernel[i]);
779   kernel=(MagickRealType **) RelinquishAlignedMemory(kernel);
780   if (status == MagickFalse)
781     sharp_image=DestroyImage(sharp_image);
782   return(sharp_image);
783 }
784 \f
785 /*
786 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
787 %                                                                             %
788 %                                                                             %
789 %                                                                             %
790 %     B l u r I m a g e                                                       %
791 %                                                                             %
792 %                                                                             %
793 %                                                                             %
794 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
795 %
796 %  BlurImage() blurs an image.  We convolve the image with a Gaussian operator
797 %  of the given radius and standard deviation (sigma).  For reasonable results,
798 %  the radius should be larger than sigma.  Use a radius of 0 and BlurImage()
799 %  selects a suitable radius for you.
800 %
801 %  The format of the BlurImage method is:
802 %
803 %      Image *BlurImage(const Image *image,const double radius,
804 %        const double sigma,ExceptionInfo *exception)
805 %
806 %  A description of each parameter follows:
807 %
808 %    o image: the image.
809 %
810 %    o radius: the radius of the Gaussian, in pixels, not counting the center
811 %      pixel.
812 %
813 %    o sigma: the standard deviation of the Gaussian, in pixels.
814 %
815 %    o exception: return any errors or warnings in this structure.
816 %
817 */
818 MagickExport Image *BlurImage(const Image *image,const double radius,
819   const double sigma,ExceptionInfo *exception)
820 {
821   char
822     geometry[MaxTextExtent];
823
824   KernelInfo
825     *kernel_info;
826
827   Image
828     *blur_image;
829
830   assert(image != (const Image *) NULL);
831   assert(image->signature == MagickSignature);
832   if (image->debug != MagickFalse)
833     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
834   assert(exception != (ExceptionInfo *) NULL);
835   assert(exception->signature == MagickSignature);
836   (void) FormatLocaleString(geometry,MaxTextExtent,
837     "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
838   kernel_info=AcquireKernelInfo(geometry);
839   if (kernel_info == (KernelInfo *) NULL)
840     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
841   blur_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
842     UndefinedCompositeOp,0.0,exception);
843   kernel_info=DestroyKernelInfo(kernel_info);
844   return(blur_image);
845 }
846 \f
847 /*
848 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
849 %                                                                             %
850 %                                                                             %
851 %                                                                             %
852 %     C o n v o l v e I m a g e                                               %
853 %                                                                             %
854 %                                                                             %
855 %                                                                             %
856 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
857 %
858 %  ConvolveImage() applies a custom convolution kernel to the image.
859 %
860 %  The format of the ConvolveImage method is:
861 %
862 %      Image *ConvolveImage(const Image *image,const KernelInfo *kernel,
863 %        ExceptionInfo *exception)
864 %
865 %  A description of each parameter follows:
866 %
867 %    o image: the image.
868 %
869 %    o kernel: the filtering kernel.
870 %
871 %    o exception: return any errors or warnings in this structure.
872 %
873 */
874 MagickExport Image *ConvolveImage(const Image *image,
875   const KernelInfo *kernel_info,ExceptionInfo *exception)
876 {
877   Image
878     *convolve_image;
879
880   convolve_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
881     UndefinedCompositeOp,0.0,exception);
882   return(convolve_image);
883 }
884 \f
885 /*
886 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
887 %                                                                             %
888 %                                                                             %
889 %                                                                             %
890 %     D e s p e c k l e I m a g e                                             %
891 %                                                                             %
892 %                                                                             %
893 %                                                                             %
894 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
895 %
896 %  DespeckleImage() reduces the speckle noise in an image while perserving the
897 %  edges of the original image.  A speckle removing filter uses a complementary %  hulling technique (raising pixels that are darker than their surrounding
898 %  neighbors, then complementarily lowering pixels that are brighter than their
899 %  surrounding neighbors) to reduce the speckle index of that image (reference
900 %  Crimmins speckle removal).
901 %
902 %  The format of the DespeckleImage method is:
903 %
904 %      Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
905 %
906 %  A description of each parameter follows:
907 %
908 %    o image: the image.
909 %
910 %    o exception: return any errors or warnings in this structure.
911 %
912 */
913
914 static void Hull(const Image *image,const ssize_t x_offset,
915   const ssize_t y_offset,const size_t columns,const size_t rows,
916   const int polarity,Quantum *restrict f,Quantum *restrict g)
917 {
918   register Quantum
919     *p,
920     *q,
921     *r,
922     *s;
923
924   ssize_t
925     y;
926
927   assert(image != (const Image *) NULL);
928   assert(image->signature == MagickSignature);
929   if (image->debug != MagickFalse)
930     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
931   assert(f != (Quantum *) NULL);
932   assert(g != (Quantum *) NULL);
933   p=f+(columns+2);
934   q=g+(columns+2);
935   r=p+(y_offset*(columns+2)+x_offset);
936 #if defined(MAGICKCORE_OPENMP_SUPPORT)
937   #pragma omp parallel for schedule(static,4) \
938     magick_threads(image,image,1,1)
939 #endif
940   for (y=0; y < (ssize_t) rows; y++)
941   {
942     MagickRealType
943       v;
944
945     register ssize_t
946       i,
947       x;
948
949     i=(2*y+1)+y*columns;
950     if (polarity > 0)
951       for (x=0; x < (ssize_t) columns; x++)
952       {
953         v=(MagickRealType) p[i];
954         if ((MagickRealType) r[i] >= (v+ScaleCharToQuantum(2)))
955           v+=ScaleCharToQuantum(1);
956         q[i]=(Quantum) v;
957         i++;
958       }
959     else
960       for (x=0; x < (ssize_t) columns; x++)
961       {
962         v=(MagickRealType) p[i];
963         if ((MagickRealType) r[i] <= (v-ScaleCharToQuantum(2)))
964           v-=ScaleCharToQuantum(1);
965         q[i]=(Quantum) v;
966         i++;
967       }
968   }
969   p=f+(columns+2);
970   q=g+(columns+2);
971   r=q+(y_offset*(columns+2)+x_offset);
972   s=q-(y_offset*(columns+2)+x_offset);
973 #if defined(MAGICKCORE_OPENMP_SUPPORT)
974   #pragma omp parallel for schedule(static,4) \
975     magick_threads(image,image,1,1)
976 #endif
977   for (y=0; y < (ssize_t) rows; y++)
978   {
979     register ssize_t
980       i,
981       x;
982
983     MagickRealType
984       v;
985
986     i=(2*y+1)+y*columns;
987     if (polarity > 0)
988       for (x=0; x < (ssize_t) columns; x++)
989       {
990         v=(MagickRealType) q[i];
991         if (((MagickRealType) s[i] >= (v+ScaleCharToQuantum(2))) &&
992             ((MagickRealType) r[i] > v))
993           v+=ScaleCharToQuantum(1);
994         p[i]=(Quantum) v;
995         i++;
996       }
997     else
998       for (x=0; x < (ssize_t) columns; x++)
999       {
1000         v=(MagickRealType) q[i];
1001         if (((MagickRealType) s[i] <= (v-ScaleCharToQuantum(2))) &&
1002             ((MagickRealType) r[i] < v))
1003           v-=ScaleCharToQuantum(1);
1004         p[i]=(Quantum) v;
1005         i++;
1006       }
1007   }
1008 }
1009
1010 MagickExport Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
1011 {
1012 #define DespeckleImageTag  "Despeckle/Image"
1013
1014   CacheView
1015     *despeckle_view,
1016     *image_view;
1017
1018   Image
1019     *despeckle_image;
1020
1021   MagickBooleanType
1022     status;
1023
1024   Quantum
1025     *restrict buffer,
1026     *restrict pixels;
1027
1028   register ssize_t
1029     i;
1030
1031   size_t
1032     length;
1033
1034   static const ssize_t
1035     X[4] = {0, 1, 1,-1},
1036     Y[4] = {1, 0, 1, 1};
1037
1038   /*
1039     Allocate despeckled image.
1040   */
1041   assert(image != (const Image *) NULL);
1042   assert(image->signature == MagickSignature);
1043   if (image->debug != MagickFalse)
1044     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1045   assert(exception != (ExceptionInfo *) NULL);
1046   assert(exception->signature == MagickSignature);
1047   despeckle_image=CloneImage(image,0,0,MagickTrue,exception);
1048   if (despeckle_image == (Image *) NULL)
1049     return((Image *) NULL);
1050   status=SetImageStorageClass(despeckle_image,DirectClass,exception);
1051   if (status == MagickFalse)
1052     {
1053       despeckle_image=DestroyImage(despeckle_image);
1054       return((Image *) NULL);
1055     }
1056   /*
1057     Allocate image buffer.
1058   */
1059   length=(size_t) ((image->columns+2)*(image->rows+2));
1060   pixels=(Quantum *) AcquireQuantumMemory(length,sizeof(*pixels));
1061   buffer=(Quantum *) AcquireQuantumMemory(length,sizeof(*buffer));
1062   if ((pixels == (Quantum *) NULL) || (buffer == (Quantum *) NULL))
1063     {
1064       if (buffer != (Quantum *) NULL)
1065         buffer=(Quantum *) RelinquishMagickMemory(buffer);
1066       if (pixels != (Quantum *) NULL)
1067         pixels=(Quantum *) RelinquishMagickMemory(pixels);
1068       despeckle_image=DestroyImage(despeckle_image);
1069       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1070     }
1071   /*
1072     Reduce speckle in the image.
1073   */
1074   status=MagickTrue;
1075   image_view=AcquireVirtualCacheView(image,exception);
1076   despeckle_view=AcquireAuthenticCacheView(despeckle_image,exception);
1077   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1078   {
1079     PixelChannel
1080        channel;
1081
1082     PixelTrait
1083       despeckle_traits,
1084       traits;
1085
1086     register ssize_t
1087       k,
1088       x;
1089
1090     ssize_t
1091       j,
1092       y;
1093
1094     if (status == MagickFalse)
1095       continue;
1096     channel=GetPixelChannelChannel(image,i);
1097     traits=GetPixelChannelTraits(image,channel);
1098     despeckle_traits=GetPixelChannelTraits(despeckle_image,channel);
1099     if ((traits == UndefinedPixelTrait) ||
1100         (despeckle_traits == UndefinedPixelTrait))
1101       continue;
1102     if ((despeckle_traits & CopyPixelTrait) != 0)
1103       continue;
1104     (void) ResetMagickMemory(pixels,0,length*sizeof(*pixels));
1105     j=(ssize_t) image->columns+2;
1106     for (y=0; y < (ssize_t) image->rows; y++)
1107     {
1108       register const Quantum
1109         *restrict p;
1110
1111       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1112       if (p == (const Quantum *) NULL)
1113         {
1114           status=MagickFalse;
1115           continue;
1116         }
1117       j++;
1118       for (x=0; x < (ssize_t) image->columns; x++)
1119       {
1120         pixels[j++]=p[i];
1121         p+=GetPixelChannels(image);
1122       }
1123       j++;
1124     }
1125     (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
1126     for (k=0; k < 4; k++)
1127     {
1128       Hull(image,X[k],Y[k],image->columns,image->rows,1,pixels,buffer);
1129       Hull(image,-X[k],-Y[k],image->columns,image->rows,1,pixels,buffer);
1130       Hull(image,-X[k],-Y[k],image->columns,image->rows,-1,pixels,buffer);
1131       Hull(image,X[k],Y[k],image->columns,image->rows,-1,pixels,buffer);
1132     }
1133     j=(ssize_t) image->columns+2;
1134     for (y=0; y < (ssize_t) image->rows; y++)
1135     {
1136       MagickBooleanType
1137         sync;
1138
1139       register Quantum
1140         *restrict q;
1141
1142       q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
1143         1,exception);
1144       if (q == (Quantum *) NULL)
1145         {
1146           status=MagickFalse;
1147           continue;
1148         }
1149       j++;
1150       for (x=0; x < (ssize_t) image->columns; x++)
1151       {
1152         SetPixelChannel(despeckle_image,channel,pixels[j++],q);
1153         q+=GetPixelChannels(despeckle_image);
1154       }
1155       sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
1156       if (sync == MagickFalse)
1157         status=MagickFalse;
1158       j++;
1159     }
1160     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1161       {
1162         MagickBooleanType
1163           proceed;
1164
1165         proceed=SetImageProgress(image,DespeckleImageTag,(MagickOffsetType) i,
1166           GetPixelChannels(image));
1167         if (proceed == MagickFalse)
1168           status=MagickFalse;
1169       }
1170   }
1171   despeckle_view=DestroyCacheView(despeckle_view);
1172   image_view=DestroyCacheView(image_view);
1173   buffer=(Quantum *) RelinquishMagickMemory(buffer);
1174   pixels=(Quantum *) RelinquishMagickMemory(pixels);
1175   despeckle_image->type=image->type;
1176   if (status == MagickFalse)
1177     despeckle_image=DestroyImage(despeckle_image);
1178   return(despeckle_image);
1179 }
1180 \f
1181 /*
1182 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1183 %                                                                             %
1184 %                                                                             %
1185 %                                                                             %
1186 %     E d g e I m a g e                                                       %
1187 %                                                                             %
1188 %                                                                             %
1189 %                                                                             %
1190 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1191 %
1192 %  EdgeImage() finds edges in an image.  Radius defines the radius of the
1193 %  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
1194 %  radius for you.
1195 %
1196 %  The format of the EdgeImage method is:
1197 %
1198 %      Image *EdgeImage(const Image *image,const double radius,
1199 %        ExceptionInfo *exception)
1200 %
1201 %  A description of each parameter follows:
1202 %
1203 %    o image: the image.
1204 %
1205 %    o radius: the radius of the pixel neighborhood.
1206 %
1207 %    o exception: return any errors or warnings in this structure.
1208 %
1209 */
1210 MagickExport Image *EdgeImage(const Image *image,const double radius,
1211   ExceptionInfo *exception)
1212 {
1213   Image
1214     *edge_image;
1215
1216   KernelInfo
1217     *kernel_info;
1218
1219   register ssize_t
1220     i;
1221
1222   size_t
1223     width;
1224
1225   assert(image != (const Image *) NULL);
1226   assert(image->signature == MagickSignature);
1227   if (image->debug != MagickFalse)
1228     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1229   assert(exception != (ExceptionInfo *) NULL);
1230   assert(exception->signature == MagickSignature);
1231   width=GetOptimalKernelWidth1D(radius,0.5);
1232   kernel_info=AcquireKernelInfo((const char *) NULL);
1233   if (kernel_info == (KernelInfo *) NULL)
1234     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1235   (void) ResetMagickMemory(kernel_info,0,sizeof(*kernel_info));
1236   kernel_info->width=width;
1237   kernel_info->height=width;
1238   kernel_info->x=(ssize_t) (kernel_info->width-1)/2;
1239   kernel_info->y=(ssize_t) (kernel_info->height-1)/2;
1240   kernel_info->signature=MagickSignature;
1241   kernel_info->values=(MagickRealType *) MagickAssumeAligned(
1242     AcquireAlignedMemory(kernel_info->width,kernel_info->height*
1243     sizeof(*kernel_info->values)));
1244   if (kernel_info->values == (MagickRealType *) NULL)
1245     {
1246       kernel_info=DestroyKernelInfo(kernel_info);
1247       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1248     }
1249   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
1250     kernel_info->values[i]=(-1.0);
1251   kernel_info->values[i/2]=(double) kernel_info->width*kernel_info->height-1.0;
1252   edge_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
1253     UndefinedCompositeOp,0.0,exception);
1254   kernel_info=DestroyKernelInfo(kernel_info);
1255   return(edge_image);
1256 }
1257 \f
1258 /*
1259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1260 %                                                                             %
1261 %                                                                             %
1262 %                                                                             %
1263 %     E m b o s s I m a g e                                                   %
1264 %                                                                             %
1265 %                                                                             %
1266 %                                                                             %
1267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1268 %
1269 %  EmbossImage() returns a grayscale image with a three-dimensional effect.
1270 %  We convolve the image with a Gaussian operator of the given radius and
1271 %  standard deviation (sigma).  For reasonable results, radius should be
1272 %  larger than sigma.  Use a radius of 0 and Emboss() selects a suitable
1273 %  radius for you.
1274 %
1275 %  The format of the EmbossImage method is:
1276 %
1277 %      Image *EmbossImage(const Image *image,const double radius,
1278 %        const double sigma,ExceptionInfo *exception)
1279 %
1280 %  A description of each parameter follows:
1281 %
1282 %    o image: the image.
1283 %
1284 %    o radius: the radius of the pixel neighborhood.
1285 %
1286 %    o sigma: the standard deviation of the Gaussian, in pixels.
1287 %
1288 %    o exception: return any errors or warnings in this structure.
1289 %
1290 */
1291 MagickExport Image *EmbossImage(const Image *image,const double radius,
1292   const double sigma,ExceptionInfo *exception)
1293 {
1294   double
1295     gamma,
1296     normalize;
1297
1298   Image
1299     *emboss_image;
1300
1301   KernelInfo
1302     *kernel_info;
1303
1304   register ssize_t
1305     i;
1306
1307   size_t
1308     width;
1309
1310   ssize_t
1311     j,
1312     k,
1313     u,
1314     v;
1315
1316   assert(image != (const Image *) NULL);
1317   assert(image->signature == MagickSignature);
1318   if (image->debug != MagickFalse)
1319     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1320   assert(exception != (ExceptionInfo *) NULL);
1321   assert(exception->signature == MagickSignature);
1322   width=GetOptimalKernelWidth1D(radius,sigma);
1323   kernel_info=AcquireKernelInfo((const char *) NULL);
1324   if (kernel_info == (KernelInfo *) NULL)
1325     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1326   kernel_info->width=width;
1327   kernel_info->height=width;
1328   kernel_info->x=(ssize_t) (width-1)/2;
1329   kernel_info->y=(ssize_t) (width-1)/2;
1330   kernel_info->values=(MagickRealType *) MagickAssumeAligned(
1331     AcquireAlignedMemory(kernel_info->width,kernel_info->width*
1332     sizeof(*kernel_info->values)));
1333   if (kernel_info->values == (MagickRealType *) NULL)
1334     {
1335       kernel_info=DestroyKernelInfo(kernel_info);
1336       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1337     }
1338   j=(ssize_t) (kernel_info->width-1)/2;
1339   k=j;
1340   i=0;
1341   for (v=(-j); v <= j; v++)
1342   {
1343     for (u=(-j); u <= j; u++)
1344     {
1345       kernel_info->values[i]=(MagickRealType) (((u < 0) || (v < 0) ? -8.0 :
1346         8.0)*exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
1347         (2.0*MagickPI*MagickSigma*MagickSigma));
1348       if (u == k)
1349         kernel_info->values[i]=v == k ? 1.0 : 0.0;
1350       i++;
1351     }
1352     k--;
1353   }
1354   normalize=0.0;
1355   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
1356     normalize+=kernel_info->values[i];
1357   gamma=PerceptibleReciprocal(normalize);
1358   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
1359     kernel_info->values[i]*=gamma;
1360   emboss_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
1361     UndefinedCompositeOp,0.0,exception);
1362   kernel_info=DestroyKernelInfo(kernel_info);
1363   if (emboss_image != (Image *) NULL)
1364     (void) EqualizeImage(emboss_image,exception);
1365   return(emboss_image);
1366 }
1367 \f
1368 /*
1369 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1370 %                                                                             %
1371 %                                                                             %
1372 %                                                                             %
1373 %     G a u s s i a n B l u r I m a g e                                       %
1374 %                                                                             %
1375 %                                                                             %
1376 %                                                                             %
1377 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1378 %
1379 %  GaussianBlurImage() blurs an image.  We convolve the image with a
1380 %  Gaussian operator of the given radius and standard deviation (sigma).
1381 %  For reasonable results, the radius should be larger than sigma.  Use a
1382 %  radius of 0 and GaussianBlurImage() selects a suitable radius for you
1383 %
1384 %  The format of the GaussianBlurImage method is:
1385 %
1386 %      Image *GaussianBlurImage(const Image *image,onst double radius,
1387 %        const double sigma,ExceptionInfo *exception)
1388 %
1389 %  A description of each parameter follows:
1390 %
1391 %    o image: the image.
1392 %
1393 %    o radius: the radius of the Gaussian, in pixels, not counting the center
1394 %      pixel.
1395 %
1396 %    o sigma: the standard deviation of the Gaussian, in pixels.
1397 %
1398 %    o exception: return any errors or warnings in this structure.
1399 %
1400 */
1401 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
1402   const double sigma,ExceptionInfo *exception)
1403 {
1404   char
1405     geometry[MaxTextExtent];
1406
1407   KernelInfo
1408     *kernel_info;
1409
1410   Image
1411     *blur_image;
1412
1413   assert(image != (const Image *) NULL);
1414   assert(image->signature == MagickSignature);
1415   if (image->debug != MagickFalse)
1416     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1417   assert(exception != (ExceptionInfo *) NULL);
1418   assert(exception->signature == MagickSignature);
1419   (void) FormatLocaleString(geometry,MaxTextExtent,"gaussian:%.20gx%.20g",
1420     radius,sigma);
1421   kernel_info=AcquireKernelInfo(geometry);
1422   if (kernel_info == (KernelInfo *) NULL)
1423     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1424   blur_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
1425     UndefinedCompositeOp,0.0,exception);
1426   kernel_info=DestroyKernelInfo(kernel_info);
1427   return(blur_image);
1428 }
1429 \f
1430 /*
1431 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1432 %                                                                             %
1433 %                                                                             %
1434 %                                                                             %
1435 %     M o t i o n B l u r I m a g e                                           %
1436 %                                                                             %
1437 %                                                                             %
1438 %                                                                             %
1439 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1440 %
1441 %  MotionBlurImage() simulates motion blur.  We convolve the image with a
1442 %  Gaussian operator of the given radius and standard deviation (sigma).
1443 %  For reasonable results, radius should be larger than sigma.  Use a
1444 %  radius of 0 and MotionBlurImage() selects a suitable radius for you.
1445 %  Angle gives the angle of the blurring motion.
1446 %
1447 %  Andrew Protano contributed this effect.
1448 %
1449 %  The format of the MotionBlurImage method is:
1450 %
1451 %    Image *MotionBlurImage(const Image *image,const double radius,
1452 %      const double sigma,const double angle,ExceptionInfo *exception)
1453 %
1454 %  A description of each parameter follows:
1455 %
1456 %    o image: the image.
1457 %
1458 %    o radius: the radius of the Gaussian, in pixels, not counting
1459 %      the center pixel.
1460 %
1461 %    o sigma: the standard deviation of the Gaussian, in pixels.
1462 %
1463 %    o angle: Apply the effect along this angle.
1464 %
1465 %    o exception: return any errors or warnings in this structure.
1466 %
1467 */
1468
1469 static MagickRealType *GetMotionBlurKernel(const size_t width,
1470   const double sigma)
1471 {
1472   MagickRealType
1473     *kernel,
1474     normalize;
1475
1476   register ssize_t
1477     i;
1478
1479   /*
1480    Generate a 1-D convolution kernel.
1481   */
1482   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
1483   kernel=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory((size_t)
1484     width,sizeof(*kernel)));
1485   if (kernel == (MagickRealType *) NULL)
1486     return(kernel);
1487   normalize=0.0;
1488   for (i=0; i < (ssize_t) width; i++)
1489   {
1490     kernel[i]=(MagickRealType) (exp((-((double) i*i)/(double) (2.0*MagickSigma*
1491       MagickSigma)))/(MagickSQ2PI*MagickSigma));
1492     normalize+=kernel[i];
1493   }
1494   for (i=0; i < (ssize_t) width; i++)
1495     kernel[i]/=normalize;
1496   return(kernel);
1497 }
1498
1499 MagickExport Image *MotionBlurImage(const Image *image,const double radius,
1500   const double sigma,const double angle,ExceptionInfo *exception)
1501 {
1502 #define BlurImageTag  "Blur/Image"
1503
1504   CacheView
1505     *blur_view,
1506     *image_view,
1507     *motion_view;
1508
1509   Image
1510     *blur_image;
1511
1512   MagickBooleanType
1513     status;
1514
1515   MagickOffsetType
1516     progress;
1517
1518   MagickRealType
1519     *kernel;
1520
1521   OffsetInfo
1522     *offset;
1523
1524   PointInfo
1525     point;
1526
1527   register ssize_t
1528     i;
1529
1530   size_t
1531     width;
1532
1533   ssize_t
1534     y;
1535
1536   assert(image != (Image *) NULL);
1537   assert(image->signature == MagickSignature);
1538   if (image->debug != MagickFalse)
1539     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1540   assert(exception != (ExceptionInfo *) NULL);
1541   width=GetOptimalKernelWidth1D(radius,sigma);
1542   kernel=GetMotionBlurKernel(width,sigma);
1543   if (kernel == (MagickRealType *) NULL)
1544     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1545   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
1546   if (offset == (OffsetInfo *) NULL)
1547     {
1548       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
1549       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1550     }
1551   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
1552   if (blur_image == (Image *) NULL)
1553     {
1554       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
1555       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
1556       return((Image *) NULL);
1557     }
1558   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
1559     {
1560       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
1561       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
1562       blur_image=DestroyImage(blur_image);
1563       return((Image *) NULL);
1564     }
1565   point.x=(double) width*sin(DegreesToRadians(angle));
1566   point.y=(double) width*cos(DegreesToRadians(angle));
1567   for (i=0; i < (ssize_t) width; i++)
1568   {
1569     offset[i].x=(ssize_t) ceil((double) (i*point.y)/hypot(point.x,point.y)-0.5);
1570     offset[i].y=(ssize_t) ceil((double) (i*point.x)/hypot(point.x,point.y)-0.5);
1571   }
1572   /*
1573     Motion blur image.
1574   */
1575   status=MagickTrue;
1576   progress=0;
1577   image_view=AcquireVirtualCacheView(image,exception);
1578   motion_view=AcquireVirtualCacheView(image,exception);
1579   blur_view=AcquireAuthenticCacheView(blur_image,exception);
1580 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1581   #pragma omp parallel for schedule(static,4) shared(progress,status) \
1582     magick_threads(image,blur_image,image->rows,1)
1583 #endif
1584   for (y=0; y < (ssize_t) image->rows; y++)
1585   {
1586     register const Quantum
1587       *restrict p;
1588
1589     register Quantum
1590       *restrict q;
1591
1592     register ssize_t
1593       x;
1594
1595     if (status == MagickFalse)
1596       continue;
1597     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1598     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
1599       exception);
1600     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
1601       {
1602         status=MagickFalse;
1603         continue;
1604       }
1605     for (x=0; x < (ssize_t) image->columns; x++)
1606     {
1607       register ssize_t
1608         i;
1609
1610       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
1611       {
1612         double
1613           alpha,
1614           gamma,
1615           pixel;
1616
1617         PixelChannel
1618           channel;
1619
1620         PixelTrait
1621           blur_traits,
1622           traits;
1623
1624         register const Quantum
1625           *restrict r;
1626
1627         register MagickRealType
1628           *restrict k;
1629
1630         register ssize_t
1631           j;
1632
1633         channel=GetPixelChannelChannel(image,i);
1634         traits=GetPixelChannelTraits(image,channel);
1635         blur_traits=GetPixelChannelTraits(blur_image,channel);
1636         if ((traits == UndefinedPixelTrait) ||
1637             (blur_traits == UndefinedPixelTrait))
1638           continue;
1639         if (((blur_traits & CopyPixelTrait) != 0) ||
1640             (GetPixelReadMask(image,p) == 0))
1641           {
1642             SetPixelChannel(blur_image,channel,p[i],q);
1643             continue;
1644           }
1645         k=kernel;
1646         pixel=0.0;
1647         if ((blur_traits & BlendPixelTrait) == 0)
1648           {
1649             for (j=0; j < (ssize_t) width; j++)
1650             {
1651               r=GetCacheViewVirtualPixels(motion_view,x+offset[j].x,y+
1652                 offset[j].y,1,1,exception);
1653               if (r == (const Quantum *) NULL)
1654                 {
1655                   status=MagickFalse;
1656                   continue;
1657                 }
1658               pixel+=(*k)*r[i];
1659               k++;
1660             }
1661             SetPixelChannel(blur_image,channel,ClampToQuantum(pixel),q);
1662             continue;
1663           }
1664         alpha=0.0;
1665         gamma=0.0;
1666         for (j=0; j < (ssize_t) width; j++)
1667         {
1668           r=GetCacheViewVirtualPixels(motion_view,x+offset[j].x,y+offset[j].y,1,
1669             1,exception);
1670           if (r == (const Quantum *) NULL)
1671             {
1672               status=MagickFalse;
1673               continue;
1674             }
1675           alpha=(double) (QuantumScale*GetPixelAlpha(image,r));
1676           pixel+=(*k)*alpha*r[i];
1677           gamma+=(*k)*alpha;
1678           k++;
1679         }
1680         gamma=PerceptibleReciprocal(gamma);
1681         SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
1682       }
1683       p+=GetPixelChannels(image);
1684       q+=GetPixelChannels(blur_image);
1685     }
1686     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
1687       status=MagickFalse;
1688     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1689       {
1690         MagickBooleanType
1691           proceed;
1692
1693 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1694         #pragma omp critical (MagickCore_MotionBlurImage)
1695 #endif
1696         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
1697         if (proceed == MagickFalse)
1698           status=MagickFalse;
1699       }
1700   }
1701   blur_view=DestroyCacheView(blur_view);
1702   motion_view=DestroyCacheView(motion_view);
1703   image_view=DestroyCacheView(image_view);
1704   kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
1705   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
1706   if (status == MagickFalse)
1707     blur_image=DestroyImage(blur_image);
1708   return(blur_image);
1709 }
1710 \f
1711 /*
1712 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1713 %                                                                             %
1714 %                                                                             %
1715 %                                                                             %
1716 %     P r e v i e w I m a g e                                                 %
1717 %                                                                             %
1718 %                                                                             %
1719 %                                                                             %
1720 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1721 %
1722 %  PreviewImage() tiles 9 thumbnails of the specified image with an image
1723 %  processing operation applied with varying parameters.  This may be helpful
1724 %  pin-pointing an appropriate parameter for a particular image processing
1725 %  operation.
1726 %
1727 %  The format of the PreviewImages method is:
1728 %
1729 %      Image *PreviewImages(const Image *image,const PreviewType preview,
1730 %        ExceptionInfo *exception)
1731 %
1732 %  A description of each parameter follows:
1733 %
1734 %    o image: the image.
1735 %
1736 %    o preview: the image processing operation.
1737 %
1738 %    o exception: return any errors or warnings in this structure.
1739 %
1740 */
1741 MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
1742   ExceptionInfo *exception)
1743 {
1744 #define NumberTiles  9
1745 #define PreviewImageTag  "Preview/Image"
1746 #define DefaultPreviewGeometry  "204x204+10+10"
1747
1748   char
1749     factor[MaxTextExtent],
1750     label[MaxTextExtent];
1751
1752   double
1753     degrees,
1754     gamma,
1755     percentage,
1756     radius,
1757     sigma,
1758     threshold;
1759
1760   extern const char
1761     DefaultTileFrame[];
1762
1763   Image
1764     *images,
1765     *montage_image,
1766     *preview_image,
1767     *thumbnail;
1768
1769   ImageInfo
1770     *preview_info;
1771
1772   MagickBooleanType
1773     proceed;
1774
1775   MontageInfo
1776     *montage_info;
1777
1778   QuantizeInfo
1779     quantize_info;
1780
1781   RectangleInfo
1782     geometry;
1783
1784   register ssize_t
1785     i,
1786     x;
1787
1788   size_t
1789     colors;
1790
1791   ssize_t
1792     y;
1793
1794   /*
1795     Open output image file.
1796   */
1797   assert(image != (Image *) NULL);
1798   assert(image->signature == MagickSignature);
1799   if (image->debug != MagickFalse)
1800     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1801   colors=2;
1802   degrees=0.0;
1803   gamma=(-0.2f);
1804   preview_info=AcquireImageInfo();
1805   SetGeometry(image,&geometry);
1806   (void) ParseMetaGeometry(DefaultPreviewGeometry,&geometry.x,&geometry.y,
1807     &geometry.width,&geometry.height);
1808   images=NewImageList();
1809   percentage=12.5;
1810   GetQuantizeInfo(&quantize_info);
1811   radius=0.0;
1812   sigma=1.0;
1813   threshold=0.0;
1814   x=0;
1815   y=0;
1816   for (i=0; i < NumberTiles; i++)
1817   {
1818     thumbnail=ThumbnailImage(image,geometry.width,geometry.height,exception);
1819     if (thumbnail == (Image *) NULL)
1820       break;
1821     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
1822       (void *) NULL);
1823     (void) SetImageProperty(thumbnail,"label",DefaultTileLabel,exception);
1824     if (i == (NumberTiles/2))
1825       {
1826         (void) QueryColorCompliance("#dfdfdf",AllCompliance,
1827           &thumbnail->matte_color,exception);
1828         AppendImageToList(&images,thumbnail);
1829         continue;
1830       }
1831     switch (preview)
1832     {
1833       case RotatePreview:
1834       {
1835         degrees+=45.0;
1836         preview_image=RotateImage(thumbnail,degrees,exception);
1837         (void) FormatLocaleString(label,MaxTextExtent,"rotate %g",degrees);
1838         break;
1839       }
1840       case ShearPreview:
1841       {
1842         degrees+=5.0;
1843         preview_image=ShearImage(thumbnail,degrees,degrees,exception);
1844         (void) FormatLocaleString(label,MaxTextExtent,"shear %gx%g",degrees,
1845           2.0*degrees);
1846         break;
1847       }
1848       case RollPreview:
1849       {
1850         x=(ssize_t) ((i+1)*thumbnail->columns)/NumberTiles;
1851         y=(ssize_t) ((i+1)*thumbnail->rows)/NumberTiles;
1852         preview_image=RollImage(thumbnail,x,y,exception);
1853         (void) FormatLocaleString(label,MaxTextExtent,"roll %+.20gx%+.20g",
1854           (double) x,(double) y);
1855         break;
1856       }
1857       case HuePreview:
1858       {
1859         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
1860         if (preview_image == (Image *) NULL)
1861           break;
1862         (void) FormatLocaleString(factor,MaxTextExtent,"100,100,%g",2.0*
1863           percentage);
1864         (void) ModulateImage(preview_image,factor,exception);
1865         (void) FormatLocaleString(label,MaxTextExtent,"modulate %s",factor);
1866         break;
1867       }
1868       case SaturationPreview:
1869       {
1870         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
1871         if (preview_image == (Image *) NULL)
1872           break;
1873         (void) FormatLocaleString(factor,MaxTextExtent,"100,%g",2.0*percentage);
1874         (void) ModulateImage(preview_image,factor,exception);
1875         (void) FormatLocaleString(label,MaxTextExtent,"modulate %s",factor);
1876         break;
1877       }
1878       case BrightnessPreview:
1879       {
1880         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
1881         if (preview_image == (Image *) NULL)
1882           break;
1883         (void) FormatLocaleString(factor,MaxTextExtent,"%g",2.0*percentage);
1884         (void) ModulateImage(preview_image,factor,exception);
1885         (void) FormatLocaleString(label,MaxTextExtent,"modulate %s",factor);
1886         break;
1887       }
1888       case GammaPreview:
1889       default:
1890       {
1891         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
1892         if (preview_image == (Image *) NULL)
1893           break;
1894         gamma+=0.4f;
1895         (void) GammaImage(preview_image,gamma,exception);
1896         (void) FormatLocaleString(label,MaxTextExtent,"gamma %g",gamma);
1897         break;
1898       }
1899       case SpiffPreview:
1900       {
1901         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
1902         if (preview_image != (Image *) NULL)
1903           for (x=0; x < i; x++)
1904             (void) ContrastImage(preview_image,MagickTrue,exception);
1905         (void) FormatLocaleString(label,MaxTextExtent,"contrast (%.20g)",
1906           (double) i+1);
1907         break;
1908       }
1909       case DullPreview:
1910       {
1911         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
1912         if (preview_image == (Image *) NULL)
1913           break;
1914         for (x=0; x < i; x++)
1915           (void) ContrastImage(preview_image,MagickFalse,exception);
1916         (void) FormatLocaleString(label,MaxTextExtent,"+contrast (%.20g)",
1917           (double) i+1);
1918         break;
1919       }
1920       case GrayscalePreview:
1921       {
1922         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
1923         if (preview_image == (Image *) NULL)
1924           break;
1925         colors<<=1;
1926         quantize_info.number_colors=colors;
1927         quantize_info.colorspace=GRAYColorspace;
1928         (void) QuantizeImage(&quantize_info,preview_image,exception);
1929         (void) FormatLocaleString(label,MaxTextExtent,
1930           "-colorspace gray -colors %.20g",(double) colors);
1931         break;
1932       }
1933       case QuantizePreview:
1934       {
1935         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
1936         if (preview_image == (Image *) NULL)
1937           break;
1938         colors<<=1;
1939         quantize_info.number_colors=colors;
1940         (void) QuantizeImage(&quantize_info,preview_image,exception);
1941         (void) FormatLocaleString(label,MaxTextExtent,"colors %.20g",(double)
1942           colors);
1943         break;
1944       }
1945       case DespecklePreview:
1946       {
1947         for (x=0; x < (i-1); x++)
1948         {
1949           preview_image=DespeckleImage(thumbnail,exception);
1950           if (preview_image == (Image *) NULL)
1951             break;
1952           thumbnail=DestroyImage(thumbnail);
1953           thumbnail=preview_image;
1954         }
1955         preview_image=DespeckleImage(thumbnail,exception);
1956         if (preview_image == (Image *) NULL)
1957           break;
1958         (void) FormatLocaleString(label,MaxTextExtent,"despeckle (%.20g)",
1959           (double) i+1);
1960         break;
1961       }
1962       case ReduceNoisePreview:
1963       {
1964         preview_image=StatisticImage(thumbnail,NonpeakStatistic,(size_t) radius,
1965           (size_t) radius,exception);
1966         (void) FormatLocaleString(label,MaxTextExtent,"noise %g",radius);
1967         break;
1968       }
1969       case AddNoisePreview:
1970       {
1971         switch ((int) i)
1972         {
1973           case 0:
1974           {
1975             (void) CopyMagickString(factor,"uniform",MaxTextExtent);
1976             break;
1977           }
1978           case 1:
1979           {
1980             (void) CopyMagickString(factor,"gaussian",MaxTextExtent);
1981             break;
1982           }
1983           case 2:
1984           {
1985             (void) CopyMagickString(factor,"multiplicative",MaxTextExtent);
1986             break;
1987           }
1988           case 3:
1989           {
1990             (void) CopyMagickString(factor,"impulse",MaxTextExtent);
1991             break;
1992           }
1993           case 4:
1994           {
1995             (void) CopyMagickString(factor,"laplacian",MaxTextExtent);
1996             break;
1997           }
1998           case 5:
1999           {
2000             (void) CopyMagickString(factor,"Poisson",MaxTextExtent);
2001             break;
2002           }
2003           default:
2004           {
2005             (void) CopyMagickString(thumbnail->magick,"NULL",MaxTextExtent);
2006             break;
2007           }
2008         }
2009         preview_image=StatisticImage(thumbnail,NonpeakStatistic,(size_t) i,
2010           (size_t) i,exception);
2011         (void) FormatLocaleString(label,MaxTextExtent,"+noise %s",factor);
2012         break;
2013       }
2014       case SharpenPreview:
2015       {
2016         preview_image=SharpenImage(thumbnail,radius,sigma,exception);
2017         (void) FormatLocaleString(label,MaxTextExtent,"sharpen %gx%g",radius,
2018           sigma);
2019         break;
2020       }
2021       case BlurPreview:
2022       {
2023         preview_image=BlurImage(thumbnail,radius,sigma,exception);
2024         (void) FormatLocaleString(label,MaxTextExtent,"blur %gx%g",radius,
2025           sigma);
2026         break;
2027       }
2028       case ThresholdPreview:
2029       {
2030         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2031         if (preview_image == (Image *) NULL)
2032           break;
2033         (void) BilevelImage(thumbnail,(double) (percentage*((double)
2034           QuantumRange+1.0))/100.0,exception);
2035         (void) FormatLocaleString(label,MaxTextExtent,"threshold %g",(double)
2036           (percentage*((double) QuantumRange+1.0))/100.0);
2037         break;
2038       }
2039       case EdgeDetectPreview:
2040       {
2041         preview_image=EdgeImage(thumbnail,radius,exception);
2042         (void) FormatLocaleString(label,MaxTextExtent,"edge %g",radius);
2043         break;
2044       }
2045       case SpreadPreview:
2046       {
2047         preview_image=SpreadImage(thumbnail,radius,thumbnail->interpolate,
2048           exception);
2049         (void) FormatLocaleString(label,MaxTextExtent,"spread %g",radius+0.5);
2050         break;
2051       }
2052       case SolarizePreview:
2053       {
2054         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2055         if (preview_image == (Image *) NULL)
2056           break;
2057         (void) SolarizeImage(preview_image,(double) QuantumRange*percentage/
2058           100.0,exception);
2059         (void) FormatLocaleString(label,MaxTextExtent,"solarize %g",
2060           (QuantumRange*percentage)/100.0);
2061         break;
2062       }
2063       case ShadePreview:
2064       {
2065         degrees+=10.0;
2066         preview_image=ShadeImage(thumbnail,MagickTrue,degrees,degrees,
2067           exception);
2068         (void) FormatLocaleString(label,MaxTextExtent,"shade %gx%g",degrees,
2069           degrees);
2070         break;
2071       }
2072       case RaisePreview:
2073       {
2074         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2075         if (preview_image == (Image *) NULL)
2076           break;
2077         geometry.width=(size_t) (2*i+2);
2078         geometry.height=(size_t) (2*i+2);
2079         geometry.x=(i-1)/2;
2080         geometry.y=(i-1)/2;
2081         (void) RaiseImage(preview_image,&geometry,MagickTrue,exception);
2082         (void) FormatLocaleString(label,MaxTextExtent,
2083           "raise %.20gx%.20g%+.20g%+.20g",(double) geometry.width,(double)
2084           geometry.height,(double) geometry.x,(double) geometry.y);
2085         break;
2086       }
2087       case SegmentPreview:
2088       {
2089         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2090         if (preview_image == (Image *) NULL)
2091           break;
2092         threshold+=0.4f;
2093         (void) SegmentImage(preview_image,sRGBColorspace,MagickFalse,threshold,
2094           threshold,exception);
2095         (void) FormatLocaleString(label,MaxTextExtent,"segment %gx%g",
2096           threshold,threshold);
2097         break;
2098       }
2099       case SwirlPreview:
2100       {
2101         preview_image=SwirlImage(thumbnail,degrees,image->interpolate,
2102           exception);
2103         (void) FormatLocaleString(label,MaxTextExtent,"swirl %g",degrees);
2104         degrees+=45.0;
2105         break;
2106       }
2107       case ImplodePreview:
2108       {
2109         degrees+=0.1f;
2110         preview_image=ImplodeImage(thumbnail,degrees,image->interpolate,
2111           exception);
2112         (void) FormatLocaleString(label,MaxTextExtent,"implode %g",degrees);
2113         break;
2114       }
2115       case WavePreview:
2116       {
2117         degrees+=5.0f;
2118         preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,
2119           image->interpolate,exception);
2120         (void) FormatLocaleString(label,MaxTextExtent,"wave %gx%g",0.5*degrees,
2121           2.0*degrees);
2122         break;
2123       }
2124       case OilPaintPreview:
2125       {
2126         preview_image=OilPaintImage(thumbnail,(double) radius,(double) sigma,
2127           exception);
2128         (void) FormatLocaleString(label,MaxTextExtent,"charcoal %gx%g",radius,
2129           sigma);
2130         break;
2131       }
2132       case CharcoalDrawingPreview:
2133       {
2134         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
2135           exception);
2136         (void) FormatLocaleString(label,MaxTextExtent,"charcoal %gx%g",radius,
2137           sigma);
2138         break;
2139       }
2140       case JPEGPreview:
2141       {
2142         char
2143           filename[MaxTextExtent];
2144
2145         int
2146           file;
2147
2148         MagickBooleanType
2149           status;
2150
2151         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
2152         if (preview_image == (Image *) NULL)
2153           break;
2154         preview_info->quality=(size_t) percentage;
2155         (void) FormatLocaleString(factor,MaxTextExtent,"%.20g",(double)
2156           preview_info->quality);
2157         file=AcquireUniqueFileResource(filename);
2158         if (file != -1)
2159           file=close(file)-1;
2160         (void) FormatLocaleString(preview_image->filename,MaxTextExtent,
2161           "jpeg:%s",filename);
2162         status=WriteImage(preview_info,preview_image,exception);
2163         if (status != MagickFalse)
2164           {
2165             Image
2166               *quality_image;
2167
2168             (void) CopyMagickString(preview_info->filename,
2169               preview_image->filename,MaxTextExtent);
2170             quality_image=ReadImage(preview_info,exception);
2171             if (quality_image != (Image *) NULL)
2172               {
2173                 preview_image=DestroyImage(preview_image);
2174                 preview_image=quality_image;
2175               }
2176           }
2177         (void) RelinquishUniqueFileResource(preview_image->filename);
2178         if ((GetBlobSize(preview_image)/1024) >= 1024)
2179           (void) FormatLocaleString(label,MaxTextExtent,"quality %s\n%gmb ",
2180             factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
2181             1024.0/1024.0);
2182         else
2183           if (GetBlobSize(preview_image) >= 1024)
2184             (void) FormatLocaleString(label,MaxTextExtent,
2185               "quality %s\n%gkb ",factor,(double) ((MagickOffsetType)
2186               GetBlobSize(preview_image))/1024.0);
2187           else
2188             (void) FormatLocaleString(label,MaxTextExtent,"quality %s\n%.20gb ",
2189               factor,(double) ((MagickOffsetType) GetBlobSize(thumbnail)));
2190         break;
2191       }
2192     }
2193     thumbnail=DestroyImage(thumbnail);
2194     percentage+=12.5;
2195     radius+=0.5;
2196     sigma+=0.25;
2197     if (preview_image == (Image *) NULL)
2198       break;
2199     (void) DeleteImageProperty(preview_image,"label");
2200     (void) SetImageProperty(preview_image,"label",label,exception);
2201     AppendImageToList(&images,preview_image);
2202     proceed=SetImageProgress(image,PreviewImageTag,(MagickOffsetType) i,
2203       NumberTiles);
2204     if (proceed == MagickFalse)
2205       break;
2206   }
2207   if (images == (Image *) NULL)
2208     {
2209       preview_info=DestroyImageInfo(preview_info);
2210       return((Image *) NULL);
2211     }
2212   /*
2213     Create the montage.
2214   */
2215   montage_info=CloneMontageInfo(preview_info,(MontageInfo *) NULL);
2216   (void) CopyMagickString(montage_info->filename,image->filename,MaxTextExtent);
2217   montage_info->shadow=MagickTrue;
2218   (void) CloneString(&montage_info->tile,"3x3");
2219   (void) CloneString(&montage_info->geometry,DefaultPreviewGeometry);
2220   (void) CloneString(&montage_info->frame,DefaultTileFrame);
2221   montage_image=MontageImages(images,montage_info,exception);
2222   montage_info=DestroyMontageInfo(montage_info);
2223   images=DestroyImageList(images);
2224   if (montage_image == (Image *) NULL)
2225     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2226   if (montage_image->montage != (char *) NULL)
2227     {
2228       /*
2229         Free image directory.
2230       */
2231       montage_image->montage=(char *) RelinquishMagickMemory(
2232         montage_image->montage);
2233       if (image->directory != (char *) NULL)
2234         montage_image->directory=(char *) RelinquishMagickMemory(
2235           montage_image->directory);
2236     }
2237   preview_info=DestroyImageInfo(preview_info);
2238   return(montage_image);
2239 }
2240 \f
2241 /*
2242 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2243 %                                                                             %
2244 %                                                                             %
2245 %                                                                             %
2246 %     R a d i a l B l u r I m a g e                                           %
2247 %                                                                             %
2248 %                                                                             %
2249 %                                                                             %
2250 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2251 %
2252 %  RadialBlurImage() applies a radial blur to the image.
2253 %
2254 %  Andrew Protano contributed this effect.
2255 %
2256 %  The format of the RadialBlurImage method is:
2257 %
2258 %    Image *RadialBlurImage(const Image *image,const double angle,
2259 %      ExceptionInfo *exception)
2260 %
2261 %  A description of each parameter follows:
2262 %
2263 %    o image: the image.
2264 %
2265 %    o angle: the angle of the radial blur.
2266 %
2267 %    o blur: the blur.
2268 %
2269 %    o exception: return any errors or warnings in this structure.
2270 %
2271 */
2272 MagickExport Image *RadialBlurImage(const Image *image,const double angle,
2273   ExceptionInfo *exception)
2274 {
2275   CacheView
2276     *blur_view,
2277     *image_view,
2278     *radial_view;
2279
2280   Image
2281     *blur_image;
2282
2283   MagickBooleanType
2284     status;
2285
2286   MagickOffsetType
2287     progress;
2288
2289   double
2290     blur_radius,
2291     *cos_theta,
2292     offset,
2293     *sin_theta,
2294     theta;
2295
2296   PointInfo
2297     blur_center;
2298
2299   register ssize_t
2300     i;
2301
2302   size_t
2303     n;
2304
2305   ssize_t
2306     y;
2307
2308   /*
2309     Allocate blur image.
2310   */
2311   assert(image != (Image *) NULL);
2312   assert(image->signature == MagickSignature);
2313   if (image->debug != MagickFalse)
2314     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2315   assert(exception != (ExceptionInfo *) NULL);
2316   assert(exception->signature == MagickSignature);
2317   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
2318   if (blur_image == (Image *) NULL)
2319     return((Image *) NULL);
2320   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
2321     {
2322       blur_image=DestroyImage(blur_image);
2323       return((Image *) NULL);
2324     }
2325   blur_center.x=(double) (image->columns-1)/2.0;
2326   blur_center.y=(double) (image->rows-1)/2.0;
2327   blur_radius=hypot(blur_center.x,blur_center.y);
2328   n=(size_t) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+2UL);
2329   theta=DegreesToRadians(angle)/(double) (n-1);
2330   cos_theta=(double *) AcquireQuantumMemory((size_t) n,
2331     sizeof(*cos_theta));
2332   sin_theta=(double *) AcquireQuantumMemory((size_t) n,
2333     sizeof(*sin_theta));
2334   if ((cos_theta == (double *) NULL) ||
2335       (sin_theta == (double *) NULL))
2336     {
2337       blur_image=DestroyImage(blur_image);
2338       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2339     }
2340   offset=theta*(double) (n-1)/2.0;
2341   for (i=0; i < (ssize_t) n; i++)
2342   {
2343     cos_theta[i]=cos((double) (theta*i-offset));
2344     sin_theta[i]=sin((double) (theta*i-offset));
2345   }
2346   /*
2347     Radial blur image.
2348   */
2349   status=MagickTrue;
2350   progress=0;
2351   image_view=AcquireVirtualCacheView(image,exception);
2352   radial_view=AcquireVirtualCacheView(image,exception);
2353   blur_view=AcquireAuthenticCacheView(blur_image,exception);
2354 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2355   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2356     magick_threads(image,blur_image,image->rows,1)
2357 #endif
2358   for (y=0; y < (ssize_t) image->rows; y++)
2359   {
2360     register const Quantum
2361       *restrict p;
2362
2363     register Quantum
2364       *restrict q;
2365
2366     register ssize_t
2367       x;
2368
2369     if (status == MagickFalse)
2370       continue;
2371     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2372     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
2373       exception);
2374     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2375       {
2376         status=MagickFalse;
2377         continue;
2378       }
2379     for (x=0; x < (ssize_t) image->columns; x++)
2380     {
2381       double
2382         radius;
2383
2384       PointInfo
2385         center;
2386
2387       register ssize_t
2388         i;
2389
2390       size_t
2391         step;
2392
2393       center.x=(double) x-blur_center.x;
2394       center.y=(double) y-blur_center.y;
2395       radius=hypot((double) center.x,center.y);
2396       if (radius == 0)
2397         step=1;
2398       else
2399         {
2400           step=(size_t) (blur_radius/radius);
2401           if (step == 0)
2402             step=1;
2403           else
2404             if (step >= n)
2405               step=n-1;
2406         }
2407       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2408       {
2409         double
2410           gamma,
2411           pixel;
2412
2413         PixelChannel
2414           channel;
2415
2416         PixelTrait
2417           blur_traits,
2418           traits;
2419
2420         register const Quantum
2421           *restrict r;
2422
2423         register ssize_t
2424           j;
2425
2426         channel=GetPixelChannelChannel(image,i);
2427         traits=GetPixelChannelTraits(image,channel);
2428         blur_traits=GetPixelChannelTraits(blur_image,channel);
2429         if ((traits == UndefinedPixelTrait) ||
2430             (blur_traits == UndefinedPixelTrait))
2431           continue;
2432         if (((blur_traits & CopyPixelTrait) != 0) ||
2433             (GetPixelReadMask(image,p) == 0))
2434           {
2435             SetPixelChannel(blur_image,channel,p[i],q);
2436             continue;
2437           }
2438         gamma=0.0;
2439         pixel=0.0;
2440         if ((blur_traits & BlendPixelTrait) == 0)
2441           {
2442             for (j=0; j < (ssize_t) n; j+=(ssize_t) step)
2443             {
2444               r=GetCacheViewVirtualPixels(radial_view, (ssize_t) (blur_center.x+
2445                 center.x*cos_theta[j]-center.y*sin_theta[j]+0.5),(ssize_t)
2446                 (blur_center.y+center.x*sin_theta[j]+center.y*cos_theta[j]+0.5),
2447                 1,1,exception);
2448               if (r == (const Quantum *) NULL)
2449                 {
2450                   status=MagickFalse;
2451                   continue;
2452                 }
2453               pixel+=r[i];
2454               gamma++;
2455             }
2456             gamma=PerceptibleReciprocal(gamma);
2457             SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
2458             continue;
2459           }
2460         for (j=0; j < (ssize_t) n; j+=(ssize_t) step)
2461         {
2462           r=GetCacheViewVirtualPixels(radial_view, (ssize_t) (blur_center.x+
2463             center.x*cos_theta[j]-center.y*sin_theta[j]+0.5),(ssize_t)
2464             (blur_center.y+center.x*sin_theta[j]+center.y*cos_theta[j]+0.5),
2465             1,1,exception);
2466           if (r == (const Quantum *) NULL)
2467             {
2468               status=MagickFalse;
2469               continue;
2470             }
2471           pixel+=GetPixelAlpha(image,r)*r[i];
2472           gamma+=GetPixelAlpha(image,r);
2473         }
2474         gamma=PerceptibleReciprocal(gamma);
2475         SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
2476       }
2477       p+=GetPixelChannels(image);
2478       q+=GetPixelChannels(blur_image);
2479     }
2480     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
2481       status=MagickFalse;
2482     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2483       {
2484         MagickBooleanType
2485           proceed;
2486
2487 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2488         #pragma omp critical (MagickCore_RadialBlurImage)
2489 #endif
2490         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
2491         if (proceed == MagickFalse)
2492           status=MagickFalse;
2493       }
2494   }
2495   blur_view=DestroyCacheView(blur_view);
2496   radial_view=DestroyCacheView(radial_view);
2497   image_view=DestroyCacheView(image_view);
2498   cos_theta=(double *) RelinquishMagickMemory(cos_theta);
2499   sin_theta=(double *) RelinquishMagickMemory(sin_theta);
2500   if (status == MagickFalse)
2501     blur_image=DestroyImage(blur_image);
2502   return(blur_image);
2503 }
2504 \f
2505 /*
2506 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2507 %                                                                             %
2508 %                                                                             %
2509 %                                                                             %
2510 %     S e l e c t i v e B l u r I m a g e                                     %
2511 %                                                                             %
2512 %                                                                             %
2513 %                                                                             %
2514 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2515 %
2516 %  SelectiveBlurImage() selectively blur pixels within a contrast threshold.
2517 %  It is similar to the unsharpen mask that sharpens everything with contrast
2518 %  above a certain threshold.
2519 %
2520 %  The format of the SelectiveBlurImage method is:
2521 %
2522 %      Image *SelectiveBlurImage(const Image *image,const double radius,
2523 %        const double sigma,const double threshold,ExceptionInfo *exception)
2524 %
2525 %  A description of each parameter follows:
2526 %
2527 %    o image: the image.
2528 %
2529 %    o radius: the radius of the Gaussian, in pixels, not counting the center
2530 %      pixel.
2531 %
2532 %    o sigma: the standard deviation of the Gaussian, in pixels.
2533 %
2534 %    o threshold: only pixels within this contrast threshold are included
2535 %      in the blur operation.
2536 %
2537 %    o exception: return any errors or warnings in this structure.
2538 %
2539 */
2540 MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
2541   const double sigma,const double threshold,ExceptionInfo *exception)
2542 {
2543 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
2544
2545   CacheView
2546     *blur_view,
2547     *image_view,
2548     *luminance_view;
2549
2550   Image
2551     *blur_image,
2552     *luminance_image;
2553
2554   MagickBooleanType
2555     status;
2556
2557   MagickOffsetType
2558     progress;
2559
2560   MagickRealType
2561     *kernel;
2562
2563   register ssize_t
2564     i;
2565
2566   size_t
2567     width;
2568
2569   ssize_t
2570     center,
2571     j,
2572     u,
2573     v,
2574     y;
2575
2576   /*
2577     Initialize blur image attributes.
2578   */
2579   assert(image != (Image *) NULL);
2580   assert(image->signature == MagickSignature);
2581   if (image->debug != MagickFalse)
2582     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2583   assert(exception != (ExceptionInfo *) NULL);
2584   assert(exception->signature == MagickSignature);
2585   width=GetOptimalKernelWidth1D(radius,sigma);
2586   kernel=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory((size_t)
2587     width,width*sizeof(*kernel)));
2588   if (kernel == (MagickRealType *) NULL)
2589     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
2590   j=(ssize_t) (width-1)/2;
2591   i=0;
2592   for (v=(-j); v <= j; v++)
2593   {
2594     for (u=(-j); u <= j; u++)
2595       kernel[i++]=(MagickRealType) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
2596         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
2597   }
2598   if (image->debug != MagickFalse)
2599     {
2600       char
2601         format[MaxTextExtent],
2602         *message;
2603
2604       register const MagickRealType
2605         *k;
2606
2607       ssize_t
2608         u,
2609         v;
2610
2611       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
2612         "  SelectiveBlurImage with %.20gx%.20g kernel:",(double) width,(double)
2613         width);
2614       message=AcquireString("");
2615       k=kernel;
2616       for (v=0; v < (ssize_t) width; v++)
2617       {
2618         *message='\0';
2619         (void) FormatLocaleString(format,MaxTextExtent,"%.20g: ",(double) v);
2620         (void) ConcatenateString(&message,format);
2621         for (u=0; u < (ssize_t) width; u++)
2622         {
2623           (void) FormatLocaleString(format,MaxTextExtent,"%+f ",(double) *k++);
2624           (void) ConcatenateString(&message,format);
2625         }
2626         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
2627       }
2628       message=DestroyString(message);
2629     }
2630   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
2631   if (blur_image == (Image *) NULL)
2632     return((Image *) NULL);
2633   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
2634     {
2635       blur_image=DestroyImage(blur_image);
2636       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2637       return((Image *) NULL);
2638     }
2639   luminance_image=CloneImage(image,0,0,MagickTrue,exception);
2640   if (luminance_image == (Image *) NULL)
2641     {
2642       blur_image=DestroyImage(blur_image);
2643       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2644       return((Image *) NULL);
2645     }
2646   status=TransformImageColorspace(luminance_image,GRAYColorspace,exception);
2647   if (status == MagickFalse)
2648     {
2649       luminance_image=DestroyImage(luminance_image);
2650       blur_image=DestroyImage(blur_image);
2651       kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2652       return((Image *) NULL);
2653     }
2654   /*
2655     Threshold blur image.
2656   */
2657   status=MagickTrue;
2658   progress=0;
2659   center=(ssize_t) (GetPixelChannels(image)*(image->columns+width)*
2660     ((width-1)/2L)+GetPixelChannels(image)*((width-1)/2L));
2661   image_view=AcquireVirtualCacheView(image,exception);
2662   luminance_view=AcquireVirtualCacheView(luminance_image,exception);
2663   blur_view=AcquireAuthenticCacheView(blur_image,exception);
2664 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2665   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2666     magick_threads(image,blur_image,image->rows,1)
2667 #endif
2668   for (y=0; y < (ssize_t) image->rows; y++)
2669   {
2670     double
2671       contrast;
2672
2673     MagickBooleanType
2674       sync;
2675
2676     register const Quantum
2677       *restrict l,
2678       *restrict p;
2679
2680     register Quantum
2681       *restrict q;
2682
2683     register ssize_t
2684       x;
2685
2686     if (status == MagickFalse)
2687       continue;
2688     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) (width-1)/2L),y-(ssize_t)
2689       ((width-1)/2L),image->columns+width,width,exception);
2690     l=GetCacheViewVirtualPixels(luminance_view,-((ssize_t) (width-1)/2L),y-
2691       (ssize_t) ((width-1)/2L),luminance_image->columns+width,width,exception);
2692     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
2693       exception);
2694     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2695       {
2696         status=MagickFalse;
2697         continue;
2698       }
2699     for (x=0; x < (ssize_t) image->columns; x++)
2700     {
2701       double
2702         intensity;
2703
2704       register ssize_t
2705         i;
2706
2707       intensity=GetPixelIntensity(image,p+center);
2708       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
2709       {
2710         double
2711           alpha,
2712           gamma,
2713           pixel;
2714
2715         PixelChannel
2716           channel;
2717
2718         PixelTrait
2719           blur_traits,
2720           traits;
2721
2722         register const MagickRealType
2723           *restrict k;
2724
2725         register const Quantum
2726           *restrict luminance_pixels,
2727           *restrict pixels;
2728
2729         register ssize_t
2730           u;
2731
2732         ssize_t
2733           v;
2734
2735         channel=GetPixelChannelChannel(image,i);
2736         traits=GetPixelChannelTraits(image,channel);
2737         blur_traits=GetPixelChannelTraits(blur_image,channel);
2738         if ((traits == UndefinedPixelTrait) ||
2739             (blur_traits == UndefinedPixelTrait))
2740           continue;
2741         if (((blur_traits & CopyPixelTrait) != 0) ||
2742             (GetPixelReadMask(image,p+center) == 0))
2743           {
2744             SetPixelChannel(blur_image,channel,p[center+i],q);
2745             continue;
2746           }
2747         k=kernel;
2748         pixel=0.0;
2749         pixels=p;
2750         luminance_pixels=l;
2751         gamma=0.0;
2752         if ((blur_traits & BlendPixelTrait) == 0)
2753           {
2754             for (v=0; v < (ssize_t) width; v++)
2755             {
2756               for (u=0; u < (ssize_t) width; u++)
2757               {
2758                 contrast=GetPixelIntensity(luminance_image,luminance_pixels)-
2759                   intensity;
2760                 if (fabs(contrast) < threshold)
2761                   {
2762                     pixel+=(*k)*pixels[i];
2763                     gamma+=(*k);
2764                   }
2765                 k++;
2766                 pixels+=GetPixelChannels(image);
2767                 luminance_pixels+=GetPixelChannels(luminance_image);
2768               }
2769               pixels+=(image->columns-1)*GetPixelChannels(image);
2770               luminance_pixels+=luminance_image->columns*
2771                 GetPixelChannels(luminance_image);
2772             }
2773             if (fabs((double) gamma) < MagickEpsilon)
2774               {
2775                 SetPixelChannel(blur_image,channel,p[center+i],q);
2776                 continue;
2777               }
2778             gamma=PerceptibleReciprocal(gamma);
2779             SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
2780             continue;
2781           }
2782         for (v=0; v < (ssize_t) width; v++)
2783         {
2784           for (u=0; u < (ssize_t) width; u++)
2785           {
2786             contrast=GetPixelIntensity(image,pixels)-intensity;
2787             if (fabs(contrast) < threshold)
2788               {
2789                 alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
2790                 pixel+=(*k)*alpha*pixels[i];
2791                 gamma+=(*k)*alpha;
2792               }
2793             k++;
2794             pixels+=GetPixelChannels(image);
2795             luminance_pixels+=GetPixelChannels(luminance_image);
2796           }
2797           pixels+=(image->columns-1)*GetPixelChannels(image);
2798           luminance_pixels+=luminance_image->columns*
2799             GetPixelChannels(luminance_image);
2800         }
2801         if (fabs((double) gamma) < MagickEpsilon)
2802           {
2803             SetPixelChannel(blur_image,channel,p[center+i],q);
2804             continue;
2805           }
2806         gamma=PerceptibleReciprocal(gamma);
2807         SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
2808       }
2809       p+=GetPixelChannels(image);
2810       l+=GetPixelChannels(luminance_image);
2811       q+=GetPixelChannels(blur_image);
2812     }
2813     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
2814     if (sync == MagickFalse)
2815       status=MagickFalse;
2816     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2817       {
2818         MagickBooleanType
2819           proceed;
2820
2821 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2822         #pragma omp critical (MagickCore_SelectiveBlurImage)
2823 #endif
2824         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress++,
2825           image->rows);
2826         if (proceed == MagickFalse)
2827           status=MagickFalse;
2828       }
2829   }
2830   blur_image->type=image->type;
2831   blur_view=DestroyCacheView(blur_view);
2832   image_view=DestroyCacheView(image_view);
2833   luminance_image=DestroyImage(luminance_image);
2834   kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
2835   if (status == MagickFalse)
2836     blur_image=DestroyImage(blur_image);
2837   return(blur_image);
2838 }
2839 \f
2840 /*
2841 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2842 %                                                                             %
2843 %                                                                             %
2844 %                                                                             %
2845 %     S h a d e I m a g e                                                     %
2846 %                                                                             %
2847 %                                                                             %
2848 %                                                                             %
2849 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2850 %
2851 %  ShadeImage() shines a distant light on an image to create a
2852 %  three-dimensional effect. You control the positioning of the light with
2853 %  azimuth and elevation; azimuth is measured in degrees off the x axis
2854 %  and elevation is measured in pixels above the Z axis.
2855 %
2856 %  The format of the ShadeImage method is:
2857 %
2858 %      Image *ShadeImage(const Image *image,const MagickBooleanType gray,
2859 %        const double azimuth,const double elevation,ExceptionInfo *exception)
2860 %
2861 %  A description of each parameter follows:
2862 %
2863 %    o image: the image.
2864 %
2865 %    o gray: A value other than zero shades the intensity of each pixel.
2866 %
2867 %    o azimuth, elevation:  Define the light source direction.
2868 %
2869 %    o exception: return any errors or warnings in this structure.
2870 %
2871 */
2872 MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
2873   const double azimuth,const double elevation,ExceptionInfo *exception)
2874 {
2875 #define ShadeImageTag  "Shade/Image"
2876
2877   CacheView
2878     *image_view,
2879     *shade_view;
2880
2881   Image
2882     *linear_image,
2883     *shade_image;
2884
2885   MagickBooleanType
2886     status;
2887
2888   MagickOffsetType
2889     progress;
2890
2891   PrimaryInfo
2892     light;
2893
2894   ssize_t
2895     y;
2896
2897   /*
2898     Initialize shaded image attributes.
2899   */
2900   assert(image != (const Image *) NULL);
2901   assert(image->signature == MagickSignature);
2902   if (image->debug != MagickFalse)
2903     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2904   assert(exception != (ExceptionInfo *) NULL);
2905   assert(exception->signature == MagickSignature);
2906   linear_image=CloneImage(image,0,0,MagickTrue,exception);
2907   shade_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
2908   if ((linear_image == (Image *) NULL) || (shade_image == (Image *) NULL))
2909     {
2910       if (linear_image != (Image *) NULL)
2911         linear_image=DestroyImage(linear_image);
2912       if (shade_image != (Image *) NULL)
2913         shade_image=DestroyImage(shade_image);
2914       return((Image *) NULL);
2915     }
2916   if (SetImageStorageClass(shade_image,DirectClass,exception) == MagickFalse)
2917     {
2918       linear_image=DestroyImage(linear_image);
2919       shade_image=DestroyImage(shade_image);
2920       return((Image *) NULL);
2921     }
2922   /*
2923     Compute the light vector.
2924   */
2925   light.x=(double) QuantumRange*cos(DegreesToRadians(azimuth))*
2926     cos(DegreesToRadians(elevation));
2927   light.y=(double) QuantumRange*sin(DegreesToRadians(azimuth))*
2928     cos(DegreesToRadians(elevation));
2929   light.z=(double) QuantumRange*sin(DegreesToRadians(elevation));
2930   /*
2931     Shade image.
2932   */
2933   status=MagickTrue;
2934   progress=0;
2935   image_view=AcquireVirtualCacheView(linear_image,exception);
2936   shade_view=AcquireAuthenticCacheView(shade_image,exception);
2937 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2938   #pragma omp parallel for schedule(static,4) shared(progress,status) \
2939     magick_threads(linear_image,shade_image,linear_image->rows,1)
2940 #endif
2941   for (y=0; y < (ssize_t) linear_image->rows; y++)
2942   {
2943     double
2944       distance,
2945       normal_distance,
2946       shade;
2947
2948     PrimaryInfo
2949       normal;
2950
2951     register const Quantum
2952       *restrict center,
2953       *restrict p,
2954       *restrict post,
2955       *restrict pre;
2956
2957     register Quantum
2958       *restrict q;
2959
2960     register ssize_t
2961       x;
2962
2963     if (status == MagickFalse)
2964       continue;
2965     p=GetCacheViewVirtualPixels(image_view,-1,y-1,linear_image->columns+2,3,
2966       exception);
2967     q=QueueCacheViewAuthenticPixels(shade_view,0,y,shade_image->columns,1,
2968       exception);
2969     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2970       {
2971         status=MagickFalse;
2972         continue;
2973       }
2974     /*
2975       Shade this row of pixels.
2976     */
2977     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
2978     pre=p+GetPixelChannels(linear_image);
2979     center=pre+(linear_image->columns+2)*GetPixelChannels(linear_image);
2980     post=center+(linear_image->columns+2)*GetPixelChannels(linear_image);
2981     for (x=0; x < (ssize_t) linear_image->columns; x++)
2982     {
2983       register ssize_t
2984         i;
2985
2986       /*
2987         Determine the surface normal and compute shading.
2988       */
2989       normal.x=(double) (
2990         GetPixelIntensity(linear_image,pre-GetPixelChannels(linear_image))+
2991         GetPixelIntensity(linear_image,center-GetPixelChannels(linear_image))+
2992         GetPixelIntensity(linear_image,post-GetPixelChannels(linear_image))-
2993         GetPixelIntensity(linear_image,pre+GetPixelChannels(linear_image))-
2994         GetPixelIntensity(linear_image,center+GetPixelChannels(linear_image))-
2995         GetPixelIntensity(linear_image,post+GetPixelChannels(linear_image)));
2996       normal.y=(double) (
2997         GetPixelIntensity(linear_image,post-GetPixelChannels(linear_image))+
2998         GetPixelIntensity(linear_image,post)+
2999         GetPixelIntensity(linear_image,post+GetPixelChannels(linear_image))-
3000         GetPixelIntensity(linear_image,pre-GetPixelChannels(linear_image))-
3001         GetPixelIntensity(linear_image,pre)-
3002         GetPixelIntensity(linear_image,pre+GetPixelChannels(linear_image)));
3003       if ((normal.x == 0.0) && (normal.y == 0.0))
3004         shade=light.z;
3005       else
3006         {
3007           shade=0.0;
3008           distance=normal.x*light.x+normal.y*light.y+normal.z*light.z;
3009           if (distance > MagickEpsilon)
3010             {
3011               normal_distance=normal.x*normal.x+normal.y*normal.y+
3012                 normal.z*normal.z;
3013               if (normal_distance > (MagickEpsilon*MagickEpsilon))
3014                 shade=distance/sqrt((double) normal_distance);
3015             }
3016         }
3017       for (i=0; i < (ssize_t) GetPixelChannels(linear_image); i++)
3018       {
3019         PixelChannel
3020           channel;
3021
3022         PixelTrait
3023           shade_traits,
3024           traits;
3025
3026         channel=GetPixelChannelChannel(linear_image,i);
3027         traits=GetPixelChannelTraits(linear_image,channel);
3028         shade_traits=GetPixelChannelTraits(shade_image,channel);
3029         if ((traits == UndefinedPixelTrait) ||
3030             (shade_traits == UndefinedPixelTrait))
3031           continue;
3032         if (((shade_traits & CopyPixelTrait) != 0) ||
3033             (GetPixelReadMask(linear_image,center) == 0))
3034           {
3035             SetPixelChannel(shade_image,channel,center[i],q);
3036             continue;
3037           }
3038         if (gray != MagickFalse)
3039           {
3040             SetPixelChannel(shade_image,channel,ClampToQuantum(shade),q);
3041             continue;
3042           }
3043         SetPixelChannel(shade_image,channel,ClampToQuantum(QuantumScale*shade*
3044           center[i]),q);
3045       }
3046       pre+=GetPixelChannels(linear_image);
3047       center+=GetPixelChannels(linear_image);
3048       post+=GetPixelChannels(linear_image);
3049       q+=GetPixelChannels(shade_image);
3050     }
3051     if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
3052       status=MagickFalse;
3053     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3054       {
3055         MagickBooleanType
3056           proceed;
3057
3058 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3059         #pragma omp critical (MagickCore_ShadeImage)
3060 #endif
3061         proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
3062         if (proceed == MagickFalse)
3063           status=MagickFalse;
3064       }
3065   }
3066   shade_view=DestroyCacheView(shade_view);
3067   image_view=DestroyCacheView(image_view);
3068   linear_image=DestroyImage(linear_image);
3069   if (status == MagickFalse)
3070     shade_image=DestroyImage(shade_image);
3071   return(shade_image);
3072 }
3073 \f
3074 /*
3075 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3076 %                                                                             %
3077 %                                                                             %
3078 %                                                                             %
3079 %     S h a r p e n I m a g e                                                 %
3080 %                                                                             %
3081 %                                                                             %
3082 %                                                                             %
3083 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3084 %
3085 %  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
3086 %  operator of the given radius and standard deviation (sigma).  For
3087 %  reasonable results, radius should be larger than sigma.  Use a radius of 0
3088 %  and SharpenImage() selects a suitable radius for you.
3089 %
3090 %  Using a separable kernel would be faster, but the negative weights cancel
3091 %  out on the corners of the kernel producing often undesirable ringing in the
3092 %  filtered result; this can be avoided by using a 2D gaussian shaped image
3093 %  sharpening kernel instead.
3094 %
3095 %  The format of the SharpenImage method is:
3096 %
3097 %    Image *SharpenImage(const Image *image,const double radius,
3098 %      const double sigma,ExceptionInfo *exception)
3099 %
3100 %  A description of each parameter follows:
3101 %
3102 %    o image: the image.
3103 %
3104 %    o radius: the radius of the Gaussian, in pixels, not counting the center
3105 %      pixel.
3106 %
3107 %    o sigma: the standard deviation of the Laplacian, in pixels.
3108 %
3109 %    o exception: return any errors or warnings in this structure.
3110 %
3111 */
3112 MagickExport Image *SharpenImage(const Image *image,const double radius,
3113   const double sigma,ExceptionInfo *exception)
3114 {
3115   double
3116     gamma,
3117     normalize;
3118
3119   Image
3120     *sharp_image;
3121
3122   KernelInfo
3123     *kernel_info;
3124
3125   register ssize_t
3126     i;
3127
3128   size_t
3129     width;
3130
3131   ssize_t
3132     j,
3133     u,
3134     v;
3135
3136   assert(image != (const Image *) NULL);
3137   assert(image->signature == MagickSignature);
3138   if (image->debug != MagickFalse)
3139     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3140   assert(exception != (ExceptionInfo *) NULL);
3141   assert(exception->signature == MagickSignature);
3142   width=GetOptimalKernelWidth2D(radius,sigma);
3143   kernel_info=AcquireKernelInfo((const char *) NULL);
3144   if (kernel_info == (KernelInfo *) NULL)
3145     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3146   (void) ResetMagickMemory(kernel_info,0,sizeof(*kernel_info));
3147   kernel_info->width=width;
3148   kernel_info->height=width;
3149   kernel_info->x=(ssize_t) (width-1)/2;
3150   kernel_info->y=(ssize_t) (width-1)/2;
3151   kernel_info->signature=MagickSignature;
3152   kernel_info->values=(MagickRealType *) MagickAssumeAligned(
3153     AcquireAlignedMemory(kernel_info->width,kernel_info->height*
3154     sizeof(*kernel_info->values)));
3155   if (kernel_info->values == (MagickRealType *) NULL)
3156     {
3157       kernel_info=DestroyKernelInfo(kernel_info);
3158       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
3159     }
3160   normalize=0.0;
3161   j=(ssize_t) (kernel_info->width-1)/2;
3162   i=0;
3163   for (v=(-j); v <= j; v++)
3164   {
3165     for (u=(-j); u <= j; u++)
3166     {
3167       kernel_info->values[i]=(MagickRealType) (-exp(-((double) u*u+v*v)/(2.0*
3168         MagickSigma*MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
3169       normalize+=kernel_info->values[i];
3170       i++;
3171     }
3172   }
3173   kernel_info->values[i/2]=(double) ((-2.0)*normalize);
3174   normalize=0.0;
3175   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
3176     normalize+=kernel_info->values[i];
3177   gamma=PerceptibleReciprocal(normalize);
3178   for (i=0; i < (ssize_t) (kernel_info->width*kernel_info->height); i++)
3179     kernel_info->values[i]*=gamma;
3180   sharp_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
3181     UndefinedCompositeOp,0.0,exception);
3182   kernel_info=DestroyKernelInfo(kernel_info);
3183   return(sharp_image);
3184 }
3185 \f
3186 /*
3187 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3188 %                                                                             %
3189 %                                                                             %
3190 %                                                                             %
3191 %     S p r e a d I m a g e                                                   %
3192 %                                                                             %
3193 %                                                                             %
3194 %                                                                             %
3195 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3196 %
3197 %  SpreadImage() is a special effects method that randomly displaces each
3198 %  pixel in a block defined by the radius parameter.
3199 %
3200 %  The format of the SpreadImage method is:
3201 %
3202 %      Image *SpreadImage(const Image *image,const double radius,
3203 %        const PixelInterpolateMethod method,ExceptionInfo *exception)
3204 %
3205 %  A description of each parameter follows:
3206 %
3207 %    o image: the image.
3208 %
3209 %    o radius:  choose a random pixel in a neighborhood of this extent.
3210 %
3211 %    o method:  the pixel interpolation method.
3212 %
3213 %    o exception: return any errors or warnings in this structure.
3214 %
3215 */
3216 MagickExport Image *SpreadImage(const Image *image,const double radius,
3217   const PixelInterpolateMethod method,ExceptionInfo *exception)
3218 {
3219 #define SpreadImageTag  "Spread/Image"
3220
3221   CacheView
3222     *image_view,
3223     *spread_view;
3224
3225   Image
3226     *spread_image;
3227
3228   MagickBooleanType
3229     status;
3230
3231   MagickOffsetType
3232     progress;
3233
3234   RandomInfo
3235     **restrict random_info;
3236
3237   size_t
3238     width;
3239
3240   ssize_t
3241     y;
3242
3243 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3244   unsigned long
3245     key;
3246 #endif
3247
3248   /*
3249     Initialize spread image attributes.
3250   */
3251   assert(image != (Image *) NULL);
3252   assert(image->signature == MagickSignature);
3253   if (image->debug != MagickFalse)
3254     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3255   assert(exception != (ExceptionInfo *) NULL);
3256   assert(exception->signature == MagickSignature);
3257   spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
3258     exception);
3259   if (spread_image == (Image *) NULL)
3260     return((Image *) NULL);
3261   if (SetImageStorageClass(spread_image,DirectClass,exception) == MagickFalse)
3262     {
3263       spread_image=DestroyImage(spread_image);
3264       return((Image *) NULL);
3265     }
3266   /*
3267     Spread image.
3268   */
3269   status=MagickTrue;
3270   progress=0;
3271   width=GetOptimalKernelWidth1D(radius,0.5);
3272   random_info=AcquireRandomInfoThreadSet();
3273   image_view=AcquireVirtualCacheView(image,exception);
3274   spread_view=AcquireAuthenticCacheView(spread_image,exception);
3275 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3276   key=GetRandomSecretKey(random_info[0]);
3277   #pragma omp parallel for schedule(static,4) shared(progress,status) \
3278     magick_threads(image,spread_image,image->rows,key == ~0UL)
3279 #endif
3280   for (y=0; y < (ssize_t) image->rows; y++)
3281   {
3282     const int
3283       id = GetOpenMPThreadId();
3284
3285     register Quantum
3286       *restrict q;
3287
3288     register ssize_t
3289       x;
3290
3291     if (status == MagickFalse)
3292       continue;
3293     q=QueueCacheViewAuthenticPixels(spread_view,0,y,spread_image->columns,1,
3294       exception);
3295     if (q == (Quantum *) NULL)
3296       {
3297         status=MagickFalse;
3298         continue;
3299       }
3300     for (x=0; x < (ssize_t) image->columns; x++)
3301     {
3302       PointInfo
3303         point;
3304
3305       point.x=GetPseudoRandomValue(random_info[id]);
3306       point.y=GetPseudoRandomValue(random_info[id]);
3307       status=InterpolatePixelChannels(image,image_view,spread_image,method,
3308         (double) x+width*point.x-0.5,(double) y+width*point.y-0.5,q,exception);
3309       q+=GetPixelChannels(spread_image);
3310     }
3311     if (SyncCacheViewAuthenticPixels(spread_view,exception) == MagickFalse)
3312       status=MagickFalse;
3313     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3314       {
3315         MagickBooleanType
3316           proceed;
3317
3318 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3319         #pragma omp critical (MagickCore_SpreadImage)
3320 #endif
3321         proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
3322         if (proceed == MagickFalse)
3323           status=MagickFalse;
3324       }
3325   }
3326   spread_view=DestroyCacheView(spread_view);
3327   image_view=DestroyCacheView(image_view);
3328   random_info=DestroyRandomInfoThreadSet(random_info);
3329   if (status == MagickFalse)
3330     spread_image=DestroyImage(spread_image);
3331   return(spread_image);
3332 }
3333 \f
3334 /*
3335 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3336 %                                                                             %
3337 %                                                                             %
3338 %                                                                             %
3339 %     U n s h a r p M a s k I m a g e                                         %
3340 %                                                                             %
3341 %                                                                             %
3342 %                                                                             %
3343 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3344 %
3345 %  UnsharpMaskImage() sharpens one or more image channels.  We convolve the
3346 %  image with a Gaussian operator of the given radius and standard deviation
3347 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
3348 %  radius of 0 and UnsharpMaskImage() selects a suitable radius for you.
3349 %
3350 %  The format of the UnsharpMaskImage method is:
3351 %
3352 %    Image *UnsharpMaskImage(const Image *image,const double radius,
3353 %      const double sigma,const double amount,const double threshold,
3354 %      ExceptionInfo *exception)
3355 %
3356 %  A description of each parameter follows:
3357 %
3358 %    o image: the image.
3359 %
3360 %    o radius: the radius of the Gaussian, in pixels, not counting the center
3361 %      pixel.
3362 %
3363 %    o sigma: the standard deviation of the Gaussian, in pixels.
3364 %
3365 %    o gain: the percentage of the difference between the original and the
3366 %      blur image that is added back into the original.
3367 %
3368 %    o threshold: the threshold in pixels needed to apply the diffence gain.
3369 %
3370 %    o exception: return any errors or warnings in this structure.
3371 %
3372 */
3373 MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
3374   const double sigma,const double gain,const double threshold,
3375   ExceptionInfo *exception)
3376 {
3377 #define SharpenImageTag  "Sharpen/Image"
3378
3379   CacheView
3380     *image_view,
3381     *unsharp_view;
3382
3383   Image
3384     *unsharp_image;
3385
3386   MagickBooleanType
3387     status;
3388
3389   MagickOffsetType
3390     progress;
3391
3392   double
3393     quantum_threshold;
3394
3395   ssize_t
3396     y;
3397
3398   assert(image != (const Image *) NULL);
3399   assert(image->signature == MagickSignature);
3400   if (image->debug != MagickFalse)
3401     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3402   assert(exception != (ExceptionInfo *) NULL);
3403   unsharp_image=BlurImage(image,radius,sigma,exception);
3404   if (unsharp_image == (Image *) NULL)
3405     return((Image *) NULL);
3406   quantum_threshold=(double) QuantumRange*threshold;
3407   /*
3408     Unsharp-mask image.
3409   */
3410   status=MagickTrue;
3411   progress=0;
3412   image_view=AcquireVirtualCacheView(image,exception);
3413   unsharp_view=AcquireAuthenticCacheView(unsharp_image,exception);
3414 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3415   #pragma omp parallel for schedule(static,4) shared(progress,status) \
3416     magick_threads(image,unsharp_image,image->rows,1)
3417 #endif
3418   for (y=0; y < (ssize_t) image->rows; y++)
3419   {
3420     register const Quantum
3421       *restrict p;
3422
3423     register Quantum
3424       *restrict q;
3425
3426     register ssize_t
3427       x;
3428
3429     if (status == MagickFalse)
3430       continue;
3431     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
3432     q=QueueCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
3433       exception);
3434     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
3435       {
3436         status=MagickFalse;
3437         continue;
3438       }
3439     for (x=0; x < (ssize_t) image->columns; x++)
3440     {
3441       register ssize_t
3442         i;
3443
3444       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
3445       {
3446         double
3447           pixel;
3448
3449         PixelChannel
3450           channel;
3451
3452         PixelTrait
3453           traits,
3454           unsharp_traits;
3455
3456         channel=GetPixelChannelChannel(image,i);
3457         traits=GetPixelChannelTraits(image,channel);
3458         unsharp_traits=GetPixelChannelTraits(unsharp_image,channel);
3459         if ((traits == UndefinedPixelTrait) ||
3460             (unsharp_traits == UndefinedPixelTrait))
3461           continue;
3462         if (((unsharp_traits & CopyPixelTrait) != 0) ||
3463             (GetPixelReadMask(image,p) != 0))
3464           {
3465             SetPixelChannel(unsharp_image,channel,p[i],q);
3466             continue;
3467           }
3468         pixel=p[i]-(double) GetPixelChannel(unsharp_image,channel,q);
3469         if (fabs(2.0*pixel) < quantum_threshold)
3470           pixel=(double) p[i];
3471         else
3472           pixel=(double) p[i]+gain*pixel;
3473         SetPixelChannel(unsharp_image,channel,ClampToQuantum(pixel),q);
3474       }
3475       p+=GetPixelChannels(image);
3476       q+=GetPixelChannels(unsharp_image);
3477     }
3478     if (SyncCacheViewAuthenticPixels(unsharp_view,exception) == MagickFalse)
3479       status=MagickFalse;
3480     if (image->progress_monitor != (MagickProgressMonitor) NULL)
3481       {
3482         MagickBooleanType
3483           proceed;
3484
3485 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3486         #pragma omp critical (MagickCore_UnsharpMaskImage)
3487 #endif
3488         proceed=SetImageProgress(image,SharpenImageTag,progress++,image->rows);
3489         if (proceed == MagickFalse)
3490           status=MagickFalse;
3491       }
3492   }
3493   unsharp_image->type=image->type;
3494   unsharp_view=DestroyCacheView(unsharp_view);
3495   image_view=DestroyCacheView(image_view);
3496   if (status == MagickFalse)
3497     unsharp_image=DestroyImage(unsharp_image);
3498   return(unsharp_image);
3499 }