]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/pixel.c
(no commit message)
[imagemagick] / MagickCore / pixel.c
index 4f32a82816be1b6c6b44caf3225ec721772dfa1e..fbf27e31d28cb584945c093aeff0390ccc994a7a 100644 (file)
 %                  MagickCore Methods to Import/Export Pixels                 %
 %                                                                             %
 %                             Software Design                                 %
-%                               John Cristy                                   %
+%                                  Cristy                                     %
 %                               October 1998                                  %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2014 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  %
 #include "MagickCore/transform.h"
 #include "MagickCore/utility.h"
 \f
-#define LogPixelChannels(image) \
-{ \
-  register ssize_t \
-    i; \
- \
-  (void) LogMagickEvent(PixelEvent,GetMagickModule(),"%s[%.20g]", \
-    image->filename,(double) image->number_channels); \
-  for (i=0; i < (ssize_t) image->number_channels; i++) \
-  { \
-    char \
-      traits[MaxTextExtent]; \
- \
-    const char \
-      *name; \
- \
-    PixelChannel \
-      channel; \
- \
-    switch (GetPixelChannelMapChannel(image,i)) \
-    { \
-      case RedPixelChannel: \
-      { \
-        name="red"; \
-        if (image->colorspace == CMYKColorspace) \
-          name="cyan"; \
-        if (image->colorspace == GRAYColorspace) \
-          name="gray"; \
-        break; \
-      } \
-      case GreenPixelChannel: \
-      { \
-        name="green"; \
-        if (image->colorspace == CMYKColorspace) \
-          name="magenta"; \
-        break; \
-      } \
-      case BluePixelChannel: \
-      { \
-        name="blue"; \
-        if (image->colorspace == CMYKColorspace) \
-          name="yellow"; \
-        break; \
-      } \
-      case BlackPixelChannel: \
-      { \
-        name="black"; \
-        if (image->storage_class == PseudoClass) \
-          name="index"; \
-        break; \
-      } \
-      case IndexPixelChannel: \
-      { \
-        name="index"; \
-        break; \
-      } \
-      case AlphaPixelChannel: \
-      { \
-        name="alpha"; \
-        break; \
-      } \
-      case MaskPixelChannel: \
-      { \
-        name="mask"; \
-        break; \
-      } \
-      case MetaPixelChannel: \
-      { \
-        name="meta"; \
-        break; \
-      } \
-      default: \
-        name="undefined"; \
-    } \
-    channel=GetPixelChannelMapChannel(image,i); \
-    *traits='\0'; \
-    if ((GetPixelChannelMapTraits(image,channel) & UpdatePixelTrait) != 0) \
-      (void) ConcatenateMagickString(traits,"update,",MaxTextExtent); \
-    if ((GetPixelChannelMapTraits(image,channel) & BlendPixelTrait) != 0) \
-      (void) ConcatenateMagickString(traits,"blend,",MaxTextExtent); \
-    if ((GetPixelChannelMapTraits(image,channel) & CopyPixelTrait) != 0) \
-      (void) ConcatenateMagickString(traits,"copy,",MaxTextExtent); \
-    if (*traits == '\0') \
-      (void) ConcatenateMagickString(traits,"undefined,",MaxTextExtent); \
-    traits[strlen(traits)-1]='\0'; \
-    (void) LogMagickEvent(PixelEvent,GetMagickModule(),"  %.20g: %s (%s)", \
-      (double) i,name,traits); \
-  } \
-}
-\f
 /*
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %                                                                             %
@@ -271,6 +182,95 @@ MagickExport PixelInfo *ClonePixelInfo(const PixelInfo *pixel)
 %                                                                             %
 %                                                                             %
 %                                                                             %
+%   D e c o d e P i x e l G a m m a                                           %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  DecodePixelGamma() applies the expansive power-law nonlinearity to the pixel.
+%
+%  The format of the DecodePixelGamma method is:
+%
+%      double DecodePixelGamma(const MagickRealType pixel)
+%
+%  A description of each parameter follows:
+%
+%    o pixel: the pixel.
+%
+*/
+
+static inline double DecodeGamma(const double x)
+{
+  div_t
+    quotient;
+
+  double
+    p,
+    term[9];
+
+  int
+    exponent;
+
+  static const double coefficient[] =  /* terms for x^(7/5), x=1.5 */
+  {
+    1.7917488588043277509,
+    0.82045614371976854984,
+    0.027694100686325412819,
+    -0.00094244335181762134018,
+    0.000064355540911469709545,
+    -5.7224404636060757485e-06,
+    5.8767669437311184313e-07,
+    -6.6139920053589721168e-08,
+    7.9323242696227458163e-09
+  };
+
+  static const double powers_of_two[] =  /* (2^x)^(7/5) */
+  {
+    1.0,
+    2.6390158215457883983,
+    6.9644045063689921093,
+    1.8379173679952558018e+01,
+    4.8502930128332728543e+01
+  };
+
+  /*
+    Compute x^2.4 == x*x^(7/5) == pow(x,2.4).
+  */
+  term[0]=1.0;
+  term[1]=4.0*frexp(x,&exponent)-3.0;
+  term[2]=2.0*term[1]*term[1]-term[0];
+  term[3]=2.0*term[1]*term[2]-term[1];
+  term[4]=2.0*term[1]*term[3]-term[2];
+  term[5]=2.0*term[1]*term[4]-term[3];
+  term[6]=2.0*term[1]*term[5]-term[4];
+  term[7]=2.0*term[1]*term[6]-term[5];
+  term[8]=2.0*term[1]*term[7]-term[6];
+  p=coefficient[0]*term[0]+coefficient[1]*term[1]+coefficient[2]*term[2]+
+    coefficient[3]*term[3]+coefficient[4]*term[4]+coefficient[5]*term[5]+
+    coefficient[6]*term[6]+coefficient[7]*term[7]+coefficient[8]*term[8];
+  quotient=div(exponent-1,5);
+  if (quotient.rem < 0)
+    {
+      quotient.quot-=1;
+      quotient.rem+=5;
+    }
+  return(x*ldexp(powers_of_two[quotient.rem]*p,7*quotient.quot));
+}
+
+MagickExport MagickRealType DecodePixelGamma(const MagickRealType pixel)
+{
+  if (pixel <= (0.0404482362771076*QuantumRange))
+    return(pixel/12.92f);
+  return((MagickRealType) (QuantumRange*DecodeGamma((double) (QuantumScale*
+    pixel+0.055)/1.055)));
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
 +   D e s t r o y P i x e l C h a n n e l M a p                               %
 %                                                                             %
 %                                                                             %
@@ -302,6 +302,102 @@ MagickExport PixelChannelMap *DestroyPixelChannelMap(
 %                                                                             %
 %                                                                             %
 %                                                                             %
++   E n c o d e P i x e l G a m m a                                           %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  EncodePixelGamma() cancels any nonlinearity in the pixel.
+%
+%  The format of the EncodePixelGamma method is:
+%
+%      MagickRealType EncodePixelGamma(const double MagickRealType)
+%
+%  A description of each parameter follows:
+%
+%    o pixel: the pixel.
+%
+*/
+
+static inline double EncodeGamma(const double x)
+{
+  div_t
+    quotient;
+
+  double
+    p,
+    term[9];
+
+  int
+    exponent;
+
+  static const double coefficient[] =  /* Chebychevi poly: x^(5/12), x=1.5 */
+  {
+    1.1758200232996901923,
+    0.16665763094889061230,
+    -0.0083154894939042125035,
+    0.00075187976780420279038,
+    -0.000083240178519391795367,
+    0.000010229209410070008679,
+    -1.3400466409860246e-06,
+    1.8333422241635376682e-07,
+    -2.5878596761348859722e-08
+  };
+
+  static const double powers_of_two[] =  /* (2^N)^(5/12) */
+  {
+    1.0,
+    1.3348398541700343678,
+    1.7817974362806785482,
+    2.3784142300054420538,
+    3.1748021039363991669,
+    4.2378523774371812394,
+    5.6568542494923805819,
+    7.5509945014535482244,
+    1.0079368399158985525e1,
+    1.3454342644059433809e1,
+    1.7959392772949968275e1,
+    2.3972913230026907883e1
+  };
+
+  /*
+    Compute x^(1/2.4) == x^(5/12) == pow(x,1.0/2.4).
+  */
+  term[0]=1.0;
+  term[1]=4.0*frexp(x,&exponent)-3.0;
+  term[2]=2.0*term[1]*term[1]-term[0];
+  term[3]=2.0*term[1]*term[2]-term[1];
+  term[4]=2.0*term[1]*term[3]-term[2];
+  term[5]=2.0*term[1]*term[4]-term[3];
+  term[6]=2.0*term[1]*term[5]-term[4];
+  term[7]=2.0*term[1]*term[6]-term[5];
+  term[8]=2.0*term[1]*term[7]-term[6];
+  p=coefficient[0]*term[0]+coefficient[1]*term[1]+coefficient[2]*term[2]+
+    coefficient[3]*term[3]+coefficient[4]*term[4]+coefficient[5]*term[5]+
+    coefficient[6]*term[6]+coefficient[7]*term[7]+coefficient[8]*term[8];
+  quotient=div(exponent-1,12);
+  if (quotient.rem < 0)
+    {
+      quotient.quot-=1;
+      quotient.rem+=12;
+    }
+  return(ldexp(powers_of_two[quotient.rem]*p,5*quotient.quot));
+}
+
+MagickExport MagickRealType EncodePixelGamma(const MagickRealType pixel)
+{
+  if (pixel <= (0.0031306684425005883*QuantumRange))
+    return(12.92f*pixel);
+  return((MagickRealType) QuantumRange*(1.055*EncodeGamma((double) QuantumScale*
+    pixel)-0.055));
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
 %   E x p o r t I m a g e P i x e l s                                         %
 %                                                                             %
 %                                                                             %
@@ -427,7 +523,6 @@ static void ExportCharPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         p=GetVirtualPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -435,7 +530,7 @@ static void ExportCharPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          *q++=ScaleQuantumToChar(GetPixelIntensity(image,p));
+          *q++=ScaleQuantumToChar(ClampToQuantum(GetPixelIntensity(image,p)));
           p+=GetPixelChannels(image);
         }
       }
@@ -546,7 +641,7 @@ static void ExportCharPixel(Image *image,const RectangleInfo *roi,
           }
           case IndexQuantum:
           {
-            *q=ScaleQuantumToChar(GetPixelIntensity(image,p));
+            *q=ScaleQuantumToChar(ClampToQuantum(GetPixelIntensity(image,p)));
             break;
           }
           default:
@@ -634,7 +729,6 @@ static void ExportDoublePixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         p=GetVirtualPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -842,7 +936,6 @@ static void ExportFloatPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         p=GetVirtualPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -1049,7 +1142,6 @@ static void ExportLongPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         p=GetVirtualPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -1057,7 +1149,7 @@ static void ExportLongPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          *q++=ScaleQuantumToLong(GetPixelIntensity(image,p));
+          *q++=ScaleQuantumToLong(ClampToQuantum(GetPixelIntensity(image,p)));
           p+=GetPixelChannels(image);
         }
       }
@@ -1168,7 +1260,7 @@ static void ExportLongPixel(Image *image,const RectangleInfo *roi,
           }
           case IndexQuantum:
           {
-            *q=ScaleQuantumToLong(GetPixelIntensity(image,p));
+            *q=ScaleQuantumToLong(ClampToQuantum(GetPixelIntensity(image,p)));
             break;
           }
           default:
@@ -1256,7 +1348,6 @@ static void ExportLongLongPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         p=GetVirtualPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -1264,7 +1355,7 @@ static void ExportLongLongPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          *q++=ScaleQuantumToLongLong(GetPixelIntensity(image,p));
+          *q++=ScaleQuantumToLongLong(ClampToQuantum(GetPixelIntensity(image,p)));
           p+=GetPixelChannels(image);
         }
       }
@@ -1375,7 +1466,8 @@ static void ExportLongLongPixel(Image *image,const RectangleInfo *roi,
           }
           case IndexQuantum:
           {
-            *q=ScaleQuantumToLongLong(GetPixelIntensity(image,p));
+            *q=ScaleQuantumToLongLong(ClampToQuantum(
+              GetPixelIntensity(image,p)));
             break;
           }
           default:
@@ -1463,7 +1555,6 @@ static void ExportQuantumPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         p=GetVirtualPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -1471,7 +1562,7 @@ static void ExportQuantumPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          *q++=GetPixelIntensity(image,p);
+          *q++=ClampToQuantum(GetPixelIntensity(image,p));
           p+=GetPixelChannels(image);
         }
       }
@@ -1582,7 +1673,7 @@ static void ExportQuantumPixel(Image *image,const RectangleInfo *roi,
           }
           case IndexQuantum:
           {
-            *q=(GetPixelIntensity(image,p));
+            *q=ClampToQuantum(GetPixelIntensity(image,p));
             break;
           }
           default:
@@ -1673,7 +1764,6 @@ static void ExportShortPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         p=GetVirtualPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -1681,7 +1771,7 @@ static void ExportShortPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          *q++=ScaleQuantumToShort(GetPixelIntensity(image,p));
+          *q++=ScaleQuantumToShort(ClampToQuantum(GetPixelIntensity(image,p)));
           p+=GetPixelChannels(image);
         }
       }
@@ -1792,7 +1882,7 @@ static void ExportShortPixel(Image *image,const RectangleInfo *roi,
           }
           case IndexQuantum:
           {
-            *q=ScaleQuantumToShort(GetPixelIntensity(image,p));
+            *q=ScaleQuantumToShort(ClampToQuantum(GetPixelIntensity(image,p)));
             break;
           }
           default:
@@ -1805,9 +1895,9 @@ static void ExportShortPixel(Image *image,const RectangleInfo *roi,
   }
 }
 
