*** empty log message ***

probtrackx.cc

@@ -15,7 +15,6 @@ using namespace NEWIMAGE;
 using namespace TRACT;
 using namespace Utilities;
 using namespace PARTICLE;
-using namespace TRACTVOLS;
 using namespace mesh;
 //using namespace NEWMAT;
 //////////////////////////

probtrackx.h

@@ -4,7 +4,6 @@
 
 #include "probtrackxOptions.h"
 #include "particle.h"
-#include "tractvols.h"
 #include "ptx_simple.h"
 #include "ptx_seedmask.h"
 #include "ptx_twomasks.h"

ptx_seedmask.cc

@@ -10,7 +10,6 @@ using namespace NEWIMAGE;
 using namespace TRACT;
 using namespace Utilities;
 using namespace PARTICLE;
-using namespace TRACTVOLS;
 using namespace mesh;
 
 

ptx_seedmask.h

@@ -8,7 +8,7 @@
 #include "meshclass/meshclass.h"
 #include "probtrackxOptions.h"
 #include "particle.h"
-#include "tractvols.h"
+
 
 
 

ptx_simple.cc

@@ -10,7 +10,6 @@ using namespace NEWIMAGE;
 using namespace TRACT;
 using namespace Utilities;
 using namespace PARTICLE;
-using namespace TRACTVOLS;
 using namespace mesh;
 
 

ptx_simple.h

@@ -8,7 +8,7 @@
 #include "meshclass/meshclass.h"
 #include "probtrackxOptions.h"
 #include "particle.h"
-#include "tractvols.h"
+
 
 
 

ptx_twomasks.cc

@@ -10,7 +10,6 @@ using namespace NEWIMAGE;
 using namespace TRACT;
 using namespace Utilities;
 using namespace PARTICLE;
-using namespace TRACTVOLS;
 using namespace mesh;
 
 

ptx_twomasks.h

@@ -8,6 +8,6 @@
 #include "meshclass/meshclass.h"
 #include "probtrackxOptions.h"
 #include "particle.h"
-#include "tractvols.h"
+
 
 void twomasks();

streamlines.cc

