Skip to content
Snippets Groups Projects
Commit 4b2fa2e7 authored by Matthew Webster's avatar Matthew Webster
Browse files

This commit was manufactured by cvs2svn to create tag 'fsl-5_0_8'.

Sprout from master 2014-02-06 17:27:54 UTC Matthew Webster <mwebster@fmrib.ox.ac.uk> 'Changes for osx 109'
Cherrypick from master 2014-07-08 16:29:20 UTC Matthew Webster <mwebster@fmrib.ox.ac.uk> 'check for bvacs/bvals .txt files':
    CUDA/bedpostx_gpu
    bedpostx
    eddy_correct
Cherrypick from master 2012-09-04 15:58:14 UTC Stamatios Sotiropoulos <stam@fmrib.ox.ac.uk> 'Hide grad_nonlin options':
    dtifitOptions.h
Cherrypick from master 2014-12-01 11:59:10 UTC Matthew Webster <mwebster@fmrib.ox.ac.uk> 'No need to have a temp file at all':
    maskdyads
Delete:
    CUDA/bedpostx_multigpu_LSF
    CUDA/bedpostx_multigpu_PBS.sh
    RLR_MAP_routines.cpp
    RLR_MAP_routines.h
    doc/fdt_bedpostx.html
    doc/fdt_bedpostx_parallel.html
    doc/fdt_biggest.html
    doc/fdt_display.html
    doc/fdt_dtifit.html
    doc/fdt_eddy.html
    doc/fdt_images/eddy.gif
    doc/fdt_images/eddy.tif
    doc/fdt_images/eddy2.gif
    doc/fdt_images/fdt_bedpost.gif
    doc/fdt_images/fdt_bedpost.tiff
    doc/fdt_images/fdt_bedpostx.gif
    doc/fdt_images/fdt_fa.gif
    doc/fdt_images/fdt_fa_big.png
    doc/fdt_images/fdt_gui.gif
    doc/fdt_images/fdt_gui.tiff
    doc/fdt_images/fdt_l1.gif
    doc/fdt_images/fdt_lines.gif
    doc/fdt_images/fdt_lines_subs.gif
    doc/fdt_images/fdt_probtrackx.gif
    doc/fdt_images/fdt_rgb.gif
    doc/fdt_images/fdt_rgb_subs.gif
    doc/fdt_images/fdt_seeds2targets_quant_eg.gif
    doc/fdt_images/fdt_seeds2targets_thal.gif
    doc/fdt_images/fdt_simple_tract3.gif
    doc/fdt_images/fdt_spherical_polars.gif
    doc/fdt_images/fdt_twomasks_tracts.gif
    doc/fdt_images/fdt_vecreg.gif
    doc/fdt_images/fdt_vectors.jpg
    doc/fdt_images/fdt_vectorsx.jpg
    doc/fdt_images/fdt_xfibres_axial.png
    doc/fdt_images/fdt_xfibres_coronal.png
    doc/fdt_images/fsl-bg.jpg
    doc/fdt_images/fsl-logo.jpg
    doc/fdt_images/fslview_dti.gif
    doc/fdt_images/fslview_x.gif
    doc/fdt_images/gui_intro.gif
    doc/fdt_images/gui_intro.tif
    doc/fdt_images/qboot.gif
    doc/fdt_pipeline.html
    doc/fdt_probtrackx.html
    doc/fdt_reg.html
    doc/fdt_surface.html
    doc/fdt_thresh.html
    doc/fdt_top.html
    doc/fdt_utils.html
    doc/fdt_vecreg.html
    doc/index.html
    nldtifit.cpp
    nldtifit.h
    rubix.cc
    rubix.h
    rubix_parallel.sh
    rubix_postproc.sh
    rubix_preproc.sh
    rubixoptions.cc
    rubixoptions.h
    rubixvox.cc
    rubixvox.h
