From 7d4aa0417a4cda59923605bed100636dd283dc1c Mon Sep 17 00:00:00 2001
From: Mark Jenkinson <mark@fmrib.ox.ac.uk>
Date: Fri, 21 Feb 2003 16:48:47 +0000
Subject: [PATCH] Initial revision

---
 doc/index.html | 100 ++++++++++++++++++++++++++++++++++++
 melodicreport  | 135 +++++++++++++++++++++++++++++++++++++++++++++++++
 out.txt        | 135 +++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 370 insertions(+)
 create mode 100644 doc/index.html
 create mode 100755 melodicreport
 create mode 100644 out.txt

diff --git a/doc/index.html b/doc/index.html
new file mode 100644
index 0000000..4899123
--- /dev/null
+++ b/doc/index.html
@@ -0,0 +1,100 @@
+<HTML><HEAD><TITLE>MELODIC</TITLE>
+</HEAD><BODY BACKGROUND="../images/fsl-bg.jpg">
+<hr><TABLE BORDER=0 WIDTH="100%"><TR>
+<TD ALIGN=CENTER><H1>MELODIC - User Guide</H1>
+(Multivariate Exploratory Linear Optimized Decomposition into Independent Components)<br>
+ICA-Based Model-Free Analysis of 4D Data<br>
+MELODIC Version 2.0
+<TD ALIGN=RIGHT><a href="../index.html"><IMG BORDER=0 SRC="../images/fsl-logo.jpg"></a>
+</TR></TABLE>
+
+<HR><H2>INTRODUCTION</H2>
+ 
+<P>MELODIC (Multivariate Exploratory Linear Optimized Decomposition
+into Independent Components) uses ICA (Independent Component Analysis)
+to decompose 4D data (for example FMRI data) into different
+interesting spatial and temporal components. It can pick out different
+activation and artefactual components without any model being
+specified. Thus the output is a new 4D data set (each volume being a
+separate spatial component) and a text file of separate time courses
+(temporal components).
+
+<P>For more detail on MELODIC and an updated journal reference, see
+the <A
+HREF="http://www.fmrib.ox.ac.uk/analysis/research/melodic/">MELODIC
+research web page</A>. If you use MELODIC in your research, please
+quote the journal reference listed there.
+
+<p>The different MELODIC programs are:
+<UL>
+<LI><a href="#Melodic">Melodic</a> - MELODIC GUI<br>
+<LI><a href="#melodic">melodic</a> - command-line MELODIC program<br>
+</UL>
+
+<A NAME="Melodic"></A><hr><H2>Melodic GUI</H2>
+
+<center><IMG SRC="./mel1.png" ALIGN=CENTER></center>
+
+<p>Press <code>Select 4D data</code> and select the input data.
+
+<p><center><IMG SRC="./mel2.png" ALIGN=CENTER></center>
+
+<p>Low-frequency drifts and motion in the data can adversely affect
+the decomposition. In most cases, you would want to motion-correct the
+data and remove these drifts first. This can be done from within the
+Melodic GUI <code>Pre-stats</code> section.
+
+<center><IMG SRC="./mel3.png" ALIGN=CENTER></center>
+
+The MELODIC section finally lets you control some of the options for
+          the decomposition.
+The default setting will most probably already be set to what you
+          would want most of the time.
+
+<p> By default, Melodic will automatically estimate the number of
+          components from the data - you can switch this
+          option off and then can specify the number of
+          components explicitly.<p>
+Melodic will also by default carry out inference on the estimated maps
+          using the mixture model approach. A threshold level of 0.5
+          in the case of alternative hypothesis testing means that a
+          voxel 'survives' thresholding as soon as the probability
+          of being in the 'active' class (as modelled by the Gamma
+          densities) exceeds the probability of being in the
+          'background' noise class. This threshold level assumes that
+          you are placing an equal loss on false-positives and
+          false-negatives. If, however, you consider e.g. false-positives as
+          being twice as bad as false-negatives you should change this value
+          to 0.66...<p>
+Now simply press &nbsp;<code>Go</code><p>
+
+Melodic will then generate the results and 
+your terminal window will tell you where to find the web report. 
+
+
+
+Each webpage shows one spatial map thresholded and rendered on top of an example functional scan
+followed by the relevant time-course of the ICA decomposition and the
+            power-spectrum of the time-course. If you click on the
+            thresholded map, you can inspect the raw IC output
+            together with probability maps and the Mixture Model
+            fit. Use <code>previous</code> and <code>next</code> to
+            view another component map.<p>
+
+<A NAME="melodic"></A><HR><H2>melodic COMMAND-LINE PROGRAM</H2>
+
+<p>Type <b>melodic --help</b> to get usage.
+  
+
+<p><HR><FONT SIZE=1>Copyright &copy; 2001-2002, University of
+Oxford. Written by C. Beckmann and <A
+HREF="http://www.fmrib.ox.ac.uk/~steve/index.html">S. Smith</A>.</FONT>
+
+<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
+<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
+<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
+<br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br><br>
+
+</BODY></HTML>
+
+
diff --git a/melodicreport b/melodicreport
new file mode 100755
index 0000000..6d6b585
--- /dev/null
+++ b/melodicreport
@@ -0,0 +1,135 @@
+#!/bin/sh
+
+# the next line restarts using wish \
+exec $FSLDIR/bin/tclsh "$0" "$@"
+
+if { [ lindex $argv 0 ] == "" } {
+    puts "Usage: melodicreport <melodic_output_directory> \[background_image \]"
+    exit
+}
+
+set background 0
+if { [ lindex $argv 1 ] != "" } {
+    set background 1
+    set bgimage [ lindex $argv 1 ]
+}
+
+set S 0.2
+catch {  [ expr $S ] } errrr 
+
+set FSLDIR $env(FSLDIR)
+
+set MD [ lindex $argv 0 ]
+catch { exec sh -c "mkdir ${MD}/report" }
+
+if { ! [ file exists ${MD}/melodic_IC.hdr ] || ! [ file exists ${MD}/melodic_IC.img ] } {
+    puts "No MELODIC output melodic_IC.hdr/.img - exiting."
+    exit
+}
+
+set REPORTFILE ${MD}/report/report.html
+set report [ open $REPORTFILE "w" ]
+fconfigure $report -buffering none
+
+puts $report "<HTML>
+
+<TITLE>MELODIC Report</TITLE>
+
+<BODY BACKGROUND=\"file:${FSLDIR}/doc/images/fsl-bg.jpg\">
+
+<hr><CENTER>
+<H1>MELODIC Report</H1>
+$REPORTFILE<br>
+[ exec date ]
+</CENTER>
+
+<hr><H2>Components:</H2><p>"
+
+set dim  [ exec sh -c "${FSLDIR}/bin/avwval ${MD}/melodic_IC dim4" ]
+set npts [ expr [ exec sh -c "cat ${MD}/melodic_mix | wc -l" ] * 1 ]
+set xscale [ expr $npts * 4 ]
+if { $xscale < 750 } {
+    set xscale 750
+}
+if { $xscale > 3000 } {
+    set xscale 3000
+}
+set xscale [ expr $xscale / 750.0 ]
+
+for { set i 1 } { $i <= $dim } { incr i 1 } {
+
+    puts $report "<a href=\"c${i}.html\">$i</a>"
+ 
+    exec sh -c "${FSLDIR}/bin/avwroi ${MD}/melodic_IC ${MD}/grot [ expr $i - 1 ] 1"
+    catch { exec sh -c "${FSLDIR}/bin/avwstats ${MD}/grot -R" } minmax
+    set absmin [ expr abs ( [ lindex $minmax 0 ] ) ]
+    set absmax [ expr abs ( [ lindex $minmax 1 ] ) ]
+    if { $absmin > $absmax } {
+	    set absmax $absmin
+    }
+    set absmax [ expr $absmax * 0.75 ]
+
+    if { $background == 0 } {
+	exec sh -c "${FSLDIR}/bin/slicer ${MD}/grot -l ${FSLDIR}/etc/luts/render3.lut -s 2 -i -$absmax $absmax -A 950 ${MD}/grot.ppm"
+    } else {
+	set llimit [ expr $absmax / 6.0 ]
+	exec sh -c "${FSLDIR}/bin/avwmaths ${MD}/grot -mul -1 ${MD}/grotinv"
+	exec sh -c "${FSLDIR}/bin/overlay 1 0 $bgimage -a ${MD}/grot $llimit $absmax ${MD}/grotinv $llimit $absmax ${MD}/grotoverlay"
+	exec sh -c "${FSLDIR}/bin/slicer ${MD}/grotoverlay -s 2 -A 950 ${MD}/grot.ppm"
+    }
+
+    exec sh -c "${FSLDIR}/bin/convert ${MD}/grot.ppm ${MD}/report/s$i.gif"
+
+    exec sh -c "cat ${MD}/melodic_mix | awk '{ print \$$i }' > ${MD}/report/t${i}.txt"
+    exec sh -c "echo \"set term pbm color ; set size $xscale,0.25 ; set title 'Timecourse (TRs)' ; plot '${MD}/report/t${i}.txt' title '' with lines\" | ${FSLDIR}/bin/gnuplot > ${MD}/grot.ppm"
+    exec sh -c "${FSLDIR}/bin/convert ${MD}/grot.ppm ${MD}/report/t$i.gif"
+
+    exec sh -c "cat ${MD}/melodic_FTmix | awk '{ print NR \" \" \$$i }' > ${MD}/report/f${i}.txt"
+    exec sh -c "echo \"set term pbm color ; set size $xscale,0.25 ; set title 'FT of Timecourse (cycles)  frequency(Hz)=cycles/($npts*TR)  period(s)=($npts*TR)/cycles' ; plot '${MD}/report/f${i}.txt' title '' with lines\" | ${FSLDIR}/bin/gnuplot > ${MD}/grot.ppm"
+    exec sh -c "${FSLDIR}/bin/convert ${MD}/grot.ppm ${MD}/report/f$i.gif"
+
+    set creport [ open ${MD}/report/c${i}.html "w" ]
+
+    puts $creport "<HTML>
+    <TITLE>MELODIC Component $i</TITLE>
+    <BODY BACKGROUND=\"file:${FSLDIR}/doc/images/fsl-bg.jpg\">
+    <hr><CENTER><H1>MELODIC Component $i</H1>
+    "
+
+    if { $i > 1 } {
+	puts $creport "<a href=\"c[ expr $i - 1 ].html\">previous</a> - "
+    }
+    
+    puts $creport "<a href=\"report.html\">index</a>"
+    
+    if { $i < $dim } {
+	puts $creport " - <a href=\"c[ expr $i + 1 ].html\">next</a>"
+    }
+
+    puts $creport "
+    <hr><p>
+
+    <IMG BORDER=0 SRC=\"s$i.gif\"><p>
+    <a href=\"t$i.txt\"><IMG BORDER=0 SRC=\"t$i.gif\"></a><p>
+    <a href=\"f$i.txt\"><IMG BORDER=0 SRC=\"f$i.gif\"></a><p>
+
+    <hr></CENTER></BODY></HTML>
+    "
+    close $creport
+}
+
+exec sh -c "/bin/rm -f ${MD}/grot*"
+
+puts $report "<HR><FONT SIZE=1>This page produced automatically by MELODIC - a part of <A HREF=\"http://www.fmrib.ox.ac.uk/fsl\">FSL - FMRIB's Software Library</A>.</FONT>
+
+</BODY></HTML>
+"
+
+close $report
+
+puts "Finished melodicreport at [ exec date ]
+To view the MELODIC report point your web browser at
+$REPORTFILE
+"
+
+exit
diff --git a/out.txt b/out.txt
new file mode 100644
index 0000000..02559be
--- /dev/null
+++ b/out.txt
@@ -0,0 +1,135 @@
+-------------------------------------
+meldata.h
+-------------------------------------
+
+
+53,55d112
+<    
+<       inline Matrix& get_pcaEVN() {return pcaEVN;}
+<       inline void set_pcaEVN(Matrix& Arg) {pcaEVN = Arg;}
+59,61c116,118
+<       
+<       inline DiagonalMatrix& get_pcaDVN() {return pcaDVN;}
+<       inline void set_pcaDVN(DiagonalMatrix& Arg) {pcaDVN = Arg;}
+---
+> 
+>       inline Matrix& get_data() {return Data;}
+>       inline void set_data(Matrix& Arg) {Data = Arg;}
+96,100c153,154
+<       //      inline Matrix& get_DataVN() {
+<       //	Matrix tmp;
+<       //      tmp = SP(Data,ones(Data.Nrows(),1)*stdDevi);
+<       //	return tmp;
+<       //      }
+---
+>       inline Matrix& get_DataVN() {return DataVN;}
+>       inline void set_DataVN(Matrix& Arg) {DataVN = Arg;}
+123,131d176
+<       inline void calc_Corr() { Corr = cov(Data.t());}
+<       inline void calc_CorrVN() { 
+< 	CorrVN = cov(SP(Data,ones(Data.Nrows(),1)*stdDevi).t());
+<       }
+< 
+<       inline SymmetricMatrix& get_Corr() {return Corr;}
+<       inline SymmetricMatrix& get_CorrVN() {return CorrVN;}
+< 
+<  
+157,162d201
+<       Matrix pcaEVN;
+<       DiagonalMatrix pcaDVN;
+< 
+<       SymmetricMatrix Corr;
+<       SymmetricMatrix CorrVN;
+<       
+178a218
+>       Matrix DataVN;
+
+
+-------------------------------------
+meldata.cc
+-------------------------------------
+
+
+103d162
+< 
+107,108c166,167
+<       //Corr = cov(Data.t());
+<       calc_Corr();
+---
+>       SymmetricMatrix Corr;
+>       Corr = cov(Data.t());
+128c187
+< 
+---
+>       
+132,145c191,192
+< 
+<     if(opts.segment.value().length()>0){
+<       Matrix Res;
+<       Res = ones(Data.Nrows(),1)*RXweight;
+<       if(opts.varnorm.value()){
+< 	Res = SP(SP(Data,ones(Data.Nrows(),1)*stdDevi),Res);
+<       }
+<       else
+< 	Res = SP(Data,Res);
+<       Data = Res;
+<     }
+< 
+<     calc_CorrVN();
+< 
+---
+>     
+>     DataVN = SP(Data,ones(Data.Nrows(),1)*stdDevi);
+147,148c194
+<       Data = SP(Data,ones(Data.Nrows(),1)*stdDevi);
+<       Corr = CorrVN;
+---
+>       Data = DataVN;
+151,153d196
+<     else
+<       calc_Corr();
+< 
+
+-------------------------------------
+melodic.cc
+-------------------------------------
+
+92c152
+< 	pcaobj.perf_pca(melodat.get_Data());
+---
+> 	pcaobj.perf_pca(melodat.get_DataVN());
+
+-------------------------------------
+melpca.cc
+-------------------------------------
+
+
+30a91,92
+>      SymmetricMatrix Corr;
+> 
+34,38c96,99
+<        if(opts.varnorm.value()){
+< 	 Res = SP(SP(melodat.get_Data(),ones(melodat.get_Data().Nrows(),1)*melodat.get_stdDevi()),Res);
+<        }
+<        else
+< 	 Res = SP(melodat.get_Data(),Res);
+---
+>        Res = SP(melodat.get_DataVN(),Res);
+>        melodat.set_DataVN(Res);
+>        Res = ones(Data.Nrows(),1)*melodat.get_RXweight();
+>        Res = SP(melodat.get_Data(),Res);
+40,41d100
+<        melodat.calc_Corr();
+<        melodat.calc_CorrVN();
+43a103,104
+>      Corr = cov(melodat.get_Data().t());
+> 
+47c108
+<      EigenValues(melodat.get_Corr(),tmpD,tmpE);
+---
+>      EigenValues(Corr,tmpD,tmpE);
+170c231
+<     Corr = cov(SP(melodat.get_Data(),ones(melodat.get_Data().Nrows(),1)*melodat.get_stdDevi()).t());
+---
+>     Corr = cov(melodat.get_DataVN().t());
+
-- 
GitLab