diff --git a/code/core/rttbBaseType.h b/code/core/rttbBaseType.h index 8bc5699..e1ef5fe 100644 --- a/code/core/rttbBaseType.h +++ b/code/core/rttbBaseType.h @@ -1,743 +1,743 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __BASE_TYPE_NEW_H #define __BASE_TYPE_NEW_H #include #include #include #include #include #include #include #include namespace rttb { const double errorConstant = 1e-5; const double reducedErrorConstant = 0.0001; typedef unsigned short UnsignedIndex1D; /*! @class UnsignedIndex2D @brief 2D index. */ class UnsignedIndex2D: public boost::numeric::ublas::vector { public: UnsignedIndex2D() : boost::numeric::ublas::vector(2) {} UnsignedIndex2D(const UnsignedIndex1D value) : boost::numeric::ublas::vector(2, value) {} const UnsignedIndex1D x() const { return (*this)(0); } const UnsignedIndex1D y() const { return (*this)(1); } }; /*! @class UnsignedIndex3D @brief 3D index. */ class UnsignedIndex3D: public boost::numeric::ublas::vector { public: UnsignedIndex3D() : boost::numeric::ublas::vector(3) {} UnsignedIndex3D(const UnsignedIndex1D value) : boost::numeric::ublas::vector(3, value) {} UnsignedIndex3D(const UnsignedIndex1D xValue, const UnsignedIndex1D yValue, const UnsignedIndex1D zValue) : boost::numeric::ublas::vector(3, xValue) { (*this)(1) = yValue; (*this)(2) = zValue; } const UnsignedIndex1D x() const { return (*this)(0); } const UnsignedIndex1D y() const { return (*this)(1); } const UnsignedIndex1D z() const { return (*this)(2); } friend bool operator==(const UnsignedIndex3D& gi1, const UnsignedIndex3D& gi2) { if (gi1.size() != gi2.size()) { return false; } for (size_t i = 0; i < gi1.size(); i++) { if (gi1(i) != gi2(i)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const UnsignedIndex3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; typedef std::list UnsignedIndexList; typedef std::string FileNameString; typedef std::string ContourGeometricTypeString; typedef double WorldCoordinate; /*! @class WorldCoordinate2D @brief 2D coordinate in real world coordinates. */ class WorldCoordinate2D: public boost::numeric::ublas::vector { public: WorldCoordinate2D() : boost::numeric::ublas::vector (2) {} WorldCoordinate2D(const WorldCoordinate value) : boost::numeric::ublas::vector(2, value) {} WorldCoordinate2D(const WorldCoordinate xValue, const WorldCoordinate yValue) : boost::numeric::ublas::vector(2, xValue) { (*this)(1) = yValue; } const WorldCoordinate x() const { return (*this)(0); } const WorldCoordinate y() const { return (*this)(1); } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y(); return ss.str(); } friend bool operator==(const WorldCoordinate2D& wc1, const WorldCoordinate2D& wc2) { if (wc1.size() != wc2.size()) { return false; } for (size_t i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } }; /*! @class WorldCoordinate3D @brief 3D coordinate in real world coordinates like in DICOM to define ImagePositionPatient. */ class WorldCoordinate3D: public boost::numeric::ublas::vector { public: WorldCoordinate3D() : boost::numeric::ublas::vector(3) {} WorldCoordinate3D(const WorldCoordinate value) : boost::numeric::ublas::vector(3, value) {} WorldCoordinate3D(const WorldCoordinate xValue, const WorldCoordinate yValue, const WorldCoordinate zValue) : boost::numeric::ublas::vector(3, xValue) { (*this)(1) = yValue; (*this)(2) = zValue; } WorldCoordinate3D(const WorldCoordinate3D& w): boost::numeric::ublas::vector(3) { (*this)(0) = w.x(); (*this)(1) = w.y(); (*this)(2) = w.z(); } const WorldCoordinate x() const { return (*this)(0); } const WorldCoordinate y() const { return (*this)(1); } const WorldCoordinate z() const { return (*this)(2); } //vector cross product (not included in boost.ublas) WorldCoordinate3D cross(WorldCoordinate3D aVector) const { WorldCoordinate3D result; WorldCoordinate x = (*this)(0); WorldCoordinate y = (*this)(1); WorldCoordinate z = (*this)(2); result(0) = y * aVector(2) - z * aVector(1); result(1) = z * aVector(0) - x * aVector(2); result(2) = x * aVector(1) - y * aVector(0); return result; } WorldCoordinate2D getXY() const { WorldCoordinate2D result; result(0) = (*this)(0); result(1) = (*this)(1); return result; } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y() << ' ' << z(); return ss.str(); } WorldCoordinate3D& operator=(const WorldCoordinate3D& wc) { (*this)(0) = wc.x(); (*this)(1) = wc.y(); (*this)(2) = wc.z(); return (*this); } WorldCoordinate3D& operator=(const boost::numeric::ublas::vector wc) { (*this)(0) = wc(0); (*this)(1) = wc(1); (*this)(2) = wc(2); return (*this); } WorldCoordinate3D operator-(const boost::numeric::ublas::vector wc) { return WorldCoordinate3D((*this)(0) - wc(0), (*this)(1) - wc(1), (*this)(2) - wc(2)); } WorldCoordinate3D operator+(const boost::numeric::ublas::vector wc) { return WorldCoordinate3D((*this)(0) + wc(0), (*this)(1) + wc(1), (*this)(2) + wc(2)); } friend bool operator==(const WorldCoordinate3D& wc1, const WorldCoordinate3D& wc2) { if (wc1.size() != wc2.size()) { return false; } for (size_t i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } bool equalsAlmost(const WorldCoordinate3D& another, double errorConstantWC = 1e-5) const { if (size() != another.size()) { return false; } double dist = norm_2(*this - another); return dist < errorConstantWC; } friend std::ostream& operator<<(std::ostream& s, const WorldCoordinate3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; /* ! @brief continuous index */ typedef WorldCoordinate3D DoubleVoxelGridIndex3D; typedef UnsignedIndex3D ImageSize; /*! @deprecated Use OrientationMatrix instead. */ typedef WorldCoordinate3D ImageOrientation; typedef double GridVolumeType; /*! @class SpacingVectorType3D @brief 3D spacing vector. @pre values of this vector may not be negative. */ class SpacingVectorType3D: public boost::numeric::ublas::vector { public: SpacingVectorType3D() : boost::numeric::ublas::vector(3) {} SpacingVectorType3D(const GridVolumeType value) : boost::numeric::ublas::vector(3, value) {} SpacingVectorType3D(const GridVolumeType xValue, const GridVolumeType yValue, const GridVolumeType zValue) : boost::numeric::ublas::vector(3, xValue) { (*this)(1) = yValue; (*this)(2) = zValue; } SpacingVectorType3D(const SpacingVectorType3D& w): boost::numeric::ublas::vector(3) { (*this)(0) = w.x(); (*this)(1) = w.y(); (*this)(2) = w.z(); } const GridVolumeType x() const { return (*this)(0); } const GridVolumeType y() const { return (*this)(1); } const GridVolumeType z() const { return (*this)(2); } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y() << ' ' << z(); return ss.str(); } SpacingVectorType3D& operator=(const SpacingVectorType3D& wc) { (*this)(0) = wc.x(); (*this)(1) = wc.y(); (*this)(2) = wc.z(); return (*this); } SpacingVectorType3D& operator=(const WorldCoordinate3D& wc) { (*this)(0) = GridVolumeType(wc.x()); (*this)(1) = GridVolumeType(wc.y()); (*this)(2) = GridVolumeType(wc.z()); return (*this); } SpacingVectorType3D& operator=(const boost::numeric::ublas::vector wc) { (*this)(0) = wc(0); (*this)(1) = wc(1); (*this)(2) = wc(2); return (*this); } friend bool operator==(const SpacingVectorType3D& wc1, const SpacingVectorType3D& wc2) { if (wc1.size() != wc2.size()) { return false; } for (size_t i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } bool equalsAlmost(const SpacingVectorType3D& another, double errorConstantSV) const { if ((*this).size() != another.size()) { return false; } double dist = norm_2(*this - another); return dist < errorConstantSV; } friend std::ostream& operator<<(std::ostream& s, const SpacingVectorType3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; /*! @class OrientationMatrix @brief Used to store image orientation information */ class OrientationMatrix : public boost::numeric::ublas::matrix { public: /*! The default constructor generates a 3x3 unit matrix */ OrientationMatrix() : boost::numeric::ublas::matrix(3, 3, 0) { for (std::size_t m = 0; m < (*this).size1(); m++) { (*this)(m, m) = 1; } } OrientationMatrix(const WorldCoordinate value) : boost::numeric::ublas::matrix(3, 3, value) {} bool equalsAlmost(const OrientationMatrix& anOrientationMatrix, double errorConstantOM) const { if (anOrientationMatrix.size1() == (*this).size1()) { if (anOrientationMatrix.size2() == (*this).size2()) { for (std::size_t m = 0; m < anOrientationMatrix.size1(); m++) { for (std::size_t n = 0; n < anOrientationMatrix.size2(); n++) { if ((std::abs((*this)(m, n) - anOrientationMatrix(m, n)) > errorConstantOM)) { return false; } } }// end element comparison } else { return false; } } else { return false; } return true; } friend bool operator==(const OrientationMatrix& om1, const OrientationMatrix& om2) { return om1.equalsAlmost(om2, 0.0); } friend std::ostream& operator<<(std::ostream& s, const OrientationMatrix& anOrientationMatrix) { s << "[ "; for (std::size_t m = 0; m < anOrientationMatrix.size1(); m++) { s << "[ "; for (std::size_t n = 0; n < anOrientationMatrix.size2(); n++) { if (n == 0) { s << anOrientationMatrix(m, n); } else { s << ", " << anOrientationMatrix(m, n); } } s << " ]"; } s << " ]"; return s; } }; /*! base for 2D and 3D VoxelIndex; Therefore required beside VoxelGridID */ typedef unsigned int GridIndexType; /*! @class VoxelGridIndex3D @brief 3D voxel grid index in a discret geometry (matrix/image). @details analogous to DICOM where ImagePositionPatient gives the position of the center of the first coordinate (0/0/0) */ class VoxelGridIndex3D: public boost::numeric::ublas::vector { public: VoxelGridIndex3D() : boost::numeric::ublas::vector(3) {} VoxelGridIndex3D(const GridIndexType value) : boost::numeric::ublas::vector(3, value) {} VoxelGridIndex3D(const GridIndexType xValue, const GridIndexType yValue, const GridIndexType zValue) : boost::numeric::ublas::vector(3, xValue) { (*this)(1) = yValue; (*this)(2) = zValue; } const GridIndexType x() const { return (*this)(0); } const GridIndexType y() const { return (*this)(1); } const GridIndexType z() const { return (*this)(2); } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y() << ' ' << z(); return ss.str(); } VoxelGridIndex3D& operator=(const UnsignedIndex3D& ui) { (*this)(0) = ui(0); (*this)(1) = ui(1); (*this)(2) = ui(2); return (*this); } friend bool operator==(const VoxelGridIndex3D& gi1, const VoxelGridIndex3D& gi2) { if (gi1.size() != gi2.size()) { return false; } for (size_t i = 0; i < gi1.size(); i++) { if (gi1(i) != gi2(i)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const VoxelGridIndex3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; /*! @class VoxelGridIndex3D @brief 2D voxel grid index. */ class VoxelGridIndex2D: public boost::numeric::ublas::vector { public: VoxelGridIndex2D() : boost::numeric::ublas::vector(2) {} VoxelGridIndex2D(const GridIndexType value) : boost::numeric::ublas::vector(2, value) {} VoxelGridIndex2D(const GridIndexType xValue, const GridIndexType yValue) : boost::numeric::ublas::vector(2, xValue) { (*this)(1) = yValue; } const GridIndexType x() const { return (*this)(0); } const GridIndexType y() const { return (*this)(1); } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y(); return ss.str(); } friend bool operator==(const VoxelGridIndex2D& gi1, const VoxelGridIndex2D& gi2) { if (gi1.size() != gi2.size()) { return false; } for (size_t i = 0; i < gi1.size(); i++) { if (gi1(i) != gi2(i)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const VoxelGridIndex2D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << " ]"; return s; } }; typedef long GridSizeType; typedef int VoxelGridID; //starts from 0 and is continuously counting all positions on the grid typedef unsigned int VoxelGridDimensionType; typedef boost::numeric::ublas::vector VoxelGridDimensionVectorType; typedef double FractionType, VoxelSizeType, DVHVoxelNumber; typedef double DoseCalcType, DoseTypeGy, GenericValueType, DoseVoxelVolumeType, VolumeType, - GridVolumeType, + GridVolumeType, PercentType, VoxelNumberType, BEDType, LQEDType; typedef std::string IDType; typedef std::string StructureLabel; struct DVHRole { enum Type { TargetVolume = 1, HealthyTissue = 2, WholeVolume = 4, UserDefined = 128 } Type; }; struct DVHType { enum Type { Differential = 1, Cumulative = 2 } Type; }; typedef std::string PatientInfoString; typedef std::string TimeString; typedef std::string FileNameType; typedef std::vector PolygonType; typedef std::vector PolygonSequenceType; typedef double IndexValueType; typedef double DoseStatisticType; typedef std::string DBStringType; typedef std::string DICOMRTFileNameString; typedef unsigned short Uint16; struct Area2D { WorldCoordinate x_begin; WorldCoordinate x_end; WorldCoordinate y_begin; WorldCoordinate y_end; VoxelGridIndex2D index_begin; VoxelGridIndex2D index_end; void Init() { x_begin = -1000000; x_end = -1000000; y_begin = -1000000; y_end = -1000000; index_begin(0) = 0; index_begin(1) = 0; index_end(0) = 0; index_end(1) = 0; } }; struct Segment2D { WorldCoordinate2D begin; WorldCoordinate2D end; }; struct Segment3D { WorldCoordinate3D begin; WorldCoordinate3D end; }; typedef int DistributionType; typedef std::string PathType; typedef std::string XMLString, StatisticsString; }//end: namespace rttb #endif diff --git a/code/core/rttbDVH.cpp b/code/core/rttbDVH.cpp index fd36886..1297e15 100644 --- a/code/core/rttbDVH.cpp +++ b/code/core/rttbDVH.cpp @@ -1,419 +1,434 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include "rttbDVH.h" #include "rttbException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace core { DVH::~DVH() {} DVH::DVH(const DataDifferentialType& aDataDifferential, const DoseTypeGy& aDeltaD, const DoseVoxelVolumeType& aDeltaV, IDType aStructureID, IDType aDoseID): _deltaD(aDeltaD), _deltaV(aDeltaV), _structureID(aStructureID), _doseID(aDoseID), _voxelizationID("NONE") { _dataDifferential.clear(); _dataDifferential = aDataDifferential; this->init(); } DVH::DVH(const DataDifferentialType& aDataDifferential, DoseTypeGy aDeltaD, DoseVoxelVolumeType aDeltaV, IDType aStructureID, IDType aDoseID, IDType aVoxelizationID): _deltaD(aDeltaD), _deltaV(aDeltaV), _structureID(aStructureID), _doseID(aDoseID), _voxelizationID(aVoxelizationID) { _dataDifferential.clear(); _dataDifferential = aDataDifferential; this->init(); } DVH::DVH(const DVH& copy) : _structureID(copy._structureID), _doseID(copy._doseID), _voxelizationID(copy._voxelizationID), _label(copy._label) { _deltaD = copy._deltaD; _deltaV = copy._deltaV; _dataDifferential.clear(); _dataDifferential = copy._dataDifferential; this->init(); } DVH& DVH::operator=(const DVH& copy) { if (this != ©) { _deltaD = copy._deltaD; _deltaV = copy._deltaV; _structureID = copy._structureID; _doseID = copy._doseID; _voxelizationID = copy._voxelizationID; _label = copy._label; _dataDifferential.clear(); _dataDifferential = copy._dataDifferential; } this->init(); return *this; } bool operator==(const DVH& aDVH, const DVH& otherDVH) { if (aDVH.getStructureID() != otherDVH.getStructureID()) { return false; } if (aDVH.getDoseID() != otherDVH.getDoseID()) { return false; } if (aDVH.getVoxelizationID() != otherDVH.getVoxelizationID()) { return false; } if (aDVH.getNumberOfVoxels() != otherDVH.getNumberOfVoxels()) { return false; } return true; } std::ostream& operator<<(std::ostream& s, const DVH& aDVH) { s << "[ " << aDVH.getStructureID() << ", " << aDVH.getDoseID() << ", " << aDVH.getVoxelizationID() << "\n " << "Number of Voxels: " << aDVH.getNumberOfVoxels() << " " << "Minimum/Maximum/Mean/Standard deviation: " << aDVH.getMinimum() << ", " << aDVH.getMaximum() << ", " << aDVH.getMean() << ", " << aDVH.getStdDeviation() << " ]"; return s; } std::deque DVH::getDataDifferential(bool relativeVolume) const { if (!relativeVolume) { return _dataDifferential; } else { return _dataDifferentialRelative; } } DoseVoxelVolumeType DVH::getDeltaV() const { return _deltaV; } DoseTypeGy DVH::getDeltaD() const { return _deltaD; } IDType DVH::getDoseID() const { return this->_doseID; } IDType DVH::getStructureID() const { return this->_structureID; } IDType DVH::getVoxelizationID() const { return this->_voxelizationID; } void DVH::setDoseID(IDType aDoseID) { _doseID = aDoseID; } void DVH::setStructureID(IDType aStrID) { _structureID = aStrID; } DoseStatisticType DVH::getMaximum() const { return _maximum; } DoseStatisticType DVH::getMinimum() const { return _minimum; } DoseStatisticType DVH::getMean() const { return _mean; } DVHVoxelNumber DVH::getNumberOfVoxels() const { return _numberOfVoxels; } DoseStatisticType DVH::getStdDeviation() const { return _stdDeviation; } DoseStatisticType DVH::getVariance() const { return _variance; } void DVH::init() { if (_deltaD == 0 || _deltaV == 0) { throw InvalidParameterException("DVH init error: neither _deltaD nor _deltaV must be zero!"); } if (this->_dataDifferential.empty()) { throw InvalidParameterException("DVH init error: data differential is empty!"); } double sum = 0; double squareSum = 0; _numberOfVoxels = 0; _maximum = 0; _minimum = 0; _dataCumulative.clear(); this->_dataCumulativeRelative.clear(); this->_dataDifferentialRelative.clear(); DataDifferentialType::iterator it; int i = 0; for (it = _dataDifferential.begin(); it != _dataDifferential.end(); ++it) { _numberOfVoxels += (*it); if ((*it) > 0) { _maximum = (i + 0.5) * this->_deltaD; } if ((_minimum == 0.0f) && ((*it) > 0)) { _minimum = (i + 0.5) * this->_deltaD; } sum += (*it) * (i + 0.5) * this->_deltaD; squareSum += (*it) * pow((i + 0.5) * this->_deltaD, 2); i++; } _mean = sum / _numberOfVoxels; for (it = _dataDifferential.begin(); it != _dataDifferential.end(); ++it) { DoseCalcType datai = ((*it) * 1.0 / _numberOfVoxels); _dataDifferentialRelative.push_back(datai); } _variance = (squareSum / _numberOfVoxels - _mean * _mean); _stdDeviation = pow(_variance, 0.5); _dataCumulative = this->calcCumulativeDVH(); } std::deque DVH::calcCumulativeDVH(bool relativeVolume) { _dataCumulative.clear(); _dataCumulativeRelative.clear(); DoseCalcType cumulativeDVHi = 0; size_t size = _dataDifferential.size(); for (size_t i = 0; i < size; i++) { cumulativeDVHi += _dataDifferential.at(size - i - 1); if (!relativeVolume) { _dataCumulative.push_front(cumulativeDVHi); } else { _dataCumulativeRelative.push_front(cumulativeDVHi / this->getNumberOfVoxels()); } } if (!relativeVolume) { return _dataCumulative; } else { return _dataCumulativeRelative; } } DoseStatisticType DVH::getMedian() const { double median_voxel = 0; int median_i = 0; for (GridIndexType i = 0; i < this->_dataDifferential.size(); i++) { if (median_voxel < (_numberOfVoxels - median_voxel)) { median_voxel += _dataDifferential[i]; median_i = i; } } double median = (median_i + 0.5) * this->_deltaD; return median; } DoseStatisticType DVH::getModal() const { double modal_voxel = 0; int modal_i = 0; for (GridIndexType i = 0; i < this->_dataDifferential.size(); i++) { if (modal_voxel < _dataDifferential[i]) { modal_voxel = _dataDifferential[i]; modal_i = i; } } double modal = (modal_i + 0.5) * this->_deltaD; return modal; } VolumeType DVH::getVx(DoseTypeGy xDoseAbsolute) { GridIndexType i = static_cast(xDoseAbsolute / _deltaD); if (i < _dataCumulative.size()) { VolumeType vx = (_dataCumulative.at(i)); vx = (vx * this->_deltaV); return vx; } else if (i < _dataCumulativeRelative.size()) { VolumeType vx = (_dataCumulativeRelative.at(i)); vx = (vx * this->_deltaV); return vx; } else { return 0; } } DoseTypeGy DVH::getDx(VolumeType xVolumeAbsolute) { GridIndexType i = 0; if (!_dataCumulative.empty()) { for (; i < _dataCumulative.size(); i++) { double volumeAbsoluteI = _dataCumulative[i] * this->_deltaV; if (xVolumeAbsolute > volumeAbsoluteI) { break; } } } else { for (; i < _dataCumulativeRelative.size(); i++) { double volumeAbsoluteI = _dataCumulativeRelative[i] * this->_deltaV; if (xVolumeAbsolute / this->getNumberOfVoxels() > volumeAbsoluteI) { break; } } } if (i <= _dataCumulative.size() && i > 0) { DoseTypeGy dx = (i - 1) * this->_deltaD; return dx; } else if (i < _dataCumulativeRelative.size() && i > 0) { DoseTypeGy dx = (i - 1) * this->_deltaD; return dx; } else { return 0; } } VolumeType DVH::getAbsoluteVolume(int relativePercent) { return (relativePercent * getNumberOfVoxels() * getDeltaV() / 100.0); } void DVH::setLabel(StructureLabel aLabel) { _label = aLabel; } StructureLabel DVH::getLabel() const { return _label; } + std::map DVH::getNormalizedDVH(bool differential) { + std::map normalizedDVH; + DataDifferentialType data = getDataDifferential(); + size_t numberOfVolumes = data.size(); + + if (!differential) { + data = calcCumulativeDVH(); + } + for (size_t i = 0; i < numberOfVolumes; i++) + { + normalizedDVH.insert(std::pair(i * getDeltaD(), data[i] * getDeltaV())); + } + return normalizedDVH; + } + }//end namespace core }//end namespace rttb diff --git a/code/core/rttbDVH.h b/code/core/rttbDVH.h index 3c3ac11..bedf9d8 100644 --- a/code/core/rttbDVH.h +++ b/code/core/rttbDVH.h @@ -1,195 +1,197 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __DVH_H #define __DVH_H #include #include +#include #include "boost/shared_ptr.hpp" #include "rttbBaseType.h" #include "RTTBCoreExports.h" namespace rttb { namespace core { /*! @class DVH @brief This is a class representing a dose volume histogram (DVH) */ class RTTBCore_EXPORT DVH { public: typedef std::deque DataDifferentialType; typedef boost::shared_ptr DVHPointer; private: /*! @brief Differential dvh data index is the dose bin, value is the voxel number (sub voxel accuracy) of the dose bin */ DataDifferentialType _dataDifferential; /*! @brief Differential dvh data relative to the total number of voxels */ DataDifferentialType _dataDifferentialRelative; /*! @brief Absolute dose value of a dose-bin in Gy */ DoseTypeGy _deltaD; /*! @brief Volume of a voxel in cm3 */ DoseVoxelVolumeType _deltaV; IDType _structureID; IDType _doseID; IDType _voxelizationID; StructureLabel _label; DoseStatisticType _maximum; DoseStatisticType _minimum; DoseStatisticType _mean; DoseStatisticType _modal; DVHVoxelNumber _numberOfVoxels; DoseStatisticType _median; DoseStatisticType _stdDeviation; DoseStatisticType _variance; DataDifferentialType _dataCumulative; DataDifferentialType _dataCumulativeRelative; /*! @brief DVH initialization The DVH is initialized and all statistical values are calculated. @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ void init(); public: ~DVH(); /*! @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ DVH(const DataDifferentialType& aDataDifferential, const DoseTypeGy& aDeltaD, const DoseVoxelVolumeType& aDeltaV, IDType aStructureID, IDType aDoseID); /*! @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ DVH(const DataDifferentialType& aDataDifferential, DoseTypeGy aDeltaD, DoseVoxelVolumeType aDeltaV, IDType aStructureID, IDType aDoseID, IDType aVoxelizationID); DVH(const DVH& copy); /*! @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ DVH& operator=(const DVH& copy); /*! equality operator DVHs are considered equal if they have the same structure, dose and voxelization ID and the same number of voxels. */ bool friend operator==(const DVH& aDVH, const DVH& otherDVH); friend std::ostream& operator<<(std::ostream& s, const DVH& aDVH); void setLabel(StructureLabel aLabel); StructureLabel getLabel() const; /*! @param relativeVolume default false-> Value is the voxel number of the dose bin; if true-> value is the relative volume % between 0 and 1, (the voxel number of this dose bin)/(number of voxels) @return Return differential data of the dvh (relative or absolute depending on the input parameter). */ DataDifferentialType getDataDifferential(bool relativeVolume = false) const; DoseVoxelVolumeType getDeltaV() const; DoseTypeGy getDeltaD() const; IDType getStructureID() const; IDType getDoseID() const; IDType getVoxelizationID() const; void setDoseID(IDType aDoseID); void setStructureID(IDType aStrID); /*! @brief Calculate number of the voxels (with sub voxel accuracy) @return Return -1 if not initialized */ DVHVoxelNumber getNumberOfVoxels() const; /*! @brief Get the maximum dose in Gy from dvh @return Return the maximum dose in Gy (i+0.5)*deltaD, i-the maximal dose-bin with volume>0 Return -1 if not initialized */ DoseStatisticType getMaximum() const; /*! @brief Get the minimum dose in Gy from dvh @return Return the minimum dose (i+0.5)*deltaD, i-the minimal dose-bin with volume>0 Return -1 if not initialized */ DoseStatisticType getMinimum() const; DoseStatisticType getMean() const; DoseStatisticType getMedian() const; DoseStatisticType getModal() const; DoseStatisticType getStdDeviation() const; DoseStatisticType getVariance() const; /*! @brief Calculate the cumulative data of dvh @param relativeVolume default false-> Value is the voxel number of the dose bin; if true-> value is the relative volume % between 0 and 1, (the voxel number of this dose bin)/(number of voxels) @return Return cumulative dvh value i-voxel number in dose-bin i */ DataDifferentialType calcCumulativeDVH(bool relativeVolume = false); /*! @brief Get Vx the volume irradiated to >= x @return Return absolute Volume in absolute cm3 Return -1 if not initialized */ VolumeType getVx(DoseTypeGy xDoseAbsolute); /*! @brief Get Dx the minimal dose delivered to x @return Return absolute dose value in Gy Return -1 if not initialized */ DoseTypeGy getDx(VolumeType xVolumeAbsolute); /*! @brief Calculate the absolute volume in cm3 @param relativePercent 0~100, the percent of the whole volume */ VolumeType getAbsoluteVolume(int relativePercent); + std::map getNormalizedDVH(bool differential = false); }; } } #endif diff --git a/testing/core/DVHTest.cpp b/testing/core/DVHTest.cpp index 1349013..feddfaa 100644 --- a/testing/core/DVHTest.cpp +++ b/testing/core/DVHTest.cpp @@ -1,191 +1,197 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVH.h" #include "DummyDoseAccessor.h" #include "DummyMaskAccessor.h" namespace rttb { namespace testing { typedef core::DVH::DataDifferentialType DataDifferentialType; /*! @brief DVHTest - test the API of DVH 1) test constructors (values as expected?) 2) test asignement 3) test set/getLabel 4) test set/get 5) test equality */ int DVHTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; //generate artificial DVH and corresponding statistical values DoseTypeGy binSize = DoseTypeGy(0.1); DoseVoxelVolumeType voxelVolume = 8; DataDifferentialType anEmptyDataDifferential; DataDifferentialType aDataDifferential; DataDifferentialType aDataDifferential2; DataDifferentialType aDataDifferentialRelative; DoseStatisticType maximum = 0; DoseStatisticType minimum = 0; double sum = 0; double squareSum = 0; DoseCalcType value = 0; DVHVoxelNumber numberOfVoxels = 0; // creat default values [0,100) for (int i = 0; i < 100; i++) { value = DoseCalcType((double(rand()) / RAND_MAX) * 1000); numberOfVoxels += value; aDataDifferential.push_back(value); aDataDifferential2.push_back(value * 10); if (value > 0) { maximum = (i + 0.5) * binSize; if (minimum == 0) { minimum = (i + 0.5) * binSize; } } sum += value * (i + 0.5) * binSize; squareSum += value * pow((i + 0.5) * binSize, 2); } DoseStatisticType mean = sum / numberOfVoxels; DoseStatisticType variance = (squareSum / numberOfVoxels - mean * mean); DoseStatisticType stdDeviation = pow(variance, 0.5); std::deque::iterator it; for (it = aDataDifferential.begin(); it != aDataDifferential.end(); ++it) { aDataDifferentialRelative.push_back((*it) / numberOfVoxels); } const IDType structureID = "myStructure"; const IDType doseID = "myDose"; const IDType voxelizationID = "myVoxelization"; //1) test default constructor (values as expected?) CHECK_THROW(core::DVH(anEmptyDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(anEmptyDataDifferential, binSize, voxelVolume, structureID, doseID)); CHECK_NO_THROW(core::DVH(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID)); CHECK_NO_THROW(core::DVH(aDataDifferential, binSize, voxelVolume, structureID, doseID)); CHECK_THROW(core::DVH(aDataDifferential, 0, voxelVolume, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(aDataDifferential, 0, voxelVolume, structureID, doseID)); CHECK_THROW(core::DVH(aDataDifferential, binSize, 0, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(aDataDifferential, binSize, 0, structureID, doseID)); //2) test asignement core::DVH myTestDVH(aDataDifferential, binSize, voxelVolume, structureID, doseID); CHECK_NO_THROW(core::DVH myOtherDVH = myTestDVH); const core::DVH myOtherDVH = myTestDVH; CHECK_NO_THROW(core::DVH aDVH(myOtherDVH)); //3) test set/getLabel core::DVH myDVH(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID); StructureLabel label = ""; CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); label = "myLabel"; CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); label = "myLabel2"; CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); //4) test set/get //IDs CHECK_EQUAL(myDVH.getStructureID(), structureID); CHECK_EQUAL(myDVH.getDoseID(), doseID); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); /*! is related to #2029*/ myDVH.setDoseID("somethingElse"); CHECK_EQUAL(myDVH.getDoseID(), "somethingElse"); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); CHECK_EQUAL(myDVH.getStructureID(), structureID); /*! is related to #2029*/ myDVH.setStructureID("somethingOther"); CHECK_EQUAL(myDVH.getDoseID(), "somethingElse"); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); CHECK_EQUAL(myDVH.getStructureID(), "somethingOther"); //dataDifferential CHECK(myDVH.getDataDifferential() == aDataDifferential); CHECK(myDVH.getDataDifferential(false) == aDataDifferential); CHECK(myDVH.getDataDifferential(true) == aDataDifferentialRelative); CHECK_EQUAL(myDVH.getNumberOfVoxels(), numberOfVoxels); CHECK_EQUAL(myDVH.getDeltaV(), voxelVolume); CHECK_EQUAL(myDVH.getDeltaD(), binSize); CHECK_EQUAL(myDVH.getMaximum(), maximum); CHECK_EQUAL(myDVH.getMinimum(), minimum); CHECK_EQUAL(myDVH.getMean(), mean); CHECK_EQUAL(myDVH.getVariance(), variance); CHECK_EQUAL(myDVH.getStdDeviation(), stdDeviation); int percentage = 20; VolumeType absVol = VolumeType(percentage * numberOfVoxels * voxelVolume / 100.0); CHECK_EQUAL(myDVH.getAbsoluteVolume(percentage), absVol); //5) test equality core::DVH myDVH2(aDataDifferential2, binSize, voxelVolume, structureID, doseID); CHECK(!(myDVH == myDVH2)); CHECK_EQUAL(myDVH, myDVH); core::DVH aDVH(myOtherDVH); CHECK_EQUAL(aDVH, myOtherDVH); + std::map normalizedDVH = myDVH.getNormalizedDVH(); + for (auto elem : normalizedDVH) + { + CHECK_EQUAL(aDataDifferential.at(round(elem.first / binSize)), (elem.second / voxelVolume)); //or +0,5 + } + RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb