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)));
}
}
}
|