diff --git a/Makefile b/Makefile
index 346add5f9edb219806e9d715ed4d96da59940416..37c61deb6156869712c4039fb2f88d85aadabc75 100644
--- a/Makefile
+++ b/Makefile
@@ -12,6 +12,7 @@ DLIBS = -lbint -lnewimage -lutils -lmiscmaths  -lnewmat -lfslio -lniftiio -lznz
 
 DTIFIT=dtifit
 PT=probtrack
+PTX=probtrackx
 FTB=find_the_biggest
 PJ=proj_thresh
 MED=medianfilter
@@ -25,6 +26,7 @@ TEST=testfile
 
 DTIFITOBJS=dtifit.o dtifitOptions.o
 PTOBJS=probtrack.o probtrackOptions.o pt_alltracts.o pt_matrix.o pt_seeds_to_targets.o pt_simple.o pt_twomasks.o pt_matrix_mesh.o
+PTXOBJS=probtrackx.o probtrackxOptions.o streamlines.o ptx_simple.o ptx_seedmask.o ptx_twomasks.o
 FTBOBJS=find_the_biggest.o
 PJOBJS=proj_thresh.o
 MEDOBJS=medianfilter.o 
@@ -46,6 +48,9 @@ RUNTCLS = Fdt
 
 all: ${XFILES} ${FXFILES} 
 
+${PTX}:		   ${PTXOBJS}
+		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${PTXOBJS} ${DLIBS}
+
 ${PT}:		   ${PTOBJS}
 		   ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${PTOBJS} ${DLIBS} 
 
diff --git a/particle.cc b/particle.cc
index cdecaf52e04542984566772b29cd3541d6319a54..b7c3f67228635256d3c47f808c09b98537542c01 100644
--- a/particle.cc
+++ b/particle.cc
@@ -1,4 +1,4 @@
-/*  Gam.cc
+/*  particle.cc
 
     Tim Behrens, FMRIB Image Analysis Group
 
@@ -7,10 +7,3 @@
 /*  CCOPYRIGHT  */
 
 #include "particle.h"
