Skip to content
Snippets Groups Projects
Commit e50138b4 authored by Mark Woolrich's avatar Mark Woolrich
Browse files

*** empty log message ***

parent dda165d0
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
AutoCorrEstimator.cc
+ 186
− 28
View file @ e50138b4
...@@ -8,12 +8,14 @@ ...@@ -8,12 +8,14 @@
#include <iostream> #include <iostream>
#include <strstream> #include <strstream>
#define WANT_STREAM
#include "AutoCorrEstimator.h" #include "AutoCorrEstimator.h"
#include "miscmaths.h" #include "miscmaths.h"
#include "Log.h" #include "Log.h"
#include "Volume.h" #include "Volume.h"
#include "histogram.h" #include "histogram.h"
#include "miscmaths.h"
using namespace NEWMAT; using namespace NEWMAT;
using namespace UTILS; using namespace UTILS;
...@@ -40,7 +42,7 @@ namespace TACO { ...@@ -40,7 +42,7 @@ namespace TACO {
{ {
dummy = 0; dummy = 0;
dmrow = 0; dmrow = 0;
mn(k) = MISCMATHS::mean(dm.Column(k)); mn(k) = MISCMATHS::mean(ColumnVector(dm.Column(k)));
dmrow.Rows(1,sizeTS) = dm.Column(k) - mn(k); dmrow.Rows(1,sizeTS) = dm.Column(k) - mn(k);
FFT(dmrow, dummy, dm_fft_real, dm_fft_imag); FFT(dmrow, dummy, dm_fft_real, dm_fft_imag);
...@@ -50,7 +52,7 @@ namespace TACO { ...@@ -50,7 +52,7 @@ namespace TACO {
} }
} }
void AutoCorrEstimator::preWhiten(ColumnVector& in, ColumnVector& ret, int i, Matrix& dmret) { void AutoCorrEstimator::preWhiten(ColumnVector& in, ColumnVector& ret, int i, Matrix& dmret, bool highfreqremovalonly) {
Tracer tr("AutoCorrEstimator::preWhiten"); Tracer tr("AutoCorrEstimator::preWhiten");
int sizeTS = xdata.getNumVolumes(); int sizeTS = xdata.getNumVolumes();
...@@ -66,6 +68,25 @@ namespace TACO { ...@@ -66,6 +68,25 @@ namespace TACO {
vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = acEst.getSeries(i).Rows(2, sizeTS/2).Reverse(); vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = acEst.getSeries(i).Rows(2, sizeTS/2).Reverse();
FFT(vrow, dummy, ac_fft_real, ac_fft_im); FFT(vrow, dummy, ac_fft_real, ac_fft_im);
float norm = ac_fft_real.SumSquare();
// Compare with raw FFT to detect high frequency artefacts:
bool violate = false;
ColumnVector violators(zeropad);
violators = 1;
for(int j = 1; j <= zeropad; j++)
{
E(j,i) = sqrt(E(j,i)/((ac_fft_real(j)*ac_fft_real(j))/norm));
// look for high frequency artefacts
if(E(j,i) > 4 && j > zeropad/4 && j < 3*zeropad/4)
{
violate = true;
violators(j) = 0;
countLargeE(j) = countLargeE(j) + 1;
}
}
// FFT x data // FFT x data
dummy = 0; dummy = 0;
...@@ -74,14 +95,25 @@ namespace TACO { ...@@ -74,14 +95,25 @@ namespace TACO {
FFT(xrow, dummy, x_fft_real, x_fft_im); FFT(xrow, dummy, x_fft_real, x_fft_im);
// inverse auto corr to give prewhitening filter: if(highfreqremovalonly)
// no DC component so set first value to 0;
ac_fft_real(1) = 0.0;
for(int j = 2; j <= zeropad; j++)
{ {
ac_fft_real(j) = 1.0/ac_fft_real(j); ac_fft_real = violators;
}
else
{
// inverse auto corr to give prewhitening filter
// no DC component so set first value to 0
ac_fft_real(1) = 0.0;
for(int j = 2; j <= zeropad; j++)
{
ac_fft_real(j) = 1.0/sqrt(ac_fft_real(j));
}
} }
// normalise ac_fft such that sum(j)(ac_fft_real)^2 = 1
ac_fft_real /= sqrt(ac_fft_real.SumSquare()/zeropad);
// filter design matrix // filter design matrix
for(int k = 1; k <= numPars; k++) for(int k = 1; k <= numPars; k++)
{ {
...@@ -90,12 +122,12 @@ namespace TACO { ...@@ -90,12 +122,12 @@ namespace TACO {
FFTI(SP(ac_fft_real, dm_fft_real), SP(ac_fft_real, dm_fft_imag), realifft, dummy); FFTI(SP(ac_fft_real, dm_fft_real), SP(ac_fft_real, dm_fft_imag), realifft, dummy);
// place result into ret: // place result into ret:
dmret.Column(k) = realifft.Rows(1,sizeTS); dmret.Column(k) = realifft.Rows(1,sizeTS) + mn(k);
float std = pow(MISCMATHS::var(dmret.Column(k)),0.5); //float std = pow(MISCMATHS::var(ColumnVector(dmret.Column(k))),0.5);
dmret.Column(k) = (dmret.Column(k)/std) + mn(k); //dmret.Column(k) = (dmret.Column(k)/std) + mn(k);
} }
// Do filtering and inverse FFT: // Do filtering of data:
FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy); FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy);
// place result into ret: // place result into ret:
...@@ -133,7 +165,7 @@ namespace TACO { ...@@ -133,7 +165,7 @@ namespace TACO {
vrow.Rows(1,sizeTS/2) = acEst.getSeries(i).Rows(1,sizeTS/2); vrow.Rows(1,sizeTS/2) = acEst.getSeries(i).Rows(1,sizeTS/2);
vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = acEst.getSeries(i).Rows(2, sizeTS/2).Reverse(); vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = acEst.getSeries(i).Rows(2, sizeTS/2).Reverse();
FFT(vrow, dummy, ac_fft_real, ac_fft_im); FFT(vrow, dummy, ac_fft_real, ac_fft_im);
// FFT x data // FFT x data
dummy = 0; dummy = 0;
xrow = 0; xrow = 0;
...@@ -144,11 +176,15 @@ namespace TACO { ...@@ -144,11 +176,15 @@ namespace TACO {
// inverse auto corr to give prewhitening filter: // inverse auto corr to give prewhitening filter:
// no DC component so set first value to 0; // no DC component so set first value to 0;
ac_fft_real(1) = 0.0; ac_fft_real(1) = 0.0;
for(int j = 2; j <= zeropad; j++) for(int j = 2; j <= zeropad; j++)
{ {
ac_fft_real(j) = 1.0/ac_fft_real(j); ac_fft_real(j) = 1.0/ac_fft_real(j);
} }
// normalise ac_fft such that sum(j)(ac_fft_real)^2 = 1
ac_fft_real /= sqrt(ac_fft_real.SumSquare()/zeropad);
// Do filtering and inverse FFT: // Do filtering and inverse FFT:
FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy); FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy);
...@@ -223,6 +259,7 @@ namespace TACO { ...@@ -223,6 +259,7 @@ namespace TACO {
} }
Log::getInstance().out("order", order); Log::getInstance().out("order", order);
countLargeE = 0;
cerr << " Completed" << endl; cerr << " Completed" << endl;
} }
...@@ -235,7 +272,7 @@ namespace TACO { ...@@ -235,7 +272,7 @@ namespace TACO {
hist.generate(); hist.generate();
float mode = hist.mode(); float mode = hist.mode();
cerr << "mode = " << mode << endl; cerr << "mode = " << mode << endl;
float sum = 0.0; float sum = 0.0;
int count = 0; int count = 0;
...@@ -255,14 +292,17 @@ namespace TACO { ...@@ -255,14 +292,17 @@ namespace TACO {
return usanthresh; return usanthresh;
} }
void AutoCorrEstimator::spatiallySmooth(const string& usanfname, const Volume& epivol, int masksize, const string& epifname, const string& susanpath) { void AutoCorrEstimator::spatiallySmooth(const string& usanfname, const Volume& epivol, int masksize, const string& epifname, const string& susanpath, int usanthresh) {
Tracer trace("AutoCorrEstimator::spatiallySmooth"); Tracer trace("AutoCorrEstimator::spatiallySmooth");
Log& logger = Log::getInstance(); Log& logger = Log::getInstance();
// Establish epi thresh to use: if(usanthresh == 0)
int usanthresh = establishUsanThresh(epivol); {
// Establish epi thresh to use:
usanthresh = establishUsanThresh(epivol);
}
// Setup external call to susan program: // Setup external call to susan program:
char callsusanstr[1000]; char callsusanstr[1000];
ostrstream osc3(callsusanstr,1000); ostrstream osc3(callsusanstr,1000);
...@@ -288,7 +328,7 @@ namespace TACO { ...@@ -288,7 +328,7 @@ namespace TACO {
cerr << callsusanstr << endl; cerr << callsusanstr << endl;
for(; i < 20; i++) for(; i < 30; i++)
{ {
// output unsmoothed estimates: // output unsmoothed estimates:
vol = acEst.getVolume(i).AsColumn()*factor; vol = acEst.getVolume(i).AsColumn()*factor;
...@@ -310,10 +350,9 @@ namespace TACO { ...@@ -310,10 +350,9 @@ namespace TACO {
char rmfilesstr[1000]; char rmfilesstr[1000];
ostrstream osc(rmfilesstr,1000); ostrstream osc(rmfilesstr,1000);
osc << "rm -rf " osc << "rm -rf "
<< logger.getDir() + "/" + postSmoothVol + "* " << logger.getDir() + "/" + postSmoothVol + "* "
<< logger.getDir() + "/" + preSmoothVol + "* " << logger.getDir() + "/" + preSmoothVol + "* "
<< logger.getDir() + "/" + epifname + "* " << logger.getDir() + "/usanSize*" << '\0';
<< logger.getDir() + "/usanSize*" << '\0';
cerr << rmfilesstr << endl; cerr << rmfilesstr << endl;
...@@ -326,7 +365,7 @@ namespace TACO { ...@@ -326,7 +365,7 @@ namespace TACO {
cerr << "Calculating raw AutoCorrs..."; cerr << "Calculating raw AutoCorrs...";
AutoCorr(xdata, acEst, zeropad); AutoCorr(xdata, acEst, E, zeropad);
cerr << " Completed" << endl; cerr << " Completed" << endl;
} }
...@@ -369,7 +408,124 @@ namespace TACO { ...@@ -369,7 +408,124 @@ namespace TACO {
cerr << " Completed" << endl; cerr << " Completed" << endl;
} }
void AutoCorrEstimator::multitaper(int M) {
Tracer tr("AutoCorrEstimator::multitaper");
cerr << "Multitapering... ";
int sizeTS = xdata.getNumVolumes();
int numTS = xdata.getNumSeries();
Matrix slepians;
getSlepians(M, sizeTS, slepians);
//Log::getInstance().out("slepians", slepians, false);
ColumnVector x(zeropad);
x = 0;
ColumnVector fft_real;
ColumnVector fft_im;
ColumnVector dummy(zeropad);
ColumnVector dummy2;
ColumnVector realifft(zeropad);
dummy = 0;
Matrix Sk(zeropad, slepians.Ncols());
acEst.ReSize(sizeTS, numTS);
acEst = 0;
for(int i = 1; i <= numTS; i++)
{
// Compute FFT for each slepian taper
for(int k = 1; k <= slepians.Ncols(); k++)
{
x.Rows(1,sizeTS) = SP(slepians.Column(k), xdata.getSeries(i));
//Log::getInstance().out("x", xdata.getSeries(i), false);
FFT(x, dummy, fft_real, fft_im);
for(int j = 1; j <= zeropad; j++)
{
// (x+iy)(x-iy) = x^2 + y^2
fft_real(j) = fft_real(j)*fft_real(j) + fft_im(j)*fft_im(j);
Sk(j,k) = fft_real(j);
}
}
// Pool multitaper FFTs
fft_im = 0;
for(int j = 1; j <= zeropad; j++)
{
fft_real(j) = MISCMATHS::mean(ColumnVector(Sk.Row(j).t()));
}
// IFFT to get autocorr
FFTI(fft_real, fft_im, realifft, dummy2);
//Log::getInstance().out("Sk", Sk, false);
//Log::getInstance().out("realifft", realifft);
//Log::getInstance().out("fftreal", fft_real);
float varx = MISCMATHS::var(ColumnVector(x.Rows(1,sizeTS)));
acEst.getSeries(i) = realifft.Rows(1,sizeTS)/varx;
}
countLargeE = 0;
cerr << "Completed" << endl;
}
void AutoCorrEstimator::getSlepians(int M, int sizeTS, Matrix& slepians) {
Tracer tr("AutoCorrEstimator::getSlepians");
slepians.ReSize(sizeTS, 2*M);
ifstream in;
char strc[10];
ostrstream osc(strc,10);
osc << sizeTS << "_" << M << '\0';
string fname("/usr/people/woolrich/parads/mt_" + string(strc));
in.open(fname.c_str(), ios::in);
if(!in)
throw Exception("Multitapering: Slepians file not found");
for(int j = 1; j <= sizeTS; j++)
{
for(int i = 1; i <= 2*M; i++)
{
in >> slepians(j,i);
}
}
in.close();
}
void AutoCorrEstimator::tukey(int M) {
Tracer tr("AutoCorrEstimator::tukey");
cerr << "Tukey M = " << M << endl;
cerr << "Tukey estimates... ";
int sizeTS = xdata.getNumVolumes();
ColumnVector window(M);
for(int j = 1; j <= M; j++)
{
window(j) = 0.5*(1+cos(M_PI*j/(float(M))));
}
for(int i = 1; i <= xdata.getNumSeries(); i++) {
acEst.SubMatrix(1,M,i,i) = SP(acEst.SubMatrix(1,M,i,i),window);
acEst.SubMatrix(M+1,sizeTS,i,i) = 0;
}
countLargeE = 0;
cerr << "Completed" << endl;
}
void AutoCorrEstimator::pava() { void AutoCorrEstimator::pava() {
Tracer tr("AutoCorrEstimator::pava"); Tracer tr("AutoCorrEstimator::pava");
...@@ -437,7 +593,7 @@ namespace TACO { ...@@ -437,7 +593,7 @@ namespace TACO {
else if(j > 2) else if(j > 2)
//if(j > 2) //if(j > 2)
{ {
int endst = j; int endst = j;
int stst = j-(int)(1+(j/8.0)); int stst = j-(int)(1+(j/8.0));
...@@ -446,11 +602,13 @@ namespace TACO { ...@@ -446,11 +602,13 @@ namespace TACO {
for(j = stst; j <= endst; j++) for(j = stst; j <= endst; j++)
{ {
acEst(j,i) = acEst(j,i)*exp(-pow((j-stst)/expwidth,exppow)); acEst(j,i) = acEst(j,i)*exp(-pow((j-stst)/float(expwidth),exppow));
} }
} }
} }
countLargeE = 0;
cerr << " Completed" << endl; cerr << " Completed" << endl;
} }
...@@ -521,7 +679,7 @@ namespace TACO { ...@@ -521,7 +679,7 @@ namespace TACO {
// Calc global Vrow: // Calc global Vrow:
for(int i = 1; i <= acEst.getNumVolumes(); i++) for(int i = 1; i <= acEst.getNumVolumes(); i++)
{ {
ret(i) = MISCMATHS::mean(acEst.getVolume(i).AsColumn()); ret(i) = MISCMATHS::mean(ColumnVector(acEst.getVolume(i).AsColumn()));
} }
} }
} }
......
Makefile
+ 8
− 3
View file @ e50138b4
...@@ -6,11 +6,13 @@ USRINCFLAGS = -I${INCDIR}/newmat -I${INCDIR}/avwio -I${INCDIR}/libprob -I${INCDI ...@@ -6,11 +6,13 @@ USRINCFLAGS = -I${INCDIR}/newmat -I${INCDIR}/avwio -I${INCDIR}/libprob -I${INCDI
LIBS = -lmiscmaths -lm -lnewmat -lavwio -lprob LIBS = -lmiscmaths -lm -lnewmat -lavwio -lprob
XFILES = film_ols film_gls quick contrast_mgr ttoz band_pass XFILES = film_ols film_gls quick contrast_mgr ttoz ftoz band_pass
OBJS = ContrastMgrOptions.o ContrastMgr.o BandPassOptions.o sigproc.o ols.o histogram.o gaussComparer.o AutoCorrEstimator.o Log.o Volume.o VolumeSeries.o glm.o t2z.o paradigm.o glim.o FilmOlsOptions.o glmrand.o FilmGlsOptions.o glimGls.o OBJS = ContrastMgrOptions.o ContrastMgr.o BandPassOptions.o sigproc.o histogram.o gaussComparer.o AutoCorrEstimator.o Log.o Volume.o VolumeSeries.o glm.o t2z.o paradigm.o glim.o FilmOlsOptions.o ols.o FilmGlsOptions.o glimGls.o base2z.o f2z.o
all: film_ols film_gls quick contrast_mgr ttoz band_pass #OPTFLAGS =
all: ${XFILES}
film_ols:${OBJS} film_ols.o film_ols:${OBJS} film_ols.o
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} film_ols.o ${LIBS} ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} film_ols.o ${LIBS}
...@@ -21,6 +23,9 @@ contrast_mgr:${OBJS} contrast_mgr.o ...@@ -21,6 +23,9 @@ contrast_mgr:${OBJS} contrast_mgr.o
ttoz:${OBJS} ttoz.o ttoz:${OBJS} ttoz.o
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} ttoz.o ${LIBS} ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} ttoz.o ${LIBS}
ftoz:${OBJS} ftoz.o
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} ftoz.o ${LIBS}
band_pass:${OBJS} band_pass.o band_pass:${OBJS} band_pass.o
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} band_pass.o ${LIBS} ${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} band_pass.o ${LIBS}
......
film_gls.cc
+ 193
− 71
View file @ e50138b4
...@@ -25,7 +25,7 @@ ...@@ -25,7 +25,7 @@
#include "paradigm.h" #include "paradigm.h"
#include "FilmGlsOptions.h" #include "FilmGlsOptions.h"
#include "glimGls.h" #include "glimGls.h"
#include <vector>
#include <string> #include <string>
#ifndef NO_NAMESPACE #ifndef NO_NAMESPACE
...@@ -94,96 +94,218 @@ int main(int argc, char *argv[]) ...@@ -94,96 +94,218 @@ int main(int argc, char *argv[])
logger.out("Gc", parad.getDesignMatrix()); logger.out("Gc", parad.getDesignMatrix());
} }
// Set up OLS GLM for non-whitened data vector<Matrix> dms;
Glim glim(x, parad.getDesignMatrix()); for(int i = 1; i <= numTS; i++)
cerr << "Computing residuals for non-whitened data... ";
const VolumeSeries& rnotw = glim.ComputeResids();
cerr << "Completed" << endl;
if(globalopts.verbose)
{ {
logger.out("rnotw", rnotw); dms.push_back(parad.getDesignMatrix());
} }
// Estimate Autocorrelations: // Setup OLS GLM for temporally filtered data:
AutoCorrEstimator acEst(rnotw); int numParams = parad.getDesignMatrix().Ncols();
GlimGls glimGls(numTS, sizeTS, numParams);
VolumeSeries residuals(sizeTS, numTS, x.getDims(), x.getPreThresholdPositions());
AutoCorrEstimator acEst(residuals);
if(globalopts.fitAutoRegressiveModel) if(!globalopts.noest)
{
acEst.fitAutoRegressiveModel();
if(globalopts.verbose)
{
AutoCorrEstimator acEstForLogging(rnotw);
acEstForLogging.calcRaw();
logger.out("rawac", acEstForLogging.getEstimates());
logger.out("autoregac", acEst.getEstimates());
}
}
else
{ {
acEst.calcRaw(); int numiters = globalopts.numiters+1;
int iters = 1;
// iters==1 is for high freq removal
if(!globalopts.highfreqremoval)
iters++;
for(; iters <= numiters; iters++)
{
cerr << "iters = " << iters << endl;
// Loop through voxels:
cerr << "Calculating residuals..." << endl;
for(int i = 1; i <= numTS; i++)
{
glimGls.setData(x.getSeries(i), dms[i-1], i);
residuals.getSeries(i) = glimGls.getResiduals();
}
cerr << "Completed" << endl;
cerr << "Estimating residual autocorrelation..." << endl;
// Estimate Autocorrelation:
acEst.calcRaw();
if(iters==1)
{
// iters==1 is for high freq removal
acEst.tukey(10);
}
else
{
if(globalopts.fitAutoRegressiveModel)
{
acEst.fitAutoRegressiveModel();
}
else if(globalopts.tukey)
{
// Smooth raw estimates:
if(globalopts.smoothACEst)
{
acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath, globalopts.epith);
}
if(globalopts.tukeysize == 0)
globalopts.tukeysize = (int)(2*sqrt(sizeTS))/2;
acEst.tukey(globalopts.tukeysize);
}
else if(globalopts.multitaper)
{
acEst.multitaper(globalopts.multitapersize);
}
else if(globalopts.pava)
{
// Smooth raw estimates:
if(globalopts.smoothACEst)
{
acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath, globalopts.epith);
}
acEst.pava();
}
if(globalopts.verbose) }
{ cerr << "Completed" << endl;
logger.out("rawac", acEst.getEstimates());
} // Loop through voxels:
cerr << "Prewhitening..." << endl;
int co = 1;
// Smooth raw estimates: for(int i = 1; i <= numTS; i++)
if(globalopts.smoothACEst) {
{ ColumnVector I(sizeTS);
acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath); I = 0;
I(1) = 1;
} ColumnVector xw(sizeTS);
ColumnVector xprew(sizeTS);
acEst.setDesignMatrix(dms[i-1]);
// Use autocorr estimate to prewhiten data:
xprew = x.getSeries(i);
Matrix designmattw;
// iters==1 is for high freq removel
acEst.preWhiten(xprew, xw, i, designmattw, bool(iters==1));
if(co > 1000)
{
co = 1;
cerr << (float)i/(float)numTS << ",";
}
else
co++;
x.getSeries(i) = xw;
dms[i-1] = designmattw;
}
cerr << "Completed" << endl;
// Apply constraints to estimate autocorr: // Add param number to "pe" to create filename:
acEst.pava(); char strc[4];
ostrstream osc(strc,4);
osc << iters - 1 << '\0';
VolumeSeries& threshac = acEst.getEstimates();
int cutoff = sizeTS/2;
if(globalopts.tukey)
cutoff = globalopts.tukeysize;
threshac = threshac.Rows(1,cutoff);
VolumeSeries::Dims dims = x.getDims();
dims.v = cutoff;
threshac.unthresholdSeries(dims,x.getPreThresholdPositions());
threshac.writeAsFloat(logger.getDir() + "/threshac" + strc);
threshac.thresholdSeries();
if(globalopts.verbose)
{
cerr << "Saving results... " << endl;
if(globalopts.verbose) ColumnVector& countLargeE = acEst.getCountLargeE();
{
logger.out("threshac", acEst.getEstimates()); logger.out(string("countLargeE") + strc, countLargeE);
residuals.unthresholdSeries(x.getDims(),x.getPreThresholdPositions());
residuals.writeAsFloat(logger.getDir() + "/res4d" + strc);
residuals.thresholdSeries();
cerr << "Completed" << endl;
}
} }
} }
else // no estimation of autocorrelations
ColumnVector I(sizeTS); {
I = 0; // do nothing
I(1) = 1; }
// Setup OLS GLM for temporally filtered data:
int numParams = parad.getDesignMatrix().Ncols();
GlimGls glimGls(numTS, sizeTS, numParams);
ColumnVector xw(sizeTS);
ColumnVector xprew(sizeTS);
acEst.setDesignMatrix(parad.getDesignMatrix());
int co = 1;
// Loop through voxels:
cerr << "Calculating params..." << endl;
// Do once more to compute real param ests:
cerr << "Computing parameter estimates..." << endl;
for(int i = 1; i <= numTS; i++) for(int i = 1; i <= numTS; i++)
{ {
// Use autocorr estimate to prewhiten data: glimGls.setData(x.getSeries(i), dms[i-1], i);
xprew = x.getSeries(i);
if(globalopts.verbose)
Matrix designmattw;
acEst.preWhiten(xprew, xw, i, designmattw);
if(co > 1000 || i == 5618 || i == 5582)
{ {
co = 1; residuals.getSeries(i) = glimGls.getResiduals();
cerr << (float)i/(float)numTS << ",";
} }
else
co++;
glimGls.setData(xw, designmattw, i);
} }
cerr << "Completed" << endl; cerr << "Completed" << endl;
// Write out necessary data: // Write out necessary data:
cerr << "Saving results... "; cerr << "Saving results... " << endl;
residuals.unthresholdSeries(x.getDims(),x.getPreThresholdPositions());
residuals.writeAsFloat(logger.getDir() + "/res4d");
// write out design matrices - a volume Series for each param
if(globalopts.verbose_dms)
{
VolumeSeries::Dims dims = x.getDims();
for(int j = 1; j <= numParams; j++)
{
char strc[4];
ostrstream osc(strc,4);
osc << j << '\0';
VolumeSeries dmsmat(sizeTS, numTS);
for(int i = 1; i <= numTS; i++)
{
dmsmat.getSeries(i) = dms[i-1].Column(j);
}
dmsmat.setDims(dims);
dmsmat.setPreThresholdPositions(x.getPreThresholdPositions());
dmsmat.unthresholdSeries();
dmsmat.writeAsFloat(logger.getDir() + "/dms" + strc);
}
// output x
x.unthresholdSeries();
x.writeAsFloat(logger.getDir() + "/wx");
}
if(globalopts.verbose)
{
VolumeSeries::Dims dims = x.getDims();
// Save E
VolumeSeries& E = acEst.getE();
dims.v = acEst.getZeroPad();
E.setDims(dims);
E.setPreThresholdPositions(x.getPreThresholdPositions());
E.unthresholdSeries();
E.writeAsFloat(logger.getDir() + "/E");
}
glimGls.Save(x.getDims(), x.getPreThresholdPositions()); glimGls.Save(x.getDims(), x.getPreThresholdPositions());
cerr << "Completed" << endl; cerr << "Completed" << endl;
} }
......
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