diff --git a/meldata.cc b/meldata.cc
index f7457755f4424cce50342c16f1ebe6bacc028b9a..582434954a3ae25f3f7c6fc819ce2f93f7a43660 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -86,7 +86,10 @@ namespace Melodic{
//variance - normalisation
if(opts.varnorm.value()){
message(" Normalising by voxel-wise variance ...");
- stdDev = varnorm(tmpData,std::min(30,tmpData.Nrows()-1),2.3);
+ if(stdDev.Storage()==0)
+ stdDev = varnorm(tmpData,std::min(30,tmpData.Nrows()-1),2.3)/numfiles;
+ else
+ stdDev += varnorm(tmpData,std::min(30,tmpData.Nrows()-1),2.3)/numfiles;
stdDevi = pow(stdDev,-1);
message(" done" << endl);
}
@@ -133,8 +136,6 @@ namespace Melodic{
{
Matrix tmp, tmpT, tmpS, tmpT2, tmpS2;
tmp = expand_dimred(mixMatrix);
- saveascii(tmp,"exp");
- saveascii(mixMatrix,"mix");
tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
tmpS = zeros(numfiles, tmp.Ncols());
krfact(tmp,tmpT,tmpS);
@@ -173,8 +174,11 @@ namespace Melodic{
opts.approach.set_T("tica");
Matrix alldat, tmpData;
-// bool tmpvarnorm = opts.varnorm.value();
-// opts.varnorm.set_T(false);
+ bool tmpvarnorm = opts.varnorm.value();
+
+ if(opts.joined_vn.value()){
+ opts.varnorm.set_T(false);
+ }
alldat = process_file(opts.inputfname.value().at(0), numfiles) / numfiles;
if(opts.debug.value())
save4D(alldat,string("preproc_dat") + num2str(1));
@@ -184,23 +188,22 @@ namespace Melodic{
save4D(tmpData /numfiles,string("preproc_dat") + num2str(ctr+1));
alldat += tmpData / numfiles;
}
-
- // opts.varnorm.set_T(tmpvarnorm);
- //variance - normalisation
- // if(opts.varnorm.value()){
- // message(" Normalising by voxel-wise variance ...");
- // stdDev = varnorm(alldat,alldat.Nrows(),3.1);
- // stdDevi = pow(stdDev,-1);
- // message(" done" << endl);
- // }
//update mask
if(opts.update_mask.value()){
message("Excluding voxels with constant value ...");
- // update_mask(Mask, alldat);
+ update_mask(Mask, alldat);
message(" done" << endl);
}
+ if(opts.joined_vn.value()){
+ //variance - normalisation
+ message(endl<<"Normalising 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);
@@ -242,11 +245,13 @@ namespace Melodic{
WM.push_back(tmp);
}
else {
+
for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
-// if(opts.varnorm.value())
-// varnorm(tmpData,stdDev);
-
+
+ if(opts.joined_vn.value()){
+ tmpData=SP(tmpData,pow(ones(tmpData.Nrows(),1)*stdDev,-1));
+ }
// whiten (separate / joint)
Matrix newWM,newDWM;
if(!opts.joined_whiten.value()){
@@ -270,9 +275,11 @@ namespace Melodic{
Data &= tmpData;
}
}
+ opts.varnorm.set_T(tmpvarnorm);
- message(" Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl);
-
+ message(endl << " Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
+ outMsize("stdDev",stdDev);
+
meanC=mean(Data,2);
if(opts.debug.value())
save4D(Data,"concat_data");
@@ -298,6 +305,8 @@ namespace Melodic{
cerr << "ERROR:: background image and data have different dimensions \n\n";
exit(2);
}
+ }else{
+ background = Mean;
}
if(!samesize(Mean,Mask)){
@@ -373,6 +382,10 @@ namespace Melodic{
// ICadjust = IC;
}
else{
+ Matrix resids = stdev(Data - mixMatrix * IC);
+ for(int ctr=1;ctr<=resids.Ncols();ctr++)
+ if(resids(1,ctr) < 0.05)
+ resids(1,ctr)=1;
stdNoisei = pow(stdev(Data - mixMatrix * IC)*
std::sqrt((float)(Data.Nrows()-1))/
std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
@@ -598,7 +611,7 @@ namespace Melodic{
}
}
return count;
- }
+ }
float MelodicData::est_resels(volume4D<float> R, volume<float> mask)
{
@@ -712,7 +725,6 @@ namespace Melodic{
char callRMstr[1000];
ostrstream osc(callRMstr,1000);
osc << "rm " << string(Mean_fname) <<"* " << '\0';
-
system(callRMstr);
message("done" << endl);
@@ -728,8 +740,8 @@ namespace Melodic{
message("done" << endl);
}
else{ //well, don't threshold then
- theMask = Mean;
- theMask = 1.0;
+ theMask = Mean;
+ theMask = 1.0;
}
}
}
@@ -764,7 +776,7 @@ namespace Melodic{
Matrix allmodes = Smodes.at(0);
for(int ctr = 1; ctr < (int)Smodes.size();++ctr)
allmodes |= Smodes.at(ctr);
- tmpscales = mean(allmodes).t();
+ tmpscales = median(allmodes).t();
} else {
// re-order wrt standard deviation of IC maps
tmpscales = stdev(IC,2);
@@ -803,7 +815,7 @@ namespace Melodic{
Matrix copeP(tmpscales), copeN(tmpscales);
Matrix max_ICs(tmpscales), min_ICs(tmpscales);
-
+
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++){
int i,j;
max_ICs(ctr_i,1) = IC.Row(ctr_i).Maximum2(i,j);
@@ -821,7 +833,7 @@ namespace Melodic{
ICstats |= copeP;
ICstats |= copeN;
}
-
+
mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
unmixMatrix = pinv(mixMatrix);
}
diff --git a/melgmix.cc b/melgmix.cc
index bbe13dcf0bdc8bf6870e8851c6b07f9516edd14d..3831a749ee56929d528d94a5b5890002776ba243 100644
--- a/melgmix.cc
+++ b/melgmix.cc
@@ -59,6 +59,8 @@ namespace Melodic{
datastdev= stdev(dat,2).AsScalar();
data=(dat - datamean)/datastdev;
+ dbgmsg(" mapdat; mean: " << datamean << " std: " <<datastdev << endl);
+
props=zeros(1,nummix);
vars=zeros(1,nummix);
means=zeros(1,nummix);
@@ -87,6 +89,7 @@ namespace Melodic{
if(epsilon >=0 && epsilon < 0.0000001)
{epsilon = std::log(float(dat.Ncols()))/1000 ;}
fixdim = fixit;
+ dbgmsg(" parameters; " << means << " : " << vars << " : " << props << endl);
}
Matrix MelGMix::threshold(const RowVector& dat, string levels){
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 78780ea30dc3572a11818150d5e465d90b161d84..95666b57a75c397d1981b62a0f40a0d773b40b06 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -121,10 +121,10 @@ namespace Melodic{
return out;
} //RowVector varnorm
- void varnorm(Matrix& in, RowVector& vars)
+ void varnorm(Matrix& in, const RowVector& vars)
{
- Matrix tmp;
- in = SP(in,pow(ones(in.Nrows(),1) * vars,-1));
+ Matrix tmp = vars;
+ in = SP(in,pow(ones(in.Nrows(),1) * tmp,-1));
}
RowVector varnorm(Matrix& in, Matrix& Corr, int dim, float level)
@@ -169,6 +169,14 @@ namespace Melodic{
return Res;
} //Matrix SP
+ Matrix corrcoef(const Matrix& in1, const Matrix& in2){
+ Matrix tmp = in1;
+ tmp |= in2;
+ Matrix out;
+ out=MISCMATHS::corrcoef(tmp,0);
+ return out.SubMatrix(1,in1.Ncols(),in1.Ncols()+1,out.Ncols());
+ }
+
Matrix calc_corr(const Matrix& in, bool econ)
{
Matrix Res;
diff --git a/melhlprfns.h b/melhlprfns.h
index eb1327e01f3477ada2c7fe7c16da2ab3128dc42e..4d95a223a35f8b149e3af491f9458bf70186f253 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -29,7 +29,7 @@ namespace Melodic{
Matrix convert_to_pbsc(Matrix& Mat);
RowVector varnorm(Matrix& in, int dim = 30, float level = 1.6);
- void varnorm(Matrix& in, RowVector& vars);
+ void varnorm(Matrix& in, const RowVector& vars);
RowVector varnorm(Matrix& in, Matrix& Corr, int dim = 30, float level = 1.6);
Matrix SP2(const Matrix& in, const Matrix& weights, bool econ = 0);
@@ -37,6 +37,7 @@ namespace Melodic{
RowVector Feta(int n1,int n2);
RowVector cumsum(const RowVector& Inp);
+ Matrix corrcoef(const Matrix& in1, const Matrix& in2);
Matrix calc_corr(const Matrix& in, bool econ = 0);
Matrix calc_corr(const Matrix& in, const Matrix& weights, bool econ = 0);
diff --git a/melica.cc b/melica.cc
index 17a3c958a2da252ca37af5c07413a657303a4dda..d264066e075e05059b8aa421d1e3727909c5348d 100644
--- a/melica.cc
+++ b/melica.cc
@@ -32,11 +32,10 @@ namespace Melodic{
// see www.cis.hut.fi/projects/ica/fastica/
//initialize matrices
- Matrix redUMM_old, redUMM_veryold;
+ Matrix redUMM_old, rank1_old;
Matrix tmpU;
//srand((unsigned int)timer(NULL));
-
redUMM = melodat.get_white()*
unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
@@ -51,12 +50,16 @@ namespace Melodic{
symm_orth(redUMM);
int itt_ctr,itt_ctr2=1,cum_itt=0,newmaxitts = opts.maxNumItt.value();
- double minAbsSin = 1.0;
+ double minAbsSin = 1.0, minAbsSin2 = 1.0;
if(opts.approach.value() == string("tica"))
opts.maxNumItt.set_T(opts.rank1interval.value());
+
+ rank1_old = melodat.get_dewhite()*redUMM;
+ rank1_old = melodat.expand_dimred(rank1_old);
+ rank1_old = krapprox(rank1_old,int(rank1_old.Nrows()/melodat.get_numfiles()));
+
do{// TICA loop
itt_ctr = 1;
- redUMM_veryold = redUMM;
do{ // da loop!!!
redUMM_old = redUMM;
//calculate IC estimates
@@ -94,37 +97,38 @@ namespace Melodic{
symm_orth(redUMM);
if(opts.approach.value() == string("tica")){
- message(" TICA:");
+ message("");
}
//termination condition : angle between old and new < epsilon
minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_old)).Minimum();
message(" Step no. " << cum_itt + itt_ctr << " change : " << minAbsSin << endl);
- if(abs(minAbsSin) < opts.epsilon.value()){ break;}
+ if((abs(minAbsSin) < opts.epsilon.value())&&
+ (opts.approach.value()!=string("tica"))){ break;}
itt_ctr++;
} while(itt_ctr < opts.maxNumItt.value());
cum_itt += itt_ctr;
itt_ctr2++;
if(opts.approach.value() == string("tica")){
- message(" Rank-1 approximation of the time courses; " <<endl;);
+ message(" Rank-1 approximation of the time courses; ");
Matrix temp(melodat.get_dewhite() * redUMM);
- outMsize(" 2nd unmmixing matrix ", temp);
temp = melodat.expand_dimred(temp);
- outMsize(" 2nd unmmixing matrix ", temp);
- temp = krapprox(temp,int(temp.Nrows()/melodat.get_numfiles()));
- outMsize(" 2nd unmmixing matrix ", temp);
- temp = melodat.reduce_dimred(temp);
- outMsize(" 2nd unmmixing matrix ", temp);
- redUMM = melodat.get_white() * temp;
- outMsize(" redUMM ", redUMM);
- minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_veryold)).Minimum();
- message(" Step no. " << cum_itt + itt_ctr << " change : " << minAbsSin << endl);
- if(abs(minAbsSin) < opts.epsilon.value()){ break;}
+ temp = krapprox(temp,int(temp.Nrows()/melodat.get_numfiles()));
+ minAbsSin2 = 1 - diag(corrcoef(temp,rank1_old)).Minimum();
+ rank1_old = temp;
+ temp = melodat.reduce_dimred(temp);
+ redUMM = melodat.get_white() * temp;
+
+ message(" change : " << abs(minAbsSin2) << endl);
+ if(abs(minAbsSin2) < 0.01 && abs(minAbsSin) < opts.epsilon.value()){ break;}
}
- } while((itt_ctr2 < 3 || itt_ctr >= opts.maxNumItt.value()) &&
- (opts.approach.value() == string("tica")) && cum_itt < newmaxitts);
+ } while(
+ (itt_ctr2 < newmaxitts/opts.maxNumItt.value()) &&
+ (opts.approach.value() == string("tica")) &&
+ cum_itt < newmaxitts);
- if(itt_ctr>=opts.maxNumItt.value()){
+ if((itt_ctr>=opts.maxNumItt.value() && (opts.approach.value()!=string("tica")))
+ || (cum_itt >= newmaxitts && opts.approach.value()==string("tica"))){
cerr << " No convergence after " << cum_itt <<" steps "<<endl;
} else {
message(" Convergence after " << cum_itt <<" steps " << endl << endl);
@@ -141,7 +145,7 @@ namespace Melodic{
if(!opts.explicitnums || opts.numICs.value()>dim){
opts.numICs.set_T(dim);
message(" Using numICs:" << opts.numICs.value() << endl);
- }
+ }
//redUMM = zeros(dim); // got to start somewhere
redUMM = melodat.get_white()*
@@ -153,7 +157,7 @@ namespace Melodic{
message(" Use columns in " << opts.guessfname.value() << " as initial values for the mixing matrix " <<endl);
guess = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
guesses = guess.Ncols();
- }
+ }
int ctrIC = 1;
int numRestart = 0;
@@ -205,12 +209,12 @@ namespace Melodic{
if((norm2(w-w_old) < 0.001*opts.epsilon.value())&&(itt_ctr>10) ||
(norm2(w+w_old) < 0.001*opts.epsilon.value())&&(itt_ctr>10)){
break;
- }
+ }
//cout << norm2(w-w_old) << " " << norm2(w+w_old) << endl;
itt_ctr++;
- } while(itt_ctr <= opts.maxNumItt.value());
+ } while(itt_ctr <= opts.maxNumItt.value());
- if(itt_ctr<opts.maxNumItt.value()){
+ if(itt_ctr<opts.maxNumItt.value()){
redUMM.Column(ctrIC) = w;
message(" estimated using " << itt_ctr << " iterations " << endl);
ctrIC++;
@@ -221,8 +225,7 @@ namespace Melodic{
<< numRestart << " attempts "
<< " giving up " << endl);
break;
- }
- else{
+ }else{
numRestart++;
message(endl <<" Estimation failed - restart "
<< numRestart << endl);
@@ -266,7 +269,7 @@ namespace Melodic{
double minAbsSin = 1.0;
Matrix Id;
Id = Identity(redUMM.Ncols());
- cerr << " nonlinearity : " << opts.nonlinearity.value() << endl;
+ //cerr << " nonlinearity : " << opts.nonlinearity.value() << endl;
do{ // da loop!!!
redUMM_old = redUMM;
@@ -454,7 +457,7 @@ namespace Melodic{
Matrix scales;
scales = stdev(melodat.get_mix());
- // cerr << " SCALES 1 " << scales << endl;
+ //cerr << " SCALES 1 " << scales << endl;
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()));
@@ -466,6 +469,7 @@ namespace Melodic{
melodat.set_IC(temp);
melodat.set_ICstats(scales);
melodat.sort();
+
message("Calculating T- and S-modes " << endl);
melodat.set_TSmode();
diff --git a/melodic.cc b/melodic.cc
index 4837a28277fe85341a67561310caa7e85a4c101b..63d3132dff8fb9260708536c958903bd4981e0c1 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -302,19 +302,20 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
message(" IC map " << ctr << " ... "<< endl;);
Matrix ICmap;
+ dbgmsg(" stdNoisei max : "<< melodat.get_stdNoisei().Maximum() <<" "<< melodat.get_stdNoisei().Minimum() << endl);
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);
-
+ }else{
+ ICmap = melodat.get_IC().Row(ctr);
+ }
string wherelog;
if(opts.genreport.value())
wherelog = report.getDir();
else
wherelog = logger.getDir();
+ dbgmsg(" ICmap max : "<< mean(ICmap,2).AsScalar() << endl);
mixmod.setup( ICmap, melodat.tempInfo,
wherelog,ctr,
melodat.get_mask(),
diff --git a/melodic.h b/melodic.h
index 949a344cceb7ec0b8b2deab2f2812693f00a1bd4..16d8516d4eb01b07d7abd9b53f16043e37d60f4e 100644
--- a/melodic.h
+++ b/melodic.h
@@ -26,6 +26,12 @@
cout.flush(); \
}
+#define dbgmsg(msg) { \
+ MelodicOptions&opt = MelodicOptions::getInstance(); \
+ if(opt.debug.value()) {\
+ cout << msg; } \
+}
+
#define outMsize(msg,Mat) { \
MelodicOptions& opt = MelodicOptions::getInstance(); \
if(opt.debug.value()) \
diff --git a/meloptions.h b/meloptions.h
index 35f4e65306b3e835ef600382e0e01175c6d5a78d..a6256aed0b081b3c83eac3725311e05cb11f380b 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -53,6 +53,7 @@ class MelodicOptions {
Option<int> pca_dim;
Option<string> pca_est;
Option<bool> joined_whiten;
+ Option<bool> joined_vn;
Option<int> numICs;
Option<string> approach;
Option<string> nonlinearity;
@@ -177,6 +178,9 @@ class MelodicOptions {
joined_whiten(string("--sep_whiten"), true,
string("switch on separate whitening"),
false, no_argument),
+ joined_vn(string("--joined_vn"), false,
+ string("switch on joined variance nomalisation"),
+ false, no_argument, false),
numICs(string("-n,--numICs"), -1,
string("numer of IC's to extract (for deflation approach)"),
false, requires_argument),
@@ -201,7 +205,7 @@ class MelodicOptions {
tsmooth(string("--spca"), false,
string("smooth the eigenvectors prior to IC decomposition"),
false, no_argument, false),
- epsilon(string("--eps,--epsilon"), 0.00005,
+ epsilon(string("--eps,--epsilon"), 0.0005,
string("minimum error change"),
false, requires_argument),
maxNumItt(string("--maxit"), 500,
@@ -210,7 +214,7 @@ class MelodicOptions {
maxRestart(string("--maxrestart"), 6,
string("maximum number of restarts\n"),
false, requires_argument),
- rank1interval(string("--rank1interval"), 5,
+ rank1interval(string("--rank1interval"), 10,
string("number of iterations between rank-1 approximation (TICA)\n"),
false, requires_argument,false),
mmthresh(string("--mmthresh"), string("0.5"),
@@ -349,6 +353,7 @@ class MelodicOptions {
options.add(pca_dim);
options.add(pca_est);
options.add(joined_whiten);
+ options.add(joined_vn);
options.add(numICs);
options.add(approach);
options.add(nonlinearity);
diff --git a/melreport.cc b/melreport.cc
index 501e8152fdd019175ccf5ba136665d24c110606d..272598e5ce088275c133439c44fa3dac237e80ab 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -43,13 +43,12 @@ namespace Melodic{
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){
+ if(ICstats.Ncols()>2&&melodat.get_numfiles()==1){
IChtml << "<p>" <<endl;
IChtml << " " << ICstats(cnum,3) << " % signal change (pos peak voxel); " << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;
}
@@ -57,8 +56,7 @@ namespace Melodic{
}
}
- volume4D<float> tempVol;
-
+ volume4D<float> tempVol;
if(mmodel.get_threshmaps().Storage()>0&&
(mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols()))
{//Output thresholded IC map
@@ -105,7 +103,6 @@ namespace Melodic{
IChtml << "<img BORDER=0 SRC=\""+mmodel.get_prefix()
+"_thresh.png\" ALT=\"MMfit\"></A><p>" << endl;
}
-
{//plot time course
IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
miscplot newplot;
@@ -174,7 +171,7 @@ namespace Melodic{
+string("f")+num2str(cnum)+".png\"></A><p>" << endl;
}//frequency plot
{//add T-mode GLM F-stats for full model fit & contrasts
- if(melodat.Tdes.Storage() > 0){
+ if(melodat.Tdes.Storage() > 0){
IChtml << " <TABLE border=1 bgcolor=ffffff cellpadding=5>" <<
"<CAPTION><EM> <b>GLM (OLS) on time series </b></EM></CAPTION>" << endl
<< "<TR valign=middle><TH ><EM>GLM β's</EM> <TH> <EM> F-test on <br> full model fit </em>";
@@ -198,7 +195,7 @@ namespace Melodic{
melodat.glmT.get_pf_fmf().Column(cnum) <<
"<BR> (uncorrected for #comp.)</TD>" << endl;
}
- if(melodat.Tcon.Storage() > 0){
+ if(melodat.Tcon.Storage() > 0){
IChtml << "<TD><TABLE border=0><TR><TD align=right>" <<endl;
for(int ctr=1; ctr <= melodat.Tcon.Nrows() ; ctr++)
IChtml << "con(" << melodat.Tcon.Row(ctr) << "): <br>" << endl;
@@ -218,8 +215,8 @@ namespace Melodic{
"<BR>" << endl;
IChtml << "</TABLE></td></tr>" << endl;
}
- IChtml << "</TABLE><p>" << endl;
- }
+ IChtml << "</TABLE><p>" << endl;
+ }
if(cnum <= (int)melodat.get_Smodes().size())
{//plot subject mode
@@ -295,7 +292,7 @@ namespace Melodic{
if(mmodel.get_threshmaps().Storage()>0&&
(mmodel.get_threshmaps().Ncols() == mmodel.get_data().Ncols())&&
(mmodel.get_threshmaps().Nrows()>1))
- /* {//Output other thresholded IC map
+ {//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());
@@ -347,10 +344,10 @@ namespace Melodic{
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 "
+ 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
@@ -473,6 +470,8 @@ namespace Melodic{
IChtml2 << "</A><p>" << endl;
}
+ RowVector dat = mmodel.get_data().Row(1);
+ if(dat.Maximum()>dat.Minimum())
{//Output GGM/GMM fit
miscplot newplot;