diff --git a/probtrackxOptions.h b/probtrackxOptions.h
index d11a53cab2b2bfedfbed0a96c68ededea80c8c96..24f39b38d27929c7b9289c35b57120ebb7dfed0b 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 6eb14f38074488f4aef9b8f59d7ef8860af30796..83e2a867f8fe186cf6c5fd8265dde68db72e3eb5 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 9e5cb183745d2b3cb0684909ebf4d57c3e433bb2..354d2817fae003499bae0bc01dca314862b819e0 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 29203e204c20eb22358cff9e71cd5ec06096bb8b..f1206e5afdb74275f414451cb97ea7bbfd2943e5 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{