医学图像配准

#include "mex.h" #include "math.h" /* 3D Bspline transformation grid function * function Vout=bspline_transform_3d_single(Ox,Oy,Oz,Vin,nx,ny,nz) * * Ox, Oy, Oz are the grid points coordinates * Vin is input image, Vout the transformed output image * nx,ny and nz the dimensions of the grid (inside the image) * * This function is an implementation of the b-spline registration * algorithm in "D. Rueckert et al. : Nonrigid Registration Using Free-Form * Deformations: Application to Breast MR Images". * * We used "Fumihiko Ino et al. : a data distrubted parallel algortihm for * nonrigid image registration" for the correct formula's, because * (most) other papers contain errors. * * Function is written by D.Kroon University of Twente (July 2008) */ // Convert 2D/3D matrix index to 1D index int mindex2(int x, int y, int sizx) { return y*sizx+x; } int mindex3(int x, int y, int z, int sizx, int sizy) { return z*sizx*sizy+y*sizx+x;} // The matlab mex function void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) { // Ox and Oy are the grid points // Zo is the input image // Zi is the transformed image // nx and ny are the number of grid points (inside the image) float *Ox,*Oy,*Oz,*ZO, *nx, *ny,*nz,*ZI; // Size of input image mwSize ZOsizex, ZOsizey, ZOsizez; const mwSize *dims; // Size of grid mwSize Osizex, Osizey, Osizez; // B-spline variables double u,v,w; int i,j,k; double Bu[4]={0,0,0,0}; double Bv[4]={0,0,0,0}; double Bw[4]={0,0,0,0}; // B-Spline loop variabels int l,m,n; // loop variable int p; // Grid distance; double dx,dy,dz; // Location of pixel which will be come the current pixel float Tlocalx; float Tlocaly; float Tlocalz; // Variables to store 1D index int indexO; int indexZI; int indexZO; // X,Y,Z coordinates of current pixel int x,y,z; // Linear interpolation variables int xBas[8], yBas[8], zBas[8]; float perc[8]; float xCom, yCom, zCom; float color[8]={0,0,0,0,0,0,0,0}; /* Check for proper number of arguments. */ if(nrhs!=7) { mexErrMsgTxt("Seven inputs are required."); } else if(nlhs!=1) { mexErrMsgTxt("One output required"); } // nsubs=mxGetNumberOfDimensions(prhs[0]); // Get the sizes of the grid dims = mxGetDimensions(prhs[0]); Osizex = dims[0]; Osizey = dims[1]; Osizez = dims[2]; // Create image matrix for the return arguments with the size of input image dims = mxGetDimensions(prhs[3]); ZOsizex = dims[0]; ZOsizey = dims[1]; ZOsizez = dims[2]; plhs[0] = mxCreateNumericArray(3, dims, mxSINGLE_CLASS, mxREAL); /* Assign pointers to each input. */ Ox=(float *)mxGetData(prhs[0]); Oy=(float *)mxGetData(prhs[1]); Oz=(float *)mxGetData(prhs[2]); ZO=(float *)mxGetData(prhs[3]); nx=(float *)mxGetData(prhs[4]); ny=(float *)mxGetData(prhs[5]); nz=(float *)mxGetData(prhs[6]); /* Assign pointer to output. */ ZI = (float *)mxGetData(plhs[0]); // Calculate grid spacing dx=(float)(ZOsizex-1)/(nx[0]-1); dy=(float)(ZOsizey-1)/(ny[0]-1); dz=(float)(ZOsizez-1)/(nz[0]-1); // Loop through all image pixel coordinates for (x=0; x<ZOsizex; x++) { for (y=0; y<ZOsizey; y++) { for (z=0; z<ZOsizez; z++) { // Calculate B-spline values u=(x/dx)-floor(x/dx); Bu[0] = pow((1-u),3)/6; Bu[1] = ( 3*pow(u,3) - 6*pow(u,2) + 4)/6; Bu[2] = (-3*pow(u,3) + 3*pow(u,2) + 3*u + 1)/6; Bu[3] = pow(u,3)/6; // Calculate B-spline values v=y/dy-floor(y/dy); Bv[0] = pow((1-v),3)/6; Bv[1] = ( 3*pow(v,3) - 6*pow(v,2) + 4)/6; Bv[2] = (-3*pow(v,3) + 3*pow(v,2) + 3*v + 1)/6; Bv[3] = pow(v,3)/6; // Calculate B-spline values w=z/dz-floor(z/dz); Bw[0] = pow((1-w),3)/6; Bw[1] = ( 3*pow(w,3) - 6*pow(w,2) + 4)/6; Bw[2] = (-3*pow(w,3) + 3*pow(w,2) + 3*w + 1)/6; Bw[3] = pow(w,3)/6; i=(int)floor(x/dx); //-1 first row against boundary artefacts j=(int)floor(y/dy); k=(int)floor(z/dz); // This part calculates the coordinates of the pixel // which will be transformed to the current x,y pixel. Tlocalx=0; Tlocaly=0; Tlocalz=0; for(l=0; l<4; l++) { for(m=0; m<4; m++) { for(n=0; n<4; n++) { if(((i+l)>=0)&&((i+l)<Osizex)&&((j+m)>=0)&&((j+m)<Osizey)&&((k+n)>=0)&&((k+n)<Osizez)) { indexO=mindex3(i+l,j+m,k+n,Osizex,Osizey); Tlocalx=Tlocalx+(float)(Bu[l]*Bv[m]*Bw[n])*Ox[indexO]; Tlocaly=Tlocaly+(float)(Bu[l]*Bv[m]*Bw[n])*Oy[indexO]; Tlocalz=Tlocalz+(float)(Bu[l]*Bv[m]*Bw[n])*Oz[indexO]; } } } } // Determine the coordinates of the pixel(s) which will be come the current pixel // (using linear interpolation) xBas[0]=(int) floor(Tlocalx); yBas[0]=(int) floor(Tlocaly); zBas[0]=(int) floor(Tlocalz); xBas[1]=xBas[0]+0; yBas[1]=yBas[0]+0; zBas[1]=zBas[0]+1; xBas[2]=xBas[0]+0; yBas[2]=yBas[0]+1; zBas[2]=zBas[0]+0; xBas[3]=xBas[0]+0; yBas[3]=yBas[0]+1; zBas[3]=zBas[0]+1; xBas[4]=xBas[0]+1; yBas[4]=yBas[0]+0; zBas[4]=zBas[0]+0; xBas[5]=xBas[0]+1; yBas[5]=yBas[0]+0; zBas[5]=zBas[0]+1; xBas[6]=xBas[0]+1; yBas[6]=yBas[0]+1; zBas[6]=zBas[0]+0; xBas[7]=xBas[0]+1; yBas[7]=yBas[0]+1; zBas[7]=zBas[0]+1; for (p=0; p<8; p++) { if(xBas[p]>=0&&xBas[p]<ZOsizex&&yBas[p]>=0&&yBas[p]<ZOsizey&&zBas[p]>=0&&zBas[p]<ZOsizez) { indexZI=mindex3(xBas[p],yBas[p],zBas[p],ZOsizex,ZOsizey); color[p]=ZO[indexZI]; } else { color[p]=0; } } // Linear interpolation constants (percentages) xCom=Tlocalx-(float)floor((double)Tlocalx); yCom=Tlocaly-(float)floor((double)Tlocaly); zCom=Tlocalz-(float)floor((double)Tlocalz); perc[0]=(1-xCom) * (1-yCom) * (1-zCom); perc[1]=(1-xCom) * (1-yCom) * zCom; perc[2]=(1-xCom) * yCom * (1-zCom); perc[3]=(1-xCom) * yCom * zCom; perc[4]=xCom * (1-yCom) * (1-zCom); perc[5]=xCom * (1-yCom) * zCom; perc[6]=xCom * yCom * (1-zCom); perc[7]=xCom * yCom * zCom; // Set the current pixel value indexZO=mindex3(x,y,z,ZOsizex,ZOsizey); ZI[indexZO]=0; for (p=0; p<8; p++) { ZI[indexZO]+=color[p]*perc[p]; } } } } }