effect.c

Go to the documentation of this file.
00001 /*
00002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00003 %                                                                             %
00004 %                                                                             %
00005 %                                                                             %
00006 %                   EEEEE  FFFFF  FFFFF  EEEEE  CCCC  TTTTT                   %
00007 %                   E      F      F      E     C        T                     %
00008 %                   EEE    FFF    FFF    EEE   C        T                     %
00009 %                   E      F      F      E     C        T                     %
00010 %                   EEEEE  F      F      EEEEE  CCCC    T                     %
00011 %                                                                             %
00012 %                                                                             %
00013 %                       MagickCore Image Effects Methods                      %
00014 %                                                                             %
00015 %                               Software Design                               %
00016 %                                 John Cristy                                 %
00017 %                                 October 1996                                %
00018 %                                                                             %
00019 %                                                                             %
00020 %  Copyright 1999-2010 ImageMagick Studio LLC, a non-profit organization      %
00021 %  dedicated to making software imaging solutions freely available.           %
00022 %                                                                             %
00023 %  You may not use this file except in compliance with the License.  You may  %
00024 %  obtain a copy of the License at                                            %
00025 %                                                                             %
00026 %    http://www.imagemagick.org/script/license.php                            %
00027 %                                                                             %
00028 %  Unless required by applicable law or agreed to in writing, software        %
00029 %  distributed under the License is distributed on an "AS IS" BASIS,          %
00030 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
00031 %  See the License for the specific language governing permissions and        %
00032 %  limitations under the License.                                             %
00033 %                                                                             %
00034 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00035 %
00036 %
00037 %
00038 */
00039 
00040 /*
00041   Include declarations.
00042 */
00043 #include "magick/studio.h"
00044 #include "magick/accelerate.h"
00045 #include "magick/blob.h"
00046 #include "magick/cache-view.h"
00047 #include "magick/color.h"
00048 #include "magick/color-private.h"
00049 #include "magick/colorspace.h"
00050 #include "magick/constitute.h"
00051 #include "magick/decorate.h"
00052 #include "magick/draw.h"
00053 #include "magick/enhance.h"
00054 #include "magick/exception.h"
00055 #include "magick/exception-private.h"
00056 #include "magick/effect.h"
00057 #include "magick/fx.h"
00058 #include "magick/gem.h"
00059 #include "magick/geometry.h"
00060 #include "magick/image-private.h"
00061 #include "magick/list.h"
00062 #include "magick/log.h"
00063 #include "magick/memory_.h"
00064 #include "magick/monitor.h"
00065 #include "magick/monitor-private.h"
00066 #include "magick/montage.h"
00067 #include "magick/morphology.h"
00068 #include "magick/paint.h"
00069 #include "magick/pixel-private.h"
00070 #include "magick/property.h"
00071 #include "magick/quantize.h"
00072 #include "magick/quantum.h"
00073 #include "magick/random_.h"
00074 #include "magick/random-private.h"
00075 #include "magick/resample.h"
00076 #include "magick/resample-private.h"
00077 #include "magick/resize.h"
00078 #include "magick/resource_.h"
00079 #include "magick/segment.h"
00080 #include "magick/shear.h"
00081 #include "magick/signature-private.h"
00082 #include "magick/string_.h"
00083 #include "magick/thread-private.h"
00084 #include "magick/transform.h"
00085 #include "magick/threshold.h"
00086 
00087 /*
00088 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00089 %                                                                             %
00090 %                                                                             %
00091 %                                                                             %
00092 %     A d a p t i v e B l u r I m a g e                                       %
00093 %                                                                             %
00094 %                                                                             %
00095 %                                                                             %
00096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00097 %
00098 %  AdaptiveBlurImage() adaptively blurs the image by blurring less
00099 %  intensely near image edges and more intensely far from edges.  We blur the
00100 %  image with a Gaussian operator of the given radius and standard deviation
00101 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
00102 %  radius of 0 and AdaptiveBlurImage() selects a suitable radius for you.
00103 %
00104 %  The format of the AdaptiveBlurImage method is:
00105 %
00106 %      Image *AdaptiveBlurImage(const Image *image,const double radius,
00107 %        const double sigma,ExceptionInfo *exception)
00108 %      Image *AdaptiveBlurImageChannel(const Image *image,
00109 %        const ChannelType channel,double radius,const double sigma,
00110 %        ExceptionInfo *exception)
00111 %
00112 %  A description of each parameter follows:
00113 %
00114 %    o image: the image.
00115 %
00116 %    o channel: the channel type.
00117 %
00118 %    o radius: the radius of the Gaussian, in pixels, not counting the center
00119 %      pixel.
00120 %
00121 %    o sigma: the standard deviation of the Laplacian, in pixels.
00122 %
00123 %    o exception: return any errors or warnings in this structure.
00124 %
00125 */
00126 
00127 MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
00128   const double sigma,ExceptionInfo *exception)
00129 {
00130   Image
00131     *blur_image;
00132 
00133   blur_image=AdaptiveBlurImageChannel(image,DefaultChannels,radius,sigma,
00134     exception);
00135   return(blur_image);
00136 }
00137 
00138 MagickExport Image *AdaptiveBlurImageChannel(const Image *image,
00139   const ChannelType channel,const double radius,const double sigma,
00140   ExceptionInfo *exception)
00141 {
00142 #define AdaptiveBlurImageTag  "Convolve/Image"
00143 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
00144 
00145   CacheView
00146     *blur_view,
00147     *edge_view,
00148     *image_view;
00149 
00150   double
00151     **kernel,
00152     normalize;
00153 
00154   Image
00155     *blur_image,
00156     *edge_image,
00157     *gaussian_image;
00158 
00159   long
00160     j,
00161     k,
00162     progress,
00163     u,
00164     v,
00165     y;
00166 
00167   MagickBooleanType
00168     status;
00169 
00170   MagickPixelPacket
00171     bias;
00172 
00173   register long
00174     i;
00175 
00176   unsigned long
00177     width;
00178 
00179   assert(image != (const Image *) NULL);
00180   assert(image->signature == MagickSignature);
00181   if (image->debug != MagickFalse)
00182     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00183   assert(exception != (ExceptionInfo *) NULL);
00184   assert(exception->signature == MagickSignature);
00185   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
00186   if (blur_image == (Image *) NULL)
00187     return((Image *) NULL);
00188   if (fabs(sigma) <= MagickEpsilon)
00189     return(blur_image);
00190   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
00191     {
00192       InheritException(exception,&blur_image->exception);
00193       blur_image=DestroyImage(blur_image);
00194       return((Image *) NULL);
00195     }
00196   /*
00197     Edge detect the image brighness channel, level, blur, and level again.
00198   */
00199   edge_image=EdgeImage(image,radius,exception);
00200   if (edge_image == (Image *) NULL)
00201     {
00202       blur_image=DestroyImage(blur_image);
00203       return((Image *) NULL);
00204     }
00205   (void) LevelImage(edge_image,"20%,95%");
00206   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
00207   if (gaussian_image != (Image *) NULL)
00208     {
00209       edge_image=DestroyImage(edge_image);
00210       edge_image=gaussian_image;
00211     }
00212   (void) LevelImage(edge_image,"10%,95%");
00213   /*
00214     Create a set of kernels from maximum (radius,sigma) to minimum.
00215   */
00216   width=GetOptimalKernelWidth2D(radius,sigma);
00217   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
00218   if (kernel == (double **) NULL)
00219     {
00220       edge_image=DestroyImage(edge_image);
00221       blur_image=DestroyImage(blur_image);
00222       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
00223     }
00224   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
00225   for (i=0; i < (long) width; i+=2)
00226   {
00227     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
00228       sizeof(**kernel));
00229     if (kernel[i] == (double *) NULL)
00230       break;
00231     normalize=0.0;
00232     j=(long) (width-i)/2;
00233     k=0;
00234     for (v=(-j); v <= j; v++)
00235     {
00236       for (u=(-j); u <= j; u++)
00237       {
00238         kernel[i][k]=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
00239           (2.0*MagickPI*MagickSigma*MagickSigma);
00240         normalize+=kernel[i][k];
00241         k++;
00242       }
00243     }
00244     if (fabs(normalize) <= MagickEpsilon)
00245       normalize=1.0;
00246     normalize=1.0/normalize;
00247     for (k=0; k < (j*j); k++)
00248       kernel[i][k]=normalize*kernel[i][k];
00249   }
00250   if (i < (long) width)
00251     {
00252       for (i-=2; i >= 0; i-=2)
00253         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
00254       kernel=(double **) RelinquishMagickMemory(kernel);
00255       edge_image=DestroyImage(edge_image);
00256       blur_image=DestroyImage(blur_image);
00257       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
00258     }
00259   /*
00260     Adaptively blur image.
00261   */
00262   status=MagickTrue;
00263   progress=0;
00264   GetMagickPixelPacket(image,&bias);
00265   SetMagickPixelPacketBias(image,&bias);
00266   image_view=AcquireCacheView(image);
00267   edge_view=AcquireCacheView(edge_image);
00268   blur_view=AcquireCacheView(blur_image);
00269 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00270   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
00271 #endif
00272   for (y=0; y < (long) blur_image->rows; y++)
00273   {
00274     register const IndexPacket
00275       *restrict indexes;
00276 
00277     register const PixelPacket
00278       *restrict p,
00279       *restrict r;
00280 
00281     register IndexPacket
00282       *restrict blur_indexes;
00283 
00284     register long
00285       x;
00286 
00287     register PixelPacket
00288       *restrict q;
00289 
00290     if (status == MagickFalse)
00291       continue;
00292     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
00293     q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
00294       exception);
00295     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
00296       {
00297         status=MagickFalse;
00298         continue;
00299       }
00300     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
00301     for (x=0; x < (long) blur_image->columns; x++)
00302     {
00303       MagickPixelPacket
00304         pixel;
00305 
00306       MagickRealType
00307         alpha,
00308         gamma;
00309 
00310       register const double
00311         *restrict k;
00312 
00313       register long
00314         i,
00315         u,
00316         v;
00317 
00318       gamma=0.0;
00319       i=(long) ceil(width*QuantumScale*PixelIntensity(r)-0.5);
00320       if (i < 0)
00321         i=0;
00322       else
00323         if (i > (long) width)
00324           i=(long) width;
00325       if ((i & 0x01) != 0)
00326         i--;
00327       p=GetCacheViewVirtualPixels(image_view,x-((long) (width-i)/2L),y-(long)
00328         ((width-i)/2L),width-i,width-i,exception);
00329       if (p == (const PixelPacket *) NULL)
00330         break;
00331       indexes=GetCacheViewVirtualIndexQueue(image_view);
00332       pixel=bias;
00333       k=kernel[i];
00334       for (v=0; v < (long) (width-i); v++)
00335       {
00336         for (u=0; u < (long) (width-i); u++)
00337         {
00338           alpha=1.0;
00339           if (((channel & OpacityChannel) != 0) &&
00340               (image->matte != MagickFalse))
00341             alpha=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(p));
00342           if ((channel & RedChannel) != 0)
00343             pixel.red+=(*k)*alpha*GetRedPixelComponent(p);
00344           if ((channel & GreenChannel) != 0)
00345             pixel.green+=(*k)*alpha*GetGreenPixelComponent(p);
00346           if ((channel & BlueChannel) != 0)
00347             pixel.blue+=(*k)*alpha*GetBluePixelComponent(p);
00348           if ((channel & OpacityChannel) != 0)
00349             pixel.opacity+=(*k)*GetOpacityPixelComponent(p);
00350           if (((channel & IndexChannel) != 0) &&
00351               (image->colorspace == CMYKColorspace))
00352             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
00353           gamma+=(*k)*alpha;
00354           k++;
00355           p++;
00356         }
00357       }
00358       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00359       if ((channel & RedChannel) != 0)
00360         q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
00361       if ((channel & GreenChannel) != 0)
00362         q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
00363       if ((channel & BlueChannel) != 0)
00364         q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
00365       if ((channel & OpacityChannel) != 0)
00366         SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
00367       if (((channel & IndexChannel) != 0) &&
00368           (image->colorspace == CMYKColorspace))
00369         blur_indexes[x]=ClampToQuantum(gamma*GetIndexPixelComponent(&pixel));
00370       q++;
00371       r++;
00372     }
00373     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
00374       status=MagickFalse;
00375     if (image->progress_monitor != (MagickProgressMonitor) NULL)
00376       {
00377         MagickBooleanType
00378           proceed;
00379 
00380 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00381   #pragma omp critical (MagickCore_AdaptiveBlurImageChannel)
00382 #endif
00383         proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress++,
00384           image->rows);
00385         if (proceed == MagickFalse)
00386           status=MagickFalse;
00387       }
00388   }
00389   blur_image->type=image->type;
00390   blur_view=DestroyCacheView(blur_view);
00391   edge_view=DestroyCacheView(edge_view);
00392   image_view=DestroyCacheView(image_view);
00393   edge_image=DestroyImage(edge_image);
00394   for (i=0; i < (long) width;  i+=2)
00395     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
00396   kernel=(double **) RelinquishMagickMemory(kernel);
00397   if (status == MagickFalse)
00398     blur_image=DestroyImage(blur_image);
00399   return(blur_image);
00400 }
00401 
00402 /*
00403 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00404 %                                                                             %
00405 %                                                                             %
00406 %                                                                             %
00407 %     A d a p t i v e S h a r p e n I m a g e                                 %
00408 %                                                                             %
00409 %                                                                             %
00410 %                                                                             %
00411 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00412 %
00413 %  AdaptiveSharpenImage() adaptively sharpens the image by sharpening more
00414 %  intensely near image edges and less intensely far from edges. We sharpen the
00415 %  image with a Gaussian operator of the given radius and standard deviation
00416 %  (sigma).  For reasonable results, radius should be larger than sigma.  Use a
00417 %  radius of 0 and AdaptiveSharpenImage() selects a suitable radius for you.
00418 %
00419 %  The format of the AdaptiveSharpenImage method is:
00420 %
00421 %      Image *AdaptiveSharpenImage(const Image *image,const double radius,
00422 %        const double sigma,ExceptionInfo *exception)
00423 %      Image *AdaptiveSharpenImageChannel(const Image *image,
00424 %        const ChannelType channel,double radius,const double sigma,
00425 %        ExceptionInfo *exception)
00426 %
00427 %  A description of each parameter follows:
00428 %
00429 %    o image: the image.
00430 %
00431 %    o channel: the channel type.
00432 %
00433 %    o radius: the radius of the Gaussian, in pixels, not counting the center
00434 %      pixel.
00435 %
00436 %    o sigma: the standard deviation of the Laplacian, in pixels.
00437 %
00438 %    o exception: return any errors or warnings in this structure.
00439 %
00440 */
00441 
00442 MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
00443   const double sigma,ExceptionInfo *exception)
00444 {
00445   Image
00446     *sharp_image;
00447 
00448   sharp_image=AdaptiveSharpenImageChannel(image,DefaultChannels,radius,sigma,
00449     exception);
00450   return(sharp_image);
00451 }
00452 
00453 MagickExport Image *AdaptiveSharpenImageChannel(const Image *image,
00454   const ChannelType channel,const double radius,const double sigma,
00455   ExceptionInfo *exception)
00456 {
00457 #define AdaptiveSharpenImageTag  "Convolve/Image"
00458 #define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
00459 
00460   CacheView
00461     *sharp_view,
00462     *edge_view,
00463     *image_view;
00464 
00465   double
00466     **kernel,
00467     normalize;
00468 
00469   Image
00470     *sharp_image,
00471     *edge_image,
00472     *gaussian_image;
00473 
00474   long
00475     j,
00476     k,
00477     progress,
00478     u,
00479     v,
00480     y;
00481 
00482   MagickBooleanType
00483     status;
00484 
00485   MagickPixelPacket
00486     bias;
00487 
00488   register long
00489     i;
00490 
00491   unsigned long
00492     width;
00493 
00494   assert(image != (const Image *) NULL);
00495   assert(image->signature == MagickSignature);
00496   if (image->debug != MagickFalse)
00497     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00498   assert(exception != (ExceptionInfo *) NULL);
00499   assert(exception->signature == MagickSignature);
00500   sharp_image=CloneImage(image,0,0,MagickTrue,exception);
00501   if (sharp_image == (Image *) NULL)
00502     return((Image *) NULL);
00503   if (fabs(sigma) <= MagickEpsilon)
00504     return(sharp_image);
00505   if (SetImageStorageClass(sharp_image,DirectClass) == MagickFalse)
00506     {
00507       InheritException(exception,&sharp_image->exception);
00508       sharp_image=DestroyImage(sharp_image);
00509       return((Image *) NULL);
00510     }
00511   /*
00512     Edge detect the image brighness channel, level, sharp, and level again.
00513   */
00514   edge_image=EdgeImage(image,radius,exception);
00515   if (edge_image == (Image *) NULL)
00516     {
00517       sharp_image=DestroyImage(sharp_image);
00518       return((Image *) NULL);
00519     }
00520   (void) LevelImage(edge_image,"20%,95%");
00521   gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
00522   if (gaussian_image != (Image *) NULL)
00523     {
00524       edge_image=DestroyImage(edge_image);
00525       edge_image=gaussian_image;
00526     }
00527   (void) LevelImage(edge_image,"10%,95%");
00528   /*
00529     Create a set of kernels from maximum (radius,sigma) to minimum.
00530   */
00531   width=GetOptimalKernelWidth2D(radius,sigma);
00532   kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
00533   if (kernel == (double **) NULL)
00534     {
00535       edge_image=DestroyImage(edge_image);
00536       sharp_image=DestroyImage(sharp_image);
00537       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
00538     }
00539   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
00540   for (i=0; i < (long) width; i+=2)
00541   {
00542     kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
00543       sizeof(**kernel));
00544     if (kernel[i] == (double *) NULL)
00545       break;
00546     normalize=0.0;
00547     j=(long) (width-i)/2;
00548     k=0;
00549     for (v=(-j); v <= j; v++)
00550     {
00551       for (u=(-j); u <= j; u++)
00552       {
00553         kernel[i][k]=(-exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
00554           (2.0*MagickPI*MagickSigma*MagickSigma));
00555         normalize+=kernel[i][k];
00556         k++;
00557       }
00558     }
00559     if (fabs(normalize) <= MagickEpsilon)
00560       normalize=1.0;
00561     normalize=1.0/normalize;
00562     for (k=0; k < (j*j); k++)
00563       kernel[i][k]=normalize*kernel[i][k];
00564   }
00565   if (i < (long) width)
00566     {
00567       for (i-=2; i >= 0; i-=2)
00568         kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
00569       kernel=(double **) RelinquishMagickMemory(kernel);
00570       edge_image=DestroyImage(edge_image);
00571       sharp_image=DestroyImage(sharp_image);
00572       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
00573     }
00574   /*
00575     Adaptively sharpen image.
00576   */
00577   status=MagickTrue;
00578   progress=0;
00579   GetMagickPixelPacket(image,&bias);
00580   SetMagickPixelPacketBias(image,&bias);
00581   image_view=AcquireCacheView(image);
00582   edge_view=AcquireCacheView(edge_image);
00583   sharp_view=AcquireCacheView(sharp_image);
00584 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00585   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
00586 #endif
00587   for (y=0; y < (long) sharp_image->rows; y++)
00588   {
00589     register const IndexPacket
00590       *restrict indexes;
00591 
00592     register const PixelPacket
00593       *restrict p,
00594       *restrict r;
00595 
00596     register IndexPacket
00597       *restrict sharp_indexes;
00598 
00599     register long
00600       x;
00601 
00602     register PixelPacket
00603       *restrict q;
00604 
00605     if (status == MagickFalse)
00606       continue;
00607     r=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns,1,exception);
00608     q=QueueCacheViewAuthenticPixels(sharp_view,0,y,sharp_image->columns,1,
00609       exception);
00610     if ((r == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
00611       {
00612         status=MagickFalse;
00613         continue;
00614       }
00615     sharp_indexes=GetCacheViewAuthenticIndexQueue(sharp_view);
00616     for (x=0; x < (long) sharp_image->columns; x++)
00617     {
00618       MagickPixelPacket
00619         pixel;
00620 
00621       MagickRealType
00622         alpha,
00623         gamma;
00624 
00625       register const double
00626         *restrict k;
00627 
00628       register long
00629         i,
00630         u,
00631         v;
00632 
00633       gamma=0.0;
00634       i=(long) ceil(width*(QuantumRange-QuantumScale*PixelIntensity(r))-0.5);
00635       if (i < 0)
00636         i=0;
00637       else
00638         if (i > (long) width)
00639           i=(long) width;
00640       if ((i & 0x01) != 0)
00641         i--;
00642       p=GetCacheViewVirtualPixels(image_view,x-((long) (width-i)/2L),y-(long)
00643         ((width-i)/2L),width-i,width-i,exception);
00644       if (p == (const PixelPacket *) NULL)
00645         break;
00646       indexes=GetCacheViewVirtualIndexQueue(image_view);
00647       k=kernel[i];
00648       pixel=bias;
00649       for (v=0; v < (long) (width-i); v++)
00650       {
00651         for (u=0; u < (long) (width-i); u++)
00652         {
00653           alpha=1.0;
00654           if (((channel & OpacityChannel) != 0) &&
00655               (image->matte != MagickFalse))
00656             alpha=(MagickRealType) (QuantumScale*GetAlphaPixelComponent(p));
00657           if ((channel & RedChannel) != 0)
00658             pixel.red+=(*k)*alpha*GetRedPixelComponent(p);
00659           if ((channel & GreenChannel) != 0)
00660             pixel.green+=(*k)*alpha*GetGreenPixelComponent(p);
00661           if ((channel & BlueChannel) != 0)
00662             pixel.blue+=(*k)*alpha*GetBluePixelComponent(p);
00663           if ((channel & OpacityChannel) != 0)
00664             pixel.opacity+=(*k)*GetOpacityPixelComponent(p);
00665           if (((channel & IndexChannel) != 0) &&
00666               (image->colorspace == CMYKColorspace))
00667             pixel.index+=(*k)*alpha*indexes[x+(width-i)*v+u];
00668           gamma+=(*k)*alpha;
00669           k++;
00670           p++;
00671         }
00672       }
00673       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00674       if ((channel & RedChannel) != 0)
00675         q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
00676       if ((channel & GreenChannel) != 0)
00677         q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
00678       if ((channel & BlueChannel) != 0)
00679         q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
00680       if ((channel & OpacityChannel) != 0)
00681         SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
00682       if (((channel & IndexChannel) != 0) &&
00683           (image->colorspace == CMYKColorspace))
00684         sharp_indexes[x]=ClampToQuantum(gamma*GetIndexPixelComponent(&pixel));
00685       q++;
00686       r++;
00687     }
00688     if (SyncCacheViewAuthenticPixels(sharp_view,exception) == MagickFalse)
00689       status=MagickFalse;
00690     if (image->progress_monitor != (MagickProgressMonitor) NULL)
00691       {
00692         MagickBooleanType
00693           proceed;
00694 
00695 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00696   #pragma omp critical (MagickCore_AdaptiveSharpenImageChannel)
00697 #endif
00698         proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress++,
00699           image->rows);
00700         if (proceed == MagickFalse)
00701           status=MagickFalse;
00702       }
00703   }
00704   sharp_image->type=image->type;
00705   sharp_view=DestroyCacheView(sharp_view);
00706   edge_view=DestroyCacheView(edge_view);
00707   image_view=DestroyCacheView(image_view);
00708   edge_image=DestroyImage(edge_image);
00709   for (i=0; i < (long) width;  i+=2)
00710     kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
00711   kernel=(double **) RelinquishMagickMemory(kernel);
00712   if (status == MagickFalse)
00713     sharp_image=DestroyImage(sharp_image);
00714   return(sharp_image);
00715 }
00716 
00717 /*
00718 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00719 %                                                                             %
00720 %                                                                             %
00721 %                                                                             %
00722 %     B l u r I m a g e                                                       %
00723 %                                                                             %
00724 %                                                                             %
00725 %                                                                             %
00726 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00727 %
00728 %  BlurImage() blurs an image.  We convolve the image with a Gaussian operator
00729 %  of the given radius and standard deviation (sigma).  For reasonable results,
00730 %  the radius should be larger than sigma.  Use a radius of 0 and BlurImage()
00731 %  selects a suitable radius for you.
00732 %
00733 %  BlurImage() differs from GaussianBlurImage() in that it uses a separable
00734 %  kernel which is faster but mathematically equivalent to the non-separable
00735 %  kernel.
00736 %
00737 %  The format of the BlurImage method is:
00738 %
00739 %      Image *BlurImage(const Image *image,const double radius,
00740 %        const double sigma,ExceptionInfo *exception)
00741 %      Image *BlurImageChannel(const Image *image,const ChannelType channel,
00742 %        const double radius,const double sigma,ExceptionInfo *exception)
00743 %
00744 %  A description of each parameter follows:
00745 %
00746 %    o image: the image.
00747 %
00748 %    o channel: the channel type.
00749 %
00750 %    o radius: the radius of the Gaussian, in pixels, not counting the center
00751 %      pixel.
00752 %
00753 %    o sigma: the standard deviation of the Gaussian, in pixels.
00754 %
00755 %    o exception: return any errors or warnings in this structure.
00756 %
00757 */
00758 
00759 MagickExport Image *BlurImage(const Image *image,const double radius,
00760   const double sigma,ExceptionInfo *exception)
00761 {
00762   Image
00763     *blur_image;
00764 
00765   blur_image=BlurImageChannel(image,DefaultChannels,radius,sigma,exception);
00766   return(blur_image);
00767 }
00768 
00769 static double *GetBlurKernel(const unsigned long width,const double sigma)
00770 {
00771   double
00772     *kernel,
00773     normalize;
00774 
00775   long
00776     j,
00777     k;
00778 
00779   register long
00780     i;
00781 
00782   /*
00783     Generate a 1-D convolution kernel.
00784   */
00785   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
00786   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
00787   if (kernel == (double *) NULL)
00788     return(0);
00789   normalize=0.0;
00790   j=(long) width/2;
00791   i=0;
00792   for (k=(-j); k <= j; k++)
00793   {
00794     kernel[i]=exp(-((double) k*k)/(2.0*MagickSigma*MagickSigma))/
00795       (MagickSQ2PI*MagickSigma);
00796     normalize+=kernel[i];
00797     i++;
00798   }
00799   for (i=0; i < (long) width; i++)
00800     kernel[i]/=normalize;
00801   return(kernel);
00802 }
00803 
00804 MagickExport Image *BlurImageChannel(const Image *image,
00805   const ChannelType channel,const double radius,const double sigma,
00806   ExceptionInfo *exception)
00807 {
00808 #define BlurImageTag  "Blur/Image"
00809 
00810   CacheView
00811     *blur_view,
00812     *image_view;
00813 
00814   double
00815     *kernel;
00816 
00817   Image
00818     *blur_image;
00819 
00820   long
00821     progress,
00822     x,
00823     y;
00824 
00825   MagickBooleanType
00826     status;
00827 
00828   MagickPixelPacket
00829     bias;
00830 
00831   register long
00832     i;
00833 
00834   unsigned long
00835     width;
00836 
00837   /*
00838     Initialize blur image attributes.
00839   */
00840   assert(image != (Image *) NULL);
00841   assert(image->signature == MagickSignature);
00842   if (image->debug != MagickFalse)
00843     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00844   assert(exception != (ExceptionInfo *) NULL);
00845   assert(exception->signature == MagickSignature);
00846   blur_image=CloneImage(image,0,0,MagickTrue,exception);
00847   if (blur_image == (Image *) NULL)
00848     return((Image *) NULL);
00849   if (fabs(sigma) <= MagickEpsilon)
00850     return(blur_image);
00851   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
00852     {
00853       InheritException(exception,&blur_image->exception);
00854       blur_image=DestroyImage(blur_image);
00855       return((Image *) NULL);
00856     }
00857   width=GetOptimalKernelWidth1D(radius,sigma);
00858   kernel=GetBlurKernel(width,sigma);
00859   if (kernel == (double *) NULL)
00860     {
00861       blur_image=DestroyImage(blur_image);
00862       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
00863     }
00864   if (image->debug != MagickFalse)
00865     {
00866       char
00867         format[MaxTextExtent],
00868         *message;
00869 
00870       register const double
00871         *k;
00872 
00873       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
00874         "  BlurImage with %ld kernel:",width);
00875       message=AcquireString("");
00876       k=kernel;
00877       for (i=0; i < (long) width; i++)
00878       {
00879         *message='\0';
00880         (void) FormatMagickString(format,MaxTextExtent,"%ld: ",i);
00881         (void) ConcatenateString(&message,format);
00882         (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
00883         (void) ConcatenateString(&message,format);
00884         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
00885       }
00886       message=DestroyString(message);
00887     }
00888   /*
00889     Blur rows.
00890   */
00891   status=MagickTrue;
00892   progress=0;
00893   GetMagickPixelPacket(image,&bias);
00894   SetMagickPixelPacketBias(image,&bias);
00895   image_view=AcquireCacheView(image);
00896   blur_view=AcquireCacheView(blur_image);
00897 #if defined(MAGICKCORE_OPENMP_SUPPORT)
00898   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
00899 #endif
00900   for (y=0; y < (long) blur_image->rows; y++)
00901   {
00902     register const IndexPacket
00903       *restrict indexes;
00904 
00905     register const PixelPacket
00906       *restrict p;
00907 
00908     register IndexPacket
00909       *restrict blur_indexes;
00910 
00911     register long
00912       x;
00913 
00914     register PixelPacket
00915       *restrict q;
00916 
00917     if (status == MagickFalse)
00918       continue;
00919     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y,image->columns+
00920       width,1,exception);
00921     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
00922       exception);
00923     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
00924       {
00925         status=MagickFalse;
00926         continue;
00927       }
00928     indexes=GetCacheViewVirtualIndexQueue(image_view);
00929     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
00930     for (x=0; x < (long) blur_image->columns; x++)
00931     {
00932       MagickPixelPacket
00933         pixel;
00934 
00935       register const double
00936         *restrict k;
00937 
00938       register const PixelPacket
00939         *restrict kernel_pixels;
00940 
00941       register long
00942         i;
00943 
00944       pixel=bias;
00945       k=kernel;
00946       kernel_pixels=p;
00947       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
00948         {
00949           for (i=0; i < (long) width; i++)
00950           {
00951             pixel.red+=(*k)*kernel_pixels->red;
00952             pixel.green+=(*k)*kernel_pixels->green;
00953             pixel.blue+=(*k)*kernel_pixels->blue;
00954             k++;
00955             kernel_pixels++;
00956           }
00957           if ((channel & RedChannel) != 0)
00958             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
00959           if ((channel & GreenChannel) != 0)
00960             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
00961           if ((channel & BlueChannel) != 0)
00962             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
00963           if ((channel & OpacityChannel) != 0)
00964             {
00965               k=kernel;
00966               kernel_pixels=p;
00967               for (i=0; i < (long) width; i++)
00968               {
00969                 pixel.opacity+=(*k)*kernel_pixels->opacity;
00970                 k++;
00971                 kernel_pixels++;
00972               }
00973               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
00974             }
00975           if (((channel & IndexChannel) != 0) &&
00976               (image->colorspace == CMYKColorspace))
00977             {
00978               register const IndexPacket
00979                 *restrict kernel_indexes;
00980 
00981               k=kernel;
00982               kernel_indexes=indexes;
00983               for (i=0; i < (long) width; i++)
00984               {
00985                 pixel.index+=(*k)*(*kernel_indexes);
00986                 k++;
00987                 kernel_indexes++;
00988               }
00989               blur_indexes[x]=ClampToQuantum(pixel.index);
00990             }
00991         }
00992       else
00993         {
00994           MagickRealType
00995             alpha,
00996             gamma;
00997 
00998           gamma=0.0;
00999           for (i=0; i < (long) width; i++)
01000           {
01001             alpha=(MagickRealType) (QuantumScale*
01002               GetAlphaPixelComponent(kernel_pixels));
01003             pixel.red+=(*k)*alpha*kernel_pixels->red;
01004             pixel.green+=(*k)*alpha*kernel_pixels->green;
01005             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
01006             gamma+=(*k)*alpha;
01007             k++;
01008             kernel_pixels++;
01009           }
01010           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
01011           if ((channel & RedChannel) != 0)
01012             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
01013           if ((channel & GreenChannel) != 0)
01014             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
01015           if ((channel & BlueChannel) != 0)
01016             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
01017           if ((channel & OpacityChannel) != 0)
01018             {
01019               k=kernel;
01020               kernel_pixels=p;
01021               for (i=0; i < (long) width; i++)
01022               {
01023                 pixel.opacity+=(*k)*kernel_pixels->opacity;
01024                 k++;
01025                 kernel_pixels++;
01026               }
01027               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
01028             }
01029           if (((channel & IndexChannel) != 0) &&
01030               (image->colorspace == CMYKColorspace))
01031             {
01032               register const IndexPacket
01033                 *restrict kernel_indexes;
01034 
01035               k=kernel;
01036               kernel_pixels=p;
01037               kernel_indexes=indexes;
01038               for (i=0; i < (long) width; i++)
01039               {
01040                 alpha=(MagickRealType) (QuantumScale*
01041                   GetAlphaPixelComponent(kernel_pixels));
01042                 pixel.index+=(*k)*alpha*(*kernel_indexes);
01043                 k++;
01044                 kernel_pixels++;
01045                 kernel_indexes++;
01046               }
01047               blur_indexes[x]=ClampToQuantum(gamma*
01048                 GetIndexPixelComponent(&pixel));
01049             }
01050         }
01051       p++;
01052       q++;
01053     }
01054     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
01055       status=MagickFalse;
01056     if (image->progress_monitor != (MagickProgressMonitor) NULL)
01057       {
01058         MagickBooleanType
01059           proceed;
01060 
01061 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01062   #pragma omp critical (MagickCore_BlurImageChannel)
01063 #endif
01064         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
01065           blur_image->columns);
01066         if (proceed == MagickFalse)
01067           status=MagickFalse;
01068       }
01069   }
01070   blur_view=DestroyCacheView(blur_view);
01071   image_view=DestroyCacheView(image_view);
01072   /*
01073     Blur columns.
01074   */
01075   image_view=AcquireCacheView(blur_image);
01076   blur_view=AcquireCacheView(blur_image);
01077 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01078   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
01079 #endif
01080   for (x=0; x < (long) blur_image->columns; x++)
01081   {
01082     register const IndexPacket
01083       *restrict indexes;
01084 
01085     register const PixelPacket
01086       *restrict p;
01087 
01088     register IndexPacket
01089       *restrict blur_indexes;
01090 
01091     register long
01092       y;
01093 
01094     register PixelPacket
01095       *restrict q;
01096 
01097     if (status == MagickFalse)
01098       continue;
01099     p=GetCacheViewVirtualPixels(image_view,x,-((long) width/2L),1,image->rows+
01100       width,exception);
01101     q=GetCacheViewAuthenticPixels(blur_view,x,0,1,blur_image->rows,exception);
01102     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
01103       {
01104         status=MagickFalse;
01105         continue;
01106       }
01107     indexes=GetCacheViewVirtualIndexQueue(image_view);
01108     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
01109     for (y=0; y < (long) blur_image->rows; y++)
01110     {
01111       MagickPixelPacket
01112         pixel;
01113 
01114       register const double
01115         *restrict k;
01116 
01117       register const PixelPacket
01118         *restrict kernel_pixels;
01119 
01120       register long
01121         i;
01122 
01123       pixel=bias;
01124       k=kernel;
01125       kernel_pixels=p;
01126       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
01127         {
01128           for (i=0; i < (long) width; i++)
01129           {
01130             pixel.red+=(*k)*kernel_pixels->red;
01131             pixel.green+=(*k)*kernel_pixels->green;
01132             pixel.blue+=(*k)*kernel_pixels->blue;
01133             k++;
01134             kernel_pixels++;
01135           }
01136           if ((channel & RedChannel) != 0)
01137             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
01138           if ((channel & GreenChannel) != 0)
01139             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
01140           if ((channel & BlueChannel) != 0)
01141             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
01142           if ((channel & OpacityChannel) != 0)
01143             {
01144               k=kernel;
01145               kernel_pixels=p;
01146               for (i=0; i < (long) width; i++)
01147               {
01148                 pixel.opacity+=(*k)*kernel_pixels->opacity;
01149                 k++;
01150                 kernel_pixels++;
01151               }
01152               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
01153             }
01154           if (((channel & IndexChannel) != 0) &&
01155               (image->colorspace == CMYKColorspace))
01156             {
01157               register const IndexPacket
01158                 *restrict kernel_indexes;
01159 
01160               k=kernel;
01161               kernel_indexes=indexes;
01162               for (i=0; i < (long) width; i++)
01163               {
01164                 pixel.index+=(*k)*(*kernel_indexes);
01165                 k++;
01166                 kernel_indexes++;
01167               }
01168               blur_indexes[y]=ClampToQuantum(pixel.index);
01169             }
01170         }
01171       else
01172         {
01173           MagickRealType
01174             alpha,
01175             gamma;
01176 
01177           gamma=0.0;
01178           for (i=0; i < (long) width; i++)
01179           {
01180             alpha=(MagickRealType) (QuantumScale*
01181               GetAlphaPixelComponent(kernel_pixels));
01182             pixel.red+=(*k)*alpha*kernel_pixels->red;
01183             pixel.green+=(*k)*alpha*kernel_pixels->green;
01184             pixel.blue+=(*k)*alpha*kernel_pixels->blue;
01185             gamma+=(*k)*alpha;
01186             k++;
01187             kernel_pixels++;
01188           }
01189           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
01190           if ((channel & RedChannel) != 0)
01191             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
01192           if ((channel & GreenChannel) != 0)
01193             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
01194           if ((channel & BlueChannel) != 0)
01195             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
01196           if ((channel & OpacityChannel) != 0)
01197             {
01198               k=kernel;
01199               kernel_pixels=p;
01200               for (i=0; i < (long) width; i++)
01201               {
01202                 pixel.opacity+=(*k)*kernel_pixels->opacity;
01203                 k++;
01204                 kernel_pixels++;
01205               }
01206               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
01207             }
01208           if (((channel & IndexChannel) != 0) &&
01209               (image->colorspace == CMYKColorspace))
01210             {
01211               register const IndexPacket
01212                 *restrict kernel_indexes;
01213 
01214               k=kernel;
01215               kernel_pixels=p;
01216               kernel_indexes=indexes;
01217               for (i=0; i < (long) width; i++)
01218               {
01219                 alpha=(MagickRealType) (QuantumScale*
01220                   GetAlphaPixelComponent(kernel_pixels));
01221                 pixel.index+=(*k)*alpha*(*kernel_indexes);
01222                 k++;
01223                 kernel_pixels++;
01224                 kernel_indexes++;
01225               }
01226               blur_indexes[y]=ClampToQuantum(gamma*
01227                 GetIndexPixelComponent(&pixel));
01228             }
01229         }
01230       p++;
01231       q++;
01232     }
01233     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
01234       status=MagickFalse;
01235     if (image->progress_monitor != (MagickProgressMonitor) NULL)
01236       {
01237         MagickBooleanType
01238           proceed;
01239 
01240 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01241   #pragma omp critical (MagickCore_BlurImageChannel)
01242 #endif
01243         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
01244           blur_image->columns);
01245         if (proceed == MagickFalse)
01246           status=MagickFalse;
01247       }
01248   }
01249   blur_view=DestroyCacheView(blur_view);
01250   image_view=DestroyCacheView(image_view);
01251   kernel=(double *) RelinquishMagickMemory(kernel);
01252   if (status == MagickFalse)
01253     blur_image=DestroyImage(blur_image);
01254   blur_image->type=image->type;
01255   return(blur_image);
01256 }
01257 
01258 /*
01259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01260 %                                                                             %
01261 %                                                                             %
01262 %                                                                             %
01263 %     C o n v o l v e I m a g e                                               %
01264 %                                                                             %
01265 %                                                                             %
01266 %                                                                             %
01267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01268 %
01269 %  ConvolveImage() applies a custom convolution kernel to the image.
01270 %
01271 %  The format of the ConvolveImage method is:
01272 %
01273 %      Image *ConvolveImage(const Image *image,const unsigned long order,
01274 %        const double *kernel,ExceptionInfo *exception)
01275 %      Image *ConvolveImageChannel(const Image *image,const ChannelType channel,
01276 %        const unsigned long order,const double *kernel,
01277 %        ExceptionInfo *exception)
01278 %
01279 %  A description of each parameter follows:
01280 %
01281 %    o image: the image.
01282 %
01283 %    o channel: the channel type.
01284 %
01285 %    o order: the number of columns and rows in the filter kernel.
01286 %
01287 %    o kernel: An array of double representing the convolution kernel.
01288 %
01289 %    o exception: return any errors or warnings in this structure.
01290 %
01291 */
01292 
01293 MagickExport Image *ConvolveImage(const Image *image,const unsigned long order,
01294   const double *kernel,ExceptionInfo *exception)
01295 {
01296   Image
01297     *convolve_image;
01298 
01299   convolve_image=ConvolveImageChannel(image,DefaultChannels,order,kernel,
01300     exception);
01301   return(convolve_image);
01302 }
01303 
01304 MagickExport Image *ConvolveImageChannel(const Image *image,
01305   const ChannelType channel,const unsigned long order,const double *kernel,
01306   ExceptionInfo *exception)
01307 {
01308 #define ConvolveImageTag  "Convolve/Image"
01309 
01310   CacheView
01311     *convolve_view,
01312     *image_view;
01313 
01314   double
01315     *normal_kernel;
01316 
01317   Image
01318     *convolve_image;
01319 
01320   long
01321     progress,
01322     y;
01323 
01324   MagickBooleanType
01325     status;
01326 
01327   MagickPixelPacket
01328     bias;
01329 
01330   MagickRealType
01331     gamma;
01332 
01333   register long
01334     i;
01335 
01336   unsigned long
01337     width;
01338 
01339   /*
01340     Initialize convolve image attributes.
01341   */
01342   assert(image != (Image *) NULL);
01343   assert(image->signature == MagickSignature);
01344   if (image->debug != MagickFalse)
01345     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01346   assert(exception != (ExceptionInfo *) NULL);
01347   assert(exception->signature == MagickSignature);
01348   width=order;
01349   if ((width % 2) == 0)
01350     ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
01351   convolve_image=CloneImage(image,0,0,MagickTrue,exception);
01352   if (convolve_image == (Image *) NULL)
01353     return((Image *) NULL);
01354   if (SetImageStorageClass(convolve_image,DirectClass) == MagickFalse)
01355     {
01356       InheritException(exception,&convolve_image->exception);
01357       convolve_image=DestroyImage(convolve_image);
01358       return((Image *) NULL);
01359     }
01360   if (image->debug != MagickFalse)
01361     {
01362       char
01363         format[MaxTextExtent],
01364         *message;
01365 
01366       long
01367         u,
01368         v;
01369 
01370       register const double
01371         *k;
01372 
01373       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
01374         "  ConvolveImage with %ldx%ld kernel:",width,width);
01375       message=AcquireString("");
01376       k=kernel;
01377       for (v=0; v < (long) width; v++)
01378       {
01379         *message='\0';
01380         (void) FormatMagickString(format,MaxTextExtent,"%ld: ",v);
01381         (void) ConcatenateString(&message,format);
01382         for (u=0; u < (long) width; u++)
01383         {
01384           (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
01385           (void) ConcatenateString(&message,format);
01386         }
01387         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
01388       }
01389       message=DestroyString(message);
01390     }
01391   /*
01392     Normalize kernel.
01393   */
01394   normal_kernel=(double *) AcquireQuantumMemory(width*width,
01395     sizeof(*normal_kernel));
01396   if (normal_kernel == (double *) NULL)
01397     {
01398       convolve_image=DestroyImage(convolve_image);
01399       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
01400     }
01401   gamma=0.0;
01402   for (i=0; i < (long) (width*width); i++)
01403     gamma+=kernel[i];
01404   gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
01405   for (i=0; i < (long) (width*width); i++)
01406     normal_kernel[i]=gamma*kernel[i];
01407   /*
01408     Convolve image.
01409   */
01410   status=MagickTrue;
01411   progress=0;
01412   GetMagickPixelPacket(image,&bias);
01413   SetMagickPixelPacketBias(image,&bias);
01414   image_view=AcquireCacheView(image);
01415   convolve_view=AcquireCacheView(convolve_image);
01416 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01417   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
01418 #endif
01419   for (y=0; y < (long) image->rows; y++)
01420   {
01421     MagickBooleanType
01422       sync;
01423 
01424     register const IndexPacket
01425       *restrict indexes;
01426 
01427     register const PixelPacket
01428       *restrict p;
01429 
01430     register IndexPacket
01431       *restrict convolve_indexes;
01432 
01433     register long
01434       x;
01435 
01436     register PixelPacket
01437       *restrict q;
01438 
01439     if (status == MagickFalse)
01440       continue;
01441     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
01442       2L),image->columns+width,width,exception);
01443     q=GetCacheViewAuthenticPixels(convolve_view,0,y,convolve_image->columns,1,
01444       exception);
01445     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
01446       {
01447         status=MagickFalse;
01448         continue;
01449       }
01450     indexes=GetCacheViewVirtualIndexQueue(image_view);
01451     convolve_indexes=GetCacheViewAuthenticIndexQueue(convolve_view);
01452     for (x=0; x < (long) image->columns; x++)
01453     {
01454       long
01455         v;
01456 
01457       MagickPixelPacket
01458         pixel;
01459 
01460       register const double
01461         *restrict k;
01462 
01463       register const PixelPacket
01464         *restrict kernel_pixels;
01465 
01466       register long
01467         u;
01468 
01469       pixel=bias;
01470       k=normal_kernel;
01471       kernel_pixels=p;
01472       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
01473         {
01474           for (v=0; v < (long) width; v++)
01475           {
01476             for (u=0; u < (long) width; u++)
01477             {
01478               pixel.red+=(*k)*kernel_pixels[u].red;
01479               pixel.green+=(*k)*kernel_pixels[u].green;
01480               pixel.blue+=(*k)*kernel_pixels[u].blue;
01481               k++;
01482             }
01483             kernel_pixels+=image->columns+width;
01484           }
01485           if ((channel & RedChannel) != 0)
01486             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
01487           if ((channel & GreenChannel) != 0)
01488             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
01489           if ((channel & BlueChannel) != 0)
01490             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
01491           if ((channel & OpacityChannel) != 0)
01492             {
01493               k=normal_kernel;
01494               kernel_pixels=p;
01495               for (v=0; v < (long) width; v++)
01496               {
01497                 for (u=0; u < (long) width; u++)
01498                 {
01499                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
01500                   k++;
01501                 }
01502                 kernel_pixels+=image->columns+width;
01503               }
01504               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
01505             }
01506           if (((channel & IndexChannel) != 0) &&
01507               (image->colorspace == CMYKColorspace))
01508             {
01509               register const IndexPacket
01510                 *restrict kernel_indexes;
01511 
01512               k=normal_kernel;
01513               kernel_indexes=indexes;
01514               for (v=0; v < (long) width; v++)
01515               {
01516                 for (u=0; u < (long) width; u++)
01517                 {
01518                   pixel.index+=(*k)*kernel_indexes[u];
01519                   k++;
01520                 }
01521                 kernel_indexes+=image->columns+width;
01522               }
01523               convolve_indexes[x]=ClampToQuantum(pixel.index);
01524             }
01525         }
01526       else
01527         {
01528           MagickRealType
01529             alpha,
01530             gamma;
01531 
01532           gamma=0.0;
01533           for (v=0; v < (long) width; v++)
01534           {
01535             for (u=0; u < (long) width; u++)
01536             {
01537               alpha=(MagickRealType) (QuantumScale*(QuantumRange-
01538                 kernel_pixels[u].opacity));
01539               pixel.red+=(*k)*alpha*kernel_pixels[u].red;
01540               pixel.green+=(*k)*alpha*kernel_pixels[u].green;
01541               pixel.blue+=(*k)*alpha*kernel_pixels[u].blue;
01542               gamma+=(*k)*alpha;
01543               k++;
01544             }
01545             kernel_pixels+=image->columns+width;
01546           }
01547           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
01548           if ((channel & RedChannel) != 0)
01549             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
01550           if ((channel & GreenChannel) != 0)
01551             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
01552           if ((channel & BlueChannel) != 0)
01553             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
01554           if ((channel & OpacityChannel) != 0)
01555             {
01556               k=normal_kernel;
01557               kernel_pixels=p;
01558               for (v=0; v < (long) width; v++)
01559               {
01560                 for (u=0; u < (long) width; u++)
01561                 {
01562                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
01563                   k++;
01564                 }
01565                 kernel_pixels+=image->columns+width;
01566               }
01567               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
01568             }
01569           if (((channel & IndexChannel) != 0) &&
01570               (image->colorspace == CMYKColorspace))
01571             {
01572               register const IndexPacket
01573                 *restrict kernel_indexes;
01574 
01575               k=normal_kernel;
01576               kernel_pixels=p;
01577               kernel_indexes=indexes;
01578               for (v=0; v < (long) width; v++)
01579               {
01580                 for (u=0; u < (long) width; u++)
01581                 {
01582                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
01583                     kernel_pixels[u].opacity));
01584                   pixel.index+=(*k)*alpha*kernel_indexes[u];
01585                   k++;
01586                 }
01587                 kernel_pixels+=image->columns+width;
01588                 kernel_indexes+=image->columns+width;
01589               }
01590               convolve_indexes[x]=ClampToQuantum(gamma*
01591                 GetIndexPixelComponent(&pixel));
01592             }
01593         }
01594       p++;
01595       q++;
01596     }
01597     sync=SyncCacheViewAuthenticPixels(convolve_view,exception);
01598     if (sync == MagickFalse)
01599       status=MagickFalse;
01600     if (image->progress_monitor != (MagickProgressMonitor) NULL)
01601       {
01602         MagickBooleanType
01603           proceed;
01604 
01605 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01606   #pragma omp critical (MagickCore_ConvolveImageChannel)
01607 #endif
01608         proceed=SetImageProgress(image,ConvolveImageTag,progress++,image->rows);
01609         if (proceed == MagickFalse)
01610           status=MagickFalse;
01611       }
01612   }
01613   convolve_image->type=image->type;
01614   convolve_view=DestroyCacheView(convolve_view);
01615   image_view=DestroyCacheView(image_view);
01616   normal_kernel=(double *) RelinquishMagickMemory(normal_kernel);
01617   if (status == MagickFalse)
01618     convolve_image=DestroyImage(convolve_image);
01619   return(convolve_image);
01620 }
01621 
01622 /*
01623 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01624 %                                                                             %
01625 %                                                                             %
01626 %                                                                             %
01627 %     D e s p e c k l e I m a g e                                             %
01628 %                                                                             %
01629 %                                                                             %
01630 %                                                                             %
01631 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01632 %
01633 %  DespeckleImage() reduces the speckle noise in an image while perserving the
01634 %  edges of the original image.
01635 %
01636 %  The format of the DespeckleImage method is:
01637 %
01638 %      Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
01639 %
01640 %  A description of each parameter follows:
01641 %
01642 %    o image: the image.
01643 %
01644 %    o exception: return any errors or warnings in this structure.
01645 %
01646 */
01647 
01648 static Quantum **DestroyPixelThreadSet(Quantum **pixels)
01649 {
01650   register long
01651     i;
01652 
01653   assert(pixels != (Quantum **) NULL);
01654   for (i=0; i < (long) GetOpenMPMaximumThreads(); i++)
01655     if (pixels[i] != (Quantum *) NULL)
01656       pixels[i]=(Quantum *) RelinquishMagickMemory(pixels[i]);
01657   pixels=(Quantum **) RelinquishAlignedMemory(pixels);
01658   return(pixels);
01659 }
01660 
01661 static Quantum **AcquirePixelThreadSet(const size_t count)
01662 {
01663   register long
01664     i;
01665 
01666   Quantum
01667     **pixels;
01668 
01669   unsigned long
01670     number_threads;
01671 
01672   number_threads=GetOpenMPMaximumThreads();
01673   pixels=(Quantum **) AcquireAlignedMemory(number_threads,sizeof(*pixels));
01674   if (pixels == (Quantum **) NULL)
01675     return((Quantum **) NULL);
01676   (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
01677   for (i=0; i < (long) number_threads; i++)
01678   {
01679     pixels[i]=(Quantum *) AcquireQuantumMemory(count,sizeof(**pixels));
01680     if (pixels[i] == (Quantum *) NULL)
01681       return(DestroyPixelThreadSet(pixels));
01682   }
01683   return(pixels);
01684 }
01685 
01686 static void Hull(const long x_offset,const long y_offset,
01687   const unsigned long columns,const unsigned long rows,Quantum *f,Quantum *g,
01688   const int polarity)
01689 {
01690   long
01691     y;
01692 
01693   MagickRealType
01694     v;
01695 
01696   register long
01697     x;
01698 
01699   register Quantum
01700     *p,
01701     *q,
01702     *r,
01703     *s;
01704 
01705   assert(f != (Quantum *) NULL);
01706   assert(g != (Quantum *) NULL);
01707   p=f+(columns+2);
01708   q=g+(columns+2);
01709   r=p+(y_offset*((long) columns+2)+x_offset);
01710   for (y=0; y < (long) rows; y++)
01711   {
01712     p++;
01713     q++;
01714     r++;
01715     if (polarity > 0)
01716       for (x=(long) columns; x != 0; x--)
01717       {
01718         v=(MagickRealType) (*p);
01719         if ((MagickRealType) *r >= (v+(MagickRealType) ScaleCharToQuantum(2)))
01720           v+=ScaleCharToQuantum(1);
01721         *q=(Quantum) v;
01722         p++;
01723         q++;
01724         r++;
01725       }
01726     else
01727       for (x=(long) columns; x != 0; x--)
01728       {
01729         v=(MagickRealType) (*p);
01730         if ((MagickRealType) *r <= (v-(MagickRealType) ScaleCharToQuantum(2)))
01731           v-=(long) ScaleCharToQuantum(1);
01732         *q=(Quantum) v;
01733         p++;
01734         q++;
01735         r++;
01736       }
01737     p++;
01738     q++;
01739     r++;
01740   }
01741   p=f+(columns+2);
01742   q=g+(columns+2);
01743   r=q+(y_offset*((long) columns+2)+x_offset);
01744   s=q-(y_offset*((long) columns+2)+x_offset);
01745   for (y=0; y < (long) rows; y++)
01746   {
01747     p++;
01748     q++;
01749     r++;
01750     s++;
01751     if (polarity > 0)
01752       for (x=(long) columns; x != 0; x--)
01753       {
01754         v=(MagickRealType) (*q);
01755         if (((MagickRealType) *s >=
01756              (v+(MagickRealType) ScaleCharToQuantum(2))) &&
01757             ((MagickRealType) *r > v))
01758           v+=ScaleCharToQuantum(1);
01759         *p=(Quantum) v;
01760         p++;
01761         q++;
01762         r++;
01763         s++;
01764       }
01765     else
01766       for (x=(long) columns; x != 0; x--)
01767       {
01768         v=(MagickRealType) (*q);
01769         if (((MagickRealType) *s <=
01770              (v-(MagickRealType) ScaleCharToQuantum(2))) &&
01771             ((MagickRealType) *r < v))
01772           v-=(MagickRealType) ScaleCharToQuantum(1);
01773         *p=(Quantum) v;
01774         p++;
01775         q++;
01776         r++;
01777         s++;
01778       }
01779     p++;
01780     q++;
01781     r++;
01782     s++;
01783   }
01784 }
01785 
01786 MagickExport Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
01787 {
01788 #define DespeckleImageTag  "Despeckle/Image"
01789 
01790   CacheView
01791     *despeckle_view,
01792     *image_view;
01793 
01794   Image
01795     *despeckle_image;
01796 
01797   long
01798     channel;
01799 
01800   MagickBooleanType
01801     status;
01802 
01803   Quantum
01804     **restrict buffers,
01805     **restrict pixels;
01806 
01807   size_t
01808     length;
01809 
01810   static const int
01811     X[4] = {0, 1, 1,-1},
01812     Y[4] = {1, 0, 1, 1};
01813 
01814   /*
01815     Allocate despeckled image.
01816   */
01817   assert(image != (const Image *) NULL);
01818   assert(image->signature == MagickSignature);
01819   if (image->debug != MagickFalse)
01820     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01821   assert(exception != (ExceptionInfo *) NULL);
01822   assert(exception->signature == MagickSignature);
01823   despeckle_image=CloneImage(image,image->columns,image->rows,MagickTrue,
01824     exception);
01825   if (despeckle_image == (Image *) NULL)
01826     return((Image *) NULL);
01827   if (SetImageStorageClass(despeckle_image,DirectClass) == MagickFalse)
01828     {
01829       InheritException(exception,&despeckle_image->exception);
01830       despeckle_image=DestroyImage(despeckle_image);
01831       return((Image *) NULL);
01832     }
01833   /*
01834     Allocate image buffers.
01835   */
01836   length=(size_t) ((image->columns+2)*(image->rows+2));
01837   pixels=AcquirePixelThreadSet(length);
01838   buffers=AcquirePixelThreadSet(length);
01839   if ((pixels == (Quantum **) NULL) || (buffers == (Quantum **) NULL))
01840     {
01841       if (buffers != (Quantum **) NULL)
01842         buffers=DestroyPixelThreadSet(buffers);
01843       if (pixels != (Quantum **) NULL)
01844         pixels=DestroyPixelThreadSet(pixels);
01845       despeckle_image=DestroyImage(despeckle_image);
01846       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
01847     }
01848   /*
01849     Reduce speckle in the image.
01850   */
01851   status=MagickTrue;
01852   image_view=AcquireCacheView(image);
01853   despeckle_view=AcquireCacheView(despeckle_image);
01854 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01855   #pragma omp parallel for schedule(dynamic,4) shared(status)
01856 #endif
01857   for (channel=0; channel <= 3; channel++)
01858   {
01859     long
01860       j,
01861       y;
01862 
01863     register long
01864       i,
01865       id,
01866       x;
01867 
01868     register Quantum
01869       *buffer,
01870       *pixel;
01871 
01872     if (status == MagickFalse)
01873       continue;
01874     id=GetOpenMPThreadId();
01875     pixel=pixels[id];
01876     (void) ResetMagickMemory(pixel,0,length*sizeof(*pixel));
01877     buffer=buffers[id];
01878     j=(long) image->columns+2;
01879     for (y=0; y < (long) image->rows; y++)
01880     {
01881       register const PixelPacket
01882         *restrict p;
01883 
01884       p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
01885       if (p == (const PixelPacket *) NULL)
01886         break;
01887       j++;
01888       for (x=0; x < (long) image->columns; x++)
01889       {
01890         switch (channel)
01891         {
01892           case 0: pixel[j]=GetRedPixelComponent(p); break;
01893           case 1: pixel[j]=GetGreenPixelComponent(p); break;
01894           case 2: pixel[j]=GetBluePixelComponent(p); break;
01895           case 3: pixel[j]=GetOpacityPixelComponent(p); break;
01896           default: break;
01897         }
01898         p++;
01899         j++;
01900       }
01901       j++;
01902     }
01903     (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
01904     for (i=0; i < 4; i++)
01905     {
01906       Hull(X[i],Y[i],image->columns,image->rows,pixel,buffer,1);
01907       Hull(-X[i],-Y[i],image->columns,image->rows,pixel,buffer,1);
01908       Hull(-X[i],-Y[i],image->columns,image->rows,pixel,buffer,-1);
01909       Hull(X[i],Y[i],image->columns,image->rows,pixel,buffer,-1);
01910     }
01911     j=(long) image->columns+2;
01912     for (y=0; y < (long) image->rows; y++)
01913     {
01914       MagickBooleanType
01915         sync;
01916 
01917       register PixelPacket
01918         *restrict q;
01919 
01920       q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
01921         1,exception);
01922       if (q == (PixelPacket *) NULL)
01923         break;
01924       j++;
01925       for (x=0; x < (long) image->columns; x++)
01926       {
01927         switch (channel)
01928         {
01929           case 0: q->red=pixel[j]; break;
01930           case 1: q->green=pixel[j]; break;
01931           case 2: q->blue=pixel[j]; break;
01932           case 3: q->opacity=pixel[j]; break;
01933           default: break;
01934         }
01935         q++;
01936         j++;
01937       }
01938       sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
01939       if (sync == MagickFalse)
01940         {
01941           status=MagickFalse;
01942           break;
01943         }
01944       j++;
01945     }
01946     if (image->progress_monitor != (MagickProgressMonitor) NULL)
01947       {
01948         MagickBooleanType
01949           proceed;
01950 
01951 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01952   #pragma omp critical (MagickCore_DespeckleImage)
01953 #endif
01954         proceed=SetImageProgress(image,DespeckleImageTag,channel,3);
01955         if (proceed == MagickFalse)
01956           status=MagickFalse;
01957       }
01958   }
01959   despeckle_view=DestroyCacheView(despeckle_view);
01960   image_view=DestroyCacheView(image_view);
01961   buffers=DestroyPixelThreadSet(buffers);
01962   pixels=DestroyPixelThreadSet(pixels);
01963   despeckle_image->type=image->type;
01964   if (status == MagickFalse)
01965     despeckle_image=DestroyImage(despeckle_image);
01966   return(despeckle_image);
01967 }
01968 
01969 /*
01970 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01971 %                                                                             %
01972 %                                                                             %
01973 %                                                                             %
01974 %     E d g e I m a g e                                                       %
01975 %                                                                             %
01976 %                                                                             %
01977 %                                                                             %
01978 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01979 %
01980 %  EdgeImage() finds edges in an image.  Radius defines the radius of the
01981 %  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
01982 %  radius for you.
01983 %
01984 %  The format of the EdgeImage method is:
01985 %
01986 %      Image *EdgeImage(const Image *image,const double radius,
01987 %        ExceptionInfo *exception)
01988 %
01989 %  A description of each parameter follows:
01990 %
01991 %    o image: the image.
01992 %
01993 %    o radius: the radius of the pixel neighborhood.
01994 %
01995 %    o exception: return any errors or warnings in this structure.
01996 %
01997 */
01998 MagickExport Image *EdgeImage(const Image *image,const double radius,
01999   ExceptionInfo *exception)
02000 {
02001   Image
02002     *edge_image;
02003 
02004   double
02005     *kernel;
02006 
02007   register long
02008     i;
02009 
02010   unsigned long
02011     width;
02012 
02013   assert(image != (const Image *) NULL);
02014   assert(image->signature == MagickSignature);
02015   if (image->debug != MagickFalse)
02016     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
02017   assert(exception != (ExceptionInfo *) NULL);
02018   assert(exception->signature == MagickSignature);
02019   width=GetOptimalKernelWidth1D(radius,0.5);
02020   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
02021   if (kernel == (double *) NULL)
02022     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
02023   for (i=0; i < (long) (width*width); i++)
02024     kernel[i]=(-1.0);
02025   kernel[i/2]=(double) (width*width-1.0);
02026   edge_image=ConvolveImage(image,width,kernel,exception);
02027   kernel=(double *) RelinquishMagickMemory(kernel);
02028   return(edge_image);
02029 }
02030 
02031 /*
02032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02033 %                                                                             %
02034 %                                                                             %
02035 %                                                                             %
02036 %     E m b o s s I m a g e                                                   %
02037 %                                                                             %
02038 %                                                                             %
02039 %                                                                             %
02040 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02041 %
02042 %  EmbossImage() returns a grayscale image with a three-dimensional effect.
02043 %  We convolve the image with a Gaussian operator of the given radius and
02044 %  standard deviation (sigma).  For reasonable results, radius should be
02045 %  larger than sigma.  Use a radius of 0 and Emboss() selects a suitable
02046 %  radius for you.
02047 %
02048 %  The format of the EmbossImage method is:
02049 %
02050 %      Image *EmbossImage(const Image *image,const double radius,
02051 %        const double sigma,ExceptionInfo *exception)
02052 %
02053 %  A description of each parameter follows:
02054 %
02055 %    o image: the image.
02056 %
02057 %    o radius: the radius of the pixel neighborhood.
02058 %
02059 %    o sigma: the standard deviation of the Gaussian, in pixels.
02060 %
02061 %    o exception: return any errors or warnings in this structure.
02062 %
02063 */
02064 MagickExport Image *EmbossImage(const Image *image,const double radius,
02065   const double sigma,ExceptionInfo *exception)
02066 {
02067   double
02068     *kernel;
02069 
02070   Image
02071     *emboss_image;
02072 
02073   long
02074     j,
02075     k,
02076     u,
02077     v;
02078 
02079   register long
02080     i;
02081 
02082   unsigned long
02083     width;
02084 
02085   assert(image != (Image *) NULL);
02086   assert(image->signature == MagickSignature);
02087   if (image->debug != MagickFalse)
02088     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
02089   assert(exception != (ExceptionInfo *) NULL);
02090   assert(exception->signature == MagickSignature);
02091   width=GetOptimalKernelWidth2D(radius,sigma);
02092   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
02093   if (kernel == (double *) NULL)
02094     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
02095   j=(long) width/2;
02096   k=j;
02097   i=0;
02098   for (v=(-j); v <= j; v++)
02099   {
02100     for (u=(-j); u <= j; u++)
02101     {
02102       kernel[i]=((u < 0) || (v < 0) ? -8.0 : 8.0)*
02103         exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
02104         (2.0*MagickPI*MagickSigma*MagickSigma);
02105       if (u != k)
02106         kernel[i]=0.0;
02107       i++;
02108     }
02109     k--;
02110   }
02111   emboss_image=ConvolveImage(image,width,kernel,exception);
02112   if (emboss_image != (Image *) NULL)
02113     (void) EqualizeImage(emboss_image);
02114   kernel=(double *) RelinquishMagickMemory(kernel);
02115   return(emboss_image);
02116 }
02117 
02118 /*
02119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02120 %                                                                             %
02121 %                                                                             %
02122 %                                                                             %
02123 %     F i l t e r I m a g e                                                   %
02124 %                                                                             %
02125 %                                                                             %
02126 %                                                                             %
02127 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02128 %
02129 %  FilterImage() applies a custom convolution kernel to the image.
02130 %
02131 %  The format of the FilterImage method is:
02132 %
02133 %      Image *FilterImage(const Image *image,const KernelInfo *kernel,
02134 %        ExceptionInfo *exception)
02135 %      Image *FilterImageChannel(const Image *image,const ChannelType channel,
02136 %        const KernelInfo *kernel,ExceptionInfo *exception)
02137 %
02138 %  A description of each parameter follows:
02139 %
02140 %    o image: the image.
02141 %
02142 %    o channel: the channel type.
02143 %
02144 %    o kernel: the filtering kernel.
02145 %
02146 %    o exception: return any errors or warnings in this structure.
02147 %
02148 */
02149 
02150 MagickExport Image *FilterImage(const Image *image,const KernelInfo *kernel,
02151   ExceptionInfo *exception)
02152 {
02153   Image
02154     *filter_image;
02155 
02156   filter_image=FilterImageChannel(image,DefaultChannels,kernel,exception);
02157   return(filter_image);
02158 }
02159 
02160 MagickExport Image *FilterImageChannel(const Image *image,
02161   const ChannelType channel,const KernelInfo *kernel,ExceptionInfo *exception)
02162 {
02163 #define FilterImageTag  "Filter/Image"
02164 
02165   CacheView
02166     *filter_view,
02167     *image_view;
02168 
02169   Image
02170     *filter_image;
02171 
02172   long
02173     progress,
02174     y;
02175 
02176   MagickBooleanType
02177     status;
02178 
02179   MagickPixelPacket
02180     bias;
02181 
02182   /*
02183     Initialize filter image attributes.
02184   */
02185   assert(image != (Image *) NULL);
02186   assert(image->signature == MagickSignature);
02187   if (image->debug != MagickFalse)
02188     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
02189   assert(exception != (ExceptionInfo *) NULL);
02190   assert(exception->signature == MagickSignature);
02191   if ((kernel->width % 2) == 0)
02192     ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
02193   filter_image=CloneImage(image,0,0,MagickTrue,exception);
02194   if (filter_image == (Image *) NULL)
02195     return((Image *) NULL);
02196   if (SetImageStorageClass(filter_image,DirectClass) == MagickFalse)
02197     {
02198       InheritException(exception,&filter_image->exception);
02199       filter_image=DestroyImage(filter_image);
02200       return((Image *) NULL);
02201     }
02202   if (image->debug != MagickFalse)
02203     {
02204       char
02205         format[MaxTextExtent],
02206         *message;
02207 
02208       long
02209         u,
02210         v;
02211 
02212       register const double
02213         *k;
02214 
02215       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
02216         "  FilterImage with %ldx%ld kernel:",kernel->width,kernel->height);
02217       message=AcquireString("");
02218       k=kernel->values;
02219       for (v=0; v < (long) kernel->height; v++)
02220       {
02221         *message='\0';
02222         (void) FormatMagickString(format,MaxTextExtent,"%ld: ",v);
02223         (void) ConcatenateString(&message,format);
02224         for (u=0; u < (long) kernel->width; u++)
02225         {
02226           (void) FormatMagickString(format,MaxTextExtent,"%g ",*k++);
02227           (void) ConcatenateString(&message,format);
02228         }
02229         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
02230       }
02231       message=DestroyString(message);
02232     }
02233   status=AccelerateConvolveImage(image,kernel,filter_image,exception);
02234   if (status == MagickTrue)
02235     return(filter_image);
02236   /*
02237     Filter image.
02238   */
02239   status=MagickTrue;
02240   progress=0;
02241   GetMagickPixelPacket(image,&bias);
02242   SetMagickPixelPacketBias(image,&bias);
02243   image_view=AcquireCacheView(image);
02244   filter_view=AcquireCacheView(filter_image);
02245 #if defined(MAGICKCORE_OPENMP_SUPPORT)
02246   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
02247 #endif
02248   for (y=0; y < (long) image->rows; y++)
02249   {
02250     MagickBooleanType
02251       sync;
02252 
02253     register const IndexPacket
02254       *restrict indexes;
02255 
02256     register const PixelPacket
02257       *restrict p;
02258 
02259     register IndexPacket
02260       *restrict filter_indexes;
02261 
02262     register long
02263       x;
02264 
02265     register PixelPacket
02266       *restrict q;
02267 
02268     if (status == MagickFalse)
02269       continue;
02270     p=GetCacheViewVirtualPixels(image_view,-((long) kernel->width/2L),
02271       y-(long) (kernel->height/2L),image->columns+kernel->width,kernel->height,
02272       exception);
02273     q=GetCacheViewAuthenticPixels(filter_view,0,y,filter_image->columns,1,
02274       exception);
02275     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
02276       {
02277         status=MagickFalse;
02278         continue;
02279       }
02280     indexes=GetCacheViewVirtualIndexQueue(image_view);
02281     filter_indexes=GetCacheViewAuthenticIndexQueue(filter_view);
02282     for (x=0; x < (long) image->columns; x++)
02283     {
02284       long
02285         v;
02286 
02287       MagickPixelPacket
02288         pixel;
02289 
02290       register const double
02291         *restrict k;
02292 
02293       register const PixelPacket
02294         *restrict kernel_pixels;
02295 
02296       register long
02297         u;
02298 
02299       pixel=bias;
02300       k=kernel->values;
02301       kernel_pixels=p;
02302       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
02303         {
02304           for (v=0; v < (long) kernel->width; v++)
02305           {
02306             for (u=0; u < (long) kernel->height; u++)
02307             {
02308               pixel.red+=(*k)*kernel_pixels[u].red;
02309               pixel.green+=(*k)*kernel_pixels[u].green;
02310               pixel.blue+=(*k)*kernel_pixels[u].blue;
02311               k++;
02312             }
02313             kernel_pixels+=image->columns+kernel->width;
02314           }
02315           if ((channel & RedChannel) != 0)
02316             SetRedPixelComponent(q,ClampRedPixelComponent(&pixel));
02317           if ((channel & GreenChannel) != 0)
02318             SetGreenPixelComponent(q,ClampGreenPixelComponent(&pixel));
02319           if ((channel & BlueChannel) != 0)
02320             SetBluePixelComponent(q,ClampBluePixelComponent(&pixel));
02321           if ((channel & OpacityChannel) != 0)
02322             {
02323               k=kernel->values;
02324               kernel_pixels=p;
02325               for (v=0; v < (long) kernel->width; v++)
02326               {
02327                 for (u=0; u < (long) kernel->height; u++)
02328                 {
02329                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
02330                   k++;
02331                 }
02332                 kernel_pixels+=image->columns+kernel->width;
02333               }
02334               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
02335             }
02336           if (((channel & IndexChannel) != 0) &&
02337               (image->colorspace == CMYKColorspace))
02338             {
02339               register const IndexPacket
02340                 *restrict kernel_indexes;
02341 
02342               k=kernel->values;
02343               kernel_indexes=indexes;
02344               for (v=0; v < (long) kernel->width; v++)
02345               {
02346                 for (u=0; u < (long) kernel->height; u++)
02347                 {
02348                   pixel.index+=(*k)*kernel_indexes[u];
02349                   k++;
02350                 }
02351                 kernel_indexes+=image->columns+kernel->width;
02352               }
02353               filter_indexes[x]=ClampToQuantum(pixel.index);
02354             }
02355         }
02356       else
02357         {
02358           MagickRealType
02359             alpha,
02360             gamma;
02361 
02362           gamma=0.0;
02363           for (v=0; v < (long) kernel->width; v++)
02364           {
02365             for (u=0; u < (long) kernel->height; u++)
02366             {
02367               alpha=(MagickRealType) (QuantumScale*(QuantumRange-
02368                 kernel_pixels[u].opacity));
02369               pixel.red+=(*k)*alpha*kernel_pixels[u].red;
02370               pixel.green+=(*k)*alpha*kernel_pixels[u].green;
02371               pixel.blue+=(*k)*alpha*kernel_pixels[u].blue;
02372               gamma+=(*k)*alpha;
02373               k++;
02374             }
02375             kernel_pixels+=image->columns+kernel->width;
02376           }
02377           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
02378           if ((channel & RedChannel) != 0)
02379             q->red=ClampToQuantum(gamma*GetRedPixelComponent(&pixel));
02380           if ((channel & GreenChannel) != 0)
02381             q->green=ClampToQuantum(gamma*GetGreenPixelComponent(&pixel));
02382           if ((channel & BlueChannel) != 0)
02383             q->blue=ClampToQuantum(gamma*GetBluePixelComponent(&pixel));
02384           if ((channel & OpacityChannel) != 0)
02385             {
02386               k=kernel->values;
02387               kernel_pixels=p;
02388               for (v=0; v < (long) kernel->width; v++)
02389               {
02390                 for (u=0; u < (long) kernel->height; u++)
02391                 {
02392                   pixel.opacity+=(*k)*kernel_pixels[u].opacity;
02393                   k++;
02394                 }
02395                 kernel_pixels+=image->columns+kernel->width;
02396               }
02397               SetOpacityPixelComponent(q,ClampOpacityPixelComponent(&pixel));
02398             }
02399           if (((channel & IndexChannel) != 0) &&
02400               (image->colorspace == CMYKColorspace))
02401             {
02402               register const IndexPacket
02403                 *restrict kernel_indexes;
02404 
02405               k=kernel->values;
02406               kernel_pixels=p;
02407               kernel_indexes=indexes;
02408               for (v=0; v < (long) kernel->width; v++)
02409               {
02410                 for (u=0; u < (long) kernel->height; u++)
02411                 {
02412                   alpha=(MagickRealType) (QuantumScale*(QuantumRange-
02413                     kernel_pixels[u].opacity));
02414                   pixel.index+=(*k)*alpha*kernel_indexes[u];
02415                   k++;
02416                 }
02417                 kernel_pixels+=image->columns+kernel->width;
02418                 kernel_indexes+=image->columns+kernel->width;
02419               }
02420               filter_indexes[x]=ClampToQuantum(gamma*
02421                 GetIndexPixelComponent(&pixel));
02422             }
02423         }
02424       p++;
02425       q++;
02426     }
02427     sync=SyncCacheViewAuthenticPixels(filter_view,exception);
02428     if (sync == MagickFalse)
02429       status=MagickFalse;
02430     if (image->progress_monitor != (MagickProgressMonitor) NULL)
02431       {
02432         MagickBooleanType
02433           proceed;
02434 
02435 #if defined(MAGICKCORE_OPENMP_SUPPORT)
02436   #pragma omp critical (MagickCore_FilterImageChannel)
02437 #endif
02438         proceed=SetImageProgress(image,FilterImageTag,progress++,image->rows);
02439         if (proceed == MagickFalse)
02440           status=MagickFalse;
02441       }
02442   }
02443   filter_image->type=image->type;
02444   filter_view=DestroyCacheView(filter_view);
02445   image_view=DestroyCacheView(image_view);
02446   if (status == MagickFalse)
02447     filter_image=DestroyImage(filter_image);
02448   return(filter_image);
02449 }
02450 
02451 /*
02452 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02453 %                                                                             %
02454 %                                                                             %
02455 %                                                                             %
02456 %     G a u s s i a n B l u r I m a g e                                       %
02457 %                                                                             %
02458 %                                                                             %
02459 %                                                                             %
02460 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02461 %
02462 %  GaussianBlurImage() blurs an image.  We convolve the image with a
02463 %  Gaussian operator of the given radius and standard deviation (sigma).
02464 %  For reasonable results, the radius should be larger than sigma.  Use a
02465 %  radius of 0 and GaussianBlurImage() selects a suitable radius for you
02466 %
02467 %  The format of the GaussianBlurImage method is:
02468 %
02469 %      Image *GaussianBlurImage(const Image *image,onst double radius,
02470 %        const double sigma,ExceptionInfo *exception)
02471 %      Image *GaussianBlurImageChannel(const Image *image,
02472 %        const ChannelType channel,const double radius,const double sigma,
02473 %        ExceptionInfo *exception)
02474 %
02475 %  A description of each parameter follows:
02476 %
02477 %    o image: the image.
02478 %
02479 %    o channel: the channel type.
02480 %
02481 %    o radius: the radius of the Gaussian, in pixels, not counting the center
02482 %      pixel.
02483 %
02484 %    o sigma: the standard deviation of the Gaussian, in pixels.
02485 %
02486 %    o exception: return any errors or warnings in this structure.
02487 %
02488 */
02489 
02490 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
02491   const double sigma,ExceptionInfo *exception)
02492 {
02493   Image
02494     *blur_image;
02495 
02496   blur_image=GaussianBlurImageChannel(image,DefaultChannels,radius,sigma,
02497     exception);
02498   return(blur_image);
02499 }
02500 
02501 MagickExport Image *GaussianBlurImageChannel(const Image *image,
02502   const ChannelType channel,const double radius,const double sigma,
02503   ExceptionInfo *exception)
02504 {
02505   double
02506     *kernel;
02507 
02508   Image
02509     *blur_image;
02510 
02511   long
02512     j,
02513     u,
02514     v;
02515 
02516   register long
02517     i;
02518 
02519   unsigned long
02520     width;
02521 
02522   assert(image != (const Image *) NULL);
02523   assert(image->signature == MagickSignature);
02524   if (image->debug != MagickFalse)
02525     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
02526   assert(exception != (ExceptionInfo *) NULL);
02527   assert(exception->signature == MagickSignature);
02528   width=GetOptimalKernelWidth2D(radius,sigma);
02529   kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
02530   if (kernel == (double *) NULL)
02531     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
02532   j=(long) width/2;
02533   i=0;
02534   for (v=(-j); v <= j; v++)
02535   {
02536     for (u=(-j); u <= j; u++)
02537       kernel[i++]=exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
02538         (2.0*MagickPI*MagickSigma*MagickSigma);
02539   }
02540   blur_image=ConvolveImageChannel(image,channel,width,kernel,exception);
02541   kernel=(double *) RelinquishMagickMemory(kernel);
02542   return(blur_image);
02543 }
02544 
02545 /*
02546 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02547 %                                                                             %
02548 %                                                                             %
02549 %                                                                             %
02550 %     M e d i a n F i l t e r I m a g e                                       %
02551 %                                                                             %
02552 %                                                                             %
02553 %                                                                             %
02554 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02555 %
02556 %  MedianFilterImage() applies a digital filter that improves the quality
02557 %  of a noisy image.  Each pixel is replaced by the median in a set of
02558 %  neighboring pixels as defined by radius.
02559 %
02560 %  The algorithm was contributed by Mike Edmonds and implements an insertion
02561 %  sort for selecting median color-channel values.  For more on this algorithm
02562 %  see "Skip Lists: A probabilistic Alternative to Balanced Trees" by William
02563 %  Pugh in the June 1990 of Communications of the ACM.
02564 %
02565 %  The format of the MedianFilterImage method is:
02566 %
02567 %      Image *MedianFilterImage(const Image *image,const double radius,
02568 %        ExceptionInfo *exception)
02569 %
02570 %  A description of each parameter follows:
02571 %
02572 %    o image: the image.
02573 %
02574 %    o radius: the radius of the pixel neighborhood.
02575 %
02576 %    o exception: return any errors or warnings in this structure.
02577 %
02578 */
02579 
02580 #define MedianListChannels  5
02581 
02582 typedef struct _MedianListNode
02583 {
02584   unsigned long
02585     next[9],
02586     count,
02587     signature;
02588 } MedianListNode;
02589 
02590 typedef struct _MedianSkipList
02591 {
02592   long
02593     level;
02594 
02595   MedianListNode
02596     *nodes;
02597 } MedianSkipList;
02598 
02599 typedef struct _MedianPixelList
02600 {
02601   unsigned long
02602     center,
02603     seed,
02604     signature;
02605 
02606   MedianSkipList
02607     lists[MedianListChannels];
02608 } MedianPixelList;
02609 
02610 static MedianPixelList *DestroyMedianPixelList(MedianPixelList *pixel_list)
02611 {
02612   register long
02613     i;
02614 
02615   if (pixel_list == (MedianPixelList *) NULL)
02616     return((MedianPixelList *) NULL);
02617   for (i=0; i < MedianListChannels; i++)
02618     if (pixel_list->lists[i].nodes != (MedianListNode *) NULL)
02619       pixel_list->lists[i].nodes=(MedianListNode *) RelinquishMagickMemory(
02620         pixel_list->lists[i].nodes);
02621   pixel_list=(MedianPixelList *) RelinquishAlignedMemory(pixel_list);
02622   return(pixel_list);
02623 }
02624 
02625 static MedianPixelList **DestroyMedianPixelListThreadSet(
02626   MedianPixelList **pixel_list)
02627 {
02628   register long
02629     i;
02630 
02631   assert(pixel_list != (MedianPixelList **) NULL);
02632   for (i=0; i < (long) GetOpenMPMaximumThreads(); i++)
02633     if (pixel_list[i] != (MedianPixelList *) NULL)
02634       pixel_list[i]=DestroyMedianPixelList(pixel_list[i]);
02635   pixel_list=(MedianPixelList **) RelinquishAlignedMemory(pixel_list);
02636   return(pixel_list);
02637 }
02638 
02639 static MedianPixelList *AcquireMedianPixelList(const unsigned long width)
02640 {
02641   MedianPixelList
02642     *pixel_list;
02643 
02644   register long
02645     i;
02646 
02647   pixel_list=(MedianPixelList *) AcquireAlignedMemory(1,sizeof(*pixel_list));
02648   if (pixel_list == (MedianPixelList *) NULL)
02649     return(pixel_list);
02650   (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
02651   pixel_list->center=width*width/2;
02652   for (i=0; i < MedianListChannels; i++)
02653   {
02654     pixel_list->lists[i].nodes=(MedianListNode *) AcquireQuantumMemory(65537UL,
02655       sizeof(*pixel_list->lists[i].nodes));
02656     if (pixel_list->lists[i].nodes == (MedianListNode *) NULL)
02657       return(DestroyMedianPixelList(pixel_list));
02658     (void) ResetMagickMemory(pixel_list->lists[i].nodes,0,65537UL*
02659       sizeof(*pixel_list->lists[i].nodes));
02660   }
02661   pixel_list->signature=MagickSignature;
02662   return(pixel_list);
02663 }
02664 
02665 static MedianPixelList **AcquireMedianPixelListThreadSet(
02666   const unsigned long width)
02667 {
02668   register long
02669     i;
02670 
02671   MedianPixelList
02672     **pixel_list;
02673 
02674   unsigned long
02675     number_threads;
02676 
02677   number_threads=GetOpenMPMaximumThreads();
02678   pixel_list=(MedianPixelList **) AcquireAlignedMemory(number_threads,
02679     sizeof(*pixel_list));
02680   if (pixel_list == (MedianPixelList **) NULL)
02681     return((MedianPixelList **) NULL);
02682   (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
02683   for (i=0; i < (long) number_threads; i++)
02684   {
02685     pixel_list[i]=AcquireMedianPixelList(width);
02686     if (pixel_list[i] == (MedianPixelList *) NULL)
02687       return(DestroyMedianPixelListThreadSet(pixel_list));
02688   }
02689   return(pixel_list);
02690 }
02691 
02692 static void AddNodeMedianPixelList(MedianPixelList *pixel_list,
02693   const long channel,const unsigned long color)
02694 {
02695   register long
02696     level;
02697 
02698   register MedianSkipList
02699     *list;
02700 
02701   unsigned long
02702     search,
02703     update[9];
02704 
02705   /*
02706     Initialize the node.
02707   */
02708   list=pixel_list->lists+channel;
02709   list->nodes[color].signature=pixel_list->signature;
02710   list->nodes[color].count=1;
02711   /*
02712     Determine where it belongs in the list.
02713   */
02714   search=65536UL;
02715   for (level=list->level; level >= 0; level--)
02716   {
02717     while (list->nodes[search].next[level] < color)
02718       search=list->nodes[search].next[level];
02719     update[level]=search;
02720   }
02721   /*
02722     Generate a pseudo-random level for this node.
02723   */
02724   for (level=0; ; level++)
02725   {
02726     pixel_list->seed=(pixel_list->seed*42893621L)+1L;
02727     if ((pixel_list->seed & 0x300) != 0x300)
02728       break;
02729   }
02730   if (level > 8)
02731     level=8;
02732   if (level > (list->level+2))
02733     level=list->level+2;
02734   /*
02735     If we're raising the list's level, link back to the root node.
02736   */
02737   while (level > list->level)
02738   {
02739     list->level++;
02740     update[list->level]=65536UL;
02741   }
02742   /*
02743     Link the node into the skip-list.
02744   */
02745   do
02746   {
02747     list->nodes[color].next[level]=list->nodes[update[level]].next[level];
02748     list->nodes[update[level]].next[level]=color;
02749   }
02750   while (level-- > 0);
02751 }
02752 
02753 static MagickPixelPacket GetMedianPixelList(MedianPixelList *pixel_list)
02754 {
02755   MagickPixelPacket
02756     pixel;
02757 
02758   register long
02759     channel;
02760 
02761   register MedianSkipList
02762     *list;
02763 
02764   unsigned long
02765     center,
02766     color,
02767     count;
02768 
02769   unsigned short
02770     channels[MedianListChannels];
02771 
02772   /*
02773     Find the median value for each of the color.
02774   */
02775   center=pixel_list->center;
02776   for (channel=0; channel < 5; channel++)
02777   {
02778     list=pixel_list->lists+channel;
02779     color=65536UL;
02780     count=0;
02781     do
02782     {
02783       color=list->nodes[color].next[0];
02784       count+=list->nodes[color].count;
02785     }
02786     while (count <= center);
02787     channels[channel]=(unsigned short) color;
02788   }
02789   GetMagickPixelPacket((const Image *) NULL,&pixel);
02790   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
02791   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
02792   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
02793   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
02794   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
02795   return(pixel);
02796 }
02797 
02798 static inline void InsertMedianPixelList(const Image *image,
02799   const PixelPacket *pixel,const IndexPacket *indexes,
02800   MedianPixelList *pixel_list)
02801 {
02802   unsigned long
02803     signature;
02804 
02805   unsigned short
02806     index;
02807 
02808   index=ScaleQuantumToShort(pixel->red);
02809   signature=pixel_list->lists[0].nodes[index].signature;
02810   if (signature == pixel_list->signature)
02811     pixel_list->lists[0].nodes[index].count++;
02812   else
02813     AddNodeMedianPixelList(pixel_list,0,index);
02814   index=ScaleQuantumToShort(pixel->green);
02815   signature=pixel_list->lists[1].nodes[index].signature;
02816   if (signature == pixel_list->signature)
02817     pixel_list->lists[1].nodes[index].count++;
02818   else
02819     AddNodeMedianPixelList(pixel_list,1,index);
02820   index=ScaleQuantumToShort(pixel->blue);
02821   signature=pixel_list->lists[2].nodes[index].signature;
02822   if (signature == pixel_list->signature)
02823     pixel_list->lists[2].nodes[index].count++;
02824   else
02825     AddNodeMedianPixelList(pixel_list,2,index);
02826   index=ScaleQuantumToShort(pixel->opacity);
02827   signature=pixel_list->lists[3].nodes[index].signature;
02828   if (signature == pixel_list->signature)
02829     pixel_list->lists[3].nodes[index].count++;
02830   else
02831     AddNodeMedianPixelList(pixel_list,3,index);
02832   if (image->colorspace == CMYKColorspace)
02833     index=ScaleQuantumToShort(*indexes);
02834   signature=pixel_list->lists[4].nodes[index].signature;
02835   if (signature == pixel_list->signature)
02836     pixel_list->lists[4].nodes[index].count++;
02837   else
02838     AddNodeMedianPixelList(pixel_list,4,index);
02839 }
02840 
02841 static void ResetMedianPixelList(MedianPixelList *pixel_list)
02842 {
02843   int
02844     level;
02845 
02846   register long
02847     channel;
02848 
02849   register MedianListNode
02850     *root;
02851 
02852   register MedianSkipList
02853     *list;
02854 
02855   /*
02856     Reset the skip-list.
02857   */
02858   for (channel=0; channel < 5; channel++)
02859   {
02860     list=pixel_list->lists+channel;
02861     root=list->nodes+65536UL;
02862     list->level=0;
02863     for (level=0; level < 9; level++)
02864       root->next[level]=65536UL;
02865   }
02866   pixel_list->seed=pixel_list->signature++;
02867 }
02868 
02869 MagickExport Image *MedianFilterImage(const Image *image,const double radius,
02870   ExceptionInfo *exception)
02871 {
02872 #define MedianFilterImageTag  "MedianFilter/Image"
02873 
02874   CacheView
02875     *image_view,
02876     *median_view;
02877 
02878   Image
02879     *median_image;
02880 
02881   long
02882     progress,
02883     y;
02884 
02885   MagickBooleanType
02886     status;
02887 
02888   MedianPixelList
02889     **restrict pixel_list;
02890 
02891   unsigned long
02892     width;
02893 
02894   /*
02895     Initialize median image attributes.
02896   */
02897   assert(image != (Image *) NULL);
02898   assert(image->signature == MagickSignature);
02899   if (image->debug != MagickFalse)
02900     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
02901   assert(exception != (ExceptionInfo *) NULL);
02902   assert(exception->signature == MagickSignature);
02903   width=GetOptimalKernelWidth2D(radius,0.5);
02904   if ((image->columns < width) || (image->rows < width))
02905     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
02906   median_image=CloneImage(image,image->columns,image->rows,MagickTrue,
02907     exception);
02908   if (median_image == (Image *) NULL)
02909     return((Image *) NULL);
02910   if (SetImageStorageClass(median_image,DirectClass) == MagickFalse)
02911     {
02912       InheritException(exception,&median_image->exception);
02913       median_image=DestroyImage(median_image);
02914       return((Image *) NULL);
02915     }
02916   pixel_list=AcquireMedianPixelListThreadSet(width);
02917   if (pixel_list == (MedianPixelList **) NULL)
02918     {
02919       median_image=DestroyImage(median_image);
02920       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
02921     }
02922   /*
02923     Median filter each image row.
02924   */
02925   status=MagickTrue;
02926   progress=0;
02927   image_view=AcquireCacheView(image);
02928   median_view=AcquireCacheView(median_image);
02929 #if defined(MAGICKCORE_OPENMP_SUPPORT)
02930   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
02931 #endif
02932   for (y=0; y < (long) median_image->rows; y++)
02933   {
02934     register const IndexPacket
02935       *restrict indexes;
02936 
02937     register const PixelPacket
02938       *restrict p;
02939 
02940     register IndexPacket
02941       *restrict median_indexes;
02942 
02943     register long
02944       id,
02945       x;
02946 
02947     register PixelPacket
02948       *restrict q;
02949 
02950     if (status == MagickFalse)
02951       continue;
02952     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
02953       2L),image->columns+width,width,exception);
02954     q=QueueCacheViewAuthenticPixels(median_view,0,y,median_image->columns,1,
02955       exception);
02956     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
02957       {
02958         status=MagickFalse;
02959         continue;
02960       }
02961     indexes=GetCacheViewVirtualIndexQueue(image_view);
02962     median_indexes=GetCacheViewAuthenticIndexQueue(median_view);
02963     id=GetOpenMPThreadId();
02964     for (x=0; x < (long) median_image->columns; x++)
02965     {
02966       MagickPixelPacket
02967         pixel;
02968 
02969       register const PixelPacket
02970         *restrict r;
02971 
02972       register const IndexPacket
02973         *restrict s;
02974 
02975       register long
02976         u,
02977         v;
02978 
02979       r=p;
02980       s=indexes+x;
02981       ResetMedianPixelList(pixel_list[id]);
02982       for (v=0; v < (long) width; v++)
02983       {
02984         for (u=0; u < (long) width; u++)
02985           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
02986         r+=image->columns+width;
02987         s+=image->columns+width;
02988       }
02989       pixel=GetMedianPixelList(pixel_list[id]);
02990       SetPixelPacket(median_image,&pixel,q,median_indexes+x);
02991       p++;
02992       q++;
02993     }
02994     if (SyncCacheViewAuthenticPixels(median_view,exception) == MagickFalse)
02995       status=MagickFalse;
02996     if (image->progress_monitor != (MagickProgressMonitor) NULL)
02997       {
02998         MagickBooleanType
02999           proceed;
03000 
03001 #if defined(MAGICKCORE_OPENMP_SUPPORT)
03002   #pragma omp critical (MagickCore_MedianFilterImage)
03003 #endif
03004         proceed=SetImageProgress(image,MedianFilterImageTag,progress++,
03005           image->rows);
03006         if (proceed == MagickFalse)
03007           status=MagickFalse;
03008       }
03009   }
03010   median_view=DestroyCacheView(median_view);
03011   image_view=DestroyCacheView(image_view);
03012   pixel_list=DestroyMedianPixelListThreadSet(pixel_list);
03013   return(median_image);
03014 }
03015 
03016 /*
03017 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03018 %                                                                             %
03019 %                                                                             %
03020 %                                                                             %
03021 %     M o t i o n B l u r I m a g e                                           %
03022 %                                                                             %
03023 %                                                                             %
03024 %                                                                             %
03025 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03026 %
03027 %  MotionBlurImage() simulates motion blur.  We convolve the image with a
03028 %  Gaussian operator of the given radius and standard deviation (sigma).
03029 %  For reasonable results, radius should be larger than sigma.  Use a
03030 %  radius of 0 and MotionBlurImage() selects a suitable radius for you.
03031 %  Angle gives the angle of the blurring motion.
03032 %
03033 %  Andrew Protano contributed this effect.
03034 %
03035 %  The format of the MotionBlurImage method is:
03036 %
03037 %    Image *MotionBlurImage(const Image *image,const double radius,
03038 %      const double sigma,const double angle,ExceptionInfo *exception)
03039 %    Image *MotionBlurImageChannel(const Image *image,const ChannelType channel,
03040 %      const double radius,const double sigma,const double angle,
03041 %      ExceptionInfo *exception)
03042 %
03043 %  A description of each parameter follows:
03044 %
03045 %    o image: the image.
03046 %
03047 %    o channel: the channel type.
03048 %
03049 %    o radius: the radius of the Gaussian, in pixels, not counting the center
03050 %    o radius: the radius of the Gaussian, in pixels, not counting
03051 %      the center pixel.
03052 %
03053 %    o sigma: the standard deviation of the Gaussian, in pixels.
03054 %
03055 %    o angle: Apply the effect along this angle.
03056 %
03057 %    o exception: return any errors or warnings in this structure.
03058 %
03059 */
03060 
03061 static double *GetMotionBlurKernel(const unsigned long width,const double sigma)
03062 {
03063   double
03064     *kernel,
03065     normalize;
03066 
03067   register long
03068     i;
03069 
03070   /*
03071    Generate a 1-D convolution kernel.
03072   */
03073   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
03074   kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
03075   if (kernel == (double *) NULL)
03076     return(kernel);
03077   normalize=0.0;
03078   for (i=0; i < (long) width; i++)
03079   {
03080     kernel[i]=exp((-((double) i*i)/(double) (2.0*MagickSigma*MagickSigma)))/
03081       (MagickSQ2PI*MagickSigma);
03082     normalize+=kernel[i];
03083   }
03084   for (i=0; i < (long) width; i++)
03085     kernel[i]/=normalize;
03086   return(kernel);
03087 }
03088 
03089 MagickExport Image *MotionBlurImage(const Image *image,const double radius,
03090   const double sigma,const double angle,ExceptionInfo *exception)
03091 {
03092   Image
03093     *motion_blur;
03094 
03095   motion_blur=MotionBlurImageChannel(image,DefaultChannels,radius,sigma,angle,
03096     exception);
03097   return(motion_blur);
03098 }
03099 
03100 MagickExport Image *MotionBlurImageChannel(const Image *image,
03101   const ChannelType channel,const double radius,const double sigma,
03102   const double angle,ExceptionInfo *exception)
03103 {
03104   CacheView
03105     *blur_view,
03106     *image_view;
03107 
03108   double
03109     *kernel;
03110 
03111   Image
03112     *blur_image;
03113 
03114   long
03115     progress,
03116     y;
03117 
03118   MagickBooleanType
03119     status;
03120 
03121   MagickPixelPacket
03122     bias;
03123 
03124   OffsetInfo
03125     *offset;
03126 
03127   PointInfo
03128     point;
03129 
03130   register long
03131     i;
03132 
03133   unsigned long
03134     width;
03135 
03136   assert(image != (Image *) NULL);
03137   assert(image->signature == MagickSignature);
03138   if (image->debug != MagickFalse)
03139     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
03140   assert(exception != (ExceptionInfo *) NULL);
03141   width=GetOptimalKernelWidth1D(radius,sigma);
03142   kernel=GetMotionBlurKernel(width,sigma);
03143   if (kernel == (double *) NULL)
03144     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
03145   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
03146   if (offset == (OffsetInfo *) NULL)
03147     {
03148       kernel=(double *) RelinquishMagickMemory(kernel);
03149       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
03150     }
03151   blur_image=CloneImage(image,0,0,MagickTrue,exception);
03152   if (blur_image == (Image *) NULL)
03153     {
03154       kernel=(double *) RelinquishMagickMemory(kernel);
03155       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
03156       return((Image *) NULL);
03157     }
03158   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
03159     {
03160       kernel=(double *) RelinquishMagickMemory(kernel);
03161       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
03162       InheritException(exception,&blur_image->exception);
03163       blur_image=DestroyImage(blur_image);
03164       return((Image *) NULL);
03165     }
03166   point.x=(double) width*sin(DegreesToRadians(angle));
03167   point.y=(double) width*cos(DegreesToRadians(angle));
03168   for (i=0; i < (long) width; i++)
03169   {
03170     offset[i].x=(long) ceil((i*point.y)/hypot(point.x,point.y)-0.5);
03171     offset[i].y=(long) ceil((i*point.x)/hypot(point.x,point.y)-0.5);
03172   }
03173   /*
03174     Motion blur image.
03175   */
03176   status=MagickTrue;
03177   progress=0;
03178   GetMagickPixelPacket(image,&bias);
03179   image_view=AcquireCacheView(image);
03180   blur_view=AcquireCacheView(blur_image);
03181 #if defined(MAGICKCORE_OPENMP_SUPPORT) && (_OPENMP > 202001)
03182   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
03183 #endif
03184   for (y=0; y < (long) image->rows; y++)
03185   {
03186     register IndexPacket
03187       *restrict blur_indexes;
03188 
03189     register long
03190       x;
03191 
03192     register PixelPacket
03193       *restrict q;
03194 
03195     if (status == MagickFalse)
03196       continue;
03197     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
03198       exception);
03199     if (q == (PixelPacket *) NULL)
03200       {
03201         status=MagickFalse;
03202         continue;
03203       }
03204     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
03205     for (x=0; x < (long) image->columns; x++)
03206     {
03207       MagickPixelPacket
03208         qixel;
03209 
03210       PixelPacket
03211         pixel;
03212 
03213       register double
03214         *restrict k;
03215 
03216       register long
03217         i;
03218 
03219       register const IndexPacket
03220         *restrict indexes;
03221 
03222       k=kernel;
03223       qixel=bias;
03224       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
03225         {
03226           for (i=0; i < (long) width; i++)
03227           {
03228             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
03229               offset[i].y,&pixel,exception);
03230             qixel.red+=(*k)*pixel.red;
03231             qixel.green+=(*k)*pixel.green;
03232             qixel.blue+=(*k)*pixel.blue;
03233             qixel.opacity+=(*k)*pixel.opacity;
03234             if (image->colorspace == CMYKColorspace)
03235               {
03236                 indexes=GetCacheViewVirtualIndexQueue(image_view);
03237                 qixel.index+=(*k)*(*indexes);
03238               }
03239             k++;
03240           }
03241           if ((channel & RedChannel) != 0)
03242             q->red=ClampToQuantum(qixel.red);
03243           if ((channel & GreenChannel) != 0)
03244             q->green=ClampToQuantum(qixel.green);
03245           if ((channel & BlueChannel) != 0)
03246             q->blue=ClampToQuantum(qixel.blue);
03247           if ((channel & OpacityChannel) != 0)
03248             q->opacity=ClampToQuantum(qixel.opacity);
03249           if (((channel & IndexChannel) != 0) &&
03250               (image->colorspace == CMYKColorspace))
03251             blur_indexes[x]=(IndexPacket) ClampToQuantum(qixel.index);
03252         }
03253       else
03254         {
03255           MagickRealType
03256             alpha,
03257             gamma;
03258 
03259           alpha=0.0;
03260           gamma=0.0;
03261           for (i=0; i < (long) width; i++)
03262           {
03263             (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
03264               offset[i].y,&pixel,exception);
03265             alpha=(MagickRealType) (QuantumScale*
03266               GetAlphaPixelComponent(&pixel));
03267             qixel.red+=(*k)*alpha*pixel.red;
03268             qixel.green+=(*k)*alpha*pixel.green;
03269             qixel.blue+=(*k)*alpha*pixel.blue;
03270             qixel.opacity+=(*k)*pixel.opacity;
03271             if (image->colorspace == CMYKColorspace)
03272               {
03273                 indexes=GetCacheViewVirtualIndexQueue(image_view);
03274                 qixel.index+=(*k)*alpha*(*indexes);
03275               }
03276             gamma+=(*k)*alpha;
03277             k++;
03278           }
03279           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
03280           if ((channel & RedChannel) != 0)
03281             q->red=ClampToQuantum(gamma*qixel.red);
03282           if ((channel & GreenChannel) != 0)
03283             q->green=ClampToQuantum(gamma*qixel.green);
03284           if ((channel & BlueChannel) != 0)
03285             q->blue=ClampToQuantum(gamma*qixel.blue);
03286           if ((channel & OpacityChannel) != 0)
03287             q->opacity=ClampToQuantum(qixel.opacity);
03288           if (((channel & IndexChannel) != 0) &&
03289               (image->colorspace == CMYKColorspace))
03290             blur_indexes[x]=(IndexPacket) ClampToQuantum(gamma*qixel.index);
03291         }
03292       q++;
03293     }
03294     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
03295       status=MagickFalse;
03296     if (image->progress_monitor != (MagickProgressMonitor) NULL)
03297       {
03298         MagickBooleanType
03299           proceed;
03300 
03301 #if defined(MAGICKCORE_OPENMP_SUPPORT) && (_OPENMP > 202001)
03302   #pragma omp critical (MagickCore_MotionBlurImageChannel)
03303 #endif
03304         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
03305         if (proceed == MagickFalse)
03306           status=MagickFalse;
03307       }
03308   }
03309   blur_view=DestroyCacheView(blur_view);
03310   image_view=DestroyCacheView(image_view);
03311   kernel=(double *) RelinquishMagickMemory(kernel);
03312   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
03313   if (status == MagickFalse)
03314     blur_image=DestroyImage(blur_image);
03315   return(blur_image);
03316 }
03317 
03318 /*
03319 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03320 %                                                                             %
03321 %                                                                             %
03322 %                                                                             %
03323 %     P r e v i e w I m a g e                                                 %
03324 %                                                                             %
03325 %                                                                             %
03326 %                                                                             %
03327 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03328 %
03329 %  PreviewImage() tiles 9 thumbnails of the specified image with an image
03330 %  processing operation applied with varying parameters.  This may be helpful
03331 %  pin-pointing an appropriate parameter for a particular image processing
03332 %  operation.
03333 %
03334 %  The format of the PreviewImages method is:
03335 %
03336 %      Image *PreviewImages(const Image *image,const PreviewType preview,
03337 %        ExceptionInfo *exception)
03338 %
03339 %  A description of each parameter follows:
03340 %
03341 %    o image: the image.
03342 %
03343 %    o preview: the image processing operation.
03344 %
03345 %    o exception: return any errors or warnings in this structure.
03346 %
03347 */
03348 MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
03349   ExceptionInfo *exception)
03350 {
03351 #define NumberTiles  9
03352 #define PreviewImageTag  "Preview/Image"
03353 #define DefaultPreviewGeometry  "204x204+10+10"
03354 
03355   char
03356     factor[MaxTextExtent],
03357     label[MaxTextExtent];
03358 
03359   double
03360     degrees,
03361     gamma,
03362     percentage,
03363     radius,
03364     sigma,
03365     threshold;
03366 
03367   Image
03368     *images,
03369     *montage_image,
03370     *preview_image,
03371     *thumbnail;
03372 
03373   ImageInfo
03374     *preview_info;
03375 
03376   long
03377     y;
03378 
03379   MagickBooleanType
03380     proceed;
03381 
03382   MontageInfo
03383     *montage_info;
03384 
03385   QuantizeInfo
03386     quantize_info;
03387 
03388   RectangleInfo
03389     geometry;
03390 
03391   register long
03392     i,
03393     x;
03394 
03395   unsigned long
03396     colors;
03397 
03398   /*
03399     Open output image file.
03400   */
03401   assert(image != (Image *) NULL);
03402   assert(image->signature == MagickSignature);
03403   if (image->debug != MagickFalse)
03404     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
03405   colors=2;
03406   degrees=0.0;
03407   gamma=(-0.2f);
03408   preview_info=AcquireImageInfo();
03409   SetGeometry(image,&geometry);
03410   (void) ParseMetaGeometry(DefaultPreviewGeometry,&geometry.x,&geometry.y,
03411     &geometry.width,&geometry.height);
03412   images=NewImageList();
03413   percentage=12.5;
03414   GetQuantizeInfo(&quantize_info);
03415   radius=0.0;
03416   sigma=1.0;
03417   threshold=0.0;
03418   x=0;
03419   y=0;
03420   for (i=0; i < NumberTiles; i++)
03421   {
03422     thumbnail=ThumbnailImage(image,geometry.width,geometry.height,exception);
03423     if (thumbnail == (Image *) NULL)
03424       break;
03425     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
03426       (void *) NULL);
03427     (void) SetImageProperty(thumbnail,"label",DefaultTileLabel);
03428     if (i == (NumberTiles/2))
03429       {
03430         (void) QueryColorDatabase("#dfdfdf",&thumbnail->matte_color,exception);
03431         AppendImageToList(&images,thumbnail);
03432         continue;
03433       }
03434     switch (preview)
03435     {
03436       case RotatePreview:
03437       {
03438         degrees+=45.0;
03439         preview_image=RotateImage(thumbnail,degrees,exception);
03440         (void) FormatMagickString(label,MaxTextExtent,"rotate %g",degrees);
03441         break;
03442       }
03443       case ShearPreview:
03444       {
03445         degrees+=5.0;
03446         preview_image=ShearImage(thumbnail,degrees,degrees,exception);
03447         (void) FormatMagickString(label,MaxTextExtent,"shear %gx%g",
03448           degrees,2.0*degrees);
03449         break;
03450       }
03451       case RollPreview:
03452       {
03453         x=(long) ((i+1)*thumbnail->columns)/NumberTiles;
03454         y=(long) ((i+1)*thumbnail->rows)/NumberTiles;
03455         preview_image=RollImage(thumbnail,x,y,exception);
03456         (void) FormatMagickString(label,MaxTextExtent,"roll %ldx%ld",x,y);
03457         break;
03458       }
03459       case HuePreview:
03460       {
03461         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03462         if (preview_image == (Image *) NULL)
03463           break;
03464         (void) FormatMagickString(factor,MaxTextExtent,"100,100,%g",
03465           2.0*percentage);
03466         (void) ModulateImage(preview_image,factor);
03467         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
03468         break;
03469       }
03470       case SaturationPreview:
03471       {
03472         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03473         if (preview_image == (Image *) NULL)
03474           break;
03475         (void) FormatMagickString(factor,MaxTextExtent,"100,%g",
03476           2.0*percentage);
03477         (void) ModulateImage(preview_image,factor);
03478         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
03479         break;
03480       }
03481       case BrightnessPreview:
03482       {
03483         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03484         if (preview_image == (Image *) NULL)
03485           break;
03486         (void) FormatMagickString(factor,MaxTextExtent,"%g",2.0*percentage);
03487         (void) ModulateImage(preview_image,factor);
03488         (void) FormatMagickString(label,MaxTextExtent,"modulate %s",factor);
03489         break;
03490       }
03491       case GammaPreview:
03492       default:
03493       {
03494         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03495         if (preview_image == (Image *) NULL)
03496           break;
03497         gamma+=0.4f;
03498         (void) GammaImageChannel(preview_image,DefaultChannels,gamma);
03499         (void) FormatMagickString(label,MaxTextExtent,"gamma %g",gamma);
03500         break;
03501       }
03502       case SpiffPreview:
03503       {
03504         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03505         if (preview_image != (Image *) NULL)
03506           for (x=0; x < i; x++)
03507             (void) ContrastImage(preview_image,MagickTrue);
03508         (void) FormatMagickString(label,MaxTextExtent,"contrast (%ld)",i+1);
03509         break;
03510       }
03511       case DullPreview:
03512       {
03513         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03514         if (preview_image == (Image *) NULL)
03515           break;
03516         for (x=0; x < i; x++)
03517           (void) ContrastImage(preview_image,MagickFalse);
03518         (void) FormatMagickString(label,MaxTextExtent,"+contrast (%ld)",i+1);
03519         break;
03520       }
03521       case GrayscalePreview:
03522       {
03523         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03524         if (preview_image == (Image *) NULL)
03525           break;
03526         colors<<=1;
03527         quantize_info.number_colors=colors;
03528         quantize_info.colorspace=GRAYColorspace;
03529         (void) QuantizeImage(&quantize_info,preview_image);
03530         (void) FormatMagickString(label,MaxTextExtent,
03531           "-colorspace gray -colors %ld",colors);
03532         break;
03533       }
03534       case QuantizePreview:
03535       {
03536         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03537         if (preview_image == (Image *) NULL)
03538           break;
03539         colors<<=1;
03540         quantize_info.number_colors=colors;
03541         (void) QuantizeImage(&quantize_info,preview_image);
03542         (void) FormatMagickString(label,MaxTextExtent,"colors %ld",colors);
03543         break;
03544       }
03545       case DespecklePreview:
03546       {
03547         for (x=0; x < (i-1); x++)
03548         {
03549           preview_image=DespeckleImage(thumbnail,exception);
03550           if (preview_image == (Image *) NULL)
03551             break;
03552           thumbnail=DestroyImage(thumbnail);
03553           thumbnail=preview_image;
03554         }
03555         preview_image=DespeckleImage(thumbnail,exception);
03556         if (preview_image == (Image *) NULL)
03557           break;
03558         (void) FormatMagickString(label,MaxTextExtent,"despeckle (%ld)",i+1);
03559         break;
03560       }
03561       case ReduceNoisePreview:
03562       {
03563         preview_image=ReduceNoiseImage(thumbnail,radius,exception);
03564         (void) FormatMagickString(label,MaxTextExtent,"noise %g",radius);
03565         break;
03566       }
03567       case AddNoisePreview:
03568       {
03569         switch ((int) i)
03570         {
03571           case 0:
03572           {
03573             (void) CopyMagickString(factor,"uniform",MaxTextExtent);
03574             break;
03575           }
03576           case 1:
03577           {
03578             (void) CopyMagickString(factor,"gaussian",MaxTextExtent);
03579             break;
03580           }
03581           case 2:
03582           {
03583             (void) CopyMagickString(factor,"multiplicative",MaxTextExtent);
03584             break;
03585           }
03586           case 3:
03587           {
03588             (void) CopyMagickString(factor,"impulse",MaxTextExtent);
03589             break;
03590           }
03591           case 4:
03592           {
03593             (void) CopyMagickString(factor,"laplacian",MaxTextExtent);
03594             break;
03595           }
03596           case 5:
03597           {
03598             (void) CopyMagickString(factor,"Poisson",MaxTextExtent);
03599             break;
03600           }
03601           default:
03602           {
03603             (void) CopyMagickString(thumbnail->magick,"NULL",MaxTextExtent);
03604             break;
03605           }
03606         }
03607         preview_image=ReduceNoiseImage(thumbnail,(double) i,exception);
03608         (void) FormatMagickString(label,MaxTextExtent,"+noise %s",factor);
03609         break;
03610       }
03611       case SharpenPreview:
03612       {
03613         preview_image=SharpenImage(thumbnail,radius,sigma,exception);
03614         (void) FormatMagickString(label,MaxTextExtent,"sharpen %gx%g",
03615           radius,sigma);
03616         break;
03617       }
03618       case BlurPreview:
03619       {
03620         preview_image=BlurImage(thumbnail,radius,sigma,exception);
03621         (void) FormatMagickString(label,MaxTextExtent,"blur %gx%g",radius,
03622           sigma);
03623         break;
03624       }
03625       case ThresholdPreview:
03626       {
03627         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03628         if (preview_image == (Image *) NULL)
03629           break;
03630         (void) BilevelImage(thumbnail,
03631           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
03632         (void) FormatMagickString(label,MaxTextExtent,"threshold %g",
03633           (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
03634         break;
03635       }
03636       case EdgeDetectPreview:
03637       {
03638         preview_image=EdgeImage(thumbnail,radius,exception);
03639         (void) FormatMagickString(label,MaxTextExtent,"edge %g",radius);
03640         break;
03641       }
03642       case SpreadPreview:
03643       {
03644         preview_image=SpreadImage(thumbnail,radius,exception);
03645         (void) FormatMagickString(label,MaxTextExtent,"spread %g",
03646           radius+0.5);
03647         break;
03648       }
03649       case SolarizePreview:
03650       {
03651         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03652         if (preview_image == (Image *) NULL)
03653           break;
03654         (void) SolarizeImage(preview_image,(double) QuantumRange*
03655           percentage/100.0);
03656         (void) FormatMagickString(label,MaxTextExtent,"solarize %g",
03657           (QuantumRange*percentage)/100.0);
03658         break;
03659       }
03660       case ShadePreview:
03661       {
03662         degrees+=10.0;
03663         preview_image=ShadeImage(thumbnail,MagickTrue,degrees,degrees,
03664           exception);
03665         (void) FormatMagickString(label,MaxTextExtent,"shade %gx%g",
03666           degrees,degrees);
03667         break;
03668       }
03669       case RaisePreview:
03670       {
03671         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03672         if (preview_image == (Image *) NULL)
03673           break;
03674         geometry.width=(unsigned long) (2*i+2);
03675         geometry.height=(unsigned long) (2*i+2);
03676         geometry.x=i/2;
03677         geometry.y=i/2;
03678         (void) RaiseImage(preview_image,&geometry,MagickTrue);
03679         (void) FormatMagickString(label,MaxTextExtent,"raise %lux%lu%+ld%+ld",
03680           geometry.width,geometry.height,geometry.x,geometry.y);
03681         break;
03682       }
03683       case SegmentPreview:
03684       {
03685         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03686         if (preview_image == (Image *) NULL)
03687           break;
03688         threshold+=0.4f;
03689         (void) SegmentImage(preview_image,RGBColorspace,MagickFalse,threshold,
03690           threshold);
03691         (void) FormatMagickString(label,MaxTextExtent,"segment %gx%g",
03692           threshold,threshold);
03693         break;
03694       }
03695       case SwirlPreview:
03696       {
03697         preview_image=SwirlImage(thumbnail,degrees,exception);
03698         (void) FormatMagickString(label,MaxTextExtent,"swirl %g",degrees);
03699         degrees+=45.0;
03700         break;
03701       }
03702       case ImplodePreview:
03703       {
03704         degrees+=0.1f;
03705         preview_image=ImplodeImage(thumbnail,degrees,exception);
03706         (void) FormatMagickString(label,MaxTextExtent,"implode %g",degrees);
03707         break;
03708       }
03709       case WavePreview:
03710       {
03711         degrees+=5.0f;
03712         preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,exception);
03713         (void) FormatMagickString(label,MaxTextExtent,"wave %gx%g",
03714           0.5*degrees,2.0*degrees);
03715         break;
03716       }
03717       case OilPaintPreview:
03718       {
03719         preview_image=OilPaintImage(thumbnail,(double) radius,exception);
03720         (void) FormatMagickString(label,MaxTextExtent,"paint %g",radius);
03721         break;
03722       }
03723       case CharcoalDrawingPreview:
03724       {
03725         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
03726           exception);
03727         (void) FormatMagickString(label,MaxTextExtent,"charcoal %gx%g",
03728           radius,sigma);
03729         break;
03730       }
03731       case JPEGPreview:
03732       {
03733         char
03734           filename[MaxTextExtent];
03735 
03736         int
03737           file;
03738 
03739         MagickBooleanType
03740           status;
03741 
03742         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
03743         if (preview_image == (Image *) NULL)
03744           break;
03745         preview_info->quality=(unsigned long) percentage;
03746         (void) FormatMagickString(factor,MaxTextExtent,"%lu",
03747           preview_info->quality);
03748         file=AcquireUniqueFileResource(filename);
03749         if (file != -1)
03750           file=close(file)-1;
03751         (void) FormatMagickString(preview_image->filename,MaxTextExtent,
03752           "jpeg:%s",filename);
03753         status=WriteImage(preview_info,preview_image);
03754         if (status != MagickFalse)
03755           {
03756             Image
03757               *quality_image;
03758 
03759             (void) CopyMagickString(preview_info->filename,
03760               preview_image->filename,MaxTextExtent);
03761             quality_image=ReadImage(preview_info,exception);
03762             if (quality_image != (Image *) NULL)
03763               {
03764                 preview_image=DestroyImage(preview_image);
03765                 preview_image=quality_image;
03766               }
03767           }
03768         (void) RelinquishUniqueFileResource(preview_image->filename);
03769         if ((GetBlobSize(preview_image)/1024) >= 1024)
03770           (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%gmb ",
03771             factor,(double) ((MagickOffsetType) GetBlobSize(preview_image))/
03772             1024.0/1024.0);
03773         else
03774           if (GetBlobSize(preview_image) >= 1024)
03775             (void) FormatMagickString(label,MaxTextExtent,
03776               "quality %s\n%gkb ",factor,(double) ((MagickOffsetType)
03777               GetBlobSize(preview_image))/1024.0);
03778           else
03779             (void) FormatMagickString(label,MaxTextExtent,"quality %s\n%lub ",
03780               factor,(unsigned long) GetBlobSize(thumbnail));
03781         break;
03782       }
03783     }
03784     thumbnail=DestroyImage(thumbnail);
03785     percentage+=12.5;
03786     radius+=0.5;
03787     sigma+=0.25;
03788     if (preview_image == (Image *) NULL)
03789       break;
03790     (void) DeleteImageProperty(preview_image,"label");
03791     (void) SetImageProperty(preview_image,"label",label);
03792     AppendImageToList(&images,preview_image);
03793     proceed=SetImageProgress(image,PreviewImageTag,i,NumberTiles);
03794     if (proceed == MagickFalse)
03795       break;
03796   }
03797   if (images == (Image *) NULL)
03798     {
03799       preview_info=DestroyImageInfo(preview_info);
03800       return((Image *) NULL);
03801     }
03802   /*
03803     Create the montage.
03804   */
03805   montage_info=CloneMontageInfo(preview_info,(MontageInfo *) NULL);
03806   (void) CopyMagickString(montage_info->filename,image->filename,MaxTextExtent);
03807   montage_info->shadow=MagickTrue;
03808   (void) CloneString(&montage_info->tile,"3x3");
03809   (void) CloneString(&montage_info->geometry,DefaultPreviewGeometry);
03810   (void) CloneString(&montage_info->frame,DefaultTileFrame);
03811   montage_image=MontageImages(images,montage_info,exception);
03812   montage_info=DestroyMontageInfo(montage_info);
03813   images=DestroyImageList(images);
03814   if (montage_image == (Image *) NULL)
03815     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
03816   if (montage_image->montage != (char *) NULL)
03817     {
03818       /*
03819         Free image directory.
03820       */
03821       montage_image->montage=(char *) RelinquishMagickMemory(
03822         montage_image->montage);
03823       if (image->directory != (char *) NULL)
03824         montage_image->directory=(char *) RelinquishMagickMemory(
03825           montage_image->directory);
03826     }
03827   preview_info=DestroyImageInfo(preview_info);
03828   return(montage_image);
03829 }
03830 
03831 /*
03832 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03833 %                                                                             %
03834 %                                                                             %
03835 %                                                                             %
03836 %     R a d i a l B l u r I m a g e                                           %
03837 %                                                                             %
03838 %                                                                             %
03839 %                                                                             %
03840 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
03841 %
03842 %  RadialBlurImage() applies a radial blur to the image.
03843 %
03844 %  Andrew Protano contributed this effect.
03845 %
03846 %  The format of the RadialBlurImage method is:
03847 %
03848 %    Image *RadialBlurImage(const Image *image,const double angle,
03849 %      ExceptionInfo *exception)
03850 %    Image *RadialBlurImageChannel(const Image *image,const ChannelType channel,
03851 %      const double angle,ExceptionInfo *exception)
03852 %
03853 %  A description of each parameter follows:
03854 %
03855 %    o image: the image.
03856 %
03857 %    o channel: the channel type.
03858 %
03859 %    o angle: the angle of the radial blur.
03860 %
03861 %    o exception: return any errors or warnings in this structure.
03862 %
03863 */
03864 
03865 MagickExport Image *RadialBlurImage(const Image *image,const double angle,
03866   ExceptionInfo *exception)
03867 {
03868   Image
03869     *blur_image;
03870 
03871   blur_image=RadialBlurImageChannel(image,DefaultChannels,angle,exception);
03872   return(blur_image);
03873 }
03874 
03875 MagickExport Image *RadialBlurImageChannel(const Image *image,
03876   const ChannelType channel,const double angle,ExceptionInfo *exception)
03877 {
03878   CacheView
03879     *blur_view,
03880     *image_view;
03881 
03882   Image
03883     *blur_image;
03884 
03885   long
03886     progress,
03887     y;
03888 
03889   MagickBooleanType
03890     status;
03891 
03892   MagickPixelPacket
03893     bias;
03894 
03895   MagickRealType
03896     blur_radius,
03897     *cos_theta,
03898     offset,
03899     *sin_theta,
03900     theta;
03901 
03902   PointInfo
03903     blur_center;
03904 
03905   register long
03906     i;
03907 
03908   unsigned long
03909     n;
03910 
03911   /*
03912     Allocate blur image.
03913   */
03914   assert(image != (Image *) NULL);
03915   assert(image->signature == MagickSignature);
03916   if (image->debug != MagickFalse)
03917     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
03918   assert(exception != (ExceptionInfo *) NULL);
03919   assert(exception->signature == MagickSignature);
03920   blur_image=CloneImage(image,0,0,MagickTrue,exception);
03921   if (blur_image == (Image *) NULL)
03922     return((Image *) NULL);
03923   if (SetImageStorageClass(blur_image,DirectClass) == MagickFalse)
03924     {
03925       InheritException(exception,&blur_image->exception);
03926       blur_image=DestroyImage(blur_image);
03927       return((Image *) NULL);
03928     }
03929   blur_center.x=(double) image->columns/2.0;
03930   blur_center.y=(double) image->rows/2.0;
03931   blur_radius=hypot(blur_center.x,blur_center.y);
03932   n=(unsigned long) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+
03933     2UL);
03934   theta=DegreesToRadians(angle)/(MagickRealType) (n-1);
03935   cos_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
03936     sizeof(*cos_theta));
03937   sin_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
03938     sizeof(*sin_theta));
03939   if ((cos_theta == (MagickRealType *) NULL) ||
03940       (sin_theta == (MagickRealType *) NULL))
03941     {
03942       blur_image=DestroyImage(blur_image);
03943       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
03944     }
03945   offset=theta*(MagickRealType) (n-1)/2.0;
03946   for (i=0; i < (long) n; i++)
03947   {
03948     cos_theta[i]=cos((double) (theta*i-offset));
03949     sin_theta[i]=sin((double) (theta*i-offset));
03950   }
03951   /*
03952     Radial blur image.
03953   */
03954   status=MagickTrue;
03955   progress=0;
03956   GetMagickPixelPacket(image,&bias);
03957   image_view=AcquireCacheView(image);
03958   blur_view=AcquireCacheView(blur_image);
03959 #if defined(MAGICKCORE_OPENMP_SUPPORT)
03960   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
03961 #endif
03962   for (y=0; y < (long) blur_image->rows; y++)
03963   {
03964     register const IndexPacket
03965       *restrict indexes;
03966 
03967     register IndexPacket
03968       *restrict blur_indexes;
03969 
03970     register long
03971       x;
03972 
03973     register PixelPacket
03974       *restrict q;
03975 
03976     if (status == MagickFalse)
03977       continue;
03978     q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
03979       exception);
03980     if (q == (PixelPacket *) NULL)
03981       {
03982         status=MagickFalse;
03983         continue;
03984       }
03985     blur_indexes=GetCacheViewAuthenticIndexQueue(blur_view);
03986     for (x=0; x < (long) blur_image->columns; x++)
03987     {
03988       MagickPixelPacket
03989         qixel;
03990 
03991       MagickRealType
03992         normalize,
03993         radius;
03994 
03995       PixelPacket
03996         pixel;
03997 
03998       PointInfo
03999         center;
04000 
04001       register long
04002         i;
04003 
04004       unsigned long
04005         step;
04006 
04007       center.x=(double) x-blur_center.x;
04008       center.y=(double) y-blur_center.y;
04009       radius=hypot((double) center.x,center.y);
04010       if (radius == 0)
04011         step=1;
04012       else
04013         {
04014           step=(unsigned long) (blur_radius/radius);
04015           if (step == 0)
04016             step=1;
04017           else
04018             if (step >= n)
04019               step=n-1;
04020         }
04021       normalize=0.0;
04022       qixel=bias;
04023       if (((channel & OpacityChannel) == 0) || (image->matte == MagickFalse))
04024         {
04025           for (i=0; i < (long) n; i+=step)
04026           {
04027             (void) GetOneCacheViewVirtualPixel(image_view,(long) (blur_center.x+
04028               center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),(long) (
04029               blur_center.y+center.x*sin_theta[i]+center.y*cos_theta[i]+0.5),
04030               &pixel,exception);
04031             qixel.red+=pixel.red;
04032             qixel.green+=pixel.green;
04033             qixel.blue+=pixel.blue;
04034             qixel.opacity+=pixel.opacity;
04035             if (image->colorspace == CMYKColorspace)
04036               {
04037                 indexes=GetCacheViewVirtualIndexQueue(image_view);
04038                 qixel.index+=(*indexes);
04039               }
04040             normalize+=1.0;
04041           }
04042           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
04043             normalize);
04044           if ((channel & RedChannel) != 0)
04045             q->red=ClampToQuantum(normalize*qixel.red);
04046           if ((channel & GreenChannel) != 0)
04047             q->green=ClampToQuantum(normalize*qixel.green);
04048           if ((channel & BlueChannel) != 0)
04049             q->blue=ClampToQuantum(normalize*qixel.blue);
04050           if ((channel & OpacityChannel) != 0)
04051             q->opacity=ClampToQuantum(normalize*qixel.opacity);
04052           if (((channel & IndexChannel) != 0) &&
04053               (image->colorspace == CMYKColorspace))
04054             blur_indexes[x]=(IndexPacket) ClampToQuantum(normalize*qixel.index);
04055         }
04056       else
04057         {
04058           MagickRealType
04059             alpha,
04060             gamma;
04061 
04062           alpha=1.0;
04063           gamma=0.0;
04064           for (i=0; i < (long) n; i+=step)
04065           {
04066             (void) GetOneCacheViewVirtualPixel(image_view,(long) (blur_center.x+
04067               center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),(long) (
04068               blur_center.y+center.x*sin_theta[i]+center.y*cos_theta[i]+0.5),
04069               &pixel,exception);
04070             alpha=(MagickRealType) (QuantumScale*
04071               GetAlphaPixelComponent(&pixel));
04072             qixel.red+=alpha*pixel.red;
04073             qixel.green+=alpha*pixel.green;
04074             qixel.blue+=alpha*pixel.blue;
04075             qixel.opacity+=pixel.opacity;
04076             if (image->colorspace == CMYKColorspace)
04077               {
04078                 indexes=GetCacheViewVirtualIndexQueue(image_view);
04079                 qixel.index+=alpha*(*indexes);
04080               }
04081             gamma+=alpha;
04082             normalize+=1.0;
04083           }
04084           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
04085           normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
04086             normalize);
04087           if ((channel & RedChannel) != 0)
04088             q->red=ClampToQuantum(gamma*qixel.red);
04089           if ((channel & GreenChannel) != 0)
04090             q->green=ClampToQuantum(gamma*qixel.green);
04091           if ((channel & BlueChannel) != 0)
04092             q->blue=ClampToQuantum(gamma*qixel.blue);
04093           if ((channel & OpacityChannel) != 0)
04094             q->opacity=ClampToQuantum(normalize*qixel.opacity);
04095           if (((channel & IndexChannel) != 0) &&
04096               (image->colorspace == CMYKColorspace))
04097             blur_indexes[x]=(IndexPacket) ClampToQuantum(gamma*qixel.index);
04098         }
04099       q++;
04100     }
04101     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
04102       status=MagickFalse;
04103     if (image->progress_monitor != (MagickProgressMonitor) NULL)
04104       {
04105         MagickBooleanType
04106           proceed;
04107 
04108 #if defined(MAGICKCORE_OPENMP_SUPPORT)
04109   #pragma omp critical (MagickCore_RadialBlurImageChannel)
04110 #endif
04111         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
04112         if (proceed == MagickFalse)
04113           status=MagickFalse;
04114       }
04115   }
04116   blur_view=DestroyCacheView(blur_view);
04117   image_view=DestroyCacheView(image_view);
04118   cos_theta=(MagickRealType *) RelinquishMagickMemory(cos_theta);
04119   sin_theta=(MagickRealType *) RelinquishMagickMemory(sin_theta);
04120   if (status == MagickFalse)
04121     blur_image=DestroyImage(blur_image);
04122   return(blur_image);
04123 }
04124 
04125 /*
04126 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
04127 %                                                                             %
04128 %                                                                             %
04129 %                                                                             %
04130 %     R e d u c e N o i s e I m a g e                                         %
04131 %                                                                             %
04132 %                                                                             %
04133 %                                                                             %
04134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
04135 %
04136 %  ReduceNoiseImage() smooths the contours of an image while still preserving
04137 %  edge information.  The algorithm works by replacing each pixel with its
04138 %  neighbor closest in value.  A neighbor is defined by radius.  Use a radius
04139 %  of 0 and ReduceNoise() selects a suitable radius for you.
04140 %
04141 %  The format of the ReduceNoiseImage method is:
04142 %
04143 %      Image *ReduceNoiseImage(const Image *image,const double radius,
04144 %        ExceptionInfo *exception)
04145 %
04146 %  A description of each parameter follows:
04147 %
04148 %    o image: the image.
04149 %
04150 %    o radius: the radius of the pixel neighborhood.
04151 %
04152 %    o exception: return any errors or warnings in this structure.
04153 %
04154 */
04155 
04156 static MagickPixelPacket GetNonpeakMedianPixelList(MedianPixelList *pixel_list)
04157 {
04158   MagickPixelPacket
04159     pixel;
04160 
04161   register long
04162     channel;
04163 
04164   register MedianSkipList
04165     *list;
04166 
04167   unsigned long
04168     center,
04169     color,
04170     count,
04171     previous,
04172     next;
04173 
04174   unsigned short
04175     channels[5];
04176 
04177   /*
04178     Finds the median value for each of the color.
04179   */
04180   center=pixel_list->center;
04181   for (channel=0; channel < 5; channel++)
04182   {
04183     list=pixel_list->lists+channel;
04184     color=65536UL;
04185     next=list->nodes[color].next[0];
04186     count=0;
04187     do
04188     {
04189       previous=color;
04190       color=next;
04191       next=list->nodes[color].next[0];
04192       count+=list->nodes[color].count;
04193     }
04194     while (count <= center);
04195     if ((previous == 65536UL) && (next != 65536UL))
04196       color=next;
04197     else
04198       if ((previous != 65536UL) && (next == 65536UL))
04199         color=previous;
04200     channels[channel]=(unsigned short) color;
04201   }
04202   GetMagickPixelPacket((const Image *) NULL,&pixel);
04203   pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
04204   pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
04205   pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
04206   pixel.opacity=(MagickRealType) ScaleShortToQuantum(channels[3]);
04207   pixel.index=(MagickRealType) ScaleShortToQuantum(channels[4]);
04208   return(pixel);
04209 }
04210 
04211 MagickExport Image *ReduceNoiseImage(const Image *image,const double radius,
04212   ExceptionInfo *exception)
04213 {
04214 #define ReduceNoiseImageTag  "ReduceNoise/Image"
04215 
04216   CacheView
04217     *image_view,
04218     *noise_view;
04219 
04220   Image
04221     *noise_image;
04222 
04223   long
04224     progress,
04225     y;
04226 
04227   MagickBooleanType
04228     status;
04229 
04230   MedianPixelList
04231     **restrict pixel_list;
04232 
04233   unsigned long
04234     width;
04235 
04236   /*
04237     Initialize noise image attributes.
04238   */
04239   assert(image != (Image *) NULL);
04240   assert(image->signature == MagickSignature);
04241   if (image->debug != MagickFalse)
04242     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
04243   assert(exception != (ExceptionInfo *) NULL);
04244   assert(exception->signature == MagickSignature);
04245   width=GetOptimalKernelWidth2D(radius,0.5);
04246   if ((image->columns < width) || (image->rows < width))
04247     ThrowImageException(OptionError,"ImageSmallerThanKernelRadius");
04248   noise_image=CloneImage(image,image->columns,image->rows,MagickTrue,
04249     exception);
04250   if (noise_image == (Image *) NULL)
04251     return((Image *) NULL);
04252   if (SetImageStorageClass(noise_image,DirectClass) == MagickFalse)
04253     {
04254       InheritException(exception,&noise_image->exception);
04255       noise_image=DestroyImage(noise_image);
04256       return((Image *) NULL);
04257     }
04258   pixel_list=AcquireMedianPixelListThreadSet(width);
04259   if (pixel_list == (MedianPixelList **) NULL)
04260     {
04261       noise_image=DestroyImage(noise_image);
04262       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
04263     }
04264   /*
04265     Reduce noise image.
04266   */
04267   status=MagickTrue;
04268   progress=0;
04269   image_view=AcquireCacheView(image);
04270   noise_view=AcquireCacheView(noise_image);
04271 #if defined(MAGICKCORE_OPENMP_SUPPORT)
04272   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
04273 #endif
04274   for (y=0; y < (long) noise_image->rows; y++)
04275   {
04276     register const IndexPacket
04277       *restrict indexes;
04278 
04279     register const PixelPacket
04280       *restrict p;
04281 
04282     register IndexPacket
04283       *restrict noise_indexes;
04284 
04285     register long
04286       id,
04287       x;
04288 
04289     register PixelPacket
04290       *restrict q;
04291 
04292     if (status == MagickFalse)
04293       continue;
04294     p=GetCacheViewVirtualPixels(image_view,-((long) width/2L),y-(long) (width/
04295       2L),image->columns+width,width,exception);
04296     q=QueueCacheViewAuthenticPixels(noise_view,0,y,noise_image->columns,1,
04297       exception);
04298     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
04299       {
04300         status=MagickFalse;
04301         continue;
04302       }
04303     indexes=GetCacheViewVirtualIndexQueue(image_view);
04304     noise_indexes=GetCacheViewAuthenticIndexQueue(noise_view);
04305     id=GetOpenMPThreadId();
04306     for (x=0; x < (long) noise_image->columns; x++)
04307     {
04308       MagickPixelPacket
04309         pixel;
04310 
04311       register const PixelPacket
04312         *restrict r;
04313 
04314       register const IndexPacket
04315         *restrict s;
04316 
04317       register long
04318         u,
04319         v;
04320 
04321       r=p;
04322       s=indexes+x;
04323       ResetMedianPixelList(pixel_list[id]);
04324       for (v=0; v < (long) width; v++)
04325       {
04326         for (u=0; u < (long) width; u++)
04327           InsertMedianPixelList(image,r+u,s+u,pixel_list[id]);
04328         r+=image->columns+width;
04329         s+=image->columns+width;
04330       }
04331       pixel=GetNonpeakMedianPixelList(pixel_list[id]);
04332       SetPixelPacket(noise_image,&pixel,q,noise_indexes+x);
04333       p++;
04334       q++;
04335     }
04336     if (SyncCacheViewAuthenticPixels(noise_view,exception) == MagickFalse)
04337       status=MagickFalse;
04338     if (image->progress_monitor != (MagickProgressMonitor) NULL)
04339       {
04340         MagickBooleanType
04341           proceed;
04342 
04343 #if defined(MAGICKCORE_OPENMP_SUPPORT)
04344   #pragma omp critical (MagickCore_ReduceNoiseImage)
04345 #endif
04346         proceed=SetImageProgress(image,ReduceNoiseImageTag,progress++,
04347