Skip to content
Snippets Groups Projects
Commit 66b5321f authored by Christian Beckmann's avatar Christian Beckmann
Browse files

fsl_sbca 1.0

parent 7ae0edbd
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fsl_sbca.cc
+ 328
− 66
View file @ 66b5321f
......@@ -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
......
melhlprfns.cc
+ 11
− 0
View file @ 66b5321f
......@@ -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)
{
......
melhlprfns.h
+ 1
− 0
View file @ 66b5321f
......@@ -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);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment