From f71525ecfc51b88ceab3bf07b2b8a70f5f5fee2d Mon Sep 17 00:00:00 2001
From: Tim Behrens <behrens@fmrib.ox.ac.uk>
Date: Fri, 10 Jun 2005 14:37:51 +0000
Subject: [PATCH] *** empty log message ***
---
Makefile | 2 +-
fibre.h | 142 +++++++++--------
xfibres.cc | 386 ++++++++++++++++++++++++++++++++++++++++++++++
xfibresoptions.cc | 56 +++++++
xfibresoptions.h | 146 ++++++++++++++++++
5 files changed, 658 insertions(+), 74 deletions(-)
create mode 100644 xfibres.cc
create mode 100644 xfibresoptions.cc
create mode 100644 xfibresoptions.h
diff --git a/Makefile b/Makefile
index fcc8a28..cb928ec 100644
--- a/Makefile
+++ b/Makefile
@@ -30,7 +30,7 @@ MEDOBJS=medianfilter.o
ROMOBJS=reord_OM.o
SAUSOBJS=sausages.o
DIFF_PVMOBJS=diff_pvm.o diff_pvmoptions.o
-XFIBOBJS=xfibres.o
+XFIBOBJS=xfibres.o xfibresoptions.o
RVOBJS=replacevols.o
MDVOBJS=make_dyadic_vectors.o
diff --git a/fibre.h b/fibre.h
index 949b56f..e0bd388 100644
--- a/fibre.h
+++ b/fibre.h
@@ -57,8 +57,8 @@ namespace FIBRE{
const Matrix& m_bvals;
public:
//constructors::
- Fibre(const float& d, const ColumnVector& alpha,
- const ColumnVector& beta, const Matrix& bvals):
+ Fibre( const ColumnVector& alpha, const ColumnVector& beta,
+ const Matrix& bvals,const float& d):
m_d(d), m_alpha(alpha), m_beta(beta), m_bvals(bvals){
m_th=0;
@@ -94,10 +94,10 @@ namespace FIBRE{
m_Signal=0;
m_Signal_old=m_Signal;
}
- Fibre(const float& d, const ColumnVector& alpha,
- const ColumnVector& beta, const Matrix& bvals,
+ Fibre(const ColumnVector& alpha,
+ const ColumnVector& beta, const Matrix& bvals, const float& d,
const float& th, const float& ph, const float& f,
- const float& lam, const int Npts) :
+ const float& lam) :
m_th(th), m_ph(ph), m_f(f), m_lam(lam), m_d(d),
m_alpha(alpha), m_beta(beta), m_bvals(bvals)
{
@@ -131,6 +131,12 @@ namespace FIBRE{
m_Signal=0;
m_Signal_old=m_Signal;
}
+ Fibre(const Fibre& rhs):
+ m_d(rhs.m_d), m_alpha(rhs.m_alpha), m_beta(rhs.m_beta), m_bvals(rhs.m_bvals){
+ *this=rhs;
+ }
+
+
~Fibre(){}
inline float get_th() const{ return m_th;}
@@ -142,8 +148,12 @@ namespace FIBRE{
inline float get_f() const{ return m_f;}
inline void set_f(const float f){ m_f=f; }
- inline ColumnVector getSignal() const{ //flitney had "inline ColumnVector&" here
- return m_Signal; // What is the difference?
+ inline float get_lam() const{ return m_lam;}
+ inline void set_lam(const float lam){ m_lam=lam; }
+
+
+ inline const ColumnVector& getSignal() const{
+ return m_Signal;
}
inline void restoreSignal() {
@@ -305,7 +315,40 @@ namespace FIBRE{
}
-
+ Fibre& operator=(const Fibre& rhs){
+ m_th=rhs.m_th;
+ m_ph=rhs.m_ph;
+ m_f=rhs.m_f;
+ 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;
+ m_lam_prop=rhs.m_lam_prop;
+ m_th_old=rhs.m_th_old;
+ m_ph_old=rhs.m_ph_old;
+ m_f_old=rhs.m_f_old;
+ m_lam_old=rhs.m_lam_old;
+ m_th_prior=rhs.m_th_prior;
+ m_ph_prior=rhs.m_ph_prior;
+ m_f_prior=rhs. m_f_prior;
+ m_lam_prior=rhs.m_lam_prior;
+ m_th_old_prior=rhs.m_th_old_prior;
+ m_ph_old_prior=rhs.m_ph_old_prior;
+ m_f_old_prior=rhs.m_f_old_prior;
+ m_lam_old_prior=rhs.m_lam_old_prior;
+ m_prior_en=rhs.m_prior_en;
+ m_old_prior_en=rhs. m_old_prior_en;
+ m_th_acc=rhs.m_th_acc;
+ m_th_rej=rhs.m_th_rej;
+ m_ph_acc=rhs.m_ph_acc;
+ m_ph_rej=rhs. m_ph_rej;
+ m_f_acc=rhs.m_f_acc;
+ m_f_rej=rhs.m_f_rej;
+ m_lam_acc=rhs.m_lam_acc;
+ m_lam_rej=rhs.m_lam_rej;
+ m_Signal=rhs.m_Signal;
+ m_Signal_old=rhs. m_Signal_old;
+ }
friend ostream& operator<<(ostream& ostr,const Fibre& p);
@@ -347,7 +390,7 @@ namespace FIBRE{
public:
Multifibre(const ColumnVector& data,const ColumnVector& alpha,
- const ColumnVector& beta, const Matrix& b, const Matrix& Amat, int N ):
+ const ColumnVector& beta, const Matrix& b, int N ):
m_data(data), m_alpha(alpha), m_beta(beta), m_bvals(b){
// initialise(Amat,N);
@@ -355,72 +398,21 @@ namespace FIBRE{
~Multifibre(){}
- /*
- void initialise(const Matrix& Amat, int N){
- //initialising
- ColumnVector logS(m_data.Nrows()),tmp(m_data.Nrows()),Dvec(7),dir(3);
- SymmetricMatrix tens; //Basser's Diffusion Tensor;
- DiagonalMatrix Dd; //eigenvalues
- Matrix Vd; //eigenvectors
- float mDd,fsquared;
- float th,ph,f,D,S0;
-
- for ( int i = 1; i <= logS.Nrows(); i++)
- {
- if(m_data(i)>0){
- logS(i)=log(m_data(i));
- }
- else{
- logS(i)=0;
- }
- }
- Dvec = -pinv(Amat)*logS;
- if( Dvec(7) > -23 ){ //23=maxlogfloat
- S0=exp(-Dvec(7));
- }
- else{
- S0=m_data.MaximumAbsoluteValue();
- }
-
- for ( int i = 1; i <= logS.Nrows(); i++)
- {
- if(S0<m_data.Sum()/m_data.Nrows()){ S0=m_data.MaximumAbsoluteValue(); }
- logS(i)=(m_data(i)/S0)>0.01 ? log((i)):log(0.01*S0);
- }
- Dvec = -pinv(Amat)*logS;
- S0=exp(-Dvec(7));
- if(S0<m_data.Sum()/S.Nrows()){ S0=m_data.Sum()/m_data.Nrows(); }
- tens = vec2tens(Dvec);
- EigenValues(tens,Dd,Vd);
- mDd = Dd.Sum()/Dd.Nrows();
- int maxind = Dd(1) > Dd(2) ? 1:2; //finding maximum eigenvalue
- maxind = Dd(maxind) > Dd(3) ? maxind:3;
- dir << Vd(1,maxind) << Vd(2,maxind) << Vd(3,maxind);
- cart2sph(dir,th,ph);
- th= mod(th,M_PI);
- ph= mod(ph,2*M_PI);
- D = Dd(maxind);
-
- float numer=1.5(*(Dd(1)-mDd)*(Dd(1)-mDd)+(Dd(2)-mDd)*(Dd(2)-mDd)+(Dd(3)-mDd)*(Dd(3)-mDd));
- float denom=(Dd(1)*Dd(1)+Dd(2)*Dd(2)+Dd(3)*Dd(3));
- if(denom>0) fsquared=numer/denom;
- else fsquared=0;
- if(fsquared>0){f=sqrt(fsquared);}
- else{f=0;}
- if(f>=0.95) f=0.95;
- if(f<=0.001) f=0.001;
-
-
- m_d=D; m_d_old=m_d;
- m_S0=S0; m_S0_old=m_S0;
-
-
-
-
-
+ const vector<Fibre>& fibres() const{
+ return m_fibres;
}
- */
+ void addfibre(const Fibre& fib){
+ m_fibres.push_back(fib);
+ }
+
+ inline float get_d() const{ return m_d;}
+ inline void set_d(const float d){ m_d=d; }
+
+
+ inline float get_S0() const{ return m_S0;}
+ inline void set_S0(const float S0){ m_S0=S0; }
+
inline bool compute_d_prior(){
m_d_old_prior=m_d_prior;
@@ -530,6 +522,10 @@ namespace FIBRE{
m_d_rej=0;
m_S0_acc=0;
m_S0_rej=0;
+ for(unsigned int f=0; f<m_fibres.size();f++){
+ m_fibres[f].update_proposals();
+ }
+
}
diff --git a/xfibres.cc b/xfibres.cc
new file mode 100644
index 0000000..b7784e4
--- /dev/null
+++ b/xfibres.cc
@@ -0,0 +1,386 @@
+/* Xfibres Diffusion Partial Volume Model
+
+ Tim Behrens - FMRIB Image Analysis Group
+
+ Copyright (C) 2005 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <strstream>
+#define WANT_STREAM
+#define WANT_MATH
+// #include "newmatap.h"
+// #include "newmatio.h"
+#include <string>
+#include <math.h>
+#include "utils/log.h"
+#include "diff_pvmoptions.h"
+#include "utils/tracer_plus.h"
+#include "miscmaths/miscprob.h"
+#include "miscmaths/miscmaths.h"
+#include "newimage/newimageall.h"
+#include "stdlib.h"
+#include "fibre.h"
+#include "xfibresoptions.h"
+
+//#include "bint/model.h"
+//#include "bint/lsmcmcmanager.h"
+//#include "bint/lslaplacemanager.h"
+
+using namespace FIBRE;
+using namespace Xfibres;
+using namespace Utilities;
+using namespace NEWMAT;
+using namespace NEWIMAGE;
+using namespace MISCMATHS;
+
+
+
+const float maxfloat=1e10;
+const float minfloat=1e-10;
+const float maxlogfloat=23;
+const float minlogfloat=-23;
+
+
+inline float min(float a,float b){
+ return a<b ? a:b;}
+inline float max(float a,float b){
+ return a>b ? a:b;}
+inline Matrix Anis()
+{
+ Matrix A(3,3);
+ A << 1 << 0 << 0
+ << 0 << 0 << 0
+ << 0 << 0 << 0;
+ return A;
+}
+
+inline Matrix Is()
+{
+ Matrix I(3,3);
+ I << 1 << 0 << 0
+ << 0 << 1 << 0
+ << 0 << 0 << 1;
+ return I;
+}
+
+inline ColumnVector Cross(const ColumnVector& A,const ColumnVector& B)
+{
+ ColumnVector res(3);
+ res << A(2)*B(3)-A(3)*B(2)
+ << A(3)*B(1)-A(1)*B(3)
+ << A(1)*B(2)-B(1)*A(2);
+ return res;
+}
+
+inline Matrix Cross(const Matrix& A,const Matrix& B)
+{
+ Matrix res(3,1);
+ res << A(2,1)*B(3,1)-A(3,1)*B(2,1)
+ << A(3,1)*B(1,1)-A(1,1)*B(3,1)
+ << A(1,1)*B(2,1)-B(1,1)*A(2,1);
+ return res;
+}
+
+float mod(float a, float b){
+ while(a>b){a=a-b;}
+ while(a<0){a=a+b;}
+ return a;
+}
+
+
+Matrix form_Amat(const Matrix& r,const Matrix& b)
+{
+ Matrix A(r.Ncols(),7);
+ Matrix tmpvec(3,1), tmpmat;
+
+ for( int i = 1; i <= r.Ncols(); i++){
+ tmpvec << r(1,i) << r(2,i) << r(3,i);
+ tmpmat = tmpvec*tmpvec.t()*b(1,i);
+ A(i,1) = tmpmat(1,1);
+ A(i,2) = 2*tmpmat(1,2);
+ A(i,3) = 2*tmpmat(1,3);
+ A(i,4) = tmpmat(2,2);
+ A(i,5) = 2*tmpmat(2,3);
+ A(i,6) = tmpmat(3,3);
+ A(i,7) = 1;
+ }
+ return A;
+}
+
+inline SymmetricMatrix vec2tens(ColumnVector& Vec){
+ SymmetricMatrix tens(3);
+ tens(1,1)=Vec(1);
+ tens(2,1)=Vec(2);
+ tens(3,1)=Vec(3);
+ tens(2,2)=Vec(4);
+ tens(3,2)=Vec(5);
+ tens(3,3)=Vec(6);
+ return tens;
+}
+
+
+
+class Samples{
+ xfibresOptions& opts;
+ Matrix m_dsamples;
+ Matrix m_S0samples;
+ vector<Matrix> m_thsamples;
+ vector<Matrix> m_phsamples;
+ vector<Matrix> m_fsamples;
+ vector<Matrix> m_lamsamples;
+public:
+
+ Samples(int nvoxels):opts(xfibresOptions::getInstance()){
+ int count=0;
+ int nsamples=0;
+
+ for(int i=0;i<=opts.njumps.value();i++){
+ count++;
+ if(count==opts.sampleevery.value()){
+ count=0;nsamples++;
+ }
+ }
+
+ m_dsamples.ReSize(nsamples,nvoxels);
+ m_S0samples.ReSize(nsamples,nvoxels);
+ for(int f=0;f<opts.nfibres.value();f++){
+ m_thsamples[f].ReSize(nsamples,nvoxels);
+ m_thsamples[f]=0;
+ m_phsamples[f].ReSize(nsamples,nvoxels);
+ m_phsamples[f]=0;
+ m_fsamples[f].ReSize(nsamples,nvoxels);
+ m_fsamples[f]=0;
+ m_lamsamples[f].ReSize(nsamples,nvoxels);
+ m_lamsamples[f]=0;
+ }
+
+ }
+
+
+ void record(Multifibre& mfib, int vox, int samp){
+ m_dsamples(samp,vox)=mfib.get_d();
+ m_S0samples(samp,vox)=mfib.get_S0();
+ for(int f=0;f<opts.nfibres.value();f++){
+ m_thsamples[f](samp,vox)=mfib.fibres()[f].get_th();
+ m_phsamples[f](samp,vox)=mfib.fibres()[f].get_ph();
+ m_fsamples[f](samp,vox)=mfib.fibres()[f].get_f();
+ m_lamsamples[f](samp,vox)=mfib.fibres()[f].get_lam();
+
+ }
+ }
+
+ void save(){
+
+ }
+
+};
+
+
+
+
+
+
+
+
+
+
+
+class xfibresVoxelManager{
+
+ xfibresOptions& opts;
+
+ Samples& m_samples;
+ int m_voxelnumber;
+ const ColumnVector& m_data;
+ const ColumnVector& m_alpha;
+ const ColumnVector& m_beta;
+ const Matrix& m_bvals;
+ Multifibre m_multifibre;
+ public:
+ xfibresVoxelManager(const ColumnVector& data,const ColumnVector& alpha,
+ const ColumnVector& beta, const Matrix& b, const Matrix& Amat,
+ Samples& samples,int voxelnumber):
+ opts(xfibresOptions::getInstance()),
+ m_samples(samples),m_voxelnumber(voxelnumber),m_data(data),
+ m_alpha(alpha), m_beta(beta), m_bvals(b),
+ m_multifibre(m_data,m_alpha,m_beta,m_bvals,opts.nfibres.value()){ }
+
+
+
+ void initialise(const Matrix& Amat){
+ //initialising
+ ColumnVector logS(m_data.Nrows()),tmp(m_data.Nrows()),Dvec(7),dir(3);
+ SymmetricMatrix tens;
+ DiagonalMatrix Dd;
+ Matrix Vd;
+ float mDd,fsquared;
+ float th,ph,f,D,S0;
+
+ for ( int i = 1; i <= logS.Nrows(); i++)
+ {
+ if(m_data(i)>0){
+ logS(i)=log(m_data(i));
+ }
+ else{
+ logS(i)=0;
+ }
+ }
+ Dvec = -pinv(Amat)*logS;
+ if( Dvec(7) > -maxlogfloat ){
+ S0=exp(-Dvec(7));
+ }
+ else{
+ S0=m_data.MaximumAbsoluteValue();
+ }
+
+ for ( int i = 1; i <= logS.Nrows(); i++)
+ {
+ if(S0<m_data.Sum()/m_data.Nrows()){ S0=m_data.MaximumAbsoluteValue(); }
+ logS(i)=(m_data(i)/S0)>0.01 ? log((i)):log(0.01*S0);
+ }
+ Dvec = -pinv(Amat)*logS;
+ S0=exp(-Dvec(7));
+ if(S0<m_data.Sum()/m_data.Nrows()){ S0=m_data.Sum()/m_data.Nrows(); }
+ tens = vec2tens(Dvec);
+ EigenValues(tens,Dd,Vd);
+ mDd = Dd.Sum()/Dd.Nrows();
+ int maxind = Dd(1) > Dd(2) ? 1:2; //finding maximum eigenvalue
+ maxind = Dd(maxind) > Dd(3) ? maxind:3;
+ dir << Vd(1,maxind) << Vd(2,maxind) << Vd(3,maxind);
+ cart2sph(dir,th,ph);
+ th= mod(th,M_PI);
+ ph= mod(ph,2*M_PI);
+ D = Dd(maxind);
+
+
+ float numer=1.5*((Dd(1)-mDd)*(Dd(1)-mDd)+(Dd(2)-mDd)*(Dd(2)-mDd)+(Dd(3)-mDd)*(Dd(3)-mDd));
+ float denom=(Dd(1)*Dd(1)+Dd(2)*Dd(2)+Dd(3)*Dd(3));
+ if(denom>0) fsquared=numer/denom;
+ else fsquared=0;
+ if(fsquared>0){f=sqrt(fsquared);}
+ else{f=0;}
+ if(f>=0.95) f=0.95;
+ if(f<=0.001) f=0.001;
+
+ if(opts.nfibres.value()>0){
+ Fibre fib(m_alpha,m_beta,m_bvals,D,th,ph,f,1);
+ m_multifibre.addfibre(fib);
+ Fibre fibun(m_alpha,m_beta,m_bvals,D);
+ for(int i=2; i<=opts.nfibres.value(); i++){
+ m_multifibre.addfibre(fibun);
+ }
+
+ }
+
+ }
+
+
+ void runmcmc(){
+ int count=0, recordcount=0,sample=0;
+ for( int i =0;i<opts.nburn.value();i++){
+ m_multifibre.jump();
+ count++;
+ if(count==opts.updateproposalevery.value()){
+ m_multifibre.update_proposals();
+ count=0;
+ }
+ }
+
+ for( int i =0;i<opts.njumps.value();i++){
+ m_multifibre.jump();
+ count++;
+ recordcount++;
+ if(recordcount==opts.sampleevery.value()){
+ m_samples.record(m_multifibre,m_voxelnumber,sample);
+ sample++;
+ recordcount=0;
+ }
+ if(count==opts.updateproposalevery.value()){
+ m_multifibre.update_proposals();
+ count=0;
+
+ }
+ }
+
+
+ }
+
+};
+
+
+
+int main(int argc, char *argv[])
+{
+ try{
+
+ // Setup logging:
+ Log& logger = LogSingleton::getInstance();
+
+ xfibresOptions& opts = xfibresOptions::getInstance();
+ opts.parse_command_line(argc,argv,logger);
+
+ srand(xfibresOptions::getInstance().seed.value());
+
+
+ Matrix datam, bvals,bvecs;
+ volume<char> mask;
+ bvals=read_ascii_matrix(opts.bvalsfile.value());
+ bvecs=read_ascii_matrix(opts.bvecsfile.value());
+
+ {//scope in whic the data exists in 4D format;
+ volume4D<float> data;
+ read_volume4D(data,opts.datafile.value());
+ read_volume(mask,opts.maskfile.value());
+ datam=data.matrix(mask);
+ }
+
+ Samples samples(datam.Ncols());
+
+
+
+ //if(opts.debuglevel.value()==1)
+ //Tracer_Plus::setrunningstackon();
+
+ //if(opts.timingon.value())
+ //Tracer_Plus::settimingon();
+
+ // read data
+
+ //VolumeSeries data;
+ //data.read(opts.datafile.value());
+ // data.writeAsFloat(LogSingleton::getInstance().appendDir("data"));
+// cout<<"done"<<endl;
+// return 0;
+ //int ntpts = data.tsize();
+ //Matrix bvecs = read_ascii_matrix(opts.bvecsfile.value());
+ //Matrix bvals = read_ascii_matrix(opts.bvalsfile.value());
+ // mask:
+ //Volume mask;
+ ///mask.read(opts.maskfile.value());
+ //mask.threshold(1e-16);
+
+ // threshold using mask:
+ //data.setPreThresholdPositions(mask.getPreThresholdPositions());
+ //data.thresholdSeries();
+ ColumnVector A,B,C;
+
+ Matrix b,Amat;
+ int n;
+
+ Multifibre(A,B,C,b,n);
+ }
+ catch(Exception& e)
+ {
+ cerr << endl << e.what() << endl;
+ }
+ catch(X_OptionError& e)
+ {
+ cerr << endl << e.what() << endl;
+ }
+
+ return 0;
+}
diff --git a/xfibresoptions.cc b/xfibresoptions.cc
new file mode 100644
index 0000000..42a0a4e
--- /dev/null
+++ b/xfibresoptions.cc
@@ -0,0 +1,56 @@
+/* xfibresoptions.cc
+
+ TIM Behrens - FMRIB Image Analysis Group
+
+ Copyright (C) 2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#define WANT_STREAM
+#define WANT_MATH
+
+#include <iostream.h>
+#include <fstream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "xfibresoptions.h"
+#include "utils/log.h"
+#include "utils/tracer_plus.h"
+
+using namespace Utilities;
+
+namespace Xfibres {
+
+ xfibresOptions* xfibresOptions::gopt = NULL;
+
+void xfibresOptions::parse_command_line(int argc, char** argv, Log& logger)
+{
+ Tracer_Plus("XfibresOptions::parse_command_line");
+
+ // do once to establish log directory name
+ for(int a = options.parse_command_line(argc, argv); a < argc; a++);
+
+ if(help.value() || ! options.check_compulsory_arguments())
+ {
+ options.usage();
+ //throw Exception("Not all of the compulsory arguments have been provided");
+ exit(2);
+ }
+ else
+ {
+ // setup logger directory
+ if(forcedir.value())
+ logger.setthenmakeDir(logdir.value());
+ else
+ logger.makeDir(logdir.value());
+
+ cout << "Log directory is: " << logger.getDir() << endl;
+
+ // do again so that options are logged
+ for(int a = 0; a < argc; a++)
+ logger.str() << argv[a] << " ";
+ logger.str() << endl << "---------------------------------------------" << endl << endl;
+
+ }
+}
+}
diff --git a/xfibresoptions.h b/xfibresoptions.h
new file mode 100644
index 0000000..57ec3e7
--- /dev/null
+++ b/xfibresoptions.h
@@ -0,0 +1,146 @@
+/* BpmOptions.h
+
+ Mark Woolrich, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2000 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#if !defined(xfibresOptions_h)
+#define xfibresOptions_h
+
+#include <string>
+#include <iostream.h>
+#include <fstream.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include "utils/options.h"
+#include "utils/log.h"
+#include "utils/tracer_plus.h"
+//#include "newmatall.h"
+using namespace Utilities;
+
+namespace Xfibres {
+
+class xfibresOptions {
+ public:
+ static xfibresOptions& getInstance();
+ ~xfibresOptions() { delete gopt; }
+
+ Option<bool> verbose;
+ Option<bool> help;
+ Option<string> logdir;
+ Option<bool> forcedir;
+ Option<string> datafile;
+ Option<string> ofile;
+ Option<string> maskfile;
+ Option<string> bvecsfile;
+ Option<string> bvalsfile;
+ Option<int> nfibres;
+ Option<int> njumps;
+ Option<int> nburn;
+ Option<int> sampleevery;
+ Option<int> updateproposalevery;
+ Option<float> seed;
+ void parse_command_line(int argc, char** argv, Log& logger);
+
+ private:
+ xfibresOptions();
+ const xfibresOptions& operator=(xfibresOptions&);
+ xfibresOptions(xfibresOptions&);
+
+ OptionParser options;
+
+ static xfibresOptions* gopt;
+
+};
+
+ inline xfibresOptions& xfibresOptions::getInstance(){
+ if(gopt == NULL)
+ gopt = new xfibresOptions();
+
+ return *gopt;
+ }
+
+ inline xfibresOptions::xfibresOptions() :
+ verbose(string("-V,--verbose"), false,
+ string("switch on diagnostic messages"),
+ false, no_argument),
+ help(string("-h,--help"), false,
+ string("display this message"),
+ false, no_argument),
+ logdir(string("--ld,--logdir"), string("logdir"),
+ string("log directory (default is logdir)"),
+ false, requires_argument),
+ forcedir(string("--forcedir"),false,string("Use the actual directory name given - i.e. don't add + to make a new directory"),false,no_argument),
+ datafile(string("--data,--datafile"), string("data"),
+ string("data file"),
+ true, requires_argument),
+ ofile(string("-o,--out"), string("dti"),
+ string("Output basename"),
+ true, requires_argument),
+ maskfile(string("-m,--mask, --maskfile"), string("nodif_brain_mask"),
+ string("mask file"),
+ true, requires_argument),
+ bvecsfile(string("-r,--bvecs"), string("bvecs"),
+ string("b vectors file"),
+ true, requires_argument),
+ bvalsfile(string("-b,--bvals"), string("bvals"),
+ string("b values file"),
+ true, requires_argument),
+ nfibres(string("--nf,--nfibres"),1,
+ string("Maximum nukmber of fibres to fit in each voxel (default 1)"),
+ false,requires_argument),
+ njumps(string("--nj,--njumps"),5000,
+ string("Num of jumps to be made by MCMC (default is 5000)"),
+ false,requires_argument),
+ nburn(string("--bi,--burnin"),1,
+ string("Num of jumps at start of MCMC to be discarded (default is 1)"),
+ false,requires_argument),
+ sampleevery(string("--se,--sampleevery"),1,
+ string("Num of jumps for each sample (MCMC) (default is 1)"),
+ false,requires_argument),
+ updateproposalevery(string("--upe,--updateproposalevery"),40,
+ string("Num of jumps for each update to the proposal density std (MCMC) (default is 40)"),
+ false,requires_argument),
+ seed(string("--seed"),0.76986654,string("seed for pseudo random number generator"),
+ false,requires_argument),
+ options("xfibres", "xfibres -k <filename>\n xfibres --verbose\n")
+ {
+
+
+ try {
+ options.add(verbose);
+ options.add(help);
+ options.add(logdir);
+ options.add(forcedir);
+ options.add(datafile);
+ options.add(ofile);
+ options.add(maskfile);
+ options.add(bvecsfile);
+ options.add(bvalsfile);
+ options.add(nfibres);
+ options.add(njumps);
+ options.add(nburn);
+ options.add(sampleevery);
+ options.add(updateproposalevery);
+ options.add(seed);
+
+ }
+ catch(X_OptionError& e) {
+ options.usage();
+ cerr << endl << e.what() << endl;
+ }
+ catch(std::exception &e) {
+ cerr << e.what() << endl;
+ }
+
+ }
+}
+
+#endif
+
+
+
+
+
--
GitLab