Added grad-nonlinearity routines

xfibres.cc

@@ -6,6 +6,10 @@
 
 /*  CCOPYRIGHT  */
 
+#ifndef EXPOSE_TREACHEROUS
+#define EXPOSE_TREACHEROUS
+#endif
+
 #include <iostream>
 #include <fstream>
 #include <iomanip>
@@ -760,6 +764,54 @@ public:
 };
 
 
+//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;
+}
+
+
+
+//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){
+  bvals_c=bvals; bvecs_c=bvecs;
+  Matrix L(3,3);  //gradient coil tensor
+  float mag;
+  L(1,1)=grad_nonlin(1);  L(1,2)=grad_nonlin(4);  L(1,3)=grad_nonlin(7);
+  L(2,1)=grad_nonlin(2);  L(2,2)=grad_nonlin(5);  L(2,3)=grad_nonlin(8);
+  L(3,1)=grad_nonlin(3);  L(3,2)=grad_nonlin(6);  L(3,3)=grad_nonlin(9);
+
+  IdentityMatrix Id(3);
+  
+  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);
+      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;
+    }
+  }
+}
+
 
 
 ////////////////////////////////////////////
@@ -790,30 +842,52 @@ int main(int argc, char *argv[])
 	bvecs(3,i)=bvecs(3,i)/tmpsum;
       }  
     }
-    
+
     volume4D<float> data;
     read_volume4D(data,opts.datafile.value());
     read_volume(mask,opts.maskfile.value());
     datam=data.matrix(mask); 
     matrix2volkey=data.matrix2volkey(mask);
     vol2matrixkey=data.vol2matrixkey(mask);
- 
-    Matrix Amat;
-    ColumnVector alpha, beta;
+    Samples samples(vol2matrixkey,matrix2volkey,datam.Ncols(),datam.Nrows());
+
+    //Read Gradient Non_linearity Maps if provided
+    volume4D<float> grad; Matrix gradm, Qform, Qform_inv;
+    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();
+    }
+
+    Matrix Amat; ColumnVector alpha, beta;
     Amat=form_Amat(bvecs,bvals);
     cart2sph(bvecs,alpha,beta);
-    Samples samples(vol2matrixkey,matrix2volkey,datam.Ncols(),datam.Nrows());
   
     if(opts.rician.value() && !opts.nonlin.value()) 
       cout<<"Rician noise model requested. Non-linear parameter initialization will be performed, overriding other initialization options!"<<endl;
-  
+
     for(int vox=1;vox<=datam.Ncols();vox++){
       cout <<vox<<"/"<<datam.Ncols()<<endl;
-      xfibresVoxelManager  vm(datam.Column(vox),alpha,beta,bvecs,bvals,samples,vox);
-      vm.initialise(Amat);
-      vm.runmcmc();
+      if (!opts.grad_file.set()){
+	xfibresVoxelManager  vm(datam.Column(vox),alpha,beta,bvecs,bvals,samples,vox);
+	vm.initialise(Amat);
+	vm.runmcmc();
+      }
+      else{ //Correct for each voxel the respective bvals/bvecs
+	Matrix Amat_c, bvals_c, bvecs_c;
+	ColumnVector alpha_c, beta_c;
+	
+	correct_bvals_bvecs(bvals,bvecs, gradm.Column(vox),Qform,Qform_inv,bvals_c,bvecs_c); //correct for gradient nonlinearities
+	Amat_c=form_Amat(bvecs_c,bvals_c);
+	cart2sph(bvecs_c,alpha_c,beta_c);
+	xfibresVoxelManager  vm(datam.Column(vox),alpha_c,beta_c,bvecs_c,bvals_c,samples,vox);
+	vm.initialise(Amat_c);
+	vm.runmcmc();
+      }
     }
-    
+
     samples.save(mask);
 
   }