/*

 * The two pass scaling function is based on:

 * Filtered Image Rescaling

 * Based on Gems III

 *  - Schumacher general filtered image rescaling

 * (pp. 414-424)

 * by Dale Schumacher

 *

 * 	Additional changes by Ray Gardener, Daylon Graphics Ltd.

 * 	December 4, 1999

 *

 * 	Ported to libgd by Pierre Joye. Support for multiple channels

 * 	added (argb for now).

 *

 * 	Initial sources code is avaibable in the Gems Source Code Packages:

 * 	http://www.acm.org/pubs/tog/GraphicsGems/GGemsIII.tar.gz

 *

 */

/*

	Summary:

		- Horizontal filter contributions are calculated on the fly,

		  as each column is mapped from src to dst image. This lets

		  us omit having to allocate a temporary full horizontal stretch

		  of the src image.

		- If none of the src pixels within a sampling region differ,

		  then the output pixel is forced to equal (any of) the source pixel.

		  This ensures that filters do not corrupt areas of constant color.

		- Filter weight contribution results, after summing, are

		  rounded to the nearest pixel color value instead of

		  being casted to ILubyte (usually an int or char). Otherwise,

		  artifacting occurs.

*/

/*

	Additional functions are available for simple rotation or up/downscaling.

	downscaling using the fixed point implementations are usually much faster

	than the existing gdImageCopyResampled while having a similar or better

	quality.

	For image rotations, the optimized versions have a lazy antialiasing for 

	the edges of the images. For a much better antialiased result, the affine

	function is recommended.

*/

/*

TODO:

 - Optimize pixel accesses and loops once we have continuous buffer

 - Add scale support for a portion only of an image (equivalent of copyresized/resampled)

 */

#ifdef HAVE_CONFIG_H

#include "config.h"

#endif /* HAVE_CONFIG_H */

#include <stdio.h>

#include <stdlib.h>

#include <string.h>

#include <math.h>

#undef NDEBUG

/* Comment out this line to enable asserts.

 * TODO: This logic really belongs in cmake and configure.

 */

#define NDEBUG 1

#include <assert.h>

#include "gd.h"

#include "gdhelpers.h"

#include "gd_intern.h"

#ifdef _MSC_VER

# pragma optimize("t", on)

# include <emmintrin.h>

#endif

static gdImagePtr gdImageScaleBilinear(gdImagePtr im, 

                                       const unsigned int new_width,

                                       const unsigned int new_height);

static gdImagePtr gdImageScaleBicubicFixed(gdImagePtr src,

                                           const unsigned int width,

                                           const unsigned int height);

static gdImagePtr gdImageScaleNearestNeighbour(gdImagePtr im,

                                               const unsigned int width,

                                               const unsigned int height);

static gdImagePtr gdImageRotateNearestNeighbour(gdImagePtr src,

                                                const float degrees,

                                                const int bgColor);

static gdImagePtr gdImageRotateGeneric(gdImagePtr src, const float degrees,

                                       const int bgColor);

#define CLAMP(x, low, high)  (((x) > (high)) ? (high) : (((x) < (low)) ? (low) : (x)))

/* only used here, let do a generic fixed point integers later if required by other

   part of GD */

typedef long gdFixed;

/* Integer to fixed point */

#define gd_itofx(x) ((x) << 8)

/* Float to fixed point */

#define gd_ftofx(x) (long)((x) * 256)

/*  Double to fixed point */

#define gd_dtofx(x) (long)((x) * 256)

/* Fixed point to integer */

#define gd_fxtoi(x) ((x) >> 8)

/* Fixed point to float */

# define gd_fxtof(x) ((float)(x) / 256)

/* Fixed point to double */

#define gd_fxtod(x) ((double)(x) / 256)

/* Multiply a fixed by a fixed */

#define gd_mulfx(x,y) (((x) * (y)) >> 8)

/* Divide a fixed by a fixed */

#define gd_divfx(x,y) (((x) << 8) / (y))

typedef struct

