From c57df158bc8b40b2a18678ebac6300ea791d96a8 Mon Sep 17 00:00:00 2001
From: Stamatios Sotiropoulos <stam@fmrib.ox.ac.uk>
Date: Fri, 21 Jun 2013 12:09:35 +0000
Subject: [PATCH] Added gradient nonlinearity correction options
---
rubix.cc | 93 +++++++++++++++++++++++++++++++++++++-------------
rubix.h | 8 ++---
rubixoptions.h | 11 ++++++
rubixvox.cc | 8 ++---
rubixvox.h | 12 +++----
5 files changed, 95 insertions(+), 37 deletions(-)
diff --git a/rubix.cc b/rubix.cc
index 33b92f8..14d22a7 100644
--- a/rubix.cc
+++ b/rubix.cc
@@ -394,32 +394,22 @@ void LRVoxelManager::initialise(){
for (int n=0; n<m_HRvoxnumber.Nrows(); n++){
if (opts.modelnum.value()==1){ //Model 1
- PVM_single_c pvm(m_dataHR[n],m_bvecsHR,m_bvalsHR,opts.nfibres.value());
+ PVM_single_c pvm(m_dataHR[n],m_bvecsHR[n],m_bvalsHR[n],opts.nfibres.value());
pvm.fit(); // this will give th,ph,f in the correct order
pvmf = pvm.get_f(); pvmth = pvm.get_th(); pvmph = pvm.get_ph();
pvmS0 = fabs(pvm.get_s0()); pvmd = pvm.get_d(); predicted_signal=pvm.get_prediction();
if(pvmd<0 || pvmd>0.01) pvmd=2e-3;
-
- // DTI dti1(m_dataHR[n],m_bvecsHR,m_bvalsHR);
- //dti1.linfit();
- //pvmS0=fabs(dti1.get_s0());
-
m_LRv.set_HRparams(n,pvmd,pvmS0,pvmth,pvmph,pvmf);
}
else{ //Model 2
- PVM_multi pvm(m_dataHR[n],m_bvecsHR,m_bvalsHR,opts.nfibres.value());
+ PVM_multi pvm(m_dataHR[n],m_bvecsHR[n],m_bvalsHR[n],opts.nfibres.value());
pvm.fit();
pvmf = pvm.get_f(); pvmth = pvm.get_th(); pvmph = pvm.get_ph(); pvmd_std=pvm.get_d_std();
pvmS0 =fabs(pvm.get_s0()); pvmd = pvm.get_d(); predicted_signal=pvm.get_prediction();
if(pvmd<0 || pvmd>0.01) pvmd=2e-3;
if(pvmd_std<0 || pvmd_std>0.01) pvmd_std=pvmd/10;
-
- // DTI dti1(m_dataHR[n],m_bvecsHR,m_bvalsHR);
- //dti1.linfit();
- //pvmS0=fabs(dti1.get_s0());
-
m_LRv.set_HRparams(n,pvmd,pvmd_std,pvmS0,pvmth,pvmph,pvmf);
}
@@ -558,6 +548,30 @@ volume<float> createHR_mask(const volume<float>& maskLR, const volume4D<float>&
}
+//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, 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);
+
+ //Correct each gradient
+ for (int l=1; l<=bvals.Ncols(); l++){
+ if (bvals(1,l)>0){ //do not correct b0s
+ 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); //mag^2 as b propto |G|^2
+ }
+ }
+}
+
////////////////////////////////////////////
// MAIN
@@ -624,32 +638,65 @@ int main(int argc, char *argv[])
HRSamples HRsampl(datamHR.Ncols(), opts.njumps.value(), opts.sampleevery.value(), opts.nfibres.value(), opts.rician.value(), opts.modelnum.value());
LRSamples LRsampl(datamLR.Ncols(), opts.njumps.value(), opts.sampleevery.value(), opts.nmodes.value(), opts.rician.value(),opts.fsumPrior.value(),opts.dPrior.value());
+ //Read Gradient Non_linearity Maps if provided
+ volume4D<float> LRgrad, HRgrad; Matrix LRgradm, HRgradm;
+ if (opts.LRgrad_file.set()){
+ read_volume4D(LRgrad,opts.LRgrad_file.value());
+ LRgradm=LRgrad.matrix(maskLR);
+ if (opts.HRgrad_file.set()){
+ read_volume4D(HRgrad,opts.HRgrad_file.value());
+ HRgradm=HRgrad.matrix(maskHR);
+ }
+ else{
+ cerr<<"HR grad_dev file should be provided as well!"<<endl;
+ return 1;
+ }
+ }
+ else{
+ if (opts.HRgrad_file.set()){
+ cerr<<"LR grad_dev file should be provided as well!"<<endl;
+ return 1;
+ }
+ }
+
//dHR.push_back(datamHR.Column(1)); dHR.push_back(datamHR.Column(2));
//dHR.push_back(datamHR.Column(3)); dHR.push_back(datamHR.Column(4));
//HRvoxnum<<1<<2<<3<<4; HRweights<<0.25<<0.25<<0.25<<0.25;
for(int vox=1;vox<=datamLR.Ncols();vox++){ //For each LR voxel
- // if (vox==19){ //vox==6
ColumnVector dLR; Matrix HRindices;
- dLR=datamLR.Column(vox); //Find the corresponding HR ones
+ dLR=datamLR.Column(vox);
+ if (opts.LRgrad_file.set()){ //Correct the respective bvals/bvecs
+ Matrix bvals_c, bvecs_c;
+ correct_bvals_bvecs(bvalsLR,bvecsLR, LRgradm.Column(vox),bvals_c,bvecs_c); //correct for gradient nonlinearities
+ bvalsLR=bvals_c; bvecsLR=bvecs_c;
+ }
+
+ //Find the corresponding HR ones
HRindices=get_HRindices((int)matrix2volkeyLR(vox,1),(int)matrix2volkeyLR(vox,2),(int)matrix2volkeyLR(vox,3),(int)xratio,(int)yratio,(int)zratio);
- //cout<<endl<<"S0LR: "<<dLR(1)<<endl<<endl; //Debugging code
-
ColumnVector HRvoxnum(HRindices.Nrows()), HRweights(HRindices.Nrows());
- vector<ColumnVector> dHR;
+ vector<ColumnVector> dHR; vector<Matrix> vbvalsHR; vector<Matrix> vbvecsHR;
- for (int n=1; n<=HRindices.Nrows(); n++){
+ for (int n=1; n<=HRindices.Nrows(); n++){ //For each HR voxel
HRweights(n)=1.0/HRindices.Nrows();
HRvoxnum(n)=vol2matrixkeyHR((int)HRindices(n,1),(int)HRindices(n,2),(int)HRindices(n,3));
ColumnVector tmp_vec=datamHR.Column((int)HRvoxnum(n));
- //cout<<(int)HRindices(n,1)<<" "<<(int)HRindices(n,2)<<" "<<(int)HRindices(n,3)<<" S0HR:"<<tmp_vec(1)<<endl;
- dHR.push_back(datamHR.Column((int)HRvoxnum(n)));
- }
+ dHR.push_back(datamHR.Column((int)HRvoxnum(n)));
+ if (!opts.HRgrad_file.set()){
+ vbvalsHR.push_back(bvalsHR);
+ vbvecsHR.push_back(bvecsHR);
+ }
+ else{ //Correct for each voxel the respective bvals/bvecs
+ Matrix bvals_c, bvecs_c;
+ correct_bvals_bvecs(bvalsHR,bvecsHR, HRgradm.Column((int)HRvoxnum(n)),bvals_c,bvecs_c); //correct for gradient nonlinearities
+ vbvalsHR.push_back(bvals_c);
+ vbvecsHR.push_back(bvecs_c);
+ }
+ }
cout <<vox<<"/"<<datamLR.Ncols()<<endl;
- LRVoxelManager vm(HRsampl,LRsampl,vox,HRvoxnum, dLR,dHR, bvecsLR, bvalsLR, bvecsHR, bvalsHR, HRweights);
+ LRVoxelManager vm(HRsampl,LRsampl,vox,HRvoxnum, dLR,dHR, bvecsLR, bvalsLR, vbvecsHR, vbvalsHR, HRweights);
vm.initialise();
vm.runmcmc();
- // }
}
HRsampl.save(maskHR);
LRsampl.save(maskLR);
diff --git a/rubix.h b/rubix.h
index 28f378f..2c6c115 100644
--- a/rubix.h
+++ b/rubix.h
@@ -197,17 +197,17 @@ namespace RUBIX{
ColumnVector m_HRvoxnumber;
LRvoxel m_LRv;
const ColumnVector& m_dataLR; //Low-Res Data for the specific LR voxel
- const vector<ColumnVector>& m_dataHR; //High-Res Data for all contained HR voxels
+ const vector<ColumnVector>& m_dataHR; //High-Res Data for all HRvoxels within a LRvoxel
const Matrix& m_bvecsLR; //bvecs at Low-Res (3 x LR_NumPoints)
const Matrix& m_bvalsLR; //bvalues at Low-Res (1 x HR_NumPoints)
- const Matrix& m_bvecsHR; //bvecs at High-Res (3 x HR_NumPoints)
- const Matrix& m_bvalsHR; //bvalues at High-Res (1 x HR_NumPoints)
+ const vector<Matrix>& m_bvecsHR; //bvecs at High-Res (HRvoxels within a LRvoxel x 3 x HR_NumPoints)
+ const vector<Matrix>& m_bvalsHR; //bvalues at High-Res (HRvoxels within a LRvoxel x 1 x HR_NumPoints)
const ColumnVector& m_HRweights; //Holds the volume fraction each HR voxel occupies out of the LR one
public:
//Constructor
LRVoxelManager(HRSamples& Hsamples, LRSamples& Lsamples, int LRvoxnum, ColumnVector& HRvoxnum,
const ColumnVector& dataLR,const vector<ColumnVector>& dataHR,
- const Matrix& bvecsLR, const Matrix& bvalsLR, const Matrix& bvecsHR, const Matrix& bvalsHR, const ColumnVector& HRweights):
+ const Matrix& bvecsLR, const Matrix& bvalsLR, const vector<Matrix>& bvecsHR, const vector<Matrix>& bvalsHR, const ColumnVector& HRweights):
opts(rubixOptions::getInstance()), m_HRsamples(Hsamples), m_LRsamples(Lsamples), m_LRvoxnumber(LRvoxnum),m_HRvoxnumber(HRvoxnum),
m_LRv(bvecsHR, bvalsHR, bvecsLR, bvalsLR, dataLR, dataHR, opts.nfibres.value(), opts.nmodes.value(), HRweights, opts.modelnum.value(), opts.fudge.value(),opts.all_ard.value(), opts.no_ard.value(),opts.kappa_ard.value(), opts.fsumPrior.value(), opts.dPrior.value(), opts.rician.value()),
m_dataLR(dataLR), m_dataHR(dataHR),m_bvecsLR(bvecsLR), m_bvalsLR(bvalsLR), m_bvecsHR(bvecsHR), m_bvalsHR(bvalsHR), m_HRweights(HRweights) { }
diff --git a/rubixoptions.h b/rubixoptions.h
index cb8654c..77fe0d5 100644
--- a/rubixoptions.h
+++ b/rubixoptions.h
@@ -54,6 +54,9 @@ namespace RUBIX{
Option<bool> fsumPrior;
Option<bool> dPrior;
Option<bool> rician;
+ FmribOption<string> LRgrad_file;
+ FmribOption<string> HRgrad_file;
+
void parse_command_line(int argc, char** argv, Log& logger);
private:
@@ -145,6 +148,12 @@ namespace RUBIX{
false,no_argument),
rician(string("--rician"),false,string("Use Rician Noise model (default is Gaussian)"),
false,no_argument),
+ LRgrad_file(string("--gLR, --gradnonlinLR"), string("grad_devLR"),
+ string("LR Gradient Nonlinearity Tensor"),
+ false, requires_argument),
+ HRgrad_file(string("--gHR, --gradnonlinHR"), string("grad_devHR"),
+ string("HR Gradient Nonlinearity Tensor"),
+ false, requires_argument),
options("RubiX v1.0", "rubix --help (for list of options)\n")
{
try {
@@ -175,6 +184,8 @@ namespace RUBIX{
options.add(fsumPrior);
options.add(dPrior);
options.add(rician);
+ options.add(LRgrad_file);
+ options.add(HRgrad_file);
}
catch(X_OptionError& e) {
options.usage();
diff --git a/rubixvox.cc b/rubixvox.cc
index fb07762..804ca3b 100644
--- a/rubixvox.cc
+++ b/rubixvox.cc
@@ -1030,7 +1030,7 @@ void LRvoxel::compute_prior(){
void LRvoxel::compute_likelihood(){
ColumnVector SLRpred, SHRpred, Diff;
SLRpred.ReSize(m_bvecsLR.Ncols());
- SHRpred.ReSize(m_bvecsHR.Ncols()); SHRpred=0;
+ SHRpred.ReSize(m_bvecsHR[0].Ncols()); SHRpred=0;
float likLR, likHR;
m_old_likelihood_en=m_likelihood_en;
@@ -1051,7 +1051,7 @@ void LRvoxel::compute_likelihood(){
SHRpred=m_dataHR[m]-SHRpred;
likHR+=log(0.5*SHRpred.SumSquare());
}
- likHR*=0.5*m_bvecsHR.Ncols();
+ likHR*=0.5*m_bvecsHR[0].Ncols();
m_likelihood_en=likLR+likHR;
}
@@ -1070,8 +1070,8 @@ void LRvoxel::compute_likelihood(){
for (unsigned int m=0; m<m_dataHR.size(); m++){
float m_tauHR=m_HRvoxels[m].get_tau();
SHRpred=m_HRvoxels[m].getSignalHR();
- likHR+=m_bvecsHR.Ncols()*log(m_tauHR);
- for (int k=1; k<=m_bvecsHR.Ncols(); k++)
+ likHR+=m_bvecsHR[0].Ncols()*log(m_tauHR);
+ for (int k=1; k<=m_bvecsHR[0].Ncols(); k++)
likHR+=m_logdataHR[m](k)-0.5*m_tauHR*(m_dataHR[m](k)*m_dataHR[m](k)+SHRpred(k)*SHRpred(k))+logIo(m_tauHR*m_dataHR[m](k)*SHRpred(k));
}
m_likelihood_en=likLR-likHR;
diff --git a/rubixvox.h b/rubixvox.h
index 572104a..9cbd5af 100644
--- a/rubixvox.h
+++ b/rubixvox.h
@@ -545,9 +545,9 @@ namespace RUBIX{
vector<ColumnVector> m_logdataHR; //Log of High-Res Data for all contained HR voxels (use it in Rician Energy)
const ColumnVector& m_dataLR; //Low-Res Data for the specific LR voxel
- const vector<ColumnVector>& m_dataHR; //High-Res Data for all contained HR voxels
- const Matrix& m_bvecsHR; //bvecs at High-Res (3 x HR_NumPoints)
- const Matrix& m_bvalsHR; //bvalues at High-Res (1 x HR_NumPoints)
+ const vector<ColumnVector>& m_dataHR; //High-Res Data for all HRvoxels within a LRvoxel
+ const vector<Matrix>& m_bvecsHR; //bvecs at High-Res (HRvoxels within a LRvoxel x 3 x HR_NumPoints)
+ const vector<Matrix>& m_bvalsHR; //bvalues at High-Res (HRvoxels within a LRvoxel x 3 x HR_NumPoints)
const Matrix& m_bvecsLR; //bvecs at Low-Res (3 x LR_NumPoints)
const Matrix& m_bvalsLR; //bvalues at Low-Res (1 x HR_NumPoints)
const int m_modelnum; //1 for single-shell, 2 for multi-shell model
@@ -564,7 +564,7 @@ namespace RUBIX{
public:
//Constructor
- LRvoxel(const Matrix& bvecsHR, const Matrix& bHR,
+ LRvoxel(const vector<Matrix>& bvecsHR, const vector<Matrix>& bHR,
const Matrix& bvecsLR, const Matrix& bLR,
const ColumnVector& dataLR, const vector<ColumnVector>& dataHR, const int N, const int Nmodes, const ColumnVector& HRweights, const int modelnum=1, const float ardfudge=1, const bool allard=false, const bool Noard=false, const bool kappa_ard=false, const bool fsumPrior_ON=false, const bool dPrior_ON=false, const bool rician=false):
m_dataLR(dataLR), m_dataHR(dataHR),m_bvecsHR(bvecsHR), m_bvalsHR(bHR), m_bvecsLR(bvecsLR), m_bvalsLR(bLR),
@@ -589,8 +589,8 @@ namespace RUBIX{
}
m_Orient_hyp_prior.ReSize(m_Nmodes,5);
- for (unsigned int m=1; m<=m_dataHR.size(); m++){ //Add HRvoxel Objects
- HRvoxel HRv(m_bvecsHR, m_bvalsHR, m_bvecsLR, m_bvalsLR, m_Orient_hyp_prior, m_mean_d, m_stdev_d, m_mean_fsum, m_stdev_fsum, m_numfibres, m_modelnum, m_ardfudge, m_fsumPrior_ON, m_dPrior_ON, m_rician);
+ for (unsigned int n=0; n<m_dataHR.size(); n++){ //Add HRvoxel Objects
+ HRvoxel HRv(m_bvecsHR[n], m_bvalsHR[n], m_bvecsLR, m_bvalsLR, m_Orient_hyp_prior, m_mean_d, m_stdev_d, m_mean_fsum, m_stdev_fsum, m_numfibres, m_modelnum, m_ardfudge, m_fsumPrior_ON, m_dPrior_ON, m_rician);
m_HRvoxels.push_back(HRv);
}
--
GitLab