*** empty log message ***

rubix.cc

+/* 
+   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;
+}

rubix.h

+/*  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

rubixoptions.cc

+/*  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;
+    }      
+  }
+
+}

rubixoptions.h

+/*  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
+
+
+
+
+