resample.c

Go to the documentation of this file.

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%           RRRR    EEEEE   SSSSS   AAA   M   M  PPPP   L      EEEEE          %
%           R   R   E       SS     A   A  MM MM  P   P  L      E              %
%           RRRR    EEE      SSS   AAAAA  M M M  PPPP   L      EEE            %
%           R R     E          SS  A   A  M   M  P      L      E              %
%           R  R    EEEEE   SSSSS  A   A  M   M  P      LLLLL  EEEEE          %
%                                                                             %
%                                                                             %
%                      MagickCore Pixel Resampling Methods                    %
%                                                                             %
%                              Software Design                                %
%                                John Cristy                                  %
%                              Anthony Thyssen                                %
%                                August 2007                                  %
%                                                                             %
%                                                                             %
%  Copyright 1999-2008 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  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    http://www.imagemagick.org/script/license.php                            %
%                                                                             %
%  Unless required by applicable law or agreed to in writing, software        %
%  distributed under the License is distributed on an "AS IS" BASIS,          %
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
%  See the License for the specific language governing permissions and        %
%  limitations under the License.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
*/

/*
  Include declarations.
*/
#include "magick/studio.h"
#include "magick/artifact.h"
#include "magick/color-private.h"
#include "magick/cache.h"
#include "magick/draw.h"
#include "magick/exception-private.h"
#include "magick/gem.h"
#include "magick/image.h"
#include "magick/image-private.h"
#include "magick/log.h"
#include "magick/memory_.h"
#include "magick/pixel-private.h"
#include "magick/quantum.h"
#include "magick/random_.h"
#include "magick/resample.h"
#include "magick/resize.h"
#include "magick/resize-private.h"
#include "magick/transform.h"
#include "magick/signature-private.h"
/*
  Typedef declarations.
*/
#define WLUT_WIDTH 1024
struct _ResampleFilter
{
  Image
    *image;

  ViewInfo
    *view;

  ExceptionInfo
    *exception;

  MagickBooleanType
    debug;

  /* Information about image being resampled */
  long
    image_area;

  InterpolatePixelMethod
    interpolate;

  VirtualPixelMethod
    virtual_pixel;

  FilterTypes
    filter;

  /* processing settings needed */
  MagickBooleanType
    limit_reached,
    do_interpolate,
    average_defined;

  MagickPixelPacket
    average_pixel;

  /* current ellipitical area being resampled around center point */
  double
    A, B, C,
    sqrtA, sqrtC, sqrtU, slope;

  /* LUT of weights for filtered average in elliptical area */
  double
    filter_lut[WLUT_WIDTH],
    support;

  unsigned long
    signature;
};

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   A c q u i r e R e s a m p l e I n f o                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AcquireResampleFilter() initializes the information resample needs do to a
%  scaled lookup of a color from an image, using area sampling.
%
%  The algorithm is based on a Elliptical Weighted Average, where the pixels
%  found in a large elliptical area is averaged together according to a
%  weighting (filter) function.  For more details see "Fundamentals of Texture
%  Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
%  1989.  Available for free from, http://www.cs.cmu.edu/~ph/
%
%  As EWA resampling (or any sort of resampling) can require a lot of
%  calculations to produce a distorted scaling of the source image for each
%  output pixel, the ResampleFilter structure generated holds that information
%  between individual image resampling.
%
%  This function will make the appropriate AcquireCacheView() calls
%  to view the image, calling functions do not need to open a cache view.
%
%  Usage Example...
%      resample_filter=AcquireResampleFilter(image,exception);
%      for (y=0; y < (long) image->rows; y++) {
%        for (x=0; x < (long) image->columns; x++) {
%          X= ....;   Y= ....;
%          ScaleResampleFilter(resample_filter, ... scaling vectors ...);
%          (void) ResamplePixelColor(resample_filter,X,Y,&pixel);
%          ... assign resampled pixel value ...
%        }
%      }
%      DestroyResampleFilter(resample_filter);
%
%  The format of the AcquireResampleFilter method is:
%
%     ResampleFilter *AcquireResampleFilter(const Image *image,
%       ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
%    o exception: return any errors or warnings in this structure.
%
*/
MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
  ExceptionInfo *exception)
{
  register ResampleFilter
    *resample_filter;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  assert(exception != (ExceptionInfo *) NULL);
  assert(exception->signature == MagickSignature);

  resample_filter=(ResampleFilter *) AcquireMagickMemory(
    sizeof(*resample_filter));
  if (resample_filter == (ResampleFilter *) NULL)
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
  (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));

  resample_filter->image=ReferenceImage((Image *) image);
  resample_filter->view=AcquireCacheView(resample_filter->image);
  resample_filter->exception=exception;

  resample_filter->debug=IsEventLogging();
  resample_filter->signature=MagickSignature;

  resample_filter->image_area = (long) resample_filter->image->columns *
    resample_filter->image->rows;
  resample_filter->average_defined = MagickFalse;

  /* initialise the resampling filter settings */
  SetResampleFilter(resample_filter, resample_filter->image->filter,
    resample_filter->image->blur);
  resample_filter->interpolate = resample_filter->image->interpolate;
  resample_filter->virtual_pixel=GetImageVirtualPixelMethod(image);

  /* init scale to a default of a unit circle */
  ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);

  return(resample_filter);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   D e s t r o y R e s a m p l e I n f o                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  DestroyResampleFilter() finalizes and cleans up the resampling
