From 90a94fa7144e4ba98386787c5e8a0d91b35fc5ec Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Mon, 11 Aug 2008 12:32:16 +0000
Subject: [PATCH] Updated for 4.1 release
---
fsl_regfilt.cc | 6 +-
ggmix.cc | 2 +-
ggmix.h | 2 +-
meldata.cc | 153 ++++++++++++++++++++++++-------------------------
meldata.h | 2 +-
melgmix.cc | 2 +-
melgmix.h | 2 +-
melhlprfns.cc | 2 +-
melhlprfns.h | 2 +-
melica.cc | 2 +-
melica.h | 2 +-
melodic.cc | 2 +-
melodic.h | 6 +-
meloptions.cc | 2 +-
meloptions.h | 4 +-
melpca.cc | 2 +-
melpca.h | 2 +-
melreport.cc | 10 ++--
melreport.h | 2 +-
test.cc | 2 +-
20 files changed, 104 insertions(+), 105 deletions(-)
diff --git a/fsl_regfilt.cc b/fsl_regfilt.cc
index 675d4cb..98f39de 100644
--- a/fsl_regfilt.cc
+++ b/fsl_regfilt.cc
@@ -1,8 +1,8 @@
/* fsl_regfilt -
- Christian Beckmann, FMRIB Image Analysis Group
+ Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 2006-2007 University of Oxford */
+ Copyright (C) 2006-2008 University of Oxford */
/* CCOPYRIGHT */
@@ -23,7 +23,7 @@ using namespace std;
// printed out as the help or usage message
string title=string("fsl_regfilt (Version 1.0)")+
- string("\n\n Copyright(c) 2007, University of Oxford (Christian F. Beckmann)\n")+
+ string("\n\n Copyright(c) 2008, 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]";
diff --git a/ggmix.cc b/ggmix.cc
index df6e4bd..3648396 100644
--- a/ggmix.cc
+++ b/ggmix.cc
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/ggmix.h b/ggmix.h
index 53af01e..8126e20 100644
--- a/ggmix.h
+++ b/ggmix.h
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/meldata.cc b/meldata.cc
index 4b2ffa3..28bf727 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
// {{{ includes/namespaces
@@ -54,8 +54,8 @@ namespace Melodic{
Matrix meanimg = convert_to_pbsc(tmpData);
meanR = meanimg.Row(1);
message(" done" << endl);
- }
- else{
+ }
+ else{
message(string(" Removing mean image ..."));
meanR = mean(tmpData);
tmpData = remmean(tmpData);
@@ -174,7 +174,7 @@ namespace Melodic{
Matrix alltcs = Smodes.at(0);
for(int ctr=1; ctr < (int)Smodes.size();ctr++)
alltcs|=Smodes.at(ctr);
- if((alltcs.Nrows()==Sdes.Nrows())&&(Sdes.Nrows()>Sdes.Ncols()))
+ if((alltcs.Nrows()==Sdes.Nrows())&&(Sdes.Nrows()>Sdes.Ncols()&&alltcs.Nrows()>2))
glmS.olsfit(alltcs,Sdes,Scon);
}
}
@@ -184,21 +184,21 @@ namespace Melodic{
numfiles = (int)opts.inputfname.value().size();
setup_misc();
- if(opts.filtermode){ // basic setup for filtering only
- Data = process_file(opts.inputfname.value().at(0));
- }
- else{
+ if(opts.filtermode){ // basic setup for filtering only
+ Data = process_file(opts.inputfname.value().at(0));
+ }
+ else{
- if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
- opts.approach.set_T("tica");
+ if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
+ opts.approach.set_T("tica");
- Matrix alldat, tmpData;
+ Matrix alldat, tmpData;
bool tmpvarnorm = opts.varnorm.value();
- if(opts.joined_vn.value()){
+ if(numfiles > 1 && opts.joined_vn.value()){
opts.varnorm.set_T(false);
}
- alldat = process_file(opts.inputfname.value().at(0), numfiles) / numfiles;
+ alldat = process_file(opts.inputfname.value().at(0), numfiles) / numfiles;
if(opts.pca_dim.value() > alldat.Nrows()-2){
cerr << "ERROR:: too many components selected \n\n";
@@ -208,51 +208,51 @@ 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);
if(opts.debug.value())
save4D(tmpData /numfiles,string("preproc_dat") + num2str(ctr+1));
if(tmpData.Ncols() == alldat.Ncols() && tmpData.Nrows() == alldat.Nrows())
- alldat += tmpData / numfiles;
+ alldat += tmpData / numfiles;
else
message("Data dimensions do not match - ignoring "+opts.inputfname.value().at(ctr) << endl);
- }
+ }
- //update mask
- if(opts.update_mask.value()){
- message("Excluding voxels with constant value ...");
- update_mask(Mask, alldat);
- message(" done" << endl);
- }
+ //update mask
+ if(opts.update_mask.value()){
+ message("Excluding voxels with constant value ...");
+ update_mask(Mask, alldat);
+ message(" done" << endl);
+ }
- if(opts.joined_vn.value()){
+ if((numfiles > 1 ) && opts.joined_vn.value() && tmpvarnorm){
//variance - normalisation
- message(endl<<"Normalising by voxel-wise variance ...");
- stdDev = varnorm(alldat,alldat.Nrows(),3.1);
- stdDevi = pow(stdDev,-1);
- message(" done" << endl);
+ message(endl<<"Normalising jointly by voxel-wise variance ...");
+ stdDev = varnorm(alldat,alldat.Nrows(),3.1);
+ stdDevi = pow(stdDev,-1);
+ message(" done" << endl);
}
if(numfiles>1)
- message(endl << "Initial data size : "<<alldat.Nrows()<<" x "<<alldat.Ncols()<<endl<<endl);
+ message(endl << "Initial data size : "<<alldat.Nrows()<<" x "<<alldat.Ncols()<<endl<<endl);
if(opts.debug.value())
save4D(alldat,"alldat");
- //estimate model order
- Matrix tmpPPCA;
- RowVector AdjEV, PercEV;
- Matrix Corr, tmpE;
- int order;
+ //estimate model order
+ Matrix tmpPPCA;
+ RowVector AdjEV, PercEV;
+ Matrix Corr, tmpE;
+ int order;
- order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
+ order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
- if(opts.pca_dim.value() == 0){
- opts.pca_dim.set_T(order);
+ if(opts.pca_dim.value() == 0){
+ opts.pca_dim.set_T(order);
PPCA=tmpPPCA;
- }
- order = opts.pca_dim.value();
- calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
+ }
+ order = opts.pca_dim.value();
+ calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
- if(opts.debug.value()){
+ if(opts.debug.value()){
outMsize("pcaE",pcaE); saveascii(pcaE,"pcaE");
outMsize("pcaD",pcaD); saveascii(pcaD,"pcaD");
outMsize("AdjEV",AdjEV); saveascii(AdjEV,"AdjEV");
@@ -266,26 +266,26 @@ namespace Melodic{
EV = AdjEV;
EVP = PercEV;
- if(numfiles == 1){
- Data = alldat;
- Matrix tmp = IdentityMatrix(Data.Nrows());
- DWM.push_back(tmp);
- WM.push_back(tmp);
- }
+ if(numfiles == 1){
+ Data = alldat;
+ Matrix tmp = IdentityMatrix(Data.Nrows());
+ DWM.push_back(tmp);
+ WM.push_back(tmp);
+ }
else {
- for(int ctr = 0; ctr < numfiles; ctr++){
+ for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
- if(opts.joined_vn.value()){
+ if(opts.joined_vn.value() && tmpvarnorm){
tmpData=SP(tmpData,pow(ones(tmpData.Nrows(),1)*stdDev,-1));
}
// whiten (separate / joint)
Matrix newWM,newDWM;
if(!opts.joined_whiten.value()){
message(" Individual whitening in a " << order << " dimensional subspace " << endl);
- std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
- calc_white(pcaE, pcaD, order, newWM, newDWM);
+ std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
+ calc_white(pcaE, pcaD, order, newWM, newDWM);
}else{
std_pca(whiteMatrix*tmpData, RXweight, Corr, pcaE, pcaD);
calc_white(pcaE, pcaD, order, newWM, newDWM);
@@ -301,19 +301,19 @@ namespace Melodic{
Data = tmpData;
else
Data &= tmpData;
- }
- }
+ }
+ }
opts.varnorm.set_T(tmpvarnorm);
- message(endl << " Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
+ message(endl << " Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
outMsize("stdDev",stdDev);
- //meanC=mean(Data,2);
+ //meanC=mean(Data,2);
if(opts.debug.value())
save4D(Data,"concat_data");
- //save the mean & mask
- save_volume(Mask,logger.appendDir("mask"));
- save_volume(Mean,logger.appendDir("mean"));
+ //save the mean & mask
+ save_volume(Mask,logger.appendDir("mask"));
+ save_volume(Mean,logger.appendDir("mean"));
}
} // void setup()
@@ -477,7 +477,7 @@ namespace Melodic{
Matrix tmp;
tmp=calc_FFT(dewhiteMatrix, opts.logPower.value());
saveascii(tmp,opts.outputfname.value() + "_FTdewhite");
- }
+ }
//Output PCA
if(opts.output_pca.value() && pcaD.Storage()>0&&pcaE.Storage()>0){
@@ -494,7 +494,7 @@ namespace Melodic{
}
message("...done" << endl);
- } //void save()
+ } //void save()
int MelodicData::remove_components()
{
@@ -528,7 +528,7 @@ namespace Melodic{
noiseMix = mixMatrix.Column(ctr);
noiseIC = IC.Row(ctr).t();
}
- else{
+ else{
cerr << endl<< "component number "<<ctr<<" does not exist" << endl;
}
@@ -614,8 +614,7 @@ namespace Melodic{
}
}
- unsigned long MelodicData::standardise(volume<float>& mask,
- volume4D<float>& R){
+ unsigned long MelodicData::standardise(volume<float>& mask, volume4D<float>& R){
unsigned long count = 0;
int M=R.tsize();
@@ -653,7 +652,7 @@ namespace Melodic{
}
}
return count;
- }
+ }
float MelodicData::est_resels(volume4D<float> R, volume<float> mask)
{
@@ -696,7 +695,7 @@ namespace Melodic{
if (usez) S2[Z] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) +
R(x, y, z-1, t)*R(x, y, z-1, t));
}
- }
+ }
float norm = 1.0/(float) N;
float v = dof; // v - degrees of freedom (nu)
@@ -806,28 +805,28 @@ namespace Melodic{
void MelodicData::sort(){
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
- numTime = mixMatrix.Nrows(), i,j;
+ numTime = mixMatrix.Nrows(), i,j;
- //flip IC maps to be positive (on average)
- //flip Subject/Session modes to be positive (on average)
- //have time courses accordingly
+ //flip IC maps to be positive (on max)
+ //flip Subject/Session modes to be positive (on average)
+ //flip time courses accordingly
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++)
if(IC.Row(ctr_i).MaximumAbsoluteValue()>IC.Row(ctr_i).Maximum()){
flipres(ctr_i);
- }
+ }
message("Sorting IC maps" << endl);
Matrix tmpscales, tmpICrow, tmpMIXcol;
- if(numfiles > 1){
- set_TSmode();
- Matrix allmodes = Smodes.at(0);
- for(int ctr = 1; ctr < (int)Smodes.size();++ctr)
- allmodes |= Smodes.at(ctr);
- tmpscales = median(allmodes).t();
- } else {
+ if(numfiles > 1){
+ set_TSmode();
+ Matrix allmodes = Smodes.at(0);
+ for(int ctr = 1; ctr < (int)Smodes.size();++ctr)
+ allmodes |= Smodes.at(ctr);
+ tmpscales = median(allmodes).t();
+ } else {
// re-order wrt standard deviation of IC maps
tmpscales = stdev(IC,2);
- }
+ }
ICstats = tmpscales;
diff --git a/meldata.h b/meldata.h
index d691e93..5944013 100644
--- a/meldata.h
+++ b/meldata.h
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melgmix.cc b/melgmix.cc
index 66d3d05..9b715bf 100644
--- a/melgmix.cc
+++ b/melgmix.cc
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melgmix.h b/melgmix.h
index cf0d966..0cf5792 100644
--- a/melgmix.h
+++ b/melgmix.h
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 5954a48..1a81a79 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melhlprfns.h b/melhlprfns.h
index 229e4a8..1d51cab 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melica.cc b/melica.cc
index cd53299..c80b17d 100644
--- a/melica.cc
+++ b/melica.cc
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melica.h b/melica.h
index 26fdc5e..af098f3 100644
--- a/melica.h
+++ b/melica.h
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melodic.cc b/melodic.cc
index fef92f1..4935a2f 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melodic.h b/melodic.h
index 2c3c009..2fcfaf5 100644
--- a/melodic.h
+++ b/melodic.h
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
@@ -40,13 +40,13 @@
namespace Melodic{
-const string version = "3.05";
+const string version = "3.09";
// The two strings below specify the title and example usage that is
// printed out as the help or usage message
const string title=string("MELODIC (Version ")+version+")"+
string("\n Multivariate Exploratory Linear Optimised Decomposition into Independent Components\n")+
- string(" Copyright(c) 2007, University of Oxford (Christian F. Beckmann)");
+ string(" Copyright(c) 2001-2008, University of Oxford (Christian F. Beckmann)");
const string usageexmpl=string(" melodic -i <filename> <options>")+
string("\n \t \t to run melodic")+
diff --git a/meloptions.cc b/meloptions.cc
index 0b92f4d..6415055 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -6,7 +6,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/meloptions.h b/meloptions.h
index ca1eec8..1db4589 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -6,7 +6,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
@@ -188,7 +188,7 @@ class MelodicOptions {
false, no_argument),
joined_vn(string("--sep_vn"), true,
string("switch off joined variance nomalisation"),
- false, no_argument, false),
+ false, no_argument),
vn_level(string("--vn_level"), float(2.3),
string("variance nomalisation threshold level (Z> value is ignored)"),
false, requires_argument, false),
diff --git a/melpca.cc b/melpca.cc
index 3819945..4c1edb9 100644
--- a/melpca.cc
+++ b/melpca.cc
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melpca.h b/melpca.h
index ce68b89..57ff0ed 100644
--- a/melpca.h
+++ b/melpca.h
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/melreport.cc b/melreport.cc
index b851f86..6110a19 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -1,11 +1,11 @@
/* MELODIC - Multivariate exploratory linear optimized decomposition into
independent components
- melodat.get_bg()
+
melreport.cc - report generation
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
@@ -50,9 +50,9 @@ namespace Melodic{
}
{//output IC stats
if(ICstats.Storage()>0&&ICstats.Nrows()>=cnum){
- IChtml << fixed << setprecision(2) << ICstats(cnum,1) << " % of explained variance";
+ IChtml << fixed << setprecision(2) << std::abs(ICstats(cnum,1)) << " % of explained variance";
if(ICstats.Ncols()>1)
- IChtml << "; " << ICstats(cnum,2) << " % of total variance";
+ IChtml << "; " << std::abs(ICstats(cnum,2)) << " % of total variance";
if(ICstats.Ncols()>2&&opts.addsigchng.value()){
IChtml << "<p>" <<endl;
IChtml << " " << ICstats(cnum,3) << " % signal change (pos peak voxel); " << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;
@@ -161,7 +161,7 @@ namespace Melodic{
if(melodat.get_numfiles()>1 && melodat.explained_var.Storage()>0
&& melodat.explained_var.Ncols()>=cnum && opts.varvals.value())
IChtml << "Rank-1 approximation of individual time courses explains "
- << melodat.explained_var(cnum) << "% of variance.<p>" << endl;
+ << std::abs(melodat.explained_var(cnum)) << "% of variance.<p>" << endl;
}//time series plot
if(!opts.pspec.value())
diff --git a/melreport.h b/melreport.h
index 9be8f82..115ea78 100644
--- a/melreport.h
+++ b/melreport.h
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
diff --git a/test.cc b/test.cc
index 33ff5a6..b51a5b4 100644
--- a/test.cc
+++ b/test.cc
@@ -1,6 +1,6 @@
/* test.cc
- Copyright (C) 1999-2007 University of Oxford */
+ Copyright (C) 1999-2008 University of Oxford */
/* CCOPYRIGHT */
--
GitLab