diff --git a/xfibres.cc b/xfibres.cc
index 24539d939f866f562bcf0f8f932a63763bae6898..c68a14c583a54aae0c34565d55328452c951d67d 100644
--- a/xfibres.cc
+++ b/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);
}