Skip to content
Snippets Groups Projects
Commit 77af356f authored by David Flitney's avatar David Flitney
Browse files

Removed Volume and VolumeSeries from here - using the one in miscmaths

parent 830f1051
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Volume.cc deleted 100644 → 0
+ 0
− 253
View file @ 830f1051
/* Volume.cc
Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#include <iostream>
#include <cstdlib>
#include "avwio/avwio.h"
#include "newmatap.h"
#include "newmatio.h"
#include "Volume.h"
using namespace NEWMAT;
namespace FILM {
void Volume::unthreshold()
{
Tracer ts("Volume::unthreshold");
int sizeVolUnthresholded = getUnthresholdSize();
int sizeVolThresholded = getVolumeSize();
// Increase this to required size and set to 0
Release();
ColumnVector X=*this;
ReSize(sizeVolUnthresholded);
ColumnVector::operator=(0);
// Loop through volume restoring non thresholded values:
for(int i = 1; i <= sizeVolThresholded; i++)
{
(*this)(int(preThresholdPositions(i)))=X(i);
}
}
void Volume::unthreshold(const VolumeSeries::Dims& pdims,const ColumnVector& in)
{
dims = pdims;
preThresholdPositions = in;
unthreshold();
}
void Volume::threshold()
{
Tracer ts("Volume::threshold");
int sizeVol = preThresholdPositions.Nrows();
ColumnVector X(sizeVol);
// Loop through pulling out non thresholded values
for(int i = 1; i <= sizeVol; i++)
{
X(i)=(*this)(int(preThresholdPositions(i)));
}
(*this) = X;
}
void Volume::threshold(float thresh)
{
Tracer ts("Volume::threshold");
int size = getVolumeSize();
int j = 0;
preThresholdPositions.ReSize(size);
float m = 0;
for(int i = 1; i <= size; i++)
{
m = (*this)(i);
if(m > thresh)
{
j++;
preThresholdPositions(j) = i;
(*this)(j) = (*this)(i);
}
}
// Discard rest:
*this = Rows(1,j);
preThresholdPositions = preThresholdPositions.Rows(1,j);
}
void Volume::writeAsFloat(const string& fname)
{
Tracer ts("Volume::writeAsFloat");
float fmin, fmax;
const ColumnVector& outputvol = *this;
AVW* OP = AvwOpen(fname.c_str(), "wc");
AvwSetDim(OP,dims.x, dims.y, dims.z, 1);
AvwSetVoxDim(OP,dims.vx, dims.vy, dims.vz, 0);
AvwSetDataType(OP, DT_FLOAT);
int sizeVol = outputvol.Nrows();
float *qv = new float[sizeVol];
fmin = outputvol(1);
fmax = outputvol(1);
for (int i=1; i<=sizeVol; i++)
{
if(outputvol(i) < fmin)
fmin = outputvol(i);
else if(outputvol(i) > fmax)
fmax = outputvol(i);
qv[i-1] = outputvol(i);
}
AvwSetMinMax(OP, (short)fmin, (short)fmax);
// fwrite(qv,sizeVol*sizeof(float),1,OP->imgfp);
//fwrite(&OP->header,sizeof(OP->header),1,OP->hdrfp);
AvwWriteVolumes(OP, qv, 1);
delete [] qv;
AvwClose(OP);
}
void Volume::writeAsInt(const string& fname)
{
Tracer ts("Volume::writeAsInt");
int fmin, fmax;
const ColumnVector& outputvol = *this;
AVW* OP = AvwOpen(fname.c_str(), "wc");
AvwSetDim(OP, dims.x, dims.y, dims.z, 1);
AvwSetVoxDim(OP, dims.vx, dims.vy, dims.vz, 0);
AvwSetDataType(OP, DT_SIGNED_SHORT);
int sizeVol = outputvol.Nrows();
short *qv = new short[sizeVol];
fmin = (int)outputvol(1);
fmax = (int)outputvol(1);
for (int i=1; i<=sizeVol; i++)
{
if(outputvol(i) < fmin)
fmin = (int)outputvol(i);
else if(outputvol(i) > fmax)
fmax = (int)outputvol(i);
qv[i-1] = (short)outputvol(i);
}
AvwSetMinMax(OP, (short)fmin, (short)(fmax+0.9999));
AvwWriteVolumes(OP, qv, 1);
delete [] qv;
AvwClose(OP);
}
void Volume::read(const string& fname)
{
Tracer ts("Volume::read");
AVW* IP = AvwOpen(fname.c_str(), "r");
ColumnVector& output = *this;
short x,y,z,v,type;
float vx,vy,vz,tr;
AvwGetDim(IP,&x,&y,&z,&v);
AvwGetVoxDim(IP,&vx,&vy,&vz,&tr);
dims.x = x; dims.y = y; dims.z = z; dims.v = v;
dims.vx = vx; dims.vy = vy; dims.vz = vz; dims.tr = tr;
size_t imagesize=x*y*z;
AvwGetDataType(IP,&type);
output.ReSize(x*y*z);
switch(type)
{
case DT_SIGNED_SHORT:
{
short* sbuffer=new short[imagesize];
AvwReadVolumes(IP,sbuffer,v);
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
output(j)=sbuffer[j-1];
}
delete[] sbuffer;
}
break;
case DT_FLOAT:
{
float* fbuffer=new float[imagesize];
AvwReadVolumes(IP,fbuffer,v);
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
output(j)=fbuffer[j-1];
}
delete[] fbuffer;
}
break;
case DT_UNSIGNED_CHAR:
{
unsigned char* cbuffer=new unsigned char[imagesize];
AvwReadVolumes(IP,cbuffer,v);
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
output(j)=cbuffer[j-1];
}
delete[] cbuffer;
}
break;
default:
perror("AvwRead: DT not supported");
}
AvwClose(IP);
return;
}
}
Volume.h deleted 100644 → 0
+ 0
− 88
View file @ 830f1051
/* Volume.h
Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#if !defined(__Volume_h)
#define __Volume_h
#include <iostream>
#include <fstream>
#define WANT_STREAM
#define WANT_MATH
#include "newmatap.h"
#include "newmatio.h"
#include "VolumeSeries.h"
#include <string>
using namespace NEWMAT;
namespace FILM {
class Volume : public ColumnVector
{
public:
Volume() : ColumnVector() {}
Volume(const VolumeSeries::Dims& pdims, const ColumnVector& in) :
ColumnVector(),
dims(pdims),
preThresholdPositions(in)
{}
Volume(int psize) : ColumnVector (psize) {}
Volume(int psize,const VolumeSeries::Dims& pdims, const ColumnVector& in) :
ColumnVector (psize),
dims(pdims),
preThresholdPositions(in){}
Volume& operator=(const Volume& vol) {
ColumnVector::operator=(vol);
preThresholdPositions = vol.preThresholdPositions;
dims = vol.dims;
return *this;
}
Volume& operator=(const ColumnVector& pvec) {
ColumnVector::operator=(pvec);
return *this;
}
Volume& operator=(float pin) {
ColumnVector::operator=(pin);
return *this;
}
Volume(const ColumnVector& pvec) {ColumnVector::operator=(pvec);}
void read(const string& fname);
void threshold(float thresh);
void threshold();
void unthreshold();
void unthreshold(const VolumeSeries::Dims& pdims, const ColumnVector& in);
void writeAsInt(const string& fname);
void writeAsFloat(const string& fname);
const VolumeSeries::Dims& getDims() const { return dims; }
void setDims(const VolumeSeries::Dims& pdims) { dims = pdims; }
int getUnthresholdSize() const { return dims.x*dims.y*dims.z; }
int getVolumeSize() const { return Nrows(); }
const ColumnVector& getPreThresholdPositions() const { return preThresholdPositions; }
void setPreThresholdPositions(const ColumnVector& in) { preThresholdPositions = in; }
protected:
VolumeSeries::Dims dims;
ColumnVector preThresholdPositions;
};
}
#endif
VolumeSeries.cc deleted 100644 → 0
+ 0
− 280
View file @ 830f1051
/* VolumeSeries.cc
Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#include <iostream>
#include <cstdlib>
#include "avwio/avwio.h"
#include "newmatap.h"
#include "newmatio.h"
#include "VolumeSeries.h"
#include "miscmaths/miscmaths.h"
using namespace NEWMAT;
namespace FILM {
void VolumeSeries::replaceMeans()
{
Tracer ts("VolumeSeries::replaceMeans");
for(int i = 1; i <= getNumSeries(); i++)
{
Column(i) = Column(i) + means(i);
}
}
void VolumeSeries::unthresholdSeries()
{
Tracer ts("VolumeSeries::unthresholdSeries");
cerr << "Unthresholding Series... ";
int numUnThresholdedSeries = getUnthresholdNumSeries();
int numThresholdedSeries = getNumSeries();
int numVolumes = getNumVolumes();
// Increase this to required size and set to 0
Release();
Matrix X=*this;
ReSize(numVolumes, numUnThresholdedSeries);
Matrix::operator=(0);
// Loop through restoring non thresholded series:
for(int i = 1; i <= numThresholdedSeries; i++)
{
Column(int(preThresholdPositions(i)))=X.Column(i);
}
cerr << "Completed" << endl;
}
void VolumeSeries::thresholdSeries()
{
Tracer ts("VolumeSeries::thresholdSeries");
cerr << "Thresholding Series... ";
int numSeries = preThresholdPositions.Nrows();
int numVolumes = getNumVolumes();
Matrix X(numVolumes, numSeries);
for(int i = 1; i <= numSeries; i++)
{
X.Column(i)=Column(int(preThresholdPositions(i)));
}
*this = X;
cerr << "Completed" << endl;
}
void VolumeSeries::thresholdSeries(float thresh, bool removeMean)
{
Tracer ts("VolumeSeries::thresholdSeries");
cerr << "Applying threshold to Series... ";
int numSeries = getNumSeries();
int j = 0;
if(removeMean)
{
means.ReSize(numSeries);
means = 0;
}
preThresholdPositions.ReSize(numSeries);
float m = 0,s = 0;
for(int i = 1; i <= numSeries; i++)
{
m = MISCMATHS::mean(ColumnVector(getSeries(i))).AsScalar();
s = MISCMATHS::var(ColumnVector(getSeries(i))).AsScalar();
/*
if(m > thresh && s == 0.0)
{
cerr << "m = " << m << endl;
cerr << "s = " << s << endl;
cerr << "i = " << i << endl;
cerr << "j = " << j+1 << endl;
}
*/
if(m > thresh && s > 1e-10)
{
j++;
preThresholdPositions(j) = i;
if(removeMean)
{
Column(j) = getSeries(i) - m;
means(i) = m;
}
else
Column(j) = getSeries(i);
}
}
// Discard rest:
*this = Columns(1,j);
preThresholdPositions = preThresholdPositions.Rows(1,j);
cerr << "Completed" << endl;
}
void VolumeSeries::removeSeriesMeans()
{
int numSeries = getNumSeries();
for(int i = 1; i <= numSeries; i++)
{
float m = MISCMATHS::mean(ColumnVector(getSeries(i))).AsScalar();
Column(i) = getSeries(i) - m;
}
}
void VolumeSeries::read(const string& fname)
{
Tracer ts("VolumeSeries::read");
cerr << "Loading image...";
AVW* IP = AvwOpen(fname.c_str(), "r");
Matrix& output = *this;
short x,y,z,v,type;
float vx,vy,vz,tr;
AvwGetDim(IP,&x,&y,&z,&v);
AvwGetVoxDim(IP,&vx,&vy,&vz,&tr);
dims.x = x; dims.y = y; dims.z = z; dims.v = v;
dims.vx = vx; dims.vy = vy; dims.vz = vz; dims.tr = tr;
size_t imagesize=x*y*z*v;
AvwGetDataType(IP,&type);
output.ReSize(v,x*y*z);
switch(type)
{
case DT_SIGNED_SHORT:
{
short* sbuffer=new short[imagesize];
AvwReadVolumes(IP, sbuffer,v);
size_t volsize = dims.x*dims.y*dims.z;
size_t volstart = 1;
for(size_t i=1; i<=(size_t)v; i++)
{
volstart = (i-1)*volsize;
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
output(i,j)=sbuffer[volstart+j-1];
}
}
delete[] sbuffer;
}
break;
case DT_FLOAT:
{
float* fbuffer=new float[imagesize];
AvwReadVolumes(IP,fbuffer,v);
size_t volsize = dims.x*dims.y*dims.z;
size_t volstart = 1;
for(size_t i=1; i<=(size_t)v; i++)
{
volstart = (i-1)*volsize;
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
output(i,j)=fbuffer[volstart+j-1];
}
}
delete[] fbuffer;
}
break;
case DT_UNSIGNED_CHAR:
{
unsigned char* cbuffer=new unsigned char[imagesize];
AvwReadVolumes(IP,cbuffer,v);
size_t volsize = dims.x*dims.y*dims.z;
size_t volstart = 1;
for(size_t i=1; i<=(size_t)v; i++)
{
volstart = (i-1)*volsize;
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
output(i,j)=cbuffer[volstart+j-1];
}
}
delete[] cbuffer;
}
break;
default:
perror("AvwRead: DT not supported");
}
cerr << " Completed" << endl;
AvwClose(IP);
return;
}
void VolumeSeries::unthresholdSeries(const Dims& pdims,const ColumnVector& in)
{
dims = pdims;
preThresholdPositions = in;
unthresholdSeries();
}
void VolumeSeries::writeAsFloat(const string& fname)
{
Tracer ts("VolumeSeries::writeAsFloat");
cerr << "Saving image...";
AVW* OP = AvwOpen(fname.c_str(), "wc");
AvwSetDim(OP,dims.x, dims.y, dims.z, dims.v);
AvwSetVoxDim(OP,dims.vx, dims.vy, dims.vz, dims.tr);
AvwSetDataType(OP, DT_FLOAT);
int volStart = 1;
int volSize = getNumSeries();
int volNum = getNumVolumes();
float *qv = new float[volSize*volNum];
for(int i = 1; i<= volNum; i++)
{
volStart = (i-1)*volSize;
for(int j = 1; j <= volSize; j++)
{
qv[volStart+j-1] = (*this)(i,j);
}
}
AvwWriteVolumes(OP, qv, volNum);
delete [] qv;
AvwClose(OP);
cerr << " Completed" << endl;
}
}
VolumeSeries.h deleted 100644 → 0
+ 0
− 126
View file @ 830f1051
/* VolumeSeries.h
Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 1999-2000 University of Oxford */
/* CCOPYRIGHT */
#include <iostream>
#include <fstream>
#define WANT_STREAM
#define WANT_MATH
#include "newmatap.h"
#include "newmatio.h"
#include <string>
using namespace NEWMAT;
namespace FILM {
#if !defined(__VolumeSeries_h)
#define __VolumeSeries_h
// Rows are volumes
// Columns are (time) series
// num Rows is size of (time) series
// num Cols is size of volumes
class VolumeSeries : public Matrix
{
public:
struct Dims
{
// Volume dimensions (no. of voxels):
int x;
int y;
int z;
int v;
// Voxel dimensions (mm)
float vx;
float vy;
float vz;
float tr;
};
VolumeSeries() : Matrix() {}
VolumeSeries(const Dims& pdims, const ColumnVector& in) :
Matrix(),
dims(pdims),
preThresholdPositions(in)
{}
VolumeSeries(int pnumVols, int pnumSeries) :
Matrix(pnumVols, pnumSeries),
means(pnumSeries){}
VolumeSeries(int pnumVols, int pnumSeries, const Dims& pdims, const ColumnVector& in) :
Matrix(pnumVols, pnumSeries),
dims(pdims),
preThresholdPositions(in),
means(pnumSeries){}
VolumeSeries& operator=(const VolumeSeries& vol) {
Matrix::operator=(vol);
preThresholdPositions = vol.preThresholdPositions;
dims = vol.dims;
means = vol.means;
return *this;
}
VolumeSeries& operator=(const Matrix& mat) {
Matrix::operator=(mat);
return *this;
}
VolumeSeries& operator=(float pin) {
Matrix::operator=(pin);
return *this;
}
VolumeSeries(const VolumeSeries& vol){operator=(vol);}
VolumeSeries(const Matrix& pmat) : Matrix(pmat) {}
void thresholdSeries(float thresh, bool removeMean);
void thresholdSeries();
void unthresholdSeries();
void unthresholdSeries(const Dims& pdims,const ColumnVector& in);
void removeSeriesMeans();
const ColumnVector& getPreThresholdPositions() const { return preThresholdPositions; }
void setPreThresholdPositions(const ColumnVector& in) { preThresholdPositions = in; }
int getNumVolumes() const { return Nrows(); }
int getNumSeries() const { return Ncols(); }
const GetSubMatrix& getSeries(int i) const { return Column(i); }
GetSubMatrix& getSeries(int i) { return Column(i); }
const GetSubMatrix& getVolume(int i) const { return Row(i); }
GetSubMatrix& getVolume(int i) { return Row(i); }
void read(const string& fname);
void writeAsInt(const string& fname);
void writeAsFloat(const string& fname);
void replaceMeans();
const Dims& getDims() const { return dims; }
void setDims(const Dims& pdims) { dims = pdims; }
int getUnthresholdNumSeries() const { return dims.x*dims.y*dims.z; }
protected:
Dims dims;
ColumnVector preThresholdPositions;
ColumnVector means;
};
#endif
}
glmrand.h
+ 2
− 2
View file @ 77af356f
......@@ -16,8 +16,8 @@
#include "newmatap.h"
#include "newmatio.h"
#include "Volume.h"
#include "VolumeSeries.h"
#include "miscmaths/volume.h"
#include "miscmaths/volumeseries.h"
using namespace NEWMAT;
......
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