diff --git a/volume.cc b/volume.cc deleted file mode 100644 index 6e7edb83367a5040a89818bb80caf457b6d435ef..0000000000000000000000000000000000000000 --- a/volume.cc +++ /dev/null @@ -1,260 +0,0 @@ -/* 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; - } -} - - - - - - - - - - - - - diff --git a/volume.h b/volume.h deleted file mode 100644 index fa11e0504bfe9ffa490f80481acb040293d66a9e..0000000000000000000000000000000000000000 --- a/volume.h +++ /dev/null @@ -1,89 +0,0 @@ -/* 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 - - - - diff --git a/volumeinfo.h b/volumeinfo.h deleted file mode 100644 index 804d5e69d45d3cd70f2b7f7eb3dc0ab8eea23f33..0000000000000000000000000000000000000000 --- a/volumeinfo.h +++ /dev/null @@ -1,42 +0,0 @@ -/* 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 - -} diff --git a/volumeseries.cc b/volumeseries.cc deleted file mode 100644 index 25872c2f6c636bb8ea3a569346844b85ae48d515..0000000000000000000000000000000000000000 --- a/volumeseries.cc +++ /dev/null @@ -1,304 +0,0 @@ -/* 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); - - } - -} - - - diff --git a/volumeseries.h b/volumeseries.h deleted file mode 100644 index 996fcef069bbbaa8c6551f671250c6197bdfeee7..0000000000000000000000000000000000000000 --- a/volumeseries.h +++ /dev/null @@ -1,122 +0,0 @@ -/* 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 - -} - - - -