/*
MEXRECON.C
        	Cmex function to do the Keller Reconstraction.

  The calling syntax is:

		R = mexrecon(Image, Tip, WaitbarFlag)

		- Image is the image matrix.
		- Tip is the matrix conteining the tip z-values.
		- WaitbarFlag 0 or 1.
		- R is the reconstracted image.
	
  Claudio   March 21, 1995

*/

#include <stdio.h>
#include <math.h>
#include "mex.h"
#include "matrix.h"

/* Input Arguments */

#define	M_IN	prhs[0]
#define M_TIP	prhs[1]
#define WF	prhs[2]


/* Output Arguments */

#define	M_OUT	plhs[0]



scan(double *pin, double *ptip, double *pout, int row, int col, int rowtip, int coltip, int flag)
{
int i, j, rt, ct;
Matrix *R, *Ti, *Tr, *D;
double *Rpr, *Tipr, *Trpr, *Dpr;
char str[48];

	R=mxCreateFull((row+rowtip), (col+coltip), REAL);
	Ti=mxCreateFull(rowtip, coltip, REAL);
	Tr=mxCreateFull(rowtip, coltip, REAL);
	D=mxCreateFull(rowtip, coltip, REAL);
	
	Rpr =(double *)mxGetPr(R);
	Tipr=(double *)mxGetPr(Ti);
	Trpr=(double *)mxGetPr(Tr);
	Dpr =(double *)mxGetPr(D);


	for(i=0; i<row; i++)
	    memcpy(Rpr+i*(col+coltip)+coltip/2+(col+coltip)*(rowtip/2), pin+i*col, sizeof(double)*col);


	for(i=0; i<row; i++)
	  {
	     if(flag == 1)
		{
	     	  sprintf(str, "progressbar(%f, 'Please wait');\0", (float)i/(row-1));
	     	  mexEvalString(str);
		}
	     for(j=0; j<col; j++)
		{
		  for(rt=0; rt<rowtip; rt++)
		    for(ct=0; ct<coltip; ct++)
			{
		 	  *(Tipr+rt*coltip+ct) = *(ptip+rt*coltip+ct) + *(pin+i*col+j);
			  *(Trpr+rt*coltip+ct) = *(Rpr+(i+rt)*(col+coltip)+j+ct);
			  *(Dpr+rt*coltip+ct)  = *(Tipr+rt*coltip+ct) - *(Trpr+rt*coltip+ct);
			  if(*(Dpr+rt*coltip+ct) < 0)
			     *(Dpr+rt*coltip+ct) = 0;

		          *(Rpr+(i+rt)*(col+coltip)+j+ct) = *(Tipr+rt*coltip+ct) - *(Dpr+rt*coltip+ct);
			}
		}
	  }

for(i=0; i<row; i++)
   memcpy(pout+i*col, Rpr+i*(col+coltip)+coltip/2+(col+coltip)*(rowtip/2), sizeof(double)*col);


mxDestroyArray(R);
mxDestroyArray(Ti);
mxDestroyArray(Tr);
mxDestroyArray(D);
return 0;
}


void mexFunction(int nlhs, Matrix *plhs[], int nrhs, Matrix *prhs[])
{
double	*pin;
double  *ptip;
double	*pout;
double  *wf;
unsigned int row,col, rowtip, coltip, wflag;


	/* Check for proper number of arguments */

	if (nrhs != 3)
		mexErrMsgTxt("MEXRECON requires two input arguments.");
	if (nlhs > 1)
		mexErrMsgTxt("MEXRECON requires one output argument.");
	


	/* Check the dimensions of M_IN */

	row = mxGetM(M_IN);
	col = mxGetN(M_IN);
	rowtip = mxGetM(M_TIP);
	coltip = mxGetN(M_TIP);


	/* Create a matrix for the return argument */

	M_OUT = mxCreateFull(row, col, REAL);


	/* Assign pointers to the various parameters */

	pout = (double *)mxGetPr(M_OUT);
	pin  = (double *)mxGetPr(M_IN);
	ptip = (double *)mxGetPr(M_TIP);
        wf   = (double *)mxGetPr(WF);

	wflag=(unsigned int)*wf;

	/* Do the actual computations in a subroutine */

	scan(pin, ptip, pout, row, col, rowtip, coltip, wflag);

return;
}