/*****************************************************************************
 *	Date: June 19, 2008
 *----------------------------------------------------------------------------
 *	This C program is based on the following paper:
 *		T. Blu, P. Thevenaz, M. Unser,
 *			"Linear Interpolation Revitalized,"
 *			IEEE Transactions on Image Processing, vol. 13, no. 5, pp. 710-719,
 * 			May 2004.
 *----------------------------------------------------------------------------
 *	Philippe Thevenaz
 *	EPFL/STI/IMT/BIG/BM.4.137
 *	Station 17
 *	CH-1015 Lausanne VD
 *  Switzerland
 *----------------------------------------------------------------------------
 *	phone (CET):	+41(21)693.51.61
 *	fax:			+41(21)693.37.01
 *	RFC-822:		philippe.thevenaz@epfl.ch
 *	X-400:			/C=ch/A=400net/P=switch/O=epfl/S=thevenaz/G=philippe/
 *	URL:			http://bigwww.epfl.ch/
 *----------------------------------------------------------------------------
 *	This file is best viewed with 4-space tabs
 * (the bars separated by dots should be aligned with those separated by tabs)
 *  |...|...|...|...|...|...|...|...|...|...|...|...|...|...|...|...|...|...|
 *	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|
 ****************************************************************************/

/*****************************************************************************
 *	System includes
 ****************************************************************************/
#include	<math.h>
#include	<stddef.h>
#include	<stdio.h>
#include	<stdlib.h>

/*****************************************************************************
 *	Other includes
 ****************************************************************************/
#include	"coeff.h"
#include	"interpol.h"
 
/*****************************************************************************
 *	Definition of extern procedures
 ****************************************************************************/
/*--------------------------------------------------------------------------*/
extern double	InterpolatedValue
				(
					float	*LinCoeff,	/* input array of shifted-linear coefficients */
					long	Width,		/* width of the image */
					long	Height,		/* height of the image */
					double	x,			/* x coordinate where to interpolate */
					double	y			/* y coordinate where to interpolate */
				)

{ /* begin InterpolatedValue */

	float	*p;
	double	xWeight[2], yWeight[2];
	double	interpolated;
	double	w;
	long	xIndex[2], yIndex[2];
	long	i, j, k;

	/* perform the shift */
	x -= BASIS_SHIFT;
	y -= BASIS_SHIFT;

	/* compute the interpolation indexes */
	i = (long)floor(x);
	j = (long)floor(y);
	for (k = 0L; k <= 1; k++) {
		xIndex[k] = i++;
		yIndex[k] = j++;
	}

	/* compute the interpolation weights */
	/* x */
	w = x - (double)xIndex[0];
	xWeight[0] = 1.0 - w;
	xWeight[1] = w;
	/* y */
	w = y - (double)yIndex[0];
	yWeight[0] = 1.0 - w;
	yWeight[1] = w;

	/* check the boundaries */
	for (k = 0L; k <= 1; k++) {
		xIndex[k] = (xIndex[k] < 0L) ? (0L)
			: ((Width <= xIndex[k]) ? (Width - 1L) : (xIndex[k]));
		yIndex[k] = (yIndex[k] < 0L) ? (0L)
			: ((Height <= yIndex[k]) ? (Height - 1L) : (yIndex[k]));
	}

	/* perform interpolation */
	interpolated = 0.0;
	for (j = 0L; j <= 1; j++) {
		p = LinCoeff + (ptrdiff_t)(yIndex[j] * Width);
		w = 0.0;
		for (i = 0L; i <= 1; i++) {
			w += xWeight[i] * p[xIndex[i]];
		}
		interpolated += yWeight[j] * w;
	}

	return(interpolated);
} /* end InterpolatedValue */