diff --git a/meldata.cc b/meldata.cc
index 18bdc4c42ca582bf0db12af13f4cbc646b05ed6e..ce24816b0eccacceb8e6079f79bd02dd94cc0fc9 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -518,16 +518,23 @@ namespace Melodic{
//update mask
if(opts.update_mask.value()){
- message("Excluding voxels with constant value ...");
+ message(endl<< "Excluding voxels with constant value ...");
update_mask(Mask, Data);
message(" done" << endl);
}
-
Matrix tmp = IdentityMatrix(Data.Nrows());
DWM.push_back(tmp);
WM.push_back(tmp);
opts.varnorm.set_T(tmpvarnorm);
+
+ if(opts.varnorm2.value()){
+ message(" Normalising by voxel-wise variance ...");
+ stdDev = varnorm(Data,std::min(30,Data.Nrows()-1),
+ opts.vn_level.value(), opts.econ.value());
+ message(" done" << endl);
+ }
+
dbgmsg(string("END: setup_migp") << endl);
}
@@ -561,7 +568,7 @@ namespace Melodic{
save_volume(Mask,logger.appendDir("mask"));
save_volume(Mean,logger.appendDir("mean"));
- dbgmsg(string("START: setup") << endl);
+ dbgmsg(string("END: setup") << endl);
} // void setup()
void MelodicData::setup_misc()
diff --git a/meloptions.h b/meloptions.h
index 65e17cd4ce5a5a248ad3056dd9f875c7af681469..44fd256ed422a36c185824e13e334991d3b7e589 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -66,6 +66,7 @@ class MelodicOptions {
Option<string> nonlinearity;
Option<bool> varnorm;
+ Option<bool> varnorm2;
Option<bool> pbsc;
Option<bool> pspec;
Option<string> segment;
@@ -195,7 +196,7 @@ class MelodicOptions {
joined_whiten(string("--sep_whiten"), false,
string("switch on separate whitening"),
false, no_argument, false),
- joined_vn(string("--sep_vn"), true,
+ joined_vn(string("--sep_vn"), false,
string("switch on group variance nomalisation (as opposed to separate VN)"),
false, no_argument, false),
dr_pca(string("--mod_pca"), true,
@@ -234,6 +235,9 @@ class MelodicOptions {
varnorm(string("--vn,--varnorm"), true,
string("switch off variance normalisation"),
false, no_argument),
+ varnorm2(string("--vn2"), true,
+ string("switch off 2nd level variance normalisation"),
+ false, no_argument, false),
pbsc(string("--pbsc"), false,
string(" switch on conversion to percent BOLD signal change"),
false, no_argument, false),
@@ -439,6 +443,7 @@ class MelodicOptions {
options.add(approach);
options.add(nonlinearity);
options.add(varnorm);
+ options.add(varnorm2);
options.add(pbsc);
options.add(pspec);
options.add(segment);
diff --git a/test.cc b/test.cc
index 73c75585b04407dd1d2fb1a0aef7c98e0dcc9a49..648574af43909f9407aa491fa2cb2483e17e7a7e 100644
--- a/test.cc
+++ b/test.cc
@@ -11,9 +11,10 @@
#include "miscmaths/miscprob.h"
#include "utils/options.h"
#include <vector>
-#include <time.h>
+#include <ctime>
#include "newimage/newimageall.h"
#include "melhlprfns.h"
+#include <iostream>
#ifdef __APPLE__
#include <mach/mach.h>
@@ -32,13 +33,16 @@
}
#endif
-
-// a simple message macro that takes care of cout
+// a simple message macro that takes care of cout and log
#define message(msg) { \
cout << msg; \
cout.flush(); \
}
+#define outMsize(msg,Mat) { \
+ cerr << " " << msg << " " <<Mat.Nrows() << " x " << Mat.Ncols() << endl; \
+}
+
using namespace MISCPLOT;
using namespace MISCMATHS;
@@ -48,7 +52,7 @@ using namespace Melodic;
// GLOBALS
-time_t sttime;
+clock_t tictime;
// The two strings below specify the title and example usage that is
@@ -63,19 +67,22 @@ time_t sttime;
//Command line Options {
Option<string> fnin(string("-i,--in"), string(""),
string("input file name (matrix 3D or 4D image)"),
- true, requires_argument);
+ false, requires_argument);
Option<string> fnmask(string("-m"), string(""),
string("mask file name "),
- true, requires_argument);
+ false, requires_argument);
Option<int> help(string("-h,--help"), 0,
string("display this help text"),
false,no_argument);
- Option<int> xdim(string("-x,--xdim"), 10,
+ Option<int> xdim(string("-x,--xdim"), 0,
string("xdim"),
false,requires_argument);
- Option<int> ydim(string("-y,--ydim"), 10,
+ Option<int> ydim(string("-y,--ydim"), 0,
string("ydim"),
false,requires_argument);
+ Option<int> econ(string("-e,--econ"), 0,
+ string("econ: how to liump stuff"),
+ false,requires_argument);
/*
}
*/
@@ -84,81 +91,106 @@ time_t sttime;
// Local functions
void tic(){
- sttime = time(NULL);
+ tictime = clock();
}
void toc(){
-
- time_t current_time;
- current_time = time(NULL) - sttime;
- cerr << endl << "TOC: " << string(ctime(¤t_time)) << endl;
+ cerr << endl << "TOC: " << float(clock()-tictime)/CLOCKS_PER_SEC << " seconds" << endl<<endl;
}
-Matrix calccorr(const Matrix& in, bool econ)
- {
- Matrix Res;
- Res = zeros(in.Nrows(),in.Nrows());
- if(econ){
- ColumnVector tmp;
- for(int ctr=1; ctr <= in.Ncols(); ctr++){
- tmp = in.Column(ctr);
- tmp = tmp - mean(tmp).AsScalar();
- Res += (tmp * tmp.t()) / in.Ncols();
- }
- }
- else
- Res = cov(in.t());
- return Res;
- } //Matrix calc_corr
+Matrix calccorr(const Matrix& in, int econ)
+ {
+ Matrix Res;
+ int nrows=in.Nrows();
+ int ncols=in.Ncols();
+ Res = zeros(nrows,nrows);
+
+ if(econ>0){
+ RowVector colmeans(ncols);
+ for (int n=1; n<=ncols; n++) {
+ colmeans(n)=0;
+ for (int m=1; m<=nrows; m++) {
+ colmeans(n)+=in(m,n);
+ }
+ colmeans(n)/=nrows;
+ }
+ int dcol = econ;
+ Matrix suba;
+
+ for(int ctr=1; ctr <= in.Ncols(); ctr+=dcol){
+ suba=in.SubMatrix(1,nrows,ctr,Min(ctr+dcol-1,ncols));
+ int scolmax = suba.Ncols();
+
+ for (int n=1; n<=scolmax; n++) {
+ double cmean=colmeans(ctr + n - 1);
+ for (int m=1; m<=nrows; m++) {
+ suba(m,n)-=cmean;
+ }
+ }
+
+ Res += suba*suba.t() / ncols;
+ }
+ }
+ else
+ Res = cov(in.t());
+ return Res;
+ } //Matrix calccorr
int do_work(int argc, char* argv[]) {
- sttime = time(NULL);
+ tic();
Matrix MatrixData;
volume<float> Mean;
-{
- volume4D<float> RawData;
- volume<float> theMask;
- toc();
- //read data
- message("Reading data file " << (string)fnin.value() << " ... ");
- read_volume4D(RawData,fnin.value());
- message(" done" << endl);
+ if(xdim.value()==0 && ydim.value()==0)
+ {
+ volume4D<float> RawData;
+ volume<float> theMask;
+ toc();
+ //read data
+ message("Reading data file " << (string)fnin.value() << " ... ");
+ read_volume4D(RawData,fnin.value());
+ message(" done" << endl);
- Mean = meanvol(RawData);
- toc();
- message("Reading mask file " << (string)fnmask.value() << " ... ");
- read_volume(theMask,fnmask.value());
+ Mean = meanvol(RawData);
+ toc();
+ message("Reading mask file " << (string)fnmask.value() << " ... ");
+ read_volume(theMask,fnmask.value());
- memmsg(" Before reshape ");
- MatrixData = RawData.matrix(theMask);
+ memmsg(" Before reshape ");
+ MatrixData = RawData.matrix(theMask);
}
- memmsg(" after reshape ");
- message(" Data size " << MatrixData.Nrows() << " x " << MatrixData.Ncols() << endl);
- toc();
- memmsg(" before remmean_econ ");
- remmean_econ(MatrixData);
- memmsg("after remmean_econ / before remmean")
- toc();
+ else{
+ Matrix data = unifrnd(xdim.value(),ydim.value());
+ outMsize("data", data);
+ tic(); calccorr(data,econ.value()); toc();
- MatrixData = remmean(MatrixData);
- memmsg(" after remmean ");
- toc();
- message(" Matrix size " << MatrixData.Nrows() << " x " << MatrixData.Ncols() << endl);
-
- memmsg(" before calc_corr ");
- toc();
- calccorr(MatrixData,FALSE);
- toc();
- memmsg(" after calc_corr ");
-
- memmsg(" before calc_corr 2");
- toc();
- calccorr(MatrixData,TRUE);
- toc();
- memmsg(" after calc_corr 2");
-
+ data = unifrnd(10,100);
+ outMsize("data", data);
+ tic(); calccorr(data,econ.value()); toc();
+
+ data = unifrnd(100,1000);
+ outMsize("data", data);
+ tic(); calccorr(data,econ.value()); toc();
+
+ data = unifrnd(100,10000);
+ outMsize("data", data);
+ tic(); calccorr(data,econ.value()); toc();
+
+ data = unifrnd(300,200000);
+ outMsize("data", data);
+ tic(); calccorr(data,econ.value()); toc();
+
+ data = unifrnd(500,20000);
+ outMsize("data", data);
+ tic(); calccorr(data,econ.value()); toc();
+
+ data = unifrnd(500,200000);
+ outMsize("data", data);
+ tic(); calccorr(data,econ.value()); toc();
+
+
+ }
return 0;
@@ -178,6 +210,7 @@ int do_work(int argc, char* argv[]) {
options.add(help);
options.add(xdim);
options.add(ydim);
+ options.add(econ);
options.parse_command_line(argc, argv);
// line below stops the program if the help was requested or