]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/effect.c
(no commit message)
[imagemagick] / MagickCore / effect.c
index 3fdfab1a2c72b1a2220c3e874dbd7bf9d6266717..ea9782d2d7a6525b0735848e81bb1246a6aed043 100644 (file)
@@ -17,7 +17,7 @@
 %                                 October 1996                                %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
 %  dedicated to making software imaging solutions freely available.           %
 %                                                                             %
 %  You may not use this file except in compliance with the License.  You may  %
@@ -49,6 +49,7 @@
 #include "MagickCore/colorspace.h"
 #include "MagickCore/constitute.h"
 #include "MagickCore/decorate.h"
+#include "MagickCore/distort.h"
 #include "MagickCore/draw.h"
 #include "MagickCore/enhance.h"
 #include "MagickCore/exception.h"
 #include "MagickCore/list.h"
 #include "MagickCore/log.h"
 #include "MagickCore/memory_.h"
+#include "MagickCore/memory-private.h"
 #include "MagickCore/monitor.h"
 #include "MagickCore/monitor-private.h"
 #include "MagickCore/montage.h"
 #include "MagickCore/morphology.h"
 #include "MagickCore/paint.h"
 #include "MagickCore/pixel-accessor.h"
+#include "MagickCore/pixel-private.h"
 #include "MagickCore/property.h"
 #include "MagickCore/quantize.h"
 #include "MagickCore/quantum.h"
@@ -81,6 +84,7 @@
 #include "MagickCore/segment.h"
 #include "MagickCore/shear.h"
 #include "MagickCore/signature-private.h"
+#include "MagickCore/statistic.h"
 #include "MagickCore/string_.h"
 #include "MagickCore/thread-private.h"
 #include "MagickCore/transform.h"
 %  The format of the AdaptiveBlurImage method is:
 %
 %      Image *AdaptiveBlurImage(const Image *image,const double radius,
-%        const double sigma,const double bias,ExceptionInfo *exception)
+%        const double sigma,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %
 %    o sigma: the standard deviation of the Laplacian, in pixels.
 %
-%    o bias: the bias.
-%
 %    o exception: return any errors or warnings in this structure.
 %
 */
@@ -167,12 +169,11 @@ MagickExport MagickBooleanType AdaptiveLevelImage(Image *image,
   return(status);
 }
 
-MagickExport Image *AdaptiveBlurImage(const Image *image,
-  const double radius,const double sigma,const double bias,
-  ExceptionInfo *exception)
+MagickExport Image *AdaptiveBlurImage(const Image *image,const double radius,
+  const double sigma,ExceptionInfo *exception)
 {
 #define AdaptiveBlurImageTag  "Convolve/Image"
-#define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
+#define MagickSigma  (fabs(sigma) < MagickEpsilon ? MagickEpsilon : sigma)
 
   CacheView
     *blur_view,
@@ -180,7 +181,6 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
     *image_view;
 
   double
-    **kernel,
     normalize;
 
   Image
@@ -194,6 +194,9 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
   MagickOffsetType
     progress;
 
+  MagickRealType
+    **kernel;
+
   register ssize_t
     i;
 
@@ -216,7 +219,7 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
   blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
   if (blur_image == (Image *) NULL)
     return((Image *) NULL);
-  if (fabs(sigma) <= MagickEpsilon)
+  if (fabs(sigma) < MagickEpsilon)
     return(blur_image);
   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
     {
@@ -233,7 +236,7 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
       return((Image *) NULL);
     }
   (void) AdaptiveLevelImage(edge_image,"20%,95%",exception);
-  gaussian_image=GaussianBlurImage(edge_image,radius,sigma,bias,exception);
+  gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
   if (gaussian_image != (Image *) NULL)
     {
       edge_image=DestroyImage(edge_image);
@@ -244,8 +247,9 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
     Create a set of kernels from maximum (radius,sigma) to minimum.
   */
   width=GetOptimalKernelWidth2D(radius,sigma);
-  kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
-  if (kernel == (double **) NULL)
+  kernel=(MagickRealType **) MagickAssumeAligned(AcquireAlignedMemory((size_t)
+    width,sizeof(*kernel)));
+  if (kernel == (MagickRealType  **) NULL)
     {
       edge_image=DestroyImage(edge_image);
       blur_image=DestroyImage(blur_image);
@@ -254,9 +258,9 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
   for (i=0; i < (ssize_t) width; i+=2)
   {
-    kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
-      sizeof(**kernel));
-    if (kernel[i] == (double *) NULL)
+    kernel[i]=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
+      (size_t) (width-i),(width-i)*sizeof(**kernel)));
+    if (kernel[i] == (MagickRealType *) NULL)
       break;
     normalize=0.0;
     j=(ssize_t) (width-i)/2;
@@ -265,23 +269,23 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
     {
       for (u=(-j); u <= j; u++)
       {
-        kernel[i][k]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
+        kernel[i][k]=(MagickRealType) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
           MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
         normalize+=kernel[i][k];
         k++;
       }
     }
-    if (fabs(normalize) <= MagickEpsilon)
-      normalize=1.0;
-    normalize=1.0/normalize;
+    if (fabs(normalize) < MagickEpsilon)
+      normalize=MagickEpsilon;
+    normalize=MagickEpsilonReciprocal(normalize);
     for (k=0; k < (j*j); k++)
       kernel[i][k]=normalize*kernel[i][k];
   }
   if (i < (ssize_t) width)
     {
       for (i-=2; i >= 0; i-=2)
-        kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
-      kernel=(double **) RelinquishMagickMemory(kernel);
+        kernel[i]=(MagickRealType *) RelinquishAlignedMemory(kernel[i]);
+      kernel=(MagickRealType **) RelinquishAlignedMemory(kernel);
       edge_image=DestroyImage(edge_image);
       blur_image=DestroyImage(blur_image);
       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
@@ -291,11 +295,12 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
-  edge_view=AcquireCacheView(edge_image);
-  blur_view=AcquireCacheView(blur_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  edge_view=AcquireVirtualCacheView(edge_image,exception);
+  blur_view=AcquireAuthenticCacheView(blur_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) blur_image->rows; y++)
   {
@@ -343,11 +348,11 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
         (ssize_t) ((width-j)/2L),width-j,width-j,exception);
       if (p == (const Quantum *) NULL)
         break;
-      center=(ssize_t) GetPixelChannels(image)*(width-j)*
-        ((width-j)/2L)+GetPixelChannels(image)*((width-j)/2);
+      center=(ssize_t) GetPixelChannels(image)*(width-j)*((width-j)/2L)+
+        GetPixelChannels(image)*((width-j)/2L);
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           alpha,
           gamma,
           pixel;
@@ -359,7 +364,7 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
           blur_traits,
           traits;
 
-        register const double
+        register const MagickRealType
           *restrict k;
 
         register const Quantum
@@ -371,20 +376,21 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
         ssize_t
           v;
 
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        blur_traits=GetPixelChannelMapTraits(blur_image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        blur_traits=GetPixelChannelTraits(blur_image,channel);
         if ((traits == UndefinedPixelTrait) ||
             (blur_traits == UndefinedPixelTrait))
           continue;
-        if ((blur_traits & CopyPixelTrait) != 0)
+        if (((blur_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
-            q[channel]=p[center+i];
+            SetPixelChannel(blur_image,channel,p[center+i],q);
             continue;
           }
         k=kernel[j];
         pixels=p;
-        pixel=bias;
+        pixel=0.0;
         gamma=0.0;
         if ((blur_traits & BlendPixelTrait) == 0)
           {
@@ -401,8 +407,8 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
                 pixels+=GetPixelChannels(image);
               }
             }
-            gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-            q[channel]=ClampToQuantum(gamma*pixel);
+            gamma=MagickEpsilonReciprocal(gamma);
+            SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
             continue;
           }
         /*
@@ -412,15 +418,15 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
         {
           for (u=0; u < (ssize_t) (width-j); u++)
           {
-            alpha=(MagickRealType) (QuantumScale*GetPixelAlpha(image,pixels));
+            alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
             pixel+=(*k)*alpha*pixels[i];
             gamma+=(*k)*alpha;
             k++;
             pixels+=GetPixelChannels(image);
           }
         }
-        gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-        q[channel]=ClampToQuantum(gamma*pixel);
+        gamma=MagickEpsilonReciprocal(gamma);
+        SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
       }
       q+=GetPixelChannels(blur_image);
       r+=GetPixelChannels(edge_image);
@@ -433,7 +439,7 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_AdaptiveBlurImage)
+        #pragma omp critical (MagickCore_AdaptiveBlurImage)
 #endif
         proceed=SetImageProgress(image,AdaptiveBlurImageTag,progress++,
           image->rows);
@@ -447,8 +453,8 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
   image_view=DestroyCacheView(image_view);
   edge_image=DestroyImage(edge_image);
   for (i=0; i < (ssize_t) width;  i+=2)
-    kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
-  kernel=(double **) RelinquishMagickMemory(kernel);
+    kernel[i]=(MagickRealType *) RelinquishAlignedMemory(kernel[i]);
+  kernel=(MagickRealType **) RelinquishAlignedMemory(kernel);
   if (status == MagickFalse)
     blur_image=DestroyImage(blur_image);
   return(blur_image);
@@ -474,7 +480,7 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
 %  The format of the AdaptiveSharpenImage method is:
 %
 %      Image *AdaptiveSharpenImage(const Image *image,const double radius,
-%        const double sigma,const double bias,ExceptionInfo *exception)
+%        const double sigma,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
@@ -485,16 +491,14 @@ MagickExport Image *AdaptiveBlurImage(const Image *image,
 %
 %    o sigma: the standard deviation of the Laplacian, in pixels.
 %
-%    o bias: the bias.
-%
 %    o exception: return any errors or warnings in this structure.
 %
 */
 MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
-  const double sigma,const double bias,ExceptionInfo *exception)
+  const double sigma,ExceptionInfo *exception)
 {
 #define AdaptiveSharpenImageTag  "Convolve/Image"
-#define MagickSigma  (fabs(sigma) <= MagickEpsilon ? 1.0 : sigma)
+#define MagickSigma  (fabs(sigma) < MagickEpsilon ? MagickEpsilon : sigma)
 
   CacheView
     *sharp_view,
@@ -502,7 +506,6 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
     *image_view;
 
   double
-    **kernel,
     normalize;
 
   Image
@@ -516,6 +519,9 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
   MagickOffsetType
     progress;
 
+  MagickRealType
+    **kernel;
+
   register ssize_t
     i;
 
@@ -538,7 +544,7 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
   sharp_image=CloneImage(image,0,0,MagickTrue,exception);
   if (sharp_image == (Image *) NULL)
     return((Image *) NULL);
-  if (fabs(sigma) <= MagickEpsilon)
+  if (fabs(sigma) < MagickEpsilon)
     return(sharp_image);
   if (SetImageStorageClass(sharp_image,DirectClass,exception) == MagickFalse)
     {
@@ -555,7 +561,7 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
       return((Image *) NULL);
     }
   (void) AdaptiveLevelImage(edge_image,"20%,95%",exception);
-  gaussian_image=GaussianBlurImage(edge_image,radius,sigma,bias,exception);
+  gaussian_image=GaussianBlurImage(edge_image,radius,sigma,exception);
   if (gaussian_image != (Image *) NULL)
     {
       edge_image=DestroyImage(edge_image);
@@ -566,8 +572,9 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
     Create a set of kernels from maximum (radius,sigma) to minimum.
   */
   width=GetOptimalKernelWidth2D(radius,sigma);
-  kernel=(double **) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
-  if (kernel == (double **) NULL)
+  kernel=(MagickRealType **) MagickAssumeAligned(AcquireAlignedMemory((size_t)
+    width,sizeof(*kernel)));
+  if (kernel == (MagickRealType **) NULL)
     {
       edge_image=DestroyImage(edge_image);
       sharp_image=DestroyImage(sharp_image);
@@ -576,9 +583,9 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
   (void) ResetMagickMemory(kernel,0,(size_t) width*sizeof(*kernel));
   for (i=0; i < (ssize_t) width; i+=2)
   {
-    kernel[i]=(double *) AcquireQuantumMemory((size_t) (width-i),(width-i)*
-      sizeof(**kernel));
-    if (kernel[i] == (double *) NULL)
+    kernel[i]=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory(
+      (size_t) (width-i),(width-i)*sizeof(**kernel)));
+    if (kernel[i] == (MagickRealType *) NULL)
       break;
     normalize=0.0;
     j=(ssize_t) (width-i)/2;
@@ -587,23 +594,23 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
     {
       for (u=(-j); u <= j; u++)
       {
-        kernel[i][k]=(double) (-exp(-((double) u*u+v*v)/(2.0*MagickSigma*
-          MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
+        kernel[i][k]=(MagickRealType) (-exp(-((double) u*u+v*v)/(2.0*
+          MagickSigma*MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
         normalize+=kernel[i][k];
         k++;
       }
     }
-    if (fabs(normalize) <= MagickEpsilon)
-      normalize=1.0;
-    normalize=1.0/normalize;
+    if (fabs(normalize) < MagickEpsilon)
+      normalize=MagickEpsilon;
+    normalize=MagickEpsilonReciprocal(normalize);
     for (k=0; k < (j*j); k++)
       kernel[i][k]=normalize*kernel[i][k];
   }
   if (i < (ssize_t) width)
     {
       for (i-=2; i >= 0; i-=2)
-        kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
-      kernel=(double **) RelinquishMagickMemory(kernel);
+        kernel[i]=(MagickRealType *) RelinquishAlignedMemory(kernel[i]);
+      kernel=(MagickRealType **) RelinquishAlignedMemory(kernel);
       edge_image=DestroyImage(edge_image);
       sharp_image=DestroyImage(sharp_image);
       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
@@ -613,11 +620,12 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
-  edge_view=AcquireCacheView(edge_image);
-  sharp_view=AcquireCacheView(sharp_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  edge_view=AcquireVirtualCacheView(edge_image,exception);
+  sharp_view=AcquireAuthenticCacheView(sharp_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) sharp_image->rows; y++)
   {
@@ -665,11 +673,11 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
         (ssize_t) ((width-j)/2L),width-j,width-j,exception);
       if (p == (const Quantum *) NULL)
         break;
-      center=(ssize_t) GetPixelChannels(image)*(width-j)*
-        ((width-j)/2L)+GetPixelChannels(image)*((width-j)/2);
-      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      center=(ssize_t) GetPixelChannels(image)*(width-j)*((width-j)/2L)+
+        GetPixelChannels(image)*((width-j)/2);
+      for (i=0; i < (ssize_t) GetPixelChannels(sharp_image); i++)
       {
-        MagickRealType
+        double
           alpha,
           gamma,
           pixel;
@@ -681,7 +689,7 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
           sharp_traits,
           traits;
 
-        register const double
+        register const MagickRealType
           *restrict k;
 
         register const Quantum
@@ -693,20 +701,21 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
         ssize_t
           v;
 
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        sharp_traits=GetPixelChannelMapTraits(sharp_image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        sharp_traits=GetPixelChannelTraits(sharp_image,channel);
         if ((traits == UndefinedPixelTrait) ||
             (sharp_traits == UndefinedPixelTrait))
           continue;
-        if ((sharp_traits & CopyPixelTrait) != 0)
+        if (((sharp_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
-            q[channel]=p[center+i];
+            SetPixelChannel(sharp_image,channel,p[center+i],q);
             continue;
           }
         k=kernel[j];
         pixels=p;
-        pixel=bias;
+        pixel=0.0;
         gamma=0.0;
         if ((sharp_traits & BlendPixelTrait) == 0)
           {
@@ -723,8 +732,8 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
                 pixels+=GetPixelChannels(image);
               }
             }
-            gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-            q[channel]=ClampToQuantum(gamma*pixel);
+            gamma=MagickEpsilonReciprocal(gamma);
+            SetPixelChannel(sharp_image,channel,ClampToQuantum(gamma*pixel),q);
             continue;
           }
         /*
@@ -734,15 +743,15 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
         {
           for (u=0; u < (ssize_t) (width-j); u++)
           {
-            alpha=(MagickRealType) (QuantumScale*GetPixelAlpha(image,pixels));
+            alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
             pixel+=(*k)*alpha*pixels[i];
             gamma+=(*k)*alpha;
             k++;
             pixels+=GetPixelChannels(image);
           }
         }
-        gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-        q[channel]=ClampToQuantum(gamma*pixel);
+        gamma=MagickEpsilonReciprocal(gamma);
+        SetPixelChannel(sharp_image,channel,ClampToQuantum(gamma*pixel),q);
       }
       q+=GetPixelChannels(sharp_image);
       r+=GetPixelChannels(edge_image);
@@ -755,7 +764,7 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_AdaptiveSharpenImage)
+        #pragma omp critical (MagickCore_AdaptiveSharpenImage)
 #endif
         proceed=SetImageProgress(image,AdaptiveSharpenImageTag,progress++,
           image->rows);
@@ -769,8 +778,8 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
   image_view=DestroyCacheView(image_view);
   edge_image=DestroyImage(edge_image);
   for (i=0; i < (ssize_t) width;  i+=2)
-    kernel[i]=(double *) RelinquishMagickMemory(kernel[i]);
-  kernel=(double **) RelinquishMagickMemory(kernel);
+    kernel[i]=(MagickRealType *) RelinquishAlignedMemory(kernel[i]);
+  kernel=(MagickRealType **) RelinquishAlignedMemory(kernel);
   if (status == MagickFalse)
     sharp_image=DestroyImage(sharp_image);
   return(sharp_image);
@@ -799,7 +808,7 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
 %  The format of the BlurImage method is:
 %
 %      Image *BlurImage(const Image *image,const double radius,
-%        const double sigma,const double bias,ExceptionInfo *exception)
+%        const double sigma,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
@@ -810,15 +819,13 @@ MagickExport Image *AdaptiveSharpenImage(const Image *image,const double radius,
 %
 %    o sigma: the standard deviation of the Gaussian, in pixels.
 %
-%    o bias: the bias.
-%
 %    o exception: return any errors or warnings in this structure.
 %
 */
 
-static double *GetBlurKernel(const size_t width,const double sigma)
+static MagickRealType *GetBlurKernel(const size_t width,const double sigma)
 {
-  double
+  MagickRealType
     *kernel,
     normalize;
 
@@ -833,16 +840,17 @@ static double *GetBlurKernel(const size_t width,const double sigma)
     Generate a 1-D convolution kernel.
   */
   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
-  kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
-  if (kernel == (double *) NULL)
+  kernel=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory((size_t)
+    width,sizeof(*kernel)));
+  if (kernel == (MagickRealType *) NULL)
     return(0);
   normalize=0.0;
   j=(ssize_t) width/2;
   i=0;
   for (k=(-j); k <= j; k++)
   {
-    kernel[i]=(double) (exp(-((double) k*k)/(2.0*MagickSigma*MagickSigma))/
-      (MagickSQ2PI*MagickSigma));
+    kernel[i]=(MagickRealType) (exp(-((double) k*k)/(2.0*MagickSigma*
+      MagickSigma))/(MagickSQ2PI*MagickSigma));
     normalize+=kernel[i];
     i++;
   }
@@ -852,7 +860,7 @@ static double *GetBlurKernel(const size_t width,const double sigma)
 }
 
 MagickExport Image *BlurImage(const Image *image,const double radius,
-  const double sigma,const double bias,ExceptionInfo *exception)
+  const double sigma,ExceptionInfo *exception)
 {
 #define BlurImageTag  "Blur/Image"
 
@@ -860,9 +868,6 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
     *blur_view,
     *image_view;
 
-  double
-    *kernel;
-
   Image
     *blur_image;
 
@@ -872,6 +877,9 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
   MagickOffsetType
     progress;
 
+  MagickRealType
+    *kernel;
+
   register ssize_t
     i;
 
@@ -879,6 +887,7 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
     width;
 
   ssize_t
+    center,
     x,
     y;
 
@@ -894,7 +903,7 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
   blur_image=CloneImage(image,0,0,MagickTrue,exception);
   if (blur_image == (Image *) NULL)
     return((Image *) NULL);
-  if (fabs(sigma) <= MagickEpsilon)
+  if (fabs(sigma) < MagickEpsilon)
     return(blur_image);
   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
     {
@@ -903,7 +912,7 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
     }
   width=GetOptimalKernelWidth1D(radius,sigma);
   kernel=GetBlurKernel(width,sigma);
-  if (kernel == (double *) NULL)
+  if (kernel == (MagickRealType *) NULL)
     {
       blur_image=DestroyImage(blur_image);
       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
@@ -914,11 +923,11 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
         format[MaxTextExtent],
         *message;
 
-      register const double
+      register const MagickRealType
         *k;
 
       (void) LogMagickEvent(TransformEvent,GetMagickModule(),
-        "  BlurImage with %.20g kernel:",(double) width);
+        "  blur image with kernel width %.20g:",(double) width);
       message=AcquireString("");
       k=kernel;
       for (i=0; i < (ssize_t) width; i++)
@@ -926,7 +935,7 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
         *message='\0';
         (void) FormatLocaleString(format,MaxTextExtent,"%.20g: ",(double) i);
         (void) ConcatenateString(&message,format);
-        (void) FormatLocaleString(format,MaxTextExtent,"%g ",*k++);
+        (void) FormatLocaleString(format,MaxTextExtent,"%g ",(double) *k++);
         (void) ConcatenateString(&message,format);
         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
       }
@@ -937,12 +946,14 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
-  blur_view=AcquireCacheView(blur_image);
+  center=(ssize_t) GetPixelChannels(image)*(width/2L);
+  image_view=AcquireVirtualCacheView(image,exception);
+  blur_view=AcquireAuthenticCacheView(blur_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
-  for (y=0; y < (ssize_t) blur_image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const Quantum
       *restrict p;
@@ -957,112 +968,85 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
       continue;
     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y,
       image->columns+width,1,exception);
-    q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
+    q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
       exception);
     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
       {
         status=MagickFalse;
         continue;
       }
-    for (x=0; x < (ssize_t) blur_image->columns; x++)
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
-      PixelInfo
-        pixel;
+      register ssize_t
+        i;
 
-      register const double
-        *restrict k;
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        double
+          alpha,
+          gamma,
+          pixel;
 
-      register const Quantum
-        *restrict kernel_pixels;
+        PixelChannel
+          channel;
 
-      register ssize_t
-        i;
+        PixelTrait
+          blur_traits,
+          traits;
 
-      pixel.red=bias;
-      pixel.green=bias;
-      pixel.blue=bias;
-      pixel.black=bias;
-      pixel.alpha=bias;
-      k=kernel;
-      kernel_pixels=p;
-      if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) == 0) ||
-          (image->matte == MagickFalse))
-        {
-          for (i=0; i < (ssize_t) width; i++)
-          {
-            pixel.red+=(*k)*GetPixelRed(image,kernel_pixels);
-            pixel.green+=(*k)*GetPixelGreen(image,kernel_pixels);
-            pixel.blue+=(*k)*GetPixelBlue(image,kernel_pixels);
-            if (image->colorspace == CMYKColorspace)
-              pixel.black+=(*k)*GetPixelBlack(image,kernel_pixels);
-            k++;
-            kernel_pixels+=GetPixelChannels(image);
-          }
-          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelRed(blur_image,ClampToQuantum(pixel.red),q);
-          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelGreen(blur_image,ClampToQuantum(pixel.green),q);
-          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelBlue(blur_image,ClampToQuantum(pixel.blue),q);
-          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-              (blur_image->colorspace == CMYKColorspace))
-            SetPixelBlack(blur_image,ClampToQuantum(pixel.black),q);
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
-            {
-              k=kernel;
-              kernel_pixels=p;
-              for (i=0; i < (ssize_t) width; i++)
-              {
-                pixel.alpha+=(*k)*GetPixelAlpha(image,kernel_pixels);
-                k++;
-                kernel_pixels+=GetPixelChannels(image);
-              }
-              SetPixelAlpha(blur_image,ClampToQuantum(pixel.alpha),q);
-            }
-        }
-      else
-        {
-          MagickRealType
-            alpha,
-            gamma;
+        register const MagickRealType
+          *restrict k;
+
+        register const Quantum
+          *restrict pixels;
+
+        register ssize_t
+          u;
 
-          gamma=0.0;
-          for (i=0; i < (ssize_t) width; i++)
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        blur_traits=GetPixelChannelTraits(blur_image,channel);
+        if ((traits == UndefinedPixelTrait) ||
+            (blur_traits == UndefinedPixelTrait))
+          continue;
+        if (((blur_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
-            alpha=(MagickRealType) (QuantumScale*
-              GetPixelAlpha(image,kernel_pixels));
-            pixel.red+=(*k)*alpha*GetPixelRed(image,kernel_pixels);
-            pixel.green+=(*k)*alpha*GetPixelGreen(image,kernel_pixels);
-            pixel.blue+=(*k)*alpha*GetPixelBlue(image,kernel_pixels);
-            if (image->colorspace == CMYKColorspace)
-              pixel.black+=(*k)*alpha*GetPixelBlack(image,kernel_pixels);
-            gamma+=(*k)*alpha;
-            k++;
-            kernel_pixels+=GetPixelChannels(image);
+            SetPixelChannel(blur_image,channel,p[center+i],q);
+            continue;
           }
-          gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelRed(blur_image,ClampToQuantum(gamma*pixel.red),q);
-          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelGreen(blur_image,ClampToQuantum(gamma*pixel.green),q);
-          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelBlue(blur_image,ClampToQuantum(gamma*pixel.blue),q);
-          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-              (blur_image->colorspace == CMYKColorspace))
-            SetPixelBlack(blur_image,ClampToQuantum(gamma*pixel.black),q);
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
+        k=kernel;
+        pixels=p;
+        pixel=0.0;
+        if ((blur_traits & BlendPixelTrait) == 0)
+          {
+            /*
+              No alpha blending.
+            */
+            for (u=0; u < (ssize_t) width; u++)
             {
-              k=kernel;
-              kernel_pixels=p;
-              for (i=0; i < (ssize_t) width; i++)
-              {
-                pixel.alpha+=(*k)*GetPixelAlpha(image,kernel_pixels);
-                k++;
-                kernel_pixels+=GetPixelChannels(image);
-              }
-              SetPixelAlpha(blur_image,ClampToQuantum(pixel.alpha),q);
+              pixel+=(*k)*pixels[i];
+              k++;
+              pixels+=GetPixelChannels(image);
             }
+            SetPixelChannel(blur_image,channel,ClampToQuantum(pixel),q);
+            continue;
+          }
+        /*
+          Alpha blending.
+        */
+        gamma=0.0;
+        for (u=0; u < (ssize_t) width; u++)
+        {
+          alpha=(double) (QuantumScale*GetPixelAlpha(image,pixels));
+          pixel+=(*k)*alpha*pixels[i];
+          gamma+=(*k)*alpha;
+          k++;
+          pixels+=GetPixelChannels(image);
         }
+        gamma=MagickEpsilonReciprocal(gamma);
+        SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
+      }
       p+=GetPixelChannels(image);
       q+=GetPixelChannels(blur_image);
     }
@@ -1074,7 +1058,7 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_BlurImage)
+        #pragma omp critical (MagickCore_BlurImage)
 #endif
         proceed=SetImageProgress(image,BlurImageTag,progress++,blur_image->rows+
           blur_image->columns);
@@ -1087,10 +1071,12 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
   /*
     Blur columns.
   */
-  image_view=AcquireCacheView(blur_image);
-  blur_view=AcquireCacheView(blur_image);
+  center=(ssize_t) GetPixelChannels(blur_image)*(width/2L);
+  image_view=AcquireVirtualCacheView(blur_image,exception);
+  blur_view=AcquireAuthenticCacheView(blur_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (x=0; x < (ssize_t) blur_image->columns; x++)
   {
@@ -1106,7 +1092,7 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
     if (status == MagickFalse)
       continue;
     p=GetCacheViewVirtualPixels(image_view,x,-((ssize_t) width/2L),1,
-      image->rows+width,exception);
+      blur_image->rows+width,exception);
     q=GetCacheViewAuthenticPixels(blur_view,x,0,1,blur_image->rows,exception);
     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
       {
@@ -1115,103 +1101,77 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
       }
     for (y=0; y < (ssize_t) blur_image->rows; y++)
     {
-      PixelInfo
-        pixel;
+      register ssize_t
+        i;
+
+      for (i=0; i < (ssize_t) GetPixelChannels(blur_image); i++)
+      {
+        double
+          alpha,
+          gamma,
+          pixel;
 
-      register const double
-        *restrict k;
+        PixelChannel
+          channel;
 
-      register const Quantum
-        *restrict kernel_pixels;
+        PixelTrait
+          blur_traits,
+          traits;
 
-      register ssize_t
-        i;
+        register const MagickRealType
+          *restrict k;
 
-      pixel.red=bias;
-      pixel.green=bias;
-      pixel.blue=bias;
-      pixel.black=bias;
-      pixel.alpha=bias;
-      k=kernel;
-      kernel_pixels=p;
-      if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) == 0) ||
-          (blur_image->matte == MagickFalse))
-        {
-          for (i=0; i < (ssize_t) width; i++)
-          {
-            pixel.red+=(*k)*GetPixelRed(blur_image,kernel_pixels);
-            pixel.green+=(*k)*GetPixelGreen(blur_image,kernel_pixels);
-            pixel.blue+=(*k)*GetPixelBlue(blur_image,kernel_pixels);
-            if (blur_image->colorspace == CMYKColorspace)
-              pixel.black+=(*k)*GetPixelBlack(blur_image,kernel_pixels);
-            k++;
-            kernel_pixels+=GetPixelChannels(blur_image);
-          }
-          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelRed(blur_image,ClampToQuantum(pixel.red),q);
-          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelGreen(blur_image,ClampToQuantum(pixel.green),q);
-          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelBlue(blur_image,ClampToQuantum(pixel.blue),q);
-          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-              (blur_image->colorspace == CMYKColorspace))
-            SetPixelBlack(blur_image,ClampToQuantum(pixel.black),q);
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
-            {
-              k=kernel;
-              kernel_pixels=p;
-              for (i=0; i < (ssize_t) width; i++)
-              {
-                pixel.alpha+=(*k)*GetPixelAlpha(blur_image,kernel_pixels);
-                k++;
-                kernel_pixels+=GetPixelChannels(blur_image);
-              }
-              SetPixelAlpha(blur_image,ClampToQuantum(pixel.alpha),q);
-            }
-        }
-      else
-        {
-          MagickRealType
-            alpha,
-            gamma;
+        register const Quantum
+          *restrict pixels;
+
+        register ssize_t
+          u;
 
-          gamma=0.0;
-          for (i=0; i < (ssize_t) width; i++)
+        channel=GetPixelChannelChannel(blur_image,i);
+        traits=GetPixelChannelTraits(blur_image,channel);
+        blur_traits=GetPixelChannelTraits(blur_image,channel);
+        if ((traits == UndefinedPixelTrait) ||
+            (blur_traits == UndefinedPixelTrait))
+          continue;
+        if (((blur_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
-            alpha=(MagickRealType) (QuantumScale*
-              GetPixelAlpha(blur_image,kernel_pixels));
-            pixel.red+=(*k)*alpha*GetPixelRed(blur_image,kernel_pixels);
-            pixel.green+=(*k)*alpha*GetPixelGreen(blur_image,kernel_pixels);
-            pixel.blue+=(*k)*alpha*GetPixelBlue(blur_image,kernel_pixels);
-            if (blur_image->colorspace == CMYKColorspace)
-              pixel.black+=(*k)*alpha*GetPixelBlack(blur_image,kernel_pixels);
-            gamma+=(*k)*alpha;
-            k++;
-            kernel_pixels+=GetPixelChannels(blur_image);
+            SetPixelChannel(blur_image,channel,p[center+i],q);
+            continue;
           }
-          gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelRed(blur_image,ClampToQuantum(gamma*pixel.red),q);
-          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelGreen(blur_image,ClampToQuantum(gamma*pixel.green),q);
-          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelBlue(blur_image,ClampToQuantum(gamma*pixel.blue),q);
-          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-              (blur_image->colorspace == CMYKColorspace))
-            SetPixelBlack(blur_image,ClampToQuantum(gamma*pixel.black),q);
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
+        k=kernel;
+        pixels=p;
+        pixel=0.0;
+        if ((blur_traits & BlendPixelTrait) == 0)
+          {
+            /*
+              No alpha blending.
+            */
+            for (u=0; u < (ssize_t) width; u++)
             {
-              k=kernel;
-              kernel_pixels=p;
-              for (i=0; i < (ssize_t) width; i++)
-              {
-                pixel.alpha+=(*k)*GetPixelAlpha(blur_image,kernel_pixels);
-                k++;
-                kernel_pixels+=GetPixelChannels(blur_image);
-              }
-              SetPixelAlpha(blur_image,ClampToQuantum(pixel.alpha),q);
+              pixel+=(*k)*pixels[i];
+              k++;
+              pixels+=GetPixelChannels(blur_image);
             }
+            SetPixelChannel(blur_image,channel,ClampToQuantum(pixel),q);
+            continue;
+          }
+        /*
+          Alpha blending.
+        */
+        gamma=0.0;
+        for (u=0; u < (ssize_t) width; u++)
+        {
+          alpha=(double) (QuantumScale*GetPixelAlpha(blur_image,
+            pixels));
+          pixel+=(*k)*alpha*pixels[i];
+          gamma+=(*k)*alpha;
+          k++;
+          pixels+=GetPixelChannels(blur_image);
         }
+        gamma=MagickEpsilonReciprocal(gamma);
+        SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
+      }
       p+=GetPixelChannels(blur_image);
       q+=GetPixelChannels(blur_image);
     }
@@ -1223,7 +1183,7 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_BlurImage)
+        #pragma omp critical (MagickCore_BlurImage)
 #endif
         proceed=SetImageProgress(blur_image,BlurImageTag,progress++,
           blur_image->rows+blur_image->columns);
@@ -1233,10 +1193,10 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
   }
   blur_view=DestroyCacheView(blur_view);
   image_view=DestroyCacheView(image_view);
-  kernel=(double *) RelinquishMagickMemory(kernel);
+  kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
+  blur_image->type=image->type;
   if (status == MagickFalse)
     blur_image=DestroyImage(blur_image);
-  blur_image->type=image->type;
   return(blur_image);
 }
 \f
