/*****************************************************************************
 *	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	<stdio.h>
#include	<stdlib.h>

/*****************************************************************************
 *	Other includes
 ****************************************************************************/
#include	"coeff.h"
#include	"interpol.h"
#include	"io.h"
 
/*****************************************************************************
 *	Defines
 ****************************************************************************/
#define		PI	((double)3.14159265358979323846264338327950288419716939937510)

/*****************************************************************************
 *	Definition of extern procedures
 ****************************************************************************/
/*--------------------------------------------------------------------------*/
extern int		main
				(
					void
				)

{ /* begin main */

	float	*ImageRasterArray, *OutputImage;
	float	*p;
	double	a11, a12, a21, a22;
	double	x0, y0, x1, y1;
	double	xOrigin, yOrigin;
	double	Angle, xShift, yShift;
	long	Width, Height;
	long	x, y;
	int		Error;

	/* access the data samples */
	Error = ReadByteImageRawData(&ImageRasterArray, &Width, &Height);
	if (Error) {
		printf("Failure to import image data\n");
		return(1);
	}

	/* ask for the transformation parameters */
	RigidBody(&Angle, &xShift, &yShift, &xOrigin, &yOrigin);

	/* allocate the output image */
	OutputImage = (float *)malloc((size_t)(Width * Height
		* (long)sizeof(float)));
	if (OutputImage == (float *)NULL) {
		free(ImageRasterArray);
		printf("Allocation of output image failed\n");
		return(1);
	}

	/* convert between a representation based on image samples */
	/* and a representation based on shifted-linear coefficients */
	Error = SamplesToCoefficients(ImageRasterArray, Width, Height);
	if (Error) {
		free(OutputImage);
		free(ImageRasterArray);
		printf("Change of basis failed\n");
		return(1);
	}

	/* prepare the geometry */
	Angle *= PI / 180.0;
	a11 = cos(Angle);
	a12 = -sin(Angle);
	a21 = sin(Angle);
	a22 = cos(Angle);
	x0 = a11 * (xShift + xOrigin) + a12 * (yShift + yOrigin);
	y0 = a21 * (xShift + xOrigin) + a22 * (yShift + yOrigin);
	xShift = xOrigin - x0;
	yShift = yOrigin - y0;

	/* visit all pixels of the output image and assign their value */
	p = OutputImage;
	for (y = 0L; y < Height; y++) {
		x0 = a12 * (double)y + xShift;
		y0 = a22 * (double)y + yShift;
		for (x = 0L; x < Width; x++) {
			x1 = x0 + a11 * (double)x;
			y1 = y0 + a21 * (double)x;
			if ((x1 <= -0.5) || (((double)Width - 0.5) <= x1)
				|| (y1 <= -0.5) || (((double)Height - 0.5) <= y1)) {
				*p++ = 0.0F;
			}
			else {
				*p++ = (float)InterpolatedValue(ImageRasterArray,
					Width, Height, x1, y1);
			}
		}
	}

	/* save output */
	Error = WriteByteImageRawData(OutputImage, Width, Height);
	if (Error) {
		free(OutputImage);
		free(ImageRasterArray);
		printf("Failure to export image data\n");
		return(1);
	}

	free(OutputImage);
	free(ImageRasterArray);
	printf("Done\n");
	return(0);
} /* end main */