From e578db22a6a4118148a4f9b1d85e6065528e9f86 Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Thu, 16 Jun 2005 13:17:33 +0000
Subject: [PATCH] *** empty log message ***
---
fibre.h | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++----
xfibres.cc | 7 ++++++
2 files changed, 71 insertions(+), 4 deletions(-)
diff --git a/fibre.h b/fibre.h
index 99fba5d..9bcf0fe 100644
--- a/fibre.h
+++ b/fibre.h
@@ -384,6 +384,13 @@ namespace FIBRE{
float m_d_old_prior;
float m_d_acc;
float m_d_rej;
+ float m_diso;
+ float m_diso_old;
+ float m_diso_prop;
+ float m_diso_prior;
+ float m_diso_old_prior;
+ float m_diso_acc;
+ float m_diso_rej;
float m_S0;
float m_S0_old;
float m_S0_prop;
@@ -429,6 +436,7 @@ namespace FIBRE{
void initialise_energies(){
compute_d_prior();
+ compute_diso_prior();
compute_S0_prior();
m_prior_en=0;
compute_prior();
@@ -440,19 +448,23 @@ namespace FIBRE{
void initialise_props(){
m_S0_prop=m_S0/10; //must have set inital values before this is called;
m_d_prop=m_d/10;
+ m_diso_prop=m_diso/10;
}
inline float get_d() const{ return m_d;}
inline void set_d(const float d){ m_d=d; }
-
+ inline float get_diso() const{ return m_diso;}
+ inline void set_diso(const float diso){ m_diso=diso; }
+
+
inline float get_S0() const{ return m_S0;}
inline void set_S0(const float S0){ m_S0=S0; }
inline bool compute_d_prior(){
m_d_old_prior=m_d_prior;
- if(m_d<0)
+ if(m_d<0||m_d>5e-3)
return true;
else{
m_d_prior=0;
@@ -460,6 +472,17 @@ namespace FIBRE{
}
}
+
+ inline bool compute_diso_prior(){
+ m_diso_old_prior=m_diso_prior;
+ if(m_diso<0||m_diso>5e-3)
+ return true;
+ else{
+ m_diso_prior=0;
+ return false;
+ }
+ }
+
inline bool compute_S0_prior(){
m_S0_old_prior=m_S0_prior;
if(m_S0<0) return true;
@@ -480,7 +503,7 @@ namespace FIBRE{
inline void compute_prior(){
m_old_prior_en=m_prior_en;
- m_prior_en=m_d_prior+m_S0_prior;
+ m_prior_en=m_d_prior+m_diso_prior+m_S0_prior;
for(unsigned int f=0;f<m_fibres.size(); f++){
m_prior_en=m_prior_en+m_fibres[f].get_prior();
}
@@ -497,7 +520,7 @@ namespace FIBRE{
}
for(int i=1;i<=pred.Nrows();i++){
- pred(i)=pred(i)+(1-fsum)*exp(-m_d*m_bvals(1,i));
+ pred(i)=pred(i)+(1-fsum)*exp(-m_diso*m_bvals(1,i));
}
pred=pred*m_S0;
@@ -549,6 +572,26 @@ namespace FIBRE{
m_d_rej++;
}
+
+ inline bool propose_diso(){
+ m_diso_old=m_diso;
+ m_diso+=normrnd().AsScalar()*m_diso_prop;
+ bool rejflag=compute_diso_prior();//inside this it stores the old prior
+ return rejflag;
+ };
+
+ inline void accept_diso(){
+ m_diso_acc++;
+ }
+
+ inline void reject_diso(){
+ m_diso=m_diso_old;
+ m_diso_prior=m_diso_old_prior;
+ m_prior_en=m_old_prior_en;
+ m_diso_rej++;
+ }
+
+
inline bool propose_S0(){
m_S0_old=m_S0;
m_S0+=normrnd().AsScalar()*m_S0_prop;
@@ -569,6 +612,7 @@ namespace FIBRE{
inline void update_proposals(){
m_d_prop*=sqrt(float(m_d_acc+1)/float(m_d_rej+1));
+ m_diso_prop*=sqrt(float(m_diso_acc+1)/float(m_diso_rej+1));
m_S0_prop*=sqrt(float(m_S0_acc+1)/float(m_S0_rej+1));
m_d_acc=0;
m_d_rej=0;
@@ -602,6 +646,22 @@ namespace FIBRE{
reject_d();
}
+ if(!propose_diso()){
+ compute_prior();
+ compute_likelihood();
+ compute_energy();
+ if(test_energy()){
+ accept_diso();
+ }
+ else{
+ reject_diso();
+ restore_energy();
+ }
+ }
+ else{
+ reject_diso();
+ }
+
if(!propose_S0()){
compute_prior();
compute_likelihood();
diff --git a/xfibres.cc b/xfibres.cc
index ca16aea..7cf8eb8 100644
--- a/xfibres.cc
+++ b/xfibres.cc
@@ -127,6 +127,7 @@ inline SymmetricMatrix vec2tens(ColumnVector& Vec){
class Samples{
xfibresOptions& opts;
Matrix m_dsamples;
+ Matrix m_disosamples;
Matrix m_S0samples;
vector<Matrix> m_thsamples;
vector<Matrix> m_phsamples;
@@ -148,6 +149,8 @@ public:
m_dsamples.ReSize(nsamples,nvoxels);
m_dsamples=0;
+ m_disosamples.ReSize(nsamples,nvoxels);
+ m_disosamples=0;
m_S0samples.ReSize(nsamples,nvoxels);
m_S0samples=0;
@@ -164,6 +167,7 @@ public:
void record(Multifibre& mfib, int vox, int samp){
m_dsamples(samp,vox)=mfib.get_d();
+ m_disosamples(samp,vox)=mfib.get_diso();
m_S0samples(samp,vox)=mfib.get_S0();
for(int f=0;f<opts.nfibres.value();f++){
m_thsamples[f](samp,vox)=mfib.fibres()[f].get_th();
@@ -186,6 +190,8 @@ public:
Log& logger = LogSingleton::getInstance();
tmp.setmatrix(m_dsamples,mask);
save_volume4D(tmp,logger.appendDir("dsamples"));
+ tmp.setmatrix(m_disosamples,mask);
+ save_volume4D(tmp,logger.appendDir("disosamples"));
tmp.setmatrix(m_S0samples,mask);
save_volume4D(tmp,logger.appendDir("S0samples"));
@@ -327,6 +333,7 @@ class xfibresVoxelManager{
if(f<=0.001) f=0.001;
if(D<=0) D=2e-3;
m_multifibre.set_d(D);
+ m_multifibre.set_diso(D);
m_multifibre.set_S0(S0);
if(opts.nfibres.value()>0){
m_multifibre.addfibre(th,ph,f,1);
--
GitLab