From d80f75561830cf35220694a425ba18d92e02d2fc Mon Sep 17 00:00:00 2001
From: Christian Beckmann <c.beckmann@donders.ru.nl>
Date: Mon, 23 Aug 2004 14:10:14 +0000
Subject: [PATCH] updated to ln(11)
---
meldata.cc | 20 ++++---
meldata.h | 7 ++-
melgmix.cc | 2 +-
melgmix.h | 12 +++-
melodic.cc | 35 ++++++-----
melodic.h | 2 +-
meloptions.cc | 1 +
meloptions.h | 5 ++
melreport.h | 2 +-
test.cc | 163 +++++++++++++++++++++++++++++++++++++++++++++++---
10 files changed, 213 insertions(+), 36 deletions(-)
diff --git a/meldata.cc b/meldata.cc
index cdcae29..c8c3b93 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -196,12 +196,12 @@ namespace Melodic{
message("Writing results to : " << endl);
//Output IC
- if((IC.Storage()>0)&&(opts.output_origIC.value()) ){
+ if((IC.Storage()>0)&&(opts.output_origIC.value())&&(after_mm==false)){
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);
+ + "_oIC"),tempInfo);
message(" " << logger.appendDir(opts.outputfname.value() + "_oIC") <<endl);
}
@@ -214,13 +214,17 @@ namespace Melodic{
//Matrix ICadjust = IC;
Matrix ICadjust;
- stdNoisei = pow(stdev(Data - mixMatrix * IC)*std::sqrt((float)(Data.Nrows()-1))/
- std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
+ if(after_mm)
+ ICadjust = IC;
+ else{
+ stdNoisei = pow(stdev(Data - mixMatrix * IC)*std::sqrt((float)(Data.Nrows()-1))/
+ std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
+
+ ColumnVector diagvals;
+ diagvals=pow(diag(unmixMatrix*unmixMatrix.t()),-0.5);
+ ICadjust = SP(IC,diagvals*stdNoisei);
+ }
- ColumnVector diagvals;
- diagvals=pow(diag(unmixMatrix*unmixMatrix.t()),-0.5);
- ICadjust = SP(IC,diagvals*stdNoisei);
-
tempVol.setmatrix(ICadjust,Mask);
//strncpy(tempInfo.header.hist.aux_file,"render3",24);
save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
diff --git a/meldata.h b/meldata.h
index f3cd6e0..491cfa9 100644
--- a/meldata.h
+++ b/meldata.h
@@ -30,7 +30,8 @@ namespace Melodic{
MelodicData(MelodicOptions &popts, Log &plogger):
opts(popts),
logger(plogger)
- {
+ {
+ after_mm = false;
}
void save();
@@ -118,6 +119,8 @@ namespace Melodic{
inline float get_resels() {return Resels;}
inline void set_resels(float& Arg) {Resels = Arg;}
+ inline void set_after_mm(bool val) {after_mm = val;}
+
Matrix smoothColumns(const Matrix &inp);
inline void flipres(int num){
@@ -162,6 +165,8 @@ namespace Melodic{
Matrix DataVN;
Matrix PPCA;
+ bool after_mm;
+
float Resels;
char Mean_fname[1000];
diff --git a/melgmix.cc b/melgmix.cc
index 9b142cd..cab5e93 100644
--- a/melgmix.cc
+++ b/melgmix.cc
@@ -730,7 +730,7 @@ namespace Melodic{
// probmap << normcdf(data,means(1),vars(1));
// }|| std::abs(mean(data,2).AsScalar()) > 1.0
- if(props(1)<0.6 ){
+ if(props(1)<0.4 ){
//set up GMM
message(" try Gaussian Mixture Model " << endl);
props=zeros(1,nummix);
diff --git a/melgmix.h b/melgmix.h
index 71891ac..b6fffae 100644
--- a/melgmix.h
+++ b/melgmix.h
@@ -53,7 +53,6 @@ namespace Melodic{
volume<float> themean, int num_mix = 3,
float eps = 0.0, bool fixdim = false);
-
void gmmfit();
void ggmfit();
@@ -123,7 +122,7 @@ namespace Melodic{
threshmaps = -threshmaps;
if(mmtype=="GGM"){
float tmp;
- tmp= means(2);means(2)= means(3);means(3)=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;
}
@@ -150,6 +149,15 @@ namespace Melodic{
threshinfo.clear();
}
+ inline void smooth_probs(float howmuch){
+ volume4D<float> tempVol;
+ tempVol.setmatrix(probmap,Mask);
+ tempVol[0]= smooth(tempVol[0],howmuch);
+ probmap = tempVol.matrix(Mask);
+ }
+
+
+
double datamean;
double datastdev;
diff --git a/melodic.cc b/melodic.cc
index f59a633..c428ab0 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -327,12 +327,6 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
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());
//re-scale spatial maps to mean 0 for nht
@@ -358,6 +352,17 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
melodat.set_IC(ctr,tmp);
}
+ if(opts.smooth_probmap.value()<0.0){
+ message(" smoothing probability map ... "<< endl);
+ mixmod.smooth_probs(0.5*(std::min(std::min(std::abs(melodat.get_mean().xdim()),std::abs(melodat.get_mean().ydim())),std::abs(melodat.get_mean().zdim()))));
+ }
+
+ if(opts.smooth_probmap.value()>0.0){
+ message(" smoothing probability map ... "<< endl);
+ mixmod.smooth_probs(opts.smooth_probmap.value());
+ }
+
+
message(" thresholding ... "<< endl);
mixmod.threshold(opts.mmthresh.value());
@@ -418,16 +423,18 @@ Matrix mmall(Log& logger, MelodicOptions& opts,
//now safe new data
// bool what = opts.verbose.value();
//opts.verbose.set_T(false);
- // melodat.save();
+
+ melodat.set_after_mm(TRUE);
+ melodat.save();
//opts.verbose.set_T(what);
- if(melodat.get_IC().Storage()>0){
- volume4D<float> tempVol;
- tempVol.setmatrix(melodat.get_IC(),melodat.get_mask());
- save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
- + "_IC"),melodat.tempInfo);
- message(endl<< endl << " Saving " << logger.appendDir(opts.outputfname.value() + "_IC") <<endl);
- }
+ //if(melodat.get_IC().Storage()>0){
+ // volume4D<float> tempVol;
+ // tempVol.setmatrix(melodat.get_IC(),melodat.get_mask());
+ // save_volume4D(tempVol,logger.appendDir(opts.outputfname.value()
+ // + "_IC"),melodat.tempInfo);
+ // message(endl<< endl << " Saving " << logger.appendDir(opts.outputfname.value() + "_IC") <<endl);
+ //}
if( bool(opts.genreport.value()) ){
message(endl << endl <<
diff --git a/melodic.h b/melodic.h
index afc7e77..78822d5 100644
--- a/melodic.h
+++ b/melodic.h
@@ -32,7 +32,7 @@
namespace Melodic{
-const string version = "ln(10)";
+const string version = "ln(11)";
}
diff --git a/meloptions.cc b/meloptions.cc
index 7087b46..eb9eb40 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -234,6 +234,7 @@ void MelodicOptions::status()
cout << " repeats = " << repeats.value() << endl;
cout << " nlconst1 = " << nlconst1.value() << endl;
cout << " nlconst2 = " << nlconst2.value() << endl;
+ cout << " smooth_probmap = " << smooth_probmap.value() << endl;
cout << " epsilon = " << epsilon.value() << endl;
cout << " threshold = " << threshold.value() << endl;
diff --git a/meloptions.h b/meloptions.h
index ecd3482..140f802 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -93,6 +93,7 @@ class MelodicOptions {
Option<int> repeats;
Option<float> nlconst1;
Option<float> nlconst2;
+ Option<float> smooth_probmap;
Option<bool> remove_meanvol;
@@ -258,6 +259,9 @@ class MelodicOptions {
nlconst2(string("--nl2,--nlconst2"), 1.0,
string("nonlinear constant 2"),
false, requires_argument, false),
+ smooth_probmap(string("--smooth_pm"), 0.0,
+ string("width of smoothing kernel for probability maps"),
+ false, requires_argument, false),
remove_meanvol(string("--keepmeanvol"), true,
string("do not subtract mean volume"),
false, no_argument, false),
@@ -331,6 +335,7 @@ class MelodicOptions {
options.add(repeats);
options.add(nlconst1);
options.add(nlconst2);
+ options.add(smooth_probmap);
options.add(remove_meanvol);
options.add(remove_endslices);
options.add(rescale_nht);
diff --git a/melreport.h b/melreport.h
index acb987b..2ab3a79 100644
--- a/melreport.h
+++ b/melreport.h
@@ -84,7 +84,7 @@ namespace Melodic{
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 2004]."<<endl;
+ report << "Analysis was carried out using MELODIC (Multivariate Exploratory Linear Decomposition into Independent Components) Version ln(11), part of FSL (FMRIB's Software Library, <A HREF=\"http://www.fmrib.ox.ac.uk/fsl/\">www.fmrib.ox.ac.uk/fsl</A>), an implementation for the estimation of a Probabilistic Independent Component Analysis model [Beckmann 2004]."<<endl;
report << "The following melodic pre-processing was applied to the input data file: "<< endl;
diff --git a/test.cc b/test.cc
index 3c8339f..bc2c81a 100644
--- a/test.cc
+++ b/test.cc
@@ -10,6 +10,8 @@
#include "miscmaths/miscprob.h"
#include <string>
#include "utils/log.h"
+#include "melgmix.h"
+#include "meloptions.h"
using namespace std;
@@ -17,7 +19,7 @@ using namespace Utilities;
using namespace NEWIMAGE;
using namespace MISCPLOT;
using namespace MISCPIC;
-
+using namespace Melodic;
Matrix calcFFT(const Matrix& Mat)
@@ -56,12 +58,153 @@ string float2str(float f, int width, int prec, int scientif)
return os.str();
}
+void usage(void)
+{
+ cout << "Usage: test infile ICfile [mixfile]" << endl;
+ exit(1);
+}
+
int main(int argc, char *argv[])
{
- Matrix data;
+ if (argc<3)
+ usage();
+
+ Matrix ICs;
+ Matrix mixMatrix;
+ Matrix fmixMatrix;
+ volumeinfo ICvolInfo;
+ volume<float> Mask;
+ volume<float> Mean;
+
+ string RAWfname;
+ RAWfname = string(argv[1]);
+ string ICfname;
+ ICfname = string(argv[2]);
+
+ string MIXfname;
+ if (argc>3)
+ MIXfname = string(argv[3]);
+
+ cerr << argc << " " << RAWfname << " " <<ICfname << " " << MIXfname << endl;
+
+ {
+ volume4D<float> RawData;
+ cout << " Reading orig. data " << RAWfname << " ... ";
+ read_volume4D(RawData,RAWfname,ICvolInfo);
+ Mean = meanvol(RawData);
+
+ float howmuch = 0.5*(std::min(std::min(std::abs(Mean.xdim()),std::abs(Mean.ydim())),std::abs(Mean.zdim())));
+
+ cout << " Smoothing by " << howmuch << endl;
+ volume<float> tmpvol = smooth(Mean,howmuch);
+ miscpic newpic;
+ char instr[10000];
+
+ sprintf(instr," ");
+ strcat(instr,"-s 2");
+ strcat(instr," -A 950 ");
+ strcat(instr,string("./res/m1.png").c_str());
+ newpic.slicer(Mean, instr, &ICvolInfo);
+
+ char instr2[10000];
+
+ sprintf(instr2," ");
+ strcat(instr2,"-s 2");
+ strcat(instr2," -A 950 ");
+ strcat(instr2,string("./res/m2.png").c_str());
+ newpic.slicer(tmpvol, instr2, &ICvolInfo);
+
+
+
+ cout << " done! " << endl;
+ }
+ {
+ volume4D<float> RawIC;
+ cout << "Reading components " << ICfname << " ... ";
+ read_volume4D(RawIC,ICfname);
+ cout << " done" << endl;
+
+ cout << "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;
+ cout << " done" << endl;
+ }
+
+ if(MIXfname.length()>0){
+ cout << "Reading mixing matrix " << MIXfname << " ... ";
+ mixMatrix = read_ascii_matrix(MIXfname);
+ if (mixMatrix.Storage()<=0) {
+ cerr <<" Please specify the mixing matrix correctly" << endl;
+ exit(2);
+ cout << " done " << endl;
+ }
+ }
+
+ cout << " ICs: " << ICs.Nrows() << " x " << ICs.Ncols() << endl;
+ if(ICs.Nrows()>0){
+ cout << " Plotting histogram for map 1" << endl;
+
+ miscplot newplot;
+ // newplot.add_label("legend1");
+ //newplot.add_label("legend2");
+ //newplot.add_ylabel(string("yll1"));
+
+ //newplot.add_xlabel(string("xl1"));
+ // newplot.set_xysize(600,600);
+
+ Matrix mu(1,3);
+ Matrix pi(1,3);
+ Matrix std(1,3);
+ mu(1,1)=0;mu(1,2)=2.10;mu(1,3)=-1.99388;
+ pi(1,1)=0.911806;pi(1,2)=0.066722;pi(1,3)=0.021472;
+ std(1,1)=1;std(1,2)=0.831590;std(1,3)=0.639299;
+ newplot.gmmfit(ICs.Row(1),mu,std,pi,string("./res/g1.png"),string("Title"));
+
+ mu(1,1)=0;mu(1,2)=2.209024;mu(1,3)=-2.016637;
+ pi(1,1)=0.873335;pi(1,2)=0.086185;pi(1,3)=0.040480;
+ std(1,1)=1;std(1,2)=1.012970;std(1,3)=0.616374;
+ newplot.gmmfit(ICs.Row(2),mu,std,pi,string("./res/g2.png"),string("Title"));
+
+ mu(1,1)=0;mu(1,2)=2.694458;mu(1,3)=-3.418763;
+ pi(1,1)=0.9779397;pi(1,2)=0.017768;pi(1,3)=0.004293;
+ std(1,1)=1;std(1,2)=2.848577;std(1,3)=0.385187;
+ newplot.gmmfit(ICs.Row(3),mu,std,pi,string("./res/g3.png"),string("Title"));
+
+ // newplot.histogram(ICs.Row(1),string("./res/h1.png"),string("Title"));
+ }
+
+ cout << endl << endl;
+ /*
+ Matrix data;
+
cerr<< "Reading mixing matrix ";
data = read_ascii_matrix(string("mixmat"));
if (data.Storage()<=0) {
@@ -77,7 +220,7 @@ int main(int argc, char *argv[])
cerr << " create new data " << M_PI <<endl;
- ColumnVector X(200);
+ ColumnVector X(1000);
for(int i=1; i<=X.Nrows(); i++){
X(i) = (float((i-1))/(X.Nrows()-1))*M_PI;
@@ -113,19 +256,23 @@ mmat.setDir(string("./res/"),string("mix.txt"));
<< "</H1>"<< endl;
- float tr = 0.0;
+ float tr = 5.0;
Matrix data2;
data2 = calcFFT(data);
-
for (int i=1; i<=data.Ncols(); i++)
{
- cerr << i << " " ;
{//plot time course
miscplot newplot;
-
+ //newplot.add_label("legend1");
+ //newplot.add_label("legend2");
+ // newplot.add_ylabel(string("yll1"));
+ //newplot.add_ylabel(string("ylabel2guiygyiasyugfuyigfuigasgfuiguiyasgfuigasuigfagifgiuygasisgfuigasgfgasgfigiagfgaggfigaif"));
+ //newplot.add_xlabel(string("xl1"));
+ //newplot.add_xlabel(string("xlabel2"));
+ //newplot.set_xysize(600,200);
if(tr>0.0)
newplot.timeseries(data.Column(i).t(),
string("./res/t")+num2str(i)+".png",
@@ -183,6 +330,6 @@ mmat.setDir(string("./res/"),string("mix.txt"));
IChtml<< "<HR><FONT SIZE=1>This page produced automatically by "
<< "FMRIB Software Library</A>.</FONT>" << endl
<< "</BODY></HTML>" << endl;
- cerr << endl;
+ cerr << endl; */
return 0;
}
--
GitLab