diff --git a/meldata.cc b/meldata.cc
index 35549c09a8f5390725718db5a8928b5275ce0c20..bf042cf309470b03635fdefcb3f99b5924cf03a2 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -83,22 +83,24 @@ namespace Melodic{
// meanR=zeros(1,Data.Ncols());
//}
}
- {//switch dimension in case temporal ICA is required
- // if(opts.temporal.value()){
- // if(opts.remove_meanvol.value()){
- // message(string("Remove mean time course for temporal ICA ... "));
- // meanC=meanR;
- // meanR=mean(Data);
- // Data=remmean(Data,2);
- // message("done"<<endl);
- // }
- // Data = Data.t();
- // }
- }
{// remove mean time course
meanC=mean(Data,2);
Data=remmean(Data,2);
}
+
+ {//switch dimension in case temporal ICA is required
+ if(opts.temporal.value()){
+ message(string("Switching dimensions for temporal ICA") << endl);
+ Data = Data.t();
+ Matrix tmp;
+ tmp = meanC;
+ meanC = meanR.t();
+ meanR = tmp.t();
+ message("Data size : " << Data.Nrows() << " x " << Data.Ncols() <<endl);
+ }
+ }
+
+
{// variance-normalize the data
DiagonalMatrix tmpD;Matrix tmpE;
message(string("Estimating data covariance ... "));
@@ -126,24 +128,35 @@ namespace Melodic{
}
}
- stdDevi = pow(stdev(Data - tmpDeWhite*WS),-1);
+ stdDevi = pow(stdev(Data - tmpDeWhite*WS),-1);
Data = Data + meanC*ones(1,Data.Ncols());
}
DataVN = SP(Data,ones(Data.Nrows(),1)*stdDevi);
+
+ if(opts.output_all.value()){
+ volume4D<float> tempVol;
+ tempVol.setmatrix(stdDevi,Mask);
+ save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
+ + "_vn_stdev"),tempInfo);
+ tempVol.setmatrix(DataVN,Mask);
+ save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
+ + "_vn"),tempInfo);
+ }
+
if(opts.varnorm.value()){
Data = DataVN;
message("done"<<endl);
}
{//remove row mean
if(opts.temporal.value()){
- Data=remmean(Data,2);
+ message(string("Removing mean image ... "));
}else{
message(string("Removing mean time course ... "));
- meanC=mean(Data,2);
- Data=remmean(Data,2);
- message("done"<<endl);
}
+ meanC=mean(Data,2);
+ Data=remmean(Data,2);
+ message("done"<<endl);
}
if(opts.segment.value().length()>0){
@@ -164,9 +177,20 @@ namespace Melodic{
//check for temporal ICA
if(opts.temporal.value()){
+ message(string("temporal ICA: transform back the data ... "));
Matrix tmpIC = mixMatrix.t();
mixMatrix=IC.t();
IC=tmpIC;
+
+ unmixMatrix=pinv(mixMatrix);
+ Data = Data.t();
+ tmpIC = meanC;
+ meanC = meanR.t();
+ meanR = tmpIC.t();
+ // whiteMatrix = whiteMatrix.t;
+ // dewhiteMatrix = dewhiteMatrix.t();
+ message(string("done") << endl);
+ opts.temporal.set_T(false); // Do not switch again!
}
message("Writing results to : " << endl);
@@ -188,19 +212,27 @@ namespace Melodic{
// read_volume4D(tempVol,opts.inputfname.value(),tempInfo);
//Matrix ICadjust = IC;
- if(opts.temporal.value()){
- Data=Data.t();
- }
-
+
Matrix ICadjust;
- stdNoisei = pow(stdev(Data - mixMatrix * IC),-1);
- ICadjust = SP(IC,ones(IC.Nrows(),1)*stdNoisei);
+ stdNoisei = pow(stdev(Data - mixMatrix * IC)*std::sqrt((float)(Data.Nrows()-1))/
+ std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
+
+ ColumnVector diagvals;
+ diagvals=pow(diag(unmixMatrix*unmixMatrix.t()),-0.5);
+ ICadjust = SP(IC,diagvals*stdNoisei);
tempVol.setmatrix(ICadjust,Mask);
//strncpy(tempInfo.header.hist.aux_file,"render3",24);
save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
+ "_IC"),tempInfo);
message(" " << logger.appendDir(opts.outputfname.value() + "_IC") <<endl);
+
+ if(opts.output_origIC.value()){
+ tempVol.setmatrix(stdNoisei,Mask);
+ save_volume4D(tempVol,logger.appendDir(string("Noise_stddev_inv")),tempInfo);
+ message(" " << logger.appendDir(string("Noise_stddev_inv")) <<endl);
+ }
+
}
//Output mixMatrix
@@ -257,19 +289,22 @@ namespace Melodic{
if(opts.output_pca.value() && pcaD.Storage()>0&&pcaE.Storage()>0){
write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_pcaE"),
pcaE);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_pcaE") <<endl);
write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_pcaD"),
(Matrix) diag(pcaD));
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_pcaD") <<endl);
+
volume4D<float> tempVol;
- //volumeinfo tempInfo;
- //read_volume4D(tempVol,opts.inputfname.value(),tempInfo);
- Matrix ICadjust = IC;
- if(opts.temporal.value()){
- Data=Data.t();
- }
-
Matrix PCAmaps;
- PCAmaps = whiteMatrix * Data;
+
+ if(whiteMatrix.Ncols()==Data.Ncols()){
+ PCAmaps = dewhiteMatrix.t();
+ }else
+ PCAmaps = whiteMatrix * Data;
+
tempVol.setmatrix(PCAmaps,Mask);
//strncpy(tempInfo.header.hist.aux_file,"render3",24);
@@ -277,11 +312,8 @@ namespace Melodic{
+ "_pca"),tempInfo);
message(" " <<
logger.appendDir(opts.outputfname.value() + "_pca") <<endl);
- message(" "<<
- logger.appendDir(opts.outputfname.value() + "_pcaE") <<endl);
- message(" "<<
- logger.appendDir(opts.outputfname.value() + "_pcaD") <<endl);
- }
+
+ }
} //void save()
int MelodicData::remove_components()
diff --git a/meldata.h b/meldata.h
index 9749f08208a529cfb9d72810207c3739aacd2c30..f3cd6e0a7c0f96dacfb156245963a2e3b14dc21b 100644
--- a/meldata.h
+++ b/meldata.h
@@ -59,6 +59,7 @@ namespace Melodic{
inline Matrix& get_IC() {return IC;}
inline void set_IC(Matrix& Arg) {IC = Arg;}
+ inline void set_IC(int ctr, Matrix& Arg) {IC.Row(ctr) = Arg;}
inline Matrix& get_white() {return whiteMatrix;}
inline void set_white(Matrix& Arg) {whiteMatrix = Arg;}
@@ -73,7 +74,10 @@ namespace Melodic{
inline void set_stdDevi(Matrix& Arg) {stdDevi = Arg;}
inline Matrix& get_mix() {return mixMatrix;}
- inline void set_mix(Matrix& Arg) {mixMatrix = Arg;}
+ inline void set_mix(Matrix& Arg)
+ {mixMatrix = Arg; unmixMatrix=pinv(mixMatrix);}
+ inline void set_mix(int ctr, Matrix& Arg)
+ {mixMatrix.Column(ctr) = Arg; unmixMatrix=pinv(mixMatrix);}
inline Matrix& get_fmix() {return mixFFT;}
inline void set_fmix(Matrix& Arg) {mixFFT = Arg;}
diff --git a/melgmix.cc b/melgmix.cc
index b586365aeb2b474d55678f77783679e2c8343423..9b142cdba96ae2de087fbd57c0dfbe03844080d8 100644
--- a/melgmix.cc
+++ b/melgmix.cc
@@ -83,7 +83,11 @@ namespace Melodic{
fitted = false;
nummix = num_mix;
numdata = dat.Ncols();
- data=dat;
+
+ //normalise data
+ datamean = mean(dat,2).AsScalar();
+ datastdev= stdev(dat,2).AsScalar();
+ data=(dat - datamean)/datastdev;
props=zeros(1,nummix);
vars=zeros(1,nummix);
@@ -181,6 +185,7 @@ namespace Melodic{
//cerr << " levles : " << levls << endl << endl;
Res = threshold(data, levls);
set_threshmaps(Res);
+
return Res;
}
diff --git a/melgmix.h b/melgmix.h
index e31d99ab412a851287cdd635ad1936fee592fe88..9adcc0eaccc0939f83c411f45bd3b42873bed617 100644
--- a/melgmix.h
+++ b/melgmix.h
@@ -64,6 +64,13 @@ namespace Melodic{
this->ggmfit();
else
this->gmmfit();
+
+ //re-insert mean and stdev
+
+ data = data*datastdev + datamean;
+ //threshmaps = threshmaps*datastdev + datamean;
+ means = means*datastdev + datamean;
+ vars = vars*datastdev*datastdev;
}
inline Matrix threshold(string levels)
@@ -143,6 +150,9 @@ namespace Melodic{
threshinfo.clear();
}
+ double datamean;
+ double datastdev;
+
private:
MelodicOptions &opts;
Log &logger; //global log file
diff --git a/melodic.cc b/melodic.cc
index 5001ba3b79b7b4f1541db2efdccb5556f7efd062..206ae552879fb9dfc757168454dc31914e503c5a 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -5,7 +5,7 @@
Christian F. Beckmann, FMRIB Image Analysis Group
- Copyright (C) 1999-2002 University of Oxford */
+ Copyright (C) 1999-2003 University of Oxford */
/* CCOPYRIGHT */
@@ -49,8 +49,7 @@ string myfloat2str(float f, int width, int prec, bool scientif)
Matrix mmall(Log& logger, MelodicOptions& opts,
MelodicData& melodat, MelodicReport& report, Matrix& probs);
-void repall(Log& logger, MelodicOptions& opts, MelodicData& melodat,
- MelodicReport& report, Matrix mmres, Matrix probs);
+
void mmonly(Log& logger, MelodicOptions& opts,
MelodicData& melodat, MelodicReport& report);
@@ -135,7 +134,17 @@ int main(int argc, char *argv[])
} while (no_conv && retry<opts.maxRestart.value() && !leaveloop);
if(!no_conv){
- //first save raw IC results
+
+ //re-normalise the decomposition to std(mix)=1
+ Matrix scales;
+ scales = stdev(melodat.get_mix());
+ Matrix tmp;
+ tmp = SP(melodat.get_mix(),ones(melodat.get_mix().Nrows(),1)*pow(scales,-1));
+ melodat.set_mix(tmp);
+ tmp = SP(melodat.get_IC(),scales.t()*ones(1,melodat.get_IC().Ncols()));
+ melodat.set_IC(tmp);
+
+ //save raw IC results
melodat.save();
Matrix pmaps;//(melodat.get_IC());
@@ -162,14 +171,12 @@ int main(int argc, char *argv[])
}
}
+
+
if( bool(opts.genreport.value()) ){
report.analysistxt();
report.PPCA_rep();
}
- //cerr << mmres.Nrows() << " x " << mmres.Ncols() << endl;
-
- // if(opts.genreport.value())
- // repall(logger,opts,melodat,report,mmres,pmaps,threshmaps);
message("finished!" << endl << endl);
} else {
@@ -271,26 +278,6 @@ void mmonly(Log& logger, MelodicOptions& opts,
mmres = mmall(logger,opts,melodat,report,pmaps);
}
-void repall(Log& logger, MelodicOptions& opts, MelodicData& melodat,
- MelodicReport& report, Matrix& mmpars, Matrix& probs)
-{
- if( bool(opts.genreport.value()) ){
-
- message(endl
- << "Creating report in " << report.getDir()
- << endl);
-
- for(int ctr=1; ctr<=probs.Nrows(); ctr++){
- message(" " << ctr);
- MelGMix mixmod(opts, logger);
-
- //load MelGMix
-
- //report.IC_rep(mixmod,ctr,melodat.get_IC().Nrows());
-
- }
- }
-}
Matrix mmall(Log& logger, MelodicOptions& opts,
MelodicData& melodat, MelodicReport& report, Matrix& pmaps)
@@ -298,8 +285,6 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
Matrix mmpars(5*melodat.get_IC().Nrows(),5);
mmpars = 0;
- //Matrix pmaps(melodat.get_IC());
-
Log stats;
if(opts.output_MMstats.value()){
@@ -309,7 +294,10 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
message(endl
<< "Running Mixture Modelling on Z-transformed IC maps ..."
<< endl);
-
+
+ ColumnVector diagvals;
+ diagvals=pow(diag(melodat.get_unmix()*melodat.get_unmix().t()),-0.5);
+
for(int ctr=1; ctr <= melodat.get_IC().Nrows(); ctr++){
MelGMix mixmod(opts, logger);
@@ -317,8 +305,9 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
Matrix ICmap;
- if(melodat.get_stdNoisei().Storage()>0)
- ICmap = SP(melodat.get_IC().Row(ctr),melodat.get_stdNoisei());
+ if(opts.varnorm.value()&&melodat.get_stdNoisei().Storage()>0){
+ ICmap = SP(melodat.get_IC().Row(ctr),diagvals(ctr)*melodat.get_stdNoisei());
+ }
else
ICmap = melodat.get_IC().Row(ctr);
@@ -345,28 +334,44 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
//else
// pmaps.Row(ctr) = zeros(1,pmaps.Ncols());
+ //re-scale spatial maps to mean 0 for nht
+
+ if(opts.rescale_nht.value()){
+ message(" re-scaling spatial maps ... "<< endl);
+ RowVector tmp;
+ tmp = mixmod.get_means();
+ float dmean = tmp(1);
+ tmp = mixmod.get_vars();
+ float dstdev = sqrt(tmp(1));
+
+ tmp = (mixmod.get_means() - dmean)/dstdev;
+ mixmod.set_means(tmp);
+ tmp = (mixmod.get_vars() / (dstdev*dstdev));
+ mixmod.set_vars(tmp);
+
+ //tmp = (mixmod.get_data() - dmean)/dstdev;
+ tmp = (ICmap - dmean)/dstdev;
+ mixmod.set_data(tmp);
+ //if(opts.varnorm.value()&&melodat.get_stdNoisei().Storage()>0)
+ // tmp = SP(tmp,pow(diagvals(ctr)*melodat.get_stdNoisei(),-1));
+
+ melodat.set_IC(ctr,tmp);
+ }
+
message(" thresholding ... "<< endl);
mixmod.threshold(opts.mmthresh.value());
- //re-orient the data
-
- //message(" done " << endl);
Matrix tmp;
tmp=(mixmod.get_threshmaps().Row(1));
float posint = SP(tmp,gt(tmp,zeros(1,tmp.Ncols()))).Sum();
float negint = -SP(tmp,lt(tmp,zeros(1,tmp.Ncols()))).Sum();
-
- //cerr << posint << " " << negint << endl;
if((posint<0.01)&&(negint<0.01)){
mixmod.clear_infstr();
- //cerr << "after infstr"<<endl;
mixmod.threshold("0.05n "+opts.mmthresh.value());
- //cerr << " back again" << endl;
posint = SP(tmp,gt(tmp,zeros(1,tmp.Ncols()))).Sum();
negint = -SP(tmp,lt(tmp,zeros(1,tmp.Ncols()))).Sum();
}
- //cerr << posint << " " << negint << endl;
if(negint>posint){//flip map
melodat.flipres(ctr);
mixmod.flipres(ctr);
@@ -378,8 +383,6 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
<< " Means : " << mixmod.get_means() << endl
<< " Vars. : " << mixmod.get_vars() << endl
<< " Prop. : " << mixmod.get_pi() << endl << endl;
- // << " Offs. : " << mixmod.get_offset() << endl << endl;
- //cerr << mixmod.get_threshmaps().Nrows() << " " << mixmod.get_threshmaps().Ncols()<< endl;
melodat.save4D(mixmod.get_threshmaps(),
string("stats/thresh_zstat")+num2str(ctr));
}
@@ -410,17 +413,27 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
message(" done" << endl);
}
}
- if( bool(opts.genreport.value()) ){
+
+ if(!opts.filtermode&&opts.filtermix.value().length()==0){
+ //now safe new data
+ // bool what = opts.verbose.value();
+ //opts.verbose.set_T(false);
+ // melodat.save();
+ //opts.verbose.set_T(what);
+
+ if(melodat.get_IC().Storage()>0){
+ volume4D<float> tempVol;
+ tempVol.setmatrix(melodat.get_IC(),melodat.get_mask());
+ save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
+ + "_IC"),melodat.tempInfo);
+ message(endl<< endl << " Saving " << logger.appendDir(opts.outputfname.value() + "_IC") <<endl);
+ }
+
+ if( bool(opts.genreport.value()) ){
message(endl << endl <<
" To view the output report point your web browser at " <<
report.getDir() + "/00index.html" << endl << endl);
- }
- if(!opts.filtermode&&opts.filtermix.value().length()==0){
- //now safe new data
- bool what = opts.verbose.value();
- opts.verbose.set_T(false);
- melodat.save();
- opts.verbose.set_T(what);
+ }
}
return mmpars;
}
diff --git a/melodic.h b/melodic.h
index c18372642341ae7e5d0142c629b2c1d5672ceb9d..afa4b1858d6e48fa44611bffe9a465713f7d528b 100644
--- a/melodic.h
+++ b/melodic.h
@@ -23,11 +23,16 @@
} \
Log& logger = LogSingleton::getInstance(); \
logger.str() << msg; \
+ cout.flush(); \
+}
+
+#define outMsize(msg,Mat) { \
+ cout << msg << " " <<Mat.Nrows() << " x " << Mat.Ncols() << endl; \
}
namespace Melodic{
-const string version = "2.0";
+const string version = "ln(10)";
}
diff --git a/meloptions.cc b/meloptions.cc
index a79c6001abb29ebbd1ae97118e31f7155a512c13..743aa23f3104e999bdc8f0871eca2bb556f9db95 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -126,6 +126,14 @@ MelodicOptions* MelodicOptions::gopt = NULL;
if (threshold.value()<=0){
use_mask.set_T(false);
}
+ if (output_all.value()){
+ output_unmix.set_T(true);
+ output_MMstats.set_T(true);
+ output_pca.set_T(true);
+ output_white.set_T(true);
+ output_origIC.set_T(true);
+ output_mean.set_T(true);
+ }
if (output_pca.value()){
output_white.set_T(true);
}
diff --git a/meloptions.h b/meloptions.h
index 17d85e1f3b2443fd936996507569596e57ecde31..1a2377550157d56ee558693073baa50e7dbccdfa 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -73,6 +73,7 @@ class MelodicOptions {
Option<float> tr;
Option<bool> logPower;
+ Option<bool> output_all;
Option<bool> output_unmix;
Option<bool> output_MMstats;
Option<bool> output_pca;
@@ -96,6 +97,7 @@ class MelodicOptions {
Option<bool> remove_meanvol;
Option<bool> remove_endslices;
+ Option<bool> rescale_nht;
Option<bool> guess_remderiv;
Option<bool> temporal;
@@ -205,6 +207,9 @@ class MelodicOptions {
logPower(string("--logPower"), false,
string("calculate log of power for frequency spectrum\n"),
false, no_argument),
+ output_all(string("--Oall"), false,
+ string("output everything"),
+ false, no_argument),
output_unmix(string("--Ounmix"), false,
string("output unmixing matrix"),
false, no_argument),
@@ -259,6 +264,9 @@ class MelodicOptions {
remove_endslices(string("--remEndslices"), false,
string("delete end slices (motion correction artefacts)"),
false, no_argument,false),
+ rescale_nht(string("--rescale_nht"), true,
+ string("switch off map rescaling after mixture-modelling"),
+ false, no_argument,false),
guess_remderiv(string("--remderiv"), false,
string("removes every second entry in paradigm file (EV derivatives)"),
false, no_argument, false),
@@ -305,6 +313,7 @@ class MelodicOptions {
options.add(genreport);
options.add(tr);
options.add(logPower);
+ options.add(output_all);
options.add(output_unmix);
options.add(output_MMstats);
options.add(output_pca);
@@ -324,6 +333,7 @@ class MelodicOptions {
options.add(nlconst2);
options.add(remove_meanvol);
options.add(remove_endslices);
+ options.add(rescale_nht);
options.add(guess_remderiv);
options.add(temporal);
}
diff --git a/melreport.cc b/melreport.cc
index a1b9c370d70b4bc1517c1584797ca85c4866f430..b39c0bc680e97926dc146c43cbd96a3c312c9c9a 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -266,7 +266,9 @@ namespace Melodic{
if(melodat.get_IC().Storage()>0)
{//Output raw IC map
- tempVol.setmatrix(melodat.get_IC().Row(cnum),
+ // tempVol.setmatrix(melodat.get_IC().Row(cnum),
+ // melodat.get_mask());
+ tempVol.setmatrix(mmodel.get_data(),
melodat.get_mask());
volume<float> map1;
volume<float> map2;
diff --git a/test.cc b/test.cc
index 51b0ee862d56dc55368463e15e0e84232b6544cc..3c8339fa099472358527d3d7bdfc2717c6b9b935 100644
--- a/test.cc
+++ b/test.cc
@@ -151,10 +151,6 @@ mmat.setDir(string("./res/"),string("mix.txt"));
{//plot frequency
miscplot newplot;
int fact = int(std::pow(10.0,int(std::log10(float(data.Nrows())))));
-
-
- cerr << " FACT : "<<fact;
-
if(tr>0.0)
newplot.timeseries(data2.Column(i).t(),