{

	double *Weights;  /* Normalized weights of neighboring pixels */

	int Left,Right;   /* Bounds of source pixels window */

} ContributionType;  /* Contirbution information for a single pixel */

typedef struct

{

	ContributionType *ContribRow; /* Row (or column) of contribution weights */

	unsigned int WindowSize,      /* Filter window size (of affecting source pixels) */

		     LineLength;      /* Length of line (no. or rows / cols) */

} LineContribType;

/* Each core filter has its own radius */

#define DEFAULT_FILTER_LINEAR               1.0f

#define DEFAULT_FILTER_BICUBIC				3.0f

#define DEFAULT_FILTER_BOX					0.5f

#define DEFAULT_FILTER_GENERALIZED_CUBIC	0.5f

#define DEFAULT_FILTER_RADIUS				1.0f

#define DEFAULT_LANCZOS8_RADIUS				8.0f

#define DEFAULT_LANCZOS3_RADIUS				3.0f

#define DEFAULT_HERMITE_RADIUS				1.0f

#define DEFAULT_BOX_RADIUS					0.5f

#define DEFAULT_TRIANGLE_RADIUS				1.0f

#define DEFAULT_BELL_RADIUS					1.5f

#define DEFAULT_CUBICSPLINE_RADIUS			2.0f

#define DEFAULT_MITCHELL_RADIUS				2.0f

#define DEFAULT_COSINE_RADIUS				1.0f

#define DEFAULT_CATMULLROM_RADIUS			2.0f

#define DEFAULT_QUADRATIC_RADIUS			1.5f

#define DEFAULT_QUADRATICBSPLINE_RADIUS		1.5f

#define DEFAULT_CUBICCONVOLUTION_RADIUS		3.0f

#define DEFAULT_GAUSSIAN_RADIUS				1.0f

#define DEFAULT_HANNING_RADIUS				1.0f

#define DEFAULT_HAMMING_RADIUS				1.0f

#define DEFAULT_SINC_RADIUS					1.0f

#define DEFAULT_WELSH_RADIUS				1.0f

static double KernelBessel_J1(const double x)

{

	double p, q;

	register long i;

	static const double

	Pone[] =

	{

		0.581199354001606143928050809e+21,

		-0.6672106568924916298020941484e+20,

		0.2316433580634002297931815435e+19,

		-0.3588817569910106050743641413e+17,

		0.2908795263834775409737601689e+15,

		-0.1322983480332126453125473247e+13,

		0.3413234182301700539091292655e+10,

		-0.4695753530642995859767162166e+7,

		0.270112271089232341485679099e+4

	},

	Qone[] =

	{

		0.11623987080032122878585294e+22,

		0.1185770712190320999837113348e+20,

		0.6092061398917521746105196863e+17,

		0.2081661221307607351240184229e+15,

		0.5243710262167649715406728642e+12,

		0.1013863514358673989967045588e+10,

		0.1501793594998585505921097578e+7,

		0.1606931573481487801970916749e+4,

		0.1e+1

	};

	p = Pone[8];

	q = Qone[8];

	for (i=7; i >= 0; i--)

	{

		p = p*x*x+Pone[i];

		q = q*x*x+Qone[i];

	}

	return (double)(p/q);

}

static double KernelBessel_P1(const double x)

{

	double p, q;

	register long i;

	static const double

	Pone[] =

	{

		0.352246649133679798341724373e+5,

		0.62758845247161281269005675e+5,

		0.313539631109159574238669888e+5,

		0.49854832060594338434500455e+4,

		0.2111529182853962382105718e+3,

		0.12571716929145341558495e+1

	},

	Qone[] =

	{

		0.352246649133679798068390431e+5,

		0.626943469593560511888833731e+5,

		0.312404063819041039923015703e+5,

		0.4930396490181088979386097e+4,

		0.2030775189134759322293574e+3,

		0.1e+1

	};

	p = Pone[5];

	q = Qone[5];

	for (i=4; i >= 0; i--)

	{

		p = p*(8.0/x)*(8.0/x)+Pone[i];

		q = q*(8.0/x)*(8.0/x)+Qone[i];

	}

	return (double)(p/q);

}

