diff --git a/fsl_sbca.cc b/fsl_sbca.cc
index 8ff23769603a489512c27bc02e46d5317763d979..e0cd66ebeff0fe7978e67cacc0060aeeb3eb5b5a 100644
--- a/fsl_sbca.cc
+++ b/fsl_sbca.cc
@@ -11,7 +11,6 @@
#include "miscmaths/miscprob.h"
#include "utils/options.h"
#include <vector>
-#include <time.h>
#include "newimage/newimageall.h"
#include "melhlprfns.h"
@@ -26,8 +25,9 @@ using namespace std;
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]";
+ string(" using either a single seed coordinate or a seed mask \n")+
+ string(" ");
+ string examples="fsl_sbca -i <input> -s <seed> -t <target> -o <basename> [options]";
//Command line Options {
Option<string> fnin(string("-i,--in"), string(""),
@@ -38,10 +38,13 @@ using namespace std;
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);
+ true, requires_argument);
Option<string> fntarget(string("-t,--target"), string(""),
string("file name of target mask(s) (3D or 4D file)"),
- false, requires_argument);
+ true, requires_argument);
+ Option<bool> regress_only(string("-r,--reg"), false,
+ string("perform time series regression rather than classification to targets"),
+ false, no_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);
@@ -54,56 +57,162 @@ using namespace std;
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,
+ Option<int> tc_order(string("--order"), 1,
string(" number of Eigenvariates (default 1)"),
false, requires_argument);
+ Option<bool> out_seeds(string("--out_seeds"), false,
+ string("output seed mask image as <basename>_seeds"),
+ false, no_argument);
+ Option<bool> out_seedmask(string("--out_seedmask"), false,
+ string("output seed mask image as <basename>_seedmask"),
+ false, no_argument);
+ Option<bool> out_ttcs(string("--out_ttcs"), false,
+ string("output target time courses as <basename>_ttc<X>.txt"),
+ false, no_argument);
+ Option<bool> out_conf(string("--out_conf"), false,
+ string("output confound time courses as <basename>_confounds.txt"),
+ false, no_argument);
+ Option<bool> out_tcorr(string("--out_tcorr"), false,
+ string("output target correlations as <basename>_tcorr.txt"),
+ false, no_argument, false);
Option<int> help(string("-h,--help"), 0,
string("display this help text"),
false,no_argument);
+ Option<bool> debug(string("--debug"), false,
+ string(" switch on debug messages"),
+ false, no_argument, false);
/*
}
*/
//Globals {
Matrix data, confounds;
volume4D<float> orig_data;
- volume<float> mask;
+ volume<float> maskS, maskT;
volumeinfo volinf;
int voxels = 0;
- Matrix seeds, corrs;
+ Matrix seeds, coords;
vector<Matrix> ttcs;
+
+ Matrix out1, out2;
/*
}
*/
////////////////////////////////////////////////////////////////////////////
// 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 save4D(Matrix what, volume<float>& msk, string fname){
+ if(debug.value())
+ cerr << "DBG: in save4D" << endl;
+ volume4D<float> tempVol;
+ tempVol.setmatrix(what,msk);
+ save_volume4D(tempVol,fname,volinf);
}
-bool isimage(Matrix what){
- if((voxels > 0)&&(what.Ncols()==voxels || what.Nrows()==voxels))
- return TRUE;
- else
- return FALSE;
+void save4D(volume<float>& in, string fname){
+ if(debug.value())
+ cerr << "DBG: in save4D" << endl;
+ volume4D<float> tempVol;
+ tempVol.addvolume(in);
+ save_volume4D(tempVol,fname,volinf);
+}
+
+ReturnMatrix create_coords(string what){
+ if(debug.value())
+ cerr << "DBG: in create_coords" << endl;
+ Matrix res;
+ ifstream fs(what.c_str());
+ if (!fs) {
+ Matrix tmp(1,3);
+ char *p;
+ char t[1024];
+ const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
+ strcpy(t, what.c_str());
+ p=strtok(t,discard);
+ tmp(1,1) = atoi(p);
+ p=strtok(NULL,discard);
+ tmp(1,2) = atoi(p);
+ p=strtok(NULL,discard);
+ tmp(1,3) = atoi(p);
+ res = tmp;
+
+ do{
+ p=strtok(NULL,discard);
+ if(p){
+ tmp(1,1) = atoi(p);
+ p=strtok(NULL,discard);
+ tmp(1,2) = atoi(p);
+ p=strtok(NULL,discard);
+ tmp(1,3) = atoi(p);
+ res &= tmp;
+ }
+ }while(p);
+ }else{
+ res = read_ascii_matrix(fs);
+ fs.close();
+ }
+
+ if(res.Ncols()!=3){
+ cerr << "ERROR: incorrect format " << what << endl;
+ exit(1);
+ }
+
+ // if(verbose.value())
+// cout << " Created seed coordinates (size: " << res.Nrows() << " x " << res.Ncols() << ")" << endl;
+ res.Release();
+
+ return res;
}
-void saveit(Matrix what, string fname){
- if(isimage(what))
- save4D(what,fname);
- else
- write_ascii_matrix(what,fname);
+void create_mask(string what){
+ if(debug.value())
+ cerr << "DBG: in create_mask" << endl;
+
+ coords = create_coords(what);
+ maskS = orig_data[0] * 0.0;
+ for(int ctr = 1; ctr <= coords.Nrows(); ctr++)
+ maskS(coords(ctr,1),coords(ctr,2),coords(ctr,3)) = 1.0;
+ maskS.binarise(1e-8);
+}
+
+void create_seeds(string what){
+ if(debug.value())
+ cerr << "DBG: in create_seeds" << endl;
+
+ if(fsl_imageexists(what)){
+ read_volume(maskS,what);
+ if(!samesize(orig_data[0],maskS)){
+ cerr << "ERROR: Seed mask image does not match input image" << endl;
+ exit(1);
+ }
+ }
+ else
+ create_mask(what);
+
+ volume4D<float> tmp_mask;
+ tmp_mask.addvolume(maskS);
+ maskS.binarise(1e-8);
+
+ Matrix scales = tmp_mask.matrix(maskS);
+ seeds = remmean(orig_data.matrix(maskS),1);
+ if(!map_bin.value())
+ seeds = SP(seeds, ones(seeds.Nrows(),1) * scales);
+
+ voxels = seeds.Ncols();
+ if(debug.value()){
+ cerr << "DBG: " << voxels << " voxels" << endl;
+ cerr << "DBG: seeds matrix is " << seeds.Nrows() << " x " << seeds.Ncols() << endl;
+ }
+
+ if(verbose.value())
+ cout << " Created seed time courses " << endl;
+
}
ReturnMatrix create_confs(string what){
- Matrix res;
+ if(debug.value())
+ cerr << "DBG: in create_confs" << endl;
+
+ Matrix res, tmp;
char *p;
char t[1024];
const char *discard = ",";
@@ -114,15 +223,25 @@ ReturnMatrix create_confs(string what){
do{
p=strtok(NULL,discard);
if(p){
- res |= remmean(read_ascii_matrix(string(p)),1);
+ tmp = read_ascii_matrix(string(p));
+ if(tmp.Nrows()!=res.Nrows()){
+ cerr << "ERROR: confound matrix" << string(p) << " is of wrong size "<< endl;
+ exit(1);
+ }
+ res |= remmean(tmp,1);
}
}while(p);
-
+
+ if(verbose.value())
+ cout << " Created confound matrix (size: " << res.Nrows() << " x " << res.Ncols() << ")" << endl;
res.Release();
return res;
}
ReturnMatrix calc_ttc(volume<float>& in){
+ if(debug.value())
+ cerr << "DBG: in calc_ttc" << endl;
+
Matrix res, tmp, scales;
volume<float> tmp1;
@@ -130,6 +249,8 @@ ReturnMatrix calc_ttc(volume<float>& in){
tmp1 = in;
tmp1.binarise(1e-8);
+ maskT += tmp1;
+
tmp2.addvolume(in);
scales = tmp2.matrix(tmp1);
tmp = remmean(orig_data.matrix(tmp1),1);
@@ -146,74 +267,208 @@ ReturnMatrix calc_ttc(volume<float>& in){
EigenValues(Corr,tmpD,res);
res = fliplr(res.Columns(res.Ncols()-tc_order.value()+1 , res.Ncols())) * std::sqrt(tmp.Nrows());
}
+
+ if(debug.value())
+ cerr << "DBG: size is " << res.Nrows() << " x " << res.Ncols() << endl;
res.Release();
return res;
}
void create_target_tcs(){
-
+ if(debug.value())
+ cerr << "DBG: in create_target_tcs" << endl;
+
volume4D<float> tmptarg;
read_volume4D(tmptarg,fntarget.value());
+ maskT = orig_data[0] * 0.0;
for(int ctr=0; ctr < tmptarg.tsize(); ctr++){
ttcs.push_back(calc_ttc(tmptarg[ctr]));
}
+
+ if(debug.value()) {
+ cerr << "DBG: " << ttcs.size() << " target matrices created " << endl;
+ }
+ if(verbose.value())
+ cout << " Created target mask time courses " << endl;
+
}
int setup(){
- if(fsl_imageexists(fnin.value())){//read data
- //input is 3D/4D vol
- read_volume4D(orig_data,fnin.value(),volinf);
+ if(debug.value())
+ cerr << "DBG: in setup" << endl;
+ if(fsl_imageexists(fnin.value())){ //read data
+ if(verbose.value())
+ cout << " Reading input file " << fnin.value() << endl;
+ read_volume4D(orig_data,fnin.value(),volinf);
+ }
+ else{
+ cerr << "ERROR: Invalid input file " << fnin.value() << endl;
+ exit(1);
}
- if(fnconf.value()>"")
- confounds = create_confs(fnconf.value());
+ create_seeds(fnseed.value());
+
+ if(!regress_only.value())
+ create_target_tcs();
+ else{
+ volume4D<float> tmptarg;
+ read_volume4D(tmptarg,fntarget.value());
+ maskT = tmptarg[0];
+ maskT.binarise(1e-8);
+ data = orig_data.matrix(maskT);
+ data = remmean(data,1);
+ }
+ 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;
- };
+ return 0;
+}
- volume4D<float> tmp_mask;
- tmp_mask.addvolume(mask);
- mask.binarise(1e-8);
+ReturnMatrix calc_tcorr(int in){
+ if(debug.value())
+ cerr << "DBG: in calc_tcorr" << endl;
- 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);
- }
+ Matrix res = zeros(1,seeds.Ncols()), partial_conf, targetcol;
- create_target_tcs();
- voxels = seeds.Ncols();
- return 0;
+ for(int ctr = 0; ctr < (int)ttcs.size(); ctr++)
+ if(ctr != in){
+ if(partial_conf.Storage() == 0)
+ partial_conf = ttcs.at(ctr);
+ else
+ partial_conf |= ttcs.at(ctr);
+ }
+
+ if(ttcs.at(in).Ncols()>1)
+ if(partial_conf.Storage()>0)
+ partial_conf = ttcs.at(in).Columns(2,ttcs.at(in).Ncols()) | partial_conf;
+ else
+ partial_conf = ttcs.at(in).Columns(2,ttcs.at(in).Ncols());
+
+ if(confounds.Storage() > 0)
+ if(partial_conf.Storage()>0)
+ partial_conf |= confounds;
+ else
+ partial_conf = confounds;
+
+ if(debug.value() && partial_conf.Storage()>0)
+ cerr << "DBG: partial_conf " << partial_conf.Nrows() << " x " << partial_conf.Ncols() << endl;
+
+ targetcol = ttcs.at(in).Column(1);
+ if(debug.value())
+ cerr << "DBG: targetcol " << targetcol.Nrows() << " x " << targetcol.Ncols() << endl;
+
+ for(int ctr = 1; ctr <= seeds.Ncols(); ctr++)
+ res(1,ctr) = Melodic::corrcoef(targetcol, seeds.Column(ctr), partial_conf).AsScalar();
+
+ res.Release();
+ return res;
}
void calc_res(){
-
+ if(debug.value())
+ cerr << "DBG: in calc_res" << endl;
+
+ out2 = zeros(1,seeds.Ncols());
+
+ if(!regress_only.value()){
+ //Target TCs exist
+ if(verbose.value())
+ cout << " Calculating partial correlation scores between seeds and targets " << endl;
+
+ Matrix tmp;
+ int tmp2;
+ out1=zeros(ttcs.size(),seeds.Ncols());
+ for(int ctr = 0 ;ctr < (int)ttcs.size(); ctr++)
+ out1.Row(ctr+1) = calc_tcorr(ctr);
+
+ for(int ctr = 1 ;ctr <= out1.Ncols(); ctr++){
+ out1.Column(ctr).Maximum1(tmp2);
+ out2(1,ctr) = tmp2;
+ }
+
+ if(debug.value()){
+ cerr << "DBG: out1 " << out1.Nrows() << " x " << out1.Ncols() << endl;
+ cerr << "DBG: out2 " << out2.Nrows() << " x " << out2.Ncols() << endl;
+ }
+ }
+ else{
+ //no Target TCs
+ if(verbose.value())
+ cout << " Calculating partial correlation maps " << endl;
+
+ out1 = zeros(seeds.Ncols(), data.Ncols());
+ if(confounds.Storage()>0){
+ data = data - confounds * pinv(confounds) * data;
+ seeds = seeds - confounds * pinv(confounds) * seeds;
+ }
+
+ if(debug.value()){
+ cerr << "DBG: seeds " << seeds.Nrows() << " x " << seeds.Ncols() << endl;
+ cerr << "DBG: data " << data.Nrows() << " x " << data.Ncols() << endl;
+ }
+
+ for(int ctr = 1 ;ctr <= seeds.Ncols(); ctr++){
+ Matrix tmp;
+ if(coords.Storage()>0){
+ tmp = orig_data.voxelts(coords(ctr,1), coords(ctr,2), coords(ctr,3));
+ volume4D<float> tmpVol;
+ tmpVol.setmatrix(out2,maskS);
+ tmpVol( coords(ctr,1), coords(ctr,2), coords(ctr,3), 0) = ctr;
+ out2 = tmpVol.matrix(maskS);
+ if(confounds.Storage()>0)
+ tmp = tmp - confounds * pinv(confounds) * tmp;
+ }
+ else{
+ tmp = seeds.Column(ctr);
+ out2(1,ctr) = ctr;
+ }
+ for(int ctr2 =1; ctr2 <= data.Ncols(); ctr2++)
+ out1(ctr,ctr2) = Melodic::corrcoef(tmp,data.Column(ctr2)).AsScalar();
+ }
+
+ if(debug.value()){
+ cerr << "DBG: out1 " << out1.Nrows() << " x " << out1.Ncols() << endl;
+ cerr << "DBG: out2 " << out2.Nrows() << " x " << out2.Ncols() << endl;
+ }
+ }
}
-void write_res(){
+void write_res(){
+ if(verbose.value())
+ cout << " Saving results " << endl;
+
+ if(debug.value())
+ cerr << "DBG: in write_res" << endl;
+
+ if(regress_only.value()){
+ save4D(out2,maskS, fnout.value()+"_index");
+ save4D(out1,maskT, fnout.value()+"_corr");
+ }
+ else{
+ save4D(out1,maskS, fnout.value()+"_corr");
+ save4D(out2,maskS, fnout.value()+"_index");
+ }
+
+ if(out_ttcs.value() && ttcs.size()>0)
+ for(int ctr = 0 ;ctr < (int)ttcs.size(); ctr++)
+ write_ascii_matrix(ttcs.at(ctr),fnout.value()+"_ttc"+num2str(ctr+1)+".txt");
+
+ if(out_conf.value() && confounds.Storage()>0)
+ write_ascii_matrix(confounds, fnout.value()+"_confounds.tx");
+ if(out_seeds.value())
+ save4D(seeds, maskS, fnout.value()+"_seeds");
+ if(out_seedmask.value())
+ save4D(maskS,fnout.value()+"_seedmask");
+
}
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;
}
@@ -227,15 +482,22 @@ int main(int argc,char *argv[]){
// 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(fnout);
options.add(fntarget);
+ options.add(regress_only);
options.add(fnconf);
options.add(map_bin);
options.add(tc_mean);
options.add(tc_order);
+ options.add(out_seeds);
+ options.add(out_seedmask);
+ options.add(out_ttcs);
+ options.add(out_conf);
+ options.add(out_tcorr);
options.add(verbose);
options.add(help);
+ options.add(debug);
options.parse_command_line(argc, argv);
// line below stops the program if the help was requested or
diff --git a/melhlprfns.cc b/melhlprfns.cc
index d160abcd32f0486eac96a8560f86e3bd48003f1d..820b857373a7dda16f8a9cd1eaa310527443f54f 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -175,7 +175,18 @@ namespace Melodic{
Matrix out;
out=MISCMATHS::corrcoef(tmp,0);
return out.SubMatrix(1,in1.Ncols(),in1.Ncols()+1,out.Ncols());
+ }
+
+ Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part){
+ Matrix tmp1 = in1, tmp2 = in2, out;
+ if(part.Storage()){
+ tmp1 = tmp1 - part * pinv(part) * tmp1;
+ tmp2 = tmp2 - part * pinv(part) * tmp2;
}
+
+ out = corrcoef(tmp1,tmp2);
+ return out;
+ }
Matrix calc_corr(const Matrix& in, bool econ)
{
diff --git a/melhlprfns.h b/melhlprfns.h
index 791b4d27c1f1a1a8a39a09bd07c56b61f8b7eab3..7b4dcbef5d24578b9708aa453071ac1593f14873 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -38,6 +38,7 @@ namespace Melodic{
RowVector cumsum(const RowVector& Inp);
Matrix corrcoef(const Matrix& in1, const Matrix& in2);
+ Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part);
Matrix calc_corr(const Matrix& in, bool econ = 0);
Matrix calc_corr(const Matrix& in, const Matrix& weights, bool econ = 0);