From 968a28a14bc8ffaf37e9f0809f783fdce721871c Mon Sep 17 00:00:00 2001
From: Saad Jbabdi <saad@fmrib.ox.ac.uk>
Date: Tue, 13 Oct 2009 12:52:33 +0000
Subject: [PATCH] added distthresh and matrix3 options
---
probtrackxOptions.h | 24 +++++++---
streamlines.cc | 109 +++++++++++++++++++++++++++++++++++++-------
streamlines.h | 14 +++++-
3 files changed, 123 insertions(+), 24 deletions(-)
diff --git a/probtrackxOptions.h b/probtrackxOptions.h
index c232363..236c2a5 100644
--- a/probtrackxOptions.h
+++ b/probtrackxOptions.h
@@ -39,7 +39,8 @@ class probtrackxOptions {
Option<bool> pathdist;
Option<bool> s2tout;
FmribOption<bool> matrix1out;
- FmribOption<bool> matrix2out;
+ Option<bool> matrix2out;
+ FmribOption<bool> matrix3out;
FmribOption<bool> maskmatrixout;
Option<string> outfile;
Option<string> rubbishfile;
@@ -53,11 +54,12 @@ class probtrackxOptions {
Option<string> waypoints;
Option<bool> network;
Option<string> meshfile;
- FmribOption<string> lrmask;
+ Option<string> lrmask;
Option<string> logdir;
Option<bool> forcedir;
Option<int> nparticles;
Option<int> nsteps;
+ Option<float> distthresh;
Option<float> c_thr;
FmribOption<float> fibthresh;
Option<float> steplength;
@@ -68,7 +70,7 @@ class probtrackxOptions {
Option<bool> modeuler;
Option<int> rseed;
Option<bool> seedcountastext;
- FmribOption<bool> splitmatrix2;
+ Option<bool> splitmatrix2;
void parse_command_line(int argc, char** argv,Log& logger);
void modecheck();
@@ -108,7 +110,7 @@ class probtrackxOptions {
string("Bet binary mask file in diffusion space"),
true, requires_argument),
seedfile(string("-x,--seed"), string("Seed"),
- string("Seed volume, or voxel, or ascii file with multiple volumes"),
+ string("Seed volume, or voxel, or ascii file with multiple volumes, or freesurfer label file"),
true, requires_argument),
mode(string("--mode"), string(""),
string("use --mode=simple for single seed voxel"),
@@ -131,6 +133,9 @@ class probtrackxOptions {
matrix2out(string("--omatrix2"), false,
string("output matrix2"),
false, no_argument),
+ matrix3out(string("--omatrix3"), false,
+ string("output matrix3 (uses the termination mask to produce NxN matrix)"),
+ false, no_argument),
maskmatrixout(string("--omaskmatrix"), false,
string("output maskmatrix"),
false, no_argument),
@@ -173,8 +178,8 @@ class probtrackxOptions {
lrmask(string("--lrmask"), string(""),
string("low resolution binary brain mask for stroring connectivity distribution in matrix2 mode"),
false, requires_argument),
- logdir(string("--dir"), string(""),
- string("Directory to put the final volumes in - code makes this directory"),
+ logdir(string("--dir"), string("logdir"),
+ string("Directory to put the final volumes in - code makes this directory - default='logdir'"),
false, requires_argument),
forcedir(string("--forcedir"), false,
string("Use the actual directory name given - i.e. don't add + to make a new directory"),
@@ -185,6 +190,9 @@ 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)"),
+ false, requires_argument),
c_thr(string("-c,--cthr"), 0.2,
string("Curvature threshold - default=0.2"),
false, requires_argument),
@@ -216,7 +224,7 @@ class probtrackxOptions {
string("Output seed-to-target counts as a text file (useful when seeding from a mesh)"),
false, no_argument),
splitmatrix2(string("--splitmatrix2"), false,
- string("split matrix 2 (in case it is too big)"),
+ string("split matrix 2 along seed dimension (in case it is too large)"),
false, no_argument),
options("probtrackx","probtrackx -s <basename> -m <maskname> -x <seedfile> -o <output> --targetmasks=<textfile>\n probtrackx --help\n")
{
@@ -244,6 +252,7 @@ class probtrackxOptions {
options.add(s2tout);
options.add(matrix1out);
options.add(matrix2out);
+ options.add(matrix3out);
options.add(maskmatrixout);
options.add(outfile);
options.add(rubbishfile);
@@ -253,6 +262,7 @@ class probtrackxOptions {
options.add(dti_to_seeds);
options.add(nparticles);
options.add(nsteps);
+ options.add(distthresh);
options.add(c_thr);
options.add(fibthresh);
options.add(steplength);
diff --git a/streamlines.cc b/streamlines.cc
index 13716c2..e639414 100644
--- a/streamlines.cc
+++ b/streamlines.cc
@@ -155,21 +155,26 @@ namespace TRACT{
// set rotation matrix
Matrix F(3,3),u(3,3),v(3,3);
DiagonalMatrix d(3);
- F << -m_Seeds_to_DTI(1,1) << m_Seeds_to_DTI(1,3) << -m_Seeds_to_DTI(1,2)
- << -m_Seeds_to_DTI(2,1) << m_Seeds_to_DTI(2,3) << -m_Seeds_to_DTI(2,2)
- << -m_Seeds_to_DTI(3,1) << m_Seeds_to_DTI(3,3) << -m_Seeds_to_DTI(3,2);
-
+
+ if(opts.meshfile.value()!=""){
+ F << -m_Seeds_to_DTI(1,1) << m_Seeds_to_DTI(1,3) << -m_Seeds_to_DTI(1,2)
+ << -m_Seeds_to_DTI(2,1) << m_Seeds_to_DTI(2,3) << -m_Seeds_to_DTI(2,2)
+ << -m_Seeds_to_DTI(3,1) << m_Seeds_to_DTI(3,3) << -m_Seeds_to_DTI(3,2);
+ }
+ else{
+ F = m_Seeds_to_DTI.SubMatrix(1,3,1,3);
+ }
SVD(F*F.t(),d,u,v);
m_rotdir.ReSize(3,3);
m_rotdir = (u*sqrt(d)*v.t()).i()*F;
-
+
}
else{
m_IsNonlinXfm = true;
FnirtFileReader ffr(opts.seeds_to_dti.value());
m_Seeds_to_DTI_warp = ffr.FieldAsNewimageVolume4D(true);
if(opts.dti_to_seeds.value()==""){
- cerr << "Error: DTI -> Seeds transform needed" << endl;
+ cerr << "TRACT::Streamliner:: DTI -> Seeds transform needed" << endl;
exit(1);
}
FnirtFileReader iffr(opts.dti_to_seeds.value());
@@ -236,7 +241,7 @@ namespace TRACT{
m_passed_flags[pf]=false; /// only keep it if this streamline went through all the masks
}
-
+ float pathlength=0;
for( int it = 1 ; it <= opts.nsteps.value()/2; it++){
if( (m_mask( round(m_part.x()), round(m_part.y()), round(m_part.z())) > 0) ){
@@ -263,6 +268,7 @@ namespace TRACT{
x=m_part.x();y=m_part.y();z=m_part.z();
xyz_dti <<x<<y<<z;
+
// now find xyz in seeds space
if(!m_IsNonlinXfm)
xyz_seeds = vox_to_vox(xyz_dti,vols.dimensions(),dim_seeds,m_DTI_to_Seeds);
@@ -291,6 +297,7 @@ namespace TRACT{
m_path.push_back(xyz_seeds);
+ //m_path.push_back(xyz_dti);
partlength++;
@@ -344,8 +351,9 @@ namespace TRACT{
}
}
-
-
+
+ // increase path length
+ pathlength += opts.steplength.value();
}
@@ -362,6 +370,8 @@ namespace TRACT{
}
if(rubbish_passed)
accept_path=false;
+ if(pathlength<opts.distthresh.value())
+ accept_path=false;
return accept_path;
}
@@ -405,6 +415,9 @@ namespace TRACT{
if(opts.matrix2out.value()){
initialise_matrix2();
}
+ if(opts.matrix3out.value()){
+ initialise_matrix3();
+ }
if(opts.maskmatrixout.value()){
initialise_maskmatrix();
}
@@ -606,6 +619,40 @@ 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.
+ 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++;
+ }
+ OUT(nstop);
+
+ m_CoordMat3.reinitialize(nstop,3,1);
+ m_ConMat3.reinitialize(nstop,nstop,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++;
+ }
+ save_volume(m_CoordMat3,logger.appendDir("coords_for_fdt_matrix3"));
+ exit(1);
+ }
void Counter::count_streamline(){
if(opts.simpleout.value()||opts.matrix1out.value()){
update_pathdist();
@@ -798,8 +845,23 @@ 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))));
+ m_ConMat3(row1,row2,0) += 1;
+ m_ConMat3(row2,row1,0) += 1;
+ }
+ }
void Counter::reset_beenhere2(const bool& forwardflag,const bool& backwardflag){
if(forwardflag){
@@ -858,6 +920,9 @@ namespace TRACT{
if(opts.matrix2out.value()){
save_matrix2();
}
+ if(opts.matrix3out.value()){
+ save_matrix3();
+ }
if(opts.maskmatrixout.value()){
save_maskmatrix();
}
@@ -983,7 +1048,11 @@ namespace TRACT{
}
}
-
+void Counter::save_matrix3(){
+ save_volume(m_ConMat3,logger.appendDir("fdt_matrix3"));
+ applycoordchange(m_CoordMat3, m_seeds.niftivox2newimagevox_mat().i());
+ save_volume(m_CoordMat3,logger.appendDir("coords_for_fdt_matrix3"));
+}
int Seedmanager::run(const float& x,const float& y,const float& z,bool onewayonly, int fibst){
ColumnVector dir(3);
dir=0;
@@ -1008,18 +1077,22 @@ namespace TRACT{
// now re-orient dir using xfm transform
ColumnVector rotdir(3);
- //m_stline.rotdir(dir,rotdir,x,y,z);
- rotdir=dir;
+ m_stline.rotdir(dir,rotdir,x,y,z);
+ //rotdir=dir;
int nlines=0;
for(int p=0;p<opts.nparticles.value();p++){
+ if(opts.randfib.value()){
+ 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();
bool forwardflag=false,backwardflag=false;
bool counted=false;
- if(!onewayonly){
+ if(!onewayonly || opts.matrix3out.value()){//always go both ways in matrix3 mode
if(m_stline.streamline(x,y,z,m_seeddims,fibst,rotdir)){ //returns whether to count the streamline or not
forwardflag=true;
m_counter.store_path();
@@ -1033,7 +1106,10 @@ namespace TRACT{
backwardflag=true;
m_counter.count_streamline();
-
+ if(counted && opts.matrix3out.value()){//if the first half is not counted, do not bother with matrix3
+ m_counter.update_matrix3();
+ }
+
if(!counted)nlines++; // the other half has is counted here
}
@@ -1043,6 +1119,7 @@ namespace TRACT{
m_counter.count_seed();
+
return nlines;
}
diff --git a/streamlines.h b/streamlines.h
index d6ca67b..031f624 100644
--- a/streamlines.h
+++ b/streamlines.h
@@ -109,7 +109,10 @@ namespace TRACT{
void rotdir(const ColumnVector& dir,ColumnVector& rotdir,const float& x,const float& y,const float& z);
+ const volume<int>& get_stop()const{return m_stop;}
+ // debug
+ //const volume<int>& get_mask_ref()const{return m_mask;}
};
@@ -148,6 +151,10 @@ namespace TRACT{
volume<int> m_beenhere2;
int m_Conrow2;
ColumnVector m_lrdim;
+
+ volume<int> m_ConMat3;
+ volume<int> m_Lookup3;
+ volume<int> m_CoordMat3;
const volume<float>& m_seeds;
ColumnVector m_seedsdim;
@@ -172,13 +179,15 @@ namespace TRACT{
void initialise_path_dist(){
m_prob.reinitialize(m_seeds.xsize(),m_seeds.ysize(),m_seeds.zsize());
copybasicproperties(m_seeds,m_prob);
+ //m_prob.reinitialize(m_stline.get_mask_ref().xsize(),m_stline.get_mask_ref().ysize(),m_stline.get_mask_ref().zsize());
+ //copybasicproperties(m_stline.get_mask_ref(),m_prob);
m_prob=0;
}
void initialise_seedcounts();
void initialise_matrix1(); //Need to make sure that initialise_path_dist is run first
-
void initialise_matrix2();
+ void initialise_matrix3();
void initialise_maskmatrix(){} //not written yet
@@ -204,6 +213,8 @@ namespace TRACT{
void update_matrix2_row(); //but run this one every streamline as with the others
void next_matrix2_row(){m_Conrow2++;}//and then run this after each voxel..
void reset_beenhere2(const bool& forwardflag,const bool& backwardflag);
+
+ void update_matrix3();
void update_maskmatrix(){} //not written yet
@@ -215,6 +226,7 @@ namespace TRACT{
void save_seedcounts();
void save_matrix1();
void save_matrix2();
+ void save_matrix3();
void save_maskmatrix(){}//not written yet
--
GitLab