From b255139e63f07daf1c28b97ac5b74a6bd1d63a3e Mon Sep 17 00:00:00 2001
From: Stamatios Sotiropoulos <stam@fmrib.ox.ac.uk>
Date: Thu, 17 May 2012 11:35:57 +0000
Subject: [PATCH] *** empty log message ***
---
rubix.cc | 645 +++++++++++++++++++++
rubix.h | 209 +++++++
rubixoptions.cc | 53 ++
rubixoptions.h | 195 +++++++
rubixvox.cc | 1480 +++++++++++++++++++++++++++++++++++++++++++++++
rubixvox.h | 729 +++++++++++++++++++++++
6 files changed, 3311 insertions(+)
create mode 100644 rubix.cc
create mode 100644 rubix.h
create mode 100644 rubixoptions.cc
create mode 100644 rubixoptions.h
create mode 100644 rubixvox.cc
create mode 100644 rubixvox.h
diff --git a/rubix.cc b/rubix.cc
new file mode 100644
index 0000000..c6a54bc
--- /dev/null
+++ b/rubix.cc
@@ -0,0 +1,645 @@
+/*
+ RubiX
+ Stamatios Sotiropoulos, Saad Jbabdi and Tim Behrens - FMRIB Image Analysis Group
+ Copyright (C) 2012 University of Oxford
+ CCOPYRIGHT
+*/
+
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#define WANT_STREAM
+#define WANT_MATH
+#include <string>
+#include <math.h>
+#include "utils/log.h"
+#include "utils/tracer_plus.h"
+#include "miscmaths/miscprob.h"
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/nonlin.h"
+#include "newimage/newimageall.h"
+#include "stdlib.h"
+#include "rubixoptions.h"
+#include "rubixvox.h"
+#include "rubix.h"
+#include "diffmodels.h"
+
+using namespace RUBIX;
+using namespace Utilities;
+using namespace NEWMAT;
+using namespace NEWIMAGE;
+using namespace MISCMATHS;
+
+
+////////////////////////////////////////////
+// HRSamples: MCMC SAMPLE STORAGE
+////////////////////////////////////////////
+
+//Store parameters for a certain sample at a certain HR voxel
+void HRSamples::record(const HRvoxel& HRv, int vox, int samp){
+ m_dsamples(samp,vox)=HRv.get_d();
+ m_S0samples(samp,vox)=HRv.get_S0();
+ if (m_rician)
+ m_tausamples(samp,vox)=HRv.get_tau();
+ if (m_modelnum==2)
+ m_d_stdsamples(samp,vox)=HRv.get_d_std();
+
+ for(int f=0; f<m_numfibres; f++){
+ m_thsamples[f](samp,vox)=HRv.fibres()[f].get_th();
+ m_phsamples[f](samp,vox)=HRv.fibres()[f].get_ph();
+ m_fsamples[f](samp,vox)=HRv.fibres()[f].get_f();
+ }
+}
+
+
+//Get the mean samples for a voxel once jumping has finished
+void HRSamples::finish_voxel(int vox){
+
+ m_mean_dsamples(vox)=m_dsamples.Column(vox).Sum()/m_nsamps;
+ m_mean_S0samples(vox)=m_S0samples.Column(vox).Sum()/m_nsamps;
+
+ if (m_rician)
+ m_mean_tausamples(vox)=m_tausamples.Column(vox).Sum()/m_nsamps;
+ if (m_modelnum==2)
+ m_mean_d_stdsamples(vox)=m_d_stdsamples.Column(vox).Sum()/m_nsamps;
+
+ for(int f=0; f<m_numfibres; f++){ //for each fibre of the voxel
+ m_mean_fsamples[f](vox)=m_fsamples[f].Column(vox).Sum()/m_nsamps;
+
+ SymmetricMatrix m_dyad(3); m_dyad=0;
+ ColumnVector m_vec(3);
+ DiagonalMatrix dyad_D; //eigenvalues
+ Matrix dyad_V; //eigenvectors
+
+ //Get the sum of dyadic tensors across samples
+ for (int n=1; n<=m_nsamps; n++){
+ float th=m_thsamples[f](n,vox);
+ float ph=m_phsamples[f](n,vox);
+ m_vec << sin(th)*cos(ph) << sin(th)*sin(ph)<<cos(th) ;
+ m_dyad << m_dyad+m_vec*m_vec.t();
+ }
+ m_dyad=m_dyad/m_nsamps;
+
+ //Eigendecompose the mean dyadic tensor to get the mean orientation across samples
+ EigenValues(m_dyad,dyad_D,dyad_V);
+ int maxeig;
+ if(dyad_D(1)>dyad_D(2)){
+ if(dyad_D(1)>dyad_D(3)) maxeig=1;
+ else maxeig=3;
+ }
+ else{
+ if(dyad_D(2)>dyad_D(3)) maxeig=2;
+ else maxeig=3;
+ }
+ m_dyadic_vectors[f](1,vox)=dyad_V(1,maxeig);
+ m_dyadic_vectors[f](2,vox)=dyad_V(2,maxeig);
+ m_dyadic_vectors[f](3,vox)=dyad_V(3,maxeig);
+ }
+}
+
+
+//save samples for all voxels
+void HRSamples::save(const volume<float>& mask){
+ volume4D<float> tmp;
+ //So that I can sort the output fibres into
+ // files ordered by fibre fractional volume..
+ vector<Matrix> thsamples_out=m_thsamples;
+ vector<Matrix> phsamples_out=m_phsamples;
+ vector<Matrix> fsamples_out=m_fsamples;
+
+ vector<Matrix> dyadic_vectors_out=m_dyadic_vectors;
+ vector<Matrix> mean_fsamples_out;
+ for(unsigned int f=0;f<m_mean_fsamples.size();f++)
+ mean_fsamples_out.push_back(m_mean_fsamples[f]);
+
+ Log& logger = LogSingleton::getInstance();
+ tmp.setmatrix(m_mean_dsamples,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),0);
+ save_volume(tmp[0],logger.appendDir("mean_dsamples"));
+
+ tmp.setmatrix(m_mean_S0samples,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),0);
+ save_volume(tmp[0],logger.appendDir("mean_S0samples"));
+
+ if (m_rician){
+ tmp.setmatrix(m_mean_tausamples,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),0);
+ save_volume(tmp[0],logger.appendDir("mean_tausamples"));
+ }
+
+ if (m_modelnum==2){
+ tmp.setmatrix(m_mean_d_stdsamples,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),0);
+ save_volume(tmp[0],logger.appendDir("mean_dstd_samples"));
+ }
+
+ //Sort the output based on mean_fsamples
+ vector<Matrix> sumf;
+ for(int f=0; f<m_numfibres; f++){
+ Matrix tmp=sum(m_fsamples[f],1);
+ sumf.push_back(tmp);
+ }
+ for(int vox=1;vox<=m_dsamples.Ncols();vox++){
+ vector<pair<float,int> > sfs;
+ pair<float,int> ftmp;
+
+ for(int f=0; f<m_numfibres; f++){
+ ftmp.first=sumf[f](1,vox);
+ ftmp.second=f;
+ sfs.push_back(ftmp);
+ }
+ sort(sfs.begin(),sfs.end());
+
+ for(int samp=1;samp<=m_dsamples.Nrows();samp++){
+ for(int f=0; f<m_numfibres; f++){;
+ thsamples_out[f](samp,vox)=m_thsamples[sfs[(sfs.size()-1)-f].second](samp,vox);
+ phsamples_out[f](samp,vox)=m_phsamples[sfs[(sfs.size()-1)-f].second](samp,vox);
+ fsamples_out[f](samp,vox)=m_fsamples[sfs[(sfs.size()-1)-f].second](samp,vox);
+ }
+ }
+
+ for(int f=0;f<m_numfibres;f++){
+ mean_fsamples_out[f](1,vox)=m_mean_fsamples[sfs[(sfs.size()-1)-f].second](vox);
+ dyadic_vectors_out[f](1,vox)=m_dyadic_vectors[sfs[(sfs.size()-1)-f].second](1,vox);
+ dyadic_vectors_out[f](2,vox)=m_dyadic_vectors[sfs[(sfs.size()-1)-f].second](2,vox);
+ dyadic_vectors_out[f](3,vox)=m_dyadic_vectors[sfs[(sfs.size()-1)-f].second](3,vox);
+ }
+ }
+
+ tmp.setmatrix(m_dsamples,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),tmp.min());
+ save_volume4D(tmp,logger.appendDir("d_samples"));
+
+ // save the sorted fibres
+ for(int f=0; f<m_numfibres; f++){
+ tmp.setmatrix(thsamples_out[f],mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),tmp.min());
+ string oname="th"+num2str(f+1)+"samples";
+ save_volume4D(tmp,logger.appendDir(oname));
+
+ tmp.setmatrix(phsamples_out[f],mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),tmp.min());
+ oname="ph"+num2str(f+1)+"samples";
+ save_volume4D(tmp,logger.appendDir(oname));
+
+ tmp.setmatrix(fsamples_out[f],mask);
+ tmp.setDisplayMaximumMinimum(1,0);
+ oname="f"+num2str(f+1)+"samples";
+ save_volume4D(tmp,logger.appendDir(oname));
+
+ tmp.setmatrix(mean_fsamples_out[f],mask);
+ tmp.setDisplayMaximumMinimum(1,0);
+ oname="mean_f"+num2str(f+1)+"samples";
+ save_volume(tmp[0],logger.appendDir(oname));
+
+ tmp.setmatrix(dyadic_vectors_out[f],mask);
+ tmp.setDisplayMaximumMinimum(1,-1);
+ oname="dyads"+num2str(f+1);
+ save_volume4D(tmp,logger.appendDir(oname));
+ }
+}
+
+
+
+
+////////////////////////////////////////////
+// LRSamples: MCMC SAMPLE STORAGE
+////////////////////////////////////////////
+//Store parameters for a certain sample at a certain LR voxel
+void LRSamples::record(const LRvoxel& LRv, int vox, int samp){
+ m_S0samples(samp,vox)=LRv.get_S0LR();
+
+ if (m_rician)
+ m_tauLRsamples(samp,vox)=LRv.get_tauLR();
+
+ m_lik_energy(samp,vox)=LRv.get_likelihood_energy();
+ m_prior_energy(samp,vox)=LRv.get_prior_energy();
+
+ for(int m=0; m<m_Nmodes; m++){
+ m_thsamples[m](samp,vox)=LRv.PModes()[m].get_th();
+ m_phsamples[m](samp,vox)=LRv.PModes()[m].get_ph();
+ m_ksamples[m](samp,vox)=1.0/LRv.PModes()[m].get_invkappa();
+ }
+}
+
+
+//Get the mean samples for a voxel once jumping has finished
+void LRSamples::finish_voxel(int vox){
+ m_mean_S0samples(vox)=m_S0samples.Column(vox).Sum()/m_nsamps;
+ if (m_rician)
+ m_mean_tausamples(vox)=m_tauLRsamples.Column(vox).Sum()/m_nsamps;
+
+ for(int m=0; m<m_Nmodes; m++){
+ m_mean_ksamples[m](vox)=m_ksamples[m].Column(vox).Sum()/m_nsamps;
+
+ DiagonalMatrix dyad_D; //eigenvalues
+ Matrix dyad_V; //eigenvectors
+ ColumnVector m_vec(3);
+ SymmetricMatrix m_dyad(3); m_dyad=0;
+
+ for (int n=1; n<=m_nsamps; n++){
+ float th=m_thsamples[m](n,vox);
+ float ph=m_phsamples[m](n,vox);
+ m_vec << sin(th)*cos(ph) << sin(th)*sin(ph)<<cos(th) ;
+ m_dyad << m_dyad+m_vec*m_vec.t();
+ }
+ m_dyad=m_dyad/m_nsamps;
+
+ //Eigendecompose the mean dyadic tensor to get the mean orientation across samples
+ EigenValues(m_dyad,dyad_D,dyad_V);
+ int maxeig;
+ if(dyad_D(1)>dyad_D(2)){
+ if(dyad_D(1)>dyad_D(3)) maxeig=1;
+ else maxeig=3;
+ }
+ else{
+ if(dyad_D(2)>dyad_D(3)) maxeig=2;
+ else maxeig=3;
+ }
+ m_dyadic_vectors[m](1,vox)=dyad_V(1,maxeig);
+ m_dyadic_vectors[m](2,vox)=dyad_V(2,maxeig);
+ m_dyadic_vectors[m](3,vox)=dyad_V(3,maxeig);
+ }
+}
+
+
+//save samples for all voxels
+void LRSamples::save(const volume<float>& mask){
+ volume4D<float> tmp;
+ //So that I can sort the output fibres into
+ // files ordered by dispersion..
+ vector<Matrix> thsamples_out=m_thsamples;
+ vector<Matrix> phsamples_out=m_phsamples;
+ vector<Matrix> ksamples_out=m_ksamples;
+
+ vector<Matrix> dyadic_vectors_out=m_dyadic_vectors;
+ vector<Matrix> mean_ksamples_out;
+ for(unsigned int f=0;f<m_mean_ksamples.size();f++)
+ mean_ksamples_out.push_back(m_mean_ksamples[f]);
+
+ Log& logger = LogSingleton::getInstance();
+ tmp.setmatrix(m_mean_S0samples,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),0);
+ save_volume(tmp[0],logger.appendDir("mean_S0LRsamples"));
+
+ tmp.setmatrix(m_lik_energy,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),0);
+ save_volume4D(tmp,logger.appendDir("En_Lik_samples"));
+
+ tmp.setmatrix(m_prior_energy,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),0);
+ save_volume4D(tmp,logger.appendDir("En_Prior_samples"));
+
+ if (m_rician){
+ tmp.setmatrix(m_mean_tausamples,mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),0);
+ save_volume(tmp[0],logger.appendDir("mean_tauLRsamples"));
+ }
+
+ //Sort the output based on mean_invksamples
+ vector<Matrix> sumk;
+ for(int f=0; f<m_Nmodes; f++){
+ Matrix tmp=sum(m_ksamples[f],1);
+ sumk.push_back(tmp);
+ }
+ for(int vox=1;vox<=m_S0samples.Ncols();vox++){
+ vector<pair<float,int> > sfs;
+ pair<float,int> ftmp;
+
+ for(int f=0; f<m_Nmodes; f++){
+ ftmp.first=sumk[f](1,vox);
+ ftmp.second=f;
+ sfs.push_back(ftmp);
+ }
+ sort(sfs.begin(),sfs.end());
+
+ for(int samp=1; samp<=m_nsamps; samp++){
+ for(int f=0; f<m_Nmodes; f++){;
+ thsamples_out[f](samp,vox)=m_thsamples[sfs[(sfs.size()-1)-f].second](samp,vox);
+ phsamples_out[f](samp,vox)=m_phsamples[sfs[(sfs.size()-1)-f].second](samp,vox);
+ ksamples_out[f](samp,vox)=m_ksamples[sfs[(sfs.size()-1)-f].second](samp,vox);
+ }
+ }
+
+ for(int f=0;f<m_Nmodes;f++){
+ mean_ksamples_out[f](1,vox)=m_mean_ksamples[sfs[(sfs.size()-1)-f].second](vox);
+ dyadic_vectors_out[f](1,vox)=m_dyadic_vectors[sfs[(sfs.size()-1)-f].second](1,vox);
+ dyadic_vectors_out[f](2,vox)=m_dyadic_vectors[sfs[(sfs.size()-1)-f].second](2,vox);
+ dyadic_vectors_out[f](3,vox)=m_dyadic_vectors[sfs[(sfs.size()-1)-f].second](3,vox);
+ }
+ }
+
+ // save the sorted fibres
+ for(int f=0; f<m_Nmodes; f++){
+ tmp.setmatrix(thsamples_out[f],mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),tmp.min());
+ string oname="Pth"+num2str(f+1)+"samples";
+ save_volume4D(tmp,logger.appendDir(oname));
+
+ tmp.setmatrix(phsamples_out[f],mask);
+ tmp.setDisplayMaximumMinimum(tmp.max(),tmp.min());
+ oname="Pph"+num2str(f+1)+"samples";
+ save_volume4D(tmp,logger.appendDir(oname));
+
+ tmp.setmatrix(ksamples_out[f],mask);
+ tmp.setDisplayMaximumMinimum(1,0);
+ oname="Pk"+num2str(f+1)+"samples";
+ save_volume4D(tmp,logger.appendDir(oname));
+
+ tmp.setmatrix(mean_ksamples_out[f],mask);
+ tmp.setDisplayMaximumMinimum(1,0);
+ oname="mean_Pk"+num2str(f+1)+"samples";
+ save_volume(tmp[0],logger.appendDir(oname));
+
+ tmp.setmatrix(dyadic_vectors_out[f],mask);
+ tmp.setDisplayMaximumMinimum(1,-1);
+ oname="Pdyads"+num2str(f+1);
+ save_volume4D(tmp,logger.appendDir(oname));
+ }
+}
+
+
+
+//////////////////////////////////////////////////////////////
+// LRVoxelManager:MCMC HANDLING for a single LR voxel
+/////////////////////////////////////////////////////////////
+
+
+void LRVoxelManager::initialise(){
+ float pvmS0, pvmd,pvmd_std;
+ ColumnVector pvmf,pvmth,pvmph,pvm2invk,pvm2th,pvm2ph,predicted_signal;
+
+ //Initialise each HR voxel using the HR data
+ float sumd=0, sumd2=0;
+ for (int n=0; n<m_HRvoxnumber.Nrows(); n++){
+ if (opts.modelnum.value()==1){ //Model 1
+
+ PVM_single_c pvm(m_dataHR[n],m_bvecsHR,m_bvalsHR,opts.nfibres.value());
+ pvm.fit(); // this will give th,ph,f in the correct order
+
+ pvmf = pvm.get_f(); pvmth = pvm.get_th(); pvmph = pvm.get_ph();
+ pvmS0 = fabs(pvm.get_s0()); pvmd = pvm.get_d(); predicted_signal=pvm.get_prediction();
+ if(pvmd<0 || pvmd>0.01) pvmd=2e-3;
+
+ // DTI dti1(m_dataHR[n],m_bvecsHR,m_bvalsHR);
+ //dti1.linfit();
+ //pvmS0=fabs(dti1.get_s0());
+
+ m_LRv.set_HRparams(n,pvmd,pvmS0,pvmth,pvmph,pvmf);
+ }
+ else{ //Model 2
+ PVM_multi pvm(m_dataHR[n],m_bvecsHR,m_bvalsHR,opts.nfibres.value());
+ pvm.fit();
+
+ pvmf = pvm.get_f(); pvmth = pvm.get_th(); pvmph = pvm.get_ph(); pvmd_std=pvm.get_d_std();
+ pvmS0 =fabs(pvm.get_s0()); pvmd = pvm.get_d(); predicted_signal=pvm.get_prediction();
+ OUT(pvmf);
+ if(pvmd<0 || pvmd>0.01) pvmd=2e-3;
+ if(pvmd_std<0 || pvmd_std>0.01) pvmd_std=pvmd/10;
+
+ // DTI dti1(m_dataHR[n],m_bvecsHR,m_bvalsHR);
+ //dti1.linfit();
+ //pvmS0=fabs(dti1.get_s0());
+
+ m_LRv.set_HRparams(n,pvmd,pvmd_std,pvmS0,pvmth,pvmph,pvmf);
+ }
+
+ if (opts.rician.value()){ //If using Rician Energy, initialize tau, using the variance of the initial fit residuals
+ ColumnVector residuals;
+ residuals=m_dataHR[n]-predicted_signal;
+ float tau=1.0/var(residuals).AsScalar();
+ m_LRv.set_tauHR(n,tau);
+ }
+ sumd+=pvmd;
+ sumd2+=pvmd*pvmd;
+ }
+
+ sumd/=m_HRvoxnumber.Nrows();
+ m_LRv.set_meand(sumd);
+ m_LRv.set_stdevd(sumd/100);//sqrt((sumd2-m_HRvoxnumber.Nrows()*sumd*sumd)/(m_HRvoxnumber.Nrows()-1.0)));
+
+ m_LRv.set_mean_fsum(0.6);
+ m_LRv.set_stdev_fsum(0.01);
+
+ //Initialise the orientation prior parameters using the LR data
+ if (opts.nmodes.value()>0){
+ PVM_single_c pvm2(m_dataLR,m_bvecsLR,m_bvalsLR,opts.nmodes.value());
+ pvm2.fit();
+
+ pvm2invk = pvm2.get_f();
+ pvm2th = pvm2.get_th();
+ pvm2ph = pvm2.get_ph();
+ for (int m=1; m<=pvm2th.Nrows(); m++)
+ pvm2invk(m)=0.02;
+
+ m_LRv.set_Priorparams(pvm2th,pvm2ph,pvm2invk);
+ }
+
+ //Initialise the rest LR params
+ //DTI dti2(m_dataLR,m_bvecsLR,m_bvalsLR);
+ //dti2.linfit();
+ PVM_single_c pvm3(m_dataLR,m_bvecsLR,m_bvalsLR,1);
+ pvm3.fit();
+ m_LRv.set_S0LR(fabs(pvm3.get_s0()));
+
+ if (opts.rician.value()){ //If using Rician Energy, initialize tau, using the variance of the initial fit residuals
+ ColumnVector residuals;
+ predicted_signal=pvm3.get_prediction();
+ residuals=m_dataLR-predicted_signal;
+ float tau=1.0/var(residuals).AsScalar();
+ m_LRv.set_tauLR(tau);
+ }
+
+ m_LRv.update_Orient_hyp_prior();
+ m_LRv.initialise_energies();
+}
+
+
+
+//Run MCMC for an LRvoxel
+void LRVoxelManager::runmcmc(){
+ int count=0, recordcount=0,sample=1;
+
+ for(int i=0; i<opts.nburn.value(); i++){ //Burn-In
+ m_LRv.jump();
+ count++;
+ if(count==opts.updateproposalevery.value()){
+ m_LRv.update_proposals();
+ count=0;
+ }
+ }
+
+ for( int i =0;i<opts.njumps.value();i++){ //Sampling
+ m_LRv.jump();
+ count++;
+ recordcount++;
+ if(recordcount==opts.sampleevery.value()){
+ for (int n=1; n<=m_HRvoxnumber.Nrows(); n++) //for each HR voxel
+ m_HRsamples.record(m_LRv.HRvoxels()[n-1],(int)m_HRvoxnumber(n),sample);
+ m_LRsamples.record(m_LRv,(int)m_LRvoxnumber,sample);
+
+ sample++;
+ recordcount=0;
+ }
+ if(count==opts.updateproposalevery.value()){
+ m_LRv.update_proposals();
+ count=0;
+ }
+ }
+
+ for (int n=1; n<=m_HRvoxnumber.Nrows(); n++) //for each HR voxel
+ m_HRsamples.finish_voxel((int)m_HRvoxnumber(n));
+ m_LRsamples.finish_voxel((int)m_LRvoxnumber);
+
+}
+
+
+//For a voxel (x,y,z) at Low_res, returns that overlapping HR voxels coordinates
+ReturnMatrix get_HRindices(const int x,const int y,const int z, const int xratio, const int yratio, const int zratio){
+ Matrix HRindices(xratio*yratio*zratio,3); //num_HR_voxels x 3 coordinates per voxel
+
+ int count=1;
+ for (int Hx=1; Hx<=xratio; Hx++)
+ for (int Hy=1; Hy<=yratio; Hy++)
+ for (int Hz=1; Hz<=zratio; Hz++){
+ HRindices(count,1)=(x+1)*xratio-Hx;
+ HRindices(count,2)=(y+1)*yratio-Hy;
+ HRindices(count,3)=(z+1)*zratio-Hz;
+ count++;
+ }
+
+ HRindices.Release();
+ return HRindices;
+}
+
+
+//Using a Low-Res mask, create an HR mask, after finding the corresponding HR voxels
+volume<float> createHR_mask(const volume<float>& maskLR, const volume4D<float>& dataHR){
+ volume<float> maskHR(dataHR.xsize(),dataHR.ysize(),dataHR.zsize());
+ copybasicproperties(dataHR,maskHR);
+
+ Matrix temp_indices;
+ maskHR=0;
+ float xratio=maskLR.xdim()/dataHR.xdim();
+ float yratio=maskLR.ydim()/dataHR.ydim();
+ float zratio=maskLR.zdim()/dataHR.zdim();
+
+ for (int i=0; i<maskLR.zsize(); i++)
+ for (int j=0; j<maskLR.ysize(); j++)
+ for (int k=0; k<maskLR.xsize(); k++)
+ if (maskLR(k,j,i)!=0){
+ temp_indices=get_HRindices(k,j,i,(int)xratio,(int)yratio,(int)zratio);
+ for (int n=1; n<=temp_indices.Nrows(); n++)
+ maskHR.value((int)temp_indices(n,1),(int)temp_indices(n,2),(int)temp_indices(n,3))=1;
+ }
+ return maskHR;
+}
+
+
+
+////////////////////////////////////////////
+// MAIN
+////////////////////////////////////////////
+
+int main(int argc, char *argv[])
+{
+ try{
+ // Setup logging:
+ Log& logger = LogSingleton::getInstance();
+ rubixOptions& opts = rubixOptions::getInstance();
+ opts.parse_command_line(argc,argv,logger);
+ srand(rubixOptions::getInstance().seed.value());
+ Matrix datamLR, bvalsLR,bvecsLR,matrix2volkeyLR;
+ Matrix datamHR, bvalsHR,bvecsHR,matrix2volkeyHR;
+ volume<float> maskLR, maskHR;
+ volume<int> vol2matrixkeyLR, vol2matrixkeyHR;
+
+ bvalsLR=read_ascii_matrix(opts.LRbvalsfile.value());
+ bvecsLR=read_ascii_matrix(opts.LRbvecsfile.value());
+ if(bvecsLR.Nrows()>3) bvecsLR=bvecsLR.t();
+ if(bvalsLR.Nrows()>1) bvalsLR=bvalsLR.t();
+ for(int i=1;i<=bvecsLR.Ncols();i++){
+ float tmpsum=sqrt(bvecsLR(1,i)*bvecsLR(1,i)+bvecsLR(2,i)*bvecsLR(2,i)+bvecsLR(3,i)*bvecsLR(3,i));
+ if(tmpsum!=0){
+ bvecsLR(1,i)=bvecsLR(1,i)/tmpsum;
+ bvecsLR(2,i)=bvecsLR(2,i)/tmpsum;
+ bvecsLR(3,i)=bvecsLR(3,i)/tmpsum;
+ }
+ }
+
+ bvalsHR=read_ascii_matrix(opts.HRbvalsfile.value());
+ bvecsHR=read_ascii_matrix(opts.HRbvecsfile.value());
+ if(bvecsHR.Nrows()>3) bvecsHR=bvecsHR.t();
+ if(bvalsHR.Nrows()>1) bvalsHR=bvalsHR.t();
+ for(int i=1;i<=bvecsHR.Ncols();i++){
+ float tmpsum=sqrt(bvecsHR(1,i)*bvecsHR(1,i)+bvecsHR(2,i)*bvecsHR(2,i)+bvecsHR(3,i)*bvecsHR(3,i));
+ if(tmpsum!=0){
+ bvecsHR(1,i)=bvecsHR(1,i)/tmpsum;
+ bvecsHR(2,i)=bvecsHR(2,i)/tmpsum;
+ bvecsHR(3,i)=bvecsHR(3,i)/tmpsum;
+ }
+ }
+
+ volume4D<float> dataLR,dataHR;
+ read_volume4D(dataLR,opts.LRdatafile.value());
+ read_volume(maskLR,opts.LRmaskfile.value());
+ datamLR=dataLR.matrix(maskLR);
+ matrix2volkeyLR=dataLR.matrix2volkey(maskLR);
+ vol2matrixkeyLR=dataLR.vol2matrixkey(maskLR);
+
+ read_volume4D(dataHR,opts.HRdatafile.value());
+ maskHR=createHR_mask(maskLR, dataHR);
+ save_volume(maskHR,logger.appendDir("HRbrain_mask")); //Generate and save an HR mask using the LR mask
+
+ datamHR=dataHR.matrix(maskHR);
+ matrix2volkeyHR=dataHR.matrix2volkey(maskHR);
+ vol2matrixkeyHR=dataHR.vol2matrixkey(maskHR);
+
+ float xratio=maskLR.xdim()/maskHR.xdim();
+ float yratio=maskLR.ydim()/maskHR.ydim();
+ float zratio=maskLR.zdim()/maskHR.zdim();
+
+ HRSamples HRsampl(datamHR.Ncols(), opts.njumps.value(), opts.sampleevery.value(), opts.nfibres.value(), opts.rician.value(), opts.modelnum.value());
+ LRSamples LRsampl(datamLR.Ncols(), opts.njumps.value(), opts.sampleevery.value(), opts.nmodes.value(), opts.rician.value());
+
+ //dHR.push_back(datamHR.Column(1)); dHR.push_back(datamHR.Column(2));
+ //dHR.push_back(datamHR.Column(3)); dHR.push_back(datamHR.Column(4));
+ //HRvoxnum<<1<<2<<3<<4; HRweights<<0.25<<0.25<<0.25<<0.25;
+
+ for(int vox=1;vox<=datamLR.Ncols();vox++){ //For each LR voxel
+ // if (vox==19){ //vox==6
+ ColumnVector dLR; Matrix HRindices;
+ dLR=datamLR.Column(vox); //Find the corresponding HR ones
+ HRindices=get_HRindices((int)matrix2volkeyLR(vox,1),(int)matrix2volkeyLR(vox,2),(int)matrix2volkeyLR(vox,3),(int)xratio,(int)yratio,(int)zratio);
+ //cout<<endl<<"S0LR: "<<dLR(1)<<endl<<endl; //Debugging code
+
+ ColumnVector HRvoxnum(HRindices.Nrows()), HRweights(HRindices.Nrows());
+ vector<ColumnVector> dHR;
+
+ for (int n=1; n<=HRindices.Nrows(); n++){
+ HRweights(n)=1.0/HRindices.Nrows();
+ HRvoxnum(n)=vol2matrixkeyHR((int)HRindices(n,1),(int)HRindices(n,2),(int)HRindices(n,3));
+ ColumnVector tmp_vec=datamHR.Column((int)HRvoxnum(n));
+ //cout<<(int)HRindices(n,1)<<" "<<(int)HRindices(n,2)<<" "<<(int)HRindices(n,3)<<" S0HR:"<<tmp_vec(1)<<endl;
+ dHR.push_back(datamHR.Column((int)HRvoxnum(n)));
+ }
+ cout <<vox<<"/"<<datamLR.Ncols()<<endl;
+ LRVoxelManager vm(HRsampl,LRsampl,vox,HRvoxnum, dLR,dHR, bvecsLR, bvalsLR, bvecsHR, bvalsHR, HRweights);
+ vm.initialise();
+ vm.runmcmc();
+ // }
+ }
+ HRsampl.save(maskHR);
+ LRsampl.save(maskLR);
+ }
+
+ catch(Exception& e){
+ cerr << endl << e.what() << endl;
+ }
+ catch(X_OptionError& e){
+ cerr << endl << e.what() << endl;
+ }
+
+ return 0;
+}
diff --git a/rubix.h b/rubix.h
new file mode 100644
index 0000000..a97af67
--- /dev/null
+++ b/rubix.h
@@ -0,0 +1,209 @@
+/* rubix.h: Classes utilized in RubiX MCMC storage and handling */
+/* Stamatios Sotiropoulos, FMRIB Analysis Group */
+/* Copyright (C) 2012 University of Oxford */
+/* CCOPYRIGHT */
+
+
+#if !defined(rubix_h)
+#define rubix_h
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#define WANT_STREAM
+#define WANT_MATH
+#include <string>
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/miscprob.h"
+#include "newimage/newimageall.h"
+#include "stdlib.h"
+
+using namespace NEWMAT;
+using namespace MISCMATHS;
+using namespace NEWIMAGE;
+
+namespace RUBIX{
+
+ ////////////////////////////////////////////
+ // MCMC SAMPLE STORAGE
+ ////////////////////////////////////////////
+ //Storing Samples for parameters of all HR voxels
+ class HRSamples{
+ Matrix m_dsamples; //Variables for storing MCMC samples of all voxels in the HRgrid
+ Matrix m_d_stdsamples;
+ Matrix m_S0samples;
+ Matrix m_tausamples;
+ vector<Matrix> m_thsamples;
+ vector<Matrix> m_phsamples;
+ vector<Matrix> m_fsamples;
+
+ //for storing means
+ RowVector m_mean_dsamples; //Storing mean_samples for all voxels in the HRgrid
+ RowVector m_mean_d_stdsamples;
+ RowVector m_mean_S0samples;
+ RowVector m_mean_tausamples;
+ vector<Matrix> m_dyadic_vectors;
+ vector<RowVector> m_mean_fsamples;
+
+ int m_nsamps;
+ const int m_njumps;
+ const int m_sample_every;
+ const int m_numfibres;
+ const bool m_rician;
+ const int m_modelnum;
+ //const string m_logdir;
+
+ public:
+ HRSamples(int nvoxels, const int njumps, const int sample_every, const int numfibres, const bool rician=false, const int modelnum=1):
+ m_njumps(njumps),m_sample_every(sample_every), m_numfibres(numfibres), m_rician(rician), m_modelnum(modelnum){
+ int count=0;
+ int nsamples=0;
+
+ for(int i=0;i<m_njumps; i++){
+ count++;
+ if(count==m_sample_every){
+ count=0;
+ nsamples++;
+ }
+ }
+ m_nsamps=nsamples;
+
+ m_dsamples.ReSize(nsamples,nvoxels); m_dsamples=0;
+ m_S0samples.ReSize(nsamples,nvoxels); m_S0samples=0;
+
+ m_mean_dsamples.ReSize(nvoxels); m_mean_dsamples=0;
+ m_mean_S0samples.ReSize(nvoxels); m_mean_S0samples=0;
+ if (m_rician){
+ m_tausamples.ReSize(nsamples,nvoxels); m_tausamples=0;
+ m_mean_tausamples.ReSize(nvoxels); m_mean_tausamples=0;
+ }
+ if (m_modelnum==2){
+ m_d_stdsamples.ReSize(nsamples,nvoxels); m_d_stdsamples=0;
+ m_mean_d_stdsamples.ReSize(nvoxels); m_mean_d_stdsamples=0;
+ }
+ Matrix tmpvecs(3,nvoxels); tmpvecs=0;
+ for(int f=0; f<m_numfibres; f++){
+ m_thsamples.push_back(m_S0samples);
+ m_phsamples.push_back(m_S0samples);
+ m_fsamples.push_back(m_S0samples);
+ m_dyadic_vectors.push_back(tmpvecs);
+ m_mean_fsamples.push_back(m_mean_S0samples);
+ }
+ }
+
+ ~HRSamples(){}
+
+ void record(const HRvoxel& HRv, int vox, int samp); //Store parameters for a certain sample at a certain HR voxel
+ void finish_voxel(int vox); //Get the mean samples for a voxel once jumping has finished
+ void save(const volume<float>& mask); //Save samples for all voxels
+ };
+
+
+
+
+ ////////////////////////////////////////////
+ // MCMC SAMPLE STORAGE
+ ////////////////////////////////////////////
+ //Storing Samples for parameters at the Low-Res level (e.g. priors, Low-res parameters)
+ class LRSamples{
+ vector<Matrix> m_thsamples; //Variables for storing MCMC samples of all voxels in the LRgrid
+ vector<Matrix> m_phsamples;
+ vector<Matrix> m_ksamples;
+ Matrix m_S0samples;
+ Matrix m_tauLRsamples;
+ Matrix m_lik_energy;
+ Matrix m_prior_energy;
+
+ //for storing means
+ vector<Matrix> m_dyadic_vectors;
+ vector<RowVector> m_mean_ksamples;
+ RowVector m_mean_S0samples;
+ RowVector m_mean_tausamples;
+
+ int m_nsamps;
+ const int m_njumps;
+ const int m_sample_every;
+ const int m_Nmodes;
+ const bool m_rician;
+ //const string m_logdir;
+
+ public:
+ LRSamples(int nvoxels, const int njumps, const int sample_every, const int Nmodes, const bool rician=false):
+ m_njumps(njumps),m_sample_every(sample_every), m_Nmodes(Nmodes), m_rician(rician){
+ int count=0;
+ int nsamples=0;
+
+ for(int i=0;i<m_njumps; i++){
+ count++;
+ if(count==m_sample_every){
+ count=0;
+ nsamples++;
+ }
+ }
+ m_nsamps=nsamples;
+
+ m_S0samples.ReSize(nsamples,nvoxels); m_S0samples=0;
+ m_lik_energy.ReSize(nsamples,nvoxels); m_lik_energy=0;
+ m_prior_energy.ReSize(nsamples,nvoxels); m_prior_energy=0;
+ m_mean_S0samples.ReSize(nvoxels); m_mean_S0samples=0;
+ if (m_rician){
+ m_tauLRsamples.ReSize(nsamples,nvoxels); m_tauLRsamples=0;
+ m_mean_tausamples.ReSize(nvoxels); m_mean_tausamples=0;
+ }
+ Matrix tmpvecs(3,nvoxels); tmpvecs=0;
+
+ for(int f=0; f<m_Nmodes; f++){
+ m_thsamples.push_back(m_S0samples);
+ m_phsamples.push_back(m_S0samples);
+ m_ksamples.push_back(m_S0samples);
+ m_dyadic_vectors.push_back(tmpvecs);
+ m_mean_ksamples.push_back(m_mean_S0samples);
+ }
+ }
+
+ ~LRSamples(){}
+
+ void record(const LRvoxel& LRv, int vox, int samp); //Store parameters for a certain sample at a certain LR voxel
+ void finish_voxel(int vox); //Get the mean samples for a voxel once jumping has finished
+ void save(const volume<float>& mask); //Save samples for all voxels
+ };
+
+
+
+ //////////////////////////////////////////////
+ // MCMC HANDLING for a single LR voxel
+ //////////////////////////////////////////////
+ class LRVoxelManager{
+ rubixOptions& opts;
+ HRSamples& m_HRsamples; //keep MCMC samples of the parameters of all voxels inferred at High-res grid
+ LRSamples& m_LRsamples; //keep MCMC samples of the parameters of all voxels inferred at Low-res grid
+ int m_LRvoxnumber;
+ ColumnVector m_HRvoxnumber;
+ LRvoxel m_LRv;
+ const ColumnVector& m_dataLR; //Low-Res Data for the specific LR voxel
+ const vector<ColumnVector>& m_dataHR; //High-Res Data for all contained HR voxels
+ const Matrix& m_bvecsLR; //bvecs at Low-Res (3 x LR_NumPoints)
+ const Matrix& m_bvalsLR; //bvalues at Low-Res (1 x HR_NumPoints)
+ const Matrix& m_bvecsHR; //bvecs at High-Res (3 x HR_NumPoints)
+ const Matrix& m_bvalsHR; //bvalues at High-Res (1 x HR_NumPoints)
+ const ColumnVector& m_HRweights; //Holds the volume fraction each HR voxel occupies out of the LR one
+ public:
+ //Constructor
+ LRVoxelManager(HRSamples& Hsamples, LRSamples& Lsamples, int LRvoxnum, ColumnVector& HRvoxnum,
+ const ColumnVector& dataLR,const vector<ColumnVector>& dataHR,
+ const Matrix& bvecsLR, const Matrix& bvalsLR, const Matrix& bvecsHR, const Matrix& bvalsHR, const ColumnVector& HRweights):
+ opts(rubixOptions::getInstance()), m_HRsamples(Hsamples), m_LRsamples(Lsamples), m_LRvoxnumber(LRvoxnum),m_HRvoxnumber(HRvoxnum),
+ m_LRv(bvecsHR, bvalsHR, bvecsLR, bvalsLR, dataLR, dataHR, opts.nfibres.value(), opts.nmodes.value(), HRweights, opts.modelnum.value(), opts.fudge.value(),opts.all_ard.value(), opts.no_ard.value(),opts.kappa_ard.value(), opts.fsumPrior.value(), opts. dPrior.value(), opts.rician.value()),
+ m_dataLR(dataLR), m_dataHR(dataHR),m_bvecsLR(bvecsLR), m_bvalsLR(bvalsLR), m_bvecsHR(bvecsHR), m_bvalsHR(bvalsHR), m_HRweights(HRweights) { }
+
+ ~LRVoxelManager() { }
+
+ void initialise(); //Initialise all parameters for an LR voxel
+ void runmcmc(); //Run MCMC for an LR voxel
+};
+
+
+}
+
+
+#endif
diff --git a/rubixoptions.cc b/rubixoptions.cc
new file mode 100644
index 0000000..c915773
--- /dev/null
+++ b/rubixoptions.cc
@@ -0,0 +1,53 @@
+/* rubixoptions.cc
+
+ Stam Sotiropoulos, FMRIB Image Analysis Group
+
+ Copyright (C) 2012 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#define WANT_STREAM
+#define WANT_MATH
+
+#include <iostream>
+#include <fstream>
+#include <stdlib.h>
+#include <stdio.h>
+#include "rubixoptions.h"
+#include "utils/log.h"
+#include "utils/tracer_plus.h"
+
+using namespace Utilities;
+
+namespace RUBIX {
+
+ rubixOptions* rubixOptions::gopt = NULL;
+
+ void rubixOptions::parse_command_line(int argc, char** argv, Log& logger){
+ Tracer_Plus("RubiXOptions::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/rubixoptions.h b/rubixoptions.h
new file mode 100644
index 0000000..cb8654c
--- /dev/null
+++ b/rubixoptions.h
@@ -0,0 +1,195 @@
+/* rubixoptions.h
+
+ Stam Sotiropoulos, FMRIB Image Analysis Group
+
+ Copyright (C) 2012 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#if !defined(rubixptions_h)
+#define rubixptions_h
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#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 RUBIX{
+
+ class rubixOptions {
+ public:
+ static rubixOptions& getInstance();
+ ~rubixOptions() { delete gopt; }
+
+ Option<bool> verbose;
+ Option<bool> help;
+ Option<string> logdir;
+ Option<bool> forcedir;
+ Option<string> LRdatafile;
+ Option<string> LRmaskfile;
+ Option<string> LRbvecsfile;
+ Option<string> LRbvalsfile;
+ Option<string> HRdatafile;
+ // Option<string> HRmaskfile;
+ Option<string> HRbvecsfile;
+ Option<string> HRbvalsfile;
+ Option<int> nfibres;
+ Option<int> nmodes; //number of modes in the orientation prior
+ Option<int> modelnum;
+ Option<float> fudge;
+ Option<int> njumps;
+ Option<int> nburn;
+ Option<int> sampleevery;
+ Option<int> updateproposalevery;
+ Option<int> seed;
+ Option<bool> no_ard;
+ Option<bool> all_ard;
+ Option<bool> kappa_ard;
+ Option<bool> fsumPrior;
+ Option<bool> dPrior;
+ Option<bool> rician;
+ void parse_command_line(int argc, char** argv, Log& logger);
+
+ private:
+ rubixOptions();
+ const rubixOptions& operator=(rubixOptions&);
+ rubixOptions(rubixOptions&);
+
+ OptionParser options;
+ static rubixOptions* gopt;
+ };
+
+ inline rubixOptions& rubixOptions::getInstance(){
+ if(gopt == NULL)
+ gopt = new rubixOptions();
+
+ return *gopt;
+ }
+
+ inline rubixOptions::rubixOptions() :
+ 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),
+ LRdatafile(string("--dLR"), string("data_LR"),
+ string("Low-Res data file"),
+ true, requires_argument),
+ LRmaskfile(string("--mLR"), string("nodif_brain_mask_LR"),
+ string("Low-Res mask file"),
+ true, requires_argument),
+ LRbvecsfile(string("--rLR"), string("bvecs_LR"),
+ string("Low-Res b vectors file"),
+ true, requires_argument),
+ LRbvalsfile(string("--bLR"), string("bvals_LR"),
+ string("Low-Res b values file"),
+ true, requires_argument),
+ HRdatafile(string("--dHR"), string("data_HR"),
+ string("High-Res data file"),
+ true, requires_argument),
+ // HRmaskfile(string("--mHR"), string("nodif_brain_mask_HR"),
+ // string("High-Res mask file"),
+ // true, requires_argument),
+ HRbvecsfile(string("--rHR"), string("bvecs_HR"),
+ string("High-Res b vectors file"),
+ true, requires_argument),
+ HRbvalsfile(string("--bHR"), string("bvals_HR"),
+ string("High-Res b values file"),
+ true, requires_argument),
+ nfibres(string("--nf,--nfibres"),1,
+ string("Maximum number of fibres to fit in each HR voxel (default 1)"),
+ false,requires_argument),
+ nmodes(string("--nM,--nmodes"),2,
+ string("Number of modes for the orientation prior (default 2)"),
+ false,requires_argument),
+ modelnum(string("--model"),1,
+ string("\tWhich model to use. 1=mono-exponential (default). 2=continous exponential"),
+ false,requires_argument),
+ fudge(string("--fudge"),1,
+ string("\tARD fudge factor"),
+ false,requires_argument),
+ njumps(string("--nj,--njumps"),1500,
+ string("Num of jumps to be made by MCMC (default is 1500)"),
+ false,requires_argument),
+ nburn(string("--bi,--burnin"),0,
+ string("Total num of jumps at start of MCMC to be discarded (default is 0)"),
+ false,requires_argument),
+ sampleevery(string("--se"),30,
+ string("\tNum of jumps for each sample (MCMC) (default is 30)"),
+ false,requires_argument),
+ updateproposalevery(string("--upe"),40,
+ string("\tNum of jumps for each update to the proposal density std (MCMC) (default is 40)"),
+ false,requires_argument),
+ seed(string("--seed"),8665904,string("\tseed for pseudo random number generator"),
+ false,requires_argument),
+ no_ard(string("--noard"),false,string("\tTurn ARD off on all fibres"),
+ false,no_argument),
+ all_ard(string("--allard"),false,string("Turn ARD on for all fibres (default: on for all but first)"),
+ false,no_argument),
+ kappa_ard(string("--ard_kappa"),false,string("Turn ARD on for the dispersion of all orientation prior modes (default: off)"),
+ false,no_argument),
+ fsumPrior(string("--fsumPrior"),false,string("Turn on prior for the fsums across an LR voxel (default: off)"),
+ false,no_argument),
+ dPrior(string("--dPrior"),false,string("Turn on prior for the diffusivity across an LR voxel (default: off)"),
+ false,no_argument),
+ rician(string("--rician"),false,string("Use Rician Noise model (default is Gaussian)"),
+ false,no_argument),
+ options("RubiX v1.0", "rubix --help (for list of options)\n")
+ {
+ try {
+ options.add(verbose);
+ options.add(help);
+ options.add(logdir);
+ options.add(forcedir);
+ options.add(LRdatafile);
+ options.add(LRmaskfile);
+ options.add(LRbvecsfile);
+ options.add(LRbvalsfile);
+ options.add(HRdatafile);
+ // options.add(HRmaskfile);
+ options.add(HRbvecsfile);
+ options.add(HRbvalsfile);
+ options.add(nfibres);
+ options.add(nmodes);
+ options.add(modelnum);
+ options.add(fudge);
+ options.add(njumps);
+ options.add(nburn);
+ options.add(sampleevery);
+ options.add(updateproposalevery);
+ options.add(seed);
+ options.add(no_ard);
+ options.add(all_ard);
+ options.add(kappa_ard);
+ options.add(fsumPrior);
+ options.add(dPrior);
+ options.add(rician);
+ }
+ catch(X_OptionError& e) {
+ options.usage();
+ cerr << endl << e.what() << endl;
+ }
+ catch(std::exception &e) {
+ cerr << e.what() << endl;
+ }
+
+ }
+}
+
+#endif
+
+
+
+
+
diff --git a/rubixvox.cc b/rubixvox.cc
new file mode 100644
index 0000000..0e85bcf
--- /dev/null
+++ b/rubixvox.cc
@@ -0,0 +1,1480 @@
+/* rubixvox.cc: Classes utilized in RubiX */
+/* Stamatios Sotiropoulos, FMRIB Analysis Group */
+/* Copyright (C) 2012 University of Oxford */
+/* CCOPYRIGHT */
+
+
+#include "rubixvox.h"
+
+namespace RUBIX{
+
+//Returns the natural log of the 0th order modified Bessel function of first kind for an argument x
+//Follows the exponential implementation of the Bessel function in Numerical Recipes, Ch. 6
+float logIo(const float x){
+ float y,b;
+
+ b=std::fabs(x);
+ if (b<3.75){
+ float a=x/3.75;
+ a*=a;
+ //Bessel function evaluation
+ y=1.0+a*(3.5156229+a*(3.0899424+a*(1.2067492+a*(0.2659732+a*(0.0360768+a*0.0045813)))));
+ y=std::log(y);
+ }
+ else{
+ float a=3.75/b;
+ //Bessel function evaluation
+ //y=(exp(b)/sqrt(b))*(0.39894228+a*(0.01328592+a*(0.00225319+a*(-0.00157565+a*(0.00916281+a*(-0.02057706+a*(0.02635537+a*(-0.01647633+a*0.00392377))))))));
+ //Logarithm of Bessel function
+ y=b+std::log((0.39894228+a*(0.01328592+a*(0.00225319+a*(-0.00157565+a*(0.00916281+a*(-0.02057706+a*(0.02635537+a*(-0.01647633+a*0.00392377))))))))/std::sqrt(b));
+ }
+ return y;
+}
+
+
+
+//////////////////////////////////////////////////////////////////////
+// RFibre: Models one anisotropic compartment in an HR voxel
+//////////////////////////////////////////////////////////////////////
+
+//Adapt the standard deviation of the proposal distributions during MCMC execution
+//to avoid over-rejection/over-acceptance of samples
+void RFibre::update_proposals(){
+ m_th_prop*=sqrt(float(m_th_acc+1)/float(m_th_rej+1));
+ m_th_prop=min(m_th_prop,maxfloat);
+ m_ph_prop*=sqrt(float(m_ph_acc+1)/float(m_ph_rej+1));
+ m_ph_prop=min(m_ph_prop,maxfloat);
+ m_f_prop*=sqrt(float(m_f_acc+1)/float(m_f_rej+1));
+ m_f_prop=min(m_f_prop,maxfloat);
+ m_th_acc=0;
+ m_th_rej=0;
+ m_ph_acc=0;
+ m_ph_rej=0;
+ m_f_acc=0;
+ m_f_rej=0;
+}
+
+
+//Call that after initializing the parameters
+void RFibre::initialise_energies(){
+ compute_th_ph_prior();
+ compute_f_prior();
+ compute_prior();
+ compute_signal();
+}
+
+
+//conditional prior for th and ph, use hyperparameters
+bool RFibre::compute_th_ph_prior(){
+ m_th_ph_old_prior=m_th_ph_prior;
+ m_th_ph_prior=0;
+ if (m_Orient_hyp_prior.Nrows()!=0){// && m_f>0.05){ //If an informative orientation prior is used
+ ColumnVector prior_vec(m_Orient_hyp_prior.Nrows());
+ for (int m=1; m<=m_Orient_hyp_prior.Nrows(); m++){ //for each mode of the prior
+ float ct=m_Orient_hyp_prior(m,5);
+ float dot=std::min(fabs(m_Orient_hyp_prior(m,1)*m_vec(1)+m_Orient_hyp_prior(m,2)*m_vec(2)+m_Orient_hyp_prior(m,3)*m_vec(3)),1.0);
+ prior_vec(m)=(1.0/m_Orient_hyp_prior(m,4))*dot*dot-ct; //This is the argument of the Watson (with the logged normalization constant incorporated)
+ //m_th_ph_prior=m_th_ph_prior+exp((1.0/m_Orient_hyp_prior(m,4))*dot*dot-ct); //compute the respective Watson pdf, this is numerically unstable
+ }
+ if (m_Orient_hyp_prior.Nrows()==1) //For one mode the prior energy is -log(exp(prior_vec(1))
+ m_th_ph_prior=-prior_vec(1);
+ else{ //For many modes, instead of calculating -log(Sum(exp(xi))) which is numerically unstable
+ float maxarg=prior_vec.Maximum(); //calculate -m-log(Sum(exp(xi-m))), with m=max(xi)
+ for (int m=1; m<=m_Orient_hyp_prior.Nrows(); m++) //For each mode of the prior
+ m_th_ph_prior=m_th_ph_prior+exp(prior_vec(m)-maxarg);
+ m_th_ph_prior=-maxarg-log(m_th_ph_prior); //Convert to -log Energy
+ }
+ //m_th_ph_prior=m_th_ph_prior-log(0.5*fabs(sin(m_th))); //Correct with the Jacobian of the transformation! We jump spherical coordinates instead of Cartesian!!
+ m_th_ph_prior=m_th_ph_prior-log(fabs(sin(m_th)/(4*M_PI))); //Correct with the Jacobian of the transformation! We jump spherical coordinates instead of Cartesian!!
+ return false; //set instant rejection flag to false
+ }
+ else{ //if no orientation prior imposed, use uniform prior
+ if (m_th==0)
+ m_th_ph_prior=0;
+ else
+ m_th_ph_prior=-log(0.5*fabs(sin(m_th)));
+ return false; //instant rejection flag
+ }
+}
+
+
+//ARD prior on f (if f_ard==true)
+bool RFibre::compute_f_prior(){
+ m_f_old_prior=m_f_prior;
+ if (m_f<=0 || m_f>=1 )
+ return true;
+ else{
+ if(!f_ard)
+ m_f_prior=0;
+ else
+ m_f_prior=std::log(m_f);
+ m_f_prior=m_ardfudge*m_f_prior;
+ return false;
+ }
+}
+
+
+//Compute Joint Prior
+void RFibre::compute_prior(){
+ m_old_prior_en=m_prior_en;
+ m_prior_en=m_th_ph_prior+m_f_prior;
+}
+
+
+//Used when orientation prior params are jumped
+void RFibre::restore_th_ph_prior(){
+ m_th_ph_prior=m_th_ph_old_prior;
+ m_prior_en=m_old_prior_en;
+}
+
+
+//Compute model predicted signal at High and Low Res, only due to the anisotropic compartment
+void RFibre::compute_signal(){
+ m_SignalHR_old=m_SignalHR;
+ m_SignalLR_old=m_SignalLR;
+
+ if(m_modelnum==1 || m_d_std<1e-5){
+ for (int i=1; i<=m_bvecsHR.Ncols(); i++){
+ float angtmp=m_vec(1)*m_bvecsHR(1,i)+m_vec(2)*m_bvecsHR(2,i)+m_vec(3)*m_bvecsHR(3,i);
+ angtmp=angtmp*angtmp;
+ m_SignalHR(i)=exp(-m_d*m_bvalsHR(1,i)*angtmp);
+ }
+
+ for (int i=1; i<=m_bvecsLR.Ncols(); i++){
+ float angtmp=m_vec(1)*m_bvecsLR(1,i)+m_vec(2)*m_bvecsLR(2,i)+m_vec(3)*m_bvecsLR(3,i);
+ angtmp=angtmp*angtmp;
+ m_SignalLR(i)=exp(-m_d*m_bvalsLR(1,i)*angtmp);
+ }
+ }
+ else if(m_modelnum==2){
+ float dbeta=m_d/(m_d_std*m_d_std);
+ float dalpha=m_d*dbeta;
+ for (int i=1; i<=m_bvecsHR.Ncols(); i++){
+ float angtmp=m_vec(1)*m_bvecsHR(1,i)+m_vec(2)*m_bvecsHR(2,i)+m_vec(3)*m_bvecsHR(3,i);
+ angtmp=angtmp*angtmp;
+ m_SignalHR(i)=exp(log(dbeta/(dbeta + m_bvalsHR(1,i)*angtmp))*dalpha);
+ }
+
+ for (int i=1; i<=m_bvecsLR.Ncols(); i++){
+ float angtmp=m_vec(1)*m_bvecsLR(1,i)+m_vec(2)*m_bvecsLR(2,i)+m_vec(3)*m_bvecsLR(3,i);
+ angtmp=angtmp*angtmp;
+ m_SignalLR(i)=exp(log(dbeta/(dbeta + m_bvalsLR(1,i)*angtmp))*dalpha);
+ }
+ }
+}
+
+
+
+bool RFibre::propose_th(){
+ m_th_old=m_th;
+ m_th+=normrnd().AsScalar()*m_th_prop;
+ m_vec_old=m_vec;
+ m_vec<<sin(m_th)*cos(m_ph) <<sin(m_th)*sin(m_ph) <<cos(m_th);
+ bool rejflag=compute_th_ph_prior();//inside this it stores the old prior
+ compute_prior();
+ compute_signal();
+ return rejflag;
+}
+
+
+void RFibre::reject_th(){
+ m_th=m_th_old;
+ m_th_ph_prior=m_th_ph_old_prior;
+ m_vec=m_vec_old;
+ m_prior_en=m_old_prior_en;
+ restoreSignals();
+ m_th_rej++;
+}
+
+
+bool RFibre::propose_ph(){
+ m_ph_old=m_ph;
+ m_ph+=normrnd().AsScalar()*m_ph_prop;
+ m_vec_old=m_vec;
+ m_vec<<sin(m_th)*cos(m_ph) <<sin(m_th)*sin(m_ph) <<cos(m_th);
+ bool rejflag=compute_th_ph_prior();//inside this it stores the old prior
+ compute_prior();
+ compute_signal();
+ return rejflag;
+}
+
+
+void RFibre::reject_ph(){
+ m_ph=m_ph_old;
+ m_th_ph_prior=m_th_ph_old_prior;
+ m_vec=m_vec_old;
+ m_prior_en=m_old_prior_en;
+ restoreSignals();
+ m_ph_rej++;
+}
+
+
+bool RFibre::propose_f(){
+ m_f_old=m_f;
+ m_f+=normrnd().AsScalar()*m_f_prop;
+ bool rejflag=compute_f_prior();
+ compute_prior();
+ return rejflag;
+}
+
+
+void RFibre::reject_f(){
+ m_f=m_f_old;
+ m_f_prior=m_f_old_prior;
+ m_prior_en=m_old_prior_en;
+ m_f_rej++;
+}
+
+
+void RFibre::report() const {
+ OUT(m_th); OUT(m_th_old); OUT(m_th_prop); OUT(m_th_acc); OUT(m_th_rej);
+ OUT(m_th_ph_prior); OUT(m_th_ph_old_prior);
+ OUT(m_ph); OUT(m_ph_old); OUT(m_ph_prop); OUT(m_ph_acc); OUT(m_ph_rej);
+ OUT(m_f); OUT(m_f_old); OUT(m_f_prop); OUT(m_f_acc); OUT(m_f_rej);
+ OUT(m_f_prior); OUT(m_f_old_prior);
+ OUT(m_prior_en); OUT(m_old_prior_en);
+}
+
+
+
+
+////////////////////////////////////////////////
+// HRvoxel
+////////////////////////////////////////////////
+
+//Initialize energies, signals and standard deviations for the proposal distributions
+//Call that after initializing the parameters
+void HRvoxel::initialise_energies_props(){
+ m_S0_prop=m_S0/10.0;
+ m_d_prop=m_d/10.0;
+ for (int f=0; f<m_numfibres; f++)
+ m_fibres[f].initialise_energies();
+ if (m_rician){
+ compute_tau_prior();
+ m_tau_prop=m_tau/2.0;
+ }
+ if (m_modelnum==2){
+ compute_d_std_prior();
+ m_d_std_prop=m_d_std/10.0;
+ }
+ compute_S0_prior();
+ compute_d_prior();
+ if (m_fsumPrior_ON)
+ compute_fsum_prior();
+ compute_prior();
+ compute_iso_signal();
+ compute_signal();
+}
+
+
+bool HRvoxel::compute_d_prior(){
+ m_d_old_prior=m_d_prior;
+ if(m_d<=0)
+ return true;
+ else{
+ //m_d_prior=0;
+
+ if (m_dPrior_ON)
+ m_d_prior=0.5*(m_d-m_mean_d)*(m_d-m_mean_d)/(m_stdev_d*m_stdev_d)+log(m_stdev_d); //Use a Gaussian centered around the neighbourhood meand
+ else{ //Use a Gamma Prior centered around 1e-3
+ float alpha=3.0; float beta=2000;
+ m_d_prior=(1.0-alpha)*log(m_d)+beta*m_d; //Gamma_prior: pow(m_d,alpha-1.0)*exp(-beta*m_d)
+ }
+ return false;
+ }
+}
+
+
+bool HRvoxel::compute_d_std_prior(){
+ m_d_std_old_prior=m_d_std_prior;
+ if(m_d_std<=0 || m_d_std>0.01)
+ return true;
+ else{
+ m_d_std_prior=0;
+ //m_d_std_prior=std::log(m_d_std);
+ return false;
+ }
+}
+
+
+bool HRvoxel::compute_S0_prior(){
+ m_S0_old_prior=m_S0_prior;
+ if(m_S0<0) return true;
+ else{
+ m_S0_prior=0;
+ return false;
+ }
+}
+
+
+bool HRvoxel::compute_tau_prior(){
+ m_tau_old_prior=m_tau_prior;
+ if(m_tau<=0) return true;
+ else{
+ m_tau_prior=0;
+ return false;
+ }
+}
+
+
+//Check if sum of volume fractions is >1
+bool HRvoxel::reject_f_sum(){
+ float fsum=0;//m_f0;
+ for(int f=0; f<m_numfibres; f++){
+ fsum+=m_fibres[f].get_f();
+ }
+ return fsum>1; //true if sum(f) > 1 and therefore, we should reject f
+}
+
+
+//Compute Joint Prior
+void HRvoxel::compute_prior(){
+ m_old_prior_en=m_prior_en;
+ m_prior_en=m_d_prior+m_S0_prior;
+ if (m_fsumPrior_ON)
+ m_prior_en=m_prior_en+m_fsum_prior;
+ if(m_rician)
+ m_prior_en=m_prior_en+m_tau_prior;
+ if(m_modelnum==2)
+ m_prior_en=m_prior_en+m_d_std_prior;
+ for(int f=0; f<m_numfibres; f++){
+ m_prior_en=m_prior_en+m_fibres[f].get_prior();
+ }
+}
+
+
+//Used to update the conditional orientation prior
+//when the prior parameters are jumped
+void HRvoxel::update_th_ph_prior(){
+ for(int f=0; f<m_numfibres; f++){
+ m_fibres[f].compute_th_ph_prior();
+ m_fibres[f].compute_prior();
+ }
+ compute_prior();
+}
+
+
+void HRvoxel::restore_th_ph_prior(){
+ for(int f=0; f<m_numfibres; f++)
+ m_fibres[f].restore_th_ph_prior();
+ m_prior_en=m_old_prior_en;
+}
+
+
+void HRvoxel::update_d_prior(){
+ compute_d_prior();
+ compute_prior();
+}
+
+
+void HRvoxel::restore_d_prior(){
+ m_d_prior=m_d_old_prior;
+ m_prior_en=m_old_prior_en;
+}
+
+
+void HRvoxel::compute_fsum_prior(){
+ m_fsum_old_prior=m_fsum_prior;
+ float fsum=0;
+ for(int f=0; f<m_numfibres; f++){
+ fsum+=m_fibres[f].get_f();
+ } //Gaussian centered around LRvox neighbourhood mean
+ m_fsum_prior=0.5*(fsum-m_mean_fsum)*(fsum-m_mean_fsum)/(m_stdev_fsum*m_stdev_fsum)+log(m_stdev_fsum);
+ //m_fsum_prior=0.5*(fsum-0.6)*(fsum-0.6)/(0.1*0.1);
+}
+
+
+void HRvoxel::update_fsum_prior(){
+ compute_fsum_prior();
+ compute_prior();
+}
+
+
+void HRvoxel::restore_fsum_prior(){
+ m_fsum_prior=m_fsum_old_prior;
+ m_prior_en=m_old_prior_en;
+}
+
+
+
+//Compute the predicted signal from the isotropic compartment only
+void HRvoxel::compute_iso_signal(){
+ m_iso_SignalHR_old=m_iso_SignalHR;
+ m_iso_SignalLR_old=m_iso_SignalLR;
+
+ if(m_modelnum==1 || m_d_std<1e-5){
+ for(int i=1; i<=m_bvecsHR.Ncols(); i++)
+ m_iso_SignalHR(i)=exp(-m_d*m_bvalsHR(1,i));
+ for(int i=1; i<=m_bvecsLR.Ncols(); i++)
+ m_iso_SignalLR(i)=exp(-m_d*m_bvalsLR(1,i));
+ }
+ else if (m_modelnum==2){
+ float dbeta=m_d/(m_d_std*m_d_std);
+ float dalpha=m_d*dbeta;
+
+ for(int i=1; i<=m_bvecsHR.Ncols(); i++)
+ m_iso_SignalHR(i)=exp(log(dbeta/(dbeta+m_bvalsHR(1,i)))*dalpha);
+ for(int i=1; i<=m_bvecsLR.Ncols(); i++)
+ m_iso_SignalLR(i)=exp(log(dbeta/(dbeta+m_bvalsLR(1,i)))*dalpha);
+ }
+}
+
+
+//Compute the total predicted signal at High and Low Res measurement points
+void HRvoxel::compute_signal(){
+ m_SignalHR_old=m_SignalHR;
+ m_SignalLR_old=m_SignalLR;
+ float fsum=0;//m_f0;
+ m_SignalHR=0; m_SignalLR=0;
+ for(int f=0; f<m_numfibres; f++){ //Signal from the anisotropic compartments
+ m_SignalHR=m_SignalHR+m_fibres[f].get_f()*m_fibres[f].getSignalHR();
+ m_SignalLR=m_SignalLR+m_fibres[f].get_f()*m_fibres[f].getSignalLR();
+ fsum+=m_fibres[f].get_f(); //Total anisotropic volume fraction
+ }
+
+ for(int i=1;i<=m_bvecsHR.Ncols();i++) //Add the signal from the isotropic compartment
+ m_SignalHR(i)=m_S0*(m_SignalHR(i)+(1-fsum)*m_iso_SignalHR(i));//+m_f0); //and multiply by S0 to get the total signal
+
+ for(int i=1;i<=m_bvecsLR.Ncols();i++) //Add the signal from the isotropic compartment
+ m_SignalLR(i)=m_S0*(m_SignalLR(i)+(1-fsum)*m_iso_SignalLR(i));//+m_f0); //and multiply by S0 to get the total signal
+}
+
+
+
+bool HRvoxel::propose_d(){
+ bool rejflag;
+ m_d_old=m_d;
+ m_d+=normrnd().AsScalar()*m_d_prop;
+ rejflag=compute_d_prior();
+ for(int f=0; f<m_numfibres; f++)
+ m_fibres[f].compute_signal();
+ compute_iso_signal();
+ return rejflag;
+}
+
+
+void HRvoxel::reject_d(){
+ m_d=m_d_old;
+ m_d_prior=m_d_old_prior;
+ for(int f=0; f<m_numfibres; f++)
+ m_fibres[f].restoreSignals();
+ m_iso_SignalHR=m_iso_SignalHR_old;
+ m_iso_SignalLR=m_iso_SignalLR_old;
+ m_d_rej++;
+}
+
+
+bool HRvoxel::propose_d_std(){
+ bool rejflag;
+ m_d_std_old=m_d_std;
+ m_d_std+=normrnd().AsScalar()*m_d_std_prop;
+ rejflag=compute_d_std_prior();
+ for(int f=0; f<m_numfibres; f++)
+ m_fibres[f].compute_signal();
+ compute_iso_signal();
+ return rejflag;
+}
+
+
+void HRvoxel::reject_d_std(){
+ m_d_std=m_d_std_old;
+ m_d_prior=m_d_old_prior;
+ for(int f=0; f<m_numfibres; f++)
+ m_fibres[f].restoreSignals();
+ m_iso_SignalHR=m_iso_SignalHR_old;
+ m_iso_SignalLR=m_iso_SignalLR_old;
+ m_d_std_rej++;
+}
+
+
+bool HRvoxel::propose_S0(){
+ m_S0_old=m_S0;
+ m_S0+=normrnd().AsScalar()*m_S0_prop;
+ bool rejflag=compute_S0_prior();//inside this it stores the old prior
+ return rejflag;
+}
+
+
+void HRvoxel::reject_S0(){
+ m_S0=m_S0_old;
+ m_S0_prior=m_S0_old_prior;
+ m_S0_rej++;
+}
+
+
+bool HRvoxel::propose_tau(){
+ m_tau_old=m_tau;
+ m_tau+=normrnd().AsScalar()*m_tau_prop;
+ bool rejflag=compute_tau_prior();//inside this it stores the old prior
+ return rejflag;
+}
+
+
+void HRvoxel::reject_tau(){
+ m_tau=m_tau_old;
+ m_tau_prior=m_tau_old_prior;
+ m_tau_rej++;
+}
+
+
+void HRvoxel::restore_prior_totsignal(){
+ m_prior_en=m_old_prior_en;
+ m_SignalLR=m_SignalLR_old;
+ m_SignalHR=m_SignalHR_old;
+}
+
+void HRvoxel::restore_prior(){
+ m_prior_en=m_old_prior_en;
+}
+
+
+//Adapt standard deviation of proposal distributions during MCMC execution
+//to avoid over-rejection/over-acceptance of MCMC samples
+void HRvoxel::update_proposals(){
+ m_d_prop*=sqrt(float(m_d_acc+1)/float(m_d_rej+1));
+ m_d_prop=min(m_d_prop,maxfloat);
+ m_S0_prop*=sqrt(float(m_S0_acc+1)/float(m_S0_rej+1));
+ m_S0_prop=min(m_S0_prop,maxfloat);
+ m_d_acc=0;
+ m_d_rej=0;
+ m_S0_acc=0;
+ m_S0_rej=0;
+ if (m_rician){
+ m_tau_prop*=sqrt(float(m_tau_acc+1)/float(m_tau_rej+1));
+ m_tau_prop=min(m_tau_prop,maxfloat);
+ m_tau_acc=0; m_tau_rej=0;
+ }
+ if(m_modelnum==2){
+ m_d_std_prop*=sqrt(float(m_d_std_acc+1)/float(m_d_std_rej+1));
+ m_d_std_prop=min(m_d_std_prop,maxfloat);
+ m_d_std_acc=0; m_d_std_rej=0;
+ }
+ for(unsigned int f=0; f<m_fibres.size();f++)
+ m_fibres[f].update_proposals();
+}
+
+
+
+void HRvoxel::report() const{
+ OUT(m_d); OUT(m_d_old); OUT(m_d_prop);
+ OUT(m_d_prior); OUT(m_d_old_prior); OUT(m_d_acc); OUT(m_d_rej);
+ OUT(m_d_std); OUT(m_d_std_old); OUT(m_d_std_prop);
+ OUT(m_d_std_prior); OUT(m_d_std_old_prior); OUT(m_d_std_acc); OUT(m_d_std_rej);
+ OUT(m_S0); OUT(m_S0_old); OUT(m_S0_prop);
+ OUT(m_S0_prior); OUT(m_S0_old_prior); OUT(m_S0_acc); OUT(m_S0_rej);
+ OUT(m_prior_en); OUT(m_old_prior_en); OUT(m_SignalHR.t()); OUT(m_SignalLR.t());
+ for (int i=0; i<m_numfibres; i++){
+ cout <<"fibre "<<i<<endl;
+ m_fibres[i].report();}
+}
+
+
+
+
+//Class that models a single mode of the orientation prior
+/////////////////////////////
+// Orient_Prior_Mode //
+/////////////////////////////
+
+//Call that after initialising the parameters
+void Orient_Prior_Mode::initialise_energies(){
+ compute_th_prior();
+ compute_ph_prior();
+ compute_invkappa_prior();
+ compute_prior();
+}
+
+
+//Uniform Prior on theta
+bool Orient_Prior_Mode::compute_th_prior(){
+ m_th_old_prior=m_th_prior;
+ if(m_th==0){ m_th_prior=0; }
+ else{
+ m_th_prior=-log(fabs(sin(m_th)/2.0));
+ }
+ return false; //instant rejection flag
+}
+
+
+//Uniform Prior on phi
+bool Orient_Prior_Mode::compute_ph_prior(){
+ m_ph_old_prior=m_ph_prior;
+ m_ph_prior=0;
+ return false;
+}
+
+
+//ARD Prior on invkappa
+bool Orient_Prior_Mode::compute_invkappa_prior(){
+ m_invkappa_old_prior=m_invkappa_prior;
+ if(m_invkappa<=0)
+ return true;
+ else{
+ if (m_kappa_ard)
+ m_invkappa_prior=log(m_invkappa);
+ else
+ m_invkappa_prior=0;
+ return false;
+ }
+}
+
+
+//Compute Joint Prior energy
+void Orient_Prior_Mode::compute_prior(){
+ m_old_prior_en=m_prior_en;
+ m_prior_en=m_th_prior+m_ph_prior+m_invkappa_prior;
+}
+
+
+//Filter out extreme values of invkappa to ensure calculations are numerically stable
+void Orient_Prior_Mode::filter_invkappa(){
+ const float MINinvk=1e-7; //1e-6
+ if (m_invkappa<MINinvk)
+ m_invkappa=MINinvk;
+}
+
+
+
+
+//Compute the normalization constant of a Watson distribution
+//with dispersion index 1/kappa. Use a saddlepoint approximation (see Kume & Wood, 2005)
+float Orient_Prior_Mode::compute_Watson_norm() {
+ float k,ct,R,t,Bm,K2,K3,K4,T;
+
+ k=1.0/m_invkappa;
+ R=sqrt(4.0*k*k+9.0-4.0*k);
+ t=-0.25*(R+3+2.0*k);
+ Bm=1.0/(R+3+2*k);
+ K2=1.0/(2.0*(k+t)*(k+t))+1.0/(t*t);
+ K3=-1.0/((k+t)*(k+t)*(k+t))-2.0/(t*t*t);
+ K4=3.0/((k+t)*(k+t)*(k+t)*(k+t))+6.0/(t*t*t*t);
+ T=K4/(8.0*K2*K2)-5.0*K3*K3/(24.0*K2*K2*K2);
+ //ct=4.0*M_PI*sqrt(Bm/R)*exp(-t+T);
+ ct=-t+T+log(4.0*M_PI*sqrt(Bm/R)); //compute the log of the normalization constant, to avoid exponential overflow
+ return ct;
+}
+
+
+bool Orient_Prior_Mode::propose_th(){
+ m_th_old=m_th;
+ m_th+=normrnd().AsScalar()*m_th_prop;
+ m_vec_old=m_vec;
+ m_vec<<sin(m_th)*cos(m_ph) <<sin(m_th)*sin(m_ph) <<cos(m_th);
+ bool rejflag=compute_th_prior();//inside this it stores the old prior
+ compute_prior();
+ return rejflag;
+}
+
+
+void Orient_Prior_Mode::reject_th(){
+ m_th=m_th_old;
+ m_th_prior=m_th_old_prior;
+ m_vec=m_vec_old;
+ m_prior_en=m_old_prior_en;
+ m_th_rej++;
+}
+
+
+bool Orient_Prior_Mode::propose_ph(){
+ m_ph_old=m_ph;
+ m_ph+=normrnd().AsScalar()*m_ph_prop;
+ m_vec_old=m_vec;
+ m_vec<<sin(m_th)*cos(m_ph) <<sin(m_th)*sin(m_ph) <<cos(m_th);
+ bool rejflag=compute_ph_prior();//inside this it stores the old prior
+ compute_prior();
+ return rejflag;
+}
+
+
+void Orient_Prior_Mode::reject_ph(){
+ m_ph=m_ph_old;
+ m_ph_prior=m_ph_old_prior;
+ m_vec=m_vec_old;
+ m_prior_en=m_old_prior_en;
+ m_ph_rej++;
+}
+
+
+bool Orient_Prior_Mode::propose_invkappa(){
+ m_invkappa_old=m_invkappa;
+ m_invkappa+=normrnd().AsScalar()*m_invkappa_prop;
+ filter_invkappa();
+ m_Watson_norm_old=m_Watson_norm;
+ m_Watson_norm=compute_Watson_norm();
+ bool rejflag=compute_invkappa_prior();
+ compute_prior();
+ return rejflag;
+}
+
+
+void Orient_Prior_Mode::reject_invkappa(){
+ m_invkappa=m_invkappa_old;
+ m_invkappa_prior=m_invkappa_old_prior;
+ m_Watson_norm=m_Watson_norm_old;
+ m_prior_en=m_old_prior_en;
+ m_invkappa_rej++;
+}
+
+
+void Orient_Prior_Mode::update_proposals(){
+ m_th_prop*=sqrt(float(m_th_acc+1)/float(m_th_rej+1));
+ m_th_prop=min(m_th_prop,maxfloat);
+ m_ph_prop*=sqrt(float(m_ph_acc+1)/float(m_ph_rej+1));
+ m_ph_prop=min(m_ph_prop,maxfloat);
+ m_invkappa_prop*=sqrt(float(m_invkappa_acc+1)/float(m_invkappa_rej+1));
+ m_invkappa_prop=min(m_invkappa_prop,maxfloat);
+ m_th_acc=0; m_th_rej=0;
+ m_ph_acc=0; m_ph_rej=0;
+ m_invkappa_acc=0; m_invkappa_rej=0;
+}
+
+
+
+
+
+////////////////////////////////////////////////
+// LRvoxel
+////////////////////////////////////////////////
+
+//Call that after initialising all parameters
+void LRvoxel::initialise_energies(){
+ m_S0LR_prop=m_S0LR/10.0;
+ compute_S0LR_prior();
+
+ if (m_dPrior_ON){
+ m_mean_d_prop=m_mean_d/10.0;
+ m_stdev_d_prop=m_stdev_d/10.0;
+ compute_meand_prior(); compute_stdevd_prior();
+ }
+ if (m_fsumPrior_ON){
+ compute_mean_fsum_prior();
+ compute_stdev_fsum_prior();
+ }
+
+ if (m_rician){
+ m_tauLR_prop=m_tauLR/2.0;
+ compute_tauLR_prior();
+ }
+ for (int m=0; m<m_Nmodes; m++)
+ m_PModes[m].initialise_energies();
+ for (unsigned int n=0; n<m_dataHR.size(); n++)
+ m_HRvoxels[n].initialise_energies_props();
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ // cout<<"Likelihood Energy:"<<m_likelihood_en<<endl;
+ //cout<<"Prior Energy:"<<m_prior_en<<endl;
+ //cout<<"Posterior Energy:"<<m_posterior_en<<endl<<endl;
+}
+
+
+//Update the matrix that keeps information
+//on the orientation prior parameters
+void LRvoxel::update_Orient_hyp_prior(){
+ for (int m=0; m<m_Nmodes; m++){
+ ColumnVector temp;
+ temp=m_PModes[m].getVec();
+ m_Orient_hyp_prior(m+1,1)=temp(1);
+ m_Orient_hyp_prior(m+1,2)=temp(2);
+ m_Orient_hyp_prior(m+1,3)=temp(3);
+ m_Orient_hyp_prior(m+1,4)=m_PModes[m].get_invkappa();
+ m_Orient_hyp_prior(m+1,5)=m_PModes[m].get_Watson_norm();
+ }
+}
+
+
+//Update the matrix that keeps information
+//on the orientation prior parameters. Update only
+//mode M (0<=M<N_modes)
+void LRvoxel::update_Orient_hyp_prior(int M){
+ ColumnVector temp;
+ temp=m_PModes[M].getVec();
+ m_Orient_hyp_prior(M+1,1)=temp(1);
+ m_Orient_hyp_prior(M+1,2)=temp(2);
+ m_Orient_hyp_prior(M+1,3)=temp(3);
+ m_Orient_hyp_prior(M+1,4)=m_PModes[M].get_invkappa();
+ m_Orient_hyp_prior(M+1,5)=m_PModes[M].get_Watson_norm();
+}
+
+
+
+//Set params for a single HRvoxel with index 0 <= n < N (model 1)
+void LRvoxel::set_HRparams(const int n, const float d, const float S0, const ColumnVector& th, const ColumnVector& ph, const ColumnVector& f){
+ m_HRvoxels[n].set_d(d);
+ m_HRvoxels[n].set_S0(S0);
+ if (m_allard && !m_Noard){
+ m_HRvoxels[n].addfibre(th(1),ph(1),f(1),true); //ARD on for the first fibre
+ }
+ else{
+ m_HRvoxels[n].addfibre(th(1),ph(1),f(1),false); //ARD off for the first fibre
+ }
+ for (int m=2; m<=m_numfibres; m++){
+ if (m_Noard){
+ m_HRvoxels[n].addfibre(th(m),ph(m),f(m),false); //ARD off for the other fibres
+ }
+ else{
+ m_HRvoxels[n].addfibre(th(m),ph(m),f(m),true); //ARD on for the other fibres
+ }
+ }
+}
+
+
+//Set params for a single HRvoxel with index 0 <= n < N (model 2)
+void LRvoxel::set_HRparams(const int n, const float d, const float d_std,const float S0, const ColumnVector& th, const ColumnVector& ph, const ColumnVector& f){
+ m_HRvoxels[n].set_d(d);
+ m_HRvoxels[n].set_d_std(d_std);
+ m_HRvoxels[n].set_S0(S0);
+ if (m_allard && !m_Noard)
+ m_HRvoxels[n].addfibre(th(1),ph(1),f(1),true); //ARD on for the first fibre
+ else
+ m_HRvoxels[n].addfibre(th(1),ph(1),f(1),false); //ARD off for the first fibre
+ for (int m=2; m<=m_numfibres; m++){
+ if (m_Noard)
+ m_HRvoxels[n].addfibre(th(m),ph(m),f(m),false); //ARD off for the other fibres
+ else
+ m_HRvoxels[n].addfibre(th(m),ph(m),f(m),true); //ARD on for the other fibres
+ }
+}
+
+
+
+//Set params for all modes of orientation priors
+void LRvoxel::set_Priorparams(const ColumnVector& th, const ColumnVector& ph, const ColumnVector& invkappa){
+ for (int m=1; m<=m_Nmodes; m++){
+ m_PModes[m-1].set_th(th(m));
+ m_PModes[m-1].set_ph(ph(m));
+ m_PModes[m-1].set_invkappa(invkappa(m));
+ }
+}
+
+
+
+bool LRvoxel::propose_S0LR(){
+ m_S0LR_old=m_S0LR;
+ m_S0LR+=normrnd().AsScalar()*m_S0LR_prop;
+ bool rejflag=compute_S0LR_prior();//inside this it stores the old prior
+ return rejflag;
+}
+
+
+void LRvoxel::reject_S0LR(){
+ m_S0LR=m_S0LR_old;
+ m_S0LR_prior=m_S0LR_old_prior;
+ m_S0LR_rej++;
+}
+
+
+bool LRvoxel::compute_S0LR_prior(){
+ m_S0LR_old_prior=m_S0LR_prior;
+ if(m_S0LR<0) return true;
+ else{
+ m_S0LR_prior=0;
+ return false;
+ }
+}
+
+
+
+bool LRvoxel::propose_meand(){
+ m_mean_d_old=m_mean_d;
+ m_mean_d+=normrnd().AsScalar()*m_mean_d_prop;
+ bool rejflag=compute_meand_prior();//inside this it stores the old prior
+ return rejflag;
+}
+
+void LRvoxel::reject_meand(){
+ m_mean_d=m_mean_d_old;
+ m_mean_d_prior=m_mean_d_old_prior;
+ m_mean_d_rej++;
+}
+
+
+bool LRvoxel::compute_meand_prior(){
+ m_mean_d_old_prior=m_mean_d_prior;
+ if(m_mean_d<=0) return true;
+ else{
+ //m_mean_d_prior=0;
+ float alpha=3.0;
+ float beta=2000; //Gamma_prior centered around 1e-3
+ m_mean_d_prior=(1.0-alpha)*log(m_mean_d)+beta*m_mean_d; //Gamma_prior: pow(m_mean_d,alpha-1.0)*exp(-beta*m_mean_d)
+
+ return false;
+ }
+}
+
+
+bool LRvoxel::propose_stdevd(){
+ m_stdev_d_old=m_stdev_d;
+ m_stdev_d+=normrnd().AsScalar()*m_stdev_d_prop;
+ bool rejflag=compute_stdevd_prior();//inside this it stores the old prior
+ return rejflag;
+}
+
+
+void LRvoxel::reject_stdevd(){
+ m_stdev_d=m_stdev_d_old;
+ m_stdev_d_prior=m_stdev_d_old_prior;
+ m_stdev_d_rej++;
+}
+
+bool LRvoxel::compute_stdevd_prior(){
+ m_stdev_d_old_prior=m_stdev_d_prior;
+ if(m_stdev_d<=0 || m_stdev_d>=0.002) return true;
+ else{
+ m_stdev_d_prior=0;
+ //m_stdev_d_prior=log(m_stdev_d); //ARD prior
+ return false;
+ }
+}
+
+
+bool LRvoxel::propose_mean_fsum(){
+ m_mean_fsum_old=m_mean_fsum;
+ m_mean_fsum+=normrnd().AsScalar()*m_mean_fsum_prop;
+ bool rejflag=compute_mean_fsum_prior();//inside this it stores the old prior
+ return rejflag;
+}
+
+
+void LRvoxel::reject_mean_fsum(){
+ m_mean_fsum=m_mean_fsum_old;
+ m_mean_fsum_prior=m_mean_fsum_old_prior;
+ m_mean_fsum_rej++;
+}
+
+
+bool LRvoxel::compute_mean_fsum_prior(){
+ m_mean_fsum_old_prior=m_mean_fsum_prior;
+ if(m_mean_fsum<=0 || m_mean_fsum>=1) return true;
+ else{
+ //m_mean_fsum_prior=0;
+ m_mean_fsum_prior=-log(m_mean_fsum)-log(1-m_mean_fsum); //Beta distribution with a=b=2: fsum^(a-1)*(1-fsum)^(b-1)
+ return false;
+ }
+}
+
+
+bool LRvoxel::propose_stdev_fsum(){
+ m_stdev_fsum_old=m_stdev_fsum;
+ m_stdev_fsum+=normrnd().AsScalar()*m_stdev_fsum_prop;
+ bool rejflag=compute_stdev_fsum_prior();//inside this it stores the old prior
+ return rejflag;
+}
+
+
+void LRvoxel::reject_stdev_fsum(){
+ m_stdev_fsum=m_stdev_fsum_old;
+ m_stdev_fsum_prior=m_stdev_fsum_old_prior;
+ m_stdev_fsum_rej++;
+}
+
+
+bool LRvoxel::compute_stdev_fsum_prior(){
+ m_stdev_fsum_old_prior=m_stdev_fsum_prior;
+ if(m_stdev_fsum<=0 || m_stdev_fsum>=0.1) return true; //Without restriction to (0,0.1), very large values (e.g. 100's) are preferred for stdev, so that the fsum_prior becomes uniform
+ else{
+ m_stdev_fsum_prior=0;
+ //m_stdev_fsum_prior=log(m_stdev_fsum); //ARD prior
+ //m_stdev_fsum_prior=-log(25*m_stdev_fsum)-log(1-5*m_stdev_fsum); //Beta distribution with a=2,b=2, defined on [0,0.2].
+ //m_stdev_fsum_prior=-2*log(1-m_stdev_fsum); //Beta distribution with a=1, b=3. Similar to an ARD, but is naturally restricted to [0,1]
+ return false;
+ }
+}
+
+
+bool LRvoxel::propose_tauLR(){
+ m_tauLR_old=m_tauLR;
+ m_tauLR+=normrnd().AsScalar()*m_tauLR_prop;
+ bool rejflag=compute_tauLR_prior();//inside this it stores the old prior
+ return rejflag;
+}
+
+
+void LRvoxel::reject_tauLR(){
+ m_tauLR=m_tauLR_old;
+ m_tauLR_prior=m_tauLR_old_prior;
+ m_tauLR_rej++;
+}
+
+
+bool LRvoxel::compute_tauLR_prior(){
+ m_tauLR_old_prior=m_tauLR_prior;
+ if(m_tauLR<=0) return true;
+ else{
+ m_tauLR_prior=0;
+ return false;
+ }
+}
+
+
+void LRvoxel::compute_prior(){
+ m_old_prior_en=m_prior_en;
+ m_prior_en=0;
+
+ m_prior_en=m_S0LR_prior;
+ if (m_rician)
+ m_prior_en+=m_tauLR_prior;
+ for (int m=0; m<m_Nmodes; m++)
+ m_prior_en+=m_PModes[m].get_prior();
+
+ for (unsigned int m=0; m<m_dataHR.size(); m++)
+ m_prior_en+=m_HRvoxels[m].get_prior();
+
+ if (m_dPrior_ON)
+ m_prior_en+=m_mean_d_prior+m_stdev_d_prior;
+ if (m_fsumPrior_ON)
+ m_prior_en+=m_mean_fsum_prior+m_stdev_fsum_prior;
+}
+
+
+void LRvoxel::compute_likelihood(){
+ ColumnVector SLRpred, SHRpred, Diff;
+ SLRpred.ReSize(m_bvecsLR.Ncols());
+ SHRpred.ReSize(m_bvecsHR.Ncols()); SHRpred=0;
+ float likLR, likHR;
+
+ m_old_likelihood_en=m_likelihood_en;
+
+ if (!m_rician){ //Gaussian Likelihood Energy Calculation
+ //likelihood of Low-Res data
+ SLRpred=0;
+ for (unsigned int m=0; m<m_dataHR.size(); m++) //add attenuations
+ SLRpred+=m_HRweights(m+1)*m_HRvoxels[m].getSignalLR()/m_HRvoxels[m].get_S0();
+ SLRpred=m_S0LR*SLRpred;
+ SLRpred=m_dataLR-SLRpred;
+ likLR=0.5*m_bvecsLR.Ncols()*log(0.5*SLRpred.SumSquare());
+
+ //likelihood of High-Res data
+ likHR=0;
+ for (unsigned int m=0; m<m_dataHR.size(); m++){
+ SHRpred=m_HRvoxels[m].getSignalHR();
+ SHRpred=m_dataHR[m]-SHRpred;
+ likHR+=log(0.5*SHRpred.SumSquare());
+ }
+ likHR*=0.5*m_bvecsHR.Ncols();
+
+ m_likelihood_en=likLR+likHR;
+ }
+ else{ //Rician Likelihood Energy Calculation
+ //likelihood of Low-Res data
+ SLRpred=0;
+ for (unsigned int m=0; m<m_dataHR.size(); m++) //add attenuations
+ SLRpred+=m_HRweights(m+1)*m_HRvoxels[m].getSignalLR()/m_HRvoxels[m].get_S0();
+ SLRpred=m_S0LR*SLRpred;
+ likLR=-m_bvecsLR.Ncols()*log(m_tauLR);
+ for (int k=1; k<=m_bvecsLR.Ncols(); k++)
+ likLR-=m_logdataLR(k)-0.5*m_tauLR*(m_dataLR(k)*m_dataLR(k)+SLRpred(k)*SLRpred(k))+logIo(m_tauLR*m_dataLR(k)*SLRpred(k));
+
+ //likelihood of High-Res data
+ likHR=0;
+ for (unsigned int m=0; m<m_dataHR.size(); m++){
+ float m_tauHR=m_HRvoxels[m].get_tau();
+ SHRpred=m_HRvoxels[m].getSignalHR();
+ likHR+=m_bvecsHR.Ncols()*log(m_tauHR);
+ for (int k=1; k<=m_bvecsHR.Ncols(); k++)
+ likHR+=m_logdataHR[m](k)-0.5*m_tauHR*(m_dataHR[m](k)*m_dataHR[m](k)+SHRpred(k)*SHRpred(k))+logIo(m_tauHR*m_dataHR[m](k)*SHRpred(k));
+ }
+ m_likelihood_en=likLR-likHR;
+ }
+}
+
+
+void LRvoxel::compute_posterior(){
+ m_old_posterior_en=m_posterior_en;
+ m_posterior_en=m_prior_en+m_likelihood_en;
+}
+
+
+//Test whether the candidate state will be accepted
+bool LRvoxel::test_energy() const{
+ float tmp=exp(m_old_posterior_en-m_posterior_en);
+ return (tmp>unifrnd().AsScalar());
+}
+
+
+void LRvoxel::restore_energies(){
+ m_prior_en=m_old_prior_en;
+ m_likelihood_en=m_old_likelihood_en;
+ m_posterior_en=m_old_posterior_en;
+}
+
+
+void LRvoxel::restore_Prior_Posterior(){
+ m_prior_en=m_old_prior_en;
+ m_posterior_en=m_old_posterior_en;
+}
+
+
+void LRvoxel::update_proposals(){
+ for(int m=0; m<m_Nmodes; m++)
+ m_PModes[m].update_proposals();
+ for(unsigned int n=0; n<m_dataHR.size(); n++)
+ m_HRvoxels[n].update_proposals();
+
+ m_S0LR_prop*=sqrt(float(m_S0LR_acc+1)/float(m_S0LR_rej+1));
+ m_S0LR_prop=min(m_S0LR_prop,maxfloat);
+ m_S0LR_acc=0; m_S0LR_rej=0;
+
+ if (m_dPrior_ON){
+ m_mean_d_prop*=sqrt(float(m_mean_d_acc+1)/float(m_mean_d_rej+1));
+ m_mean_d_prop=min(m_mean_d_prop,maxfloat);
+ m_mean_d_acc=0; m_mean_d_rej=0;
+
+ m_stdev_d_prop*=sqrt(float(m_stdev_d_acc+1)/float(m_stdev_d_rej+1));
+ m_stdev_d_prop=min(m_stdev_d_prop,maxfloat);
+ m_stdev_d_acc=0; m_stdev_d_rej=0;
+ }
+
+ if (m_fsumPrior_ON){
+ m_mean_fsum_prop*=sqrt(float(m_mean_fsum_acc+1)/float(m_mean_fsum_rej+1));
+ m_mean_fsum_prop=min(m_mean_fsum_prop,maxfloat);
+ m_mean_fsum_acc=0; m_mean_fsum_rej=0;
+
+ m_stdev_fsum_prop*=sqrt(float(m_stdev_fsum_acc+1)/float(m_stdev_fsum_rej+1));
+ m_stdev_fsum_prop=min(m_stdev_fsum_prop,maxfloat);
+ m_stdev_fsum_acc=0; m_stdev_fsum_rej=0;
+ }
+
+ if (m_rician){
+ m_tauLR_prop*=sqrt(float(m_tauLR_acc+1)/float(m_tauLR_rej+1));
+ m_tauLR_prop=min(m_tauLR_prop,maxfloat);
+ m_tauLR_acc=0; m_tauLR_rej=0;
+ }
+}
+
+
+//Single MCMC iteration with all parameters jumped
+void LRvoxel::jump(){
+
+ //Jump LRvoxel parameters first
+
+ if(!propose_S0LR()){ //Try S0LR
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy()){
+ accept_S0LR();
+ }
+ else{
+ restore_energies();
+ reject_S0LR();
+ }
+ }
+ else
+ reject_S0LR();
+
+ if (m_dPrior_ON){
+ if(!propose_meand()){ //Try mean_d
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update nuisance_prior for each HR voxel
+ m_HRvoxels[n].update_d_prior();
+
+ compute_prior();
+ compute_posterior();
+ if(test_energy()){
+ accept_meand();
+ }
+ else{
+ restore_Prior_Posterior();
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //restore nuisance_prior_energy for each HR voxel
+ m_HRvoxels[n].restore_d_prior();
+ reject_meand();
+ }
+ }
+ else
+ reject_meand();
+
+ if(!propose_stdevd()){ //Try stdev_d
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update nuisance_prior energy for each HR voxel
+ m_HRvoxels[n].update_d_prior();
+
+ compute_prior();
+ compute_posterior();
+ if(test_energy()){
+ accept_stdevd();
+ }
+ else{
+ restore_Prior_Posterior();
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //restore nuisance_prior_energy for each HR voxel
+ m_HRvoxels[n].restore_d_prior();
+ reject_stdevd();
+ }
+ }
+ else
+ reject_stdevd();
+ }
+
+ if (m_fsumPrior_ON){
+ if(!propose_mean_fsum()){ //Try mean_fsum
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update nuisance_prior for each HR voxel
+ m_HRvoxels[n].update_fsum_prior();
+
+ compute_prior();
+ compute_posterior();
+ if(test_energy()){
+ accept_mean_fsum();
+ }
+ else{
+ restore_Prior_Posterior();
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //restore nuisance_prior_energy for each HR voxel
+ m_HRvoxels[n].restore_fsum_prior();
+ reject_mean_fsum();
+ }
+ }
+ else
+ reject_mean_fsum();
+
+ if(!propose_stdev_fsum()){ //Try stdev_fsum
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update nuisance_prior energy for each HR voxel
+ m_HRvoxels[n].update_fsum_prior();
+
+ compute_prior();
+ compute_posterior();
+ if(test_energy()){
+ accept_stdev_fsum();
+ }
+ else{
+ restore_Prior_Posterior();
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //restore nuisance_prior_energy for each HR voxel
+ m_HRvoxels[n].restore_fsum_prior();
+ reject_stdev_fsum();
+ }
+ }
+ else
+ reject_stdev_fsum();
+ }
+
+ if (m_rician){
+ if(!propose_tauLR()){ //Try tauLR
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy()){
+ accept_tauLR();
+ }
+ else{
+ restore_energies();
+ reject_tauLR();
+ }
+ }
+ else
+ reject_tauLR();
+ }
+
+ //For each mode of the Orientation Prior
+ for(int m=0; m<m_Nmodes; m++){
+
+ if(!m_PModes[m].propose_th()){ //Try theta
+ update_Orient_hyp_prior(m);
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update th_ph_prior for each HR voxel
+ m_HRvoxels[n].update_th_ph_prior();
+
+ compute_prior();
+ //compute_likelihood();
+ compute_posterior();
+ if(test_energy()){
+ m_PModes[m].accept_th();
+ }
+ else{
+ //restore_energies();
+ restore_Prior_Posterior();
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update th_ph_prior for each HR voxel
+ m_HRvoxels[n].restore_th_ph_prior();
+ m_PModes[m].reject_th();
+ update_Orient_hyp_prior(m);
+ }
+ }
+ else
+ m_PModes[m].reject_th();
+
+
+ if(!m_PModes[m].propose_ph()){ //Try phi
+ update_Orient_hyp_prior(m);
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update th_ph_prior for each HR voxel
+ m_HRvoxels[n].update_th_ph_prior();
+
+ compute_prior();
+ //compute_likelihood();
+ compute_posterior();
+ if(test_energy()){
+ m_PModes[m].accept_ph();
+ }
+ else{
+ //restore_energies();
+ restore_Prior_Posterior();
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update th_ph_prior for each HR voxel
+ m_HRvoxels[n].restore_th_ph_prior();
+ m_PModes[m].reject_ph();
+ update_Orient_hyp_prior(m);
+ }
+ }
+ else
+ m_PModes[m].reject_ph();
+
+
+ if(!m_PModes[m].propose_invkappa()){ //Try invkappa
+ update_Orient_hyp_prior(m);
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update th_ph_prior for each HR voxel
+ m_HRvoxels[n].update_th_ph_prior();
+
+ compute_prior();
+ //compute_likelihood();
+ compute_posterior();
+ if(test_energy()){
+ m_PModes[m].accept_invkappa();
+ }
+ else{
+ //restore_energies();
+ restore_Prior_Posterior();
+ for (unsigned int n=0; n<m_dataHR.size(); n++) //update th_ph_prior for each HR voxel
+ m_HRvoxels[n].restore_th_ph_prior();
+ m_PModes[m].reject_invkappa();
+ update_Orient_hyp_prior(m);
+ }
+ }
+ else
+ m_PModes[m].reject_invkappa();
+ } //end Orientation Prior params
+
+
+ //For each HRvoxel
+ for (unsigned int n=0; n<m_dataHR.size(); n++){
+
+ if(!m_HRvoxels[n].propose_d()){ //Try d
+ m_HRvoxels[n].compute_prior();
+ m_HRvoxels[n].compute_signal();
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy()){
+ m_HRvoxels[n].accept_d();
+ }
+ else{
+ restore_energies();
+ m_HRvoxels[n].restore_prior_totsignal();
+ m_HRvoxels[n].reject_d();
+ }
+ }
+ else
+ m_HRvoxels[n].reject_d();
+
+
+ if (m_modelnum==2){
+ if(!m_HRvoxels[n].propose_d_std()){ //Try d_std
+ m_HRvoxels[n].compute_prior();
+ m_HRvoxels[n].compute_signal();
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy()){
+ m_HRvoxels[n].accept_d_std();
+ }
+ else{
+ restore_energies();
+ m_HRvoxels[n].restore_prior_totsignal();
+ m_HRvoxels[n].reject_d_std();
+ }
+ }
+ else
+ m_HRvoxels[n].reject_d_std();
+ }
+
+ if(!m_HRvoxels[n].propose_S0()){ //Try S0
+ m_HRvoxels[n].compute_prior();
+ m_HRvoxels[n].compute_signal();
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy())
+ m_HRvoxels[n].accept_S0();
+ else{
+ restore_energies();
+ m_HRvoxels[n].restore_prior_totsignal();
+ m_HRvoxels[n].reject_S0();
+ }
+ }
+ else
+ m_HRvoxels[n].reject_S0();
+
+ if (m_rician){
+ if(!m_HRvoxels[n].propose_tau()){ //Try tau_HR
+ m_HRvoxels[n].compute_prior();
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy())
+ m_HRvoxels[n].accept_tau();
+ else{
+ restore_energies();
+ m_HRvoxels[n].restore_prior();
+ m_HRvoxels[n].reject_tau();
+ }
+ }
+ else
+ m_HRvoxels[n].reject_tau();
+ }
+
+ for(int f=0; f<m_numfibres; f++){ //For each fibre in the HR voxel
+ if(!m_HRvoxels[n].propose_th(f)){ //Try theta
+ m_HRvoxels[n].compute_prior();
+ m_HRvoxels[n].compute_signal();
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy())
+ m_HRvoxels[n].accept_th(f);
+ else{
+ restore_energies();
+ m_HRvoxels[n].restore_prior_totsignal();
+ m_HRvoxels[n].reject_th(f);
+ }
+ }
+ else
+ m_HRvoxels[n].reject_th(f);
+
+
+ if(!m_HRvoxels[n].propose_ph(f)){ //Try phi
+ m_HRvoxels[n].compute_prior();
+ m_HRvoxels[n].compute_signal();
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy())
+ m_HRvoxels[n].accept_ph(f);
+ else{
+ restore_energies();
+ m_HRvoxels[n].restore_prior_totsignal();
+ m_HRvoxels[n].reject_ph(f);
+ }
+ }
+ else
+ m_HRvoxels[n].reject_ph(f);
+
+
+ if(!m_HRvoxels[n].propose_f(f)){ //Try f
+ if(!m_HRvoxels[n].reject_f_sum()){
+ m_HRvoxels[n].compute_fsum_prior();
+ m_HRvoxels[n].compute_prior();
+ m_HRvoxels[n].compute_signal();
+ compute_prior();
+ compute_likelihood();
+ compute_posterior();
+ if(test_energy())
+ m_HRvoxels[n].accept_f(f);
+ else{
+ restore_energies();
+ m_HRvoxels[n].restore_prior_totsignal();
+ m_HRvoxels[n].restore_fsum_prior();
+ m_HRvoxels[n].reject_f(f);
+ }
+ }
+ else //else for rejectin fsum>1
+ m_HRvoxels[n].reject_f(f);
+ }
+ else //else for rejecting rejflag returned from propose_f()
+ m_HRvoxels[n].reject_f(f);
+ } //end anisotropic compartment parameters
+
+ } //end HR voxel parameters
+}
+
+} //end namespace
diff --git a/rubixvox.h b/rubixvox.h
new file mode 100644
index 0000000..572104a
--- /dev/null
+++ b/rubixvox.h
@@ -0,0 +1,729 @@
+/* rubixvox.h: Classes utilized in RubiX */
+
+/* Stam Sotiropoulos, FMRIB Analysis Group */
+
+/* Copyright (C) 2012 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#if !defined(rubixvox_h)
+#define rubixvox_h
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#define WANT_STREAM
+#define WANT_MATH
+#include <string>
+#include <math.h>
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/miscprob.h"
+#include "stdlib.h"
+#include "utils/log.h"
+
+using namespace std;
+using namespace NEWMAT;
+using namespace MISCMATHS;
+
+namespace RUBIX{
+ const float maxfloat=1e10;
+ const float minfloat=1e-10;
+ const float maxlogfloat=23;
+ const float minlogfloat=-23;
+
+ //////////////////////////////////////////////////////////////////////
+ // RFibre: Models one anisotropic compartment in an HR voxel
+ //////////////////////////////////////////////////////////////////////
+ class RFibre{
+ float m_th; //Current/candidate MCMC state
+ float m_ph;
+ float m_f;
+ ColumnVector m_vec; //Holds the current/candidate orientation in cartesian coordinates
+ ColumnVector m_vec_old;
+ float m_th_prop; //Standard deviation for Gaussian proposal distributions of parameters
+ float m_ph_prop;
+ float m_f_prop;
+ float m_th_old; //Last accepted value. If a sample is rejected this value is restored
+ float m_ph_old;
+ float m_f_old;
+ float m_th_ph_prior; //Priors for the model parameters
+ float m_f_prior;
+ float m_th_ph_old_prior;
+ float m_f_old_prior;
+ float m_prior_en; //Joint Prior
+ float m_old_prior_en;
+ int m_th_acc;
+ int m_th_rej;
+ int m_ph_acc;
+ int m_ph_rej;
+ int m_f_acc;
+ int m_f_rej;
+ bool f_ard; //By default ARD is on, on the volume fraction f
+ ColumnVector m_SignalHR; //Vector that stores the predicted signal from the anisotropic compartment during the candidate/current MCMC state at High-Res measurement points
+ ColumnVector m_SignalHR_old;
+ ColumnVector m_SignalLR; //Vector that stores the predicted signal from the anisotropic compartment during the candidate/current MCMC state at Low-Res measurement points
+ ColumnVector m_SignalLR_old;
+
+ const Matrix& m_Orient_hyp_prior;//Matrix Nmodes x 5 that contains the hyperparameters for the orientation prior
+ //columns 1-3 contains the (x,y,z) coordinates for the mode, 4th column contains the invkappa value, 5th the Watson normalization constant
+ const float m_ardfudge;
+ const float& m_d;
+ const float& m_d_std;
+ const Matrix& m_bvecsHR; //bvecs at High-Res (3 x HR_NumPoints)
+ const Matrix& m_bvalsHR; //bvalues at High-Res (1 x HR_NumPoints)
+ const Matrix& m_bvecsLR; //bvecs at Low-Res (3 x LR_NumPoints)
+ const Matrix& m_bvalsLR; //bvalues at Low-Res (1 x HR_NumPoints)
+ const int m_modelnum; //1 for single-shell, 2 for multi-shell model
+
+ public:
+ //constructor
+ RFibre(const float th, const float ph,const float f, const Matrix& bvecsHR, const Matrix& bvalsHR,
+ const Matrix& bvecsLR, const Matrix& bvalsLR, const Matrix& Orient_hyp_prior,
+ const float& d, const float& d_std, const bool ard=true, const int modelnum=1,const float ardfudge=1):
+ m_th(th), m_ph(ph), m_f(f), f_ard(ard), m_Orient_hyp_prior(Orient_hyp_prior), m_ardfudge(ardfudge), m_d(d),
+ m_d_std(d_std),m_bvecsHR(bvecsHR), m_bvalsHR(bvalsHR), m_bvecsLR(bvecsLR), m_bvalsLR(bvalsLR), m_modelnum(modelnum){
+ m_th_old=m_th; m_ph_old=m_ph;
+ m_vec.ReSize(3);
+ m_vec<<sin(m_th)*cos(m_ph) <<sin(m_th)*sin(m_ph) <<cos(m_th);
+ m_vec_old=m_vec;
+ m_f_old=m_f;
+ m_th_prop=0.2; m_ph_prop=0.2; m_f_prop=m_f/5.0;
+
+ m_th_ph_prior=0; m_th_ph_old_prior=0;
+ m_f_prior=0; m_f_old_prior=0;
+
+ m_SignalHR.ReSize(m_bvecsHR.Ncols()); m_SignalHR=0; m_SignalHR_old=m_SignalHR;
+ m_SignalLR.ReSize(m_bvecsLR.Ncols()); m_SignalLR=0; m_SignalLR_old=m_SignalLR;
+
+ m_th_acc=0; m_th_rej=0;
+ m_ph_acc=0; m_ph_rej=0;
+ m_f_acc=0; m_f_rej=0;
+ }
+
+ ~RFibre(){}
+
+ inline float get_th() const{ return m_th;}
+ inline float get_ph() const{ return m_ph;}
+ inline const ColumnVector& getVec() const { return m_vec; }
+ inline float get_f() const{ return m_f;}
+ //inline void set_th(const float th){ m_th=th; m_vec<<sin(m_th)*cos(m_ph) <<sin(m_th)*sin(m_ph) <<cos(m_th); }
+ //inline void set_ph(const float ph){ m_ph=ph; m_vec<<sin(m_th)*cos(m_ph) <<sin(m_th)*sin(m_ph) <<cos(m_th); }
+ //inline void set_f(const float f){ m_f=f; }
+
+ inline const ColumnVector& getSignalHR() const{ return m_SignalHR; }
+ inline const ColumnVector& getSignalLR() const{ return m_SignalLR; }
+ inline void restoreSignals() { m_SignalHR=m_SignalHR_old; m_SignalLR=m_SignalLR_old; }
+ inline float get_prior() const{ return m_prior_en;}
+
+ void initialise_energies(); //Initialize energies, signals
+ void update_proposals();
+ bool compute_th_ph_prior();
+ bool compute_f_prior();
+ void compute_prior();
+ void compute_signal(); //Compute model predicted signal, only due to the anisotropic compartment
+ void restore_th_ph_prior();
+
+ bool propose_th();
+ inline void accept_th(){ m_th_acc++; }
+ void reject_th();
+
+ bool propose_ph();
+ inline void accept_ph(){ m_ph_acc++; }
+ void reject_ph();
+
+ bool propose_f();
+ inline void accept_f(){ m_f_acc++; }
+ void reject_f();
+
+ void report() const;
+
+ RFibre& operator=(const RFibre& rhs){
+ m_th=rhs.m_th; m_ph=rhs.m_ph; m_f=rhs.m_f;
+ m_vec=rhs.m_vec; m_vec_old=rhs.m_vec_old;
+ m_th_prop=rhs.m_th_prop; m_ph_prop=rhs.m_ph_prop; m_f_prop=rhs.m_f_prop;
+ m_th_old=rhs.m_th_old; m_ph_old=rhs.m_ph_old; m_f_old=rhs.m_f_old;
+ m_th_ph_prior=rhs.m_th_ph_prior; m_f_prior=rhs.m_f_prior;
+ m_th_ph_old_prior=rhs.m_th_ph_old_prior; m_f_old_prior=rhs.m_f_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;
+ f_ard=rhs.f_ard;
+ m_SignalHR=rhs.m_SignalHR; m_SignalHR_old=rhs.m_SignalHR_old;
+ m_SignalLR=rhs.m_SignalLR; m_SignalLR_old=rhs.m_SignalLR_old;
+ return *this;
+ }
+
+ };
+
+
+
+ //////////////////////////////
+ // HRvoxel //
+ //////////////////////////////
+ class HRvoxel{
+ vector<RFibre> m_fibres;
+ float m_d;
+ float m_d_old;
+ float m_d_prop;
+ float m_d_prior;
+ float m_d_old_prior;
+ float m_d_acc;
+ float m_d_rej;
+
+ float m_d_std;
+ float m_d_std_old;
+ float m_d_std_prop;
+ float m_d_std_prior;
+ float m_d_std_old_prior;
+ float m_d_std_acc;
+ float m_d_std_rej;
+
+ float m_S0;
+ float m_S0_old;
+ float m_S0_prop;
+ float m_S0_prior;
+ float m_S0_old_prior;
+ float m_S0_acc;
+ float m_S0_rej;
+
+ float m_tau; //Noise at a High-res voxel (precision)
+ float m_tau_old;
+ float m_tau_prop;
+ float m_tau_prior;
+ float m_tau_old_prior;
+ float m_tau_acc;
+ float m_tau_rej;
+
+ const float& m_mean_d; //hyperparameters for d prior
+ const float& m_stdev_d;
+ const float& m_mean_fsum; //hyperparameters for fsum prior
+ const float& m_stdev_fsum;
+ float m_fsum_prior;
+ float m_fsum_old_prior;
+
+ float m_prior_en; //Joint Prior
+ float m_old_prior_en;
+ ColumnVector m_iso_SignalHR; //Vector that stores the predicted signal from the isotropic compartment only
+ ColumnVector m_iso_SignalHR_old;
+ ColumnVector m_iso_SignalLR; //Vector that stores the predicted signal from the isotropic compartment only
+ ColumnVector m_iso_SignalLR_old;
+
+ ColumnVector m_SignalHR; //Vector that stores the total predicted signal from the specific voxel at High-Res measurement points
+ ColumnVector m_SignalHR_old;
+ ColumnVector m_SignalLR; //Vector that stores the total predicted signal from the specific voxel at Low-Res measurement points
+ ColumnVector m_SignalLR_old;
+
+ const Matrix& m_Orient_hyp_prior;//Matrix Nmodes x 5 that contains the hyperparameters for the orientation prior
+
+ const Matrix& m_bvecsHR; //bvecs at High-Res (3 x HR_NumPoints)
+ const Matrix& m_bvalsHR; //bvalues at High-Res (1 x HR_NumPoints)
+ const Matrix& m_bvecsLR; //bvecs at Low-Res (3 x LR_NumPoints)
+ const Matrix& m_bvalsLR; //bvalues at Low-Res (1 x HR_NumPoints)
+ const int m_modelnum; //1 for single-shell, 2 for multi-shell model
+ const int m_numfibres; //number of fibres in this HR voxel
+ const float m_ardfudge;
+ const bool m_fsumPrior_ON; //Flag that indicates whether a prior on the fsum will be imposed (based on the LRvox neighbourhood mean)
+ const bool m_dPrior_ON; //Flag that indicates whether a prior on diffusivity will be imposed (based on the LRvox neighbourhood mean)
+ const bool m_rician; //indicates whether Rician Noise model is used
+
+ //const ColumnVector& LRvox_inter; //array with indices on LR voxels intersected by this HR
+
+ public:
+ //Constructor
+ HRvoxel(const Matrix& bvecsHR, const Matrix& bHR,
+ const Matrix& bvecsLR, const Matrix& bLR,
+ const Matrix& Orient_hyp_prior,
+ const float& mean_d, const float& stdev_d, const float& mean_fsum, const float& stdev_fsum,
+ const int N, const int modelnum=1,const float ardfudge=1, const bool fsumPrior_ON=false, const bool dPrior_ON=false, const bool rician=false):
+ m_mean_d(mean_d), m_stdev_d(stdev_d), m_mean_fsum(mean_fsum), m_stdev_fsum(stdev_fsum),
+ m_Orient_hyp_prior(Orient_hyp_prior),
+ m_bvecsHR(bvecsHR), m_bvalsHR(bHR), m_bvecsLR(bvecsLR), m_bvalsLR(bLR), m_modelnum(modelnum), m_numfibres(N), m_ardfudge(ardfudge), m_fsumPrior_ON(fsumPrior_ON), m_dPrior_ON(dPrior_ON), m_rician(rician) {
+
+ //Initialize vectors that keep the signal from the isotropic compartment
+ m_iso_SignalHR.ReSize(m_bvecsHR.Ncols()); m_iso_SignalHR=0; m_iso_SignalHR_old=m_iso_SignalHR;
+ m_SignalHR.ReSize(m_bvecsHR.Ncols()); m_SignalHR=0; m_SignalHR_old=m_SignalHR;
+ m_iso_SignalLR.ReSize(m_bvecsLR.Ncols()); m_iso_SignalLR=0; m_iso_SignalLR_old=m_iso_SignalLR;
+ m_SignalLR.ReSize(m_bvecsLR.Ncols()); m_SignalLR=0; m_SignalLR_old=m_SignalLR;
+
+ m_d_acc=0; m_d_rej=0; m_d_prior=0; m_d_old_prior=0;
+ m_d_std_acc=0; m_d_std_rej=0; m_d_std_prior=0; m_d_std_old_prior=0;
+ m_S0_acc=0; m_S0_rej=0; m_S0_prior=0; m_S0_old_prior=0;
+ m_d=0; m_d_old=0; m_d_std=0; m_d_std_old=0; m_S0=0; m_S0_old=0;
+ m_tau_acc=0; m_tau_rej=0; m_tau_prior=0; m_tau_old_prior=0; m_tau=0; m_tau_old=0;
+ m_prior_en=0; m_old_prior_en=0;
+ m_fsum_prior=0; m_fsum_old_prior=0;
+ }
+
+ //Destructor
+ ~HRvoxel(){}
+
+
+ const vector<RFibre>& fibres() const{ return m_fibres; }
+ void addfibre(const float th, const float ph, const float f,const bool use_ard=true) {
+ RFibre fib(th,ph,f,m_bvecsHR,m_bvalsHR,m_bvecsLR,m_bvalsLR,m_Orient_hyp_prior,m_d,m_d_std,use_ard,m_modelnum,m_ardfudge);
+ m_fibres.push_back(fib);
+ }
+
+ void initialise_energies_props(); //Initialize energies, signals and standard deviations for the proposal distributions
+ inline float get_d() const{ return m_d;}
+ inline void set_d(const float d){ m_d=d; }
+ inline float get_d_std() const{ return m_d_std; }
+ inline void set_d_std(const float d_std){ m_d_std=d_std; }
+ inline float get_S0() const{ return m_S0; }
+ inline void set_S0(const float S0){ m_S0=S0; }
+ inline float get_tau() const{ return m_tau; }
+ inline void set_tau(const float tau){ m_tau=tau; }
+ inline float get_prior() const { return m_prior_en; }
+
+ inline const ColumnVector& getSignalHR() const{ return m_SignalHR; }
+ inline const ColumnVector& getSignalLR() const{ return m_SignalLR; }
+
+ bool compute_d_prior();
+ bool compute_d_std_prior();
+ bool compute_S0_prior();
+ bool compute_tau_prior();
+ bool reject_f_sum(); //Check if sum of volume fractions is >1
+
+ void update_d_prior();
+ void restore_d_prior();
+
+ void compute_fsum_prior();
+ void update_fsum_prior();
+ void restore_fsum_prior();
+
+ void compute_prior(); //Compute Joint Prior energy
+ void compute_iso_signal(); //Compute the predicted signal from the isotropic compartment only
+ void compute_signal(); //Compute the total predicted signal
+
+ bool propose_d();
+ void accept_d(){ m_d_acc++; }
+ void reject_d();
+
+ bool propose_d_std();
+ void accept_d_std(){ m_d_std_acc++; }
+ void reject_d_std();
+
+ bool propose_S0();
+ void accept_S0(){ m_S0_acc++; }
+ void reject_S0();
+
+ bool propose_tau();
+ void accept_tau(){ m_tau_acc++; }
+ void reject_tau();
+
+ bool propose_th(int n){ return m_fibres[n].propose_th(); } //n is the fibre index from 0 to N-1
+ void accept_th(int n) { m_fibres[n].accept_th(); }
+ void reject_th(int n) { m_fibres[n].reject_th(); }
+ bool propose_ph(int n){ return m_fibres[n].propose_ph(); }
+ void accept_ph(int n) { m_fibres[n].accept_ph(); }
+ void reject_ph(int n) { m_fibres[n].reject_ph(); }
+ bool propose_f(int n) { return m_fibres[n].propose_f(); }
+ void accept_f(int n) { m_fibres[n].accept_f(); }
+ void reject_f(int n) { m_fibres[n].reject_f(); }
+
+ void restore_prior_totsignal();
+ void restore_prior();
+
+ void update_th_ph_prior(); //Used to update the conditional orientation prior when the prior parameters are jumped
+ void restore_th_ph_prior();
+
+ void update_proposals(); //Adapt standard deviation of proposal distributions during MCMC execution
+ void report() const;
+
+ HRvoxel& operator=(const HRvoxel& rhs){
+ m_fibres=rhs.m_fibres;
+ m_d=rhs.m_d; m_d_old=rhs.m_d_old; m_d_prop=rhs.m_d_prop;
+ m_d_prior=rhs.m_d_prior; m_d_old_prior=rhs.m_d_old_prior;
+ m_d_acc=rhs.m_d_acc; m_d_rej=rhs.m_d_rej;
+ m_d_std=rhs.m_d_std; m_d_std_old=rhs.m_d_std_old; m_d_std_prop=rhs.m_d_std_prop;
+ m_d_std_prior=rhs.m_d_std_prior; m_d_std_old_prior=rhs.m_d_std_old_prior;
+ m_d_std_acc=rhs.m_d_std_acc; m_d_std_rej=rhs.m_d_std_rej;
+
+ m_fsum_prior=rhs.m_fsum_prior; m_fsum_old_prior=rhs.m_fsum_old_prior;
+ m_S0=rhs.m_S0; m_S0_old=rhs.m_S0_old; m_S0_prop=rhs.m_S0_prop;
+ m_S0_prior=rhs.m_S0_prior; m_S0_old_prior=rhs.m_S0_old_prior;
+ m_S0_acc=rhs.m_S0_acc; m_S0_rej=rhs.m_S0_rej;
+ m_tau=rhs.m_tau; m_tau_old=rhs.m_tau_old; m_tau_prop=rhs.m_tau_prop;
+ m_tau_prior=rhs.m_tau_prior; m_tau_old_prior=rhs.m_tau_old_prior;
+ m_tau_acc=rhs.m_tau_acc; m_tau_rej=rhs.m_tau_rej;
+ m_prior_en=rhs.m_prior_en; m_old_prior_en=rhs.m_old_prior_en;
+ m_iso_SignalHR=rhs.m_iso_SignalHR; m_iso_SignalHR_old=rhs.m_iso_SignalHR_old;
+ m_iso_SignalLR=rhs.m_iso_SignalLR; m_iso_SignalLR_old=rhs.m_iso_SignalLR_old;
+ m_SignalHR=rhs.m_SignalHR; m_SignalHR_old=rhs.m_SignalHR_old;
+ m_SignalLR=rhs.m_SignalLR; m_SignalLR_old=rhs.m_SignalLR_old;
+ return *this;
+ }
+
+ };
+
+
+
+
+ /////////////////////////////////////////////////////////////////////////////
+ // Orient_Prior_Mode: Models a single mode of the orientation prior ///
+ // The orientation prior is a sum of Watson distributions centred ///
+ // around different modes ///
+ /////////////////////////////////////////////////////////////////////////////
+ class Orient_Prior_Mode{
+ float m_th;
+ float m_th_old;
+ float m_th_prop;
+ float m_th_prior;
+ float m_th_old_prior;
+ float m_th_acc;
+ float m_th_rej;
+
+ float m_ph;
+ float m_ph_old;
+ float m_ph_prop;
+ float m_ph_prior;
+ float m_ph_old_prior;
+ float m_ph_acc;
+ float m_ph_rej;
+
+ float m_invkappa; //Dispersion Index of a Watson distribution (1/kappa actually)
+ float m_invkappa_old;
+ float m_invkappa_prop;
+ float m_invkappa_prior;
+ float m_invkappa_old_prior;
+ float m_invkappa_acc;
+ float m_invkappa_rej;
+
+ float m_Watson_norm; //this is a function of current m_kappa
+ float m_Watson_norm_old;
+ ColumnVector m_vec;
+ ColumnVector m_vec_old;
+ float m_prior_en; //Joint HyperPrior
+ float m_old_prior_en;
+
+ const bool m_kappa_ard; //Flag for setting ARD on the dispersion index
+
+ public:
+ //Constructor
+ Orient_Prior_Mode(bool kappa_ard): m_kappa_ard(kappa_ard) {
+ m_th=M_PI/2; m_th_old=m_th;
+ m_ph=0; m_ph_old=m_ph;
+ m_vec.ReSize(3);
+ m_vec<<sin(m_th)*cos(m_ph) <<sin(m_th)*sin(m_ph) <<cos(m_th);
+ m_vec_old=m_vec;
+ m_invkappa=0.02; m_invkappa_old=m_invkappa;
+ m_Watson_norm=compute_Watson_norm(); m_Watson_norm_old=m_Watson_norm;
+ m_th_prop=0.2; m_ph_prop=0.2; m_invkappa_prop=0.2;
+
+ m_th_prior=0; m_th_old_prior=0;
+ m_ph_prior=0; m_ph_old_prior=0;
+ m_invkappa_prior=0; m_invkappa_old_prior=0;
+ m_prior_en=0; m_old_prior_en=0;
+
+ m_th_acc=0; m_th_rej=0;
+ m_ph_acc=0; m_ph_rej=0;
+ m_invkappa_acc=0; m_invkappa_rej=0;
+ }
+
+ //Destructor
+ ~Orient_Prior_Mode(){}
+
+ inline float get_th() const{ return m_th;}
+ inline float get_ph() const{ return m_ph;}
+ inline float get_invkappa() const{ return m_invkappa;}
+ inline void set_th(const float th) { m_th=th; m_vec<< sin(m_th)*cos(m_ph) << sin(m_th)*sin(m_ph) << cos(m_th);}
+ inline void set_ph(const float ph) { m_ph=ph; m_vec<< sin(m_th)*cos(m_ph) << sin(m_th)*sin(m_ph) << cos(m_th);}
+ inline void set_invkappa(const float invkappa) { m_invkappa=invkappa; filter_invkappa(); m_Watson_norm=compute_Watson_norm(); }
+ void filter_invkappa(); //Filter out extreme values of invkappa to ensure calculations are numerically stable
+
+ inline float get_prior() const{ return m_prior_en;}
+ inline const ColumnVector& getVec() const { return m_vec; }
+ inline float get_Watson_norm() const{ return m_Watson_norm; }
+ float compute_Watson_norm();
+
+ void initialise_energies();
+ bool compute_th_prior();
+ bool compute_ph_prior();
+ bool compute_invkappa_prior();
+ void compute_prior(); //Compute Joint Prior energy
+ void update_proposals();
+
+ bool propose_th();
+ inline void accept_th(){ m_th_acc++; }
+ void reject_th();
+
+ bool propose_ph();
+ inline void accept_ph(){ m_ph_acc++; }
+ void reject_ph();
+
+ bool propose_invkappa();
+ inline void accept_invkappa(){ m_invkappa_acc++; }
+ void reject_invkappa();
+
+ Orient_Prior_Mode& operator=(const Orient_Prior_Mode& rhs){
+ m_th=rhs.m_th;
+ m_th_old=rhs.m_th_old;
+ m_th_prop=rhs.m_th_prop;
+ m_th_prior=rhs.m_th_prior; m_th_old_prior=rhs.m_th_old_prior;
+ m_th_acc=rhs.m_th_acc; m_th_rej=rhs.m_th_rej;
+ m_ph=rhs.m_ph; m_ph_old=rhs.m_ph_old; m_ph_prop=rhs.m_ph_prop;
+ m_ph_prior=rhs.m_ph_prior; m_ph_old_prior=rhs.m_ph_old_prior;
+ m_ph_acc=rhs.m_ph_acc; m_ph_rej=rhs.m_ph_rej;
+ m_invkappa=rhs.m_invkappa; m_invkappa_old=rhs.m_invkappa_old; m_invkappa_prop=rhs.m_invkappa_prop;
+ m_invkappa_prior=rhs.m_invkappa_prior; m_invkappa_old_prior=rhs.m_invkappa_old_prior;
+ m_invkappa_acc=rhs.m_invkappa_acc; m_invkappa_rej=rhs.m_invkappa_rej;
+ m_Watson_norm=rhs.m_Watson_norm; m_Watson_norm_old=rhs.m_Watson_norm_old;
+ m_vec=rhs.m_vec; m_vec_old=rhs.m_vec_old;
+ m_prior_en=rhs.m_prior_en; m_old_prior_en=rhs.m_old_prior_en;
+ return *this;
+ }
+ };
+
+
+
+ //////////////////////////////
+ // LRvoxel //
+ //////////////////////////////
+ class LRvoxel{
+ vector<HRvoxel> m_HRvoxels;
+ vector<Orient_Prior_Mode> m_PModes;
+
+ float m_tauLR; //models noise at Low-res
+ float m_tauLR_old;
+ float m_tauLR_prop;
+ float m_tauLR_prior;
+ float m_tauLR_old_prior;
+ float m_tauLR_acc;
+ float m_tauLR_rej;
+
+ float m_S0LR; //models S0 intensity at the Low-Res acquisition
+ float m_S0LR_old;
+ float m_S0LR_prop;
+ float m_S0LR_prior;
+ float m_S0LR_old_prior;
+ float m_S0LR_acc;
+ float m_S0LR_rej;
+ Matrix m_Orient_hyp_prior; //Matrix Nmodes x 5 that contains the hyperparameters for the orientation prior
+ //columns 1-3 contains the (x,y,z) coordinates for the mode, 4th column contains the invkappa value, 5th the Watson normalization constant
+
+ float m_mean_d; //models mean of d hyperprior for the High-Res voxels intersected by the LR_voxel
+ float m_mean_d_old;
+ float m_mean_d_prop;
+ float m_mean_d_prior;
+ float m_mean_d_old_prior;
+ float m_mean_d_acc;
+ float m_mean_d_rej;
+
+ float m_stdev_d; //models std_dev of d hyperprior for the High-Res voxels intersected by the LR_voxel
+ float m_stdev_d_old;
+ float m_stdev_d_prop;
+ float m_stdev_d_prior;
+ float m_stdev_d_old_prior;
+ float m_stdev_d_acc;
+ float m_stdev_d_rej;
+
+ float m_mean_fsum; //models mean of fsum hyperprior for the High-Res voxels intersected by the LR_voxel
+ float m_mean_fsum_old;
+ float m_mean_fsum_prop;
+ float m_mean_fsum_prior;
+ float m_mean_fsum_old_prior;
+ float m_mean_fsum_acc;
+ float m_mean_fsum_rej;
+
+ float m_stdev_fsum; //models std_dev of fsum hyperprior for the High-Res voxels intersected by the LR_voxel
+ float m_stdev_fsum_old;
+ float m_stdev_fsum_prop;
+ float m_stdev_fsum_prior;
+ float m_stdev_fsum_old_prior;
+ float m_stdev_fsum_acc;
+ float m_stdev_fsum_rej;
+
+ float m_prior_en; //Joint Prior Energy
+ float m_old_prior_en;
+ float m_likelihood_en; //Likelihood Energy
+ float m_old_likelihood_en;
+ float m_posterior_en; //Posterior Energy
+ float m_old_posterior_en;
+
+ ColumnVector m_logdataLR; //Log of Low-Res Data for the specific LR voxel (use it in Rician Energy)
+ vector<ColumnVector> m_logdataHR; //Log of High-Res Data for all contained HR voxels (use it in Rician Energy)
+
+ const ColumnVector& m_dataLR; //Low-Res Data for the specific LR voxel
+ const vector<ColumnVector>& m_dataHR; //High-Res Data for all contained HR voxels
+ const Matrix& m_bvecsHR; //bvecs at High-Res (3 x HR_NumPoints)
+ const Matrix& m_bvalsHR; //bvalues at High-Res (1 x HR_NumPoints)
+ const Matrix& m_bvecsLR; //bvecs at Low-Res (3 x LR_NumPoints)
+ const Matrix& m_bvalsLR; //bvalues at Low-Res (1 x HR_NumPoints)
+ const int m_modelnum; //1 for single-shell, 2 for multi-shell model
+ const int m_numfibres; //Number of fibres in each HR voxel
+ const int m_Nmodes; //Number of modes for the Orientation Prior
+ const float m_ardfudge;
+ const bool m_allard; //Flag for setting ARD on for all fibres in all HR voxels
+ const bool m_Noard; //Flag for setting ARD off for all fibres in all HR voxels
+ const bool m_kappa_ard; //Flag for setting ARD on the dispersion index
+ const bool m_fsumPrior_ON; //Flag for setting on a prior on fsums across intersected HR voxels
+ const bool m_dPrior_ON; //Flag for setting on a prior on the diffusivity across intersected HR voxels
+ const bool m_rician; //Flag for using a Rician noise model
+ const ColumnVector& m_HRweights;//Holds the volume fraction each HR voxel occupies out of the LR one
+
+ public:
+ //Constructor
+ LRvoxel(const Matrix& bvecsHR, const Matrix& bHR,
+ const Matrix& bvecsLR, const Matrix& bLR,
+ const ColumnVector& dataLR, const vector<ColumnVector>& dataHR, const int N, const int Nmodes, const ColumnVector& HRweights, const int modelnum=1, const float ardfudge=1, const bool allard=false, const bool Noard=false, const bool kappa_ard=false, const bool fsumPrior_ON=false, const bool dPrior_ON=false, const bool rician=false):
+ m_dataLR(dataLR), m_dataHR(dataHR),m_bvecsHR(bvecsHR), m_bvalsHR(bHR), m_bvecsLR(bvecsLR), m_bvalsLR(bLR),
+ m_modelnum(modelnum), m_numfibres(N), m_Nmodes(Nmodes), m_ardfudge(ardfudge), m_allard(allard), m_Noard(Noard), m_kappa_ard(kappa_ard), m_fsumPrior_ON(fsumPrior_ON), m_dPrior_ON(dPrior_ON), m_rician(rician), m_HRweights(HRweights) {
+
+ m_S0LR=0; m_S0LR_old=0; m_S0LR_prior=0; m_S0LR_old_prior=0; m_S0LR_acc=0; m_S0LR_rej=0; m_S0LR_prop=0.2;
+ m_tauLR=0; m_tauLR_old=0; m_tauLR_prior=0; m_tauLR_old_prior=0; m_tauLR_acc=0; m_tauLR_rej=0; m_tauLR_prop=0.2;
+
+ m_mean_d=0; m_mean_d_old=0; m_mean_d_prior=0; m_mean_d_old_prior=0; m_mean_d_acc=0; m_mean_d_rej=0; m_mean_d_prop=0.001;
+ m_stdev_d=0; m_stdev_d_old=0; m_stdev_d_prior=0; m_stdev_d_old_prior=0; m_stdev_d_acc=0; m_stdev_d_rej=0; m_stdev_d_prop=0.001;
+
+ m_mean_fsum=0; m_mean_fsum_old=0; m_mean_fsum_prior=0; m_mean_fsum_old_prior=0; m_mean_fsum_acc=0; m_mean_fsum_rej=0; m_mean_fsum_prop=0.1;
+ m_stdev_fsum=0; m_stdev_fsum_old=0; m_stdev_fsum_prior=0; m_stdev_fsum_old_prior=0; m_stdev_fsum_acc=0; m_stdev_fsum_rej=0; m_stdev_fsum_prop=0.1;
+
+ m_prior_en=0; m_old_prior_en=0;
+ m_likelihood_en=0; m_old_likelihood_en=0;
+ m_posterior_en=0; m_old_posterior_en=0;
+
+ for (int m=1; m<=m_Nmodes; m++){ //Add Modes for the orientation Prior
+ Orient_Prior_Mode pMod(m_kappa_ard);
+ m_PModes.push_back(pMod);
+ }
+ m_Orient_hyp_prior.ReSize(m_Nmodes,5);
+
+ for (unsigned int m=1; m<=m_dataHR.size(); m++){ //Add HRvoxel Objects
+ HRvoxel HRv(m_bvecsHR, m_bvalsHR, m_bvecsLR, m_bvalsLR, m_Orient_hyp_prior, m_mean_d, m_stdev_d, m_mean_fsum, m_stdev_fsum, m_numfibres, m_modelnum, m_ardfudge, m_fsumPrior_ON, m_dPrior_ON, m_rician);
+
+ m_HRvoxels.push_back(HRv);
+ }
+
+ // ofstream myfile; //Debugging code
+ //myfile.open("/Users/stam/Rubix_data/Energies.txt", ios::trunc);
+ //myfile.close();
+
+ if (m_rician){//Store the log of the data for energy calculations
+ m_logdataLR=dataLR; m_logdataHR=dataHR;
+ for (int m=1; m<=dataLR.Nrows(); m++){
+ if (dataLR(m)!=0)
+ m_logdataLR(m)=log(dataLR(m));
+ else
+ m_logdataLR(m)=minlogfloat;
+ }
+ for (unsigned int n=0; n<dataHR.size(); n++)
+ for (int m=1; m<=dataHR[n].Nrows(); m++){
+ if (dataHR[n](m)!=0)
+ m_logdataHR[n](m)=log(dataHR[n](m));
+ else
+ m_logdataHR[n](m)=minlogfloat;
+ }
+ }
+ }
+
+ //Destructor
+ ~LRvoxel(){ }
+
+ const vector<HRvoxel>& HRvoxels() const {return m_HRvoxels;}
+ const vector<Orient_Prior_Mode>& PModes() const {return m_PModes;}
+
+ void update_Orient_hyp_prior(); //Update the matrix that keeps information on the orientation prior parameters
+ void update_Orient_hyp_prior(int M); //Update the entry only for a specific Mode 0<=M<N_modes
+ void set_HRparams(const int n, const float d, const float S0, const ColumnVector& th, const ColumnVector& ph, const ColumnVector& f); //Set params for a single HR voxel
+ void set_HRparams(const int n, const float d, const float d_std,const float S0, const ColumnVector& th, const ColumnVector& ph, const ColumnVector& f); //Set params for a single HR voxel
+ void set_Priorparams(const ColumnVector& th, const ColumnVector& ph, const ColumnVector& invkappa); //Set params for all modes of orientation priors
+ float get_S0LR() const { return m_S0LR; }
+ void set_S0LR(const float S0) { m_S0LR=S0; }
+ float get_tauLR() const { return m_tauLR; }
+ void set_tauLR(const float tau) { m_tauLR=tau; }
+ void set_tauHR(const int n, const float tau) { m_HRvoxels[n].set_tau(tau); }
+ float get_likelihood_energy() const { return m_likelihood_en; }
+ float get_prior_energy() const { return m_prior_en; }
+ void initialise_energies();
+
+ bool propose_S0LR();
+ void accept_S0LR(){ m_S0LR_acc++; }
+ void reject_S0LR();
+ bool compute_S0LR_prior();
+
+ bool propose_tauLR();
+ void accept_tauLR(){ m_tauLR_acc++; }
+ void reject_tauLR();
+ bool compute_tauLR_prior();
+
+ void compute_prior();
+ void compute_likelihood();
+ void compute_posterior();
+ bool test_energy() const;
+ void restore_energies();
+ void restore_Prior_Posterior();
+ void update_proposals();
+ void jump(); //Single MCMC iteration with all parameters jumped
+
+
+ bool propose_meand();
+ void accept_meand(){ m_mean_d_acc++; }
+ void reject_meand();
+ bool compute_meand_prior();
+ bool propose_stdevd();
+ void accept_stdevd(){ m_stdev_d_acc++; }
+ void reject_stdevd();
+ bool compute_stdevd_prior();
+ void set_meand(const float meand) { m_mean_d=meand; }
+ void set_stdevd(const float stdevd) { m_stdev_d=stdevd; }
+ float get_meand() const{ return m_mean_d;}
+ float get_stdevd() const{ return m_stdev_d;}
+
+ bool propose_mean_fsum();
+ void accept_mean_fsum(){ m_mean_fsum_acc++; }
+ void reject_mean_fsum();
+ bool compute_mean_fsum_prior();
+ bool propose_stdev_fsum();
+ void accept_stdev_fsum(){ m_stdev_fsum_acc++; }
+ void reject_stdev_fsum();
+ bool compute_stdev_fsum_prior();
+ void set_mean_fsum(const float fsum) { m_mean_fsum=fsum; }
+ void set_stdev_fsum(const float stdev_fsum) { m_stdev_fsum=stdev_fsum; }
+ float get_mean_fsum() const{ return m_mean_fsum;}
+ float get_stdev_fsum() const{ return m_stdev_fsum;}
+
+
+ LRvoxel& operator=(const LRvoxel& rhs){
+ m_HRvoxels=rhs.m_HRvoxels; m_PModes=rhs.m_PModes;
+ m_tauLR=rhs.m_tauLR; m_tauLR_old=rhs.m_tauLR_old;
+ m_tauLR_prop=rhs.m_tauLR_prop;
+ m_tauLR_prior=rhs.m_tauLR_prior; m_tauLR_old_prior=rhs.m_tauLR_old_prior;
+ m_tauLR_acc=rhs.m_tauLR_acc; m_tauLR_rej=rhs.m_tauLR_rej;
+ m_S0LR=rhs.m_S0LR; m_S0LR_old=rhs.m_S0LR_old;
+ m_S0LR_prop=rhs.m_S0LR_prop;
+ m_S0LR_prior=rhs.m_S0LR_prior; m_S0LR_old_prior=rhs.m_S0LR_old_prior;
+ m_S0LR_acc=rhs.m_S0LR_acc; m_S0LR_rej=rhs.m_S0LR_rej;
+ m_Orient_hyp_prior=rhs.m_Orient_hyp_prior;
+ m_prior_en=rhs.m_prior_en; m_old_prior_en=rhs.m_old_prior_en;
+ m_likelihood_en=rhs.m_likelihood_en; m_old_likelihood_en=rhs.m_old_likelihood_en;
+ m_posterior_en=rhs.m_posterior_en; m_old_posterior_en=rhs.m_old_posterior_en;
+ m_logdataLR=rhs.m_logdataLR; m_logdataHR=rhs.m_logdataHR;
+
+ m_mean_d=rhs.m_mean_d; m_mean_d_old=rhs.m_mean_d_old;
+ m_mean_d_prop=rhs.m_mean_d_prop;
+ m_mean_d_prior=rhs.m_mean_d_prior; m_mean_d_old_prior=rhs.m_mean_d_old_prior;
+ m_mean_d_acc=rhs.m_mean_d_acc; m_mean_d_rej=rhs.m_mean_d_rej;
+ m_stdev_d=rhs.m_stdev_d; m_stdev_d_old=rhs.m_stdev_d_old;
+ m_stdev_d_prop=rhs.m_stdev_d_prop;
+ m_stdev_d_prior=rhs.m_stdev_d_prior; m_stdev_d_old_prior=rhs.m_stdev_d_old_prior;
+ m_stdev_d_acc=rhs.m_stdev_d_acc; m_stdev_d_rej=rhs.m_stdev_d_rej;
+
+ m_mean_fsum=rhs.m_mean_fsum; m_mean_fsum_old=rhs.m_mean_fsum_old;
+ m_mean_fsum_prop=rhs.m_mean_fsum_prop;
+ m_mean_fsum_prior=rhs.m_mean_fsum_prior; m_mean_fsum_old_prior=rhs.m_mean_fsum_old_prior;
+ m_mean_fsum_acc=rhs.m_mean_fsum_acc; m_mean_fsum_rej=rhs.m_mean_fsum_rej;
+ m_stdev_fsum=rhs.m_stdev_fsum; m_stdev_fsum_old=rhs.m_stdev_fsum_old;
+ m_stdev_fsum_prop=rhs.m_stdev_fsum_prop;
+ m_stdev_fsum_prior=rhs.m_stdev_fsum_prior; m_stdev_fsum_old_prior=rhs.m_stdev_fsum_old_prior;
+ m_stdev_fsum_acc=rhs.m_stdev_fsum_acc; m_stdev_fsum_rej=rhs.m_stdev_fsum_rej;
+
+ return *this;
+ }
+ };
+
+
+
+}
+
+#endif
--
GitLab