From 79991dcf54961e60b8742b63cc366d503cf92a40 Mon Sep 17 00:00:00 2001
From: Mark Jenkinson <mark@fmrib.ox.ac.uk>
Date: Fri, 21 Feb 2003 16:22:32 +0000
Subject: [PATCH] Initial revision
---
meldata.cc | 500 ++++++++++++++++++++++++++++++++++++++++
meldata.h | 172 ++++++++++++++
melgmix.h | 195 ++++++++++++++++
melica.cc | 408 +++++++++++++++++++++++++++++++++
melica.h | 65 ++++++
melodic.cc | 426 ++++++++++++++++++++++++++++++++++
melodic.h | 34 +++
meloptions.cc | 247 ++++++++++++++++++++
meloptions.h | 344 ++++++++++++++++++++++++++++
melpca.cc | 586 +++++++++++++++++++++++++++++++++++++++++++++++
melpca.h | 58 +++++
melreport.cc | 614 ++++++++++++++++++++++++++++++++++++++++++++++++++
melreport.h | 248 ++++++++++++++++++++
test.cc | 19 ++
14 files changed, 3916 insertions(+)
create mode 100644 meldata.cc
create mode 100644 meldata.h
create mode 100644 melgmix.h
create mode 100644 melica.cc
create mode 100644 melica.h
create mode 100644 melodic.cc
create mode 100644 melodic.h
create mode 100644 meloptions.cc
create mode 100644 meloptions.h
create mode 100644 melpca.cc
create mode 100644 melpca.h
create mode 100644 melreport.cc
create mode 100644 melreport.h
create mode 100644 test.cc
diff --git a/meldata.cc b/meldata.cc
new file mode 100644
index 0000000..01da593
--- /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 0000000..8eadddf
--- /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 0000000..2b4755e
--- /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 0000000..1f355fb
--- /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 0000000..8cbce76
--- /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 0000000..13e6a1a
--- /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 0000000..c183726
--- /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 0000000..96e29c5
--- /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 0000000..8242c0d
--- /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 0000000..66a5f42
--- /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 0000000..5fe727f
--- /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 0000000..87fb7d8
--- /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 0000000..eb6b051
--- /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 0000000..2640659
--- /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;
+}
--
GitLab