From 56439f0dd336ee38fb32418b2f202ddcef88c45e Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Mon, 23 Jul 2007 12:32:04 +0000
Subject: [PATCH] added fsl_glm and fsl_regfilt
---
Makefile | 15 ++-
fsl_glm.cc | 275 +++++++++++++++++++++++++++++++++++++++++
fsl_regfilt.cc | 230 ++++++++++++++++++++++++++++++++++
meldata.cc | 23 +++-
meldata.h | 3 +-
melhlprfns.cc | 29 ++---
melhlprfns.h | 3 +-
melica.cc | 3 +-
meloptions.cc | 2 +-
meloptions.h | 7 +-
melreport.cc | 121 +++++++++++++++---
test.cc | 326 +++++++++++++++++++++++++++++++++++++++++--------
12 files changed, 941 insertions(+), 96 deletions(-)
create mode 100644 fsl_glm.cc
create mode 100644 fsl_regfilt.cc
diff --git a/Makefile b/Makefile
index 3919431..6cc7f4a 100644
--- a/Makefile
+++ b/Makefile
@@ -16,17 +16,21 @@ TEST_OBJS = melhlprfns.o test.o
GGMIX_OBJS = ggmix.o
+FSL_GLM_OBJS = melhlprfns.o fsl_glm.o
+
+FSL_REGFILT_OBJS = melhlprfns.o fsl_regfilt.o
+
MELODIC_OBJS = meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o melreport.o melodic.o
TESTXFILES = test
-XFILES = melodic
+XFILES = fsl_glm fsl_regfilt melodic
RUNTCLS = Melodic
SCRIPTS = melodicreport
-all: ggmix melodic
+all: ggmix fsl_regfilt fsl_glm melodic
ggmix: ${GGMIX_OBJS}
${AR} -r libggmix.a ${GGMIX_OBJS}
@@ -34,8 +38,13 @@ ggmix: ${GGMIX_OBJS}
melodic: ${MELODIC_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${MELODIC_OBJS} ${LIBS}
-
test: ${TEST_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${TEST_OBJS} ${LIBS}
+fsl_glm: ${FSL_GLM_OBJS}
+ $(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_GLM_OBJS} ${LIBS}
+
+fsl_regfilt: ${FSL_REGFILT_OBJS}
+ $(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_REGFILT_OBJS} ${LIBS}
+
diff --git a/fsl_glm.cc b/fsl_glm.cc
new file mode 100644
index 0000000..e7ce0d0
--- /dev/null
+++ b/fsl_glm.cc
@@ -0,0 +1,275 @@
+/* fsl_glm -
+
+ Christian Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 2006-2007 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "libvis/miscplot.h"
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/miscprob.h"
+#include "utils/options.h"
+#include <vector>
+#include "newimage/newimageall.h"
+#include "melhlprfns.h"
+
+using namespace MISCPLOT;
+using namespace MISCMATHS;
+using namespace Utilities;
+using namespace std;
+
+// The two strings below specify the title and example usage that is
+// printed out as the help or usage message
+
+ string title=string("fsl_glm (Version 1.0)")+
+ string("\nCopyright(c) 2007, University of Oxford (Christian F. Beckmann)\n")+
+ string(" \n Simple GLM usign ordinary least-squares regression on\n")+
+ string(" time courses and/or 3D/4D imges against time courses \n")+
+ string(" or 3D/4D images\n\n");
+ string examples="fsl_glm -i <input> -d <design> [options]";
+
+//Command line Options {
+ Option<string> fnin(string("-i,--in"), string(""),
+ string(" input file name (matrix 3D or 4D image)"),
+ true, requires_argument);
+ Option<string> fnout(string("-o,--out"), string(""),
+ string(" output file name for GLM parameter estimates"),
+ true, requires_argument);
+ Option<string> fndesign(string("-d,--design"), string(""),
+ string("file name of the GLM design matrix (time courses or spatial maps)"),
+ true, requires_argument);
+ Option<string> fnmask(string("-m,--mask"), string(""),
+ string("mask image file name"),
+ false, requires_argument);
+ Option<string> fncontrasts(string("-c,--contrasts"), string(""),
+ string("matrix of t-statistics contrasts"),
+ false, requires_argument);
+ Option<string> fnftest(string("-f,--ftests"), string(""),
+ string("matrix of F-tests on contrasts"),
+ false, requires_argument);
+ Option<int> dofset(string("--dof"),0,
+ string(" set degrees-of-freedom explicitly"),
+ false, requires_argument);
+ Option<bool> perfvn(string("--vn"),FALSE,
+ string(" perfrom variance-normalisation on data"),
+ false, requires_argument);
+ Option<int> help(string("-h,--help"), 0,
+ string("display this help text"),
+ false,no_argument);
+ // Output options
+ Option<string> outcope(string("--out_cope"),string(""),
+ string("output COPEs"),
+ false, requires_argument);
+ Option<string> outz(string("--out_z"),string(""),
+ string(" output Z-stats"),
+ false, requires_argument);
+ Option<string> outt(string("--out_t"),string(""),
+ string(" output t-stats"),
+ false, requires_argument);
+ Option<string> outp(string("--out_p"),string(""),
+ string(" output p-values of Z-stats"),
+ false, requires_argument);
+ Option<string> outf(string("--out_f"),string(""),
+ string(" output F-value of full model fit"),
+ false, requires_argument);
+ Option<string> outpf(string("--out_pf"),string(""),
+ string("output p-value for full model fit"),
+ false, requires_argument);
+ Option<string> outres(string("--out_res"),string(""),
+ string("output residuals"),
+ false, requires_argument);
+ Option<string> outvarcb(string("--out_varcb"),string(""),
+ string("output variance of COPEs"),
+ false, requires_argument);
+ Option<string> outsigsq(string("--out_sigsq"),string(""),
+ string("output residual noise variance sigma-square"),
+ false, requires_argument);
+ Option<string> outdata(string("--out_data"),string(""),
+ string("output data"),
+ false, requires_argument);
+ Option<string> outvnscales(string("--out_vnscales"),string(""),
+ string("output scaling factors for variance normalisation"),
+ false, requires_argument);
+ /*
+}
+*/
+//Globals {
+ Melodic::basicGLM glm;
+ int voxels = 0;
+ Matrix data;
+ Matrix design;
+ Matrix contrasts;
+ Matrix fcontrasts;
+ Matrix meanR;
+ RowVector vnscales;
+ volume<float> mask;
+ volumeinfo volinf; /*
+}
+*/
+////////////////////////////////////////////////////////////////////////////
+
+// Local functions
+void save4D(Matrix what, string fname){
+ if(what.Ncols()==data.Ncols()||what.Nrows()==data.Nrows()){
+ volume4D<float> tempVol;
+ if(what.Nrows()>what.Ncols())
+ tempVol.setmatrix(what.t(),mask);
+ else
+ tempVol.setmatrix(what,mask);
+ save_volume4D(tempVol,fname,volinf);
+ }
+}
+
+bool isimage(Matrix what){
+ if((voxels > 0)&&(what.Ncols()==voxels || what.Nrows()==voxels))
+ return TRUE;
+ else
+ return FALSE;
+}
+
+void saveit(Matrix what, string fname){
+ if(isimage(what))
+ save4D(what,fname);
+ else
+ write_ascii_matrix(what,fname);
+}
+
+int setup(){
+ if(fsl_imageexists(fnin.value())){//read data
+ //input is 3D/4D vol
+ volume4D<float> tmpdata;
+ read_volume4D(tmpdata,fnin.value(),volinf);
+
+ // create mask
+ if(fnmask.value()>""){
+ read_volume(mask,fnmask.value());
+ if(!samesize(tmpdata[0],mask)){
+ cerr << "ERROR: Mask image does not match input image" << endl;
+ return 1;
+ };
+ }else{
+ mask = tmpdata[0]*0.0+1.0;
+ }
+
+ data = tmpdata.matrix(mask);
+ voxels = data.Ncols();
+ }
+ else
+ data = read_ascii_matrix(fnin.value());
+
+ if(fsl_imageexists(fndesign.value())){//read design
+ volume4D<float> tmpdata;
+ read_volume4D(tmpdata,fndesign.value());
+ if(!samesize(tmpdata[0],mask)){
+ cerr << "ERROR: GLM design does not match input image in size" << endl;
+ return 1;
+ }
+ design = tmpdata.matrix(mask).t();
+ data = data.t();
+ }else{
+ design = read_ascii_matrix(fndesign.value());
+ }
+
+ meanR=mean(data,1);
+ data = remmean(data,1);
+ design = remmean(design,1);
+ if(perfvn.value())
+ vnscales = Melodic::varnorm(data);
+ if(fncontrasts.value()>""){//read contrast
+ contrasts = read_ascii_matrix(fncontrasts.value());
+ if(!(contrasts.Ncols()==design.Ncols())){
+ cerr << "ERROR: contrast matrix GLM design does not match GLM design" << endl;
+ return 1;
+ }
+ }else{
+ contrasts = Identity(design.Ncols());
+ contrasts &= -1.0 * contrasts;
+ }
+ return 0;
+}
+
+void write_res(){
+ if(fnout.value()>"")
+ saveit(glm.get_beta(),fnout.value());
+ if(outcope.value()>"")
+ saveit(glm.get_cbeta(),outcope.value());
+ if(outz.value()>"")
+ saveit(glm.get_z(),outz.value());
+ if(outt.value()>"")
+ saveit(glm.get_t(),outt.value());
+ if(outp.value()>"")
+ saveit(glm.get_p(),outp.value());
+ if(outf.value()>"")
+ saveit(glm.get_f_fmf(),outf.value());
+ if(outpf.value()>"")
+ saveit(glm.get_pf_fmf(),outpf.value());
+ if(outres.value()>"")
+ saveit(glm.get_residu(),outres.value());
+ if(outvarcb.value()>"")
+ saveit(glm.get_varcb(),outvarcb.value());
+ if(outsigsq.value()>"")
+ saveit(glm.get_sigsq(),outsigsq.value());
+ if(outdata.value()>"")
+ saveit(data,outdata.value());
+ if(outvnscales.value()>"")
+ saveit(vnscales,outvnscales.value());
+}
+
+int do_work(int argc, char* argv[]) {
+ if(setup())
+ exit(1);
+
+ glm.olsfit(data,design,contrasts,dofset.value());
+ write_res();
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+int main(int argc,char *argv[]){
+ Tracer tr("main");
+ OptionParser options(title, examples);
+ try{
+ // must include all wanted options here (the order determines how
+ // the help message is printed)
+ options.add(fnin);
+ options.add(fnout);
+ options.add(fndesign);
+ options.add(fnmask);
+ options.add(fncontrasts);
+ options.add(fnftest);
+ options.add(dofset);
+ options.add(perfvn);
+ options.add(help);
+ options.add(outcope);
+ options.add(outz);
+ options.add(outt);
+ options.add(outp);
+ options.add(outf);
+ options.add(outpf);
+ options.add(outres);
+ options.add(outvarcb);
+ options.add(outsigsq);
+ options.add(outdata);
+ options.add(outvnscales);
+ options.parse_command_line(argc, argv);
+
+ // line below stops the program if the help was requested or
+ // a compulsory option was not set
+ if ( (help.value()) || (!options.check_compulsory_arguments(true)) ){
+ options.usage();
+ exit(EXIT_FAILURE);
+ }else{
+ // Call the local functions
+ return do_work(argc,argv);
+ }
+ }catch(X_OptionError& e) {
+ options.usage();
+ cerr << endl << e.what() << endl;
+ exit(EXIT_FAILURE);
+ }catch(std::exception &e) {
+ cerr << e.what() << endl;
+ }
+ }
+
diff --git a/fsl_regfilt.cc b/fsl_regfilt.cc
new file mode 100644
index 0000000..e8f4c1c
--- /dev/null
+++ b/fsl_regfilt.cc
@@ -0,0 +1,230 @@
+/* fsl_regfilt -
+
+ Christian Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 2006-2007 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "libvis/miscplot.h"
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/miscprob.h"
+#include "utils/options.h"
+#include <vector>
+#include "newimage/newimageall.h"
+#include "melhlprfns.h"
+
+using namespace MISCPLOT;
+using namespace MISCMATHS;
+using namespace Utilities;
+using namespace std;
+
+// The two strings below specify the title and example usage that is
+// printed out as the help or usage message
+
+ string title=string("fsl_regfilt (Version 1.0)")+
+ string("\nCopyright(c) 2007, University of Oxford (Christian F. Beckmann)\n")+
+ string(" Data filtering by regressing out part of a design matrix\n \n")+
+ string(" using simple OLS regression on 4D images\n\n");
+ string examples="fsl_regfilt -i <input> -d <design> -f -o <out> [options]";
+
+//Command line Options {
+ Option<string> fnin(string("-i,--in"), string(""),
+ string(" input file name (4D image)"),
+ true, requires_argument);
+ Option<string> fnout(string("-o,--out"), string(""),
+ string(" output file name for the filtered data"),
+ true, requires_argument);
+ Option<string> fndesign(string("-d,--design"), string(""),
+ string("file name of the GLM design matrix (time courses)"),
+ true, requires_argument);
+ Option<string> fnmask(string("-m,--mask"), string(""),
+ string("mask image file name"),
+ false, requires_argument);
+ Option<string> filter(string("-f,--filter"),string(""),
+ string("filter out part of the regression model"),
+ true, requires_argument);
+ Option<bool> perfvn(string("--vn"),FALSE,
+ string(" perfrom variance-normalisation on data"),
+ false, requires_argument);
+ Option<int> help(string("-h,--help"), 0,
+ string("display this help text"),
+ false,no_argument);
+ // Output options
+ Option<string> outdata(string("--out_data"),string(""),
+ string("output data"),
+ false, requires_argument);
+ Option<string> outvnscales(string("--out_vnscales"),string(""),
+ string("output scaling factors for variance normalisation"),
+ false, requires_argument);
+ /*
+}
+*/
+//Globals {
+ int voxels = 0;
+ Matrix data;
+ Matrix design;
+ Matrix meanR;
+ RowVector vnscales;
+ volume<float> mask;
+ volumeinfo volinf; /*
+}
+*/
+////////////////////////////////////////////////////////////////////////////
+
+// Local functions
+void save4D(Matrix what, string fname){
+ if(what.Ncols()==data.Ncols()||what.Nrows()==data.Nrows()){
+ volume4D<float> tempVol;
+ if(what.Nrows()>what.Ncols())
+ tempVol.setmatrix(what.t(),mask);
+ else
+ tempVol.setmatrix(what,mask);
+ save_volume4D(tempVol,fname,volinf);
+ }
+}
+
+bool isimage(Matrix what){
+ if((voxels > 0)&&(what.Ncols()==voxels || what.Nrows()==voxels))
+ return TRUE;
+ else
+ return FALSE;
+}
+
+void saveit(Matrix what, string fname){
+ if(isimage(what))
+ save4D(what,fname);
+ else
+ write_ascii_matrix(what,fname);
+}
+
+int dofilter(){
+ if(!isimage(data)){
+ cerr << "ERROR: need to specify 4D input to use filtering" << endl;
+ return 1;
+ }
+ Matrix unmixMatrix = pinv(design);
+ Matrix maps = unmixMatrix * data;
+
+ Matrix noisedes;
+ Matrix noisemaps;
+
+ int ctr=0;
+ char *p;
+ char t[1024];
+ const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
+
+ strcpy(t, filter.value().c_str());
+ p=strtok(t,discard);
+ ctr = atoi(p);
+ if(ctr>0 && ctr<=design.Ncols()){
+ noisedes = design.Column(ctr);
+ noisemaps = maps.Row(ctr).t();
+ }
+
+ do{
+ p=strtok(NULL,discard);
+ if(p){
+ ctr = atoi(p);
+ if(ctr>0 && ctr<=design.Ncols()){
+ noisedes |= design.Column(ctr);
+ noisemaps |= maps.Row(ctr).t();
+ }
+ }
+ }while(p);
+ Matrix newData;
+ newData = data - noisedes * noisemaps.t();
+ newData = newData + ones(newData.Nrows(),1)*meanR;
+
+ save4D(newData,fnout.value());
+ return 0;
+}
+
+int setup(){
+ if(fsl_imageexists(fnin.value())){//read data
+ //input is 3D/4D vol
+ volume4D<float> tmpdata;
+ read_volume4D(tmpdata,fnin.value(),volinf);
+
+ // create mask
+ if(fnmask.value()>""){
+ read_volume(mask,fnmask.value());
+ if(!samesize(tmpdata[0],mask)){
+ cerr << "ERROR: Mask image does not match input image" << endl;
+ return 1;
+ };
+ }else{
+ mask = tmpdata[0]*0.0+1.0;
+ }
+
+ data = tmpdata.matrix(mask);
+ voxels = data.Ncols();
+ }else{
+ cerr << "ERROR: cannot read input image " << fnin.value()<<endl;
+ return 1;
+ }
+
+ design = read_ascii_matrix(fndesign.value());
+
+ meanR=mean(data,1);
+ data = remmean(data,1);
+ design = remmean(design,1);
+ if(perfvn.value())
+ vnscales = Melodic::varnorm(data);
+ return 0;
+}
+
+void write_res(){
+ if(outdata.value()>"")
+ saveit(data,outdata.value());
+ if(outvnscales.value()>"")
+ saveit(vnscales,outvnscales.value());
+}
+
+int do_work(int argc, char* argv[]) {
+ if(setup())
+ exit(1);
+
+ if(dofilter())
+ exit(1);
+ write_res();
+ return 0;
+}
+
+////////////////////////////////////////////////////////////////////////////
+
+int main(int argc,char *argv[]){
+ Tracer tr("main");
+ OptionParser options(title, examples);
+ try{
+ // must include all wanted options here (the order determines how
+ // the help message is printed)
+ options.add(fnin);
+ options.add(fnout);
+ options.add(fndesign);
+ options.add(fnmask);
+ options.add(filter);
+ options.add(perfvn);
+ options.add(help);
+ options.add(outdata);
+ options.add(outvnscales);
+ options.parse_command_line(argc, argv);
+
+ // line below stops the program if the help was requested or
+ // a compulsory option was not set
+ if ( (help.value()) || (!options.check_compulsory_arguments(true)) ){
+ options.usage();
+ exit(EXIT_FAILURE);
+ }else{
+ // Call the local functions
+ return do_work(argc,argv);
+ }
+ }catch(X_OptionError& e) {
+ options.usage();
+ cerr << endl << e.what() << endl;
+ exit(EXIT_FAILURE);
+ }catch(std::exception &e) {
+ cerr << e.what() << endl;
+ }
+ }
+
diff --git a/meldata.cc b/meldata.cc
index 5a2f43c..89fffd3 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -140,9 +140,27 @@ namespace Melodic{
for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
tmpT2 << tmpT.Column(ctr);
tmpS2 << tmpS.Column(ctr);
+ if(mean(tmpS2,1).AsScalar()<0){
+ tmpT2*=-1.0;
+ tmpS2*=-1.0;
+ }
add_Tmodes(tmpT2);
add_Smodes(tmpS2);
}
+
+ //add GLM OLS fit
+ if(Tdes.Storage()){
+ Matrix alltcs = Tmodes.at(0);
+ for(int ctr=1; ctr < (int)Tmodes.size();ctr++)
+ alltcs|=Tmodes.at(ctr);
+ glmT.olsfit(alltcs,Tdes,Tcon);
+ }
+ if(Sdes.Storage()){
+ Matrix alltcs = Smodes.at(0);
+ for(int ctr=1; ctr < (int)Smodes.size();ctr++)
+ alltcs|=Smodes.at(ctr);
+ glmS.olsfit(alltcs,Sdes,Scon);
+ }
}
void MelodicData::setup()
@@ -668,10 +686,13 @@ namespace Melodic{
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
numTime = mixMatrix.Nrows(), i,j;
+ //flip IC maps to be positive (on average)
+ //flip Subject/Session modes to be positive (on average)
+ //have time courses accordingly
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++)
if(IC.Row(ctr_i).Sum()<0)
flipres(ctr_i);
-
+
// re-order wrt standard deviation of IC maps
message("Sorting IC maps" << endl);
Matrix tmpscales, tmpICrow, tmpMIXcol;
diff --git a/meldata.h b/meldata.h
index 6231469..4ebec20 100644
--- a/meldata.h
+++ b/meldata.h
@@ -183,6 +183,7 @@ namespace Melodic{
volumeinfo tempInfo;
vector<Matrix> DWM, WM;
basicGLM glmT, glmS;
+ Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;
private:
MelodicOptions &opts;
@@ -202,8 +203,6 @@ namespace Melodic{
Matrix mixFFT;
Matrix IC;
Matrix ICstats;
- Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;
-
vector<Matrix> Tmodes;
vector<Matrix> Smodes;
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 6c88188..05a73d4 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -499,7 +499,6 @@ namespace Melodic{
return res;
} //int ppca_dim
-
int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which)
{
ColumnVector PPCA;
@@ -783,7 +782,7 @@ namespace Melodic{
return res;
} //Matrix gen_arCorr
- Matrix basicGLM::olsfit(const Matrix& data, const Matrix& design,
+ void basicGLM::olsfit(const Matrix& data, const Matrix& design,
const Matrix& contrasts, int DOFadjust)
{
beta = zeros(design.Ncols(),1);
@@ -793,28 +792,30 @@ namespace Melodic{
if(data.Nrows()==design.Nrows()){
Matrix dat = remmean(data,1);
- Matrix pinvdes = pinv(design);
+ Matrix tmp = design.t()*design;
+ Matrix pinvdes = tmp.i()*design.t();
beta = pinvdes * dat;
residu = dat - design*beta;
- Matrix R = Identity(design.Nrows()) - design * pinvdes;
- Matrix R2 = R*R;
- float tR = R.Trace();
- sigsq = sum(SP(residu,residu))/tR;
+// Matrix R = Identity(design.Nrows()) - design * pinvdes;
+// Matrix R2 = R*R;
+// float tR = R.Trace();
+// sigsq = sum(SP(residu,residu))/tR;
- dof = (int)(tR*tR/R2.Trace() - DOFadjust);
+// dof = (int)(tR*tR/R2.Trace() - DOFadjust);
+ dof = design.Nrows() - design.Ncols();
+ sigsq = sum(SP(residu,residu))/dof;
float fact = float(design.Nrows() - 1 - design.Ncols()) / design.Ncols();
f_fmf = SP(var(design*beta),pow(var(residu),-1))* fact;
- pf_fmf = f_fmf.Row(1); pf_fmf &= pf_fmf;
+ pf_fmf = f_fmf.Row(1);
for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),
int(design.Nrows() -1 -design.Ncols()),f_fmf.Column(ctr1).AsScalar());
- pf_fmf.Row(2) = pf_fmf.Row(1) * pf_fmf.Ncols();
- if(contrasts.Storage()>0 && contrasts.Nrows()==beta.Ncols()){
- cbeta = contrasts.t()*beta;
- Matrix tmp = contrasts.t()*pinvdes*pinvdes.t()*contrasts;
+ if(contrasts.Storage()>0 && contrasts.Ncols()==beta.Nrows()){
+ cbeta = contrasts*beta;
+ Matrix tmp = contrasts*pinvdes*pinvdes.t()*contrasts.t();
varcb = diag(tmp)*sigsq;
t = SP(cbeta,pow(varcb,-0.5));
z = t; p=t;
@@ -827,7 +828,7 @@ namespace Melodic{
}
}
}
- return beta;
+
}
diff --git a/melhlprfns.h b/melhlprfns.h
index e0334cf..6f8d5a4 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -80,7 +80,7 @@ namespace Melodic{
//destructor
~basicGLM(){}
- Matrix olsfit(const Matrix& data, const Matrix& design,
+ void olsfit(const Matrix& data, const Matrix& design,
const Matrix& contrasts, int DOFadjust = 0);
inline Matrix& get_t(){return t;}
@@ -89,6 +89,7 @@ namespace Melodic{
inline Matrix& get_f_fmf(){return f_fmf;}
inline Matrix& get_pf_fmf(){return pf_fmf;}
inline Matrix& get_cbeta(){return cbeta;}
+ inline Matrix& get_beta(){return beta;}
inline Matrix& get_varcb(){return varcb;}
inline Matrix& get_sigsq(){return sigsq;}
inline Matrix& get_residu(){return residu;}
diff --git a/melica.cc b/melica.cc
index 253c620..ee4fb3f 100644
--- a/melica.cc
+++ b/melica.cc
@@ -102,7 +102,7 @@ namespace Melodic{
cum_itt += itt_ctr;
itt_ctr2++;
if(opts.approach.value() == string("tica")){
- message(" Rank-1 approximation of the time courses; ");
+ message(" Rank-1 approximation of the time courses; " <<endl;);
Matrix temp(melodat.get_dewhite() * redUMM);
outMsize(" 2nd unmmixing matrix ", temp);
temp = melodat.expand_dimred(temp);
@@ -459,7 +459,6 @@ namespace Melodic{
melodat.set_IC(temp);
melodat.set_ICstats(scales);
melodat.sort();
-
melodat.set_TSmode();
}
diff --git a/meloptions.cc b/meloptions.cc
index c534053..5742c17 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -157,7 +157,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
}
//in the case of indirect inputs, create the vector of input names here
- if( inputindirect.value()){
+ if(!fsl_imageexists(inputfname.value().at(0))){
std::vector< string > tmpfnames;
ifstream fs(inputfname.value().at(0).c_str());
string cline;
diff --git a/meloptions.h b/meloptions.h
index d1c11fa..81d7bc5 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -46,7 +46,6 @@ class MelodicOptions {
Option<string> maskfname;
Option<bool> use_mask;
- Option<bool> inputindirect;
Option<bool> update_mask;
Option<bool> perf_bet;
Option<float> threshold;
@@ -147,7 +146,7 @@ class MelodicOptions {
string("output directory name\n"),
false, requires_argument),
inputfname(string("-i,--in"), std::vector<string>(),
- string("input file names (either single file name or comma-separated list)"),
+ string("input file names (either single file name or comma-separated list or text file)"),
true, requires_argument),
outputfname(string("-O,--out"), string("melodic"),
string("output file name"),
@@ -158,9 +157,6 @@ class MelodicOptions {
use_mask(string("--nomask"), true,
string("switch off masking"),
false, no_argument),
- inputindirect(string("--filelist"), false,
- string("input file contains list of input file names"),
- false, no_argument),
update_mask(string("--update_mask"), true,
string("switch off mask updating"),
false, no_argument),
@@ -351,7 +347,6 @@ class MelodicOptions {
options.add(guessfname);
options.add(maskfname);
options.add(use_mask);
- options.add(inputindirect);
options.add(update_mask);
options.add(perf_bet);
options.add(threshold);
diff --git a/melreport.cc b/melreport.cc
index 3fe7994..3597b78 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -108,23 +108,17 @@ namespace Melodic{
}
{//plot time course
- IChtml << "<H3> Temporal mode <p>" << endl <<endl;
+ IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
miscplot newplot;
Matrix tmptc = melodat.get_Tmodes(cnum-1).t();
-
+
//add GLM OLS fit
- /* basicGLM glm;
- if(melodat.Tdes.Storage()){
- Matrix betas = glm.olsfit(tmptc.t(),melodat.Tdes,melodat.Tcon).t();
- tmptc &= betas*melodat.Tdes.t();
+ if(melodat.Tdes.Storage() > 0){
+ tmptc &= melodat.glmT.get_beta().Column(cnum).t() * melodat.Tdes.t();
newplot.add_label(string("IC ")+num2str(cnum)+" time course");
newplot.add_label("full model fit");
-
-cerr << endl << endl <<
-glm.get_f_fmf() << endl<<
-glm.get_pf_fmf() << endl << endl;
}
-*/
+
if(opts.tr.value()>0.0)
newplot.timeseries(tmptc,
report.appendDir(string("t")+num2str(cnum)+".png"),
@@ -140,10 +134,9 @@ glm.get_pf_fmf() << endl << endl;
IChtml << "<A HREF=\"" << string("t")
+num2str(cnum)+".txt" << "\"> ";
IChtml << "<img BORDER=0 SRC=\""
- +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
- }//time series plot
-
- {//plot frequency
+ +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
+ }//time series plot
+ {//plot frequency
miscplot newplot;
RowVector empty(1);
empty = 0.0;
@@ -181,17 +174,74 @@ glm.get_pf_fmf() << endl << endl;
IChtml << "<img BORDER=0 SRC=\""
+string("f")+num2str(cnum)+".png\"></A><p>" << endl;
}//frequency plot
-
+ {//add T-mode GLM F-stats for full model fit & contrasts
+ if(melodat.Tdes.Storage() > 0){
+ IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
+ "<CAPTION><EM> <b>GLM (OLS) on time series </b></EM></CAPTION>" << endl
+ << "<TR valign=middle><TH ><EM>GLM β's</EM> <TH> <EM> F-test on <br> full model fit </em>";
+ if(melodat.Tcon.Storage() > 0)
+ IChtml << "<TH ><EM>Contrasts</EM>"<<endl;
+ IChtml << "<TR><TD><TABLE border=0><TR><TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Tdes.Ncols();ctr++)
+ IChtml << " PE(" <<num2str(ctr)+"): <br>" << endl;
+ IChtml << "<TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Tdes.Ncols();ctr++)
+ IChtml << melodat.glmT.get_beta().Column(cnum).Row(ctr) << "<br>" <<endl;
+ IChtml << "</TABLE>" <<
+ " <TD align=center> F = "<< melodat.glmT.get_f_fmf().Column(cnum) <<
+ " <BR> dof1 = " << melodat.Tdes.Ncols() << "; dof2 = "
+ << melodat.glmT.get_dof() << "<BR>" <<endl;
+ if(melodat.glmT.get_pf_fmf().Column(cnum).AsScalar() < 0.05)
+ IChtml << "<b> p < " << melodat.glmT.get_pf_fmf().Column(cnum) <<
+ "<BR> (uncorrected for #comp.)<b></TD>" << endl;
+ else
+ IChtml << " p < " <<
+ melodat.glmT.get_pf_fmf().Column(cnum) <<
+ "<BR> (uncorrected for #comp.)</TD>" << endl;
+ }
+ if(melodat.Tcon.Storage() > 0){
+ IChtml << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
+ for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
+ IChtml << "con(" << melodat.Tcon.Row(ctr) << "): <br>" << endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
+ IChtml <<" z = <BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
+ IChtml << melodat.glmT.get_z().Column(cnum).Row(ctr) <<";<BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
+ if(melodat.glmT.get_p().Column(cnum).Row(ctr).AsScalar() < 0.05)
+ IChtml << "<b> p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
+ "</b><BR>" << endl;
+ else
+ IChtml << " p < " << melodat.glmT.get_p().Column(cnum).Row(ctr) <<
+ "<BR>" << endl;
+ IChtml << "</TABLE></td></tr>" << endl;
+ }
+ IChtml << "</TABLE><p>" << endl;
+ }
+
if(cnum <= (int)melodat.get_Smodes().size())
{//plot subject mode
Matrix smode;
smode = melodat.get_Smodes(cnum-1);
+
if(smode.Nrows() > 1){
miscplot newplot;
+
+ //add GLM OLS fit
+ if(melodat.Sdes.Storage() > 0){
+ smode |= melodat.Sdes * melodat.glmS.get_beta().Column(cnum);
+ newplot.add_label(string("IC ")+num2str(cnum)+" subject/session-mode");
+ newplot.add_label("full model fit");
+ }
newplot.setscatter(smode,5);
newplot.timeseries(smode.t(),
report.appendDir(string("s")+num2str(cnum)+".png"),
string("Subject/Session mode"));
+ newplot.set_xysize(120,200);
+ newplot.set_minmaxscale(1.1);
newplot.boxplot(smode,
report.appendDir(string("b")+num2str(cnum)+".png"),
string("Subject/Session mode"));
@@ -206,8 +256,43 @@ glm.get_pf_fmf() << endl << endl;
IChtml << "<img BORDER=0 SRC=\""
+string("b")+num2str(cnum)+".png\"></A><p>" << endl;
}
- }//subject mode plot
-
+ {//add S-mode GLM F-stats for full model fit & contrasts
+ if(melodat.Sdes.Storage() > 0){
+ IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
+ "<CAPTION><EM> <b>GLM (OLS) on subject/session-mode </b></EM></CAPTION>" << endl
+ << "<TR valign=middle><TH colspan=2>Betas <TH> <EM> F-test on <br> full model fit </em>";
+ if(melodat.Scon.Storage() > 0)
+ IChtml << "<TH colspan=3><EM>Contrasts</EM>"<<endl;
+ IChtml << "<TR><TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
+ IChtml << " β(" <<num2str(ctr)+"): <br>" << endl;
+ IChtml << "<TD align=right>" << endl;
+ for(int ctr=1;ctr <= melodat.Sdes.Ncols();ctr++)
+ IChtml << melodat.glmS.get_beta().Column(cnum).Row(ctr) << "<br>" <<endl;
+ IChtml <<
+ " <TD align=center> F = "<< melodat.glmS.get_f_fmf().Column(cnum) <<
+ " <BR> dof1 = " << melodat.Sdes.Ncols() << "; dof2 = "
+ << melodat.glmS.get_dof() << "<BR> p < " <<
+ melodat.glmS.get_pf_fmf().Column(cnum) << "</TD>" << endl;
+ }
+ if(melodat.Scon.Storage() > 0){
+ IChtml << "<td>" <<endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml << "con(" << melodat.Scon.Row(ctr) << ") <br>" << endl;
+ IChtml << "<td align=center>" << endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml << " z = " << melodat.glmS.get_z().Column(cnum).Row(ctr) <<
+ "<BR>" <<endl;
+ IChtml << "<td align=right>" << endl;
+ for(int ctr=1; ctr <= melodat.Scon.Nrows() ; ctr++)
+ IChtml << " p < " << melodat.glmS.get_p().Column(cnum).Row(ctr) <<
+ "<BR>" << endl;
+ IChtml << "</td></tr>" << endl;
+ }
+ IChtml << "</TABLE><p>" << endl;
+ }
+ }//subject mode plot
+
if(mmodel.get_threshmaps().Storage()>0&&
(mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&&
(mmodel.get_threshmaps().Nrows()>1))
diff --git a/test.cc b/test.cc
index 1077016..9423bad 100644
--- a/test.cc
+++ b/test.cc
@@ -9,83 +9,313 @@
#include "miscmaths/miscprob.h"
#include "utils/options.h"
#include <vector>
-
+#include "newimage/newimageall.h"
+#include "melhlprfns.h"
+
using namespace MISCPLOT;
using namespace MISCMATHS;
using namespace Utilities;
using namespace std;
// The two strings below specify the title and example usage that is
-// printed out as the help or usage message
-
-string title="test (Version 1.0)\nCopyright(c) 2007, University of Oxford (Christian F. Beckmann)";
-string examples="test int";
+// printed out as the help or usage message
-Option<bool> help(string("--help"), false,
- string(" display this message"),
- false, no_argument);
-Option<int> num(string("--num"), 1,
- string("number of iterations"),
- false, requires_argument);
-
-int nonoptarg;
+ string title=string("fsl_glm (Version 1.0)")+
+ string("\nCopyright(c) 2007, University of Oxford (Christian F. Beckmann)\n")+
+ string(" \n Simple GLM usign ordinary least-squares regression on\n")+
+ string(" time courses and/or 3D/4D volumes\n\n");
+ string examples="fsl_glm <input> -d <design> [options]";
+//Command line Options {
+ Option<string> fnin(string("-i,--in"), string(""),
+ string("input file name (matrix 3D or 4D image)"),
+ true, requires_argument);
+ Option<string> fnout(string("-o,--out"), string(""),
+ string(""),
+ true, requires_argument);
+ Option<string> fndesign(string("-d,--design"), string(""),
+ string("file name of the GLM design matrix (time courses or spatial maps)"),
+ true, requires_argument);
+ Option<string> fnmask(string("-m,--mask"), string(""),
+ string("mask image"),
+ false, requires_argument);
+ Option<string> fncontrasts(string("-c,--contrasts"), string(""),
+ string("matrix of t-statistics contrasts"),
+ false, requires_argument);
+ Option<string> fnftest(string("-f,--ftests"), string(""),
+ string("matrix of F-tests on contrasts"),
+ false, requires_argument);
+ Option<int> dofset(string("--dof"),0,
+ string("set degrees-of-freedom explicitly"),
+ false, requires_argument);
+ Option<string> filter(string("--filter"),string(""),
+ string("filter out part of the regression model"),
+ false, requires_argument);
+ Option<bool> perfvn(string("--vn"),FALSE,
+ string("perfrom variance-normalisation on data"),
+ false, requires_argument);
+ Option<int> help(string("-h,--help"), 0,
+ string("display this help text"),
+ false,no_argument);
+ // Output options
+ Option<string> outcope(string("--out_cope"),string(""),
+ string("output COPEs"),
+ false, requires_argument);
+ Option<string> outz(string("--out_z"),string(""),
+ string("output Z-stats"),
+ false, requires_argument);
+ Option<string> outt(string("--out_t"),string(""),
+ string("output t-stats"),
+ false, requires_argument);
+ Option<string> outp(string("--out_p"),string(""),
+ string("output p-values of Z-stats"),
+ false, requires_argument);
+ Option<string> outf(string("--out_f"),string(""),
+ string("output F-value of full model fit"),
+ false, requires_argument);
+ Option<string> outpf(string("--out_pf"),string(""),
+ string("output p-value for full model fit"),
+ false, requires_argument);
+ Option<string> outfilt(string("--out_filt"),string(""),
+ string("output filtered data"),
+ false, requires_argument);
+ Option<string> outres(string("--out_res"),string(""),
+ string("output residuals"),
+ false, requires_argument);
+ Option<string> outvarcb(string("--out_varcb"),string(""),
+ string("output variance of COPEs"),
+ false, requires_argument);
+ Option<string> outsigsq(string("--out_sigsq"),string(""),
+ string("output residual noise variance sigma-square"),
+ false, requires_argument);
+ Option<string> outdata(string("--out_data"),string(""),
+ string("output data"),
+ false, requires_argument);
+ Option<string> outvnscales(string("--out_vnscales"),string(""),
+ string("output scaling factors for variance normalisation"),
+ false, requires_argument);
+ /*
+}
+*/
+//Globals {
+ Melodic::basicGLM glm;
+ bool isimage = FALSE;
+ Matrix data;
+ Matrix design;
+ Matrix contrasts;
+ Matrix fcontrasts;
+ Matrix meanR;
+ RowVector vnscales;
+ volume<float> mask;
+ volumeinfo volinf; /*
+}
+*/
////////////////////////////////////////////////////////////////////////////
// Local functions
+void save4D(Matrix what, string fname){
+ if(what.Ncols()==data.Ncols()||what.Nrows()==data.Nrows()){
+ volume4D<float> tempVol;
+ if(what.Nrows()>what.Ncols())
+ tempVol.setmatrix(what.t(),mask);
+ else
+ tempVol.setmatrix(what,mask);
+ save_volume4D(tempVol,fname,volinf);
+ }
+}
+
+void saveit(Matrix what, string fname){
+ if(isimage)
+ save4D(what,fname);
+ else
+ write_ascii_matrix(what,fname);
+}
+
+int dofilter(){
+ if(!isimage){
+ cerr << "ERROR: need to specify 4D input to use filtering" << endl;
+ return 1;
+ }
+ Matrix unmixMatrix = pinv(design);
+ Matrix maps = unmixMatrix * data;
+
+ Matrix noisedes;
+ Matrix noisemaps;
+
+ int ctr=0;
+ char *p;
+ char t[1024];
+ const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
+
+ strcpy(t, filter.value().c_str());
+ p=strtok(t,discard);
+ ctr = atoi(p);
+ if(ctr>0 && ctr<=design.Ncols()){
+ noisedes = design.Column(ctr);
+ noisemaps = maps.Row(ctr).t();
+ }
+
+ do{
+ p=strtok(NULL,discard);
+ if(p){
+ ctr = atoi(p);
+ if(ctr>0 && ctr<=design.Ncols()){
+ noisedes |= design.Column(ctr);
+ noisemaps |= maps.Row(ctr).t();
+ }
+ }
+ }while(p);
+ Matrix newData;
+ newData = data - noisedes * noisemaps.t();
+ newData = newData + ones(newData.Nrows(),1)*meanR;
+
+ save4D(newData,outfilt.value());
+ return 0;
+}
+
+int setup(){
+ if(fsl_imageexists(fnin.value())){//read data
+ //input is 3D/4D vol
+ isimage = TRUE;
+ volume4D<float> tmpdata;
+ read_volume4D(tmpdata,fnin.value(),volinf);
+
+ // create mask
+ if(fnmask.value()>""){
+ read_volume(mask,fnmask.value());
+ if(!samesize(tmpdata[0],mask)){
+ cerr << "ERROR: Mask image does not match input image" << endl;
+ return 1;
+ };
+ }else{
+ mask = tmpdata[0]*0.0+1.0;
+ }
+
+ data = tmpdata.matrix(mask);
+ }
+ else
+ data = read_ascii_matrix(fnin.value());
+
+ if(fsl_imageexists(fndesign.value())){//read design
+ volume4D<float> tmpdata;
+ read_volume4D(tmpdata,fndesign.value());
+ if(!samesize(tmpdata[0],mask)){
+ cerr << "ERROR: GLM design does not match input image in size" << endl;
+ return 1;
+ }
+ design = tmpdata.matrix(mask).t();
+ data = data.t();
+ isimage = FALSE;
+ }else{
+ design = read_ascii_matrix(fndesign.value());
+ }
+
+ meanR=mean(data,1);
+ data = remmean(data,1);
+ design = remmean(design,1);
+ if(perfvn.value())
+ vnscales = Melodic::varnorm(data);
+ if(fncontrasts.value()>""){//read contrast
+ contrasts = read_ascii_matrix(fncontrasts.value());
+ if(!(contrasts.Ncols()==design.Ncols())){
+ cerr << "ERROR: contrast matrix GLM design does not match GLM design" << endl;
+ return 1;
+ }
+ }else{
+ contrasts = Identity(design.Ncols());
+ contrasts &= -1.0 * contrasts;
+ }
+ return 0;
+}
+
+void write_res(){
+ if(fnout.value()>"")
+ saveit(glm.get_beta(),fnout.value());
+ if(outcope.value()>"")
+ saveit(glm.get_cbeta(),outcope.value());
+ if(outz.value()>"")
+ saveit(glm.get_z(),outz.value());
+ if(outt.value()>"")
+ saveit(glm.get_t(),outt.value());
+ if(outp.value()>"")
+ saveit(glm.get_p(),outp.value());
+ if(outf.value()>"")
+ saveit(glm.get_f_fmf(),outf.value());
+ if(outpf.value()>"")
+ saveit(glm.get_pf_fmf(),outpf.value());
+ if(outres.value()>"")
+ saveit(glm.get_residu(),outres.value());
+ if(outvarcb.value()>"")
+ saveit(glm.get_varcb(),outvarcb.value());
+ if(outsigsq.value()>"")
+ saveit(glm.get_sigsq(),outsigsq.value());
+ if(outdata.value()>"")
+ saveit(data,outdata.value());
+ if(outvnscales.value()>"")
+ saveit(vnscales,outvnscales.value());
+}
+
int do_work(int argc, char* argv[]) {
-
- Matrix mat;
+ if(setup())
+ exit(1);
- cout << "BLAH " << num.value() << endl;
- mat=normrnd(300,1);
- miscplot::Timeseries(mat.t(),string("test0.png"),string("TEST"));
-
- for (int i=1; i <= (int)num.value();i++){
- cout << "Processing " << i << endl;
- miscplot newplot;
- newplot.GDCglobals_set();
- mat = normrnd(300,3)+2;
- Matrix col;
- col = mat;
- newplot.add_bpdata(col);
- // newplot.add_bpdata(col);
- newplot.boxplot(string("test")+num2str(i)+string(".png"),string("TEST"));
- }
- return 0;
+ if(filter.value()>"")
+ if(dofilter())
+ exit(1);
+ else{
+ glm.olsfit(data,design,contrasts,dofset.value());
+ write_res();
+ }
+ return 0;
}
////////////////////////////////////////////////////////////////////////////
int main(int argc,char *argv[]){
-
Tracer tr("main");
OptionParser options(title, examples);
-
- try {
+ try{
// must include all wanted options here (the order determines how
// the help message is printed)
- options.add(help);
- options.add(num);
+ options.add(fnin);
+ options.add(fnout);
+ options.add(fndesign);
+ options.add(fnmask);
+ options.add(fncontrasts);
+ options.add(fnftest);
+ options.add(dofset);
+ options.add(filter);
+ options.add(perfvn);
+ options.add(help);
+ options.add(outcope);
+ options.add(outz);
+ options.add(outt);
+ options.add(outp);
+ options.add(outf);
+ options.add(outpf);
+ options.add(outfilt);
+ options.add(outres);
+ options.add(outvarcb);
+ options.add(outsigsq);
+ options.add(outdata);
+ options.add(outvnscales);
options.parse_command_line(argc, argv);
// line below stops the program if the help was requested or
// a compulsory option was not set
- if ( (help.value()) || (!options.check_compulsory_arguments(true)) )
- {
+ if ( (help.value()) || (!options.check_compulsory_arguments(true)) ){
+ options.usage();
+ exit(EXIT_FAILURE);
+ }else{
+ // Call the local functions
+ return do_work(argc,argv);
+ }
+ }catch(X_OptionError& e) {
options.usage();
- exit(EXIT_FAILURE);
- }
-
- } catch(X_OptionError& e) {
- options.usage();
- cerr << endl << e.what() << endl;
+ cerr << endl << e.what() << endl;
exit(EXIT_FAILURE);
- } catch(std::exception &e) {
+ }catch(std::exception &e) {
cerr << e.what() << endl;
}
-
- // Call the local functions
- return do_work(argc,argv);
}
--
GitLab