static double KernelBessel_Q1(const double x)

{

	double p, q;

	register long i;

	static const double

	Pone[] =

	{

		0.3511751914303552822533318e+3,

		0.7210391804904475039280863e+3,

		0.4259873011654442389886993e+3,

		0.831898957673850827325226e+2,

		0.45681716295512267064405e+1,

		0.3532840052740123642735e-1

	},

	Qone[] =

	{

		0.74917374171809127714519505e+4,

		0.154141773392650970499848051e+5,

		0.91522317015169922705904727e+4,

		0.18111867005523513506724158e+4,

		0.1038187585462133728776636e+3,

		0.1e+1

	};

	p = Pone[5];

	q = Qone[5];

	for (i=4; i >= 0; i--)

	{

		p = p*(8.0/x)*(8.0/x)+Pone[i];

		q = q*(8.0/x)*(8.0/x)+Qone[i];

	}

	return (double)(p/q);

}

static double KernelBessel_Order1(double x)

{

	double p, q;

	if (x == 0.0)

		return (0.0f);

	p = x;

	if (x < 0.0)

		x=(-x);

	if (x < 8.0)

		return (p*KernelBessel_J1(x));

	q = (double)sqrt(2.0f/(M_PI*x))*(double)(KernelBessel_P1(x)*(1.0f/sqrt(2.0f)*(sin(x)-cos(x)))-8.0f/x*KernelBessel_Q1(x)*

		(-1.0f/sqrt(2.0f)*(sin(x)+cos(x))));

	if (p < 0.0f)

		q = (-q);

	return (q);

}

static double filter_bessel(const double x)

{

	if (x == 0.0f)

		return (double)(M_PI/4.0f);

	return (KernelBessel_Order1((double)M_PI*x)/(2.0f*x));

}

static double filter_blackman(const double x)

{

	return (0.42f+0.5f*(double)cos(M_PI*x)+0.08f*(double)cos(2.0f*M_PI*x));

}

double filter_linear(const double x) {

	double ax = fabs(x);

	if (ax < 1.0f) {

		return (1.0f - ax);

	}

	return 0.0f;

}

/**

 * Bicubic interpolation kernel (a=-1):

  \verbatim

          /

         | 1-2|t|**2+|t|**3          , if |t| < 1

  h(t) = | 4-8|t|+5|t|**2-|t|**3     , if 1<=|t|<2

         | 0                         , otherwise

          \

  \endverbatim

 * ***bd*** 2.2004

 */

static double filter_bicubic(const double t)

{

	const double abs_t = (double)fabs(t);

	const double abs_t_sq = abs_t * abs_t;

	if (abs_t<1) return 1-2*abs_t_sq+abs_t_sq*abs_t;

	if (abs_t<2) return 4 - 8*abs_t +5*abs_t_sq - abs_t_sq*abs_t;

	return 0;

}

/**

 * Generalized cubic kernel (for a=-1 it is the same as BicubicKernel):

  \verbatim

          /

         | (a+2)|t|**3 - (a+3)|t|**2 + 1     , |t| <= 1

  h(t) = | a|t|**3 - 5a|t|**2 + 8a|t| - 4a   , 1 < |t| <= 2

         | 0                                 , otherwise

          \

  \endverbatim

 * Often used values for a are -1 and -1/2.

 */

static double filter_generalized_cubic(const double t)

{

	const double a = -DEFAULT_FILTER_GENERALIZED_CUBIC;

	double abs_t = (double)fabs(t);

	double abs_t_sq = abs_t * abs_t;

	if (abs_t < 1) return (a + 2) * abs_t_sq * abs_t - (a + 3) * abs_t_sq + 1;

	if (abs_t < 2) return a * abs_t_sq * abs_t - 5 * a * abs_t_sq + 8 * a * abs_t - 4 * a;

	return 0;

}

