diff --git a/Makefile b/Makefile
index 744fb651452b52d246487b1c5311cb894aaac5d9..940ec2d3af7611a8a6685e100c8914e1c33ca7e4 100644
--- a/Makefile
+++ b/Makefile
@@ -16,6 +16,8 @@ TEST_OBJS = test.o
GGMIX_OBJS = ggmix.o
+FSL_DR_OBJS = melhlprfns.o fsl_dualreg.o
+
FSL_GLM_OBJS = melhlprfns.o fsl_glm.o
FSL_SBCA_OBJS = melhlprfns.o fsl_sbca.o
@@ -32,7 +34,7 @@ MELODIC_OBJS = meloptions.o melhlprfns.o melgmix.o meldata.o melpca.o melica.o
TESTXFILES = test
-XFILES = fsl_glm fsl_sbca fsl_mvlm fsl_regfilt fsl_schurprod fsl_bwtf melodic
+XFILES = fsl_dualreg fsl_glm fsl_sbca fsl_mvlm fsl_regfilt fsl_schurprod fsl_bwtf melodic
OTHERFILES =
@@ -63,7 +65,7 @@ RUNTCLS = Melodic
SCRIPTS = melodicreport dual_regression
-all: ggmix fsl_regfilt fsl_glm melodic ${OTHERFILES}
+all: ggmix fsl_regfilt fsl_glm fsl_dr melodic ${OTHERFILES}
ggmix: ${GGMIX_OBJS}
${AR} -r libggmix.a ${GGMIX_OBJS}
@@ -74,6 +76,9 @@ melodic: ${MELODIC_OBJS}
test: ${TEST_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${TEST_OBJS} ${LIBS}
+fsl_dr: ${FSL_DR_OBJS}
+ $(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_DR_OBJS} ${LIBS}
+
fsl_glm: ${FSL_GLM_OBJS}
$(CXX) ${CXXFLAGS} ${LDFLAGS} -o $@ ${FSL_GLM_OBJS} ${LIBS}
diff --git a/meldata.cc b/meldata.cc
index fe93c7fa12a94b06112982f7b21bbf9092783ccb..1eba93ed671ce380b969247ad223101212b0f4b8 100644
--- a/meldata.cc
+++ b/meldata.cc
@@ -79,7 +79,7 @@ namespace Melodic{
if(opts.remove_meantc.value()){
remmean(tmpData,meanC,2);
}
-
+
//convert to power spectra
if(opts.pspec.value()){
message(" Converting data to powerspectra ...");
@@ -115,18 +115,42 @@ namespace Melodic{
message(" done" << endl);
}
- //convert to instacorrs
- if(opts.insta_fn.value()>""){
- Matrix vscales = pow(stdev(tmpData,1),-1);
- varnorm(tmpData,vscales);
-
- Matrix tmpTC = tmpData * insta_mask.t();
- varnorm(tmpTC,pow(stdev(tmpTC),-1));
-
- for(int ctr=1; ctr <=tmpData.Ncols();ctr++)
- tmpData.Column(ctr) = SP(tmpData.Column(ctr),tmpTC);
+ //INSTACORRS
+ Matrix tmpTC;
+ tmpTC = tmpData * insta_maps.t();
+
+ if(opts.insta_fn.value().length()>0){
+
+ if(opts.insta_idx.value()<1 || opts.insta_idx.value()>tmpTC.Ncols()){
+ cerr << "ERROR:: INSTACORR index is wrong \n\n";
+ exit(2);
+ }
+ Matrix tmpRef = tmpTC.Column(opts.insta_idx.value());
+ if(opts.insta_idx.value()>1){
+ // swap columns
+ tmpTC.Column(opts.insta_idx.value()) << tmpTC.Column(1);
+ tmpTC.Column(1) << tmpRef;
+ }
+
+ if(opts.insta_partial.value() && tmpTC.Ncols()>1){
+ // partal correlations
+ Matrix tmpConf = tmpTC.Columns(2,tmpTC.Ncols());
+ tmpData -= tmpConf * (pinv(tmpConf) * tmpData);
+ tmpRef -= tmpConf * (pinv(tmpConf) * tmpRef);
+ }
+
+ if(opts.insta_varnorm.value()){
+ Matrix vscales = pow(stdev(tmpData,1),-1);
+ varnorm(tmpData,vscales);
+ varnorm(tmpRef,pow(stdev(tmpRef,1),-1));
+ }
+
+ // Shur product
+ SP4(tmpData,tmpRef);
}
+ //END INSTACORRS
+
tmpData.Release();
dbgmsg(string("END: process_file") << endl);
@@ -620,19 +644,6 @@ namespace Melodic{
create_RXweight();
}
- //set up instacorr mask image
- if(opts.insta_fn.value()>""){
- dbgmsg(string(" Setting up instacorr mask") << endl);
- volume4D<float> tmp_im;
- read_volume4D(tmp_im,opts.insta_fn.value());
-
- if(!samesize(Mean,tmp_im[0])){
- cerr << "ERROR:: instacorr mask and data have different voxel dimensions \n\n";
- exit(2);
- }
- insta_mask = tmp_im.matrix(Mask);
- }
-
//seed the random number generator
double tmptime = time(NULL);
if ( opts.seed.value() != -1 ) {
@@ -680,6 +691,20 @@ namespace Melodic{
}
remmean(Tdes);
+ //INSTACORRS
+ if(opts.insta_fn.value().length()>0){
+ message(" Reading in " << opts.insta_fn.value()
+ << " for instantaneous correlation analysis" <<endl);
+ volume4D<float> tmp_im;
+ read_volume4D(tmp_im,opts.insta_fn.value());
+
+ if(!samesize(Mean,tmp_im[0])){
+ cerr << "ERROR:: instacorr mask and data have different voxel dimensions \n\n";
+ exit(2);
+ }
+ insta_maps = tmp_im.matrix(Mask);
+ }
+
dbgmsg(string("END: setup_misc") << endl);
}
diff --git a/meldata.h b/meldata.h
index f1ba9d23b04e7f82e44801742a2b72cd7211c342..9a04ac24be3b7e2bef58a7c476d947def7e62a67 100644
--- a/meldata.h
+++ b/meldata.h
@@ -241,7 +241,7 @@ namespace Melodic{
volume<float> Mask;
volume<float> Mean;
volume<float> background;
- Matrix insta_mask;
+ Matrix insta_maps;
Matrix Data;
Matrix PPCA;
diff --git a/melhlprfns.cc b/melhlprfns.cc
index 5f7965d200e66a30af022a8a7e764a202acf7873..c662e17f27525f3ba12c16266f3c4e291ec9e6ee 100644
--- a/melhlprfns.cc
+++ b/melhlprfns.cc
@@ -171,12 +171,14 @@ namespace Melodic{
Matrix Res;
Res = in;
if(econ>0){
- ColumnVector tmp;
- for(int ctr=1; ctr <= in.Ncols(); ctr++){
- tmp = in.Column(ctr);
- tmp = tmp * weights(1,ctr);
- Res.Column(ctr) = tmp;
- }
+ if(weights.Ncols() == in.Ncols()){
+ ColumnVector tmp;
+ for(int ctr=1; ctr <= in.Ncols(); ctr++){
+ tmp = in.Column(ctr);
+ tmp = tmp * weights(1,ctr);
+ Res.Column(ctr) = tmp;
+ }
+ }
}
else{
Res = ones(in.Nrows(),1)*weights.Row(1);
@@ -192,6 +194,15 @@ namespace Melodic{
}
}
+ void SP4(Matrix& in, const Matrix& weights){
+ RowVector tmp;
+ for(int ctr=1; ctr <= in.Nrows(); ctr++){
+ tmp = in.Row(ctr);
+ tmp = tmp * weights(ctr,1);
+ in.Row(ctr) << tmp;
+ }
+ }
+
Matrix corrcoef(const Matrix& in1, const Matrix& in2){
Matrix tmp = in1;
tmp |= in2;
diff --git a/melhlprfns.h b/melhlprfns.h
index 0f06f333241ba65d133e869aa5e172dcbd6504f6..1ff66df9b639d2c80db736b777c8309f087d8e06 100644
--- a/melhlprfns.h
+++ b/melhlprfns.h
@@ -50,6 +50,7 @@ namespace Melodic{
Matrix SP2(const Matrix& in, const Matrix& weights, int econ = 20000);
void SP3(Matrix& in, const Matrix& weights);
+ void SP4(Matrix& in, const Matrix& weights);
RowVector Feta(int n1,int n2);
RowVector cumsum(const RowVector& Inp);
diff --git a/meloptions.cc b/meloptions.cc
index e7dc519aa158f9e293f1074fe30589161b07e1ae..24c9af7fd5325f55c204cb377b11cb3313c0718a 100644
--- a/meloptions.cc
+++ b/meloptions.cc
@@ -157,6 +157,10 @@ MelodicOptions* MelodicOptions::gopt = NULL;
}
if (insta_fn.value() > ""){
varnorm.set_T(false);
+ cout << " inc_fn = " << insta_fn.value() << endl;
+ cout << " inc_idx = " << insta_idx.value() << endl;
+ cout << " inc_partial = " << insta_partial.value() << endl;
+ cout << " inc_varnorm = " << insta_varnorm.value() << endl;
}
//in the case of indirect inputs, create the vector of input names here
diff --git a/meloptions.h b/meloptions.h
index e2e9be6b9c3feb001eff19689bd25b97fb6e051b..020884ff7e5e23fcfb264d37b188d321c6644ab9 100644
--- a/meloptions.h
+++ b/meloptions.h
@@ -125,7 +125,11 @@ class MelodicOptions {
Option<float> nlconst1;
Option<float> nlconst2;
Option<float> smooth_probmap;
+
Option<string> insta_fn;
+ Option<int> insta_idx;
+ Option<bool> insta_partial;
+ Option<bool> insta_varnorm;
Option<bool> remove_meanvol;
Option<bool> remove_meantc;
@@ -392,9 +396,18 @@ class MelodicOptions {
smooth_probmap(string("--smooth_pm"), 0.0,
string("width of smoothing kernel for probability maps"),
false, requires_argument, false),
- insta_fn(string("--insta_fn"), string(""),
- string(" mask file name for instacorr calculation"),
+ insta_fn(string("--inc_fn"), string(""),
+ string(" file name for instacorr volumes"),
false, requires_argument, false),
+ insta_idx(string("--inc_idx"), 1,
+ string(" index for instacorr calculation"),
+ false, requires_argument, false),
+ insta_partial(string("--inc_partial"), true,
+ string(" switch off partial correlation analysis for instacorrs"),
+ false, requires_argument, false),
+ insta_varnorm(string("--inc_varnorm"), true,
+ string(" switch off varnorm for instantaneous correlations"),
+ false, no_argument, false),
remove_meanvol(string("--keep_meanvol"), true,
string("do not subtract mean volume"),
false, no_argument, false),
@@ -501,6 +514,9 @@ class MelodicOptions {
options.add(nlconst2);
options.add(smooth_probmap);
options.add(insta_fn);
+ options.add(insta_idx);
+ options.add(insta_partial);
+ options.add(insta_varnorm);
options.add(remove_meanvol);
options.add(remove_meantc);
options.add(remove_endslices);
diff --git a/test.cc b/test.cc
index 648574af43909f9407aa491fa2cb2483e17e7a7e..bc0fb918f2217522e8e5f8d25f6772b86ec58428 100644
--- a/test.cc
+++ b/test.cc
@@ -136,7 +136,18 @@ Matrix calccorr(const Matrix& in, int econ)
return Res;
} //Matrix calccorr
+void convert_to_mat(volume4D<float>& in, volume<float>& mask, Matrix& out)
+{
+ out = in[0].vec(mask).t();
+ in.deletevolume(0);
+ while(in.tsize()>0){
+ out &= in[0].vec(mask).t();
+ in.deletevolume(0);
+ }
+}
+
int do_work(int argc, char* argv[]) {
+
tic();
Matrix MatrixData;
@@ -151,14 +162,24 @@ int do_work(int argc, char* argv[]) {
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());
+ MatrixData = RawData.matrix();
+
+ memmsg("start");
+ Matrix res;
+ convert_to_mat(RawData,theMask, res);
+
memmsg(" Before reshape ");
- MatrixData = RawData.matrix(theMask);
+
+ write_ascii_matrix("mat1", res);
+ write_ascii_matrix("mat2", MatrixData);
+
+
}
else{
Matrix data = unifrnd(xdim.value(),ydim.value());