diff --git a/fsl_sbca.cc b/fsl_sbca.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8ff23769603a489512c27bc02e46d5317763d979
--- /dev/null
+++ b/fsl_sbca.cc
@@ -0,0 +1,258 @@
+/* fsl_glm -
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 2008 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "libvis/miscplot.h"
+#include "miscmaths/miscmaths.h"
+#include "miscmaths/miscprob.h"
+#include "utils/options.h"
+#include <vector>
+#include <time.h>
+#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_sbca (Version 1.0)")+
+ string("\nCopyright(c) 2008, University of Oxford (Christian F. Beckmann)\n")+
+ string(" \n Performs seed-based correlation analysis on FMRI data\n")+
+ string(" using either a single seed coordinate or a seed mask ");
+ string examples="fsl_sbca -i <input> -o <basename> [options]";
+
+//Command line Options {
+ Option<string> fnin(string("-i,--in"), string(""),
+ string(" input file name (4D image file)"),
+ true, requires_argument);
+ Option<string> fnout(string("-o,--out"), string(""),
+ string("output file base name"),
+ true, requires_argument);
+ Option<string> fnseed(string("-s,--seed"), string(""),
+ string("seed voxel coordinate or file name of seed mask (3D file)"),
+ false, requires_argument);
+ Option<string> fntarget(string("-t,--target"), string(""),
+ string("file name of target mask(s) (3D or 4D file)"),
+ false, requires_argument);
+ Option<string> fnconf(string("--conf"), string(""),
+ string(" file name (or comma-separated list of file name) for confound ascii txt files"),
+ false, requires_argument);
+ Option<bool> map_bin(string("--bin"), false,
+ string(" binarise spatial maps prior to calculation of time courses"),
+ false, no_argument);
+ Option<bool> verbose(string("-v,--verbose"), false,
+ string("switch on diagnostic messages"),
+ false, no_argument);
+ Option<bool> tc_mean(string("--mean"), false,
+ string(" use mean instead of Eigenvariates for calculation of time courses"),
+ false, no_argument);
+ Option<int> tc_order(string("--mean"), 1,
+ string(" number of Eigenvariates (default 1)"),
+ false, requires_argument);
+ Option<int> help(string("-h,--help"), 0,
+ string("display this help text"),
+ false,no_argument);
+ /*
+}
+*/
+//Globals {
+ Matrix data, confounds;
+ volume4D<float> orig_data;
+ volume<float> mask;
+ volumeinfo volinf;
+ int voxels = 0;
+ Matrix seeds, corrs;
+ vector<Matrix> ttcs;
+ /*
+}
+*/
+////////////////////////////////////////////////////////////////////////////
+
+// 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);
+}
+
+ReturnMatrix create_confs(string what){
+ Matrix res;
+ char *p;
+ char t[1024];
+ const char *discard = ",";
+
+ strcpy(t, what.c_str());
+ p=strtok(t,discard);
+ res = remmean(read_ascii_matrix(string(p)),1);
+ do{
+ p=strtok(NULL,discard);
+ if(p){
+ res |= remmean(read_ascii_matrix(string(p)),1);
+ }
+ }while(p);
+
+ res.Release();
+ return res;
+}
+
+ReturnMatrix calc_ttc(volume<float>& in){
+ Matrix res, tmp, scales;
+
+ volume<float> tmp1;
+ volume4D<float> tmp2;
+
+ tmp1 = in;
+ tmp1.binarise(1e-8);
+ tmp2.addvolume(in);
+ scales = tmp2.matrix(tmp1);
+ tmp = remmean(orig_data.matrix(tmp1),1);
+
+ if(!map_bin.value())
+ tmp = SP(tmp, ones(tmp.Nrows(),1) * scales);
+
+ if(tc_mean.value())
+ res = mean(tmp,2);
+ else{
+ SymmetricMatrix Corr;
+ Corr << tmp * tmp.t() / tmp.Ncols();
+ DiagonalMatrix tmpD;
+ EigenValues(Corr,tmpD,res);
+ res = fliplr(res.Columns(res.Ncols()-tc_order.value()+1 , res.Ncols())) * std::sqrt(tmp.Nrows());
+ }
+ res.Release();
+ return res;
+}
+
+void create_target_tcs(){
+
+ volume4D<float> tmptarg;
+ read_volume4D(tmptarg,fntarget.value());
+
+ for(int ctr=0; ctr < tmptarg.tsize(); ctr++){
+ ttcs.push_back(calc_ttc(tmptarg[ctr]));
+ }
+}
+
+int setup(){
+ if(fsl_imageexists(fnin.value())){//read data
+ //input is 3D/4D vol
+ read_volume4D(orig_data,fnin.value(),volinf);
+ }
+
+ if(fnconf.value()>"")
+ confounds = create_confs(fnconf.value());
+
+
+ if(fnseed.value()>""){
+ read_volume(mask,fnseed.value());
+ if(!samesize(orig_data[0],mask)){
+ cerr << "ERROR: Seed mask image does not match input image" << endl;
+ return 1;
+ };
+
+ volume4D<float> tmp_mask;
+ tmp_mask.addvolume(mask);
+ mask.binarise(1e-8);
+
+ Matrix scales = tmp_mask.matrix(mask);
+ seeds = remmean(orig_data.matrix(mask),1);
+ if(!map_bin.value())
+ seeds = SP(seeds, ones(seeds.Nrows(),1) * scales);
+ }
+
+ create_target_tcs();
+ voxels = seeds.Ncols();
+ return 0;
+}
+
+void calc_res(){
+
+}
+
+void write_res(){
+
+}
+
+int do_work(int argc, char* argv[]) {
+
+ double tmptime = time(NULL);
+ srand((unsigned int) tmptime);
+
+cerr << (unsigned int) tmptime << endl << endl;
+
+cerr << unifrnd(2,2) << endl;
+exit(1);
+ if(setup())
+ exit(1);
+
+ calc_res();
+
+ 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(fnseed);
+ options.add(fntarget);
+ options.add(fnconf);
+ options.add(map_bin);
+ options.add(tc_mean);
+ options.add(tc_order);
+ options.add(verbose);
+ options.add(help);
+ 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;
+ }
+}
+