From 57b767823b9cbedd92e4d1a4482e8bc1d67155db Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Mon, 28 Sep 2009 16:38:05 +0000
Subject: [PATCH] changed randfib functionality
/bin/bash: i: command not found
---
probtrackxOptions.h | 6 ++---
streamlines.cc | 11 +++------
tractvolsx.h | 56 +++++++++++++++++++++++++++++++++++++++++++--
3 files changed, 60 insertions(+), 13 deletions(-)
diff --git a/probtrackxOptions.h b/probtrackxOptions.h
index fd9ee8d..c232363 100644
--- a/probtrackxOptions.h
+++ b/probtrackxOptions.h
@@ -63,7 +63,7 @@ class probtrackxOptions {
Option<float> steplength;
Option<bool> loopcheck;
Option<bool> usef;
- Option<bool> randfib;
+ Option<int> randfib;
Option<int> fibst;
Option<bool> modeuler;
Option<int> rseed;
@@ -200,8 +200,8 @@ class probtrackxOptions {
usef(string("-f,--usef"), false,
string("Use anisotropy to constrain tracking"),
false, no_argument),
- randfib(string("--randfib"), false,
- string("Select randomly from one of the fibres"),
+ randfib(string("--randfib"), 0,
+ string("Set to 1 to randomly sample initial fibres. Set to 2 to sample in proportion to f"),
false, no_argument),
fibst(string("--fibst"),1,
string("Force a starting fibre for tracking - default=1, i.e. first fibre orientation"),
diff --git a/streamlines.cc b/streamlines.cc
index 3709304..13716c2 100644
--- a/streamlines.cc
+++ b/streamlines.cc
@@ -1001,14 +1001,9 @@ namespace TRACT{
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);
- //if(tmp>0.5)
- //fibst=0;
- //else
- //fibst=1;// fix this for > 2 fibres
- }
+ //TB moved randfib option inside tractvols.h 28/10/2009
+ // This means that we have access to fsamples when figuring out fibst
+ // so we can choose to seed in proportion to f in that voxel.
}
// now re-orient dir using xfm transform
diff --git a/tractvolsx.h b/tractvolsx.h
index f9cabd3..624c17f 100644
--- a/tractvolsx.h
+++ b/tractvolsx.h
@@ -35,7 +35,7 @@ namespace TRACTVOLSX{
public:
//constructors::
- Tractvolsx(const bool& usefin=false):opts(probtrackxOptions::getInstance()),init_sample(true),fibst(1),usef(usefin){}
+ Tractvolsx(const bool& usefin=false):opts(probtrackxOptions::getInstance()),init_sample(true),fibst(0),usef(usefin){}
Tractvolsx():opts(probtrackxOptions::getInstance()){}
~Tractvolsx(){
for(unsigned int m=0;m<thsamples.size();m++)
@@ -155,10 +155,62 @@ namespace TRACTVOLSX{
float dotmax=0,dottmp=0;
int fibind=0;
if(thsamples.size()>1){//more than 1 fibre
- if(init_sample){//go for the specified option on the first jump
+
+ if(init_sample){//go for the specified fibre on the first jump or generate at random
+
+ if(opts.randfib.value()==1){//this generates startfib at random (except for fibres where f<fibthresh)
+ vector<int> fibvec;
+ for(unsigned int fib=0;fib<thsamples.size();fib++){
+ float ft=(*fsamples[fib])(int(newx),int(newy),int(newz),int(samp));
+ if(ft>opts.fibthresh.value()){
+ fibvec.push_back(fib);
+ }
+ }
+
+ if(fibvec.size()==0){
+ fibst=0;
+ }
+ else{
+ float rtmp=rand()/RAND_MAX * float(fibvec.size()-1);
+ fibst = fibvec[ (int)round(rtmp) ];
+ }
+
+ }
+
+
+ else if(opts.randfib.value()==2){ //this generates startfib with probability proportional to f (except for fibres where f<fibthresh).
+ //this chooses at random but in proportion to fsamples.
+ float fsumtmp=0;
+ for(unsigned int fib=0;fib<thsamples.size();fib++){
+ float ft=(*fsamples[fib])(int(newx),int(newy),int(newz),int(samp));
+ if(ft>opts.fibthresh.value()){
+ fsumtmp+=ft; //count total weight of f in this voxel.
+ }
+ }
+
+ if(fsumtmp==0){
+ fibst=0;
+ }
+ else{
+ float fsumtmp2=0;
+ int fib=0;
+ float rtmp=rand()/RAND_MAX;
+
+ while( fsumtmp2<rtmp){
+ float ft=(*fsamples[fib])(int(newx),int(newy),int(newz),int(samp));
+ if(ft>opts.fibthresh.value()){
+ fsumtmp2+=(ft/fsumtmp);
+ }
+ fibst=fib;
+ fib++;
+ }
+
+ }
+
theta=(*thsamples[fibst])(int(newx),int(newy),int(newz),int(samp));
phi=(*phsamples[fibst])(int(newx),int(newy),int(newz),int(samp));
init_sample=false;
+ }
}
else{
if((fabs(prefer_x)+fabs(prefer_y)+fabs(prefer_z))==0){
--
GitLab