#ifdef FUNCTION_NOT_USED_YET

/* CubicSpline filter, default radius 2 */

static double filter_cubic_spline(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

	if (x < 1.0 ) {

		const double x2 = x*x;

		return (0.5 * x2 * x - x2 + 2.0 / 3.0);

	}

	if (x < 2.0) {

		return (pow(2.0 - x, 3.0)/6.0);

	}

	return 0;

}

#endif

#ifdef FUNCTION_NOT_USED_YET

/* CubicConvolution filter, default radius 3 */

static double filter_cubic_convolution(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

	const double x2 = x1 * x1;

	const double x2_x = x2 * x;

	if (x <= 1.0) return ((4.0 / 3.0)* x2_x - (7.0 / 3.0) * x2 + 1.0);

	if (x <= 2.0) return (- (7.0 / 12.0) * x2_x + 3 * x2 - (59.0 / 12.0) * x + 2.5);

	if (x <= 3.0) return ( (1.0/12.0) * x2_x - (2.0 / 3.0) * x2 + 1.75 * x - 1.5);

	return 0;

}

#endif

static double filter_box(double x) {

	if (x < - DEFAULT_FILTER_BOX)

		return 0.0f;

	if (x < DEFAULT_FILTER_BOX)

		return 1.0f;

	return 0.0f;

}

static double filter_catmullrom(const double x)

{

	if (x < -2.0)

		return(0.0f);

	if (x < -1.0)

		return(0.5f*(4.0f+x*(8.0f+x*(5.0f+x))));

	if (x < 0.0)

		return(0.5f*(2.0f+x*x*(-5.0f-3.0f*x)));

	if (x < 1.0)

		return(0.5f*(2.0f+x*x*(-5.0f+3.0f*x)));

	if (x < 2.0)

		return(0.5f*(4.0f+x*(-8.0f+x*(5.0f-x))));

	return(0.0f);

}

#ifdef FUNCTION_NOT_USED_YET

static double filter_filter(double t)

{

	/* f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1 */

	if(t < 0.0) t = -t;

	if(t < 1.0) return((2.0 * t - 3.0) * t * t + 1.0);

	return(0.0);

}

#endif

#ifdef FUNCTION_NOT_USED_YET

/* Lanczos8 filter, default radius 8 */

static double filter_lanczos8(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

#define R DEFAULT_LANCZOS8_RADIUS

	if ( x == 0.0) return 1;

	if ( x < R) {

		return R * sin(x*M_PI) * sin(x * M_PI/ R) / (x * M_PI * x * M_PI);

	}

	return 0.0;

#undef R

}

#endif

#ifdef FUNCTION_NOT_USED_YET

/* Lanczos3 filter, default radius 3 */

static double filter_lanczos3(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

#define R DEFAULT_LANCZOS3_RADIUS

	if ( x == 0.0) return 1;

	if ( x < R)

	{

		return R * sin(x*M_PI) * sin(x * M_PI / R) / (x * M_PI * x * M_PI);

	}

	return 0.0;

#undef R

}

#endif

/* Hermite filter, default radius 1 */

static double filter_hermite(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

	if (x < 1.0) return ((2.0 * x - 3) * x * x + 1.0 );

	return 0.0;

}

/* Trangle filter, default radius 1 */

static double filter_triangle(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

	if (x < 1.0) return (1.0 - x);

	return 0.0;

}

/* Bell filter, default radius 1.5 */

static double filter_bell(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

	if (x < 0.5) return (0.75 - x*x);

	if (x < 1.5) return (0.5 * pow(x - 1.5, 2.0));

	return 0.0;

}

/* Mitchell filter, default radius 2.0 */

static double filter_mitchell(const double x)

