diff --git a/Makefile b/Makefile
index 3919431c916f8e8248d3ad2dafeedb7d96667acc..6cc7f4aa8847eb624fa3601a5d3828238be170fc 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 0000000000000000000000000000000000000000..e7ce0d0df26f1ae32dc33c59655b4f86b9f6f121
--- /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 0000000000000000000000000000000000000000..e8f4c1c7ce37ba7bf3308ce8a5ce535e65c46e4e
--- /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 5a2f43c5b384861245241d05c24c8c79b43436d8..89fffd3a959852ad18fb01150032253f828d48f3 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 62314691cf01bcaa2c9fe231af3b568f030314c3..4ebec202466067c58cb8462a830c382093f51240 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 6c881885d69cb2fc0fee31248f6907542cd3f61c..05a73d49b4c90ff4d28d1893efe333eaeaf4e3ea 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 e0334cfd994e7464004fdb4a4de412c5738fb884..6f8d5a41d48be7aa98a72330315b57fc96425e4a 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 253c6200f954b3cc94460930ecf72db0a2082624..ee4fb3f16f2048c82ac3b46d0d43b49692be7f26 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 c534053d436666034d0cbe431a88c4487ebc0480..5742c173ee9a56c679398c5431d5cc0944fd419c 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 d1c11fa8e2b63ae8779906c78e27f520759d01e3..81d7bc5894274f2009f1a03184c93f4db8bf8372 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 3fe7994f5ba8546f3ffaaf977ece6bed257d8045..3597b780ec713e5a59d319fcbb2f446d31b81bbf 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 &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())
 	    {//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 << " &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&&
 	 			(mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&&
 	 			(mmodel.get_threshmaps().Nrows()>1))
diff --git a/test.cc b/test.cc
index 1077016e79348046d051618d191489d555afcff1..9423bad2c8e89b22413e85a777e8405d1c6efa55 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);
 	}