Changes in grad-nonlin functions

CUDA/xfibres_gpu.cu

@@ -31,8 +31,6 @@ void xfibres_gpu(	//INPUT
 			const Matrix			bvecs,
 			const Matrix			bvals,
 			const Matrix	 		gradm, 
-			const Matrix 			Qform, 
-			const Matrix 			Qform_inv,
 			const NEWIMAGE::volume<int> 	vol2matrixkey,
 			const NEWMAT::Matrix		matrix2volkey,
 			const NEWIMAGE::volume<float>	mask,
@@ -64,7 +62,7 @@ void xfibres_gpu(	//INPUT
 		vector<Matrix> bvecs_vec, bvals_vec;
 
 		///// FIT /////
-		prepare_data_gpu_FIT(datam,bvecs,bvals,gradm,Qform,Qform_inv,datam_vec, bvecs_vec, bvals_vec, datam_host, bvecs_host,  bvals_host, alpha_host, beta_host, params_host, tau_host);	
+		prepare_data_gpu_FIT(datam,bvecs,bvals,gradm,datam_vec, bvecs_vec, bvals_vec, datam_host, bvecs_host,  bvals_host, alpha_host, beta_host, params_host, tau_host);	
 
 		thrust::device_vector<double> datam_gpu=datam_host;
 		thrust::device_vector<double> bvecs_gpu=bvecs_host;
@@ -125,7 +123,7 @@ void xfibres_gpu(	//INPUT
 
 // Correct bvals/bvecs accounting for Gradient Nonlinearities
 // ColumnVector grad_nonlin has 9 entries, corresponding to the 3 components of each of the x,y and z gradient deviation
-void correct_bvals_bvecs(const Matrix& bvals,const Matrix& bvecs, const ColumnVector& grad_nonlin, const Matrix& Qform, const Matrix& Qform_inv, Matrix& bvals_c, Matrix& bvecs_c){
+void correct_bvals_bvecs(const Matrix& bvals,const Matrix& bvecs, const ColumnVector& grad_nonlin, Matrix& bvals_c, Matrix& bvecs_c){
   	bvals_c=bvals; bvecs_c=bvecs;
   	Matrix L(3,3);  //gradient coil tensor
   	float mag;
@@ -137,22 +135,11 @@ void correct_bvals_bvecs(const Matrix& bvals,const Matrix& bvecs, const ColumnVe
   
   	for (int l=1; l<=bvals.Ncols(); l++){
     		if (bvals(1,l)>0){ //do not correct b0s
-      			//Rotate bvecs to scanner's coordinate system
-      			ColumnVector bvec_tmp(3);
-      			bvec_tmp=Qform*bvecs.Column(l);
-      			bvec_tmp(1)=-bvec_tmp(1); //Sign-flip X coordinate
-
-      			//Correct for grad-nonlin in scanner's coordinate system
-     		 	bvecs_c.Column(l)=(Id+L)*bvec_tmp;//bvecs.Column(l);
+     		 	bvecs_c.Column(l)=(Id+L)*bvecs.Column(l);
       			mag=sqrt(bvecs_c(1,l)*bvecs_c(1,l)+bvecs_c(2,l)*bvecs_c(2,l)+bvecs_c(3,l)*bvecs_c(3,l));
       			if (mag!=0)
 				bvecs_c.Column(l)=bvecs_c.Column(l)/mag;
-      			bvals_c(1,l)=mag*mag*bvals(1,l);
-      			bvec_tmp=bvecs_c.Column(l);
-
-      			//Rotate corrected bvecs back to voxel coordinate system
-      			bvec_tmp(1)=-bvec_tmp(1); //Sign-flip X coordinate
-      			bvecs_c.Column(l)=Qform_inv*bvec_tmp;
+      			bvals_c(1,l)=mag*mag*bvals(1,l);//mag^2 as b propto |G|^2
     		}
   	}
 }
@@ -335,8 +322,6 @@ void prepare_data_gpu_FIT(	//INPUT
 				const Matrix				bvecs,
 				const Matrix				bvals,
 				const Matrix	 			gradm, 
-				const Matrix 				Qform, 
-				const Matrix 				Qform_inv,
 				//OUTPUT
 				vector<ColumnVector>&			datam_vec,
 				vector<Matrix>&				bvecs_vec,
@@ -373,7 +358,7 @@ void prepare_data_gpu_FIT(	//INPUT
 
 	if (opts.grad_file.set()){
 		for(int vox=0;vox<nvox;vox++){
-			correct_bvals_bvecs(bvals,bvecs, gradm.Column(vox+1),Qform,Qform_inv,bvals_vec[vox],bvecs_vec[vox]); //correct for gradient nonlinearities
+			correct_bvals_bvecs(bvals,bvecs, gradm.Column(vox+1),bvals_vec[vox],bvecs_vec[vox]); //correct for gradient nonlinearities
  			MISCMATHS::cart2sph(bvecs_vec[vox],alpha,beta);
 			
 			for(int dir=0;dir<NDIRECTIONS;dir++){

CUDA/xfibres_gpu.cuh

@@ -35,8 +35,6 @@ void prepare_data_gpu_FIT(	//INPUT
 				const Matrix				bvecs,
 				const Matrix				bvals,
 				const Matrix	 			gradm, 
-				const Matrix 				Qform, 
-				const Matrix 				Qform_inv,
 				//OUTPUT
 				vector<ColumnVector>&			datam_vec,
 				vector<Matrix>&				bvecs_vec,

CUDA/xfibres_gpu.h

@@ -17,8 +17,6 @@ void xfibres_gpu(	//INPUT
 			const NEWMAT::Matrix			bvals,
 			const NEWMAT::Matrix			bvecs,
 			const NEWMAT::Matrix	 		gradm, 
-			const NEWMAT::Matrix 			Qform, 
-			const NEWMAT::Matrix 			Qform_inv,
 			const NEWIMAGE::volume<int> 		vol2matrixkey,
 			const NEWMAT::Matrix			matrix2volkey,
 			const NEWIMAGE::volume<float>		mask,

xfibres_gpu.cc

@@ -6,10 +6,6 @@
 
 /*  CCOPYRIGHT  */
 
-#ifndef EXPOSE_TREACHEROUS
-#define EXPOSE_TREACHEROUS
-#endif
-
 #include "CUDA/xfibres_gpu.h"
 #include "xfibresoptions.h"
 #include "newmat.h"
@@ -22,9 +18,6 @@
 
 using namespace Xfibres;
 
-
-void Return_Qform(const Matrix& qform_mat, Matrix& QMat, const float xdim, const float ydim, const float zdim);
-
 //////////////////////////////////////////////////////////
 //       XFIBRES CPU PART. IT CALLS TO GPU PART
 //////////////////////////////////////////////////////////
@@ -114,13 +107,10 @@ int xfibres_gpu(char *subjdir,int slice, int nextras, char** extras)
     	//Samples samples(vol2matrixkey,matrix2volkey,datam.Ncols(),datam.Nrows());	//in xfibres_gpu.cu
 
 	//Read Gradient Non_linearity Maps if provided
-    	NEWIMAGE::volume4D<float> grad; Matrix gradm, Qform, Qform_inv;
+    	NEWIMAGE::volume4D<float> grad; Matrix gradm;
     	if (opts.grad_file.set()){
       		read_volume4D(grad,opts.grad_file.value());
       		gradm=grad.matrix(mask);
-      		//Get the scale-free Qform matrix to rotate bvecs back to scanner's coordinate system 
-      		Return_Qform(data.qform_mat(), Qform, data.xdim(), data.ydim(), data.zdim());
-      		Qform_inv=Qform.i();
     	}
 	
 	myfile << "Number of Voxels: " << datam.Ncols() << endl;  
@@ -129,7 +119,7 @@ int xfibres_gpu(char *subjdir,int slice, int nextras, char** extras)
     	if(opts.rician.value() && !opts.nonlin.value()) 
       		cout<<"Rician noise model requested. Non-linear parameter initialization will be performed, overriding other initialization options!"<<endl;
 
-	xfibres_gpu(datam,bvecs,bvals,gradm,Qform,Qform_inv,vol2matrixkey,matrix2volkey,mask,slice,subjdir);
+	xfibres_gpu(datam,bvecs,bvals,gradm,vol2matrixkey,matrix2volkey,mask,slice,subjdir);
 
 	gettimeofday(&t2,NULL);
     	time=timeval_diff(&t2,&t1); 
@@ -139,15 +129,3 @@ int xfibres_gpu(char *subjdir,int slice, int nextras, char** extras)
   	return 0;
 }
 
-
-
-//Get the scale-free Qform matrix to rotate bvecs back to scanner's coordinate system 
-//After applying this matrix, make sure to sign flip the x coordinate of the resulted bvecs!
-//If you apply the inverse of the matrix, sign flip the x coordinate of the input bvecs
-void Return_Qform(const Matrix& qform_mat, Matrix& QMat, const float xdim, const float ydim, const float zdim){ 
-  	Matrix QMat_tmp; DiagonalMatrix Scale(3); 
-  	QMat_tmp=qform_mat.SubMatrix(1,3,1,3);
-  	Scale(1)=xdim; Scale(2)=ydim; Scale(3)=zdim;
-  	QMat_tmp=Scale.i()*QMat_tmp;
-  	QMat=QMat_tmp;
-}