diff --git a/meldata.cc b/meldata.cc
index 89fffd3a959852ad18fb01150032253f828d48f3..ef4421f16f1af7a740a6135999134d2a7532a994 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -203,7 +203,7 @@ namespace Melodic{
calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
}
- if(numfiles < 2){
+ if(numfiles == 1){
Data = alldat;
Matrix tmp = Identity(Data.Nrows());
DWM.push_back(tmp);
@@ -246,20 +246,18 @@ namespace Melodic{
//initialize Mean
read_volume(Mean,opts.inputfname.value().at(0),tempInfo);
- //save first image
- tmpnam(Mean_fname); // generate a tmp name
- save_volume(Mean,Mean_fname);
-
//create mask
create_mask(Mask);
- // clean /tmp
- char callRMstr[1000];
- ostrstream osc(callRMstr,1000);
- osc << "rm " << string(Mean_fname) <<"* " << '\0';
-
- system(callRMstr);
-
+ //setup background image
+ if(opts.bgimage.value()>""){
+ read_volume(background,opts.bgimage.value());
+ if(!samesize(Mean,background)){
+ cerr << "ERROR:: background image and data have different dimensions \n\n";
+ exit(2);
+ }
+ }
+
if(!samesize(Mean,Mask)){
cerr << "ERROR:: mask and data have different dimensions \n\n";
exit(2);
@@ -633,6 +631,10 @@ namespace Melodic{
else{
if(opts.perf_bet.value() && opts.use_mask.value()){ //use BET
message("Create mask ... ");
+ //save first image
+ tmpnam(Mean_fname); // generate a tmp name
+ save_volume(Mean,Mean_fname);
+
// set up all strings
string BET_outputfname = string(Mean_fname)+"_brain";
@@ -652,6 +654,14 @@ namespace Melodic{
// read back the Mask file
read_volume(theMask,Mask_fname);
+
+ // clean /tmp
+ char callRMstr[1000];
+ ostrstream osc(callRMstr,1000);
+ osc << "rm " << string(Mean_fname) <<"* " << '\0';
+
+ system(callRMstr);
+
message("done" << endl);
}
else{
diff --git a/meldata.h b/meldata.h
index 4ebec202466067c58cb8462a830c382093f51240..623db9234bbc0c1ff8a2703f0dd2bdff115fcf1c 100644
--- a/meldata.h
+++ b/meldata.h
@@ -133,6 +133,14 @@ namespace Melodic{
inline volume<float>& get_mean() {return Mean;}
inline void set_mean(volume<float>& Arg) {Mean = Arg;}
+ inline volume<float>& get_bg() {
+ if(opts.bgimage.value()>"")
+ return background;
+ else
+ return Mean;
+ }
+ inline void set_bg(volume<float>& Arg) {background = Arg;}
+
inline Matrix& get_Data() {return Data;}
inline void set_Data(Matrix& Arg) {Data = Arg;}
@@ -143,13 +151,16 @@ namespace Melodic{
inline void set_ICstats(Matrix& Arg) {ICstats = Arg;}
inline Matrix& get_EVP() {return EVP;}
- inline void set_EVP(Matrix& Arg) {EVP = Arg;}
+ inline void set_EVP(Matrix& Arg) {if(EVP.Storage()==0)
+ EVP = Arg;}
inline Matrix& get_EV() {return EV;}
- inline void set_EV(Matrix& Arg) {EV = Arg;}
+ inline void set_EV(Matrix& Arg) {if(EV.Storage()==0)
+ EV = Arg;}
inline Matrix& get_PPCA() {return PPCA;}
- inline void set_PPCA(Matrix& Arg) {PPCA = Arg;}
+ inline void set_PPCA(Matrix& Arg) {if(PPCA.Storage()==0)
+ PPCA = Arg;}
inline Matrix& get_stdNoisei() {return stdNoisei;}
inline void set_stdNoisei(Matrix& Arg) {stdNoisei = Arg;}
@@ -211,6 +222,8 @@ namespace Melodic{
Matrix stdNoisei;
volume<float> Mask;
volume<float> Mean;
+ volume<float> background;
+
Matrix Data;
Matrix PPCA;
Matrix jointCC;
diff --git a/melica.cc b/melica.cc
index ee4fb3f16f2048c82ac3b46d0d43b49692be7f26..1379340ef85398412b09fde334bc47847316f45c 100644
--- a/melica.cc
+++ b/melica.cc
@@ -115,7 +115,7 @@ namespace Melodic{
outMsize(" redUMM ", redUMM);
}
} while((itt_ctr2 < 3 || itt_ctr >= opts.maxNumItt.value()) &&
- (opts.approach.value() == string("tica")) && cum_itt < 25*opts.maxNumItt.value());
+ (opts.approach.value() == string("tica")) && cum_itt < 100*opts.maxNumItt.value());
if(itt_ctr>=opts.maxNumItt.value()){
cerr << " No convergence after " << cum_itt <<" steps "<<endl;
diff --git a/meloptions.h b/meloptions.h
index 81d7bc5894274f2009f1a03184c93f4db8bf8372..3ca7a3fd8f4024ca38a2b5f7e358a3987d0fe763 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -75,6 +75,7 @@ class MelodicOptions {
Option<string> filter;
Option<bool> genreport;
+ Option<string> bgimage;
Option<float> tr;
Option<bool> logPower;
@@ -172,7 +173,7 @@ class MelodicOptions {
pca_est(string("--dimest"), string("lap"),
string("use specific dim. estimation technique: lap, bic, mdl, aic, mean (default: lap)"),
false, requires_argument),
- joined_whiten(string("--separate_whiten"), true,
+ joined_whiten(string("--sep_whiten"), true,
string("switch on separate whitening"),
false, no_argument),
numICs(string("-n,--numICs"), -1,
@@ -188,10 +189,10 @@ class MelodicOptions {
string("switch off variance normalisation"),
false, no_argument),
pbsc(string("--pbsc"), false,
- string("switch off conversion to percent BOLD signal change"),
+ string(" switch off conversion to percent BOLD signal change"),
false, no_argument),
pspec(string("--pspec"), false,
- string("switch on conversion to powerspectra"),
+ string(" switch on conversion to powerspectra"),
false, no_argument),
segment(string("--covarweight"), string(""),
string("voxel-wise weights for the covariance matrix (e.g. segmentation information)"),
@@ -226,12 +227,15 @@ class MelodicOptions {
smodename(string("--smode"), string(""),
string("\tmatrix of session modes for report generation"),
false, requires_argument),
- filter(string("-f,--filter"), string(""),
- string("component numbers to remove\n "),
- false, requires_argument),
+ filter(string("-f,--filter"), string(""),
+ string("component numbers to remove\n "),
+ false, requires_argument),
genreport(string("--report"), false,
string("generate Melodic web report"),
false, no_argument),
+ bgimage(string("--bgimage"), string(""),
+ string("specify background image for report (default: mean image)\n "),
+ false, requires_argument),
tr(string("--tr"), 0.0,
string("\tTR in seconds"),
false, requires_argument),
@@ -239,25 +243,25 @@ class MelodicOptions {
string("calculate log of power for frequency spectrum\n"),
false, no_argument),
fn_Tdesign(string("--Tdes"), string(""),
- string("design matrix across time-domain"),
+ string(" design matrix across time-domain"),
false, requires_argument),
fn_Tcon(string("--Tcon"), string(""),
- string("t-contrast matrix across time-domain"),
+ string(" t-contrast matrix across time-domain"),
false, requires_argument),
fn_TconF(string("--Tconf"), string(""),
- string("F-contrast matrix across time-domain"),
+ string(" F-contrast matrix across time-domain"),
false, requires_argument),
fn_Sdesign(string("--Sdes"), string(""),
- string("design matrix across subject-domain"),
+ string(" design matrix across subject-domain"),
false, requires_argument),
fn_Scon(string("--Scon"), string(""),
- string("t-contrast matrix across subject-domain"),
+ 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),
+ string(" F-contrast matrix across subject-domain"),
+ false, requires_argument,false),
output_all(string("--Oall"), false,
- string("output everything"),
+ string(" output everything"),
false, no_argument),
output_unmix(string("--Ounmix"), false,
string("output unmixing matrix"),
@@ -290,7 +294,7 @@ class MelodicOptions {
string("prints this help message"),
false, no_argument),
debug(string("--debug"), false,
- string("switch on debug messages"),
+ string(" switch on debug messages"),
false, no_argument),
guessfname(string("-g,--guess"), string(""),
string("file name of guess of mixing matrix"),
@@ -332,16 +336,16 @@ class MelodicOptions {
options(string(""),
string(" melodic -i <filename> <options>")+
string("\n \t \t to run melodic")+
- string("\n melodic -i <filename> --mix=melodic_mix")+
- string(" --filter=\"string of component numbers\"")+
- string("\n \t \t to remove estimated ICs from input")+
+// string("\n melodic -i <filename> --mix=melodic_mix")+
+// string(" --filter=\"string of component numbers\"")+
+// string("\n \t \t to remove estimated ICs from input")+
string("\n melodic -i <filename> --ICs=melodic_IC")+
string(" --mix=melodic_mix <options>")+
string("\n \t \t to run Mixture Model based inference on estimated ICs")+
string("\n melodic --help "))
{
try {
- options.add(logdir);
+ options.add(logdir);
options.add(inputfname);
options.add(outputfname);
options.add(guessfname);
@@ -372,6 +376,7 @@ class MelodicOptions {
options.add(smodename);
options.add(filter);
options.add(genreport);
+ options.add(bgimage);
options.add(tr);
options.add(logPower);
options.add(fn_Tdesign);
diff --git a/melreport.cc b/melreport.cc
index 3597b780ec713e5a59d319fcbb2f446d31b81bbf..a61db4082fbb9790c3f2496e312addbd6209ad61 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -1,6 +1,6 @@
/* MELODIC - Multivariate exploratory linear optimized decomposition into
independent components
-
+ melodat.get_bg()
melreport.cc - report generation
Christian F. Beckmann, FMRIB Image Analysis Group
@@ -79,9 +79,9 @@ namespace Melodic{
float map2max = std::min((map2 + binarise(tempVol[0],float(0.0),
tempVol[0].max()) * map2.min()).max(),float(-0.01));
- newpic.overlay(newvol, melodat.get_mean(), map1, map2,
- melodat.get_mean().percentile(0.02),
- melodat.get_mean().percentile(0.98),
+ newpic.overlay(newvol, melodat.get_bg(), map1, map2,
+ melodat.get_bg().percentile(0.02),
+ melodat.get_bg().percentile(0.98),
map1min, map1max, map2max, map2min,
0, 0, &melodat.tempInfo);
char instr[10000];
@@ -320,9 +320,9 @@ namespace Melodic{
// << 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),
+ newpic.overlay(newvol, melodat.get_bg(), map1, map2,
+ melodat.get_bg().percentile(0.02),
+ melodat.get_bg().percentile(0.98),
map1min, map1max, map2max, map2min,
0, 0, &melodat.tempInfo);
@@ -411,7 +411,7 @@ namespace Melodic{
//float map2max = (map2 + binarise(tempVol[0],float(0.0),
// tempVol[0].max()) * map2.min()).robustmax();
- newpic.overlay(newvol, melodat.get_mean(), map1, map2,
+ newpic.overlay(newvol, melodat.get_bg(), map1, map2,
float(0.0),
float(0.0),
float(0.01), map1max, float(-0.01), map2min,
@@ -448,9 +448,9 @@ namespace Melodic{
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),
+ newpic.overlay(newvol, melodat.get_bg(), map, map,
+ melodat.get_bg().percentile(0.02),
+ melodat.get_bg().percentile(0.98),
float(0.1), float(1.0), float(0.0), float(0.0),
0, 0, &melodat.tempInfo);
@@ -582,7 +582,7 @@ namespace Melodic{
//float map2max = (map2 + binarise(tempVol[0],float(0.0),
// tempVol[0].max()) * map2.min()).robustmax();
- newpic.overlay(newvol, melodat.get_mean(), map1, map2,
+ newpic.overlay(newvol, melodat.get_bg(), map1, map2,
float(0.0),
float(0.0),
float(0.01), map1max, float(-0.01), map2min,
diff --git a/melreport.h b/melreport.h
index 87567f24c5f5ca07baee35ab458a4f94b5257c71..68662e91bd1e8f8331cb3299d17e18de66e82167 100644
--- a/melreport.h
+++ b/melreport.h
@@ -42,7 +42,7 @@ namespace Melodic{
logger(plogger){
if( bool(opts.genreport.value()) ){
const time_t tmptime = time(NULL);
-
+
report.makeDir(logger.appendDir("report"),"00index.html");
report << "<HTML>" << endl << endl
<< "<TITLE>MELODIC Report</TITLE>"<< endl <<endl
diff --git a/test.cc b/test.cc
index 7fc130b4c9c7643fbb2f0f30358711685c5e3dc1..d0e01658e7abd71331c42045f78a26a7d469f374 100644
--- a/test.cc
+++ b/test.cc
@@ -81,6 +81,14 @@ int do_work(int argc, char* argv[]) {
cerr << weights << endl <<endl << Rows << endl;
cerr << SP(weights,pow(Rows,-1)) << endl;
+
+ vector<int> vi;
+ for(int ctr=1;ctr<=5;ctr++)
+ vi.push_back(ctr);
+
+ for(int ctr=0;ctr<5;ctr++)
+ cerr << " vi["<<ctr<<"] = " <<vi.at(ctr)<<endl;
+
return 0;
}