Dear SPM-ers,

This is a question for the (C-)programmers in the user group. I'm having
strange effects when I try to carry out a trilinear interpolation.
I'm working on a program that transforms dynamic ECAT files, having
calculated the transformation with SPM routines. The transformation must
be the same for every frame, so the quickest way to do this is to build
a 3d mesh grid once, and interpolate all frames on that grid. 
The interp3.m interpolation function in Matlab is very slow, so I
decided to write a quick C/mexsol file to carry out the interpolation.
The interpolation now takes about 9 seconds per frame, instead of the
previous 45 seconds. But the images produced by the Matlab routine are
much smoother than those produced by my own routine. Another effect is a
high background signal.
Does anybody  know what I'm doing wrong? Are there any steps I've
forgotten (like for the boundaries of the volume)?
Included are the source file of my interpolation routine, a screen shot
of a transformed file interpolated by my own routine, and a screen shot
of a transformed file interpolated by Matlab's interp3 function.

Thanks a lot for any suggestions.

Best regards,

-- 

  Alle Meije Wink

  kamer Y0203	
  PET-centrum Academisch Ziekenhuis Groningen
  Hanzeplein 1
  9713 GZ Groningen

  tel: +31 50 361 34 39  
  email: [log in to unmask]

--



/*
	File		azgTLinterpolate.c

	Author		Alle Meije Wink

	Description
	  Matlab MEX function to interpolate 3D data on a mesh grid

	Calling syntax		
	  Vi = azgTLinterpolate(V,Xi,Yi,Zi)

	Input
	  V	3D data set
	  Xi	x co-ordinates of the mesh grid
	  Yi	y co-ordinates of the mesh grid
	  Zi	z co-ordinates of the mesh grid
	
	Output
	  Vi	the interpolated volume
	
	Remarks 
	  No dimension checks!
	
	See also
	  interp3

	PET Center, Groningen.
	University Hospital Groningen.
*/

static char Sccsid[] = "@(#)azgReadFrame.c	1.5   09/29/99";

void TLI_core();

#include <stdio.h>
#include <math.h>
#include "mex.h"
#include "/home/sources/Matlab5/toolbox/mexsol/cti_matrix.h"

#define offsetxyz(x,y,z,xysize,xsize) z*xysize+x*xsize+y

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

{ 
  double *Vi,*V,*Xi,*Yi,*Zi;
  const int  *dim_array;         

  /* Check for a valid number of arguments */
  if (nrhs != 4) 
     mexErrMsgTxt("azgReadFrame: 4 input arguments required");
  if (nlhs != 1) 
     mexErrMsgTxt("azgReadFrame: 1 output argument required");


  /* Get pointers to input arguments */
  V  = mxGetPr(prhs[0]);
  Xi = mxGetPr(prhs[1]);
  Yi = mxGetPr(prhs[2]);
  Zi = mxGetPr(prhs[3]);

  /* Get pointer to output argument. */
  dim_array=mxGetDimensions(prhs[0]);

  plhs[0] = mxCreateNumericArray(3,dim_array,mxDOUBLE_CLASS,mxREAL); 
  Vi = mxGetPr(plhs[0]);

  /* Carry out the interpolation */
  TLI_core(Vi,V,Xi,Yi,Zi,dim_array);
}


void TLI_core(ivol,vol,xco,yco,zco,dims)

/* 	Tri-Linear data interpolation		*/
/* 	Warning!	i <=> y			*/
/*			j <=> x			*/
/*						*/
/* 	(see first call of offfsetxyz)		*/

double		*ivol,*vol,*xco,*yco,*zco;		
int		*dims;
{
  int i,j,k;
  long offset;
  int xmin,xmax,ymin,ymax,zmin,zmax;
  double abovex,abovey,abovez,belowx,belowy,belowz;

  /* number of points in a plane */
  long plane=dims[0]*dims[1];

  /* x dimension of a plane */
  long xdim=dims[0];

  /* interpolate every point in the new volume */
  for (i=0;i<dims[0];i++)
      for (j=0;j<dims[1];j++)
	  for (k=0;k<dims[2];k++) 	  
	  {
	      /* calculate the offset for point (i,j,k) in all of the co-ordinate arrays */
	      offset=offsetxyz(j,i,k,plane,xdim);

	      /* check if x,y,z are inside the range */
	      if ( ((*(xco+offset)>0) && (*(xco+offset)<dims[0])) &&
	           ((*(yco+offset)>0) && (*(yco+offset)<dims[1])) &&
	           ((*(zco+offset)>0) && (*(zco+offset)<dims[2])) )

	      {
		 /* Use co-ordinates from the Matlab array */
		 /* Subtract 1: C programmers start at 0!  */

		 xmin=floor(*(xco+offset))-1;
		 xmax=xmin+1;
		 ymin=floor(*(yco+offset))-1;
		 ymax=ymin+1;
		 zmin=floor(*(zco+offset))-1;
		 zmax=zmin+1;

		 /* Compute the distance to the nearest integer co-ordinates */
		 /* Also compensate for the previous subtraction	     */

		 belowx=*(xco+offset)-xmin+1;
		 abovex=1-belowx;
		 belowy=*(yco+offset)-ymin+1;
		 abovey=1-belowy;
		 belowz=*(zco+offset)-zmin+1;
		 abovez=1-belowz;		 

		 *(ivol+offset) = belowx*belowy*belowz*(*(vol+offsetxyz(xmin,ymin,zmin,plane,xdim)))+
			       	  belowx*belowy*abovez*(*(vol+offsetxyz(xmin,ymin,zmax,plane,xdim)))+
			       	  belowx*abovey*belowz*(*(vol+offsetxyz(xmin,ymax,zmin,plane,xdim)))+
			       	  belowx*abovey*abovez*(*(vol+offsetxyz(xmin,ymax,zmax,plane,xdim)))+
			       	  abovex*belowy*belowz*(*(vol+offsetxyz(xmax,ymin,zmin,plane,xdim)))+
			       	  abovex*belowy*abovez*(*(vol+offsetxyz(xmax,ymin,zmax,plane,xdim)))+
			       	  abovex*abovey*belowz*(*(vol+offsetxyz(xmax,ymax,zmin,plane,xdim)))+
			       	  abovex*abovey*abovez*(*(vol+offsetxyz(xmax,ymax,zmax,plane,xdim)));
	      }

	  }   

}