Skip to content
Snippets Groups Projects
Commit 34c76008 authored by Christian Beckmann's avatar Christian Beckmann
Browse files

Refs updated

parent 0e6ee78b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
fsl_glm.cc
+ 5
− 3
View file @ 34c76008
...@@ -27,14 +27,14 @@ using namespace std; ...@@ -27,14 +27,14 @@ using namespace std;
string(" \n Simple GLM usign ordinary least-squares regression on\n")+ string(" \n Simple GLM usign ordinary least-squares regression on\n")+
string(" time courses and/or 3D/4D imges against time courses \n")+ string(" time courses and/or 3D/4D imges against time courses \n")+
string(" or 3D/4D images\n\n"); string(" or 3D/4D images\n\n");
string examples="fsl_glm -i <input> -d <design> [options]"; string examples="fsl_glm -i <input> -d <design> -o <output> [options]";
//Command line Options { //Command line Options {
Option<string> fnin(string("-i,--in"), string(""), Option<string> fnin(string("-i,--in"), string(""),
string(" input file name (matrix 3D or 4D image)"), string(" input file name (matrix 3D or 4D image)"),
true, requires_argument); true, requires_argument);
Option<string> fnout(string("-o,--out"), string(""), Option<string> fnout(string("-o,--out"), string(""),
string(" output file name for GLM parameter estimates"), string("output file name for GLM parameter estimates"),
true, requires_argument); true, requires_argument);
Option<string> fndesign(string("-d,--design"), string(""), Option<string> fndesign(string("-d,--design"), string(""),
string("file name of the GLM design matrix (time courses or spatial maps)"), string("file name of the GLM design matrix (time courses or spatial maps)"),
...@@ -44,7 +44,7 @@ using namespace std; ...@@ -44,7 +44,7 @@ using namespace std;
false, requires_argument); false, requires_argument);
Option<string> fncontrasts(string("-c,--contrasts"), string(""), Option<string> fncontrasts(string("-c,--contrasts"), string(""),
string("matrix of t-statistics contrasts"), string("matrix of t-statistics contrasts"),
false, requires_argument); false, requires_argument, false);
Option<string> fnftest(string("-f,--ftests"), string(""), Option<string> fnftest(string("-f,--ftests"), string(""),
string("matrix of F-tests on contrasts"), string("matrix of F-tests on contrasts"),
false, requires_argument); false, requires_argument);
...@@ -131,6 +131,8 @@ bool isimage(Matrix what){ ...@@ -131,6 +131,8 @@ bool isimage(Matrix what){
void saveit(Matrix what, string fname){ void saveit(Matrix what, string fname){
if(isimage(what)) if(isimage(what))
save4D(what,fname); save4D(what,fname);
else if(fsl_imageexists(fndesign.value()))
write_ascii_matrix(what.t(),fname);
else else
write_ascii_matrix(what,fname); write_ascii_matrix(what,fname);
} }
......
meldata.cc
+ 34
− 14
View file @ 34c76008
...@@ -45,7 +45,7 @@ namespace Melodic{ ...@@ -45,7 +45,7 @@ namespace Melodic{
tmpData = RawData.matrix(Mask); tmpData = RawData.matrix(Mask);
//estimate smoothness //estimate smoothness
if(Resels == 0) if((Resels == 0)&&(!opts.filtermode))
Resels = est_resels(RawData,Mask); Resels = est_resels(RawData,Mask);
//convert to percent BOLD signal change //convert to percent BOLD signal change
...@@ -62,7 +62,7 @@ namespace Melodic{ ...@@ -62,7 +62,7 @@ namespace Melodic{
message(" done" << endl); message(" done" << endl);
} }
meanC = mean(tmpData,2); // meanC = mean(tmpData,2);
// tmpData = remmean(tmpData,2); // tmpData = remmean(tmpData,2);
//convert to power spectra //convert to power spectra
...@@ -134,15 +134,24 @@ namespace Melodic{ ...@@ -134,15 +134,24 @@ namespace Melodic{
void MelodicData::set_TSmode() void MelodicData::set_TSmode()
{ {
Matrix tmp, tmpT, tmpS, tmpT2, tmpS2; Matrix tmp, tmpT, tmpS, tmpT2, tmpS2, tmpT3;
tmp = expand_dimred(mixMatrix); tmp = expand_dimred(mixMatrix);
tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols()); tmpT = zeros(tmp.Nrows()/numfiles, tmp.Ncols());
tmpS = zeros(numfiles, tmp.Ncols()); tmpS = zeros(numfiles, tmp.Ncols());
krfact(tmp,tmpT,tmpS); explained_var = krfact(tmp,tmpT,tmpS);
outMsize("tmp",tmp);
outMsize("tmpT",tmpT);
outMsize("tmpS",tmpS);
Tmodes.clear(); Smodes.clear(); Tmodes.clear(); Smodes.clear();
for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){ for(int ctr = 1; ctr <= tmp.Ncols(); ctr++){
tmpT3 << reshape(tmp.Column(ctr),tmpT.Nrows(),numfiles);
outMsize("tmpT3", tmpT3);
tmpT2 << tmpT.Column(ctr); tmpT2 << tmpT.Column(ctr);
tmpS2 << tmpS.Column(ctr); tmpS2 << tmpS.Column(ctr);
tmpT3 << SP(tmpT3,pow(ones(tmpT3.Nrows(),1)*tmpS2.t(),-1));
if(numfiles>1)
tmpT2 |= tmpT3;
if(mean(tmpS2,1).AsScalar()<0){ if(mean(tmpS2,1).AsScalar()<0){
tmpT2*=-1.0; tmpT2*=-1.0;
tmpS2*=-1.0; tmpS2*=-1.0;
...@@ -153,9 +162,9 @@ namespace Melodic{ ...@@ -153,9 +162,9 @@ namespace Melodic{
//add GLM OLS fit //add GLM OLS fit
if(Tdes.Storage()){ if(Tdes.Storage()){
Matrix alltcs = Tmodes.at(0); Matrix alltcs = Tmodes.at(0).Column(1);
for(int ctr=1; ctr < (int)Tmodes.size();ctr++) for(int ctr=1; ctr < (int)Tmodes.size();ctr++)
alltcs|=Tmodes.at(ctr); alltcs|=Tmodes.at(ctr).Column(1);
if((alltcs.Nrows()==Tdes.Nrows())&&(Tdes.Nrows()>Tdes.Ncols())) if((alltcs.Nrows()==Tdes.Nrows())&&(Tdes.Nrows()>Tdes.Ncols()))
glmT.olsfit(alltcs,Tdes,Tcon); glmT.olsfit(alltcs,Tdes,Tcon);
} }
...@@ -170,8 +179,14 @@ namespace Melodic{ ...@@ -170,8 +179,14 @@ namespace Melodic{
void MelodicData::setup() void MelodicData::setup()
{ {
setup_misc();
numfiles = (int)opts.inputfname.value().size(); numfiles = (int)opts.inputfname.value().size();
setup_misc();
if(opts.filtermode){ // basic setup for filtering only
Data = process_file(opts.inputfname.value().at(0));
}
else{
if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm"))) if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
opts.approach.set_T("tica"); opts.approach.set_T("tica");
...@@ -288,15 +303,15 @@ namespace Melodic{ ...@@ -288,15 +303,15 @@ namespace Melodic{
message(endl << " Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl); message(endl << " Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
outMsize("stdDev",stdDev); outMsize("stdDev",stdDev);
meanC=mean(Data,2); //meanC=mean(Data,2);
if(opts.debug.value()) if(opts.debug.value())
save4D(Data,"concat_data"); save4D(Data,"concat_data");
//save the mean & mask //save the mean & mask
save_volume(Mask,logger.appendDir("mask")); save_volume(Mask,logger.appendDir("mask"));
save_volume(Mean,logger.appendDir("mean")); save_volume(Mean,logger.appendDir("mean"));
}
} // void setup() } // void setup()
void MelodicData::setup_misc() void MelodicData::setup_misc()
{ {
...@@ -348,6 +363,13 @@ namespace Melodic{ ...@@ -348,6 +363,13 @@ namespace Melodic{
if(opts.fn_SconF.value().length()>0) if(opts.fn_SconF.value().length()>0)
SconF = read_ascii_matrix(opts.fn_SconF.value()); SconF = read_ascii_matrix(opts.fn_SconF.value());
if(numfiles>1 && Sdes.Storage() == 0){
Sdes = ones(numfiles,1);
if(Scon.Storage() == 0){
Scon = ones(1,1);
Scon &= -1*Scon;
}
}
Tdes = remmean(Tdes,1); Tdes = remmean(Tdes,1);
} }
...@@ -378,7 +400,6 @@ namespace Melodic{ ...@@ -378,7 +400,6 @@ namespace Melodic{
if((IC.Storage()>0)&&(opts.output_origIC.value())&&(after_mm==false)) if((IC.Storage()>0)&&(opts.output_origIC.value())&&(after_mm==false))
save4D(IC,opts.outputfname.value() + "_oIC"); save4D(IC,opts.outputfname.value() + "_oIC");
//Output IC -- adjusted for noise //Output IC -- adjusted for noise
if(IC.Storage()>0){ if(IC.Storage()>0){
volume4D<float> tempVol; volume4D<float> tempVol;
...@@ -412,7 +433,7 @@ namespace Melodic{ ...@@ -412,7 +433,7 @@ namespace Melodic{
//Output mixMatrix //Output mixMatrix
if(mixMatrix.Storage()>0){ if(mixMatrix.Storage()>0){
saveascii(mixMatrix, opts.outputfname.value() + "_mix"); saveascii(expand_mix(), opts.outputfname.value() + "_mix");
mixFFT=calc_FFT(expand_mix(), opts.logPower.value()); mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
saveascii(mixFFT,opts.outputfname.value() + "_FTmix"); saveascii(mixFFT,opts.outputfname.value() + "_FTmix");
} }
...@@ -521,8 +542,7 @@ namespace Melodic{ ...@@ -521,8 +542,7 @@ namespace Melodic{
outMsize("meanC",meanC); outMsize("meanC",meanC);
newData = Data - noiseMix * noiseIC.t(); newData = Data - noiseMix * noiseIC.t();
if(meanC.Storage()>0)
newData = newData + meanC*ones(1,newData.Ncols());
if(meanR.Storage()>0) if(meanR.Storage()>0)
newData = newData + ones(newData.Nrows(),1)*meanR; newData = newData + ones(newData.Nrows(),1)*meanR;
......
meldata.h
+ 3
− 2
View file @ 34c76008
...@@ -36,7 +36,7 @@ namespace Melodic{ ...@@ -36,7 +36,7 @@ namespace Melodic{
void save(); void save();
Matrix process_file(string fname, int numfiles); Matrix process_file(string fname, int numfiles = 1);
inline void save4D(Matrix what, string fname){ inline void save4D(Matrix what, string fname){
volume4D<float> tempVol; volume4D<float> tempVol;
...@@ -193,7 +193,8 @@ namespace Melodic{ ...@@ -193,7 +193,8 @@ namespace Melodic{
volumeinfo tempInfo; volumeinfo tempInfo;
vector<Matrix> DWM, WM; vector<Matrix> DWM, WM;
basicGLM glmT, glmS; basicGLM glmT, glmS;
Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF; Matrix Tdes, Tcon, TconF, Sdes, Scon, SconF;
RowVector explained_var;
private: private:
MelodicOptions &opts; MelodicOptions &opts;
......
melhlprfns.cc
+ 49
− 31
View file @ 34c76008
...@@ -316,7 +316,7 @@ namespace Melodic{ ...@@ -316,7 +316,7 @@ namespace Melodic{
evecs = C * Evc; evecs = C * Evc;
} //void em_pca } //void em_pca
void rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim) float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim)
{ {
Matrix Corr, Evecs, tmpWM, tmpDWM, tmp; Matrix Corr, Evecs, tmpWM, tmpDWM, tmp;
RowVector Evals; RowVector Evals;
...@@ -325,29 +325,36 @@ namespace Melodic{ ...@@ -325,29 +325,36 @@ namespace Melodic{
tmp = tmpWM * Mat.t(); tmp = tmpWM * Mat.t();
cols = tmp.t(); cols = tmp.t();
rows << tmpDWM; rows << tmpDWM;
float res;
Evals=fliplr(Evals);
res = sum(Evals.Columns(1,dim),2).AsScalar()/sum(Evals,2).AsScalar()*100;
return res;
} // rankapprox } // rankapprox
void krfact(const Matrix& Mat, Matrix& cols, Matrix& rows) RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows)
{ {
Matrix out; Matrix out; RowVector res(Mat.Ncols());
for(int ctr1 = 1; ctr1 <= Mat.Ncols(); ctr1++) for(int ctr1 = 1; ctr1 <= Mat.Ncols(); ctr1++){
{ Matrix tmpVals(cols.Nrows(),rows.Nrows());
Matrix tmpVals(cols.Nrows(),rows.Nrows()); for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++)
for(int ctr2 = 1; ctr2 <= rows.Nrows(); ctr2++) tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() *
tmpVals.Column(ctr2) << Mat.SubMatrix(cols.Nrows() * (ctr2 - 1) + 1,cols.Nrows()*ctr2 ,ctr1,ctr1); (ctr2 - 1) + 1,cols.Nrows()*ctr2 ,ctr1,ctr1);
Matrix tmpcols, tmprows; Matrix tmpcols, tmprows;
rankapprox(tmpVals, tmpcols, tmprows); res(ctr1) =rankapprox(tmpVals, tmpcols, tmprows);
cols.Column(ctr1) = tmpcols; cols.Column(ctr1) = tmpcols;
rows.Column(ctr1) = tmprows; rows.Column(ctr1) = tmprows;
} }
return res;
} // krfact } // krfact
void krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows) RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows)
{ {
RowVector res;
cols = zeros(colnum,Mat.Ncols()); cols = zeros(colnum,Mat.Ncols());
rows = zeros(int(Mat.Nrows() / colnum),Mat.Ncols()); rows = zeros(int(Mat.Nrows() / colnum),Mat.Ncols());
krfact(Mat,cols,rows); res = krfact(Mat,cols,rows);
return res;
} // krfact } // krfact
Matrix krprod(const Matrix& cols, const Matrix& rows) Matrix krprod(const Matrix& cols, const Matrix& rows)
...@@ -558,35 +565,46 @@ namespace Melodic{ ...@@ -558,35 +565,46 @@ namespace Melodic{
int res = 0; int res = 0;
ColumnVector PPCA; ColumnVector PPCA;
if(which == string("aut"))
if(int(estimators(2)) < int(estimators(1)) && int(estimators(2)) > 15){
res=int(estimators(2));
PPCA << PPCA2.Column(3);
}else{
res = int(estimators(1));
PPCA << PPCA2.Column(2);
}
if(which == string("lap")){ if(which == string("lap")){
res = int(estimators(1)); res = int(estimators(1));
PPCA << PPCA2.Column(2); PPCA << PPCA2.Column(2);
} }
if(which == string("bic")){ if(which == string("bic")){
res = int(estimators(2)); res = int(estimators(2));
PPCA << PPCA2.Column(2); PPCA << PPCA2.Column(3);
} }
if(which == string("mdl")){ if(which == string("mdl")){
res = int(estimators(3)); res = int(estimators(3));
PPCA << PPCA2.Column(4); PPCA << PPCA2.Column(4);
} }
if(which == string("rrn")){
res = int(estimators(4));
PPCA << PPCA2.Column(5);
}
if(which == string("aic")){ if(which == string("aic")){
res = int(estimators(5)); res = int(estimators(5));
PPCA << PPCA2.Column(6); PPCA << PPCA2.Column(6);
} }
if(res==0){//median estimator if(which == string("median")){
PPCA = PPCA2.Column(2); RowVector unsorted(estimators);
SortAscending(unsorted);
for(int ctr=1; ctr<=PPCA2.Nrows(); ctr++){ ctr_i=1;
RowVector tmp = PPCA2.SubMatrix(ctr,ctr,2,6); res=int(unsorted(3));
PPCA(ctr) = float(tmp.Sum()/5); while(res != int(estimators(ctr_i)))
} ctr_i++;
PPCA << PPCA2.Column(ctr_i);
ctr_i = 1; }
while((PPCA(ctr_i) < PPCA(ctr_i+1))&&(ctr_i<maxEV)){ if(res==0 || which == string("mean")){//median estimator
res=ctr_i+1;ctr_i++; PPCA = mean(PPCA2.Columns(2,6),2);
} res=int(mean(estimators,2).AsScalar());
} }
dim = res; dim = res;
...@@ -812,10 +830,10 @@ namespace Melodic{ ...@@ -812,10 +830,10 @@ namespace Melodic{
beta = pinvdes * dat; beta = pinvdes * dat;
residu = dat - design*beta; residu = dat - design*beta;
dof = design.Nrows() - design.Ncols(); dof = design.Nrows() - design.Ncols()-1;
sigsq = sum(SP(residu,residu))/dof; sigsq = sum(SP(residu,residu))/dof;
float fact = float(design.Nrows() - 1 - design.Ncols()) / design.Ncols(); float fact = float(dof) / design.Ncols();
f_fmf = SP(sum(SP(design*beta,design*beta)),pow(sum(SP(residu,residu)),-1)) * fact; f_fmf = SP(sum(SP(design*beta,design*beta)),pow(sum(SP(residu,residu)),-1)) * fact;
pf_fmf = f_fmf.Row(1); pf_fmf = f_fmf.Row(1);
......
melhlprfns.h
+ 3
− 3
View file @ 34c76008
...@@ -50,9 +50,9 @@ namespace Melodic{ ...@@ -50,9 +50,9 @@ namespace Melodic{
void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20); void em_pca(const Matrix& Mat, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20); void em_pca(const Matrix& Mat, Matrix& guess, Matrix& evecs, RowVector& evals, int num_pc = 1, int iter = 20);
void rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim = 1); float rankapprox(const Matrix& Mat, Matrix& cols, Matrix& rows, int dim = 1);
void krfact(const Matrix& Mat, Matrix& cols, Matrix& rows); RowVector krfact(const Matrix& Mat, Matrix& cols, Matrix& rows);
void krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows); RowVector krfact(const Matrix& Mat, int colnum, Matrix& cols, Matrix& rows);
Matrix krprod(const Matrix& cols, const Matrix& rows); Matrix krprod(const Matrix& cols, const Matrix& rows);
Matrix krapprox(const Matrix& Mat, int size_col, int dim = 1); Matrix krapprox(const Matrix& Mat, int size_col, int dim = 1);
......
melodic.cc
+ 10
− 13
View file @ 34c76008
...@@ -71,13 +71,12 @@ int main(int argc, char *argv[]){ ...@@ -71,13 +71,12 @@ int main(int argc, char *argv[]){
if (opts.filtermode || opts.filtermix.value().length()>0){ if (opts.filtermode || opts.filtermix.value().length()>0){
if(opts.filtermode){ // just filter out some noise from a previous run if(opts.filtermode){ // just filter out some noise from a previous run
melodat.setup(); melodat.setup();
if(opts.debug.value())
message(" Denoising data setup completed "<< endl);
melodat.remove_components(); melodat.remove_components();
} } else
else
mmonly(logger,opts,melodat,report); mmonly(logger,opts,melodat,report);
} } else { // standard PICA now
else
{ // standard PICA now
int retry = 0; int retry = 0;
bool no_conv; bool no_conv;
bool leaveloop = false; bool leaveloop = false;
...@@ -104,12 +103,10 @@ int main(int argc, char *argv[]){ ...@@ -104,12 +103,10 @@ int main(int argc, char *argv[]){
opts.approach.set_T("defl"); opts.approach.set_T("defl");
message(endl << "Restarting MELODIC using deflation approach" message(endl << "Restarting MELODIC using deflation approach"
<< endl << endl); << endl << endl);
} }else{
else{
leaveloop = true; leaveloop = true;
} }
} }else{
else{
if(no_conv){ if(no_conv){
retry++; retry++;
if(opts.pca_dim.value()-retry*opts.retrystep > if(opts.pca_dim.value()-retry*opts.retrystep >
...@@ -125,10 +122,10 @@ int main(int argc, char *argv[]){ ...@@ -125,10 +122,10 @@ int main(int argc, char *argv[]){
} }
} }
if(!leaveloop){ if(!leaveloop){
message(endl << "Restarting MELODIC using -d " message(endl << "Restarting MELODIC using -d "
<< opts.pca_dim.value() << opts.pca_dim.value()
<< endl << endl); << endl << endl);
} }
} }
} }
} while (no_conv && retry<opts.maxRestart.value() && !leaveloop); } while (no_conv && retry<opts.maxRestart.value() && !leaveloop);
......
meloptions.cc
+ 5
− 2
View file @ 34c76008
...@@ -32,7 +32,7 @@ MelodicOptions* MelodicOptions::gopt = NULL; ...@@ -32,7 +32,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
version = p_version; version = p_version;
filtermode = false; filtermode = false;
explicitnums = false; explicitnums = false;
logfname = string("melodic.log"); logfname = string("log.txt");
// work out the path to the $FSLDIR/bin directory // work out the path to the $FSLDIR/bin directory
if(getenv("FSLDIR")!=0){ if(getenv("FSLDIR")!=0){
...@@ -112,8 +112,11 @@ MelodicOptions* MelodicOptions::gopt = NULL; ...@@ -112,8 +112,11 @@ MelodicOptions* MelodicOptions::gopt = NULL;
print_usage(argc,argv); print_usage(argc,argv);
exit(2); exit(2);
} else { } else {
temporal.set_T(false);
filtermode = true; filtermode = true;
varnorm.set_T(false); varnorm.set_T(false);
pbsc.set_T(false);
cerr << "WARNING: melodic denoising is deprecated, please use fsl_regfilt instead!" <<endl;
} }
} }
if (threshold.value()<=0){ if (threshold.value()<=0){
...@@ -182,7 +185,7 @@ MelodicOptions* MelodicOptions::gopt = NULL; ...@@ -182,7 +185,7 @@ MelodicOptions* MelodicOptions::gopt = NULL;
inputfname.set_T(tmpfnames); inputfname.set_T(tmpfnames);
//create melodic directory name //create melodic directory name
if(logdir.value()==string("melodic.log")){ if(logdir.value()==string("log.txt")){
logdir.set_T(string(inputfname.value().at(0)+".ica")); logdir.set_T(string(inputfname.value().at(0)+".ica"));
logger.makeDir(logdir.value(),logfname); logger.makeDir(logdir.value(),logfname);
} else{ } else{
......
meloptions.h
+ 14
− 9
View file @ 34c76008
...@@ -81,6 +81,7 @@ class MelodicOptions { ...@@ -81,6 +81,7 @@ class MelodicOptions {
Option<string> bgimage; Option<string> bgimage;
Option<float> tr; Option<float> tr;
Option<bool> logPower; Option<bool> logPower;
Option<bool> addsigchng;
Option<string> fn_Tdesign; Option<string> fn_Tdesign;
Option<string> fn_Tcon; Option<string> fn_Tcon;
...@@ -106,8 +107,8 @@ class MelodicOptions { ...@@ -106,8 +107,8 @@ class MelodicOptions {
Option<string> guessfname; Option<string> guessfname;
Option<string> paradigmfname; Option<string> paradigmfname;
Option<int> dummy; Option<int> dummy;
Option<int> repeats; Option<int> repeats;
Option<float> nlconst1; Option<float> nlconst1;
Option<float> nlconst2; Option<float> nlconst2;
Option<float> smooth_probmap; Option<float> smooth_probmap;
...@@ -147,7 +148,7 @@ class MelodicOptions { ...@@ -147,7 +148,7 @@ class MelodicOptions {
} }
inline MelodicOptions::MelodicOptions() : inline MelodicOptions::MelodicOptions() :
logdir(string("-o,--outdir"), string("melodic.log"), logdir(string("-o,--outdir"), string("log.txt"),
string("output directory name\n"), string("output directory name\n"),
false, requires_argument), false, requires_argument),
inputfname(string("-i,--in"), std::vector<string>(), inputfname(string("-i,--in"), std::vector<string>(),
...@@ -180,8 +181,8 @@ class MelodicOptions { ...@@ -180,8 +181,8 @@ class MelodicOptions {
joined_whiten(string("--sep_whiten"), true, joined_whiten(string("--sep_whiten"), true,
string("switch on separate whitening"), string("switch on separate whitening"),
false, no_argument), false, no_argument),
joined_vn(string("--joined_vn"), false, joined_vn(string("--sep_vn"), true,
string("switch on joined variance nomalisation"), string("switch off joined variance nomalisation"),
false, no_argument, false), false, no_argument, false),
numICs(string("-n,--numICs"), -1, numICs(string("-n,--numICs"), -1,
string("numer of IC's to extract (for deflation approach)"), string("numer of IC's to extract (for deflation approach)"),
...@@ -200,7 +201,7 @@ class MelodicOptions { ...@@ -200,7 +201,7 @@ class MelodicOptions {
false, no_argument), false, no_argument),
pspec(string("--pspec"), false, pspec(string("--pspec"), false,
string(" switch on conversion to powerspectra"), string(" switch on conversion to powerspectra"),
false, no_argument), false, no_argument, false),
segment(string("--covarweight"), string(""), segment(string("--covarweight"), string(""),
string("voxel-wise weights for the covariance matrix (e.g. segmentation information)"), string("voxel-wise weights for the covariance matrix (e.g. segmentation information)"),
false, requires_argument), false, requires_argument),
...@@ -249,9 +250,12 @@ class MelodicOptions { ...@@ -249,9 +250,12 @@ class MelodicOptions {
tr(string("--tr"), 0.0, tr(string("--tr"), 0.0,
string("\tTR in seconds"), string("\tTR in seconds"),
false, requires_argument), false, requires_argument),
logPower(string("--logPower"), false, logPower(string("--logPower"), false,
string("calculate log of power for frequency spectrum\n"), string("calculate log of power for frequency spectrum\n"),
false, no_argument), false, no_argument),
addsigchng(string("--sigchng"), false,
string("add signal change estimates to report pages\n"),
false, no_argument, false),
fn_Tdesign(string("--Tdes"), string(""), fn_Tdesign(string("--Tdes"), string(""),
string(" design matrix across time-domain"), string(" design matrix across time-domain"),
false, requires_argument), false, requires_argument),
...@@ -382,6 +386,7 @@ class MelodicOptions { ...@@ -382,6 +386,7 @@ class MelodicOptions {
options.add(bgimage); options.add(bgimage);
options.add(tr); options.add(tr);
options.add(logPower); options.add(logPower);
options.add(addsigchng);
options.add(fn_Tdesign); options.add(fn_Tdesign);
options.add(fn_Tcon); options.add(fn_Tcon);
options.add(fn_TconF); options.add(fn_TconF);
......
melreport.cc
+ 34
− 8
View file @ 34c76008
...@@ -53,7 +53,7 @@ namespace Melodic{ ...@@ -53,7 +53,7 @@ namespace Melodic{
IChtml << fixed << setprecision(2) << ICstats(cnum,1) << " % of explained variance"; IChtml << fixed << setprecision(2) << ICstats(cnum,1) << " % of explained variance";
if(ICstats.Ncols()>1) if(ICstats.Ncols()>1)
IChtml << "; &nbsp; &nbsp; " << ICstats(cnum,2) << " % of total variance"; IChtml << "; &nbsp; &nbsp; " << ICstats(cnum,2) << " % of total variance";
if(ICstats.Ncols()>2&&melodat.get_numfiles()==1){ if(ICstats.Ncols()>2&&opts.addsigchng.value()){
IChtml << "<p>" <<endl; IChtml << "<p>" <<endl;
IChtml << " &nbsp; &nbsp; " << ICstats(cnum,3) << " % signal change (pos peak voxel); &nbsp; &nbsp;" << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ; IChtml << " &nbsp; &nbsp; " << ICstats(cnum,3) << " % signal change (pos peak voxel); &nbsp; &nbsp;" << ICstats(cnum,4) << "% signal change (peak neg. voxel)" << endl ;
} }
...@@ -112,16 +112,27 @@ namespace Melodic{ ...@@ -112,16 +112,27 @@ namespace Melodic{
{//plot time course {//plot time course
IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl; IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
miscplot newplot; miscplot newplot;
Matrix tmptc = melodat.get_Tmodes(cnum-1).t(); Matrix tmptc = melodat.get_Tmodes(cnum-1).Column(1).t();
newplot.col_replace(0,0xFF0000);
newplot.add_label(string("IC ")+num2str(cnum)+" time course");
//add GLM OLS fit //add GLM OLS fit
if(melodat.Tdes.Storage()>0 && if(melodat.Tdes.Storage()>0 &&
melodat.glmT.get_beta().Nrows() == melodat.Tdes.Ncols()){ melodat.glmT.get_beta().Nrows() == melodat.Tdes.Ncols()){
tmptc &= melodat.glmT.get_beta().Column(cnum).t() * melodat.Tdes.t(); tmptc &= melodat.glmT.get_beta().Column(cnum).t() * melodat.Tdes.t();
newplot.add_label(string("IC ")+num2str(cnum)+" time course");
newplot.add_label("full model fit"); newplot.add_label("full model fit");
} }
//add deviation around time course
if(melodat.get_Tmodes(cnum-1).Ncols()>1){
Matrix tmp = stdev(melodat.get_Tmodes(cnum-1).Columns(2,melodat.get_Tmodes(cnum-1).Ncols()).t(),1);
tmptc &= melodat.get_Tmodes(cnum-1).Column(1).t()+tmp;
tmptc &= melodat.get_Tmodes(cnum-1).Column(1).t()-tmp;
newplot.add_label("std error across subjects");
newplot.col_replace(tmptc.Nrows()-1,0x808080);
newplot.col_replace(tmptc.Nrows()-2,0x808080);
}
if(opts.tr.value()>0.0) if(opts.tr.value()>0.0)
newplot.add_xlabel(string("Time (seconds); TR = ")+ newplot.add_xlabel(string("Time (seconds); TR = ")+
float2str(opts.tr.value(),0,2,0)+" s"); float2str(opts.tr.value(),0,2,0)+" s");
...@@ -137,12 +148,19 @@ namespace Melodic{ ...@@ -137,12 +148,19 @@ namespace Melodic{
string("Timecourse No. ")+num2str(cnum), string("Timecourse No. ")+num2str(cnum),
opts.tr.value(),150,12,1,false); opts.tr.value(),150,12,1,false);
if(melodat.get_Tmodes(cnum-1).Ncols()>1)
tmptc &= melodat.get_Tmodes(cnum-1).Columns(2,melodat.get_Tmodes(cnum-1).Ncols()).t();
write_ascii_matrix(report.appendDir(string("t") write_ascii_matrix(report.appendDir(string("t")
+num2str(cnum)+".txt"),tmptc.t()); +num2str(cnum)+".txt"),tmptc.t());
IChtml << "<A HREF=\"" << string("t") IChtml << "<A HREF=\"" << string("t")
+num2str(cnum)+".txt" << "\"> "; +num2str(cnum)+".txt" << "\"> ";
IChtml << "<img BORDER=0 SRC=\"" IChtml << "<img BORDER=0 SRC=\""
+string("t")+num2str(cnum)+".png\"></A><p>" << endl; +string("t")+num2str(cnum)+".png\"></A><p>" << endl;
if(melodat.get_numfiles()>1 && melodat.explained_var.Storage()>0
&& melodat.explained_var.Ncols()>=cnum)
IChtml << "Rank-1 approximation of individual time courses explains "
<< melodat.explained_var(cnum) << "% of variance.<p>" << endl;
}//time series plot }//time series plot
if(!opts.pspec.value()) if(!opts.pspec.value())
...@@ -156,7 +174,7 @@ namespace Melodic{ ...@@ -156,7 +174,7 @@ namespace Melodic{
else else
newplot.add_ylabel(string("Power")); newplot.add_ylabel(string("Power"));
Matrix fmixtc = calc_FFT(melodat.get_Tmodes(cnum-1), opts.logPower.value()); Matrix fmixtc = calc_FFT(melodat.get_Tmodes(cnum-1).Column(1), opts.logPower.value());
newplot.set_Ylabel_fmt("%.0f"); newplot.set_Ylabel_fmt("%.0f");
newplot.set_yrange(0.0,1.02*fmixtc.Maximum()); newplot.set_yrange(0.0,1.02*fmixtc.Maximum());
...@@ -239,11 +257,13 @@ namespace Melodic{ ...@@ -239,11 +257,13 @@ namespace Melodic{
if(cnum <= (int)melodat.get_Smodes().size()) if(cnum <= (int)melodat.get_Smodes().size())
{//plot subject mode {//plot subject mode
Matrix smode; Matrix smode;
smode = melodat.get_Smodes(cnum-1); smode = melodat.get_Smodes(cnum-1);
if(smode.Nrows() > 1){ if(smode.Nrows() > 1){
miscplot newplot; IChtml << "<hr><H3> Sessions/Subjects mode </H3><p>" << endl <<endl;
miscplot newplot;
//add GLM OLS fit //add GLM OLS fit
//newplot.setscatter(smode,2); //newplot.setscatter(smode,2);
...@@ -260,7 +280,7 @@ namespace Melodic{ ...@@ -260,7 +280,7 @@ namespace Melodic{
// newplot.set_xysize(smode.Nrows()*80,150); // newplot.set_xysize(smode.Nrows()*80,150);
newplot.timeseries(smode.t(), newplot.timeseries(smode.t(),
report.appendDir(string("s")+num2str(cnum)+".png"), report.appendDir(string("s")+num2str(cnum)+".png"),
string("Subject/Session mode")); string("Subject/Session mode No. ") + num2str(cnum));
newplot.clear_xlabel(); newplot.clear_xlabel();
newplot.clear_labels(); newplot.clear_labels();
newplot.set_xysize(120,200); newplot.set_xysize(120,200);
...@@ -722,7 +742,7 @@ namespace Melodic{ ...@@ -722,7 +742,7 @@ namespace Melodic{
void MelodicReport::PPCA_rep(){ void MelodicReport::PPCA_rep(){
{//plot time course {//plot time course
report << "<p><hr><b>PCA estimates </b> <br>" << endl; report << "<p><hr><b>PCA estimates </b> <p>" << endl;
Matrix what; Matrix what;
miscplot newplot; miscplot newplot;
...@@ -753,12 +773,18 @@ namespace Melodic{ ...@@ -753,12 +773,18 @@ namespace Melodic{
} }
void MelodicReport::Smode_rep(){ void MelodicReport::Smode_rep(){
report << "<b>TICA Subject/Session modes </b> <p>" << endl; report << "<p><hr><b>TICA Subject/Session modes </b> <br>" << endl;
miscplot newplot; miscplot newplot;
report << "Boxplots show the relative response amplitudes across the "
<< "session/subject domain (" << melodat.get_numfiles()
<< " input files). Components have been sorted in decreasing order of "
<< " the median response per component. <br><br>";
Matrix allmodes = melodat.get_Smodes().at(0); Matrix allmodes = melodat.get_Smodes().at(0);
for(int ctr = 1; ctr < (int)melodat.get_Smodes().size();++ctr) for(int ctr = 1; ctr < (int)melodat.get_Smodes().size();++ctr)
allmodes |= melodat.get_Smodes().at(ctr); allmodes |= melodat.get_Smodes().at(ctr);
newplot.add_xlabel("Component No.");
newplot.add_ylabel("");
newplot.set_xysize(100+30*allmodes.Ncols(),300); newplot.set_xysize(100+30*allmodes.Ncols(),300);
newplot.boxplot(allmodes,report.appendDir(string("bp_Smodes.png")), newplot.boxplot(allmodes,report.appendDir(string("bp_Smodes.png")),
string("Subject/Session modes")); string("Subject/Session modes"));
......
melreport.h
+ 88
− 27
View file @ 34c76008
...@@ -82,12 +82,12 @@ namespace Melodic{ ...@@ -82,12 +82,12 @@ namespace Melodic{
} }
report << "<IFRAME height=80px width=100% src=nav.html frameborder=0></IFRAME><p>"<< endl; report << "<IFRAME height=80px width=100% src=nav.html frameborder=0></IFRAME><p>"<< endl;
loghtml << "<IFRAME height=100px width=100% src=nav.html frameborder=0></IFRAME><p>" loghtml << "<IFRAME height=100px width=100% src=nav.html frameborder=0></IFRAME><p>"
<<"<pre>../melodic.log</pre>" <<endl; <<"<IFRAME width=100% height=100% src=\"../log.txt\"></IFRAME>" <<endl;
navigator <<"<CENTER><TABLE BORDER=0><TR>" << endl navigator <<"<CENTER><TABLE BORDER=0><TR>" << endl
<<"<TD ALIGN=CENTER WIDTH=100%><FONT SIZE=-1>"<<endl <<"<TD ALIGN=CENTER WIDTH=100%><FONT SIZE=-1>"<<endl
<<"<A HREF=\"00index.html\" target=\"_top\">Main</A>&nbsp;-&nbsp;"; <<"<A HREF=\"00index.html\" target=\"_top\">Main</A>&nbsp;-&nbsp;";
// if(opts.guireport.value()=="") if(opts.guireport.value()=="")
// navigator << "<A HREF=\"log.html\" target=\"_top\">Log</A>&nbsp;-&nbsp;"; navigator << "<A HREF=\"log.html\" target=\"_top\">Log</A>&nbsp;-&nbsp;";
navigator <<"Components: "; navigator <<"Components: ";
navigator.flush(); navigator.flush();
} }
...@@ -113,9 +113,23 @@ namespace Melodic{ ...@@ -113,9 +113,23 @@ namespace Melodic{
if( bool(opts.genreport.value()) ){ if( bool(opts.genreport.value()) ){
report << "<b>Analysis methods</b> <br>"<<endl; report << "<b>Analysis methods</b> <br>"<<endl;
report << "Analysis was carried out using MELODIC (Multivariate Exploratory Linear Decomposition into Independent Components) Version "<< version <<", part of FSL (FMRIB's Software Library, <A HREF=\"http://www.fmrib.ox.ac.uk/fsl/\">www.fmrib.ox.ac.uk/fsl</A>), an implementation for the estimation of a Probabilistic Independent Component Analysis model [Beckmann 2004]."<<endl; report << "Analysis was carried out using ";
report << "The following melodic pre-processing was applied to the input data file: "<< endl; if(opts.approach.value() != string("tica"))
report << "Probabilistic Independent Component Analysis"
<<" [Beckmann 2004] as implemented in "<<endl;
else
report << "Tensorial Independent Component Analysis "
<<"[Beckmann 2005] as implemented in "<< endl;
report << " MELODIC (Multivariate"
<<" Exploratory Linear Decomposition into Independent Components)"
<<" Version "<< version <<", part of FSL (FMRIB's Software"
<<" Library, <A HREF=\"http://www.fmrib.ox.ac.uk/fsl/\">"
<<"www.fmrib.ox.ac.uk/fsl</A>).<br>";
report << "The following data pre-processing was applied to"
<<" the input data: "<< endl;
if(opts.use_mask.value()) if(opts.use_mask.value())
report << " masking of non-brain voxels;"; report << " masking of non-brain voxels;";
...@@ -124,35 +138,52 @@ namespace Melodic{ ...@@ -124,35 +138,52 @@ namespace Melodic{
if(opts.varnorm.value()) if(opts.varnorm.value())
report << " normalisation of the voxel-wise variance; "; report << " normalisation of the voxel-wise variance; ";
if(opts.pbsc.value())
report << " conversion to %BOLD signal change; ";
report << "<br>"<<endl; report << "<br>"<<endl;
report << " Pre-processed data was whitened and projected into a " report << " Pre-processed data were whitened and projected into a "
<< melodat.get_mix().Ncols()<< "-dimensional subspace using "; << melodat.get_mix().Ncols()<< "-dimensional subspace using ";
if(melodat.get_PPCA().Storage()>0){ if(melodat.get_PPCA().Storage()>0){
report << "probabilistic Principal Component Analysis where the number of dimensions was estimated using "; report << "probabilistic Principal Component Analysis where the"
<<" number of dimensions was estimated using ";
if(opts.pca_est.value() == string("lap")) if(opts.pca_est.value() == string("lap"))
report << "the Laplace approximation to the Bayesian evidence of the model order [Minka 2000, Beckmann 2004]. " << endl; report << "the Laplace approximation to the Bayesian"
<<" evidence of the model order [Minka 2000, Beckmann 2004]. " << endl;
else else
if(opts.pca_est.value() == string("bic")) if(opts.pca_est.value() == string("bic"))
report << "the <em> Bayesian Information Criterion</em> (BIC) [Kass 1993]. " << endl; report << "the <em> Bayesian Information Criterion</em>"
<<" (BIC) [Kass 1993]. " << endl;
else else
if(opts.pca_est.value() == string("mdl")) if(opts.pca_est.value() == string("mdl"))
report << "<em> Minimum Description Length</em> (MDL) [Rissanen 1978]. " << endl; report << "<em> Minimum Description Length</em> (MDL)"
<<" [Rissanen 1978]. " << endl;
else else
if(opts.pca_est.value() == string("aic")) if(opts.pca_est.value() == string("aic"))
report << "the <em> Akaike Information Criterion</em> (AIC) [Akaike 1969]. " << endl; report << "the <em> Akaike Information Criterion</em>"
<<" (AIC) [Akaike 1969]. " << endl;
else else
report << "a committee of approximations to Bayesian the model order [Beckmann 2004]. " << endl; report << " approximations to Bayesian the"
<<" model order [Beckmann 2004]. " << endl;
} }
else else
report << "Principal Component Analysis. "; report << "Principal Component Analysis. ";
report << " The whitened observations were decomposed into a set of time-courses and spatial maps by optimising for non-Gaussian spatial source distributions using a fixed-point iteration technique [Hyv&auml;rinen 1999]. " << endl; report << " <BR>The whitened observations were decomposed into "
report << "Estimated Component maps were divided by the standard deviation of the residual noise"; <<" sets of vectors which describe signal variation across"
<<" the temporal domain (time-courses)";
if(opts.approach.value() == string("tica") ||
opts.approach.value() == string("concat"))
report << ", the session/subject domain ";
report <<" and across the spatial domain (maps) by optimising for"
<<" non-Gaussian spatial source distributions using a"
<<" fixed-point iteration technique [Hyv&auml;rinen 1999]. " << endl;
report << "Estimated Component maps were divided by the standard"
<<" deviation of the residual noise";
if(opts.perf_mm.value()) if(opts.perf_mm.value())
report << " and thresholded by fitting a mixture model to the histogram of intensity values [Beckmann 2004]. <p>" << endl; report << " and thresholded by fitting a mixture model "
<<"to the histogram of intensity values [Beckmann 2004]. <p>" << endl;
else else
report <<".<p>" << endl; report <<".<p>" << endl;
...@@ -164,26 +195,56 @@ namespace Melodic{ ...@@ -164,26 +195,56 @@ namespace Melodic{
if( bool(opts.genreport.value()) ){ if( bool(opts.genreport.value()) ){
report << "<b>References</b> <br>"<<endl; report << "<b>References</b> <br>"<<endl;
report << "[Hyv&auml;rinen 1999] A. Hyv&auml;rinen. Fast and Robust Fixed-Point Algorithms for Independent Component Analysis. IEEE Transactions on Neural Networks 10(3):626-634, 1999.<br> " << endl; report << "[Hyv&auml;rinen 1999] A. Hyv&auml;rinen. Fast and"
report << "[Beckmann 2004] C.F. Beckmann and S.M. Smith. Probabilistic Independent Component Analysis for Functional Magnetic Resonance Imaging. IEEE Transactions on Medical Imaging 23(2):137-152 2004. <br>" << endl; <<" Robust Fixed-Point Algorithms for Independent Component"
<<" Analysis. IEEE Transactions on Neural Networks 10(3):"
<<"626-634, 1999.<br> " << endl;
report << "[Beckmann 2004] C.F. Beckmann and S.M. Smith."
<<" Probabilistic Independent Component Analysis for Functional"
<<" Magnetic Resonance Imaging. IEEE Transactions on Medical"
<<" Imaging 23(2):137-152 2004. <br>" << endl;
if(opts.approach.value() == string("tica") ||
opts.approach.value() == string("concat") )
report << "[Beckmann 2005] C.F. Beckmann and S.M. Smith."
<<" Tensorial extensions of independent component analysis"
<< " for multisubject FMRI analysis. Neuroimage "
<< " 25(1):294-311 2005. <br>";
if(melodat.get_PPCA().Storage()>0){ if(melodat.get_PPCA().Storage()>0){
report << "[Everson 2000] R. Everson and S. Roberts. Inferring the eigenvalues of covariance matrices from limited, noisy data. IEEE Trans Signal Processing, 48(7):2083-2091, 2000<br>"<<endl; report << "[Everson 2000] R. Everson and S. Roberts."
report << "[Tipping 1999] M.E. Tipping and C.M.Bishop. Probabilistic Principal component analysis. J Royal Statistical Society B, 61(3), 1999. <br>" << endl; <<" Inferring the eigenvalues of covariance matrices from"
/* report << "[Beckmann 2001] C.F. Beckmann, J.A. Noble and S.M. Smith. Investigating the intrinsic dimensionality of FMRI data for ICA. In Seventh Int. Conf. on Functional Mapping of the Human Brain, 2001. <br>" << endl;*/ <<" limited, noisy data. IEEE Trans Signal Processing,"
<<" 48(7):2083-2091, 2000<br>"<<endl;
report << "[Tipping 1999] M.E. Tipping and C.M.Bishop."
<<" Probabilistic Principal component analysis. J Royal"
<<" Statistical Society B, 61(3), 1999. <br>" << endl;
report << "[Beckmann 2001] C.F. Beckmann, J.A. Noble and"
<<" S.M. Smith. Investigating the intrinsic dimensionality"
<<" of FMRI data for ICA. In Seventh Int. Conf. on Functional"
<<" Mapping of the Human Brain, 2001. <br>" << endl;
if(opts.pca_est.value() == string("lap")) if(opts.pca_est.value() == string("lap"))
report << "[Minka 2000] T. Minka. Automatic choice of dimensionality for PCA. Technical Report 514, MIT Media Lab Vision and Modeling Group, 2000. <BR>"<< endl; report << "[Minka 2000] T. Minka. Automatic choice of"
<<" dimensionality for PCA. Technical Report 514, MIT"
<<" Media Lab Vision and Modeling Group, 2000. <BR>"<< endl;
else else
if(opts.pca_est.value() == string("bic")) if(opts.pca_est.value() == string("bic"))
report << "[Kass 1995] R.E. Kass and A. E. Raftery. Bayes factors. Journal of the American Statistical Association, 90:733-795, 1995 <br>" << endl; report << "[Kass 1995] R.E. Kass and A. E. Raftery. Bayes"
<<" factors. Journal of the American Statistical"
<<" Association, 90:733-795, 1995 <br>" << endl;
else else
if(opts.pca_est.value() == string("mdl")) if(opts.pca_est.value() == string("mdl"))
report << "[Rissanen 1978]. J. Rissanen. Modelling by shortest data description. Automatica, 14:465-471, 1978. <br>" << endl; report << "[Rissanen 1978]. J. Rissanen. Modelling by"
<<" shortest data description. Automatica,"
<<" 14:465-471, 1978. <br>" << endl;
else else
if(opts.pca_est.value() == string("aic")) if(opts.pca_est.value() == string("aic"))
report << "[Akaike 1974]. H. Akaike. A new look at statistical model identification. IEEE Transactions on Automatic Control, 19:716-723, 1974. <br>" << endl; report << "[Akaike 1974]. H. Akaike. A new look at"
<<" statistical model identification. IEEE Transactions"
<<" on Automatic Control, 19:716-723, 1974. <br>" << endl;
else else
report << "[Minka 2000]. T. Minka. Automatic choice of dimensionality for PCA. Technical Report 514, MIT Media Lab Vision and Modeling Group, 2000. <BR>" << endl; report << "[Minka 2000]. T. Minka. Automatic choice of"
<<" dimensionality for PCA. Technical Report 514, MIT"
<<" Media Lab Vision and Modeling Group, 2000. <BR>" << endl;
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment