diff --git a/meldata.cc b/meldata.cc
new file mode 100644
index 0000000000000000000000000000000000000000..01da5935cd6776a05525911fcd4aae8489c07aa0
--- /dev/null
+++ b/meldata.cc
@@ -0,0 +1,500 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ meldata.cc - data handler / container class
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "newimageall.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melodic.h"
+#include "log.h"
+#include <time.h>
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+
+ void MelodicData::setup()
+ {
+ {
+ volume4D<float> RawData;
+ message("Reading data file " << opts.inputfname.value() << " ... ");
+ read_volume4D(RawData,opts.inputfname.value(),tempInfo);
+ message(" done" << endl);
+ for(int ctr=1; ctr<=opts.dummy.value(); ctr++){
+ RawData.deletevolume(ctr);
+ }
+
+ // calculate a Mean image and save it
+ Mean = meanvol(RawData);
+ tmpnam(Mean_fname); // generate a tmp name
+ save_volume(Mean,Mean_fname);
+ create_mask(RawData, Mask);
+ if(Mask.xsize()==RawData.xsize() &&
+ Mask.ysize()==RawData.ysize() &&
+ Mask.zsize()==RawData.zsize())
+ Data = RawData.matrix(Mask);
+ else{
+ cerr << "ERROR:: mask and data have different dimensions \n\n";
+ exit(2);
+ }
+
+
+ // clean /tmp
+ char callRMstr[1000];
+ ostrstream osc(callRMstr,1000);
+ osc << "rm " << string(Mean_fname) <<"* " << '\0';
+ system(callRMstr);
+
+ //mask out constant voxels
+ message("Excluding voxels with constant value " << endl);
+ Matrix DStDev=stdev(Data);
+ volume4D<float> tmpMask;
+ tmpMask.setmatrix(DStDev,Mask);
+
+ float tMmax;
+ volume<float> tmpMask2;
+ tmpMask2 = tmpMask[0];
+ tMmax = tmpMask2.max();
+ double st_mean = DStDev.Sum()/DStDev.Ncols();
+ double st_std = stdev(DStDev.t()).AsScalar();
+
+ Mask = binarise(tmpMask2,(float) max((float) st_mean-3*st_std,
+ (float) 0.01*st_mean),tMmax);
+ Data = RawData.matrix(Mask);
+ }
+
+ message("Data size : " << Data.Nrows() << " x " << Data.Ncols() <<endl);
+
+ {// remove mean volume
+ // if((opts.remove_meanvol.value()||opts.varnorm.value())){
+ message(string("Removing mean image ... "));
+ meanR=mean(Data);
+ Data=remmean(Data);
+ message("done" << endl);
+ //}else{
+ // meanR=zeros(1,Data.Ncols());
+ //}
+ }
+ {//switch dimension in case temporal ICA is required
+ // if(opts.temporal.value()){
+ // if(opts.remove_meanvol.value()){
+ // message(string("Remove mean time course for temporal ICA ... "));
+ // meanC=meanR;
+ // meanR=mean(Data);
+ // Data=remmean(Data,2);
+ // message("done"<<endl);
+ // }
+ // Data = Data.t();
+ // }
+ }
+ {// remove mean time course
+ meanC=mean(Data,2);
+ Data=remmean(Data,2);
+ }
+ {// variance-normalize the data
+ DiagonalMatrix tmpD;Matrix tmpE;
+ message(string("Estimating data covariance ... "));
+ SymmetricMatrix Corr;
+ Corr = cov(Data.t());
+ EigenValues(Corr,tmpD,tmpE);
+
+ Matrix RE; DiagonalMatrix RD;
+ //RE = tmpE.Columns(2,Corr.Ncols());
+ RE = tmpE;
+ //RD << abs(tmpD.SymSubMatrix(2,Corr.Ncols()));
+ RD << abs(tmpD);
+ Matrix tmpWhite;Matrix tmpDeWhite;
+ tmpWhite = sqrt(abs(RD.i()))*RE.t();
+ tmpDeWhite = RE*sqrt(RD);
+ message("done"<<endl);
+ if(opts.varnorm.value())
+ message(string("Perform variance-normalisation ... "));
+ Matrix WS;
+ WS = tmpWhite * Data;
+ for(int ctr1 =1; ctr1<=WS.Nrows(); ctr1++){
+ for(int ctr2 =1; ctr2<=WS.Ncols(); ctr2++){
+ if(abs(WS(ctr1,ctr2))<3.1)
+ WS(ctr1,ctr2)=0.0;
+ }
+ }
+
+ stdDevi = pow(stdev(Data - tmpDeWhite*WS),-1);
+ Data = Data + meanC*ones(1,Data.Ncols());
+ }
+
+ DataVN = SP(Data,ones(Data.Nrows(),1)*stdDevi);
+ if(opts.varnorm.value()){
+ Data = DataVN;
+ message("done"<<endl);
+ }
+ {//remove row mean
+ if(opts.temporal.value()){
+ Data=remmean(Data,2);
+ }else{
+ message(string("Removing mean time course ... "));
+ meanC=mean(Data,2);
+ Data=remmean(Data,2);
+ message("done"<<endl);
+ }
+ }
+
+ if(opts.segment.value().length()>0){
+ create_RXweight();
+ }
+
+ //save the mask
+ save_volume(Mask,logger.appendDir("mask"));
+
+ //seed the random number generator
+ double tmptime = time(NULL);
+ srand((unsigned int) tmptime/Data.Ncols()*Data.Nrows());
+
+ } // void setup()
+
+ void MelodicData::save()
+ {
+
+ //check for temporal ICA
+ if(opts.temporal.value()){
+ Matrix tmpIC = mixMatrix.t();
+ mixMatrix=IC.t();
+ IC=tmpIC;
+ }
+
+ message("Writing results to : " << endl);
+
+ //Output IC
+ if((IC.Storage()>0)&&(opts.output_origIC.value()) ){
+ volume4D<float> tempVol;
+ tempVol.setmatrix(IC,Mask);
+ //strncpy(tempInfo.header.hist.aux_file,"render3",24);
+ save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
+ + "_origIC"),tempInfo);
+ message(" " << logger.appendDir(opts.outputfname.value() + "_oIC") <<endl);
+ }
+
+ //Output IC -- adjusted for noise
+ if(IC.Storage()>0){
+ volume4D<float> tempVol;
+ // volumeinfo tempInfo;
+ // read_volume4D(tempVol,opts.inputfname.value(),tempInfo);
+
+ //Matrix ICadjust = IC;
+ if(opts.temporal.value()){
+ Data=Data.t();
+ }
+
+ Matrix ICadjust;
+ stdNoisei = pow(stdev(Data - mixMatrix * IC),-1);
+ ICadjust = SP(IC,ones(IC.Nrows(),1)*stdNoisei);
+
+ tempVol.setmatrix(ICadjust,Mask);
+ //strncpy(tempInfo.header.hist.aux_file,"render3",24);
+ save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
+ + "_IC"),tempInfo);
+ message(" " << logger.appendDir(opts.outputfname.value() + "_IC") <<endl);
+ }
+
+ //Output mixMatrix
+ if(mixMatrix.Storage()>0){
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_mix"),
+ mixMatrix);
+ mixFFT=calc_FFT(mixMatrix);
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_FTmix"),
+ mixFFT);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_mix") <<endl);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_FTmix") <<endl);
+ }
+
+ //Output unmixMatrix
+ if(opts.output_unmix.value() && unmixMatrix.Storage()>0){
+ write_ascii_matrix(opts.outputfname.value() + "_unmix",unmixMatrix);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_unmix") <<endl);
+ }
+
+ //Output Mask
+ message(" "<< logger.appendDir("mask") <<endl);
+
+ //Output mean
+ if(opts.output_mean.value() && meanC.Storage()>0 && meanR.Storage()>0){
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_meanR"),
+ meanR);
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_meanC"),
+ meanC);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_meanR") <<endl);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_meanC") <<endl);
+ }
+
+ //Output white
+ if(opts.output_white.value() && whiteMatrix.Storage()>0&&
+ dewhiteMatrix.Storage()>0){
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_white"),
+ whiteMatrix);
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_dewhite"),dewhiteMatrix);
+ Matrix tmp;
+ tmp=calc_FFT(dewhiteMatrix);
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_FTdewhite"),tmp);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_white") <<endl);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_dewhite") <<endl);
+ }
+
+ //Output PCA
+ if(opts.output_pca.value() && pcaD.Storage()>0&&pcaE.Storage()>0){
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_pcaE"),
+ pcaE);
+ write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_pcaD"),
+ (Matrix) diag(pcaD));
+ volume4D<float> tempVol;
+ //volumeinfo tempInfo;
+ //read_volume4D(tempVol,opts.inputfname.value(),tempInfo);
+
+ Matrix ICadjust = IC;
+ if(opts.temporal.value()){
+ Data=Data.t();
+ }
+
+ Matrix PCAmaps;
+ PCAmaps = whiteMatrix * Data;
+
+ tempVol.setmatrix(PCAmaps,Mask);
+ //strncpy(tempInfo.header.hist.aux_file,"render3",24);
+ save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
+ + "_pca"),tempInfo);
+ message(" " <<
+ logger.appendDir(opts.outputfname.value() + "_pca") <<endl);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_pcaE") <<endl);
+ message(" "<<
+ logger.appendDir(opts.outputfname.value() + "_pcaD") <<endl);
+ }
+ } //void save()
+
+ int MelodicData::remove_components()
+ {
+ message("Reading " << opts.filtermix.value() << endl)
+ mixMatrix = read_ascii_matrix(opts.filtermix.value());
+ if (mixMatrix.Storage()<=0) {
+ cerr <<" Please specify the mixing matrix correctly" << endl;
+ exit(2);
+ }
+
+ unmixMatrix = pinv(mixMatrix);
+ IC = unmixMatrix * Data;
+
+ string tmpstr;
+ tmpstr = opts.filter.value();
+
+ Matrix noiseMix;
+ Matrix noiseIC;
+
+ int ctr=0;
+ char *p;
+ char t[1024];
+ const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
+
+ message("Filtering the data...");
+ strcpy(t, tmpstr.c_str());
+ p=strtok(t,discard);
+ ctr = atoi(p);
+ if(ctr>0 && ctr<=mixMatrix.Ncols()){
+ message(" "<< ctr );
+ noiseMix = mixMatrix.Column(ctr);
+ noiseIC = IC.Row(ctr).t();
+ }else{
+ cerr << endl<< "component number "<<ctr<<" does not exist" << endl;
+ }
+
+ do{
+ p=strtok(NULL,discard);
+ if(p){
+ ctr = atoi(p);
+
+ if(ctr>0 && ctr<=mixMatrix.Ncols()){
+ message(" "<<ctr);
+ noiseMix |= mixMatrix.Column(ctr);
+ noiseIC |= IC.Row(ctr).t();
+ }
+ else{
+ cerr << endl<< "component number "<<ctr<<" does not exist" << endl;
+ }
+ }
+ }while(p);
+ message(endl);
+ Matrix newData;
+ newData = Data - noiseMix * noiseIC.t();
+
+ //cerr << newData.Nrows() << " x " << newData.Ncols() << endl;
+ //cerr << meanC.Nrows() << " x " << meanC.Ncols() << endl;
+ //cerr << meanR.Nrows() << " x " << meanR.Ncols() << endl;
+ newData = newData + meanC*ones(1,newData.Ncols());
+ newData = newData + ones(newData.Nrows(),1)*meanR;
+
+ volume4D<float> tmp;
+ read_volume4D(tmp,opts.inputfname.value());
+ tmp.setmatrix(newData,Mask);
+ save_volume4D(tmp,logger.appendDir(opts.outputfname.value() + "_ICAfiltered"));
+
+ return 0;
+ } // int remove_components()
+
+ void MelodicData::create_RXweight()
+ {
+ message("Reading the weights for the covariance R_X from file "<< opts.segment.value() << endl);
+
+ volume4D<float> tmpRX;
+ read_volume4D(tmpRX,opts.segment.value());
+ RXweight = tmpRX.matrix(Mask);
+ }
+
+ void MelodicData::create_mask(volume4D<float> &theData,
+ volume<float> &theMask)
+ {
+ if(opts.use_mask.value() && opts.maskfname.value().size()>0){ // mask provided
+ read_volume(theMask,opts.maskfname.value());
+ message("Mask provided : " << opts.maskfname.value()<<endl);
+ }
+ else{
+ if(opts.perf_bet.value() && opts.use_mask.value()){ //use BET
+ message("Create mask ... ");
+ // set up all strings
+ string BET_outputfname = string(Mean_fname)+"_brain";
+
+ string BET_path = opts.binpath + "bet";
+ string BET_optarg = "-m -f 0.4"; // see man bet
+ string Mask_fname = BET_outputfname+"_mask.hdr";
+
+ // Setup external call to BET:
+
+ char callBETstr[1000];
+ ostrstream osc(callBETstr,1000);
+ osc << BET_path << " " << Mean_fname << " "
+ << BET_outputfname << " " << BET_optarg << " > /dev/null " << '\0';
+
+ message(" Calling BET: " << callBETstr << endl);
+ system(callBETstr);
+
+ // read back the Mask file
+ read_volume(theMask,Mask_fname);
+
+ message("done" << endl);
+ }
+ else{
+ if(opts.use_mask.value()){ //just threshold the Mean
+ message("Create mask ... ");
+ float Mmin, Mmax, Mtmp;
+ Mmin = Mean.min(); Mmax = Mean.max();
+ theMask = binarise(Mean,Mmin + opts.threshold.value()* (Mmax-Mmin),Mmax);
+ Mtmp = Mmin + opts.threshold.value()* (Mmax-Mmin);
+ message("done" << endl);
+ }
+ else{ //well, don't threshold then
+ theMask = Mean;
+ theMask = 1.0;
+ }
+ }
+ }
+ if(opts.remove_endslices.value()){
+ // just in case mc introduced something nasty
+ message(" Deleting end slices" << endl);
+ for(int ctr1=theMask.miny(); ctr1<=theMask.maxy(); ctr1++){
+ for(int ctr2=theMask.minx(); ctr2<=theMask.maxx(); ctr2++){
+ theMask(ctr2,ctr1,Mask.minz()) = 0.0;
+ theMask(ctr2,ctr1,Mask.maxz()) = 0.0;
+ }
+ }
+ }
+ } //void create_mask()
+
+
+ void MelodicData::status(const string &txt)
+ {
+ cout << "MelodicData Object " << txt << endl;
+ if(Data.Storage()>0){cout << "Data: " << Data.Nrows() <<"x" << Data.Ncols() << endl;}else{cout << "Data empty " <<endl;}
+ if(pcaE.Storage()>0){cout << "pcaE: " << pcaE.Nrows() <<"x" << pcaE.Ncols() << endl;}else{cout << "pcaE empty " <<endl;}
+ if(pcaD.Storage()>0){cout << "pcaD: " << pcaD.Nrows() <<"x" << pcaD.Ncols() << endl;}else{cout << "pcaD empty " <<endl;}
+ if(whiteMatrix.Storage()>0){cout << "white: " << whiteMatrix.Nrows() <<"x" << whiteMatrix.Ncols() << endl;}else{cout << "white empty " <<endl;}
+ if(dewhiteMatrix.Storage()>0){cout << "dewhite: " << dewhiteMatrix.Nrows() <<"x" << dewhiteMatrix.Ncols() << endl;}else{cout << "dewhite empty " <<endl;}
+ if(mixMatrix.Storage()>0){cout << "mix: " << mixMatrix.Nrows() <<"x" << mixMatrix.Ncols() << endl;}else{cout << "mix empty " <<endl;}
+ if(unmixMatrix.Storage()>0){cout << "unmix: " << unmixMatrix.Nrows() <<"x" << unmixMatrix.Ncols() << endl;}else{cout << "unmix empty " <<endl;}
+ if(IC.Storage()>0){cout << "IC: " << IC.Nrows() <<"x" << IC.Ncols() << endl;}else{cout << "IC empty " <<endl;}
+
+ } //void status()
+
+ Matrix MelodicData::calc_FFT(const Matrix& Mat)
+ {
+ Matrix res;
+ for(int ctr=1; ctr <= Mat.Ncols(); ctr++)
+ {
+ ColumnVector tmpCol;
+ tmpCol=Mat.Column(ctr);
+ ColumnVector FtmpCol_real;
+ ColumnVector FtmpCol_imag;
+ ColumnVector tmpPow;
+ if(tmpCol.Nrows()%2 != 0){
+ Matrix empty(1,1); empty=0;
+ tmpCol &= empty;}
+ RealFFT(tmpCol,FtmpCol_real,FtmpCol_imag);
+ tmpPow = pow(FtmpCol_real,2)+pow(FtmpCol_imag,2);
+ tmpPow = tmpPow.Rows(2,tmpPow.Nrows());
+ if(opts.logPower.value()) tmpPow = log(tmpPow);
+ if(res.Storage()==0){res= tmpPow;}else{res|=tmpPow;}
+ }
+ return res;
+ } //Matrix calc_FFT()
+
+ Matrix MelodicData::smoothColumns(const Matrix &inp)
+ {
+ Matrix temp(inp);
+ int ctr1 = temp.Nrows();
+ Matrix temp2(temp);
+ temp2=0;
+
+ temp = temp.Row(4) & temp.Row(3) & temp.Row(2) & temp & temp.Row(ctr1-1)
+ & temp.Row(ctr1-2) &temp.Row(ctr1-3);
+
+ double kern[] ={0.0045 , 0.055, 0.25, 0.4, 0.25, 0.055, 0.0045};
+ double fac = 0.9090909;
+
+ // Matrix FFTinp;
+ //FFTinp = calc_FFT(inp);
+ //write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_FT"),
+ // FFTinp); FFTinp = calc_FFT(inp);
+ //write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_Sinp"),
+ // inp);
+
+ for(int cc=1;cc<=temp2.Ncols();cc++){
+ // double all = FFTinp.Column(cc).Rows(1,30).SumAbsoluteValue() / FFTinp.Column(cc).SumAbsoluteValue();
+ // if( all > 0.5){
+ for(int cr=1;cr<=temp2.Nrows();cr++){
+ temp2(cr,cc) = fac*( kern[0] * temp(cr,cc) + kern[1] * temp(cr+1,cc) +
+ kern[2] * temp(cr+2,cc) + kern[3] * temp(cr+3,cc) +
+ kern[4] * temp(cr+4,cc) + kern[5] * temp(cr+5,cc) +
+ kern[6] * temp(cr+6,cc));
+ }//}
+ //else{
+ // for(int cr=1;cr<=temp2.Nrows();cr++){
+ // temp2(cr,cc) = temp(cr+3,cc);
+ // }
+ //}
+ }
+//write_ascii_matrix(logger.appendDir(opts.outputfname.value() + "_Sout"),
+// temp2);
+ return temp2;
+ }
+}
+
diff --git a/meldata.h b/meldata.h
new file mode 100644
index 0000000000000000000000000000000000000000..8eadddf3d3ecadedfc5ad4216bbda6152631beb7
--- /dev/null
+++ b/meldata.h
@@ -0,0 +1,172 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ meldata.h - data container class
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+
+#ifndef __MELODICDATA_h
+#define __MELODICDATA_h
+
+#include "newimageall.h"
+#include "log.h"
+#include "meloptions.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+
+ class MelodicData
+ {
+ public:
+
+ //constructor
+ MelodicData(MelodicOptions &popts, Log &plogger):
+ opts(popts),
+ logger(plogger)
+ {
+ }
+
+ void save();
+
+ inline void save4D(Matrix what, string fname){
+ volume4D<float> tempVol;
+ tempVol.setmatrix(what,Mask);
+ save_volume4D(tempVol,logger.appendDir(fname),tempInfo);
+ }
+
+ inline void saveascii(Matrix what, string fname){
+ write_ascii_matrix(logger.appendDir(fname),what); }
+
+ int remove_components();
+ void setup();
+ void status(const string &txt);
+
+ inline Matrix& get_pcaE() {return pcaE;}
+ inline void set_pcaE(Matrix& Arg) {pcaE = Arg;}
+
+ inline DiagonalMatrix& get_pcaD() {return pcaD;}
+ inline void set_pcaD(DiagonalMatrix& Arg) {pcaD = Arg;}
+
+ inline Matrix& get_data() {return Data;}
+ inline void set_data(Matrix& Arg) {Data = Arg;}
+
+ inline Matrix& get_IC() {return IC;}
+ inline void set_IC(Matrix& Arg) {IC = Arg;}
+
+ inline Matrix& get_white() {return whiteMatrix;}
+ inline void set_white(Matrix& Arg) {whiteMatrix = Arg;}
+
+ inline Matrix& get_dewhite() {return dewhiteMatrix;}
+ inline void set_dewhite(Matrix& Arg) {dewhiteMatrix = Arg;}
+
+ inline Matrix& get_meanC() {return meanC;}
+ inline Matrix& get_meanR() {return meanR;}
+
+ inline Matrix& get_stdDevi() {return stdDevi;}
+ inline void set_stdDevi(Matrix& Arg) {stdDevi = Arg;}
+
+ inline Matrix& get_mix() {return mixMatrix;}
+ inline void set_mix(Matrix& Arg) {mixMatrix = Arg;}
+
+ inline Matrix& get_fmix() {return mixFFT;}
+ inline void set_fmix(Matrix& Arg) {mixFFT = Arg;}
+
+ inline Matrix& get_unmix() {return unmixMatrix;}
+ inline void set_unmix(Matrix& Arg) {unmixMatrix = Arg;}
+
+ inline volume<float>& get_mask() {return Mask;}
+ inline void set_mask(volume<float>& Arg) {Mask = Arg;}
+
+ inline volume<float>& get_mean() {return Mean;}
+ inline void set_mean(volume<float>& Arg) {Mean = Arg;}
+
+ inline Matrix& get_Data() {return Data;}
+ inline void set_Data(Matrix& Arg) {Data = Arg;}
+
+ inline Matrix& get_DataVN() {return DataVN;}
+ inline void set_DataVN(Matrix& Arg) {DataVN = Arg;}
+
+ inline Matrix& get_RXweight() {return RXweight;}
+ inline void set_RXweight(Matrix& Arg) {RXweight = Arg;}
+
+ inline Matrix& get_EVP() {return EVP;}
+ inline void set_EVP(Matrix& Arg) {EVP = Arg;}
+
+ inline Matrix& get_EV() {return EV;}
+ inline void set_EV(Matrix& Arg) {EV = Arg;}
+
+ inline Matrix& get_PPCA() {return PPCA;}
+ inline void set_PPCA(Matrix& Arg) {PPCA = Arg;}
+
+ inline Matrix& get_stdNoisei() {return stdNoisei;}
+ inline void set_stdNoisei(Matrix& Arg) {stdNoisei = Arg;}
+
+ inline int data_dim() {return Data.Nrows();}
+ inline int data_samples() {return Data.Ncols();}
+
+ inline float get_resels() {return Resels;}
+ inline void set_resels(float& Arg) {Resels = Arg;}
+
+ Matrix smoothColumns(const Matrix &inp);
+
+ inline void flipres(int num){
+ IC.Row(num) = -IC.Row(num);
+ mixMatrix.Column(num) = -mixMatrix.Column(num);
+ mixFFT=calc_FFT(mixMatrix);
+ unmixMatrix = pinv(mixMatrix);
+ }
+
+
+ inline void varnorm(){
+ Data = SP(Data,ones(Data.Nrows(),1)*stdDevi);
+ }
+
+ volumeinfo tempInfo;
+
+ Matrix calc_FFT(const Matrix& Mat);
+
+ private:
+ MelodicOptions &opts;
+ Log &logger;
+ // static MelodicData* melodat;
+
+ Matrix pcaE;
+ DiagonalMatrix pcaD;
+ Matrix whiteMatrix;
+ Matrix dewhiteMatrix;
+ Matrix meanC;
+ Matrix meanR;
+ Matrix stdDevi;
+ Matrix RXweight;
+ Matrix mixMatrix;
+ Matrix unmixMatrix;
+ Matrix mixFFT;
+ Matrix IC;
+ Matrix EVP;
+ Matrix EV;
+ Matrix stdNoisei;
+ volume<float> Mask;
+ volume<float> Mean;
+ Matrix Data;
+ Matrix DataVN;
+ Matrix PPCA;
+
+ float Resels;
+
+ char Mean_fname[1000];
+
+ void create_mask(volume4D<float> &theData, volume<float> &theMask);
+ // void remove_comp(Matrix& Mat, const int numIC);
+ void create_RXweight();
+
+ };
+}
+
+#endif
diff --git a/melgmix.h b/melgmix.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b4755e2a5665ca58011550db6b7a325f751927e
--- /dev/null
+++ b/melgmix.h
@@ -0,0 +1,195 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melgmix.h - class for Gaussian/Gamma Mixture Model
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#ifndef __MELGMIX_h
+#define __MELGMIX_h
+
+#include "newimageall.h"
+#include "log.h"
+#include "melodic.h"
+#include "options.h"
+#include "meloptions.h"
+//#include "melreport.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+
+ class MelGMix
+ {
+ public:
+
+ //constructor
+ /* MelGMix(MelodicOptions &popts, Log &plogger): */
+/* opts(popts), */
+/* logger(plogger), */
+/* mainhtml() */
+/* { */
+/* } */
+
+ MelGMix(MelodicOptions &popts, Log &plogger):
+ opts(popts),
+ logger(plogger)
+ {
+ }
+
+ ~MelGMix() {
+ //mainhtml << endl << "<hr></CENTER></BODY></HTML>" << endl;
+ }
+
+ void save();
+
+ void setup(const RowVector& dat, volumeinfo inf, const string dirname,
+ int here, volume<float> themask,
+ volume<float> themean, int num_mix = 3,
+ float eps = 0.0, bool fixdim = false);
+
+
+ void gmmfit();
+ void ggmfit();
+
+ inline void fit(string mtype = string("GGM"))
+ {
+ mmtype = mtype;
+ if(mmtype==string("GGM"))
+ this->ggmfit();
+ else
+ this->gmmfit();
+ }
+
+ inline Matrix threshold(string levels)
+ {return this->threshold(data, levels);}
+ inline Matrix threshold(RowVector& levels)
+ {return this->threshold(data, levels);}
+ Matrix threshold(const RowVector& dat, Matrix& levels);
+ Matrix threshold(const RowVector& dat, string levels);
+
+ void status(const string &txt);
+
+ inline RowVector& get_means() {return means;}
+ inline void set_means(RowVector& Arg) {means = Arg;}
+
+ inline RowVector& get_vars() {return vars;}
+ inline void set_vars(RowVector& Arg) {vars = Arg;}
+
+ inline RowVector& get_pi() {return props;}
+ inline void set_pi(RowVector& Arg) {props = Arg;}
+
+ inline RowVector& get_data() {return data;}
+ inline void set_data(RowVector& Arg) {data = Arg;}
+
+ inline RowVector& get_prob() {return probmap;}
+
+ inline float get_eps() {return epsilon;}
+ inline void set_eps(float Arg) {epsilon = Arg;}
+
+ inline Matrix& get_threshmaps() {return threshmaps;}
+ inline void set_threshmaps(Matrix& Arg) {threshmaps = Arg;}
+
+ inline bool isfitted(){return fitted;}
+
+ inline int mixtures(){return nummix;}
+
+ inline string get_type() { return mmtype;}
+ inline void set_type(string Arg) { mmtype = Arg;}
+
+ inline string get_prefix() { return prefix;}
+ inline void set_prefix(string Arg) { prefix = Arg;}
+
+ inline RowVector get_probmap() {return probmap;}
+
+ inline float get_offset() {return offset;}
+ inline void set_offset(float Arg) {offset = Arg;}
+
+ inline void flipres(int num){
+ means = -means;
+ data = -data;
+ threshmaps = -threshmaps;
+ if(mmtype=="GGM"){
+ float tmp;
+ tmp= means(2);means(2)= means(3);means(3)=tmp;
+ tmp=vars(2);vars(2)=vars(3);vars(3)=tmp;
+ tmp=props(2);props(2)=props(3);props(3)=tmp;
+ }
+ }
+
+ void create_rep();
+
+ inline void add_infstr(string what){
+ threshinfo.push_back(what);
+ }
+
+ inline string get_infstr(int num){
+ if((threshinfo.size()<(unsigned int)(num-1))||(num<1))
+ return string("");
+ else
+ return threshinfo[num-1];
+ }
+
+ inline int size_infstr(){
+ return threshinfo.size();
+ }
+
+ inline void clear_infstr(){
+ threshinfo.clear();
+ }
+
+ private:
+ MelodicOptions &opts;
+ Log &logger; //global log file
+
+ //Log mainhtml;
+
+ void gmmupdate();
+ float gmmevidence();
+ void gmmreducemm();
+ void add_params(Matrix& mu, Matrix& sig, Matrix& pi,
+ float logLH, float MDL, float Evi, bool advance = false);
+ void get_params(int index, Matrix& mu, Matrix& sig, Matrix& pi,
+ float logLH, float MDL, float Evi);
+
+ Matrix Params;
+ Matrix threshmaps;
+
+ RowVector means;
+ RowVector vars;
+ RowVector props;
+ RowVector data;
+ RowVector probmap;
+
+ volume<float> Mean;
+ volume<float> Mask;
+
+ float epsilon;
+ float logprobY;
+ float MDL;
+ float Evi;
+ float offset;
+
+ int nummix;
+ int numdata;
+ int cnumber;
+
+ bool fitted;
+ bool fixdim;
+
+ string prefix;
+ string mmtype;
+ string dirname;
+
+ volumeinfo bginfo;
+ vector<string> threshinfo;
+
+ };
+}
+
+#endif
diff --git a/melica.cc b/melica.cc
new file mode 100644
index 0000000000000000000000000000000000000000..1f355fb3dbeb0f91400cb6d3751e474d7a54b3e2
--- /dev/null
+++ b/melica.cc
@@ -0,0 +1,408 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melica.cc - ICA estimation
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <stdlib.h>
+#include "newimageall.h"
+#include "log.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melodic.h"
+#include "newmatap.h"
+#include "newmatio.h"
+#include "melica.h"
+#include "melpca.h"
+#include "miscprob.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+
+ void MelodicICA::ica_fastica_symm(const Matrix &Data)
+ {
+ // based on Aapo Hyvärinen's fastica method
+ // see www.cis.hut.fi/projects/ica/fastica/
+
+ //initialize matrices
+ Matrix redUMM_old;
+ Matrix tmpU;
+
+ //srand((unsigned int)timer(NULL));
+
+ redUMM = melodat.get_white()*
+ unifrnd(melodat.get_white().Ncols(),dim); // got to start somewhere
+
+ if(opts.guessfname.value().size()>1){
+ message(" Use columns in " << opts.guessfname.value()
+ << " as initial values for the mixing matrix " <<endl);
+ Matrix guess ;
+ guess = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
+ redUMM.Columns(1,guess.Ncols()) = guess;
+ }
+
+ symm_orth(redUMM);
+
+ int itt_ctr=1;
+ double minAbsSin = 1.0;
+ do{ // da loop!!!
+
+ redUMM_old = redUMM;
+
+ //calculate IC estimates
+ tmpU = Data.t() * redUMM;
+
+ //update redUMM depending on nonlinearity
+ if(opts.nonlinearity.value()=="pow4"){
+ redUMM = (Data * pow(tmpU,3.0)) / samples - 3 * redUMM;
+ }
+ if(opts.nonlinearity.value()=="pow3"){
+ tmpU /= opts.nlconst1.value();
+ redUMM = 3 * (Data * pow(tmpU,2.0)) / samples -
+ (SP(ones(dim,1)*sum(tmpU,1),redUMM))/ samples;
+ }
+ if(opts.nonlinearity.value()=="tanh"){
+ Matrix hyptanh;
+ hyptanh = tanh(opts.nlconst1.value()*tmpU);
+ redUMM = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
+ sum(1-pow(hyptanh,2),1),redUMM))/samples;
+ }
+ if(opts.nonlinearity.value()=="gauss"){
+ Matrix tmpUsq;
+ Matrix tmpU2;
+ Matrix gauss;
+ Matrix dgauss;
+ tmpUsq = pow(tmpU,2);
+ tmpU2 = exp(-(opts.nlconst2.value()/2) * tmpUsq);
+ gauss = SP(tmpU,tmpU2);
+ dgauss = SP(1-opts.nlconst2.value()*tmpUsq,tmpU2);
+ redUMM = (Data * gauss - SP(ones(dim,1)*
+ sum(dgauss,1),redUMM))/samples;
+ }
+
+ // orthogonalize the unmixing-matrix
+ symm_orth(redUMM);
+
+ //termination condition : angle between old and new < epsilon
+ minAbsSin = 1 - diag(abs(redUMM.t()*redUMM_old)).Minimum();
+ message(" Step no. " << itt_ctr << " change : " << minAbsSin << endl);
+ if(abs(minAbsSin) < opts.epsilon.value()){ break;}
+
+ itt_ctr++;
+ } while(itt_ctr < opts.maxNumItt.value());
+
+ if(itt_ctr>=opts.maxNumItt.value()){
+ cerr << " No convergence after " << itt_ctr <<" steps "<<endl;
+ } else {
+ message(" Convergence after " << itt_ctr <<" steps " << endl << endl);
+ no_convergence = false;
+ {Matrix temp(melodat.get_dewhite() * redUMM);
+ melodat.set_mix(temp);}
+ {Matrix temp(redUMM.t()*melodat.get_white());
+ melodat.set_unmix(temp);}
+ }
+ }
+
+
+ void MelodicICA::ica_fastica_defl(const Matrix &Data)
+ {
+ 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()*
+ unifrnd(melodat.get_white().Ncols(),opts.numICs.value());
+ int guesses=0;
+ if(opts.guessfname.value().size()>1){
+ message(" Use columns in " << opts.guessfname.value() << " as initial values for the mixing matrix " <<endl);
+ Matrix guess;
+ guess = melodat.get_white()*read_ascii_matrix(opts.guessfname.value());
+ guesses = guess.Ncols();
+ redUMM.Columns(1,guesses) = guess;
+ }
+
+ int ctrIC = 1;
+ int numRestart = 0;
+
+ while(ctrIC<=opts.numICs.value()){
+
+ message(" Extracting IC " << ctrIC << " ... ");
+ ColumnVector w;
+ ColumnVector w_old;
+ ColumnVector tmpU;
+ if(ctrIC <= guesses){
+ w = redUMM.Column(ctrIC);}
+ else{
+ w = unifrnd(dim,1);}
+
+ w = w - redUMM * redUMM.t() * w;
+ w = w / norm2(w);
+ int itt_ctr = 1;
+
+ do{
+ w_old = w;
+
+ tmpU = Data.t() * w;
+
+ if(opts.nonlinearity.value()=="pow4"){
+ w = (Data * pow(tmpU,3.0)) / samples - 3 * w;
+ }
+ if(opts.nonlinearity.value()=="tanh"){
+ ColumnVector hyptanh;
+ hyptanh = tanh(opts.nlconst1.value()*tmpU);
+ w = (Data * hyptanh - opts.nlconst1.value()*SP(ones(dim,1)*
+ sum(1-pow(hyptanh,2),1),w))/samples;
+ }
+ if(opts.nonlinearity.value()=="pow3"){
+ tmpU /= opts.nlconst1.value();
+ w = 3*(Data * pow(tmpU,2.0)) / samples - 2*(SP(ones(dim,1)*
+ sum(tmpU,1),w))/samples;
+ }
+ if(opts.nonlinearity.value()=="gauss"){
+ ColumnVector tmpUsq;
+ ColumnVector tmpU2;
+ ColumnVector gauss;
+ ColumnVector dgauss;
+ tmpUsq = pow(tmpU,2);
+ tmpU2 = exp(-(opts.nlconst2.value()/2) * tmpUsq);
+ gauss = SP(tmpU,tmpU2);
+ dgauss = SP(1-opts.nlconst2.value()*tmpUsq,tmpU2);
+ w = (Data * gauss - SP(ones(dim,1)*
+ sum(dgauss,1),w))/samples;
+ }
+
+ // orthogonalize w
+ w = w - redUMM * redUMM.t() * w;
+ w = w / norm2(w);
+
+ //termination condition : angle between old and new < epsilon
+ if(norm2(w-w_old) < opts.epsilon.value() ||
+ norm2(w+w_old) < opts.epsilon.value()){break;}
+ //cout << norm2(w-w_old) << " " << norm2(w+w_old) << endl;
+ itt_ctr++;
+ } while(itt_ctr <= opts.maxNumItt.value());
+
+ if(itt_ctr<opts.maxNumItt.value()){
+ redUMM.Column(ctrIC) = w;
+ message(" estimated using " << itt_ctr << " iterations " << endl);
+ ctrIC++;
+ numRestart = 0;
+ }
+ else{
+ if(numRestart > opts.maxRestart.value()){
+ message(endl << " Estimation failed after "
+ << numRestart << " attempts "
+ << " giving up " << endl);
+ break;
+ }
+ else{
+ numRestart++;
+ message(endl <<" Estimation failed - restart "
+ << numRestart << endl);
+ }
+ }
+ }
+ if(numRestart < opts.maxRestart.value()){
+ no_convergence = false;
+ {Matrix temp(melodat.get_dewhite() * redUMM);
+ melodat.set_mix(temp);}
+ {Matrix temp(redUMM.t()*melodat.get_white());
+ melodat.set_unmix(temp);}
+ }
+ }
+
+
+ void MelodicICA::ica_maxent(const Matrix &Data)
+ {
+ cerr << " Not fully implemented yet " << endl;
+ // cout << "maxent" <<endl;
+ }
+
+
+ void MelodicICA::ica_jade(const Matrix &Data)
+ {
+ int dim_sym = (int) (dim*(dim+1))/2;
+ int num_CM = dim_sym;
+ Matrix CM;
+ Matrix R; R = Identity(dim);
+ Matrix Qij; Qij = zeros(dim);
+ Matrix Xim;
+ Matrix Xjm;
+ Matrix scale; scale = ones(dim,1)/samples;
+
+ for (int im =1; im <= dim; im++){
+ Xim = Data.Row(im);
+write_ascii_matrix("Xim",Data.Row(1));
+ //Qij = SP((scale * pow(Xim,2)),Data) * Data.t();//- R - 2*R.Column(im)*R.Column(im).t();
+ Qij = (pow(Xim,2)) * Data.t();//- R - 2*R.Column(im)*R.Column(im).t();
+ if(im==1){CM = Qij; write_ascii_matrix("CM",CM);exit(2);}else{CM |= Qij;}
+ for (int jm = 1; jm < im; jm++){
+ Xjm = Data.Row(jm);
+ Qij = (SP((scale * SP(Xim,Xjm)),Data) * Data.t() - R.Column(im)*R.Column(jm).t()
+ - R.Column(jm)*R.Column(im).t());
+ Qij *= sqrt(2);
+ CM |= Qij;
+ }
+ }
+
+ write_ascii_matrix("CM",CM);
+
+ Matrix redUMM; redUMM = Identity(dim);
+
+ bool exitloop = false;
+ int ctr_itt = 0;
+ int ctr_updates = 0;
+ Matrix Givens; Givens = zeros(2,num_CM);
+ Matrix Givens_ip; Givens_ip = zeros(2);
+ Matrix Givens_ro; Givens_ro = zeros(2);
+ double det_on, det_off;
+ double cos_theta, sin_theta, theta;
+
+ while(!exitloop && ctr_itt <= opts.maxNumItt.value()){
+ ctr_itt++;
+ cout << "loop" <<endl;
+ for(int ctr_p = 1; ctr_p < dim; ctr_p++){
+ for(int ctr_q = ctr_p+1; ctr_q <= dim; ctr_q++){
+
+ for(int ctr_i = 0; ctr_i < num_CM; ctr_i++){
+ int Ip = ctr_p + ctr_i * dim;
+ int Iq = ctr_q + ctr_i * dim;
+ Givens(1,ctr_i + 1) = CM(ctr_p,Ip) - CM(ctr_q,Iq);
+ Givens(2,ctr_i + 1) = CM(ctr_p,Iq) - CM(ctr_q,Ip);
+ }
+
+ Givens_ip = Givens * Givens.t();
+ det_on = Givens_ip(1,1) - Givens_ip(2,2);
+ det_off = Givens_ip(1,2) + Givens_ip(2,1);
+ theta = 0.5 * atan2(det_off, det_on + sqrt(det_on*det_on + det_off*det_off));
+
+ cout << theta << endl;
+
+ if(abs(theta) > opts.epsilon.value()){
+ ctr_updates++;
+ message(" Step No. "<< ctr_itt << " change : " << theta << endl);
+
+ //create Givens rotation matrix
+ cos_theta = cos(theta);
+ sin_theta = sin(theta);
+ Givens_ro(1,1) = cos_theta;
+ Givens_ro(1,2) = -sin_theta;
+ Givens_ro(2,1) = sin_theta;
+ Givens_ro(2,2) = cos_theta;
+
+ //update 2 entries of redUMM
+ {Matrix tmp_redUMM;
+ tmp_redUMM = redUMM.Column(ctr_p);
+ tmp_redUMM |= redUMM.Column(ctr_q);
+ tmp_redUMM *= Givens_ro;
+ redUMM.Column(ctr_p) = tmp_redUMM.Column(1);
+ redUMM.Column(ctr_q) = tmp_redUMM.Column(2);}
+
+ //update Cumulant matrix
+ {Matrix tmp_CM;
+ tmp_CM = CM.Row(ctr_p);
+ tmp_CM &= CM.Row(ctr_q);
+ tmp_CM = Givens_ro.t() * tmp_CM;
+ CM.Row(ctr_p) = tmp_CM.Row(1);
+ CM.Row(ctr_q) = tmp_CM.Row(2);}
+
+ //update Cumulant matrices
+ for(int ctr_i = 0; ctr_i < num_CM; ctr_i++){
+ int Ip = ctr_p + ctr_i * dim;
+ int Iq = ctr_q + ctr_i * dim;
+ CM.Column(Ip) = cos_theta*CM.Column(Ip)+sin_theta*CM.Column(Iq);
+ CM.Column(Iq) = cos_theta*CM.Column(Iq)-sin_theta*CM.Column(Ip);
+ }
+ }
+ else{
+ exitloop = true;
+ }
+ }
+ }
+ }//while loop
+ if(ctr_itt > opts.maxNumItt.value()){
+ cerr << " No convergence after " << ctr_itt <<" steps "<<endl;
+ } else {
+ message(" Convergence after " << ctr_itt <<" steps " << endl << endl);
+ no_convergence = false;
+ {Matrix temp(melodat.get_dewhite() * redUMM);
+ melodat.set_mix(temp);}
+ {Matrix temp(redUMM.t()*melodat.get_white());
+ melodat.set_unmix(temp);
+ }
+ }
+ }
+
+ Matrix MelodicICA::sign(const Matrix &Inp)
+ {
+ Matrix Res = Inp;
+ Res = 1;
+ for(int ctr_i = 1; ctr_i <= Inp.Ncols(); ctr_i++){
+ for(int ctr_j = 1; ctr_j <= Inp.Nrows(); ctr_j++){
+ if(Inp(ctr_j,ctr_i)<0){Res(ctr_j,ctr_i)=-1;}
+ }
+ }
+ return Res;
+ }
+
+ void MelodicICA::sort()
+ {
+ int numComp = melodat.get_mix().Ncols();
+
+ Matrix tmpIC = melodat.get_IC();
+ Matrix tmpA = melodat.get_mix();
+
+ for(int ctr_i = 1; ctr_i <= numComp; ctr_i++){
+ if(tmpIC.Row(ctr_i).Sum()>0){
+ tmpIC.Row(ctr_i) = -tmpIC.Row(ctr_i);
+ tmpA.Column(ctr_i) = -tmpA.Column(ctr_i);
+ };
+
+ }
+ melodat.set_mix(tmpA);
+ melodat.set_IC(tmpIC);
+ Matrix tmpW = pinv(tmpA);
+ melodat.set_unmix(tmpW);
+ }
+
+ void MelodicICA::perf_ica(const Matrix &Data)
+ {
+ message("Starting ICA estimation using " << opts.approach.value()
+ << endl << endl);
+ dim = Data.Nrows();
+ samples = Data.Ncols();
+
+ no_convergence = true;
+ //switch to the chosen method
+
+ // cout << endl << "Dim = " << dim << endl << "Samples = " << samples << endl;
+ if(opts.approach.value()==string("symm")){
+ ica_fastica_symm(Data);}
+ if(opts.approach.value()==string("defl")){
+ ica_fastica_defl(Data);}
+ if(opts.approach.value()==string("jade")){
+ 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());
+ // Add the mean time course again
+ temp += melodat.get_unmix()*melodat.get_meanC()*ones(1,temp.Ncols());
+
+ melodat.set_IC(temp);
+ sort();
+ }
+ }
+}
+
+
diff --git a/melica.h b/melica.h
new file mode 100644
index 0000000000000000000000000000000000000000..8cbce7669968416c2e9afb195fe05e510264feba
--- /dev/null
+++ b/melica.h
@@ -0,0 +1,65 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melica.h - ICA estimation
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#ifndef __MELODICICA_h
+#define __MELODICICA_h
+#include "newimageall.h"
+#include "log.h"
+#include "melpca.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melodic.h"
+#include "newmatap.h"
+#include "newmatio.h"
+#include "melreport.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+
+ class MelodicICA
+ {
+ public:
+ MelodicICA(MelodicData &pmelodat, MelodicOptions &popts, Log &plogger, MelodicReport &preport):
+ melodat(pmelodat),
+ opts(popts),
+ logger(plogger),
+ report(preport)
+ {
+ }
+
+ void perf_ica(const Matrix &Data);
+
+ bool no_convergence;
+
+ private:
+ MelodicData &melodat;
+ MelodicOptions &opts;
+ Log &logger;
+ MelodicReport &report;
+
+ int dim;
+ int samples;
+ Matrix redUMM;
+
+ void ica_fastica_symm(const Matrix &Data);
+ void ica_fastica_defl(const Matrix &Data);
+ void ica_maxent(const Matrix &Data);
+ void ica_jade(const Matrix &Data);
+ Matrix randm(const int dim1, const int dim2);
+ void sort();
+ Matrix sign(const Matrix &Inp);
+
+ };
+}
+
+#endif
diff --git a/melodic.cc b/melodic.cc
new file mode 100644
index 0000000000000000000000000000000000000000..13e6a1a3d0376b0aa31a80a9badd25565a9c2e62
--- /dev/null
+++ b/melodic.cc
@@ -0,0 +1,426 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melodic.cc - main program file
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <iostream>
+#include "newmatap.h"
+#include "newmatio.h"
+#include "newimageall.h"
+#include "miscmaths.h"
+#include "miscprob.h"
+#include <string>
+#include <math.h>
+#include "options.h"
+#include "log.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melpca.h"
+#include "melica.h"
+#include "melodic.h"
+#include "melreport.h"
+#include "melgmix.h"
+
+using namespace Utilities;
+using namespace NEWMAT;
+using namespace NEWIMAGE;
+using namespace Melodic;
+using namespace MISCPLOT;
+
+string myfloat2str(float f, int width, int prec, bool scientif)
+ {
+ ostrstream os;
+ int redw = int(std::abs(std::log10(std::abs(f))))+1;
+ if(width>0)
+ os.width(width);
+ if(scientif)
+ os.setf(ios::scientific);
+ os.precision(redw+std::abs(prec));
+ os.setf(ios::internal, ios::adjustfield);
+ os << f << '\0';
+ return os.str();
+ }
+
+Matrix mmall(Log& logger, MelodicOptions& opts,
+ MelodicData& melodat, MelodicReport& report, Matrix& probs);
+void repall(Log& logger, MelodicOptions& opts, MelodicData& melodat,
+ MelodicReport& report, Matrix mmres, Matrix probs);
+void mmonly(Log& logger, MelodicOptions& opts,
+ MelodicData& melodat, MelodicReport& report);
+
+int main(int argc, char *argv[])
+{
+
+ try{
+ // Setup logging:
+ Log& logger = LogSingleton::getInstance();
+
+ // parse command line - will output arguments to logfile
+ MelodicOptions& opts = MelodicOptions::getInstance();
+ opts.parse_command_line(argc, argv, logger, Melodic::version);
+
+ //set up data object
+ MelodicData melodat(opts,logger);
+
+ 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();
+ }
+ else
+ mmonly(logger,opts,melodat,report);
+ }
+ else
+ { // standard PICA now
+ int retry = 0;
+ bool no_conv;
+ bool leaveloop = false;
+
+ melodat.setup();
+
+ do{
+ //do PCA pre-processing
+ MelodicPCA pcaobj(melodat,opts,logger,report);
+ pcaobj.perf_pca(melodat.get_DataVN());
+ pcaobj.perf_white(melodat.get_Data());
+
+ //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){
+ message(endl << "Restarting MELODIC using -d "
+ << opts.pca_dim.value()
+ << endl << endl);
+ }
+ }
+ }
+ } while (no_conv && retry<opts.maxRestart.value() && !leaveloop);
+
+ if(!no_conv){
+ //first 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();
+ }
+ //cerr << mmres.Nrows() << " x " << mmres.Ncols() << endl;
+
+ // if(opts.genreport.value())
+ // repall(logger,opts,melodat,report,mmres,pmaps,threshmaps);
+
+ message("finished!" << endl << endl);
+ } else {
+ message(endl <<"No convergence -- giving up " << endl <<
+ "please contact fsl@fmrib.ox.ac.uk " << endl);
+ }
+ }
+ }
+ catch(Exception e)
+ {
+ cerr << endl << e.what() << endl;
+ }
+ catch(X_OptionError& e)
+ {
+ cerr << endl << e.what() << endl;
+ }
+
+ return 0;
+}
+
+void mmonly(Log& logger, MelodicOptions& opts,
+ MelodicData& melodat, MelodicReport& report){
+
+ Matrix ICs;
+ Matrix mixMatrix;
+ Matrix fmixMatrix;
+ volumeinfo ICvolInfo;
+ volume<float> Mask;
+ volume<float> Mean;
+
+ {
+ volume4D<float> RawData;
+ message("Reading data file " << opts.inputfname.value() << " ... ");
+ read_volume4D(RawData,opts.inputfname.value(),ICvolInfo);
+ message(" done" << endl);
+ Mean = meanvol(RawData);
+ }
+
+ {
+ volume4D<float> RawIC;
+ message("Reading components " << opts.ICsfname.value() << " ... ");
+ read_volume4D(RawIC,opts.ICsfname.value());
+ message(" done" << endl);
+
+ message("Creating mask ... ");
+ Mask = binarise(RawIC[0],float(RawIC[0].min()),float(RawIC[0].max()));
+
+ ICs = RawIC.matrix(Mask);
+ if(ICs.Nrows()>1){
+ Matrix DStDev=stdev(ICs);
+
+ volume4D<float> tmpMask;
+ tmpMask.setmatrix(DStDev,Mask);
+
+ float tMmax;
+ volume<float> tmpMask2;
+ tmpMask2 = tmpMask[0];
+ tMmax = tmpMask2.max();
+ double st_mean = DStDev.Sum()/DStDev.Ncols();
+ double st_std = stdev(DStDev.t()).AsScalar();
+
+ Mask = binarise(tmpMask2,(float) max((float) st_mean-3*st_std,
+ (float) 0.01*st_mean),tMmax);
+ ICs = RawIC.matrix(Mask);
+ }
+ else{
+ Mask = binarise(RawIC[0],float(0.001),float(RawIC[0].max()))
+ + binarise(RawIC[0],float(RawIC[0].min()),float(-0.001));
+ ICs = RawIC.matrix(Mask);
+ }
+
+ //cerr << "ICs : " << ICs.Ncols() << ICs.Nrows() << endl;
+ message(" done" << endl);
+ }
+
+ message("Reading mixing matrix " << opts.filtermix.value());
+ mixMatrix = read_ascii_matrix(opts.filtermix.value());
+ if (mixMatrix.Storage()<=0) {
+ cerr <<" Please specify the mixing matrix correctly" << endl;
+ exit(2);
+ }
+ message(" done" << endl);
+
+ melodat.tempInfo = ICvolInfo;
+ melodat.set_mask(Mask);
+ melodat.set_mean(Mean);
+ melodat.set_IC(ICs);
+ melodat.set_mix(mixMatrix);
+ fmixMatrix = melodat.calc_FFT(mixMatrix);
+ melodat.set_fmix(fmixMatrix);
+ fmixMatrix = pinv(mixMatrix);
+ melodat.set_unmix(fmixMatrix);
+
+ // write_ascii_matrix("ICs",ICs);
+
+ Matrix mmres;
+ Matrix pmaps;//(ICs);
+ if(opts.perf_mm.value())
+ mmres = mmall(logger,opts,melodat,report,pmaps);
+}
+
+void repall(Log& logger, MelodicOptions& opts, MelodicData& melodat,
+ MelodicReport& report, Matrix& mmpars, Matrix& probs)
+{
+ if( bool(opts.genreport.value()) ){
+
+ message(endl
+ << "Creating report in " << report.getDir()
+ << endl);
+
+ for(int ctr=1; ctr<=probs.Nrows(); ctr++){
+ message(" " << ctr);
+ MelGMix mixmod(opts, logger);
+
+ //load MelGMix
+
+ //report.IC_rep(mixmod,ctr,melodat.get_IC().Nrows());
+
+ }
+ }
+}
+
+Matrix mmall(Log& logger, MelodicOptions& opts,
+ MelodicData& melodat, MelodicReport& report, Matrix& pmaps)
+{
+
+ Matrix mmpars(5*melodat.get_IC().Nrows(),5);
+ mmpars = 0;
+ //Matrix pmaps(melodat.get_IC());
+
+ Log stats;
+
+ if(opts.output_MMstats.value()){
+ stats.makeDir(logger.appendDir("stats"),"stats.log");
+ }
+
+ message(endl
+ << "Running Mixture Modelling on Z-transformed IC maps ..."
+ << endl);
+
+ for(int ctr=1; ctr <= melodat.get_IC().Nrows(); ctr++){
+ MelGMix mixmod(opts, logger);
+
+ message(" IC map " << ctr << " ... "<< endl;);
+
+ Matrix ICmap;
+
+ if(melodat.get_stdNoisei().Storage()>0)
+ ICmap = SP(melodat.get_IC().Row(ctr),melodat.get_stdNoisei());
+ else
+ ICmap = melodat.get_IC().Row(ctr);
+
+ string wherelog;
+ if(opts.genreport.value())
+ wherelog = report.getDir();
+ else
+ wherelog = logger.getDir();
+
+ mixmod.setup( ICmap, melodat.tempInfo,
+ wherelog,ctr,
+ melodat.get_mask(),
+ melodat.get_mean(),3);
+ mixmod.fit("GGM");
+
+ if(opts.output_MMstats.value()){
+ melodat.save4D(mixmod.get_probmap(),
+ string("stats/probmap_")+num2str(ctr));
+ }
+ // save probability map
+ //if((mixmod.get_probmap().Storage()>0)&&
+ // (mixmod.get_probmap().Ncols() == pmaps.Ncols()))
+ // pmaps.Row(ctr) = mixmod.get_probmap();
+ //else
+ // pmaps.Row(ctr) = zeros(1,pmaps.Ncols());
+
+ message(" thresholding ... "<< endl);
+ mixmod.threshold(opts.mmthresh.value());
+
+ //re-orient the data
+
+ //message(" done " << endl);
+ Matrix tmp;
+ tmp=(mixmod.get_threshmaps().Row(1));
+ float posint = SP(tmp,gt(tmp,zeros(1,tmp.Ncols()))).Sum();
+ float negint = -SP(tmp,lt(tmp,zeros(1,tmp.Ncols()))).Sum();
+
+ //cerr << posint << " " << negint << endl;
+
+ if((posint<0.01)&&(negint<0.01)){
+ mixmod.clear_infstr();
+ //cerr << "after infstr"<<endl;
+ mixmod.threshold("0.05n "+opts.mmthresh.value());
+ //cerr << " back again" << endl;
+ posint = SP(tmp,gt(tmp,zeros(1,tmp.Ncols()))).Sum();
+ negint = -SP(tmp,lt(tmp,zeros(1,tmp.Ncols()))).Sum();
+ }
+ //cerr << posint << " " << negint << endl;
+ if(negint>posint){//flip map
+ melodat.flipres(ctr);
+ mixmod.flipres(ctr);
+ }
+
+ //save mixture model stats
+ if(opts.output_MMstats.value()){
+ stats << " IC " << num2str(ctr) << " " << mixmod.get_type() << endl
+ << " Means : " << mixmod.get_means() << endl
+ << " Vars. : " << mixmod.get_vars() << endl
+ << " Prop. : " << mixmod.get_pi() << endl << endl;
+ // << " Offs. : " << mixmod.get_offset() << endl << endl;
+ //cerr << mixmod.get_threshmaps().Nrows() << " " << mixmod.get_threshmaps().Ncols()<< endl;
+ melodat.save4D(mixmod.get_threshmaps(),
+ string("stats/thresh_zstat")+num2str(ctr));
+ }
+
+ //save mmpars
+ // mmpars((ctr-1)*5+1,1) = ctr;
+// if(mixmod.get_type()=="GGM")
+// mmpars((ctr-1)*5+1,2) = 1.0;
+// else
+// mmpars((ctr-1)*5+1,2) = 0.0;
+// mmpars((ctr-1)*5+1,2) = mixmod.get_means().Ncols();
+// tmp = mixmod.get_means();
+// for(int ctr2=1;ctr2<=mixmod.get_means().Ncols();ctr2++)
+// mmpars((ctr-1)*5+2,ctr2) = tmp(1,ctr2);
+// tmp = mixmod.get_vars();
+// for(int ctr2=1;ctr2<=mixmod.get_vars().Ncols();ctr2++)
+// mmpars((ctr-1)*5+3,ctr2) = tmp(1,ctr2);
+// tmp = mixmod.get_pi();
+// for(int ctr2=1;ctr2<=mixmod.get_pi().Ncols();ctr2++)
+// mmpars((ctr-1)*5+4,ctr2) = tmp(1,ctr2);
+// mmpars((ctr-1)*5+5,1) = mixmod.get_offset();
+
+
+
+ if( bool(opts.genreport.value()) ){
+ message(" creating report page ... ");
+ report.IC_rep(mixmod,ctr,melodat.get_IC().Nrows());
+ message(" done" << endl);
+ }
+ }
+ if( bool(opts.genreport.value()) ){
+ message(endl << endl <<
+ " To view the output report point your web browser at " <<
+ report.getDir() + "/00index.html" << endl << endl);
+ }
+ if(!opts.filtermode&&opts.filtermix.value().length()==0){
+ //now safe new data
+ bool what = opts.verbose.value();
+ opts.verbose.set_T(false);
+ melodat.save();
+ opts.verbose.set_T(what);
+ }
+ return mmpars;
+}
diff --git a/melodic.h b/melodic.h
new file mode 100644
index 0000000000000000000000000000000000000000..c18372642341ae7e5d0142c629b2c1d5672ceb9d
--- /dev/null
+++ b/melodic.h
@@ -0,0 +1,34 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melodic.h - main program header
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#ifndef __MELODIC_h
+#define __MELODIC_h
+
+#include<strstream>
+
+// a simple message macro that takes care of cout and log
+#define message(msg) { \
+ MelodicOptions& opt = MelodicOptions::getInstance(); \
+ if(opt.verbose.value()) \
+ { \
+ cout << msg; \
+ } \
+ Log& logger = LogSingleton::getInstance(); \
+ logger.str() << msg; \
+}
+
+namespace Melodic{
+
+const string version = "2.0";
+
+}
+
+#endif
diff --git a/meloptions.cc b/meloptions.cc
new file mode 100644
index 0000000000000000000000000000000000000000..96e29c5ec07a37496225f0cdce54cb08446ef713
--- /dev/null
+++ b/meloptions.cc
@@ -0,0 +1,247 @@
+
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ meloptions.cc - class for command line options
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include <iostream>
+#include <fstream>
+#include <stdlib.h>
+#include <stdio.h>
+#include "log.h"
+#include "meloptions.h"
+#include "newimageall.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic {
+
+MelodicOptions* MelodicOptions::gopt = NULL;
+
+ void MelodicOptions::parse_command_line(int argc, char** argv, Log& logger, const string &p_version)
+{
+ //set version number and some other stuff
+ version = p_version;
+ filtermode = false;
+ explicitnums = false;
+ logfname = string("melodic.log");
+
+ // work out the path to the $FSLDIR/bin directory
+ if(getenv("FSLDIR")!=0){
+ binpath = (string) getenv("FSLDIR") + "/bin/";
+ }
+ else{
+ binpath = argv[0];
+ binpath = binpath.substr(0,binpath.length()-7);
+ }
+
+ // parse once to establish log directory name
+ for(int a = options.parse_command_line(argc, argv); a < argc; a++);
+
+ // act on simple command line arguments
+ if(help.value())
+ {
+ print_usage(argc, argv);
+ exit(0);
+ }
+ if(vers.value())
+ {
+ print_version();
+ cout << endl;
+ exit(0);
+ }
+ if(copyright.value())
+ {
+ print_copyright();
+ exit(0);
+ }
+ 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;
+ exit(2);
+ }
+
+ // check for invalid values
+ if (inputfname.value().size()<1) {
+ cerr << "ERROR:: Input volume not found\n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (approach.value() != "symm" && approach.value() != "defl" &&
+ approach.value() != "jade" && approach.value() != "maxent"){
+ cerr << "ERROR:: unknown approach \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (nonlinearity.value() != "pow3" && nonlinearity.value() != "pow4" && nonlinearity.value() != "tanh"
+ && nonlinearity.value() != "gauss" ){
+ cerr << "ERROR:: unknown nonlinearity \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (maxNumItt.value() < 1){
+ cerr << "ERROR:: maxItt too small \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (epsilon.value() < 0.000000001){
+ cerr << "ERROR:: epsilon too small \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (epsilon.value() >= 0.01){
+ cerr << "ERROR:: epsilon too large \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (nlconst1.value() <= 0){
+ cerr << "ERROR:: nlconst1 negative \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (!remove_meanvol.value()){
+ varnorm.set_T(false);
+ }
+ if (filter.value().length()>0){
+ if (filtermix.value().length()<0){
+ cerr << "ERROR:: no mixing matrix for filtering (use --mix='filename') \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ } else {
+ filtermode = true;
+ varnorm.set_T(false);
+ }
+ }
+ if (threshold.value()<=0){
+ use_mask.set_T(false);
+ }
+ if (output_pca.value()){
+ output_white.set_T(true);
+ }
+ if(threshold.value()>=1){
+ threshold.set_T(threshold.value()/100);
+ if(threshold.value()>=1){
+ cerr << "ERROR:: threshold level not a percentage value \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ }
+ if (nlconst2.value() <= 0){
+ cerr << "ERROR:: nlconst2 negative \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (dummy.value() < 0){
+ cerr << "ERROR:: negative dummy value \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (repeats.value() < 1){
+ cerr << "ERROR:: repeats < 1 \n\n";
+ print_usage(argc,argv);
+ exit(2);
+ }
+ if (numICs.value() > 0){
+ explicitnums = true;
+ }
+
+ //transform inputfname to its basename
+ string basename = inputfname.value();
+ make_basename(basename);
+ inputfname.set_T(basename);
+
+ //create melodic directory name
+ if(logdir.value()==string("melodic.log")){
+ logdir.set_T(string(inputfname.value()+".ica"));
+ logger.makeDir(logdir.value(),logfname);
+ }
+ else{
+ // setup logger directory
+ system(("mkdir "+ logdir.value() + " 2>/dev/null").c_str());
+ logger.setDir(logdir.value(),logfname);
+ }
+ message(endl << "Melodic Version " << version << endl << endl);
+
+ // parse again so that options are logged
+ for(int a = 0; a < argc; a++)
+ logger.str() << argv[a] << " ";
+ logger.str() << endl << "---------------------------------------------" << endl << endl;
+
+ message("Melodic results will be in " << logger.getDir() << endl << endl);
+}
+
+void MelodicOptions::print_usage(int argc, char *argv[])
+{
+ print_version();
+ options.usage();
+ /* cout << "Usage: " << argv[0] << " ";
+
+ Have own usage output here
+ */
+}
+
+void MelodicOptions::print_version()
+{
+ cout << "MELODIC Version " << version;
+}
+
+void MelodicOptions::print_copyright()
+{
+ print_version();
+ cout << endl << "Copyright (C) 1999-2002 University of Oxford " << endl;
+}
+
+void MelodicOptions::status()
+{
+ cout << " version = " << version << endl;
+
+ cout << " logdir = " << logdir.value() << endl;
+ cout << " inputfname = " << inputfname.value() << endl;
+ cout << " outputfname = " << outputfname.value() << endl;
+ cout << " guessfname = " << guessfname.value() << endl;
+ cout << " maskfname = " << maskfname.value() << endl;
+ cout << " paradigmfname = " << paradigmfname.value() << endl;
+ cout << " nonlinearity = " << nonlinearity.value() << endl;
+ cout << " approach = " << approach.value() << endl;
+
+ cout << " pca_dim = " << pca_dim.value() << endl;
+ cout << " segment = " << segment.value() << endl;
+ cout << " numICs = " << numICs.value() << endl;
+ cout << " maxNumItt = " << maxNumItt.value() << endl;
+ cout << " maxRestart = " << maxRestart.value() << endl;
+ cout << " dummy = " << dummy.value() << endl;
+ cout << " repeats = " << repeats.value() << endl;
+ cout << " nlconst1 = " << nlconst1.value() << endl;
+ cout << " nlconst2 = " << nlconst2.value() << endl;
+ cout << " epsilon = " << epsilon.value() << endl;
+ cout << " threshold = " << threshold.value() << endl;
+
+ cout << " help = " << help.value() << endl;
+ cout << " verbose = " << verbose.value() << endl;
+ cout << " vers = " << vers.value() << endl;
+ cout << " copyright = " << copyright.value() << endl;
+ cout << " perf_bet = " << perf_bet.value() << endl;
+ cout << " remove_meanvol = " << remove_meanvol.value() << endl;
+ cout << " remove_endslices = " << remove_endslices.value() << endl;
+ cout << " output_pca = " << output_pca.value() << endl;
+ cout << " output_white = " << output_white.value() << endl;
+ cout << " output_mean = " << output_mean.value() << endl;
+ cout << " use_mask = " << use_mask.value() << endl;
+ cout << " output_unmix = " << output_unmix.value() << endl;
+ cout << " guess_remderiv = " << guess_remderiv.value() << endl;
+ cout << " filter = " << filter.value() << endl;
+ cout << " logPower = " << logPower.value() << endl;
+}
+
+}
diff --git a/meloptions.h b/meloptions.h
new file mode 100644
index 0000000000000000000000000000000000000000..8242c0db4c58ffd44a83bf0001d76dadbcd79643
--- /dev/null
+++ b/meloptions.h
@@ -0,0 +1,344 @@
+
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ meloptions.h - class for command line options
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+
+#ifndef __MELODICOPTIONS_h
+#define __MELODICOPTIONS_h
+
+#include <string>
+#include <strstream>
+#include <iostream>
+#include <fstream>
+#include <stdlib.h>
+#include <stdio.h>
+#include "options.h"
+#include "log.h"
+#include "melodic.h"
+
+
+using namespace Utilities;
+
+namespace Melodic {
+
+class MelodicOptions {
+ public:
+ static MelodicOptions& getInstance();
+ ~MelodicOptions() { delete gopt; }
+
+ string version;
+ string binpath;
+ string logfname;
+ bool filtermode;
+ bool explicitnums;
+
+ Option<string> logdir;
+ Option<string> inputfname;
+
+ Option<string> outputfname;
+
+ Option<string> maskfname;
+ Option<bool> use_mask;
+ Option<bool> perf_bet;
+ Option<float> threshold;
+
+ Option<int> pca_dim;
+ Option<string> pca_est;
+ Option<int> numICs;
+ Option<string> approach;
+ Option<string> nonlinearity;
+
+ Option<bool> varnorm;
+ Option<string> segment;
+ Option<bool> tsmooth;
+ Option<float> epsilon;
+ Option<int> maxNumItt;
+ Option<int> maxRestart;
+
+ Option<string> mmthresh;
+ Option<bool> perf_mm;
+ Option<string> ICsfname;
+ Option<string> filtermix;
+ Option<string> filter;
+
+ Option<bool> genreport;
+ Option<float> tr;
+ Option<bool> logPower;
+
+ Option<bool> output_unmix;
+ Option<bool> output_MMstats;
+ Option<bool> output_pca;
+ Option<bool> output_white;
+ Option<bool> output_origIC;
+ Option<bool> output_mean;
+
+ Option<bool> verbose;
+ Option<bool> vers;
+ Option<bool> copyright;
+ Option<bool> help;
+
+ Option<string> guessfname;
+ Option<string> paradigmfname;
+
+ Option<int> dummy;
+ Option<int> repeats;
+ Option<float> nlconst1;
+ Option<float> nlconst2;
+
+
+ Option<bool> remove_meanvol;
+ Option<bool> remove_endslices;
+
+ Option<bool> guess_remderiv;
+ Option<bool> temporal;
+
+ int retrystep;
+
+ void parse_command_line(int argc, char** argv, Log& logger, const string &p_version);
+
+ private:
+ MelodicOptions();
+ const MelodicOptions& operator=(MelodicOptions&);
+ MelodicOptions(MelodicOptions&);
+
+ OptionParser options;
+
+ static MelodicOptions* gopt;
+
+ void print_usage(int argc, char *argv[]);
+ void print_version();
+ void print_copyright();
+ void status();
+};
+
+ inline MelodicOptions& MelodicOptions::getInstance(){
+ if(gopt == NULL)
+ gopt = new MelodicOptions();
+
+ return *gopt;
+ }
+
+ inline MelodicOptions::MelodicOptions() :
+ logdir(string("-o,--outdir"), string("melodic.log"),
+ string("output directory name\n"),
+ false, requires_argument),
+ inputfname(string("-i,--in"), string(""),
+ string("input file name"),
+ true, requires_argument),
+ outputfname(string("-O,--out"), string("melodic"),
+ string("output file name"),
+ false, requires_argument,false),
+ maskfname(string("-m,--mask"), string(""),
+ string("file name of mask for thresholding"),
+ false, requires_argument),
+ use_mask(string("--nomask"), true,
+ string("switch off masking"),
+ false, no_argument),
+ perf_bet(string("--nobet"), true,
+ string("\tswitch off BET"),
+ false, no_argument),
+ threshold(string("--bgthreshold"), 0.00000001,
+ string("brain / non-brain threshold (only if --nobet selected)\n"),
+ false, requires_argument),
+ pca_dim(string("-d,--dim"), 0,
+ string("dimensionality reduction into #num dimensions (default: automatic estimation)"),
+ false, requires_argument),
+ pca_est(string("--dimest"), string("lap"),
+ string("use specific dim. estimation technique: lap, bic, mdl, aic, mean (default: lap)"),
+ false, requires_argument),
+ numICs(string("-n,--numICs"), -1,
+ string("numer of IC's to extract (for deflation approach)"),
+ false, requires_argument),
+ approach(string("-a,--approach"), string("symm"),
+ string("approach for ICA estimation: defl, symm"),
+ false, requires_argument),
+ nonlinearity(string("--nl"), string("pow3"),
+ string("\tnonlinearity: gauss, tanh, pow3, pow4"),
+ false, requires_argument),
+ varnorm(string("--vn,--varnorm"), true,
+ string("switch off variance normalisation"),
+ false, no_argument),
+ segment(string("--covarweight"), string(""),
+ string("voxel-wise weights for the covariance matrix (e.g. segmentation information)"),
+ false, requires_argument),
+ tsmooth(string("--spca"), false,
+ string("smooth the eigenvectors prior to IC decomposition"),
+ false, no_argument, false),
+ epsilon(string("--eps,--epsilon"), 0.00005,
+ string("minimum error change"),
+ false, requires_argument),
+ maxNumItt(string("--maxit"), 500,
+ string("\tmaximum number of iterations before restart"),
+ false, requires_argument),
+ maxRestart(string("--maxrestart"), 6,
+ string("maximum number of restarts\n"),
+ false, requires_argument),
+ mmthresh(string("--mmthresh"), string("0.5"),
+ string("threshold for Mixture Model based inference"),
+ false, requires_argument),
+ perf_mm(string("--no_mm"), true,
+ string("\tswitch off mixture modelling on IC maps\n "),
+ false, no_argument),
+ ICsfname(string("--ICs"), string(""),
+ string("\tfilename of the IC components file for mixture modelling"),
+ false, requires_argument),
+ filtermix(string("--mix"), string(""),
+ string("\tmixing matrix for mixture modelling / filtering"),
+ 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),
+ tr(string("--tr"), 0.0,
+ string("\tTR in seconds"),
+ false, requires_argument),
+ logPower(string("--logPower"), false,
+ string("calculate log of power for frequency spectrum\n"),
+ false, no_argument),
+ output_unmix(string("--Ounmix"), false,
+ string("output unmixing matrix"),
+ false, no_argument),
+ output_MMstats(string("--Ostats"), false,
+ string("output thresholded maps and probability maps"),
+ false, no_argument),
+ output_pca(string("--Opca"), false,
+ string("\toutput PCA results"),
+ false, no_argument),
+ output_white(string("--Owhite"), false,
+ string("output whitening/dewhitening matrices"),
+ false, no_argument),
+ output_origIC(string("--Oorig"), false,
+ string("\toutput the original ICs"),
+ false, no_argument),
+ output_mean(string("--Omean"), false,
+ string("\toutput mean volume\n"),
+ false, no_argument),
+ verbose(string("-v,--verbose"), false,
+ string("switch on diagnostic messages"),
+ false, no_argument),
+ vers(string("-V,--version"), false,
+ string("prints version information"),
+ false, no_argument),
+ copyright(string("--copyright"), false,
+ string("prints copyright information"),
+ false, no_argument),
+ help(string("-h,--help"), false,
+ string("prints this help message"),
+ false, no_argument),
+ guessfname(string("-g,--guess"), string(""),
+ string("file name of guess of mixing matrix"),
+ false, requires_argument, false),
+ paradigmfname(string("-p,--paradigm"), string(""),
+ string("file name of FEAT paradigm file"),
+ false, requires_argument, false),
+ dummy(string("--dummy"), 0,
+ string("number of dummy volumes"),
+ false, requires_argument,false),
+ repeats(string("--repeats"), 1,
+ string("number of repeats (multistart)"),
+ false, requires_argument, false),
+ nlconst1(string("--nl1,--nlconst1"), 1.0,
+ string("nonlinear constant 1"),
+ false, requires_argument, false),
+ nlconst2(string("--nl2,--nlconst2"), 1.0,
+ string("nonlinear constant 2"),
+ false, requires_argument, false),
+ remove_meanvol(string("--keepmeanvol"), true,
+ string("do not subtract mean volume"),
+ false, no_argument, false),
+ remove_endslices(string("--remEndslices"), false,
+ string("delete end slices (motion correction artefacts)"),
+ false, no_argument,false),
+ guess_remderiv(string("--remderiv"), false,
+ string("removes every second entry in paradigm file (EV derivatives)"),
+ false, no_argument, false),
+ temporal(string("--temporal"), false,
+ string("perform temporal ICA"),
+ false, no_argument, false),
+ retrystep(3),
+ 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> --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(inputfname);
+ options.add(outputfname);
+ options.add(guessfname);
+ options.add(maskfname);
+ options.add(use_mask);
+ options.add(perf_bet);
+ options.add(threshold);
+ options.add(pca_dim);
+ options.add(pca_est);
+ options.add(numICs);
+ options.add(approach);
+ options.add(nonlinearity);
+ options.add(varnorm);
+ options.add(segment);
+ options.add(tsmooth);
+ options.add(epsilon);
+ options.add(maxNumItt);
+ options.add(maxRestart);
+ options.add(mmthresh);
+ options.add(perf_mm);
+ options.add(ICsfname);
+ options.add(filtermix);
+ options.add(filter);
+ options.add(genreport);
+ options.add(tr);
+ options.add(logPower);
+ options.add(output_unmix);
+ options.add(output_MMstats);
+ options.add(output_pca);
+ options.add(output_white);
+ options.add(output_origIC);
+ options.add(output_mean);
+ options.add(verbose);
+ options.add(vers);
+ options.add(copyright);
+ options.add(help);
+
+ options.add(guessfname);
+ options.add(paradigmfname);
+ options.add(dummy);
+ options.add(repeats);
+ options.add(nlconst1);
+ options.add(nlconst2);
+ options.add(remove_meanvol);
+ options.add(remove_endslices);
+ options.add(guess_remderiv);
+ options.add(temporal);
+ }
+ catch(X_OptionError& e) {
+ options.usage();
+ cerr << endl << e.what() << endl;
+ }
+ catch(std::exception &e) {
+ cerr << e.what() << endl;
+ }
+
+ }
+}
+
+#endif
+
+
+
diff --git a/melpca.cc b/melpca.cc
new file mode 100644
index 0000000000000000000000000000000000000000..66a5f4277be06707cf9c8bd4bb1b90f284d2f5ef
--- /dev/null
+++ b/melpca.cc
@@ -0,0 +1,586 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melpca.cc - PCA and whitening
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-200 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "newimageall.h"
+#include "log.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melodic.h"
+#include "newmatap.h"
+#include "newmatio.h"
+#include "melpca.h"
+#include "libprob.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+
+ void MelodicPCA::perf_pca(const Matrix &Data)
+ {
+ message("Starting PCA ... ");
+
+ SymmetricMatrix Corr;
+
+ if(opts.segment.value().length()>0){
+ Matrix Res;
+ Res = ones(Data.Nrows(),1)*melodat.get_RXweight();
+ Res = SP(melodat.get_Data(),Res);
+ Corr = cov(Res.t());
+ }
+ else{
+ Corr = cov(melodat.get_Data().t());
+ }
+
+ Matrix tmpE;
+ DiagonalMatrix tmpD;
+
+ EigenValues(Corr,tmpD,tmpE);
+
+ if(opts.tsmooth.value()){
+ message(" temporal smoothing of Eigenvectors " << endl);
+ tmpE=melodat.smoothColumns(tmpE);
+ }
+
+ melodat.set_pcaE(tmpE);
+ melodat.set_pcaD(tmpD);
+
+ RowVector AdjEV;
+ AdjEV = tmpD.AsRow().Reverse();
+ SortDescending(AdjEV);
+
+ RowVector PercEV(AdjEV);
+ PercEV = cumsum(AdjEV / sum(AdjEV,2).AsScalar());
+ write_ascii_matrix(logger.appendDir("eigenvalues_percent"),PercEV);
+ melodat.set_EVP(PercEV);
+ AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar());
+ melodat.set_EV((AdjEV));
+ message("done" << endl);
+ }
+
+ void MelodicPCA::perf_white(const Matrix &Data)
+ {
+ Matrix RE;
+ DiagonalMatrix RD;
+ Matrix tmpWhite;
+ Matrix tmpDeWhite;
+
+ int N = melodat.get_pcaE().Ncols();
+ if(opts.pca_dim.value() > N){
+ message("dimensionality too large - using -dim " << N <<
+ " instead " << endl);
+ opts.pca_dim.set_T(N);
+ }
+ if(opts.pca_dim.value() < 0){
+ if(opts.remove_meanvol.value()){
+ opts.pca_dim.set_T(N-2);
+ }else{
+ opts.pca_dim.set_T(N-1);
+ }
+ }
+ if(opts.pca_dim.value() ==0){
+ opts.pca_dim.set_T(pcadim());
+ if(melodat.get_Data().Nrows()<20){
+ opts.pca_dim.set_T(N-2);
+ message("too few data points for full estimation, using "
+ << opts.pca_dim.value() << " instead"<< endl);
+ }
+ }
+ if(opts.approach.value()==string("jade") && opts.pca_dim.value() > 30){
+ message("dimensionality too large for jade estimation - using --dim 30 instead" << endl);
+ opts.pca_dim.set_T(30);
+ }
+ message("Start whitening using "<< opts.pca_dim.value()<<" dimensions ... " << endl);
+ RowVector tmpEVP;
+ tmpEVP << melodat.get_EVP();
+ float varp = 1.0;
+ if(opts.pca_dim.value() <= tmpEVP.Ncols()){
+ varp = tmpEVP(opts.pca_dim.value());
+ }
+ message(" retaining "<< varp*100 <<" percent of the variability " << endl);
+
+ RE = melodat.get_pcaE().Columns(N-opts.pca_dim.value()+1,N);
+ RD << abs(melodat.get_pcaD().SymSubMatrix(N-opts.pca_dim.value()+1,N));
+
+ tmpWhite = sqrt(abs(RD.i()))*RE.t();
+ tmpDeWhite = RE*sqrt(RD);
+ melodat.set_white(tmpWhite);
+ melodat.set_dewhite(tmpDeWhite);
+ message(" ... done"<< endl << endl);
+ }
+
+ int MelodicPCA::pcadim()
+ {
+ message("Estimating the number of sources from the data (PPCA) ..." << endl);
+
+ //First, estimate the smoothness of the data;
+ // set up all strings
+
+ string SM_path = opts.binpath + "smoothest";
+ string Mask_fname = logger.appendDir("mask");
+
+ if(opts.segment.value().length()>0){
+ Mask_fname = opts.segment.value();
+ }
+
+ // Setup external call to smoothest:
+ char callSMOOTHESTstr[1000];
+ ostrstream osc(callSMOOTHESTstr,1000);
+ osc << SM_path << " -d " << melodat.data_dim()
+ << " -r " << opts.inputfname.value() << " -m "
+ << Mask_fname << " > " << logger.appendDir("smoothest") << '\0';
+
+ message(" Calling Smoothest: " << callSMOOTHESTstr << endl);
+ system(callSMOOTHESTstr);
+
+ //read back the results
+ ifstream in;
+ string str;
+ float Resel = 1.0;
+
+ in.open(logger.appendDir("smoothest").c_str(), ios::in);
+ if(in>0){
+ for(int ctr=1; ctr<7; ctr++){ in >> str;}
+ in.close();
+ if(str!="nan"){
+ Resel = atof(str.c_str());
+ }
+ }
+
+ //cerr << " Resels : " << Resel << endl << endl;
+
+ melodat.set_resels(Resel);
+
+ Matrix PPCAest;
+
+ // if(!opts.varnorm.value()){
+ SymmetricMatrix Corr;
+ if(opts.segment.value().length()>0){
+ Matrix Res;
+ Res = ones(melodat.get_DataVN().Nrows(),1)*melodat.get_RXweight();
+ Res = SP(melodat.get_DataVN(),Res);
+ Corr = cov(Res.t());
+ }
+ else{
+ Corr = cov(melodat.get_DataVN().t());
+ }
+ DiagonalMatrix tmpD;
+ Matrix tmpE;
+ EigenValues(Corr,tmpD,tmpE);
+ // }
+
+ RowVector AdjEV;
+ AdjEV << tmpD.AsRow();
+ AdjEV = AdjEV.Columns(3,AdjEV.Ncols());
+ AdjEV = AdjEV.Reverse();
+
+ RowVector CircleLaw;
+ int NumVox = (int) floor(melodat.data_samples()/(2.5*Resel));
+
+
+ CircleLaw = Feta(int(AdjEV.Ncols()), NumVox);
+
+ for(int ctr=1;ctr<=CircleLaw.Ncols(); ctr++){
+ if(CircleLaw(ctr)<5*10e-10){CircleLaw(ctr) = 5*10e-10;}
+ }
+
+ //write_ascii_matrix(logger.appendDir("tmpA1"),AdjEV);
+ //AdjEV = AdjEV.Columns(2,AdjEV.Ncols());
+ //write_ascii_matrix(logger.appendDir("tmpA2"),AdjEV);
+
+ //cerr<< AdjEV.Nrows() << " x " << AdjEV.Ncols() << endl;
+
+ //cerr<< CircleLaw.Nrows() << " x " << CircleLaw.Ncols() << endl;
+
+ float slope;
+ slope = CircleLaw.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
+ int(AdjEV.Ncols()/4)).Sum() /
+ AdjEV.Columns(int(AdjEV.Ncols()/4),AdjEV.Ncols() -
+ int(AdjEV.Ncols()/4)).Sum();
+
+ //CircleLaw = slope * (CircleLaw-1) + 1;
+
+ // write_ascii_matrix(logger.appendDir("claw"),CircleLaw.Columns(1,AdjEV.Ncols()));
+
+ RowVector PercEV(AdjEV);
+ PercEV = cumsum(AdjEV / sum(AdjEV,2).AsScalar());
+
+ // write_ascii_matrix(logger.appendDir("ev"),AdjEV);
+
+ //cerr << int(AdjEV.Ncols()) << " " << NumVox << " " << slope << endl;
+ AdjEV << SP(AdjEV,pow(CircleLaw.Columns(1,AdjEV.Ncols()),-1));
+ // cerr << "recalculated" << endl;
+
+ SortDescending(AdjEV);
+ int maxEV = 1;
+ float threshold = 0.98;
+ for(int ctr_i = 1; ctr_i < PercEV.Ncols(); ctr_i++){
+ if((PercEV(ctr_i)<threshold)&&(PercEV(ctr_i+1)>=threshold)){maxEV=ctr_i;}
+ }
+
+ if(maxEV<3){maxEV=PercEV.Ncols()/2;}
+ RowVector NewEV;
+ Matrix temp1;
+ temp1 = abs(AdjEV);
+ NewEV << temp1.SubMatrix(1,1,1,maxEV);
+ PPCAest = ppca_est(NewEV, NumVox);
+ RowVector estimators(5);
+ estimators = 1.0;
+
+ Matrix PPCA2(PPCAest);
+
+ for(int ctr=1; ctr<=PPCA2.Ncols(); ctr++){
+ PPCA2.Column(ctr) = (PPCA2.Column(ctr) -
+ min(PPCA2.Column(ctr)).AsScalar()) /
+ ( max(PPCA2.Column(ctr)).AsScalar() -
+ min(PPCA2.Column(ctr)).AsScalar());
+ }
+
+ int ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,2) < PPCAest(ctr_i+1,2))&&(ctr_i<maxEV))
+ {estimators(1)=ctr_i+1;ctr_i++;}
+ ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,3) < PPCAest(ctr_i+1,3))&&(ctr_i<maxEV))
+ {estimators(2)=ctr_i+1;ctr_i++;}
+ ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,4) < PPCAest(ctr_i+1,4))&&(ctr_i<maxEV))
+ {estimators(3)=ctr_i+1;ctr_i++;}
+ ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,5) < PPCAest(ctr_i+1,5))&&(ctr_i<maxEV))
+ {estimators(4)=ctr_i+1;ctr_i++;}
+ ctr_i = 1;
+ while((ctr_i< PPCAest.Nrows()-1)&&
+ (PPCAest(ctr_i,6) < PPCAest(ctr_i+1,6))&&(ctr_i<maxEV))
+ {estimators(5)=ctr_i+1;ctr_i++;}
+
+ int res = 0;
+ ColumnVector PPCA;
+
+ if(opts.pca_est.value() == string("lap")){
+ res = int(estimators(1));
+ PPCA << PPCA2.Column(2);
+ }
+
+ if(opts.pca_est.value() == string("bic")){
+ res = int(estimators(2));
+ PPCA << PPCA2.Column(2);
+ }
+ if(opts.pca_est.value() == string("mdl")){
+ res = int(estimators(3));
+ PPCA << PPCA2.Column(4);
+ }
+ if(opts.pca_est.value() == string("aic")){
+ res = int(estimators(5));
+ PPCA << PPCA2.Column(6);
+ }
+ if(res==0){//median estimator
+ PPCA = PPCA2.Column(2);
+
+ for(int ctr=1; ctr<=PPCA2.Nrows(); ctr++){
+ RowVector tmp = PPCA2.SubMatrix(ctr,ctr,2,6);
+// SortAscending(tmp);
+// float themean = float(tmp.Sum()/5);
+// if(std::abs(int(tmp(2)-themean)) < std::abs(int(tmp(3)-themean)))
+// PPCA(ctr) = tmp(2);
+// else
+// PPCA(ctr) = tmp(3);
+
+ PPCA(ctr) = float(tmp.Sum()/5);
+ }
+
+ ctr_i = 1;
+ while((PPCA(ctr_i) < PPCA(ctr_i+1))&&(ctr_i<maxEV)){
+ res=ctr_i+1;ctr_i++;
+ }
+ }
+
+ // cerr << estimators << " " << res << endl;
+
+ //write_ascii_matrix(logger.appendDir("PPCA2"),PPCA2);
+ AdjEV = (AdjEV - min(AdjEV).AsScalar())/(max(AdjEV).AsScalar() - min(AdjEV).AsScalar());
+
+
+ write_ascii_matrix(logger.appendDir("PPCA"),PPCA);
+ write_ascii_matrix(logger.appendDir("eigenvalues_adjusted"),AdjEV.t());
+ write_ascii_matrix(logger.appendDir("eigenvalues_percent"),PercEV.t());
+
+ melodat.set_EVP(PercEV);
+ melodat.set_EV(AdjEV);
+ melodat.set_PPCA(PPCA);
+
+ //PPCA << sum(PPCAest.Columns(2,6),2);
+
+ //ctr_i = 1;
+
+ //while((PPCA(ctr_i) < PPCA(ctr_i+1))&&(ctr_i<maxEV)){
+ // res=ctr_i+1;ctr_i++;
+ //}
+
+ //res = int(sum(estimators,2).AsScalar()/5);
+
+ // res = int(estimators(1)); // Laplace approximation
+ // SortAscending(estimators);
+ // if(std::abs(int(estimators(2))-res) < std::abs(int(estimators(3))-res))
+ // res = int(estimators(2));
+ // else
+ // res = int(estimators(3));
+
+
+ //write_ascii_matrix(logger.appendDir("PPCA"),PPCAest);
+ //write_ascii_matrix(logger.appendDir("dimest"),estimators);
+
+ return res;
+ }
+
+RowVector MelodicPCA::Feta(int n1, int n2)
+ {
+ float nu = (float) n1/n2;
+ float bm = pow((1-sqrt(nu)),2.0);
+ float bp = pow((1+sqrt(nu)),2.0);
+
+ //cerr << "nu, bm, bp " << nu << " " <<bm << " " << bp << endl;
+
+ float lrange = 0.9*bm;
+ float urange = 1.1*bp;
+
+ // int dummy;
+
+ RowVector eta(30*n1);
+ float rangestepsize = (urange - lrange) / eta.Ncols();
+ for(int ctr_i = 0; ctr_i < eta.Ncols(); ctr_i++){
+ eta(ctr_i+1) = lrange + rangestepsize * (ctr_i);
+ }
+
+ RowVector teta(10*n1);
+ teta = 0;
+ float stepsize = (bp - bm) / teta.Ncols();
+ for(int ctr_i = 0; ctr_i < teta.Ncols(); ctr_i++){
+ teta(ctr_i+1) = stepsize*(ctr_i);
+ }
+
+ //cerr << teta(1)+bm << " " << teta(1000)+bm << endl;
+ //cerr << eta(1)<< " " << eta(eta.Ncols())<< endl;
+
+ //cerr << "BP1" << endl;
+
+ //write_ascii_matrix(logger.appendDir("teta"),teta.t());
+
+
+ // RowVector tmp1(teta);
+ // tmp1 = teta + bm;
+
+ // cerr << "tmp1" << endl;
+ // tmp1 = pow(2*M_PI*nu*(tmp1),-1);
+ // cerr << "tmp1" << endl;
+
+ // RowVector tmp2(teta);
+ // cerr << "tmp2" << endl;
+ // tmp2 = SP(teta, bp-bm-teta);
+ // cerr << "tmp2" << endl;
+ // tmp2=abs(tmp2);
+ // cerr << "tmp2" << endl;
+
+
+ RowVector feta(teta);
+ feta = SP(pow(2*M_PI*nu*(teta + bm),-1), pow(SP(teta, bp-bm-teta),0.5));
+
+ //Matrix location;
+ teta = teta + bm;
+
+ //cerr << "teta : " << teta.Nrows() << " x " << teta.Ncols() << endl;
+ //cerr << "eta : " << eta.Nrows() << " x " << eta.Ncols() << endl;
+ //cerr << "feta : " << feta.Nrows() << " x " << feta.Ncols() << endl;
+ //c/err << "vor location (input)" << endl;
+ //cin >> dummy;
+ //location = SP(teta.t()*ones(1,eta.Ncols()),pow(ones(teta.Ncols(),1)*eta,-1));
+ //cerr << "nach location (input)" << endl;
+ //cin >> dummy;
+ //cerr << " weiter " << endl;
+
+ //for(int ctr_i = 1; ctr_i <= location.Ncols(); ctr_i++){
+ // for(int ctr_j = 1; ctr_j <= location.Nrows(); ctr_j++){
+ // if(location(ctr_j,ctr_i)<1){location(ctr_j,ctr_i)=1;}
+ // else {location(ctr_j,ctr_i)=0;}
+ // }
+ // }
+ //write_ascii_matrix(logger.appendDir("location"),location);
+ // write_ascii_matrix(logger.appendDir("teta"),teta);
+ //write_ascii_matrix(logger.appendDir("eta"),eta);
+ //write_ascii_matrix(logger.appendDir("feta"),feta);
+
+
+ //RowVector claw;
+ // claw = n1*(1-sum(SP(stepsize*feta.t()*ones(1,eta.Ncols()),location),1).AsRow());
+
+ RowVector claw(eta);
+ claw = 0;
+ for(int ctr_i = 1; ctr_i <= eta.Ncols(); ctr_i++){
+ double tmpval = 0.0;
+ for(int ctr_j = 1; ctr_j <= teta.Ncols(); ctr_j++){
+ if(( double(teta(ctr_j))/double(eta(ctr_i)) )<1)
+ tmpval += feta(ctr_j);
+ }
+ claw(ctr_i) = n1*(1-stepsize*tmpval);
+ }
+
+ //write_ascii_matrix(logger.appendDir("claw"),claw);
+ //cerr << "BP1" << endl;
+ RowVector Res(n1); //invert the CDF
+ //cerr << "n1=" << n1 << endl;
+ for(int ctr_i = 1; ctr_i < eta.Ncols(); ctr_i++){
+ if(floor(claw(ctr_i))>floor(claw(ctr_i+1))){
+ // cerr << int(floor(claw(ctr_i))) << " ";
+ Res(int(floor(claw(ctr_i)))) = eta(ctr_i);
+ }
+ }
+
+ //cerr << endl;
+ // cerr << " Done with loop " << int(floor(tmp4b(ctr_i))) << endl;
+ //write_ascii_matrix(logger.appendDir("claw-dstn"),Res);
+ return Res;
+ }
+
+
+ RowVector MelodicPCA::cumsum(const RowVector& Inp)
+ {
+ UpperTriangularMatrix UT(Inp.Ncols());
+ UT=1.0;
+ RowVector Res;
+ Res = Inp * UT;
+ return Res;
+ }
+
+
+ Matrix MelodicPCA::ppca_est(const RowVector& eigenvalues, const int N)
+ {
+ RowVector logLambda(eigenvalues);
+ logLambda = log(eigenvalues);
+
+ int d = logLambda.Ncols();
+
+ RowVector k(d);
+ for(int ctr_i = 1; ctr_i <=d; ctr_i++){
+ k(ctr_i)=ctr_i;
+ }
+
+ RowVector m(d);
+ m=d*k-0.5*SP(k,k+1);
+
+ RowVector loggam(d);
+ loggam = 0.5*k.Reverse();
+ for(int ctr_i = 1; ctr_i <=d; ctr_i++){
+ loggam(ctr_i)=lgam(loggam(ctr_i));
+ }
+ loggam = cumsum(loggam);
+
+ RowVector l_probU(d);
+ l_probU = -log(2)*k + loggam - cumsum(0.5*log(M_PI)*k.Reverse());
+
+ RowVector tmp1;
+ tmp1 = -cumsum(eigenvalues).Reverse()+sum(eigenvalues,2).AsScalar();
+ tmp1(1) = 0.95*tmp1(2);
+ tmp1=tmp1.Reverse();
+
+ RowVector tmp2;
+ tmp2 = -cumsum(logLambda).Reverse()+sum(logLambda,2).AsScalar();
+ tmp2(1)=tmp2(2);
+ tmp2=tmp2.Reverse();
+
+ RowVector tmp3;
+ tmp3 = d-k;
+ tmp3(d) = 1.0;
+
+ RowVector tmp4;
+ tmp4 = SP(tmp1,pow(tmp3,-1));
+ for(int ctr_i = 1; ctr_i <=d; ctr_i++){
+ if(tmp4(ctr_i)<0.01){tmp4(ctr_i)=0.01;}
+ if(tmp3(ctr_i)<0.01){tmp3(ctr_i)=0.01;}
+ if(tmp1(ctr_i)<0.01){tmp1(ctr_i)=0.01;}
+ }
+
+ RowVector l_nu;
+ l_nu = SP(-N/2*(d-k),log(tmp4));
+ l_nu(d) = 0;
+
+ RowVector l_lam;
+ l_lam = -(N/2)*cumsum(logLambda);
+
+ RowVector l_lhood;
+ l_lhood = SP(pow(tmp3,-1),tmp2)-log(SP(pow(tmp3,-1),tmp1));
+
+ Matrix t1,t2, t3;
+ UpperTriangularMatrix triu(d);
+ triu = 1.0;
+ for(int ctr_i = 1; ctr_i <= triu.Ncols(); ctr_i++){
+ triu(ctr_i,ctr_i)=0.0;
+ }
+ t1 = (ones(d,1) * eigenvalues);
+ t1 = SP(triu,t1.t() - t1);
+ t2 = pow(tmp4,-1).t()*ones(1,d);
+ t3 = ones(d,1)*pow(eigenvalues,-1);
+ t2 = SP(triu, t2.t()-t3.t());
+ for(int ctr_i = 1; ctr_i <= t1.Ncols(); ctr_i++){
+ for(int ctr_j = 1; ctr_j <= t1.Nrows(); ctr_j++){
+ if(t1(ctr_j,ctr_i)<=0){t1(ctr_j,ctr_i)=1;}
+ }
+ }
+ for(int ctr_i = 1; ctr_i <= t2.Ncols(); ctr_i++){
+ for(int ctr_j = 1; ctr_j <= t2.Nrows(); ctr_j++){
+ if(t2(ctr_j,ctr_i)<=0){t2(ctr_j,ctr_i)=1;}
+ }
+ }
+ t1 = cumsum(sum(log(t1),2).AsRow());
+ t2 = cumsum(sum(log(t2),2).AsRow());
+
+ RowVector l_Az(d);
+ l_Az << (t1+t2);
+
+ RowVector l_lap;
+ l_lap = l_probU + l_nu +l_Az + l_lam + 0.5*log(2*M_PI)*(m+k)-0.5*log(N)*k;
+
+ RowVector l_BIC;
+ l_BIC = l_lam + l_nu - 0.5*log(N)*(m+k);
+
+ RowVector l_RRN;
+ l_RRN = -0.5*N*SP(k,log(SP(cumsum(eigenvalues),pow(k,-1))))+l_nu;
+
+ RowVector l_AIC;
+ l_AIC = -2*N*SP(tmp3,l_lhood)+ 2*(1+d*k+0.5*(k-1));
+ l_AIC = -l_AIC;
+
+ RowVector l_MDL;
+ l_MDL = -N*SP(tmp3,l_lhood)+ 0.5*(1+d*k+0.5*(k-1))*log(N);
+ l_MDL = -l_MDL;
+
+ Matrix Res;
+
+ Res = eigenvalues.t();
+ Res |= l_lap.t();
+ Res |= l_BIC.t();
+ Res |= l_MDL.t();
+ Res |= l_RRN.t();
+ Res |= l_AIC.t();
+
+
+ return Res;
+
+ }
+
+
+
+
+}
+
+
diff --git a/melpca.h b/melpca.h
new file mode 100644
index 0000000000000000000000000000000000000000..5fe727f6f7311a65aaeeeebf2d1770896701401c
--- /dev/null
+++ b/melpca.h
@@ -0,0 +1,58 @@
+ /* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melpca.h - PCA and whitening
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#ifndef __MELODICPCA_h
+#define __MELODICPCA_h
+#include "newimageall.h"
+#include "log.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melodic.h"
+//#include "melreport.h"
+#include "newmatap.h"
+#include "newmatio.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+
+namespace Melodic{
+
+ class MelodicReport;
+
+ class MelodicPCA
+ {
+ public:
+ MelodicPCA(MelodicData &pmelodat, MelodicOptions &popts, Log &plogger, MelodicReport &preport):
+ melodat(pmelodat),
+ opts(popts),
+ logger(plogger),
+ report(preport)
+ {
+ }
+
+ void perf_pca(const Matrix &Data);
+ void perf_white(const Matrix &Data);
+
+ private:
+ MelodicData &melodat;
+ MelodicOptions &opts;
+ Log &logger;
+ MelodicReport &report;
+
+ int pcadim();
+ RowVector cumsum(const RowVector& Inp);
+ Matrix ppca_est(const RowVector& eigenvalues, const int N);
+
+ RowVector Feta(int n1,int n2);
+ };
+}
+
+#endif
diff --git a/melreport.cc b/melreport.cc
new file mode 100644
index 0000000000000000000000000000000000000000..87fb7d82e2407f85a90aa6724a9cd902a71226f3
--- /dev/null
+++ b/melreport.cc
@@ -0,0 +1,614 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melreport.cc - report generation
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+#include "newimageall.h"
+#include "log.h"
+#include "melreport.h"
+
+namespace Melodic{
+
+ void MelodicReport::IC_rep(MelGMix &mmodel, int cnum, int dim)
+ {
+ if( bool(opts.genreport.value()) ){
+ addlink(mmodel.get_prefix()+".html",num2str(cnum));
+ IChtml.setDir(report.getDir(),mmodel.get_prefix()+".html");
+
+ {//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;
+ }
+
+
+ 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));
+
+
+ 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.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);
+
+ }
+ // }
+// 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;
+
+
+ {//plot time course
+ miscplot newplot;
+
+ if(opts.tr.value()>0.0)
+ newplot.timeseries(melodat.get_mix().Column(cnum).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_mix().Column(cnum).t(),
+ report.appendDir(string("t")+num2str(cnum)+".png"),
+ string("Timecourse (in TRs)"));
+ write_ascii_matrix(report.appendDir(string("t")
+ +num2str(cnum)+".txt"),
+ melodat.get_mix().Column(cnum));
+ 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_mix().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_mix().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_mix().Nrows())
+ +"* TR); period(s)=("
+ +num2str(melodat.get_mix().Nrows())
+ +"* TR)/cycles"));
+ write_ascii_matrix(report.appendDir(string("f")
+ +num2str(cnum)+".txt"),
+ melodat.get_mix().Column(cnum));
+ IChtml << "<A HREF=\"" << string("f")
+ +num2str(cnum)+".txt" << "\"> ";
+ IChtml << "<img BORDER=0 SRC=\""
+ +string("f")+num2str(cnum)+".png\"></A><p>" << endl;
+ }//frequency 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
+
+ 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();
+
+ //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);
+
+
+ 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;
+ } //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;
+
+ 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());
+ 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);
+
+ 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);
+ }
+ IChtml2 << "<a href=\""+mmodel.get_prefix()+".html\">";
+ IChtml2 << "<img BORDER=0 SRC=\""+ mmodel.get_prefix()+
+ ".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());
+
+ 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"));
+
+ 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;
+ }
+
+ {//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));
+
+ }
+ // 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;
+ } //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;
+ // }
+ }
+
+ { //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;
+ } //finish IC2 page
+ }
+ }
+
+
+ void MelodicReport::IC_simplerep(string prefix, int cnum, int dim)
+ {
+ if( bool(opts.genreport.value()) ){
+ addlink(prefix+".html",num2str(cnum));
+ IChtml.setDir(report.getDir(),prefix+".html");
+
+ {//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;
+ }
+
+
+ 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();
+
+ 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];
+
+ //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);
+ }
+
+ IChtml << "<img BORDER=0 SRC=\""+ prefix+
+ ".png\"><p>" << endl;
+
+ {//plot time course
+ miscplot newplot;
+
+ if(opts.tr.value()>0.0)
+ newplot.timeseries(melodat.get_mix().Column(cnum).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_mix().Column(cnum).t(),
+ report.appendDir(string("t")+
+ num2str(cnum)+".png"),
+ string("Timecourse (in TRs)"));
+ write_ascii_matrix(report.appendDir(string("t")
+ +num2str(cnum)+".txt"),
+ melodat.get_mix().Column(cnum));
+ 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_mix().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_mix().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_mix().Nrows())
+ +"* TR); period(s)=("
+ +num2str(melodat.get_mix().Nrows())
+ +"* TR)/cycles"));
+ write_ascii_matrix(report.appendDir(string("f")
+ +num2str(cnum)+".txt"),
+ melodat.get_mix().Column(cnum));
+ 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;
+ } //finish IC page
+ }
+ }
+
+ void MelodicReport::PPCA_rep(){
+
+ {//plot time course
+ report << "<p> <H3>PCA estimates </H3> <p>" << endl;
+
+ Matrix what;
+ miscplot newplot;
+
+ what = melodat.get_EV();
+ what &= melodat.get_EVP();
+
+ newplot.add_label("ordered Eigenvalues");
+ 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");
+ }
+
+ newplot.timeseries(what,
+ 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/melreport.h b/melreport.h
new file mode 100644
index 0000000000000000000000000000000000000000..eb6b0512ba205f05af51e98f47b993924c50da46
--- /dev/null
+++ b/melreport.h
@@ -0,0 +1,248 @@
+/* MELODIC - Multivariate exploratory linear optimized decomposition into
+ independent components
+
+ melreport.h - report generation
+
+ Christian F. Beckmann, FMRIB Image Analysis Group
+
+ Copyright (C) 1999-2002 University of Oxford */
+
+/* CCOPYRIGHT */
+
+
+#ifndef __MELODICREPORT_h
+#define __MELODICREPORT_h
+#include "newimageall.h"
+#include "log.h"
+#include "melpca.h"
+#include "meloptions.h"
+#include "meldata.h"
+#include "melgmix.h"
+#include "melodic.h"
+#include "newmatap.h"
+#include "newmatio.h"
+#include <time.h>
+#include <strstream>
+#include "miscplot.h"
+#include "miscpic.h"
+#include "options.h"
+
+using namespace Utilities;
+using namespace NEWIMAGE;
+using namespace MISCPLOT;
+using namespace MISCPIC;
+
+
+namespace Melodic{
+
+ class MelodicReport
+ {
+ public:
+ MelodicReport(MelodicData &pmelodat, MelodicOptions &popts, Log &plogger):
+ melodat(pmelodat),
+ opts(popts),
+ 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
+ << "<BODY BACKGROUND=\"file:" << getenv("FSLDIR")
+ << "/doc/images/fsl-bg.jpg\">" << endl
+ << endl << "<hr><CENTER> " << endl
+ << "<H1>MELODIC Report</H1>"
+ << endl
+ << report.getDir() << "/" << report.getLogFileName()
+ << "<br>" << endl
+ << ctime(&tmptime) << "<br>" << endl;
+ /* ifstream reptest(string("../report.log").c_str()); */
+/* if(!reptest) */
+/* report << "<A HREF=\"/" << logger.appendDir("report.log") */
+/* << "\">report.log</A><br>" << endl; */
+/* else */
+/* report << "<A HREF=\"" << logger.appendDir("melodic.log") */
+/* << "\">melodic.log</A><br>" << endl; */
+/* reptest.close(); */
+ report << "</CENTER>" << endl << endl
+ << "<hr><H2>Components:</H2> <p>" << endl << endl;
+ }
+ }
+
+ ~MelodicReport(){
+ if( bool(opts.genreport.value()) ){
+ report << "<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;
+ }
+ }
+
+ inline void analysistxt(){
+ if( bool(opts.genreport.value()) ){
+
+ report << "<hr><h2>Analysis methods</h2> <p>"<<endl;
+ report << "Analysis was carried out using MELODIC (Multivariate Exploratory Linear Decomposition into Independent Components) Version 2.00, 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 2002]."<<endl;
+
+ report << "The following melodic pre-processing was applied to the input data file: "<< endl;
+
+ if(opts.use_mask.value())
+ report << " masking of non-brain voxels;";
+
+ report << " voxel-wise de-meaning of the data;" << endl;
+
+ if(opts.varnorm.value())
+ report << " normalisation of the voxel-wise variance; ";
+
+ report << "<br>"<<endl;
+
+ report << " Pre-processed data was whitened and projected into a "
+ << melodat.get_mix().Ncols()<< "-dimensional subspace using ";
+ if(melodat.get_PPCA().Storage()>0){
+ report << "probabilistic Principal Component Analysis where the number of dimensions was estimated using ";
+ if(opts.pca_est.value() == string("lap"))
+ report << "the Laplace approximation to the Bayesian evidence of the model order [Beckmann 2001, Minka 2000]. " << endl;
+ else
+ if(opts.pca_est.value() == string("bic"))
+ report << "the <em> Bayesian Information Criterion</em> (BIC) [Kass 1993]. " << endl;
+ else
+ if(opts.pca_est.value() == string("mdl"))
+ report << "<em> Minimum Description Length</em> (MDL) [Rissanen 1978]. " << endl;
+ else
+ if(opts.pca_est.value() == string("aic"))
+ report << "the <em> Akaike Information Criterion</em> (AIC) [Akaike 1969]. " << endl;
+ else
+ report << "a committee of approximations to Bayesian the model order [Beckmann 2001, Minka 2000]. " << endl;
+
+ }
+ else
+ 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ärinen 1999]. " << endl;
+
+ report << "Estimated Component maps were divided by the standard deviation of the residual noise";
+
+ if(opts.perf_mm.value())
+ report << " and thresholded by fitting a mixture model to the histogram of intensity values [Beckmann 2002]. <p>" << endl;
+ else
+ report <<".<p>" << endl;
+
+ refstxt();
+ }
+ }
+
+ inline void refstxt(){
+
+ if( bool(opts.genreport.value()) ){
+ report << "<h3>References</h3> <p>"<<endl;
+
+ report << "[Hyvärinen 1999] A. Hyvärinen. Fast and Robust Fixed-Point Algorithms for Independent Component Analysis. IEEE Transactions on Neural Networks 10(3):626-634, 1999.<br> " << endl;
+ report << "[Beckmann 2002] C.F. Beckmann and S.M. Smith. Probabilistic Independent Component Analysis for FMRI. <A HREF=\"http://www.fmrib.ox.ac.uk/analysis/techrep/\">Technical Report TR02CB1</A>, FMRIB Analysis Group, 2002. <br>" << endl;
+
+ /* if(opts.perf_mm.value()){ */
+/* report << "[Bullmore 1996] E. Bullmore <em>et. al.</em> Statistical methods of estimation and inference for functional MR image analysis. Magnetic Resonance in Medicine, 35(2):261-177, 1996. <br>" << endl; */
+
+/* } */
+
+ 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 << "[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"))
+ 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
+ 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;
+ else
+ if(opts.pca_est.value() == string("mdl"))
+ report << "[Rissanen 1978]. J. Rissanen. Modelling by shortest data description. Automatica, 14:465-471, 1978. <br>" << endl;
+ else
+ 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;
+ 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;
+
+ }
+ }
+ }
+
+ inline void addtxt(string what){
+ if( bool(opts.genreport.value()) ){
+ report << what << endl;
+ }
+ }
+
+ inline void addpar(string what){
+ if( bool(opts.genreport.value()) ){
+ report << "<p>" << what << endl;
+ }
+ }
+
+ inline void addlink(string where, string what){
+ if( bool(opts.genreport.value()) ){
+ report << "<A HREF=\"" << where << "\"> " << what << "</A> ";
+ }
+ }
+
+ inline void addpic(string what, string link = ""){
+ if( bool(opts.genreport.value()) ){
+ if( link.length() > 0)
+ report << "<A HREF=\"" << link << "\"> ";
+
+ report << "<img BORDER=0 SRC=\"" << what<< ".png\"><p>";
+ if( link.length() > 0)
+ report << "</A> ";
+ }
+ }
+
+ inline string getDir(){
+ return report.getDir();
+ }
+
+ void IC_rep(MelGMix &mmodel, int cnum, int dim);
+ void IC_simplerep(string prefix, int cnum, int dim);
+
+ void PPCA_rep();
+
+ private:
+ MelodicData &melodat;
+ MelodicOptions &opts;
+ Log &logger;
+ Log report;
+
+ Log IChtml;
+ Log IChtml2;
+
+ void IC_rep_det(MelGMix &mmodel, int cnum, int dim);
+
+ string int2str(int n)
+ {
+ ostrstream os;
+ // os.fill(' ');
+ // os.width(width);
+ os.setf(ios::internal, ios::adjustfield);
+ os << n << '\0';
+ return os.str();
+ }
+
+
+ string float2str(float f, int width, int prec, int scientif)
+ {
+ ostrstream os;
+ int redw = int(std::abs(std::log10(std::abs(f))))+1;
+ if(width>0)
+ os.width(width);
+ if(scientif>0)
+ os.setf(ios::scientific);
+ os.precision(redw+std::abs(prec));
+ os.setf(ios::internal, ios::adjustfield);
+ os << f << '\0';
+ return os.str();
+ }
+ };
+
+}
+#endif
diff --git a/test.cc b/test.cc
new file mode 100644
index 0000000000000000000000000000000000000000..26406598591ed445a9ccf69bf32975a322d4ed48
--- /dev/null
+++ b/test.cc
@@ -0,0 +1,19 @@
+#include <iostream.h>
+#include <iomanip>
+#include "miscplot.h"
+#include "miscpic.h"
+
+
+
+using namespace std;
+
+template<class T>
+int foo(const T& bar){
+ return 1;
+}
+
+int main(int argc, char *argv[])
+{
+ cerr << foo(int(5));
+ return 1;
+}