From a7f5560431ae278bc70931ab391e548db3cafd06 Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Thu, 14 Jul 2005 10:06:11 +0000
Subject: [PATCH] *** empty log message ***
---
Makefile | 5 ++
particle.cc | 9 +---
particle.h | 47 +++++++++++++----
probtrackx.cc | 2 +-
probtrackx.h | 2 +-
ptx_seedmask.cc | 2 +-
ptx_seedmask.h | 2 +-
ptx_simple.cc | 2 +-
ptx_simple.h | 2 +-
ptx_twomasks.cc | 2 +-
ptx_twomasks.h | 2 +-
streamlines.cc | 2 +-
streamlines.h | 12 ++---
tractvols.h | 135 +++++++++++++++++++++++++++++++++++++++---------
14 files changed, 171 insertions(+), 55 deletions(-)
diff --git a/Makefile b/Makefile
index 346add5..37c61de 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 cdecaf5..b7c3f67 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 f30696d..6ea2118 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 dfffb8c..4c9cf19 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 856a6d0..1e2fb56 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 e755cf6..9539738 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 9fc135a..dbc4968 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 3425efb..fd47364 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 3f8590e..eb57c1f 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 1abcde0..6010297 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 626d000..40fce4a 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 c32eab2..aa03ccc 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 d2869a6..c7bf349 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 4eb9da7..b7a9b27 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;
}
};
--
GitLab