/*****************************************************************************
 *	Date: November 28th 2008  /  modified by P. Chainais
 *----------------------------------------------------------------------------
 *	This C program is based on the following paper:
 *		P. Thevenaz, T. Blu, M. Unser, "Interpolation Revisited,"
 *		IEEE Transactions on Medical Imaging,
 *		vol. 19, no. 7, pp. 739-758, July 2000.
 *----------------------------------------------------------------------------
 *	EPFL/STI/IOA/LIB/BM.4.137
 *	Philippe Thevenaz
 *	Station 17
 *	CH-1015 Lausanne VD
 *----------------------------------------------------------------------------
 *	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 below should be aligned)
 *	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|	|
 *  |...|...|...|...|...|...|...|...|...|...|...|...|...|...|...|...|...|...|
 ****************************************************************************/

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

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

/*****************************************************************************
 *	Definition of extern procedures
 ****************************************************************************/
/*--------------------------------------------------------------------------*/
#include "mex.h"
#include "matrix.h"

/********* CONVENTION REPERE IMAGE facon MATLAB ************/
/* --------------> Oy 
 * |
 * |
 * |
 * |
 * v
 *   Ox
 ***************************************/

void mexFunction(int nb_sorties,      mxArray *sorties[],
                 int nb_entrees, const mxArray *entrees[]) 

{ /* begin main */

	double	*ImageRasterArray, *CoeffMatrix ;
	double  *OutputImage;
	double	*p;
	double  *xOut, *yOut ; /* listes des abscisses et ordonnees ou on interpole */
/*	double	a11, a12, a21, a22;
	double	x0, y0, x1, y1;
	double	xOrigin, yOrigin;
	double	Angle, xShift, yShift; */
	long	Width, Height;
	long	Nx, Ny ;  /* taille de xOut et yOut */
	long	SplineDegree;
	long	x, y;
/*	int		Masking; */
	long	i,j ;
	int		Error;

/* lecture et verification des entrees */
	ImageRasterArray = mxGetPr(entrees[0]) ;
	Width = mxGetN(entrees[0]) ;
	Height = mxGetM(entrees[0]) ;		
	
	xOut = (double *)mxGetPr(entrees[1]) ;
	Nx = mxGetN(entrees[1]) ;
	yOut = (double *)mxGetPr(entrees[2]) ;
	Ny = mxGetN(entrees[2]) ;
	
	if (nb_entrees>3) SplineDegree = (long)*mxGetPr(entrees[3]);
	else SplineDegree=3;

/*	Masking  = 0 ; */ /*(int)*mxGetPr(entrees[4]);*/
/*	Angle = 0 ;
	xShift = 0 ;
	yShift = 0 ;
	xOrigin = 0 ;
	yOrigin = 0 ;
*/
	
	/* allocate output image */
	/* OutputImage = (float *)malloc((size_t)(Width * Height * (long)sizeof(float))); */
	sorties[0] =  mxCreateDoubleMatrix(Nx,Ny,mxREAL) ;
	OutputImage = mxGetPr(sorties[0]) ;
	
			/* mexPrintf("Allocation : Done\n"); */

	/* (float *)malloc((size_t)(Nx * Ny * (long)sizeof(float))); */
	if (OutputImage == (double *)NULL) {
		/* free(ImageRasterArray); */
		mexErrMsgTxt("Allocation of output image failed\n");
		/* return(1); */
	}

	/* convert between a representation based on image samples */
	/* and a representation based on image B-spline coefficients */
    CoeffMatrix = (double *)malloc((size_t)(Width * Height * (long)sizeof(double)));
	Error = SamplesToCoefficients(ImageRasterArray, CoeffMatrix, Width, Height, SplineDegree);
	if (Error) {
		/* mxDestroyArray(OutputImage); */
		mexErrMsgTxt("Change of basis failed\n");
		/* return(1); */
	}
	
		/* mexPrintf("Samples to coeff : Done\n"); */

	/* 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 */
	for (j=0 ; j<Ny ; j++)
	{
		for (i=0 ; i<Nx ; i++)
		{
			/* mexPrintf(" %d %d\n",i,j) ; */
            *OutputImage++ = (double)InterpolatedValue(CoeffMatrix, Width, Height, xOut[i], yOut[j], SplineDegree); 
            /* *OutputImage++ = CoeffMatrix[i+j*Height] ; */
		}
	}
		
	/* 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 (Masking) {
				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, SplineDegree);
				}
			}
			else {
				*p++ = (float)InterpolatedValue(ImageRasterArray, Width, Height,
					x1, y1, SplineDegree);
			}
		}
	}
	*/

	free(CoeffMatrix);
	/* mexPrintf("Done\n"); */
	/* return(0); */
} /* end main */