{

#define KM_B (1.0f/3.0f)

#define KM_C (1.0f/3.0f)

#define KM_P0 ((  6.0f - 2.0f * KM_B ) / 6.0f)

#define KM_P2 ((-18.0f + 12.0f * KM_B + 6.0f * KM_C) / 6.0f)

#define KM_P3 (( 12.0f - 9.0f  * KM_B - 6.0f * KM_C) / 6.0f)

#define KM_Q0 ((  8.0f * KM_B + 24.0f * KM_C) / 6.0f)

#define KM_Q1 ((-12.0f * KM_B - 48.0f * KM_C) / 6.0f)

#define KM_Q2 ((  6.0f * KM_B + 30.0f * KM_C) / 6.0f)

#define KM_Q3 (( -1.0f * KM_B -  6.0f * KM_C) / 6.0f)

	if (x < -2.0)

		return(0.0f);

	if (x < -1.0)

		return(KM_Q0-x*(KM_Q1-x*(KM_Q2-x*KM_Q3)));

	if (x < 0.0f)

		return(KM_P0+x*x*(KM_P2-x*KM_P3));

	if (x < 1.0f)

		return(KM_P0+x*x*(KM_P2+x*KM_P3));

	if (x < 2.0f)

		return(KM_Q0+x*(KM_Q1+x*(KM_Q2+x*KM_Q3)));

	return(0.0f);

}

#ifdef FUNCTION_NOT_USED_YET

/* Cosine filter, default radius 1 */

static double filter_cosine(const double x)

{

	if ((x >= -1.0) && (x <= 1.0)) return ((cos(x * M_PI) + 1.0)/2.0);

	return 0;

}

#endif

/* Quadratic filter, default radius 1.5 */

static double filter_quadratic(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

	if (x <= 0.5) return (- 2.0 * x * x + 1);

	if (x <= 1.5) return (x * x - 2.5* x + 1.5);

	return 0.0;

}

static double filter_bspline(const double x)

{

	if (x>2.0f) {

		return 0.0f;

	} else {

		double a, b, c, d;

		/* Was calculated anyway cause the "if((x-1.0f) < 0)" */

		const double xm1 = x - 1.0f;

		const double xp1 = x + 1.0f;

		const double xp2 = x + 2.0f;

		if ((xp2) <= 0.0f) a = 0.0f; else a = xp2*xp2*xp2;

		if ((xp1) <= 0.0f) b = 0.0f; else b = xp1*xp1*xp1;

		if (x <= 0) c = 0.0f; else c = x*x*x;

		if ((xm1) <= 0.0f) d = 0.0f; else d = xm1*xm1*xm1;

		return (0.16666666666666666667f * (a - (4.0f * b) + (6.0f * c) - (4.0f * d)));

	}

}

#ifdef FUNCTION_NOT_USED_YET

/* QuadraticBSpline filter, default radius 1.5 */

static double filter_quadratic_bspline(const double x1)

{

	const double x = x1 < 0.0 ? -x1 : x1;

	if (x <= 0.5) return (- x * x + 0.75);

	if (x <= 1.5) return (0.5 * x * x - 1.5 * x + 1.125);

	return 0.0;

}

#endif

static double filter_gaussian(const double x)

{

	/* return(exp((double) (-2.0 * x * x)) * sqrt(2.0 / M_PI)); */

	return (double)(exp(-2.0f * x * x) * 0.79788456080287f);

}

static double filter_hanning(const double x)

{

	/* A Cosine windowing function */

	return(0.5 + 0.5 * cos(M_PI * x));

}

static double filter_hamming(const double x)

{

	/* should be

	(0.54+0.46*cos(M_PI*(double) x));

	but this approximation is sufficient */

	if (x < -1.0f)

		return 0.0f;

	if (x < 0.0f)

		return 0.92f*(-2.0f*x-3.0f)*x*x+1.0f;

	if (x < 1.0f)

		return 0.92f*(2.0f*x-3.0f)*x*x+1.0f;

	return 0.0f;

}

static double filter_power(const double x)

{

	const double a = 2.0f;

	if (fabs(x)>1) return 0.0f;

	return (1.0f - (double)fabs(pow(x,a)));

}

static double filter_sinc(const double x)

