*** empty log message ***

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();

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);