@@ -1270,220 +1230,299 @@ MagickExport Image *BlurImage(const Image *image,const double radius,
 MagickExport Image *ConvolveImage(const Image *image,
   const KernelInfo *kernel_info,ExceptionInfo *exception)
 {
-#define ConvolveImageTag  "Convolve/Image"
+  return(MorphologyImage(image,CorrelateMorphology,1,kernel_info,exception));
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%     D e s p e c k l e I m a g e                                             %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  DespeckleImage() reduces the speckle noise in an image while perserving the
+%  edges of the original image.  A speckle removing filter uses a complementary %  hulling technique (raising pixels that are darker than their surrounding
+%  neighbors, then complementarily lowering pixels that are brighter than their
+%  surrounding neighbors) to reduce the speckle index of that image (reference
+%  Crimmins speckle removal).
+%
+%  The format of the DespeckleImage method is:
+%
+%      Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+
+static void Hull(const Image *image,const ssize_t x_offset,
+  const ssize_t y_offset,const size_t columns,const size_t rows,
+  const int polarity,Quantum *restrict f,Quantum *restrict g)
+{
+  register Quantum
+    *p,
+    *q,
+    *r,
+    *s;
+
+  ssize_t
+    y;
+
+  assert(f != (Quantum *) NULL);
+  assert(g != (Quantum *) NULL);
+  p=f+(columns+2);
+  q=g+(columns+2);
+  r=p+(y_offset*(columns+2)+x_offset);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static) \
+    dynamic_number_threads(image,columns,rows,1)
+#endif
+  for (y=0; y < (ssize_t) rows; y++)
+  {
+    register ssize_t
+      i,
+      x;
+
+    SignedQuantum
+      v;
+
+    i=(2*y+1)+y*columns;
+    if (polarity > 0)
+      for (x=0; x < (ssize_t) columns; x++)
+      {
+        v=(SignedQuantum) p[i];
+        if ((SignedQuantum) r[i] >= (v+ScaleCharToQuantum(2)))
+          v+=ScaleCharToQuantum(1);
+        q[i]=(Quantum) v;
+        i++;
+      }
+    else
+      for (x=0; x < (ssize_t) columns; x++)
+      {
+        v=(SignedQuantum) p[i];
+        if ((SignedQuantum) r[i] <= (v-ScaleCharToQuantum(2)))
+          v-=ScaleCharToQuantum(1);
+        q[i]=(Quantum) v;
+        i++;
+      }
+  }
+  p=f+(columns+2);
+  q=g+(columns+2);
+  r=q+(y_offset*(columns+2)+x_offset);
+  s=q-(y_offset*(columns+2)+x_offset);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static) \
+    dynamic_number_threads(image,columns,rows,1)
+#endif
+  for (y=0; y < (ssize_t) rows; y++)
+  {
+    register ssize_t
+      i,
+      x;
+
+    SignedQuantum
+      v;
+
+    i=(2*y+1)+y*columns;
+    if (polarity > 0)
+      for (x=0; x < (ssize_t) columns; x++)
+      {
+        v=(SignedQuantum) q[i];
+        if (((SignedQuantum) s[i] >= (v+ScaleCharToQuantum(2))) &&
+            ((SignedQuantum) r[i] > v))
+          v+=ScaleCharToQuantum(1);
+        p[i]=(Quantum) v;
+        i++;
+      }
+    else
+      for (x=0; x < (ssize_t) columns; x++)
+      {
+        v=(SignedQuantum) q[i];
+        if (((SignedQuantum) s[i] <= (v-ScaleCharToQuantum(2))) &&
+            ((SignedQuantum) r[i] < v))
+          v-=ScaleCharToQuantum(1);
+        p[i]=(Quantum) v;
+        i++;
+      }
+  }
+}
+
+MagickExport Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
+{
+#define DespeckleImageTag  "Despeckle/Image"
 
   CacheView
-    *convolve_view,
+    *despeckle_view,
     *image_view;
 
   Image
-    *convolve_image;
+    *despeckle_image;
 
   MagickBooleanType
     status;
 
-  MagickOffsetType
-    progress;
+  Quantum
+    *restrict buffer,
+    *restrict pixels;
 
-  ssize_t
-    center,
-    y;
+  register ssize_t
+    i;
+
+  size_t
+    length;
+
+  static const ssize_t
+    X[4] = {0, 1, 1,-1},
+    Y[4] = {1, 0, 1, 1};
 
   /*
-    Initialize convolve image attributes.
+    Allocate despeckled image.
   */
-  assert(image != (Image *) NULL);
+  assert(image != (const Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
-  if ((kernel_info->width % 2) == 0)
-    ThrowImageException(OptionError,"KernelWidthMustBeAnOddNumber");
-  convolve_image=CloneImage(image,image->columns,image->rows,MagickTrue,
-    exception);
-  if (convolve_image == (Image *) NULL)
+  despeckle_image=CloneImage(image,0,0,MagickTrue,exception);
+  if (despeckle_image == (Image *) NULL)
     return((Image *) NULL);
-  if (SetImageStorageClass(convolve_image,DirectClass,exception) == MagickFalse)
+  status=SetImageStorageClass(despeckle_image,DirectClass,exception);
+  if (status == MagickFalse)
     {
-      convolve_image=DestroyImage(convolve_image);
+      despeckle_image=DestroyImage(despeckle_image);
       return((Image *) NULL);
     }
-  if (image->debug != MagickFalse)
+  /*
+    Allocate image buffer.
+  */
+  length=(size_t) ((image->columns+2)*(image->rows+2));
+  pixels=(Quantum *) AcquireQuantumMemory(length,sizeof(*pixels));
+  buffer=(Quantum *) AcquireQuantumMemory(length,sizeof(*buffer));
+  if ((pixels == (Quantum *) NULL) || (buffer == (Quantum *) NULL))
     {
-      char
-        format[MaxTextExtent],
-        *message;
-
-      register const double
-        *k;
-
-      register ssize_t
-        u;
-
-      ssize_t
-        v;
-
-      (void) LogMagickEvent(TransformEvent,GetMagickModule(),
-        "  ConvolveImage with %.20gx%.20g kernel:",(double) kernel_info->width,
-        (double) kernel_info->height);
-      message=AcquireString("");
-      k=kernel_info->values;
-      for (v=0; v < (ssize_t) kernel_info->width; v++)
-      {
-        *message='\0';
-        (void) FormatLocaleString(format,MaxTextExtent,"%.20g: ",(double) v);
-        (void) ConcatenateString(&message,format);
-        for (u=0; u < (ssize_t) kernel_info->height; u++)
-        {
-          (void) FormatLocaleString(format,MaxTextExtent,"%g ",*k++);
-          (void) ConcatenateString(&message,format);
-        }
-        (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
-      }
-      message=DestroyString(message);
+      if (buffer != (Quantum *) NULL)
+        buffer=(Quantum *) RelinquishMagickMemory(buffer);
+      if (pixels != (Quantum *) NULL)
+        pixels=(Quantum *) RelinquishMagickMemory(pixels);
+      despeckle_image=DestroyImage(despeckle_image);
+      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
     }
   /*
-    Convolve image.
+    Reduce speckle in the image.
   */
-  center=(ssize_t) GetPixelChannels(image)*(image->columns+kernel_info->width)*
-    (kernel_info->height/2L)+GetPixelChannels(image)*(kernel_info->width/2);
   status=MagickTrue;
-  progress=0;
-  image_view=AcquireCacheView(image);
-  convolve_view=AcquireCacheView(convolve_image);
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
-#endif
-  for (y=0; y < (ssize_t) image->rows; y++)
+  image_view=AcquireVirtualCacheView(image,exception);
+  despeckle_view=AcquireAuthenticCacheView(despeckle_image,exception);
+  for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
   {
-    register const Quantum
-      *restrict p;
+    PixelChannel
+       channel;
 
-    register Quantum
-      *restrict q;
+    PixelTrait
+      despeckle_traits,
+      traits;
 
     register ssize_t
+      k,
       x;
 
+    ssize_t
+      j,
+      y;
+
     if (status == MagickFalse)
       continue;
-    p=GetCacheViewVirtualPixels(image_view,-((ssize_t) kernel_info->width/2L),y-
-      (ssize_t) (kernel_info->height/2L),image->columns+kernel_info->width,
-      kernel_info->height,exception);
-    q=QueueCacheViewAuthenticPixels(convolve_view,0,y,convolve_image->columns,1,
-      exception);
-    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
-      {
-        status=MagickFalse;
-        continue;
-      }
-    for (x=0; x < (ssize_t) image->columns; x++)
+    channel=GetPixelChannelChannel(image,i);
+    traits=GetPixelChannelTraits(image,channel);
+    despeckle_traits=GetPixelChannelTraits(despeckle_image,channel);
+    if ((traits == UndefinedPixelTrait) ||
+        (despeckle_traits == UndefinedPixelTrait))
+      continue;
+    if ((despeckle_traits & CopyPixelTrait) != 0)
+      continue;
+    (void) ResetMagickMemory(pixels,0,length*sizeof(*pixels));
+    j=(ssize_t) image->columns+2;
+    for (y=0; y < (ssize_t) image->rows; y++)
     {
-      register ssize_t
-        i;
+      register const Quantum
+        *restrict p;
 
-      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+      if (p == (const Quantum *) NULL)
+        {
+          status=MagickFalse;
+          continue;
+        }
+      j++;
+      for (x=0; x < (ssize_t) image->columns; x++)
       {
-        MagickRealType
-          alpha,
-          gamma,
-          pixel;
-
-        PixelChannel
-          channel;
-
-        PixelTrait
-          convolve_traits,
-          traits;
-
-        register const double
-          *restrict k;
-
-        register const Quantum
-          *restrict pixels;
-
-        register ssize_t
-          u;
+        pixels[j++]=p[i];
+        p+=GetPixelChannels(image);
+      }
+      j++;
+    }
+    (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
+    for (k=0; k < 4; k++)
+    {
+      Hull(image,X[k],Y[k],image->columns,image->rows,1,pixels,buffer);
+      Hull(image,-X[k],-Y[k],image->columns,image->rows,1,pixels,buffer);
+      Hull(image,-X[k],-Y[k],image->columns,image->rows,-1,pixels,buffer);
+      Hull(image,X[k],Y[k],image->columns,image->rows,-1,pixels,buffer);
+    }
+    j=(ssize_t) image->columns+2;
+    for (y=0; y < (ssize_t) image->rows; y++)
+    {
+      MagickBooleanType
+        sync;
 
-        ssize_t
-          v;
+      register Quantum
+        *restrict q;
 
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        convolve_traits=GetPixelChannelMapTraits(convolve_image,channel);
-        if ((traits == UndefinedPixelTrait) ||
-            (convolve_traits == UndefinedPixelTrait))
-          continue;
-        if ((convolve_traits & CopyPixelTrait) != 0)
-          {
-            q[channel]=p[center+i];
-            continue;
-          }
-        k=kernel_info->values;
-        pixels=p;
-        pixel=kernel_info->bias;
-        if ((convolve_traits & BlendPixelTrait) == 0)
-          {
-            /*
-              No alpha blending.
-            */
-            for (v=0; v < (ssize_t) kernel_info->height; v++)
-            {
-              for (u=0; u < (ssize_t) kernel_info->width; u++)
-              {
-                pixel+=(*k)*pixels[i];
-                k++;
-                pixels+=GetPixelChannels(image);
-              }
-              pixels+=image->columns*GetPixelChannels(image);
-            }
-            q[channel]=ClampToQuantum(pixel);
-            continue;
-          }
-        /*
-          Alpha blending.
-        */
-        gamma=0.0;
-        for (v=0; v < (ssize_t) kernel_info->height; v++)
+      q=QueueCacheViewAuthenticPixels(despeckle_view,0,y,
+        despeckle_image->columns,1,exception);
+      if (q == (Quantum *) NULL)
         {
-          for (u=0; u < (ssize_t) kernel_info->width; u++)
-          {
-            alpha=(MagickRealType) (QuantumScale*GetPixelAlpha(image,pixels));
-            pixel+=(*k)*alpha*pixels[i];
-            gamma+=(*k)*alpha;
-            k++;
-            pixels+=GetPixelChannels(image);
-          }
-          pixels+=image->columns*GetPixelChannels(image);
+          status=MagickFalse;
+          continue;
         }
-        gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-        q[channel]=ClampToQuantum(gamma*pixel);
+      j++;
+      for (x=0; x < (ssize_t) image->columns; x++)
+      {
+        SetPixelChannel(despeckle_image,channel,pixels[j++],q);
+        q+=GetPixelChannels(despeckle_image);
       }
-      p+=GetPixelChannels(image);
-      q+=GetPixelChannels(convolve_image);
+      sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
+      if (sync == MagickFalse)
+        status=MagickFalse;
+      j++;
     }
-    if (SyncCacheViewAuthenticPixels(convolve_view,exception) == MagickFalse)
-      status=MagickFalse;
     if (image->progress_monitor != (MagickProgressMonitor) NULL)
       {
         MagickBooleanType
           proceed;
 
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_ConvolveImage)
-#endif
-        proceed=SetImageProgress(image,ConvolveImageTag,progress++,image->rows);
+        proceed=SetImageProgress(image,DespeckleImageTag,(MagickOffsetType) i,
+          GetPixelChannels(image));
         if (proceed == MagickFalse)
           status=MagickFalse;
       }
   }
-  convolve_image->type=image->type;
-  convolve_view=DestroyCacheView(convolve_view);
+  despeckle_view=DestroyCacheView(despeckle_view);
   image_view=DestroyCacheView(image_view);
+  buffer=(Quantum *) RelinquishMagickMemory(buffer);
+  pixels=(Quantum *) RelinquishMagickMemory(pixels);
+  despeckle_image->type=image->type;
   if (status == MagickFalse)
-    convolve_image=DestroyImage(convolve_image);
-  return(convolve_image);
+    despeckle_image=DestroyImage(despeckle_image);
+  return(despeckle_image);
 }
 \f
 /*
@@ -1491,375 +1530,71 @@ MagickExport Image *ConvolveImage(const Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%     D e s p e c k l e I m a g e                                             %
+%     E d g e I m a g e                                                       %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  DespeckleImage() reduces the speckle noise in an image while perserving the
-%  edges of the original image.
+%  EdgeImage() finds edges in an image.  Radius defines the radius of the
+%  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
+%  radius for you.
 %
-%  The format of the DespeckleImage method is:
+%  The format of the EdgeImage method is:
 %
-%      Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
+%      Image *EdgeImage(const Image *image,const double radius,
+%        const double sigma,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %    o image: the image.
 %
+%    o radius: the radius of the pixel neighborhood.
+%
+%    o sigma: the standard deviation of the Gaussian, in pixels.
+%
 %    o exception: return any errors or warnings in this structure.
 %
 */
-
-static void Hull(const ssize_t x_offset,const ssize_t y_offset,
-  const size_t columns,const size_t rows,Quantum *f,Quantum *g,
-  const int polarity)
-{
-  MagickRealType
-    v;
-
-  register Quantum
-    *p,
-    *q,
-    *r,
-    *s;
-
-  register ssize_t
-    x;
-
-  ssize_t
-    y;
-
-  assert(f != (Quantum *) NULL);
-  assert(g != (Quantum *) NULL);
-  p=f+(columns+2);
-  q=g+(columns+2);
-  r=p+(y_offset*((ssize_t) columns+2)+x_offset);
-  for (y=0; y < (ssize_t) rows; y++)
-  {
-    p++;
-    q++;
-    r++;
-    if (polarity > 0)
-      for (x=(ssize_t) columns; x != 0; x--)
-      {
-        v=(MagickRealType) (*p);
-        if ((MagickRealType) *r >= (v+(MagickRealType) ScaleCharToQuantum(2)))
-          v+=ScaleCharToQuantum(1);
-        *q=(Quantum) v;
-        p++;
-        q++;
-        r++;
-      }
-    else
-      for (x=(ssize_t) columns; x != 0; x--)
-      {
-        v=(MagickRealType) (*p);
-        if ((MagickRealType) *r <= (v-(MagickRealType) ScaleCharToQuantum(2)))
-          v-=(ssize_t) ScaleCharToQuantum(1);
-        *q=(Quantum) v;
-        p++;
-        q++;
-        r++;
-      }
-    p++;
-    q++;
-    r++;
-  }
-  p=f+(columns+2);
-  q=g+(columns+2);
-  r=q+(y_offset*((ssize_t) columns+2)+x_offset);
-  s=q-(y_offset*((ssize_t) columns+2)+x_offset);
-  for (y=0; y < (ssize_t) rows; y++)
-  {
-    p++;
-    q++;
-    r++;
-    s++;
-    if (polarity > 0)
-      for (x=(ssize_t) columns; x != 0; x--)
-      {
-        v=(MagickRealType) (*q);
-        if (((MagickRealType) *s >=
-             (v+(MagickRealType) ScaleCharToQuantum(2))) &&
-            ((MagickRealType) *r > v))
-          v+=ScaleCharToQuantum(1);
-        *p=(Quantum) v;
-        p++;
-        q++;
-        r++;
-        s++;
-      }
-    else
-      for (x=(ssize_t) columns; x != 0; x--)
-      {
-        v=(MagickRealType) (*q);
-        if (((MagickRealType) *s <=
-             (v-(MagickRealType) ScaleCharToQuantum(2))) &&
-            ((MagickRealType) *r < v))
-          v-=(MagickRealType) ScaleCharToQuantum(1);
-        *p=(Quantum) v;
-        p++;
-        q++;
-        r++;
-        s++;
-      }
-    p++;
-    q++;
-    r++;
-    s++;
-  }
-}
-
-MagickExport Image *DespeckleImage(const Image *image,ExceptionInfo *exception)
+MagickExport Image *EdgeImage(const Image *image,const double radius,
+  const double sigma,ExceptionInfo *exception)
 {
-#define DespeckleImageTag  "Despeckle/Image"
-
-  CacheView
-    *despeckle_view,
-    *image_view;
-
   Image
-    *despeckle_image;
+    *edge_image;
 
-  MagickBooleanType
-    status;
+  KernelInfo
+    *kernel_info;
 
   register ssize_t
     i;
 
-  Quantum
-    *restrict buffers,
-    *restrict pixels;
-
   size_t
-    length,
-    number_channels;
+    width;
 
-  static const ssize_t
-    X[4] = {0, 1, 1,-1},
-    Y[4] = {1, 0, 1, 1};
+  ssize_t
+    j,
+    u,
+    v;
 
-  /*
-    Allocate despeckled image.
-  */
   assert(image != (const Image *) NULL);
   assert(image->signature == MagickSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
-  despeckle_image=CloneImage(image,image->columns,image->rows,MagickTrue,
-    exception);
-  if (despeckle_image == (Image *) NULL)
-    return((Image *) NULL);
-  if (SetImageStorageClass(despeckle_image,DirectClass,exception) == MagickFalse)
+  width=GetOptimalKernelWidth1D(radius,sigma);
+  kernel_info=AcquireKernelInfo((const char *) NULL);
+  if (kernel_info == (KernelInfo *) NULL)
+    ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
+  kernel_info->width=width;
+  kernel_info->height=width;
+  kernel_info->values=(MagickRealType *) MagickAssumeAligned(
+    AcquireAlignedMemory(kernel_info->width,kernel_info->width*
+    sizeof(*kernel_info->values)));
+  if (kernel_info->values == (MagickRealType *) NULL)
     {
-      despeckle_image=DestroyImage(despeckle_image);
-      return((Image *) NULL);
-    }
-  /*
-    Allocate image buffers.
-  */
-  length=(size_t) ((image->columns+2)*(image->rows+2));
-  pixels=(Quantum *) AcquireQuantumMemory(length,2*sizeof(*pixels));
-  buffers=(Quantum *) AcquireQuantumMemory(length,2*sizeof(*pixels));
-  if ((pixels == (Quantum *) NULL) || (buffers == (Quantum *) NULL))
-    {
-      if (buffers != (Quantum *) NULL)
-        buffers=(Quantum *) RelinquishMagickMemory(buffers);
-      if (pixels != (Quantum *) NULL)
-        pixels=(Quantum *) RelinquishMagickMemory(pixels);
-      despeckle_image=DestroyImage(despeckle_image);
-      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-    }
-  /*
-    Reduce speckle in the image.
-  */
-  status=MagickTrue;
-  number_channels=(size_t) (image->colorspace == CMYKColorspace ? 5 : 4);
-  image_view=AcquireCacheView(image);
-  despeckle_view=AcquireCacheView(despeckle_image);
-  for (i=0; i < (ssize_t) number_channels; i++)
-  {
-    register Quantum
-      *buffer,
-      *pixel;
-
-    register ssize_t
-      k,
-      x;
-
-    ssize_t
-      j,
-      y;
-
-    if (status == MagickFalse)
-      continue;
-    pixel=pixels;
-    (void) ResetMagickMemory(pixel,0,length*sizeof(*pixel));
-    buffer=buffers;
-    j=(ssize_t) image->columns+2;
-    for (y=0; y < (ssize_t) image->rows; y++)
-    {
-      register const Quantum
-        *restrict p;
-
-      p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
-      if (p == (const Quantum *) NULL)
-        break;
-      j++;
-      for (x=0; x < (ssize_t) image->columns; x++)
-      {
-        switch (i)
-        {
-          case 0: pixel[j]=GetPixelRed(image,p); break;
-          case 1: pixel[j]=GetPixelGreen(image,p); break;
-          case 2: pixel[j]=GetPixelBlue(image,p); break;
-          case 3: pixel[j]=GetPixelAlpha(image,p); break;
-          case 4: pixel[j]=GetPixelBlack(image,p); break;
-          default: break;
-        }
-        p+=GetPixelChannels(image);
-        j++;
-      }
-      j++;
-    }
-    (void) ResetMagickMemory(buffer,0,length*sizeof(*buffer));
-    for (k=0; k < 4; k++)
-    {
-      Hull(X[k],Y[k],image->columns,image->rows,pixel,buffer,1);
-      Hull(-X[k],-Y[k],image->columns,image->rows,pixel,buffer,1);
-      Hull(-X[k],-Y[k],image->columns,image->rows,pixel,buffer,-1);
-      Hull(X[k],Y[k],image->columns,image->rows,pixel,buffer,-1);
-    }
-    j=(ssize_t) image->columns+2;
-    for (y=0; y < (ssize_t) image->rows; y++)
-    {
-      MagickBooleanType
-        sync;
-
-      register Quantum
-        *restrict q;
-
-      q=GetCacheViewAuthenticPixels(despeckle_view,0,y,despeckle_image->columns,
-        1,exception);
-      if (q == (Quantum *) NULL)
-        break;
-      j++;
-      for (x=0; x < (ssize_t) image->columns; x++)
-      {
-        switch (i)
-        {
-          case 0: SetPixelRed(despeckle_image,pixel[j],q); break;
-          case 1: SetPixelGreen(despeckle_image,pixel[j],q); break;
-          case 2: SetPixelBlue(despeckle_image,pixel[j],q); break;
-          case 3: SetPixelAlpha(despeckle_image,pixel[j],q); break;
-          case 4: SetPixelBlack(despeckle_image,pixel[j],q); break;
-          default: break;
-        }
-        q+=GetPixelChannels(despeckle_image);
-        j++;
-      }
-      sync=SyncCacheViewAuthenticPixels(despeckle_view,exception);
-      if (sync == MagickFalse)
-        {
-          status=MagickFalse;
-          break;
-        }
-      j++;
-    }
-    if (image->progress_monitor != (MagickProgressMonitor) NULL)
-      {
-        MagickBooleanType
-          proceed;
-
-        proceed=SetImageProgress(image,DespeckleImageTag,(MagickOffsetType) i,
-          number_channels);
-        if (proceed == MagickFalse)
-          status=MagickFalse;
-      }
-  }
-  despeckle_view=DestroyCacheView(despeckle_view);
-  image_view=DestroyCacheView(image_view);
-  buffers=(Quantum *) RelinquishMagickMemory(buffers);
-  pixels=(Quantum *) RelinquishMagickMemory(pixels);
-  despeckle_image->type=image->type;
-  if (status == MagickFalse)
-    despeckle_image=DestroyImage(despeckle_image);
-  return(despeckle_image);
-}
-\f
-/*
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%     E d g e I m a g e                                                       %
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%  EdgeImage() finds edges in an image.  Radius defines the radius of the
-%  convolution filter.  Use a radius of 0 and EdgeImage() selects a suitable
-%  radius for you.
-%
-%  The format of the EdgeImage method is:
-%
-%      Image *EdgeImage(const Image *image,const double radius,
-%        const double sigma,ExceptionInfo *exception)
-%
-%  A description of each parameter follows:
-%
-%    o image: the image.
-%
-%    o radius: the radius of the pixel neighborhood.
-%
-%    o sigma: the standard deviation of the Gaussian, in pixels.
-%
-%    o exception: return any errors or warnings in this structure.
-%
-*/
-MagickExport Image *EdgeImage(const Image *image,const double radius,
-  const double sigma,ExceptionInfo *exception)
-{
-  Image
-    *edge_image;
-
-  KernelInfo
-    *kernel_info;
-
-  register ssize_t
-    i;
-
-  size_t
-    width;
-
-  ssize_t
-    j,
-    u,
-    v;
-
-  assert(image != (const Image *) NULL);
-  assert(image->signature == MagickSignature);
-  if (image->debug != MagickFalse)
-    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  assert(exception != (ExceptionInfo *) NULL);
-  assert(exception->signature == MagickSignature);
-  width=GetOptimalKernelWidth2D(radius,sigma);
-  kernel_info=AcquireKernelInfo((const char *) NULL);
-  if (kernel_info == (KernelInfo *) NULL)
-    ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-  kernel_info->width=width;
-  kernel_info->height=width;
-  kernel_info->values=(double *) AcquireAlignedMemory(kernel_info->width,
-    kernel_info->width*sizeof(*kernel_info->values));
-  if (kernel_info->values == (double *) NULL)
-    {
-      kernel_info=DestroyKernelInfo(kernel_info);
-      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
+      kernel_info=DestroyKernelInfo(kernel_info);
+      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
     }
   j=(ssize_t) kernel_info->width/2;
   i=0;
@@ -1867,13 +1602,14 @@ MagickExport Image *EdgeImage(const Image *image,const double radius,
   {
     for (u=(-j); u <= j; u++)
     {
-      kernel_info->values[i]=(-1.0);
+      kernel_info->values[i]=(MagickRealType) (-1.0);
       i++;
     }
   }
-  kernel_info->values[i/2]=(double) (width*width-1.0);
-  kernel_info->bias=image->bias;
+  kernel_info->values[i/2]=(MagickRealType) (width*width-1.0);
   edge_image=ConvolveImage(image,kernel_info,exception);
+  if (edge_image != (Image *) NULL)
+    (void) ClampImage(edge_image,exception);
   kernel_info=DestroyKernelInfo(kernel_info);
   return(edge_image);
 }
@@ -1938,15 +1674,16 @@ MagickExport Image *EmbossImage(const Image *image,const double radius,
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
-  width=GetOptimalKernelWidth2D(radius,sigma);
+  width=GetOptimalKernelWidth1D(radius,sigma);
   kernel_info=AcquireKernelInfo((const char *) NULL);
   if (kernel_info == (KernelInfo *) NULL)
     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
   kernel_info->width=width;
   kernel_info->height=width;
-  kernel_info->values=(double *) AcquireAlignedMemory(kernel_info->width,
-    kernel_info->width*sizeof(*kernel_info->values));
-  if (kernel_info->values == (double *) NULL)
+  kernel_info->values=(MagickRealType *) MagickAssumeAligned(
+    AcquireAlignedMemory(kernel_info->width,kernel_info->width*
+    sizeof(*kernel_info->values)));
+  if (kernel_info->values == (MagickRealType *) NULL)
     {
       kernel_info=DestroyKernelInfo(kernel_info);
       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
@@ -1958,8 +1695,8 @@ MagickExport Image *EmbossImage(const Image *image,const double radius,
   {
     for (u=(-j); u <= j; u++)
     {
-      kernel_info->values[i]=(double) (((u < 0) || (v < 0) ? -8.0 : 8.0)*
-        exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
+      kernel_info->values[i]=(MagickRealType) (((u < 0) || (v < 0) ? -8.0 :
+        8.0)*exp(-((double) u*u+v*v)/(2.0*MagickSigma*MagickSigma))/
         (2.0*MagickPI*MagickSigma*MagickSigma));
       if (u != k)
         kernel_info->values[i]=0.0;
@@ -1967,7 +1704,6 @@ MagickExport Image *EmbossImage(const Image *image,const double radius,
     }
     k--;
   }
-  kernel_info->bias=image->bias;
   emboss_image=ConvolveImage(image,kernel_info,exception);
   kernel_info=DestroyKernelInfo(kernel_info);
   if (emboss_image != (Image *) NULL)
@@ -1994,7 +1730,7 @@ MagickExport Image *EmbossImage(const Image *image,const double radius,
 %  The format of the GaussianBlurImage method is:
 %
 %      Image *GaussianBlurImage(const Image *image,onst double radius,
-%        const double sigma,const double bias,ExceptionInfo *exception)
+%        const double sigma,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
@@ -2005,13 +1741,11 @@ MagickExport Image *EmbossImage(const Image *image,const double radius,
 %
 %    o sigma: the standard deviation of the Gaussian, in pixels.
 %
-%    o bias: the bias.
-%
 %    o exception: return any errors or warnings in this structure.
 %
 */
 MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
-  const double sigma,const double bias,ExceptionInfo *exception)
+  const double sigma,ExceptionInfo *exception)
 {
   Image
     *blur_image;
@@ -2043,11 +1777,11 @@ MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
   (void) ResetMagickMemory(kernel_info,0,sizeof(*kernel_info));
   kernel_info->width=width;
   kernel_info->height=width;
-  kernel_info->bias=bias;
   kernel_info->signature=MagickSignature;
-  kernel_info->values=(double *) AcquireAlignedMemory(kernel_info->width,
-    kernel_info->width*sizeof(*kernel_info->values));
-  if (kernel_info->values == (double *) NULL)
+  kernel_info->values=(MagickRealType *) MagickAssumeAligned(
+    AcquireAlignedMemory(kernel_info->width,kernel_info->width*
+    sizeof(*kernel_info->values)));
+  if (kernel_info->values == (MagickRealType *) NULL)
     {
       kernel_info=DestroyKernelInfo(kernel_info);
       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
@@ -2058,7 +1792,7 @@ MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
   {
     for (u=(-j); u <= j; u++)
     {
-      kernel_info->values[i]=(double) (exp(-((double) u*u+v*v)/(2.0*
+      kernel_info->values[i]=(MagickRealType) (exp(-((double) u*u+v*v)/(2.0*
         MagickSigma*MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
       i++;
     }
@@ -2107,9 +1841,10 @@ MagickExport Image *GaussianBlurImage(const Image *image,const double radius,
 %
 */
 
-static double *GetMotionBlurKernel(const size_t width,const double sigma)
+static MagickRealType *GetMotionBlurKernel(const size_t width,
+  const double sigma)
 {
-  double
+  MagickRealType
     *kernel,
     normalize;
 
@@ -2120,13 +1855,14 @@ static double *GetMotionBlurKernel(const size_t width,const double sigma)
    Generate a 1-D convolution kernel.
   */
   (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
-  kernel=(double *) AcquireQuantumMemory((size_t) width,sizeof(*kernel));
-  if (kernel == (double *) NULL)
+  kernel=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory((size_t)
+    width,sizeof(*kernel)));
+  if (kernel == (MagickRealType *) NULL)
     return(kernel);
   normalize=0.0;
   for (i=0; i < (ssize_t) width; i++)
   {
-    kernel[i]=(double) (exp((-((double) i*i)/(double) (2.0*MagickSigma*
+    kernel[i]=(MagickRealType) (exp((-((double) i*i)/(double) (2.0*MagickSigma*
       MagickSigma)))/(MagickSQ2PI*MagickSigma));
     normalize+=kernel[i];
   }
@@ -2135,16 +1871,13 @@ static double *GetMotionBlurKernel(const size_t width,const double sigma)
   return(kernel);
 }
 
-MagickExport Image *MotionBlurImage(const Image *image,
-  const double radius,const double sigma,const double angle,
-  ExceptionInfo *exception)
+MagickExport Image *MotionBlurImage(const Image *image,const double radius,
+  const double sigma,const double angle,ExceptionInfo *exception)
 {
   CacheView
     *blur_view,
-    *image_view;
-
-  double
-    *kernel;
+    *image_view,
+    *motion_view;
 
   Image
     *blur_image;
@@ -2155,8 +1888,8 @@ MagickExport Image *MotionBlurImage(const Image *image,
   MagickOffsetType
     progress;
 
-  PixelInfo
-    bias;
+  MagickRealType
+    *kernel;
 
   OffsetInfo
     *offset;
@@ -2180,24 +1913,24 @@ MagickExport Image *MotionBlurImage(const Image *image,
   assert(exception != (ExceptionInfo *) NULL);
   width=GetOptimalKernelWidth1D(radius,sigma);
   kernel=GetMotionBlurKernel(width,sigma);
-  if (kernel == (double *) NULL)
+  if (kernel == (MagickRealType *) NULL)
     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
   offset=(OffsetInfo *) AcquireQuantumMemory(width,sizeof(*offset));
   if (offset == (OffsetInfo *) NULL)
     {
-      kernel=(double *) RelinquishMagickMemory(kernel);
+      kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
     }
-  blur_image=CloneImage(image,0,0,MagickTrue,exception);
+  blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
   if (blur_image == (Image *) NULL)
     {
-      kernel=(double *) RelinquishMagickMemory(kernel);
+      kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
       return((Image *) NULL);
     }
   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
     {
-      kernel=(double *) RelinquishMagickMemory(kernel);
+      kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
       offset=(OffsetInfo *) RelinquishMagickMemory(offset);
       blur_image=DestroyImage(blur_image);
       return((Image *) NULL);
@@ -2214,14 +1947,18 @@ MagickExport Image *MotionBlurImage(const Image *image,
   */
   status=MagickTrue;
   progress=0;
-  GetPixelInfo(image,&bias);
-  image_view=AcquireCacheView(image);
-  blur_view=AcquireCacheView(blur_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  motion_view=AcquireVirtualCacheView(image,exception);
+  blur_view=AcquireAuthenticCacheView(blur_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
+    register const Quantum
+      *restrict p;
+
     register Quantum
       *restrict q;
 
@@ -2230,100 +1967,93 @@ MagickExport Image *MotionBlurImage(const Image *image,
 
     if (status == MagickFalse)
       continue;
-    q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+    q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
       exception);
-    if (q == (Quantum *) NULL)
+    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
       {
         status=MagickFalse;
         continue;
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
-      PixelInfo
-        qixel;
+      register ssize_t
+        i;
+
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        double
+          alpha,
+          gamma,
+          pixel;
+
+        PixelChannel
+          channel;
 
-      PixelPacket
-        pixel;
+        PixelTrait
+          blur_traits,
+          traits;
 
-      register double
-        *restrict k;
+        register const Quantum
+          *restrict r;
 
-      register ssize_t
-        i;
+        register MagickRealType
+          *restrict k;
 
-      k=kernel;
-      qixel=bias;
-      if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) == 0) || (image->matte == MagickFalse))
-        {
-          for (i=0; i < (ssize_t) width; i++)
+        register ssize_t
+          j;
+
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        blur_traits=GetPixelChannelTraits(blur_image,channel);
+        if ((traits == UndefinedPixelTrait) ||
+            (blur_traits == UndefinedPixelTrait))
+          continue;
+        if (((blur_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
-            (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
-              offset[i].y,&pixel,exception);
-            qixel.red+=(*k)*pixel.red;
-            qixel.green+=(*k)*pixel.green;
-            qixel.blue+=(*k)*pixel.blue;
-            qixel.alpha+=(*k)*pixel.alpha;
-            if (image->colorspace == CMYKColorspace)
-              qixel.black+=(*k)*pixel.black;
-            k++;
+            SetPixelChannel(blur_image,channel,p[i],q);
+            continue;
           }
-          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelRed(blur_image,
-              ClampToQuantum(qixel.red),q);
-          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelGreen(blur_image,
-              ClampToQuantum(qixel.green),q);
-          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelBlue(blur_image,
-              ClampToQuantum(qixel.blue),q);
-          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-              (image->colorspace == CMYKColorspace))
-            SetPixelBlack(blur_image,
-              ClampToQuantum(qixel.black),q);
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelAlpha(blur_image,
-              ClampToQuantum(qixel.alpha),q);
-        }
-      else
-        {
-          MagickRealType
-            alpha,
-            gamma;
-
-          alpha=0.0;
-          gamma=0.0;
-          for (i=0; i < (ssize_t) width; i++)
+        k=kernel;
+        pixel=0.0;
+        if ((blur_traits & BlendPixelTrait) == 0)
           {
-            (void) GetOneCacheViewVirtualPixel(image_view,x+offset[i].x,y+
-              offset[i].y,&pixel,exception);
-            alpha=(MagickRealType) (QuantumScale*pixel.alpha);
-            qixel.red+=(*k)*alpha*pixel.red;
-            qixel.green+=(*k)*alpha*pixel.green;
-            qixel.blue+=(*k)*alpha*pixel.blue;
-            qixel.alpha+=(*k)*pixel.alpha;
-            if (image->colorspace == CMYKColorspace)
-              qixel.black+=(*k)*alpha*pixel.black;
-            gamma+=(*k)*alpha;
-            k++;
+            for (j=0; j < (ssize_t) width; j++)
+            {
+              r=GetCacheViewVirtualPixels(motion_view,x+offset[j].x,y+
+                offset[j].y,1,1,exception);
+              if (r == (const Quantum *) NULL)
+                {
+                  status=MagickFalse;
+                  continue;
+                }
+              pixel+=(*k)*r[i];
+              k++;
+            }
+            SetPixelChannel(blur_image,channel,ClampToQuantum(pixel),q);
+            continue;
           }
-          gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelRed(blur_image,
-              ClampToQuantum(gamma*qixel.red),q);
-          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelGreen(blur_image,
-              ClampToQuantum(gamma*qixel.green),q);
-          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelBlue(blur_image,
-              ClampToQuantum(gamma*qixel.blue),q);
-          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-              (image->colorspace == CMYKColorspace))
-            SetPixelBlack(blur_image,
-              ClampToQuantum(gamma*qixel.black),q);
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelAlpha(blur_image,
-              ClampToQuantum(qixel.alpha),q);
+        alpha=0.0;
+        gamma=0.0;
+        for (j=0; j < (ssize_t) width; j++)
+        {
+          r=GetCacheViewVirtualPixels(motion_view,x+offset[j].x,y+offset[j].y,1,
+            1,exception);
+          if (r == (const Quantum *) NULL)
+            {
+              status=MagickFalse;
+              continue;
+            }
+          alpha=(double) (QuantumScale*GetPixelAlpha(image,r));
+          pixel+=(*k)*alpha*r[i];
+          gamma+=(*k)*alpha;
+          k++;
         }
+        gamma=MagickEpsilonReciprocal(gamma);
+        SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
+      }
+      p+=GetPixelChannels(image);
       q+=GetPixelChannels(blur_image);
     }
     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
@@ -2334,7 +2064,7 @@ MagickExport Image *MotionBlurImage(const Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_MotionBlurImage)
+        #pragma omp critical (MagickCore_MotionBlurImage)
 #endif
         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
@@ -2342,8 +2072,9 @@ MagickExport Image *MotionBlurImage(const Image *image,
       }
   }
   blur_view=DestroyCacheView(blur_view);
+  motion_view=DestroyCacheView(motion_view);
   image_view=DestroyCacheView(image_view);
-  kernel=(double *) RelinquishMagickMemory(kernel);
+  kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
   offset=(OffsetInfo *) RelinquishMagickMemory(offset);
   if (status == MagickFalse)
     blur_image=DestroyImage(blur_image);
@@ -2399,6 +2130,9 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
     sigma,
     threshold;
 
+  extern const char
+    DefaultTileFrame[];
+
   Image
     *images,
     *montage_image,
@@ -2459,10 +2193,11 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
       break;
     (void) SetImageProgressMonitor(thumbnail,(MagickProgressMonitor) NULL,
       (void *) NULL);
-    (void) SetImageProperty(thumbnail,"label",DefaultTileLabel);
+    (void) SetImageProperty(thumbnail,"label",DefaultTileLabel,exception);
     if (i == (NumberTiles/2))
       {
-        (void) QueryColorDatabase("#dfdfdf",&thumbnail->matte_color,exception);
+        (void) QueryColorCompliance("#dfdfdf",AllCompliance,
+          &thumbnail->matte_color,exception);
         AppendImageToList(&images,thumbnail);
         continue;
       }
@@ -2652,15 +2387,14 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
       }
       case SharpenPreview:
       {
-        preview_image=SharpenImage(thumbnail,radius,sigma,image->bias,
-          exception);
+        preview_image=SharpenImage(thumbnail,radius,sigma,exception);
         (void) FormatLocaleString(label,MaxTextExtent,"sharpen %gx%g",
           radius,sigma);
         break;
       }
       case BlurPreview:
       {
-        preview_image=BlurImage(thumbnail,radius,sigma,image->bias,exception);
+        preview_image=BlurImage(thumbnail,radius,sigma,exception);
         (void) FormatLocaleString(label,MaxTextExtent,"blur %gx%g",radius,
           sigma);
         break;
@@ -2670,10 +2404,10 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
         preview_image=CloneImage(thumbnail,0,0,MagickTrue,exception);
         if (preview_image == (Image *) NULL)
           break;
-        (void) BilevelImage(thumbnail,
-          (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
+        (void) BilevelImage(thumbnail,(double) (percentage*((double)
+          QuantumRange+1.0))/100.0,exception);
         (void) FormatLocaleString(label,MaxTextExtent,"threshold %g",
-          (double) (percentage*((MagickRealType) QuantumRange+1.0))/100.0);
+          (double) (percentage*((double) QuantumRange+1.0))/100.0);
         break;
       }
       case EdgeDetectPreview:
@@ -2684,7 +2418,8 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
       }
       case SpreadPreview:
       {
-        preview_image=SpreadImage(thumbnail,radius,exception);
+        preview_image=SpreadImage(thumbnail,radius,thumbnail->interpolate,
+          exception);
         (void) FormatLocaleString(label,MaxTextExtent,"spread %g",
           radius+0.5);
         break;
@@ -2730,7 +2465,7 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
         if (preview_image == (Image *) NULL)
           break;
         threshold+=0.4f;
-        (void) SegmentImage(preview_image,RGBColorspace,MagickFalse,threshold,
+        (void) SegmentImage(preview_image,sRGBColorspace,MagickFalse,threshold,
           threshold,exception);
         (void) FormatLocaleString(label,MaxTextExtent,"segment %gx%g",
           threshold,threshold);
@@ -2738,7 +2473,8 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
       }
       case SwirlPreview:
       {
-        preview_image=SwirlImage(thumbnail,degrees,exception);
+        preview_image=SwirlImage(thumbnail,degrees,image->interpolate,
+          exception);
         (void) FormatLocaleString(label,MaxTextExtent,"swirl %g",degrees);
         degrees+=45.0;
         break;
@@ -2746,14 +2482,16 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
       case ImplodePreview:
       {
         degrees+=0.1f;
-        preview_image=ImplodeImage(thumbnail,degrees,exception);
+        preview_image=ImplodeImage(thumbnail,degrees,image->interpolate,
+          exception);
         (void) FormatLocaleString(label,MaxTextExtent,"implode %g",degrees);
         break;
       }
       case WavePreview:
       {
         degrees+=5.0f;
-        preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,exception);
+        preview_image=WaveImage(thumbnail,0.5*degrees,2.0*degrees,
+          image->interpolate,exception);
         (void) FormatLocaleString(label,MaxTextExtent,"wave %gx%g",
           0.5*degrees,2.0*degrees);
         break;
@@ -2769,7 +2507,7 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
       case CharcoalDrawingPreview:
       {
         preview_image=CharcoalImage(thumbnail,(double) radius,(double) sigma,
-          image->bias,exception);
+          exception);
         (void) FormatLocaleString(label,MaxTextExtent,"charcoal %gx%g",
           radius,sigma);
         break;
@@ -2834,7 +2572,7 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
     if (preview_image == (Image *) NULL)
       break;
     (void) DeleteImageProperty(preview_image,"label");
-    (void) SetImageProperty(preview_image,"label",label);
+    (void) SetImageProperty(preview_image,"label",label,exception);
     AppendImageToList(&images,preview_image);
     proceed=SetImageProgress(image,PreviewImageTag,(MagickOffsetType) i,
       NumberTiles);
@@ -2901,15 +2639,18 @@ MagickExport Image *PreviewImage(const Image *image,const PreviewType preview,
 %
 %    o angle: the angle of the radial blur.
 %
+%    o blur: the blur.
+%
 %    o exception: return any errors or warnings in this structure.
 %
 */
-MagickExport Image *RadialBlurImage(const Image *image,
-  const double angle,ExceptionInfo *exception)
+MagickExport Image *RadialBlurImage(const Image *image,const double angle,
+  ExceptionInfo *exception)
 {
   CacheView
     *blur_view,
-    *image_view;
+    *image_view,
+    *radial_view;
 
   Image
     *blur_image;
@@ -2920,10 +2661,7 @@ MagickExport Image *RadialBlurImage(const Image *image,
   MagickOffsetType
     progress;
 
-  PixelInfo
-    bias;
-
-  MagickRealType
+  double
     blur_radius,
     *cos_theta,
     offset,
@@ -2951,7 +2689,7 @@ MagickExport Image *RadialBlurImage(const Image *image,
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
-  blur_image=CloneImage(image,0,0,MagickTrue,exception);
+  blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
   if (blur_image == (Image *) NULL)
     return((Image *) NULL);
   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
@@ -2963,18 +2701,18 @@ MagickExport Image *RadialBlurImage(const Image *image,
   blur_center.y=(double) image->rows/2.0;
   blur_radius=hypot(blur_center.x,blur_center.y);
   n=(size_t) fabs(4.0*DegreesToRadians(angle)*sqrt((double) blur_radius)+2UL);
-  theta=DegreesToRadians(angle)/(MagickRealType) (n-1);
-  cos_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
+  theta=DegreesToRadians(angle)/(double) (n-1);
+  cos_theta=(double *) AcquireQuantumMemory((size_t) n,
     sizeof(*cos_theta));
-  sin_theta=(MagickRealType *) AcquireQuantumMemory((size_t) n,
+  sin_theta=(double *) AcquireQuantumMemory((size_t) n,
     sizeof(*sin_theta));
-  if ((cos_theta == (MagickRealType *) NULL) ||
-      (sin_theta == (MagickRealType *) NULL))
+  if ((cos_theta == (double *) NULL) ||
+      (sin_theta == (double *) NULL))
     {
       blur_image=DestroyImage(blur_image);
       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
     }
-  offset=theta*(MagickRealType) (n-1)/2.0;
+  offset=theta*(double) (n-1)/2.0;
   for (i=0; i < (ssize_t) n; i++)
   {
     cos_theta[i]=cos((double) (theta*i-offset));
@@ -2985,14 +2723,18 @@ MagickExport Image *RadialBlurImage(const Image *image,
   */
   status=MagickTrue;
   progress=0;
-  GetPixelInfo(image,&bias);
-  image_view=AcquireCacheView(image);
-  blur_view=AcquireCacheView(blur_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  radial_view=AcquireVirtualCacheView(image,exception);
+  blur_view=AcquireAuthenticCacheView(blur_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
-  for (y=0; y < (ssize_t) blur_image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
+    register const Quantum
+      *restrict p;
+
     register Quantum
       *restrict q;
 
@@ -3001,25 +2743,19 @@ MagickExport Image *RadialBlurImage(const Image *image,
 
     if (status == MagickFalse)
       continue;
-    q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+    q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
       exception);
-    if (q == (Quantum *) NULL)
+    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
       {
         status=MagickFalse;
         continue;
       }
-    for (x=0; x < (ssize_t) blur_image->columns; x++)
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
-      PixelInfo
-        qixel;
-
-      MagickRealType
-        normalize,
+      double
         radius;
 
-      PixelPacket
-        pixel;
-
       PointInfo
         center;
 
@@ -3043,87 +2779,77 @@ MagickExport Image *RadialBlurImage(const Image *image,
             if (step >= n)
               step=n-1;
         }
-      normalize=0.0;
-      qixel=bias;
-      if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) == 0) || (image->matte == MagickFalse))
-        {
-          for (i=0; i < (ssize_t) n; i+=(ssize_t) step)
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        double
+          gamma,
+          pixel;
+
+        PixelChannel
+          channel;
+
+        PixelTrait
+          blur_traits,
+          traits;
+
+        register const Quantum
+          *restrict r;
+
+        register ssize_t
+          j;
+
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        blur_traits=GetPixelChannelTraits(blur_image,channel);
+        if ((traits == UndefinedPixelTrait) ||
+            (blur_traits == UndefinedPixelTrait))
+          continue;
+        if (((blur_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
-            (void) GetOneCacheViewVirtualPixel(image_view,(ssize_t)
-              (blur_center.x+center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),
-              (ssize_t) (blur_center.y+center.x*sin_theta[i]+center.y*
-              cos_theta[i]+0.5),&pixel,exception);
-            qixel.red+=pixel.red;
-            qixel.green+=pixel.green;
-            qixel.blue+=pixel.blue;
-            if (image->colorspace == CMYKColorspace)
-              qixel.black+=pixel.black;
-            qixel.alpha+=pixel.alpha;
-            normalize+=1.0;
+            SetPixelChannel(blur_image,channel,p[i],q);
+            continue;
           }
-          normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
-            normalize);
-          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelRed(blur_image,
-              ClampToQuantum(normalize*qixel.red),q);
-          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelGreen(blur_image,
-              ClampToQuantum(normalize*qixel.green),q);
-          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelBlue(blur_image,
-              ClampToQuantum(normalize*qixel.blue),q);
-          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-              (image->colorspace == CMYKColorspace))
-            SetPixelBlack(blur_image,
-              ClampToQuantum(normalize*qixel.black),q);
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelAlpha(blur_image,
-              ClampToQuantum(normalize*qixel.alpha),q);
-        }
-      else
-        {
-          MagickRealType
-            alpha,
-            gamma;
-
-          alpha=1.0;
-          gamma=0.0;
-          for (i=0; i < (ssize_t) n; i+=(ssize_t) step)
+        gamma=0.0;
+        pixel=0.0;
+        if ((blur_traits & BlendPixelTrait) == 0)
           {
-            (void) GetOneCacheViewVirtualPixel(image_view,(ssize_t)
-              (blur_center.x+center.x*cos_theta[i]-center.y*sin_theta[i]+0.5),
-              (ssize_t) (blur_center.y+center.x*sin_theta[i]+center.y*
-              cos_theta[i]+0.5),&pixel,exception);
-            alpha=(MagickRealType) (QuantumScale*pixel.alpha);
-            qixel.red+=alpha*pixel.red;
-            qixel.green+=alpha*pixel.green;
-            qixel.blue+=alpha*pixel.blue;
-            qixel.alpha+=pixel.alpha;
-            if (image->colorspace == CMYKColorspace)
-              qixel.black+=alpha*pixel.black;
-            gamma+=alpha;
-            normalize+=1.0;
+            for (j=0; j < (ssize_t) n; j+=(ssize_t) step)
+            {
+              r=GetCacheViewVirtualPixels(radial_view, (ssize_t) (blur_center.x+
+                center.x*cos_theta[j]-center.y*sin_theta[j]+0.5),(ssize_t)
+                (blur_center.y+center.x*sin_theta[j]+center.y*cos_theta[j]+0.5),
+                1,1,exception);
+              if (r == (const Quantum *) NULL)
+                {
+                  status=MagickFalse;
+                  continue;
+                }
+              pixel+=r[i];
+              gamma++;
+            }
+            gamma=MagickEpsilonReciprocal(gamma);
+            SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
+            continue;
           }
-          gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-          normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 :
-            normalize);
-          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelRed(blur_image,
-              ClampToQuantum(gamma*qixel.red),q);
-          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelGreen(blur_image,
-              ClampToQuantum(gamma*qixel.green),q);
-          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelBlue(blur_image,
-              ClampToQuantum(gamma*qixel.blue),q);
-          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-              (image->colorspace == CMYKColorspace))
-            SetPixelBlack(blur_image,
-              ClampToQuantum(gamma*qixel.black),q);
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
-            SetPixelAlpha(blur_image,
-              ClampToQuantum(normalize*qixel.alpha),q);
+        for (j=0; j < (ssize_t) n; j+=(ssize_t) step)
+        {
+          r=GetCacheViewVirtualPixels(radial_view, (ssize_t) (blur_center.x+
+            center.x*cos_theta[j]-center.y*sin_theta[j]+0.5),(ssize_t)
+            (blur_center.y+center.x*sin_theta[j]+center.y*cos_theta[j]+0.5),
+            1,1,exception);
+          if (r == (const Quantum *) NULL)
+            {
+              status=MagickFalse;
+              continue;
+            }
+          pixel+=GetPixelAlpha(image,r)*r[i];
+          gamma+=GetPixelAlpha(image,r);
         }
+        gamma=MagickEpsilonReciprocal(gamma);
+        SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
+      }
+      p+=GetPixelChannels(image);
       q+=GetPixelChannels(blur_image);
     }
     if (SyncCacheViewAuthenticPixels(blur_view,exception) == MagickFalse)
@@ -3134,7 +2860,7 @@ MagickExport Image *RadialBlurImage(const Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_RadialBlurImage)
+        #pragma omp critical (MagickCore_RadialBlurImage)
 #endif
         proceed=SetImageProgress(image,BlurImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
@@ -3142,9 +2868,10 @@ MagickExport Image *RadialBlurImage(const Image *image,
       }
   }
   blur_view=DestroyCacheView(blur_view);
+  radial_view=DestroyCacheView(radial_view);
   image_view=DestroyCacheView(image_view);
-  cos_theta=(MagickRealType *) RelinquishMagickMemory(cos_theta);
-  sin_theta=(MagickRealType *) RelinquishMagickMemory(sin_theta);
+  cos_theta=(double *) RelinquishMagickMemory(cos_theta);
+  sin_theta=(double *) RelinquishMagickMemory(sin_theta);
   if (status == MagickFalse)
     blur_image=DestroyImage(blur_image);
   return(blur_image);
@@ -3185,21 +2912,19 @@ MagickExport Image *RadialBlurImage(const Image *image,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-MagickExport Image *SelectiveBlurImage(const Image *image,
-  const double radius,const double sigma,const double threshold,
-  ExceptionInfo *exception)
+MagickExport Image *SelectiveBlurImage(const Image *image,const double radius,
+  const double sigma,const double threshold,ExceptionInfo *exception)
 {
 #define SelectiveBlurImageTag  "SelectiveBlur/Image"
 
   CacheView
     *blur_view,
-    *image_view;
-
-  double
-    *kernel;
+    *image_view,
+    *luminance_view;
 
   Image
-    *blur_image;
+    *blur_image,
+    *luminance_image;
 
   MagickBooleanType
     status;
@@ -3207,8 +2932,8 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
   MagickOffsetType
     progress;
 
-  PixelInfo
-    bias;
+  MagickRealType
+    *kernel;
 
   register ssize_t
     i;
@@ -3217,6 +2942,7 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
     width;
 
   ssize_t
+    center,
     j,
     u,
     v,
@@ -3232,15 +2958,16 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
   width=GetOptimalKernelWidth1D(radius,sigma);
-  kernel=(double *) AcquireQuantumMemory((size_t) width,width*sizeof(*kernel));
-  if (kernel == (double *) NULL)
+  kernel=(MagickRealType *) MagickAssumeAligned(AcquireAlignedMemory((size_t)
+    width,width*sizeof(*kernel)));
+  if (kernel == (MagickRealType *) NULL)
     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
   j=(ssize_t) width/2;
   i=0;
   for (v=(-j); v <= j; v++)
   {
     for (u=(-j); u <= j; u++)
-      kernel[i++]=(double) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
+      kernel[i++]=(MagickRealType) (exp(-((double) u*u+v*v)/(2.0*MagickSigma*
         MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
   }
   if (image->debug != MagickFalse)
@@ -3249,7 +2976,7 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
         format[MaxTextExtent],
         *message;
 
-      register const double
+      register const MagickRealType
         *k;
 
       ssize_t
@@ -3268,19 +2995,35 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
         (void) ConcatenateString(&message,format);
         for (u=0; u < (ssize_t) width; u++)
         {
-          (void) FormatLocaleString(format,MaxTextExtent,"%+f ",*k++);
+          (void) FormatLocaleString(format,MaxTextExtent,"%+f ",(double) *k++);
           (void) ConcatenateString(&message,format);
         }
         (void) LogMagickEvent(TransformEvent,GetMagickModule(),"%s",message);
       }
       message=DestroyString(message);
     }
-  blur_image=CloneImage(image,0,0,MagickTrue,exception);
+  blur_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
   if (blur_image == (Image *) NULL)
     return((Image *) NULL);
   if (SetImageStorageClass(blur_image,DirectClass,exception) == MagickFalse)
     {
       blur_image=DestroyImage(blur_image);
+      kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
+      return((Image *) NULL);
+    }
+  luminance_image=CloneImage(image,0,0,MagickTrue,exception);
+  if (luminance_image == (Image *) NULL)
+    {
+      blur_image=DestroyImage(blur_image);
+      kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
+      return((Image *) NULL);
+    }
+  status=TransformImageColorspace(luminance_image,GRAYColorspace,exception);
+  if (status == MagickFalse)
+    {
+      luminance_image=DestroyImage(luminance_image);
+      blur_image=DestroyImage(blur_image);
+      kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
       return((Image *) NULL);
     }
   /*
@@ -3288,12 +3031,14 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
   */
   status=MagickTrue;
   progress=0;
-  GetPixelInfo(image,&bias);
-  SetPixelInfoBias(image,&bias);
-  image_view=AcquireCacheView(image);
-  blur_view=AcquireCacheView(blur_image);
+  center=(ssize_t) (GetPixelChannels(image)*(image->columns+width)*(width/2L)+
+    GetPixelChannels(image)*(width/2L));
+  image_view=AcquireVirtualCacheView(image,exception);
+  luminance_view=AcquireVirtualCacheView(luminance_image,exception);
+  blur_view=AcquireAuthenticCacheView(blur_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -3303,10 +3048,8 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
     MagickBooleanType
       sync;
 
-    MagickRealType
-      gamma;
-
     register const Quantum
+      *restrict l,
       *restrict p;
 
     register Quantum
@@ -3319,7 +3062,9 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
       continue;
     p=GetCacheViewVirtualPixels(image_view,-((ssize_t) width/2L),y-(ssize_t)
       (width/2L),image->columns+width,width,exception);
-    q=GetCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
+    l=GetCacheViewVirtualPixels(luminance_view,-((ssize_t) width/2L),y-(ssize_t)
+      (width/2L),luminance_image->columns+width,width,exception);
+    q=QueueCacheViewAuthenticPixels(blur_view,0,y,blur_image->columns,1,
       exception);
     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
       {
@@ -3328,166 +3073,117 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
-      PixelInfo
-        pixel;
-
-      register const double
-        *restrict k;
+      double
+        intensity;
 
       register ssize_t
-        u;
+        i;
 
-      ssize_t
-        j,
-        v;
+      intensity=GetPixelIntensity(image,p+center);
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
+      {
+        double
+          alpha,
+          gamma,
+          pixel;
 
-      pixel=bias;
-      k=kernel;
-      gamma=0.0;
-      j=0;
-      if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) == 0) || (image->matte == MagickFalse))
-        {
-          for (v=0; v < (ssize_t) width; v++)
+        PixelChannel
+          channel;
+
+        PixelTrait
+          blur_traits,
+          traits;
+
+        register const MagickRealType
+          *restrict k;
+
+        register const Quantum
+          *restrict luminance_pixels,
+          *restrict pixels;
+
+        register ssize_t
+          u;
+
+        ssize_t
+          v;
+
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        blur_traits=GetPixelChannelTraits(blur_image,channel);
+        if ((traits == UndefinedPixelTrait) ||
+            (blur_traits == UndefinedPixelTrait))
+          continue;
+        if (((blur_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
-            for (u=0; u < (ssize_t) width; u++)
-            {
-              contrast=GetPixelIntensity(image,p+(u+j)*GetPixelChannels(image))-
-                (double) GetPixelIntensity(blur_image,q);
-              if (fabs(contrast) < threshold)
-                {
-                  pixel.red+=(*k)*
-                    GetPixelRed(image,p+(u+j)*GetPixelChannels(image));
-                  pixel.green+=(*k)*
-                    GetPixelGreen(image,p+(u+j)*GetPixelChannels(image));
-                  pixel.blue+=(*k)*
-                    GetPixelBlue(image,p+(u+j)*GetPixelChannels(image));
-                  if (image->colorspace == CMYKColorspace)
-                    pixel.black+=(*k)*
-                      GetPixelBlack(image,p+(u+j)*GetPixelChannels(image));
-                  gamma+=(*k);
-                  k++;
-                }
-            }
-            j+=(ssize_t) (image->columns+width);
+            SetPixelChannel(blur_image,channel,p[center+i],q);
+            continue;
           }
-          if (gamma != 0.0)
-            {
-              gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-              if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-                SetPixelRed(blur_image,ClampToQuantum(gamma*pixel.red),q);
-              if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-                SetPixelGreen(blur_image,ClampToQuantum(gamma*pixel.green),q);
-              if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-                SetPixelBlue(blur_image,ClampToQuantum(gamma*pixel.blue),q);
-              if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-                  (image->colorspace == CMYKColorspace))
-                SetPixelBlack(blur_image,ClampToQuantum(gamma*pixel.black),q);
-            }
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
+        k=kernel;
+        pixel=0.0;
+        pixels=p;
+        luminance_pixels=l;
+        gamma=0.0;
+        if ((blur_traits & BlendPixelTrait) == 0)
+          {
+            for (v=0; v < (ssize_t) width; v++)
             {
-              gamma=0.0;
-              j=0;
-              for (v=0; v < (ssize_t) width; v++)
+              for (u=0; u < (ssize_t) width; u++)
               {
-                for (u=0; u < (ssize_t) width; u++)
-                {
-                  contrast=GetPixelIntensity(image,p+(u+j)*
-                    GetPixelChannels(image))-(double)
-                    GetPixelIntensity(blur_image,q);
-                  if (fabs(contrast) < threshold)
-                    {
-                      pixel.alpha+=(*k)*
-                        GetPixelAlpha(image,p+(u+j)*GetPixelChannels(image));
-                      gamma+=(*k);
-                      k++;
-                    }
-                }
-                j+=(ssize_t) (image->columns+width);
+                contrast=GetPixelIntensity(luminance_image,luminance_pixels)-
+                  intensity;
+                if (fabs(contrast) < threshold)
+                  {
+                    pixel+=(*k)*pixels[i];
+                    gamma+=(*k);
+                  }
+                k++;
+                pixels+=GetPixelChannels(image);
+                luminance_pixels+=GetPixelChannels(luminance_image);
               }
-              if (gamma != 0.0)
-                {
-                  gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
-                    gamma);
-                  SetPixelAlpha(blur_image,ClampToQuantum(gamma*pixel.alpha),q);
-                }
+              pixels+=image->columns*GetPixelChannels(image);
+              luminance_pixels+=luminance_image->columns*
+                GetPixelChannels(luminance_image);
             }
-        }
-      else
+            if (fabs((double) gamma) < MagickEpsilon)
+              {
+                SetPixelChannel(blur_image,channel,p[center+i],q);
+                continue;
+              }
+            gamma=MagickEpsilonReciprocal(gamma);
+            SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
+            continue;
+          }
+        for (v=0; v < (ssize_t) width; v++)
         {
-          MagickRealType
-            alpha;
-
-          for (v=0; v < (ssize_t) width; v++)
+          for (u=0; u < (ssize_t) width; u++)
           {
-            for (u=0; u < (ssize_t) width; u++)
-            {
-              contrast=GetPixelIntensity(image,p+(u+j)*
-                GetPixelChannels(image))-(double)
-                GetPixelIntensity(blur_image,q);
-              if (fabs(contrast) < threshold)
-                {
-                  alpha=(MagickRealType) (QuantumScale*
-                    GetPixelAlpha(image,p+(u+j)*GetPixelChannels(image)));
-                  pixel.red+=(*k)*alpha*
-                    GetPixelRed(image,p+(u+j)*GetPixelChannels(image));
-                  pixel.green+=(*k)*alpha*GetPixelGreen(image,p+(u+j)*
-                    GetPixelChannels(image));
-                  pixel.blue+=(*k)*alpha*GetPixelBlue(image,p+(u+j)*
-                    GetPixelChannels(image));
-                  pixel.alpha+=(*k)*GetPixelAlpha(image,p+(u+j)*
-                    GetPixelChannels(image));
-                  if (image->colorspace == CMYKColorspace)
-                    pixel.black+=(*k)*GetPixelBlack(image,p+(u+j)*
-                      GetPixelChannels(image));
-                  gamma+=(*k)*alpha;
-                  k++;
-                }
-            }
-            j+=(ssize_t) (image->columns+width);
-          }
-          if (gamma != 0.0)
-            {
-              gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-              if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-                SetPixelRed(blur_image,ClampToQuantum(gamma*pixel.red),q);
-              if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-                SetPixelGreen(blur_image,ClampToQuantum(gamma*pixel.green),q);
-              if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-                SetPixelBlue(blur_image,ClampToQuantum(gamma*pixel.blue),q);
-              if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-                  (image->colorspace == CMYKColorspace))
-                SetPixelBlack(blur_image,ClampToQuantum(gamma*pixel.black),q);
-            }
-          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0)
-            {
-              gamma=0.0;
-              j=0;
-              for (v=0; v < (ssize_t) width; v++)
+            contrast=GetPixelIntensity(image,pixels)-intensity;
+            if (fabs(contrast) < threshold)
               {
-                for (u=0; u < (ssize_t) width; u++)
-                {
-                  contrast=GetPixelIntensity(image,p+(u+j)*
-                    GetPixelChannels(image))-(double)
-                    GetPixelIntensity(blur_image,q);
-                  if (fabs(contrast) < threshold)
-                    {
-                      pixel.alpha+=(*k)*
-                        GetPixelAlpha(image,p+(u+j)*GetPixelChannels(image));
-                      gamma+=(*k);
-                      k++;
-                    }
-                }
-                j+=(ssize_t) (image->columns+width);
+                alpha=(double) (QuantumScale*
+                  GetPixelAlpha(image,pixels));
+                pixel+=(*k)*alpha*pixels[i];
+                gamma+=(*k)*alpha;
               }
-              if (gamma != 0.0)
-                {
-                  gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 :
-                    gamma);
-                  SetPixelAlpha(blur_image,ClampToQuantum(pixel.alpha),q);
-                }
-            }
+            k++;
+            pixels+=GetPixelChannels(image);
+            luminance_pixels+=GetPixelChannels(luminance_image);
+          }
+          pixels+=image->columns*GetPixelChannels(image);
+          luminance_pixels+=luminance_image->columns*
+            GetPixelChannels(luminance_image);
         }
+        if (fabs((double) gamma) < MagickEpsilon)
+          {
+            SetPixelChannel(blur_image,channel,p[center+i],q);
+            continue;
+          }
+        gamma=MagickEpsilonReciprocal(gamma);
+        SetPixelChannel(blur_image,channel,ClampToQuantum(gamma*pixel),q);
+      }
       p+=GetPixelChannels(image);
+      l+=GetPixelChannels(luminance_image);
       q+=GetPixelChannels(blur_image);
     }
     sync=SyncCacheViewAuthenticPixels(blur_view,exception);
@@ -3499,7 +3195,7 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_SelectiveBlurImage)
+        #pragma omp critical (MagickCore_SelectiveBlurImage)
 #endif
         proceed=SetImageProgress(image,SelectiveBlurImageTag,progress++,
           image->rows);
@@ -3510,7 +3206,8 @@ MagickExport Image *SelectiveBlurImage(const Image *image,
   blur_image->type=image->type;
   blur_view=DestroyCacheView(blur_view);
   image_view=DestroyCacheView(image_view);
-  kernel=(double *) RelinquishMagickMemory(kernel);
+  luminance_image=DestroyImage(luminance_image);
+  kernel=(MagickRealType *) RelinquishAlignedMemory(kernel);
   if (status == MagickFalse)
     blur_image=DestroyImage(blur_image);
   return(blur_image);
@@ -3602,14 +3299,15 @@ MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
-  shade_view=AcquireCacheView(shade_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  shade_view=AcquireAuthenticCacheView(shade_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
-    MagickRealType
+    double
       distance,
       normal_distance,
       shade;
@@ -3618,10 +3316,10 @@ MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
       normal;
 
     register const Quantum
+      *restrict center,
       *restrict p,
-      *restrict s0,
-      *restrict s1,
-      *restrict s2;
+      *restrict post,
+      *restrict pre;
 
     register Quantum
       *restrict q;
@@ -3643,26 +3341,28 @@ MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
       Shade this row of pixels.
     */
     normal.z=2.0*(double) QuantumRange;  /* constant Z of surface normal */
-    s0=p+GetPixelChannels(image);
-    s1=s0+(image->columns+2)*GetPixelChannels(image);
-    s2=s1+(image->columns+2)*GetPixelChannels(image);
+    pre=p+GetPixelChannels(image);
+    center=pre+(image->columns+2)*GetPixelChannels(image);
+    post=center+(image->columns+2)*GetPixelChannels(image);
     for (x=0; x < (ssize_t) image->columns; x++)
     {
+      register ssize_t
+        i;
+
       /*
         Determine the surface normal and compute shading.
       */
-      normal.x=(double) (GetPixelIntensity(image,s0-GetPixelChannels(image))+
-        GetPixelIntensity(image,s1-GetPixelChannels(image))+
-        GetPixelIntensity(image,s2-GetPixelChannels(image))-
-        GetPixelIntensity(image,s0+GetPixelChannels(image))-
-        GetPixelIntensity(image,s1+GetPixelChannels(image))-
-        GetPixelIntensity(image,s2+GetPixelChannels(image)));
-      normal.y=(double) (GetPixelIntensity(image,s2-GetPixelChannels(image))+
-        GetPixelIntensity(image,s2)+
-        GetPixelIntensity(image,s2+GetPixelChannels(image))-
-        GetPixelIntensity(image,s0-GetPixelChannels(image))-
-        GetPixelIntensity(image,s0)-
-        GetPixelIntensity(image,s0+GetPixelChannels(image)));
+      normal.x=(double) (GetPixelIntensity(image,pre-GetPixelChannels(image))+
+        GetPixelIntensity(image,center-GetPixelChannels(image))+
+        GetPixelIntensity(image,post-GetPixelChannels(image))-
+        GetPixelIntensity(image,pre+GetPixelChannels(image))-
+        GetPixelIntensity(image,center+GetPixelChannels(image))-
+        GetPixelIntensity(image,post+GetPixelChannels(image)));
+      normal.y=(double) (GetPixelIntensity(image,post-GetPixelChannels(image))+
+        GetPixelIntensity(image,post)+GetPixelIntensity(image,post+
+        GetPixelChannels(image))-GetPixelIntensity(image,pre-
+        GetPixelChannels(image))-GetPixelIntensity(image,pre)-
+        GetPixelIntensity(image,pre+GetPixelChannels(image)));
       if ((normal.x == 0.0) && (normal.y == 0.0))
         shade=light.z;
       else
@@ -3677,947 +3377,207 @@ MagickExport Image *ShadeImage(const Image *image,const MagickBooleanType gray,
                 shade=distance/sqrt((double) normal_distance);
             }
         }
-      if (gray != MagickFalse)
-        {
-          SetPixelRed(shade_image,ClampToQuantum(shade),q);
-          SetPixelGreen(shade_image,ClampToQuantum(shade),q);
-          SetPixelBlue(shade_image,ClampToQuantum(shade),q);
-        }
-      else
-        {
-          SetPixelRed(shade_image,ClampToQuantum(QuantumScale*shade*
-            GetPixelRed(image,s1)),q);
-          SetPixelGreen(shade_image,ClampToQuantum(QuantumScale*shade*
-            GetPixelGreen(image,s1)),q);
-          SetPixelBlue(shade_image,ClampToQuantum(QuantumScale*shade*
-            GetPixelBlue(image,s1)),q);
-        }
-      SetPixelAlpha(shade_image,GetPixelAlpha(image,s1),q);
-      s0+=GetPixelChannels(image);
-      s1+=GetPixelChannels(image);
-      s2+=GetPixelChannels(image);
-      q+=GetPixelChannels(shade_image);
-    }
-    if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
-      status=MagickFalse;
-    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickBooleanType
-          proceed;
-
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_ShadeImage)
-#endif
-        proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
-        if (proceed == MagickFalse)
-          status=MagickFalse;
-      }
-  }
-  shade_view=DestroyCacheView(shade_view);
-  image_view=DestroyCacheView(image_view);
-  if (status == MagickFalse)
-    shade_image=DestroyImage(shade_image);
-  return(shade_image);
-}
-\f
-/*
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%     S h a r p e n I m a g e                                                 %
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
-%  operator of the given radius and standard deviation (sigma).  For
-%  reasonable results, radius should be larger than sigma.  Use a radius of 0
-%  and SharpenImage() selects a suitable radius for you.
-%
-%  Using a separable kernel would be faster, but the negative weights cancel
-%  out on the corners of the kernel producing often undesirable ringing in the
-%  filtered result; this can be avoided by using a 2D gaussian shaped image
-%  sharpening kernel instead.
-%
-%  The format of the SharpenImage method is:
-%
-%    Image *SharpenImage(const Image *image,const double radius,
-%      const double sigma,const double bias,ExceptionInfo *exception)
-%
-%  A description of each parameter follows:
-%
-%    o image: the image.
-%
-%    o radius: the radius of the Gaussian, in pixels, not counting the center
-%      pixel.
-%
-%    o sigma: the standard deviation of the Laplacian, in pixels.
-%
-%    o bias: bias.
-%
-%    o exception: return any errors or warnings in this structure.
-%
-*/
-MagickExport Image *SharpenImage(const Image *image,const double radius,
-  const double sigma,const double bias,ExceptionInfo *exception)
-{
-  double
-    normalize;
-
-  Image
-    *sharp_image;
-
-  KernelInfo
-    *kernel_info;
-
-  register ssize_t
-    i;
-
-  size_t
-    width;
-
-  ssize_t
-    j,
-    u,
-    v;
+        PixelChannel
+          channel;
 
-  assert(image != (const Image *) NULL);
-  assert(image->signature == MagickSignature);
-  if (image->debug != MagickFalse)
-    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  assert(exception != (ExceptionInfo *) NULL);
-  assert(exception->signature == MagickSignature);
-  width=GetOptimalKernelWidth2D(radius,sigma);
-  kernel_info=AcquireKernelInfo((const char *) NULL);
-  if (kernel_info == (KernelInfo *) NULL)
-    ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-  (void) ResetMagickMemory(kernel_info,0,sizeof(*kernel_info));
-  kernel_info->width=width;
-  kernel_info->height=width;
-  kernel_info->bias=bias;
-  kernel_info->signature=MagickSignature;
-  kernel_info->values=(double *) AcquireAlignedMemory(kernel_info->width,
-    kernel_info->width*sizeof(*kernel_info->values));
-  if (kernel_info->values == (double *) NULL)
-    {
-      kernel_info=DestroyKernelInfo(kernel_info);
-      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-    }
-  normalize=0.0;
-  j=(ssize_t) kernel_info->width/2;
-  i=0;
-  for (v=(-j); v <= j; v++)
-  {
-    for (u=(-j); u <= j; u++)
-    {
-      kernel_info->values[i]=(double) (-exp(-((double) u*u+v*v)/(2.0*
-        MagickSigma*MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
-      normalize+=kernel_info->values[i];
-      i++;
-    }
-  }
-  kernel_info->values[i/2]=(double) ((-2.0)*normalize);
-  sharp_image=ConvolveImage(image,kernel_info,exception);
-  kernel_info=DestroyKernelInfo(kernel_info);
-  return(sharp_image);
-}
-\f
-/*
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%     S p r e a d I m a g e                                                   %
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%  SpreadImage() is a special effects method that randomly displaces each
-%  pixel in a block defined by the radius parameter.
-%
-%  The format of the SpreadImage method is:
-%
-%      Image *SpreadImage(const Image *image,const double radius,
-%        ExceptionInfo *exception)
-%
-%  A description of each parameter follows:
-%
-%    o image: the image.
-%
-%    o radius:  Choose a random pixel in a neighborhood of this extent.
-%
-%    o exception: return any errors or warnings in this structure.
-%
-*/
-MagickExport Image *SpreadImage(const Image *image,const double radius,
-  ExceptionInfo *exception)
-{
-#define SpreadImageTag  "Spread/Image"
-
-  CacheView
-    *image_view,
-    *spread_view;
-
-  Image
-    *spread_image;
-
-  MagickBooleanType
-    status;
-
-  MagickOffsetType
-    progress;
-
-  PixelInfo
-    bias;
-
-  RandomInfo
-    **restrict random_info;
-
-  size_t
-    width;
-
-  ssize_t
-    y;
-
-  /*
-    Initialize spread image attributes.
-  */
-  assert(image != (Image *) NULL);
-  assert(image->signature == MagickSignature);
-  if (image->debug != MagickFalse)
-    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  assert(exception != (ExceptionInfo *) NULL);
-  assert(exception->signature == MagickSignature);
-  spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
-    exception);
-  if (spread_image == (Image *) NULL)
-    return((Image *) NULL);
-  if (SetImageStorageClass(spread_image,DirectClass,exception) == MagickFalse)
-    {
-      spread_image=DestroyImage(spread_image);
-      return((Image *) NULL);
-    }
-  /*
-    Spread image.
-  */
-  status=MagickTrue;
-  progress=0;
-  GetPixelInfo(spread_image,&bias);
-  width=GetOptimalKernelWidth1D(radius,0.5);
-  random_info=AcquireRandomInfoThreadSet();
-  image_view=AcquireCacheView(image);
-  spread_view=AcquireCacheView(spread_image);
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status) omp_throttle(1)
-#endif
-  for (y=0; y < (ssize_t) spread_image->rows; y++)
-  {
-    const int
-      id = GetOpenMPThreadId();
-
-    PixelInfo
-      pixel;
-
-    register Quantum
-      *restrict q;
-
-    register ssize_t
-      x;
-
-    if (status == MagickFalse)
-      continue;
-    q=QueueCacheViewAuthenticPixels(spread_view,0,y,spread_image->columns,1,
-      exception);
-    if (q == (Quantum *) NULL)
-      {
-        status=MagickFalse;
-        continue;
-      }
-    pixel=bias;
-    for (x=0; x < (ssize_t) spread_image->columns; x++)
-    {
-      (void) InterpolatePixelInfo(image,image_view,
-        UndefinedInterpolatePixel,(double) x+width*(GetPseudoRandomValue(
-        random_info[id])-0.5),(double) y+width*(GetPseudoRandomValue(
-        random_info[id])-0.5),&pixel,exception);
-      SetPixelPixelInfo(spread_image,&pixel,q);
-      q+=GetPixelChannels(spread_image);
-    }
-    if (SyncCacheViewAuthenticPixels(spread_view,exception) == MagickFalse)
-      status=MagickFalse;
-    if (image->progress_monitor != (MagickProgressMonitor) NULL)
-      {
-        MagickBooleanType
-          proceed;
-
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_SpreadImage)
-#endif
-        proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
-        if (proceed == MagickFalse)
-          status=MagickFalse;
-      }
-  }
-  spread_view=DestroyCacheView(spread_view);
-  image_view=DestroyCacheView(image_view);
-  random_info=DestroyRandomInfoThreadSet(random_info);
-  return(spread_image);
-}
-\f
-/*
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%     S t a t i s t i c I m a g e                                             %
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%  StatisticImage() makes each pixel the min / max / median / mode / etc. of
-%  the neighborhood of the specified width and height.
-%
-%  The format of the StatisticImage method is:
-%
-%      Image *StatisticImage(const Image *image,const StatisticType type,
-%        const size_t width,const size_t height,ExceptionInfo *exception)
-%
-%  A description of each parameter follows:
-%
-%    o image: the image.
-%
-%    o type: the statistic type (median, mode, etc.).
-%
-%    o width: the width of the pixel neighborhood.
-%
-%    o height: the height of the pixel neighborhood.
-%
-%    o exception: return any errors or warnings in this structure.
-%
-*/
-
-#define ListChannels  5
-
-typedef struct _ListNode
-{
-  size_t
-    next[9],
-    count,
-    signature;
-} ListNode;
-
-typedef struct _SkipList
-{
-  ssize_t
-    level;
-
-  ListNode
-    *nodes;
-} SkipList;
-
-typedef struct _PixelList
-{
-  size_t
-    length,
-    seed,
-    signature;
-
-  SkipList
-    lists[ListChannels];
-} PixelList;
-
-static PixelList *DestroyPixelList(PixelList *pixel_list)
-{
-  register ssize_t
-    i;
-
-  if (pixel_list == (PixelList *) NULL)
-    return((PixelList *) NULL);
-  for (i=0; i < ListChannels; i++)
-    if (pixel_list->lists[i].nodes != (ListNode *) NULL)
-      pixel_list->lists[i].nodes=(ListNode *) RelinquishMagickMemory(
-        pixel_list->lists[i].nodes);
-  pixel_list=(PixelList *) RelinquishMagickMemory(pixel_list);
-  return(pixel_list);
-}
-
-static PixelList **DestroyPixelListThreadSet(PixelList **pixel_list)
-{
-  register ssize_t
-    i;
-
-  assert(pixel_list != (PixelList **) NULL);
-  for (i=0; i < (ssize_t) GetOpenMPMaximumThreads(); i++)
-    if (pixel_list[i] != (PixelList *) NULL)
-      pixel_list[i]=DestroyPixelList(pixel_list[i]);
-  pixel_list=(PixelList **) RelinquishMagickMemory(pixel_list);
-  return(pixel_list);
-}
-
-static PixelList *AcquirePixelList(const size_t width,const size_t height)
-{
-  PixelList
-    *pixel_list;
-
-  register ssize_t
-    i;
-
-  pixel_list=(PixelList *) AcquireMagickMemory(sizeof(*pixel_list));
-  if (pixel_list == (PixelList *) NULL)
-    return(pixel_list);
-  (void) ResetMagickMemory((void *) pixel_list,0,sizeof(*pixel_list));
-  pixel_list->length=width*height;
-  for (i=0; i < ListChannels; i++)
-  {
-    pixel_list->lists[i].nodes=(ListNode *) AcquireQuantumMemory(65537UL,
-      sizeof(*pixel_list->lists[i].nodes));
-    if (pixel_list->lists[i].nodes == (ListNode *) NULL)
-      return(DestroyPixelList(pixel_list));
-    (void) ResetMagickMemory(pixel_list->lists[i].nodes,0,65537UL*
-      sizeof(*pixel_list->lists[i].nodes));
-  }
-  pixel_list->signature=MagickSignature;
-  return(pixel_list);
-}
-
-static PixelList **AcquirePixelListThreadSet(const size_t width,
-  const size_t height)
-{
-  PixelList
-    **pixel_list;
-
-  register ssize_t
-    i;
-
-  size_t
-    number_threads;
-
-  number_threads=GetOpenMPMaximumThreads();
-  pixel_list=(PixelList **) AcquireQuantumMemory(number_threads,
-    sizeof(*pixel_list));
-  if (pixel_list == (PixelList **) NULL)
-    return((PixelList **) NULL);
-  (void) ResetMagickMemory(pixel_list,0,number_threads*sizeof(*pixel_list));
-  for (i=0; i < (ssize_t) number_threads; i++)
-  {
-    pixel_list[i]=AcquirePixelList(width,height);
-    if (pixel_list[i] == (PixelList *) NULL)
-      return(DestroyPixelListThreadSet(pixel_list));
-  }
-  return(pixel_list);
-}
-
-static void AddNodePixelList(PixelList *pixel_list,const ssize_t channel,
-  const size_t color)
-{
-  register SkipList
-    *list;
-
-  register ssize_t
-    level;
-
-  size_t
-    search,
-    update[9];
-
-  /*
-    Initialize the node.
-  */
-  list=pixel_list->lists+channel;
-  list->nodes[color].signature=pixel_list->signature;
-  list->nodes[color].count=1;
-  /*
-    Determine where it belongs in the list.
-  */
-  search=65536UL;
-  for (level=list->level; level >= 0; level--)
-  {
-    while (list->nodes[search].next[level] < color)
-      search=list->nodes[search].next[level];
-    update[level]=search;
-  }
-  /*
-    Generate a pseudo-random level for this node.
-  */
-  for (level=0; ; level++)
-  {
-    pixel_list->seed=(pixel_list->seed*42893621L)+1L;
-    if ((pixel_list->seed & 0x300) != 0x300)
-      break;
-  }
-  if (level > 8)
-    level=8;
-  if (level > (list->level+2))
-    level=list->level+2;
-  /*
-    If we're raising the list's level, link back to the root node.
-  */
-  while (level > list->level)
-  {
-    list->level++;
-    update[list->level]=65536UL;
-  }
-  /*
-    Link the node into the skip-list.
-  */
-  do
-  {
-    list->nodes[color].next[level]=list->nodes[update[level]].next[level];
-    list->nodes[update[level]].next[level]=color;
-  } while (level-- > 0);
-}
-
-static PixelInfo GetMaximumPixelList(PixelList *pixel_list)
-{
-  PixelInfo
-    pixel;
-
-  register SkipList
-    *list;
-
-  register ssize_t
-    channel;
-
-  size_t
-    color,
-    maximum;
-
-  ssize_t
-    count;
-
-  unsigned short
-    channels[ListChannels];
-
-  /*
-    Find the maximum value for each of the color.
-  */
-  for (channel=0; channel < 5; channel++)
-  {
-    list=pixel_list->lists+channel;
-    color=65536L;
-    count=0;
-    maximum=list->nodes[color].next[0];
-    do
-    {
-      color=list->nodes[color].next[0];
-      if (color > maximum)
-        maximum=color;
-      count+=list->nodes[color].count;
-    } while (count < (ssize_t) pixel_list->length);
-    channels[channel]=(unsigned short) maximum;
-  }
-  GetPixelInfo((const Image *) NULL,&pixel);
-  pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
-  pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
-  pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
-  pixel.alpha=(MagickRealType) ScaleShortToQuantum(channels[3]);
-  pixel.black=(MagickRealType) ScaleShortToQuantum(channels[4]);
-  return(pixel);
-}
-
-static PixelInfo GetMeanPixelList(PixelList *pixel_list)
-{
-  PixelInfo
-    pixel;
-
-  MagickRealType
-    sum;
-
-  register SkipList
-    *list;
-
-  register ssize_t
-    channel;
-
-  size_t
-    color;
-
-  ssize_t
-    count;
-
-  unsigned short
-    channels[ListChannels];
-
-  /*
-    Find the mean value for each of the color.
-  */
-  for (channel=0; channel < 5; channel++)
-  {
-    list=pixel_list->lists+channel;
-    color=65536L;
-    count=0;
-    sum=0.0;
-    do
-    {
-      color=list->nodes[color].next[0];
-      sum+=(MagickRealType) list->nodes[color].count*color;
-      count+=list->nodes[color].count;
-    } while (count < (ssize_t) pixel_list->length);
-    sum/=pixel_list->length;
-    channels[channel]=(unsigned short) sum;
-  }
-  GetPixelInfo((const Image *) NULL,&pixel);
-  pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
-  pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
-  pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
-  pixel.black=(MagickRealType) ScaleShortToQuantum(channels[4]);
-  pixel.alpha=(MagickRealType) ScaleShortToQuantum(channels[3]);
-  return(pixel);
-}
-
-static PixelInfo GetMedianPixelList(PixelList *pixel_list)
-{
-  PixelInfo
-    pixel;
-
-  register SkipList
-    *list;
-
-  register ssize_t
-    channel;
-
-  size_t
-    color;
-
-  ssize_t
-    count;
-
-  unsigned short
-    channels[ListChannels];
-
-  /*
-    Find the median value for each of the color.
-  */
-  for (channel=0; channel < 5; channel++)
-  {
-    list=pixel_list->lists+channel;
-    color=65536L;
-    count=0;
-    do
-    {
-      color=list->nodes[color].next[0];
-      count+=list->nodes[color].count;
-    } while (count <= (ssize_t) (pixel_list->length >> 1));
-    channels[channel]=(unsigned short) color;
-  }
-  GetPixelInfo((const Image *) NULL,&pixel);
-  pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
-  pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
-  pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
-  pixel.black=(MagickRealType) ScaleShortToQuantum(channels[4]);
-  pixel.alpha=(MagickRealType) ScaleShortToQuantum(channels[3]);
-  return(pixel);
-}
-
-static PixelInfo GetMinimumPixelList(PixelList *pixel_list)
-{
-  PixelInfo
-    pixel;
-
-  register SkipList
-    *list;
-
-  register ssize_t
-    channel;
-
-  size_t
-    color,
-    minimum;
-
-  ssize_t
-    count;
-
-  unsigned short
-    channels[ListChannels];
-
-  /*
-    Find the minimum value for each of the color.
-  */
-  for (channel=0; channel < 5; channel++)
-  {
-    list=pixel_list->lists+channel;
-    count=0;
-    color=65536UL;
-    minimum=list->nodes[color].next[0];
-    do
-    {
-      color=list->nodes[color].next[0];
-      if (color < minimum)
-        minimum=color;
-      count+=list->nodes[color].count;
-    } while (count < (ssize_t) pixel_list->length);
-    channels[channel]=(unsigned short) minimum;
-  }
-  GetPixelInfo((const Image *) NULL,&pixel);
-  pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
-  pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
-  pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
-  pixel.black=(MagickRealType) ScaleShortToQuantum(channels[4]);
-  pixel.alpha=(MagickRealType) ScaleShortToQuantum(channels[3]);
-  return(pixel);
-}
-
-static PixelInfo GetModePixelList(PixelList *pixel_list)
-{
-  PixelInfo
-    pixel;
-
-  register SkipList
-    *list;
-
-  register ssize_t
-    channel;
-
-  size_t
-    color,
-    max_count,
-    mode;
-
-  ssize_t
-    count;
-
-  unsigned short
-    channels[5];
-
-  /*
-    Make each pixel the 'predominant color' of the specified neighborhood.
-  */
-  for (channel=0; channel < 5; channel++)
-  {
-    list=pixel_list->lists+channel;
-    color=65536L;
-    mode=color;
-    max_count=list->nodes[mode].count;
-    count=0;
-    do
-    {
-      color=list->nodes[color].next[0];
-      if (list->nodes[color].count > max_count)
-        {
-          mode=color;
-          max_count=list->nodes[mode].count;
-        }
-      count+=list->nodes[color].count;
-    } while (count < (ssize_t) pixel_list->length);
-    channels[channel]=(unsigned short) mode;
-  }
-  GetPixelInfo((const Image *) NULL,&pixel);
-  pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
-  pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
-  pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
-  pixel.black=(MagickRealType) ScaleShortToQuantum(channels[4]);
-  pixel.alpha=(MagickRealType) ScaleShortToQuantum(channels[3]);
-  return(pixel);
-}
-
-static PixelInfo GetNonpeakPixelList(PixelList *pixel_list)
-{
-  PixelInfo
-    pixel;
-
-  register SkipList
-    *list;
-
-  register ssize_t
-    channel;
-
-  size_t
-    color,
-    next,
-    previous;
-
-  ssize_t
-    count;
+        PixelTrait
+          shade_traits,
+          traits;
 
-  unsigned short
-    channels[5];
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        shade_traits=GetPixelChannelTraits(shade_image,channel);
+        if ((traits == UndefinedPixelTrait) ||
+            (shade_traits == UndefinedPixelTrait))
+          continue;
+        if (((shade_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
+          {
+            SetPixelChannel(shade_image,channel,center[i],q);
+            continue;
+          }
+        if (gray != MagickFalse)
+          {
+            SetPixelChannel(shade_image,channel,ClampToQuantum(shade),q);
+            continue;
+          }
+        SetPixelChannel(shade_image,channel,ClampToQuantum(QuantumScale*shade*
+          center[i]),q);
+      }
+      pre+=GetPixelChannels(image);
+      center+=GetPixelChannels(image);
+      post+=GetPixelChannels(image);
+      q+=GetPixelChannels(shade_image);
+    }
+    if (SyncCacheViewAuthenticPixels(shade_view,exception) == MagickFalse)
+      status=MagickFalse;
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      {
+        MagickBooleanType
+          proceed;
 
-  /*
-    Finds the non peak value for each of the colors.
-  */
-  for (channel=0; channel < 5; channel++)
-  {
-    list=pixel_list->lists+channel;
-    color=65536L;
-    next=list->nodes[color].next[0];
-    count=0;
-    do
-    {
-      previous=color;
-      color=next;
-      next=list->nodes[color].next[0];
-      count+=list->nodes[color].count;
-    } while (count <= (ssize_t) (pixel_list->length >> 1));
-    if ((previous == 65536UL) && (next != 65536UL))
-      color=next;
-    else
-      if ((previous != 65536UL) && (next == 65536UL))
-        color=previous;
-    channels[channel]=(unsigned short) color;
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+        #pragma omp critical (MagickCore_ShadeImage)
+#endif
+        proceed=SetImageProgress(image,ShadeImageTag,progress++,image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
   }
-  GetPixelInfo((const Image *) NULL,&pixel);
-  pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
-  pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
-  pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
-  pixel.alpha=(MagickRealType) ScaleShortToQuantum(channels[3]);
-  pixel.black=(MagickRealType) ScaleShortToQuantum(channels[4]);
-  return(pixel);
+  shade_view=DestroyCacheView(shade_view);
+  image_view=DestroyCacheView(image_view);
+  if (status == MagickFalse)
+    shade_image=DestroyImage(shade_image);
+  return(shade_image);
 }
-
-static PixelInfo GetStandardDeviationPixelList(PixelList *pixel_list)
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%     S h a r p e n I m a g e                                                 %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  SharpenImage() sharpens the image.  We convolve the image with a Gaussian
+%  operator of the given radius and standard deviation (sigma).  For
+%  reasonable results, radius should be larger than sigma.  Use a radius of 0
+%  and SharpenImage() selects a suitable radius for you.
+%
+%  Using a separable kernel would be faster, but the negative weights cancel
+%  out on the corners of the kernel producing often undesirable ringing in the
+%  filtered result; this can be avoided by using a 2D gaussian shaped image
+%  sharpening kernel instead.
+%
+%  The format of the SharpenImage method is:
+%
+%    Image *SharpenImage(const Image *image,const double radius,
+%      const double sigma,ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o radius: the radius of the Gaussian, in pixels, not counting the center
+%      pixel.
+%
+%    o sigma: the standard deviation of the Laplacian, in pixels.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+MagickExport Image *SharpenImage(const Image *image,const double radius,
+  const double sigma,ExceptionInfo *exception)
 {
-  PixelInfo
-    pixel;
+  double
+    normalize;
 
-  MagickRealType
-    sum,
-    sum_squared;
+  Image
+    *sharp_image;
 
-  register SkipList
-    *list;
+  KernelInfo
+    *kernel_info;
 
   register ssize_t
-    channel;
+    i;
 
   size_t
-    color;
+    width;
 
   ssize_t
-    count;
-
-  unsigned short
-    channels[ListChannels];
+    j,
+    u,
+    v;
 
-  /*
-    Find the standard-deviation value for each of the color.
-  */
-  for (channel=0; channel < 5; channel++)
-  {
-    list=pixel_list->lists+channel;
-    color=65536L;
-    count=0;
-    sum=0.0;
-    sum_squared=0.0;
-    do
+  assert(image != (const Image *) NULL);
+  assert(image->signature == MagickSignature);
+  if (image->debug != MagickFalse)
+    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  assert(exception != (ExceptionInfo *) NULL);
+  assert(exception->signature == MagickSignature);
+  width=GetOptimalKernelWidth2D(radius,sigma);
+  kernel_info=AcquireKernelInfo((const char *) NULL);
+  if (kernel_info == (KernelInfo *) NULL)
+    ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
+  (void) ResetMagickMemory(kernel_info,0,sizeof(*kernel_info));
+  kernel_info->width=width;
+  kernel_info->height=width;
+  kernel_info->signature=MagickSignature;
+  kernel_info->values=(MagickRealType *) MagickAssumeAligned(
+    AcquireAlignedMemory(kernel_info->width,kernel_info->width*
+    sizeof(*kernel_info->values)));
+  if (kernel_info->values == (MagickRealType *) NULL)
     {
-      register ssize_t
-        i;
-
-      color=list->nodes[color].next[0];
-      sum+=(MagickRealType) list->nodes[color].count*color;
-      for (i=0; i < (ssize_t) list->nodes[color].count; i++)
-        sum_squared+=((MagickRealType) color)*((MagickRealType) color);
-      count+=list->nodes[color].count;
-    } while (count < (ssize_t) pixel_list->length);
-    sum/=pixel_list->length;
-    sum_squared/=pixel_list->length;
-    channels[channel]=(unsigned short) sqrt(sum_squared-(sum*sum));
-  }
-  GetPixelInfo((const Image *) NULL,&pixel);
-  pixel.red=(MagickRealType) ScaleShortToQuantum(channels[0]);
-  pixel.green=(MagickRealType) ScaleShortToQuantum(channels[1]);
-  pixel.blue=(MagickRealType) ScaleShortToQuantum(channels[2]);
-  pixel.alpha=(MagickRealType) ScaleShortToQuantum(channels[3]);
-  pixel.black=(MagickRealType) ScaleShortToQuantum(channels[4]);
-  return(pixel);
-}
-
-static inline void InsertPixelList(const Image *image,const Quantum *pixel,
-  PixelList *pixel_list)
-{
-  size_t
-    signature;
-
-  unsigned short
-    index;
-
-  index=ScaleQuantumToShort(GetPixelRed(image,pixel));
-  signature=pixel_list->lists[0].nodes[index].signature;
-  if (signature == pixel_list->signature)
-    pixel_list->lists[0].nodes[index].count++;
-  else
-    AddNodePixelList(pixel_list,0,index);
-  index=ScaleQuantumToShort(GetPixelGreen(image,pixel));
-  signature=pixel_list->lists[1].nodes[index].signature;
-  if (signature == pixel_list->signature)
-    pixel_list->lists[1].nodes[index].count++;
-  else
-    AddNodePixelList(pixel_list,1,index);
-  index=ScaleQuantumToShort(GetPixelBlue(image,pixel));
-  signature=pixel_list->lists[2].nodes[index].signature;
-  if (signature == pixel_list->signature)
-    pixel_list->lists[2].nodes[index].count++;
-  else
-    AddNodePixelList(pixel_list,2,index);
-  index=ScaleQuantumToShort(GetPixelAlpha(image,pixel));
-  signature=pixel_list->lists[3].nodes[index].signature;
-  if (signature == pixel_list->signature)
-    pixel_list->lists[3].nodes[index].count++;
-  else
-    AddNodePixelList(pixel_list,3,index);
-  if (image->colorspace == CMYKColorspace)
-    index=ScaleQuantumToShort(GetPixelBlack(image,pixel));
-  signature=pixel_list->lists[4].nodes[index].signature;
-  if (signature == pixel_list->signature)
-    pixel_list->lists[4].nodes[index].count++;
-  else
-    AddNodePixelList(pixel_list,4,index);
-}
-
-static inline MagickRealType MagickAbsoluteValue(const MagickRealType x)
-{
-  if (x < 0)
-    return(-x);
-  return(x);
-}
-
-static void ResetPixelList(PixelList *pixel_list)
-{
-  int
-    level;
-
-  register ListNode
-    *root;
-
-  register SkipList
-    *list;
-
-  register ssize_t
-    channel;
-
-  /*
-    Reset the skip-list.
-  */
-  for (channel=0; channel < 5; channel++)
+      kernel_info=DestroyKernelInfo(kernel_info);
+      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
+    }
+  normalize=0.0;
+  j=(ssize_t) kernel_info->width/2;
+  i=0;
+  for (v=(-j); v <= j; v++)
   {
-    list=pixel_list->lists+channel;
-    root=list->nodes+65536UL;
-    list->level=0;
-    for (level=0; level < 9; level++)
-      root->next[level]=65536UL;
+    for (u=(-j); u <= j; u++)
+    {
+      kernel_info->values[i]=(MagickRealType) (-exp(-((double) u*u+v*v)/(2.0*
+        MagickSigma*MagickSigma))/(2.0*MagickPI*MagickSigma*MagickSigma));
+      normalize+=kernel_info->values[i];
+      i++;
+    }
   }
-  pixel_list->seed=pixel_list->signature++;
+  kernel_info->values[i/2]=(double) ((-2.0)*normalize);
+  sharp_image=ConvolveImage(image,kernel_info,exception);
+  if (sharp_image != (Image *) NULL)
+    (void) ClampImage(sharp_image,exception);
+  kernel_info=DestroyKernelInfo(kernel_info);
+  return(sharp_image);
 }
-
-MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
-  const size_t width,const size_t height,ExceptionInfo *exception)
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%     S p r e a d I m a g e                                                   %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  SpreadImage() is a special effects method that randomly displaces each
+%  pixel in a block defined by the radius parameter.
+%
+%  The format of the SpreadImage method is:
+%
+%      Image *SpreadImage(const Image *image,const double radius,
+%        const PixelInterpolateMethod method,ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o radius:  choose a random pixel in a neighborhood of this extent.
+%
+%    o method:  the pixel interpolation method.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+MagickExport Image *SpreadImage(const Image *image,const double radius,
+  const PixelInterpolateMethod method,ExceptionInfo *exception)
 {
-#define StatisticWidth \
-  (width == 0 ? GetOptimalKernelWidth2D((double) width,0.5) : width)
-#define StatisticHeight \
-  (height == 0 ? GetOptimalKernelWidth2D((double) height,0.5) : height)
-#define StatisticImageTag  "Statistic/Image"
+#define SpreadImageTag  "Spread/Image"
 
   CacheView
     *image_view,
-    *statistic_view;
+    *spread_view;
 
   Image
-    *statistic_image;
+    *spread_image;
 
   MagickBooleanType
     status;
@@ -4625,14 +3585,22 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
   MagickOffsetType
     progress;
 
-  PixelList
-    **restrict pixel_list;
+  RandomInfo
+    **restrict random_info;
+
+  size_t
+    width;
 
   ssize_t
     y;
 
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  unsigned long
+    key;
+#endif
+
   /*
-    Initialize statistics image attributes.
+    Initialize spread image attributes.
   */
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
@@ -4640,32 +3608,32 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
-  statistic_image=CloneImage(image,image->columns,image->rows,MagickTrue,
+  spread_image=CloneImage(image,image->columns,image->rows,MagickTrue,
     exception);
-  if (statistic_image == (Image *) NULL)
+  if (spread_image == (Image *) NULL)
     return((Image *) NULL);
-  if (SetImageStorageClass(statistic_image,DirectClass,exception) == MagickFalse)
+  if (SetImageStorageClass(spread_image,DirectClass,exception) == MagickFalse)
     {
-      statistic_image=DestroyImage(statistic_image);
+      spread_image=DestroyImage(spread_image);
       return((Image *) NULL);
     }
-  pixel_list=AcquirePixelListThreadSet(StatisticWidth,StatisticHeight);
-  if (pixel_list == (PixelList **) NULL)
-    {
-      statistic_image=DestroyImage(statistic_image);
-      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-    }
   /*
-    Make each pixel the min / max / median / mode / etc. of the neighborhood.
+    Spread image.
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
-  statistic_view=AcquireCacheView(statistic_image);
+  width=GetOptimalKernelWidth1D(radius,0.5);
+  random_info=AcquireRandomInfoThreadSet();
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  key=GetRandomSecretKey(random_info[0]);
+#endif
+  image_view=AcquireVirtualCacheView(image,exception);
+  spread_view=AcquireAuthenticCacheView(spread_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,key == ~0UL)
 #endif
-  for (y=0; y < (ssize_t) statistic_image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
     const int
       id = GetOpenMPThreadId();
@@ -4681,109 +3649,26 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
 
     if (status == MagickFalse)
       continue;
-    p=GetCacheViewVirtualPixels(image_view,-((ssize_t) StatisticWidth/2L),y-
-      (ssize_t) (StatisticHeight/2L),image->columns+StatisticWidth,
-      StatisticHeight,exception);
-    q=QueueCacheViewAuthenticPixels(statistic_view,0,y,statistic_image->columns,      1,exception);
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+    q=QueueCacheViewAuthenticPixels(spread_view,0,y,spread_image->columns,1,
+      exception);
     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
       {
         status=MagickFalse;
         continue;
       }
-    for (x=0; x < (ssize_t) statistic_image->columns; x++)
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
-      PixelInfo
-        pixel;
-
-      register const Quantum
-        *restrict r;
-
-      register ssize_t
-        u,
-        v;
+      PointInfo
+        point;
 
-      r=p;
-      ResetPixelList(pixel_list[id]);
-      for (v=0; v < (ssize_t) StatisticHeight; v++)
-      {
-        for (u=0; u < (ssize_t) StatisticWidth; u++)
-          InsertPixelList(image,r+u*GetPixelChannels(image),pixel_list[id]);
-        r+=(image->columns+StatisticWidth)*GetPixelChannels(image);
-      }
-      GetPixelInfo(image,&pixel);
-      SetPixelInfo(image,p+(StatisticWidth*StatisticHeight/2)*
-        GetPixelChannels(image),&pixel);
-      switch (type)
-      {
-        case GradientStatistic:
-        {
-          PixelInfo
-            maximum,
-            minimum;
-
-          minimum=GetMinimumPixelList(pixel_list[id]);
-          maximum=GetMaximumPixelList(pixel_list[id]);
-          pixel.red=MagickAbsoluteValue(maximum.red-minimum.red);
-          pixel.green=MagickAbsoluteValue(maximum.green-minimum.green);
-          pixel.blue=MagickAbsoluteValue(maximum.blue-minimum.blue);
-          pixel.alpha=MagickAbsoluteValue(maximum.alpha-minimum.alpha);
-          if (image->colorspace == CMYKColorspace)
-            pixel.black=MagickAbsoluteValue(maximum.black-minimum.black);
-          break;
-        }
-        case MaximumStatistic:
-        {
-          pixel=GetMaximumPixelList(pixel_list[id]);
-          break;
-        }
-        case MeanStatistic:
-        {
-          pixel=GetMeanPixelList(pixel_list[id]);
-          break;
-        }
-        case MedianStatistic:
-        default:
-        {
-          pixel=GetMedianPixelList(pixel_list[id]);
-          break;
-        }
-        case MinimumStatistic:
-        {
-          pixel=GetMinimumPixelList(pixel_list[id]);
-          break;
-        }
-        case ModeStatistic:
-        {
-          pixel=GetModePixelList(pixel_list[id]);
-          break;
-        }
-        case NonpeakStatistic:
-        {
-          pixel=GetNonpeakPixelList(pixel_list[id]);
-          break;
-        }
-        case StandardDeviationStatistic:
-        {
-          pixel=GetStandardDeviationPixelList(pixel_list[id]);
-          break;
-        }
-      }
-      if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
-        SetPixelRed(statistic_image,ClampToQuantum(pixel.red),q);
-      if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
-        SetPixelGreen(statistic_image,ClampToQuantum(pixel.green),q);
-      if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
-        SetPixelBlue(statistic_image,ClampToQuantum(pixel.blue),q);
-      if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
-          (image->colorspace == CMYKColorspace))
-        SetPixelBlack(statistic_image,ClampToQuantum(pixel.black),q);
-      if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
-          (image->matte != MagickFalse))
-        SetPixelAlpha(statistic_image,ClampToQuantum(pixel.alpha),q);
-      p+=GetPixelChannels(image);
-      q+=GetPixelChannels(statistic_image);
+      point.x=GetPseudoRandomValue(random_info[id]);
+      point.y=GetPseudoRandomValue(random_info[id]);
+      status=InterpolatePixelChannels(image,image_view,spread_image,method,
+        (double) x+width*point.x-0.5,(double) y+width*point.y-0.5,q,exception);
+      q+=GetPixelChannels(spread_image);
     }
-    if (SyncCacheViewAuthenticPixels(statistic_view,exception) == MagickFalse)
+    if (SyncCacheViewAuthenticPixels(spread_view,exception) == MagickFalse)
       status=MagickFalse;
     if (image->progress_monitor != (MagickProgressMonitor) NULL)
       {
@@ -4791,18 +3676,17 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_StatisticImage)
+        #pragma omp critical (MagickCore_SpreadImage)
 #endif
-        proceed=SetImageProgress(image,StatisticImageTag,progress++,
-          image->rows);
+        proceed=SetImageProgress(image,SpreadImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
           status=MagickFalse;
       }
   }
-  statistic_view=DestroyCacheView(statistic_view);
+  spread_view=DestroyCacheView(spread_view);
   image_view=DestroyCacheView(image_view);
-  pixel_list=DestroyPixelListThreadSet(pixel_list);
-  return(statistic_image);
+  random_info=DestroyRandomInfoThreadSet(random_info);
+  return(spread_image);
 }
 \f
 /*
@@ -4844,9 +3728,9 @@ MagickExport Image *StatisticImage(const Image *image,const StatisticType type,
 %    o exception: return any errors or warnings in this structure.
 %
 */
-MagickExport Image *UnsharpMaskImage(const Image *image,
-  const double radius,const double sigma,const double amount,
-  const double threshold,ExceptionInfo *exception)
+MagickExport Image *UnsharpMaskImage(const Image *image,const double radius,
+  const double sigma,const double amount,const double threshold,
+  ExceptionInfo *exception)
 {
 #define SharpenImageTag  "Sharpen/Image"
 
@@ -4863,7 +3747,7 @@ MagickExport Image *UnsharpMaskImage(const Image *image,
   MagickOffsetType
     progress;
 
-  MagickRealType
+  double
     quantum_threshold;
 
   ssize_t
@@ -4874,19 +3758,20 @@ MagickExport Image *UnsharpMaskImage(const Image *image,
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   assert(exception != (ExceptionInfo *) NULL);
-  unsharp_image=BlurImage(image,radius,sigma,image->bias,exception);
+  unsharp_image=BlurImage(image,radius,sigma,exception);
   if (unsharp_image == (Image *) NULL)
     return((Image *) NULL);
-  quantum_threshold=(MagickRealType) QuantumRange*threshold;
+  quantum_threshold=(double) QuantumRange*threshold;
   /*
     Unsharp-mask image.
   */
   status=MagickTrue;
   progress=0;
-  image_view=AcquireCacheView(image);
-  unsharp_view=AcquireCacheView(unsharp_image);
+  image_view=AcquireVirtualCacheView(image,exception);
+  unsharp_view=AcquireAuthenticCacheView(unsharp_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+  #pragma omp parallel for schedule(static,4) shared(progress,status) \
+    dynamic_number_threads(image,image->columns,image->rows,1)
 #endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
@@ -4902,7 +3787,7 @@ MagickExport Image *UnsharpMaskImage(const Image *image,
     if (status == MagickFalse)
       continue;
     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
-    q=GetCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
+    q=QueueCacheViewAuthenticPixels(unsharp_view,0,y,unsharp_image->columns,1,
       exception);
     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
       {
@@ -4916,7 +3801,7 @@ MagickExport Image *UnsharpMaskImage(const Image *image,
 
       for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
       {
-        MagickRealType
+        double
           pixel;
 
         PixelChannel
@@ -4926,23 +3811,24 @@ MagickExport Image *UnsharpMaskImage(const Image *image,
           traits,
           unsharp_traits;
 
-        traits=GetPixelChannelMapTraits(image,(PixelChannel) i);
-        channel=GetPixelChannelMapChannel(image,(PixelChannel) i);
-        unsharp_traits=GetPixelChannelMapTraits(unsharp_image,channel);
+        channel=GetPixelChannelChannel(image,i);
+        traits=GetPixelChannelTraits(image,channel);
+        unsharp_traits=GetPixelChannelTraits(unsharp_image,channel);
         if ((traits == UndefinedPixelTrait) ||
             (unsharp_traits == UndefinedPixelTrait))
           continue;
-        if ((unsharp_traits & CopyPixelTrait) != 0)
+        if (((unsharp_traits & CopyPixelTrait) != 0) ||
+            (GetPixelMask(image,p) != 0))
           {
-            q[channel]=p[i];
+            SetPixelChannel(unsharp_image,channel,p[i],q);
             continue;
           }
-        pixel=p[i]-(MagickRealType) q[channel];
+        pixel=p[i]-(double) GetPixelChannel(unsharp_image,channel,q);
         if (fabs(2.0*pixel) < quantum_threshold)
-          pixel=(MagickRealType) p[i];
+          pixel=(double) p[i];
         else
-          pixel=(MagickRealType) p[i]+amount*pixel;
-        q[channel]=ClampToQuantum(pixel);
+          pixel=(double) p[i]+amount*pixel;
+        SetPixelChannel(unsharp_image,channel,ClampToQuantum(pixel),q);
       }
       p+=GetPixelChannels(image);
       q+=GetPixelChannels(unsharp_image);
@@ -4955,7 +3841,7 @@ MagickExport Image *UnsharpMaskImage(const Image *image,
           proceed;
 
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_UnsharpMaskImage)
+        #pragma omp critical (MagickCore_UnsharpMaskImage)
 #endif
         proceed=SetImageProgress(image,SharpenImageTag,progress++,image->rows);
         if (proceed == MagickFalse)
@@ -4964,6 +3850,8 @@ MagickExport Image *UnsharpMaskImage(const Image *image,
   }
   unsharp_image->type=image->type;
   unsharp_view=DestroyCacheView(unsharp_view);
+  if (unsharp_image != (Image *) NULL)
+    (void) ClampImage(unsharp_image,exception);
   image_view=DestroyCacheView(image_view);
   if (status == MagickFalse)
     unsharp_image=DestroyImage(unsharp_image);