@@ -73,7 +73,7 @@ namespace TRACT{
   
   bool Streamliner::streamline(const float& x_init,const float& y_init,const float& z_init, const ColumnVector& dim_seeds,const int& fibst){ 
     
-    //fibst tells tractvols which fibre to start with if there are more than one..
+    //fibst tells tractvolsx which fibre to start with if there are more than one..
     //x_init etc. are in seed space...
     vols.reset(fibst);
     m_x_s_init=x_init; //seed x position in voxels

streamlines.h

@@ -4,12 +4,12 @@
 #include "meshclass/meshclass.h"
 #include "probtrackxOptions.h"
 #include "particle.h"
-#include "tractvols.h"
+#include "tractvolsx.h"
 
 using namespace std;
 using namespace NEWIMAGE;
 using namespace Utilities;
-using namespace TRACTVOLS;
+using namespace TRACTVOLSX;
 using namespace mesh;
 using namespace PARTICLE;
 
@@ -18,7 +18,7 @@ namespace TRACT{
 
 
   class Streamliner{
-    //Everything in DTI space is done INSIDE this class and lower level classes (particle and tractvols)
+    //Everything in DTI space is done INSIDE this class and lower level classes (particle and tractvolsx)
     //This class communicates with higher level classes in Seed voxels.
     //
     probtrackxOptions& opts;
@@ -33,7 +33,7 @@ namespace TRACT{
     vector<bool> m_passed_flags;
     vector<bool> m_own_waymasks;
     Matrix m_Seeds_to_DTI;
-    TractVols vols;
+    Tractvolsx vols;
     float m_lcrat;
     float m_x_s_init;
     float m_y_s_init;

tractvols.h

-/*  tractvolsX.h
+/*  tractvols.h
 
     Tim Behrens, FMRIB Image Analysis Group
 
     Copyright (C) 2004 University of Oxford  */
 
-/*  CCOPYRIGHT  */
+/*  Part of FSL - FMRIB's Software Library
+    http://www.fmrib.ox.ac.uk/fsl
+    fsl@fmrib.ox.ac.uk
+    
+    Developed at FMRIB (Oxford Centre for Functional Magnetic Resonance
+    Imaging of the Brain), Department of Clinical Neurology, Oxford
+    University, Oxford, UK
+    
+    
+    LICENCE
+    
+    FMRIB Software Library, Release 3.2beta (c) 2004, The University of
+    Oxford (the "Software")
+    
+    The Software remains the property of the University of Oxford ("the
+    University").
+    
+    The Software is distributed "AS IS" under this Licence solely for
+    non-commercial use in the hope that it will be useful, but in order
+    that the University as a charitable foundation protects its assets for
+    the benefit of its educational and research purposes, the University
+    makes clear that no condition is made or to be implied, nor is any
+    warranty given or to be implied, as to the accuracy of the Software,
+    or that it will be suitable for any particular purpose or for use
+    under any specific conditions. Furthermore, the University disclaims
+    all responsibility for the use which is made of the Software. It
+    further disclaims any liability for the outcomes arising from using
+    the Software.
+    
+    The Licensee agrees to indemnify the University and hold the
+    University harmless from and against any and all claims, damages and
+    liabilities asserted by third parties (including claims for
+    negligence) which arise directly or indirectly from the use of the
+    Software or the sale of any products based on the Software.
+    
+    No part of the Software may be reproduced, modified, transmitted or
+    transferred in any form or by any means, electronic or mechanical,
+    without the express permission of the University. The permission of
+    the University is not required if the said reproduction, modification,
+    transmission or transference is done without financial return, the
+    conditions of this Licence are imposed upon the receiver of the
+    product, and all original and amended source code is included in any
+    transmitted product. You may be held legally responsible for any
+    copyright infringement that is caused or encouraged by your failure to
+    abide by these terms and conditions.
+    
+    You are not permitted under this Licence to use this Software
+    commercially. Use for which any financial return is received shall be
+    defined as commercial use, and includes (1) integration of all or part
+    of the source code or the Software into a product for sale or license
+    by or on behalf of Licensee to third parties or (2) use of the
+    Software or any derivative of it for research with the final aim of
+    developing software products for sale or license to a third party or
+    (3) use of the Software or any derivative of it for research with the
+    final aim of developing non-software products for sale or license to a
+    third party, or (4) use of the Software to provide any service to an
+    external organisation for which payment is received. If you are
+    interested in using the Software commercially, please contact Isis
+    Innovation Limited ("Isis"), the technology transfer company of the
+    University, to negotiate a licence. Contact details are:
+    innovation@isis.ox.ac.uk quoting reference DE/1112. */
 
 #ifndef __TRACTVOLS_H_
 #define __TRACTVOLS_H_
@@ -16,93 +76,36 @@
 #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:
-      probtrackxOptions& opts;
-      vector<volume4D<float>* > thsamples;
-      vector<volume4D<float>* > phsamples;
-      vector<volume4D<float>* > fsamples;
-      bool init_sample;
-      int fibst;
+      volume4D<float> thsamples;
+      volume4D<float> phsamples;
+      volume4D<float> fsamples;
       bool usef;
     public:
       //constructors::
-      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];
-      }
+      TractVols(const bool& usefin=false):usef(usefin){}
+      ~TractVols(){}
       
-      
-      void reset(const int& fibst_in){
-	init_sample=true;
-	fibst=fibst_in;
-      }
       //Initialise
       void initialise(const string& basename){
-	
-
-	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;
+	read_volume4D(thsamples,basename+"_thsamples");
+	read_volume4D(phsamples,basename+"_phsamples");
+	if(usef)
+	  read_volume4D(fsamples,basename+"_fsamples");
       }
       
-      
-      ColumnVector sample(const float& x,const float& y,const float&z,const float& r_x,const float& r_y,const float& r_z){
+
+      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);
 	
 	////////Probabilistic interpolation
 	int cx =(int) ceil(x),fx=(int) floor(x);
@@ -150,49 +153,17 @@ namespace TRACTVOLS{
 	
 	
 	float samp=rand(); samp/=RAND_MAX;
-	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;
-		}
+	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);
 		
-	      }
-	      
-	    }
-	    
-	    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 phi = phsamples(int(newx),int(newy),int(newz),int(samp));
+	float theta = thsamples(int(newx),int(newy),int(newz),int(samp));
 	
 	float f;
 	
 	if(usef){
-	  f = (*fsamples[fibind])(int(newx),int(newy),int(newz),int(samp));
+	  f = fsamples(int(newx),int(newy),int(newz),int(samp));
 	}
 	else
 	  f=1;
@@ -205,7 +176,7 @@ namespace TRACTVOLS{
 
       ColumnVector dimensions(){
 	ColumnVector dims(3);
-	dims << (*thsamples[0]).xdim() <<(*thsamples[0]).ydim() << (*thsamples[0]).zdim();
+	dims << thsamples.xdim() <<thsamples.ydim() << thsamples.zdim();
 	return dims;
       }
     };

tractvolsx.h

+/*  tractvolsx.h
+
+    Tim Behrens, FMRIB Image Analysis Group
+
+    Copyright (C) 2004 University of Oxford  */
+
+/*  CCOPYRIGHT  */
+
+#ifndef __TRACTVOLSX_H_
+#define __TRACTVOLSX_H_
+
+/////////////////////////////////////////////////////////
+//         Class TractVolsx                             //
+/////////////////////////////////////////////////////////
+
+#include "newimage/newimageall.h"
+#include <iostream>
+#include "stdlib.h"
+#include "probtrackxOptions.h"
+using namespace std;
+using namespace NEWIMAGE;
+using namespace TRACT;
+
+namespace TRACTVOLSX{
+  class Tractvolsx
+    {
+    private:
+      probtrackxOptions& opts;
+      vector<volume4D<float>* > thsamples;
+      vector<volume4D<float>* > phsamples;
+      vector<volume4D<float>* > fsamples;
+      bool init_sample;
+      int fibst;
+      bool usef;
+    public:
+      //constructors::
+      Tractvolsx(const bool& usefin=false):opts(probtrackxOptions::getInstance()),init_sample(true),fibst(1),usef(usefin){}
+      Tractvolsx():opts(probtrackxOptions::getInstance()){}
+      ~Tractvolsx(){
+	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){
+	
+
+	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,const float& r_x,const float& r_y,const float& r_z){
+	
+	////////Probabilistic interpolation
+	int cx =(int) ceil(x),fx=(int) floor(x);
+	int cy =(int) ceil(y),fy=(int) floor(y);
+	int cz =(int) ceil(z),fz=(int) floor(z);
+	
+	//cerr<<x<<" "<<y<<" "<<z<<" "<<cx<<" "<<cy<<" "<<cz<<" "<<fx<<" "<<fy<<" "<<fz<<endl;
+	float pcx,pcy,pcz;
+	if(cx==fx)
+	  pcx=1;
+	else
+	  pcx=(x-fx)/(cx-fx);
+	
+	if(cy==fy)
+	  pcy=1;
+	else
+	  pcy=(y-fy)/(cy-fy);
+	
+	if(cz==fz)
+	  pcz=1;
+	else
+	  pcz=(z-fz)/(cz-fz);
+	
+	///////new xyz values from probabilistic interpolation
+	int newx,newy,newz; 
+	float tmp=rand(); tmp/=RAND_MAX;
+	if(tmp>pcx)
+	  newx=fx;
+	else
+	  newx=cx;
+	
+	tmp=rand(); tmp/=RAND_MAX;
+	if(tmp>pcy)
+	  newy=fy;
+	else
+	  newy=cy;
+	
+	tmp=rand(); tmp/=RAND_MAX;
+	if(tmp>pcz)
+	  newz=fz;
+	else
+	  newz=cz;
+ 
+	ColumnVector th_ph_f(3);	
+	
+	
+	float samp=rand(); samp/=RAND_MAX;
+	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;
+		}
+		
+	      }
+	      
+	    }
+	    
+	    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[fibind])(int(newx),int(newy),int(newz),int(samp));
+	}
+	else
+	  f=1;
+	
+	th_ph_f(1)=theta;
+	th_ph_f(2)=phi;
+	th_ph_f(3)=f;
+	return th_ph_f;
+      }
+
+      ColumnVector dimensions(){
+	ColumnVector dims(3);
+	dims << (*thsamples[0]).xdim() <<(*thsamples[0]).ydim() << (*thsamples[0]).zdim();
+	return dims;
+      }
+    };
+}
+
+#endif
+
+
+