-MagickExport MagickBooleanType ExportImagePixels(Image *image,
-  const ssize_t x,const ssize_t y,const size_t width,const size_t height,
-  const char *map,const StorageType type,void *pixels,ExceptionInfo *exception)
+MagickExport MagickBooleanType ExportImagePixels(Image *image,const ssize_t x,
+  const ssize_t y,const size_t width,const size_t height,const char *map,
+  const StorageType type,void *pixels,ExceptionInfo *exception)
 {
   QuantumType
     *quantum_map;
@@ -1830,7 +1920,7 @@ MagickExport MagickBooleanType ExportImagePixels(Image *image,
   if (quantum_map == (QuantumType *) NULL)
     {
       (void) ThrowMagickException(exception,GetMagickModule(),
-        ResourceLimitError,"MemoryAllocationFailed","'%s'",image->filename);
+        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
       return(MagickFalse);
     }
   for (i=0; i < (ssize_t) length; i++)
@@ -1857,7 +1947,7 @@ MagickExport MagickBooleanType ExportImagePixels(Image *image,
           break;
         quantum_map=(QuantumType *) RelinquishMagickMemory(quantum_map);
         (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
-          "ColorSeparatedImageRequired","'%s'",map);
+          "ColorSeparatedImageRequired","`%s'",map);
         return(MagickFalse);
       }
       case 'g':
@@ -1880,7 +1970,7 @@ MagickExport MagickBooleanType ExportImagePixels(Image *image,
           break;
         quantum_map=(QuantumType *) RelinquishMagickMemory(quantum_map);
         (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
-          "ColorSeparatedImageRequired","'%s'",map);
+          "ColorSeparatedImageRequired","`%s'",map);
         return(MagickFalse);
       }
       case 'M':
@@ -1891,7 +1981,7 @@ MagickExport MagickBooleanType ExportImagePixels(Image *image,
           break;
         quantum_map=(QuantumType *) RelinquishMagickMemory(quantum_map);
         (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
-          "ColorSeparatedImageRequired","'%s'",map);
+          "ColorSeparatedImageRequired","`%s'",map);
         return(MagickFalse);
       }
       case 'o':
@@ -1920,14 +2010,14 @@ MagickExport MagickBooleanType ExportImagePixels(Image *image,
           break;
         quantum_map=(QuantumType *) RelinquishMagickMemory(quantum_map);
         (void) ThrowMagickException(exception,GetMagickModule(),ImageError,
-          "ColorSeparatedImageRequired","'%s'",map);
+          "ColorSeparatedImageRequired","`%s'",map);
         return(MagickFalse);
       }
       default:
       {
         quantum_map=(QuantumType *) RelinquishMagickMemory(quantum_map);
         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-          "UnrecognizedPixelMap","'%s'",map);
+          "UnrecognizedPixelMap","`%s'",map);
         return(MagickFalse);
       }
     }
@@ -1977,7 +2067,7 @@ MagickExport MagickBooleanType ExportImagePixels(Image *image,
     {
       quantum_map=(QuantumType *) RelinquishMagickMemory(quantum_map);
       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-        "UnrecognizedPixelMap","'%s'",map);
+        "UnrecognizedPixelMap","`%s'",map);
       break;
     }
   }
@@ -2004,7 +2094,7 @@ MagickExport MagickBooleanType ExportImagePixels(Image *image,
 %
 %  A description of each parameter follows:
 %
-%    o image: the image.
+%    o image: the image. (optional - may be NULL)
 %
 %    o pixel: Specifies a pointer to a PixelInfo structure.
 %
@@ -2013,20 +2103,20 @@ MagickExport void GetPixelInfo(const Image *image,PixelInfo *pixel)
 {
   pixel->storage_class=DirectClass;
   pixel->colorspace=sRGBColorspace;
-  pixel->matte=MagickFalse;
+  pixel->alpha_trait=UndefinedPixelTrait;
   pixel->fuzz=0.0;
   pixel->depth=MAGICKCORE_QUANTUM_DEPTH;
   pixel->red=0.0;
   pixel->green=0.0;
   pixel->blue=0.0;
   pixel->black=0.0;
-  pixel->alpha=(MagickRealType) OpaqueAlpha;
+  pixel->alpha=(double) OpaqueAlpha;
   pixel->index=0.0;
   if (image == (const Image *) NULL)
     return;
   pixel->storage_class=image->storage_class;
   pixel->colorspace=image->colorspace;
-  pixel->matte=image->matte;
+  pixel->alpha_trait=image->alpha_trait;
   pixel->depth=image->depth;
   pixel->fuzz=image->fuzz;
 }
@@ -2036,6 +2126,153 @@ MagickExport void GetPixelInfo(const Image *image,PixelInfo *pixel)
 %                                                                             %
 %                                                                             %
 %                                                                             %
+%   G e t P i x e l I n t e n s i t y                                         %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  GetPixelIntensity() returns a single sample intensity value from the red,
+%  green, and blue components of a pixel based on the selected method:
+%
+%    Rec601Luma       0.298839R' + 0.586811G' + 0.114350B'
+%    Rec601Luminance  0.298839R + 0.586811G + 0.114350B
+%    Rec709Luma       0.212656R' + 0.715158G' + 0.072186B'
+%    Rec709Luminance  0.212656R + 0.715158G + 0.072186B
+%    Brightness       max(R', G', B')
+%    Lightness        (min(R', G', B') + max(R', G', B')) / 2.0
+% 
+%    MS               (R^2 + G^2 + B^2) / 3.0
+%    RMS              sqrt((R^2 + G^2 + B^2) / 3.0
+%    Average          (R + G + B') / 3.0
+%
+%  The format of the GetPixelIntensity method is:
+%
+%      MagickRealType GetPixelIntensity(const Image *image,
+%        const Quantum *pixel)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o pixel: Specifies a pointer to a Quantum structure.
+%
+*/
+
+static inline MagickRealType MagickMax(const MagickRealType x,
+  const MagickRealType y)
+{
+  if (x > y)
+    return(x);
+  return(y);
+}
+
+static inline MagickRealType MagickMin(const MagickRealType x,
+  const MagickRealType y)
+{
+  if (x < y)
+    return(x);
+  return(y);
+}
+
+MagickExport MagickRealType GetPixelIntensity(const Image *restrict image,
+  const Quantum *restrict pixel)
+{
+  MagickRealType
+    blue,
+    green,
+    red,
+    intensity;
+
+  if (image->colorspace == GRAYColorspace)
+    return((MagickRealType) GetPixelGray(image,pixel));
+  red=(MagickRealType) GetPixelRed(image,pixel);
+  green=(MagickRealType) GetPixelGreen(image,pixel);
+  blue=(MagickRealType) GetPixelBlue(image,pixel);
+  switch (image->intensity)
+  {
+    case AveragePixelIntensityMethod:
+    {
+      intensity=(red+green+blue)/3.0;
+      break;
+    }
+    case BrightnessPixelIntensityMethod:
+    {
+      intensity=MagickMax(MagickMax(red,green),blue);
+      break;
+    }
+    case LightnessPixelIntensityMethod:
+    {
+      intensity=(MagickMin(MagickMin(red,green),blue)+
+        MagickMax(MagickMax(red,green),blue))/2.0;
+      break;
+    }
+    case MSPixelIntensityMethod:
+    {
+      intensity=(MagickRealType) (((double) red*red+green*green+blue*blue)/
+        (3.0*QuantumRange));
+      break;
+    }
+    case Rec601LumaPixelIntensityMethod:
+    {
+      if (image->colorspace == RGBColorspace)
+        {
+          red=EncodePixelGamma(red);
+          green=EncodePixelGamma(green);
+          blue=EncodePixelGamma(blue);
+        }
+      intensity=0.298839*red+0.586811*green+0.114350*blue;
+      break;
+    }
+    case Rec601LuminancePixelIntensityMethod:
+    {
+      if (image->colorspace == sRGBColorspace)
+        {
+          red=DecodePixelGamma(red);
+          green=DecodePixelGamma(green);
+          blue=DecodePixelGamma(blue);
+        }
+      intensity=0.298839*red+0.586811*green+0.114350*blue;
+      break;
+    }
+    case Rec709LumaPixelIntensityMethod:
+    default:
+    {
+      if (image->colorspace == RGBColorspace)
+        {
+          red=EncodePixelGamma(red);
+          green=EncodePixelGamma(green);
+          blue=EncodePixelGamma(blue);
+        }
+      intensity=0.212656*red+0.715158*green+0.072186*blue;
+      break;
+    }
+    case Rec709LuminancePixelIntensityMethod:
+    {
+      if (image->colorspace == sRGBColorspace)
+        {
+          red=DecodePixelGamma(red);
+          green=DecodePixelGamma(green);
+          blue=DecodePixelGamma(blue);
+        }
+      intensity=0.212656*red+0.715158*green+0.072186*blue;
+      break;
+    }
+    case RMSPixelIntensityMethod:
+    {
+      intensity=(MagickRealType) (sqrt((double) red*red+green*green+blue*blue)/
+        sqrt(3.0));
+      break;
+    }
+  }
+  return(intensity);
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
 %   I m p o r t I m a g e P i x e l s                                         %
 %                                                                             %
 %                                                                             %
@@ -2188,7 +2425,6 @@ static void ImportCharPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         q=GetAuthenticPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -2377,14 +2613,11 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2402,17 +2635,13 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelAlpha(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelAlpha(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2430,14 +2659,11 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           p++;
           q+=GetPixelChannels(image);
@@ -2449,7 +2675,6 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         q=GetAuthenticPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -2457,8 +2682,7 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelGray(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGray(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2476,14 +2700,11 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2501,17 +2722,13 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelAlpha(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelAlpha(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2529,14 +2746,11 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2563,46 +2777,39 @@ static void ImportDoublePixel(Image *image,const RectangleInfo *roi,
           case RedQuantum:
           case CyanQuantum:
           {
-            SetPixelRed(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case GreenQuantum:
           case MagentaQuantum:
           {
-            SetPixelGreen(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case BlueQuantum:
           case YellowQuantum:
           {
-            SetPixelBlue(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case AlphaQuantum:
           {
-            SetPixelAlpha(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelAlpha(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case OpacityQuantum:
           {
-            SetPixelAlpha(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelAlpha(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case BlackQuantum:
           {
-            SetPixelBlack(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelBlack(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case IndexQuantum:
           {
-            SetPixelGray(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelGray(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           default:
@@ -2646,14 +2853,11 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2671,17 +2875,13 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelAlpha(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelAlpha(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2699,14 +2899,11 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           p++;
           q+=GetPixelChannels(image);
@@ -2718,7 +2915,6 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         q=GetAuthenticPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -2726,8 +2922,7 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelGray(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGray(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2745,14 +2940,11 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2770,17 +2962,13 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelRed(image,ClampToQuantum((MagickRealType)
-            QuantumRange*(*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelAlpha(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelAlpha(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2798,14 +2986,11 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
           break;
         for (x=0; x < (ssize_t) roi->width; x++)
         {
-          SetPixelRed(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelGreen(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
-          SetPixelBlue(image,ClampToQuantum((MagickRealType) QuantumRange*
-            (*p)),q);
+          SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
           p++;
           q+=GetPixelChannels(image);
         }
@@ -2832,46 +3017,39 @@ static void ImportFloatPixel(Image *image,const RectangleInfo *roi,
           case RedQuantum:
           case CyanQuantum:
           {
-            SetPixelRed(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelRed(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case GreenQuantum:
           case MagentaQuantum:
           {
-            SetPixelGreen(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelGreen(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case BlueQuantum:
           case YellowQuantum:
           {
-            SetPixelBlue(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelBlue(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case AlphaQuantum:
           {
-            SetPixelAlpha(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelAlpha(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case OpacityQuantum:
           {
-            SetPixelAlpha(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelAlpha(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case BlackQuantum:
           {
-            SetPixelBlack(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelBlack(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           case IndexQuantum:
           {
-            SetPixelGray(image,ClampToQuantum((MagickRealType)
-              QuantumRange*(*p)),q);
+            SetPixelGray(image,ClampToQuantum(QuantumRange*(*p)),q);
             break;
           }
           default:
@@ -2967,7 +3145,6 @@ static void ImportLongPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         q=GetAuthenticPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -3188,7 +3365,6 @@ static void ImportLongLongPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         q=GetAuthenticPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -3409,7 +3585,6 @@ static void ImportQuantumPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         q=GetAuthenticPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -3630,7 +3805,6 @@ static void ImportShortPixel(Image *image,const RectangleInfo *roi,
     }
   if (LocaleCompare(map,"I") == 0)
     {
-      (void) SetImageColorspace(image,GRAYColorspace,exception);
       for (y=0; y < (ssize_t) roi->height; y++)
       {
         q=GetAuthenticPixels(image,roi->x,roi->y+y,roi->width,1,exception);
@@ -3806,7 +3980,7 @@ MagickExport MagickBooleanType ImportImagePixels(Image *image,const ssize_t x,
       case 'A':
       {
         quantum_map[i]=AlphaQuantum;
-        image->matte=MagickTrue;
+        image->alpha_trait=BlendPixelTrait;
         break;
       }
       case 'B':
@@ -3839,6 +4013,7 @@ MagickExport MagickBooleanType ImportImagePixels(Image *image,const ssize_t x,
       case 'i':
       {
         quantum_map[i]=IndexQuantum;
+        (void) SetImageColorspace(image,GRAYColorspace,exception);
         break;
       }
       case 'm':
@@ -3852,7 +4027,7 @@ MagickExport MagickBooleanType ImportImagePixels(Image *image,const ssize_t x,
       case 'o':
       {
         quantum_map[i]=OpacityQuantum;
-        image->matte=MagickTrue;
+        image->alpha_trait=BlendPixelTrait;
         break;
       }
       case 'P':
@@ -3878,7 +4053,7 @@ MagickExport MagickBooleanType ImportImagePixels(Image *image,const ssize_t x,
       {
         quantum_map=(QuantumType *) RelinquishMagickMemory(quantum_map);
         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-          "UnrecognizedPixelMap","'%s'",map);
+          "UnrecognizedPixelMap","`%s'",map);
         return(MagickFalse);
       }
     }
@@ -3933,7 +4108,7 @@ MagickExport MagickBooleanType ImportImagePixels(Image *image,const ssize_t x,
     {
       quantum_map=(QuantumType *) RelinquishMagickMemory(quantum_map);
       (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-        "UnrecognizedPixelMap","'%s'",map);
+        "UnrecognizedStorageType","`%d'",type);
       break;
     }
   }
@@ -3963,6 +4138,101 @@ MagickExport MagickBooleanType ImportImagePixels(Image *image,const ssize_t x,
 %    o image: the image.
 %
 */
+
+static void LogPixelChannels(const Image *image)
+{
+  register ssize_t
+    i;
+
+  (void) LogMagickEvent(PixelEvent,GetMagickModule(),"%s[%.20g]",
+    image->filename,(double) image->number_channels);
+  for (i=0; i < (ssize_t) image->number_channels; i++)
+  {
+    char
+      traits[MaxTextExtent];
+
+    const char
+      *name;
+
+    PixelChannel
+      channel;
+
+    switch (GetPixelChannelChannel(image,i))
+    {
+      case RedPixelChannel:
+      {
+        name="red";
+        if (image->colorspace == CMYKColorspace)
+          name="cyan";
+        if (image->colorspace == GRAYColorspace)
+          name="gray";
+        break;
+      }
+      case GreenPixelChannel:
+      {
+        name="green";
+        if (image->colorspace == CMYKColorspace)
+          name="magenta";
+        break;
+      }
+      case BluePixelChannel:
+      {
+        name="blue";
+        if (image->colorspace == CMYKColorspace)
+          name="yellow";
+        break;
+      }
+      case BlackPixelChannel:
+      {
+        name="black";
+        if (image->storage_class == PseudoClass)
+          name="index";
+        break;
+      }
+      case IndexPixelChannel:
+      {
+        name="index";
+        break;
+      }
+      case AlphaPixelChannel:
+      {
+        name="alpha";
+        break;
+      }
+      case ReadMaskPixelChannel:
+      {
+        name="read-mask";
+        break;
+      }
+      case WriteMaskPixelChannel:
+      {
+        name="write-mask";
+        break;
+      }
+      case MetaPixelChannel:
+      {
+        name="meta";
+        break;
+      }
+      default:
+        name="undefined";
+    }
+    channel=GetPixelChannelChannel(image,i);
+    *traits='\0';
+    if ((GetPixelChannelTraits(image,channel) & UpdatePixelTrait) != 0)
+      (void) ConcatenateMagickString(traits,"update,",MaxTextExtent);
+    if ((GetPixelChannelTraits(image,channel) & BlendPixelTrait) != 0)
+      (void) ConcatenateMagickString(traits,"blend,",MaxTextExtent);
+    if ((GetPixelChannelTraits(image,channel) & CopyPixelTrait) != 0)
+      (void) ConcatenateMagickString(traits,"copy,",MaxTextExtent);
+    if (*traits == '\0')
+      (void) ConcatenateMagickString(traits,"undefined,",MaxTextExtent);
+    traits[strlen(traits)-1]='\0';
+    (void) LogMagickEvent(PixelEvent,GetMagickModule(),"  %.20g: %s (%s)",
+      (double) i,name,traits);
+  }
+}
+
 MagickExport void InitializePixelChannelMap(Image *image)
 {
   PixelTrait
@@ -3979,37 +4249,39 @@ MagickExport void InitializePixelChannelMap(Image *image)
   (void) ResetMagickMemory(image->channel_map,0,MaxPixelChannels*
     sizeof(*image->channel_map));
   trait=UpdatePixelTrait;
-  if (image->matte != MagickFalse)
+  if (image->alpha_trait == BlendPixelTrait)
     trait=(PixelTrait) (trait | BlendPixelTrait);
   n=0;
   if (image->colorspace == GRAYColorspace)
     {
-      SetPixelChannelMap(image,BluePixelChannel,trait,n);
-      SetPixelChannelMap(image,GreenPixelChannel,trait,n);
-      SetPixelChannelMap(image,RedPixelChannel,trait,n++);
+      SetPixelChannelAttributes(image,BluePixelChannel,trait,n);
+      SetPixelChannelAttributes(image,GreenPixelChannel,trait,n);
+      SetPixelChannelAttributes(image,RedPixelChannel,trait,n++);
     }
   else
     {
-      SetPixelChannelMap(image,RedPixelChannel,trait,n++);
-      SetPixelChannelMap(image,GreenPixelChannel,trait,n++);
-      SetPixelChannelMap(image,BluePixelChannel,trait,n++);
+      SetPixelChannelAttributes(image,RedPixelChannel,trait,n++);
+      SetPixelChannelAttributes(image,GreenPixelChannel,trait,n++);
+      SetPixelChannelAttributes(image,BluePixelChannel,trait,n++);
     }
   if (image->colorspace == CMYKColorspace)
-    SetPixelChannelMap(image,BlackPixelChannel,trait,n++);
-  if (image->matte != MagickFalse)
-    SetPixelChannelMap(image,AlphaPixelChannel,CopyPixelTrait,n++);
+    SetPixelChannelAttributes(image,BlackPixelChannel,trait,n++);
+  if (image->alpha_trait != UndefinedPixelTrait)
+    SetPixelChannelAttributes(image,AlphaPixelChannel,CopyPixelTrait,n++);
   if (image->storage_class == PseudoClass)
-    SetPixelChannelMap(image,IndexPixelChannel,CopyPixelTrait,n++);
-  if (image->mask != MagickFalse)
-    SetPixelChannelMap(image,MaskPixelChannel,CopyPixelTrait,n++);
+    SetPixelChannelAttributes(image,IndexPixelChannel,CopyPixelTrait,n++);
+  if (image->read_mask != MagickFalse)
+    SetPixelChannelAttributes(image,ReadMaskPixelChannel,CopyPixelTrait,n++);
+  if (image->write_mask != MagickFalse)
+    SetPixelChannelAttributes(image,WriteMaskPixelChannel,CopyPixelTrait,n++);
   assert((n+image->number_meta_channels) < MaxPixelChannels);
   for (i=0; i < (ssize_t) image->number_meta_channels; i++)
-    SetPixelChannelMap(image,(PixelChannel) (MetaPixelChannel+i),CopyPixelTrait,
-      n++);
+    SetPixelChannelAttributes(image,(PixelChannel) (MetaPixelChannel+i),
+      CopyPixelTrait,n++);
   image->number_channels=(size_t) n;
   if (image->debug != MagickFalse)
     LogPixelChannels(image);
-  (void) SetPixelChannelMask(image,image->channel_mask);
+  (void) SetImageChannelMask(image,image->channel_mask);
 }
 \f
 /*
@@ -4054,59 +4326,48 @@ MagickExport void InitializePixelChannelMap(Image *image)
 %
 */
 
-static inline double MagickMax(const MagickRealType x,const MagickRealType y)
+static inline void CatromWeights(const double x,double (*weights)[4])
 {
-  if (x > y)
-    return(x);
-  return(y);
-}
-
-static inline void CatromWeights(const MagickRealType x,
-  MagickRealType (*weights)[4])
-{
-  /*
-    Nicolas Robidoux' 10 flops (4* + 5- + 1+) refactoring of the
-    computation of the standard four 1D Catmull-Rom weights. The
-    sampling location is assumed between the second and third input
-    pixel locations, and x is the position relative to the second
-    input pixel location. Formulas originally derived for the VIPS
-    (Virtual Image Processing System) library.
-  */
-  MagickRealType
+  double
     alpha,
     beta,
     gamma;
 
-  alpha=(MagickRealType) 1.0-x;
-  beta=(MagickRealType) (-0.5)*x*alpha;
+  /*
+    Nicolas Robidoux' 10 flops (4* + 5- + 1+) refactoring of the computation
+    of the standard four 1D Catmull-Rom weights. The sampling location is
+    assumed between the second and third input pixel locations, and x is the
+    position relative to the second input pixel location. Formulas originally
+    derived for the VIPS (Virtual Image Processing System) library.
+  */
+  alpha=(double) 1.0-x;
+  beta=(double) (-0.5)*x*alpha;
   (*weights)[0]=alpha*beta;
   (*weights)[3]=x*beta;
   /*
-    The following computation of the inner weights from the outer ones
-    works for all Keys cubics.
+    The following computation of the inner weights from the outer ones work
+    for all Keys cubics.
   */
   gamma=(*weights)[3]-(*weights)[0];
   (*weights)[1]=alpha-(*weights)[0]+gamma;
   (*weights)[2]=x-(*weights)[3]-gamma;
 }
 
-static inline void SplineWeights(const MagickRealType x,
-  MagickRealType (*weights)[4])
+static inline void SplineWeights(const double x,double (*weights)[4])
 {
-  /*
-    Nicolas Robidoux' 12 flops (6* + 5- + 1+) refactoring of the
-    computation of the standard four 1D cubic B-spline smoothing
-    weights. The sampling location is assumed between the second and
-    third input pixel locations, and x is the position relative to the
-    second input pixel location.
-  */
-  MagickRealType
+  double
     alpha,
     beta;
 
-  alpha=(MagickRealType) 1.0-x;
-  (*weights)[3]=(MagickRealType) (1.0/6.0)*x*x*x;
-  (*weights)[0]=(MagickRealType) (1.0/6.0)*alpha*alpha*alpha;
+  /*
+    Nicolas Robidoux' 12 flops (6* + 5- + 1+) refactoring of the computation
+    of the standard four 1D cubic B-spline smoothing weights. The sampling
+    location is assumed between the second and third input pixel locations,
+    and x is the position relative to the second input pixel location.
+  */
+  alpha=(double) 1.0-x;
+  (*weights)[3]=(double) (1.0/6.0)*x*x*x;
+  (*weights)[0]=(double) (1.0/6.0)*alpha*alpha*alpha;
   beta=(*weights)[3]-(*weights)[0];
   (*weights)[1]=alpha-(*weights)[0]+beta;
   (*weights)[2]=x-(*weights)[3]-beta;
@@ -4119,7 +4380,7 @@ static inline double MeshInterpolate(const PointInfo *delta,const double p,
 }
 
 /*
-static inline ssize_t NearestNeighbor(const MagickRealType x)
+static inline ssize_t NearestNeighbor(const double x)
 {
   if (x >= 0.0)
     return((ssize_t) (x+0.5));
@@ -4132,14 +4393,17 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
   const PixelInterpolateMethod method,const double x,const double y,
   double *pixel,ExceptionInfo *exception)
 {
-  MagickBooleanType
-    status;
-
-  MagickRealType
+  double
     alpha[16],
     gamma,
     pixels[16];
 
+  MagickBooleanType
+    status;
+
+  PixelInterpolateMethod
+    interpolate;
+
   PixelTrait
     traits;
 
@@ -4153,16 +4417,13 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
     x_offset,
     y_offset;
 
-  PixelInterpolateMethod
-    interpolate;
-
   assert(image != (Image *) NULL);
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
   assert(image_view != (CacheView *) NULL);
   status=MagickTrue;
   *pixel=0.0;
-  traits=GetPixelChannelMapTraits(image,channel);
+  traits=GetPixelChannelTraits(image,channel);
   x_offset=(ssize_t) floor(x);
   y_offset=(ssize_t) floor(y);
   interpolate = method;
@@ -4170,39 +4431,40 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
     interpolate = image->interpolate;
   switch (interpolate)
   {
-    case AverageInterpolatePixel:        /* nearest 4 neighbours */
-    case Average9InterpolatePixel:       /* nearest 9 neighbours */
-    case Average16InterpolatePixel:      /* nearest 16 neighbours */
+    case AverageInterpolatePixel:  /* nearest 4 neighbours */
+    case Average9InterpolatePixel:  /* nearest 9 neighbours */
+    case Average16InterpolatePixel:  /* nearest 16 neighbours */
     {
-      size_t
-        count=2; /* size of the area to average - default nearest 4 */
+      ssize_t
+        count;
 
+      count=2;  /* size of the area to average - default nearest 4 */
       if (interpolate == Average9InterpolatePixel)
         {
           count=3;
           x_offset=(ssize_t) (floor(x+0.5)-1);
           y_offset=(ssize_t) (floor(y+0.5)-1);
         }
-      else if (interpolate == Average16InterpolatePixel)
-        {
-          count=4;
-          x_offset--;
-          y_offset--;
-        }
-      p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset,count,count,
-        exception);
+      else
+        if (interpolate == Average16InterpolatePixel)
+          {
+            count=4;
+            x_offset--;
+            y_offset--;
+          }
+      p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset,(size_t) count,
+        (size_t) count,exception);
       if (p == (const Quantum *) NULL)
         {
           status=MagickFalse;
           break;
         }
-
-      count*=count;   /* Number of pixels to Average */
+      count*=count;  /* Number of pixels to average */
       if ((traits & BlendPixelTrait) == 0)
         for (i=0; i < (ssize_t) count; i++)
         {
           alpha[i]=1.0;
-          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
+          pixels[i]=(double) p[i*GetPixelChannels(image)+channel];
         }
       else
         for (i=0; i < (ssize_t) count; i++)
@@ -4213,7 +4475,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
         }
       for (i=0; i < (ssize_t) count; i++)
       {
-        gamma=MagickEpsilonReciprocal(alpha[i])/count;
+        gamma=PerceptibleReciprocal(alpha[i])/count;
         *pixel+=gamma*pixels[i];
       }
       break;
@@ -4235,7 +4497,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
         for (i=0; i < 4; i++)
         {
           alpha[i]=1.0;
-          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
+          pixels[i]=(double) p[i*GetPixelChannels(image)+channel];
         }
       else
         for (i=0; i < 4; i++)
@@ -4250,14 +4512,69 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
       epsilon.y=1.0-delta.y;
       gamma=((epsilon.y*(epsilon.x*alpha[0]+delta.x*alpha[1])+delta.y*
         (epsilon.x*alpha[2]+delta.x*alpha[3])));
-      gamma=MagickEpsilonReciprocal(gamma);
+      gamma=PerceptibleReciprocal(gamma);
       *pixel=gamma*(epsilon.y*(epsilon.x*pixels[0]+delta.x*pixels[1])+delta.y*
         (epsilon.x*pixels[2]+delta.x*pixels[3]));
       break;
     }
+    case BlendInterpolatePixel:
+    {
+      p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset,2,2,exception);
+      if (p == (const Quantum *) NULL)
+        {
+          status=MagickFalse;
+          break;
+        }
+      if ((traits & BlendPixelTrait) == 0)
+        for (i=0; i < 4; i++)
+        {
+          alpha[i]=1.0;
+          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
+        }
+      else
+        for (i=0; i < 4; i++)
+        {
+          alpha[i]=QuantumScale*GetPixelAlpha(image,p+i*
+            GetPixelChannels(image));
+          pixels[i]=alpha[i]*p[i*GetPixelChannels(image)+channel];
+        }
+      gamma=1.0;  /* number of pixels blended together (its variable) */
+      for (i=0; i <= 1L; i++) {
+        if ((y-y_offset) >= 0.75)
+          {
+            alpha[i]=alpha[i+2];  /* take right pixels */
+            pixels[i]=pixels[i+2];
+          }
+        else
+          if ((y-y_offset) > 0.25)
+            {
+              gamma=2.0;  /* blend both pixels in row */
+              alpha[i]+=alpha[i+2];  /* add up alpha weights */
+              pixels[i]+=pixels[i+2];
+            }
+      }
+      if ((x-x_offset) >= 0.75)
+        {
+          alpha[0]=alpha[1];  /* take bottom row blend */
+          pixels[0]=pixels[1];
+        }
+      else
+        if ((x-x_offset) > 0.25)
+          {
+            gamma*=2.0;  /* blend both rows */
+            alpha[0]+=alpha[1];  /* add up alpha weights */
+            pixels[0]+=pixels[1];
+          }
+      if (channel != AlphaPixelChannel)
+        gamma=PerceptibleReciprocal(alpha[0]);  /* (color) 1/alpha_weights */
+      else
+        gamma=PerceptibleReciprocal(gamma);  /* (alpha) 1/number_of_pixels */
+      *pixel=gamma*pixels[0];
+      break;
+    }
     case CatromInterpolatePixel:
     {
-      MagickRealType
+      double
         cx[4],
         cy[4];
 
@@ -4272,7 +4589,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
         for (i=0; i < 16; i++)
         {
           alpha[i]=1.0;
-          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
+          pixels[i]=(double) p[i*GetPixelChannels(image)+channel];
         }
       else
         for (i=0; i < 16; i++)
@@ -4281,10 +4598,10 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
             GetPixelChannels(image));
           pixels[i]=alpha[i]*p[i*GetPixelChannels(image)+channel];
         }
-      CatromWeights((MagickRealType) (x-x_offset),&cx);
-      CatromWeights((MagickRealType) (y-y_offset),&cy);
-      gamma=(channel == AlphaPixelChannel ? (MagickRealType) 1.0 :
-        MagickEpsilonReciprocal(cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*
+      CatromWeights((double) (x-x_offset),&cx);
+      CatromWeights((double) (y-y_offset),&cy);
+      gamma=(channel == AlphaPixelChannel ? (double) 1.0 :
+        PerceptibleReciprocal(cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*
         alpha[2]+cx[3]*alpha[3])+cy[1]*(cx[0]*alpha[4]+cx[1]*alpha[5]+cx[2]*
         alpha[6]+cx[3]*alpha[7])+cy[2]*(cx[0]*alpha[8]+cx[1]*alpha[9]+cx[2]*
         alpha[10]+cx[3]*alpha[11])+cy[3]*(cx[0]*alpha[12]+cx[1]*alpha[13]+
@@ -4296,45 +4613,6 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
         pixels[13]+cx[2]*pixels[14]+cx[3]*pixels[15]));
       break;
     }
-#if 0
-    /* deprecated useless and very slow interpolator */
-    case FilterInterpolatePixel:
-    {
-      CacheView
-        *filter_view;
-
-      Image
-        *excerpt_image,
-        *filter_image;
-
-      RectangleInfo
-        geometry;
-
-      geometry.width=4L;
-      geometry.height=4L;
-      geometry.x=x_offset-1;
-      geometry.y=y_offset-1;
-      excerpt_image=ExcerptImage(image,&geometry,exception);
-      if (excerpt_image == (Image *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      filter_image=ResizeImage(excerpt_image,1,1,image->filter,exception);
-      excerpt_image=DestroyImage(excerpt_image);
-      if (filter_image == (Image *) NULL)
-        break;
-      filter_view=AcquireVirtualCacheView(filter_image,exception);
-      p=GetCacheViewVirtualPixels(filter_view,0,0,1,1,exception);
-      if (p == (const Quantum *) NULL)
-        status=MagickFalse;
-      else
-        *pixel=(double) GetPixelChannel(image,channel,p);
-      filter_view=DestroyCacheView(filter_view);
-      filter_image=DestroyImage(filter_image);
-      break;
-    }
-#endif
     case IntegerInterpolatePixel:
     {
       p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset,1,1,exception);
@@ -4375,7 +4653,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
         for (i=0; i < 4; i++)
         {
           alpha[i]=1.0;
-          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
+          pixels[i]=(double) p[i*GetPixelChannels(image)+channel];
         }
       else
         for (i=0; i < 4; i++)
@@ -4386,10 +4664,10 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
         }
       delta.x=x-x_offset;
       delta.y=y-y_offset;
-      luminance.x=GetPixelLuminance(image,p)-(double)
-        GetPixelLuminance(image,p+3*GetPixelChannels(image));
-      luminance.y=GetPixelLuminance(image,p+GetPixelChannels(image))-(double)
-        GetPixelLuminance(image,p+2*GetPixelChannels(image));
+      luminance.x=GetPixelLuma(image,p)-(double)
+        GetPixelLuma(image,p+3*GetPixelChannels(image));
+      luminance.y=GetPixelLuma(image,p+GetPixelChannels(image))-(double)
+        GetPixelLuma(image,p+2*GetPixelChannels(image));
       if (fabs(luminance.x) < fabs(luminance.y))
         {
           /*
@@ -4402,7 +4680,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
               */
               delta.y=1.0-delta.y;
               gamma=MeshInterpolate(&delta,alpha[2],alpha[3],alpha[0]);
-              gamma=MagickEpsilonReciprocal(gamma);
+              gamma=PerceptibleReciprocal(gamma);
               *pixel=gamma*MeshInterpolate(&delta,pixels[2],pixels[3],
                 pixels[0]);
             }
@@ -4413,7 +4691,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
               */
               delta.x=1.0-delta.x;
               gamma=MeshInterpolate(&delta,alpha[1],alpha[0],alpha[3]);
-              gamma=MagickEpsilonReciprocal(gamma);
+              gamma=PerceptibleReciprocal(gamma);
               *pixel=gamma*MeshInterpolate(&delta,pixels[1],pixels[0],
                 pixels[3]);
             }
@@ -4429,7 +4707,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
                 Top-left triangle (pixel: 0, diagonal: 1-2).
               */
               gamma=MeshInterpolate(&delta,alpha[0],alpha[1],alpha[2]);
-              gamma=MagickEpsilonReciprocal(gamma);
+              gamma=PerceptibleReciprocal(gamma);
               *pixel=gamma*MeshInterpolate(&delta,pixels[0],pixels[1],
                 pixels[2]);
             }
@@ -4441,7 +4719,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
               delta.x=1.0-delta.x;
               delta.y=1.0-delta.y;
               gamma=MeshInterpolate(&delta,alpha[3],alpha[2],alpha[1]);
-              gamma=MagickEpsilonReciprocal(gamma);
+              gamma=PerceptibleReciprocal(gamma);
               *pixel=gamma*MeshInterpolate(&delta,pixels[3],pixels[2],
                 pixels[1]);
             }
@@ -4450,7 +4728,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
     }
     case SplineInterpolatePixel:
     {
-      MagickRealType
+      double
         cx[4],
         cy[4];
 
@@ -4465,7 +4743,7 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
         for (i=0; i < 16; i++)
         {
           alpha[i]=1.0;
-          pixels[i]=(MagickRealType) p[i*GetPixelChannels(image)+channel];
+          pixels[i]=(double) p[i*GetPixelChannels(image)+channel];
         }
       else
         for (i=0; i < 16; i++)
@@ -4474,10 +4752,10 @@ MagickExport MagickBooleanType InterpolatePixelChannel(const Image *image,
             GetPixelChannels(image));
           pixels[i]=alpha[i]*p[i*GetPixelChannels(image)+channel];
         }
-      SplineWeights((MagickRealType) (x-x_offset),&cx);
-      SplineWeights((MagickRealType) (y-y_offset),&cy);
-      gamma=(channel == AlphaPixelChannel ? (MagickRealType) 1.0 :
-        MagickEpsilonReciprocal(cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*
+      SplineWeights((double) (x-x_offset),&cx);
+      SplineWeights((double) (y-y_offset),&cy);
+      gamma=(channel == AlphaPixelChannel ? (double) 1.0 :
+        PerceptibleReciprocal(cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*
         alpha[2]+cx[3]*alpha[3])+cy[1]*(cx[0]*alpha[4]+cx[1]*alpha[5]+cx[2]*
         alpha[6]+cx[3]*alpha[7])+cy[2]*(cx[0]*alpha[8]+cx[1]*alpha[9]+cx[2]*
         alpha[10]+cx[3]*alpha[11])+cy[3]*(cx[0]*alpha[12]+cx[1]*alpha[13]+
@@ -4543,18 +4821,11 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
   MagickBooleanType
     status;
 
-  MagickRealType
+  double
     alpha[16],
     gamma,
     pixels[16];
 
-  PixelChannel
-    channel;
-
-  PixelTrait
-    destination_traits,
-    traits;
-
   register const Quantum
     *p;
 
@@ -4580,33 +4851,35 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
     interpolate = source->interpolate;
   switch (interpolate)
   {
-    case AverageInterpolatePixel:        /* nearest 4 neighbours */
-    case Average9InterpolatePixel:       /* nearest 9 neighbours */
-    case Average16InterpolatePixel:      /* nearest 16 neighbours */
+    case AverageInterpolatePixel:  /* nearest 4 neighbours */
+    case Average9InterpolatePixel:  /* nearest 9 neighbours */
+    case Average16InterpolatePixel:  /* nearest 16 neighbours */
     {
-      size_t
-        count=2; /* size of the area to average - default nearest 4 */
+      ssize_t
+        count;
 
+      count=2;  /* size of the area to average - default nearest 4 */
       if (interpolate == Average9InterpolatePixel)
         {
           count=3;
           x_offset=(ssize_t) (floor(x+0.5)-1);
           y_offset=(ssize_t) (floor(y+0.5)-1);
         }
-      else if (interpolate == Average16InterpolatePixel)
-        {
-          count=4;
-          x_offset--;
-          y_offset--;
-        }
-      p=GetCacheViewVirtualPixels(source_view,x_offset,y_offset,count,count,
-        exception);
+      else
+        if (interpolate == Average16InterpolatePixel)
+          {
+            count=4;
+            x_offset--;
+            y_offset--;
+          }
+      p=GetCacheViewVirtualPixels(source_view,x_offset,y_offset,(size_t) count,
+        (size_t) count,exception);
       if (p == (const Quantum *) NULL)
         {
           status=MagickFalse;
           break;
         }
-      count*=count;   /* Number of pixels to Average */
+      count*=count;  /* Number of pixels to average */
       for (i=0; i < (ssize_t) GetPixelChannels(source); i++)
       {
         double
@@ -4615,14 +4888,15 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
         register ssize_t
           j;
 
-        channel=GetPixelChannelMapChannel(source,i);
-        traits=GetPixelChannelMapTraits(source,channel);
-        destination_traits=GetPixelChannelMapTraits(destination,channel);
+        PixelChannel channel=GetPixelChannelChannel(source,i);
+        PixelTrait traits=GetPixelChannelTraits(source,channel);
+        PixelTrait destination_traits=GetPixelChannelTraits(destination,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
             (destination_traits == UndefinedPixelTrait))
           continue;
         for (j=0; j < (ssize_t) count; j++)
-          pixels[j]=(MagickRealType) p[j*GetPixelChannels(source)+i];
+          pixels[j]=(double) p[j*GetPixelChannels(source)+i];
         sum=0.0;
         if ((traits & BlendPixelTrait) == 0)
           {
@@ -4637,7 +4911,7 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
           alpha[j]=QuantumScale*GetPixelAlpha(source,p+j*
             GetPixelChannels(source));
           pixels[j]*=alpha[j];
-          gamma=MagickEpsilonReciprocal(alpha[j]);
+          gamma=PerceptibleReciprocal(alpha[j]);
           sum+=gamma*pixels[j];
         }
         sum/=count;
@@ -4660,9 +4934,10 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
           delta,
           epsilon;
 
-        channel=GetPixelChannelMapChannel(source,i);
-        traits=GetPixelChannelMapTraits(source,channel);
-        destination_traits=GetPixelChannelMapTraits(destination,channel);
+        PixelChannel channel=GetPixelChannelChannel(source,i);
+        PixelTrait traits=GetPixelChannelTraits(source,channel);
+        PixelTrait destination_traits=GetPixelChannelTraits(destination,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
             (destination_traits == UndefinedPixelTrait))
           continue;
@@ -4670,14 +4945,14 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
         delta.y=y-y_offset;
         epsilon.x=1.0-delta.x;
         epsilon.y=1.0-delta.y;
-        pixels[0]=(MagickRealType) p[i];
-        pixels[1]=(MagickRealType) p[GetPixelChannels(source)+i];
-        pixels[2]=(MagickRealType) p[2*GetPixelChannels(source)+i];
-        pixels[3]=(MagickRealType) p[3*GetPixelChannels(source)+i];
+        pixels[0]=(double) p[i];
+        pixels[1]=(double) p[GetPixelChannels(source)+i];
+        pixels[2]=(double) p[2*GetPixelChannels(source)+i];
+        pixels[3]=(double) p[3*GetPixelChannels(source)+i];
         if ((traits & BlendPixelTrait) == 0)
           {
             gamma=((epsilon.y*(epsilon.x+delta.x)+delta.y*(epsilon.x+delta.x)));
-            gamma=MagickEpsilonReciprocal(gamma);
+            gamma=PerceptibleReciprocal(gamma);
             SetPixelChannel(destination,channel,ClampToQuantum(gamma*(epsilon.y*
               (epsilon.x*pixels[0]+delta.x*pixels[1])+delta.y*(epsilon.x*
               pixels[2]+delta.x*pixels[3]))),pixel);
@@ -4695,16 +4970,86 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
         pixels[3]*=alpha[3];
         gamma=((epsilon.y*(epsilon.x*alpha[0]+delta.x*alpha[1])+delta.y*
           (epsilon.x*alpha[2]+delta.x*alpha[3])));
-        gamma=MagickEpsilonReciprocal(gamma);
+        gamma=PerceptibleReciprocal(gamma);
         SetPixelChannel(destination,channel,ClampToQuantum(gamma*(epsilon.y*
           (epsilon.x*pixels[0]+delta.x*pixels[1])+delta.y*(epsilon.x*pixels[2]+
           delta.x*pixels[3]))),pixel);
       }
       break;
     }
+    case BlendInterpolatePixel:
+    {
+      p=GetCacheViewVirtualPixels(source_view,x_offset,y_offset,2,2,exception);
+      if (p == (const Quantum *) NULL)
+        {
+          status=MagickFalse;
+          break;
+        }
+      for (i=0; i < (ssize_t) GetPixelChannels(source); i++)
+      {
+        register ssize_t
+          j;
+
+        PixelChannel channel=GetPixelChannelChannel(source,i);
+        PixelTrait traits=GetPixelChannelTraits(source,channel);
+        PixelTrait destination_traits=GetPixelChannelTraits(destination,
+          channel);
+        if ((traits == UndefinedPixelTrait) ||
+            (destination_traits == UndefinedPixelTrait))
+          continue;
+        if ((traits & BlendPixelTrait) == 0)
+          for (j=0; j < 4; j++)
+          {
+            alpha[j]=1.0;
+            pixels[j]=(MagickRealType) p[j*GetPixelChannels(source)+channel];
+          }
+        else
+          for (j=0; j < 4; j++)
+          {
+            alpha[j]=QuantumScale*GetPixelAlpha(source,p+j*
+              GetPixelChannels(source));
+            pixels[j]=alpha[j]*p[j*GetPixelChannels(source)+channel];
+          }
+        gamma=1.0;  /* number of pixels blended together (its variable) */
+        for (j=0; j <= 1L; j++)
+        {
+          if ((y-y_offset) >= 0.75)
+            {
+              alpha[j]=alpha[j+2];  /* take right pixels */
+              pixels[j]=pixels[j+2];
+            }
+          else
+            if ((y-y_offset) > 0.25)
+              {
+                gamma=2.0;  /* blend both pixels in row */
+                alpha[j]+=alpha[j+2];  /* add up alpha weights */
+                pixels[j]+=pixels[j+2];
+              }
+        }
+        if ((x-x_offset) >= 0.75)
+          {
+            alpha[0]=alpha[1];  /* take bottom row blend */
+            pixels[0]=pixels[1];
+          }
+        else
+           if ((x-x_offset) > 0.25)
+             {
+               gamma*=2.0;  /* blend both rows */
+               alpha[0]+=alpha[1];  /* add up alpha weights */
+               pixels[0]+=pixels[1];
+             }
+        if ((traits & BlendPixelTrait) == 0)
+          gamma=PerceptibleReciprocal(alpha[0]);  /* (color) 1/alpha_weights */
+        else
+          gamma=PerceptibleReciprocal(gamma);  /* (alpha) 1/number_of_pixels */
+        SetPixelChannel(destination,channel,ClampToQuantum(gamma*pixels[0]),
+          pixel);
+      }
+      break;
+    }
     case CatromInterpolatePixel:
     {
-      MagickRealType
+      double
         cx[4],
         cy[4];
 
@@ -4720,9 +5065,10 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
         register ssize_t
           j;
 
-        channel=GetPixelChannelMapChannel(source,i);
-        traits=GetPixelChannelMapTraits(source,channel);
-        destination_traits=GetPixelChannelMapTraits(destination,channel);
+        PixelChannel channel=GetPixelChannelChannel(source,i);
+        PixelTrait traits=GetPixelChannelTraits(source,channel);
+        PixelTrait destination_traits=GetPixelChannelTraits(destination,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
             (destination_traits == UndefinedPixelTrait))
           continue;
@@ -4730,7 +5076,7 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
           for (j=0; j < 16; j++)
           {
             alpha[j]=1.0;
-            pixels[j]=(MagickRealType) p[j*GetPixelChannels(source)+i];
+            pixels[j]=(double) p[j*GetPixelChannels(source)+i];
           }
         else
           for (j=0; j < 16; j++)
@@ -4739,10 +5085,10 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
               GetPixelChannels(source));
             pixels[j]=alpha[j]*p[j*GetPixelChannels(source)+i];
           }
-        CatromWeights((MagickRealType) (x-x_offset),&cx);
-        CatromWeights((MagickRealType) (y-y_offset),&cy);
-        gamma=((traits & BlendPixelTrait) ? (MagickRealType) (1.0) :
-          MagickEpsilonReciprocal(cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*
+        CatromWeights((double) (x-x_offset),&cx);
+        CatromWeights((double) (y-y_offset),&cy);
+        gamma=((traits & BlendPixelTrait) ? (double) (1.0) :
+          PerceptibleReciprocal(cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*
           alpha[2]+cx[3]*alpha[3])+cy[1]*(cx[0]*alpha[4]+cx[1]*alpha[5]+cx[2]*
           alpha[6]+cx[3]*alpha[7])+cy[2]*(cx[0]*alpha[8]+cx[1]*alpha[9]+cx[2]*
           alpha[10]+cx[3]*alpha[11])+cy[3]*(cx[0]*alpha[12]+cx[1]*alpha[13]+
@@ -4756,56 +5102,6 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
       }
       break;
     }
-#if 0
-    /* deprecated useless and very slow interpolator */
-    case FilterInterpolatePixel:
-    {
-      for (i=0; i < (ssize_t) GetPixelChannels(source); i++)
-      {
-        CacheView
-          *filter_view;
-
-        Image
-          *excerpt_source,
-          *filter_source;
-
-        RectangleInfo
-          geometry;
-
-        channel=GetPixelChannelMapChannel(source,i);
-        traits=GetPixelChannelMapTraits(source,channel);
-        destination_traits=GetPixelChannelMapTraits(destination,channel);
-        if ((traits == UndefinedPixelTrait) ||
-            (destination_traits == UndefinedPixelTrait))
-          continue;
-        geometry.width=4L;
-        geometry.height=4L;
-        geometry.x=x_offset-1;
-        geometry.y=y_offset-1;
-        excerpt_source=ExcerptImage(source,&geometry,exception);
-        if (excerpt_source == (Image *) NULL)
-          {
-            status=MagickFalse;
-            continue;
-          }
-        filter_source=ResizeImage(excerpt_source,1,1,source->filter,exception);
-        excerpt_source=DestroyImage(excerpt_source);
-        if (filter_source == (Image *) NULL)
-          continue;
-        filter_view=AcquireVirtualCacheView(filter_source,exception);
-        p=GetCacheViewVirtualPixels(filter_view,0,0,1,1,exception);
-        if (p == (const Quantum *) NULL)
-          status=MagickFalse;
-        else
-          {
-            SetPixelChannel(destination,channel,p[i],pixel);
-          }
-        filter_view=DestroyCacheView(filter_view);
-        filter_source=DestroyImage(filter_source);
-      }
-      break;
-    }
-#endif
     case IntegerInterpolatePixel:
     {
       p=GetCacheViewVirtualPixels(source_view,x_offset,y_offset,1,1,exception);
@@ -4816,9 +5112,10 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
         }
       for (i=0; i < (ssize_t) GetPixelChannels(source); i++)
       {
-        channel=GetPixelChannelMapChannel(source,i);
-        traits=GetPixelChannelMapTraits(source,channel);
-        destination_traits=GetPixelChannelMapTraits(destination,channel);
+        PixelChannel channel=GetPixelChannelChannel(source,i);
+        PixelTrait traits=GetPixelChannelTraits(source,channel);
+        PixelTrait destination_traits=GetPixelChannelTraits(destination,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
             (destination_traits == UndefinedPixelTrait))
           continue;
@@ -4838,9 +5135,10 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
         }
       for (i=0; i < (ssize_t) GetPixelChannels(source); i++)
       {
-        channel=GetPixelChannelMapChannel(source,i);
-        traits=GetPixelChannelMapTraits(source,channel);
-        destination_traits=GetPixelChannelMapTraits(destination,channel);
+        PixelChannel channel=GetPixelChannelChannel(source,i);
+        PixelTrait traits=GetPixelChannelTraits(source,channel);
+        PixelTrait destination_traits=GetPixelChannelTraits(destination,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
             (destination_traits == UndefinedPixelTrait))
           continue;
@@ -4862,16 +5160,17 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
           delta,
           luminance;
 
-        channel=GetPixelChannelMapChannel(source,i);
-        traits=GetPixelChannelMapTraits(source,channel);
-        destination_traits=GetPixelChannelMapTraits(destination,channel);
+        PixelChannel channel=GetPixelChannelChannel(source,i);
+        PixelTrait traits=GetPixelChannelTraits(source,channel);
+        PixelTrait destination_traits=GetPixelChannelTraits(destination,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
             (destination_traits == UndefinedPixelTrait))
           continue;
-        pixels[0]=(MagickRealType) p[i];
-        pixels[1]=(MagickRealType) p[GetPixelChannels(source)+i];
-        pixels[2]=(MagickRealType) p[2*GetPixelChannels(source)+i];
-        pixels[3]=(MagickRealType) p[3*GetPixelChannels(source)+i];
+        pixels[0]=(double) p[i];
+        pixels[1]=(double) p[GetPixelChannels(source)+i];
+        pixels[2]=(double) p[2*GetPixelChannels(source)+i];
+        pixels[3]=(double) p[3*GetPixelChannels(source)+i];
         if ((traits & BlendPixelTrait) == 0)
           {
             alpha[0]=1.0;
@@ -4891,12 +5190,11 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
           }
         delta.x=x-x_offset;
         delta.y=y-y_offset;
-        luminance.x=fabs((double)(
-              GetPixelLuminance(source,p)
-               -GetPixelLuminance(source,p+3*GetPixelChannels(source))));
-        luminance.y=fabs((double)(
-              GetPixelLuminance(source,p+GetPixelChannels(source))
-               -GetPixelLuminance(source,p+2*GetPixelChannels(source))));
+        luminance.x=fabs((double) (GetPixelLuma(source,p)-
+          GetPixelLuma(source,p+3*GetPixelChannels(source))));
+        luminance.y=fabs((double) (GetPixelLuma(source,p+
+          GetPixelChannels(source))-GetPixelLuma(source,p+2*
+          GetPixelChannels(source))));
         if (luminance.x < luminance.y)
           {
             /*
@@ -4909,7 +5207,7 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
                 */
                 delta.y=1.0-delta.y;
                 gamma=MeshInterpolate(&delta,alpha[2],alpha[3],alpha[0]);
-                gamma=MagickEpsilonReciprocal(gamma);
+                gamma=PerceptibleReciprocal(gamma);
                 SetPixelChannel(destination,channel,ClampToQuantum(gamma*
                   MeshInterpolate(&delta,pixels[2],pixels[3],pixels[0])),pixel);
               }
@@ -4920,7 +5218,7 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
                 */
                 delta.x=1.0-delta.x;
                 gamma=MeshInterpolate(&delta,alpha[1],alpha[0],alpha[3]);
-                gamma=MagickEpsilonReciprocal(gamma);
+                gamma=PerceptibleReciprocal(gamma);
                 SetPixelChannel(destination,channel,ClampToQuantum(gamma*
                   MeshInterpolate(&delta,pixels[1],pixels[0],pixels[3])),pixel);
               }
@@ -4936,7 +5234,7 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
                   Top-left triangle (pixel: 0, diagonal: 1-2).
                 */
                 gamma=MeshInterpolate(&delta,alpha[0],alpha[1],alpha[2]);
-                gamma=MagickEpsilonReciprocal(gamma);
+                gamma=PerceptibleReciprocal(gamma);
                 SetPixelChannel(destination,channel,ClampToQuantum(gamma*
                   MeshInterpolate(&delta,pixels[0],pixels[1],pixels[2])),pixel);
               }
@@ -4948,7 +5246,7 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
                 delta.x=1.0-delta.x;
                 delta.y=1.0-delta.y;
                 gamma=MeshInterpolate(&delta,alpha[3],alpha[2],alpha[1]);
-                gamma=MagickEpsilonReciprocal(gamma);
+                gamma=PerceptibleReciprocal(gamma);
                 SetPixelChannel(destination,channel,ClampToQuantum(gamma*
                   MeshInterpolate(&delta,pixels[3],pixels[2],pixels[1])),pixel);
               }
@@ -4958,7 +5256,7 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
     }
     case SplineInterpolatePixel:
     {
-      MagickRealType
+      double
         cx[4],
         cy[4];
 
@@ -4974,9 +5272,10 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
         register ssize_t
           j;
 
-        channel=GetPixelChannelMapChannel(source,i);
-        traits=GetPixelChannelMapTraits(source,channel);
-        destination_traits=GetPixelChannelMapTraits(destination,channel);
+        PixelChannel channel=GetPixelChannelChannel(source,i);
+        PixelTrait traits=GetPixelChannelTraits(source,channel);
+        PixelTrait destination_traits=GetPixelChannelTraits(destination,
+          channel);
         if ((traits == UndefinedPixelTrait) ||
             (destination_traits == UndefinedPixelTrait))
           continue;
@@ -4984,7 +5283,7 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
           for (j=0; j < 16; j++)
           {
             alpha[j]=1.0;
-            pixels[j]=(MagickRealType) p[j*GetPixelChannels(source)+i];
+            pixels[j]=(double) p[j*GetPixelChannels(source)+i];
           }
         else
           for (j=0; j < 16; j++)
@@ -4993,10 +5292,10 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
               GetPixelChannels(source));
             pixels[j]=alpha[j]*p[j*GetPixelChannels(source)+i];
           }
-        SplineWeights((MagickRealType) (x-x_offset),&cx);
-        SplineWeights((MagickRealType) (y-y_offset),&cy);
-        gamma=((traits & BlendPixelTrait) ? (MagickRealType) (1.0) :
-          MagickEpsilonReciprocal(cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*
+        SplineWeights((double) (x-x_offset),&cx);
+        SplineWeights((double) (y-y_offset),&cy);
+        gamma=((traits & BlendPixelTrait) ? (double) (1.0) :
+          PerceptibleReciprocal(cy[0]*(cx[0]*alpha[0]+cx[1]*alpha[1]+cx[2]*
           alpha[2]+cx[3]*alpha[3])+cy[1]*(cx[0]*alpha[4]+cx[1]*alpha[5]+cx[2]*
           alpha[6]+cx[3]*alpha[7])+cy[2]*(cx[0]*alpha[8]+cx[1]*alpha[9]+cx[2]*
           alpha[10]+cx[3]*alpha[11])+cy[3]*(cx[0]*alpha[12]+cx[1]*alpha[13]+
@@ -5055,18 +5354,18 @@ MagickExport MagickBooleanType InterpolatePixelChannels(const Image *source,
 */
 
 static inline void AlphaBlendPixelInfo(const Image *image,
-  const Quantum *pixel,PixelInfo *pixel_info,MagickRealType *alpha)
+  const Quantum *pixel,PixelInfo *pixel_info,double *alpha)
 {
-  if (image->matte == MagickFalse)
+  if (image->alpha_trait != BlendPixelTrait)
     {
       *alpha=1.0;
-      pixel_info->red=(MagickRealType) GetPixelRed(image,pixel);
-      pixel_info->green=(MagickRealType) GetPixelGreen(image,pixel);
-      pixel_info->blue=(MagickRealType) GetPixelBlue(image,pixel);
+      pixel_info->red=(double) GetPixelRed(image,pixel);
+      pixel_info->green=(double) GetPixelGreen(image,pixel);
+      pixel_info->blue=(double) GetPixelBlue(image,pixel);
       pixel_info->black=0.0;
       if (image->colorspace == CMYKColorspace)
-        pixel_info->black=(MagickRealType) GetPixelBlack(image,pixel);
-      pixel_info->alpha=(MagickRealType) GetPixelAlpha(image,pixel);
+        pixel_info->black=(double) GetPixelBlack(image,pixel);
+      pixel_info->alpha=(double) GetPixelAlpha(image,pixel);
       return;
     }
   *alpha=QuantumScale*GetPixelAlpha(image,pixel);
@@ -5076,7 +5375,7 @@ static inline void AlphaBlendPixelInfo(const Image *image,
   pixel_info->black=0.0;
   if (image->colorspace == CMYKColorspace)
     pixel_info->black=(*alpha*GetPixelBlack(image,pixel));
-  pixel_info->alpha=(MagickRealType) GetPixelAlpha(image,pixel);
+  pixel_info->alpha=(double) GetPixelAlpha(image,pixel);
 }
 
 MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
@@ -5086,7 +5385,7 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
   MagickBooleanType
     status;
 
-  MagickRealType
+  double
     alpha[16],
     gamma;
 
@@ -5117,13 +5416,14 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
     interpolate = image->interpolate;
   switch (interpolate)
   {
-    case AverageInterpolatePixel:        /* nearest 4 neighbours */
-    case Average9InterpolatePixel:       /* nearest 9 neighbours */
-    case Average16InterpolatePixel:      /* nearest 16 neighbours */
+    case AverageInterpolatePixel:  /* nearest 4 neighbours */
+    case Average9InterpolatePixel:  /* nearest 9 neighbours */
+    case Average16InterpolatePixel:  /* nearest 16 neighbours */
     {
-      size_t
-        count=2; /* size of the area to average - default nearest 4 */
+      ssize_t
+        count;
 
+      count=2;  /* size of the area to average - default nearest 4 */
       if (interpolate == Average9InterpolatePixel)
         {
           count=3;
@@ -5136,8 +5436,8 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
           x_offset--;
           y_offset--;
         }
-      p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset,count,count,
-        exception);
+      p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset,(size_t) count,
+        (size_t) count,exception);
       if (p == (const Quantum *) NULL)
         {
           status=MagickFalse;
@@ -5148,29 +5448,29 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
       pixel->blue=0.0;
       pixel->black=0.0;
       pixel->alpha=0.0;
-      count*=count;         /* number of pixels - square of size */
+      count*=count;  /* number of pixels - square of size */
       for (i=0; i < (ssize_t) count; i++)
       {
         AlphaBlendPixelInfo(image,p,pixels,alpha);
-        gamma=MagickEpsilonReciprocal(alpha[0]);
-        pixel->red   += gamma*pixels[0].red;
-        pixel->green += gamma*pixels[0].green;
-        pixel->blue  += gamma*pixels[0].blue;
-        pixel->black += gamma*pixels[0].black;
-        pixel->alpha +=       pixels[0].alpha;
+        gamma=PerceptibleReciprocal(alpha[0]);
+        pixel->red+=gamma*pixels[0].red;
+        pixel->green+=gamma*pixels[0].green;
+        pixel->blue+=gamma*pixels[0].blue;
+        pixel->black+=gamma*pixels[0].black;
+        pixel->alpha+=pixels[0].alpha;
         p += GetPixelChannels(image);
       }
       gamma=1.0/count;   /* average weighting of each pixel in area */
-      pixel->red   *= gamma;
-      pixel->green *= gamma;
-      pixel->blue  *= gamma;
-      pixel->black *= gamma;
-      pixel->alpha *= gamma;
+      pixel->red*=gamma;
+      pixel->green*=gamma;
+      pixel->blue*=gamma;
+      pixel->black*=gamma;
+      pixel->alpha*=gamma;
       break;
     }
     case BackgroundInterpolatePixel:
     {
-      *pixel = image->background_color;  /* Copy PixelInfo Structure  */
+      *pixel=image->background_color;  /* Copy PixelInfo Structure  */
       break;
     }
     case BilinearInterpolatePixel:
@@ -5194,7 +5494,7 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
       epsilon.y=1.0-delta.y;
       gamma=((epsilon.y*(epsilon.x*alpha[0]+delta.x*alpha[1])+delta.y*
         (epsilon.x*alpha[2]+delta.x*alpha[3])));
-      gamma=MagickEpsilonReciprocal(gamma);
+      gamma=PerceptibleReciprocal(gamma);
       pixel->red=gamma*(epsilon.y*(epsilon.x*pixels[0].red+delta.x*
         pixels[1].red)+delta.y*(epsilon.x*pixels[2].red+delta.x*pixels[3].red));
       pixel->green=gamma*(epsilon.y*(epsilon.x*pixels[0].green+delta.x*
@@ -5208,7 +5508,7 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
           pixels[1].black)+delta.y*(epsilon.x*pixels[2].black+delta.x*
           pixels[3].black));
       gamma=((epsilon.y*(epsilon.x+delta.x)+delta.y*(epsilon.x+delta.x)));
-      gamma=MagickEpsilonReciprocal(gamma);
+      gamma=PerceptibleReciprocal(gamma);
       pixel->alpha=(epsilon.y*(epsilon.x*pixels[0].alpha+delta.x*
         pixels[1].alpha)+delta.y*(epsilon.x*pixels[2].alpha+delta.x*
         pixels[3].alpha));
@@ -5224,47 +5524,54 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
         }
       for (i=0; i < 4L; i++)
         AlphaBlendPixelInfo(image,p+i*GetPixelChannels(image),pixels+i,alpha+i);
-      gamma=1.0;       /* number of pixels blended together */
-      for (i=0; i <= 1L; i++) {
-        if ( y-y_offset >= 0.75 ) {
-          alpha[i]  = alpha[i+2];
-          pixels[i] = pixels[i+2];
-        }
-        else if ( y-y_offset > 0.25 ) {
-          gamma = 2.0;             /* each y pixels have been blended */
-          alpha[i]        += alpha[i+2];  /* add up alpha weights */
-          pixels[i].red   += pixels[i+2].red;
-          pixels[i].green += pixels[i+2].green;
-          pixels[i].blue  += pixels[i+2].blue;
-          pixels[i].black += pixels[i+2].black;
-          pixels[i].alpha += pixels[i+2].alpha;
-        }
-      }
-      if ( x-x_offset >= 0.75 ) {
-        alpha[0]  = alpha[1];
-        pixels[0] = pixels[1];
-      }
-      else if ( x-x_offset > 0.25 ) {
-        gamma *= 2.0;          /* double number of pixels blended */
-        alpha[0]        += alpha[1];  /* add up alpha weights */
-        pixels[0].red   += pixels[1].red;
-        pixels[0].green += pixels[1].green;
-        pixels[0].blue  += pixels[1].blue;
-        pixels[0].black += pixels[1].black;
-        pixels[0].alpha += pixels[1].alpha;
-      }
-      gamma = 1.0/gamma;
-      alpha[0]=MagickEpsilonReciprocal(alpha[0]);
-      pixel->red   = alpha[0]*pixels[0].red;
-      pixel->green = alpha[0]*pixels[0].green; /* divide by sum of alpha */
-      pixel->blue  = alpha[0]*pixels[0].blue;
-      pixel->black = alpha[0]*pixels[0].black;
-      pixel->alpha =    gamma*pixels[0].alpha; /* divide by number of pixels */
+      gamma=1.0;  /* number of pixels blended together (its variable) */
+      for (i=0; i <= 1L; i++)
+      {
+        if ((y-y_offset) >= 0.75)
+          {
+            alpha[i]=alpha[i+2];  /* take right pixels */
+            pixels[i]=pixels[i+2];
+          }
+        else
+          if ((y-y_offset) > 0.25)
+            {
+              gamma=2.0;  /* blend both pixels in row */
+              alpha[i]+=alpha[i+2];  /* add up alpha weights */
+              pixels[i].red+=pixels[i+2].red;
+              pixels[i].green+=pixels[i+2].green;
+              pixels[i].blue+=pixels[i+2].blue;
+              pixels[i].black+=pixels[i+2].black;
+              pixels[i].alpha+=pixels[i+2].alpha;
+            }
+      }
+      if ((x-x_offset) >= 0.75)
+        {
+          alpha[0]=alpha[1];
+          pixels[0]=pixels[1];
+        }
+      else
+        if ((x-x_offset) > 0.25)
+          {
+            gamma*=2.0;  /* blend both rows */
+            alpha[0]+= alpha[1];  /* add up alpha weights */
+            pixels[0].red+=pixels[1].red;
+            pixels[0].green+=pixels[1].green;
+            pixels[0].blue+=pixels[1].blue;
+            pixels[0].black+=pixels[1].black;
+            pixels[0].alpha+=pixels[1].alpha;
+          }
+      gamma=1.0/gamma;
+      alpha[0]=PerceptibleReciprocal(alpha[0]);
+      pixel->red=alpha[0]*pixels[0].red;
+      pixel->green=alpha[0]*pixels[0].green;  /* divide by sum of alpha */
+      pixel->blue=alpha[0]*pixels[0].blue;
+      pixel->black=alpha[0]*pixels[0].black;
+      pixel->alpha=gamma*pixels[0].alpha;   /* divide by number of pixels */
       break;
     }
     case CatromInterpolatePixel:
     {
-      MagickRealType
+      double
         cx[4],
         cy[4];
 
@@ -5277,88 +5584,43 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
         }
       for (i=0; i < 16L; i++)
         AlphaBlendPixelInfo(image,p+i*GetPixelChannels(image),pixels+i,alpha+i);
-      CatromWeights((MagickRealType) (x-x_offset),&cx);
-      CatromWeights((MagickRealType) (y-y_offset),&cy);
-      pixel->red=(cy[0]*(cx[0]*pixels[0].red+cx[1]*
-        pixels[1].red+cx[2]*pixels[2].red+cx[3]*
-        pixels[3].red)+cy[1]*(cx[0]*pixels[4].red+cx[1]*
-        pixels[5].red+cx[2]*pixels[6].red+cx[3]*
-        pixels[7].red)+cy[2]*(cx[0]*pixels[8].red+cx[1]*
-        pixels[9].red+cx[2]*pixels[10].red+cx[3]*
-        pixels[11].red)+cy[3]*(cx[0]*pixels[12].red+cx[1]*
-        pixels[13].red+cx[2]*pixels[14].red+cx[3]*pixels[15].red));
-      pixel->green=(cy[0]*(cx[0]*pixels[0].green+cx[1]*
-        pixels[1].green+cx[2]*pixels[2].green+cx[3]*
-        pixels[3].green)+cy[1]*(cx[0]*pixels[4].green+cx[1]*
-        pixels[5].green+cx[2]*pixels[6].green+cx[3]*
-        pixels[7].green)+cy[2]*(cx[0]*pixels[8].green+cx[1]*
-        pixels[9].green+cx[2]*pixels[10].green+cx[3]*
-        pixels[11].green)+cy[3]*(cx[0]*pixels[12].green+cx[1]*
-        pixels[13].green+cx[2]*pixels[14].green+cx[3]*pixels[15].green));
-      pixel->blue=(cy[0]*(cx[0]*pixels[0].blue+cx[1]*
-        pixels[1].blue+cx[2]*pixels[2].blue+cx[3]*
-        pixels[3].blue)+cy[1]*(cx[0]*pixels[4].blue+cx[1]*
-        pixels[5].blue+cx[2]*pixels[6].blue+cx[3]*
-        pixels[7].blue)+cy[2]*(cx[0]*pixels[8].blue+cx[1]*
-        pixels[9].blue+cx[2]*pixels[10].blue+cx[3]*
-        pixels[11].blue)+cy[3]*(cx[0]*pixels[12].blue+cx[1]*
-        pixels[13].blue+cx[2]*pixels[14].blue+cx[3]*pixels[15].blue));
+      CatromWeights((double) (x-x_offset),&cx);
+      CatromWeights((double) (y-y_offset),&cy);
+      pixel->red=(cy[0]*(cx[0]*pixels[0].red+cx[1]*pixels[1].red+cx[2]*
+        pixels[2].red+cx[3]*pixels[3].red)+cy[1]*(cx[0]*pixels[4].red+cx[1]*
+        pixels[5].red+cx[2]*pixels[6].red+cx[3]*pixels[7].red)+cy[2]*(cx[0]*
+        pixels[8].red+cx[1]*pixels[9].red+cx[2]*pixels[10].red+cx[3]*
+        pixels[11].red)+cy[3]*(cx[0]*pixels[12].red+cx[1]*pixels[13].red+cx[2]*
+        pixels[14].red+cx[3]*pixels[15].red));
+      pixel->green=(cy[0]*(cx[0]*pixels[0].green+cx[1]*pixels[1].green+cx[2]*
+        pixels[2].green+cx[3]*pixels[3].green)+cy[1]*(cx[0]*pixels[4].green+
+        cx[1]*pixels[5].green+cx[2]*pixels[6].green+cx[3]*pixels[7].green)+
+        cy[2]*(cx[0]*pixels[8].green+cx[1]*pixels[9].green+cx[2]*
+        pixels[10].green+cx[3]*pixels[11].green)+cy[3]*(cx[0]*
+        pixels[12].green+cx[1]*pixels[13].green+cx[2]*pixels[14].green+cx[3]*
+        pixels[15].green));
+      pixel->blue=(cy[0]*(cx[0]*pixels[0].blue+cx[1]*pixels[1].blue+cx[2]*
+        pixels[2].blue+cx[3]*pixels[3].blue)+cy[1]*(cx[0]*pixels[4].blue+cx[1]*
+        pixels[5].blue+cx[2]*pixels[6].blue+cx[3]*pixels[7].blue)+cy[2]*(cx[0]*
+        pixels[8].blue+cx[1]*pixels[9].blue+cx[2]*pixels[10].blue+cx[3]*
+        pixels[11].blue)+cy[3]*(cx[0]*pixels[12].blue+cx[1]*pixels[13].blue+
+        cx[2]*pixels[14].blue+cx[3]*pixels[15].blue));
       if (image->colorspace == CMYKColorspace)
-        pixel->black=(cy[0]*(cx[0]*pixels[0].black+cx[1]*
-          pixels[1].black+cx[2]*pixels[2].black+cx[3]*
-          pixels[3].black)+cy[1]*(cx[0]*pixels[4].black+cx[1]*
-          pixels[5].black+cx[2]*pixels[6].black+cx[3]*
-          pixels[7].black)+cy[2]*(cx[0]*pixels[8].black+cx[1]*
-          pixels[9].black+cx[2]*pixels[10].black+cx[3]*
-          pixels[11].black)+cy[3]*(cx[0]*pixels[12].black+cx[1]*
-          pixels[13].black+cx[2]*pixels[14].black+cx[3]*pixels[15].black));
-      pixel->alpha=(cy[0]*(cx[0]*pixels[0].alpha+cx[1]*
-        pixels[1].alpha+cx[2]*pixels[2].alpha+cx[3]*
-        pixels[3].alpha)+cy[1]*(cx[0]*pixels[4].alpha+cx[1]*
-        pixels[5].alpha+cx[2]*pixels[6].alpha+cx[3]*
-        pixels[7].alpha)+cy[2]*(cx[0]*pixels[8].alpha+cx[1]*
-        pixels[9].alpha+cx[2]*pixels[10].alpha+cx[3]*
-        pixels[11].alpha)+cy[3]*(cx[0]*pixels[12].alpha+cx[1]*
-        pixels[13].alpha+cx[2]*pixels[14].alpha+cx[3]*pixels[15].alpha));
-      break;
-    }
-#if 0
-    /* deprecated useless and very slow interpolator */
-    case FilterInterpolatePixel:
-    {
-      CacheView
-        *filter_view;
-
-      Image
-        *excerpt_image,
-        *filter_image;
-
-      RectangleInfo
-        geometry;
-
-      geometry.width=4L;
-      geometry.height=4L;
-      geometry.x=x_offset-1;
-      geometry.y=y_offset-1;
-      excerpt_image=ExcerptImage(image,&geometry,exception);
-      if (excerpt_image == (Image *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      filter_image=ResizeImage(excerpt_image,1,1,image->filter,exception);
-      excerpt_image=DestroyImage(excerpt_image);
-      if (filter_image == (Image *) NULL)
-        break;
-      filter_view=AcquireVirtualCacheView(filter_image,exception);
-      p=GetCacheViewVirtualPixels(filter_view,0,0,1,1,exception);
-      if (p != (const Quantum *) NULL)
-        GetPixelInfoPixel(image,p,pixel);
-      filter_view=DestroyCacheView(filter_view);
-      filter_image=DestroyImage(filter_image);
+        pixel->black=(cy[0]*(cx[0]*pixels[0].black+cx[1]*pixels[1].black+cx[2]*
+          pixels[2].black+cx[3]*pixels[3].black)+cy[1]*(cx[0]*pixels[4].black+
+          cx[1]*pixels[5].black+cx[2]*pixels[6].black+cx[3]*pixels[7].black)+
+          cy[2]*(cx[0]*pixels[8].black+cx[1]*pixels[9].black+cx[2]*
+          pixels[10].black+cx[3]*pixels[11].black)+cy[3]*(cx[0]*
+          pixels[12].black+cx[1]*pixels[13].black+cx[2]*pixels[14].black+cx[3]*
+          pixels[15].black));
+      pixel->alpha=(cy[0]*(cx[0]*pixels[0].alpha+cx[1]*pixels[1].alpha+cx[2]*
+        pixels[2].alpha+cx[3]*pixels[3].alpha)+cy[1]*(cx[0]*pixels[4].alpha+
+        cx[1]*pixels[5].alpha+cx[2]*pixels[6].alpha+cx[3]*pixels[7].alpha)+
+        cy[2]*(cx[0]*pixels[8].alpha+cx[1]*pixels[9].alpha+cx[2]*
+        pixels[10].alpha+cx[3]*pixels[11].alpha)+cy[3]*(cx[0]*pixels[12].alpha+
+        cx[1]*pixels[13].alpha+cx[2]*pixels[14].alpha+cx[3]*pixels[15].alpha));
       break;
     }
-#endif
     case IntegerInterpolatePixel:
     {
       p=GetCacheViewVirtualPixels(image_view,x_offset,y_offset,1,1,exception);
@@ -5384,10 +5646,10 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
         }
       delta.x=x-x_offset;
       delta.y=y-y_offset;
-      luminance.x=GetPixelLuminance(image,p)-(double)
-        GetPixelLuminance(image,p+3*GetPixelChannels(image));
-      luminance.y=GetPixelLuminance(image,p+GetPixelChannels(image))-(double)
-        GetPixelLuminance(image,p+2*GetPixelChannels(image));
+      luminance.x=GetPixelLuma(image,p)-(double)
+        GetPixelLuma(image,p+3*GetPixelChannels(image));
+      luminance.y=GetPixelLuma(image,p+GetPixelChannels(image))-(double)
+        GetPixelLuma(image,p+2*GetPixelChannels(image));
       AlphaBlendPixelInfo(image,p,pixels+0,alpha+0);
       AlphaBlendPixelInfo(image,p+GetPixelChannels(image),pixels+1,alpha+1);
       AlphaBlendPixelInfo(image,p+2*GetPixelChannels(image),pixels+2,alpha+2);
@@ -5404,7 +5666,7 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
               */
               delta.y=1.0-delta.y;
               gamma=MeshInterpolate(&delta,alpha[2],alpha[3],alpha[0]);
-              gamma=MagickEpsilonReciprocal(gamma);
+              gamma=PerceptibleReciprocal(gamma);
               pixel->red=gamma*MeshInterpolate(&delta,pixels[2].red,
                 pixels[3].red,pixels[0].red);
               pixel->green=gamma*MeshInterpolate(&delta,pixels[2].green,
@@ -5425,7 +5687,7 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
               */
               delta.x=1.0-delta.x;
               gamma=MeshInterpolate(&delta,alpha[1],alpha[0],alpha[3]);
-              gamma=MagickEpsilonReciprocal(gamma);
+              gamma=PerceptibleReciprocal(gamma);
               pixel->red=gamma*MeshInterpolate(&delta,pixels[1].red,
                 pixels[0].red,pixels[3].red);
               pixel->green=gamma*MeshInterpolate(&delta,pixels[1].green,
@@ -5451,7 +5713,7 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
                 Top-left triangle (pixel: 0, diagonal: 1-2).
               */
               gamma=MeshInterpolate(&delta,alpha[0],alpha[1],alpha[2]);
-              gamma=MagickEpsilonReciprocal(gamma);
+              gamma=PerceptibleReciprocal(gamma);
               pixel->red=gamma*MeshInterpolate(&delta,pixels[0].red,
                 pixels[1].red,pixels[2].red);
               pixel->green=gamma*MeshInterpolate(&delta,pixels[0].green,
@@ -5473,7 +5735,7 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
               delta.x=1.0-delta.x;
               delta.y=1.0-delta.y;
               gamma=MeshInterpolate(&delta,alpha[3],alpha[2],alpha[1]);
-              gamma=MagickEpsilonReciprocal(gamma);
+              gamma=PerceptibleReciprocal(gamma);
               pixel->red=gamma*MeshInterpolate(&delta,pixels[3].red,
                 pixels[2].red,pixels[1].red);
               pixel->green=gamma*MeshInterpolate(&delta,pixels[3].green,
@@ -5505,7 +5767,7 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
     }
     case SplineInterpolatePixel:
     {
-      MagickRealType
+      double
         cx[4],
         cy[4];
 
@@ -5518,49 +5780,40 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
         }
       for (i=0; i < 16L; i++)
         AlphaBlendPixelInfo(image,p+i*GetPixelChannels(image),pixels+i,alpha+i);
-      SplineWeights((MagickRealType) (x-x_offset),&cx);
-      SplineWeights((MagickRealType) (y-y_offset),&cy);
-      pixel->red=(cy[0]*(cx[0]*pixels[0].red+cx[1]*
-        pixels[1].red+cx[2]*pixels[2].red+cx[3]*
-        pixels[3].red)+cy[1]*(cx[0]*pixels[4].red+cx[1]*
-        pixels[5].red+cx[2]*pixels[6].red+cx[3]*
-        pixels[7].red)+cy[2]*(cx[0]*pixels[8].red+cx[1]*
-        pixels[9].red+cx[2]*pixels[10].red+cx[3]*
-        pixels[11].red)+cy[3]*(cx[0]*pixels[12].red+cx[1]*
-        pixels[13].red+cx[2]*pixels[14].red+cx[3]*pixels[15].red));
-      pixel->green=(cy[0]*(cx[0]*pixels[0].green+cx[1]*
-        pixels[1].green+cx[2]*pixels[2].green+cx[3]*
-        pixels[3].green)+cy[1]*(cx[0]*pixels[4].green+cx[1]*
-        pixels[5].green+cx[2]*pixels[6].green+cx[3]*
-        pixels[7].green)+cy[2]*(cx[0]*pixels[8].green+cx[1]*
-        pixels[9].green+cx[2]*pixels[10].green+cx[3]*
-        pixels[11].green)+cy[3]*(cx[0]*pixels[12].green+cx[1]*
-        pixels[13].green+cx[2]*pixels[14].green+cx[3]*pixels[15].green));
-      pixel->blue=(cy[0]*(cx[0]*pixels[0].blue+cx[1]*
-        pixels[1].blue+cx[2]*pixels[2].blue+cx[3]*
-        pixels[3].blue)+cy[1]*(cx[0]*pixels[4].blue+cx[1]*
-        pixels[5].blue+cx[2]*pixels[6].blue+cx[3]*
-        pixels[7].blue)+cy[2]*(cx[0]*pixels[8].blue+cx[1]*
-        pixels[9].blue+cx[2]*pixels[10].blue+cx[3]*
-        pixels[11].blue)+cy[3]*(cx[0]*pixels[12].blue+cx[1]*
-        pixels[13].blue+cx[2]*pixels[14].blue+cx[3]*pixels[15].blue));
+      SplineWeights((double) (x-x_offset),&cx);
+      SplineWeights((double) (y-y_offset),&cy);
+      pixel->red=(cy[0]*(cx[0]*pixels[0].red+cx[1]*pixels[1].red+cx[2]*
+        pixels[2].red+cx[3]*pixels[3].red)+cy[1]*(cx[0]*pixels[4].red+cx[1]*
+        pixels[5].red+cx[2]*pixels[6].red+cx[3]*pixels[7].red)+cy[2]*(cx[0]*
+        pixels[8].red+cx[1]*pixels[9].red+cx[2]*pixels[10].red+cx[3]*
+        pixels[11].red)+cy[3]*(cx[0]*pixels[12].red+cx[1]*pixels[13].red+cx[2]*
+        pixels[14].red+cx[3]*pixels[15].red));
+      pixel->green=(cy[0]*(cx[0]*pixels[0].green+cx[1]*pixels[1].green+cx[2]*
+        pixels[2].green+cx[3]*pixels[3].green)+cy[1]*(cx[0]*pixels[4].green+
+        cx[1]*pixels[5].green+cx[2]*pixels[6].green+cx[3]*pixels[7].green)+
+        cy[2]*(cx[0]*pixels[8].green+cx[1]*pixels[9].green+cx[2]*
+        pixels[10].green+cx[3]*pixels[11].green)+cy[3]*(cx[0]*pixels[12].green+
+        cx[1]*pixels[13].green+cx[2]*pixels[14].green+cx[3]*pixels[15].green));
+      pixel->blue=(cy[0]*(cx[0]*pixels[0].blue+cx[1]*pixels[1].blue+cx[2]*
+        pixels[2].blue+cx[3]*pixels[3].blue)+cy[1]*(cx[0]*pixels[4].blue+cx[1]*
+        pixels[5].blue+cx[2]*pixels[6].blue+cx[3]*pixels[7].blue)+cy[2]*(cx[0]*
+        pixels[8].blue+cx[1]*pixels[9].blue+cx[2]*pixels[10].blue+cx[3]*
+        pixels[11].blue)+cy[3]*(cx[0]*pixels[12].blue+cx[1]*pixels[13].blue+
+        cx[2]*pixels[14].blue+cx[3]*pixels[15].blue));
       if (image->colorspace == CMYKColorspace)
-        pixel->black=(cy[0]*(cx[0]*pixels[0].black+cx[1]*
-          pixels[1].black+cx[2]*pixels[2].black+cx[3]*
-          pixels[3].black)+cy[1]*(cx[0]*pixels[4].black+cx[1]*
-          pixels[5].black+cx[2]*pixels[6].black+cx[3]*
-          pixels[7].black)+cy[2]*(cx[0]*pixels[8].black+cx[1]*
-          pixels[9].black+cx[2]*pixels[10].black+cx[3]*
-          pixels[11].black)+cy[3]*(cx[0]*pixels[12].black+cx[1]*
-          pixels[13].black+cx[2]*pixels[14].black+cx[3]*pixels[15].black));
-      pixel->alpha=(cy[0]*(cx[0]*pixels[0].alpha+cx[1]*
-        pixels[1].alpha+cx[2]*pixels[2].alpha+cx[3]*
-        pixels[3].alpha)+cy[1]*(cx[0]*pixels[4].alpha+cx[1]*
-        pixels[5].alpha+cx[2]*pixels[6].alpha+cx[3]*
-        pixels[7].alpha)+cy[2]*(cx[0]*pixels[8].alpha+cx[1]*
-        pixels[9].alpha+cx[2]*pixels[10].alpha+cx[3]*
-        pixels[11].alpha)+cy[3]*(cx[0]*pixels[12].alpha+cx[1]*
-        pixels[13].alpha+cx[2]*pixels[14].alpha+cx[3]*pixels[15].alpha));
+        pixel->black=(cy[0]*(cx[0]*pixels[0].black+cx[1]*pixels[1].black+cx[2]*
+          pixels[2].black+cx[3]*pixels[3].black)+cy[1]*(cx[0]*pixels[4].black+
+          cx[1]*pixels[5].black+cx[2]*pixels[6].black+cx[3]*pixels[7].black)+
+          cy[2]*(cx[0]*pixels[8].black+cx[1]*pixels[9].black+cx[2]*
+          pixels[10].black+cx[3]*pixels[11].black)+cy[3]*(cx[0]*
+          pixels[12].black+cx[1]*pixels[13].black+cx[2]*pixels[14].black+cx[3]*
+          pixels[15].black));
+      pixel->alpha=(cy[0]*(cx[0]*pixels[0].alpha+cx[1]*pixels[1].alpha+cx[2]*
+        pixels[2].alpha+cx[3]*pixels[3].alpha)+cy[1]*(cx[0]*pixels[4].alpha+
+        cx[1]*pixels[5].alpha+cx[2]*pixels[6].alpha+cx[3]*pixels[7].alpha)+
+        cy[2]*(cx[0]*pixels[8].alpha+cx[1]*pixels[9].alpha+cx[2]*
+        pixels[10].alpha+cx[3]*pixels[11].alpha)+cy[3]*(cx[0]*pixels[12].alpha+
+        cx[1]*pixels[13].alpha+cx[2]*pixels[14].alpha+cx[3]*pixels[15].alpha));
       break;
     }
   }
@@ -5601,25 +5854,24 @@ MagickExport MagickBooleanType InterpolatePixelInfo(const Image *image,
 MagickExport MagickBooleanType IsFuzzyEquivalencePixel(const Image *source,
   const Quantum *p,const Image *destination,const Quantum *q)
 {
-  MagickRealType
+  double
     fuzz,
     pixel;
 
-  register MagickRealType
+  register double
     distance,
     scale;
 
-  fuzz=MagickMax(source->fuzz,(MagickRealType) MagickSQ1_2)*MagickMax(
+  fuzz=(double) MagickMax(source->fuzz,(MagickRealType) MagickSQ1_2)*MagickMax(
     destination->fuzz,(MagickRealType) MagickSQ1_2);
   scale=1.0;
   distance=0.0;
-  if (source->matte != MagickFalse)
+  if (source->alpha_trait == BlendPixelTrait)
     {
       /*
         Transparencies are involved - set alpha distance
       */
-      pixel=GetPixelAlpha(source,p)-(MagickRealType)
-        GetPixelAlpha(destination,q);
+      pixel=GetPixelAlpha(source,p)-(double) GetPixelAlpha(destination,q);
       distance=pixel*pixel;
       if (distance > fuzz)
         return(MagickFalse);
@@ -5637,7 +5889,7 @@ MagickExport MagickBooleanType IsFuzzyEquivalencePixel(const Image *source,
   */
   distance*=3.0;  /* rescale appropriately */
   fuzz*=3.0;
-  pixel=GetPixelRed(source,p)-(MagickRealType) GetPixelRed(destination,q);
+  pixel=GetPixelRed(source,p)-(double) GetPixelRed(destination,q);
   if ((source->colorspace == HSLColorspace) ||
       (source->colorspace == HSBColorspace) ||
       (source->colorspace == HWBColorspace))
@@ -5653,11 +5905,11 @@ MagickExport MagickBooleanType IsFuzzyEquivalencePixel(const Image *source,
   distance+=scale*pixel*pixel;
   if (distance > fuzz)
     return(MagickFalse);
-  pixel=GetPixelGreen(source,p)-(MagickRealType) GetPixelGreen(destination,q);
+  pixel=GetPixelGreen(source,p)-(double) GetPixelGreen(destination,q);
   distance+=scale*pixel*pixel;
   if (distance > fuzz)
     return(MagickFalse);
-  pixel=GetPixelBlue(source,p)-(MagickRealType) GetPixelBlue(destination,q);
+  pixel=GetPixelBlue(source,p)-(double) GetPixelBlue(destination,q);
   distance+=scale*pixel*pixel;
   if (distance > fuzz)
     return(MagickFalse);
@@ -5711,34 +5963,35 @@ MagickExport MagickBooleanType IsFuzzyEquivalencePixel(const Image *source,
 MagickExport MagickBooleanType IsFuzzyEquivalencePixelInfo(const PixelInfo *p,
   const PixelInfo *q)
 {
-  MagickRealType
+  double
     fuzz,
     pixel;
 
-  register MagickRealType
+  register double
     scale,
     distance;
 
   if ((p->fuzz == 0.0) && (q->fuzz == 0.0))
     return(IsPixelInfoEquivalent(p,q));
   if (p->fuzz == 0.0)
-    fuzz=MagickMax(q->fuzz,(MagickRealType) MagickSQ1_2)*MagickMax(q->fuzz,
-      (MagickRealType) MagickSQ1_2);
+    fuzz=MagickMax(q->fuzz,(MagickRealType) MagickSQ1_2)*
+      MagickMax(q->fuzz,(MagickRealType) MagickSQ1_2);
   else if (q->fuzz == 0.0)
-    fuzz=MagickMax(p->fuzz,(MagickRealType) MagickSQ1_2)*MagickMax(p->fuzz,
-      (MagickRealType) MagickSQ1_2);
+    fuzz=MagickMax(p->fuzz,(MagickRealType) MagickSQ1_2)*
+      MagickMax(p->fuzz,(MagickRealType) MagickSQ1_2);
   else
-    fuzz=MagickMax(p->fuzz,(MagickRealType) MagickSQ1_2)*MagickMax(q->fuzz,
-      (MagickRealType) MagickSQ1_2);
+    fuzz=MagickMax(p->fuzz,(MagickRealType) MagickSQ1_2)*
+      MagickMax(q->fuzz,(MagickRealType) MagickSQ1_2);
   scale=1.0;
   distance=0.0;
-  if ((p->matte != MagickFalse) || (q->matte != MagickFalse))
+  if ((p->alpha_trait == BlendPixelTrait) ||
+      (q->alpha_trait == BlendPixelTrait))
     {
       /*
         Transparencies are involved - set alpha distance.
       */
-      pixel=(p->matte != MagickFalse ? p->alpha : OpaqueAlpha)-
-        (q->matte != MagickFalse ? q->alpha : OpaqueAlpha);
+      pixel=(p->alpha_trait == BlendPixelTrait ? p->alpha : OpaqueAlpha)-
+        (q->alpha_trait == BlendPixelTrait ? q->alpha : OpaqueAlpha);
       distance=pixel*pixel;
       if (distance > fuzz)
         return(MagickFalse);
@@ -5746,9 +5999,9 @@ MagickExport MagickBooleanType IsFuzzyEquivalencePixelInfo(const PixelInfo *p,
         Generate a alpha scaling factor to generate a 4D cone on colorspace.
         If one color is transparent, distance has no color component.
       */
-      if (p->matte != MagickFalse)
+      if (p->alpha_trait == BlendPixelTrait)
         scale=(QuantumScale*p->alpha);
-      if (q->matte != MagickFalse)
+      if (q->alpha_trait == BlendPixelTrait)
         scale*=(QuantumScale*q->alpha);
       if (scale <= MagickEpsilon )
         return(MagickTrue);
@@ -5762,8 +6015,8 @@ MagickExport MagickBooleanType IsFuzzyEquivalencePixelInfo(const PixelInfo *p,
       distance+=pixel*pixel*scale;
       if (distance > fuzz)
         return(MagickFalse);
-      scale*=(MagickRealType) (QuantumScale*(QuantumRange-p->black));
-      scale*=(MagickRealType) (QuantumScale*(QuantumRange-q->black));
+      scale*=(double) (QuantumScale*(QuantumRange-p->black));
+      scale*=(double) (QuantumScale*(QuantumRange-q->black));
     }
   /*
     RGB or CMY color cube.
@@ -5775,9 +6028,9 @@ MagickExport MagickBooleanType IsFuzzyEquivalencePixelInfo(const PixelInfo *p,
       (p->colorspace == HWBColorspace))
     {
       /*
-        This calculates a arc distance for hue-- it should be a vector angle
-        of 'S'/'W' length with 'L'/'B' forming appropriate cones.  In other
-        words this is a hack - Anthony.
+        This calculates a arc distance for hue-- it should be a vector
+        angle of 'S'/'W' length with 'L'/'B' forming appropriate cones.
+        In other words this is a hack - Anthony.
       */
       if (fabs((double) pixel) > (QuantumRange/2))
         pixel-=QuantumRange;
@@ -5802,18 +6055,18 @@ MagickExport MagickBooleanType IsFuzzyEquivalencePixelInfo(const PixelInfo *p,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   S e t P i x e l C h a n n e l M a p M a s k                               %
+%   S e t P i x e l C h a n n e l M a s k                                     %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  SetPixelChannelMapMask() sets the pixel channel map from the specified
-%  channel mask.
+%  SetPixelChannelMask() sets the pixel channel map from the specified channel
+%  mask.
 %
-%  The format of the SetPixelChannelMapMask method is:
+%  The format of the SetPixelChannelMask method is:
 %
-%      void SetPixelChannelMapMask(Image *image,const ChannelType channel_mask)
+%      void SetPixelChannelMask(Image *image,const ChannelType channel_mask)
 %
 %  A description of each parameter follows:
 %
@@ -5822,7 +6075,7 @@ MagickExport MagickBooleanType IsFuzzyEquivalencePixelInfo(const PixelInfo *p,
 %    o channel_mask: the channel mask.
 %
 */
-MagickExport void SetPixelChannelMapMask(Image *image,
+MagickExport void SetPixelChannelMask(Image *image,
   const ChannelType channel_mask)
 {
 #define GetChannelBit(mask,bit)  (((size_t) (mask) >> (size_t) (bit)) & 0x01)
@@ -5836,19 +6089,18 @@ MagickExport void SetPixelChannelMapMask(Image *image,
   image->channel_mask=channel_mask;
   for (i=0; i < (ssize_t) GetPixelChannels(image); i++)
   {
-    PixelChannel
-      channel;
-
-    channel=GetPixelChannelMapChannel(image,i);
-    SetPixelChannelMapTraits(image,channel,
+    PixelChannel channel=GetPixelChannelChannel(image,i);
+    SetPixelChannelTraits(image,channel,
       GetChannelBit(channel_mask,channel) == 0 ? CopyPixelTrait :
-      image->matte == MagickFalse || (channel == AlphaPixelChannel) ?
-      UpdatePixelTrait : (PixelTrait) (UpdatePixelTrait | BlendPixelTrait));
+      image->alpha_trait != BlendPixelTrait || (channel == AlphaPixelChannel) ?
+      UpdatePixelTrait : (PixelTrait) (UpdatePixelTrait | image->alpha_trait));
   }
   if (image->storage_class == PseudoClass)
-    SetPixelChannelMapTraits(image,IndexPixelChannel,CopyPixelTrait);
-  if (image->mask != MagickFalse)
-    SetPixelChannelMapTraits(image,MaskPixelChannel,CopyPixelTrait);
+    SetPixelChannelTraits(image,IndexPixelChannel,CopyPixelTrait);
+  if (image->read_mask != MagickFalse)
+    SetPixelChannelTraits(image,ReadMaskPixelChannel,CopyPixelTrait);
+  if (image->write_mask != MagickFalse)
+    SetPixelChannelTraits(image,WriteMaskPixelChannel,CopyPixelTrait);
   if (image->debug != MagickFalse)
     LogPixelChannels(image);
 }
@@ -5858,44 +6110,6 @@ MagickExport void SetPixelChannelMapMask(Image *image,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   S e t P i x e l C h a n n e l M a s k                                     %
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%  SetPixelChannelMask() sets the pixel channel mask from the specified channel
-%  mask.
-%
-%  The format of the SetPixelChannelMask method is:
-%
-%      ChannelType SetPixelChannelMask(Image *image,
-%        const ChannelType channel_mask)
-%
-%  A description of each parameter follows:
-%
-%    o image: the image.
-%
-%    o channel_mask: the channel mask.
-%
-*/
-MagickExport ChannelType SetPixelChannelMask(Image *image,
-  const ChannelType channel_mask)
-{
-  ChannelType
-    mask;
-
-  mask=image->channel_mask;
-  image->channel_mask=channel_mask;
-  SetPixelChannelMapMask(image,channel_mask);
-  return(mask);
-}
-\f
-/*
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                             %
-%                                                                             %
-%                                                                             %
 %   S e t P i x e l M e t a C h a n n e l s                                   %
 %                                                                             %
 %                                                                             %