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

MIGP v 0.1

parent 5e1f5f83
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
meldata.cc
+ 104
− 67
View file @ 33a7ba6d
......@@ -15,6 +15,7 @@
#include "melodic.h"
#include "utils/log.h"
#include <time.h>
#include <algorithm>
#include "miscmaths/miscprob.h"
#include "melhlprfns.h"
......@@ -30,34 +31,28 @@ namespace Melodic{
{
volume4D<float> RawData;
if(opts.debug.value())
memmsg(" before reading file "<< fname);
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);
memmsg(" after reading file "<< fname);
del_vols(RawData,opts.dummy.value());
Mean += meanvol(RawData)/numfiles;
//estimate smoothness
if(opts.debug.value())
memmsg(" before est smoothness ");
memmsg(" before est smoothness ");
if((Resels == 0)&&(!opts.filtermode))
Resels = est_resels(RawData,Mask);
if(opts.debug.value())
memmsg(" after smoothness ");
memmsg(" after smoothness ");
//reshape
if(opts.debug.value())
memmsg(" before reshape ");
memmsg(" before reshape ");
tmpData = RawData.matrix(Mask);
if(opts.debug.value())
memmsg(" after reshape ");
memmsg(" after reshape ");
}
//convert to percent BOLD signal change
......@@ -71,11 +66,9 @@ namespace Melodic{
if(opts.remove_meanvol.value())
{
message(string(" Removing mean image ..."));
if(opts.debug.value())
memmsg(" before remmean ");
memmsg(" before remmean ");
remmean(tmpData,meanR,1);
if(opts.debug.value())
memmsg(" after remmean ");
memmsg(" after remmean ");
message(" done" << endl);
}
else meanR=ones(1,tmpData.Ncols());
......@@ -104,8 +97,7 @@ namespace Melodic{
}
//variance - normalisation
if(opts.debug.value())
memmsg(" before VN ");
memmsg(" before VN ");
if(opts.varnorm.value()){
message(" Normalising by voxel-wise variance ...");
if(stdDev.Storage()==0)
......@@ -115,8 +107,7 @@ 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 ");
memmsg(" in VN ");
message(" done" << endl);
}
......@@ -233,20 +224,8 @@ namespace Melodic{
dbgmsg(string("END: set_TSmode"));
}
void MelodicData::setup()
{
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));
}
else{
if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
opts.approach.set_T("tica");
void MelodicData::setup_classic()
{
Matrix alldat, tmpData;
bool tmpvarnorm = opts.varnorm.value();
......@@ -255,16 +234,16 @@ 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 ");
memmsg(" after process_file ");
if(opts.pca_dim.value() > alldat.Nrows()-2){
cerr << "ERROR:: too many components selected \n\n";
exit(2);
}
if(opts.debug.value())
save4D(alldat,string("preproc_dat") + num2str(1));
for(int ctr = 1; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles) / numfiles;
if(opts.debug.value())
......@@ -276,7 +255,7 @@ namespace Melodic{
cerr << "ERROR:: data dimensions do not match, TICA not possible \n\n";
exit(2);
}
if(tmpData.Ncols() == alldat.Ncols()){
int mindim = min(alldat.Nrows(),tmpData.Nrows());
alldat = alldat.Rows(1,mindim);
......@@ -285,7 +264,7 @@ namespace Melodic{
}
else
message("Data dimensions do not match - ignoring "+opts.inputfname.value().at(ctr) << endl);
}
}
}
//update mask
......@@ -302,7 +281,7 @@ namespace Melodic{
stdDevi = pow(stdDev,-1);
message(" done" << endl);
}
if(numfiles>1)
message(endl << "Initial data size : "<<alldat.Nrows()<<" x "<<alldat.Ncols()<<endl<<endl);
......@@ -333,7 +312,7 @@ namespace Melodic{
Matrix tmpPscales;
tmpPscales = param.t() * alldat;
paramS = stdev(tmpPscales.t());
calc_white(pcaE, pcaD, order, param, paramS, whiteMatrix, dewhiteMatrix);
}else
calc_white(pcaE, pcaD, order, whiteMatrix, dewhiteMatrix);
......@@ -358,11 +337,11 @@ namespace Melodic{
WM.push_back(tmp);
}
else {
dbgmsg("Multi-Subject ICA");
for(int ctr = 0; ctr < numfiles; ctr++){
dbgmsg("Multi-Subject ICA");
for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles);
if(opts.joined_vn.value() && tmpvarnorm){
dbgmsg("tmpData normalisation"<< endl);
dbgmsg("stdDev " << stdDev(1,2)<< endl);
......@@ -399,29 +378,95 @@ namespace Melodic{
DWM.push_back(newDWM);
WM.push_back(newWM);
tmpData = newWM * tmpData;
//concatenate Data
if(Data.Storage() == 0)
Data = tmpData;
else
Data &= tmpData;
}
}
}
opts.varnorm.set_T(tmpvarnorm);
opts.varnorm.set_T(tmpvarnorm);
}
message(endl << " Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
outMsize("stdDev",stdDev);
//meanC=mean(Data,2);
void MelodicData::setup_migp()
{
dbgmsg("starting MIGP");
//permute input vector if desired std::random_shuffle ( opts.inputfname.value().begin(), opts.inputfname.value().end() );
Matrix tmpData;
bool tmpvarnorm = opts.varnorm.value();
if(numfiles > 1 && opts.joined_vn.value()){
opts.varnorm.set_T(false);
}
for(int ctr = 1; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles) / numfiles;
if(opts.debug.value())
save4D(Data,"concat_data");
//save the mean & mask
save_volume(Mask,logger.appendDir("mask"));
save_volume(Mean,logger.appendDir("mean"));
save4D(tmpData,string("preproc_dat") + num2str(ctr+1));
if(Data.Storage()==0)
Data = tmpData;
else
Data &= tmpData;
//reduce dim down to manageable level
if(Data.Nrows() > opts.migpN.value()){
message(" Reducing data matrix to a " << opt.migpN.value() << " dimensional subspace " << endl);
Matrix pcaE, Corr;
RowVector pcaD;
std_pca(Data, RXweight, Corr, pcaE, pcaD);
pcaE = pcaE.Columns(pcaE.Ncols()-opts.migpN.value()+1,pcaE.Ncols());
Data = pcaE.t() * Data;
}
}
//update mask
if(opts.update_mask.value()){
message("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);
}
void MelodicData::setup()
{
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));
}
else{
if((numfiles > 1) && (opts.approach.value()==string("defl") || opts.approach.value()==string("symm")))
opts.approach.set_T("concat");
if(opts.migp.value())
setup_migp();
else
setup_classic();
}
message(endl << " Data size : "<<Data.Nrows()<<" x "<<Data.Ncols()<<endl<<endl);
outMsize("stdDev",stdDev);
//meanC=mean(Data,2);
if(opts.debug.value())
save4D(Data,"concat_data");
//save the mean & mask
save_volume(Mask,logger.appendDir("mask"));
save_volume(Mean,logger.appendDir("mean"));
} // void setup()
void MelodicData::setup_misc()
{
......@@ -503,7 +548,6 @@ namespace Melodic{
remmean(Tdes);
}
void MelodicData::save()
{
......@@ -560,7 +604,7 @@ namespace Melodic{
}
if(opts.output_origIC.value())
save4D(stdNoisei,string("Noise_stddev_inv"));
save4D(stdNoisei,string("Noise__inv"));
}
//Output T- & S-modes
......@@ -623,7 +667,6 @@ namespace Melodic{
message("...done" << endl);
} //void save()
int MelodicData::remove_components()
{
message("Reading " << opts.filtermix.value() << endl)
......@@ -698,7 +741,6 @@ namespace Melodic{
return 0;
} // int remove_components()
void MelodicData::create_RXweight()
{
message("Reading the weights for the covariance R_X from file "<< opts.segment.value() << endl);
......@@ -708,7 +750,6 @@ namespace Melodic{
RXweight = tmpRX.matrix(Mask);
}
void MelodicData::est_smoothness()
{
if(Resels == 0){
......@@ -745,7 +786,6 @@ namespace Melodic{
}
}
unsigned long MelodicData::standardise(volume<float>& mask, volume4D<float>& R)
{
unsigned long count = 0;
......@@ -787,7 +827,6 @@ namespace Melodic{
return count;
}
float MelodicData::est_resels(volume4D<float> R, volume<float> mask)
{
message(" Estimating data smoothness ... ");
......@@ -863,7 +902,6 @@ namespace Melodic{
return resels;
}
void MelodicData::create_mask(volume<float>& theMask)
{
if(opts.use_mask.value() && opts.maskfname.value().size()>0){ // mask provided
......@@ -938,7 +976,6 @@ namespace Melodic{
}
} //void create_mask()
void MelodicData::sort()
{
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
......
meldata.h
+ 4
− 0
View file @ 33a7ba6d
......@@ -56,7 +56,11 @@ namespace Melodic{
}
int remove_components();
void setup_classic();
void setup_migp();
void setup();
void status(const string &txt);
inline Matrix& get_pcaE() {return pcaE;}
......
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