blob: ae998b6e3fe0200476fe52d1ce181c4edf406379 [file] [log] [blame]
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% %
% 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-2010 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;
CacheView
*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;
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;
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,
epsilon;
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);
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 (resample_filter->image->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:
{
Image
*excerpt_image,
*filter_image;
MagickPixelPacket
pixels[1];
RectangleInfo
geometry;
CacheView
*filter_view;
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];
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;
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];
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;
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);
status=MagickTrue;
GetMagickPixelPacket(resample_filter->image,pixel);
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
*/
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_filter->sqrtA > (double) resample_filter->image->rows )
)
hit++;
break;
case HorizontalTileVirtualPixelMethod:
if ( v0 + resample_filter->sqrtA < 0.0
|| v0 - resample_filter->sqrtA > (double) resample_filter->image->rows
)
hit++; /* outside the horizontally tiled images. */
break;
case VerticalTileVirtualPixelMethod:
if ( u0 + resample_filter->sqrtC < 0.0
|| u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
)
hit++; /* outside the vertically tiled images. */
break;
case DitherVirtualPixelMethod:
if ( ( u0 + resample_filter->sqrtC < -32.0 && v0 + resample_filter->sqrtA < -32.0 )
|| ( u0 + resample_filter->sqrtC < -32.0
&& v0 - resample_filter->sqrtA > (double) resample_filter->image->rows+32.0 )
|| ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns+32.0
&& v0 + resample_filter->sqrtA < -32.0 )
|| ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns+32.0
&& v0 - resample_filter->sqrtA > (double) resample_filter->image->rows+32.0 )
)
hit++;
break;
case TileVirtualPixelMethod:
case MirrorVirtualPixelMethod:
case RandomVirtualPixelMethod:
case HorizontalTileEdgeVirtualPixelMethod:
case VerticalTileEdgeVirtualPixelMethod:
case CheckerTileVirtualPixelMethod:
/* resampling of area is always needed - no VP limits */
break;
}
if ( hit ) {
/* whole area is a solid color -- just return that color */
status=InterpolateResampleFilter(resample_filter,IntegerInterpolatePixel,
u0,v0,pixel);
return(status);
}
/*
Scaling limits reached, return an 'averaged' result.
*/
if ( resample_filter->limit_reached ) {
switch ( resample_filter->virtual_pixel ) {
/* This is always handled by the above, so no need.
case BackgroundVirtualPixelMethod:
case ConstantVirtualPixelMethod:
case TransparentVirtualPixelMethod:
case GrayVirtualPixelMethod,
case WhiteVirtualPixelMethod
case MaskVirtualPixelMethod:
*/
case UndefinedVirtualPixelMethod:
case EdgeVirtualPixelMethod:
case DitherVirtualPixelMethod:
case HorizontalTileEdgeVirtualPixelMethod:
case VerticalTileEdgeVirtualPixelMethod:
/* We need an average edge pixel, for the right edge!
How should I calculate an average edge color?
Just returning an averaged neighbourhood,
works well in general, but falls down for TileEdge methods.
This needs to be done properly!!!!!!
*/
status=InterpolateResampleFilter(resample_filter,
AverageInterpolatePixel,u0,v0,pixel);
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);
break;
case TileVirtualPixelMethod:
case MirrorVirtualPixelMethod:
case RandomVirtualPixelMethod:
case CheckerTileVirtualPixelMethod:
default:
/* generate a average color of the WHOLE image */
if ( resample_filter->average_defined == MagickFalse ) {
Image
*average_image;
CacheView
*average_view;
GetMagickPixelPacket(resample_filter->image,
(MagickPixelPacket *)&(resample_filter->average_pixel));
resample_filter->average_defined = MagickTrue;
/* Try to get an averaged pixel color of whole image */
average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
resample_filter->exception);
if (average_image == (Image *) NULL)
{
*pixel=resample_filter->average_pixel; /* FAILED */
break;
}
average_view=AcquireCacheView(average_image);
pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
resample_filter->exception);
if (pixels == (const PixelPacket *) NULL) {
average_view=DestroyCacheView(average_view);
average_image=DestroyImage(average_image);
*pixel=resample_filter->average_pixel; /* FAILED */
break;
}
indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
SetMagickPixelPacket(resample_filter->image,pixels,indexes,
&(resample_filter->average_pixel));
average_view=DestroyCacheView(average_view);
average_image=DestroyImage(average_image);
#if 0
/* CheckerTile should average the image with background color */
//if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod ) {
#if 0
resample_filter->average_pixel.red =
( resample_filter->average_pixel.red +
resample_filter->image->background_color.red ) /2;
resample_filter->average_pixel.green =
( resample_filter->average_pixel.green +
resample_filter->image->background_color.green ) /2;
resample_filter->average_pixel.blue =
( resample_filter->average_pixel.blue +
resample_filter->image->background_color.blue ) /2;
resample_filter->average_pixel.matte =
( resample_filter->average_pixel.matte +
resample_filter->image->background_color.matte ) /2;
resample_filter->average_pixel.black =
( resample_filter->average_pixel.black +
resample_filter->image->background_color.black ) /2;
#else
resample_filter->average_pixel =
resample_filter->image->background_color;
#endif
}
#endif
}
*pixel=resample_filter->average_pixel;
break;
}
return(status);
}
/*
Initialize weighted average data collection
*/
hit = 0;
divisor_c = 0.0;
divisor_m = 0.0;
pixel->red = pixel->green = pixel->blue = 0.0;
if (resample_filter->image->matte != MagickFalse) pixel->opacity = 0.0;
if (resample_filter->image->colorspace == CMYKColorspace) pixel->index = 0.0;
/*
Determine the parellelogram bounding box fitted to the ellipse
centered at u0,v0. This area is bounding by the lines...
v = +/- sqrt(A)
u = -By/2A +/- sqrt(F/A)
Which has been pre-calculated above.
*/
v1 = (long)(v0 - resample_filter->sqrtA); /* range of scan lines */
v2 = (long)(v0 + resample_filter->sqrtA + 1);
u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->sqrtU; /* start of scanline for v=v1 */
uw = (long)(2*resample_filter->sqrtU)+1; /* width of parallelogram */
/*
Do weighted resampling of all pixels, within the scaled ellipse,
bound by a Parellelogram fitted to the ellipse.
*/
DDQ = 2*resample_filter->A;
for( v=v1; v<=v2; v++, u1+=resample_filter->slope ) {
u = (long)u1; /* first pixel in scanline ( floor(u1) ) */
U = (double)u-u0; /* location of that pixel, relative to u0,v0 */
V = (double)v-v0;
/* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
Q = U*(resample_filter->A*U + resample_filter->B*V) + resample_filter->C*V*V;
DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
/* get the scanline of pixels for this v */
pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(unsigned long) uw,
1,resample_filter->exception);
if (pixels == (const PixelPacket *) NULL)
return(MagickFalse);
indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
/* count up the weighted pixel colors */
for( u=0; u<uw; u++ ) {
/* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
if ( Q < (double)WLUT_WIDTH ) {
weight = resample_filter->filter_lut[(int)Q];
pixel->opacity += weight*pixels->opacity;
divisor_m += weight;
if (resample_filter->image->matte != MagickFalse)
weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
pixel->red += weight*pixels->red;
pixel->green += weight*pixels->green;
pixel->blue += weight*pixels->blue;
if (resample_filter->image->colorspace == CMYKColorspace)
pixel->index += weight*(*indexes);
divisor_c += weight;
hit++;
}
pixels++;
indexes++;
Q += DQ;
DQ += DDQ;
}
}
/*
Result sanity check -- this should NOT happen
*/
if ( hit < 4 || divisor_c < 1.0 ) {
/* not enough pixels in resampling, resort to direct interpolation */
status=InterpolateResampleFilter(resample_filter,
resample_filter->interpolate,u0,v0,pixel);
return status;
}
/*
Finialize results of resampling
*/
divisor_m = 1.0/divisor_m;
pixel->opacity = (MagickRealType) RoundToQuantum(divisor_m*pixel->opacity);
divisor_c = 1.0/divisor_c;
pixel->red = (MagickRealType) RoundToQuantum(divisor_c*pixel->red);
pixel->green = (MagickRealType) RoundToQuantum(divisor_c*pixel->green);
pixel->blue = (MagickRealType) RoundToQuantum(divisor_c*pixel->blue);
if (resample_filter->image->colorspace == CMYKColorspace)
pixel->index = (MagickRealType) RoundToQuantum(divisor_c*pixel->index);
return(MagickTrue);
}
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% %
% S c a l e R e s a m p l e F i l t e r %
% %
% %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% ScaleResampleFilter() does all the calculations needed to resample an image
% at a specific scale, defined by two scaling vectors. This not using
% a orthogonal scaling, but two distorted scaling vectors, to allow the
% generation of a angled ellipse.
%
% As only two deritive scaling vectors are used the center of the ellipse
% must be the center of the lookup. That is any curvature that the
% distortion may produce is discounted.
%
% The input vectors are produced by either finding the derivitives of the
% distortion function, or the partial derivitives from a distortion mapping.
% They do not need to be the orthogonal dx,dy scaling vectors, but can be
% calculated from other derivatives. For example you could use dr,da/r
% polar coordinate vector scaling vectors
%
% If u,v = DistortEquation(x,y)
% Then the scaling vectors dx,dy (in u,v space) are the derivitives...
% du/dx, dv/dx and du/dy, dv/dy
% If the scaling is only othogonally aligned then...
% dv/dx = 0 and du/dy = 0
% Producing an othogonally alligned ellipse for the area to be resampled.
%
% Note that scaling vectors are different to argument order. Argument order
% is the general order the deritives are extracted from the distortion
% equations, EG: U(x,y), V(x,y). Caution is advised if you are trying to
% define the ellipse directly from scaling vectors.
%
% The format of the ScaleResampleFilter method is:
%
% void ScaleResampleFilter(const ResampleFilter *resample_filter,
% const double dux,const double duy,const double dvx,const double dvy)
%
% A description of each parameter follows:
%
% o resample_filter: the resampling resample_filterrmation defining the
% image being resampled
%
% o dux,duy,dvx,dvy:
% The partial derivitives or scaling vectors for resampling.
% dx = du/dx, dv/dx and dy = du/dy, dv/dy
%
% The values are used to define the size and angle of the
% elliptical resampling area, centered on the lookup point.
%
*/
MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
const double dux,const double duy,const double dvx,const double dvy)
{
double A,B,C,F, area;
assert(resample_filter != (ResampleFilter *) NULL);
assert(resample_filter->signature == MagickSignature);
resample_filter->limit_reached = MagickFalse;
resample_filter->do_interpolate = MagickFalse;
/* A 'point' filter forces use of interpolation instead of area sampling */
if ( resample_filter->filter == PointFilter ) {
resample_filter->do_interpolate = MagickTrue;
return;
}
/* Find Ellipse Coefficents such that
A*u^2 + B*u*v + C*v^2 = F
With u,v relative to point around which we are resampling.
And the given scaling dx,dy vectors in u,v space
du/dx,dv/dx and du/dy,dv/dy
*/
#if 0
/* Direct conversions of derivatives to elliptical coefficients
No scaling will result in F == 1.0 and a unit circle.
*/
A = dvx*dvx+dvy*dvy;
B = (dux*dvx+duy*dvy)*-2.0;
C = dux*dux+duy*duy;
F = dux*dvy+duy*dvx;
F *= F;
#define F_UNITY 1.0
#else
/* This Paul Heckbert's recomended "Higher Quality EWA" formula, from page
60 in his thesis, which adds a unit circle to the elliptical area so are
to do both Reconstruction and Prefiltering of the pixels in the
resampling. It also means it is likely to have at least 4 pixels within
the area of the ellipse, for weighted averaging.
No scaling will result if F == 4.0 and a circle of radius 2.0
*/
A = dvx*dvx+dvy*dvy+1;
B = (dux*dvx+duy*dvy)*-2.0;
C = dux*dux+duy*duy+1;
F = A*C - B*B/4;
#define F_UNITY 4.0
#endif
/* DEBUGGING OUTPUT */
#if 0
fprintf(stderr, "dux=%lf; dvx=%lf; duy=%lf; dvy%lf;\n",
dux, dvx, duy, dvy);
fprintf(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
#endif
#if 0
/* Figure out the Ellipses Major and Minor Axis, and other info.
This information currently not needed at this time, but may be
needed later for better limit determination.
*/
{ double alpha, beta, gamma, Major, Minor;
double Eccentricity, Ellipse_Area, Ellipse_angle;
double max_horizontal_cross_section, max_vertical_cross_section;
alpha = A+C;
beta = A-C;
gamma = sqrt(beta*beta + B*B );
if ( alpha - gamma <= MagickEpsilon )
Major = MagickHuge;
else
Major = sqrt(2*F/(alpha - gamma));
Minor = sqrt(2*F/(alpha + gamma));
fprintf(stderr, "\tMajor=%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, "\tAngle=%lf Area=%lf\n",
RadiansToDegrees(Ellipse_angle), Ellipse_Area );
/* Ellipse limits */
/* orthogonal rectangle - improved ellipse */
max_horizontal_orthogonal = sqrt(A); /* = sqrt(4*A*F/(4*A*C-B*B)) */
max_vertical_orthogonal = sqrt(C); /* = sqrt(4*C*F/(4*A*C-B*B)) */
/* parallelogram bounds -- what we are using */
max_horizontal_cross_section = sqrt(F/A);
max_vertical_cross_section = sqrt(F/C);
}
#endif
/* Is default elliptical area, too small? Image being magnified?
Switch to doing pure 'point' interpolation of the pixel.
That is turn off EWA Resampling.
*/
if ( F <= F_UNITY ) {
resample_filter->do_interpolate = MagickTrue;
return;
}
/* If F is impossibly large, we may as well not bother doing any
* form of resampling, as you risk an infinite resampled area.
*/
if ( F > MagickHuge ) {
resample_filter->limit_reached = MagickTrue;
return;
}
/* Othogonal bounds of the ellipse */
resample_filter->sqrtA = sqrt(A)+1.0; /* Vertical Orthogonal Limit */
resample_filter->sqrtC = sqrt(C)+1.0; /* Horizontal Orthogonal Limit */
/* Horizontally aligned Parallelogram fitted to ellipse */
resample_filter->sqrtU = sqrt(F/A)+1.0; /* Parallelogram Width */
resample_filter->slope = -B/(2*A); /* Slope of the parallelogram */
/* The size of the area of the parallelogram we will be sampling */
area = 4 * resample_filter->sqrtA * resample_filter->sqrtU;
/* Absolute limit on the area to be resampled
* This limit needs more work, as it gets too slow for
* larger images involved with tiled views of the horizon. */
if ( area > 20.0*resample_filter->image_area ) {
resample_filter->limit_reached = MagickTrue;
return;
}
/* Scale ellipse formula to directly fit the Filter Lookup Table */
{ register double scale;
scale = (double)WLUT_WIDTH/F;
resample_filter->A = A*scale;
resample_filter->B = B*scale;
resample_filter->C = C*scale;
/* ..ple_filter->F = WLUT_WIDTH; -- hardcoded */
}
}
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% %
% S e t R e s a m p l e F i l t e r %
% %
% %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SetResampleFilter() set the resampling filter lookup table based on a
% specific filter. Note that the filter is used as a radial filter not as a
% two pass othogonally aligned resampling filter.
%
% The default Filter, is Gaussian, which is the standard filter used by the
% original paper on the Elliptical Weighted Everage Algorithm. However other
% filters can also be used.
%
% The format of the SetResampleFilter method is:
%
% void SetResampleFilter(ResampleFilter *resample_filter,
% const FilterTypes filter,const double blur)
%
% A description of each parameter follows:
%
% o resample_filter: resampling resample_filterrmation structure
%
% o filter: the resize filter for elliptical weighting LUT
%
% o blur: filter blur factor (radial scaling) for elliptical weighting LUT
%
*/
MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
const FilterTypes filter,const double blur)
{
register int
Q;
double
r_scale;
ResizeFilter
*resize_filter;
assert(resample_filter != (ResampleFilter *) NULL);
assert(resample_filter->signature == MagickSignature);
resample_filter->filter = filter;
/* Scale radius so it equals 1.0, at edge of ellipse when a
default blurring factor of 1.0 is used.
Note that these filters are being used as a radial filter, not as
an othoginally alligned filter. How this effects results is still
being worked out.
Future: Direct use of teh resize filters in "resize.c" to set the lookup
table, based on the filters working support window.
*/
r_scale = sqrt(1.0/(double)WLUT_WIDTH)/blur;
r_scale *= 2; /* for 2 pixel radius of Improved Elliptical Formula */
switch ( filter ) {
case PointFilter:
/* This equivelent to turning off the EWA algroithm.
Only Interpolated lookup will be used. */
break;
default:
/*
Fill the LUT with a 1D resize filter function
But make the Sinc/Bessel tapered window 2.0
I also normalize the result so the filter is 1.0
*/
resize_filter = AcquireResizeFilter(resample_filter->image,filter,
(MagickRealType)1.0,MagickTrue,resample_filter->exception);
if (resize_filter != (ResizeFilter *) NULL) {
resample_filter->support = GetResizeFilterSupport(resize_filter);
resample_filter->support /= blur; /* taken into account above */
resample_filter->support *= resample_filter->support;
resample_filter->support *= (double)WLUT_WIDTH/4;
if ( resample_filter->support >= (double)WLUT_WIDTH )
resample_filter->support = (double)WLUT_WIDTH; /* hack */
for(Q=0; Q<WLUT_WIDTH; Q++)
if ( (double) Q < resample_filter->support )
resample_filter->filter_lut[Q] = (double)
GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
else
resample_filter->filter_lut[Q] = 0.0;
resize_filter = DestroyResizeFilter(resize_filter);
break;
}
else {
(void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
ModuleError, "UnableToSetFilteringValue",
"Fall back to default EWA gaussian filter");
}
/* FALLTHRU - on exception */
/*case GaussianFilter:*/
case UndefinedFilter:
/*
Create Normal Gaussian 2D Filter Weighted Lookup Table.
A normal EWA guassual lookup would use exp(Q*ALPHA)
where Q = distantce squared from 0.0 (center) to 1.0 (edge)
and ALPHA = -4.0*ln(2.0) ==> -2.77258872223978123767
However the table is of length 1024, and equates to a radius of 2px
thus needs to be scaled by ALPHA*4/1024 and any blur factor squared
*/
/*r_scale = -2.77258872223978123767*4/WLUT_WIDTH/blur/blur;*/
r_scale = -2.77258872223978123767/WLUT_WIDTH/blur/blur;
for(Q=0; Q<WLUT_WIDTH; Q++)
resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
resample_filter->support = WLUT_WIDTH;
break;
}
if (GetImageArtifact(resample_filter->image,"resample:verbose")
!= (const char *) NULL)
/* Debug output of the filter weighting LUT
Gnuplot the LUT with hoizontal adjusted to 'r' using...
plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
THe filter values is normalized for comparision
*/
for(Q=0; Q<WLUT_WIDTH; Q++)
printf("%lf\n", resample_filter->filter_lut[Q]
/resample_filter->filter_lut[0] );
return;
}
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% %
% S e t R e s a m p l e F i l t e r I n t e r p o l a t e M e t h o d %
% %
% %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SetResampleFilterInterpolateMethod() changes the interpolation method
% associated with the specified resample filter.
%
% The format of the SetResampleFilterInterpolateMethod method is:
%
% MagickBooleanType SetResampleFilterInterpolateMethod(
% ResampleFilter *resample_filter,const InterpolateMethod method)
%
% A description of each parameter follows:
%
% o resample_filter: the resample filter.
%
% o method: the interpolation method.
%
*/
MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
ResampleFilter *resample_filter,const InterpolatePixelMethod method)
{
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->interpolate=method;
return(MagickTrue);
}
/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% %
% %
% 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 %
% %
% %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
% associated with the specified resample filter.
%
% The format of the SetResampleFilterVirtualPixelMethod method is:
%
% MagickBooleanType SetResampleFilterVirtualPixelMethod(
% ResampleFilter *resample_filter,const VirtualPixelMethod method)
%
% A description of each parameter follows:
%
% o resample_filter: the resample filter.
%
% o method: the virtual pixel method.
%
*/
MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
ResampleFilter *resample_filter,const VirtualPixelMethod method)
{
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->virtual_pixel=method;
(void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
return(MagickTrue);
}