parent d3b3343a
No related branches found
No related merge requests found
/* 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_sumfsamples;
Matrix m_meandsamples;
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;
RowVector m_mean_sumfsamples;
RowVector m_mean_meandsamples;
int m_nsamps;
const int m_njumps;
const int m_sample_every;
const int m_Nmodes;
const bool m_rician;
const bool m_fsumPrior;
const bool m_dPrior;
//const string m_logdir;
public:
LRSamples(int nvoxels, const int njumps, const int sample_every, const int Nmodes, const bool rician=false, const bool fsumPrior=false, const bool dPrior=false):
m_njumps(njumps),m_sample_every(sample_every), m_Nmodes(Nmodes), m_rician(rician), m_fsumPrior(fsumPrior), m_dPrior(dPrior){
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;
}
if (m_fsumPrior){
m_sumfsamples.ReSize(nsamples,nvoxels); m_sumfsamples=0;
m_mean_sumfsamples.ReSize(nvoxels); m_mean_sumfsamples=0;
}
if (m_dPrior){
m_meandsamples.ReSize(nsamples,nvoxels); m_meandsamples=0;
m_mean_meandsamples.ReSize(nvoxels); m_mean_meandsamples=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 HRvoxels within a LRvoxel
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 vector<Matrix>& m_bvecsHR; //bvecs at High-Res (HRvoxels within a LRvoxel x 3 x HR_NumPoints)
const vector<Matrix>& m_bvalsHR; //bvalues at High-Res (HRvoxels within a LRvoxel x 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 vector<Matrix>& bvecsHR, const vector<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(), opts.noS0jump.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
This diff is collapsed.
#!/bin/sh
subjdir=$1
filterflag=$2
out_dir=${subjdir}.rubiX
maskLR=nodif_brain_maskLR
numfib=`${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_0000/f*samples* | wc -w | awk '{print $1}'`
fib=1
while [ $fib -le $numfib ]
do
${FSLDIR}/bin/fslmerge -z ${out_dir}/merged_th${fib}samples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/th${fib}samples*`
${FSLDIR}/bin/fslmerge -z ${out_dir}/merged_ph${fib}samples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/ph${fib}samples*`
${FSLDIR}/bin/fslmerge -z ${out_dir}/merged_f${fib}samples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/f${fib}samples*`
${FSLDIR}/bin/fslmaths ${out_dir}/merged_f${fib}samples -Tmean ${out_dir}/mean_f${fib}samples
${FSLDIR}/bin/make_dyadic_vectors ${out_dir}/merged_th${fib}samples ${out_dir}/merged_ph${fib}samples ${out_dir}/mean_f${fib}samples ${out_dir}/dyads${fib} 95
if [ $fib -ge 2 ];then
${FSLDIR}/bin/maskdyads ${out_dir}/dyads${fib} ${out_dir}/mean_f${fib}samples
fi
fib=$(($fib + 1))
done
if [ `${FSLDIR}/bin/imtest ${out_dir}/mean_f1samples` -eq 1 ];then
${FSLDIR}/bin/fslmaths ${out_dir}/mean_f1samples -mul 0 ${out_dir}/mean_fsumsamples
fib=1
while [ $fib -le $numfib ]
do
${FSLDIR}/bin/fslmaths ${out_dir}/mean_fsumsamples -add ${out_dir}/mean_f${fib}samples ${out_dir}/mean_fsumsamples
fib=$(($fib + 1))
done
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/mean_dsamples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/mean_dsamples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/mean_dsamples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/mean_S0samples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/mean_S0samples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/mean_S0samples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/mean_tausamples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/mean_tausamples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/mean_tausamples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/mean_tau_LRsamples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/mean_tau_LRsamples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/mean_tau_LRsamples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/mean_S0_LRsamples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/mean_S0_LRsamples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/mean_S0_LRsamples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/En_Lik_LRsamples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/En_Lik_LRsamples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/En_Lik_LRsamples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/En_Prior_LRsamples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/En_Prior_LRsamples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/En_Prior_LRsamples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/mean_dstd_samples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/mean_dstd_samples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/mean_dstd_samples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/HRbrain_mask` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/HRbrain_mask `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/HRbrain_mask*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/mean_fsumPriorMode_LRsamples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/mean_fsumPriorMode_LRsamples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/mean_fsumPriorMode_LRsamples*`
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/diff_slices/data_slice_0000/mean_dPriorMode_LRsamples` -eq 1 ];then
${FSLDIR}/bin/fslmerge -z ${out_dir}/mean_dPriorMode_LRsamples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/mean_dPriorMode_LRsamples*`
fi
numModes=`${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_0000/Pk*samples* | wc -w | awk '{print $1}'`
mode=1
while [ $mode -le $numModes ]
do
${FSLDIR}/bin/fslmerge -z ${out_dir}/merged_Pth${mode}samples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/Pth${mode}samples*`
${FSLDIR}/bin/fslmerge -z ${out_dir}/merged_Pph${mode}samples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/Pph${mode}samples*`
${FSLDIR}/bin/fslmerge -z ${out_dir}/merged_Pk${mode}samples `${FSLDIR}/bin/imglob ${out_dir}/diff_slices/data_slice_*/Pk${mode}samples*`
${FSLDIR}/bin/fslmaths ${out_dir}/merged_Pk${mode}samples -Tmean ${out_dir}/mean_Pk${mode}samples
${FSLDIR}/bin/make_dyadic_vectors ${out_dir}/merged_Pth${mode}samples ${out_dir}/merged_Pph${mode}samples ${out_dir}/mean_Pk${mode}samples ${out_dir}/Pdyads${mode}
mode=$(($mode + 1))
done
echo Removing intermediate files
if [ `${FSLDIR}/bin/imtest ${out_dir}/merged_th1samples` -eq 1 ];then
if [ `${FSLDIR}/bin/imtest ${out_dir}/merged_ph1samples` -eq 1 ];then
if [ `${FSLDIR}/bin/imtest ${out_dir}/merged_f1samples` -eq 1 ];then
rm -rf ${out_dir}/diff_slices
rm -f ${out_dir}/dataLR_slice_*
rm -f ${out_dir}/dataHR_newslice_*
rm -f ${out_dir}/nodif_brain_maskLR_slice_*
if [ `${FSLDIR}/bin/imtest ${out_dir}/grad_devLR_slice_0000` -eq 1 ];then
rm -f ${out_dir}/grad_devLR_slice_*
fi
if [ `${FSLDIR}/bin/imtest ${out_dir}/grad_devHR_newslice_0000` -eq 1 ];then
rm -f ${out_dir}/grad_devHR_newslice_*
fi
fi
fi
fi
echo Done
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment