From 8bdecf30f5550a04c60935a0064e014f77d135e5 Mon Sep 17 00:00:00 2001
From: Saad Jbabdi <saad@fmrib.ox.ac.uk>
Date: Tue, 21 Aug 2007 17:00:38 +0000
Subject: [PATCH] meshes work now for probtrackx, and streamlines go in one
direction only
---
ptx_meshmask.cc | 58 +++++++++++++++++++++++++------------------------
streamlines.cc | 25 +++++++++++++++------
streamlines.h | 3 ++-
3 files changed, 50 insertions(+), 36 deletions(-)
diff --git a/ptx_meshmask.cc b/ptx_meshmask.cc
index ade7b6d..2fd96ba 100644
--- a/ptx_meshmask.cc
+++ b/ptx_meshmask.cc
@@ -18,7 +18,6 @@ void meshmask()
{
probtrackxOptions& opts =probtrackxOptions::getInstance();
-
// load seed mesh
Mesh mseeds;
mseeds.load(opts.meshfile.value());
@@ -33,7 +32,6 @@ void meshmask()
seeds=0;
-
////////////////////////////////
// Log& logger = LogSingleton::getInstance();
Streamliner stline;
@@ -41,40 +39,44 @@ void meshmask()
counter.initialise();
Seedmanager seedmanager(counter);
- Matrix Seeds_to_DTI;
- if(opts.seeds_to_dti.value()!=""){
- read_ascii_matrix(Seeds_to_DTI,opts.seeds_to_dti.value());
- }
- else{
- Seeds_to_DTI=Identity(4);
- }
- ColumnVector mni_origin(3),fs_coord_mm(3),xyz_dti;
- mni_origin << 92 << 128 << 37;
+
+ ColumnVector fs_coord_mm(4),xyz_vox,seeddim(3);
+ seeddim << seeds.xdim() << seeds.ydim() << seeds.zdim();
+ ColumnVector dir(3);
+ int keeptotal=0;
+
for(vector<Mpoint*>::iterator i = mseeds._points.begin();i!=mseeds._points.end();i++){
if((*i)->get_value() > 0){
- fs_coord_mm<<(*i)->get_coord().X<<(*i)->get_coord().Y<<(*i)->get_coord().Z;
- xyz_dti=mni_to_imgvox(fs_coord_mm,mni_origin,Seeds_to_DTI,counter.get_streamline().get_vols().dimensions());
- float x=xyz_dti(1);float y=xyz_dti(2);float z=xyz_dti(3);
-
- seeds(round(x),round(y),round(z)) = 1;
- }
- }
+ fs_coord_mm<<(*i)->get_coord().X<<(*i)->get_coord().Y<<(*i)->get_coord().Z << 1.0;
+ xyz_vox = seeds.qform_mat().i()*fs_coord_mm;
+ float x=xyz_vox(1);float y=xyz_vox(2);float z=xyz_vox(3);
- int keeptotal=0;
- for(int z=0;z<seeds.zsize();z++){
- cout <<"sl "<<z<<endl;
- for(int y=0;y<seeds.ysize();y++){
- for(int x=0;x<seeds.xsize();x++){
- if(seeds(x,y,z)>0){
- cout <<"run"<<endl;
- keeptotal += seedmanager.run(x,y,z);
- }
- }
+ seeds(round(x),round(y),round(z)) = 1;
+
+ cout <<"run"<<endl;
+ dir << (*i)->local_normal().X << (*i)->local_normal().Y << (*i)->local_normal().Z;
+ keeptotal += seedmanager.run(round(x),round(y),round(z),true,-1,dir);
+
}
}
+ //return;
+
+// for(int z=0;z<seeds.zsize();z++){
+// cout <<"sl "<<z<<endl;
+// for(int y=0;y<seeds.ysize();y++){
+// for(int x=0;x<seeds.xsize();x++){
+// if(seeds(x,y,z)>0){
+// cout <<"run"<<endl;
+// dir << (*i)->local_normal().X << (*i)->local_normal().Y << (*i)->local_normal().Z;
+// keeptotal += seedmanager.run(x,y,z,true,-1,dir);
+// }
+// }
+// }
+// }
+
counter.save_total(keeptotal);
counter.save();
diff --git a/streamlines.cc b/streamlines.cc
index dc8d685..59f38d9 100644
--- a/streamlines.cc
+++ b/streamlines.cc
@@ -74,7 +74,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){
+ bool Streamliner::streamline(const float& x_init,const float& y_init,const float& z_init, const ColumnVector& dim_seeds,const int& fibst,const ColumnVector& dir){
//fibst tells tractvolsx which fibre to start with if there are more than one..
//x_init etc. are in seed space...
@@ -92,6 +92,10 @@ namespace TRACT{
m_path.clear();
x=xst;y=yst;z=zst;
m_part.change_xyz(x,y,z);
+ if(opts.meshfile.value()!="")
+ 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;
@@ -100,7 +104,9 @@ namespace TRACT{
for(unsigned int pf=0;pf<m_passed_flags.size();pf++) {
m_passed_flags[pf]=false; /// only keep it if this streamline went through all the masks
}
-
+
+ Matrix DTI_to_Seeds(4,4);
+ DTI_to_Seeds = m_Seeds_to_DTI.i();
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) ){
///////////////////////////////////
@@ -126,7 +132,7 @@ namespace TRACT{
x=m_part.x();y=m_part.y();z=m_part.z();
xyz_dti <<x<<y<<z;
- xyz_seeds=vox_to_vox(xyz_dti,vols.dimensions(),dim_seeds,m_Seeds_to_DTI.i());
+ xyz_seeds=vox_to_vox(xyz_dti,vols.dimensions(),dim_seeds,DTI_to_Seeds);
int x_s =(int)round((float)xyz_seeds(1));
int y_s =(int)round((float)xyz_seeds(2));
int z_s =(int)round((float)xyz_seeds(3));
@@ -614,8 +620,13 @@ namespace TRACT{
}
- // this function now returns the total number of pathways that survived a streamlining (SJ)
int Seedmanager::run(const float& x,const float& y,const float& z,bool onewayonly, int fibst){
+ ColumnVector dir(3);
+ dir=0;
+ return run(x,y,z,onewayonly,fibst,dir);
+ }
+ // this function now returns the total number of pathways that survived a streamlining (SJ)
+ int Seedmanager::run(const float& x,const float& y,const float& z,bool onewayonly, int fibst,const ColumnVector& dir){
//onewayonly for mesh things..
cout <<x<<" "<<y<<" "<<z<<endl;
if(fibst == -1){
@@ -628,7 +639,7 @@ namespace TRACT{
//fibst=0;
//else
//fibst=1;// fix this for > 2 fibres
- }
+ }
int nlines=0;
for(int p=0;p<opts.nparticles.value();p++){
@@ -638,7 +649,7 @@ namespace TRACT{
m_stline.reset();
bool forwardflag=false,backwardflag=false;
if(!onewayonly){
- if(m_stline.streamline(x,y,z,m_seeddims,fibst)){ //returns whether to count the streamline or not
+ if(m_stline.streamline(x,y,z,m_seeddims,fibst,dir)){ //returns whether to count the streamline or not
forwardflag=true;
m_counter.store_path();
m_counter.count_streamline();
@@ -646,7 +657,7 @@ namespace TRACT{
}
m_stline.reverse();
}
- if(m_stline.streamline(x,y,z,m_seeddims,fibst)){
+ if(m_stline.streamline(x,y,z,m_seeddims,fibst,dir)){
backwardflag=true;
m_counter.count_streamline();
//nlines++; //count twice ?
diff --git a/streamlines.h b/streamlines.h
index 0558f57..eb74f3e 100644
--- a/streamlines.h
+++ b/streamlines.h
@@ -82,7 +82,7 @@ namespace TRACT{
inline void reverse(){
m_part.restart_reverse();
}
- bool streamline(const float& x_init,const float& y_init, const float& z_init,const ColumnVector& dim_seeds,const int& fibst);
+ bool streamline(const float& x_init,const float& y_init, const float& z_init,const ColumnVector& dim_seeds,const int& fibst,const ColumnVector& dir);
};
@@ -206,6 +206,7 @@ namespace TRACT{
m_seeddims<<m_seeds.xdim()<<m_seeds.ydim()<<m_seeds.zdim();
}
int run(const float& x,const float& y,const float& z,bool onewayonly=false, int fibst=-1);
+ int run(const float& x,const float& y,const float& z,bool onewayonly, int fibst,const ColumnVector& dir);
};
}
--
GitLab