From a33476e37f85ef52d5f6ccd7e6694a6e4c3b7dfd Mon Sep 17 00:00:00 2001
From: Stamatios Sotiropoulos <stam@fmrib.ox.ac.uk>
Date: Mon, 10 Oct 2011 11:29:35 +0000
Subject: [PATCH] Include noise floor for model2
---
diffmodels.cc | 74 +++++++++++++++++++++++++++++++++++++++------------
diffmodels.h | 12 +++++++--
xfibres.cc | 45 ++++++++++++++++++-------------
3 files changed, 94 insertions(+), 37 deletions(-)
diff --git a/diffmodels.cc b/diffmodels.cc
index cb5c250..09a35a1 100644
--- a/diffmodels.cc
+++ b/diffmodels.cc
@@ -1156,7 +1156,7 @@ boost::shared_ptr<BFMatrix> PVM_single::hess(const NEWMAT::ColumnVector& p,boost
void PVM_multi::fit(){
// initialise with simple pvm
- PVM_single pvm1(Y,bvecs,bvals,nfib);
+ PVM_single pvm1(Y,bvecs,bvals,nfib,m_include_f0);
pvm1.fit();
float _a,_b;
@@ -1168,10 +1168,13 @@ void PVM_multi::fit(){
start(2) = _a;
start(3) = _b;
for(int i=1,ii=4;i<=nfib;i++,ii+=3){
- start(ii) = pvm1.get_f(i);
+ start(ii) = f2x(pvm1.get_f(i)); //Is this a bug?
start(ii+1) = pvm1.get_th(i);
start(ii+2) = pvm1.get_ph(i);
}
+ if (m_include_f0)
+ start(nparams)=f2x(pvm1.get_f0());
+
// do the fit
NonlinParam lmpar(start.Nrows(),NL_LM);
@@ -1192,10 +1195,12 @@ void PVM_multi::fit(){
m_th(k) = final_par(i+1);
m_ph(k) = final_par(i+2);
}
+ if (m_include_f0)
+ m_f0=x2f(final_par(nparams));
sort();
fix_fsum();
-
}
+
void PVM_multi::sort(){
vector< pair<float,int> > fvals(nfib);
ColumnVector ftmp(nfib),thtmp(nfib),phtmp(nfib);
@@ -1211,13 +1216,19 @@ void PVM_multi::sort(){
m_ph(i) = phtmp(fvals[ii].second);
}
}
+
+
void PVM_multi::fix_fsum(){
float sumf=0;
+ if (m_include_f0)
+ sumf=m_f0;
for(int i=1;i<=nfib;i++){
sumf+=m_f(i);
if(sumf>=1){for(int j=i;j<=nfib;j++)m_f(j)=0;break;}
}
}
+
+
ReturnMatrix PVM_multi::get_prediction()const{
ColumnVector pred(npts);
ColumnVector p(nparams);
@@ -1230,10 +1241,14 @@ ReturnMatrix PVM_multi::get_prediction()const{
p(kk+1) = m_th(k);
p(kk+2) = m_ph(k);
}
+ if (m_include_f0)
+ p(nparams)=f2x(m_f0);
pred = forwardModel(p);
pred.Release();
return pred;
}
+
+
NEWMAT::ReturnMatrix PVM_multi::forwardModel(const NEWMAT::ColumnVector& p)const{
ColumnVector pred(npts);
pred = 0;
@@ -1258,11 +1273,17 @@ NEWMAT::ReturnMatrix PVM_multi::forwardModel(const NEWMAT::ColumnVector& p)const
for(int k=1;k<=nfib;k++){
val += fs(k)*anisoterm(i,_a,_b,x.Row(k).t());
}
- pred(i) = p(1)*((1-sumf)*isoterm(i,_a,_b)+val);
+ if (m_include_f0){
+ float temp_f0=x2f(p(nparams));
+ pred(i) = std::abs(p(1))*(temp_f0+(1-sumf-temp_f0)*isoterm(i,_a,_b)+val);
+ }
+ else
+ pred(i) = std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+val);
}
pred.Release();
return pred;
}
+
double PVM_multi::cf(const NEWMAT::ColumnVector& p)const{
//cout<<"CF"<<endl;
//OUT(p.t());
@@ -1287,8 +1308,13 @@ double PVM_multi::cf(const NEWMAT::ColumnVector& p)const{
err = 0.0;
for(int k=1;k<=nfib;k++){
err += fs(k)*anisoterm(i,_a,_b,x.Row(k).t());
+ }
+ if (m_include_f0){
+ float temp_f0=x2f(p(nparams));
+ err = (std::abs(p(1))*(temp_f0+(1-sumf-temp_f0)*isoterm(i,_a,_b)+err) - Y(i));
}
- err = (std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+err) - Y(i));
+ else
+ err = (std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+err) - Y(i));
cfv += err*err;
}
//OUT(cfv);
@@ -1338,12 +1364,21 @@ NEWMAT::ReturnMatrix PVM_multi::grad(const NEWMAT::ColumnVector& p)const{
// ph
J(i,kk+2) = std::abs(p(1))*fs(k)*anisoterm_ph(i,_a,_b,xx,p(kk+1),p(kk+2));
}
- sig = std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+sig);
+ if (m_include_f0){
+ float temp_f0=x2f(p(nparams));
+ //derivative with respect to f0
+ J(i,nparams)= std::abs(p(1))*(1-isoterm(i,_a,_b))*two_pi*sign(p(nparams))*1/(1+p(nparams)*p(nparams));
+ sig=std::abs(p(1))*(temp_f0+(1-sumf-temp_f0)*isoterm(i,_a,_b)+sig);
+ J(i,2) += (p(2)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf-temp_f0)*isoterm_a(i,_a,_b);
+ J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf-temp_f0)*isoterm_b(i,_a,_b);
+ }
+ else{
+ sig = std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+sig);
+ J(i,2) += (p(2)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_a(i,_a,_b);
+ J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_b(i,_a,_b);
+ }
diff(i) = sig - Y(i);
- // other stuff for derivatives
J(i,1) = (p(1)>0?1.0:-1.0)*sig/p(1);
- J(i,2) += (p(2)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_a(i,_a,_b);
- J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_b(i,_a,_b);
}
gradv = 2*J.t()*diff;
@@ -1378,7 +1413,6 @@ boost::shared_ptr<BFMatrix> PVM_multi::hess(const NEWMAT::ColumnVector& p,boost:
}
////////////////////////////////////
Matrix J(npts,nparams);
- ColumnVector diff(npts);
float sig;
for(int i=1;i<=Y.Nrows();i++){
sig = 0;
@@ -1401,16 +1435,22 @@ boost::shared_ptr<BFMatrix> PVM_multi::hess(const NEWMAT::ColumnVector& p,boost:
// ph
J(i,kk+2) = std::abs(p(1))*fs(k)*anisoterm_ph(i,_a,_b,xx,p(kk+1),p(kk+2));
}
- sig = std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+sig);
- diff(i) = sig - Y(i);
- // other stuff for derivatives
+ if (m_include_f0){
+ float temp_f0=x2f(p(nparams));
+ //derivative with respect to f0
+ J(i,nparams)= std::abs(p(1))*(1-isoterm(i,_a,_b))*two_pi*sign(p(nparams))*1/(1+p(nparams)*p(nparams));
+ sig=std::abs(p(1))*(temp_f0+(1-sumf-temp_f0)*isoterm(i,_a,_b)+sig);
+ J(i,2) += (p(2)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf-temp_f0)*isoterm_a(i,_a,_b);
+ J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf-temp_f0)*isoterm_b(i,_a,_b);
+ }
+ else{
+ sig = std::abs(p(1))*((1-sumf)*isoterm(i,_a,_b)+sig);
+ J(i,2) += (p(2)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_a(i,_a,_b);
+ J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_b(i,_a,_b);
+ }
J(i,1) = (p(1)>0?1.0:-1.0)*sig/p(1);
- J(i,2) += (p(2)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_a(i,_a,_b);
- J(i,3) += (p(3)>0?1.0:-1.0)*std::abs(p(1))*(1-sumf)*isoterm_b(i,_a,_b);
-
}
-
for (int i=1; i<=p.Nrows(); i++){
for (int j=i; j<=p.Nrows(); j++){
sig = 0.0;
diff --git a/diffmodels.h b/diffmodels.h
index fab457f..0b79afa 100644
--- a/diffmodels.h
+++ b/diffmodels.h
@@ -461,9 +461,12 @@ class PVM_multi : public PVM, public NonlinCF {
public:
PVM_multi(const ColumnVector& iY,
const Matrix& ibvecs, const Matrix& ibvals,
- const int& nfibres):PVM(iY,ibvecs,ibvals,nfibres){
+ const int& nfibres, bool incl_f0=false):PVM(iY,ibvecs,ibvals,nfibres), m_include_f0(incl_f0){
- nparams = nfib*3 + 3;
+ if (m_include_f0)
+ nparams = nfib*3 + 4;
+ else
+ nparams = nfib*3 + 3;
m_f.ReSize(nfib);
m_th.ReSize(nfib);
@@ -506,10 +509,13 @@ public:
cout << "TH" << ii << " :" << p(i+1) << endl;
cout << "PH" << ii << " :" << p(i+2) << endl;
}
+ if (m_include_f0)
+ cout << "f0 :" << x2f(p(nparams)) << endl;
}
float get_s0()const{return m_s0;}
float get_d()const{return m_d;}
+ float get_f0()const{return m_f0;}
float get_d_std()const{return m_d_std;}
ColumnVector get_f()const{return m_f;}
ColumnVector get_th()const{return m_th;}
@@ -536,10 +542,12 @@ private:
int nparams;
float m_s0;
float m_d;
+ float m_f0;
float m_d_std;
ColumnVector m_f;
ColumnVector m_th;
ColumnVector m_ph;
+ const bool m_include_f0; //Indicate whether f0 will be used in the model (an unattenuated signal compartment)
};
diff --git a/xfibres.cc b/xfibres.cc
index ab9fe40..93fb4b9 100644
--- a/xfibres.cc
+++ b/xfibres.cc
@@ -643,27 +643,36 @@ public:
else{
//////////////////////////////////////////////////////
// model 2 : non-mono-exponential
- if (opts.f0.value())
- m_multifibre.set_f0(0.001); //Need to include f0 in the non-linear initialization of model2 as well!! For now set it initially to almost zero
-
- PVM_multi pvm(m_data,m_bvecs,m_bvals,opts.nfibres.value());
+ float pvmS0, pvmd, pvmd_std, pvmf0=0.001;
+ ColumnVector pvmf,pvmth,pvmph;
+
+ PVM_multi pvm(m_data,m_bvecs,m_bvals,opts.nfibres.value(),opts.f0.value());
pvm.fit();
- m_multifibre.set_S0(pvm.get_s0());
- if(pvm.get_d()>=0)
- m_multifibre.set_d(pvm.get_d());
- else
- m_multifibre.set_d(2e-3);
- if(pvm.get_d_std()>=0)
- m_multifibre.set_d_std(pvm.get_d());
- else
- m_multifibre.set_d(2e-3);
+ pvmf = pvm.get_f(); pvmth = pvm.get_th(); pvmph = pvm.get_ph(); pvmd_std=pvm.get_d_std();
+ pvmS0 = pvm.get_s0(); pvmd = pvm.get_d(); predicted_signal=pvm.get_prediction();
+
+ if (opts.f0.value()){
+ pvmf0=pvm.get_f0();
+ //If the full model gives values that are implausible, or we are in a CSF voxel (f1<0.05)
+ //then fit a model without the f0 and drive f0_init to almost zero
+ if ((opts.nfibres.value()>0 && pvmf(1)<0.05) || pvmd>0.007 || pvmf0>0.4){
+ PVM_multi pvm2(m_data,m_bvecs,m_bvals,opts.nfibres.value(),false);
+ pvm2.fit();
+ pvmf0=0.001; pvmS0=pvm2.get_s0(); pvmd=pvm2.get_d(); pvmd_std=pvm2.get_d_std();
+ pvmf = pvm2.get_f(); pvmth = pvm2.get_th(); pvmph = pvm2.get_ph();
+ predicted_signal=pvm2.get_prediction();
+ }
+ m_multifibre.set_f0(pvmf0);
+ }
- ColumnVector pvmf,pvmth,pvmph;
- pvmf = pvm.get_f();
- pvmth = pvm.get_th();
- pvmph = pvm.get_ph();
+ if(pvmd<0) pvmd=2e-3;
+ if(pvmd_std<0) pvmd_std=2e-3;
+ m_multifibre.set_S0(pvmS0);
+ m_multifibre.set_d(pvmd);
+ m_multifibre.set_d_std(pvmd_std);
+
if(opts.nfibres.value()>0){
m_multifibre.addfibre(pvmth(1),
pvmph(1),
@@ -677,7 +686,7 @@ public:
}
}
- residuals=m_data-pvm.get_prediction();
+ residuals=m_data-predicted_signal;
}
residuals.Release();
return residuals;
--
GitLab