{

	/* X-scaled Sinc(x) function. */

	if (x == 0.0) return(1.0);

	return (sin(M_PI * (double) x) / (M_PI * (double) x));

}

#ifdef FUNCTION_NOT_USED_YET

static double filter_welsh(const double x)

{

	/* Welsh parabolic windowing filter */

	if (x <  1.0)

		return(1 - x*x);

	return(0.0);

}

#endif

#if defined(_MSC_VER) && !defined(inline)

# define inline __inline

#endif 

/* keep it for future usage for affine copy over an existing image, targetting fix for 2.2.2 */

#ifdef FUNCTION_NOT_USED_YET

/* Copied from upstream's libgd */

static inline int _color_blend (const int dst, const int src)

{

	const int src_alpha = gdTrueColorGetAlpha(src);

	if( src_alpha == gdAlphaOpaque ) {

		return src;

	} else {

		const int dst_alpha = gdTrueColorGetAlpha(dst);

		if( src_alpha == gdAlphaTransparent ) return dst;

		if( dst_alpha == gdAlphaTransparent ) {

			return src;

		} else {

			register int alpha, red, green, blue;

			const int src_weight = gdAlphaTransparent - src_alpha;

			const int dst_weight = (gdAlphaTransparent - dst_alpha) * src_alpha / gdAlphaMax;

			const int tot_weight = src_weight + dst_weight;

			alpha = src_alpha * dst_alpha / gdAlphaMax;

			red = (gdTrueColorGetRed(src) * src_weight

				   + gdTrueColorGetRed(dst) * dst_weight) / tot_weight;

			green = (gdTrueColorGetGreen(src) * src_weight

				   + gdTrueColorGetGreen(dst) * dst_weight) / tot_weight;

			blue = (gdTrueColorGetBlue(src) * src_weight

				   + gdTrueColorGetBlue(dst) * dst_weight) / tot_weight;

			return ((alpha << 24) + (red << 16) + (green << 8) + blue);

		}

	}

}

static inline int _setEdgePixel(const gdImagePtr src, unsigned int x, unsigned int y, gdFixed coverage, const int bgColor) 

{

	const gdFixed f_127 = gd_itofx(127);

	register int c = src->tpixels[y][x];

	c = c | (( (int) (gd_fxtof(gd_mulfx(coverage, f_127)) + 50.5f)) << 24);

	return _color_blend(bgColor, c);

}

#endif

static inline int getPixelOverflowTC(gdImagePtr im, const int x, const int y, const int bgColor)

{

	if (gdImageBoundsSafe(im, x, y)) {

		const int c = im->tpixels[y][x];

		if (c == im->transparent) {

			return bgColor == -1 ? gdTrueColorAlpha(0, 0, 0, 127) : bgColor;

		}

		return c;

	} else {

		return bgColor;

	}

}

#define colorIndex2RGBA(c) gdTrueColorAlpha(im->red[(c)], im->green[(c)], im->blue[(c)], im->alpha[(c)])

#define colorIndex2RGBcustomA(c, a) gdTrueColorAlpha(im->red[(c)], im->green[(c)], im->blue[(c)], im->alpha[(a)])

static inline int getPixelOverflowPalette(gdImagePtr im, const int x, const int y, const int bgColor)

{

	if (gdImageBoundsSafe(im, x, y)) {

		const int c = im->pixels[y][x];

		if (c == im->transparent) {

			return bgColor == -1 ? gdTrueColorAlpha(0, 0, 0, 127) : bgColor;

		}

		return colorIndex2RGBA(c);

	} else {

		return bgColor;

	}

}

static int getPixelInterpolateWeight(gdImagePtr im, const double x, const double y, const int bgColor)

