diff --git a/meldata.cc b/meldata.cc
index 79dd63804ec3abd5115330fe4b2b8abc41bcb422..dfc440f7bcaced63194621955aa6f1cc9cabcf0c 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -181,10 +181,10 @@ namespace Melodic{
}
else{
order = opts.pca_dim.value();
- std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
+ std_pca(alldat, RXweight, Corr, pcaE, pcaD);
calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
}
-
+
if(numfiles < 2){
Data = alldat;
Matrix tmp = Identity(Data.Nrows());
@@ -193,7 +193,7 @@ namespace Melodic{
} else {
for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
-
+
// whiten (separate / joint)
if(!opts.joined_whiten.value()){
std_pca(tmpData, RXweight, Corr, pcaE, pcaD);
@@ -248,6 +248,7 @@ namespace Melodic{
char callRMstr[1000];
ostrstream osc(callRMstr,1000);
osc << "rm " << string(Mean_fname) <<"* " << '\0';
+
system(callRMstr);
if(!samesize(Mean,Mask)){
@@ -266,9 +267,24 @@ namespace Melodic{
//seed the random number generator
double tmptime = time(NULL);
srand((unsigned int) tmptime);
- }
-
+ //read in post-proc design matrices etc
+ if(opts.fn_Tdesign.value().length()>0)
+ Tdes = read_vest(opts.fn_Tdesign.value());
+ if(opts.fn_Sdesign.value().length()>0)
+ Sdes = read_vest(opts.fn_Sdesign.value());
+ if(opts.fn_Tcon.value().length()>0)
+ Tcon = read_vest(opts.fn_Tcon.value());
+ if(opts.fn_Scon.value().length()>0)
+ Scon = read_vest(opts.fn_Scon.value());
+ if(opts.fn_TconF.value().length()>0)
+ TconF = read_vest(opts.fn_TconF.value());
+ if(opts.fn_SconF.value().length()>0)
+ SconF = read_vest(opts.fn_SconF.value());
+
+ Tdes = remmean(Tdes,1);
+ Sdes = remmean(Sdes,1);
+ }
// {{{ Save
void MelodicData::save()
@@ -334,7 +350,7 @@ namespace Melodic{
//Output mixMatrix
if(mixMatrix.Storage()>0){
saveascii(mixMatrix, opts.outputfname.value() + "_mix");
- mixFFT=calc_FFT(mixMatrix, opts.logPower.value());
+ mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
saveascii(mixFFT,opts.outputfname.value() + "_FTmix");
}
@@ -746,7 +762,7 @@ namespace Melodic{
ICstats |= copeN;
}
- mixFFT=calc_FFT(mixMatrix, opts.logPower.value());
+ mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
unmixMatrix = pinv(mixMatrix);
//if(ICstats.Storage()>0){cout << "ICstats: " << ICstats.Nrows() <<"x" << ICstats.Ncols() << endl;}else{cout << "ICstats empty " <<endl;}
diff --git a/meldata.h b/meldata.h
index f860cff7fbca084cf0eccc9561138eb701941ba1..cb9eae7adc876e6ff10010a430031fccff5e1e8e 100644
--- a/meldata.h
+++ b/meldata.h
@@ -29,32 +29,31 @@ namespace Melodic{
//constructor
MelodicData(MelodicOptions &popts, Log &plogger):
- opts(popts),
- logger(plogger)
- {
- after_mm = false;
- Resels = 0;
- }
+ opts(popts),logger(plogger)
+ {
+ after_mm = false;
+ Resels = 0;
+ }
void save();
Matrix process_file(string fname, int numfiles);
inline void save4D(Matrix what, string fname){
- volume4D<float> tempVol;
- tempVol.setmatrix(what,Mask);
- save_volume4D(tempVol,logger.appendDir(fname),tempInfo);
- message(" " << logger.appendDir(fname) << endl);
+ volume4D<float> tempVol;
+ tempVol.setmatrix(what,Mask);
+ save_volume4D(tempVol,logger.appendDir(fname),tempInfo);
+ message(" " << logger.appendDir(fname) << endl);
}
inline void saveascii(Matrix what, string fname){
- write_ascii_matrix(logger.appendDir(fname),what);
- message(" " << logger.appendDir(fname) << endl);
+ write_ascii_matrix(logger.appendDir(fname),what);
+ message(" " << logger.appendDir(fname) << endl);
}
inline void savebinary(Matrix what, string fname){
- write_binary_matrix(what,logger.appendDir(fname));
- message(" " << logger.appendDir(fname) << endl);
+ write_binary_matrix(what,logger.appendDir(fname));
+ message(" " << logger.appendDir(fname) << endl);
}
int remove_components();
@@ -78,20 +77,20 @@ namespace Melodic{
inline Matrix& get_Smodes(int what) {return Smodes.at(what);}
inline void add_Smodes(Matrix& Arg) {Smodes.push_back(Arg);}
inline void save_Smodes(){
- Matrix tmp = Smodes.at(0);
- for(unsigned int ctr = 1; ctr < Smodes.size(); ctr++)
- tmp |= Smodes.at(ctr);
- saveascii(tmp,opts.outputfname.value() + "_Smodes");
+ Matrix tmp = Smodes.at(0);
+ for(unsigned int ctr = 1; ctr < Smodes.size(); ctr++)
+ tmp |= Smodes.at(ctr);
+ saveascii(tmp,opts.outputfname.value() + "_Smodes");
}
inline vector<Matrix>& get_Tmodes() {return Tmodes;}
inline Matrix& get_Tmodes(int what) {return Tmodes.at(what);}
inline void add_Tmodes(Matrix& Arg) {Tmodes.push_back(Arg);}
inline void save_Tmodes(){
- Matrix tmp = Tmodes.at(0);
- for(unsigned int ctr = 1; ctr < Tmodes.size(); ctr++)
- tmp |= Tmodes.at(ctr);
- saveascii(tmp,opts.outputfname.value() + "_Tmodes");
+ Matrix tmp = Tmodes.at(0);
+ for(unsigned int ctr = 1; ctr < Tmodes.size(); ctr++)
+ tmp |= Tmodes.at(ctr);
+ saveascii(tmp,opts.outputfname.value() + "_Tmodes");
}
void set_TSmode();
@@ -111,12 +110,12 @@ namespace Melodic{
inline Matrix& get_mix() {return mixMatrix;}
inline void set_mix(Matrix& Arg) {
- mixMatrix = Arg;
- if (Tmodes.size() < 1)
- for (int ctr = 1; ctr <= Arg.Ncols(); ctr++){
- Matrix tmp = Arg.Column(ctr);
- add_Tmodes(tmp);
- }
+ mixMatrix = Arg;
+ if (Tmodes.size() < 1)
+ for (int ctr = 1; ctr <= Arg.Ncols(); ctr++){
+ Matrix tmp = Arg.Column(ctr);
+ add_Tmodes(tmp);
+ }
}
Matrix expand_mix();
@@ -167,24 +166,24 @@ namespace Melodic{
inline void set_after_mm(bool val) {after_mm = val;}
inline void flipres(int num){
- IC.Row(num) = -IC.Row(num);
- mixMatrix.Column(num) = -mixMatrix.Column(num);
- mixFFT=calc_FFT(mixMatrix);
- unmixMatrix = pinv(mixMatrix);
+ IC.Row(num) = -IC.Row(num);
+ mixMatrix.Column(num) = -mixMatrix.Column(num);
+ mixFFT=calc_FFT(mixMatrix);
+ unmixMatrix = pinv(mixMatrix);
if(ICstats.Storage()>0&&ICstats.Ncols()>3){
- double tmp;
- tmp = ICstats(num,3);
- ICstats(num,3) = -1.0*ICstats(num,4);
- ICstats(num,4) = -1.0*tmp;
- }
+ double tmp;
+ tmp = ICstats(num,3);
+ ICstats(num,3) = -1.0*ICstats(num,4);
+ ICstats(num,4) = -1.0*tmp;
+ }
}
void sort();
volumeinfo tempInfo;
-
vector<Matrix> DWM, WM;
+ basicGLM glmT, glmS;
private:
MelodicOptions &opts;
@@ -204,6 +203,7 @@ namespace Melodic{
Matrix mixFFT;
Matrix IC;
Matrix ICstats;
+ Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;
vector<Matrix> Tmodes;
vector<Matrix> Smodes;
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 6d2c371b86f4283e84b4d1f6510eec6e74f9eab5..6c881885d69cb2fc0fee31248f6907542cd3f61c 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -12,6 +12,8 @@
#include "melhlprfns.h"
#include "libprob.h"
#include "miscmaths/miscprob.h"
+#include "miscmaths/t2z.h"
+#include "miscmaths/f2z.h"
namespace Melodic{
@@ -781,4 +783,52 @@ namespace Melodic{
return res;
} //Matrix gen_arCorr
+ Matrix basicGLM::olsfit(const Matrix& data, const Matrix& design,
+ const Matrix& contrasts, int DOFadjust)
+ {
+ beta = zeros(design.Ncols(),1);
+ residu = zeros(1); sigsq = -1.0*ones(1); varcb = -1.0*ones(1);
+ t = zeros(1); z = zeros(1); p=-1.0*ones(1);
+ dof = (int)-1; cbeta = -1.0*ones(1);
+
+ if(data.Nrows()==design.Nrows()){
+ Matrix dat = remmean(data,1);
+ Matrix pinvdes = pinv(design);
+
+ beta = pinvdes * dat;
+ residu = dat - design*beta;
+ Matrix R = Identity(design.Nrows()) - design * pinvdes;
+ Matrix R2 = R*R;
+ float tR = R.Trace();
+ sigsq = sum(SP(residu,residu))/tR;
+
+ dof = (int)(tR*tR/R2.Trace() - DOFadjust);
+
+ float fact = float(design.Nrows() - 1 - design.Ncols()) / design.Ncols();
+ f_fmf = SP(var(design*beta),pow(var(residu),-1))* fact;
+ pf_fmf = f_fmf.Row(1); pf_fmf &= pf_fmf;
+ for(int ctr1=1;ctr1<=f_fmf.Ncols();ctr1++)
+ pf_fmf(1,ctr1) = 1.0-MISCMATHS::fdtr(design.Ncols(),
+ int(design.Nrows() -1 -design.Ncols()),f_fmf.Column(ctr1).AsScalar());
+ pf_fmf.Row(2) = pf_fmf.Row(1) * pf_fmf.Ncols();
+
+ if(contrasts.Storage()>0 && contrasts.Nrows()==beta.Ncols()){
+ cbeta = contrasts.t()*beta;
+ Matrix tmp = contrasts.t()*pinvdes*pinvdes.t()*contrasts;
+ varcb = diag(tmp)*sigsq;
+ t = SP(cbeta,pow(varcb,-0.5));
+ z = t; p=t;
+ for(int ctr1=1;ctr1<=t.Ncols();ctr1++){
+ ColumnVector tmp = t.Column(ctr1);
+ T2z::ComputeZStats(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
+ z.Column(ctr1) << tmp;
+ T2z::ComputePs(varcb.Column(ctr1),cbeta.Column(ctr1),dof, tmp);
+ p.Column(ctr1) << exp(tmp);
+ }
+ }
+ }
+ return beta;
+ }
+
+
}
diff --git a/melhlprfns.h b/melhlprfns.h
index 6c18ae2f42058c46c8a08d43808c8660c3011c48..58f457568d78f7a2fdd6064eb0580fc805c3cb6e 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -70,7 +70,44 @@ namespace Melodic{
ColumnVector gen_ar(const ColumnVector& in, int maxorder = 1);
Matrix gen_ar(const Matrix& in, int maxorder);
Matrix gen_arCorr(const Matrix& in, int maxorder);
-
+
+ class basicGLM
+ {
+ public:
+
+ //constructor
+ basicGLM(){}
+
+ //destructor
+ ~basicGLM(){}
+
+ Matrix olsfit(const Matrix& data, const Matrix& design,
+ const Matrix& contrasts, int DOFadjust = 0);
+
+ inline Matrix& get_t(){return t;}
+ inline Matrix& get_z(){return z;}
+ inline Matrix& get_p(){return p;}
+ inline Matrix& get_f_fmf(){return f_fmf;}
+ inline Matrix& get_pf_fmf(){return pf_fmf;}
+ inline Matrix& get_cbeta(){return cbeta;}
+ inline Matrix& get_varcb(){return varcb;}
+ inline Matrix& get_sigsq(){return sigsq;}
+ inline Matrix& get_residu(){return residu;}
+ inline int get_dof(){return dof;}
+
+ private:
+ Matrix beta;
+ Matrix residu;
+ Matrix sigsq;
+ Matrix varcb;
+ Matrix cbeta;
+ Matrix f_fmf, pf_fmf;
+ int dof;
+ Matrix t;
+ Matrix z;
+ Matrix p;
+ };
+// Matrix glm_ols(const Matrix& dat, const Matrix& design);
}
#endif
diff --git a/melica.cc b/melica.cc
index f7bf5aa23cafa71260abe16c0571c6f03b419b4a..aa4e181e8124502945f3608e2b2cfba584ba873d 100644
--- a/melica.cc
+++ b/melica.cc
@@ -499,7 +499,6 @@ write_ascii_matrix("Xim",Data.Row(1));
ica_jade(Data);}
if(opts.approach.value()==string("maxent")){
ica_maxent(Data);}
-
if(!no_convergence){//calculate the IC
Matrix temp(melodat.get_unmix()*melodat.get_Data());
@@ -511,18 +510,20 @@ write_ascii_matrix("Xim",Data.Row(1));
scales = stdev(melodat.get_mix());
// cerr << " SCALES 1 " << scales << endl;
- Matrix tmp;
+ Matrix tmp, tmp2;
tmp = SP(melodat.get_mix(),ones(melodat.get_mix().Nrows(),1)*pow(scales,-1));
temp = SP(temp,scales.t()*ones(1,temp.Ncols()));
scales = scales.t();
melodat.set_mix(tmp);
+
melodat.set_IC(temp);
melodat.set_ICstats(scales);
melodat.sort();
melodat.set_TSmode();
+
}
}
}
diff --git a/melodic.cc b/melodic.cc
index 8a9ad7a074439e636e2d95b72f5a16804c13afbb..c0f178658d6c70a97e5cc3fb453698046d63c900 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -56,7 +56,6 @@ void mmonly(Log& logger, MelodicOptions& opts,
int main(int argc, char *argv[])
{
-
try{
// Setup logging:
Log& logger = LogSingleton::getInstance();
@@ -67,16 +66,16 @@ int main(int argc, char *argv[])
//set up data object
MelodicData melodat(opts,logger);
-
+ //set up report object
MelodicReport report(melodat,opts,logger);
if (opts.filtermode || opts.filtermix.value().length()>0){
if(opts.filtermode){ // just filter out some noise from a previous run
- melodat.setup();
- melodat.remove_components();
+ melodat.setup();
+ melodat.remove_components();
}
else
- mmonly(logger,opts,melodat,report);
+ mmonly(logger,opts,melodat,report);
}
else
{ // standard PICA now
@@ -84,105 +83,100 @@ int main(int argc, char *argv[])
bool no_conv;
bool leaveloop = false;
- // melodat.setup2();
melodat.setup();
do{
- //do PCA pre-processing
- MelodicPCA pcaobj(melodat,opts,logger,report);
- pcaobj.perf_pca();
- pcaobj.perf_white();
-
- //do ICA
- MelodicICA icaobj(melodat,opts,logger,report);
- icaobj.perf_ica(melodat.get_white()*melodat.get_Data());
+ //do PCA pre-processing
+ MelodicPCA pcaobj(melodat,opts,logger,report);
+ pcaobj.perf_pca();
+ pcaobj.perf_white();
+
+ //do ICA
+ MelodicICA icaobj(melodat,opts,logger,report);
+ icaobj.perf_ica(melodat.get_white()*melodat.get_Data());
- no_conv = icaobj.no_convergence;
-
- opts.maxNumItt.set_T(500);
- if((opts.approach.value()=="symm")&&
- (retry > std::min(opts.retrystep,3))){
- if(no_conv){
- retry++;
- opts.approach.set_T("defl");
- message(endl << "Restarting MELODIC using deflation approach"
- << endl << endl);
- }
- else{
- leaveloop = true;
- }
- }
- else{
- if(no_conv){
- retry++;
- if(opts.pca_dim.value()-retry*opts.retrystep >
- 0.1*melodat.data_dim()){
- opts.pca_dim.set_T(opts.pca_dim.value()-retry*opts.retrystep);
- }
- else{
- if(opts.pca_dim.value()-retry*opts.retrystep < melodat.data_dim()){
- opts.pca_dim.set_T(opts.pca_dim.value()+retry*opts.retrystep);
- }else{
- leaveloop = true; //stupid, but break does not compile
- //on all platforms
- }
- }
- if(!leaveloop){
+ no_conv = icaobj.no_convergence;
+
+ opts.maxNumItt.set_T(500);
+ if((opts.approach.value()=="symm")&&
+ (retry > std::min(opts.retrystep,3))){
+ if(no_conv){
+ retry++;
+ opts.approach.set_T("defl");
+ message(endl << "Restarting MELODIC using deflation approach"
+ << endl << endl);
+ }
+ else{
+ leaveloop = true;
+ }
+ }
+ else{
+ if(no_conv){
+ retry++;
+ if(opts.pca_dim.value()-retry*opts.retrystep >
+ 0.1*melodat.data_dim()){
+ opts.pca_dim.set_T(opts.pca_dim.value()-retry*opts.retrystep);
+ }
+ else{
+ if(opts.pca_dim.value()-retry*opts.retrystep < melodat.data_dim()){
+ opts.pca_dim.set_T(opts.pca_dim.value()+retry*opts.retrystep);
+ }else{
+ leaveloop = true; //stupid, but break does not compile
+ //on all platforms
+ }
+ }
+ if(!leaveloop){
message(endl << "Restarting MELODIC using -d "
<< opts.pca_dim.value()
<< endl << endl);
}
- }
- }
+ }
+ }
} while (no_conv && retry<opts.maxRestart.value() && !leaveloop);
if(!no_conv){
- //save raw IC results
- melodat.save();
-
- Matrix pmaps;//(melodat.get_IC());
- Matrix mmres;
-
- if(opts.perf_mm.value())
- mmres = mmall(logger,opts,melodat,report,pmaps);
- else{
- if( bool(opts.genreport.value()) ){
- message(endl
- << "Creating web report in " << report.getDir()
- << " " << endl);
- for(int ctr=1; ctr<= melodat.get_IC().Nrows(); ctr++){
- string prefix = "IC_"+num2str(ctr);
- message(" " << ctr);
- report.IC_simplerep(prefix,ctr,melodat.get_IC().Nrows());
- }
-
-
-
- message(endl << endl <<
- " To view the output report point your web browser at " <<
- report.getDir() + "/00index.html" << endl<< endl);
- }
- }
-
- if( bool(opts.genreport.value()) ){
- report.analysistxt();
- report.PPCA_rep();
- }
-
- message("finished!" << endl << endl);
- } else {
- message(endl <<"No convergence -- giving up " << endl << endl);
+ //save raw IC results
+ melodat.save();
+ Matrix pmaps;//(melodat.get_IC());
+ Matrix mmres;
+
+ if(opts.perf_mm.value())
+ mmres = mmall(logger,opts,melodat,report,pmaps);
+ else{
+ if( bool(opts.genreport.value()) ){
+ message(endl
+ << "Creating web report in " << report.getDir()
+ << " " << endl);
+ for(int ctr=1; ctr<= melodat.get_IC().Nrows(); ctr++){
+ string prefix = "IC_"+num2str(ctr);
+ message(" " << ctr);
+ report.IC_simplerep(prefix,ctr,melodat.get_IC().Nrows());
+ }
+
+ message(endl << endl <<
+ " To view the output report point your web browser at " <<
+ report.getDir() + "/00index.html" << endl<< endl);
+ }
+ }
+
+ if( bool(opts.genreport.value()) ){
+ report.analysistxt();
+ report.PPCA_rep();
+ }
+
+ message("finished!" << endl << endl);
+ }
+ else {
+ message(endl <<"No convergence -- giving up " << endl << endl);
}
}
}
- catch(Exception e)
- {
+ catch(Exception e) {
cerr << endl << e.what() << endl;
- }
- catch(X_OptionError& e)
- {
+ }
+ catch(X_OptionError& e) {
cerr << endl << e.what() << endl;
- }
+ }
return 0;
}
@@ -271,7 +265,7 @@ void mmonly(Log& logger, MelodicOptions& opts,
melodat.set_mean(Mean);
melodat.set_IC(ICs);
melodat.set_mix(mixMatrix);
- fmixMatrix = calc_FFT(mixMatrix, opts.logPower.value());
+ fmixMatrix = calc_FFT(melodat.expand_mix(), opts.logPower.value());
melodat.set_fmix(fmixMatrix);
fmixMatrix = pinv(mixMatrix);
melodat.set_unmix(fmixMatrix);
diff --git a/melodic.h b/melodic.h
index 029b2659dd16edaa78be25b5ead3d3c3b3cbe8dc..5669b86039beb6db7c52208361af1468a2ff66df 100644
--- a/melodic.h
+++ b/melodic.h
@@ -34,7 +34,7 @@
namespace Melodic{
-const string version = "3.0 alpha 1";
+const string version = "3.0 beta";
}
diff --git a/meloptions.cc b/meloptions.cc
index 594c216eba0b888d9a4f38e7c3cd2819f835d578..a6ff2a93bc8e3470cd8deae7112cb3bda24d6c08 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -64,11 +64,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
}
if(! options.check_compulsory_arguments())
{
- print_version();
- cout << endl;
- cout << "Usage: melodic <options> -i <filename>" << endl <<
- "Please specify input file name correctly or try melodic --help"
- << endl << endl;
+ print_usage(argc, argv);
exit(2);
}
@@ -173,8 +169,8 @@ MelodicOptions* MelodicOptions::gopt = NULL;
while (!fs.eof()) {
getline(fs,cline);
if(cline.length()>0)
- tmpfnames.push_back(cline);
- }
+ tmpfnames.push_back(cline);
+ }
fs.close();
inputfname.set_T(tmpfnames);
}
@@ -211,7 +207,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
void MelodicOptions::print_usage(int argc, char *argv[])
{
- print_version();
+ print_copyright();
options.usage();
/* cout << "Usage: " << argv[0] << " ";
@@ -221,14 +217,14 @@ void MelodicOptions::print_usage(int argc, char *argv[])
void MelodicOptions::print_version()
{
- cout << "MELODIC Version " << version;
+ cout << endl <<"MELODIC Version " << version;
cout.flush();
}
void MelodicOptions::print_copyright()
{
print_version();
- cout << endl << "Copyright (C) 1999-2002 University of Oxford " << endl;
+ cout << endl << "Copyright (C) 1999-2007 University of Oxford " << endl;
cout.flush();
}
diff --git a/meloptions.h b/meloptions.h
index 52c4eacda8785a6b454f7172e5f6b39057a5135e..d91c18c4969f97a50bd214cab87f9ba85048c78d 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -80,6 +80,13 @@ class MelodicOptions {
Option<float> tr;
Option<bool> logPower;
+ Option<string> fn_Tdesign;
+ Option<string> fn_Tcon;
+ Option<string> fn_TconF;
+ Option<string> fn_Sdesign;
+ Option<string> fn_Scon;
+ Option<string> fn_SconF;
+
Option<bool> output_all;
Option<bool> output_unmix;
Option<bool> output_MMstats;
@@ -103,7 +110,6 @@ class MelodicOptions {
Option<float> nlconst2;
Option<float> smooth_probmap;
-
Option<bool> remove_meanvol;
Option<bool> remove_endslices;
Option<bool> rescale_nht;
@@ -186,7 +192,7 @@ class MelodicOptions {
varnorm(string("--vn,--varnorm"), true,
string("switch off variance normalisation"),
false, no_argument),
- pbsc(string("--pbsc"), true,
+ pbsc(string("--pbsc"), false,
string("switch off conversion to percent BOLD signal change"),
false, no_argument),
pspec(string("--pspec"), false,
@@ -207,7 +213,7 @@ class MelodicOptions {
maxRestart(string("--maxrestart"), 6,
string("maximum number of restarts\n"),
false, requires_argument),
- rank1interval(string("--rank1interval"), 30,
+ rank1interval(string("--rank1interval"), 5,
string("number of iterations between rank-1 approximation (TICA)\n"),
false, requires_argument,false),
mmthresh(string("--mmthresh"), string("0.5"),
@@ -237,6 +243,24 @@ class MelodicOptions {
logPower(string("--logPower"), false,
string("calculate log of power for frequency spectrum\n"),
false, no_argument),
+ fn_Tdesign(string("--Tdes"), string(""),
+ string("design matrix across time-domain"),
+ false, requires_argument),
+ fn_Tcon(string("--Tcon"), string(""),
+ string("t-contrast matrix across time-domain"),
+ false, requires_argument),
+ fn_TconF(string("--Tconf"), string(""),
+ string("F-contrast matrix across time-domain"),
+ false, requires_argument),
+ fn_Sdesign(string("--Sdes"), string(""),
+ string("design matrix across subject-domain"),
+ false, requires_argument),
+ fn_Scon(string("--Scon"), string(""),
+ string("t-contrast matrix across subject-domain"),
+ false, requires_argument),
+ fn_SconF(string("--Sconf"), string(""),
+ string("F-contrast matrix across subject-domain"),
+ false, requires_argument),
output_all(string("--Oall"), false,
string("output everything"),
false, no_argument),
@@ -356,6 +380,12 @@ class MelodicOptions {
options.add(genreport);
options.add(tr);
options.add(logPower);
+ options.add(fn_Tdesign);
+ options.add(fn_Tcon);
+ options.add(fn_TconF);
+ options.add(fn_Sdesign);
+ options.add(fn_Scon);
+ options.add(fn_SconF);
options.add(output_all);
options.add(output_unmix);
options.add(output_MMstats);
diff --git a/melreport.cc b/melreport.cc
index 5315ed151c3c925e1a11ac4c593827f7d2e33e27..7c987696af26ed7a8f41074f60b1b7a8a5fca230 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -1,5 +1,5 @@
/* MELODIC - Multivariate exploratory linear optimized decomposition into
- independent components
+ independent components
melreport.cc - report generation
@@ -13,465 +13,440 @@
#include "utils/log.h"
#include "melreport.h"
#include "melhlprfns.h"
+#include "miscmaths/miscprob.h"
namespace Melodic{
- void MelodicReport::IC_rep(MelGMix &mmodel, int cnum, int dim, Matrix ICstats)
- {
- if( bool(opts.genreport.value()) ){
- addlink(mmodel.get_prefix()+".html",num2str(cnum));
- IChtml.setDir(report.getDir(),mmodel.get_prefix()+".html");
-
- {//start IC page
+ void MelodicReport::IC_rep(MelGMix &mmodel, int cnum, int dim, Matrix ICstats)
+ {
- IChtml << "<HTML> " << endl
- << "<TITLE>MELODIC Component " << num2str(cnum)
- << "</TITLE>" << endl
- << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
- << "/doc/images/fsl-bg.jpg\">" << endl
- << "<hr><CENTER><H1>MELODIC Component " << num2str(cnum)
- << "</H1>"<< endl;
-
- if(cnum>1)
- IChtml << "<a href=\"" << string("IC_")+num2str(cnum-1)
- <<".html\">previous</a> - ";
-
- IChtml << "<a href=\"00index.html\"> index </a>" ;
+ if( bool(opts.genreport.value()) ){
+ addlink(mmodel.get_prefix()+".html",num2str(cnum));
+ IChtml.setDir(report.getDir(),mmodel.get_prefix()+".html");
- if(cnum<dim)
- IChtml << " - <a href=\"" << string("IC_")+num2str(cnum+1)
- <<".html\">next</a><p>";
-
- IChtml << "<hr><p>" << endl;
- }
+ {//start IC page
+ IChtml << "<HTML> " << endl
+ << "<TITLE>MELODIC Component " << num2str(cnum)
+ << "</TITLE>" << endl
+ << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
+ << "/doc/images/fsl-bg.jpg\">" << endl
+ << "<hr><CENTER><H1>MELODIC Component " << num2str(cnum)
+ << "</H1>"<< endl;
+
+ if(cnum>1)
+ IChtml << "<a href=\"" << string("IC_")+num2str(cnum-1)
+ <<".html\">previous</a> - ";
+ IChtml << "<a href=\"00index.html\"> index </a>" ;
+
+ if(cnum<dim)
+ IChtml << " - <a href=\"" << string("IC_")+num2str(cnum+1)
+ <<".html\">next</a><p>";
+
+ IChtml << "<hr><p>" << endl;
+ }
{//output IC stats
- if(ICstats.Storage()>0&&ICstats.Nrows()>=cnum){
- IChtml << ICstats(cnum,1) << " % of explained variance";
- if(ICstats.Ncols()>1)
- IChtml << "; " << ICstats(cnum,2) << " % of total variance";
- if(ICstats.Ncols()>2){
- IChtml << "<p>" <<endl;
- IChtml << " " << ICstats(cnum,3) << " % signal change (pos peak voxel); " << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;}
- IChtml << "<hr><p>" << endl;
- }
+ if(ICstats.Storage()>0&&ICstats.Nrows()>=cnum){
+ IChtml << ICstats(cnum,1) << " % of explained variance";
+ if(ICstats.Ncols()>1)
+ IChtml << "; " << ICstats(cnum,2) << " % of total variance";
+ if(ICstats.Ncols()>2){
+ IChtml << "<p>" <<endl;
+ IChtml << " " << ICstats(cnum,3) << " % signal change (pos peak voxel); " << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;
+ }
+ IChtml << "<hr><p>" << endl;
+ }
}
-
+
volume4D<float> tempVol;
if(mmodel.get_threshmaps().Storage()>0&&
- (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols()))
- {//Output thresholded IC map
-
-
- tempVol.setmatrix(mmodel.get_threshmaps().Row(1),melodat.get_mask());
- volume<float> map1;
- volume<float> map2;
- map1 = threshold(tempVol[0],float(0.0), tempVol[0].max());
- map2 = threshold(tempVol[0],tempVol[0].min(), float(0.0));
- // save_volume(map,report.appendDir(mmodel.get_prefix()+"_thresh")
- // ,melodat.tempInfo);
- volume<float> newvol;
- // for(int ctr=1; ctr<500; ctr++){
-// cerr << ctr << " ";
- miscpic newpic;
-
- float map1min = std::max((map1 + binarise(tempVol[0],tempVol[0].min(),
- float(0.0)) * map1.max()).min(),float(0.01));
- float map1max = std::max(map1.max(),float(0.01));
- float map2min = std::min(map2.min(),float(-0.01));
- float map2max = std::min((map2 + binarise(tempVol[0],float(0.0),
- tempVol[0].max()) * map2.min()).max(),float(-0.01));
+ (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols()))
+ {//Output thresholded IC map
+ tempVol.setmatrix(mmodel.get_threshmaps().Row(1),melodat.get_mask());
+ volume<float> map1;
+ volume<float> map2;
+ map1 = threshold(tempVol[0],float(0.0), tempVol[0].max());
+ map2 = threshold(tempVol[0],tempVol[0].min(), float(0.0));
+ volume<float> newvol;
- newpic.overlay(newvol, melodat.get_mean(), map1, map2,
- melodat.get_mean().percentile(0.02),
- melodat.get_mean().percentile(0.98),
- map1min, map1max, map2max, map2min,
- 0, 0, &melodat.tempInfo);
+ miscpic newpic;
- char instr[10000];
+ float map1min = std::max((map1 + binarise(tempVol[0],tempVol[0].min(),
+ float(0.0)) * map1.max()).min(),float(0.01));
+ float map1max = std::max(map1.max(),float(0.01));
+ float map2min = std::min(map2.min(),float(-0.01));
+ float map2max = std::min((map2 + binarise(tempVol[0],float(0.0),
+ tempVol[0].max()) * map2.min()).max(),float(-0.01));
- //save_volume(newvol,report.appendDir(mmodel.get_prefix()+"rendered"),
- // melodat.tempInfo);
- sprintf(instr," ");
- strcat(instr," -s 2");
- strcat(instr," -A 950 ");
- strcat(instr,string(report.appendDir(mmodel.get_prefix()+
- "_thresh.png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - thresholded IC map ")+
- mmodel.get_infstr(1));
- newpic.set_cbar(string("ysb"));
- //cerr << instr << endl;
-
-
+ newpic.overlay(newvol, melodat.get_mean(), map1, map2,
+ melodat.get_mean().percentile(0.02),
+ melodat.get_mean().percentile(0.98),
+ map1min, map1max, map2max, map2min,
+ 0, 0, &melodat.tempInfo);
+ char instr[10000];
- if(map1.max()-map1.min()>0.01)
- newpic.slicer(newvol, instr, &melodat.tempInfo);
- else
- newpic.slicer(map1, instr, &melodat.tempInfo);
-
- }
- // }
-// exit(2);
-
- IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
- IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
- +"_thresh.png\" ALT=\"MMfit\"></A><p>" << endl;
+ //save_volume(newvol,report.appendDir(mmodel.get_prefix()+"rendered"),
+ // melodat.tempInfo);
+ sprintf(instr," ");
+ strcat(instr," -s 2");
+ strcat(instr," -A 950 ");
+ strcat(instr,string(report.appendDir(mmodel.get_prefix()+
+ "_thresh.png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - thresholded IC map ") + mmodel.get_infstr(1));
+ newpic.set_cbar(string("ysb"));
+ //cerr << instr << endl;
+
+ if(map1.max()-map1.min()>0.01)
+ newpic.slicer(newvol, instr, &melodat.tempInfo);
+ else
+ newpic.slicer(map1, instr, &melodat.tempInfo);
+ IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
+ IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
+ +"_thresh.png\" ALT=\"MMfit\"></A><p>" << endl;
+ }
-
{//plot time course
- IChtml << "<H3> Temporal mode <p>" << endl <<endl;
- miscplot newplot;
-
- if(opts.tr.value()>0.0)
- newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
- report.appendDir(string("t")+num2str(cnum)+".png"),
- string("Timecourse (in seconds); TR = ")+
- float2str(opts.tr.value(),0,2,0)+" s",
- opts.tr.value(),150,4,1);
- else
- newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
- report.appendDir(string("t")+num2str(cnum)+".png"),
- string("Timecourse (in TRs)"));
- write_ascii_matrix(report.appendDir(string("t")
- +num2str(cnum)+".txt"),
- melodat.get_Tmodes(cnum-1));
- IChtml << "<A HREF=\"" << string("t")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
+ IChtml << "<H3> Temporal mode <p>" << endl <<endl;
+ miscplot newplot;
+ Matrix tmptc = melodat.get_Tmodes(cnum-1).t();
+
+ //add GLM OLS fit
+ /* basicGLM glm;
+ if(melodat.Tdes.Storage()){
+ Matrix betas = glm.olsfit(tmptc.t(),melodat.Tdes,melodat.Tcon).t();
+ tmptc &= betas*melodat.Tdes.t();
+ newplot.add_label(string("IC ")+num2str(cnum)+" time course");
+ newplot.add_label("full model fit");
+
+cerr << endl << endl <<
+glm.get_f_fmf() << endl<<
+glm.get_pf_fmf() << endl << endl;
+ }
+*/
+ if(opts.tr.value()>0.0)
+ newplot.timeseries(tmptc,
+ report.appendDir(string("t")+num2str(cnum)+".png"),
+ string("Timecourse (in seconds); TR = ")+
+ float2str(opts.tr.value(),0,2,0)+" s",
+ opts.tr.value(),150,4,1);
+ else
+ newplot.timeseries(tmptc,
+ report.appendDir(string("t")+num2str(cnum)+".png"),
+ string("Timecourse (in TRs)"));
+ write_ascii_matrix(report.appendDir(string("t")
+ +num2str(cnum)+".txt"),tmptc.t());
+ IChtml << "<A HREF=\"" << string("t")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
}//time series plot
-
- {//plot frequency
- miscplot newplot;
- RowVector empty(1);
- empty = 0.0;
- int fact = int(std::pow(10.0,int(std::log10(float(melodat.get_Tmodes(0).Nrows())))));
- if(opts.tr.value()>0.0)
- newplot.timeseries(empty | melodat.get_fmix().Column(cnum).t(),
- report.appendDir(string("f")+
- num2str(cnum)+".png"),
- string("FFT of timecourse (in Hz / ") +num2str(fact)+")",
- fact/(opts.tr.value()*melodat.get_Tmodes(0).Nrows()), 150,
- 0,2);
- else
- newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
- report.appendDir(string("f")+
- num2str(cnum)+".png"),
- string(string("FFT of timecourse (in cycles); ")
- +"frequency(Hz)=cycles/("
- +num2str(melodat.get_Tmodes(0).Nrows())
- +"* TR); period(s)=("
- +num2str(melodat.get_Tmodes(0).Nrows())
- +"* TR)/cycles"));
- write_ascii_matrix(report.appendDir(string("f")
- +num2str(cnum)+".txt"),
- melodat.get_fmix().Column(cnum));
- IChtml << "<A HREF=\"" << string("f")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
+ {//plot frequency
+ miscplot newplot;
+ RowVector empty(1);
+ empty = 0.0;
+ int fact = int(std::pow(10.0,int(std::log10(float(melodat.get_Tmodes(0).Nrows())))));
+ if(opts.logPower.value())
+ newplot.add_ylabel(string("log-Power"));
+ else
+ newplot.add_ylabel(string("Power"));
+
+ Matrix fmixtc = calc_FFT(melodat.get_Tmodes(cnum-1), opts.logPower.value());
+
+ if(opts.tr.value()>0.0){
+ newplot.add_xlabel(string("Frequency (in Hz / ")+num2str(fact)+ " )");
+ newplot.timeseries(empty | fmixtc.t(),
+ report.appendDir(string("f")+
+ num2str(cnum)+".png"),
+ string("Powerspectrum of timecourse"),
+ fact/(opts.tr.value()*melodat.get_Tmodes(0).Nrows()), 150,0,2);
+ }else{
+ newplot.add_xlabel(string("Frequency (in cycles); ")
+ +"frequency(Hz)=cycles/("
+ +num2str(melodat.get_Tmodes(0).Nrows())
+ +"* TR); period(s)=("
+ +num2str(melodat.get_Tmodes(0).Nrows())
+ +"* TR)/cycles"
+ );
+ newplot.timeseries(fmixtc.t(),
+ report.appendDir(string("f")+num2str(cnum)+".png"),
+ string("Powerspectrum of timecourse"));
+ }
+ write_ascii_matrix(report.appendDir(string("f")
+ +num2str(cnum)+".txt"), fmixtc);
+ IChtml << "<A HREF=\"" << string("f")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
}//frequency plot
if(cnum <= (int)melodat.get_Smodes().size())
- {//plot subject mode
- Matrix smode;
- smode = melodat.get_Smodes(cnum-1);
- if(smode.Nrows() > 1){
- miscplot newplot;
- newplot.timeseries(smode.t(),
- report.appendDir(string("s")+num2str(cnum)+".png"),
- string("Subject/Session mode"));
- write_ascii_matrix(report.appendDir(string("s")
- +num2str(cnum)+".txt"),
- smode);
- IChtml << "<A HREF=\"" << string("s")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("s")+num2str(cnum)+".png\"></A><p>" << endl;
- }
- }//subject mode plot
-
+ {//plot subject mode
+ Matrix smode;
+ smode = melodat.get_Smodes(cnum-1);
+ if(smode.Nrows() > 1){
+ miscplot newplot;
+ newplot.setscatter(smode,5);
+ newplot.timeseries(smode.t(),
+ report.appendDir(string("s")+num2str(cnum)+".png"),
+ string("Subject/Session mode"));
+ newplot.boxplot(smode,
+ report.appendDir(string("b")+num2str(cnum)+".png"),
+ string("Subject/Session mode"));
+ write_ascii_matrix(report.appendDir(string("s")
+ +num2str(cnum)+".txt"), smode);
+ IChtml << "<A HREF=\"" << string("s")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("s")+num2str(cnum)+".png\"></A>" << endl;
+ IChtml << "<A HREF=\"" << string("s")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("b")+num2str(cnum)+".png\"></A><p>" << endl;
+ }
+ }//subject mode plot
+
if(mmodel.get_threshmaps().Storage()>0&&
- (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&&
- (mmodel.get_threshmaps().Nrows()>1))
- {//Output other thresholded IC map
+ (mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&&
+ (mmodel.get_threshmaps().Nrows()>1))
+ {//Output other thresholded IC map
- for(int tctr=2; tctr<=mmodel.get_threshmaps().Nrows(); tctr++){
- tempVol.setmatrix(mmodel.get_threshmaps().Row(tctr),melodat.get_mask());
- volume<float> map1;
- volume<float> map2;
- map1 = threshold(tempVol[0],float(0.0), tempVol[0].max());
- map2 = threshold(tempVol[0],tempVol[0].min(), float(0.0));
- // save_volume(map,report.appendDir(mmodel.get_prefix()+"_thresh")
- // ,melodat.tempInfo);
- volume<float> newvol;
- miscpic newpic;
-
- float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
- float(0.0)) * map1.max()).min();
- float map1max = map1.max();
- float map2min = map2.min();
- float map2max = (map2 + binarise(tempVol[0],float(0.0),
- tempVol[0].max()) * map2.min()).max();
+ for(int tctr=2; tctr<=mmodel.get_threshmaps().Nrows(); tctr++){
+ tempVol.setmatrix(mmodel.get_threshmaps().Row(tctr),melodat.get_mask());
+ volume<float> map1;
+ volume<float> map2;
+ map1 = threshold(tempVol[0],float(0.0), tempVol[0].max());
+ map2 = threshold(tempVol[0],tempVol[0].min(), float(0.0));
+ // save_volume(map,report.appendDir(mmodel.get_prefix()+"_thresh")
+ // ,melodat.tempInfo);
+ volume<float> newvol;
+ miscpic newpic;
- //cerr << endl << map1min << " " << map1max << endl
- // << map2min << " " << map2max << endl;
-
- // if(map1.max()-map1.min()>0.01)
- newpic.overlay(newvol, melodat.get_mean(), map1, map2,
- melodat.get_mean().percentile(0.02),
- melodat.get_mean().percentile(0.98),
- map1min, map1max, map2max, map2min,
- 0, 0, &melodat.tempInfo);
+ float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
+ float(0.0)) * map1.max()).min();
+ float map1max = map1.max();
+ float map2min = map2.min();
+ float map2max = (map2 + binarise(tempVol[0],float(0.0),
+ tempVol[0].max()) * map2.min()).max();
-
- char instr[10000];
+ //cerr << endl << map1min << " " << map1max << endl
+ // << map2min << " " << map2max << endl;
- //save_volume(newvol,report.appendDir(mmodel.get_prefix()+"rendered"),
- // melodat.tempInfo);
- sprintf(instr," ");
- strcat(instr," -s 2");
- strcat(instr," -A 950 ");
- strcat(instr,string(report.appendDir(mmodel.get_prefix()+"_thresh"+
- num2str(tctr)+".png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - thresholded IC map ("+
- num2str(tctr)+") ")+
- mmodel.get_infstr(tctr));
- newpic.set_cbar(string("ysb"));
- //cerr << instr << endl;
- newpic.slicer(newvol, instr, &melodat.tempInfo);
-
- IC_rep_det(mmodel, cnum, dim);
- IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
- IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
- +"_thresh"+num2str(tctr)+".png\" ALT=\"MMfit\"></A><p>" << endl;
- }
- }
-
-
+ // if(map1.max()-map1.min()>0.01)
+ newpic.overlay(newvol, melodat.get_mean(), map1, map2,
+ melodat.get_mean().percentile(0.02),
+ melodat.get_mean().percentile(0.98),
+ map1min, map1max, map2max, map2min,
+ 0, 0, &melodat.tempInfo);
+
+ char instr[10000];
+
+ //save_volume(newvol,report.appendDir(mmodel.get_prefix()+"rendered"),
+ // melodat.tempInfo);
+ sprintf(instr," ");
+ strcat(instr," -s 2");
+ strcat(instr," -A 950 ");
+ strcat(instr,string(report.appendDir(mmodel.get_prefix()+"_thresh"+
+ num2str(tctr)+".png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - thresholded IC map ("+
+ num2str(tctr)+") ")+
+ mmodel.get_infstr(tctr));
+ newpic.set_cbar(string("ysb"));
+ //cerr << instr << endl;
+ newpic.slicer(newvol, instr, &melodat.tempInfo);
+
+ IC_rep_det(mmodel, cnum, dim);
+ IChtml << "<a href=\""+mmodel.get_prefix()+"_MM.html\">";
+ IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
+ +"_thresh"+num2str(tctr)+".png\" ALT=\"MMfit\"></A><p>" << endl;
+ }
+ }
+
{ //finish IC page
- IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
- << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
- << "FMRIB Software Library</A>.</FONT>" << endl
- << "</BODY></HTML>" << endl;
+ IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
+ << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "FMRIB Software Library</A>.</FONT>" << endl
+ << "</BODY></HTML>" << endl;
} //finish IC page
IC_rep_det(mmodel, cnum, dim);
- }
- }
+ }
+ }
void MelodicReport::IC_rep_det(MelGMix &mmodel, int cnum, int dim)
{
if( bool(opts.genreport.value()) ){
{//start IC2 page
- IChtml2.setDir(report.getDir(),mmodel.get_prefix()+"_MM.html");
- IChtml2 << "<HTML> " << endl
- << "<TITLE>Component " << num2str(cnum)
- << " Mixture Model fit </TITLE>" << endl
- << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
- << "/doc/images/fsl-bg.jpg\">" << endl
- << "<hr><CENTER><H1>Component " << num2str(cnum)
- << " Mixture Model fit </H1>"<< endl;
+ IChtml2.setDir(report.getDir(),mmodel.get_prefix()+"_MM.html");
+ IChtml2 << "<HTML> " << endl
+ << "<TITLE>Component " << num2str(cnum)
+ << " Mixture Model fit </TITLE>" << endl
+ << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
+ << "/doc/images/fsl-bg.jpg\">" << endl
+ << "<hr><CENTER><H1>Component " << num2str(cnum)
+ << " Mixture Model fit </H1>"<< endl;
- if(cnum>1)
- IChtml2 << "<a href=\"" << string("IC_")+num2str(cnum-1)
- <<"_MM.html\">previous</a> - ";
-
- IChtml2 << "<a href=\""+ mmodel.get_prefix() +
- ".html\"> up to IC report </a>";
-
- if(cnum<dim)
- IChtml2 << " - <a href=\"" << string("IC_")+num2str(cnum+1)
- <<"_MM.html\">next</a><p>";
- IChtml2 << "<hr><p>" << endl;
+ if(cnum>1)
+ IChtml2 << "<a href=\"" << string("IC_")+num2str(cnum-1)
+ <<"_MM.html\">previous</a> - ";
+ IChtml2 << "<a href=\""+ mmodel.get_prefix() +
+ ".html\"> up to IC report </a>";
+ if(cnum<dim)
+ IChtml2 << " - <a href=\"" << string("IC_")+num2str(cnum+1)
+ <<"_MM.html\">next</a><p>";
+ IChtml2 << "<hr><p>" << endl;
}
volume4D<float> tempVol;
if(melodat.get_IC().Storage()>0)
- {//Output raw IC map
-
- // tempVol.setmatrix(melodat.get_IC().Row(cnum),
- // melodat.get_mask());
- tempVol.setmatrix(mmodel.get_data(),
- melodat.get_mask());
- volume<float> map1;
- volume<float> map2;
- map1 = threshold(tempVol[0],float(0.0),
- tempVol[0].max());
- map2 = threshold(tempVol[0],tempVol[0].min(),
- float(-0.0));
-
- volume<float> newvol;
- miscpic newpic;
-
- // float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
- // float(0.0)) * map1.max()).robustmin();
- float map1max = map1.percentile(0.99);
- float map2min = map2.percentile(0.01);
- //float map2max = (map2 + binarise(tempVol[0],float(0.0),
- // tempVol[0].max()) * map2.min()).robustmax();
+ {//Output raw IC map
+
+ // tempVol.setmatrix(melodat.get_IC().Row(cnum),
+ // melodat.get_mask());
+ tempVol.setmatrix(mmodel.get_data(),
+ melodat.get_mask());
+ volume<float> map1;
+ volume<float> map2;
+ map1 = threshold(tempVol[0],float(0.0),
+ tempVol[0].max());
+ map2 = threshold(tempVol[0],tempVol[0].min(),
+ float(-0.0));
+
+ volume<float> newvol;
+ miscpic newpic;
+
+ // float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
+ // float(0.0)) * map1.max()).robustmin();
+ float map1max = map1.percentile(0.99);
+ float map2min = map2.percentile(0.01);
+ //float map2max = (map2 + binarise(tempVol[0],float(0.0),
+ // tempVol[0].max()) * map2.min()).robustmax();
- newpic.overlay(newvol, melodat.get_mean(), map1, map2,
- float(0.0),
- float(0.0),
- float(0.01), map1max, float(-0.01), map2min,
- 0, 0, &melodat.tempInfo);
+ newpic.overlay(newvol, melodat.get_mean(), map1, map2,
+ float(0.0),
+ float(0.0),
+ float(0.01), map1max, float(-0.01), map2min,
+ 0, 0, &melodat.tempInfo);
- char instr[10000];
+ char instr[10000];
- //save_volume(newvol,report.appendDir(mmodel.get_prefix()+"rendered"),
- // melodat.tempInfo);
- sprintf(instr," ");
- strcat(instr," -s 2");
- strcat(instr," -A 950 ");
- strcat(instr,string(report.appendDir(mmodel.get_prefix()+
- ".png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - raw Z transformed IC map (1 - 99 percentile)"));
- newpic.set_cbar(string("ysb"));
-
- newpic.slicer(newvol, instr, &melodat.tempInfo);
- }
+ //save_volume(newvol,report.appendDir(mmodel.get_prefix()+"rendered"),
+ // melodat.tempInfo);
+ sprintf(instr," ");
+ strcat(instr," -s 2");
+ strcat(instr," -A 950 ");
+ strcat(instr,string(report.appendDir(mmodel.get_prefix()+
+ ".png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - raw Z transformed IC map (1 - 99 percentile)"));
+ newpic.set_cbar(string("ysb"));
+
+ newpic.slicer(newvol, instr, &melodat.tempInfo);
+ }
IChtml2 << "<a href=\""+mmodel.get_prefix()+".html\">";
IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
- ".png\"><A><p>" << endl;
-
-
+ ".png\"><A><p>" << endl;
if(mmodel.get_probmap().Storage()>0&&
- (mmodel.get_probmap().Ncols() == mmodel.get_data().Ncols())&&
- (mmodel.get_probmap().Nrows() == mmodel.get_data().Nrows()))
- {//Output probmap
- tempVol.setmatrix(mmodel.get_probmap(),melodat.get_mask());
+ (mmodel.get_probmap().Ncols() == mmodel.get_data().Ncols())&&
+ (mmodel.get_probmap().Nrows() == mmodel.get_data().Nrows()))
+ {//Output probmap
+ tempVol.setmatrix(mmodel.get_probmap(),melodat.get_mask());
- volume<float> map;
- map = tempVol[0];
+ volume<float> map;
+ map = tempVol[0];
- volume<float> newvol;
- miscpic newpic;
-
- newpic.overlay(newvol, melodat.get_mean(), map, map,
- melodat.get_mean().percentile(0.02),
- melodat.get_mean().percentile(0.98),
- float(0.1), float(1.0), float(0.0), float(0.0),
- 0, 0, &melodat.tempInfo);
-
- //newpic.set_minmax(float(0.0),float(0.0),float(0.0),
- // float(1.0),float(0.0),float(0.0));
-
- char instr[10000];
-
- sprintf(instr," ");
- strcat(instr,"-l render1 -s 2");
- strcat(instr," -A 950 ");
- strcat(instr,string(report.appendDir(mmodel.get_prefix()+
- "_prob.png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - Mixture Model probability map"));
+ volume<float> newvol;
+ miscpic newpic;
+
+ newpic.overlay(newvol, melodat.get_mean(), map, map,
+ melodat.get_mean().percentile(0.02),
+ melodat.get_mean().percentile(0.98),
+ float(0.1), float(1.0), float(0.0), float(0.0),
+ 0, 0, &melodat.tempInfo);
+
+ char instr[10000];
+
+ sprintf(instr," ");
+ strcat(instr,"-l render1 -s 2");
+ strcat(instr," -A 950 ");
+ strcat(instr,string(report.appendDir(mmodel.get_prefix()+
+ "_prob.png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - Mixture Model probability map"));
- newpic.set_cbar(string("y"));
- newpic.slicer(newvol, instr, &melodat.tempInfo);
-
-// {
-// Log MMdetails;
-// MMdetails.setDir(report.getDir(),mmodel.get_prefix()+"_MM.txt");
-// MMdetails << mmodel.get_prefix() << " Mixture Model fit "
-// << endl << endl;
-// MMdetails << " Means : " << mmodel.get_means() << endl
-// << " Vars : " << mmodel.get_vars() << endl
-// << " Prop. : " << mmodel.get_pi() << endl;
-// // if(mmodel.get_type()=="GGM"){
-// // MMdetails << endl << " Gamma offset : "
-// // << mmodel.get_offset() << endl;
-// // }
-// }
- // IChtml2 << "<A HREF=\"" +mmodel.get_prefix()+"_MM.txt\"> ";
- IChtml2 << "<a href=\""+mmodel.get_prefix()+".html\">";
- IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
- "_prob.png\">" << endl;
- // IChtml2 << "</A>";
- IChtml2 << "</A><p>" << endl;
- }
+ newpic.set_cbar(string("y"));
+ newpic.slicer(newvol, instr, &melodat.tempInfo);
+
+ IChtml2 << "<a href=\""+mmodel.get_prefix()+".html\">";
+ IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
+ "_prob.png\">" << endl;
+
+ IChtml2 << "</A><p>" << endl;
+ }
{//Output GGM/GMM fit
- miscplot newplot;
-
-// cerr << endl << endl;
-// cerr << mmodel.get_means() << endl;
-// cerr << mmodel.get_vars() << endl;
-// cerr << mmodel.get_pi() << endl;
-
- if(mmodel.get_type()=="GGM"){
- newplot.add_label("IC map histogram");
- newplot.add_label("full GGM fit");
- newplot.add_label("background Gaussian");
- newplot.add_label("Gamma distributions");
- newplot.gmmfit(mmodel.get_data().Row(1),
- mmodel.get_means(),
- mmodel.get_vars(),
- mmodel.get_pi(),
- report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
- string(mmodel.get_prefix() +
- " GGM("+num2str(mmodel.mixtures())+") fit"),
- true, float(0.0), float(2.0));
-
- // newplot.histogram(mmodel.get_data().Row(1),
- // report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
- // string(mmodel.get_prefix() +
- // " GGM("+num2str(mmodel.mixtures())+") fit"));
-
- //, mmodel.get_offset());
- }
- else{
- newplot.add_label("IC map histogram");
- newplot.add_label("full GMM fit");
- newplot.add_label("individual Gaussians");
- newplot.gmmfit(mmodel.get_data().Row(1),
- mmodel.get_means(),
- mmodel.get_vars(),
- mmodel.get_pi(),
- report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
- string(mmodel.get_prefix() +
- " GMM("+num2str(mmodel.mixtures())+") fit"),
- false, float(0.0), float(2.0));
+ miscplot newplot;
+
+ if(mmodel.get_type()=="GGM"){
+ newplot.add_label("IC map histogram");
+ newplot.add_label("full GGM fit");
+ newplot.add_label("background Gaussian");
+ newplot.add_label("Gamma distributions");
+ newplot.gmmfit(mmodel.get_data().Row(1),
+ mmodel.get_means(),
+ mmodel.get_vars(),
+ mmodel.get_pi(),
+ report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
+ string(mmodel.get_prefix() +
+ " GGM("+num2str(mmodel.mixtures())+") fit"),
+ true, float(0.0), float(2.0));
- }
- // cerr << "finish HTML2 " << endl;
- IChtml2 << "<A HREF=\"" +mmodel.get_prefix()+"_MMfit_detail.png\"> ";
- IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
- "_MMfit.png\"></A><p>" << endl;
+ }
+ else{
+ newplot.add_label("IC map histogram");
+ newplot.add_label("full GMM fit");
+ newplot.add_label("individual Gaussians");
+ newplot.gmmfit(mmodel.get_data().Row(1),
+ mmodel.get_means(),
+ mmodel.get_vars(),
+ mmodel.get_pi(),
+ report.appendDir(mmodel.get_prefix()+"_MMfit.png"),
+ string(mmodel.get_prefix() +
+ " GMM("+num2str(mmodel.mixtures())+") fit"),
+ false, float(0.0), float(2.0));
+
+ }
+
+ // IChtml2 << "<A HREF=\"" +mmodel.get_prefix()+"_MMfit_detail.png\"> ";
+ IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
+ "_MMfit.png\"><p>" << endl;
} //GGM/GMM plot
- {
- IChtml2 << "<br> " << mmodel.get_prefix()
- << " Mixture Model fit <br>" << endl
- << "<br> Means : " << mmodel.get_means() << endl
- << "<br> Vars : " << mmodel.get_vars() << endl
- << "<br> Prop. : " << mmodel.get_pi() << endl;
- // if(mmodel.get_type()=="GGM"){
- // IChtml2 << "<br> Gamma offset : "
- // << mmodel.get_offset() << "<br><p>"<< endl;
- // }
- }
+ {//MM parameters
+ IChtml2 << "<br> " << mmodel.get_prefix()
+ << " Mixture Model fit <br>" << endl
+ << "<br> Means : " << mmodel.get_means() << endl
+ << "<br> Vars : " << mmodel.get_vars() << endl
+ << "<br> Prop. : " << mmodel.get_pi() << endl;
+ }
{ //finish IC2 page
- IChtml2<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
- << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
- << "FMRIB Software Library</A>.</FONT>" << endl
- << "</BODY></HTML>" << endl;
+ IChtml2<< "<HR><FONT SIZE=1>This page produced automatically by "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
+ << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "FMRIB Software Library</A>.</FONT>" << endl
+ << "</BODY></HTML>" << endl;
} //finish IC2 page
}
}
-
void MelodicReport::IC_simplerep(string prefix, int cnum, int dim)
{
if( bool(opts.genreport.value()) ){
@@ -480,140 +455,138 @@ namespace Melodic{
{//start IC page
- IChtml << "<HTML> " << endl
- << "<TITLE>MELODIC Component " << num2str(cnum)
- << "</TITLE>" << endl
- << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
- << "/doc/images/fsl-bg.jpg\">" << endl
- << "<hr><CENTER><H1>MELODIC Component " << num2str(cnum)
- << "</H1>"<< endl;
+ IChtml << "<HTML> " << endl
+ << "<TITLE>MELODIC Component " << num2str(cnum)
+ << "</TITLE>" << endl
+ << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
+ << "/doc/images/fsl-bg.jpg\">" << endl
+ << "<hr><CENTER><H1>MELODIC Component " << num2str(cnum)
+ << "</H1>"<< endl;
- if(cnum>1)
- IChtml << "<a href=\"" << string("IC_")+num2str(cnum-1)
- <<".html\">previous</a> - ";
+ if(cnum>1)
+ IChtml << "<a href=\"" << string("IC_")+num2str(cnum-1)
+ <<".html\">previous</a> - ";
- IChtml << "<a href=\"00index.html\"> index </a>" ;
+ IChtml << "<a href=\"00index.html\"> index </a>" ;
- if(cnum<dim)
- IChtml << " - <a href=\"" << string("IC_")+num2str(cnum+1)
- <<".html\">next</a><p>";
+ if(cnum<dim)
+ IChtml << " - <a href=\"" << string("IC_")+num2str(cnum+1)
+ <<".html\">next</a><p>";
- IChtml << "<hr><p>" << endl;
+ IChtml << "<hr><p>" << endl;
}
-
volume4D<float> tempVol;
if(melodat.get_IC().Storage()>0)
- {//Output raw IC map
-
- tempVol.setmatrix(melodat.get_IC().Row(cnum),
- melodat.get_mask());
- volume<float> map1;
- volume<float> map2;
- map1 = threshold(tempVol[0],float(0.0),
- tempVol[0].max());
- map2 = threshold(tempVol[0],tempVol[0].min(),
- float(-0.0));
-
- volume<float> newvol;
- miscpic newpic;
-
- // float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
- // float(0.0)) * map1.max()).robustmin();
- float map1max = map1.percentile(0.99);
- float map2min = map2.percentile(0.01);
- //float map2max = (map2 + binarise(tempVol[0],float(0.0),
- // tempVol[0].max()) * map2.min()).robustmax();
+ {//Output raw IC map
+
+ tempVol.setmatrix(melodat.get_IC().Row(cnum),
+ melodat.get_mask());
+ volume<float> map1;
+ volume<float> map2;
+ map1 = threshold(tempVol[0],float(0.0),
+ tempVol[0].max());
+ map2 = threshold(tempVol[0],tempVol[0].min(),
+ float(-0.0));
+
+ volume<float> newvol;
+ miscpic newpic;
+
+ // float map1min = (map1 + binarise(tempVol[0],tempVol[0].min(),
+ // float(0.0)) * map1.max()).robustmin();
+ float map1max = map1.percentile(0.99);
+ float map2min = map2.percentile(0.01);
+ //float map2max = (map2 + binarise(tempVol[0],float(0.0),
+ // tempVol[0].max()) * map2.min()).robustmax();
- newpic.overlay(newvol, melodat.get_mean(), map1, map2,
- float(0.0),
- float(0.0),
- float(0.01), map1max, float(-0.01), map2min,
- 0, 0, &melodat.tempInfo);
+ newpic.overlay(newvol, melodat.get_mean(), map1, map2,
+ float(0.0),
+ float(0.0),
+ float(0.01), map1max, float(-0.01), map2min,
+ 0, 0, &melodat.tempInfo);
- char instr[10000];
+ char instr[10000];
- //save_volume(newvol,report.appendDir(prefix+"rendered"),
- // melodat.tempInfo);
- sprintf(instr," ");
- strcat(instr," -s 2");
- strcat(instr," -A 950 ");
- strcat(instr,string(report.appendDir(prefix+
- ".png")).c_str());
- newpic.set_title(string("Component No. "+num2str(cnum)+
- " - raw Z transformed IC map (1 - 99 percentile)"));
- newpic.set_cbar(string("ysb"));
+ //save_volume(newvol,report.appendDir(prefix+"rendered"),
+ // melodat.tempInfo);
+ sprintf(instr," ");
+ strcat(instr," -s 2");
+ strcat(instr," -A 950 ");
+ strcat(instr,string(report.appendDir(prefix+
+ ".png")).c_str());
+ newpic.set_title(string("Component No. "+num2str(cnum)+
+ " - raw Z transformed IC map (1 - 99 percentile)"));
+ newpic.set_cbar(string("ysb"));
- newpic.slicer(newvol, instr, &melodat.tempInfo);
- }
+ newpic.slicer(newvol, instr, &melodat.tempInfo);
+ }
IChtml << "<img BORDER=0 SRC=\""+ prefix+
- ".png\"><p>" << endl;
+ ".png\"><p>" << endl;
{//plot time course
- miscplot newplot;
+ miscplot newplot;
- if(opts.tr.value()>0.0)
- newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
- report.appendDir(string("t")+
- num2str(cnum)+".png"),
- string("Timecourse (in seconds); TR = ")+
- float2str(opts.tr.value(),0,2,0)+" s",
- opts.tr.value(),150,4,1);
- else
- newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
- report.appendDir(string("t")+
- num2str(cnum)+".png"),
- string("Timecourse (in TRs)"));
- write_ascii_matrix(report.appendDir(string("t")
- +num2str(cnum)+".txt"),
- melodat.get_Tmodes(cnum-1));
- IChtml << "<A HREF=\"" << string("t")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
+ if(opts.tr.value()>0.0)
+ newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
+ report.appendDir(string("t")+
+ num2str(cnum)+".png"),
+ string("Timecourse (in seconds); TR = ")+
+ float2str(opts.tr.value(),0,2,0)+" s",
+ opts.tr.value(),150,4,1);
+ else
+ newplot.timeseries(melodat.get_Tmodes(cnum-1).t(),
+ report.appendDir(string("t")+
+ num2str(cnum)+".png"),
+ string("Timecourse (in TRs)"));
+ write_ascii_matrix(report.appendDir(string("t")
+ +num2str(cnum)+".txt"),
+ melodat.get_Tmodes(cnum-1));
+ IChtml << "<A HREF=\"" << string("t")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
}//time series plot
{//plot frequency
- miscplot newplot;
- int fact = int(std::pow(10.0,
- int(std::log10(float(melodat.get_Tmodes(0).Nrows())))));
+ miscplot newplot;
+ int fact = int(std::pow(10.0,
+ int(std::log10(float(melodat.get_Tmodes(0).Nrows())))));
- if(opts.tr.value()>0.0)
- newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
- report.appendDir(string("f")+
- num2str(cnum)+".png"),
- string("FFT of timecourse (in Hz / ") +
- num2str(fact)+")",
- fact/(opts.tr.value()*melodat.get_Tmodes(0).Nrows()),
- 150,
- 0,2);
- else
- newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
- report.appendDir(string("f")+
- num2str(cnum)+".png"),
- string(string("FFT of timecourse (in cycles); ")
+ if(opts.tr.value()>0.0)
+ newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
+ report.appendDir(string("f")+
+ num2str(cnum)+".png"),
+ string("FFT of timecourse (in Hz / ") +
+ num2str(fact)+")",
+ fact/(opts.tr.value()*melodat.get_Tmodes(0).Nrows()),
+ 150,0,2);
+ else
+ newplot.timeseries(melodat.get_fmix().Column(cnum).t(),
+ report.appendDir(string("f")+
+ num2str(cnum)+".png"),
+ string(string("FFT of timecourse (in cycles); ")
+"frequency(Hz)=cycles/("
+num2str(melodat.get_Tmodes(0).Nrows())
+"* TR); period(s)=("
+num2str(melodat.get_Tmodes(0).Nrows())
+"* TR)/cycles"));
- write_ascii_matrix(report.appendDir(string("f")
- +num2str(cnum)+".txt"),
- melodat.get_Tmodes(cnum-1));
- IChtml << "<A HREF=\"" << string("f")
- +num2str(cnum)+".txt" << "\"> ";
- IChtml << "<img BORDER=0 SRC=\""
- +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
+ write_ascii_matrix(report.appendDir(string("f")
+ +num2str(cnum)+".txt"),
+ melodat.get_Tmodes(cnum-1));
+ IChtml << "<A HREF=\"" << string("f")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
}//frequency plot
{ //finish IC page
- IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
- << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
- << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
- << "FMRIB Software Library</A>.</FONT>" << endl
- << "</BODY></HTML>" << endl;
+ IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
+ << "<A HREF=\"http://www.fmrib.ox.ac.uk/fsl/melodic/index.html\"> MELODIC </A>"
+ << " - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - "
+ << "FMRIB Software Library</A>.</FONT>" << endl
+ << "</BODY></HTML>" << endl;
} //finish IC page
}
}
@@ -633,15 +606,15 @@ namespace Melodic{
newplot.add_label("% of expl. variance");
if(melodat.get_PPCA().Storage()>0){
- what = what.Columns(1,melodat.get_PPCA().Nrows());
- what &= melodat.get_PPCA().t();
- newplot.add_label("dim. estimate");
+ what = what.Columns(1,melodat.get_PPCA().Nrows());
+ what &= melodat.get_PPCA().t();
+ newplot.add_label("dim. estimate");
}
newplot.timeseries(what,
- report.appendDir("EVplot.png"),
- string("Eigenspectrum Analysis"),
- 0,450,4,0);
+ report.appendDir("EVplot.png"),
+ string("Eigenspectrum Analysis"),
+ 0,450,4,0);
report << "<img BORDER=0 SRC=\"EVplot.png\"><p>" << endl;
}//time series plot
diff --git a/test.cc b/test.cc
index 91ce76f37f0fc1a8c9e0253f23319b8c909107cb..1c2553466686da4e07e1afa019ecc9089949b98c 100644
--- a/test.cc
+++ b/test.cc
@@ -37,18 +37,23 @@ int do_work(int argc, char* argv[])
{
Matrix mat;
- mat = normrnd(300,(int)num.value())+2;
cout << "BLAH " << num.value() << endl;
- for (int i=1; i <= (int)num.value();i++){
+ mat=normrnd(300,1);
+ miscplot::Timeseries(mat.t(),string("test0.png"),string("TEST"));
+
+ for (int i=1; i <= (int)num.value();i++){
cout << "Processing " << i << endl;
- miscplot newplot;
- ColumnVector col;
- col = mat.Column(i);
- newplot.add_bpdata(col);
- newplot.boxplot(string("test")+num2str(i)+string(".png"),string("TEST"));
+ miscplot newplot;
+ newplot.GDCglobals_set();
+ mat = normrnd(300,3)+2;
+ Matrix col;
+ col = mat;
+ newplot.add_bpdata(col);
+ // newplot.add_bpdata(col);
+ newplot.boxplot(string("test")+num2str(i)+string(".png"),string("TEST"));
}
- return 0;
+ return 0;
}
////////////////////////////////////////////////////////////////////////////