Skip to content
Snippets Groups Projects
Commit 56439f0d authored by Christian Beckmann's avatar Christian Beckmann
Browse files

added fsl_glm and fsl_regfilt

parent ee86ea29
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Makefile
+ 12
− 3
View file @ 56439f0d
...@@ -16,17 +16,21 @@ TEST_OBJS = melhlprfns.o test.o ...@@ -16,17 +16,21 @@ TEST_OBJS = melhlprfns.o test.o
GGMIX_OBJS = ggmix.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 MELODIC_OBJS = meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o melreport.o melodic.o
TESTXFILES = test TESTXFILES = test
XFILES = melodic XFILES = fsl_glm fsl_regfilt melodic
RUNTCLS = Melodic RUNTCLS = Melodic
SCRIPTS = melodicreport SCRIPTS = melodicreport
all: ggmix melodic all: ggmix fsl_regfilt fsl_glm melodic
ggmix: ${GGMIX_OBJS} ggmix: ${GGMIX_OBJS}
${AR} -r libggmix.a ${GGMIX_OBJS} ${AR} -r libggmix.a ${GGMIX_OBJS}
...@@ -34,8 +38,13 @@ ggmix: ${GGMIX_OBJS} ...@@ -34,8 +38,13 @@ ggmix: ${GGMIX_OBJS}
melodic: ${MELODIC_OBJS} melodic: ${MELODIC_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${MELODIC_OBJS} ${LIBS} $(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${MELODIC_OBJS} ${LIBS}
test: ${TEST_OBJS} test: ${TEST_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${TEST_OBJS} ${LIBS} $(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}
fsl_glm.cc 0 → 100644
+ 275
− 0
View file @ 56439f0d
/* 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;
}
}
fsl_regfilt.cc 0 → 100644
+ 230
− 0
View file @ 56439f0d
/* 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;
}
}
meldata.cc
+ 22
− 1
View file @ 56439f0d
...@@ -140,9 +140,27 @@ namespace Melodic{ ...@@ -140,9 +140,27 @@ namespace Melodic{
for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){ for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
tmpT2 << tmpT.Column(ctr); tmpT2 << tmpT.Column(ctr);
tmpS2 << tmpS.Column(ctr); tmpS2 << tmpS.Column(ctr);
if(mean(tmpS2,1).AsScalar()<0){
tmpT2*=-1.0;
tmpS2*=-1.0;
}
add_Tmodes(tmpT2); add_Tmodes(tmpT2);
add_Smodes(tmpS2); 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() void MelodicData::setup()
...@@ -668,10 +686,13 @@ namespace Melodic{ ...@@ -668,10 +686,13 @@ namespace Melodic{
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(), int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
numTime = mixMatrix.Nrows(), i,j; 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++) for(int ctr_i = 1; ctr_i <= numComp; ctr_i++)
if(IC.Row(ctr_i).Sum()<0) if(IC.Row(ctr_i).Sum()<0)
flipres(ctr_i); flipres(ctr_i);
// re-order wrt standard deviation of IC maps // re-order wrt standard deviation of IC maps
message("Sorting IC maps" << endl); message("Sorting IC maps" << endl);
Matrix tmpscales, tmpICrow, tmpMIXcol; Matrix tmpscales, tmpICrow, tmpMIXcol;
......
meldata.h
+ 1
− 2
View file @ 56439f0d
...@@ -183,6 +183,7 @@ namespace Melodic{ ...@@ -183,6 +183,7 @@ namespace Melodic{
volumeinfo tempInfo; volumeinfo tempInfo;
vector<Matrix> DWM, WM; vector<Matrix> DWM, WM;
basicGLM glmT, glmS; basicGLM glmT, glmS;
Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;
private: private:
MelodicOptions &opts; MelodicOptions &opts;
...@@ -202,8 +203,6 @@ namespace Melodic{ ...@@ -202,8 +203,6 @@ namespace Melodic{
Matrix mixFFT; Matrix mixFFT;
Matrix IC; Matrix IC;
Matrix ICstats; Matrix ICstats;
Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;
vector<Matrix> Tmodes; vector<Matrix> Tmodes;
vector<Matrix> Smodes; vector<Matrix> Smodes;
......
melhlprfns.cc
+ 15
− 14
View file @ 56439f0d
...@@ -499,7 +499,6 @@ namespace Melodic{ ...@@ -499,7 +499,6 @@ namespace Melodic{
return res; return res;
} //int ppca_dim } //int ppca_dim
int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which) int ppca_dim(const Matrix& in, const Matrix& weights, float resels, string which)
{ {
ColumnVector PPCA; ColumnVector PPCA;
...@@ -783,7 +782,7 @@ namespace Melodic{ ...@@ -783,7 +782,7 @@ namespace Melodic{
return res; return res;
} //Matrix gen_arCorr } //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) const Matrix& contrasts, int DOFadjust)
{ {
beta = zeros(design.Ncols(),1); beta = zeros(design.Ncols(),1);
...@@ -793,28 +792,30 @@ namespace Melodic{ ...@@ -793,28 +792,30 @@ namespace Melodic{
if(data.Nrows()==design.Nrows()){ if(data.Nrows()==design.Nrows()){
Matrix dat = remmean(data,1); Matrix dat = remmean(data,1);
Matrix pinvdes = pinv(design); Matrix tmp = design.t()*design;
Matrix pinvdes = tmp.i()*design.t();
beta = pinvdes * dat; beta = pinvdes * dat;
residu = dat - design*beta; residu = dat - design*beta;
Matrix R = Identity(design.Nrows()) - design * pinvdes; // Matrix R = Identity(design.Nrows()) - design * pinvdes;
Matrix R2 = R*R; // Matrix R2 = R*R;
float tR = R.Trace(); // float tR = R.Trace();
sigsq = sum(SP(residu,residu))/tR; // 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(); float fact = float(design.Nrows() - 1 - design.Ncols()) / design.Ncols();
f_fmf = SP(var(design*beta),pow(var(residu),-1))* fact; 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++) for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(), pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),
int(design.Nrows() -1 -design.Ncols()),f_fmf.Column(ctr1).AsScalar()); 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()){ if(contrasts.Storage()>0 && contrasts.Ncols()==beta.Nrows()){
cbeta = contrasts.t()*beta; cbeta = contrasts*beta;
Matrix tmp = contrasts.t()*pinvdes*pinvdes.t()*contrasts; Matrix tmp = contrasts*pinvdes*pinvdes.t()*contrasts.t();
varcb = diag(tmp)*sigsq; varcb = diag(tmp)*sigsq;
t = SP(cbeta,pow(varcb,-0.5)); t = SP(cbeta,pow(varcb,-0.5));
z = t; p=t; z = t; p=t;
...@@ -827,7 +828,7 @@ namespace Melodic{ ...@@ -827,7 +828,7 @@ namespace Melodic{
} }
} }
} }
return beta;
} }
......
melhlprfns.h
+ 2
− 1
View file @ 56439f0d
...@@ -80,7 +80,7 @@ namespace Melodic{ ...@@ -80,7 +80,7 @@ namespace Melodic{
//destructor //destructor
~basicGLM(){} ~basicGLM(){}
Matrix olsfit(const Matrix& data, const Matrix& design, void olsfit(const Matrix& data, const Matrix& design,
const Matrix& contrasts, int DOFadjust = 0); const Matrix& contrasts, int DOFadjust = 0);
inline Matrix& get_t(){return t;} inline Matrix& get_t(){return t;}
...@@ -89,6 +89,7 @@ namespace Melodic{ ...@@ -89,6 +89,7 @@ namespace Melodic{
inline Matrix& get_f_fmf(){return f_fmf;} inline Matrix& get_f_fmf(){return f_fmf;}
inline Matrix& get_pf_fmf(){return pf_fmf;} inline Matrix& get_pf_fmf(){return pf_fmf;}
inline Matrix& get_cbeta(){return cbeta;} inline Matrix& get_cbeta(){return cbeta;}
inline Matrix& get_beta(){return beta;}
inline Matrix& get_varcb(){return varcb;} inline Matrix& get_varcb(){return varcb;}
inline Matrix& get_sigsq(){return sigsq;} inline Matrix& get_sigsq(){return sigsq;}
inline Matrix& get_residu(){return residu;} inline Matrix& get_residu(){return residu;}
......
melica.cc
+ 1
− 2
View file @ 56439f0d
...@@ -102,7 +102,7 @@ namespace Melodic{ ...@@ -102,7 +102,7 @@ namespace Melodic{
cum_itt += itt_ctr; cum_itt += itt_ctr;
itt_ctr2++; itt_ctr2++;
if(opts.approach.value() == string("tica")){ 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); Matrix temp(melodat.get_dewhite() * redUMM);
outMsize(" 2nd unmmixing matrix ", temp); outMsize(" 2nd unmmixing matrix ", temp);
temp = melodat.expand_dimred(temp); temp = melodat.expand_dimred(temp);
...@@ -459,7 +459,6 @@ namespace Melodic{ ...@@ -459,7 +459,6 @@ namespace Melodic{
melodat.set_IC(temp); melodat.set_IC(temp);
melodat.set_ICstats(scales); melodat.set_ICstats(scales);
melodat.sort(); melodat.sort();
melodat.set_TSmode(); melodat.set_TSmode();
} }
......
meloptions.cc
+ 1
− 1
View file @ 56439f0d
...@@ -157,7 +157,7 @@ MelodicOptions* MelodicOptions::gopt = NULL; ...@@ -157,7 +157,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
} }
//in the case of indirect inputs, create the vector of input names here //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; std::vector< string > tmpfnames;
ifstream fs(inputfname.value().at(0).c_str()); ifstream fs(inputfname.value().at(0).c_str());
string cline; string cline;
......
meloptions.h
+ 1
− 6
View file @ 56439f0d
...@@ -46,7 +46,6 @@ class MelodicOptions { ...@@ -46,7 +46,6 @@ class MelodicOptions {
Option<string> maskfname; Option<string> maskfname;
Option<bool> use_mask; Option<bool> use_mask;
Option<bool> inputindirect;
Option<bool> update_mask; Option<bool> update_mask;
Option<bool> perf_bet; Option<bool> perf_bet;
Option<float> threshold; Option<float> threshold;
...@@ -147,7 +146,7 @@ class MelodicOptions { ...@@ -147,7 +146,7 @@ class MelodicOptions {
string("output directory name\n"), string("output directory name\n"),
false, requires_argument), false, requires_argument),
inputfname(string("-i,--in"), std::vector<string>(), 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), true, requires_argument),
outputfname(string("-O,--out"), string("melodic"), outputfname(string("-O,--out"), string("melodic"),
string("output file name"), string("output file name"),
...@@ -158,9 +157,6 @@ class MelodicOptions { ...@@ -158,9 +157,6 @@ class MelodicOptions {
use_mask(string("--nomask"), true, use_mask(string("--nomask"), true,
string("switch off masking"), string("switch off masking"),
false, no_argument), 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, update_mask(string("--update_mask"), true,
string("switch off mask updating"), string("switch off mask updating"),
false, no_argument), false, no_argument),
...@@ -351,7 +347,6 @@ class MelodicOptions { ...@@ -351,7 +347,6 @@ class MelodicOptions {
options.add(guessfname); options.add(guessfname);
options.add(maskfname); options.add(maskfname);
options.add(use_mask); options.add(use_mask);
options.add(inputindirect);
options.add(update_mask); options.add(update_mask);
options.add(perf_bet); options.add(perf_bet);
options.add(threshold); options.add(threshold);
......
melreport.cc
+ 103
− 18
View file @ 56439f0d
...@@ -108,23 +108,17 @@ namespace Melodic{ ...@@ -108,23 +108,17 @@ namespace Melodic{
} }
{//plot time course {//plot time course
IChtml << "<H3> Temporal mode <p>" << endl <<endl; IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
miscplot newplot; miscplot newplot;
Matrix tmptc = melodat.get_Tmodes(cnum-1).t(); Matrix tmptc = melodat.get_Tmodes(cnum-1).t();
//add GLM OLS fit //add GLM OLS fit
/* basicGLM glm; if(melodat.Tdes.Storage() > 0){
if(melodat.Tdes.Storage()){ tmptc &= melodat.glmT.get_beta().Column(cnum).t() * melodat.Tdes.t();
Matrix betas = glm.olsfit(tmptc.t(),melodat.Tdes,melodat.Tcon).t();
tmptc &= betas*melodat.Tdes.t();
newplot.add_label(string("IC ")+num2str(cnum)+" time course"); newplot.add_label(string("IC ")+num2str(cnum)+" time course");
newplot.add_label("full model fit"); 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) if(opts.tr.value()>0.0)
newplot.timeseries(tmptc, newplot.timeseries(tmptc,
report.appendDir(string("t")+num2str(cnum)+".png"), report.appendDir(string("t")+num2str(cnum)+".png"),
...@@ -140,10 +134,9 @@ glm.get_pf_fmf() << endl << endl; ...@@ -140,10 +134,9 @@ glm.get_pf_fmf() << endl << endl;
IChtml << "<A HREF=\"" << string("t") IChtml << "<A HREF=\"" << string("t")
+num2str(cnum)+".txt" << "\"> "; +num2str(cnum)+".txt" << "\"> ";
IChtml << "<img BORDER=0 SRC=\"" IChtml << "<img BORDER=0 SRC=\""
+string("t")+num2str(cnum)+".png\"></A><p>" << endl; +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
}//time series plot }//time series plot
{//plot frequency
{//plot frequency
miscplot newplot; miscplot newplot;
RowVector empty(1); RowVector empty(1);
empty = 0.0; empty = 0.0;
...@@ -181,17 +174,74 @@ glm.get_pf_fmf() << endl << endl; ...@@ -181,17 +174,74 @@ glm.get_pf_fmf() << endl << endl;
IChtml << "<img BORDER=0 SRC=\"" IChtml << "<img BORDER=0 SRC=\""
+string("f")+num2str(cnum)+".png\"></A><p>" << endl; +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
}//frequency plot }//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 &beta;'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()) if(cnum <= (int)melodat.get_Smodes().size())
{//plot subject mode {//plot subject mode
Matrix smode; Matrix smode;
smode = melodat.get_Smodes(cnum-1); smode = melodat.get_Smodes(cnum-1);
if(smode.Nrows() > 1){ if(smode.Nrows() > 1){
miscplot newplot; 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.setscatter(smode,5);
newplot.timeseries(smode.t(), newplot.timeseries(smode.t(),
report.appendDir(string("s")+num2str(cnum)+".png"), report.appendDir(string("s")+num2str(cnum)+".png"),
string("Subject/Session mode")); string("Subject/Session mode"));
newplot.set_xysize(120,200);
newplot.set_minmaxscale(1.1);
newplot.boxplot(smode, newplot.boxplot(smode,
report.appendDir(string("b")+num2str(cnum)+".png"), report.appendDir(string("b")+num2str(cnum)+".png"),
string("Subject/Session mode")); string("Subject/Session mode"));
...@@ -206,8 +256,43 @@ glm.get_pf_fmf() << endl << endl; ...@@ -206,8 +256,43 @@ glm.get_pf_fmf() << endl << endl;
IChtml << "<img BORDER=0 SRC=\"" IChtml << "<img BORDER=0 SRC=\""
+string("b")+num2str(cnum)+".png\"></A><p>" << endl; +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 << " &beta;(" <<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&& if(mmodel.get_threshmaps().Storage()>0&&
(mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&& (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&&
(mmodel.get_threshmaps().Nrows()>1)) (mmodel.get_threshmaps().Nrows()>1))
......
test.cc
+ 278
− 48
View file @ 56439f0d
...@@ -9,83 +9,313 @@ ...@@ -9,83 +9,313 @@
#include "miscmaths/miscprob.h" #include "miscmaths/miscprob.h"
#include "utils/options.h" #include "utils/options.h"
#include <vector> #include <vector>
#include "newimage/newimageall.h"
#include "melhlprfns.h"
using namespace MISCPLOT; using namespace MISCPLOT;
using namespace MISCMATHS; using namespace MISCMATHS;
using namespace Utilities; using namespace Utilities;
using namespace std; using namespace std;
// The two strings below specify the title and example usage that is // The two strings below specify the title and example usage that is
// printed out as the help or usage message // 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";
Option<bool> help(string("--help"), false, string title=string("fsl_glm (Version 1.0)")+
string(" display this message"), string("\nCopyright(c) 2007, University of Oxford (Christian F. Beckmann)\n")+
false, no_argument); string(" \n Simple GLM usign ordinary least-squares regression on\n")+
Option<int> num(string("--num"), 1, string(" time courses and/or 3D/4D volumes\n\n");
string("number of iterations"), string examples="fsl_glm <input> -d <design> [options]";
false, requires_argument);
int nonoptarg;
//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 // 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[]) { int do_work(int argc, char* argv[]) {
if(setup())
Matrix mat; exit(1);
cout << "BLAH " << num.value() << endl; if(filter.value()>"")
mat=normrnd(300,1); if(dofilter())
miscplot::Timeseries(mat.t(),string("test0.png"),string("TEST")); exit(1);
else{
for (int i=1; i <= (int)num.value();i++){ glm.olsfit(data,design,contrasts,dofset.value());
cout << "Processing " << i << endl; write_res();
miscplot newplot; }
newplot.GDCglobals_set(); return 0;
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;
} }
//////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////
int main(int argc,char *argv[]){ int main(int argc,char *argv[]){
Tracer tr("main"); Tracer tr("main");
OptionParser options(title, examples); OptionParser options(title, examples);
try{
try {
// must include all wanted options here (the order determines how // must include all wanted options here (the order determines how
// the help message is printed) // the help message is printed)
options.add(help); options.add(fnin);
options.add(num); 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); options.parse_command_line(argc, argv);
// line below stops the program if the help was requested or // line below stops the program if the help was requested or
// a compulsory option was not set // 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(); options.usage();
exit(EXIT_FAILURE); cerr << endl << e.what() << endl;
}
} catch(X_OptionError& e) {
options.usage();
cerr << endl << e.what() << endl;
exit(EXIT_FAILURE); exit(EXIT_FAILURE);
} catch(std::exception &e) { }catch(std::exception &e) {
cerr << e.what() << endl; cerr << e.what() << endl;
} }
// Call the local functions
return do_work(argc,argv);
} }
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment