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

removed junk files

parent 809d3868
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
AutoCorrEstimator.bak2 deleted 100755 → 0
+ 0
− 536
View file @ 809d3868
#include <iostream>
#include <strstream>
// Written by Mark Woolrich 6/1/99
#include "AutoCorrEstimator.h"
#include "miscmaths.h"
#include "Log.h"
#include "Volume.h"
#include "histogram.h"
using namespace NEWMAT;
using namespace UTILS;
namespace TACO {
void AutoCorrEstimator::setDesignMatrix(const Matrix& dm) {
Tracer tr("AutoCorrEstimator::setDesignMatrix");
int sizeTS = xdata.getNumVolumes();
int numPars = dm.Ncols();
dminFFTReal.ReSize(zeropad, numPars);
dminFFTImag.ReSize(zeropad, numPars);
ColumnVector dmrow;
dmrow.ReSize(zeropad);
ColumnVector dm_fft_real, dm_fft_imag;
ColumnVector dummy(zeropad);
ColumnVector realifft(zeropad);
// FFT design matrix
for(int k = 1; k <= numPars; k++)
{
dummy = 0;
dmrow = 0;
mn(k) = MISCMATHS::mean(dm.Column(k));
dmrow.Rows(1,sizeTS) = dm.Column(k) - mn(k);
FFT(dmrow, dummy, dm_fft_real, dm_fft_imag);
dminFFTImag.Column(k) = dm_fft_imag;
dminFFTReal.Column(k) = dm_fft_real;
}
}
void AutoCorrEstimator::preWhiten(ColumnVector& in, ColumnVector& ret, int i, Matrix& dmret) {
Tracer tr("AutoCorrEstimator::preWhiten");
int sizeTS = xdata.getNumVolumes();
int numPars = dminFFTReal.getNumSeries();
ret.ReSize(sizeTS);
dmret.ReSize(sizeTS, numPars);
// FFT auto corr estimate
dummy = 0;
vrow = 0;
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;
xrow.Rows(1,sizeTS) = in;
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++)
{
ac_fft_real(j) = 1.0/ac_fft_real(j);
}
// filter design matrix
for(int k = 1; k <= numPars; k++)
{
dm_fft_real = dminFFTReal.getSeries(k);
dm_fft_imag = dminFFTImag.getSeries(k);
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);
}
// Do filtering and inverse FFT:
FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy);
// place result into ret:
ret = realifft.Rows(1,sizeTS);
}
void AutoCorrEstimator::preWhiten(VolumeSeries& in, VolumeSeries& ret)
{
Tracer tr("AutoCorrEstimator::preWhiten");
cerr << "Prewhitening... ";
int sizeTS = xdata.getNumVolumes();
int numTS = xdata.getNumSeries();
ret.ReSize(sizeTS, numTS);
// make sure p_vrow is cyclic (even function)
ColumnVector vrow, xrow;
vrow.ReSize(zeropad);
xrow.ReSize(zeropad);
ColumnVector x_fft_real, ac_fft_real;
ColumnVector x_fft_im, ac_fft_im;
ColumnVector dummy(zeropad);
ColumnVector realifft(zeropad);
int co = 1;
for(int i = 1; i <= numTS; i++)
{
// FFT auto corr estimate
dummy = 0;
vrow = 0;
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;
xrow.Rows(1,sizeTS/2) = in.getSeries(i).Rows(1,sizeTS/2);
xrow.Rows(zeropad - sizeTS/2 + 2, zeropad) = in.getSeries(i).Rows(2, sizeTS/2).Reverse();
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++)
{
ac_fft_real(j) = 1.0/ac_fft_real(j);
}
// Do filtering and inverse FFT:
FFTI(SP(ac_fft_real, x_fft_real), SP(ac_fft_real, x_fft_im), realifft, dummy);
// place result into ret:
ret.Column(i) = realifft.Rows(1,sizeTS);
if(co > 100)
{
co = 1;
cerr << (float)i/(float)numTS << ",";
}
else
co++;
}
cerr << " Completed" << endl;
}
void AutoCorrEstimator::fitAutoRegressiveModel()
{
Tracer trace("AutoCorrEstimator::fitAutoRegressiveModel");
cerr << "Fitting autoregressive model..." << endl;
const int maxorder = 10;
const int minorder = 1;
int sizeTS = xdata.getNumVolumes();
int numTS = xdata.getNumSeries();
// setup temp variables
ColumnVector x(sizeTS);
ColumnVector order(numTS);
ColumnVector betas(maxorder);
acEst.ReSize(sizeTS, numTS);
acEst = 0;
int co = 1;
for(int i = 1; i <= numTS; i++)
{
x = xdata.getSeries(i).AsColumn();
order(i) = SIGPROC::Pacf(x, minorder, maxorder, betas);
if(order(i) != -1)
{
// Calculate auto corr:
ColumnVector Krow(sizeTS);
Krow = 0;
Krow(sizeTS) = 1;
Krow.Rows(sizeTS-order(i), sizeTS-1) = -betas.Reverse();
Matrix Kinv(sizeTS, sizeTS);
Kinv = 0;
for(int j = 1; j <= sizeTS; j++)
{
Kinv.SubMatrix(j,j,1,j) = Krow.Rows(sizeTS-j+1,sizeTS).t();
}
// Kinv now becomes V:
Kinv = (Kinv.t()*Kinv).i();
acEst.SubMatrix(1,sizeTS/2+1,i,i) = (Kinv.SubMatrix(sizeTS/2, sizeTS/2, sizeTS/2, sizeTS)/Kinv.MaximumAbsoluteValue()).AsColumn();
if(co > 200)
{
co = 1;
cerr << (float)i/(float)numTS << ",";
}
else
co++;
}
}
Log::getInstance().out("order", order);
cerr << " Completed" << endl;
}
int AutoCorrEstimator::establishUsanThresh(const Volume& epivol)
{
int usanthresh = 100;
int num = epivol.getVolumeSize();
Histogram hist(epivol, num/200);
hist.generate();
float mode = hist.mode();
cerr << "mode = " << mode << endl;
float sum = 0.0;
int count = 0;
// Work out standard deviation from mode for values greater than mode:
for(int i = 1; i <= num; i++) {
if(epivol(i) > mode) {
sum = sum + (epivol(i) - mode)*(epivol(i) - mode);
count++;
}
}
int sig = (int)pow(sum/num, 0.5);
cerr << "sig = " << sig << endl;
usanthresh = sig/3;
return usanthresh;
}
void AutoCorrEstimator::spatiallySmooth(const string& usanfname, const Volume& epivol, int masksize, const string& epifname, const string& susanpath) {
Tracer trace("AutoCorrEstimator::spatiallySmooth");
Log& logger = Log::getInstance();
// Establish epi thresh to use:
int usanthresh = establishUsanThresh(epivol);
// Setup external call to susan program:
char callsusanstr[1000];
ostrstream osc3(callsusanstr,1000);
string preSmoothVol = "preSmoothVol";
string postSmoothVol = "postSmoothVol";
osc3 << susanpath << " "
<< logger.getDir() << "/" << preSmoothVol << " 1 "
<< logger.getDir() << "/" << postSmoothVol << " "
<< masksize << " 3D 0 1 " << usanfname << " " << usanthresh << " "
<< logger.getDir() << "/" << "usanSize" << '\0';
// Loop through first third of volumes
// assume the rest are zero
int factor = 10000;
// Setup volume for reading and writing volumes:
Volume vol(acEst.getNumSeries());
vol.setDims(xdata.getDims());
vol.setPreThresholdPositions(xdata.getPreThresholdPositions());
acEst.setDims(xdata.getDims());
acEst.setPreThresholdPositions(xdata.getPreThresholdPositions());
int i = 2;
acEst.unthresholdSeries();
cerr << "Spatially smoothing auto corr estimates" << endl;
cerr << callsusanstr << endl;
for(; i < 20; i++)
{
// output unsmoothed estimates:
vol = acEst.getVolume(i).AsColumn()*factor;
vol.writeAsInt(logger.getDir() + "/" + preSmoothVol);
// call susan:
system(callsusanstr);
// read in smoothed volume:
vol.read(logger.getDir() + "/" + postSmoothVol);
acEst.getVolume(i) = vol.AsRow()/factor;
cerr << ".";
}
cerr << endl;
// Clear unwanted written files
char rmfilesstr[1000];
ostrstream osc(rmfilesstr,1000);
osc << "rm -rf "
<< logger.getDir() + "/" + postSmoothVol + "* "
<< logger.getDir() + "/" + preSmoothVol + "* "
<< logger.getDir() + "/" + epifname + "* "
<< logger.getDir() + "/usanSize*" << '\0';
cerr << rmfilesstr << endl;
system(rmfilesstr);
cerr << "Completed" << endl;
acEst.thresholdSeries();
}
void AutoCorrEstimator::calcRaw() {
cerr << "Calculating raw AutoCorrs...";
Matrix fts;
AutoCorr(xdata.t(), acEst, fts, zeropad);
acEst = acEst.t();
cerr << " Completed" << endl;
}
void AutoCorrEstimator::filter(const ColumnVector& filterFFT) {
Tracer tr("AutoCorrEstimator::filter");
cerr << "Combining temporal filtering effects with AutoCorr estimates... ";
// This function adjusts the autocorrelations as if the
// xdata has been filtered by the passed in filterFFT
// DOES NOT filter the xdata itself
ColumnVector vrow;
// make sure p_vrow is cyclic (even function)
vrow.ReSize(zeropad);
ColumnVector fft_real;
ColumnVector fft_im;
ColumnVector dummy(zeropad);
ColumnVector realifft(zeropad);
int sizeTS = xdata.getNumVolumes();
for(int i = 1; i <= xdata.getNumSeries(); i++)
{
dummy = 0;
vrow = 0;
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, fft_real, fft_im);
FFTI(SP(fft_real, filterFFT), dummy, realifft, dummy);
// place result into acEst:
acEst.Column(i) = realifft.Rows(1,sizeTS)/realifft(1);
}
cerr << " Completed" << endl;
}
void AutoCorrEstimator::pava() {
Tracer tr("AutoCorrEstimator::pava");
cerr << "Using New PAVA on AutoCorr estimates... ";
for(int i = 1; i <= xdata.getNumSeries(); i++) {
int sizeTS = xdata.getNumVolumes();
int stopat = (int)sizeTS/2;
// 5% point of distribution of autocorr about zero
const float th = (-1/sizeTS)+(2/sqrt(sizeTS));
ColumnVector values = acEst.Column(i);
ColumnVector zero(1);
zero = 0;
values = values.Rows(1,stopat) & zero;
ColumnVector gm(stopat + 1);
for(int j = 1; j <= stopat + 1; gm(j) = j++);
ColumnVector weights(stopat+1);
weights = 1;
bool anyviolators = true;
while(anyviolators) {
anyviolators = false;
for(int k = 2; k <= values.Nrows(); k++) {
if(values(k) > values(k-1)) {
anyviolators = true;
values(k-1) = (values(k-1)*weights(k-1) + values(k)*weights(k))/(weights(k-1) + weights(k));
values = values.Rows(1,k-1) & values.Rows(k+1,values.Nrows());
weights(k-1) = weights(k) + weights(k-1);
weights = weights.Rows(1,k-1) & weights.Rows(k+1,weights.Nrows());
for(int j = 1; j <= stopat + 1; j++) {
if(gm(j) >= k)
gm(j) = gm(j)-1;
}
break;
}
}
}
acEst.Column(i) = 0.0;
int j=1;
for(; j <= stopat; j++) {
acEst(j,i) = values(gm(j));
if(acEst(j,i) <= 0.0)
{
acEst(j,i) = 0.0;
break;
}
}
if(acEst(2,i) < th/2)
{
acEst.SubMatrix(3,stopat,i,i) = 0;
}
else if(j > 2)
{
int endst = j;
int stst = j-(int)(1+(j/8.0));
const int expwidth = MISCMATHS::Max((endst - stst)/2,1);
const int exppow = 2;
for(j = stst; j <= endst; j++)
{
acEst(j,i) = acEst(j,i)*exp(-pow((j-stst)/expwidth,exppow));
}
}
}
cerr << " Completed" << endl;
}
void AutoCorrEstimator::applyConstraints() {
Tracer tr("AutoCorrEstimator::applyConstraints");
cerr << "Applying constraints to AutoCorr estimates... ";
for(int i = 1; i <= xdata.getNumSeries(); i++)
{
int sizeTS = xdata.getNumVolumes();
int j = 3;
int stopat = (int)sizeTS/4;
// found1 is last valid value above threshold
int found1 = stopat;
// 5% point of distribution of autocorr about zero
const float thresh = (-1/sizeTS)+(2/sqrt(sizeTS));
acEst(2,i) = (acEst(2,i)+ acEst(3,i))/2;
if(acEst(2,i) < 0)
{
acEst(2,i) = 0;
}
else
{
float grad = 0.0;
while(j <= stopat && j < found1 + 2)
{
grad = ((acEst(j,i) + acEst(j+1,i))/2 - acEst(j-1,i))/1.5;
if(grad < 0)
acEst(j,i) = grad + acEst(j-1,i);
else
acEst(j,i) = acEst(j-1,i);
// look for threshold
if(acEst(j,i) < thresh/3.0 && found1 == stopat)
{
found1 = j;
}
if(acEst(j,i) < 0)
{
acEst(j,i) = 0;
}
j++;
}
}
// set rest to zero:
for(; j <= sizeTS; j++)
{
acEst(j,i) = 0;
}
}
cerr << "Completed" << endl;
}
void AutoCorrEstimator::getMeanEstimate(ColumnVector& ret)
{
Tracer tr("AutoCorrEstimator::getMeanEstimate");
ret.ReSize(acEst.getNumVolumes());
// Calc global Vrow:
for(int i = 1; i <= acEst.getNumVolumes(); i++)
{
ret(i) = MISCMATHS::mean(acEst.getVolume(i).AsColumn());
}
}
}
Logbak.cc deleted 100644 → 0
+ 0
− 53
View file @ 809d3868
/* Log.cc
Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#include "Log.h"
using namespace UTILS;
namespace UTILS {
Log* Log::logger = NULL;
void Log::establishDir(const string& name)
{
dir = name;
// make directory to place results into:
// keep adding "+" until directory is made:
while(true)
{
int ret = system(("mkdir "+ dir).c_str());
if(ret == 0)
{
break;
}
dir = dir + "+";
}
cerr << "Results are in directory " + dir << endl;
// setup logfile
logfileout.open((dir + "/" + logfilename).c_str(), ios::out);
}
void Log::setDir(const string& name)
{
dir = name;
cerr << "Results are in directory " + dir << endl;
// setup logfile
logfileout.open((dir + "/" + logfilename).c_str(), ios::out);
}
}
Logbak.h deleted 100644 → 0
+ 0
− 111
View file @ 809d3868
/* Log.h
Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <math.h>
#include "newmatap.h"
#include "newmatio.h"
using namespace NEWMAT;
namespace UTILS{
#if !defined(__Log_h)
#define __Log_h
class Log
{
public:
static Log& getInstance();
~Log() { delete logger; }
void establishDir(const string& name);
void setDir(const string& name);
void setLogFile(const string& name) {logfilename = name;}
const string& getDir() const { return dir; }
void out(const string& p_fname, const Matrix& p_mat, bool p_addMatrixString = true);
void out(const string& p_fname, const RowVector& p_mat);
void out(const string& p_fname, const ColumnVector& p_mat);
ofstream& str() { return logfileout; }
private:
Log() {}
const Log& operator=(Log&);
Log(Log&);
static Log* logger;
string dir;
ofstream logfileout;
string logfilename;
};
inline void Log::out(const string& p_fname, const Matrix& p_mat, bool p_addMatrixString)
{
ofstream out;
out.open((dir + "/" + p_fname).c_str(), ios::out);
out.setf(ios::scientific, ios::floatfield);
if(p_addMatrixString)
out << "/Matrix" << endl;
for(int i=1; i<=p_mat.Nrows(); i++)
{
for(int j=1; j<=p_mat.Ncols(); j++)
{
out << p_mat(i,j) << " ";
}
out << endl;
}
out.close();
}
inline void Log::out(const string& p_fname, const ColumnVector& p_mat)
{
ofstream out;
out.open((dir + "/" + p_fname).c_str(), ios::out);
out.setf(ios::scientific, ios::floatfield);
for(int j=1; j<=p_mat.Nrows(); j++)
{
out << p_mat(j);
out << endl;
}
out.close();
}
inline void Log::out(const string& p_fname, const RowVector& p_mat)
{
ofstream out;
out.open((dir + "/" + p_fname).c_str(), ios::out);
out.setf(ios::scientific, ios::floatfield);
for(int j=1; j<=p_mat.Ncols(); j++)
{
out << p_mat(j) << " ";
}
out << endl;
out.close();
}
inline Log& Log::getInstance(){
if(logger == NULL)
logger = new Log();
return *logger;
}
#endif
}
deb deleted 100755 → 0
+ 0
− 1
View file @ 809d3868
1.2246468e-16
logfile deleted 100755 → 0
+ 0
− 0
View file @ 809d3868
tmp deleted 100755 → 0
+ 0
− 2
View file @ 809d3868
x = sin(pi)
save -ASCII deb x
\ No newline at end of file
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