{

	/* Closest pixel <= (xf,yf) */

	int sx = (int)(x);

	int sy = (int)(y);

	const double xf = x - (double)sx;

	const double yf = y - (double)sy;

	const double nxf = (double) 1.0 - xf;

	const double nyf = (double) 1.0 - yf;

	const double m1 = xf * yf;

	const double m2 = nxf * yf;

	const double m3 = xf * nyf;

	const double m4 = nxf * nyf;

	/* get color values of neighbouring pixels */

	const int c1 = im->trueColor == 1 ? getPixelOverflowTC(im, sx, sy, bgColor)         : getPixelOverflowPalette(im, sx, sy, bgColor);

	const int c2 = im->trueColor == 1 ? getPixelOverflowTC(im, sx - 1, sy, bgColor)     : getPixelOverflowPalette(im, sx - 1, sy, bgColor);

	const int c3 = im->trueColor == 1 ? getPixelOverflowTC(im, sx, sy - 1, bgColor)     : getPixelOverflowPalette(im, sx, sy - 1, bgColor);

	const int c4 = im->trueColor == 1 ? getPixelOverflowTC(im, sx - 1, sy - 1, bgColor) : getPixelOverflowPalette(im, sx, sy - 1, bgColor);

	int r, g, b, a;

	if (x < 0) sx--;

	if (y < 0) sy--;

	/* component-wise summing-up of color values */

	if (im->trueColor) {

		r = (int)(m1*gdTrueColorGetRed(c1)   + m2*gdTrueColorGetRed(c2)   + m3*gdTrueColorGetRed(c3)   + m4*gdTrueColorGetRed(c4));

		g = (int)(m1*gdTrueColorGetGreen(c1) + m2*gdTrueColorGetGreen(c2) + m3*gdTrueColorGetGreen(c3) + m4*gdTrueColorGetGreen(c4));

		b = (int)(m1*gdTrueColorGetBlue(c1)  + m2*gdTrueColorGetBlue(c2)  + m3*gdTrueColorGetBlue(c3)  + m4*gdTrueColorGetBlue(c4));

		a = (int)(m1*gdTrueColorGetAlpha(c1) + m2*gdTrueColorGetAlpha(c2) + m3*gdTrueColorGetAlpha(c3) + m4*gdTrueColorGetAlpha(c4));

	} else {

		r = (int)(m1*im->red[(c1)]   + m2*im->red[(c2)]   + m3*im->red[(c3)]   + m4*im->red[(c4)]);

		g = (int)(m1*im->green[(c1)] + m2*im->green[(c2)] + m3*im->green[(c3)] + m4*im->green[(c4)]);

		b = (int)(m1*im->blue[(c1)]  + m2*im->blue[(c2)]  + m3*im->blue[(c3)]  + m4*im->blue[(c4)]);

		a = (int)(m1*im->alpha[(c1)] + m2*im->alpha[(c2)] + m3*im->alpha[(c3)] + m4*im->alpha[(c4)]);

	}

	r = CLAMP(r, 0, 255);

	g = CLAMP(g, 0, 255);

	b = CLAMP(b, 0, 255);

	a = CLAMP(a, 0, gdAlphaMax);

	return gdTrueColorAlpha(r, g, b, a);

}

/**

 * InternalFunction: getPixelInterpolated

 *  Returns the interpolated color value using the default interpolation

 *  method. The returned color is always in the ARGB format (truecolor).

 *

 * Parameters:

 * 	im - Image to set the default interpolation method

 *  y - X value of the ideal position

 *  y - Y value of the ideal position

 *  method - Interpolation method <gdInterpolationMethod>

 *

 * Returns:

 *  GD_TRUE if the affine is rectilinear or GD_FALSE

 *

 * See also:

 *  <gdSetInterpolationMethod>

 */

int getPixelInterpolated(gdImagePtr im, const double x, const double y, const int bgColor)

