Skip to content
Snippets Groups Projects
Commit c43ac9d0 authored by Matthew Webster's avatar Matthew Webster
Browse files

removed old volume

parent 8482373a
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
− 260
View file @ 8482373a
/* volume.cc
Mark Woolrich, FMRIB Image Analysis Group
Copyright (C) 2002 University of Oxford */
/* CCOPYRIGHT */
#include <iostream>
#include <cstdlib>
#include "fslio/fslio.h"
#include "newmatap.h"
#include "newmatio.h"
#include "volume.h"
#include "utils/time_tracer.h"
using namespace NEWMAT;
using namespace Utilities;
namespace MISCMATHS {
void Volume::unthreshold()
{
Time_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 VolumeInfo& pvolinfo,const ColumnVector& in)
{
volinfo = pvolinfo;
preThresholdPositions = in;
unthreshold();
}
void Volume::threshold()
{
Time_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)
{
Time_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)
{
Time_Tracer ts(string("Volume::writeAsFloat" + fname).c_str());
const ColumnVector& outputvol = *this;
FSLIO* OP = FslOpen(fname.c_str(), "wb");
FslCloneHeader(OP,volinfo.miscinfo);
FslSetCalMinMax(OP,Minimum(),Maximum());
FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, 1);
FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, 0);
FslSetDataType(OP, DT_FLOAT);
FslSetIntent(OP, volinfo.intent_code, volinfo.intent_p1, volinfo.intent_p2,
volinfo.intent_p3);
int sizeVol = outputvol.Nrows();
float *qv = new float[sizeVol];
for (int i=1; i<=sizeVol; i++)
{
qv[i-1] = outputvol(i);
}
FslWriteHeader(OP);
FslWriteVolumes(OP, qv, 1);
delete [] qv;
FslClose(OP);
}
void Volume::writeAsInt(const string& fname)
{
Time_Tracer ts("Volume::writeAsInt");
const ColumnVector& outputvol = *this;
FSLIO* OP = FslOpen(fname.c_str(), "wb");
FslCloneHeader(OP,volinfo.miscinfo);
FslSetCalMinMax(OP,Minimum(),Maximum());
FslSetDim(OP, volinfo.x, volinfo.y, volinfo.z, 1);
FslSetVoxDim(OP, volinfo.vx, volinfo.vy, volinfo.vz, 0);
FslSetDataType(OP, DT_SIGNED_SHORT);
FslSetIntent(OP, volinfo.intent_code, volinfo.intent_p1, volinfo.intent_p2,
volinfo.intent_p3);
int sizeVol = outputvol.Nrows();
short *qv = new short[sizeVol];
for (int i=1; i<=sizeVol; i++)
{
qv[i-1] = (short)outputvol(i);
}
FslWriteHeader(OP);
FslWriteVolumes(OP, qv, 1);
delete [] qv;
FslClose(OP);
}
void Volume::read(const string& fname)
{
Time_Tracer ts(string("Volume::read" + fname).c_str());
FSLIO* IP = FslOpen(fname.c_str(), "rb");
ColumnVector& output = *this;
short x,y,z,v,type;
float vx,vy,vz,tr;
float slope, intercept;
int doscaling;
FslGetDim(IP,&x,&y,&z,&v);
FslGetVoxDim(IP,&vx,&vy,&vz,&tr);
FslGetIntent(IP, &(volinfo.intent_code), &(volinfo.intent_p1), &(volinfo.intent_p2),
&(volinfo.intent_p3));
doscaling = FslGetIntensityScaling(IP,&slope,&intercept);
volinfo.x = x; volinfo.y = y; volinfo.z = z; volinfo.v = v;
volinfo.vx = vx; volinfo.vy = vy; volinfo.vz = vz; volinfo.tr = tr;
volinfo.miscinfo = FslInit();
FslCloneHeader(volinfo.miscinfo,IP);
size_t imagesize=x*y*z;
FslGetDataType(IP,&type);
output.ReSize(x*y*z);
switch(type)
{
case DT_SIGNED_SHORT:
{
short* sbuffer=new short[imagesize];
FslReadVolumes(IP,sbuffer,v);
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
if (doscaling==0) { output(j)=sbuffer[j-1]; }
else { output(j)=(slope * sbuffer[j-1]) + intercept; }
}
delete[] sbuffer;
}
break;
case DT_FLOAT:
{
float* fbuffer=new float[imagesize];
FslReadVolumes(IP,fbuffer,v);
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
if (doscaling==0) { output(j)=fbuffer[j-1]; }
else { output(j)=(slope * fbuffer[j-1]) + intercept; }
}
delete[] fbuffer;
}
break;
case DT_UNSIGNED_CHAR:
{
unsigned char* cbuffer=new unsigned char[imagesize];
FslReadVolumes(IP,cbuffer,v);
for(size_t j = 1; j<=(size_t)x*y*z; j++)
{
if (doscaling==0) { output(j)=cbuffer[j-1]; }
else { output(j)=(slope * cbuffer[j-1]) + intercept; }
}
delete[] cbuffer;
}
break;
default:
perror("FslRead: DT not supported");
}
FslClose(IP);
return;
}
}
volume.h deleted 100644 → 0
+ 0
− 89
View file @ 8482373a
/* volume.h
Mark Woolrich - FMRIB Image Analysis Group
Copyright (C) 2002 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 "volumeinfo.h"
#include <string>
using namespace NEWMAT;
namespace MISCMATHS {
class Volume : public ColumnVector
{
public:
Volume() : ColumnVector() {}
Volume(const VolumeInfo& pvolinfo, const ColumnVector& in) :
ColumnVector(),
volinfo(pvolinfo),
preThresholdPositions(in)
{}
Volume(int psize) : ColumnVector (psize) {}
Volume(int psize,const VolumeInfo& pvolinfo, const ColumnVector& in) :
ColumnVector (psize),
volinfo(pvolinfo),
preThresholdPositions(in){}
Volume& operator=(const Volume& vol) {
ColumnVector::operator=(vol);
preThresholdPositions = vol.preThresholdPositions;
volinfo = vol.volinfo;
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 std::string& fname);
void threshold(float thresh);
void threshold();
void unthreshold();
void unthreshold(const VolumeInfo& pvolinfo, const ColumnVector& in);
void writeAsInt(const std::string& fname);
void writeAsFloat(const std::string& fname);
const VolumeInfo& getInfo() const { return volinfo; }
void setInfo(const VolumeInfo& pvolinfo) { volinfo = pvolinfo; }
int getUnthresholdSize() const { return volinfo.x*volinfo.y*volinfo.z; }
int getVolumeSize() const { return Nrows(); }
const ColumnVector& getPreThresholdPositions() const { return preThresholdPositions; }
void setPreThresholdPositions(const ColumnVector& in) { preThresholdPositions = in; }
protected:
VolumeInfo volinfo;
ColumnVector preThresholdPositions;
};
}
#endif
volumeinfo.h deleted 100644 → 0
+ 0
− 42
View file @ 8482373a
/* volumeinfo.h
Mark Woolrich - FMRIB Image Analysis Group
Copyright (C) 2002 University of Oxford */
/* CCOPYRIGHT */
#include "fslio/fslio.h"
namespace MISCMATHS {
#if !defined(__VolumeInfo_h)
#define __VolumeInfo_h
struct VolumeInfo
{
// Volume dimensions (no. of voxels):
int x;
int y;
int z;
int v;
// Voxel dimensions (mm)
float vx;
float vy;
float vz;
float tr;
// Intent codes and parameters
short intent_code;
float intent_p1;
float intent_p2;
float intent_p3;
FSLIO* miscinfo;
};
#endif
}
volumeseries.cc deleted 100644 → 0
+ 0
− 304
View file @ 8482373a
/* volumeseries.cc
Mark Woolrich - FMRIB Image Analysis Group
Copyright (C) 2002 University of Oxford */
/* CCOPYRIGHT */
#include <iostream>
#include <cstdlib>
#include "fslio/fslio.h"
#include "newmatap.h"
#include "newmatio.h"
#include "volumeseries.h"
#include "miscmaths.h"
#include "utils/time_tracer.h"
using namespace NEWMAT;
using namespace Utilities;
namespace MISCMATHS {
void VolumeSeries::replaceMeans()
{
Time_Tracer ts("VolumeSeries::replaceMeans");
for(int i = 1; i <= getNumSeries(); i++)
{
Column(i) = Column(i) + means(i);
}
}
void VolumeSeries::unthresholdSeries()
{
Time_Tracer ts("VolumeSeries::unthresholdSeries");
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);
}
}
void VolumeSeries::thresholdSeries()
{
Time_Tracer ts("VolumeSeries::thresholdSeries");
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;
}
void VolumeSeries::thresholdSeries(float thresh, bool removeMean)
{
Time_Tracer ts("VolumeSeries::thresholdSeries");
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);
}
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)
{
Time_Tracer ts(string("VolumeSeries::read-" + fname).c_str());
FSLIO* IP = FslOpen(fname.c_str(), "rb");
Matrix& output = *this;
short x,y,z,v,type;
float vx,vy,vz,tr;
float slope, intercept;
int doscaling;
FslGetDim(IP,&x,&y,&z,&v);
FslGetVoxDim(IP,&vx,&vy,&vz,&tr);
FslGetIntent(IP, &(volinfo.intent_code), &(volinfo.intent_p1), &(volinfo.intent_p2),
&(volinfo.intent_p3));
doscaling = FslGetIntensityScaling(IP,&slope,&intercept);
volinfo.x = x; volinfo.y = y; volinfo.z = z; volinfo.v = v;
volinfo.vx = vx; volinfo.vy = vy; volinfo.vz = vz; volinfo.tr = tr;
volinfo.miscinfo = FslInit();
FslCloneHeader(volinfo.miscinfo,IP);
size_t imagesize=x*y*z;
FslGetDataType(IP,&type);
output.ReSize(v,x*y*z);
for(size_t i=1; i<=(size_t)v; i++)
{
switch(type)
{
case DT_SIGNED_SHORT:
{
short* sbuffer=new short[imagesize];
FslReadVolumes(IP, sbuffer,1);
for(size_t j = 1; j<=imagesize; j++)
{
if (doscaling==0) { output(i,j)=sbuffer[j-1]; }
else { output(i,j)=(slope * sbuffer[j-1]) + intercept;}
}
delete[] sbuffer;
}
break;
case DT_FLOAT:
{
float* fbuffer=new float[imagesize];
FslReadVolumes(IP,fbuffer,1);
for(size_t j = 1; j<=imagesize; j++)
{
if (doscaling==0) { output(i,j)=fbuffer[j-1]; }
else { output(i,j)=(slope * fbuffer[j-1]) + intercept;}
}
delete[] fbuffer;
}
break;
case DT_UNSIGNED_CHAR:
{
unsigned char* cbuffer=new unsigned char[imagesize];
FslReadVolumes(IP,cbuffer,1);
for(size_t j = 1; j<=imagesize; j++)
{
if (doscaling==0) { output(i,j)=cbuffer[j-1]; }
else { output(i,j)=(slope * cbuffer[j-1]) + intercept;}
}
delete[] cbuffer;
}
break;
default:
perror("FslRead: DT not supported");
}
}
FslClose(IP);
return;
}
void VolumeSeries::unthresholdSeries(const VolumeInfo& pvolinfo,const ColumnVector& in)
{
volinfo = pvolinfo;
preThresholdPositions = in;
unthresholdSeries();
}
void VolumeSeries::writeAsFloat(const string& fname)
{
Time_Tracer ts(string("VolumeSeries::writeAsFloat" + fname).c_str());
FSLIO* OP = FslOpen(fname.c_str(), "wb");
FslCloneHeader(OP,volinfo.miscinfo);
FslSetCalMinMax(OP,Minimum(),Maximum());
FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, volinfo.v);
FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, volinfo.tr);
FslSetDataType(OP, DT_FLOAT);
FslSetIntent(OP, volinfo.intent_code, volinfo.intent_p1, volinfo.intent_p2,
volinfo.intent_p3);
int volSize = getNumSeries();
int volNum = getNumVolumes();
// This should assert() but I'm too scared to add that to production code!
if ( (volinfo.x * volinfo.y * volinfo.z != volSize) || (volinfo.v != volNum) )
cerr << "WARNING: Dimensions are inconsistent in VolumeSeries::writeAsFloat()\n"
<< "x, y, z, v = " << volinfo.x << ", " << volinfo.y << ", " << volinfo.z
<< ", " << volinfo.v << "\nvolSize, volNum = " << volSize << ", "
<< volNum << "\nThis is probably NOT what you intended!!" << endl;
FslWriteHeader(OP);
float *qv = new float[volSize];
for(int i = 1; i<= volNum; i++)
{
for(int j = 1; j <= volSize; j++)
qv[j-1] = (*this)(i,j);
FslWriteVolumes(OP, qv, 1);
}
delete [] qv;
FslClose(OP);
}
void VolumeSeries::writeThresholdedSeriesAsFloat(const VolumeInfo& pvolinfo,const ColumnVector& in,const string& fname)
{
volinfo = pvolinfo;
preThresholdPositions = in;
Time_Tracer ts(string("VolumeSeries::writeThresholdedSeriesAsFloat" + fname).c_str());
FSLIO* OP = FslOpen(fname.c_str(), "wb");
FslCloneHeader(OP,volinfo.miscinfo);
FslSetCalMinMax(OP,Minimum(),Maximum());
FslSetDim(OP,volinfo.x, volinfo.y, volinfo.z, volinfo.v);
FslSetVoxDim(OP,volinfo.vx, volinfo.vy, volinfo.vz, volinfo.tr);
FslSetDataType(OP, DT_FLOAT);
FslSetIntent(OP, volinfo.intent_code, volinfo.intent_p1, volinfo.intent_p2, volinfo.intent_p3);
int volSize = getUnthresholdNumSeries();
int numThresholdedSeries = getNumSeries();
int volNum = getNumVolumes();
FslWriteHeader(OP);
float *qv = new float[volSize];
for(int i = 1; i<= volNum; i++)
{
for(int j = 1; j <= volSize; j++)
qv[j-1]=0;
for(int j = 1; j <= numThresholdedSeries; j++)
qv[int(preThresholdPositions(j))-1] = (*this)(i,j);
FslWriteVolumes(OP, qv, 1);
}
delete [] qv;
FslClose(OP);
}
}
volumeseries.h deleted 100644 → 0
+ 0
− 122
View file @ 8482373a
/* volumeseries.h
Mark Woolrich - FMRIB Image Analysis Group
Copyright (C) 2002 University of Oxford */
/* CCOPYRIGHT */
#include <iostream>
#include <fstream>
#define WANT_STREAM
#define WANT_MATH
#include "newmatap.h"
#include "newmatio.h"
#include "volumeinfo.h"
#include <string>
using namespace NEWMAT;
namespace MISCMATHS {
#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:
VolumeSeries() : Matrix() {}
VolumeSeries(const VolumeInfo& pvolinfo, const ColumnVector& in) :
Matrix(),
volinfo(pvolinfo),
preThresholdPositions(in)
{}
VolumeSeries(int pnumVols, int pnumSeries) :
Matrix(pnumVols, pnumSeries),
means(pnumSeries){}
VolumeSeries(int pnumVols, int pnumSeries, const VolumeInfo& pvolinfo, const ColumnVector& in) :
Matrix(pnumVols, pnumSeries),
volinfo(pvolinfo),
preThresholdPositions(in),
means(pnumSeries){}
VolumeSeries& operator=(const VolumeSeries& vol) {
Matrix::operator=(vol);
preThresholdPositions = vol.preThresholdPositions;
volinfo = vol.volinfo;
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 VolumeInfo& pvolinfo,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(); }
int nvoxels() const { return Ncols(); }
int tsize() const { return Nrows(); }
const ReturnMatrix getSeries(int i) const { ColumnVector tmp = Column(i);tmp.Release();return tmp; }
ReturnMatrix getSeries(int i) { ColumnVector tmp = Column(i);tmp.Release();return tmp; }
void setSeries(const ColumnVector& in, int i) { Column(i) = in; }
const ReturnMatrix getVolume(int i) const { RowVector tmp = Row(i);tmp.Release();return tmp; }
ReturnMatrix getVolume(int i) { RowVector tmp = Row(i);tmp.Release();return tmp; }
void setVolume(const ColumnVector& in, int i) { Row(i) = in.AsRow(); }
void setVolume(const RowVector& in, int i) { Row(i) = in; }
void read(const std::string& fname);
void writeAsInt(const std::string& fname);
void writeAsFloat(const std::string& fname);
void writeThresholdedSeriesAsFloat(const VolumeInfo& pvolinfo,const ColumnVector& in,const std::string& fname);
void replaceMeans();
const VolumeInfo& getInfo() const { return volinfo; }
void setInfo(const VolumeInfo& pvolinfo) { volinfo = pvolinfo; }
int getUnthresholdNumSeries() const { return volinfo.x*volinfo.y*volinfo.z; }
protected:
VolumeInfo volinfo;
ColumnVector preThresholdPositions;
ColumnVector means;
};
#endif
}
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