From 15e6a9b44c47a2792b2bc139e63d63659f7a249e Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Tue, 13 Sep 2011 15:00:18 +0000
Subject: [PATCH] added steve\'s ad-hoc denoising to fsl_regfilt
---
fsl_regfilt.cc | 273 +++++++++++++++++++++++++++++++++++++++----------
meldata.cc | 67 +++++++-----
melica.cc | 2 +-
meloptions.h | 8 +-
4 files changed, 264 insertions(+), 86 deletions(-)
diff --git a/fsl_regfilt.cc b/fsl_regfilt.cc
index 063b124..3e00500 100644
--- a/fsl_regfilt.cc
+++ b/fsl_regfilt.cc
@@ -2,8 +2,8 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 2006-2008 University of Oxford */
-
+ Copyright (C) 2006-2011 University of Oxford / Christian F. Beckmann */
+
/* Part of FSL - FMRIB's Software Library
http://www.fmrib.ox.ac.uk/fsl
fsl@fmrib.ox.ac.uk
@@ -82,11 +82,11 @@ 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_regfilt (Version 1.0)")+
- string("\n\n Copyright(c) 2008, University of Oxford (Christian F. Beckmann)\n")+
+ string title=string("fsl_regfilt (Version 1.2)")+
+ string("\n\n Copyright(c) 2011, University of Oxford (Christian F. Beckmann)\n")+
string(" Data de-noising by regressing out part of a design matrix\n")+
string(" using simple OLS regression on 4D images");
- string examples="fsl_regfilt -i <input> -d <design> -f <components> -o <out> [options]";
+ string examples="fsl_regfilt -i <input> -d <design> -f <component numbers or filter threshold> -o <out> [options]";
//Command line Options {
Option<string> fnin(string("-i,--in"), string(""),
@@ -100,37 +100,67 @@ using namespace std;
true, requires_argument);
Option<string> fnmask(string("-m,--mask"), string(""),
string("mask image file name"),
+ false, requires_argument);
+ Option<string> filter(string("-f,--filter"),string(""),
+ string("filter out part of the regression model, e.g. -f \"1,2,3\" "),
false, requires_argument);
- Option<string> filter(string("-f,--filter"),string(""),
- string("filter out part of the regression model, e.g. -f \"1,2,3\""),
- true, requires_argument);
+ Option<bool> freqfilt(string("-F,--freqfilt"),false,
+ string("filter out components based on high vs. low frequency content "),
+ false, no_argument);
+ Option<bool> freq_ic(string("--freq_ic"),true,
+ string("switch off IC Z-stats filtering as part of frequency filtering"),
+ false, no_argument);
+ Option<float> freq_ic_smooth(string("--freq_ic_smooth"),5.0,
+ string("smoothing width for IC Z-stats filtering as part of frequency filtering"),
+ false, no_argument);
+ Option<float> freqthresh(string("--fthresh"),0.15,
+ string("frequency threshold ratio - default: 0.15"),
+ false,requires_argument);
+ Option<float> freqthresh2(string("--fthresh2"),0.02,
+ string("frequency filter score threshold - default: 0.02"),
+ false,requires_argument);
Option<bool> verbose(string("-v"),FALSE,
string(" switch on diagnostic messages"),
false, no_argument);
+ Option<bool> aggressive(string("-a"),FALSE,
+ string(" switch on aggressive filtering (full instead of partial regression)"),
+ false, no_argument);
Option<bool> perfvn(string("--vn"),FALSE,
string(" perfrom variance-normalisation on data"),
false, no_argument);
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);
// Output options
Option<string> outdata(string("--out_data"),string(""),
- string("output data"),
+ string("output file name for pre-processed data (prior to denoising)"),
+ false, requires_argument);
+ Option<string> outmix(string("--out_mix"),string(""),
+ string("output file name for new mixing matrix"),
false, requires_argument);
Option<string> outvnscales(string("--out_vnscales"),string(""),
- string("output scaling factors for variance normalisation"),
+ string("output file name for scaling factors from variance normalisation"),
false, requires_argument);
/*
}
*/
//Globals {
int voxels = 0;
+ float TR;
Matrix data;
Matrix design;
- Matrix meanR;
+ Matrix fdesign;
+ Matrix meanR, meanC;
+ Matrix newData, newMix;
RowVector vnscales;
volume<float> mask;
- volume<float> Mean; /*
+ volume<float> Mean;
+ vector<int> comps, ind;
+ vector<int>::iterator it;
+ /*
}
*/
////////////////////////////////////////////////////////////////////////////
@@ -161,50 +191,164 @@ void saveit(Matrix what, string fname){
write_ascii_matrix(what,fname);
}
-int dofilter(){
- if(!isimage(data)){
- cerr << "ERROR: need to specify 4D input to use filtering" << endl;
+Matrix smooth_map(Matrix what, float howmuch){
+ volume4D<float> tempVol;
+ tempVol.setmatrix(what,mask);
+ tempVol= smooth(tempVol,howmuch);
+ Matrix out;
+ out = tempVol.matrix(mask);
+ return out;
+}
+
+int parse_filterstring(){
+ int ctr=0;
+ char *p;
+ char t[1024];
+ const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
+
+ strcpy(t, filter.value().c_str());
+ p=strtok(t,discard);
+ ctr = atoi(p);
+ if(ctr>0 && ctr<=design.Ncols())
+ comps.push_back(ctr);
+ do{
+ p=strtok(NULL,discard);
+ if(p){
+ ctr = atoi(p);
+ if(ctr>0 && ctr<=design.Ncols())
+ comps.push_back(ctr);
+ }
+ }while(p);
+ return 0;
+}
+
+int calc_freqindex(){
+ if(debug.value()) cerr << " In calc_freqindex " << endl;
+
+ fdesign = Melodic::calc_FFT(design);
+ if(debug.value()) cerr << " fdesign: " << fdesign.Nrows() << " x " << fdesign.Ncols() << endl;
+
+ int Nps = fdesign.Nrows();
+ float MAXf = 1/(2*TR);
+ float Nthresh = ceil(Nps * freqthresh.value()/MAXf);
+ if(debug.value()) cerr << " Nps: " << Nps << " MAXf: " << MAXf << " Nthresh: " << Nthresh << endl;
+
+ Matrix sum_ratio;
+ sum_ratio = SP(sum(fdesign.Rows(1,Nthresh),1),pow(sum(sum(fdesign.Rows(Nthresh+1,Nps))),-1));
+ sum_ratio /= (float)sum_ratio.MaximumAbsoluteValue();
+ if(debug.value()) cerr << " sum_ratio: " << sum_ratio << endl;
+
+ if(freq_ic.value()){
+ Matrix scores = zeros(1,design.Ncols());
+ {
+ Matrix ICs, noisestddev, stdNoisei,unmixMatrix;
+
+ unmixMatrix = pinv(design);
+ ICs = unmixMatrix * data;
+ noisestddev = stdev(data-design*ICs);
+ stdNoisei = pow(noisestddev*
+ std::sqrt((float)(data.Nrows()-1))/
+ std::sqrt((float)(data.Nrows()-ICs.Nrows())),-1);
+
+ ColumnVector diagvals;
+ diagvals = pow(diag( unmixMatrix*unmixMatrix.t()),-0.5);
+
+ ICs=smooth_map(SP(ICs,diagvals*stdNoisei),freq_ic_smooth.value());
+ ICs= SP(ICs,ones(ICs.Nrows(),1)*meanR);
+ volume4D<float> tempVol;
+ tempVol.setmatrix(ICs,mask);
+ tempVol.threshold(0.0);
+ for(int ctr = 0; ctr < design.Ncols(); ctr++ )
+ scores(1,ctr+1) = tempVol[ctr].percentile(0.99,mask);
+ scores-=scores.MinimumAbsoluteValue();
+ scores/=scores.MaximumAbsoluteValue();
+ if(debug.value()) cerr << " initial scores: " << scores << endl;
+ }
+ scores = SP(scores,sum_ratio);
+ scores /= scores.Maximum();
+ if(debug.value()) cerr << " scores: " << scores << endl;
+ for(int ctr = 1; ctr <= design.Ncols(); ctr++ )
+ if(scores(1,ctr) < freqthresh2.value())
+ comps.push_back(ctr);
+ }
+ return 0;
+}
+
+int get_comp(){
+ if(filter.value().length()>0 && parse_filterstring())
+ return 1;
+
+ if(freqfilt.value() && calc_freqindex())
return 1;
+
+ //sort and remove duplicates
+ sort (comps.begin(), comps.end());
+ it = unique (comps.begin(), comps.end());
+ comps.resize( it - comps.begin() );
+
+ if(debug.value()){
+ for (it=comps.begin(); it!=comps.end(); ++it)
+ cout << " " << *it;
+ cout << endl;
}
+
+ return 0;
+}
+
+int dofilter(){
+
if(verbose.value())
cout << " Calculating maps " << endl;
- Matrix unmixMatrix = pinv(design);
- Matrix maps = unmixMatrix * data;
-
- Matrix noisedes;
- Matrix noisemaps;
-
- int ctr=0;
- char *p;
- char t[1024];
- const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
-
- strcpy(t, filter.value().c_str());
- p=strtok(t,discard);
- ctr = atoi(p);
- if(ctr>0 && ctr<=design.Ncols()){
- noisedes = design.Column(ctr);
- noisemaps = maps.Row(ctr).t();
- }
-
- do{
- p=strtok(NULL,discard);
- if(p){
- ctr = atoi(p);
- if(ctr>0 && ctr<=design.Ncols()){
- noisedes |= design.Column(ctr);
- noisemaps |= maps.Row(ctr).t();
- }
- }
- }while(p);
- Matrix newData;
+
+ Matrix unmixMatrix = pinv(design);
+ Matrix maps = unmixMatrix * data;
+
+ Matrix noisedes;
+ Matrix noisemaps;
+
+ noisedes = design.Column(comps.at(0));
+ noisemaps = maps.Row(comps.at(0)).t();
+
+ for(int ctr = 1; ctr < (int)comps.size();++ctr){
+ noisedes |= design.Column(comps.at(ctr));
+ noisemaps |= maps.Row(comps.at(ctr)).t();
+ }
+ if(debug.value()) cerr << " noisedes " << noisedes.Nrows() << " x " << noisedes.Ncols() << endl;
+
if(verbose.value())
cout << " Calculating filtered data " << endl;
- newData = data - noisedes * noisemaps.t();
- newData = newData + ones(newData.Nrows(),1)*meanR;
-
- save4D(newData,fnout.value());
- return 0;
+
+ if(aggressive.value())
+ newData = data - noisedes * (pinv(noisedes)*data);
+ else
+ newData = data - noisedes * noisemaps.t();
+
+ if(perfvn.value())
+ newData = SP(newData,ones(newData.Nrows(),1)*vnscales);
+ newData = newData + ones(newData.Nrows(),1)*meanR;
+
+ for(int ctr = 1; ctr <= design.Ncols();++ctr)
+ ind.push_back(ctr);
+ for(int ctr = 0; ctr < (int)comps.size();++ctr)
+ it=remove(ind.begin(),ind.end(),comps.at(ctr));
+ ind.resize(design.Ncols()-comps.size());
+ if(debug.value()){
+ for (it=ind.begin(); it!=ind.end(); ++it)
+ cout << " " << *it;
+ cout << endl;
+ }
+
+ if(ind.size()>0){
+ newMix=design.Column(ind.at(0));
+ for(int ctr = 1; ctr < (int)ind.size();++ctr)
+ newMix |= design.Column(ind.at(ctr));
+ newMix = newMix - noisedes * (pinv(noisedes)*newMix);
+
+ if(debug.value())
+ cerr << " newMix " << newMix.Nrows() << " x " << newMix.Ncols() << endl;
+ }
+
+ return 0;
}
int setup(){
@@ -212,6 +356,7 @@ int setup(){
//input is 3D/4D vol
volume4D<float> tmpdata;
read_volume4D(tmpdata,fnin.value());
+ TR=tmpdata.TR();
// create mask
if(fnmask.value()>""){
@@ -238,28 +383,42 @@ int setup(){
cerr << "ERROR: cannot read input image " << fnin.value()<<endl;
return 1;
}
-
+
design = read_ascii_matrix(fndesign.value());
+ if(!isimage(data)){
+ cerr << "ERROR: need to specify 4D input to use filtering" << endl;
+ return 1;
+ }
+
meanR=mean(data,1);
data = remmean(data,1);
+ meanC=mean(design,1);
design = remmean(design,1);
if(perfvn.value())
vnscales = Melodic::varnorm(data);
+
+ if(debug.value()) cerr << " data: " << data.Nrows() << " x " << data.Ncols() << endl;
+ if(debug.value()) cerr << " design: " << design.Nrows() << " x " << design.Ncols() << endl;
+
return 0;
}
void write_res(){
+ saveit(newData,fnout.value());
if(outdata.value()>"")
saveit(data,outdata.value());
if(outvnscales.value()>"")
saveit(vnscales,outvnscales.value());
+ if(outmix.value()>"" && newMix.Storage()>0)
+ saveit(newMix,outmix.value());
}
int do_work(int argc, char* argv[]) {
- if(setup())
+ if(setup())
+ exit(1);
+ if(get_comp())
exit(1);
-
if(dofilter())
exit(1);
write_res();
@@ -279,10 +438,18 @@ int main(int argc,char *argv[]){
options.add(fndesign);
options.add(fnmask);
options.add(filter);
+ options.add(freqfilt);
+ options.add(freq_ic);
+ options.add(freq_ic_smooth);
+ options.add(freqthresh);
+ options.add(freqthresh2);
options.add(perfvn);
options.add(verbose);
+ options.add(aggressive);
options.add(help);
+ options.add(debug);
options.add(outdata);
+ options.add(outmix);
options.add(outvnscales);
options.parse_command_line(argc, argv);
diff --git a/meldata.cc b/meldata.cc
index 44fd9e6..3298324 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -202,7 +202,7 @@ namespace Melodic{
void MelodicData::set_TSmode()
{
- Matrix tmp, tmpT, tmpS, tmpT2, tmpS2, tmpT3;
+ Matrix tmp, tmpT, tmpS, tmpT2, tmpS2, tmpT3;
tmp = expand_dimred(mixMatrix);
tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
tmpS = zeros(numfiles, tmp.Ncols());
@@ -210,24 +210,25 @@ namespace Melodic{
outMsize("tmp",tmp);
outMsize("tmpT",tmpT);
outMsize("tmpS",tmpS);
+ if(opts.approach.value()==string("tica")){
+ Tmodes.clear(); Smodes.clear();
+ for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
+ tmpT3 << reshape(tmp.Column(ctr),tmpT.Nrows(),numfiles);
+ outMsize("tmpT3", tmpT3);
+ tmpT2 << tmpT.Column(ctr);
+ tmpS2 << tmpS.Column(ctr);
+ tmpT3 << SP(tmpT3,pow(ones(tmpT3.Nrows(),1)*tmpS2.t(),-1));
+ if(numfiles>1)
+ tmpT2 |= tmpT3;
+ if(mean(tmpS2,1).AsScalar()<0){
+ tmpT2*=-1.0;
+ tmpS2*=-1.0;
+ }
+ add_Tmodes(tmpT2);
+ add_Smodes(tmpS2);
+ }
+ }
- Tmodes.clear(); Smodes.clear();
- for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
- tmpT3 << reshape(tmp.Column(ctr),tmpT.Nrows(),numfiles);
- outMsize("tmpT3", tmpT3);
- tmpT2 << tmpT.Column(ctr);
- tmpS2 << tmpS.Column(ctr);
- tmpT3 << SP(tmpT3,pow(ones(tmpT3.Nrows(),1)*tmpS2.t(),-1));
- if(numfiles>1)
- tmpT2 |= tmpT3;
- if(mean(tmpS2,1).AsScalar()<0){
- tmpT2*=-1.0;
- tmpS2*=-1.0;
- }
- add_Tmodes(tmpT2);
- add_Smodes(tmpS2);
- }
-
//add GLM OLS fit
if(Tdes.Storage()){
Matrix alltcs = Tmodes.at(0).Column(1);
@@ -274,14 +275,22 @@ namespace Melodic{
if(opts.debug.value())
save4D(alldat,string("preproc_dat") + num2str(1));
for(int ctr = 1; ctr < numfiles; ctr++){
- tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
+ tmpData = process_file(opts.inputfname.value().at(ctr), numfiles) / numfiles;
if(opts.debug.value())
- save4D(tmpData /numfiles,string("preproc_dat") + num2str(ctr+1));
+ save4D(tmpData,string("preproc_dat") + num2str(ctr+1));
if(tmpData.Ncols() == alldat.Ncols() && tmpData.Nrows() == alldat.Nrows())
- alldat += tmpData / numfiles;
- else
- message("Data dimensions do not match - ignoring "+opts.inputfname.value().at(ctr) << endl);
- }
+ alldat += tmpData;
+ else{
+ if(tmpData.Ncols() == alldat.Ncols()){
+ int mindim = min(alldat.Nrows(),tmpData.Nrows());
+ alldat = alldat.Rows(1,mindim);
+ tmpData = tmpData.Rows(1,mindim);
+ alldat += tmpData;
+ }
+ else
+ message("Data dimensions do not match - ignoring "+opts.inputfname.value().at(ctr) << endl);
+ }
+ }
//update mask
if(opts.update_mask.value()){
@@ -308,16 +317,18 @@ namespace Melodic{
RowVector AdjEV, PercEV;
Matrix Corr, tmpE;
int order;
-
+ cerr << "here1" << endl;
order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
if (opts.paradigmfname.value().length()>0)
order += param.Ncols();
-
+ cerr << "here2" << endl;
if(opts.pca_dim.value() == 0){
opts.pca_dim.set_T(order);
PPCA=tmpPPCA;
}
order = opts.pca_dim.value();
+ if(opts.debug.value())
+ message(endl << "Model order : "<<order<<endl<<endl);
if (opts.paradigmfname.value().length()>0){
Matrix tmpPscales;
@@ -349,7 +360,7 @@ namespace Melodic{
WM.push_back(tmp);
}
else {
-
+ cerr << "here" << endl;
for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
@@ -930,7 +941,7 @@ namespace Melodic{
}
message("Sorting IC maps" << endl);
Matrix tmpscales, tmpICrow, tmpMIXcol;
- if(numfiles > 1){
+ if(numfiles > 1 && opts.approach.value()==string("tica")){
set_TSmode();
Matrix allmodes = Smodes.at(0);
for(int ctr = 1; ctr < (int)Smodes.size();++ctr)
diff --git a/melica.cc b/melica.cc
index e4b945d..6bfa57e 100644
--- a/melica.cc
+++ b/melica.cc
@@ -536,7 +536,7 @@ namespace Melodic{
melodat.set_ICstats(scales);
melodat.sort();
- message("Calculating T- and S-modes " << endl);
+ message("Calculating T- and S-modes " << endl);
melodat.set_TSmode();
}
diff --git a/meloptions.h b/meloptions.h
index 9db9fcc..6f5035f 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -272,8 +272,8 @@ class MelodicOptions {
string("switch off variance normalisation"),
false, no_argument),
pbsc(string("--pbsc"), false,
- string(" switch off conversion to percent BOLD signal change"),
- false, no_argument),
+ string(" switch on conversion to percent BOLD signal change"),
+ false, no_argument, false),
pspec(string("--pspec"), false,
string(" switch on conversion to powerspectra"),
false, no_argument, false),
@@ -413,7 +413,7 @@ class MelodicOptions {
string("number of repeats (multistart)"),
false, requires_argument, false),
seed(string("--seed"), -1,
- string("integer seed for melodic"),
+ string("integer seed for random number generator within melodic"),
false, requires_argument, false),
nlconst1(string("--nl1,--nlconst1"), 1.0,
string("nonlinear constant 1"),
@@ -438,7 +438,7 @@ class MelodicOptions {
false, no_argument,false),
guess_remderiv(string("--remove_deriv"), false,
string("removes every second entry in paradigm file (EV derivatives)"),
- false, no_argument),
+ false, no_argument, false),
temporal(string("--temporal"), false,
string("perform temporal ICA"),
false, no_argument, false),
--
GitLab