Skip to content
Snippets Groups Projects

Rf meanvol

Merged Rf meanvol
rf_meanvol into master
Merged Matthew Webster requested to merge rf_meanvol into master
Compare
and
2 files
+ 192
− 196
Compare changes
  • Side-by-side
  • Inline
Files
2
meldata.cc
+ 191
− 195
/* MELODIC - Multivariate exploratory linear optimized decomposition into
/* MELODIC - Multivariate exploratory linear optimized decomposition into
independent components
meldata.cc - data handler / container class
Christian F. Beckmann, FMRIB Analysis Group
Copyright (C) 1999-2013 University of Oxford */
/* CCOPYRIGHT */
@@ -12,7 +12,7 @@
#include "newimage/newimageall.h"
#include "meloptions.h"
#include "meldata.h"
#include "melodic.h"
#include "melodic.h"
#include "utils/log.h"
#include <time.h>
#include <algorithm>
@@ -21,15 +21,15 @@
using namespace cifti;
using namespace Utilities;
using namespace NEWIMAGE;
namespace Melodic{
// {{{ Setup
ReturnMatrix MelodicData::process_file(string fname, int numfiles)
{
dbgmsg(string("START: process_file") << endl);
dbgmsg(string("START: process_file") << endl);
Matrix tmpData;
if ( !opts.readCIFTI.value() ) //Process NIFTI
{
@@ -42,9 +42,9 @@ namespace Melodic{
read_volume4D(RawData,fname);
message(" done" << endl);
memmsg(" after reading file "<< fname);
del_vols(RawData,opts.dummy.value());
Mean += meanvol(RawData)/numfiles;
//estimate smoothness
@@ -52,11 +52,11 @@ namespace Melodic{
if((Resels == 0)&&(!opts.filtermode))
Resels = est_resels(RawData,Mask);
memmsg(" after smoothness ");
//reshape
memmsg(" before reshape ");
tmpData = RawData.matrix(Mask);
memmsg(" after reshape ");
memmsg(" after reshape ");
} else { //Read in Cifti
inputCifti.openFile(fname+".nii");
const vector<int64_t>& dims = inputCifti.getDimensions();
@@ -64,33 +64,33 @@ namespace Melodic{
vector<float> scratchRow(dims[0]);//read/write a row at a time
for (int64_t row=0;row<dims[1];row++) {
inputCifti.getRow(scratchRow.data(),row);
for (int64_t col=0;col<dims[0];col++)
for (int64_t col=0;col<dims[0];col++)
tmpData(col+1,row+1)=scratchRow[col];
}
Resels=1;
}
// If a time series model design was specified, check
// If a time series model design was specified, check
// that the data dimensions match the model dimensions
if (Tdes.Storage() && (tmpData.Nrows() != Tdes.Nrows())) {
cerr << "ERROR: " << fname << " " <<
"- data dimensions (" << tmpData.Nrows() << ") " <<
cerr << "ERROR: " << fname << " " <<
"- data dimensions (" << tmpData.Nrows() << ") " <<
"do not match model dimensions (" << Tdes.Nrows() << ")" << endl;
exit(2);
}
}
//convert to percent BOLD signal change
if(opts.pbsc.value()){
message(" Converting data to percent BOLD signal change ...");
Matrix meanimg = convert_to_pbsc(tmpData);
meanR = meanimg.Row(1);
message(" done" << endl);
}
}
else{
if(opts.remove_meanvol.value())
{
{
message(string(" Removing mean image ..."));
memmsg(" before remmean ");
remmean(tmpData,meanR,1);
@@ -103,14 +103,14 @@ namespace Melodic{
if(opts.remove_meantc.value()){
remmean(tmpData,meanC,2);
}
//convert to power spectra
if(opts.pspec.value()){
message(" Converting data to powerspectra ...");
tmpData = calc_FFT(tmpData);
message(" done" << endl);
}
//switch dimension in case temporal ICA is required
if(opts.temporal.value()){
message(string(" Switching dimensions for temporal ICA") << endl);
@@ -121,20 +121,20 @@ namespace Melodic{
meanR = tmp.t();
message(" Data size : " << Data.Nrows() << " x " << Data.Ncols() <<endl);
}
//variance - normalisation
if(opts.varnorm.value()){
memmsg(" before VN ");
message(" Normalising by voxel-wise variance ...");
message(" Normalising by voxel-wise variance ...");
outMsize("stdDev",stdDev);
// if(stdDev.Storage()==0)
stdDev = varnorm(tmpData,std::min(30,tmpData.Nrows()-1),
opts.vn_level.value(), opts.econ.value());
// else
// else
// stdDev += varnorm(tmpData,std::min(30,tmpData.Nrows()-1),
// opts.vn_level.value(), opts.econ.value())/numfiles;
stdDevi = pow(stdDev,-1);
stdDevi = pow(stdDev,-1);
memmsg(" in VN ");
message(" done" << endl);
}
@@ -143,17 +143,17 @@ namespace Melodic{
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);
}
tmpData.Release();
dbgmsg(string("END: process_file") << endl);
dbgmsg(string("END: process_file") << endl);
return tmpData;
}
@@ -204,7 +204,7 @@ namespace Melodic{
outMsize("tmpT",tmpT);
outMsize("tmpS",tmpS);
dbgmsg(string(" approach ") << opts.approach.value() << endl);
dbgmsg(string(" approach ") << opts.approach.value() << endl);
if(opts.approach.value()==string("tica")){
message("Calculating T- and S-modes " << endl);
@@ -264,7 +264,7 @@ namespace Melodic{
void MelodicData::dual_regression()
{
dbgmsg(string("START: dual_regression") << endl);
dbgmsg(string("START: dual_regression") << endl);
Tmodes.clear();
Smodes.clear();
@@ -292,10 +292,10 @@ namespace Melodic{
alltcs=s1;
else
alltcs&=s1;
// output DR
if(opts.dr_out.value()){
dbgmsg(string("START: dual_regression output") << endl);
write_ascii_matrix(drO.appendDir("dr_stage1_subject"+num2str(ctr,4)+".txt"),s1);
//des_norm
@@ -312,7 +312,7 @@ namespace Melodic{
tmpcont << alltcs.Column(ctr);
add_Tmodes(tmpcont);
}
for(int ctrC = 1; ctrC <=IC.Nrows(); ctrC++){
Matrix tmpall = zeros(numfiles,IC.Ncols());
string fnout = string("dr/dr_stage2_ic"+num2str(ctrC-1,4));
@@ -321,31 +321,31 @@ namespace Melodic{
dbgmsg(fnout << endl << fnin << endl);
volume4D<float> vol;
read_volumeROI(vol,fnin,0,0,0,ctrC-1,-1,-1,-1,ctrC-1);
Matrix tmp2 = vol.matrix(Mask);
tmpall.Row(ctrS+1) << vol.matrix(Mask);
}
save4D(tmpall,fnout);
}
opts.varnorm.set_T(tmpvarnorm);
dbgmsg(string("END: dual_regression") << endl);
opts.varnorm.set_T(tmpvarnorm);
dbgmsg(string("END: dual_regression") << endl);
}
void MelodicData::set_TSmode()
{
dbgmsg(string("START: set_TSmode")<< endl);
dbgmsg(string("START: set_TSmode")<< endl);
if(opts.dr.value())
dual_regression();
else
set_TSmode_depr();
dbgmsg(string("END: set_TSmode")<< endl);
dbgmsg(string("END: set_TSmode")<< endl);
}
void MelodicData::setup_classic()
{
dbgmsg(string("START: setup_classic") << endl);
dbgmsg(string("START: setup_classic") << endl);
Matrix alldat, tmpData;
bool tmpvarnorm = opts.varnorm.value();
@@ -363,36 +363,36 @@ namespace Melodic{
for(int ctr = 1; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(ctr), numfiles) / numfiles;
if(tmpData.Ncols() == alldat.Ncols() && tmpData.Nrows() == alldat.Nrows())
alldat = alldat + tmpData;
alldat = alldat + tmpData;
else{
if(opts.approach.value()==string("tica")){
cerr << "ERROR:: data dimensions do not match, TICA not possible \n\n";
exit(2);
exit(2);
}
if(tmpData.Ncols() == alldat.Ncols()){
int mindim = min(alldat.Nrows(),tmpData.Nrows());
alldat = alldat.Rows(1,mindim);
tmpData = tmpData.Rows(1,mindim);
alldat += tmpData;
}
else
alldat += tmpData;
}
else
message("Data dimensions do not match - ignoring "+opts.inputfname.value().at(ctr) << endl);
}
}
}
//update mask
if(opts.update_mask.value()){
message("Excluding voxels with constant value ...");
update_mask(Mask, alldat);
update_mask(Mask, alldat);
message(" done" << endl);
}
if((numfiles > 1 ) && opts.joined_vn.value() && tmpvarnorm){
if((numfiles > 1 ) && opts.joined_vn.value() && tmpvarnorm){
//variance - normalisation
message(endl<<"Normalising jointly by voxel-wise variance ...");
message(endl<<"Normalising jointly by voxel-wise variance ...");
stdDev = varnorm(alldat,alldat.Nrows(),opts.vn_level.value(),opts.econ.value());
stdDevi = pow(stdDev,-1);
stdDevi = pow(stdDev,-1);
message(" done" << endl);
}
@@ -408,7 +408,7 @@ namespace Melodic{
SymmetricMatrix Corr;
int order;
order = ppca_dim(remmean(alldat,2), RXweight, tmpPPCA, AdjEV, PercEV, Corr, pcaE, pcaD, Resels, opts.pca_est.value());
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();
@@ -450,7 +450,7 @@ namespace Melodic{
Matrix tmp = IdentityMatrix(Data.Nrows());
DWM.push_back(tmp);
WM.push_back(tmp);
}
}
else {
dbgmsg("Multi-Subject ICA");
@@ -465,15 +465,15 @@ namespace Melodic{
SP3(tmpData,pow(stdDev,-1));
}
// whiten (separate / joint)
Matrix newWM,newDWM;
if(!opts.joined_whiten.value()){
Matrix newWM,newDWM;
if(!opts.joined_whiten.value()){
message(" Individual whitening in a " << order << " dimensional subspace " << endl);
std_pca(tmpData, RXweight, Corr, pcaE, pcaD, opts.econ.value());
calc_white(pcaE, pcaD, order, newWM, newDWM);
}else{
if(!opts.dr_pca.value()){
std_pca(whiteMatrix*tmpData, RXweight, Corr, pcaE, pcaD, opts.econ.value());
calc_white(pcaE, pcaD, order, newWM, newDWM);
calc_white(pcaE, pcaD, order, newWM, newDWM);
newDWM=(dewhiteMatrix*newDWM);
newWM=(newWM*whiteMatrix);
}
@@ -484,9 +484,9 @@ namespace Melodic{
tmp1 = whiteMatrix * alldat;
remmean(tmp1,2);
tmp1 *= tmpData.t();
tmp2 = MISCMATHS::pinv(tmp1.t()).t();
tmp2 = MISCMATHS::pinv(tmp1.t()).t();
std_pca(tmp1 * tmpData, RXweight, Corr, pcaE, pcaD, opts.econ.value());
calc_white(pcaE, pcaD, order, newWM, newDWM);
calc_white(pcaE, pcaD, order, newWM, newDWM);
newDWM=(tmp2*newDWM);
newWM=(newWM * tmp1);
}
@@ -503,38 +503,38 @@ namespace Melodic{
}
}
opts.varnorm.set_T(tmpvarnorm);
dbgmsg(string("END: setup_classic") << endl);
dbgmsg(string("END: setup_classic") << endl);
}
void MelodicData::setup_migp()
{
dbgmsg(string("START: setup_migp") << endl);
dbgmsg(string("START: setup_migp") << endl);
std::vector<int> myctr;
for (int i=0; i< numfiles ; ++i) myctr.push_back(i);
for (int i=0; i< numfiles ; ++i) myctr.push_back(i);
if(opts.migp_shuffle.value()){
message("Randomising input file order" << endl);
std::random_shuffle ( myctr.begin(), myctr.end() );
}
Matrix tmpData;
bool tmpvarnorm = opts.varnorm.value();
if(numfiles > 1 && opts.joined_vn.value()){
opts.varnorm.set_T(false);
}
for(int ctr = 0; ctr < numfiles; ctr++){
tmpData = process_file(opts.inputfname.value().at(myctr.at(ctr)), numfiles) / numfiles;
if (opts.migpN.value()==0){
opts.migpN.set_T(2*tmpData.Nrows()-1);
}
if(opts.debug.value())
save4D(tmpData,string("preproc_dat") + num2str(ctr+1));
if(Data.Storage()==0)
Data = tmpData;
else
@@ -549,16 +549,16 @@ namespace Melodic{
RowVector pcaD;
std_pca(Data, RXweight, Corr, pcaE, pcaD, opts.econ.value());
pcaE = pcaE.Columns(pcaE.Ncols()-opts.migpN.value()+1,pcaE.Ncols());
Data = pcaE.t() * Data;
Data = pcaE.t() * Data;
}
outMsize("Data", Data);
}
//update mask
if(opts.update_mask.value()){
message(endl<< "Excluding voxels with constant value ...");
update_mask(Mask, Data);
update_mask(Mask, Data);
message(" done" << endl);
}
@@ -568,30 +568,30 @@ namespace Melodic{
opts.varnorm.set_T(tmpvarnorm);
if(opts.varnorm2.value()){
message(" Normalising by voxel-wise variance ...");
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);
dbgmsg(string("END: setup_migp") << endl);
}
void MelodicData::setup()
{
dbgmsg(string("START: setup") << endl);
{
dbgmsg(string("START: setup") << endl);
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.set_T("symm");
opts.approach.set_T("symm");
if(opts.deflation.value())
opts.approach.set_T("defl");
opts.migp.set_T(false);
@@ -601,25 +601,25 @@ namespace Melodic{
if( opts.approach.value()==string("concat") && opts.migp.value() )
setup_migp();
else
setup_classic();
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");
save4D(Data,"concat_data");
//save the mean & mask
if ( !opts.readCIFTI.value() ) {
save_volume(Mask,logger.appendDir("mask"));
save_volume(Mean,logger.appendDir("mean"));
}
dbgmsg(string("END: setup") << endl);
dbgmsg(string("END: setup") << endl);
} // void setup()
void MelodicData::setup_misc()
{
dbgmsg(string("START: setup_misc") << endl);
dbgmsg(string("START: setup_misc") << endl);
if (!opts.readCIFTI.value()) {
//initialize Mean
read_volumeROI(Mean,opts.inputfname.value().at(0),-1,-1,-1,0,-1,-1,-1,0);
@@ -637,7 +637,7 @@ namespace Melodic{
}else{
background = Mean;
}
if(!samesize(Mean,Mask,3)){
cerr << "ERROR:: mask and data have different dimensions \n\n";
exit(2);
@@ -645,7 +645,7 @@ namespace Melodic{
//reset mean
Mean *= 0;
//set up weighting
if(opts.segment.value().length()>0){
create_RXweight();
@@ -656,32 +656,32 @@ namespace Melodic{
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);
}
insta_mask = tmp_im.matrix(Mask);
}
}
//seed the random number generator
double tmptime = time(NULL);
if ( opts.seed.value() != -1 ) {
tmptime = opts.seed.value();
tmptime = opts.seed.value();
}
srand((unsigned int) tmptime);
if(opts.paradigmfname.value().length()>0){
message(" Use columns in " << opts.paradigmfname.value()
message(" Use columns in " << opts.paradigmfname.value()
<< " for PCA initialisation" <<endl);
param = read_ascii_matrix(opts.paradigmfname.value());
Matrix tmpPU, tmpPV;
DiagonalMatrix tmpPD;
SVD(param, tmpPD, tmpPU, tmpPV);
param = tmpPU;
opts.pca_dim.set_T(std::max(opts.pca_dim.value(), param.Ncols()+3));
opts.pca_dim.set_T(std::max(opts.pca_dim.value(), param.Ncols()+3));
if(opts.debug.value()){
outMsize("Paradigm",param); saveascii(param,"param");
}
@@ -702,7 +702,7 @@ namespace Melodic{
if(opts.fn_SconF.value().length()>0)
SconF = read_ascii_matrix(opts.fn_SconF.value());
// Check that the number of input
// Check that the number of input
// files matches the session design
if (Sdes.Storage()) {
if (Sdes.Nrows() != numfiles) {
@@ -710,11 +710,11 @@ namespace Melodic{
"does not match subject/session design (" << Sdes.Nrows() << ")" << endl;
exit(2);
}
}
// Or create a default session design
}
// Or create a default session design
// if one was not specified
else if(numfiles>1){
else if(numfiles>1){
Sdes = ones(numfiles,1);
if(Scon.Storage() == 0){
Scon = ones(1,1);
@@ -722,13 +722,13 @@ namespace Melodic{
}
}
remmean(Tdes);
dbgmsg(string("END: setup_misc") << endl);
dbgmsg(string("END: setup_misc") << endl);
}
void MelodicData::save()
{
{
//check for temporal ICA
if(opts.temporal.value()){
@@ -747,24 +747,24 @@ namespace Melodic{
message(string("done") << endl);
opts.temporal.set_T(false); // Do not switch again!
}
message(endl << "Writing results to : " << endl);
//Output IC
//Output IC
if((IC.Storage()>0)&&(opts.output_origIC.value())&&(after_mm==false))
save4D(IC,opts.outputfname.value() + "_oIC");
//Output IC -- adjusted for noise
//Output IC -- adjusted for noise
if(IC.Storage()>0){
volume4D<float> tempVol;
volume4D<float> tempVol;
//Matrix ICadjust;
if(after_mm){
save4D(IC,opts.outputfname.value() + "_IC");
// ICadjust = IC;
}
}
else{
Matrix resids = stdev(Data - mixMatrix * IC);
for(int ctr=1;ctr<=resids.Ncols();ctr++)
if(resids(1,ctr) < 0.05)
@@ -775,19 +775,19 @@ namespace Melodic{
stdNoisei = pow(resids*
std::sqrt((float)(Data.Nrows()-1))/
std::sqrt((float)(Data.Nrows()-IC.Nrows())),-1);
ColumnVector diagvals;
diagvals=pow(diag(unmixMatrix*unmixMatrix.t()),-0.5);
save4D(SP(IC,diagvals*stdNoisei),opts.outputfname.value() + "_IC");
}
if(opts.output_origIC.value())
save4D(stdNoisei,string("Noise__inv"));
}
//Output T- & S-modes
save_Tmodes();
save_Smodes();
@@ -795,17 +795,17 @@ namespace Melodic{
if(mixMatrix.Storage()>0){
saveascii(expand_mix(), opts.outputfname.value() + "_mix");
mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
saveascii(mixFFT,opts.outputfname.value() + "_FTmix");
saveascii(mixFFT,opts.outputfname.value() + "_FTmix");
}
//Output PPCA
if(PPCA.Storage()>0)
saveascii(PPCA, opts.outputfname.value() + "_PPCA");
//Output ICstats
if(ICstats.Storage()>0)
saveascii(ICstats,opts.outputfname.value() + "_ICstats");
saveascii(ICstats,opts.outputfname.value() + "_ICstats");
//Output unmixMatrix
if(opts.output_unmix.value() && unmixMatrix.Storage()>0)
saveascii(unmixMatrix,opts.outputfname.value() + "_unmix");
@@ -833,7 +833,7 @@ namespace Melodic{
if(opts.output_pca.value() && pcaD.Storage()>0&&pcaE.Storage()>0){
saveascii(pcaE,opts.outputfname.value() + "_pcaE");
saveascii((Matrix) diag(pcaD),opts.outputfname.value() + "_pcaD");
Matrix PCAmaps;
if(whiteMatrix.Ncols()==Data.Ncols())
PCAmaps = dewhiteMatrix.t();
@@ -842,19 +842,19 @@ namespace Melodic{
save4D(PCAmaps,opts.outputfname.value() + "_pca");
}
message("...done" << endl);
} //void save()
int MelodicData::remove_components()
{
message("Reading " << opts.filtermix.value() << endl)
{
message("Reading " << opts.filtermix.value() << endl)
mixMatrix = read_ascii_matrix(opts.filtermix.value());
if (mixMatrix.Storage()<=0) {
cerr <<" Please specify the mixing matrix correctly" << endl;
exit(2);
}
unmixMatrix = pinv(mixMatrix);
IC = unmixMatrix * Data;
@@ -864,11 +864,11 @@ namespace Melodic{
Matrix noiseMix;
Matrix noiseIC;
int ctr=0;
int ctr=0;
char *p;
char t[1024];
const char *discard = ", [];{(})abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ~!@#$%^&*_-=+|\':><./?";
message("Filtering the data...");
strcpy(t, tmpstr.c_str());
p=strtok(t,discard);
@@ -876,17 +876,17 @@ namespace Melodic{
if(ctr>0 && ctr<=mixMatrix.Ncols()){
message(" "<< ctr );
noiseMix = mixMatrix.Column(ctr);
noiseIC = IC.Row(ctr).t();
noiseIC = IC.Row(ctr).t();
}
else{
cerr << endl<< "component number "<<ctr<<" does not exist" << endl;
}
do{
p=strtok(NULL,discard);
if(p){
ctr = atoi(p);
if(ctr>0 && ctr<=mixMatrix.Ncols()){
message(" "<<ctr);
noiseMix |= mixMatrix.Column(ctr);
@@ -911,23 +911,23 @@ namespace Melodic{
if(meanR.Storage()>0)
newData = newData + ones(newData.Nrows(),1)*meanR;
volume4D<float> tmp;
read_volume4D(tmp,opts.inputfname.value().at(0));
read_volume4D(tmp,opts.inputfname.value().at(0));
tmp.setmatrix(newData,Mask);
save_volume4D(tmp,logger.appendDir(opts.outputfname.value() + "_ICAfiltered"));
save_volume4D(tmp,logger.appendDir(opts.outputfname.value() + "_ICAfiltered"));
return 0;
} // int remove_components()
void MelodicData::create_RXweight()
{
message("Reading the weights for the covariance R_X from file "<< opts.segment.value() << endl);
volume4D<float> tmpRX;
read_volume4D(tmpRX,opts.segment.value());
RXweight = tmpRX.matrix(Mask);
}
}
void MelodicData::est_smoothness()
{
@@ -937,15 +937,15 @@ namespace Melodic{
if(opts.segment.value().length()>0){
Mask_fname = opts.segment.value();
}
}
// Setup external call to smoothest:
char callSMOOTHESTstr[1000];
ostrstream osc(callSMOOTHESTstr,1000);
osc << SM_path << " -d " << data_dim()
<< " -r " << opts.inputfname.value().at(0) << " -m "
<< " -r " << opts.inputfname.value().at(0) << " -m "
<< Mask_fname << " > " << logger.appendDir("smoothest") << '\0';
message(" Calling Smoothest: " << callSMOOTHESTstr << endl);
system(callSMOOTHESTstr);
@@ -953,9 +953,9 @@ namespace Melodic{
ifstream in;
string str;
Resels = 1.0;
in.open(logger.appendDir("smoothest").c_str(), ios::in);
if(in>0){
if(in.good()){
for(int ctr=1; ctr<7; ctr++)
in >> str;
in.close();
@@ -969,25 +969,25 @@ namespace Melodic{
{
unsigned long count = 0;
int M=R.tsize();
for (int z=mask.minz(); z<=mask.maxz(); z++) {
for (int y=mask.miny(); y<=mask.maxy(); y++) {
for (int x=mask.minx(); x<=mask.maxx(); x++) {
if( mask(x,y,z) > 0.5) {
count ++;
if( M > 2 ) {
// For each voxel
// For each voxel
// calculate mean and standard deviation...
double Sx = 0.0, SSx = 0.0;
double Sx = 0.0, SSx = 0.0;
for ( int t = 0; t < M; t++ ) {
float R_it = R(x,y,z,t);
Sx += R_it;
SSx += (R_it)*(R_it);
}
float mean = Sx / M;
float sdsq = (SSx - ((Sx)*(Sx) / M)) / (M - 1) ;
if (sdsq<=0) {
// trap for differences between mask and invalid data
mask(x,y,z)=0;
@@ -997,12 +997,12 @@ namespace Melodic{
for ( unsigned short t = 0; t < M; t++ ) {
R(x,y,z,t) = (R(x,y,z,t) - mean) / sqrt(sdsq);
}
}
}
}
}
}
}
}
}
return count;
}
@@ -1028,47 +1028,47 @@ namespace Melodic{
for ( unsigned short x = 1; x < R.xsize() ; x++ )
// Sum over N
if( (mask(x, y, z)>0.5) &&
(mask(x-1, y, z)>0.5) &&
(mask(x, y-1, z)>0.5) &&
(mask(x-1, y, z)>0.5) &&
(mask(x, y-1, z)>0.5) &&
( (!usez) || (mask(x, y, z-1)>0.5) ) ) {
N++;
for ( unsigned short t = 0; t < R.tsize(); t++ ) {
// Sum over M
SSminus[X] += R(x, y, z, t) * R(x-1, y, z, t);
SSminus[Y] += R(x, y, z, t) * R(x, y-1, z, t);
if (usez) SSminus[Z] += R(x, y, z, t) * R(x, y, z-1, t);
S2[X] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) +
S2[X] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) +
R(x-1, y, z, t)*R(x-1, y, z, t));
S2[Y] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) +
S2[Y] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) +
R(x, y-1, z, t)*R(x, y-1, z, t));
if (usez) S2[Z] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) +
if (usez) S2[Z] += 0.5 * (R(x, y, z, t)*R(x, y, z, t) +
R(x, y, z-1, t)*R(x, y, z-1, t));
}
}
double norm = 1.0/(double) N;
double v = dof; // v - degrees of freedom (nu)
double v = dof; // v - degrees of freedom (nu)
if(R.tsize() > 1) {
norm = (v - 2) / ((v - 1) * N * R.tsize());
}
// for extreme smoothness
if (SSminus[X]>=0.99999999*S2[X])
SSminus[X]=0.99999*S2[X];
if (SSminus[Y]>=0.99999999*S2[Y])
// for extreme smoothness
if (SSminus[X]>=0.99999999*S2[X])
SSminus[X]=0.99999*S2[X];
if (SSminus[Y]>=0.99999999*S2[Y])
SSminus[Y]=0.99999*S2[Y];
if (usez)
if (SSminus[Z]>=0.99999999*S2[Z])
if (usez)
if (SSminus[Z]>=0.99999999*S2[Z])
SSminus[Z]=0.99999*S2[Z];
// Convert to sigma squared
double sigmasq[3] = {0,0,0};
sigmasq[X] = -1.0 / (4 * log(fabs(SSminus[X]/S2[X])));
sigmasq[Y] = -1.0 / (4 * log(fabs(SSminus[Y]/S2[Y])));
if (usez) { sigmasq[Z] = -1.0 / (4 * log(fabs(SSminus[Z]/S2[Z]))); }
// Convert to full width half maximum
double FWHM[3] = {0,0,0};
FWHM[X] = sqrt(8 * log(2) * sigmasq[X]);
@@ -1083,7 +1083,7 @@ namespace Melodic{
void MelodicData::create_mask(volume<float>& theMask)
{
if(opts.use_mask.value() && opts.maskfname.value().size()>0){ // mask provided
if(opts.use_mask.value() && opts.maskfname.value().size()>0){ // mask provided
read_volume(theMask,opts.maskfname.value());
message("Mask provided : " << opts.maskfname.value()<<endl<<endl);
}
@@ -1092,7 +1092,7 @@ namespace Melodic{
message("Create mask ... ");
//save first image
tmpnam(Mean_fname); // generate a tmp name
save_volume(Mean,Mean_fname);
save_volume(Mean,Mean_fname);
// set up all strings
string BET_outputfname = string(Mean_fname)+"_brain";
@@ -1105,34 +1105,34 @@ namespace Melodic{
// char callBETstr[1000];
// ostrstream betosc(callBETstr,1000);
// betosc << BET_path << " " << Mean_fname << " "
// betosc << BET_path << " " << Mean_fname << " "
// << BET_outputfname << " " << BET_optarg << " > /dev/null " << '\0';
// message(" Calling BET: " << callBETstr << endl);
// system(callBETstr);
string tmpstr = BET_path + string(" ") +
Mean_fname + string(" ") + BET_outputfname + string(" ") +
string tmpstr = BET_path + string(" ") +
Mean_fname + string(" ") + BET_outputfname + string(" ") +
BET_optarg + string(" > /dev/null ");
system(tmpstr.c_str());
// read back the Mask file
// read back the Mask file
read_volume(theMask,Mask_fname);
// clean /tmp
char callRMstr[1000];
ostrstream osc(callRMstr,1000);
osc << "rm " << string(Mean_fname) <<"* " << '\0';
system(callRMstr);
message("done" << endl);
}
}
else{
if(opts.use_mask.value()){ //just threshold the Mean
message("Create mask ... ");
float Mmin, Mmax, Mtmp;
Mmin = Mean.min(); Mmax = Mean.max();
theMask = binarise(Mean,Mmin + opts.threshold.value()*
theMask = binarise(Mean,Mmin + opts.threshold.value()*
(Mmax-Mmin),Mmax);
Mtmp = Mmin + opts.threshold.value()* (Mmax-Mmin);
message("done" << endl);
@@ -1143,11 +1143,11 @@ namespace Melodic{
}
}
}
if(opts.remove_endslices.value()){
if(opts.remove_endslices.value()){
// just in case mc introduced something nasty
message(" Deleting end slices" << endl);
for(int ctr1=theMask.miny(); ctr1<=theMask.maxy(); ctr1++){
for(int ctr2=theMask.minx(); ctr2<=theMask.maxx(); ctr2++){
for(int ctr2=theMask.minx(); ctr2<=theMask.maxx(); ctr2++){
theMask(ctr2,ctr1,Mask.minz()) = 0.0;
theMask(ctr2,ctr1,Mask.maxz()) = 0.0;
}
@@ -1157,7 +1157,7 @@ namespace Melodic{
void MelodicData::sort()
{
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
int numComp = mixMatrix.Ncols(), numVox = IC.Ncols(),
numTime = mixMatrix.Nrows(), i,j;
//flip IC maps to be positive (on max)
@@ -1166,21 +1166,21 @@ namespace Melodic{
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++)
if(IC.Row(ctr_i).MaximumAbsoluteValue()>IC.Row(ctr_i).Maximum()){
flipres(ctr_i);
flipres(ctr_i);
}
message("Sorting IC maps" << endl);
message("Sorting IC maps" << endl);
Matrix tmpscales, tmpICrow, tmpMIXcol;
if(numfiles > 1 && opts.approach.value()==string("tica")){
set_TSmode();
Matrix allmodes = Smodes.at(0);
for(int ctr = 1; ctr < (int)Smodes.size();++ctr)
allmodes |= Smodes.at(ctr);
tmpscales = median(allmodes).t();
tmpscales = median(allmodes).t();
} else {
// re-order wrt standard deviation of IC maps
// re-order wrt standard deviation of IC maps
tmpscales = stdev(IC,2);
}
ICstats = tmpscales;
double max_val, min_val = tmpscales.Minimum()-1;
@@ -1188,34 +1188,34 @@ namespace Melodic{
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++){
max_val = tmpscales.Maximum2(i,j);
ICstats(ctr_i,1)=max_val;
tmpICrow = IC.Row(ctr_i);
tmpMIXcol = mixMatrix.Column(ctr_i);
IC.SubMatrix(ctr_i,ctr_i,1,numVox) = IC.SubMatrix(i,i,1,numVox);
mixMatrix.SubMatrix(1,numTime,ctr_i,ctr_i) =
mixMatrix.SubMatrix(1,numTime,ctr_i,ctr_i) =
mixMatrix.SubMatrix(1,numTime,i,i);
IC.SubMatrix(i,i,1,numVox) = tmpICrow.SubMatrix(1,1,1,numVox);
mixMatrix.SubMatrix(1,numTime,i,i) = tmpMIXcol.SubMatrix(1,numTime,1,1);
tmpscales(i,1)=tmpscales(ctr_i,1);
tmpscales(ctr_i,1)=min_val;
}
ICstats /= ICstats.Column(1).Sum();
ICstats *= 100;
if(EVP.Storage()>0){
tmpscales = ICstats.Column(1).AsMatrix(ICstats.Nrows(),1) * EVP(1,numComp);
ICstats |= tmpscales;
}
if(Data.Storage()>0&&stdDev.Storage()>0){
Matrix copeP(tmpscales), copeN(tmpscales);
Matrix max_ICs(tmpscales), min_ICs(tmpscales);
for(int ctr_i = 1; ctr_i <= numComp; ctr_i++){
int i,j;
max_ICs(ctr_i,1) = IC.Row(ctr_i).Maximum2(i,j);
@@ -1233,7 +1233,7 @@ namespace Melodic{
ICstats |= copeP;
ICstats |= copeN;
}
mixFFT=calc_FFT(expand_mix(), opts.logPower.value());
unmixMatrix = pinv(mixMatrix);
}
@@ -1249,11 +1249,7 @@ namespace Melodic{
if(mixMatrix.Storage()>0){cout << "mix: " << mixMatrix.Nrows() <<"x" << mixMatrix.Ncols() << endl;}else{cout << "mix empty " <<endl;}
if(unmixMatrix.Storage()>0){cout << "unmix: " << unmixMatrix.Nrows() <<"x" << unmixMatrix.Ncols() << endl;}else{cout << "unmix empty " <<endl;}
if(IC.Storage()>0){cout << "IC: " << IC.Nrows() <<"x" << IC.Ncols() << endl;}else{cout << "IC empty " <<endl;}
} //void status()
}