Skip to content
Snippets Groups Projects
Commit bc86c836 authored by Christian Beckmann's avatar Christian Beckmann
Browse files

ready for 5_0_4

parent e2a18321
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
meldata.cc
+ 10
− 3
View file @ bc86c836
......@@ -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()
......
meloptions.h
+ 6
− 1
View file @ bc86c836
......@@ -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);
......
test.cc
+ 101
− 68
View file @ bc86c836
......@@ -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(&current_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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment