diff --git a/probtrackxOptions.h b/probtrackxOptions.h
index 236c2a5831f1be338ea491c87407d0951e3069f1..d11a53cab2b2bfedfbed0a96c68ededea80c8c96 100644
--- a/probtrackxOptions.h
+++ b/probtrackxOptions.h
@@ -41,6 +41,7 @@ class probtrackxOptions {
FmribOption<bool> matrix1out;
Option<bool> matrix2out;
FmribOption<bool> matrix3out;
+ FmribOption<string> maskmatrix3;
FmribOption<bool> maskmatrixout;
Option<string> outfile;
Option<string> rubbishfile;
@@ -134,8 +135,11 @@ class probtrackxOptions {
string("output matrix2"),
false, no_argument),
matrix3out(string("--omatrix3"), false,
- string("output matrix3 (uses the termination mask to produce NxN matrix)"),
+ string("output matrix3 (NxN connectivity matrix)"),
false, no_argument),
+ maskmatrix3(string("--mask3"), "",
+ string("mask used for NxN connectivity matrix"),
+ false, requires_argument),
maskmatrixout(string("--omaskmatrix"), false,
string("output maskmatrix"),
false, no_argument),
@@ -190,8 +194,8 @@ class probtrackxOptions {
nsteps(string("-S,--nsteps"), 2000,
string("Number of steps per sample - default=2000"),
false, requires_argument),
- distthresh(string("--distthresh"), 1000,
- string("Discards samples shorter than this threshold (in mm - default=1000)"),
+ distthresh(string("--distthresh"), 0,
+ string("Discards samples shorter than this threshold (in mm - default=0)"),
false, requires_argument),
c_thr(string("-c,--cthr"), 0.2,
string("Curvature threshold - default=0.2"),
@@ -253,6 +257,7 @@ class probtrackxOptions {
options.add(matrix1out);
options.add(matrix2out);
options.add(matrix3out);
+ options.add(maskmatrix3);
options.add(maskmatrixout);
options.add(outfile);
options.add(rubbishfile);
diff --git a/streamlines.cc b/streamlines.cc
index 1c5ebca57b8f04d48b0e17fcd88d96bca705ac4f..6eb14f38074488f4aef9b8f59d7ef8860af30796 100644
--- a/streamlines.cc
+++ b/streamlines.cc
@@ -195,7 +195,6 @@ namespace TRACT{
m_y_s_init=0;
m_z_s_init=0;
-
}
@@ -230,8 +229,6 @@ namespace TRACT{
m_part.set_dir(dir(1),dir(2),dir(3));//Set the start dir so that we track inwards from cortex
}
-
- int partlength=0;
bool rubbish_passed=false;
bool stop_flag=false;
//bool has_goneout=false;
@@ -300,6 +297,17 @@ namespace TRACT{
//m_path.push_back(xyz_dti);
pathlength += opts.steplength.value();
+ // //////////////////////////////
+ // update coordinates for matrix3
+ if(opts.matrix3out.value()){
+ if( m_mask3(x_s,y_s,z_s)!=0 ){
+ if( m_beenhere3(x_s,y_s,z_s)==0 ){
+ m_beenhere3(x_s,y_s,z_s)=1;
+ m_inmask3.push_back(xyz_seeds);
+ }
+ }
+ }
+ // //////////////////////////////
if(opts.rubbishfile.value()!=""){
@@ -324,7 +332,7 @@ namespace TRACT{
th_ph_f=vols.sample(m_part.x(),m_part.y(),m_part.z(),m_part.rx(),m_part.ry(),m_part.rz(),pref_x,pref_y,pref_z);
}
-
+
tmp2=rand(); tmp2/=RAND_MAX;
if(th_ph_f(3)>tmp2){
if(!m_part.check_dir(th_ph_f(1),th_ph_f(2),opts.c_thr.value())){
@@ -617,39 +625,43 @@ namespace TRACT{
}
- // the following will use the termination mask voxels to initialise
- // and NxN matrix. This matrix will store the number of samples from
- // each seed voxel that have made it to the termination mask from one
- // side and from the other.
+ // the following will use the voxels of a given mask to initialise
+ // an NxN matrix. This matrix will store the number of samples from
+ // each seed voxel that have made it to each pair of voxels in the mask
void Counter::initialise_matrix3(){
- const volume<int>& stop = m_stline.get_stop();
- int nstop=0;
- for(int Wz=stop.minz();Wz<=stop.maxz();Wz++)
- for(int Wy=stop.miny();Wy<=stop.maxy();Wy++)
- for(int Wx=stop.minx();Wx<=stop.maxx();Wx++){
- if(stop(Wx,Wy,Wz)==0)continue;
- nstop++;
+ volume<int>& m_mask3 = m_stline.get_mask3();
+ volume<int>& m_beenhere3 = m_stline.get_beenhere3();
+ vector<ColumnVector> m_inmask3 = m_stline.get_inmask3();
+
+ read_volume(m_mask3,opts.maskmatrix3.value());
+ m_beenhere3.reinitialize(m_mask3.xsize(),m_mask3.ysize(),m_mask3.zsize());
+ m_beenhere3=0;
+
+ int nmask3=0;
+ for(int Wz=m_mask3.minz();Wz<=m_mask3.maxz();Wz++)
+ for(int Wy=m_mask3.miny();Wy<=m_mask3.maxy();Wy++)
+ for(int Wx=m_mask3.minx();Wx<=m_mask3.maxx();Wx++){
+ if(m_mask3(Wx,Wy,Wz)==0)continue;
+ nmask3++;
}
- OUT(nstop);
- m_CoordMat3.reinitialize(nstop,3,1);
- m_ConMat3.reinitialize(nstop,nstop,1);
+ m_CoordMat3.reinitialize(nmask3,3,1);
+ m_ConMat3.reinitialize(nmask3,nmask3,1);
m_ConMat3=0;
- m_Lookup3.reinitialize(stop.xsize(),stop.ysize(),stop.zsize());
-
- nstop=0;
- for(int Wz=stop.minz();Wz<=stop.maxz();Wz++)
- for(int Wy=stop.miny();Wy<=stop.maxy();Wy++)
- for(int Wx=stop.minx();Wx<=stop.maxx();Wx++){
- if(stop(Wx,Wy,Wz)==0)continue;
- m_CoordMat3(nstop,0,0) = Wx;
- m_CoordMat3(nstop,1,0) = Wy;
- m_CoordMat3(nstop,2,0) = Wz;
- m_Lookup3(Wx,Wy,Wz) = nstop;
- nstop++;
+ m_Lookup3.reinitialize(m_mask3.xsize(),m_mask3.ysize(),m_mask3.zsize());
+
+ nmask3=0;
+ for(int Wz=m_mask3.minz();Wz<=m_mask3.maxz();Wz++)
+ for(int Wy=m_mask3.miny();Wy<=m_mask3.maxy();Wy++)
+ for(int Wx=m_mask3.minx();Wx<=m_mask3.maxx();Wx++){
+ if(m_mask3(Wx,Wy,Wz)==0)continue;
+ m_CoordMat3(nmask3,0,0) = Wx;
+ m_CoordMat3(nmask3,1,0) = Wy;
+ m_CoordMat3(nmask3,2,0) = Wz;
+ m_Lookup3(Wx,Wy,Wz) = nmask3;
+ nmask3++;
}
save_volume(m_CoordMat3,logger.appendDir("coords_for_fdt_matrix3"));
- exit(1);
}
void Counter::count_streamline(){
if(opts.simpleout.value()||opts.matrix1out.value()){
@@ -845,20 +857,26 @@ namespace TRACT{
}
void Counter::update_matrix3(){
- ColumnVector endpoint1(3),endpoint2(3);
- const vector<ColumnVector>& path=m_stline.get_path_ref();
- const volume<int>& stop = m_stline.get_stop();
- if(m_path.size()==0 || path.size()==0)return;
- endpoint1 = m_path.back();
- endpoint2 = path.back();
-
- if( stop(round(float(endpoint1(1))),round(float(endpoint1(2))),round(float(endpoint1(3))))!=0 )
- if( stop(round(float(endpoint2(1))),round(float(endpoint2(2))),round(float(endpoint2(3))))!=0 ){
- int row1 = m_Lookup3(round(float(endpoint1(1))),round(float(endpoint1(2))),round(float(endpoint1(3))));
- int row2 = m_Lookup3(round(float(endpoint2(1))),round(float(endpoint2(2))),round(float(endpoint2(3))));
+ volume<int> m_beenhere3 = m_stline.get_beenhere3();
+ vector<ColumnVector> m_inmask3 = m_stline.get_inmask3();
+ if(m_inmask3.size()==0)return;
+
+ for(unsigned int i=0;i<m_inmask3.size();i++){
+ for(unsigned int j=i+1;j<m_inmask3.size();j++){
+ int row1 = m_Lookup3((int)round(float(m_inmask3[i](1))),(int)round(float(m_inmask3[i](2))),(int)round(float(m_inmask3[i](3))));
+ int row2 = m_Lookup3((int)round(float(m_inmask3[j](1))),(int)round(float(m_inmask3[j](2))),(int)round(float(m_inmask3[j](3))));
m_ConMat3(row1,row2,0) += 1;
m_ConMat3(row2,row1,0) += 1;
}
+ }
+
+ // cleanup
+ for(unsigned int i=0;i<m_inmask3.size();i++){
+ m_beenhere3((int)round(float(m_inmask3[i](1))),
+ (int)round(float(m_inmask3[i](2))),
+ (int)round(float(m_inmask3[i](3))))=0;
+ }
+ m_inmask3.clear();
}
void Counter::reset_beenhere2(const bool& forwardflag,const bool& backwardflag){
@@ -1104,7 +1122,7 @@ void Counter::save_matrix3(){
backwardflag=true;
m_counter.count_streamline();
- if(counted && opts.matrix3out.value()){//if the first half is not counted, do not bother with matrix3
+ if(opts.matrix3out.value()){
m_counter.update_matrix3();
}
diff --git a/streamlines.h b/streamlines.h
index 031f624e2de9ad0de2b18071005070e4f84b0cae..9e5cb183745d2b3cb0684909ebf4d57c3e433bb2 100644
--- a/streamlines.h
+++ b/streamlines.h
@@ -59,6 +59,11 @@ namespace TRACT{
float m_y_s_init;
float m_z_s_init;
+ // Streamliner needs to know about matrix3
+ volume<int> m_mask3;
+ volume<int> m_beenhere3;
+ vector<ColumnVector> m_inmask3;
+
// we need this class to know about seed space
const volume<float>& m_seeds;
@@ -113,6 +118,12 @@ namespace TRACT{
// debug
//const volume<int>& get_mask_ref()const{return m_mask;}
+
+ // streamliner needs to be able to tell Counter about matrix3 stuff
+ volume<int>& get_beenhere3(){return m_beenhere3;}
+ volume<int>& get_mask3() {return m_mask3;}
+ vector<ColumnVector>& get_inmask3() {return m_inmask3;}
+
};
@@ -152,13 +163,13 @@ namespace TRACT{
int m_Conrow2;
ColumnVector m_lrdim;
- volume<int> m_ConMat3;
+ volume<int> m_ConMat3;
volume<int> m_Lookup3;
volume<int> m_CoordMat3;
const volume<float>& m_seeds;
ColumnVector m_seedsdim;
- const Streamliner& m_stline;
+ Streamliner& m_stline;
Streamliner& m_nonconst_stline;
public:
diff --git a/tractvolsx.h b/tractvolsx.h
index 624c17fcde32754d691d4d166f4a72d1bdd4a2cc..29203e204c20eb22358cff9e71cd5ec06096bb8b 100644
--- a/tractvolsx.h
+++ b/tractvolsx.h
@@ -43,7 +43,7 @@ namespace TRACTVOLSX{
for(unsigned int m=0;m<phsamples.size();m++)
delete phsamples[m];
for(unsigned int m=0;m<fsamples.size();m++)
- delete fsamples[m];
+ delete fsamples[m];
}
inline int nfibres()const{return (int)thsamples.size();}
@@ -68,10 +68,8 @@ namespace TRACTVOLSX{
cout<<"4"<<endl;
phsamples.push_back(tmpphptr);
cout<<"5"<<endl;
- if(usef){
- read_volume4D(*tmpfptr,basename+"_fsamples");
- fsamples.push_back(tmpfptr);
- }
+ read_volume4D(*tmpfptr,basename+"_fsamples");
+ fsamples.push_back(tmpfptr);
cout<<"6"<<endl;
}
else{
@@ -130,26 +128,26 @@ namespace TRACTVOLSX{
///////new xyz values from probabilistic interpolation
int newx,newy,newz;
- float tmp=rand(); tmp/=RAND_MAX;
+ float tmp=rand(); tmp/=float(RAND_MAX);
if(tmp>pcx)
newx=fx;
else
newx=cx;
- tmp=rand(); tmp/=RAND_MAX;
+ tmp=rand(); tmp/=float(RAND_MAX);
if(tmp>pcy)
newy=fy;
else
newy=cy;
- tmp=rand(); tmp/=RAND_MAX;
+ tmp=rand(); tmp/=float(RAND_MAX);
if(tmp>pcz)
newz=fz;
else
newz=cz;
ColumnVector th_ph_f(3);
- float samp=rand(); samp/=RAND_MAX;
+ float samp=rand(); samp/=float(RAND_MAX);
samp=round(samp*((*thsamples[0]).tsize()-1));
float theta=0,phi=0;
float dotmax=0,dottmp=0;
@@ -171,7 +169,7 @@ namespace TRACTVOLSX{
fibst=0;
}
else{
- float rtmp=rand()/RAND_MAX * float(fibvec.size()-1);
+ float rtmp=rand()/float(RAND_MAX) * float(fibvec.size()-1);
fibst = fibvec[ (int)round(rtmp) ];
}
@@ -194,7 +192,7 @@ namespace TRACTVOLSX{
else{
float fsumtmp2=0;
int fib=0;
- float rtmp=rand()/RAND_MAX;
+ float rtmp=rand()/float(RAND_MAX);
while( fsumtmp2<rtmp){
float ft=(*fsamples[fib])(int(newx),int(newy),int(newz),int(samp));
@@ -207,10 +205,11 @@ namespace TRACTVOLSX{
}
- 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;
+ 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{
if((fabs(prefer_x)+fabs(prefer_y)+fabs(prefer_z))==0){