/*
  mex file: pchirp2ia.c                     Steve Mann and Shawn Becker, 1993 
  resamples an image according to a p-vector of length 8.
  The p vector is calculated by one of:
        a) dechirp2.m (matlab function) or
	b) rechirp2.c (mex file)
  (this function is seldom called by itself)

  must specify interp order and antialias method
*/

#include <math.h>
#include <stdio.h>
#include "mex.h"
#include "dst_to_src_pchirp2.h"
#include "img_resample.h"

/* Input arguments are image, the vector of pchirp2 parameters, order, method */
#define P_IN prhs[0]
#define	IMAGE_IN prhs[1]
#define	INTERP_ORDER prhs[2]
#define	ANTIALIAS_METHOD prhs[3]
#define	NOTANUMBER prhs[4]
#define	M_OUT prhs[5]
#define	N_OUT prhs[6]
#define	OUTPUT_SCALING prhs[7]

/* output argument */
#define	OUTPUT	plhs[0]

#define	max(A, B)	((A) > (B) ? (A) : (B))
#define	min(A, B)	((A) < (B) ? (A) : (B))
#define pi 3.14159265
static int in_vdim, in_hdim;
static int out_vdim, out_hdim;
static int out_index, in_index;     /* one dim indices from rasterized.' 2d */
static double p1,p2,p3,p4,p5,p6,p7,p8;  /* pchirp parameters */


static
#ifdef __STDC__
void pchirp(
	double	image_out[],
	double	p_in[],
	double	image_in[],
	unsigned int M, unsigned int N,
	int interp_order, int antialias_method,
	double background_greylevel,
        unsigned int M_out, unsigned int N_out,
        unsigned int output_scaling
	)
#else
pchirp(image_out,
       p_in,
       image_in,
       M,N,
       interp_order,antialias_method,
       background_greylevel,
       M_out, N_out,
       output_scaling
      )
double	image_out[];
double	p_in[];
double	image_in[];
unsigned int M;
unsigned int N;
int interp_order;
int antialias_method;
double background_greylevel;
unsigned int M_out;
unsigned int N_out;
unsigned int output_scaling;
#endif
{
        float std_dev=1.0/20.0; /* Gaussian method now separate function*/
	in_vdim = (int) M;  /* convert unsigned int to int */
	in_hdim = (int) N;
	out_vdim = (int) M_out;
	out_hdim = (int) N_out;
	p1 = p_in[0];  /* set global variables that will be seen by */
	p2 = p_in[1];  /* the function img_resample_dst_to_src_once_per_pixel */
	p3 = p_in[2];
	p4 = p_in[3];
	p5 = p_in[4];
	p6 = p_in[5];
	p7 = p_in[6];
	p8 = p_in[7];
	/* print out the ``p-vector'' */
	/* printf("p=%g,%g, %g,%g, %g,%g, %g,%g\n",p1,p2,p3,p4,p5,p6,p7,p8); */
        dst_to_src_once_per_image(in_vdim,in_hdim,out_vdim,out_hdim,
				  p1,p2,p3,p4,p5,p6,p7,p8);
        img_resample(image_in,
	             in_vdim, in_vdim, in_hdim, 
                     interp_order, antialias_method,
                     image_out,
		     out_vdim, out_vdim, out_hdim,
		     background_greylevel,
                     output_scaling);

        return;         /* ?????? why is this return here??? not in book */
}


/* gateway routine */
#ifdef __STDC__
void mexFunction(
	int		nlhs,
	Matrix	*plhs[],
	int		nrhs,
	Matrix	*prhs[]
	)
