From eebba49123a759843a598647a5cd09e66f429c2c Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Thu, 22 Oct 2009 08:28:43 +0000
Subject: [PATCH] Fixed randfib, fibst bugs. Updated options text to make it
 clearer what these things do..

---
 probtrackxOptions.h |  4 ++--
 streamlines.cc      | 14 ++++++++++----
 streamlines.h       |  1 +
 tractvolsx.h        | 15 +++++----------
 4 files changed, 18 insertions(+), 16 deletions(-)

diff --git a/probtrackxOptions.h b/probtrackxOptions.h
index d11a53c..24f39b3 100644
--- a/probtrackxOptions.h
+++ b/probtrackxOptions.h
@@ -213,10 +213,10 @@ class probtrackxOptions {
 	 string("Use anisotropy to constrain tracking"), 
 	 false, no_argument),
   randfib(string("--randfib"), 0, 
-	 string("Set to 1 to randomly sample initial fibres. Set to 2 to sample in proportion to f"), 
+	 string("Default 0. Set to 1 to randomly sample initial fibres (with f > fibthresh). \n                        Set to 2 to sample in proportion fibres (with f>fibthresh) to f. \n                        Set to 3 to sample ALL populations at random (even if f<fibthresh)"), 
 	 false, no_argument),
   fibst(string("--fibst"),1, 
-	 string("Force a starting fibre for tracking - default=1, i.e. first fibre orientation"), 
+	 string("Force a starting fibre for tracking - default=1, i.e. first fibre orientation. Only works if randfib==0"), 
 	 false, requires_argument),
   modeuler(string("--modeuler"), false, 
 	   string("Use modified euler streamlining"), 
diff --git a/streamlines.cc b/streamlines.cc
index 6eb14f3..83e2a86 100644
--- a/streamlines.cc
+++ b/streamlines.cc
@@ -1098,14 +1098,20 @@ void Counter::save_matrix3(){
 
     int nlines=0;
     for(int p=0;p<opts.nparticles.value();p++){
-      if(opts.randfib.value()){
+      
+      if(opts.randfib.value()>2){ 
+	//This bit of code just does random sampling from all fibre populations - even those whose f value is less than f-thresh. 
+	//3 other possibilities - randfib==0 -> use fibst (default first fibre but can be set)
+	// randfib==1 - random sampling of fibres bigger than fthresh
+	// randfib==2 random sampling of fibres bigger than fthresh in proporthion to their f-values. 
 	float tmp=rand()/float(RAND_MAX) * float(m_stline.nfibres()-1);
 	fibst = (int)round(tmp);
-      } 
+      }
+      
       if(opts.verbose.value()>1)
 	logger.setLogFile("particle"+num2str(p));
-      
-      m_stline.reset();
+   
+      m_stline.reset(); //This now includes a vols.reset() in order to get fibst right. 
       bool forwardflag=false,backwardflag=false;
       bool counted=false;
       if(!onewayonly || opts.matrix3out.value()){//always go both ways in matrix3 mode
diff --git a/streamlines.h b/streamlines.h
index 9e5cb18..354d281 100644
--- a/streamlines.h
+++ b/streamlines.h
@@ -106,6 +106,7 @@ namespace TRACT{
     inline vector<ColumnVector> get_path() const{return m_path;}
     inline void reset(){
       m_part.reset();
+      vols.reset(opts.fibst.value());
     }
     inline void reverse(){
       m_part.restart_reverse();
diff --git a/tractvolsx.h b/tractvolsx.h
index 29203e2..f1206e5 100644
--- a/tractvolsx.h
+++ b/tractvolsx.h
@@ -153,9 +153,8 @@ namespace TRACTVOLSX{
 	float dotmax=0,dottmp=0;
 	int fibind=0;
 	if(thsamples.size()>1){//more than 1 fibre
-
+	  
 	  if(init_sample){//go for the specified fibre on the first jump or generate at random
-	    
 	    if(opts.randfib.value()==1){//this generates startfib at random (except for fibres where f<fibthresh)
 	      vector<int> fibvec;
 	      for(unsigned int fib=0;fib<thsamples.size();fib++){	    
@@ -172,10 +171,7 @@ namespace TRACTVOLSX{
 		float rtmp=rand()/float(RAND_MAX) * float(fibvec.size()-1);
 		fibst = fibvec[ (int)round(rtmp) ];	      
 	      }
-	      
 	    }
-
-
 	    else if(opts.randfib.value()==2){ //this generates startfib with probability proportional to f (except for fibres where f<fibthresh). 
 	      //this chooses at random but in proportion to fsamples. 
 	      float fsumtmp=0;
@@ -204,11 +200,10 @@ namespace TRACTVOLSX{
 		}
 		
 	      }
-	      
-	      theta=(*thsamples[fibst])(int(newx),int(newy),int(newz),int(samp));
-	      phi=(*phsamples[fibst])(int(newx),int(newy),int(newz),int(samp));
-	      
-	    }
+	    }  
+
+	    theta=(*thsamples[fibst])(int(newx),int(newy),int(newz),int(samp));
+	    phi=(*phsamples[fibst])(int(newx),int(newy),int(newz),int(samp));
 	    init_sample=false;
 	  }
 	  else{
-- 
GitLab