diff --git a/meldata.cc b/meldata.cc
index eb72dfbe8c4cd57fea5e9824ff2a168735fc4b86..bc4f70a5f3e381155e384390c43b46dade6a7d4c 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -24,25 +24,41 @@ using namespace NEWIMAGE;
namespace Melodic{
// {{{ Setup
- Matrix MelodicData::process_file(string fname, int numfiles)
+ ReturnMatrix MelodicData::process_file(string fname, int numfiles)
{
- volume4D<float> RawData;
+ Matrix tmpData;
+ {
+ volume4D<float> RawData;
- //read data
- message("Reading data file " << fname << " ... ");
- read_volume4D(RawData,fname);
- message(" done" << endl);
- del_vols(RawData,opts.dummy.value());
+ if(opts.debug.value())
+ memmsg(" before reading file "<< fname);
+
+ //read data
+ message("Reading data file " << fname << " ... ");
+ read_volume4D(RawData,fname);
+ message(" done" << endl);
+ if(opts.debug.value())
+ memmsg(" after reading file "<< fname);
+
+ del_vols(RawData,opts.dummy.value());
- Mean += meanvol(RawData)/numfiles;
+ Mean += meanvol(RawData)/numfiles;
+
+ //estimate smoothness
+ if(opts.debug.value())
+ memmsg(" before est smoothness ");
+ if((Resels == 0)&&(!opts.filtermode))
+ Resels = est_resels(RawData,Mask);
+ if(opts.debug.value())
+ memmsg(" after smoothness ");
- //reshape
- Matrix tmpData;
- tmpData = RawData.matrix(Mask);
-
- //estimate smoothness
- if((Resels == 0)&&(!opts.filtermode))
- Resels = est_resels(RawData,Mask);
+ //reshape
+ if(opts.debug.value())
+ memmsg(" before reshape ");
+ tmpData = RawData.matrix(Mask);
+ if(opts.debug.value())
+ memmsg(" after reshape ");
+ }
//convert to percent BOLD signal change
if(opts.pbsc.value()){
@@ -55,16 +71,18 @@ namespace Melodic{
if(opts.remove_meanvol.value())
{
message(string(" Removing mean image ..."));
- meanR = mean(tmpData);
- tmpData = remmean(tmpData);
+ if(opts.debug.value())
+ memmsg(" before remmean ");
+ remmean(tmpData,meanR,1);
+ if(opts.debug.value())
+ memmsg(" after remmean ");
message(" done" << endl);
}
else meanR=ones(1,tmpData.Ncols());
}
if(opts.remove_meantc.value()){
- meanC = mean(tmpData,2);
- tmpData = remmean(tmpData,2);
+ remmean(tmpData,meanC,2);
}
//convert to power spectra
@@ -86,6 +104,8 @@ namespace Melodic{
}
//variance - normalisation
+ if(opts.debug.value())
+ memmsg(" before VN ");
if(opts.varnorm.value()){
message(" Normalising by voxel-wise variance ...");
if(stdDev.Storage()==0)
@@ -95,9 +115,12 @@ namespace Melodic{
stdDev += varnorm(tmpData,std::min(30,tmpData.Nrows()-1),
opts.vn_level.value())/numfiles;
stdDevi = pow(stdDev,-1);
+ if(opts.debug.value())
+ memmsg(" in VN ");
message(" done" << endl);
}
+ tmpData.Release();
return tmpData;
}
@@ -152,7 +175,7 @@ namespace Melodic{
dbgmsg(string(" approach ") << opts.approach.value() << endl);
- if(opts.approach.value()!=string("concat")){
+ if(opts.approach.value()==string("tica")){
message("Calculating T- and S-modes " << endl);
explained_var = krfact(tmp,tmpT,tmpS);
Tmodes.clear(); Smodes.clear();
@@ -180,6 +203,7 @@ namespace Melodic{
add_Tmodes(tmpT3);
}
}
+
if(opts.approach.value()!=string("concat")){
//add GLM OLS fit
dbgmsg(string(" GLM fitting ") << endl);
@@ -213,7 +237,9 @@ namespace Melodic{
{
numfiles = (int)opts.inputfname.value().size();
setup_misc();
-
+ if(opts.debug.value())
+ memmsg(" after setup_misc ");
+
if(opts.filtermode){ // basic setup for filtering only
Data = process_file(opts.inputfname.value().at(0));
}
@@ -229,7 +255,9 @@ namespace Melodic{
opts.varnorm.set_T(false);
}
alldat = process_file(opts.inputfname.value().at(0), numfiles) / numfiles;
-
+ if(opts.debug.value())
+ memmsg(" after process_file ");
+
if(opts.pca_dim.value() > alldat.Nrows()-2){
cerr << "ERROR:: too many components selected \n\n";
exit(2);
@@ -285,11 +313,11 @@ namespace Melodic{
RowVector AdjEV, PercEV;
Matrix Corr, tmpE;
int order;
- // cerr << "here1" << endl;
+
order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
if (opts.paradigmfname.value().length()>0)
order += param.Ncols();
- // cerr << "here2" << endl;
+
if(opts.pca_dim.value() == 0){
opts.pca_dim.set_T(order);
PPCA=tmpPPCA;
@@ -332,7 +360,6 @@ namespace Melodic{
WM.push_back(tmp);
}
else {
- //cerr << "here" << endl;
for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
@@ -468,7 +495,7 @@ namespace Melodic{
Scon &= -1*Scon;
}
}
- Tdes = remmean(Tdes,1);
+ remmean(Tdes);
}
void MelodicData::save()
diff --git a/meldata.h b/meldata.h
index 902183482fc464742a8de7c1666456e69442a808..12629db521e52c8a5b16752e3c3631445c25b634 100644
--- a/meldata.h
+++ b/meldata.h
@@ -36,7 +36,7 @@ namespace Melodic{
void save();
- Matrix process_file(string fname, int numfiles = 1);
+ ReturnMatrix process_file(string fname, int numfiles = 1);
inline void save4D(Matrix what, string fname){
volume4D<float> tempVol;
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 9035a8f762041784106671af05a1b90976880410..0758bc36316406bdd4e1cc5ff46b7664dc4d58f6 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -22,7 +22,6 @@ namespace Melodic{
Matrix DStDev=stdev(Data);
volume4D<float> tmpMask, RawData;
tmpMask.setmatrix(DStDev,mask);
- RawData.setmatrix(Data,mask);
float tMmax;
volume<float> tmpMask2;
@@ -32,9 +31,23 @@ namespace Melodic{
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);
+ volume4D<float> newmask;
+ newmask = binarise(tmpMask2,(float) max((float) st_mean-3*st_std,(float) 0.01*st_mean),tMmax);
- Data = RawData.matrix(mask);
+ Matrix newmaskM,newData;
+ newmaskM = newmask.matrix(mask);
+ int N = Data.Nrows();
+
+ if(int(newmaskM.Row(1).SumAbsoluteValue() + 0.3) < newmaskM.Ncols()){
+ RawData.setmatrix(Data.Row(1),mask);
+ newData = RawData.matrix(newmask[0]);
+ for(int r=2; r <= N; r++){
+ RawData.setmatrix(Data.Row(r),mask);
+ newData &= RawData.matrix(newmask[0]);
+ }
+ Data = newData;
+ mask = newmask[0];
+ }
}
void del_vols(volume4D<float>& in, int howmany)
@@ -123,12 +136,14 @@ namespace Melodic{
void varnorm(Matrix& in, const RowVector& vars)
{
- Matrix tmp = vars;
- in = SP(in,pow(ones(in.Nrows(),1) * tmp,-1));
+ for(int ctr=1; ctr <=in.Nrows();ctr++){
+ in.Row(ctr) = SD(in.Row(ctr),vars);
+ }
}
RowVector varnorm(Matrix& in, Matrix& Corr, int dim, float level)
{
+
Matrix tmpE, white, dewhite;
RowVector tmpD, tmpD2;
@@ -147,6 +162,7 @@ namespace Melodic{
in.Column(ctr1) = 0.0*in.Column(ctr1);
}
varnorm(in,tmpD);
+
return tmpD;
} //RowVector varnorm
diff --git a/melhlprfns.h b/melhlprfns.h
index 516445193d0be4d93c541ca1ec65c73378fe453d..8e16fd1b6963280fd332443f91265932a3cfecef 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -16,6 +16,22 @@
#include "newmatap.h"
#include "newmatio.h"
+ #ifdef __APPLE__
+ #include <mach/mach.h>
+ #define mmsg(msg) { \
+ struct task_basic_info t_info; \
+ mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
+ if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
+ { \
+ cout << " MEM: " << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
+ } \
+ }
+ #else
+ #define mmsg(msg) { \
+ cout << msg; \
+ }
+ #endif
+
using namespace NEWIMAGE;
namespace Melodic{
@@ -39,8 +55,8 @@ namespace Melodic{
Matrix corrcoef(const Matrix& in1, const Matrix& in2);
Matrix corrcoef(const Matrix& in1, const Matrix& in2, const Matrix& part);
- Matrix calc_corr(const Matrix& in, bool econ = 0);
- Matrix calc_corr(const Matrix& in, const Matrix& weights, bool econ = 0);
+ Matrix calc_corr(const Matrix& in, bool econ = 1);
+ Matrix calc_corr(const Matrix& in, const Matrix& weights, bool econ = 1);
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& param, Matrix& paramS, Matrix& white, Matrix& dewhite);
float calc_white(const Matrix& tmpE, const RowVector& tmpD, const RowVector& PercEV, int dim, Matrix& white, Matrix& dewhite);
diff --git a/melodic.cc b/melodic.cc
index b78b5ff644485c5b39ee7d22ad68292c7656c55a..d111afed03169f277c33e04a6966c192549d0301 100644
--- a/melodic.cc
+++ b/melodic.cc
@@ -311,7 +311,7 @@ Matrix mmall(Log& logger, MelodicOptions& opts,MelodicData& melodat, MelodicRepo
if(opts.varnorm.value()&&melodat.get_stdNoisei().Storage()>0){
ICmap = SP(melodat.get_IC().Row(ctr),diagvals(ctr)*melodat.get_stdNoisei());
}else{
- ICmap = melodat.get_IC().Row(ctr);
+ ICmap = melodat.get_IC().Row(ctr);
}
string wherelog;
if(opts.genreport.value())
diff --git a/melodic.h b/melodic.h
index 1f6c80a2ba804fa55013e5556fc059b4b95449ca..126042356b9ca829d383a475ef31df4167bc9acc 100644
--- a/melodic.h
+++ b/melodic.h
@@ -14,6 +14,25 @@
#include <strstream>
+#ifdef __APPLE__
+#include <mach/mach.h>
+#define memmsg(msg) { \
+ struct task_basic_info t_info; \
+ mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
+ if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
+ { \
+ cout << " MEM: " << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
+ cout.flush(); \
+ } \
+}
+#else
+#define memmsg(msg) { \
+ cout << msg; \
+}
+#endif
+
+
+
// a simple message macro that takes care of cout and log
#define message(msg) { \
MelodicOptions& opt = MelodicOptions::getInstance(); \
@@ -40,7 +59,7 @@
namespace Melodic{
-const string version = "3.10";
+const string version = "3.12";
// The two strings below specify the title and example usage that is
// printed out as the help or usage message
diff --git a/melreport.cc b/melreport.cc
index a1bc2f640c05fcb56e277a05b11231c2d3d73ae6..76fe5c610b787dd827f7a82f20201de15a6170ee 100644
--- a/melreport.cc
+++ b/melreport.cc
@@ -107,9 +107,11 @@ namespace Melodic{
}
{//plot time course
- IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
- miscplot newplot;
+ IChtml << "<H3> Temporal mode </H3><p>" << endl <<endl;
+ miscplot newplot;
+
Matrix tmptc = melodat.get_Tmodes(cnum-1).Column(1).t();
+
newplot.col_replace(0,0xFF0000);
newplot.add_label(string("IC ")+num2str(cnum)+" time course");
@@ -141,7 +143,8 @@ namespace Melodic{
tmptc.Row(1).Minimum()),tmptc.Row(1).Maximum()+
0.05*(tmptc.Row(1).Maximum()-tmptc.Row(1).Minimum()));
newplot.grid_swapdefault();
- newplot.timeseries(tmptc,
+
+ newplot.timeseries(tmptc,
report.appendDir(string("t")+num2str(cnum)+".png"),
string("Timecourse No. ")+num2str(cnum),
opts.tr.value(),150,12,1,false);
diff --git a/test.cc b/test.cc
index 3a13d0e9c3754428e942a1851e3af31a3fcedc45..29e422c532d6cb8adc5335d43a7593454792bb3b 100644
--- a/test.cc
+++ b/test.cc
@@ -13,6 +13,31 @@
#include "newimage/newimageall.h"
#include "melhlprfns.h"
+#ifdef __APPLE__
+#include <mach/mach.h>
+#define memmsg(msg) { \
+ struct task_basic_info t_info; \
+ mach_msg_type_number_t t_info_count = TASK_BASIC_INFO_COUNT; \
+ if (KERN_SUCCESS == task_info(mach_task_self(), TASK_BASIC_INFO, (task_info_t) &t_info, &t_info_count)) \
+ { \
+ cout << msg << " res: " << t_info.resident_size/1000000 << " virt: " << t_info.virtual_size/1000000 << "\n"; \
+ cout.flush(); \
+ } \
+}
+#else
+#define memmsg(msg) { \
+ cout << msg; \
+}
+#endif
+
+
+// a simple message macro that takes care of cout
+#define message(msg) { \
+ cout << msg; \
+ cout.flush(); \
+}
+
+
using namespace MISCPLOT;
using namespace MISCMATHS;
using namespace Utilities;
@@ -49,11 +74,36 @@ using namespace std;
int do_work(int argc, char* argv[]) {
- Matrix test;
+ Matrix MatrixData;
+ volume<float> Mean;
+
+ {
+ volume4D<float> RawData;
+
+ //read data
+ message("Reading data file " << (string)fnin.value() << " ... ");
+ read_volume4D(RawData,fnin.value());
+ message(" done" << endl);
+
+ Mean = meanvol(RawData);
+
+ memmsg(" Before reshape ");
+ MatrixData = RawData.matrix();
+ }
+ memmsg(" after reshape ");
+ message(" Data size " << MatrixData.Nrows() << " x " << MatrixData.Ncols() << endl);
+
+ memmsg(" before remmean_econ ");
+ remmean_econ(MatrixData);
+ memmsg("after remmean_econ / before remmean")
+
+ MatrixData = remmean(MatrixData);
+
+ memmsg(" after remmean ");
+ message(" Mean size " << MatrixData.Nrows() << " x " << MatrixData.Ncols() << endl);
+
+ message(" Mean size " << MatrixData.Nrows() << " x " << MatrixData.Ncols() << endl);
- cerr << " X: "<< xdim.value() << ", Y: "<< ydim.value() << endl;
- test = zeros(xdim.value(),ydim.value());
- cerr << "Created matrix of size " << test.Nrows() << " x " << test.Ncols() << endl;
return 0;
}