/*****************************************************************************
 *	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"
 
/*****************************************************************************
 *	Declaration of static procedures
 ****************************************************************************/
/*--------------------------------------------------------------------------*/
static void		GetColumn
				(
					float	*Image,		/* input image array */
					long	Width,		/* width of the image */
					long	x,			/* x coordinate of the selected line */
					double	Line[],		/* output linear array */
					long	Height		/* length of the line and height of the image */
				);

/*--------------------------------------------------------------------------*/
static void		GetRow
				(
					float	*Image,		/* input image array */
					long	y,			/* y coordinate of the selected line */
					double	Line[],		/* output linear array */
					long	Width		/* length of the line and width of the image */
				);

/*--------------------------------------------------------------------------*/
static void		PutColumn
				(
					float	*Image,		/* output image array */
					long	Width,		/* width of the image */
					long	x,			/* x coordinate of the selected line */
					double	Line[],		/* input linear array */
					long	Height		/* length of the line and height of the image */
				);

/*--------------------------------------------------------------------------*/
static void		PutRow
				(
					float	*Image,		/* output image array */
					long	y,			/* y coordinate of the selected line */
					double	Line[],		/* input linear array */
					long	Width		/* length of the line and width of the image */
				);

/*--------------------------------------------------------------------------*/
static void		ShiftedBasisCoefficients1D
				(
					double	c[],		/* in-place processing; 1D array */
					long	DataLength	/* length of the 1D array */
				);

/*--------------------------------------------------------------------------*/
static double	ShiftedBasisFirstCoefficient1D
				(
					double	c[]			/* input 1D array */
				);

/*****************************************************************************
 *	Definition of static procedures
 ****************************************************************************/
/*--------------------------------------------------------------------------*/
static void		GetColumn
				(
					float	*Image,		/* input image array */
					long	Width,		/* width of the image */
					long	x,			/* x coordinate of the selected line */
					double	Line[],		/* output linear array */
					long	Height		/* length of the line */
				)

{ /* begin GetColumn */

	long	y;

	Image = Image + (ptrdiff_t)x;
	for (y = 0L; (y < Height); y++) {
		Line[y] = (double)*Image;
		Image += (ptrdiff_t)Width;
	}
} /* end GetColumn */

/*--------------------------------------------------------------------------*/
static void		GetRow
				(
					float	*Image,		/* input image array */
					long	y,			/* y coordinate of the selected line */
					double	Line[],		/* output linear array */
					long	Width		/* length of the line */
				)

{ /* begin GetRow */

	long	x;

	Image = Image + (ptrdiff_t)(y * Width);
	for (x = 0L; (x < Width); x++) {
		Line[x] = (double)*Image++;
	}
} /* end GetRow */

/*--------------------------------------------------------------------------*/
static void		PutColumn
				(
					float	*Image,		/* output image array */
					long	Width,		/* width of the image */
					long	x,			/* x coordinate of the selected line */
					double	Line[],		/* input linear array */
					long	Height		/* length of the line and height of the image */
				)

{ /* begin PutColumn */

	long	y;

	Image = Image + (ptrdiff_t)x;
	for (y = 0L; (y < Height); y++) {
		*Image = (float)Line[y];
		Image += (ptrdiff_t)Width;
	}
} /* end PutColumn */

/*--------------------------------------------------------------------------*/
static void		PutRow
				(
					float	*Image,		/* output image array */
					long	y,			/* y coordinate of the selected line */
					double	Line[],		/* input linear array */
					long	Width		/* length of the line and width of the image */
				)

{ /* begin PutRow */

	long	x;

	Image = Image + (ptrdiff_t)(y * Width);
	for (x = 0L; (x < Width); x++) {
		*Image++ = (float)Line[x];
	}
} /* end PutRow */

/*--------------------------------------------------------------------------*/
static void		ShiftedBasisCoefficients1D
				(
					double	c[],		/* in-place processing; 1D array */
					long	DataLength	/* length of the 1D array */
				)

{ /* begin ShiftedBasisCoefficients1D */

	double	Factor = 1.0 / (1.0 - BASIS_SHIFT);
	double	Pole = BASIS_SHIFT / (BASIS_SHIFT - 1.0);
	long	n;

	c[0] = ShiftedBasisFirstCoefficient1D(c);
	for (n = 1L; (n < DataLength); n++) {
		c[n] = Factor * c[n] + Pole * c[n - 1L];
	}
} /* end ShiftedBasisCoefficients1D */

/*--------------------------------------------------------------------------*/
static double	ShiftedBasisFirstCoefficient1D
				(
					double	c[]			/* input 1D array */
				)

{ /* begin ShiftedBasisFirstCoefficient1D */

	return(c[0]);
} /* end ShiftedBasisFirstCoefficient1D */

/*****************************************************************************
 *	Definition of extern procedures
 ****************************************************************************/
/*--------------------------------------------------------------------------*/
extern int		SamplesToCoefficients
				(
					float	*Image,		/* in-place processing */
					long	Width,		/* width of the image */
					long	Height		/* height of the image */
				)

{ /* begin SamplesToCoefficients */

	double	*Line;
	long	x, y;

	/* convert the image samples into interpolation coefficients */
	/* in-place separable process, along x */
	Line = (double *)malloc((size_t)(Width * (long)sizeof(double)));
	if (Line == (double *)NULL) {
		printf("Row allocation failed\n");
		return(1);
	}
	for (y = 0L; y < Height; y++) {
		GetRow(Image, y, Line, Width);
		ShiftedBasisCoefficients1D(Line, Width);
		PutRow(Image, y, Line, Width);
	}
	free(Line);

	/* in-place separable process, along y */
	Line = (double *)malloc((size_t)(Height * (long)sizeof(double)));
	if (Line == (double *)NULL) {
		printf("Column allocation failed\n");
		return(1);
	}
	for (x = 0L; x < Width; x++) {
		GetColumn(Image, Width, x, Line, Height);
		ShiftedBasisCoefficients1D(Line, Height);
		PutColumn(Image, Width, x, Line, Height);
	}
	free(Line);

	return(0);
} /* end SamplesToCoefficients */