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;
	  } 
}