%  resample_filter as returned by AcquireResampleFilter(), freeing any memory
%  or other information as needed.
%
%  The format of the DestroyResampleFilter method is:
%
%      ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
%
%  A description of each parameter follows:
%
%    o resample_filter: resampling information structure
%
*/
MagickExport ResampleFilter *DestroyResampleFilter(
  ResampleFilter *resample_filter)
{
  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->view=DestroyCacheView(resample_filter->view);
  resample_filter->image=DestroyImage(resample_filter->image);
  resample_filter->signature=(~MagickSignature);
  resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
  return(resample_filter);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   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 long NearestNeighbor(MagickRealType x)
{
  if (x >= 0.0)
    return((long) (x+0.5));
  return((long) (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 long
    i;

  assert(resample_filter != (ResampleFilter *) NULL);
  assert(resample_filter->signature == MagickSignature);
  status=MagickTrue;
  switch (method)
  {
    case AverageInterpolatePixel:
    {
      MagickPixelPacket
        pixels[16];

      MagickRealType
        alpha[16],
        gamma;

      GetMagickPixelPacket(resample_filter->image,pixel);
      p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
        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->image->matte != MagickFalse)
          {
            alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
            pixels[i].red*=alpha[i];
            pixels[i].green*=alpha[i];
            pixels[i].blue*=alpha[i];
            if (resample_filter->image->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 (resample_filter->image->colorspace == CMYKColorspace)
          pixel->index+=gamma*0.0625*pixels[i].index;
        p++;
      }
      break;
    }
    case BicubicInterpolatePixel:
    {
      MagickPixelPacket
        pixels[16],
        u[4];

      MagickRealType
        alpha[16];

      PointInfo
        delta;

      GetMagickPixelPacket(resample_filter->image,pixel);
      p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
        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->image->matte != MagickFalse)
          {
            alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
            pixels[i].red*=alpha[i];
            pixels[i].green*=alpha[i];
            pixels[i].blue*=alpha[i];
            if (resample_filter->image->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;

      p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
        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->image->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 (resample_filter->image->colorspace == CMYKColorspace)
            pixels[i].index*=alpha[i];
        }
      delta.x=x-floor(x);
      delta.y=y-floor(y);
      gamma=(((1.0-delta.y)*((1.0-delta.x)*alpha[0]+delta.x*alpha[1])+delta.y*
        ((1.0-delta.x)*alpha[2]+delta.x*alpha[3])));
      gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
      pixel->red=gamma*((1.0-delta.y)*((1.0-delta.x)*pixels[0].red+delta.x*
        pixels[1].red)+delta.y*((1.0-delta.x)*pixels[2].red+delta.x*
        pixels[3].red));
      pixel->green=gamma*((1.0-delta.y)*((1.0-delta.x)*pixels[0].green+delta.x*
        pixels[1].green)+delta.y*((1.0-delta.x)*pixels[2].green+
        delta.x*pixels[3].green));
      pixel->blue=gamma*((1.0-delta.y)*((1.0-delta.x)*pixels[0].blue+delta.x*
        pixels[1].blue)+delta.y*((1.0-delta.x)*pixels[2].blue+delta.x*
        pixels[3].blue));
      pixel->opacity=((1.0-delta.y)*((1.0-delta.x)*pixels[0].opacity+delta.x*
        pixels[1].opacity)+delta.y*((1.0-delta.x)*pixels[2].opacity+delta.x*
        pixels[3].opacity));
      if (resample_filter->image->colorspace == CMYKColorspace)
        pixel->index=gamma*((1.0-delta.y)*((1.0-delta.x)*pixels[0].index+
          delta.x*pixels[1].index)+delta.y*((1.0-delta.x)*pixels[2].index+
          delta.x*pixels[3].index));
      break;
    }
    case FilterInterpolatePixel:
    {
      Image
        *excerpt_image,
        *filter_image;

      MagickPixelPacket
        pixels[1];

      RectangleInfo
        geometry;

      ViewInfo
        *filter_view;

      GetMagickPixelPacket(resample_filter->image,pixel);
      geometry.width=4L;
      geometry.height=4L;
      geometry.x=(long) floor(x)-1L;
      geometry.y=(long) 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];

      GetMagickPixelPacket(resample_filter->image,pixel);
      p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
        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;

      GetMagickPixelPacket(resample_filter->image,pixel);
      p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
        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->image->matte != MagickFalse)
          {
            alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
            pixels[i].red*=alpha[i];
            pixels[i].green*=alpha[i];
            pixels[i].blue*=alpha[i];
            if (resample_filter->image->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 (resample_filter->image->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 (resample_filter->image->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 (resample_filter->image->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 (resample_filter->image->colorspace == CMYKColorspace)
                pixel->index=gamma*MeshInterpolate(&delta,pixels[3].index,
                  pixels[2].index,pixels[1].index);
            }
        }
      break;
    }
    case NearestNeighborInterpolatePixel:
    {
      MagickPixelPacket
        pixels[1];

      GetMagickPixelPacket(resample_filter->image,pixel);
      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:
    {
      long
        j,
        n;

      MagickPixelPacket
        pixels[16];

      MagickRealType
        alpha[16],
        dx,
        dy,
        gamma;

      PointInfo
        delta;

      GetMagickPixelPacket(resample_filter->image,pixel);
      p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
        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->image->matte != MagickFalse)
            {
              alpha[n]=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
              pixels[n].red*=alpha[n];
              pixels[n].green*=alpha[n];
              pixels[n].blue*=alpha[n];
              if (resample_filter->image->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->image->matte != MagickFalse)
            pixel->opacity+=dx*dy*pixels[n].opacity;
          if (resample_filter->image->colorspace == CMYKColorspace)
            pixel->index+=gamma*dx*dy*pixels[n].index;
          n++;
          p++;
        }
      }
      break;
    }
  }
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   R e s a m p l e P i x e l C o l o r                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ResamplePixelColor() samples the pixel values surrounding the location
%  given using an elliptical weighted average, at the scale previously
%  calculated, and in the most efficent manner possible for the
%  VirtualPixelMethod setting.
%
%  The format of the ResamplePixelColor method is:
%
%     MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
%       const double u0,const double v0,MagickPixelPacket *pixel)
%
%  A description of each parameter follows:
%
%    o resample_filter: the resample filter.
%
%    o u0,v0: A double representing the center of the area to resample,
%        The distortion transformed transformed x,y coordinate.
%
%    o pixel: the resampled pixel is returned here.
%
*/
MagickExport MagickBooleanType ResamplePixelColor(
  ResampleFilter *resample_filter,const double u0,const double v0,
  MagickPixelPacket *pixel)
{
  MagickBooleanType
    status;

  long u,v, uw,v1,v2, hit;
  double u1;
  double U,V,Q,DQ,DDQ;
  double divisor_c,divisor_m;
  register double weight;
  register const PixelPacket *pixels;
  register const IndexPacket *indexes;
  assert(resample_filter != (ResampleFilter *) NULL);
  assert(resample_filter->signature == MagickSignature);

  if ( resample_filter->do_interpolate ) {
    status=InterpolateResampleFilter(resample_filter,
      resample_filter->interpolate,u0,v0,pixel);
    return(status);
  }

  /*
    Does resample area Miss the image?
    And is that area a simple solid color - then return that color
  */
  status=MagickTrue;
  GetMagickPixelPacket(resample_filter->image,pixel);
  hit = 0;
  switch ( resample_filter->virtual_pixel ) {
    case BackgroundVirtualPixelMethod:
    case ConstantVirtualPixelMethod:
    case TransparentVirtualPixelMethod:
    case BlackVirtualPixelMethod:
    case GrayVirtualPixelMethod:
    case WhiteVirtualPixelMethod:
    case MaskVirtualPixelMethod:
      if ( resample_filter->limit_reached
           || u0 + resample_filter->sqrtC < 0.0
           || u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
           || v0 + resample_filter->sqrtA < 0.0
           || v0 - resample_filter->sqrtA > (double) resample_filter->image->rows
           )
        hit++;
      break;

    case UndefinedVirtualPixelMethod:
    case EdgeVirtualPixelMethod:
      if (    ( u0 + resample_filter->sqrtC < 0.0 && v0 + resample_filter->sqrtA < 0.0 )
           || ( u0 + resample_filter->sqrtC < 0.0
                && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows )
           || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
                && v0 + resample_filter->sqrtA < 0.0 )
           || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
                && v0 - resample_