From 9c00435d02e4ccb874a27cd4df3793c66db7a668 Mon Sep 17 00:00:00 2001
From: Matthew Webster <mwebster@fmrib.ox.ac.uk>
Date: Fri, 9 Dec 2011 14:54:29 +0000
Subject: [PATCH] Confounds
---
fsl_glm.cc | 44 ++++++++++++++++++++++++++++++++++++--------
melhlprfns.cc | 7 ++++---
melhlprfns.h | 2 +-
3 files changed, 41 insertions(+), 12 deletions(-)
diff --git a/fsl_glm.cc b/fsl_glm.cc
index f6801cc..3d97261 100644
--- a/fsl_glm.cc
+++ b/fsl_glm.cc
@@ -48,7 +48,7 @@ using namespace std;
Option<string> fnftest(string("-f,--ftests"), string(""),
string("matrix of F-tests on contrasts"),
false, requires_argument,false);
- Option<int> dofset(string("--dof"),0,
+ Option<int> dofset(string("--dof"),-1,
string(" set degrees-of-freedom explicitly"),
false, requires_argument);
Option<bool> normdes(string("--des_norm"),FALSE,
@@ -103,6 +103,9 @@ using namespace std;
Option<string> outvnscales(string("--out_vnscales"),string(""),
string("output file name for scaling factors for variance normalisation"),
false, requires_argument);
+ Option<vector<string> > voxelwiseConfounds(string("--vxf"), vector<string>(),
+ string("\tlist of 4D images containing voxelwise EVs (list order corresponds to numbers in vxl option). caution BETA option."),
+ false, requires_argument);
/*
}
*/
@@ -147,7 +150,7 @@ void saveit(Matrix what, string fname){
write_ascii_matrix(what,fname);
}
-int setup(){
+int setup(int &dof){
if(fsl_imageexists(fnin.value())){//read data
//input is 3D/4D vol
volume4D<float> tmpdata;
@@ -199,10 +202,34 @@ int setup(){
design = read_ascii_matrix(fndesign.value());
}
+ dof=(int)ols_dof(design);
+
+ if ( voxelwiseConfounds.set() ) {
+ vector<Matrix> confounds;
+ confounds.resize(voxelwiseConfounds.value().size());
+ cerr << confounds.size() << endl;
+ volume4D<float> input;
+ for(unsigned int i=0; i< confounds.size(); i++) {
+ read_volume4D(input,voxelwiseConfounds.value().at(i));
+ if ( mask.nvoxels() )
+ confounds.at(i)=input.matrix(mask);
+ else
+ confounds.at(i)=input.matrix();
+ }
+ for(int voxel=1;voxel<=data.Ncols();voxel++) {
+ Matrix confound(confounds.at(0).Column(voxel) );
+ for(unsigned int i=1; i< confounds.size(); i++)
+ confound|=confounds.at(i).Column(voxel);
+ data.Column(voxel)=(IdentityMatrix(confound.Nrows())-confound*pinv(confound))*data.Column(voxel);
+ }
+ dof-=2;
+ }
+
if(perf_demean.value()){
if(debug.value())
cout << "De-meaning the data matrix" << endl;
data = remmean(data,1);
+ dof-=1;
}
if(normdat.value()){
if(debug.value())
@@ -262,12 +289,12 @@ void write_res(){
}
int do_work(int argc, char* argv[]) {
- if(setup())
- exit(1);
-
- glm.olsfit(data,design,contrasts,dofset.value());
- write_res();
- return 0;
+ int dof(-1);
+ if(setup(dof))
+ exit(1);
+ glm.olsfit(data,design,contrasts,dof);
+ write_res();
+ return 0;
}
////////////////////////////////////////////////////////////////////////////
@@ -302,6 +329,7 @@ int main(int argc,char *argv[]){
options.add(outsigsq);
options.add(outdata);
options.add(outvnscales);
+ options.add(voxelwiseConfounds);
options.parse_command_line(argc, argv);
// line below stops the program if the help was requested or
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 0ce17c0..7d5a013 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -884,7 +884,7 @@ namespace Melodic{
} //Matrix gen_arCorr
void basicGLM::olsfit(const Matrix& data, const Matrix& design,
- const Matrix& contrasts, int DOFadjust)
+ const Matrix& contrasts, int requestedDOF)
{
beta = zeros(design.Ncols(),1);
residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1);
@@ -898,8 +898,9 @@ namespace Melodic{
beta = pinvdes * dat;
residu = dat - design*beta;
-
- dof = design.Nrows() - design.Ncols()-1;
+ dof = ols_dof(design);
+ if ( requestedDOF>0)
+ dof = requestedDOF;
sigsq = sum(SP(residu,residu))/dof;
float fact = float(dof) / design.Ncols();
diff --git a/melhlprfns.h b/melhlprfns.h
index 7b4dcbe..5164451 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -86,7 +86,7 @@ namespace Melodic{
~basicGLM(){}
void olsfit(const Matrix& data, const Matrix& design,
- const Matrix& contrasts, int DOFadjust = 0);
+ const Matrix& contrasts, int DOFadjust = -1);
inline Matrix& get_t(){return t;}
inline Matrix& get_z(){return z;}
--
GitLab