]> granicus.if.org Git - imagemagick/blobdiff - magick/resample.c
Horizon validity (anti-aliased) added to Plane2Cylinder
[imagemagick] / magick / resample.c
index 537715da47b085009f9ce4c0c348f232e909b0fb..755b154de48213494a0a771c14aaf574501d7b1d 100644 (file)
@@ -52,6 +52,7 @@
 #include "magick/log.h"
 #include "magick/magick.h"
 #include "magick/memory_.h"
+#include "magick/pixel.h"
 #include "magick/pixel-private.h"
 #include "magick/quantum.h"
 #include "magick/random_.h"
@@ -92,12 +93,6 @@ struct _ResampleFilter
   Image
     *image;
 
-  MagickBooleanType
-    matte;
-
-  ColorspaceType
-    colorspace;
-
   ExceptionInfo
     *exception;
 
@@ -226,7 +221,6 @@ MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
 
   resample_filter->exception=exception;
   resample_filter->image=ReferenceImage((Image *) image);
-  resample_filter->matte=image->matte;
   resample_filter->view=AcquireCacheView(resample_filter->image);
 
   resample_filter->debug=IsEventLogging();
@@ -293,575 +287,6 @@ MagickExport ResampleFilter *DestroyResampleFilter(
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   I n t e r p o l a t e R e s a m p l e F i l t e r                         %
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%  InterpolateResampleFilter() applies bi-linear or tri-linear interpolation
-%  between a floating point coordinate and the pixels surrounding that
-%  coordinate.  No pixel area resampling, or scaling of the result is
-%  performed.
-%
-%  The format of the InterpolateResampleFilter method is:
-%
-%      MagickBooleanType InterpolateResampleFilter(
-%        ResampleInfo *resample_filter,const InterpolatePixelMethod method,
-%        const double x,const double y,MagickPixelPacket *pixel)
-%
-%  A description of each parameter follows:
-%
-%    o resample_filter: the resample filter.
-%
-%    o method: the pixel clor interpolation method.
-%
-%    o x,y: A double representing the current (x,y) position of the pixel.
-%
-%    o pixel: return the interpolated pixel here.
-%
-*/
-
-static inline double MagickMax(const double x,const double y)
-{
-  if (x > y)
-    return(x);
-  return(y);
-}
-
-static void BicubicInterpolate(const MagickPixelPacket *pixels,const double dx,
-  MagickPixelPacket *pixel)
-{
-  MagickRealType
-    dx2,
-    p,
-    q,
-    r,
-    s;
-
-  dx2=dx*dx;
-  p=(pixels[3].red-pixels[2].red)-(pixels[0].red-pixels[1].red);
-  q=(pixels[0].red-pixels[1].red)-p;
-  r=pixels[2].red-pixels[0].red;
-  s=pixels[1].red;
-  pixel->red=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
-  p=(pixels[3].green-pixels[2].green)-(pixels[0].green-pixels[1].green);
-  q=(pixels[0].green-pixels[1].green)-p;
-  r=pixels[2].green-pixels[0].green;
-  s=pixels[1].green;
-  pixel->green=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
-  p=(pixels[3].blue-pixels[2].blue)-(pixels[0].blue-pixels[1].blue);
-  q=(pixels[0].blue-pixels[1].blue)-p;
-  r=pixels[2].blue-pixels[0].blue;
-  s=pixels[1].blue;
-  pixel->blue=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
-  p=(pixels[3].opacity-pixels[2].opacity)-(pixels[0].opacity-pixels[1].opacity);
-  q=(pixels[0].opacity-pixels[1].opacity)-p;
-  r=pixels[2].opacity-pixels[0].opacity;
-  s=pixels[1].opacity;
-  pixel->opacity=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
-  if (pixel->colorspace == CMYKColorspace)
-    {
-      p=(pixels[3].index-pixels[2].index)-(pixels[0].index-pixels[1].index);
-      q=(pixels[0].index-pixels[1].index)-p;
-      r=pixels[2].index-pixels[0].index;
-      s=pixels[1].index;
-      pixel->index=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
-    }
-}
-
-static inline MagickRealType CubicWeightingFunction(const MagickRealType x)
-{
-  MagickRealType
-    alpha,
-    gamma;
-
-  alpha=MagickMax(x+2.0,0.0);
-  gamma=1.0*alpha*alpha*alpha;
-  alpha=MagickMax(x+1.0,0.0);
-  gamma-=4.0*alpha*alpha*alpha;
-  alpha=MagickMax(x+0.0,0.0);
-  gamma+=6.0*alpha*alpha*alpha;
-  alpha=MagickMax(x-1.0,0.0);
-  gamma-=4.0*alpha*alpha*alpha;
-  return(gamma/6.0);
-}
-
-static inline double MeshInterpolate(const PointInfo *delta,const double p,
-  const double x,const double y)
-{
-  return(delta->x*x+delta->y*y+(1.0-delta->x-delta->y)*p);
-}
-
-static inline ssize_t NearestNeighbor(MagickRealType x)
-{
-  if (x >= 0.0)
-    return((ssize_t) (x+0.5));
-  return((ssize_t) (x-0.5));
-}
-
-static MagickBooleanType InterpolateResampleFilter(
-  ResampleFilter *resample_filter,const InterpolatePixelMethod method,
-  const double x,const double y,MagickPixelPacket *pixel)
-{
-  MagickBooleanType
-    status;
-
-  register const IndexPacket
-    *indexes;
-
-  register const PixelPacket
-    *p;
-
-  register ssize_t
-    i;
-
-  assert(resample_filter != (ResampleFilter *) NULL);
-  assert(resample_filter->signature == MagickSignature);
-  status=MagickTrue;
-
-  switch (method)
-  {
-    case AverageInterpolatePixel:
-    {
-      MagickPixelPacket
-        pixels[16];
-
-      MagickRealType
-        alpha[16],
-        gamma;
-
-      p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
-        (ssize_t) floor(y)-1,4,4,resample_filter->exception);
-      if (p == (const PixelPacket *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
-      for (i=0; i < 16L; i++)
-      {
-        GetMagickPixelPacket(resample_filter->image,pixels+i);
-        SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
-        alpha[i]=1.0;
-        if (resample_filter->matte != MagickFalse)
-          {
-            alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
-            pixels[i].red*=alpha[i];
-            pixels[i].green*=alpha[i];
-            pixels[i].blue*=alpha[i];
-            if (pixels[i].colorspace == CMYKColorspace)
-              pixels[i].index*=alpha[i];
-          }
-        gamma=alpha[i];
-        gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-        pixel->red+=gamma*0.0625*pixels[i].red;
-        pixel->green+=gamma*0.0625*pixels[i].green;
-        pixel->blue+=gamma*0.0625*pixels[i].blue;
-        pixel->opacity+=0.0625*pixels[i].opacity;
-        if (pixel->colorspace == CMYKColorspace)
-          pixel->index+=gamma*0.0625*pixels[i].index;
-        p++;
-      }
-      break;
-    }
-    case BicubicInterpolatePixel:
-    {
-      MagickPixelPacket
-        pixels[16],
-        u[4];
-
-      MagickRealType
-        alpha[16];
-
-      PointInfo
-        delta;
-
-      p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
-        (ssize_t) floor(y)-1,4,4,resample_filter->exception);
-      if (p == (const PixelPacket *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
-      for (i=0; i < 16L; i++)
-      {
-        GetMagickPixelPacket(resample_filter->image,pixels+i);
-        SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
-        alpha[i]=1.0;
-        if (resample_filter->matte != MagickFalse)
-          {
-            alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
-            pixels[i].red*=alpha[i];
-            pixels[i].green*=alpha[i];
-            pixels[i].blue*=alpha[i];
-            if (pixels[i].colorspace == CMYKColorspace)
-              pixels[i].index*=alpha[i];
-          }
-        p++;
-      }
-      delta.x=x-floor(x);
-      for (i=0; i < 4L; i++)
-        BicubicInterpolate(pixels+4*i,delta.x,u+i);
-      delta.y=y-floor(y);
-      BicubicInterpolate(u,delta.y,pixel);
-      break;
-    }
-    case BilinearInterpolatePixel:
-    default:
-    {
-      MagickPixelPacket
-        pixels[4];
-
-      MagickRealType
-        alpha[4],
-        gamma;
-
-      PointInfo
-        delta,
-        epsilon;
-
-      p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
-        (ssize_t) floor(y),2,2,resample_filter->exception);
-      if (p == (const PixelPacket *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
-      for (i=0; i < 4L; i++)
-      {
-        pixels[i].red=(MagickRealType) p[i].red;
-        pixels[i].green=(MagickRealType) p[i].green;
-        pixels[i].blue=(MagickRealType) p[i].blue;
-        pixels[i].opacity=(MagickRealType) p[i].opacity;
-        alpha[i]=1.0;
-      }
-      if (resample_filter->matte != MagickFalse)
-        for (i=0; i < 4L; i++)
-        {
-          alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p[i].opacity);
-          pixels[i].red*=alpha[i];
-          pixels[i].green*=alpha[i];
-          pixels[i].blue*=alpha[i];
-        }
-      if (indexes != (IndexPacket *) NULL)
-        for (i=0; i < 4L; i++)
-        {
-          pixels[i].index=(MagickRealType) indexes[i];
-          if (pixels[i].colorspace == CMYKColorspace)
-            pixels[i].index*=alpha[i];
-        }
-      delta.x=x-floor(x);
-      delta.y=y-floor(y);
-      epsilon.x=1.0-delta.x;
-      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=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : 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*
-        pixels[1].green)+delta.y*(epsilon.x*pixels[2].green+delta.x*
-        pixels[3].green));
-      pixel->blue=gamma*(epsilon.y*(epsilon.x*pixels[0].blue+delta.x*
-        pixels[1].blue)+delta.y*(epsilon.x*pixels[2].blue+delta.x*
-        pixels[3].blue));
-      pixel->opacity=(epsilon.y*(epsilon.x*pixels[0].opacity+delta.x*
-        pixels[1].opacity)+delta.y*(epsilon.x*pixels[2].opacity+delta.x*
-        pixels[3].opacity));
-      if (pixel->colorspace == CMYKColorspace)
-        pixel->index=gamma*(epsilon.y*(epsilon.x*pixels[0].index+delta.x*
-          pixels[1].index)+delta.y*(epsilon.x*pixels[2].index+delta.x*
-          pixels[3].index));
-      break;
-    }
-    case FilterInterpolatePixel:
-    {
-      CacheView
-        *filter_view;
-
-      Image
-        *excerpt_image,
-        *filter_image;
-
-      MagickPixelPacket
-        pixels[1];
-
-      RectangleInfo
-        geometry;
-
-      geometry.width=4L;
-      geometry.height=4L;
-      geometry.x=(ssize_t) floor(x)-1L;
-      geometry.y=(ssize_t) floor(y)-1L;
-      excerpt_image=ExcerptImage(resample_filter->image,&geometry,
-        resample_filter->exception);
-      if (excerpt_image == (Image *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      filter_image=ResizeImage(excerpt_image,1,1,resample_filter->image->filter,
-        resample_filter->image->blur,resample_filter->exception);
-      excerpt_image=DestroyImage(excerpt_image);
-      if (filter_image == (Image *) NULL)
-        break;
-      filter_view=AcquireCacheView(filter_image);
-      p=GetCacheViewVirtualPixels(filter_view,0,0,1,1,
-        resample_filter->exception);
-      if (p != (const PixelPacket *) NULL)
-        {
-          indexes=GetVirtualIndexQueue(filter_image);
-          GetMagickPixelPacket(resample_filter->image,pixels);
-          SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
-        }
-      filter_view=DestroyCacheView(filter_view);
-      filter_image=DestroyImage(filter_image);
-      break;
-    }
-    case IntegerInterpolatePixel:
-    {
-      MagickPixelPacket
-        pixels[1];
-
-      p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
-        (ssize_t) floor(y),1,1,resample_filter->exception);
-      if (p == (const PixelPacket *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
-      GetMagickPixelPacket(resample_filter->image,pixels);
-      SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
-      break;
-    }
-    case MeshInterpolatePixel:
-    {
-      MagickPixelPacket
-        pixels[4];
-
-      MagickRealType
-        alpha[4],
-        gamma;
-
-      PointInfo
-        delta,
-        luminance;
-
-      p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
-        (ssize_t) floor(y),2,2,resample_filter->exception);
-      if (p == (const PixelPacket *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
-      for (i=0; i < 4L; i++)
-      {
-        GetMagickPixelPacket(resample_filter->image,pixels+i);
-        SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
-        alpha[i]=1.0;
-        if (resample_filter->matte != MagickFalse)
-          {
-            alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
-            pixels[i].red*=alpha[i];
-            pixels[i].green*=alpha[i];
-            pixels[i].blue*=alpha[i];
-            if (pixels[i].colorspace == CMYKColorspace)
-              pixels[i].index*=alpha[i];
-          }
-        p++;
-      }
-      delta.x=x-floor(x);
-      delta.y=y-floor(y);
-      luminance.x=MagickPixelLuminance(pixels+0)-MagickPixelLuminance(pixels+3);
-      luminance.y=MagickPixelLuminance(pixels+1)-MagickPixelLuminance(pixels+2);
-      if (fabs(luminance.x) < fabs(luminance.y))
-        {
-          /*
-            Diagonal 0-3 NW-SE.
-          */
-          if (delta.x <= delta.y)
-            {
-              /*
-                Bottom-left triangle  (pixel:2, diagonal: 0-3).
-              */
-              delta.y=1.0-delta.y;
-              gamma=MeshInterpolate(&delta,alpha[2],alpha[3],alpha[0]);
-              gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-              pixel->red=gamma*MeshInterpolate(&delta,pixels[2].red,
-                pixels[3].red,pixels[0].red);
-              pixel->green=gamma*MeshInterpolate(&delta,pixels[2].green,
-                pixels[3].green,pixels[0].green);
-              pixel->blue=gamma*MeshInterpolate(&delta,pixels[2].blue,
-                pixels[3].blue,pixels[0].blue);
-              pixel->opacity=gamma*MeshInterpolate(&delta,pixels[2].opacity,
-                pixels[3].opacity,pixels[0].opacity);
-              if (pixel->colorspace == CMYKColorspace)
-                pixel->index=gamma*MeshInterpolate(&delta,pixels[2].index,
-                  pixels[3].index,pixels[0].index);
-            }
-          else
-            {
-              /*
-                Top-right triangle (pixel:1, diagonal: 0-3).
-              */
-              delta.x=1.0-delta.x;
-              gamma=MeshInterpolate(&delta,alpha[1],alpha[0],alpha[3]);
-              gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-              pixel->red=gamma*MeshInterpolate(&delta,pixels[1].red,
-                pixels[0].red,pixels[3].red);
-              pixel->green=gamma*MeshInterpolate(&delta,pixels[1].green,
-                pixels[0].green,pixels[3].green);
-              pixel->blue=gamma*MeshInterpolate(&delta,pixels[1].blue,
-                pixels[0].blue,pixels[3].blue);
-              pixel->opacity=gamma*MeshInterpolate(&delta,pixels[1].opacity,
-                pixels[0].opacity,pixels[3].opacity);
-              if (pixel->colorspace == CMYKColorspace)
-                pixel->index=gamma*MeshInterpolate(&delta,pixels[1].index,
-                  pixels[0].index,pixels[3].index);
-            }
-        }
-      else
-        {
-          /*
-            Diagonal 1-2 NE-SW.
-          */
-          if (delta.x <= (1.0-delta.y))
-            {
-              /*
-                Top-left triangle (pixel 0, diagonal: 1-2).
-              */
-              gamma=MeshInterpolate(&delta,alpha[0],alpha[1],alpha[2]);
-              gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-              pixel->red=gamma*MeshInterpolate(&delta,pixels[0].red,
-                pixels[1].red,pixels[2].red);
-              pixel->green=gamma*MeshInterpolate(&delta,pixels[0].green,
-                pixels[1].green,pixels[2].green);
-              pixel->blue=gamma*MeshInterpolate(&delta,pixels[0].blue,
-                pixels[1].blue,pixels[2].blue);
-              pixel->opacity=gamma*MeshInterpolate(&delta,pixels[0].opacity,
-                pixels[1].opacity,pixels[2].opacity);
-              if (pixel->colorspace == CMYKColorspace)
-                pixel->index=gamma*MeshInterpolate(&delta,pixels[0].index,
-                  pixels[1].index,pixels[2].index);
-            }
-          else
-            {
-              /*
-                Bottom-right triangle (pixel: 3, diagonal: 1-2).
-              */
-              delta.x=1.0-delta.x;
-              delta.y=1.0-delta.y;
-              gamma=MeshInterpolate(&delta,alpha[3],alpha[2],alpha[1]);
-              gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-              pixel->red=gamma*MeshInterpolate(&delta,pixels[3].red,
-                pixels[2].red,pixels[1].red);
-              pixel->green=gamma*MeshInterpolate(&delta,pixels[3].green,
-                pixels[2].green,pixels[1].green);
-              pixel->blue=gamma*MeshInterpolate(&delta,pixels[3].blue,
-                pixels[2].blue,pixels[1].blue);
-              pixel->opacity=gamma*MeshInterpolate(&delta,pixels[3].opacity,
-                pixels[2].opacity,pixels[1].opacity);
-              if (pixel->colorspace == CMYKColorspace)
-                pixel->index=gamma*MeshInterpolate(&delta,pixels[3].index,
-                  pixels[2].index,pixels[1].index);
-            }
-        }
-      break;
-    }
-    case NearestNeighborInterpolatePixel:
-    {
-      MagickPixelPacket
-        pixels[1];
-
-      p=GetCacheViewVirtualPixels(resample_filter->view,NearestNeighbor(x),
-        NearestNeighbor(y),1,1,resample_filter->exception);
-      if (p == (const PixelPacket *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
-      GetMagickPixelPacket(resample_filter->image,pixels);
-      SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
-      break;
-    }
-    case SplineInterpolatePixel:
-    {
-      MagickPixelPacket
-        pixels[16];
-
-      MagickRealType
-        alpha[16],
-        dx,
-        dy,
-        gamma;
-
-      PointInfo
-        delta;
-
-      ssize_t
-        j,
-        n;
-
-      p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
-        (ssize_t) floor(y)-1,4,4,resample_filter->exception);
-      if (p == (const PixelPacket *) NULL)
-        {
-          status=MagickFalse;
-          break;
-        }
-      indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
-      n=0;
-      delta.x=x-floor(x);
-      delta.y=y-floor(y);
-      for (i=(-1); i < 3L; i++)
-      {
-        dy=CubicWeightingFunction((MagickRealType) i-delta.y);
-        for (j=(-1); j < 3L; j++)
-        {
-          GetMagickPixelPacket(resample_filter->image,pixels+n);
-          SetMagickPixelPacket(resample_filter->image,p,indexes+n,pixels+n);
-          alpha[n]=1.0;
-          if (resample_filter->matte != MagickFalse)
-            {
-              alpha[n]=QuantumScale*((MagickRealType)
-                GetAlphaPixelComponent(p));
-              pixels[n].red*=alpha[n];
-              pixels[n].green*=alpha[n];
-              pixels[n].blue*=alpha[n];
-              if (pixels[i].colorspace == CMYKColorspace)
-                pixels[n].index*=alpha[n];
-            }
-          dx=CubicWeightingFunction(delta.x-(MagickRealType) j);
-          gamma=alpha[n];
-          gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-          pixel->red+=gamma*dx*dy*pixels[n].red;
-          pixel->green+=gamma*dx*dy*pixels[n].green;
-          pixel->blue+=gamma*dx*dy*pixels[n].blue;
-          if (resample_filter->matte != MagickFalse)
-            pixel->opacity+=dx*dy*pixels[n].opacity;
-          if (pixel->colorspace == CMYKColorspace)
-            pixel->index+=gamma*dx*dy*pixels[n].index;
-          n++;
-          p++;
-        }
-      }
-      break;
-    }
-  }
-  return(status);
-}
-\f
-/*
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                             %
-%                                                                             %
-%                                                                             %
 %   R e s a m p l e P i x e l C o l o r                                       %
 %                                                                             %
 %                                                                             %
@@ -906,15 +331,16 @@ MagickExport MagickBooleanType ResamplePixelColor(
   assert(resample_filter->signature == MagickSignature);
 
   status=MagickTrue;
-  GetMagickPixelPacket(resample_filter->image,pixel);
+  /* GetMagickPixelPacket(resample_filter->image,pixel); */
   if ( resample_filter->do_interpolate ) {
-    status=InterpolateResampleFilter(resample_filter,
-      resample_filter->interpolate,u0,v0,pixel);
+    status=InterpolateMagickPixelPacket(resample_filter->image,
+      resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
+      resample_filter->exception);
     return(status);
   }
 
 #if DEBUG_ELLIPSE
-  fprintf(stderr, "u0=%lf; v0=%lf;\n", u0, v0);
+  (void) FormatLocaleFile(stderr, "u0=%lf; v0=%lf;\n", u0, v0);
 #endif
 
   /*
@@ -985,8 +411,9 @@ MagickExport MagickBooleanType ResamplePixelColor(
   }
   if ( hit ) {
     /* whole area is a solid color -- just return that color */
-    status=InterpolateResampleFilter(resample_filter,IntegerInterpolatePixel,
-      u0,v0,pixel);
+    status=InterpolateMagickPixelPacket(resample_filter->image,
+      resample_filter->view,IntegerInterpolatePixel,u0,v0,pixel,
+      resample_filter->exception);
     return(status);
   }
 
@@ -1014,14 +441,16 @@ MagickExport MagickBooleanType ResamplePixelColor(
            works well in general, but falls down for TileEdge methods.
            This needs to be done properly!!!!!!
         */
-        status=InterpolateResampleFilter(resample_filter,
-          AverageInterpolatePixel,u0,v0,pixel);
+        status=InterpolateMagickPixelPacket(resample_filter->image,
+          resample_filter->view,AverageInterpolatePixel,u0,v0,pixel,
+          resample_filter->exception);
         break;
       case HorizontalTileVirtualPixelMethod:
       case VerticalTileVirtualPixelMethod:
         /* just return the background pixel - Is there more direct way? */
-        status=InterpolateResampleFilter(resample_filter,
-           IntegerInterpolatePixel,(double)-1,(double)-1,pixel);
+        status=InterpolateMagickPixelPacket(resample_filter->image,
+          resample_filter->view,IntegerInterpolatePixel,-1.0,-1.0,pixel,
+          resample_filter->exception);
         break;
       case TileVirtualPixelMethod:
       case MirrorVirtualPixelMethod:
@@ -1107,7 +536,7 @@ MagickExport MagickBooleanType ResamplePixelColor(
   divisor_c = 0.0;
   divisor_m = 0.0;
   pixel->red = pixel->green = pixel->blue = 0.0;
-  if (resample_filter->matte != MagickFalse) pixel->opacity = 0.0;
+  if (pixel->matte != MagickFalse) pixel->opacity = 0.0;
   if (pixel->colorspace == CMYKColorspace) pixel->index = 0.0;
 
   /*
@@ -1122,8 +551,8 @@ MagickExport MagickBooleanType ResamplePixelColor(
   uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
 
 #if DEBUG_ELLIPSE
-  fprintf(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
-  fprintf(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
+  (void) FormatLocaleFile(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
+  (void) FormatLocaleFile(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
 #else
 # define DEBUG_HIT_MISS 0 /* only valid if DEBUG_ELLIPSE is enabled */
 #endif
@@ -1136,7 +565,7 @@ MagickExport MagickBooleanType ResamplePixelColor(
   for( v=v1; v<=v2;  v++ ) {
 #if DEBUG_HIT_MISS
     long uu = ceil(u1);   /* actual pixel location (for debug only) */
-    fprintf(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
+    (void) FormatLocaleFile(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
 #endif
     u = (ssize_t)ceil(u1);        /* first pixel in scanline */
     u1 += resample_filter->slope; /* start of next scan line */
@@ -1173,7 +602,7 @@ MagickExport MagickBooleanType ResamplePixelColor(
         pixel->opacity  += weight*pixels->opacity;
         divisor_m += weight;
 
-        if (resample_filter->matte != MagickFalse)
+        if (pixel->matte != MagickFalse)
           weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
         pixel->red   += weight*pixels->red;
         pixel->green += weight*pixels->green;
@@ -1185,14 +614,14 @@ MagickExport MagickBooleanType ResamplePixelColor(
         hit++;
 #if DEBUG_HIT_MISS
         /* mark the pixel according to hit/miss of the ellipse */
-        fprintf(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
+        (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
-        fprintf(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
+        (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
       } else {
-        fprintf(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
+        (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
-        fprintf(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
+        (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
       }
       uu++;
@@ -1206,7 +635,7 @@ MagickExport MagickBooleanType ResamplePixelColor(
     }
   }
 #if DEBUG_ELLIPSE
-  fprintf(stderr, "Hit=%ld;  Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
+  (void) FormatLocaleFile(stderr, "Hit=%ld;  Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
 #endif
 
   /*
@@ -1218,8 +647,9 @@ MagickExport MagickBooleanType ResamplePixelColor(
     pixel->opacity = pixel->red = pixel->green = pixel->blue = 0;
     pixel->red = QuantumRange; /* show pixels for which EWA fails */
 #else
-    status=InterpolateResampleFilter(resample_filter,
-      resample_filter->interpolate,u0,v0,pixel);
+    status=InterpolateMagickPixelPacket(resample_filter->image,
+      resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
+      resample_filter->exception);
 #endif
     return status;
   }
@@ -1298,7 +728,7 @@ static inline void ClampUpAxes(const double dux,
    * position in input space and [x,y].)
    */
   /*
-   * Outputs:
+   * Output:
    *
    * major_mag is the half-length of the major axis of the "new"
    * ellipse.
@@ -1317,8 +747,17 @@ static inline void ClampUpAxes(const double dux,
    *
    * Unit vectors are useful for computing projections, in particular,
    * to compute the distance between a point in output space and the
-   * center (of a disk) from the position of the corresponding point
-   * in input space.
+   * center of a unit disk in output space, using the position of the
+   * corresponding point [s,t] in input space. Following the clamping,
+   * the square of this distance is
+   *
+   * ( ( s * major_unit_x + t * major_unit_y ) / major_mag )^2
+   * +
+   * ( ( s * minor_unit_x + t * minor_unit_y ) / minor_mag )^2
+   *
+   * If such distances will be computed for many [s,t]'s, it makes
+   * sense to actually compute the reciprocal of major_mag and
+   * minor_mag and multiply them by the above unit lengths.
    *
    * Now, if you want to modify the input pair of tangent vectors so
    * that it defines the modified ellipse, all you have to do is set
@@ -1376,7 +815,7 @@ static inline void ClampUpAxes(const double dux,
    * pulling pixel locations in the output image back to locations in
    * the input image. Jinv is
    *
-   *   Jinb = [ a, b ] = [ dx/dX, dx/dY ]
+   *   Jinv = [ a, b ] = [ dx/dX, dx/dY ]
    *          [ c, d ]   [ dy/dX, dy/dY ]
    *
    * Note: Jinv can be computed from J with the following matrix
@@ -1395,7 +834,7 @@ static inline void ClampUpAxes(const double dux,
    *
    * be an SVD decomposition of Jinv. (The SVD is not unique, but the
    * final ellipse does not depend on the particular SVD.)
-   * 
+   *
    * We could clamp up the entries of the diagonal matrix Sigma so
    * that they are at least 1, and then set
    *
@@ -1597,8 +1036,8 @@ MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
     return; /* EWA turned off - nothing to do */
 
 #if DEBUG_ELLIPSE
-  fprintf(stderr, "# -----\n" );
-  fprintf(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy=%lf;\n",
+  (void) FormatLocaleFile(stderr, "# -----\n" );
+  (void) FormatLocaleFile(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy=%lf;\n",
        dux, dvx, duy, dvy);
 #endif
 
@@ -1628,7 +1067,7 @@ MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
   major_x *= major_mag;  major_y *= major_mag;
   minor_x *= minor_mag;  minor_y *= minor_mag;
 #if DEBUG_ELLIPSE
-  fprintf(stderr, "major_x=%lf; major_y=%lf;  minor_x=%lf; minor_y=%lf;\n",
+  (void) FormatLocaleFile(stderr, "major_x=%lf; major_y=%lf;  minor_x=%lf; minor_y=%lf;\n",
         major_x, major_y, minor_x, minor_y);
 #endif
   A = major_y*major_y+minor_y*minor_y;
@@ -1667,7 +1106,7 @@ MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
 #endif
 
 #if DEBUG_ELLIPSE
-  fprintf(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
+  (void) FormatLocaleFile(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
 
   /* Figure out the various information directly about the ellipse.
      This information currently not needed at this time, but may be
@@ -1688,14 +1127,14 @@ MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
       Major = sqrt(2*F/(alpha - gamma));
     Minor = sqrt(2*F/(alpha + gamma));
 
-    fprintf(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
+    (void) FormatLocaleFile(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
 
     /* other information about ellipse include... */
     Eccentricity = Major/Minor;
     Ellipse_Area = MagickPI*Major*Minor;
     Ellipse_Angle = atan2(B, A-C);
 
-    fprintf(stderr, "# Angle=%lf   Area=%lf\n",
+    (void) FormatLocaleFile(stderr, "# Angle=%lf   Area=%lf\n",
          RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
   }
 #endif
@@ -1719,15 +1158,15 @@ MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
   F *= resample_filter->support;
 
   /* Orthogonal bounds of the ellipse */
-  resample_filter->Ulimit = sqrt(4*C*F/(4*A*C-B*B));
-  resample_filter->Vlimit = sqrt(4*A*F/(4*A*C-B*B));
+  resample_filter->Ulimit = sqrt(C*F/(A*C-0.25*B*B));
+  resample_filter->Vlimit = sqrt(A*F/(A*C-0.25*B*B));
 
   /* Horizontally aligned parallelogram fitted to Ellipse */
   resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
-  resample_filter->slope = -B/(2*A); /* Reciprocal slope of the parallelogram */
+  resample_filter->slope = -B/(2.0*A); /* Reciprocal slope of the parallelogram */
 
 #if DEBUG_ELLIPSE
-  fprintf(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
+  (void) FormatLocaleFile(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
            resample_filter->Ulimit, resample_filter->Vlimit,
            resample_filter->Uwidth, resample_filter->slope );
 #endif
@@ -1962,44 +1401,6 @@ MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%   S e t R e s a m p l e F i l t e r M a t t e                               %
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%  SetResampleFilterColorspace() sets the resample filter matte.
-%
-%  The format of the SetResampleFilterColorspace method is:
-%
-%      MagickBooleanType SetResampleFilterColorspace(
-%        ResampleFilter *resample_filter,const MagickBooleanType matte)
-%
-%  A description of each parameter follows:
-%
-%    o resample_filter: the resample filter.
-%
-%    o matte: the filter matte.
-%
-*/
-MagickExport MagickBooleanType SetResampleFilterMatte(
-  ResampleFilter *resample_filter,const MagickBooleanType matte)
-{
-  assert(resample_filter != (ResampleFilter *) NULL);
-  assert(resample_filter->signature == MagickSignature);
-  assert(resample_filter->image != (Image *) NULL);
-  if (resample_filter->debug != MagickFalse)
-    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
-      resample_filter->image->filename);
-  resample_filter->matte=matte;
-  return(MagickTrue);
-}
-\f
-/*
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                             %
-%                                                                             %
-%                                                                             %
 %   S e t R e s a m p l e F i l t e r V i r t u a l P i x e l M e t h o d     %
 %                                                                             %
 %                                                                             %