{

	const int xi=(int)(x);

	const int yi=(int)(y);

	int yii;

	int i;

	double kernel, kernel_cache_y;

	double kernel_x[12], kernel_y[4];

	double new_r = 0.0f, new_g = 0.0f, new_b = 0.0f, new_a = 0.0f;

	/* These methods use special implementations */

	if (im->interpolation_id == GD_NEAREST_NEIGHBOUR) {

		return -1;

	}

	if (im->interpolation_id == GD_WEIGHTED4) {

		return getPixelInterpolateWeight(im, x, y, bgColor);

	}

	if (im->interpolation_id == GD_NEAREST_NEIGHBOUR) {

		if (im->trueColor == 1) {

			return getPixelOverflowTC(im, xi, yi, bgColor);

		} else {

			return getPixelOverflowPalette(im, xi, yi, bgColor);

		}

	}

	if (im->interpolation) {

		for (i=0; i<4; i++) {

			kernel_x[i] = (double) im->interpolation((double)(xi+i-1-x));

			kernel_y[i] = (double) im->interpolation((double)(yi+i-1-y));

		}

	} else {

		return -1;

	}

	/*

	 * TODO: use the known fast rgba multiplication implementation once

	 * the new formats are in place

	 */

	for (yii = yi-1; yii < yi+3; yii++) {

		int xii;

		kernel_cache_y = kernel_y[yii-(yi-1)];

		if (im->trueColor) {

			for (xii=xi-1; xii

				const int rgbs = getPixelOverflowTC(im, xii, yii, bgColor);

				kernel = kernel_cache_y * kernel_x[xii-(xi-1)];

				new_r += kernel * gdTrueColorGetRed(rgbs);

				new_g += kernel * gdTrueColorGetGreen(rgbs);

				new_b += kernel * gdTrueColorGetBlue(rgbs);

				new_a += kernel * gdTrueColorGetAlpha(rgbs);

			}

		} else {

			for (xii=xi-1; xii

				const int rgbs = getPixelOverflowPalette(im, xii, yii, bgColor);

				kernel = kernel_cache_y * kernel_x[xii-(xi-1)];

				new_r += kernel * gdTrueColorGetRed(rgbs);

				new_g += kernel * gdTrueColorGetGreen(rgbs);

				new_b += kernel * gdTrueColorGetBlue(rgbs);

				new_a += kernel * gdTrueColorGetAlpha(rgbs);

			}

		}

	}

	new_r = CLAMP(new_r, 0, 255);

	new_g = CLAMP(new_g, 0, 255);

	new_b = CLAMP(new_b, 0, 255);

	new_a = CLAMP(new_a, 0, gdAlphaMax);

	return gdTrueColorAlpha(((int)new_r), ((int)new_g), ((int)new_b), ((int)new_a));

}

static inline LineContribType * _gdContributionsAlloc(unsigned int line_length, unsigned int windows_size)

{

	unsigned int u = 0;

	LineContribType *res;

	size_t weights_size;

	if (overflow2(windows_size, sizeof(double))) {

		return NULL;

	} else {

		weights_size = windows_size * sizeof(double);

	}

	res = (LineContribType *) gdMalloc(sizeof(LineContribType));

	if (!res) {

		return NULL;

	}

	res->WindowSize = windows_size;

	res->LineLength = line_length;

	if (overflow2(line_length, sizeof(ContributionType))) {

		gdFree(res);

		return NULL;

	}

	res->ContribRow = (ContributionType *) gdMalloc(line_length * sizeof(ContributionType));

	if (res->ContribRow == NULL) {

		gdFree(res);

		return NULL;

	}

	for (u = 0 ; u < line_length ; u++) {

		res->ContribRow[u].Weights = (double *) gdMalloc(weights_size);

		if (res->ContribRow[u].Weights == NULL) {

			unsigned int i;

			for (i=0;i

				gdFree(res->ContribRow[i].Weights);

			}

			gdFree(res->ContribRow);

			gdFree(res);

			return NULL;

		}

	}

	return res;

}

static inline void _gdContributionsFree(LineContribType * p)

{

	unsigned int u;

	for (u = 0; u < p->LineLength; u++)  {

		gdFree(p->ContribRow[u].Weights);

	}

	gdFree(p->ContribRow);

	gdFree(p);

}

static inline LineContribType *_gdContributionsCalc(unsigned int line_size, unsigned int src_size, double scale_d,  const interpolation_method pFilter)