diff --git a/probtrackxOptions.h b/probtrackxOptions.h
index f38519e1f01394b2ba918024f1ad0f1f04695867..9327d4c8eaa622ac32fa3bb6fc68f7129c83387b 100644
--- a/probtrackxOptions.h
+++ b/probtrackxOptions.h
@@ -46,6 +46,7 @@ class probtrackxOptions {
Option<string> stopfile;
Option<string> prefdirfile;
Option<string> seeds_to_dti;
+ Option<string> dti_to_seeds;
FmribOption<string> skipmask;
Option<string> seedref;
Option<string> mask2;
@@ -146,7 +147,11 @@ class probtrackxOptions {
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),
+ string("Transform taking seed space to DTI space (either FLIRT matrix or FNIRT warpfield) - default is identity"),
+ false, requires_argument),
+ dti_to_seeds(string("--invxfm"), string(""),
+ string("Transform taking DTI space to seed space (compulsory when using a warpfield for seeds_to_dti)"),
+ false, requires_argument),
skipmask(string("--no_integrity"), string(""),
string("no explanation needed"),
false, requires_argument),
@@ -245,6 +250,7 @@ class probtrackxOptions {
options.add(stopfile);
options.add(prefdirfile);
options.add(seeds_to_dti);
+ options.add(dti_to_seeds);
options.add(nparticles);
options.add(nsteps);
options.add(c_thr);
diff --git a/ptx_meshmask.cc b/ptx_meshmask.cc
index 8bd9850717d83a835fe071261348e139cc2e4b79..a86b0f5c1cfeaa9b6bb02a22aa2bc6a9993ec02b 100644
--- a/ptx_meshmask.cc
+++ b/ptx_meshmask.cc
@@ -36,7 +36,7 @@ void meshmask()
////////////////////////////////
// Log& logger = LogSingleton::getInstance();
- Streamliner stline;
+ Streamliner stline(seeds);
Counter counter(seeds,stline);
counter.initialise();
Seedmanager seedmanager(counter);
diff --git a/ptx_nmasks.cc b/ptx_nmasks.cc
index 67a98a28b64240cb42dcc9bce47e52397e8634b6..c01c08231d470c275b77990f427412146b91d17e 100644
--- a/ptx_nmasks.cc
+++ b/ptx_nmasks.cc
@@ -38,7 +38,7 @@ void nmasks()
seeds.push_back(tmpvol);
}
- Streamliner stline;
+ Streamliner stline(seeds[0]);
Counter counter(seeds[0],stline);
counter.initialise();
Seedmanager seedmanager(counter);
diff --git a/ptx_seedmask.cc b/ptx_seedmask.cc
index 0db11c91a4c2fc2ded39b13834629a4d3c0eeaa9..f98871ef01b85dda2fbc0bcaa86d99d60aa3fb68 100644
--- a/ptx_seedmask.cc
+++ b/ptx_seedmask.cc
@@ -22,7 +22,7 @@ void seedmask()
// Log& logger = LogSingleton::getInstance();
volume<float> seeds;
read_volume(seeds,opts.seedfile.value());
- Streamliner stline;
+ Streamliner stline(seeds);
Counter counter(seeds,stline);
counter.initialise();
Seedmanager seedmanager(counter);
diff --git a/ptx_simple.cc b/ptx_simple.cc
index 58b23c44b99506f0b66748f83dbf5660863f7be5..258598aa8164e23cbc3e5b4aef463f43f11fd4cc 100644
--- a/ptx_simple.cc
+++ b/ptx_simple.cc
@@ -30,7 +30,7 @@ void track(){
read_volume(seedref,opts.maskfile.value());
}
- Streamliner stline;
+ Streamliner stline(seedref);
Counter counter(seedref,stline);
counter.initialise();
Seedmanager seedmanager(counter);
diff --git a/ptx_twomasks.cc b/ptx_twomasks.cc
index d3b05a9b33192b443651cf9846468cb98b3a75fd..59f4c9152d7c61de495c9ba33cb8c34e538dd3dc 100644
--- a/ptx_twomasks.cc
+++ b/ptx_twomasks.cc
@@ -30,7 +30,7 @@ void twomasks()
seeds2=NEWIMAGE::abs(seeds2);
seeds2.binarise(0,seeds2.max()+1,exclusive);
- Streamliner stline;
+ Streamliner stline(seeds);
Counter counter(seeds,stline);
counter.initialise();
Seedmanager seedmanager(counter);
diff --git a/streamlines.cc b/streamlines.cc
index b218252ca0e5e934046f1b40645ee9e3dcb233c2..69dac58b61f33b80b8d855be5f634b30e52d27b8 100644
--- a/streamlines.cc
+++ b/streamlines.cc
@@ -1,5 +1,6 @@
#include "streamlines.h"
-
+#include "warpfns/fnirt_file_reader.h"
+#include "warpfns/warpfns.h"
namespace TRACT{
@@ -38,8 +39,8 @@ namespace TRACT{
}
- Streamliner::Streamliner():opts(probtrackxOptions::getInstance()),logger(LogSingleton::getInstance()),
- vols(opts.usef.value()){
+ Streamliner::Streamliner(const volume<float>& seeds):opts(probtrackxOptions::getInstance()),logger(LogSingleton::getInstance()),
+ vols(opts.usef.value()),m_seeds(seeds){
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);
@@ -77,13 +78,32 @@ namespace TRACT{
m_passed_flags.push_back(false);
m_own_waymasks.push_back(true);
}
- }
- if(opts.seeds_to_dti.value()!=""){
- m_Seeds_to_DTI = read_ascii_matrix(opts.seeds_to_dti.value());
}
- else{
- m_Seeds_to_DTI=IdentityMatrix(4);
+
+ // Allow for either matrix transform (12dof affine) or nonlinear (warpfield)
+ m_Seeds_to_DTI = IdentityMatrix(4);
+ m_DTI_to_Seeds = IdentityMatrix(4);
+
+ m_IsNonlinXfm = false;
+ if(opts.seeds_to_dti.value()!=""){
+ if(!fsl_imageexists(opts.seeds_to_dti.value())){
+ m_Seeds_to_DTI = read_ascii_matrix(opts.seeds_to_dti.value());
+ m_DTI_to_Seeds = m_Seeds_to_DTI.i();
+ }
+ 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;
+ exit(1);
+ }
+ FnirtFileReader iffr(opts.dti_to_seeds.value());
+ m_DTI_to_Seeds_warp = iffr.FieldAsNewimageVolume4D(true);
+ }
}
+
+
vols.initialise(opts.basename.value());
m_path.reserve(opts.nparticles.value());
m_x_s_init=0;
@@ -118,8 +138,12 @@ namespace TRACT{
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);
+ // find xyz in dti space
+ if(!m_IsNonlinXfm)
+ xyz_dti = vox_to_vox(xyz_seeds,dim_seeds,vols.dimensions(),m_Seeds_to_DTI);
+ else{
+ xyz_dti = NewimageCoord2NewimageCoord(m_DTI_to_Seeds_warp,false,m_seeds,m_mask,xyz_seeds);
+ }
xst=xyz_dti(1);yst=xyz_dti(2);zst=xyz_dti(3);
m_path.clear();
@@ -142,8 +166,6 @@ namespace TRACT{
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) ){
///////////////////////////////////
@@ -169,7 +191,14 @@ 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,DTI_to_Seeds);
+ // now find xyz in seeds space
+ if(!m_IsNonlinXfm)
+ xyz_seeds = vox_to_vox(xyz_dti,vols.dimensions(),dim_seeds,m_DTI_to_Seeds);
+ else{
+ xyz_seeds = NewimageCoord2NewimageCoord(m_Seeds_to_DTI_warp,false,m_mask,m_seeds,xyz_dti);
+ }
+
+
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));
@@ -812,11 +841,11 @@ namespace TRACT{
if(opts.fibst.set()){
fibst=opts.fibst.value()-1;
OUT(fibst);
- }
+ }
else{
if(fibst == -1){
fibst=0;//m_seeds(int(round(x)),int(round(y)),int(round(z)))-1;//fibre to start with is taken from seed volume..
- }
+ }
if(opts.randfib.value()){
float tmp=rand()/RAND_MAX * float(m_stline.nfibres()-1);
fibst = (int)round(tmp);
diff --git a/streamlines.h b/streamlines.h
index 39748fc8cba34a58d749ca9da17be11bb1d9ac75..a17e16705fc823de55963b9ecce0b17ecd9d9ef0 100644
--- a/streamlines.h
+++ b/streamlines.h
@@ -1,4 +1,7 @@
#include <fstream>
+#ifndef EXPOSE_TREACHEROUS
+#define EXPOSE_TREACHEROUS
+#endif
#include "newimage/newimageall.h"
#include "utils/log.h"
#include "meshclass/meshclass.h"
@@ -6,6 +9,9 @@
#include "particle.h"
#include "tractvolsx.h"
+
+
+
using namespace std;
using namespace NEWIMAGE;
using namespace Utilities;
@@ -35,15 +41,23 @@ namespace TRACT{
vector<bool> m_passed_flags;
vector<bool> m_own_waymasks;
Matrix m_Seeds_to_DTI;
+ Matrix m_DTI_to_Seeds;
+ volume4D<float> m_Seeds_to_DTI_warp;
+ volume4D<float> m_DTI_to_Seeds_warp;
+ bool m_IsNonlinXfm;
Matrix m_rotdir;
Tractvolsx vols;
float m_lcrat;
float m_x_s_init;
float m_y_s_init;
float m_z_s_init;
+
+ // we need this class to know about seed space
+ const volume<float>& m_seeds;
+
public:
//Constructors
- Streamliner();
+ Streamliner(const volume<float>&);
~Streamliner(){
for(unsigned int i=0;i<m_waymasks.size();i++)
if(m_own_waymasks[i]) delete m_waymasks[i];