initial checkin

fsl_sbca.cc

+/*  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;
+	  } 
+}
+