-
-
-
-
-
-
-
diff --git a/particle.h b/particle.h
index f30696dbd07ff6ceea87793f3b0da25a8eedb1a7..6ea2118a3d5d080dde75b1eb11dfd56dfb3b75fb 100644
--- a/particle.h
+++ b/particle.h
@@ -46,25 +46,51 @@ namespace PARTICLE{
       int m_jumpsign;
     public:
       //constructors::
-      Particle(const float& xin=0,const float& yin=0,
-	       const float& zin=0,const float& rxin=0,
-	       const float& ryin=0,const float &rzin=0,
-	       const float& steplengthin=0.5,
-	       const float& xdimin=2,
-	       const float& ydimin=2,
-	       const float& zdimin=2,
+      Particle(const float& xin,const float& yin,
+	       const float& zin,const float& rxin,
+	       const float& ryin,const float &rzin,
+	       const float& steplengthin,
+	       const float& xdimin,
+	       const float& ydimin,
+	       const float& zdimin,
 	       const bool& hasjumpedin=false,
 	       const bool& simdiffin=false) : 
 	m_x(xin), m_y(yin), m_z(zin), m_rx(rxin), 
 	m_ry(ryin),m_rz(rzin),m_steplength(steplengthin),
 	m_xdim(xdimin),m_ydim(ydimin),m_zdim(zdimin),
 	m_has_jumped(hasjumpedin),m_simdiff(false){}
+      Particle(){}
       ~Particle(){}
-
+      
+      //initialise
+      void initialise(const float& xin=0,const float& yin=0,
+		      const float& zin=0,const float& rxin=0,
+		      const float& ryin=0,const float &rzin=0,
+		 const float& steplengthin=0.5,
+		 const float& xdimin=2,
+		 const float& ydimin=2,
+		 const float& zdimin=2,
+		 const bool& hasjumpedin=false,
+		 const bool& simdiffin=false){
+	m_x=xin;
+	m_y=yin;
+	m_z=zin;
+	m_rx=rxin; 
+	m_ry=ryin;
+	m_rz=rzin;
+	m_steplength=steplengthin;
+	m_xdim=xdimin;
+	m_ydim=ydimin;
+	m_zdim=zdimin;
+	m_has_jumped=hasjumpedin;
+	m_simdiff=simdiffin;
+      }
+      
+      
       //return values
       const float& x() const { return m_x; }
       float x() { return m_x; }
-  
+      
       const float& y() const { return m_y; }
       float y() { return m_y; }
   
@@ -210,6 +236,9 @@ namespace PARTICLE{
     return ostr;
   }
 
+  
+
+
 }
 
 #endif
diff --git a/probtrackx.cc b/probtrackx.cc
index dfffb8c78569064d386f8d8dfbb8315f867a08e0..4c9cf199a9fae8e3f8fc821c5062e7424d47cedd 100644
--- a/probtrackx.cc
+++ b/probtrackx.cc
@@ -24,7 +24,7 @@ using namespace mesh;
 
 
 int main ( int argc, char **argv ){
-  probtrackOptions& opts =probtrackOptions::getInstance();
+  probtrackxOptions& opts =probtrackxOptions::getInstance();
   Log& logger = LogSingleton::getInstance();
   opts.parse_command_line(argc,argv,logger);
   srand(opts.rseed.value());
diff --git a/probtrackx.h b/probtrackx.h
index 856a6d0882952e4ab4096cce177a9be2eaeb5a5d..1e2fb56138b6ec69c630b135ea065dcb817ccc74 100644
--- a/probtrackx.h
+++ b/probtrackx.h
@@ -2,7 +2,7 @@
 
 /*  CCOPYRIGHT  */
 
-#include "probtrackOptions.h"
+#include "probtrackxOptions.h"
 #include "particle.h"
 #include "tractvols.h"
 #include "ptx_simple.h"
diff --git a/ptx_seedmask.cc b/ptx_seedmask.cc
index e755cf6d245cee7798f1f40ce86724de89f28144..9539738b06bf16baca008f14b98427cf44fcb6c2 100644
--- a/ptx_seedmask.cc
+++ b/ptx_seedmask.cc
@@ -17,7 +17,7 @@ using namespace mesh;
 
 void seedmask()
 { 
-  probtrackOptions& opts =probtrackOptions::getInstance();
+  probtrackxOptions& opts =probtrackxOptions::getInstance();
 
   ////////////////////////////////
   //  Log& logger = LogSingleton::getInstance();
diff --git a/ptx_seedmask.h b/ptx_seedmask.h
index 9fc135aa47ed1a0ee0cbedadefbc3296413732db..dbc49682f1ddeea7276d637b19a7794bf096ebd7 100644
--- a/ptx_seedmask.h
+++ b/ptx_seedmask.h
@@ -6,7 +6,7 @@
 #include "newimage/newimageall.h"
 #include "utils/log.h"
 #include "meshclass/meshclass.h"
-#include "probtrackOptions.h"
+#include "probtrackxOptions.h"
 #include "particle.h"
 #include "tractvols.h"
 
diff --git a/ptx_simple.cc b/ptx_simple.cc
index 3425efb33a70adcd2cbbc392beb556b7e7ead60d..fd47364c23f6f51fcd7e48ba8a936be643029336 100644
--- a/ptx_simple.cc
+++ b/ptx_simple.cc
@@ -15,7 +15,7 @@ using namespace mesh;
 
 
 void track(){
-  probtrackOptions& opts =probtrackOptions::getInstance();
+  probtrackxOptions& opts =probtrackxOptions::getInstance();
   
   ////////////////////////////
   Log& logger = LogSingleton::getInstance();
diff --git a/ptx_simple.h b/ptx_simple.h
index 3f8590eb1b851fdafcba2948b79934763b1a3c6d..eb57c1fa92ec7b9f5c4a503fa10327c9c7bc21e9 100644
--- a/ptx_simple.h
+++ b/ptx_simple.h
@@ -6,7 +6,7 @@
 #include "newimage/newimageall.h"
 #include "utils/log.h"
 #include "meshclass/meshclass.h"
-#include "probtrackOptions.h"
+#include "probtrackxOptions.h"
 #include "particle.h"
 #include "tractvols.h"
 
diff --git a/ptx_twomasks.cc b/ptx_twomasks.cc
index 1abcde0bc705540ced33e4eba424401f85606980..6010297599238dbebce6fda1f4b39c85f5687503 100644
--- a/ptx_twomasks.cc
+++ b/ptx_twomasks.cc
@@ -17,7 +17,7 @@ using namespace mesh;
 
 void twomasks()
 { 
-  probtrackOptions& opts =probtrackOptions::getInstance();
+  probtrackxOptions& opts =probtrackxOptions::getInstance();
 
   ////////////////////////////////
   //  Log& logger = LogSingleton::getInstance();
diff --git a/ptx_twomasks.h b/ptx_twomasks.h
index 626d00023a2d2d0c6a93c06901e0154e828705ef..40fce4a2dceccb57883fb46d4a4a3bd1427b4ef7 100644
--- a/ptx_twomasks.h
+++ b/ptx_twomasks.h
@@ -6,7 +6,7 @@
 #include "newimage/newimageall.h"
 #include "utils/log.h"
 #include "meshclass/meshclass.h"
-#include "probtrackOptions.h"
+#include "probtrackxOptions.h"
 #include "particle.h"
 #include "tractvols.h"
 
diff --git a/streamlines.cc b/streamlines.cc
index c32eab29e49e6cda17e20160b82539596afbab01..aa03cccb714fe0bfd1f6f86e29a210afd3fb2fe4 100644
--- a/streamlines.cc
+++ b/streamlines.cc
@@ -23,7 +23,7 @@ namespace TRACT{
   }
   
   
-  Streamliner::Streamliner():opts(probtrackOptions::getInstance()),logger(LogSingleton::getInstance()),
+  Streamliner::Streamliner():opts(probtrackxOptions::getInstance()),logger(LogSingleton::getInstance()),
 			     vols(opts.usef.value()){
     
     read_volume(m_mask,opts.maskfile.value());
diff --git a/streamlines.h b/streamlines.h
index d2869a61e8e155a25b81009434eb49516fdd7a4b..c7bf3497e2027e10e289e193332b331e6f20aba4 100644
--- a/streamlines.h
+++ b/streamlines.h
@@ -2,7 +2,7 @@
 #include "newimage/newimageall.h"
 #include "utils/log.h"
 #include "meshclass/meshclass.h"
-#include "probtrackOptions.h"
+#include "probtrackxOptions.h"
 #include "particle.h"
 #include "tractvols.h"
 
@@ -21,7 +21,7 @@ namespace TRACT{
     //Everything in DTI space is done INSIDE this class and lower level classes (particle and tractvols)
     //This class communicates with higher level classes in Seed voxels.
     //
-    probtrackOptions& opts;
+    probtrackxOptions& opts;
     Log& logger;
     Particle m_part;
     vector<ColumnVector> m_path;
@@ -80,7 +80,7 @@ namespace TRACT{
 
 
   class Counter{
-    probtrackOptions& opts;
+    probtrackxOptions& opts;
     Log& logger;
     volume<int> m_prob;
     volume<int> m_beenhere;
@@ -113,7 +113,7 @@ namespace TRACT{
     Streamliner& m_nonconst_stline;
     
   public:
-    Counter(const volume<int>& seeds,Streamliner& stline):opts(probtrackOptions::getInstance()),
+    Counter(const volume<int>& seeds,Streamliner& stline):opts(probtrackxOptions::getInstance()),
 							  logger(LogSingleton::getInstance()),
 							  m_seeds(seeds),m_stline(stline),
 							  m_nonconst_stline(stline){
@@ -180,14 +180,14 @@ namespace TRACT{
   };
   
   class Seedmanager{
-    probtrackOptions& opts;
+    probtrackxOptions& opts;
     Log& logger;
     Counter& m_counter;    
     Streamliner& m_stline;
     const volume<int>& m_seeds;
     ColumnVector m_seeddims;
   public:
-    Seedmanager(Counter& counter):opts(probtrackOptions::getInstance()),
+    Seedmanager(Counter& counter):opts(probtrackxOptions::getInstance()),
 				  logger(LogSingleton::getInstance()),
 				  m_counter(counter),
 				  m_stline(m_counter.get_nonconst_streamline()),
diff --git a/tractvols.h b/tractvols.h
index 4eb9da744fa2a351ae064eeef918ba64d1391617..b7a9b2755a333e59daf149670b12992d8d8c2f5f 100644
--- a/tractvols.h
+++ b/tractvols.h
@@ -1,4 +1,4 @@
-/*  tractvols.h
+/*  tractvolsX.h
 
     Tim Behrens, FMRIB Image Analysis Group
 
@@ -16,36 +16,93 @@
 #include "newimage/newimageall.h"
 #include <iostream>
 #include "stdlib.h"
-
+#include "probtrackxOptions.h"
 using namespace std;
 using namespace NEWIMAGE;
+using namespace TRACT;
 
 namespace TRACTVOLS{
   class TractVols
     {
     private:
-      volume4D<float> thsamples;
-      volume4D<float> phsamples;
-      volume4D<float> fsamples;
+      probtrackxOptions& opts;
+      vector<volume4D<float>* > thsamples;
+      vector<volume4D<float>* > phsamples;
+      vector<volume4D<float>* > fsamples;
+      bool init_sample;
+      int fibst;
       bool usef;
     public:
       //constructors::
-      TractVols(const bool& usefin=false):usef(usefin){}
-      ~TractVols(){}
+      TractVols(const bool& usefin=false):opts(probtrackxOptions::getInstance()),init_sample(true),fibst(1),usef(usefin){}
+      TractVols():opts(probtrackxOptions::getInstance()){}
+      ~TractVols(){
+	for(unsigned int m=0;m<thsamples.size();m++)
+	  delete thsamples[m]; //ask flitney, do you just delete the ptr??
+	for(unsigned int m=0;m<phsamples.size();m++)
+	  delete phsamples[m];
+	for(unsigned int m=0;m<fsamples.size();m++)
+	  delete fsamples[m];
+      }
       
+      
+      void reset(const int& fibst_in){
+	init_sample=true;
+	fibst=fibst_in;
+      }
       //Initialise
       void initialise(const string& basename){
-	read_volume4D(thsamples,basename+"_thsamples");
-	read_volume4D(phsamples,basename+"_phsamples");
-	if(usef)
-	  read_volume4D(fsamples,basename+"_fsamples");
+	
+
+	if(fsl_imageexists(basename+"_thsamples")){
+	  volume4D<float> *tmpthptr= new volume4D<float>;
+	  volume4D<float> *tmpphptr= new volume4D<float>;
+	  volume4D<float> *tmpfptr= new volume4D<float>;
+	  cout<<"1"<<endl;
+	  read_volume4D(*tmpthptr,basename+"_thsamples");
+	  cout<<"2"<<endl;
+	  thsamples.push_back(tmpthptr);
+	  cout<<"3"<<endl;
+	  read_volume4D(*tmpphptr,basename+"_phsamples");
+	  cout<<"4"<<endl;
+	  phsamples.push_back(tmpphptr);
+	  cout<<"5"<<endl;
+	  if(usef){
+	    read_volume4D(*tmpfptr,basename+"_fsamples");
+	    fsamples.push_back(tmpfptr);
+	  }
+	  cout<<"6"<<endl;
+	}
+	else{
+	  int fib=1;
+	  bool fib_existed=true;
+	  while(fib_existed){
+	    if(fsl_imageexists(basename+"_th"+num2str(fib)+"samples")){
+	      volume4D<float> *tmpthptr= new volume4D<float>;
+	      volume4D<float> *tmpphptr= new volume4D<float>;
+	      volume4D<float> *tmpfptr= new volume4D<float>;
+	      cout<<fib<<"_1"<<endl;
+	      read_volume4D(*tmpthptr,basename+"_th"+num2str(fib)+"samples");
+	      thsamples.push_back(tmpthptr);
+	      cout<<fib<<"_2"<<endl;
+	      read_volume4D(*tmpphptr,basename+"_ph"+num2str(fib)+"samples");
+	      phsamples.push_back(tmpphptr);
+	      cout<<fib<<"_3"<<endl;
+	      read_volume4D(*tmpfptr,basename+"_f"+num2str(fib)+"samples");
+	      fsamples.push_back(tmpfptr);
+	      fib++;
+	    }
+	    else{
+	      fib_existed=false;
+	    }
+	  }
+	  
+	}
+	cout<<"7"<<endl;
       }
       
-
-      ColumnVector sample(const float& x,const float& y,const float&z){
-	// 	int r_x=(int) round(x);
-	// 	int r_y=(int) round(y);
-	// 	int r_z=(int) round(z);
+      
+      ColumnVector sample(const float& x,const float& y,const float&z,const float& r_x,const float& r_y,const float& r_z){
 	
 	////////Probabilistic interpolation
 	int cx =(int) ceil(x),fx=(int) floor(x);
@@ -93,17 +150,49 @@ namespace TRACTVOLS{
 	
 	
 	float samp=rand(); samp/=RAND_MAX;
-	samp=round(samp*(thsamples.tsize()-1));
-	//float phi = phsamples(r_x,r_y,r_z,samp);
-	//float theta = thsamples(r_x,r_y,r_z,samp);
+	samp=round(samp*((*thsamples[0]).tsize()-1));
+	float theta=0,phi=0;
+	float dotmax=0,dottmp=0;
+	int fibind=0;
+	if(thsamples.size()>1){//more than 1 fibre
+	  if(init_sample){//go for the specified option on the first jump
+	    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{
+	    for(unsigned int fib=0;fib<thsamples.size();fib++){
+	      if((*fsamples[fib])(int(newx),int(newy),int(newz),int(samp))>opts.fibthresh.value()){
+		float phtmp=(*phsamples[fib])(int(newx),int(newy),int(newz),int(samp));
+		float thtmp=(*thsamples[fib])(int(newx),int(newy),int(newz),int(samp));
+		dottmp=fabs(sin(thtmp)*cos(phtmp)*r_x + sin(thtmp)*sin(phtmp)*r_y + cos(thtmp)*r_z);
+		if(dottmp>dotmax){
+		  dotmax=dottmp;
+		  theta=thtmp;
+		  phi=phtmp;
+		  fibind=fib;
+		}
 		
-	float phi = phsamples(int(newx),int(newy),int(newz),int(samp));
-	float theta = thsamples(int(newx),int(newy),int(newz),int(samp));
+	      }
+	      
+	    }
+	    
+	    if(dotmax==0){
+	      theta=(*thsamples[0])(int(newx),int(newy),int(newz),int(samp));
+	      phi=(*phsamples[0])(int(newx),int(newy),int(newz),int(samp));
+	    }
+	  }
+	}
+	else{
+	  theta=(*thsamples[0])(int(newx),int(newy),int(newz),int(samp));
+	  phi=(*phsamples[0])(int(newx),int(newy),int(newz),int(samp));
+	}
+
 	
 	float f;
 	
 	if(usef){
-	  f = fsamples(int(newx),int(newy),int(newz),int(samp));
+	  f = (*fsamples[fibind])(int(newx),int(newy),int(newz),int(samp));
 	}
 	else
 	  f=1;
@@ -116,7 +205,7 @@ namespace TRACTVOLS{
 
       ColumnVector dimensions(){
 	ColumnVector dims(3);
-	dims << thsamples.xdim() <<thsamples.ydim() << thsamples.zdim();
+	dims << (*thsamples[0]).xdim() <<(*thsamples[0]).ydim() << (*thsamples[0]).zdim();
 	return dims;
       }
     };