#else
mexFunction(nlhs, plhs, nrhs, prhs)
int nlhs, nrhs;
Matrix *plhs[], *prhs[];
#endif
{
	double	*image_out;
	double	*p_in;
	double  *i_in;
	double  *a_in;
        int interp_order;
        int antialias_method;
	double  *notanumber_ptr;
	double   notanumber;
	double	*image_in;
	unsigned int	M, N;
        double  *M_out_ptr, *N_out_ptr;
	unsigned int	M_out,N_out;
        double  *output_scaling_ptr;
        unsigned int output_scaling;

	/* Check for proper number of arguments */
        
        /* fprintf(stderr,"pchirp2ia: in gateway routine\n"); */

	if (nrhs != 8) {
		mexErrMsgTxt("pchirp2ia requires 8 input arguments:");
		mexErrMsgTxt("       p, image, interp_order, antialias_method");
		mexErrMsgTxt("       notanumber (background grey level)");
		mexErrMsgTxt("       M_out, N_out (desired output size)");
		mexErrMsgTxt("       output_scaling (0 for off; default = 1)");
	} else if (nlhs > 1) {
		mexErrMsgTxt("pchirp2ia requires one output argument.");
	}

	/* Check the dimensions of P_IN;  can be 8 X 1 or 1 X 8. */
	M = mxGetM(P_IN);
	N = mxGetN(P_IN);
	if (!mxIsNumeric(P_IN) || mxIsComplex(P_IN) || 
		!mxIsFull(P_IN)  || !mxIsDouble(P_IN) ||
		(max(M,N) != 8) || (min(M,N) != 1)) {
	  	  mexErrMsgTxt("pchirp requires p to be an 8 element vector.");
	  }
	/* Check the dimensions of INTERP_ORDER; must be 1x1 */
	M = mxGetM(INTERP_ORDER);
	N = mxGetN(INTERP_ORDER);
	if (!mxIsNumeric(INTERP_ORDER) || mxIsComplex(INTERP_ORDER) ||
          !mxIsFull(INTERP_ORDER)  || !mxIsDouble(INTERP_ORDER) ||
	  (max(M,N) != 1) || (min(M,N) != 1)) {
    	    mexErrMsgTxt("pchirp requires that interp_order be real scalar.");
	  }
	/* Check the dimensions of ANTIALIAS_METHOD; must be 1x1 */
	M = mxGetM(ANTIALIAS_METHOD);
	N = mxGetN(ANTIALIAS_METHOD);
	if (!mxIsNumeric(ANTIALIAS_METHOD) || mxIsComplex(ANTIALIAS_METHOD) ||
          !mxIsFull(ANTIALIAS_METHOD)  || !mxIsDouble(ANTIALIAS_METHOD) ||
	  (max(M,N) != 1) || (min(M,N) != 1)) {
    	    mexErrMsgTxt("pchirp requires that antialias_method be scalar.");
	  }
        /* Check the dimensions of NOTANUMBER; must be 1x1 */
        M = mxGetM(NOTANUMBER);
        N = mxGetN(NOTANUMBER);
        if (!mxIsNumeric(NOTANUMBER) || mxIsComplex(NOTANUMBER) ||
          !mxIsFull(NOTANUMBER)  || !mxIsDouble(NOTANUMBER) ||
          (max(M,N) != 1) || (min(M,N) != 1)) {
            mexErrMsgTxt("pchirp requires mynan be real scalar.");
	  }
        /* Check the dimensions of M_OUT; must be 1x1 */
        M = mxGetM(M_OUT);
        N = mxGetN(M_OUT);
        if (!mxIsNumeric(M_OUT) || mxIsComplex(M_OUT) ||
          !mxIsFull(M_OUT)  || !mxIsDouble(M_OUT) ||
          (max(M,N) != 1) || (min(M,N) != 1)) {
            mexErrMsgTxt("pchirp requires that M_out be real scalar.");
	  }
        /* Check the dimensions of N_OUT; must be 1x1 */
        M = mxGetM(N_OUT);
        N = mxGetN(N_OUT);
        if (!mxIsNumeric(N_OUT) || mxIsComplex(N_OUT) ||
          !mxIsFull(N_OUT)  || !mxIsDouble(N_OUT) ||
          (max(M,N) != 1) || (min(M,N) != 1)) {
            mexErrMsgTxt("pchirp requires that N_out be real scalar.");
	  }
        /* Check the dimensions of OUTPUT_SCALING; must be 1x1 */
        M = mxGetM(OUTPUT_SCALING);
        N = mxGetN(OUTPUT_SCALING);
        if (!mxIsNumeric(OUTPUT_SCALING) || mxIsComplex(OUTPUT_SCALING) ||
          !mxIsFull(OUTPUT_SCALING)  || !mxIsDouble(OUTPUT_SCALING) ||
          (max(M,N) != 1) || (min(M,N) != 1)) {
            mexErrMsgTxt("pchirp requires that output_scaling be real scalar.");
	  }
	/* get input image size; reuse variables M and N */
	M = mxGetM(IMAGE_IN);
	N = mxGetN(IMAGE_IN);

	/* OUTPUT = mxCreateFull(M, N, REAL); */

	/* Assign pointers to the various parameters */
	/* image_out = mxGetPr(OUTPUT); */
	image_in = mxGetPr(IMAGE_IN);
	p_in = mxGetPr(P_IN);
	i_in = mxGetPr(INTERP_ORDER);
	a_in = mxGetPr(ANTIALIAS_METHOD);
	interp_order = i_in[0];
	antialias_method = a_in[0];
	notanumber_ptr = mxGetPr(NOTANUMBER);
        notanumber = notanumber_ptr[0];
        M_out_ptr = mxGetPr(M_OUT);
        N_out_ptr = mxGetPr(N_OUT);
        M_out = M_out_ptr[0]; /* appears to cast automatically from double to unsigned int */
        N_out = (unsigned int) N_out_ptr[0]; /* try casting explicitly */
        fprintf(stderr,"pchirp2ia: M_out=%d, N_out=%d\n",M_out,N_out);
        output_scaling_ptr = mxGetPr(OUTPUT_SCALING);
        output_scaling = (unsigned int) output_scaling_ptr[0];
        fprintf(stderr,"pchirp2ia: output_scaling = %d\n",output_scaling);
        fprintf(stderr,"pchirp2ia: about to mxCreateFull M_out by N_out size\n");

        /* Create a matrix for the return image */
        OUTPUT = mxCreateFull(M_out, N_out, REAL);
	image_out = mxGetPr(OUTPUT);

	/* Do the actual computations in a subroutine */

	pchirp(image_out,
	       p_in,
	       image_in,
	       M,N,
	       (int) interp_order, (int) antialias_method,
	       notanumber,           /* BACKGROUND for undefined values */
               M_out, N_out,
               output_scaling
	      );
	return;
}