diff --git a/AutoCorrEstimator.cc b/AutoCorrEstimator.cc
index f7fdf927ef22e9b0bfa9b48e8a6696ba2ed83e79..972c451829cfcbc9d9e1057846843adf03ac89d7 100644
--- a/AutoCorrEstimator.cc
+++ b/AutoCorrEstimator.cc
@@ -8,12 +8,14 @@
#include <iostream>
#include <strstream>
+#define WANT_STREAM
#include "AutoCorrEstimator.h"
#include "miscmaths.h"
#include "Log.h"
#include "Volume.h"
#include "histogram.h"
+#include "miscmaths.h"
using namespace NEWMAT;
using namespace UTILS;
@@ -40,7 +42,7 @@ namespace TACO {
{
dummy = 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);
FFT(dmrow, dummy, dm_fft_real, dm_fft_imag);
@@ -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");
int sizeTS = xdata.getNumVolumes();
@@ -66,6 +68,25 @@ namespace TACO {
vrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = acEst.getSeries(i).Rows(2, sizeTS/2).Reverse();
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
dummy = 0;
@@ -74,14 +95,25 @@ namespace TACO {
FFT(xrow, dummy, x_fft_real, x_fft_im);
- // 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++)
+ if(highfreqremovalonly)
{
- 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
for(int k = 1; k <= numPars; k++)
{
@@ -90,12 +122,12 @@ namespace TACO {
FFTI(SP(ac_fft_real, dm_fft_real), SP(ac_fft_real, dm_fft_imag), realifft, dummy);
// place result into ret:
- dmret.Column(k) = realifft.Rows(1,sizeTS);
- float std = pow(MISCMATHS::var(dmret.Column(k)),0.5);
- dmret.Column(k) = (dmret.Column(k)/std) + mn(k);
+ dmret.Column(k) = realifft.Rows(1,sizeTS) + mn(k);
+ //float std = pow(MISCMATHS::var(ColumnVector(dmret.Column(k))),0.5);
+ //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);
// place result into ret:
@@ -133,7 +165,7 @@ namespace TACO {
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();
FFT(vrow, dummy, ac_fft_real, ac_fft_im);
-
+
// FFT x data
dummy = 0;
xrow = 0;
@@ -144,11 +176,15 @@ namespace TACO {
// 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/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:
FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy);
@@ -223,6 +259,7 @@ namespace TACO {
}
Log::getInstance().out("order", order);
+ countLargeE = 0;
cerr << " Completed" << endl;
}
@@ -235,7 +272,7 @@ namespace TACO {
hist.generate();
float mode = hist.mode();
cerr << "mode = " << mode << endl;
-
+
float sum = 0.0;
int count = 0;
@@ -255,14 +292,17 @@ namespace TACO {
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");
Log& logger = Log::getInstance();
- // Establish epi thresh to use:
- int usanthresh = establishUsanThresh(epivol);
-
+ if(usanthresh == 0)
+ {
+ // Establish epi thresh to use:
+ usanthresh = establishUsanThresh(epivol);
+ }
+
// Setup external call to susan program:
char callsusanstr[1000];
ostrstream osc3(callsusanstr,1000);
@@ -288,7 +328,7 @@ namespace TACO {
cerr << callsusanstr << endl;
- for(; i < 20; i++)
+ for(; i < 30; i++)
{
// output unsmoothed estimates:
vol = acEst.getVolume(i).AsColumn()*factor;
@@ -310,10 +350,9 @@ namespace TACO {
char rmfilesstr[1000];
ostrstream osc(rmfilesstr,1000);
osc << "rm -rf "
- << logger.getDir() + "/" + postSmoothVol + "* "
- << logger.getDir() + "/" + preSmoothVol + "* "
- << logger.getDir() + "/" + epifname + "* "
- << logger.getDir() + "/usanSize*" << '\0';
+ << logger.getDir() + "/" + postSmoothVol + "* "
+ << logger.getDir() + "/" + preSmoothVol + "* "
+ << logger.getDir() + "/usanSize*" << '\0';
cerr << rmfilesstr << endl;
@@ -326,7 +365,7 @@ namespace TACO {
cerr << "Calculating raw AutoCorrs...";
- AutoCorr(xdata, acEst, zeropad);
+ AutoCorr(xdata, acEst, E, zeropad);
cerr << " Completed" << endl;
}
@@ -369,7 +408,124 @@ namespace TACO {
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() {
Tracer tr("AutoCorrEstimator::pava");
@@ -437,7 +593,7 @@ namespace TACO {
else if(j > 2)
//if(j > 2)
- {
+ {
int endst = j;
int stst = j-(int)(1+(j/8.0));
@@ -446,11 +602,13 @@ namespace TACO {
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;
}
@@ -521,7 +679,7 @@ namespace TACO {
// Calc global Vrow:
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()));
}
}
}
diff --git a/Makefile b/Makefile
index 004fad0da88a1d7074864194c1af2c9a8c79195e..b6b492fc9bf1abc81a9228c4752fd7fd16e1eb78 100644
--- a/Makefile
+++ b/Makefile
@@ -6,11 +6,13 @@ USRINCFLAGS = -I${INCDIR}/newmat -I${INCDIR}/avwio -I${INCDIR}/libprob -I${INCDI
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
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} film_ols.o ${LIBS}
@@ -21,6 +23,9 @@ contrast_mgr:${OBJS} contrast_mgr.o
ttoz:${OBJS} ttoz.o
${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
${CXX} ${CXXFLAGS} ${LDFLAGS} -o $@ ${OBJS} band_pass.o ${LIBS}
diff --git a/film_gls.cc b/film_gls.cc
index d1b675635413ebc514999de772160831b47f0143..15c13ad14f1bf104471a1b215ad973527a95da0a 100644
--- a/film_gls.cc
+++ b/film_gls.cc
@@ -25,7 +25,7 @@
#include "paradigm.h"
#include "FilmGlsOptions.h"
#include "glimGls.h"
-
+#include <vector>
#include <string>
#ifndef NO_NAMESPACE
@@ -94,96 +94,218 @@ int main(int argc, char *argv[])
logger.out("Gc", parad.getDesignMatrix());
}
- // Set up OLS GLM for non-whitened data
- Glim glim(x, parad.getDesignMatrix());
-
- cerr << "Computing residuals for non-whitened data... ";
- const VolumeSeries& rnotw = glim.ComputeResids();
- cerr << "Completed" << endl;
-
- if(globalopts.verbose)
+ vector<Matrix> dms;
+ for(int i = 1; i <= numTS; i++)
{
- logger.out("rnotw", rnotw);
+ dms.push_back(parad.getDesignMatrix());
}
- // Estimate Autocorrelations:
- AutoCorrEstimator acEst(rnotw);
+ // Setup OLS GLM for temporally filtered data:
+ int numParams = parad.getDesignMatrix().Ncols();
+ GlimGls glimGls(numTS, sizeTS, numParams);
+
+ VolumeSeries residuals(sizeTS, numTS, x.getDims(), x.getPreThresholdPositions());
+ AutoCorrEstimator acEst(residuals);
- if(globalopts.fitAutoRegressiveModel)
- {
- acEst.fitAutoRegressiveModel();
- if(globalopts.verbose)
- {
- AutoCorrEstimator acEstForLogging(rnotw);
- acEstForLogging.calcRaw();
- logger.out("rawac", acEstForLogging.getEstimates());
- logger.out("autoregac", acEst.getEstimates());
- }
- }
- else
+ if(!globalopts.noest)
{
- 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)
- {
- logger.out("rawac", acEst.getEstimates());
- }
+ }
+ cerr << "Completed" << endl;
+
+ // Loop through voxels:
+ cerr << "Prewhitening..." << endl;
+ int co = 1;
- // Smooth raw estimates:
- if(globalopts.smoothACEst)
- {
- acEst.spatiallySmooth(logger.getDir() + "/" + globalopts.epifname, epivol, globalopts.ms, globalopts.epifname, globalopts.susanpath);
+ for(int i = 1; i <= numTS; i++)
+ {
+ ColumnVector I(sizeTS);
+ 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:
- acEst.pava();
+ // Add param number to "pe" to create filename:
+ 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)
- {
- logger.out("threshac", acEst.getEstimates());
+ ColumnVector& countLargeE = acEst.getCountLargeE();
+
+ logger.out(string("countLargeE") + strc, countLargeE);
+
+ residuals.unthresholdSeries(x.getDims(),x.getPreThresholdPositions());
+ residuals.writeAsFloat(logger.getDir() + "/res4d" + strc);
+ residuals.thresholdSeries();
+ cerr << "Completed" << endl;
+ }
}
}
-
- ColumnVector I(sizeTS);
- I = 0;
- 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;
+ else // no estimation of autocorrelations
+ {
+ // do nothing
+ }
+ // Do once more to compute real param ests:
+ cerr << "Computing parameter estimates..." << endl;
for(int i = 1; i <= numTS; i++)
- {
- // Use autocorr estimate to prewhiten data:
- xprew = x.getSeries(i);
-
- Matrix designmattw;
- acEst.preWhiten(xprew, xw, i, designmattw);
-
- if(co > 1000 || i == 5618 || i == 5582)
+ {
+ glimGls.setData(x.getSeries(i), dms[i-1], i);
+
+ if(globalopts.verbose)
{
- co = 1;
- cerr << (float)i/(float)numTS << ",";
+ residuals.getSeries(i) = glimGls.getResiduals();
}
- else
- co++;
-
- glimGls.setData(xw, designmattw, i);
}
cerr << "Completed" << endl;
// 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());
cerr << "Completed" << endl;
}