interp2bili.c

本文件采用Matlab软件

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

/*




  Bilinear interpolation 


  Usage
  -----

  Z     = interp2bili(map , X , Y)

  Inputs
  -------

  map               map (nR x nC x m)
  X                 X coordinates (1 x N)
  Y                 Y coordinates (1 x N)

  Ouputs
  -------

  Z                 Interpolated values at (X , Y)



  Compile with:
  ------------


  mex -f mexopts_intel10amd.bat -output interp2bili.dll interp2bili.c


  Author          S閎astien PARIS ([e-mail]) 
  -------

  
*/


double my_nan(void) 
{
	double zero = 0.0;
	
	return zero/zero;
}

#define NAN my_nan()


/*-------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------*/
/*------------------------ Headers   --------------------------------------*/
/*-------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------*/


void interp2bili(double *, double * , double * , double * , int  , int , int , int);


/*-------------------------------------------------------------------------*/
/*-------------------------------------------------------------------------*/



void mexFunction( int nlhs, mxArray *plhs[] , int nrhs, const mxArray *prhs[] )

{
	
	
	double *Z , *xi , *yi , *Zi ;
	
	int *dimsZ=NULL , *dimsxi=NULL , *dimsyi=NULL , *dimsZi;
	
	int i , j , n , m , zz = 1 , zxi = 1 , zyi = 1 , numdimsZ;
	
	int numdimsxi , numdimsyi , numdimsZi;
	
	
	
	
	/*---------------------------------------------------------------*/
	/*---------------------- PARSE INPUT ----------------------------*/	
	/*---------------------------------------------------------------*/
	
	if(nrhs != 3)
		
	{
		
		mexErrMsgTxt("3 inputs are required for interp2bili");
		
	}
	
	/*----------- Input 1 ---------------*/
	
    
    Z        = mxGetPr(prhs[0]);
    
    numdimsZ = mxGetNumberOfDimensions(prhs[0]);
    
	dimsZ    = mxGetDimensions(prhs[0]);
	
    if (numdimsZ < 1)
		
	{
		mexErrMsgTxt("First input must at least a Matrix (n x m x s1,....,sn)");
	}
	
	for (i = 2 ; i < numdimsZ ; i++)
	{
		zz *= dimsZ[i];
	}
	
	n  = dimsZ[0];
	
	m  = dimsZ[1];
	
	
	/*----------- Input 2 ---------------*/
	
	xi        = mxGetPr(prhs[1]);
    
    numdimsxi = mxGetNumberOfDimensions(prhs[1]);
    
	dimsxi    = mxGetDimensions(prhs[1]);
	
	if (numdimsxi == 0)
		
	{
		
		mexErrMsgTxt("Second input must not be empty");
	}
	
    for (i = 0 ; i<numdimsxi ; i++)
	{
		zxi *= dimsxi[i];
	}
	
	
	
	/*----------- Input 3 ---------------*/
	
	
	yi        = mxGetPr(prhs[2]);
    
    numdimsyi = mxGetNumberOfDimensions(prhs[2]);
    
	dimsyi    = mxGetDimensions(prhs[2]);
	
	if (numdimsyi == 0)
		
	{  
		mexErrMsgTxt("Third input must not be empty");
	}
	
	for (i = 0 ; i<numdimsyi ; i++)
	{
		zyi *= dimsyi[i];
	}
	
    if (zxi != zyi)
		
	{
		mexErrMsgTxt("xi and yi must have the same size");
	}
	
	
	
	/*----------- Output 1 ---------------*/
	
	numdimsZi      = numdimsxi + (numdimsZ - 2);
	
	dimsZi         = (int *)mxMalloc(numdimsZi*sizeof(int));
	
	for (i = 0 ; i < numdimsxi ; i++)
	{
		dimsZi[i] = dimsxi[i];
	}	 
	
	j = 2;
	
	for (i = numdimsxi  ; i < numdimsZi ; i++)
	{
		
		dimsZi[i] = dimsZ[j];
		
		j        += 1;
		
	}
	
	
	plhs[0] = mxCreateNumericArray(numdimsZi, dimsZi, mxDOUBLE_CLASS, mxREAL);
	
	Zi      = mxGetPr(plhs[0]);
	
	
	/*----------- Computation of the quadratic form ---------------*/
	
	interp2bili(Z , xi , yi , Zi , n , m , zz , zxi);
	
	/*--------------------------------------------------------------*/
	
	mxFree(dimsZi);
	
	
}


/* ----------------------------------------------------------------------------------------------- */



void interp2bili(double *Z, double *xi , double *yi , double *Zi , int n  , int m , int zz  , int zxi)

{
	
	unsigned int i , v , nm = n*m , ndx , r , u;

	int floors , floort;
	
	double s , t , NotANumber = NAN;
	
	for (i = 0 ; i < zxi ; i++)
		
	{
		
		/* Out-Ranging Values */
		
		s = xi[i];
		
		t = yi[i];
		
		if (s < 1 || s > m || t < 1 || t > n)
			
		{
			for (v = 0 ; v < zz ; v++)
				
			{
				Zi[i + v*zxi] = NotANumber;				
			}
		}
		
		/* Interpolate values */
		
		
		else 
			
		{
			
          floors              = floor(s);

		  floort              = floor(t);

		  ndx                 = floort + (floors - 1)*n - 1;
			
		  s                   = (s - floors);
			
		  t                   = (t - floort);
			
			for (v = 0 ; v < zz ; v++)
				
			{
				r            = v*nm;

				u            = v*zxi;
				
				Zi[i + u]    = (Z[ndx + r]*(1 - t) + Z[ndx + r + 1]*t)*(1 - s) + ( Z[ndx + r + n]*(1 - t) + Z[ndx + n + 1 + r]*t )*s;
			}
			
			
		}
	}
	
}