From df7f66f705aa371bdfcd5440d29e014b5ac15c43 Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Fri, 27 Aug 2004 12:48:49 +0000
Subject: [PATCH] *** empty log message ***
---
probtrack.cc | 4 +-
probtrackOptions.cc | 10 +--
pt_twomasks.cc | 191 ++++++++++++++++++++++++--------------------
pt_twomasks.h | 2 +-
4 files changed, 113 insertions(+), 94 deletions(-)
diff --git a/probtrack.cc b/probtrack.cc
index a3de30f..f13babc 100644
--- a/probtrack.cc
+++ b/probtrack.cc
@@ -35,8 +35,8 @@ int main ( int argc, char **argv ){
alltracts();
else if(opts.mode.value()=="twomasks_symm")
twomasks_symm();
- else if(opts.mode.value()=="twomasks_asymm")
- twomasks_asymm();
+ else if(opts.mode.value()=="waypoints")
+ waypoints();
else if(opts.mode.value()=="matrix1")
matrix1();
else if(opts.mode.value()=="matrix2"){
diff --git a/probtrackOptions.cc b/probtrackOptions.cc
index 0834d12..91ee23e 100644
--- a/probtrackOptions.cc
+++ b/probtrackOptions.cc
@@ -104,17 +104,17 @@ probtrackOptions* probtrackOptions::gopt = NULL;
check=false;
}
}
- if(mode.value()=="twomasks_asymm"){
+ if(mode.value()=="waypoints"){
if(logdir.value()==""){
- mesg+="You must set an output directory name in twomasks_asymm mode: --dir\n";
+ mesg+="You must set an output directory name in waypoints mode: --dir\n";
check=false;
}
if(mask2.value()==""){
- mesg+="You must set a second mask in twomasks_symm mode: --mask2\n";
+ mesg+="You must set a waypoint mask or list of waypoint masks in waypoints mode: --mask2\n";
check=false;
}
if(outfile.value()==""){
- mesg+="You must set output name in twomasks_asymm mode: -o\n";
+ mesg+="You must set output name in waypoints mode: -o\n";
check=false;
}
}
@@ -175,7 +175,7 @@ probtrackOptions* probtrackOptions::gopt = NULL;
void probtrackOptions::modehelp()
{
- cout<<"tracking mode - Options are:\n\n simple (default)\n Input is text file defining start point.\n Output is dti space volume with connectivity values at each voxel.\n Note-everything in simple mode is in\n diffusion space.\n\n seeds_to_targets\n Inputs are a volume with ones at all seed points\n (requires a matrix taking it to DTI space (--xfm))\n and a text file with the name of a single target mask\n on each new line (--targetmasks).\n Output is a single volume for each target mask where the value of each\n volume within the seed mask corresponds to the number of particles \n seeded from this voxel reaching this target mask.\n\n seedmask\n Input is a volume with ones at all seed points\n (requires a matrix taking it to DTI space (--xfm))\n Output is seed space volume with connectivity values at each voxel\n summed from every seed.\n\n twomasks_symm\n Input is a binary volume for first mask (-x)\n also requires a binary volume for second mask (--mask2)\n (requires a matrix taking both seeds to DTI space (--xfm))\n Output is the sum of all paths from the ones which pass\n through the twos and vice-versa i.e. only paths\n which pass through both masks are retained.\n\n twomasks_asymm\n Exactly as above but tracking is asymmetric\n i.e. paths are seeded from the first mask and only\n paths which pass through the second mask (--mask2) are kept.\n\n matrix1\n Input is a volume with ones at all seed points\n (requires a matrix taking it to DTI space (--xfm))\n Output is an avw format nseeds x nseeds matrix \n where element ij contains the number of particles leaving seed \n voxel i and passing through seed voxel j. \n Second output is an nsees x 3 matrix containing the \n anatomical locations in Seed space of each matrix node.\n\n matrix2\n Inputs are a volume with ones at all seed points\n (requires a matrix taking it to DTI space (--xfm))\n and a low resolution binary mask in register with the\n seed volume (--lrmask)\n Output is a matrix of size nseeds x nlowres \n which contains the connectivity distributions from every seed voxel.\n Auxiliary outputs are index matrices for each axis \n of this matrix containing the anatomical location of each\n node in seed space (x-axis) and low res space (y-axis)\n\n maskmatrix\n Input is volume with integers in each seed location\n Integer values distinguish different clusters\n (e.g. output of Feat clustering)\n (requires a matrix taking it to DTI space (--xfm))\n Output is two matrices whose ij^th element is the\n mean and max number of particles seeded from cluster i which pass \n through cluster j\n"<<endl;
+ cout<<"tracking mode - Options are:\n\n simple (default)\n Input is text file defining start point.\n Output is dti space volume with connectivity values at each voxel.\n Note-everything in simple mode is in\n diffusion space.\n\n seeds_to_targets\n Inputs are a volume with ones at all seed points\n (requires a matrix taking it to DTI space (--xfm))\n and a text file with the name of a single target mask\n on each new line (--targetmasks).\n Output is a single volume for each target mask where the value of each\n volume within the seed mask corresponds to the number of particles \n seeded from this voxel reaching this target mask.\n\n seedmask\n Input is a volume with ones at all seed points\n (requires a matrix taking it to DTI space (--xfm))\n Output is seed space volume with connectivity values at each voxel\n summed from every seed.\n\n twomasks_symm\n Input is a binary volume for first mask (-x)\n also requires a binary volume for second mask (--mask2)\n (requires a matrix taking both seeds to DTI space (--xfm))\n Output is the sum of all paths the from first mask which pass\n through the second and vice-versa i.e. only paths\n which pass through both masks are retained.\n\n waypoints\n Input is a binary volume for seed mask (-x)\n also requires a binary waypoint mask or ascii list\n of waypoint masks (--mask2)\n (requires a matrix taking both seeds to DTI space (--xfm))\n Output is a volume containing all the paths from the\n seedmask which pass through ALL of the waypoint masks. \n\n matrix1\n Input is a volume with ones at all seed points\n (requires a matrix taking it to DTI space (--xfm))\n Output is an avw format nseeds x nseeds matrix \n where element ij contains the number of particles leaving seed \n voxel i and passing through seed voxel j. \n Second output is an nsees x 3 matrix containing the \n anatomical locations in Seed space of each matrix node.\n\n matrix2\n Inputs are a volume with ones at all seed points\n (requires a matrix taking it to DTI space (--xfm))\n and a low resolution binary mask in register with the\n seed volume (--lrmask)\n Output is a matrix of size nseeds x nlowres \n which contains the connectivity distributions from every seed voxel.\n Auxiliary outputs are index matrices for each axis \n of this matrix containing the anatomical location of each\n node in seed space (x-axis) and low res space (y-axis)\n\n maskmatrix\n Input is volume with integers in each seed location\n Integer values distinguish different clusters\n (e.g. output of Feat clustering)\n (requires a matrix taking it to DTI space (--xfm))\n Output is two matrices whose ij^th element is the\n mean and max number of particles seeded from cluster i which pass \n through cluster j\n"<<endl;
}
diff --git a/pt_twomasks.cc b/pt_twomasks.cc
index d00a875..8b7ed47 100644
--- a/pt_twomasks.cc
+++ b/pt_twomasks.cc
@@ -1,4 +1,5 @@
#include "pt_twomasks.h"
+#include "pt_seeds_to_targets.h"
using namespace std;
using namespace NEWIMAGE;
@@ -12,7 +13,7 @@ void twomasks_symm(){
probtrackOptions& opts = probtrackOptions::getInstance();
Log& logger = LogSingleton::getInstance();
-
+
float xst,yst,zst,x,y,z;
int nparticles=opts.nparticles.value();
int nsteps=opts.nsteps.value();
@@ -27,7 +28,7 @@ void twomasks_symm(){
loopcheck.reinitialize(int(ceil(mask.xsize()/lcrat)+1),int(ceil(mask.ysize()/lcrat)+1),int(ceil(mask.zsize()/lcrat)+1),3);
loopcheck=0;
}
-
+
volume<int> Seeds,mask2;
read_volume(Seeds,opts.seedfile.value());
@@ -93,18 +94,16 @@ void twomasks_symm(){
xyz_dti=vox_to_vox(xyz_seeds,dim_seeds,vols.dimensions(),Seeds_to_DTI);
xst=xyz_dti(1);yst=xyz_dti(2);zst=xyz_dti(3);
Particle part(0,0,0,0,0,0,opts.steplength.value(),mask.xdim(),mask.ydim(),mask.zdim(),false);
-
+
for( int p = 0; p < nparticles ; p++ ){
- bool keepflag=false;
- ColumnVector partlength(2);
- partlength=0;
for(int direc=1;direc<=2;direc++){
x=xst;y=yst;z=zst;
part.change_xyz(x,y,z);
if(direc==2){
part.restart_reverse(); //go again in the opposite direction
}
-
+ int partlength=0;
+ bool keepflag=false; /// only keep it if this direction went through the masks
for( int it = 1 ; it <= nsteps/2; it++){
if( (mask( round(part.x()), round(part.y()), round(part.z())) > 0) ){
@@ -137,10 +136,10 @@ void twomasks_symm(){
if(opts.rubbishfile.value()!=""){
if(RUBBISH(x_s,y_s,z_s)>0) break;
}
- path(it+(direc-1)*nsteps/2,1)=x_s;
- path(it+(direc-1)*nsteps/2,2)=y_s;
- path(it+(direc-1)*nsteps/2,3)=z_s;
- partlength(direc)+=1;
+ path(it,1)=x_s;
+ path(it,2)=y_s;
+ path(it,3)=z_s;
+ partlength+=1;
if( Seeds(x_s,y_s,z_s)==(maskno*-1 + 3) ){ //mult by -1 and add 3 makes 1 into 2 and 2 into 1
keepflag=true;
}
@@ -188,26 +187,25 @@ void twomasks_symm(){
loopcheck=0;
}
- }
+
// Do the counting after a particle has finished
// And only if keepflag is set
- if(keepflag){
-
- for(int direc=1;direc<=2;direc++){
- for(int s=1;s<=int(partlength(direc));s++){
- int x_s=int(path(s+(direc-1)*nsteps/2,1));
- int y_s=int(path(s+(direc-1)*nsteps/2,2));
- int z_s=int(path(s+(direc-1)*nsteps/2,3));
+ if(keepflag){
+ for(int s=1;s<=partlength;s++){
+ int x_s=int(path(s,1));
+ int y_s=int(path(s,2));
+ int z_s=int(path(s,3));
if(beenhere(x_s,y_s,z_s)==0){
prob(x_s,y_s,z_s)+=1;
beenhere(x_s,y_s,z_s)=1;
}
- }
+ }
+
+ indexset(beenhere,path,0);
}
-
- indexset(beenhere,path,0);
- }
- part.reset();
+ } //close direction loop
+
+ part.reset();
@@ -225,11 +223,10 @@ void twomasks_symm(){
}
-void twomasks_asymm(){
-
+void waypoints(){
probtrackOptions& opts = probtrackOptions::getInstance();
Log& logger = LogSingleton::getInstance();
-
+
float xst,yst,zst,x,y,z;
int nparticles=opts.nparticles.value();
int nsteps=opts.nsteps.value();
@@ -244,17 +241,27 @@ void twomasks_asymm(){
loopcheck.reinitialize(int(ceil(mask.xsize()/lcrat)+1),int(ceil(mask.ysize()/lcrat)+1),int(ceil(mask.zsize()/lcrat)+1),3);
loopcheck=0;
}
-
volume<int> Seeds,mask2;
read_volume(Seeds,opts.seedfile.value());
- if(opts.mask2.value()!=""){
- read_volume(mask2,opts.mask2.value());
- Seeds=Seeds+(2*mask2);
- }
-
-
+
+ vector<string> masknames;
+ if(fsl_imageexists(opts.mask2.value())){
+ masknames.push_back( opts.mask2.value() );
+ }
+ else{
+ read_masks(masknames,opts.mask2.value());
+ }
+ vector<volume<int> > waymasks;
+ volume<int> tmpway;
+ vector<bool> passed_flags; //store whether the current sample has passed through each of the waypoints
+ for( unsigned int m = 0; m < masknames.size(); m++ ){
+ read_volume(tmpway,masknames[m]);
+ waymasks.push_back(tmpway);
+ passed_flags.push_back(false);
+ }
+
volume<int> RUBBISH;
if(opts.rubbishfile.value()!=""){
read_volume(RUBBISH,opts.rubbishfile.value());
@@ -265,30 +272,33 @@ void twomasks_asymm(){
beenhere=Seeds;
beenhere=0;
prob=beenhere;
- // int numseeds=0;
- // for(int Wz=Seeds.minz();Wz<=Seeds.maxz();Wz++){
- // for(int Wy=Seeds.miny();Wy<=Seeds.maxy();Wy++){
- // for(int Wx=Seeds.minx();Wx<=Seeds.maxx();Wx++){
- // if(Seeds.value(Wx,Wy,Wz)>0){
- // numseeds++;
- // }
- // }
- // }
- // }
-
+// int numseeds=0;
+// for(int Wz=Seeds.minz();Wz<=Seeds.maxz();Wz++){
+// for(int Wy=Seeds.miny();Wy<=Seeds.maxy();Wy++){
+// for(int Wx=Seeds.minx();Wx<=Seeds.maxx();Wx++){
+// if(Seeds.value(Wx,Wy,Wz)>0){
+// numseeds++;
+// }
+// }
+// }
+// }
+
cout<<"Loading MCMC volumes"<<endl;
TractVols vols(opts.usef.value());
vols.initialise(opts.basename.value());
-
+
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);
}
-
+
+
+
Matrix path(nsteps,3);
path=1;
@@ -299,8 +309,7 @@ void twomasks_asymm(){
cout<<Sz<<endl;
for(int Sy=Seeds.miny();Sy<=Seeds.maxy();Sy++){
for(int Sx=Seeds.minx();Sx<=Seeds.maxx();Sx++){
- if(Seeds(Sx,Sy,Sz)==1){
- int maskno=Seeds(Sx,Sy,Sz);
+ if(Seeds(Sx,Sy,Sz)>0){
ColumnVector xyz_seeds(3),dim_seeds(3),xyz_dti;
xyz_seeds << Sx << Sy << Sz;
dim_seeds <<Seeds.xdim()<<Seeds.ydim()<<Seeds.zdim();
@@ -308,18 +317,20 @@ void twomasks_asymm(){
xyz_dti=vox_to_vox(xyz_seeds,dim_seeds,vols.dimensions(),Seeds_to_DTI);
xst=xyz_dti(1);yst=xyz_dti(2);zst=xyz_dti(3);
Particle part(0,0,0,0,0,0,opts.steplength.value(),mask.xdim(),mask.ydim(),mask.zdim(),false);
-
+
for( int p = 0; p < nparticles ; p++ ){
- bool keepflag=false;
- ColumnVector partlength(2);
- partlength=0;
+
for(int direc=1;direc<=2;direc++){
x=xst;y=yst;z=zst;
part.change_xyz(x,y,z);
if(direc==2){
part.restart_reverse(); //go again in the opposite direction
}
-
+ int partlength=0;
+ for(unsigned int pf=0;pf<passed_flags.size();pf++) {
+ passed_flags[pf]=false; /// only keep it if this direction went through all the masks
+ }
+
for( int it = 1 ; it <= nsteps/2; it++){
if( (mask( round(part.x()), round(part.y()), round(part.z())) > 0) ){
@@ -341,9 +352,8 @@ void twomasks_asymm(){
loopcheck((int)round(part.x()/lcrat),(int)round(part.y()/lcrat),(int)round(part.z()/lcrat),2)=part.rz();
}
-
-
-
+
+
x=part.x();y=part.y();z=part.z();
xyz_dti <<x<<y<<z;
xyz_seeds=vox_to_vox(xyz_dti,vols.dimensions(),dim_seeds,Seeds_to_DTI.i());
@@ -353,25 +363,27 @@ void twomasks_asymm(){
if(opts.rubbishfile.value()!=""){
if(RUBBISH(x_s,y_s,z_s)>0) break;
}
-
- path(it+(direc-1)*nsteps/2,1)=x_s;
- path(it+(direc-1)*nsteps/2,2)=y_s;
- path(it+(direc-1)*nsteps/2,3)=z_s;
- partlength(direc)+=1;
- if( Seeds(x_s,y_s,z_s)==(maskno*-1 + 3) ){ //mult by -1 and add 3 makes 1 into 2 and 2 into 1
- keepflag=true;
+ path(it,1)=x_s;
+ path(it,2)=y_s;
+ path(it,3)=z_s;
+ partlength+=1;
+
+ //update every passed_flag
+ for( unsigned int wm=0;wm<waymasks.size();wm++ ){
+ if( waymasks[wm](x_s,y_s,z_s)>0 )
+ passed_flags[wm]=true;
}
-
-
-
+
+
if(opts.skipmask.value() == ""){
th_ph_f=vols.sample(part.x(),part.y(),part.z());
}
else{
- if(skipmask((int)round(part.x()),(int)round(part.y()),(int)round(part.z()))==0)
+ if(skipmask(x_s,y_s,z_s)==0)
th_ph_f=vols.sample(part.x(),part.y(),part.z());
}
-
+
+
tmp2=rand(); tmp2/=RAND_MAX;
if(th_ph_f(3)>tmp2){
if(!part.check_dir(th_ph_f(1),th_ph_f(2),opts.c_thr.value())){
@@ -389,7 +401,12 @@ void twomasks_asymm(){
test_th_ph_f=vols.sample(part.testx(),part.testy(),part.testz());
part.jump(test_th_ph_f(1),test_th_ph_f(2));
}
- }
+ }
+
+
+
+
+
}
}
@@ -402,29 +419,31 @@ void twomasks_asymm(){
loopcheck=0;
}
- }
+
// Do the counting after a particle has finished
- // And only if keepflag is set
- if(keepflag){
-
- for(int direc=1;direc<=2;direc++){
- for(int s=1;s<=int(partlength(direc));s++){
- int x_s=int(path(s+(direc-1)*nsteps/2,1));
- int y_s=int(path(s+(direc-1)*nsteps/2,2));
- int z_s=int(path(s+(direc-1)*nsteps/2,3));
+ // And only if all passed_flags are set
+ bool keepflag=true;
+ for(unsigned int pf=0;pf<passed_flags.size();pf++){
+ keepflag=(keepflag && passed_flags[pf]);
+ }
+
+ if(keepflag){
+ for(int s=1;s<=partlength;s++){
+ int x_s=int(path(s,1));
+ int y_s=int(path(s,2));
+ int z_s=int(path(s,3));
if(beenhere(x_s,y_s,z_s)==0){
prob(x_s,y_s,z_s)+=1;
beenhere(x_s,y_s,z_s)=1;
}
- }
+ }
+ indexset(beenhere,path,0);
}
-
- indexset(beenhere,path,0);
- }
- part.reset();
+ } //close direction loop
-
-
+ part.reset();
+
+
} // Close Particle Number Loop
diff --git a/pt_twomasks.h b/pt_twomasks.h
index b550e20..c8f24f9 100644
--- a/pt_twomasks.h
+++ b/pt_twomasks.h
@@ -7,4 +7,4 @@
#include "tractvols.h"
void twomasks_symm();
-void twomasks_asymm();
+void waypoints();
--
GitLab