From 855016b09e78c13b5fec406f54efa62d8eb638ce Mon Sep 17 00:00:00 2001
From: Saad Jbabdi <saad@fmrib.ox.ac.uk>
Date: Wed, 2 Apr 2008 15:16:46 +0000
Subject: [PATCH] changes for Jill
---
probtrackxOptions.h | 5 +++++
streamlines.cc | 16 ++++++++++++----
streamlines.h | 1 +
tractvolsx.h | 11 ++++++-----
4 files changed, 24 insertions(+), 9 deletions(-)
diff --git a/probtrackxOptions.h b/probtrackxOptions.h
index beaacc1..0c19c10 100644
--- a/probtrackxOptions.h
+++ b/probtrackxOptions.h
@@ -44,6 +44,7 @@ class probtrackxOptions {
Option<string> outfile;
Option<string> rubbishfile;
Option<string> stopfile;
+ Option<string> prefdirfile;
Option<string> seeds_to_dti;
FmribOption<string> skipmask;
Option<string> seedref;
@@ -138,6 +139,9 @@ class probtrackxOptions {
stopfile(string("--stop"), string(""),
string("Stop tracking at locations given by this mask file"),
false, requires_argument),
+ prefdirfile(string("--prefdir"), string(""),
+ string("prefered orientation preset in a 4D mask"),
+ false, requires_argument),
seeds_to_dti(string("--xfm"), string(""),
string("Transform Matrix taking seed space to DTI space default is to use the identity"),false, requires_argument),
skipmask(string("--no_integrity"), string(""),
@@ -230,6 +234,7 @@ class probtrackxOptions {
options.add(outfile);
options.add(rubbishfile);
options.add(stopfile);
+ options.add(prefdirfile);
options.add(seeds_to_dti);
options.add(nparticles);
options.add(nsteps);
diff --git a/streamlines.cc b/streamlines.cc
index bbdae30..16d0caf 100644
--- a/streamlines.cc
+++ b/streamlines.cc
@@ -40,6 +40,9 @@ namespace TRACT{
if(opts.stopfile.value()!=""){
read_volume(m_stop,opts.stopfile.value());
}
+ if(opts.prefdirfile.value()!=""){
+ read_volume4D(m_prefdir,opts.prefdirfile.value());
+ }
vector<string> masknames;
if(opts.waypoints.value()!=""){
@@ -137,7 +140,12 @@ namespace TRACT{
int y_s =(int)round((float)xyz_seeds(2));
int z_s =(int)round((float)xyz_seeds(3));
-
+ float pref_x=0,pref_y=0,pref_z=0;
+ if(opts.prefdirfile.value()!=""){
+ pref_x = m_prefdir(xyz_seeds(1),xyz_seeds(2),xyz_seeds(3),0);
+ pref_y = m_prefdir(xyz_seeds(1),xyz_seeds(2),xyz_seeds(3),1);
+ pref_z = m_prefdir(xyz_seeds(1),xyz_seeds(2),xyz_seeds(3),2);
+ }
//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 ) {
@@ -166,11 +174,11 @@ namespace TRACT{
}
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());
+ 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);
}
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());
+ 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);
}
@@ -188,7 +196,7 @@ namespace TRACT{
{
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());
+ test_th_ph_f=vols.sample(m_part.testx(),m_part.testy(),m_part.testz(),m_part.rx(),m_part.ry(),m_part.rz(),pref_x,pref_y,pref_z);
m_part.jump(test_th_ph_f(1),test_th_ph_f(2));
}
}
diff --git a/streamlines.h b/streamlines.h
index faca837..8f4548e 100644
--- a/streamlines.h
+++ b/streamlines.h
@@ -29,6 +29,7 @@ namespace TRACT{
volume<int> m_skipmask;
volume<int> m_rubbish;
volume<int> m_stop;
+ volume4D<float> m_prefdir;
volume4D<float> m_loopcheck;
vector<volume<float>* > m_waymasks;
vector<bool> m_passed_flags;
diff --git a/tractvolsx.h b/tractvolsx.h
index c8b56a7..f4cf8f3 100644
--- a/tractvolsx.h
+++ b/tractvolsx.h
@@ -103,7 +103,8 @@ namespace TRACTVOLSX{
}
- ColumnVector sample(const float& x,const float& y,const float&z,const float& r_x,const float& r_y,const float& r_z){
+ ColumnVector sample(const float& x,const float& y,const float&z,const float& r_x,const float& r_y,const float& r_z,
+ float& prefer_x=0,float& prefer_y=0,float& prefer_z=0){
////////Probabilistic interpolation
int cx =(int) ceil(x),fx=(int) floor(x);
@@ -148,8 +149,6 @@ namespace TRACTVOLSX{
newz=cz;
ColumnVector th_ph_f(3);
-
-
float samp=rand(); samp/=RAND_MAX;
samp=round(samp*((*thsamples[0]).tsize()-1));
float theta=0,phi=0;
@@ -162,12 +161,14 @@ namespace TRACTVOLSX{
init_sample=false;
}
else{
-
+ if(fabs(prefer_x)+fabs(prefer_y)+fabs(prefer_z)==0){
+ prefer_x=r_x;prefer_y=r_y;prefer_z=r_z;
+ }
for(unsigned int fib=0;fib<thsamples.size();fib++){
if((*fsamples[fib])(int(newx),int(newy),int(newz),int(samp))>opts.fibthresh.value()){
float phtmp=(*phsamples[fib])(int(newx),int(newy),int(newz),int(samp));
float thtmp=(*thsamples[fib])(int(newx),int(newy),int(newz),int(samp));
- dottmp=fabs(sin(thtmp)*cos(phtmp)*r_x + sin(thtmp)*sin(phtmp)*r_y + cos(thtmp)*r_z);
+ dottmp=fabs(sin(thtmp)*cos(phtmp)*prefer_x + sin(thtmp)*sin(phtmp)*prefer_y + cos(thtmp)*prefer_z);
if(dottmp>dotmax){
dotmax=dottmp;
theta=thtmp;
--
GitLab