From 83e34a10ab40c5681b0ccd48dfc21a115296d6e5 Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Thu, 14 Jul 2005 09:31:35 +0000
Subject: [PATCH] *** empty log message ***
---
fibre.h | 19 +-
probtrackx.cc | 81 +++++++
probtrackx.h | 10 +
ptx_seedmask.cc | 45 ++++
ptx_seedmask.h | 16 ++
ptx_simple.cc | 51 +++++
ptx_simple.h | 15 ++
ptx_twomasks.cc | 61 ++++++
ptx_twomasks.h | 13 ++
streamlines.cc | 572 ++++++++++++++++++++++++++++++++++++++++++++++++
streamlines.h | 201 +++++++++++++++++
11 files changed, 1074 insertions(+), 10 deletions(-)
create mode 100644 probtrackx.cc
create mode 100644 probtrackx.h
create mode 100644 ptx_seedmask.cc
create mode 100644 ptx_seedmask.h
create mode 100644 ptx_simple.cc
create mode 100644 ptx_simple.h
create mode 100644 ptx_twomasks.cc
create mode 100644 ptx_twomasks.h
create mode 100644 streamlines.cc
create mode 100644 streamlines.h
diff --git a/fibre.h b/fibre.h
index 6109d12..87c98a6 100644
--- a/fibre.h
+++ b/fibre.h
@@ -108,7 +108,6 @@ namespace FIBRE{
m_th(th), m_ph(ph), m_f(f), m_lam(lam),m_lam_jump(lam_jump), m_d(d),
m_alpha(alpha), m_beta(beta), m_bvals(bvals)
{
-
m_th_old=m_th;
m_ph_old=m_ph;
m_f_old=m_f;
@@ -117,7 +116,7 @@ namespace FIBRE{
m_ph_prop=0.2;
m_f_prop=0.2;
m_lam_prop=1;
-
+
m_th_prior=0;
compute_th_prior();
@@ -313,12 +312,12 @@ namespace FIBRE{
inline bool propose_lam(){
if(m_lam_jump){
- m_lam_old=m_lam;
- m_lam+=normrnd().AsScalar()*m_lam_prop;
- bool rejflag=compute_lam_prior();
- compute_f_prior();
- compute_prior();
- return rejflag;
+ m_lam_old=m_lam;
+ m_lam+=normrnd().AsScalar()*m_lam_prop;
+ bool rejflag=compute_lam_prior();
+ compute_f_prior();
+ compute_prior();
+ return rejflag;
}
else {return true;}
};
@@ -339,7 +338,7 @@ namespace FIBRE{
m_th=rhs.m_th;
m_ph=rhs.m_ph;
m_f=rhs.m_f;
- m_lam=rhs. m_lam;
+ m_lam=rhs.m_lam;
m_th_prop=rhs. m_th_prop;
m_ph_prop=rhs.m_ph_prop;
m_f_prop=rhs.m_f_prop;
@@ -693,7 +692,7 @@ namespace FIBRE{
m_fibres[f].reject_f();
}
-
+
if(!m_fibres[f].propose_lam()){
compute_prior();
compute_energy();
diff --git a/probtrackx.cc b/probtrackx.cc
new file mode 100644
index 0000000..dfffb8c
--- /dev/null
+++ b/probtrackx.cc
@@ -0,0 +1,81 @@
+/* Copyright (C) 2004 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <iostream>
+#include <fstream>
+#include "newimage/newimageall.h"
+#include "utils/log.h"
+#include "meshclass/meshclass.h"
+#include "probtrackx.h"
+
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace TRACT;
+using namespace Utilities;
+using namespace PARTICLE;
+using namespace TRACTVOLS;
+using namespace mesh;
+//using namespace NEWMAT;
+//////////////////////////
+/////////////////////////
+
+
+
+int main ( int argc, char **argv ){
+ probtrackOptions& opts =probtrackOptions::getInstance();
+ Log& logger = LogSingleton::getInstance();
+ opts.parse_command_line(argc,argv,logger);
+ srand(opts.rseed.value());
+
+
+ if(opts.verbose.value()>0){
+ opts.status();
+ }
+ if(opts.mode.value()=="simple")
+ track();
+ else if(fsl_imageexists(opts.mask2.value())){ twomasks();}
+ else if(fsl_imageexists(opts.seedfile.value())){ seedmask();}
+
+ // else if(opts.mode.value()=="seeds_to_targets")
+ // seeds_to_targets();
+ // else if(opts.mode.value()=="seedmask")
+ // alltracts();
+ // else if(opts.mode.value()=="twomasks_symm")
+ // twomasks_symm();
+ // else if(opts.mode.value()=="waypoints")
+ // waypoints();
+ // else if(opts.mode.value()=="matrix1")
+ // matrix1();
+ // else if(opts.mode.value()=="matrix2"){
+ // if(opts.meshfile.value()=="")
+ // matrix2();
+ // else
+ // mesh_matrix2();
+ // }
+ // else if(opts.mode.value()=="maskmatrix")
+ // maskmatrix();
+ // else if(opts.mode.value()=="meshlengths")
+ // mesh_lengths();
+ else{
+ cout <<"Invalid setting for option mode -- try setting mode=help"<<endl;
+ }
+
+ return 0;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/probtrackx.h b/probtrackx.h
new file mode 100644
index 0000000..856a6d0
--- /dev/null
+++ b/probtrackx.h
@@ -0,0 +1,10 @@
+/* Copyright (C) 2004 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "probtrackOptions.h"
+#include "particle.h"
+#include "tractvols.h"
+#include "ptx_simple.h"
+#include "ptx_seedmask.h"
+#include "ptx_twomasks.h"
diff --git a/ptx_seedmask.cc b/ptx_seedmask.cc
new file mode 100644
index 0000000..e755cf6
--- /dev/null
+++ b/ptx_seedmask.cc
@@ -0,0 +1,45 @@
+/* Copyright (C) 2004 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "ptx_seedmask.h"
+#include "streamlines.h"
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace TRACT;
+using namespace Utilities;
+using namespace PARTICLE;
+using namespace TRACTVOLS;
+using namespace mesh;
+
+
+
+void seedmask()
+{
+ probtrackOptions& opts =probtrackOptions::getInstance();
+
+ ////////////////////////////////
+ // Log& logger = LogSingleton::getInstance();
+ volume<int> seeds;
+ read_volume(seeds,opts.seedfile.value());
+ Streamliner stline;
+ Counter counter(seeds,stline);
+ counter.initialise();
+ Seedmanager seedmanager(counter);
+
+ for(int z=0;z<seeds.zsize();z++){
+ for(int y=0;y<seeds.ysize();y++){
+ for(int x=0;x<seeds.xsize();x++){
+ if(seeds(x,y,z)>0){
+ seedmanager.run(x,y,z);
+ }
+ }
+ }
+ }
+
+ counter.save();
+
+}
+
+
diff --git a/ptx_seedmask.h b/ptx_seedmask.h
new file mode 100644
index 0000000..9fc135a
--- /dev/null
+++ b/ptx_seedmask.h
@@ -0,0 +1,16 @@
+/* Copyright (C) 2004 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <fstream>
+#include "newimage/newimageall.h"
+#include "utils/log.h"
+#include "meshclass/meshclass.h"
+#include "probtrackOptions.h"
+#include "particle.h"
+#include "tractvols.h"
+
+
+
+
+void seedmask();
diff --git a/ptx_simple.cc b/ptx_simple.cc
new file mode 100644
index 0000000..3425efb
--- /dev/null
+++ b/ptx_simple.cc
@@ -0,0 +1,51 @@
+/* Copyright (C) 2004 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "ptx_simple.h"
+#include "streamlines.h"
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace TRACT;
+using namespace Utilities;
+using namespace PARTICLE;
+using namespace TRACTVOLS;
+using namespace mesh;
+
+
+void track(){
+ probtrackOptions& opts =probtrackOptions::getInstance();
+
+ ////////////////////////////
+ Log& logger = LogSingleton::getInstance();
+ if(opts.verbose.value()>1){
+ logger.makeDir("particles","particle0",true,false);
+ }
+
+ volume<int> seedref;
+ if(opts.seedref.value()!=""){
+ read_volume(seedref,opts.seedref.value());
+ }
+ else{
+ read_volume(seedref,opts.maskfile.value());
+ }
+
+ Streamliner stline;
+ Counter counter(seedref,stline);
+ counter.initialise();
+ Seedmanager seedmanager(counter);
+
+
+ Matrix Seeds = read_ascii_matrix(opts.seedfile.value());
+ for(int SN=1; SN<=Seeds.Nrows();SN++){
+ float xst=Seeds(SN,1);
+ float yst=Seeds(SN,2);
+ float zst=Seeds(SN,3);
+ seedmanager.run(xst,yst,zst);
+ string add=num2str(Seeds(SN,1))+(string)"_"+num2str(Seeds(SN,2))+(string)"_"+num2str(Seeds(SN,3));
+
+ counter.save_pathdist(add);
+ counter.reset_prob();
+ } //Close Seed number Loop
+}
diff --git a/ptx_simple.h b/ptx_simple.h
new file mode 100644
index 0000000..3f8590e
--- /dev/null
+++ b/ptx_simple.h
@@ -0,0 +1,15 @@
+/* Copyright (C) 2004 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <fstream>
+#include "newimage/newimageall.h"
+#include "utils/log.h"
+#include "meshclass/meshclass.h"
+#include "probtrackOptions.h"
+#include "particle.h"
+#include "tractvols.h"
+
+
+
+void track();
diff --git a/ptx_twomasks.cc b/ptx_twomasks.cc
new file mode 100644
index 0000000..1abcde0
--- /dev/null
+++ b/ptx_twomasks.cc
@@ -0,0 +1,61 @@
+/* Copyright (C) 2004 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "ptx_twomasks.h"
+#include "streamlines.h"
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace TRACT;
+using namespace Utilities;
+using namespace PARTICLE;
+using namespace TRACTVOLS;
+using namespace mesh;
+
+
+
+void twomasks()
+{
+ probtrackOptions& opts =probtrackOptions::getInstance();
+
+ ////////////////////////////////
+ // Log& logger = LogSingleton::getInstance();
+ volume<int> seeds,seeds2;
+ read_volume(seeds,opts.seedfile.value());
+ read_volume(seeds2,opts.mask2.value());
+ Streamliner stline;
+ Counter counter(seeds,stline);
+ counter.initialise();
+ Seedmanager seedmanager(counter);
+
+ stline.add_waymask(seeds2);
+ for(int z=0;z<seeds.zsize();z++){
+ for(int y=0;y<seeds.ysize();y++){
+ for(int x=0;x<seeds.xsize();x++){
+ if(seeds(x,y,z)>0){
+ seedmanager.run(x,y,z);
+ }
+ }
+ }
+ }
+ stline.pop_waymasks();
+
+
+ stline.add_waymask(seeds);
+ for(int z=0;z<seeds2.zsize();z++){
+ for(int y=0;y<seeds2.ysize();y++){
+ for(int x=0;x<seeds2.xsize();x++){
+ if(seeds2(x,y,z)>0){
+ seedmanager.run(x,y,z);
+ }
+ }
+ }
+ }
+
+
+ counter.save();
+
+}
+
+
diff --git a/ptx_twomasks.h b/ptx_twomasks.h
new file mode 100644
index 0000000..626d000
--- /dev/null
+++ b/ptx_twomasks.h
@@ -0,0 +1,13 @@
+/* Copyright (C) 2004 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <fstream>
+#include "newimage/newimageall.h"
+#include "utils/log.h"
+#include "meshclass/meshclass.h"
+#include "probtrackOptions.h"
+#include "particle.h"
+#include "tractvols.h"
+
+void twomasks();
diff --git a/streamlines.cc b/streamlines.cc
new file mode 100644
index 0000000..c32eab2
--- /dev/null
+++ b/streamlines.cc
@@ -0,0 +1,572 @@
+#include "streamlines.h"
+
+
+
+namespace TRACT{
+
+
+
+ void read_masks(vector<string>& masks, const string& filename){
+ ifstream fs(filename.c_str());
+ string tmp;
+ if(fs){
+ fs>>tmp;
+ while(!fs.eof()){
+ masks.push_back(tmp);
+ fs>>tmp;
+ }
+ }
+ else{
+ cerr<<filename<<" does not exist"<<endl;
+ exit(0);
+ }
+ }
+
+
+ Streamliner::Streamliner():opts(probtrackOptions::getInstance()),logger(LogSingleton::getInstance()),
+ vols(opts.usef.value()){
+
+ read_volume(m_mask,opts.maskfile.value());
+ m_part.initialise(0,0,0,0,0,0,opts.steplength.value(),m_mask.xdim(),m_mask.ydim(),m_mask.zdim(),false);
+ if(opts.skipmask.value()!="") read_volume(m_skipmask,opts.skipmask.value());
+ m_lcrat=5;
+ if(opts.loopcheck.value()){
+ m_loopcheck.reinitialize(int(ceil(m_mask.xsize()/m_lcrat)+1),int(ceil(m_mask.ysize()/m_lcrat)+1),int(ceil(m_mask.zsize()/m_lcrat)+1),3);
+ m_loopcheck=0;
+ }
+ if(opts.rubbishfile.value()!=""){
+ read_volume(m_rubbish,opts.rubbishfile.value());
+ }
+
+ vector<string> masknames;
+ if(opts.waypoints.value()!=""){
+ if(fsl_imageexists(opts.waypoints.value())){
+ masknames.push_back( opts.waypoints.value() );
+ }
+ else{
+ read_masks(masknames,opts.waypoints.value());
+ }
+
+ for( unsigned int m = 0; m < masknames.size(); m++ ){
+ volume<int>* tmpptr =new volume<int>;
+ if(opts.verbose.value()>0)
+ cout<<masknames[m]<<endl;
+ read_volume(*tmpptr,masknames[m]);
+ m_waymasks.push_back(tmpptr);
+ m_passed_flags.push_back(false);
+ m_own_waymasks.push_back(true);
+ }
+ }
+ if(opts.seeds_to_dti.value()!=""){
+ read_ascii_matrix(m_Seeds_to_DTI,opts.seeds_to_dti.value());
+ }
+ else{
+ m_Seeds_to_DTI=Identity(4);
+ }
+ vols.initialise(opts.basename.value());
+ m_path.reserve(opts.nparticles.value());
+ m_x_s_init=0;
+ m_y_s_init=0;
+ m_z_s_init=0;
+ }
+
+
+ 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..
+ //x_init etc. are in seed space...
+ vols.reset(fibst);
+ m_x_s_init=x_init; //seed x position in voxels
+ m_y_s_init=y_init; // and y
+ m_z_s_init=z_init; // and z
+ ColumnVector xyz_seeds(3);
+ xyz_seeds<<x_init<<y_init<<z_init;
+ ColumnVector xyz_dti;
+ ColumnVector th_ph_f;
+ float xst,yst,zst,x,y,z,tmp2;
+ xyz_dti=vox_to_vox(xyz_seeds,dim_seeds,vols.dimensions(),m_Seeds_to_DTI);
+ xst=xyz_dti(1);yst=xyz_dti(2);zst=xyz_dti(3);
+ m_path.clear();
+ x=xst;y=yst;z=zst;
+ m_part.change_xyz(x,y,z);
+ int partlength=0;
+ //NB - this only goes in one direction!!
+ 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
+ }
+
+ 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) ){
+
+ ///////////////////////////////////
+ //loopchecking
+ ///////////////////////////////////
+ if(opts.loopcheck.value()){
+ float oldrx=m_loopcheck((int)round(m_part.x()/m_lcrat),(int)round(m_part.y()/m_lcrat),(int)round(m_part.z()/m_lcrat),0);
+ float oldry=m_loopcheck((int)round(m_part.x()/m_lcrat),(int)round(m_part.y()/m_lcrat),(int)round(m_part.z()/m_lcrat),1);
+ float oldrz=m_loopcheck((int)round(m_part.x()/m_lcrat),(int)round(m_part.y()/m_lcrat),(int)round(m_part.z()/m_lcrat),2);
+ if(m_part.rx()*oldrx+m_part.ry()*oldry+m_part.rz()*oldrz<0)
+ {
+ break;
+ }
+
+ m_loopcheck((int)round(m_part.x()/m_lcrat),(int)round(m_part.y()/m_lcrat),(int)round(m_part.z()/m_lcrat),0)=m_part.rx();
+ m_loopcheck((int)round(m_part.x()/m_lcrat),(int)round(m_part.y()/m_lcrat),(int)round(m_part.z()/m_lcrat),1)=m_part.ry();
+ m_loopcheck((int)round(m_part.x()/m_lcrat),(int)round(m_part.y()/m_lcrat),(int)round(m_part.z()/m_lcrat),2)=m_part.rz();
+
+ }
+
+ if(opts.verbose.value()>1)
+ logger<<m_part;
+
+ 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());
+ 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));
+ if(opts.rubbishfile.value()!=""){
+ if(m_rubbish(x_s,y_s,z_s)>0) break;
+ }
+
+ m_path.push_back(xyz_seeds);
+ // m_path(it,1)=x_s;
+ // m_path(it,2)=y_s;
+ // m_path(it,3)=z_s;
+ partlength++;
+
+ //update every passed_flag
+ for( unsigned int wm=0;wm<m_waymasks.size();wm++ ){
+ if( (*m_waymasks[wm])(x_s,y_s,z_s)>0 )
+ m_passed_flags[wm]=true;
+ }
+
+
+ if(opts.skipmask.value() == ""){
+ th_ph_f=vols.sample(m_part.x(),m_part.y(),m_part.z(),m_part.rx(),m_part.ry(),m_part.rz());
+ }
+ else{
+ if(m_skipmask(x_s,y_s,z_s)==0)
+ th_ph_f=vols.sample(m_part.x(),m_part.y(),m_part.z(),m_part.rx(),m_part.ry(),m_part.rz());
+ }
+
+
+ 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())){
+ break;
+ }
+
+ if((th_ph_f(1)!=0&&th_ph_f(2)!=0)){
+ if( (m_mask( round(m_part.x()), round(m_part.y()), round(m_part.z())) != 0) ){
+ if(!opts.modeuler.value())
+ m_part.jump(th_ph_f(1),th_ph_f(2));
+ else
+ {
+ ColumnVector test_th_ph_f;
+ m_part.testjump(th_ph_f(1),th_ph_f(2));
+ test_th_ph_f=vols.sample(m_part.testx(),m_part.testy(),m_part.testz(),m_part.rx(),m_part.ry(),m_part.rz());
+ m_part.jump(test_th_ph_f(1),test_th_ph_f(2));
+ }
+ }
+
+
+ }
+ }
+
+
+ }
+
+ } // Close Step Number Loop
+ if(opts.loopcheck.value()){
+ m_loopcheck=0;
+ }
+ bool accept_path=true;
+ if(m_passed_flags.size()!=0){
+ for(unsigned int i=0; i<m_passed_flags.size();i++)
+ if(!m_passed_flags[i])
+ accept_path=false;
+ }
+ return accept_path;
+ }
+
+
+ void Counter::initialise(){
+ if(opts.simpleout.value()||opts.matrix1out.value()){
+ initialise_path_dist();
+ }
+ if(opts.s2tout.value()){
+ initialise_seedcounts();
+ }
+ if(opts.matrix1out.value()){
+ initialise_matrix1();
+ }
+ if(opts.matrix2out.value()){
+ initialise_matrix2();
+ }
+ if(opts.maskmatrixout.value()){
+ initialise_maskmatrix();
+ }
+ }
+
+ void Counter::initialise_seedcounts(){
+
+ volume<int> tmp;
+ //vector<int> tmpvec;
+ read_masks(m_targetmasknames,opts.targetfile.value());
+ m_targflags.resize(m_targetmasknames.size(),0);
+ //m_particle_numbers.resize(m_targetmasknames.size());
+ //tmpvec.reserve(opts.nparticles.value());
+ cout<<"Number of masks "<<m_targetmasknames.size()<<endl;
+ //are they initialised to zero?
+ for(unsigned int m=0;m<m_targetmasknames.size();m++){
+ read_volume(tmp,m_targetmasknames[m]);
+ m_targetmasks.push_back(tmp);
+ tmp=0;
+ m_seedcounts.push_back(tmp);
+ //m_particle_numbers.push_back(tmpvec);
+ }
+ }
+
+
+ void Counter::initialise_matrix1(){
+ m_Conrow=0;
+ int numseeds=0;
+ for(int Wz=m_seeds.minz();Wz<=m_seeds.maxz();Wz++)
+ for(int Wy=m_seeds.miny();Wy<=m_seeds.maxy();Wy++)
+ for(int Wx=m_seeds.minx();Wx<=m_seeds.maxx();Wx++)
+ if(m_seeds.value(Wx,Wy,Wz)>0)
+ numseeds++;
+
+ m_ConMat.reinitialize(numseeds,numseeds,1);
+ m_CoordMat.reinitialize(numseeds,3,1);
+ int myrow=0;
+
+ for(int Wz=m_seeds.minz();Wz<=m_seeds.maxz();Wz++){
+ for(int Wy=m_seeds.miny();Wy<=m_seeds.maxy();Wy++){
+ for(int Wx=m_seeds.minx();Wx<=m_seeds.maxx();Wx++){
+ if(m_seeds(Wx,Wy,Wz)>0){
+ m_CoordMat(myrow,0,0)=Wx;
+ m_CoordMat(myrow,1,0)=Wy;
+ m_CoordMat(myrow,2,0)=Wz;
+ myrow++;
+ }
+ }
+ }
+ }
+
+ }
+
+ void Counter::initialise_matrix2(){
+ m_Conrow2=0;
+ read_volume(m_lrmask,opts.lrmask.value());
+ m_beenhere2.reinitialize(m_lrmask.xsize(),m_lrmask.ysize(),m_lrmask.zsize());
+ m_lrdim.ReSize(3);
+ m_lrdim<<m_lrmask.xdim()<<m_lrmask.ydim()<<m_lrmask.zdim();
+ int numseeds=0,numnz=0;
+ for(int Wz=m_seeds.minz();Wz<=m_seeds.maxz();Wz++)
+ for(int Wy=m_seeds.miny();Wy<=m_seeds.maxy();Wy++)
+ for(int Wx=m_seeds.minx();Wx<=m_seeds.maxx();Wx++)
+ if(m_seeds.value(Wx,Wy,Wz)>0)
+ numseeds++;
+
+ for(int Wz=m_lrmask.minz();Wz<=m_lrmask.maxz();Wz++)
+ for(int Wy=m_lrmask.miny();Wy<=m_lrmask.maxy();Wy++)
+ for(int Wx=m_lrmask.minx();Wx<=m_lrmask.maxx();Wx++)
+ if(m_lrmask.value(Wx,Wy,Wz)>0)
+ numnz++;
+
+
+ if(numnz> pow(2,(float)sizeof(short)*8-1)-1){
+ cerr<<"Output matrix too big for AVW - stopping."<<endl;
+ cerr<<" Remember - you can store your tracts in "<<endl;
+ cerr<<" low res even if you want your seeds in high res"<<endl;
+ cerr<<" Just subsample the structural space mask"<<endl;
+ cerr<<" Although, it must stay in line with the seeds"<<endl;
+ exit(-1);
+ }
+ m_ConMat2.reinitialize(numseeds,numnz,1);
+ m_CoordMat2.reinitialize(numseeds,3,1);
+ m_CoordMat_tract2.reinitialize(numnz,3,1);
+
+ Matrix tempy(numnz,1);
+ for(int i=1;i<=numnz;i++){tempy(i,1)=i-1;}
+ m_lookup2.addvolume(m_lrmask);
+ m_lookup2.setmatrix(tempy.t(),m_lrmask);
+
+ int mytrow=0;
+ for(int Wz=m_lrmask.minz();Wz<=m_lrmask.maxz();Wz++)
+ for(int Wy=m_lrmask.miny();Wy<=m_lrmask.maxy();Wy++)
+ for(int Wx=m_lrmask.minx();Wx<=m_lrmask.maxx();Wx++)
+ if(m_lrmask(Wx,Wy,Wz)>0){
+ m_CoordMat_tract2(mytrow,0,0)=Wx;
+ m_CoordMat_tract2(mytrow,1,0)=Wy;
+ m_CoordMat_tract2(mytrow,2,0)=Wz;
+ mytrow++;
+ }
+
+ int myrow=0;
+ for(int Wz=m_seeds.minz();Wz<=m_seeds.maxz();Wz++)
+ for(int Wy=m_seeds.miny();Wy<=m_seeds.maxy();Wy++)
+ for(int Wx=m_seeds.minx();Wx<=m_seeds.maxx();Wx++)
+ if(m_seeds(Wx,Wy,Wz)>0){
+ m_CoordMat2(myrow,0,0)=Wx;
+ m_CoordMat2(myrow,1,0)=Wy;
+ m_CoordMat2(myrow,2,0)=Wz;
+ myrow++;
+ }
+
+
+ }
+
+ void Counter::count_streamline(){
+ if(opts.simpleout.value()||opts.matrix1out.value()){
+ update_pathdist();
+ }
+ if(opts.s2tout.value()){
+ update_seedcounts();
+ }
+ if(opts.matrix2out.value()){
+ update_matrix2_row();
+ }
+ if(opts.maskmatrixout.value()){
+ update_maskmatrix();
+ }
+ }
+
+ void Counter::count_seed(){
+ if(opts.matrix1out.value()){
+ update_matrix1();
+ }
+ if(opts.matrix2out.value()){
+ next_matrix2_row();
+ }
+ }
+
+
+ void Counter::clear_streamline(const bool& forwardflag,const bool& backwardflag){
+ if(opts.simpleout.value()||opts.matrix1out.value()){
+ reset_beenhere(forwardflag,backwardflag);
+ }
+ if(opts.s2tout.value()){
+ reset_targetflags();
+ }
+ if(opts.matrix2out.value()){
+ reset_beenhere2(forwardflag,backwardflag);
+ }
+ if(opts.maskmatrixout.value()){
+ //Do whatever it it you have to do!!
+ }
+ }
+
+ void Counter::update_pathdist(){
+ const vector<ColumnVector>& path=m_stline.get_path_ref();
+ for(unsigned int i=0;i<path.size();i++){
+ int x_s=int(round(float(path[i](1)))),y_s=int(round(float(path[i](2)))),z_s=int(round(float(path[i](3))));
+ if(m_beenhere(x_s,y_s,z_s)==0){
+ m_prob(x_s,y_s,z_s)+=1;
+ m_beenhere(x_s,y_s,z_s)=1;
+ }
+ }
+
+ }
+
+ void Counter::reset_beenhere(const bool& forwardflag,const bool& backwardflag){
+ if(forwardflag){
+ for(unsigned int i=0;i<m_path.size();i++){
+ int x_s=int(round(float(m_path[i](1)))),y_s=int(round(float(m_path[i](2)))),z_s=int(round(float(m_path[i](3))));
+ m_beenhere(x_s,y_s,z_s)=0;
+ }
+ }
+ if(backwardflag){
+ const vector<ColumnVector>& path=m_stline.get_path_ref();
+ for(unsigned int i=0;i<path.size();i++){
+ int x_s=int(round(float(path[i](1)))),y_s=int(round(float(path[i](2)))),z_s=int(round(float(path[i](3))));
+ m_beenhere(x_s,y_s,z_s)=0;
+ }
+ }
+ }
+
+
+ void Counter::update_seedcounts(){
+ const vector<ColumnVector>& path=m_stline.get_path_ref();
+ int xseedvox=int(round(m_stline.get_x_seed()));
+ int yseedvox=int(round(m_stline.get_y_seed()));
+ int zseedvox=int(round(m_stline.get_z_seed()));
+ for(unsigned int i=0;i<path.size();i++){
+ int x_s=int(round(float(path[i](1)))),y_s=int(round(float(path[i](2)))),z_s=int(round(float(path[i](3))));
+ for(unsigned int m=0;m<m_targetmasknames.size();m++){
+ if(m_targetmasks[m](x_s,y_s,z_s)>0 && m_targflags[m]==0){
+ m_seedcounts[m](xseedvox,yseedvox,zseedvox)=m_seedcounts[m](xseedvox,yseedvox,zseedvox)+1;
+ m_targflags[m]=1;
+ //m_particle_numbers[m].push_back(particle_number);
+ }
+ }
+ }
+
+ }
+
+
+
+ void Counter::update_matrix1(){
+ //after each particle, update_pathdist(), only run this after each voxel
+ int Concol=0;
+ for(int Wz=m_prob.minz();Wz<=m_prob.maxz();Wz++){
+ for(int Wy=m_prob.miny();Wy<=m_prob.maxy();Wy++){
+ for(int Wx=m_prob.minx();Wx<=m_prob.maxx();Wx++){
+ if(m_seeds(Wx,Wy,Wz)>0){
+ if(m_prob(Wx,Wy,Wz)>0){
+ m_ConMat(m_Conrow,Concol,0)=m_prob(Wx,Wy,Wz);
+ }
+ Concol++;
+ }
+ m_prob(Wx,Wy,Wz)=0;
+
+ }
+ }
+ }
+
+ m_Conrow++;
+ }
+
+ void Counter::update_matrix2_row(){
+ //run this one every streamline - not every voxel..
+ const vector<ColumnVector>& path=m_stline.get_path_ref();
+ for(unsigned int i=0;i<path.size();i++){
+ ColumnVector xyz_seeds=path[i];
+ ColumnVector xyz_lr=vox_to_vox(xyz_seeds,m_seedsdim,m_lrdim,m_I);
+ int x_lr=int(round(float(xyz_lr(1)))),y_lr=int(round(float(xyz_lr(2)))),z_lr=int(round(float(xyz_lr(3))));
+ int Concol2=m_lookup2(x_lr,y_lr,z_lr,0);
+ if(Concol2!=0){
+ if(m_beenhere2(x_lr,y_lr,z_lr)==0){
+ m_ConMat2(m_Conrow2,Concol2,0)+=1;
+ m_beenhere2(x_lr,y_lr,z_lr)=1;
+ }
+ }
+
+ }
+
+ }
+
+
+
+ void Counter::reset_beenhere2(const bool& forwardflag,const bool& backwardflag){
+ if(forwardflag){
+ for(unsigned int i=0;i<m_path.size();i++){
+ ColumnVector xyz_seeds=m_path[i];
+ ColumnVector xyz_lr=vox_to_vox(xyz_seeds,m_seedsdim,m_lrdim,m_I);
+ int x_lr=int(round(float(xyz_lr(1)))),y_lr=int(round(float(xyz_lr(2)))),z_lr=int(round(float(xyz_lr(3))));
+ m_beenhere2(x_lr,y_lr,z_lr)=0;
+ }
+ }
+ if(backwardflag){
+ const vector<ColumnVector>& path=m_stline.get_path_ref();
+ for(unsigned int i=0;i<path.size();i++){
+ ColumnVector xyz_seeds=path[i];
+ ColumnVector xyz_lr=vox_to_vox(xyz_seeds,m_seedsdim,m_lrdim,m_I);
+ int x_lr=int(round(float(xyz_lr(1)))),y_lr=int(round(float(xyz_lr(2)))),z_lr=int(round(float(xyz_lr(3))));
+ m_beenhere2(x_lr,y_lr,z_lr)=0;
+ }
+ }
+
+ }
+
+ void Counter::save(){
+ if(opts.simpleout.value()){
+ save_pathdist();
+ }
+ if(opts.s2tout.value()){
+ save_seedcounts();
+ }
+ if(opts.matrix1out.value()){
+ save_matrix1();
+ }
+ if(opts.matrix2out.value()){
+ save_matrix2();
+ }
+ if(opts.maskmatrixout.value()){
+ save_maskmatrix();
+ }
+
+ }
+
+ void Counter::save_pathdist(){
+ save_volume(m_prob,logger.appendDir(opts.outfile.value()));
+ }
+
+ void Counter::save_pathdist(string add){ //for simple mode
+ string thisout=opts.outfile.value();
+ make_basename(thisout);
+ thisout+=add;
+ save_volume(m_prob,thisout);
+ }
+
+ void Counter::save_seedcounts(){
+ for(unsigned int m=0;m<m_targetmasknames.size();m++){
+ string tmpname=m_targetmasknames[m];
+
+ int pos=tmpname.find("/",0);
+ int lastpos=pos;
+
+ while(pos>=0){
+ lastpos=pos;
+ pos=tmpname.find("/",pos);
+ // replace / with _
+ tmpname[pos]='_';
+ }
+
+ //only take things after the last pos
+ tmpname=tmpname.substr(lastpos+1,tmpname.length()-lastpos-1);
+
+ save_volume(m_seedcounts[m],logger.appendDir("seeds_to_"+tmpname));
+ }
+ }
+
+ void Counter::save_matrix1(){
+ save_volume(m_ConMat,logger.appendDir("fdt_matrix1"));
+ save_volume(m_CoordMat,logger.appendDir("coords_for_fdt_matrix1"));
+ }
+
+ void Counter::save_matrix2(){
+ save_volume(m_ConMat2,logger.appendDir("fdt_matrix2"));
+ save_volume(m_CoordMat2,logger.appendDir("coords_for_fdt_matrix2"));
+ save_volume(m_CoordMat_tract2,logger.appendDir("tract_space_coords_for_fdt_matrix2"));
+ save_volume4D(m_lookup2,logger.appendDir("lookup_tractspace_fdt_matrix2"));
+
+ }
+
+ void Seedmanager::run(const float& x,const float& y,const float& z,bool onewayonly){
+ //onewayonly for mesh things..
+ int fibst=m_seeds(int(round(x)),int(round(y)),int(round(z)))-1;//fibre to start with is taken from seed volume..
+ for(int p=0;p<opts.nparticles.value();p++){
+ if(opts.verbose.value()>1)
+ logger.setLogFile("particle"+num2str(p));
+
+ m_stline.reset();
+ bool forwardflag=false,backwardflag=false;
+ if(!onewayonly){
+ if(m_stline.streamline(x,y,z,m_seeddims,fibst)){
+ forwardflag=true;
+ m_counter.store_path();
+ m_counter.count_streamline();
+ }
+ m_stline.reverse();
+ }
+ if(m_stline.streamline(x,y,z,m_seeddims,fibst)){
+ backwardflag=true;
+ m_counter.count_streamline();
+ }
+
+ m_counter.clear_streamline(forwardflag,backwardflag);
+ }
+ m_counter.count_seed();
+
+
+ }
+
+
+
+}
+
+
+
diff --git a/streamlines.h b/streamlines.h
new file mode 100644
index 0000000..d2869a6
--- /dev/null
+++ b/streamlines.h
@@ -0,0 +1,201 @@
+#include <fstream>
+#include "newimage/newimageall.h"
+#include "utils/log.h"
+#include "meshclass/meshclass.h"
+#include "probtrackOptions.h"
+#include "particle.h"
+#include "tractvols.h"
+
+using namespace std;
+using namespace NEWIMAGE;
+using namespace Utilities;
+using namespace TRACTVOLS;
+using namespace mesh;
+using namespace PARTICLE;
+
+namespace TRACT{
+ void read_masks(vector<string>& masks, const string& filename);
+
+
+ class Streamliner{
+ //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;
+ Log& logger;
+ Particle m_part;
+ vector<ColumnVector> m_path;
+ volume<int> m_mask;
+ volume<int> m_skipmask;
+ volume<int> m_rubbish;
+ volume4D<float> m_loopcheck;
+ vector<volume<int>* > m_waymasks;
+ vector<bool> m_passed_flags;
+ vector<bool> m_own_waymasks;
+ Matrix m_Seeds_to_DTI;
+ TractVols vols;
+ float m_lcrat;
+ float m_x_s_init;
+ float m_y_s_init;
+ float m_z_s_init;
+ public:
+ //Constructors
+ Streamliner();
+ ~Streamliner(){
+ for(unsigned int i=0;i<m_waymasks.size();i++)
+ if(m_own_waymasks[i]) delete m_waymasks[i];
+ }
+ void add_waymask(volume<int>& myway,const bool& ownership=false){
+ //make sure that the waymask to add will not be deleted before
+ //this streamliner goes out of scope!!
+ m_waymasks.push_back(&myway);
+ m_own_waymasks.push_back(ownership);
+ m_passed_flags.push_back(false);
+ }
+ void pop_waymasks(){
+ volume<int>* tmpptr=m_waymasks[m_waymasks.size()-1];
+ m_waymasks.pop_back();
+ m_passed_flags.pop_back();
+ if(m_own_waymasks[m_own_waymasks.size()-1]){
+ delete tmpptr;
+ }
+ m_own_waymasks.pop_back();
+
+ }
+
+ inline const float get_x_seed() const {return m_x_s_init;}
+ inline const float get_y_seed() const {return m_y_s_init;}
+ inline const float get_z_seed() const {return m_z_s_init;}
+ inline const vector<ColumnVector>& get_path_ref() const{return m_path;}
+ inline vector<ColumnVector> get_path() const{return m_path;}
+ inline void reset(){
+ m_part.reset();
+ }
+ 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);
+
+ };
+
+
+ class Counter{
+ probtrackOptions& opts;
+ Log& logger;
+ volume<int> m_prob;
+ volume<int> m_beenhere;
+ Matrix m_I;
+ vector<ColumnVector> m_path;
+
+ vector<volume<int> > m_seedcounts;
+ vector<volume<int> > m_targetmasks;
+ vector<string> m_targetmasknames;
+ vector<int> m_targflags;
+ //vector<vector<int> > m_particle_numbers;
+
+
+ volume<int> m_ConMat;
+ volume<int> m_CoordMat;
+ int m_Conrow;
+
+ volume<int> m_ConMat2;
+ volume<int> m_CoordMat2;
+ volume<int> m_CoordMat_tract2;
+ volume<int> m_lrmask;
+ volume4D<int> m_lookup2;
+ volume<int> m_beenhere2;
+ int m_Conrow2;
+ ColumnVector m_lrdim;
+
+ const volume<int>& m_seeds;
+ ColumnVector m_seedsdim;
+ const Streamliner& m_stline;
+ Streamliner& m_nonconst_stline;
+
+ public:
+ Counter(const volume<int>& seeds,Streamliner& stline):opts(probtrackOptions::getInstance()),
+ logger(LogSingleton::getInstance()),
+ m_seeds(seeds),m_stline(stline),
+ m_nonconst_stline(stline){
+ //are they initialised to zero?
+ m_beenhere.reinitialize(m_seeds.xsize(),m_seeds.ysize(),m_seeds.zsize());
+ m_seedsdim.ReSize(3);
+ m_seedsdim << m_seeds.xdim() <<m_seeds.ydim() <<m_seeds.zdim();
+ m_I=Identity(4);
+
+ }
+
+ void initialise();
+
+ void initialise_path_dist(){
+ m_prob.reinitialize(m_seeds.xsize(),m_seeds.ysize(),m_seeds.zsize());
+ }
+ void initialise_seedcounts();
+
+ void initialise_matrix1(); //Need to make sure that initialise_path_dist is run first
+
+ void initialise_matrix2();
+
+ void initialise_maskmatrix(){} //not written yet
+
+ inline void store_path(){ m_path=m_stline.get_path();}
+
+ void count_streamline();
+ void count_seed();
+ void clear_streamline(const bool& forwardflag,const bool& backwardflag);
+
+
+ void update_pathdist();
+ void reset_beenhere(const bool& forwardflag,const bool& backwardflag);
+
+ void reset_prob(){m_prob=0;}
+ void update_seedcounts();
+ void reset_targetflags(){
+ for(unsigned int i=0;i<m_targflags.size();i++) m_targflags[i]=0;
+ }
+
+
+ void update_matrix1(); //update path_dist after each streamline, only run this after each voxel!!
+
+ 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_maskmatrix(){} //not written yet
+
+ void save();
+ void save_pathdist();
+ void save_pathdist(string add);
+ void save_seedcounts();
+ void save_matrix1();
+ void save_matrix2();
+ void save_maskmatrix(){}//not written yet
+
+
+ inline const Streamliner& get_streamline() const {return m_stline;}
+ inline Streamliner& get_nonconst_streamline() const {return m_nonconst_stline;}
+ inline const volume<int>& get_seeds() const {return m_seeds;}
+
+
+ };
+
+ class Seedmanager{
+ probtrackOptions& opts;
+ Log& logger;
+ Counter& m_counter;
+ Streamliner& m_stline;
+ const volume<int>& m_seeds;
+ ColumnVector m_seeddims;
+ public:
+ Seedmanager(Counter& counter):opts(probtrackOptions::getInstance()),
+ logger(LogSingleton::getInstance()),
+ m_counter(counter),
+ m_stline(m_counter.get_nonconst_streamline()),
+ m_seeds(m_counter.get_seeds()){
+ m_seeddims.ReSize(3);
+ m_seeddims<<m_seeds.xdim()<<m_seeds.ydim()<<m_seeds.zdim();
+ }
+ void run(const float& x,const float& y,const float& z,bool onewayonly=false);
+ };
+
+}
--
GitLab