diff --git a/code/core/rttbBaseType.h b/code/core/rttbBaseType.h index b338a3e..643a4f8 100644 --- a/code/core/rttbBaseType.h +++ b/code/core/rttbBaseType.h @@ -1,742 +1,742 @@ // ----------------------------------------------------------------------- // 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 #include #include namespace rttb { const double errorConstant = 1e-5; 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 (int i = 0; i < gi1.size(); i++) + for (unsigned int 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 (int i = 0; i < wc1.size(); i++) + for (unsigned int i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } }; /*! @class WorldCoordinate3D @brief 3D coordinate in real world coordinates. */ 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 (int i = 0; i < wc1.size(); i++) + for (unsigned int i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } bool equalsAlmost(const WorldCoordinate3D& another, double errorConstant = 1e-5) const { if (size() != another.size()) { return false; } double dist = norm_2(*this - another); return dist < errorConstant; } friend std::ostream& operator<<(std::ostream& s, const WorldCoordinate3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; 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 (int i = 0; i < wc1.size(); i++) + for (unsigned int i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } bool equalsAlmost(const SpacingVectorType3D& another, double errorConstant) const { if ((*this).size() != another.size()) { return false; } double dist = norm_2(*this - another); return dist < errorConstant; } 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 errorConstant) 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 ((abs((*this)(m, n) - anOrientationMatrix(m, n)) > errorConstant)) { 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. */ 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 (int i = 0; i < gi1.size(); i++) + for (unsigned int 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 (int i = 0; i < gi1.size(); i++) + for (unsigned int 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 int VoxelGridDimensionType; typedef boost::numeric::ublas::vector VoxelGridDimensionVectorType; typedef double FractionType, VoxelSizeType, DVHVoxelNumber; typedef double DoseCalcType, DoseTypeGy, GenericValueType, DoseVoxelVolumeType, VolumeType, GridVolumeType, 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 VirtuosFileNameString, 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 74c0ac0..ecc4228 100644 --- a/code/core/rttbDVH.cpp +++ b/code/core/rttbDVH.cpp @@ -1,418 +1,418 @@ // ----------------------------------------------------------------------- // 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() << ", " << aDVH.getNumberOfVoxels() << " ]"; 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 (int i = 0; i < size; i++) + for (unsigned int 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; } }//end namespace core }//end namespace rttb diff --git a/code/core/rttbDVH.h b/code/core/rttbDVH.h index 5c49c66..f9d9a82 100644 --- a/code/core/rttbDVH.h +++ b/code/core/rttbDVH.h @@ -1,191 +1,194 @@ // ----------------------------------------------------------------------- // 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 "rttbBaseType.h" #include "rttbStructure.h" namespace rttb { namespace core { class Structure; /*! @class DVH @brief This is a class representing a dose volume histogram (DVH) */ class DVH { public: typedef std::deque DataDifferentialType; typedef boost::shared_ptr DVHPointer; private: - IDType _structureID; - IDType _doseID; - IDType _voxelizationID; /*! @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 Volume of a voxel in cm3 - */ - DoseVoxelVolumeType _deltaV; /*! @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 initialisation 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 abolute 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); }; } } #endif diff --git a/code/core/rttbDVHSet.cpp b/code/core/rttbDVHSet.cpp index 2410755..4db072a 100644 --- a/code/core/rttbDVHSet.cpp +++ b/code/core/rttbDVHSet.cpp @@ -1,223 +1,223 @@ // ----------------------------------------------------------------------- // 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 "rttbDVHSet.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace core { DVHSet::DVHSet(IDType aStructureSetID, IDType aDoseID) : _structureSetID(aStructureSetID), _doseID(aDoseID) {} DVHSet::DVHSet(DVHSetType aDVHTVSet, DVHSetType aDVHHTSet, IDType aStructureSetID, IDType aDoseID) { _dvhHTSet = aDVHHTSet; _dvhTVSet = aDVHTVSet; this->_structureSetID = aStructureSetID; this->_doseID = aDoseID; } DVHSet::DVHSet(DVHSetType aDVHTVSet, DVHSetType aDVHHTSet, DVHSetType aDVHWVSet, IDType aStructureSetID, IDType aDoseID) { _dvhHTSet = aDVHHTSet; _dvhTVSet = aDVHTVSet; _dvhWVSet = aDVHWVSet; this->_structureSetID = aStructureSetID; this->_doseID = aDoseID; } std::size_t DVHSet::size() const { return _dvhHTSet.size() + _dvhTVSet.size() + _dvhWVSet.size(); } void DVHSet::setStrSetID(IDType aStrSetID) { _structureSetID = aStrSetID; } void DVHSet::setDoseID(IDType aDoseID) { _doseID = aDoseID; } IDType DVHSet::getStrSetID() const { return _structureSetID; } IDType DVHSet::getDoseID() const { return _doseID; } DVH* DVHSet::getDVH(IDType structureID) { DVHSetType::iterator itTV = _dvhTVSet.begin(); for (; itTV != _dvhTVSet.end(); ++itTV) { if ((*(itTV)).getStructureID() == structureID) { return &(*itTV); } } DVHSetType::iterator itHT = _dvhHTSet.begin(); for (; itHT != _dvhHTSet.end(); ++itHT) { if ((*(itHT)).getStructureID() == structureID) { return &(*itHT); } } DVHSetType::iterator itWV = _dvhWVSet.begin(); for (; itWV != _dvhWVSet.end(); ++itWV) { if ((*(itWV)).getStructureID() == structureID) { return &(*itWV); } } std::cout << "No DVH with the structure id: " << structureID << " was found!" << std::endl; return NULL; } void DVHSet::insert(DVH& aDvh, DVHRole aDVHRole) { if (aDVHRole.Type == DVHRole::TargetVolume) { _dvhTVSet.push_back(aDvh); } else if (aDVHRole.Type == DVHRole::HealthyTissue) { _dvhHTSet.push_back(aDvh); } else if (aDVHRole.Type == DVHRole::WholeVolume) { _dvhWVSet.push_back(aDvh); } else { throw core::InvalidParameterException("aDVHType must be TV or HT or WV!"); } } VolumeType DVHSet::getWholeVolume(DoseTypeGy aDoseAbsolute) const { VolumeType volume = this->getHealthyTissueVolume(aDoseAbsolute) + this->getTargetVolume( aDoseAbsolute); return volume; } VolumeType DVHSet::getHealthyTissueVolume(DoseTypeGy aDoseAbsolute) const { DVHSetType::const_iterator itHT = _dvhHTSet.begin(); VolumeType volume = 0; while (itHT != _dvhHTSet.end()) { VolumeType testVol = 0; DVH dvh = *(itHT); testVol = dvh.getVx(aDoseAbsolute); if (testVol >= 0) { volume += testVol; } ++itHT; } return volume; } VolumeType DVHSet::getTargetVolume(DoseTypeGy aDoseAbsolute) const { DVHSetType::const_iterator itTV = _dvhTVSet.begin(); VolumeType volume = 0; while (itTV != _dvhTVSet.end()) { VolumeType testVol = 0; DVH dvh = *(itTV); testVol = dvh.getVx(aDoseAbsolute); if (testVol >= 0) { volume += testVol; } ++itTV; } return volume; } bool operator==(const DVHSet& aDVHSet, const DVHSet& otherDVHSet) { if (aDVHSet.getStrSetID() != otherDVHSet.getStrSetID()) { return false; } if (aDVHSet.getDoseID() != otherDVHSet.getDoseID()) { return false; } if (aDVHSet.size() != otherDVHSet.size()) { return false; } return true; } std::ostream& operator<<(std::ostream& s, const DVHSet& aDVHSet) { s << "[ " << aDVHSet.getStrSetID() << ", " << aDVHSet.getDoseID() << " ]"; return s; } std::ostream& operator<<(std::ostream& s, const DVHSet::DVHSetType& aDVHSet) { s << "[ "; - for (int i = 0; i < aDVHSet.size(); i++) + for (unsigned int i = 0; i < aDVHSet.size(); i++) { s << aDVHSet.at(i); } s << " ]"; return s; } } } diff --git a/code/core/rttbGeometricInfo.h b/code/core/rttbGeometricInfo.h index f802ce7..b3df594 100644 --- a/code/core/rttbGeometricInfo.h +++ b/code/core/rttbGeometricInfo.h @@ -1,190 +1,190 @@ // ----------------------------------------------------------------------- // 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 __GEOMETRIC_INFO_NEW_H #define __GEOMETRIC_INFO_NEW_H #include #include #include #include #include "rttbBaseType.h" namespace rttb { namespace core { /*! @brief GeometricInfo objects contain all the information required for transformations between voxel grid coordinates and world coordinates. Corresponding converter functions are also available. @note ITK Pixel Indexing used (http://www.itk.org/Doxygen45/html/classitk_1_1Image.html): The Index type reverses the order so that with Index[0] = col, Index[1] = row, Index[2] = slice. */ class GeometricInfo { private: WorldCoordinate3D _imagePositionPatient; OrientationMatrix _orientationMatrix; OrientationMatrix _invertedOrientationMatrix; SpacingVectorType3D _spacing; VoxelGridDimensionType _numberOfColumns; VoxelGridDimensionType _numberOfRows; VoxelGridDimensionType _numberOfFrames; /* @brief Matrix inversion routine. Uses lu_factorize and lu_substitute in uBLAS to invert a matrix http://savingyoutime.wordpress.com/2009/09/21/c-matrix-inversion-boostublas/ */ bool computeInvertOrientation(); public: /*! @brief Constructor, initializes orientation matrix, spacing vector and patient position with zeros. */ - GeometricInfo() : _orientationMatrix(0), _spacing(0), _imagePositionPatient(0), _numberOfFrames(0), - _numberOfRows(0), - _numberOfColumns(0) {} + GeometricInfo() : _imagePositionPatient(0),_orientationMatrix(0), _spacing(0), _numberOfColumns(0), + _numberOfRows(0), _numberOfFrames(0) + {} void setSpacing(const SpacingVectorType3D& aSpacingVector); const SpacingVectorType3D& getSpacing() const; void setImagePositionPatient(const WorldCoordinate3D& aImagePositionPatient); const WorldCoordinate3D& getImagePositionPatient() const; void setOrientationMatrix(const OrientationMatrix& anOrientationMatrix); const OrientationMatrix getOrientationMatrix() const { return _orientationMatrix; }; /*! @brief Returns the 1st row of the OrientationMatrix. @deprecated please use getOrientationMatrix() (x,0) instead*/ const ImageOrientation getImageOrientationRow() const; /*! @brief Returns the 2nd row of the OrientationMatrix. @deprecated please use getOrientationMatrix() (x,1) instead*/ const ImageOrientation getImageOrientationColumn() const; void setPixelSpacingRow(const GridVolumeType aValue); /*! @brief Returns the spacing in the x-dimension (rows) of the data grid. @deprecated please use getSpacing() (0) instead*/ const GridVolumeType getPixelSpacingRow() const; void setPixelSpacingColumn(const GridVolumeType aValue); /*! @brief Returns the spacing in the y-dimension (columns) of the data grid. @deprecated please use getSpacing() (1) instead*/ const GridVolumeType getPixelSpacingColumn() const; void setSliceThickness(const GridVolumeType aValue); /*! @brief Returns the spacing in the z-dimension (slices) of the data grid. @deprecated please use getSpacing() (2) instead*/ const GridVolumeType getSliceThickness() const; void setImageSize(const ImageSize& aSize); const ImageSize getImageSize() const; void setNumColumns(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumColumns() const; void setNumRows(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumRows() const; void setNumSlices(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumSlices() const; /*! @brief determines equality of two GeometricInfo objects. */ friend bool operator==(const GeometricInfo& gInfo, const GeometricInfo& gInfo1); bool equalsAlmost(const GeometricInfo& another, double errorConstant = 1e-5) const; /*! @brief convert world coordinates to voxel grid index. The conversion of values is done even if the target index is not inside the given voxel grid (return false). If the target is inside the grid return true. */ bool worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, VoxelGridIndex3D& aIndex) const; /*! @brief convert world coordinates to double geometry coordinate. The world coordinate of the image position patient (center of the first voxel) will be convert to the double voxel (0.0, 0.0, 0.0) The conversion of values is done even if the target index is not inside the given voxel grid (return false). If the target is inside the grid return true. */ bool worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, DoubleVoxelGridIndex3D& aIndex) const; /*! @brief convert double geometry coordinate to world coordinates. The double voxel index (0.0, 0.0, 0.0) will be convert to the world coordinate of the image postion patient (center of the first voxel) The conversion of values is done even if the target is not inside the given voxel grid (return false). If the target is inside the voxel grid return true. */ bool geometryCoordinateToWorldCoordinate(const DoubleVoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const; /*! @brief convert int voxel grid index to world coordinates. The conversion of values is done even if the target is not inside the given voxel grid (return false). If the target is inside the voxel grid return true. */ bool indexToWorldCoordinate(const VoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const; /*! @brief check if a given voxel grid index is inside the given voxel grid.*/ bool isInside(const VoxelGridIndex3D& aIndex) const; /*! @brief check if a given world coordinate is inside the given voxel grid.*/ bool isInside(const WorldCoordinate3D& aWorldCoordinate); const GridSizeType getNumberOfVoxels() const; bool convert(const VoxelGridID& gridID, VoxelGridIndex3D& gridIndex) const; bool convert(const VoxelGridIndex3D& gridIndex, VoxelGridID& gridID) const; /*! @brief test if given ID is inside current dose grid */ bool validID(const VoxelGridID aID) const; /*! @brief test if given index is inside current dose grid */ bool validIndex(const VoxelGridIndex3D& aIndex) const; /*!@ brief generates string stream representation of the GeometricInfo object. */ friend std::ostream& operator<<(std::ostream& s, const GeometricInfo& anGeometricInfo); }; } } -#endif \ No newline at end of file +#endif diff --git a/code/indices/rttbCoverageIndex.cpp b/code/indices/rttbCoverageIndex.cpp index 426086f..f79ec24 100644 --- a/code/indices/rttbCoverageIndex.cpp +++ b/code/indices/rttbCoverageIndex.cpp @@ -1,80 +1,79 @@ // ----------------------------------------------------------------------- // 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 "rttbCoverageIndex.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace indices { CoverageIndex::CoverageIndex(DVHSetPtr dvhSet, DoseTypeGy aDoseReference) : DvhBasedDoseIndex(dvhSet, aDoseReference) { init(); } bool CoverageIndex::calcIndex() { VolumeType TV = _dvhSet->getTargetVolume(0); if (TV != 0) { _value = _dvhSet->getTargetVolume(this->_doseReference) / TV; } else { throw core::InvalidParameterException("DVH Set invalid: Target volume should not be 0!"); } return true; } IndexValueType CoverageIndex::getValueAt(core::DVHSet::IndexType tvIndex) { std::vector dvhTVSet = this->_dvhSet->getTargetVolumeSet(); - VolumeType Vref = _dvhSet->getWholeVolume(_doseReference); if (tvIndex >= dvhTVSet.size()) { rttbExceptionMacro(core::InvalidParameterException, << "tvIndex invalid: it should be <" << dvhTVSet.size() << "!"); } core::DVH dvh = dvhTVSet.at(tvIndex); VolumeType TV = dvh.getVx(0); if (TV == 0) { throw core::InvalidParameterException("DVH invalid: Volume of tvIndex should not be 0!"); } IndexValueType value = dvh.getVx(_doseReference) / TV; //the irradiation factor of i-th treated volume return value; } }//end namespace indices }//end namespace rttb diff --git a/code/indices/rttbHomogeneityIndex.cpp b/code/indices/rttbHomogeneityIndex.cpp index 556c24d..fd6ec79 100644 --- a/code/indices/rttbHomogeneityIndex.cpp +++ b/code/indices/rttbHomogeneityIndex.cpp @@ -1,102 +1,103 @@ // ----------------------------------------------------------------------- // 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 "rttbHomogeneityIndex.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace indices { HomogeneityIndex::HomogeneityIndex(DVHSetPtr dvhSet, DoseTypeGy aDoseReference) : DvhBasedDoseIndex(dvhSet, aDoseReference) { init(); } bool HomogeneityIndex::calcIndex() { double max = 0; - double min; + double min = std::numeric_limits::max(); std::vector dvhTVSet = this->_dvhSet->getTargetVolumeSet(); std::vector::iterator it; for (it = dvhTVSet.begin(); it != dvhTVSet.end(); ++it) { core::DVH dvh = *(it); if (it == dvhTVSet.begin()) { min = dvh.getMinimum(); } if (dvh.getMaximum() > max) { max = dvh.getMaximum(); } if (dvh.getMinimum() < min) { min = dvh.getMinimum(); } } if (this->getDoseReference() != 0) { _value = (max - min) / this->getDoseReference(); } else { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } return true; } IndexValueType HomogeneityIndex::getValueAt(core::DVHSet::IndexType tvIndex) { std::vector dvhTVSet = this->_dvhSet->getTargetVolumeSet(); if (tvIndex >= dvhTVSet.size()) { rttbExceptionMacro(core::InvalidParameterException, << "tvIndex invalid: it should be <" << dvhTVSet.size() << "!"); } core::DVH dvh = dvhTVSet.at(tvIndex); if (this->getDoseReference() <= 0) { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } return (dvh.getMaximum() - dvh.getMinimum()) / this->getDoseReference(); } }//end namespace indices -}//end namespace rttb \ No newline at end of file +}//end namespace rttb diff --git a/code/io/helax/rttbDicomHelaxDoseAccessor.cpp b/code/io/helax/rttbDicomHelaxDoseAccessor.cpp index 850c9c0..af46a09 100644 --- a/code/io/helax/rttbDicomHelaxDoseAccessor.cpp +++ b/code/io/helax/rttbDicomHelaxDoseAccessor.cpp @@ -1,328 +1,328 @@ // ----------------------------------------------------------------------- // 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 "boost/filesystem/operations.hpp" #include "boost/filesystem/path.hpp" #include "boost/progress.hpp" #include #include "rttbDicomHelaxDoseAccessor.h" #include "rttbInvalidDoseException.h" #include "rttbDcmrtException.h" #include "rttbIndexOutOfBoundsException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace io { namespace helax { DicomHelaxDoseAccessor::~DicomHelaxDoseAccessor() { } DicomHelaxDoseAccessor::DicomHelaxDoseAccessor(std::vector aDICOMRTDoseVector) { - for (int i = 0; i < aDICOMRTDoseVector.size(); i++) + for (size_t i = 0; i < aDICOMRTDoseVector.size(); i++) { _doseVector.push_back(aDICOMRTDoseVector.at(i)); } this->begin(); } bool DicomHelaxDoseAccessor::begin() { if (_doseVector.size() == 0) { throw core::InvalidParameterException(" The size of aDICOMRTDoseVector is 0!"); } assembleGeometricInfo(); _doseData.clear(); OFString doseGridScalingStr; _doseVector.at(0)->getDoseGridScaling( doseGridScalingStr);//get the first dose grid scaling as _doseGridScaling try { _doseGridScaling = boost::lexical_cast(doseGridScalingStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Dose grid scaling not readable or = 0!") ; } - for (int i = 0; i < _doseVector.size(); i++) + for (size_t i = 0; i < _doseVector.size(); i++) { DRTDoseIODPtr dose = _doseVector.at(i); OFString currentDoseGridScalingStr; dose->getDoseGridScaling(currentDoseGridScalingStr); double currentDoseGridScaling; try { currentDoseGridScaling = boost::lexical_cast(currentDoseGridScalingStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Dose grid scaling not readable or = 0!") ; } OFCondition status; DcmFileFormat fileformat; DcmItem doseitem; status = dose->write(doseitem); if (status.good()) { unsigned long count; const Uint16* pixelData; status = doseitem.findAndGetUint16Array(DcmTagKey(0x7fe0, 0x0010), pixelData, &count); if (status.good()) { - for (unsigned int i = 0; - i < static_cast(_geoInfo.getNumColumns()*_geoInfo.getNumRows()); i++) + for (unsigned int j = 0; + j < static_cast(_geoInfo.getNumColumns()*_geoInfo.getNumRows()); j++) { - this->_doseData.push_back(pixelData[i]*currentDoseGridScaling / + this->_doseData.push_back(pixelData[j]*currentDoseGridScaling / _doseGridScaling); //recalculate dose data } } else { throw dicom::DcmrtException("Read Pixel Data (7FE0,0010) failed!"); } } else { throw dicom::DcmrtException("Read DICOM-RT Dose file failed!"); } } return true; } bool DicomHelaxDoseAccessor::assembleGeometricInfo() { DRTDoseIODPtr dose = _doseVector.at(0); Uint16 temp = 0; dose->getColumns(temp); _geoInfo.setNumColumns(temp); temp = 0; dose->getRows(temp); _geoInfo.setNumRows(temp); OFString numberOfFramesStr; dose->getNumberOfFrames(numberOfFramesStr); if (!numberOfFramesStr.empty()) { _geoInfo.setNumSlices(boost::lexical_cast(numberOfFramesStr.c_str())); } else { _geoInfo.setNumSlices((VoxelGridDimensionType)_doseVector.size()); } if (_geoInfo.getNumColumns() == 0 || _geoInfo.getNumRows() == 0 || _geoInfo.getNumSlices() == 0) { throw core::InvalidDoseException("Empty dicom dose!") ; } OFString imageOrientationRowX; dose->getImageOrientationPatient(imageOrientationRowX, 0); OFString imageOrientationRowY; dose->getImageOrientationPatient(imageOrientationRowY, 1); OFString imageOrientationRowZ; dose->getImageOrientationPatient(imageOrientationRowZ, 2); OFString imageOrientationColumnX; dose->getImageOrientationPatient(imageOrientationColumnX, 3); OFString imageOrientationColumnY; dose->getImageOrientationPatient(imageOrientationColumnY, 4); OFString imageOrientationColumnZ; dose->getImageOrientationPatient(imageOrientationColumnZ, 5); WorldCoordinate3D imageOrientationRow; WorldCoordinate3D imageOrientationColumn; try { imageOrientationRow(0) = boost::lexical_cast(imageOrientationRowX.c_str()); imageOrientationRow(1) = boost::lexical_cast(imageOrientationRowY.c_str()); imageOrientationRow(2) = boost::lexical_cast(imageOrientationRowZ.c_str()); imageOrientationColumn(0) = boost::lexical_cast(imageOrientationColumnX.c_str()); imageOrientationColumn(1) = boost::lexical_cast(imageOrientationColumnY.c_str()); imageOrientationColumn(2) = boost::lexical_cast(imageOrientationColumnZ.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("boost::lexical_cast ImageOrientation failed! Can not read image orientation X/Y/Z!") ; } OrientationMatrix orientation; orientation(0, 0) = imageOrientationRow.x(); orientation(1, 0) = imageOrientationRow.y(); orientation(2, 0) = imageOrientationRow.z(); orientation(0, 1) = imageOrientationColumn.x(); orientation(1, 1) = imageOrientationColumn.y(); orientation(2, 1) = imageOrientationColumn.z(); WorldCoordinate3D perpendicular = imageOrientationRow.cross(imageOrientationColumn); orientation(0, 2) = perpendicular.x(); orientation(1, 2) = perpendicular.y(); orientation(2, 2) = perpendicular.z(); _geoInfo.setOrientationMatrix(orientation); OFString imagePositionX; dose->getImagePositionPatient(imagePositionX, 0); OFString imagePositionY; dose->getImagePositionPatient(imagePositionY, 1); OFString imagePositionZ; dose->getImagePositionPatient(imagePositionZ, 2); WorldCoordinate3D imagePositionPatient; try { imagePositionPatient(0) = boost::lexical_cast(imagePositionX.c_str()); imagePositionPatient(1) = boost::lexical_cast(imagePositionY.c_str()); imagePositionPatient(2) = boost::lexical_cast(imagePositionZ.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("boost::lexical_cast ImagePosition failed! Can not read image position X/Y/Z!") ; } _geoInfo.setImagePositionPatient(imagePositionPatient); SpacingVectorType3D spacingVector; OFString pixelSpacingRowStr; dose->getPixelSpacing(pixelSpacingRowStr, 0); OFString pixelSpacingColumnStr; dose->getPixelSpacing(pixelSpacingColumnStr, 1); try { spacingVector(1) = boost::lexical_cast(pixelSpacingRowStr.c_str()); spacingVector(0) = boost::lexical_cast(pixelSpacingColumnStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Can not read Pixel Spacing Row/Column!") ; } _geoInfo.setSpacing(spacingVector); if (_geoInfo.getPixelSpacingRow() == 0 || _geoInfo.getPixelSpacingColumn() == 0) { throw core::InvalidDoseException("Pixel spacing not readable or = 0!"); } OFString sliceThicknessStr; dose->getSliceThickness(sliceThicknessStr); try { spacingVector(2) = boost::lexical_cast(sliceThicknessStr.c_str()); } catch (boost::bad_lexical_cast&) { spacingVector(2) = 0 ;//if no information about slice thickness, set to 0 and calculate it using z coordinate difference between 1. and 2. dose } if (spacingVector(2) == 0) { if (_doseVector.size() > 1) { DRTDoseIODPtr dose2 = _doseVector.at(1);//get the 2. dose OFString imagePositionZ2; dose2->getImagePositionPatient(imagePositionZ2, 2); try { spacingVector(2) = boost::lexical_cast(imagePositionZ2.c_str()) - imagePositionPatient( 2); //caculate slicethickness } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Can not read image position Z of the 2. dose!"); } } else { std::cerr << "sliceThickness == 0! It will be replaced with pixelSpacingRow=" << _geoInfo.getPixelSpacingRow() << "!" << std::endl; spacingVector(2) = spacingVector(0); } } _geoInfo.setSpacing(spacingVector); return true; } GenericValueType DicomHelaxDoseAccessor::getValueAt(const VoxelGridID aID) const { return _doseData.at(aID) * _doseGridScaling; } GenericValueType DicomHelaxDoseAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getValueAt(aVoxelGridID); } else { return -1; } } } } } diff --git a/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp b/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp index 228dc8d..1882187 100644 --- a/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp +++ b/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp @@ -1,94 +1,94 @@ // ----------------------------------------------------------------------- // 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 "boost/filesystem/operations.hpp" #include "boost/filesystem/path.hpp" #include "boost/progress.hpp" #include "rttbDicomHelaxFileDoseAccessorGenerator.h" #include "rttbDicomHelaxDoseAccessor.h" #include "rttbInvalidDoseException.h" #include "rttbDcmrtException.h" #include "rttbIndexOutOfBoundsException.h" #include "rttbInvalidParameterException.h" #include "rttbDicomFileReaderHelper.h" namespace rttb { namespace io { namespace helax { DicomHelaxFileDoseAccessorGenerator::~DicomHelaxFileDoseAccessorGenerator() {} DicomHelaxFileDoseAccessorGenerator::DicomHelaxFileDoseAccessorGenerator( FileNameType aDICOMRTDoseDirName) { _doseDirName = aDICOMRTDoseDirName; } core::DoseAccessorGeneratorInterface::DoseAccessorPointer DicomHelaxFileDoseAccessorGenerator::generateDoseAccessor() { rttb::io::dicom::Modality doseModality = {rttb::io::dicom::Modality::RTDOSE}; std::vector fileVector = rttb::io::dicom::getFileNamesWithSameUID(_doseDirName, doseModality); OFCondition status; DcmFileFormat fileformat; std::vector doseVector; - for (int i = 0; i < fileVector.size(); i++) + for (size_t i = 0; i < fileVector.size(); i++) { DRTDoseIODPtr dose = boost::make_shared(); status = fileformat.loadFile(fileVector.at(i).c_str()); if (!status.good()) { throw core::InvalidDoseException("Error: load dose fileformat.loadFile failed!"); } status = dose->read(*fileformat.getDataset()); if (!status.good()) { throw core::InvalidDoseException("Error: read DRTDoseIOD failed!"); } doseVector.push_back(dose); } _doseAccessor = boost::make_shared(doseVector); return _doseAccessor; } } } } diff --git a/code/masks/boost/rttbBoostMaskAccessor.h b/code/masks/boost/rttbBoostMaskAccessor.h index 96f1f42..96cc8d6 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/boost/rttbBoostMaskAccessor.h @@ -1,128 +1,128 @@ // ----------------------------------------------------------------------- // 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 __BOOST_MASK_ACCESSOR__H #define __BOOST_MASK_ACCESSOR__H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include "rttbGenericDoseIterator.h" #include "rttbStructure.h" #include namespace rttb { namespace masks { namespace boost { /*! @class BoostMaskAccessor * @brief Using the voxelization based on boost::geometry and generate the mask accessor. * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. * (A structure without self interseciton means all contours of the structure have no self intersection, and * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored * and the mask will be calculated, however, it may cause errors in the mask results. */ class BoostMaskAccessor: public core::MaskAccessorInterface { public: typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::Structure::StructTypePointer StructTypePointer; private: - + StructTypePointer _spStructure; core::GeometricInfo _geoInfo; /*! vector containing list of mask voxels*/ MaskVoxelListPointer _spRelevantVoxelVector; - StructTypePointer _spStructure; + IDType _maskUID; bool _strict; public: /*! @brief Constructor with a structure pointer and a geometric info pointer * @param aStructurePointer smart pointer of the structure * @param aGeometricInfoPtr smart pointer of the geometricinfo of the dose * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. * @exception InvalidParameterException thrown if strict is true and the structure has self intersections */ BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, bool strict = true); /*! @brief destructor*/ ~BoostMaskAccessor(); /*! @brief voxelization of the given structures using boost algorithms*/ void updateMask(); /*! @brief get vector containing all relevant voxels that are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(); /*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold); /*!@brief determine how a given voxel on the dose grid is masked * @param aID ID of the voxel in grid. * @param voxel Reference to the voxel. * @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined. * The relevant volume fraction will be set to zero. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const; /*!@brief determine how a given voxel on the dose grid is masked * @param aIndex 3d index of the voxel in grid. * @param voxel Reference to the voxel. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const; /*! @brief give access to GeometricInfo*/ const core::GeometricInfo& getGeometricInfo() const; /* @ brief is true if dose is on a homogeneous grid * @remark Inhomogeneous grids are not supported at the moment, but if they will be supported in the future the interface does not need to change.*/ bool isGridHomogeneous() const { return true; }; IDType getMaskUID() const { return _maskUID; }; }; } } } #endif diff --git a/code/masks/rttbVOIindexIdentifier.cpp b/code/masks/rttbVOIindexIdentifier.cpp index 5b9f874..be15378 100644 --- a/code/masks/rttbVOIindexIdentifier.cpp +++ b/code/masks/rttbVOIindexIdentifier.cpp @@ -1,103 +1,103 @@ // ----------------------------------------------------------------------- // 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: // @date $Date: // @author $Author: */ #include "rttbVOIindexIdentifier.h" #include #include "rttbExceptionMacros.h" namespace rttb { namespace masks { const std::vector VOIindexIdentifier::getIndicesByVoiRegex( StructSetTypePointer spStructSet, const std::string& nameAsRegEx) { if (!spStructSet) { rttbDefaultExceptionStaticMacro("spStructSet is NULL"); } std::vector voiLabelList; std::vector resultVOiIndices; boost::regex e(nameAsRegEx); - for (int i = 0; i < spStructSet->getNumberOfStructures(); i++) + for (size_t i = 0; i < spStructSet->getNumberOfStructures(); i++) { voiLabelList.push_back(spStructSet->getStructure(i)->getLabel()); std::string s = spStructSet->getStructure(i)->getLabel(); if (boost::regex_match(s, e)) { resultVOiIndices.push_back(i); } } return resultVOiIndices; } const unsigned int VOIindexIdentifier::getIndexByVoiName(StructSetTypePointer spStructSet, const std::string& name) { if (!spStructSet) { rttbDefaultExceptionStaticMacro("spStructSet is NULL"); } for (unsigned int i = 0; i < spStructSet->getNumberOfStructures(); i++) { if (spStructSet->getStructure(i)->getLabel() == name) { return i; } } rttbDefaultExceptionStaticMacro("no VOI was found with the given name"); } const std::string VOIindexIdentifier::getVoiNameByIndex(StructSetTypePointer spStructSet, unsigned int index) { if (!spStructSet) { rttbDefaultExceptionStaticMacro("spStructSet is NULL!"); } if (index >= spStructSet->getNumberOfStructures()) { rttbDefaultExceptionStaticMacro("invalid index, voiLabelList out of range"); } return spStructSet->getStructure(index)->getLabel(); } } } diff --git a/code/models/rttbIntegration.cpp b/code/models/rttbIntegration.cpp index 52e780f..d811c3f 100644 --- a/code/models/rttbIntegration.cpp +++ b/code/models/rttbIntegration.cpp @@ -1,192 +1,192 @@ // ----------------------------------------------------------------------- // 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 "rttbIntegration.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace models { double tcpModelFunctor::calculate(double x) const { if (x == 0) { x = 1e-30; } return tcpFunction(a + (1 - x) / x, params) / (x * x); } double LkbModelFunctor::calculate(double x) const { if (x == 0) { x = 1e-30; } return lkbFunction(b - (1 - x) / x) / (x * x); } double tcpFunction(double x, const TcpParams& tcp_params) { double alphaVariance = tcp_params.alphaVariance; if (alphaVariance == 0) { alphaVariance = 1e-30; } double f = exp(-pow((x - tcp_params.alphaMean) / alphaVariance, 2) / 2); - for (int i = 0; i < tcp_params.volumeVector.size(); ++i) + for (unsigned int i = 0; i < tcp_params.volumeVector.size(); ++i) { double tmp, tmp1, tmp2, tmp3; tmp1 = exp(-x * tcp_params.bedVector.at(i)); tmp2 = -(tcp_params.rho) * tcp_params.volumeVector.at(i); tmp3 = tmp2 * tmp1; tmp = exp(tmp3); if (tmp != 1) { f = f * tmp; } } return f; } double integrateTCP(double a, const TcpParams& params) { double aNew = 1e-30; double bNew = 1.0; tcpModelFunctor BMFunction; BMFunction.params = params; BMFunction.a = a; return iterativeIntegration(BMFunction, aNew, bNew); } double lkbFunction(double x) { double tmp = -pow(x, 2) / 2; double step = exp(tmp); return step; } double integrateLKB(double b) { double aNew = 1e-30; double bNew = 1.0; LkbModelFunctor BMFunction; BMFunction.b = b; return iterativeIntegration(BMFunction, aNew, bNew); } //returns the nth stage of refinement of the extended trapezoidal rule template integrationType trapzd(const FunctorType& BMfunction, integrationType a, integrationType b, int stepNum) { static integrationType result; if (stepNum == 1) { result = 0.5 * (b - a) * (BMfunction.calculate(a) + BMfunction.calculate(b)); } else { integrationType x, tnm, sum, del; int it, j; for (it = 1, j = 1; j < stepNum - 1; j++) { it <<= 1; } tnm = it; del = (b - a) / tnm; x = a + 0.5 * del; for (sum = 0.0, j = 0; j < it; j++, x += del) { sum += BMfunction.calculate(x); } result = 0.5 * (result + (b - a) * sum / tnm); } return result; } template integrationType iterativeIntegration(const FunctorType& BMfunction, integrationType a, integrationType b) { integrationType ost = 0.0; integrationType os = 0.0; int maxSteps = 50; double eps = 1e-6; int i = 1; for (; i <= maxSteps; ++i) { integrationType st = trapzd(BMfunction, a, b, i); integrationType s = (4.0 * st - ost) / 3.0; if (i > 5) { if (fabs(s - os) < eps * fabs(os) || (s == 0.0 && os == 0.0)) { return s; } } os = s; ost = st; } //too many iterations, this should never be reachable! throw rttb::core::InvalidParameterException("Integral calculation failed: too many iterations! "); } } -} \ No newline at end of file +} diff --git a/testing/core/CreateTestStructure.cpp b/testing/core/CreateTestStructure.cpp index 60ee7c1..79df5dc 100644 --- a/testing/core/CreateTestStructure.cpp +++ b/testing/core/CreateTestStructure.cpp @@ -1,670 +1,670 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include "CreateTestStructure.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace testing { CreateTestStructure::~CreateTestStructure() {} CreateTestStructure::CreateTestStructure(const core::GeometricInfo& aGeoInfo) { _geoInfo = aGeoInfo; } PolygonType CreateTestStructure::createPolygonLeftUpper(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; - for (int i = 0; i < aVoxelVector.size(); i++) + for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); polygon.push_back(p1); std::cout << "(" << p1.x() << "," << p1.y() << "," << p1.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonCenter(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; - for (int i = 0; i < aVoxelVector.size(); i++) + for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() / 2; p(1) = p1.y() + delta.y() / 2; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonBetweenUpperLeftAndCenter( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; - for (int i = 0; i < aVoxelVector.size(); i++) + for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() / 4; p(1) = p1.y() + delta.y() / 4; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonBetweenLowerRightAndCenter( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; - for (int i = 0; i < aVoxelVector.size(); i++) + for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() * 0.75; p(1) = p1.y() + delta.y() * 0.75; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonLeftEdgeMiddle(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; - for (int i = 0; i < aVoxelVector.size(); i++) + for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x(); p(1) = p1.y() + delta.y() * 0.5; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonUpperCenter(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; - for (int i = 0; i < aVoxelVector.size(); i++) + for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() * 0.5; p(1) = p1.y(); p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonIntermediatePoints(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; - for (int i = 0; i < aVoxelVector.size(); i++) + for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; VoxelGridIndex3D voxelIndex2; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; if (i < (aVoxelVector.size() - 1)) { voxelIndex2(0) = (aVoxelVector.at(i + 1)).x(); voxelIndex2(1) = (aVoxelVector.at(i + 1)).y(); voxelIndex2(2) = sliceNumber; } else { voxelIndex2(0) = (aVoxelVector.at(0)).x(); voxelIndex2(1) = (aVoxelVector.at(0)).y(); voxelIndex2(2) = sliceNumber; } WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p2; _geoInfo.indexToWorldCoordinate(voxelIndex2, p2); SpacingVectorType3D delta2 = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x() + delta.x() * 0.75; wp1(1) = p1.y() + delta.y() * 0.75; wp1(2) = p1.z(); WorldCoordinate3D wp2; wp2(0) = p2.x() + delta.x() * 0.75; wp2(1) = p2.y() + delta.y() * 0.75; wp2(2) = p2.z(); polygon.push_back(wp1); double diffX = (wp2.x() - wp1.x()) / 1000.0; double diffY = (wp2.y() - wp1.y()) / 1000.0; WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; for (int i = 0 ; i < 1000 ; i++) { wp_intermediate(0) = wp1.x() + i * diffX; wp_intermediate(1) = wp1.y() + i * diffY; polygon.push_back(wp_intermediate); } } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonCircle(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); double radius = 2 * delta.x(); double frac_radius = (radius * 0.001); double correct_y = (delta.x() / delta.y()); for (int i = 0 ; i <= 1000 ; i++) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() + frac_radius * i; wp_intermediate(1) = wp1.y() - (y_offset * correct_y) ; polygon.push_back(wp_intermediate); } for (int i = 1000 ; i >= 0 ; i--) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))) ; wp_intermediate(0) = wp1.x() + frac_radius * i; wp_intermediate(1) = wp1.y() + (y_offset * correct_y); polygon.push_back(wp_intermediate); } for (int i = 0 ; i <= 1000 ; i++) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() - frac_radius * i; wp_intermediate(1) = wp1.y() + y_offset * correct_y ; polygon.push_back(wp_intermediate); } for (int i = 1000 ; i >= 0 ; i--) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() - frac_radius * i ; wp_intermediate(1) = wp1.y() - (y_offset * correct_y); polygon.push_back(wp_intermediate); } } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSeveralSectionsInsideOneVoxelA( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x() + (delta.x() * 0.25); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.25); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.75); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.75); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.25); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.25); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.5); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.75); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.75); wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSelfTouchingA(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y(); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSelfTouchingB(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y(); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } }//testing }//rttb diff --git a/testing/core/DummyStructure.cpp b/testing/core/DummyStructure.cpp index 4fef460..f8dd8dc 100644 --- a/testing/core/DummyStructure.cpp +++ b/testing/core/DummyStructure.cpp @@ -1,656 +1,655 @@ // ----------------------------------------------------------------------- // 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 "DummyStructure.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace testing { DummyStructure::~DummyStructure() {} DummyStructure::DummyStructure(const core::GeometricInfo& aGeoInfo) { _geoInfo = aGeoInfo; } core::Structure DummyStructure::CreateRectangularStructureCentered(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 1); VoxelGridIndex2D another_i2(5, 1); VoxelGridIndex2D another_i3(5, 5); VoxelGridIndex2D another_i4(2, 5); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure_rectangular_centered = core::Structure(another_polySeq); return test_structure_rectangular_centered; } core::Structure DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacement( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(5, 1); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(5, 7); VoxelGridIndex2D another_i4(2, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonBetweenUpperLeftAndCenter( another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure - DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRight( - GridIndexType zPlane) + DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRight() { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(5, 1); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(5, 7); VoxelGridIndex2D another_i4(2, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonSequenceType another_polySeq; core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRightCounterClock( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(5, 7); VoxelGridIndex2D another_i3(8, 4); VoxelGridIndex2D another_i4(5, 1); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonBetweenLowerRightAndCenter( another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRightCounterClockIntermediatePoints( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(5, 7); VoxelGridIndex2D another_i3(8, 4); VoxelGridIndex2D another_i4(5, 1); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonIntermediatePoints(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureSeveralSeperateSectionsInsideOneVoxel( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 2); another_voxelVector.push_back(another_i1); PolygonType another_polygon1 = another_cts.createStructureSeveralSectionsInsideOneVoxelA( another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureSelfTouchingA(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 2); another_voxelVector.push_back(another_i1); PolygonType another_polygon1 = another_cts.createStructureSelfTouchingA(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureIntersectingTwoPolygonsInDifferentSlices( GridIndexType zPlane1, GridIndexType zPlane2) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); CreateTestStructure one_more_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(8, 6); VoxelGridIndex2D another_i4(2, 6); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector , zPlane1); std::vector one_more_voxelVector; VoxelGridIndex2D one_more_i1(3, 5); VoxelGridIndex2D one_more_i2(9, 5); VoxelGridIndex2D one_more_i3(9, 7); VoxelGridIndex2D one_more_i4(3, 7); one_more_voxelVector.push_back(one_more_i1); one_more_voxelVector.push_back(one_more_i2); one_more_voxelVector.push_back(one_more_i3); one_more_voxelVector.push_back(one_more_i4); PolygonType another_polygon2 = one_more_cts.createPolygonCenter(one_more_voxelVector , zPlane2); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); another_polySeq.push_back(another_polygon2); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureIntersectingTwoPolygons(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); CreateTestStructure one_more_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(8, 6); VoxelGridIndex2D another_i4(2, 6); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector , zPlane); std::vector one_more_voxelVector; VoxelGridIndex2D one_more_i1(3, 5); VoxelGridIndex2D one_more_i2(9, 5); VoxelGridIndex2D one_more_i3(9, 7); VoxelGridIndex2D one_more_i4(3, 7); one_more_voxelVector.push_back(one_more_i1); one_more_voxelVector.push_back(one_more_i2); one_more_voxelVector.push_back(one_more_i3); one_more_voxelVector.push_back(one_more_i4); PolygonType another_polygon2 = one_more_cts.createPolygonCenter(one_more_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); another_polySeq.push_back(another_polygon2); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureIntersecting(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(2, 6); VoxelGridIndex2D another_i4(8, 6); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouches(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 4); VoxelGridIndex2D another_i2(2, 8); VoxelGridIndex2D another_i3(4, 8); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); PolygonType another_polygon1 = another_cts.createPolygonUpperCenter(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesRotatedPointDoubeled( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i4(2, 4); VoxelGridIndex2D another_i2(4, 4); VoxelGridIndex2D another_i3(3, 8); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); PolygonType another_polygon1 = another_cts.createPolygonUpperCenter(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesCounterClockRotatedOnePointFivePi( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 5); VoxelGridIndex2D another_i2(3, 5); VoxelGridIndex2D another_i3(7, 6); VoxelGridIndex2D another_i4(7, 6); VoxelGridIndex2D another_i5(7, 4); VoxelGridIndex2D another_i6(7, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftEdgeMiddle(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesCounterClockRotatedQuaterPi( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(7, 5); VoxelGridIndex2D another_i2(7, 5); VoxelGridIndex2D another_i3(3, 6); VoxelGridIndex2D another_i4(3, 6); VoxelGridIndex2D another_i5(3, 4); VoxelGridIndex2D another_i6(3, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftEdgeMiddle(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesRotatedQuaterPi( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 4); VoxelGridIndex2D another_i2(3, 4); VoxelGridIndex2D another_i3(3, 6); VoxelGridIndex2D another_i4(3, 6); VoxelGridIndex2D another_i5(7, 5); VoxelGridIndex2D another_i6(7, 5); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftEdgeMiddle(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureCircle(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(4, 4); another_voxelVector.push_back(another_i1); PolygonType another_polygon1 = another_cts.createPolygonCircle(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesUpperRight( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 4); VoxelGridIndex2D another_i3(4, 5); VoxelGridIndex2D another_i5(5, 3); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i5); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesLowerRight( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 2); VoxelGridIndex2D another_i2(2, 2); VoxelGridIndex2D another_i3(3, 1); VoxelGridIndex2D another_i4(3, 1); VoxelGridIndex2D another_i5(4, 3); VoxelGridIndex2D another_i6(4, 3); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesLowerLeft( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 3); VoxelGridIndex2D another_i2(3, 3); VoxelGridIndex2D another_i3(4, 4); VoxelGridIndex2D another_i4(4, 4); VoxelGridIndex2D another_i5(2, 5); VoxelGridIndex2D another_i6(2, 5); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesUpperLeft( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(5, 5); VoxelGridIndex2D another_i2(5, 5); VoxelGridIndex2D another_i3(4, 6); VoxelGridIndex2D another_i4(4, 6); VoxelGridIndex2D another_i5(3, 4); VoxelGridIndex2D another_i6(3, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } void DummyStructure::ShowTestStructure(core::Structure aStructure) { WorldCoordinate3D aPoint(0); PolygonSequenceType strVector = aStructure.getStructureVector(); - for (unsigned int struct_index = 0 ; struct_index < strVector.size() ; struct_index++) + for (size_t struct_index = 0 ; struct_index < strVector.size() ; struct_index++) { - for (int point_index = 0 ; point_index < strVector.at(struct_index).size() ; point_index++) + for (size_t point_index = 0 ; point_index < strVector.at(struct_index).size() ; point_index++) { aPoint = strVector.at(struct_index).at(point_index); std::cout << " aPoint.x " << aPoint.x() << std::endl; std::cout << " aPoint.y " << aPoint.y() << std::endl; std::cout << " aPoint.z " << aPoint.z() << std::endl; } } } core::Structure DummyStructure::CreateRectangularStructureUpperLeftRotated(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(5, 1); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(5, 7); VoxelGridIndex2D another_i4(2, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } }//testing }//rttb diff --git a/testing/core/DummyStructure.h b/testing/core/DummyStructure.h index bdc8894..c59ea11 100644 --- a/testing/core/DummyStructure.h +++ b/testing/core/DummyStructure.h @@ -1,101 +1,100 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include "rttbStructure.h" #include "CreateTestStructure.h" #include "rttbBaseType.h" #include "rttbGeometricInfo.h" namespace rttb { namespace testing { /*! @class DummyStructure @brief generate simple geometric testing structures. The maximal x coordinate used is 9 and the maximal y coordinate is 8. Make sure the geometricInfo corresponds to a sufficiently large data grid. @see CreateTestStructures */ class DummyStructure { private: core::GeometricInfo _geoInfo; public: ~DummyStructure(); DummyStructure(const core::GeometricInfo& aGeoInfo); const core::GeometricInfo& getGeometricInfo() { return _geoInfo; }; core::Structure CreateRectangularStructureCentered(GridIndexType zPlane); core::Structure CreateTestStructureCircle(GridIndexType zPlane); core::Structure CreateRectangularStructureUpperLeftRotated(GridIndexType zPlane); core::Structure CreateTestStructureSeveralSeperateSectionsInsideOneVoxel(GridIndexType zPlane); core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacement( GridIndexType zPlane); - core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRight( - GridIndexType zPlane); + core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRight(); core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRightCounterClock( GridIndexType zPlane); core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRightCounterClockIntermediatePoints( GridIndexType zPlane); core::Structure CreateTestStructureSelfTouchingA(GridIndexType zPlane); core::Structure CreateTestStructureIntersecting(GridIndexType zPlane); core::Structure CreateTestStructureIntersectingTwoPolygons(GridIndexType zPlane); core::Structure CreateTestStructureIntersectingTwoPolygonsInDifferentSlices(GridIndexType zPlane1, GridIndexType zPlane2); core::Structure CreateTestStructureInsideInsideTouches(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesRotatedQuaterPi(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesCounterClockRotatedQuaterPi( GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesCounterClockRotatedOnePointFivePi( GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesRotatedPointDoubeled(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesUpperLeft(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesLowerLeft(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesLowerRight(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesUpperRight(GridIndexType zPlane); void ShowTestStructure(core::Structure aStructure); }; }//testing }//rttb diff --git a/testing/core/StructureTest.cpp b/testing/core/StructureTest.cpp index 42b550a..493da11 100644 --- a/testing/core/StructureTest.cpp +++ b/testing/core/StructureTest.cpp @@ -1,138 +1,138 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbStructure.h" #include "rttbInvalidParameterException.h" #include "DummyStructure.h" #include "DummyDoseAccessor.h" namespace rttb { namespace testing { /*! @brief StructureTest - tests the API for Structure 1) constructors 2) get/setXX */ int StructureTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; boost::shared_ptr spTestDoseAccessor = boost::make_shared(); DummyStructure myStructGenerator(spTestDoseAccessor->getGeometricInfo()); //1) constructors CHECK_NO_THROW(core::Structure()); core::Structure emptyTestStruct; CHECK_EQUAL("None", emptyTestStruct.getLabel()); CHECK_NO_THROW(emptyTestStruct.getUID()); GridIndexType zPlane = 4; core::Structure rect = myStructGenerator.CreateRectangularStructureCentered(zPlane); CHECK_NO_THROW(core::Structure(rect.getStructureVector())); core::Structure rect2 = core::Structure(rect.getStructureVector()); CHECK_EQUAL(rect.getLabel(), rect2.getLabel()); CHECK(rect.getUID() != rect2.getUID()); PolygonSequenceType rectVec = rect.getStructureVector(); PolygonSequenceType rect2Vec = rect2.getStructureVector(); CHECK_EQUAL(rectVec.size(), rect2Vec.size()); PolygonSequenceType::iterator it = rectVec.begin(); PolygonSequenceType::iterator it2 = rect2Vec.begin(); for (; it != rectVec.end(); it++) { CHECK_EQUAL(it->size(), it2->size()); PolygonType::iterator pit = it->begin(); PolygonType::iterator pit2 = it2->begin(); for (; pit != it->end(); pit++) { CHECK_EQUAL(*(pit), *(pit2)); pit2++; } it2++; } CHECK_NO_THROW(core::Structure rect3 = rect); core::Structure rect3 = rect; CHECK_EQUAL(rect.getLabel(), rect3.getLabel()); CHECK_EQUAL(rect.getUID(), rect3.getUID()); PolygonSequenceType rect3Vec = rect3.getStructureVector(); CHECK_EQUAL(rectVec.size(), rect3Vec.size()); it = rectVec.begin(); PolygonSequenceType::iterator it3 = rect3Vec.begin(); for (; it != rectVec.end(); it++) { CHECK_EQUAL(it->size(), it3->size()); PolygonType::iterator pit = it->begin(); PolygonType::iterator pit3 = it3->begin(); for (; pit != it->end(); pit++) { CHECK_EQUAL(*(pit), *(pit3)); pit3++; } it3++; } //2) get/setXX CHECK_EQUAL("None", emptyTestStruct.getLabel()); CHECK_NO_THROW(emptyTestStruct.setLabel("NEW Label")); CHECK_EQUAL("NEW Label", emptyTestStruct.getLabel()); CHECK_NO_THROW(emptyTestStruct.getUID()); CHECK_NO_THROW(emptyTestStruct.setUID("1.2.345.67.8.9")); CHECK_EQUAL("1.2.345.67.8.9", emptyTestStruct.getUID()); CHECK((emptyTestStruct.getStructureVector()).empty()); CHECK_EQUAL(0, emptyTestStruct.getNumberOfEndpoints()); CHECK_EQUAL(4, rect.getNumberOfEndpoints()); CHECK_EQUAL(rect.getNumberOfEndpoints(), rect2.getNumberOfEndpoints()); core::Structure circ = myStructGenerator.CreateTestStructureCircle(zPlane); CHECK_EQUAL(4004, circ.getNumberOfEndpoints()); core::Structure multiPoly = myStructGenerator.CreateTestStructureIntersectingTwoPolygonsInDifferentSlices(zPlane, zPlane + 1); CHECK_EQUAL(8, multiPoly.getNumberOfEndpoints()); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing -}//rttb \ No newline at end of file +}//rttb diff --git a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolDifferentCommandsTest.cpp b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolDifferentCommandsTest.cpp index e5ec12b..09cb68d 100644 --- a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolDifferentCommandsTest.cpp +++ b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolDifferentCommandsTest.cpp @@ -1,110 +1,110 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include #include #include /*! @brief VoxelizerToolTest3. Test the output, multipleStructs and the booleanVoxelization parameter. */ namespace rttb { namespace testing { //path to the current running directory. VoxelizerTool is in the same directory (Debug/Release) extern const char* _callingAppPath; int VoxelizerToolDifferentCommandsTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string voxelizerToolExe; std::string tempDirectory; std::string structFile; std::string referenceFile; if (argc > 4) { voxelizerToolExe = argv[1]; tempDirectory = argv[2]; structFile = argv[3]; referenceFile = argv[4]; } std::vector commands; commands.push_back("\"Niere.*\" -m -o Test.hdr"); commands.push_back("\"Leber\" -o Boolean.hdr -v"); std::vector filenames; filenames.push_back("Test_Niere li"); filenames.push_back("Test_Niere re"); filenames.push_back("Boolean_Leber"); boost::filesystem::path callingPath(_callingAppPath); std::string voxelizerToolExeWithPath = callingPath.parent_path().string() + "/" + voxelizerToolExe; std::string baseCommand = voxelizerToolExeWithPath; baseCommand += " -s " + structFile; baseCommand += " -r " + referenceFile; baseCommand += " -e "; - for (int i = 0; i < commands.size(); i++) + for (size_t i = 0; i < commands.size(); i++) { std::string command = baseCommand + commands.at(i); int returnValue = system(command.c_str()); std::cout << "Command line call: " + command << std::endl; CHECK_EQUAL(returnValue, 0); } std::string helpCommand = voxelizerToolExeWithPath + " -h"; int returnValue = system(helpCommand.c_str()); std::cout << "Command line call: " + helpCommand << std::endl; CHECK_EQUAL(returnValue, 0); - for (int i = 0; i < filenames.size(); i++) + for (size_t i = 0; i < filenames.size(); i++) { const std::string HDRfileName = tempDirectory + "/" + filenames.at(i) + ".hdr"; boost::filesystem::path HDRFile(HDRfileName); const std::string IMGfileName = tempDirectory + "/" + filenames.at(i) + ".img"; boost::filesystem::path IMGFile(IMGfileName); CHECK_EQUAL(boost::filesystem::exists(HDRFile), true); CHECK_EQUAL(boost::filesystem::exists(IMGFile), true); if (boost::filesystem::exists(IMGFile)) { boost::filesystem::remove(IMGFile); } if (boost::filesystem::exists(HDRFile)) { boost::filesystem::remove(HDRFile); } } RETURN_AND_REPORT_TEST_SUCCESS; } } } diff --git a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelValueTest.cpp b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelValueTest.cpp index c03d58a..2db1965 100644 --- a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelValueTest.cpp +++ b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelValueTest.cpp @@ -1,134 +1,134 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include #include "itkImage.h" #include "itkImageFileReader.h" #include #include /*! @brief VoxelizerToolTest5. Search the coordinate at the Image and return the Voxel(Pixel) value. */ namespace rttb { namespace testing { //path to the current running directory. VoxelizerTool is in the same directory (Debug/Release) extern const char* _callingAppPath; int VoxelizerToolVoxelValue(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; typedef itk::Image< double, 3 > ImageType; typedef itk::ImageFileReader ReaderType; std::string voxelizerToolExe; std::string tempDirectory; std::string structFile; std::string referenceFile; std::string structName; if (argc > 5) { voxelizerToolExe = argv[1]; tempDirectory = argv[2]; structFile = argv[3]; referenceFile = argv[4]; structName = argv[5]; } boost::filesystem::path callingPath(_callingAppPath); std::string voxelizerToolExeWithPath = callingPath.parent_path().string() + "/" + voxelizerToolExe; std::string command = voxelizerToolExeWithPath; command += " -s " + structFile; command += " -r " + referenceFile; command += " -e " + structName; int returnValue = system(command.c_str()); CHECK_EQUAL(returnValue, 0); //image values taken in Mevislab //Index inside ImageType::IndexType voxelInside1 = {{20, 30, 30}}; ImageType::IndexType voxelInside2 = {{30, 10, 40}}; //Outside index ImageType::IndexType voxelOutside1 = {{40, 30, 30}}; ImageType::IndexType voxelOutside2 = {{10, 40, 30}}; //Border index ImageType::IndexType voxelBorder1 = {{12, 23, 27}}; ImageType::IndexType voxelBorder2 = {{34, 21, 31}}; std::vector voxelIndices = boost::assign::list_of(voxelInside1)(voxelInside2)( voxelOutside1)( voxelOutside2)(voxelBorder1)(voxelBorder2); std::vector expectedVoxelValues = boost::assign::list_of(1.0)(1.0)(0.0)(0.0)( 0.265865)(0.819613); std::string filenameHDRWithVoxelization = tempDirectory + "/out_" + structName + ".hdr"; std::string filenameIMGWithVoxelization = tempDirectory + "/out_" + structName + ".img"; CHECK(boost::filesystem::exists(filenameHDRWithVoxelization)); CHECK(boost::filesystem::exists(filenameIMGWithVoxelization)); if (boost::filesystem::exists(filenameHDRWithVoxelization)) { ReaderType::Pointer reader = ReaderType::New(); reader->SetFileName(filenameHDRWithVoxelization); reader->Update(); ReaderType::OutputImageType::ConstPointer image = reader->GetOutput(); - for (int i = 0; i < voxelIndices.size(); i++) + for (size_t i = 0; i < voxelIndices.size(); i++) { ImageType::PixelType voxelValue = image->GetPixel(voxelIndices.at(i)); ImageType::PixelType expectedVoxelValue = expectedVoxelValues.at(i); if (expectedVoxelValue == 0.0 || expectedVoxelValue == 1.0) { CHECK_EQUAL(voxelValue, expectedVoxelValue); } else { CHECK_CLOSE(voxelValue, expectedVoxelValue, 1e-4); } } if (boost::filesystem::exists(filenameHDRWithVoxelization)) { boost::filesystem::remove(filenameHDRWithVoxelization); } if (boost::filesystem::exists(filenameIMGWithVoxelization)) { boost::filesystem::remove(filenameIMGWithVoxelization); } } RETURN_AND_REPORT_TEST_SUCCESS; } } } diff --git a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerAllStructsTest.cpp b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerAllStructsTest.cpp index 6ce2031..96880da 100644 --- a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerAllStructsTest.cpp +++ b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerAllStructsTest.cpp @@ -1,111 +1,111 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include #include /*! @brief VoxelizerToolTest. Tests a selection of structs. */ namespace rttb { namespace testing { //path to the current running directory. VoxelizerTool is in the same directory (Debug/Release) extern const char* _callingAppPath; int VoxelizerToolVoxelizerAllStructsTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string voxelizerToolExe; std::string tempDirectory; std::string structFile; std::string referenceFile; if (argc > 4) { voxelizerToolExe = argv[1]; tempDirectory = argv[2]; structFile = argv[3]; referenceFile = argv[4]; } boost::filesystem::path callingPath(_callingAppPath); std::string voxelizerToolExeWithPath = callingPath.parent_path().string() + "/" + voxelizerToolExe; std::vector structNames; structNames.push_back("Niere.*"); structNames.push_back("Magen/DD"); structNames.push_back("PTV"); //structure names will be used for file naming, BUT '.' in the end will be cropped and '/' will be replaced by '_'. Thus, the different filenames. std::vector filenames; filenames.push_back("Niere re"); filenames.push_back("Magen_DD"); filenames.push_back("PTV"); std::string baseCommand = voxelizerToolExeWithPath; baseCommand += " -s " + structFile; baseCommand += " -r " + referenceFile; baseCommand += " -e \""; - for (int i = 0; i < structNames.size(); i++) + for (size_t i = 0; i < structNames.size(); i++) { std::string command = baseCommand + structNames.at(i) + "\""; std::cout << "Command line call: " + command << std::endl; int returnValue = system(command.c_str()); CHECK_EQUAL(returnValue, 0); const std::string HDRfileName = tempDirectory + "/out_" + filenames.at(i) + ".hdr"; boost::filesystem::path HDRFile(tempDirectory); HDRFile /= "out_" + filenames.at(i) + ".hdr"; const std::string IMGfileName = tempDirectory + "/out_" + filenames.at(i) + ".img"; boost::filesystem::path IMGFile(tempDirectory); IMGFile /= "out_" + filenames.at(i) + ".img"; CHECK_EQUAL( boost::filesystem::exists(HDRFile), true); CHECK_EQUAL( boost::filesystem::exists(IMGFile), true); if (boost::filesystem::exists(IMGFile)) { boost::filesystem::remove(IMGFile); } if (boost::filesystem::exists(HDRFile)) { boost::filesystem::remove(HDRFile); } } RETURN_AND_REPORT_TEST_SUCCESS; } } } diff --git a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerBoostLegacyTest.cpp b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerBoostLegacyTest.cpp index 51ff1ad..aad393b 100644 --- a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerBoostLegacyTest.cpp +++ b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerBoostLegacyTest.cpp @@ -1,104 +1,104 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include #include #include /*! @brief VoxelizerToolTest5. Test the paramter boost and legacy Voxelization. */ namespace rttb { namespace testing { //path to the current running directory. VoxelizerTool is in the same directory (Debug/Release) extern const char* _callingAppPath; int VoxelizerToolVoxelizerBoostLegacy(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string voxelizerToolExe; std::string tempDirectory; std::string structFile; std::string referenceFile; if (argc > 4) { voxelizerToolExe = argv[1]; tempDirectory = argv[2]; structFile = argv[3]; referenceFile = argv[4]; } std::vector commands; commands.push_back("PTV -o Legacy.hdr -l"); commands.push_back("PTV -o Boost.hdr -b"); std::vector filenames; filenames.push_back("Boost_PTV"); filenames.push_back("Legacy_PTV"); boost::filesystem::path callingPath(_callingAppPath); std::string voxelizerToolExeWithPath = callingPath.parent_path().string() + "/" + voxelizerToolExe; std::string baseCommand = voxelizerToolExeWithPath; baseCommand += " -s " + structFile; baseCommand += " -r " + referenceFile; baseCommand += " -e "; - for (int i = 0; i < commands.size(); i++) + for (size_t i = 0; i < commands.size(); i++) { std::string command = baseCommand + commands.at(i); int returnValue = system(command.c_str()); std::cout << "Command line call: " + command << std::endl; CHECK_EQUAL(returnValue, 0); } - for (int i = 0; i < filenames.size(); i++) + for (size_t i = 0; i < filenames.size(); i++) { const std::string HDRfileName = tempDirectory + "/" + filenames.at(i) + ".hdr"; boost::filesystem::path HDRFile(HDRfileName); const std::string IMGfileName = tempDirectory + "/" + filenames.at(i) + ".img"; boost::filesystem::path IMGFile(IMGfileName); CHECK_EQUAL(boost::filesystem::exists(HDRFile), true); CHECK_EQUAL(boost::filesystem::exists(IMGFile), true); if (boost::filesystem::exists(IMGFile)) { boost::filesystem::remove(IMGFile); } if (boost::filesystem::exists(HDRFile)) { boost::filesystem::remove(HDRFile); } } RETURN_AND_REPORT_TEST_SUCCESS; } } -} \ No newline at end of file +} diff --git a/testing/examples/DVHCalculatorComparisonTest.cpp b/testing/examples/DVHCalculatorComparisonTest.cpp index 8fa933b..3c7c539 100644 --- a/testing/examples/DVHCalculatorComparisonTest.cpp +++ b/testing/examples/DVHCalculatorComparisonTest.cpp @@ -1,365 +1,365 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVHCalculator.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbOTBMaskAccessor.h" #include "rttbDVHTxtFileReader.h" namespace rttb { namespace testing { /*! @brief DVHCalculatorTest. Test if calculation in new architecture returns similar results to the original implementation. WARNING: The values for comparison need to be adjusted if the input files are changed! This test can be used to get more detailed information, but it will always fail, because differences in voxelization accuracy, especially the ones caused by the change from double to float precission will not cause considerable deviations in the structure sizes, which correspond to considerable differences in the calculated DVHs. Even in double precission differences up to 0.005 between values from old and new implementation can occure. */ int DVHCalculatorComparisonTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: structure file name // 2: dose1 file name // 3: dose2 file name // 4: dose3 file name std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string RTDOSE2_FILENAME; std::string RTDOSE3_FILENAME; std::string COMPARISON_DVH_FOLDER; if (argc > 1) { RTSTRUCT_FILENAME = argv[1]; } if (argc > 2) { RTDOSE_FILENAME = argv[2]; } if (argc > 3) { RTDOSE2_FILENAME = argv[3]; } if (argc > 4) { RTDOSE3_FILENAME = argv[4]; } if (argc > 5) { COMPARISON_DVH_FOLDER = argv[5]; } OFCondition status; DcmFileFormat fileformat; /* read dicom-rt dose */ io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create a vector of MaskAccessors (one for each structure) StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); std::vector rtStructSetMaskAccessorVec; ::DRTDoseIOD rtdose2; io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator2(RTDOSE2_FILENAME.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); ::DRTDoseIOD rtdose3; io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator3(RTDOSE3_FILENAME.c_str()); DoseAccessorPointer doseAccessor3(doseAccessorGenerator3.generateDoseAccessor()); double maxDifference = 0; double difference = 0; double minDifference = 1000; clock_t start(clock()); if (rtStructureSet->getNumberOfStructures() > 0) { - for (int j = 0; j < rtStructureSet->getNumberOfStructures(); j++) + for (int j = 0; j < static_cast(rtStructureSet->getNumberOfStructures()); j++) { //create MaskAccessor ::boost::shared_ptr spOTBMaskAccessor = ::boost::make_shared(rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); spOTBMaskAccessor->updateMask(); MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor); //create corresponding MaskedDoseIterator ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(spMaskAccessor, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); //store MaskAccessor rtStructSetMaskAccessorVec.push_back(spMaskAccessor); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); std::string dvhFileName = "dvh1"; std::string label = (rtStructureSet->getStructure(j))->getLabel(); dvhFileName.append(label); if (dvhFileName.find("/") != std::string::npos) { dvhFileName.replace(dvhFileName.find("/"), 1, ""); } std::cout << "=== Dose 1: " << label << "===" << std::endl; rttb::io::other::DVHTxtFileReader dvhOriginalReader = rttb::io::other::DVHTxtFileReader( COMPARISON_DVH_FOLDER + dvhFileName + ".txt"); rttb::core::DVH dvhOrig = *(dvhOriginalReader.generateDVH()); rttb::core::DVH::DataDifferentialType dvhOrigData = dvhOrig.getDataDifferential(); rttb::core::DVH dvh = *(calc.generateDVH()); rttb::core::DVH::DataDifferentialType dvhData = dvh.getDataDifferential(); CHECK_EQUAL(dvhOrig, dvh); CHECK_CLOSE(dvhOrig.getMaximum(), dvh.getMaximum(), errorConstant); CHECK_CLOSE(dvhOrig.getMinimum(), dvh.getMinimum(), errorConstant); CHECK_CLOSE(dvhOrig.getMean(), dvh.getMean(), errorConstant); rttb::core::DVH::DataDifferentialType::iterator it; rttb::core::DVH::DataDifferentialType::iterator itOrig; itOrig = dvhOrigData.begin(); for (it = dvhData.begin(); it != dvhData.end() || itOrig != dvhOrigData.end(); ++it, ++itOrig) { CHECK_CLOSE(*(itOrig), *(it), errorConstant); std::cout << std::setprecision(20) << "difference: " << abs(*(itOrig) - * (it)) << std::endl; difference = abs(*(itOrig) - * (it)); if (difference < minDifference) { minDifference = difference; } if (difference > maxDifference) { maxDifference = difference; } } } } clock_t finish(clock()); std::cout << std::setprecision(20) << "max(difference): " << maxDifference << std::endl; std::cout << std::setprecision(20) << "min(difference): " << minDifference << std::endl; std::cout << "DVH Calculation time: " << finish - start << " ms" << std::endl; maxDifference = 0; minDifference = 1000; clock_t start2(clock()); - for (int j = 0; j < rtStructSetMaskAccessorVec.size(); j++) + for (size_t j = 0; j < rtStructSetMaskAccessorVec.size(); j++) { //create corresponding MaskedDoseIterator ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(rtStructSetMaskAccessorVec.at(j), doseAccessor2); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor2->getUID()); std::string dvhFileName = "dvh2"; std::string label = (rtStructureSet->getStructure(j))->getLabel(); dvhFileName.append(label); if (dvhFileName.find("/") != std::string::npos) { dvhFileName.replace(dvhFileName.find("/"), 1, ""); } std::cout << "=== Dose 2: " << label << "===" << std::endl; rttb::io::other::DVHTxtFileReader dvhOriginalReader = rttb::io::other::DVHTxtFileReader( COMPARISON_DVH_FOLDER + dvhFileName + ".txt"); rttb::core::DVH dvhOrig = *(dvhOriginalReader.generateDVH()); rttb::core::DVH::DataDifferentialType dvhOrigData = dvhOrig.getDataDifferential(); rttb::core::DVH dvh = *(calc.generateDVH()); rttb::core::DVH::DataDifferentialType dvhData = dvh.getDataDifferential(); CHECK_EQUAL(dvhOrig, dvh); CHECK_CLOSE(dvhOrig.getMaximum(), dvh.getMaximum(), errorConstant); CHECK_CLOSE(dvhOrig.getMinimum(), dvh.getMinimum(), errorConstant); CHECK_CLOSE(dvhOrig.getMean(), dvh.getMean(), errorConstant); rttb::core::DVH::DataDifferentialType::iterator it; rttb::core::DVH::DataDifferentialType::iterator itOrig; itOrig = dvhOrigData.begin(); for (it = dvhData.begin(); it != dvhData.end() || itOrig != dvhOrigData.end(); ++it, ++itOrig) { CHECK_CLOSE(*(itOrig), *(it), errorConstant); std::cout << std::setprecision(20) << "difference: " << abs(*(itOrig) - * (it)) << std::endl; difference = abs(*(itOrig) - * (it)); if (difference > 10) { std::cout << "large difference: " << abs(*(itOrig) - * (it)) << " = " << *(itOrig) << " - " << * (it) << std::endl; } if (difference < minDifference) { minDifference = difference; } if (difference > maxDifference) { maxDifference = difference; } } } clock_t finish2(clock()); std::cout << std::setprecision(20) << "max(difference): " << maxDifference << std::endl; std::cout << std::setprecision(20) << "min(difference): " << minDifference << std::endl; std::cout << "Reset dose 2, DVH Calculation time: " << finish2 - start2 << " ms" << std::endl; maxDifference = 0; minDifference = 1000; clock_t start3(clock()); - for (int j = 0; j < rtStructSetMaskAccessorVec.size(); j++) + for (size_t j = 0; j < rtStructSetMaskAccessorVec.size(); j++) { //create corresponding MaskedDoseIterator ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(rtStructSetMaskAccessorVec.at(j), doseAccessor3); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor3->getUID()); std::string dvhFileName = "dvh3"; std::string label = (rtStructureSet->getStructure(j))->getLabel(); dvhFileName.append(label); if (dvhFileName.find("/") != std::string::npos) { dvhFileName.replace(dvhFileName.find("/"), 1, ""); } std::cout << "=== Dose 3: " << label << "===" << std::endl; rttb::io::other::DVHTxtFileReader dvhOriginalReader = rttb::io::other::DVHTxtFileReader( COMPARISON_DVH_FOLDER + dvhFileName + ".txt"); rttb::core::DVH dvhOrig = *(dvhOriginalReader.generateDVH()); rttb::core::DVH::DataDifferentialType dvhOrigData = dvhOrig.getDataDifferential(); rttb::core::DVH dvh = *(calc.generateDVH()); rttb::core::DVH::DataDifferentialType dvhData = dvh.getDataDifferential(); CHECK_EQUAL(dvhOrig, dvh); CHECK_CLOSE(dvhOrig.getMaximum(), dvh.getMaximum(), errorConstant); CHECK_CLOSE(dvhOrig.getMinimum(), dvh.getMinimum(), errorConstant); CHECK_CLOSE(dvhOrig.getMean(), dvh.getMean(), errorConstant); rttb::core::DVH::DataDifferentialType::iterator it; rttb::core::DVH::DataDifferentialType::iterator itOrig; itOrig = dvhOrigData.begin(); for (it = dvhData.begin(); it != dvhData.end() || itOrig != dvhOrigData.end(); ++it, ++itOrig) { CHECK_CLOSE(*(itOrig), *(it), errorConstant); std::cout << std::setprecision(20) << "difference: " << abs(*(itOrig) - * (it)) << std::endl; if (difference > 10) { std::cout << "large difference: " << abs(*(itOrig) - * (it)) << " = " << *(itOrig) << " - " << * (it) << std::endl; } difference = abs(*(itOrig) - * (it)); if (difference < minDifference) { minDifference = difference; } if (difference > maxDifference) { maxDifference = difference; } } } clock_t finish3(clock()); std::cout << std::setprecision(20) << "max(difference): " << maxDifference << std::endl; std::cout << std::setprecision(20) << "min(difference): " << minDifference << std::endl; std::cout << "Reset dose 3, DVH Calculation time: " << finish3 - start3 << " ms" << std::endl; RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/examples/DVHCalculatorExampleTest.cpp b/testing/examples/DVHCalculatorExampleTest.cpp index 6685f09..89ce459 100644 --- a/testing/examples/DVHCalculatorExampleTest.cpp +++ b/testing/examples/DVHCalculatorExampleTest.cpp @@ -1,305 +1,305 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVHCalculator.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbOTBMaskAccessor.h" namespace rttb { namespace testing { /*! @brief DVHCalculatorTest. Test if calculation in new architecture returns similar results to the original implementation. WARNING: The values for comparison need to be adjusted if the input files are changed! */ int DVHCalculatorExampleTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::DVHCalculator::MaskedDoseIteratorPointer MaskedDoseIteratorPointer; typedef masks::legacy::OTBMaskAccessor::StructTypePointer StructTypePointer; typedef core::DVH::DVHPointer DVHPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: structure file name // 2: dose1 file name // 3: dose2 file name // 4: dose3 file name std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string RTDOSE2_FILENAME; std::string RTDOSE3_FILENAME; if (argc > 1) { RTSTRUCT_FILENAME = argv[1]; } if (argc > 2) { RTDOSE_FILENAME = argv[2]; } if (argc > 3) { RTDOSE2_FILENAME = argv[3]; } if (argc > 4) { RTDOSE3_FILENAME = argv[4]; } OFCondition status; DcmFileFormat fileformat; // read dicom-rt dose io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create a vector of MaskAccessors (one for each structure) StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); //storeage for mask accessors to reduce time spent on voxelization (perform only once) std::vector rtStructSetMaskAccessorVec; ::DRTDoseIOD rtdose2; io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator2(RTDOSE2_FILENAME.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); ::DRTDoseIOD rtdose3; io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator3(RTDOSE3_FILENAME.c_str()); DoseAccessorPointer doseAccessor3(doseAccessorGenerator3.generateDoseAccessor()); //start evaluation clock_t start(clock()); if (rtStructureSet->getNumberOfStructures() > 0) { - for (int j = 0; j < rtStructureSet->getNumberOfStructures(); j++) + for (int j = 0; j < static_cast(rtStructureSet->getNumberOfStructures()); j++) { std::cout << rtStructureSet->getStructure(j)->getLabel() << std::endl; //create MaskAccessor for each structure ::boost::shared_ptr spOTBMaskAccessor = ::boost::make_shared(rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); spOTBMaskAccessor->updateMask(); MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor); //create corresponding MaskedDoseIterator ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(spMaskAccessor, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); //store MaskAccessor for each structure (later reuse) rtStructSetMaskAccessorVec.push_back(spMaskAccessor); //calculate DVH rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvh = *(calc.generateDVH()); //DEBUG OUTPUT std::cout << "=== Dose 1 Structure " << j << "===" << std::endl; std::cout << std::setprecision(20) << "max: " << dvh.getMaximum() << std::endl; std::cout << std::setprecision(20) << "min: " << dvh.getMinimum() << std::endl; std::cout << std::setprecision(20) << "mean: " << dvh.getMean() << std::endl; std::cout << std::setprecision(20) << "median: " << dvh.getMedian() << std::endl; std::cout << std::setprecision(20) << "modal: " << dvh.getModal() << std::endl; std::cout << std::setprecision(20) << "std: " << dvh.getStdDeviation() << std::endl; std::cout << std::setprecision(20) << "var: " << dvh.getVariance() << std::endl; std::cout << std::setprecision(20) << "numV: " << dvh.getNumberOfVoxels() << std::endl; //compare explicit values for some results. //expected values were generated from the original implementation if (j == 0) { CHECK_CLOSE(5.6658345820895519e-005, dvh.getMaximum(), 1e-3); CHECK_CLOSE(5.6658345820895519e-005, dvh.getMinimum(), errorConstant); CHECK_CLOSE(5.6658345820895525e-005, dvh.getMean(), errorConstant); CHECK_CLOSE(5.6658345820895519e-005, dvh.getMedian(), errorConstant); CHECK_CLOSE(5.6658345820895519e-005, dvh.getModal(), errorConstant); CHECK_CLOSE(-4.1359030627651384e-025, dvh.getVariance(), errorConstant); CHECK_CLOSE(89230.858052685173, dvh.getNumberOfVoxels(), 2e-1); } if (j == 1) { CHECK_CLOSE(-1.2407709188295415e-024, dvh.getVariance(), errorConstant); CHECK_CLOSE(595.30645355716683, dvh.getNumberOfVoxels(), 1e-3); } if (j == 2) { CHECK_CLOSE(-4.1359030627651384e-025, dvh.getVariance(), errorConstant); CHECK_CLOSE(1269.9125811291087, dvh.getNumberOfVoxels(), 1e-3); } if (j == 9) { CHECK_CLOSE(5.6658345820895519e-005, dvh.getMaximum(), 1e-3); CHECK_CLOSE(5.6658345820895519e-005, dvh.getMinimum(), errorConstant); CHECK_CLOSE(5.6658345820895519e-005, dvh.getMean(), errorConstant); CHECK_CLOSE(5.6658345820895519e-005, dvh.getMedian(), errorConstant); CHECK_CLOSE(5.6658345820895519e-005, dvh.getModal(), errorConstant); CHECK_CLOSE(0, dvh.getStdDeviation(), errorConstant); CHECK_CLOSE(0, dvh.getVariance(), errorConstant); CHECK_CLOSE(1328.4918116279605, dvh.getNumberOfVoxels(), 1e-3); } } } clock_t finish(clock()); std::cout << "DVH Calculation time: " << finish - start << " ms" << std::endl; clock_t start2(clock()); - for (int j = 0; j < rtStructSetMaskAccessorVec.size(); j++) + for (size_t j = 0; j < rtStructSetMaskAccessorVec.size(); j++) { //create corresponding MaskedDoseIterator ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(rtStructSetMaskAccessorVec.at(j), doseAccessor2); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); //calculate DVH rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor2->getUID()); rttb::core::DVH dvh = *(calc.generateDVH()); if (j == 0) { CHECK_CLOSE(1.8423272053074631, dvh.getMaximum(), 1e-1); CHECK_CLOSE(0.0069001018925373145, dvh.getMinimum(), errorConstant); CHECK_CLOSE(0.5534586388640208, dvh.getMean(), errorConstant); CHECK_CLOSE(0.42090621544477619, dvh.getMedian(), errorConstant); CHECK_CLOSE(0.075901120817910464, dvh.getModal(), errorConstant); CHECK_CLOSE(0.44688344565881616, dvh.getStdDeviation(), 1e-4); CHECK_CLOSE(0.19970481400389611, dvh.getVariance(), 1e-4); CHECK_CLOSE(89230.858052685187, dvh.getNumberOfVoxels(), 1e-1); } if (j == 4) { CHECK_CLOSE(1.6264736515373135, dvh.getMaximum(), 1e-3); CHECK_CLOSE(0.10433981915522389, dvh.getMinimum(), errorConstant); CHECK_CLOSE(0.70820073085773427, dvh.getMean(), errorConstant); CHECK_CLOSE(0.71810346124477609, dvh.getMedian(), errorConstant); CHECK_CLOSE(0.23936782041492538, dvh.getModal(), errorConstant); CHECK_CLOSE(0.36355099006842068, dvh.getStdDeviation(), errorConstant); CHECK_CLOSE(0.13216932237972889, dvh.getVariance(), errorConstant); CHECK_CLOSE(2299.7105030728999, dvh.getNumberOfVoxels(), 1e-3); } } clock_t finish2(clock()); std::cout << "Reset dose 2, DVH Calculation time: " << finish2 - start2 << " ms" << std::endl; clock_t start3(clock()); - for (int j = 0; j < rtStructSetMaskAccessorVec.size(); j++) + for (size_t j = 0; j < rtStructSetMaskAccessorVec.size(); j++) { //create corresponding MaskedDoseIterator ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(rtStructSetMaskAccessorVec.at(j), doseAccessor3); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); //calculate DVH rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor3->getUID()); rttb::core::DVH dvh = *(calc.generateDVH()); if (j == 1) { CHECK_CLOSE(0.0010765085705970151, dvh.getMaximum(), 1e-3); CHECK_CLOSE(0.00087641404074626872, dvh.getMinimum(), errorConstant); CHECK_CLOSE(0.0009788401527774486, dvh.getMean(), errorConstant); CHECK_CLOSE(0.00098846697746268666, dvh.getMedian(), errorConstant); CHECK_CLOSE(0.00098846697746268666, dvh.getModal(), errorConstant); CHECK_CLOSE(3.2977969280849566e-005, dvh.getStdDeviation(), errorConstant); CHECK_CLOSE(1.0875464578886574e-009, dvh.getVariance(), errorConstant); CHECK_CLOSE(595.30645355716683, dvh.getNumberOfVoxels(), 1e-3); } if (j == 6) { CHECK_CLOSE(0.0016589942782835824, dvh.getMaximum(), 1e-3); CHECK_CLOSE(0.00027960577723880602, dvh.getMinimum(), errorConstant); CHECK_CLOSE(0.0010389077643351956, dvh.getMean(), errorConstant); CHECK_CLOSE(0.0011246365706716419, dvh.getMedian(), errorConstant); CHECK_CLOSE(0.0013856019627611941, dvh.getModal(), errorConstant); CHECK_CLOSE(0.00036431958148461669, dvh.getStdDeviation(), errorConstant); CHECK_CLOSE(1.3272875745312625e-007, dvh.getVariance(), errorConstant); CHECK_CLOSE(8034.8878045085003, dvh.getNumberOfVoxels(), 1e-2); } } clock_t finish3(clock()); std::cout << "Reset dose 3, DVH Calculation time: " << finish3 - start3 << " ms" << std::endl; RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/examples/RTBioModelExampleTest.cpp b/testing/examples/RTBioModelExampleTest.cpp index c292641..f5950b9 100644 --- a/testing/examples/RTBioModelExampleTest.cpp +++ b/testing/examples/RTBioModelExampleTest.cpp @@ -1,416 +1,415 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include "litCheckMacros.h" #include "rttbBioModel.h" #include "rttbDVHTxtFileReader.h" #include "rttbDVH.h" #include "rttbTCPLQModel.h" #include "rttbNTCPLKBModel.h" #include "rttbNTCPRSModel.h" #include "rttbBioModelScatterPlots.h" #include "rttbBioModelCurve.h" #include "rttbDvhBasedModels.h" #include "rttbDoseIteratorInterface.h" namespace rttb { namespace testing { /*! @brief RTBioModelTest. TCP calculated using a DVH PTV and LQ Model. NTCP tested using 3 Normal Tissue DVHs and LKB/RS Model. Test if calculation in new architecture returns similar results to the original implementation. WARNING: The values for comparison need to be adjusted if the input files are changed! */ int RTBioModelExampleTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; typedef rttb::models::CurveDataType CurveDataType; typedef std::multimap > ScatterPlotType; typedef core::DVH::DVHPointer DVHPointer; //increased accuracy requires double values in the calculation (rttbBaseType.h) double toleranceEUD = 1e-5; double tolerance = 1e-7; //ARGUMENTS: 1: ptv dvh file name // 2: normal tissue 1 dvh file name // 3: normal tissue 2 dvh file name // 4: normal tissue 3 dvh file name //...........5: Virtuos MPM_LR_ah dvh lung file name //...........6: Virtuos MPM_LR_ah dvh target file name std::string DVH_FILENAME_PTV; std::string DVH_FILENAME_NT1; std::string DVH_FILENAME_NT2; std::string DVH_FILENAME_NT3; std::string DVH_FILENAME_TV_TEST; std::string DVH_Virtuos_Target; std::string DVH_Virtuos_Lung; if (argc > 1) { DVH_FILENAME_PTV = argv[1]; } if (argc > 2) { DVH_FILENAME_NT1 = argv[2]; } if (argc > 3) { DVH_FILENAME_NT2 = argv[3]; } if (argc > 4) { DVH_FILENAME_NT3 = argv[4]; } if (argc > 5) { DVH_FILENAME_TV_TEST = argv[5]; } if (argc > 6) { DVH_Virtuos_Lung = argv[6]; } if (argc > 7) { DVH_Virtuos_Target = argv[7]; } //DVH PTV rttb::io::other::DVHTxtFileReader dvhReader = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_PTV); DVHPointer dvhPtr = dvhReader.generateDVH(); CHECK_CLOSE(6.04759613161786830000e+001, models::getEUD(dvhPtr, 10), toleranceEUD); rttb::io::other::DVHTxtFileReader dvhReader_test_tv = rttb::io::other::DVHTxtFileReader( DVH_FILENAME_TV_TEST); DVHPointer dvh_test_tv = dvhReader_test_tv.generateDVH(); //test TCP LQ Model models::BioModelParamType alpha = 0.35; models::BioModelParamType beta = 0.023333333333333; models::BioModelParamType roh = 10000000; int numFractions = 2; DoseTypeGy normalizationDose = 68; rttb::models::TCPLQModel tcplq = rttb::models::TCPLQModel(dvhPtr, alpha, beta, roh, numFractions); CHECK_EQUAL(alpha, tcplq.getAlphaMean()); CHECK_EQUAL(alpha / beta, tcplq.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); CHECK_NO_THROW(tcplq.init()); if (tcplq.init()) { CHECK_CLOSE(1.00497232941856940000e-127, tcplq.getValue(), tolerance); } CurveDataType curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); CurveDataType::iterator it; for (it = curve.begin(); it != curve.end(); it++) { if ((*it).first < 72) { CHECK_EQUAL(0, (*it).second); } else if ((*it).first > 150) { CHECK((*it).second > 0.9); } } models::BioModelParamType alphaBeta = 10; tcplq.setParameters(alpha, alphaBeta, roh, 0.08); CHECK_EQUAL(alpha, tcplq.getAlphaMean()); CHECK_EQUAL(alphaBeta, tcplq.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); if (tcplq.init()) { CHECK_CLOSE(1.84e-005, tcplq.getValue(), tolerance); } normalizationDose = 40; curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); alpha = 1; alphaBeta = 14.5; tcplq.setAlpha(alpha); tcplq.setAlphaBeta(alphaBeta); tcplq.setRho(roh); CHECK_EQUAL(alpha, tcplq.getAlphaMean()); CHECK_EQUAL(alphaBeta, tcplq.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); if (tcplq.init()) { CHECK_CLOSE(0.954885, tcplq.getValue(), toleranceEUD); } alpha = 0.9; alphaBeta = 1; tcplq.setAlpha(alpha); tcplq.setAlphaBeta(alphaBeta); tcplq.setRho(roh); CHECK_EQUAL(alpha, tcplq.getAlphaMean()); CHECK_EQUAL(alphaBeta, tcplq.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); if (tcplq.init()) { CHECK_EQUAL(1, tcplq.getValue()); } //TCP LQ Test alpha = 0.3; beta = 0.03; roh = 10000000; numFractions = 20; rttb::models::TCPLQModel tcplq_test = rttb::models::TCPLQModel(dvh_test_tv, alpha, beta, roh, numFractions); CHECK_EQUAL(alpha, tcplq_test.getAlphaMean()); CHECK_EQUAL(alpha / beta, tcplq_test.getAlphaBeta()); CHECK_EQUAL(roh, tcplq_test.getRho()); CHECK_NO_THROW(tcplq_test.init()); if (tcplq_test.init()) { CHECK_CLOSE(9.79050278878883180000e-001, tcplq_test.getValue(), tolerance); } normalizationDose = 60; curve = models::getCurveDoseVSBioModel(tcplq_test, normalizationDose); //DVH HT 1 rttb::io::other::DVHTxtFileReader dvhReader2 = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_NT1); DVHPointer dvhPtr2 = dvhReader2.generateDVH(); CHECK_CLOSE(1.07920836034015810000e+001, models::getEUD(dvhPtr2, 10), toleranceEUD); //test RTNTCPLKBModel rttb::models::NTCPLKBModel lkb = rttb::models::NTCPLKBModel(); models::BioModelParamType aVal = 10; models::BioModelParamType mVal = 0.16; models::BioModelParamType d50Val = 55; CHECK_EQUAL(0, lkb.getA()); CHECK_EQUAL(0, lkb.getM()); CHECK_EQUAL(0, lkb.getD50()); lkb.setDVH(dvhPtr2); CHECK_EQUAL(dvhPtr2, lkb.getDVH()); lkb.setA(aVal); CHECK_EQUAL(aVal, lkb.getA()); lkb.setM(mVal); CHECK_EQUAL(mVal, lkb.getM()); lkb.setD50(d50Val); CHECK_EQUAL(d50Val, lkb.getD50()); CHECK_NO_THROW(lkb.init()); if (lkb.init()) { CHECK_CLOSE(2.53523522831366570000e-007, lkb.getValue(), tolerance); } //test RTNTCPRSModel rttb::models::NTCPRSModel rs = rttb::models::NTCPRSModel(); models::BioModelParamType gammaVal = 1.7; models::BioModelParamType sVal = 1; CHECK_EQUAL(0, rs.getGamma()); CHECK_EQUAL(0, rs.getS()); CHECK_EQUAL(0, rs.getD50()); rs.setDVH(dvhPtr2); CHECK_EQUAL(dvhPtr2, rs.getDVH()); rs.setD50(d50Val); CHECK_EQUAL(d50Val, rs.getD50()); rs.setGamma(gammaVal); CHECK_EQUAL(gammaVal, rs.getGamma()); rs.setS(sVal); CHECK_EQUAL(sVal, rs.getS()); CHECK_NO_THROW(rs.init()); if (rs.init()) { CHECK_CLOSE(3.70385888626145740000e-009, rs.getValue(), tolerance); } //DVH HT 2 rttb::io::other::DVHTxtFileReader dvhReader3 = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_NT2); DVHPointer dvhPtr3 = dvhReader3.generateDVH(); CHECK_CLOSE(1.26287047025885110000e+001, models::getEUD(dvhPtr3, 10), toleranceEUD); //test RTNTCPLKBModel aVal = 10; mVal = 0.16; d50Val = 55; lkb.setDVH(dvhPtr3); CHECK_EQUAL(dvhPtr3, lkb.getDVH()); lkb.setA(aVal); CHECK_EQUAL(aVal, lkb.getA()); lkb.setM(mVal); CHECK_EQUAL(mVal, lkb.getM()); lkb.setD50(d50Val); CHECK_EQUAL(d50Val, lkb.getD50()); if (lkb.init()) { CHECK_CLOSE(7.36294657754956700000e-007, lkb.getValue(), tolerance); } //test RTNTCPRSModel rs = rttb::models::NTCPRSModel(); gammaVal = 1.7; sVal = 1; CHECK_EQUAL(0, rs.getGamma()); CHECK_EQUAL(0, rs.getS()); CHECK_EQUAL(0, rs.getD50()); rs.setDVH(dvhPtr3); CHECK_EQUAL(dvhPtr3, rs.getDVH()); rs.setD50(d50Val); CHECK_EQUAL(d50Val, rs.getD50()); rs.setGamma(gammaVal); CHECK_EQUAL(gammaVal, rs.getGamma()); rs.setS(sVal); CHECK_EQUAL(sVal, rs.getS()); if (rs.init()) { CHECK_CLOSE(1.76778795490939440000e-007, rs.getValue(), tolerance); } //DVH HT 3 rttb::io::other::DVHTxtFileReader dvhReader4 = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_NT3); DVHPointer dvhPtr4 = dvhReader4.generateDVH(); CHECK_CLOSE(2.18212982041056310000e+001, models::getEUD(dvhPtr4, 10), toleranceEUD); //test RTNTCPLKBModel aVal = 10; mVal = 0.16; d50Val = 55; lkb.setDVH(dvhPtr4); CHECK_EQUAL(dvhPtr4, lkb.getDVH()); lkb.setA(aVal); CHECK_EQUAL(aVal, lkb.getA()); lkb.setM(mVal); CHECK_EQUAL(mVal, lkb.getM()); lkb.setD50(d50Val); CHECK_EQUAL(d50Val, lkb.getD50()); if (lkb.init()) { CHECK_CLOSE(8.15234192641929420000e-005, lkb.getValue(), tolerance); } //test RTNTCPRSModel rs = rttb::models::NTCPRSModel(); gammaVal = 1.7; sVal = 1; CHECK_EQUAL(0, rs.getGamma()); CHECK_EQUAL(0, rs.getS()); CHECK_EQUAL(0, rs.getD50()); rs.setDVH(dvhPtr4); CHECK_EQUAL(dvhPtr4, rs.getDVH()); rs.setD50(d50Val); CHECK_EQUAL(d50Val, rs.getD50()); rs.setGamma(gammaVal); CHECK_EQUAL(gammaVal, rs.getGamma()); rs.setS(sVal); CHECK_EQUAL(sVal, rs.getS()); if (rs.init()) { CHECK_CLOSE(2.02607985020919480000e-004, rs.getValue(), tolerance); } //test using Virtuos Pleuramesotheliom MPM_LR_ah //DVH PTV rttb::io::other::DVHTxtFileReader dR_Target = rttb::io::other::DVHTxtFileReader(DVH_Virtuos_Target); DVHPointer dvhPtrTarget = dR_Target.generateDVH(); rttb::io::other::DVHTxtFileReader dR_Lung = rttb::io::other::DVHTxtFileReader(DVH_Virtuos_Lung); DVHPointer dvhPtrLung = dR_Lung.generateDVH(); //test TCP LQ Model models::BioModelParamType alphaMean = 0.34; models::BioModelParamType alphaVarianz = 0.02; models::BioModelParamType alpha_beta = 28; models::BioModelParamType rho = 1200; int numFractionsVirtuos = 27; rttb::models::TCPLQModel tcplqVirtuos = rttb::models::TCPLQModel(dvhPtrTarget, rho, numFractionsVirtuos, alpha_beta, alphaMean, alphaVarianz); if (tcplqVirtuos.init()) { CHECK_CLOSE(0.8894, tcplqVirtuos.getValue(), 1e-4); } models::BioModelParamType d50Mean = 20; - models::BioModelParamType d50Varianz = 2; models::BioModelParamType m = 0.36; models::BioModelParamType a = 1.06; rttb::models::NTCPLKBModel lkbVirtuos = rttb::models::NTCPLKBModel(dvhPtrLung, d50Mean, m, a); if (lkbVirtuos.init()) { CHECK_CLOSE(0.0397, lkbVirtuos.getValue(), 1e-4); } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing -}//rttb \ No newline at end of file +}//rttb diff --git a/testing/examples/RTBioModelScatterPlotExampleTest.cpp b/testing/examples/RTBioModelScatterPlotExampleTest.cpp index 4cefe2f..aa96469 100644 --- a/testing/examples/RTBioModelScatterPlotExampleTest.cpp +++ b/testing/examples/RTBioModelScatterPlotExampleTest.cpp @@ -1,480 +1,479 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include "litCheckMacros.h" #include "rttbBioModel.h" #include "rttbDVHTxtFileReader.h" #include "rttbDVH.h" #include "rttbTCPLQModel.h" #include "rttbNTCPLKBModel.h" #include "rttbNTCPRSModel.h" #include "rttbBioModelScatterPlots.h" #include "rttbBioModelCurve.h" #include "rttbDvhBasedModels.h" #include "../models/rttbScatterTester.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace testing { /*! @brief RTBioModelScatterPlotExampleTest. calculating Curves and Scatterplots for TCP and NTCP models. The values on curve and scatterplot need to be similar for similar dose values. The range of difference is given by the variance used to generate the scatter. WARNING: The values for comparison need to be adjusted if the input files are changed! */ int RTBioModelScatterPlotExampleTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; typedef rttb::models::CurveDataType CurveDataType; typedef rttb::models::ScatterPlotType ScatterPlotType; typedef core::DVH::DVHPointer DVHPointer; //increased accuracy requires double values in the calculation (rttbBaseType.h) double toleranceEUD = 1e-5; - double tolerance = 1e-7; //ARGUMENTS: 1: ptv dvh file name // 2: normal tissue 1 dvh file name // 3: TV dvh file name std::string DVH_FILENAME_PTV; std::string DVH_FILENAME_NT1; std::string DVH_FILENAME_TV_TEST; if (argc > 1) { DVH_FILENAME_PTV = argv[1]; } if (argc > 2) { DVH_FILENAME_NT1 = argv[2]; } if (argc > 3) { DVH_FILENAME_TV_TEST = argv[3]; } //DVH PTV rttb::io::other::DVHTxtFileReader dvhReader = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_PTV); DVHPointer dvhPtr = dvhReader.generateDVH(); rttb::io::other::DVHTxtFileReader dvhReader_test_tv = rttb::io::other::DVHTxtFileReader( DVH_FILENAME_TV_TEST); DVHPointer dvh_test_tv = dvhReader_test_tv.generateDVH(); //test TCP LQ Model models::BioModelParamType alpha = 0.35; models::BioModelParamType beta = 0.023333333333333; models::BioModelParamType roh = 10000000; int numFractions = 2; DoseTypeGy normalizationDose = 68; rttb::models::TCPLQModel tcplq = rttb::models::TCPLQModel(dvhPtr, alpha, beta, roh, numFractions); CHECK_NO_THROW(tcplq.init()); CurveDataType curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq, 0, alpha, 0, normalizationDose, 100)); ScatterPlotType tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 0, alpha, 0, normalizationDose, 100); CHECK_EQUAL(100, tcpScatter.size()); ScatterTester scatterCompare(curve, tcpScatter); CHECK_TESTER(scatterCompare); //test also with other parameter tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 3, roh, 0, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); CHECK_TESTER(scatterCompare); std::vector paramIdVec; models::BioModel::ParamVectorType meanVec; models::BioModel::ParamVectorType meanVecTest; meanVecTest.push_back(alpha); models::BioModel::ParamVectorType varianceVec; //"alphaMean":0,"alphaVariance":1,"alpha_beta":2, "rho":3 paramIdVec.push_back(0); meanVec.push_back(tcplq.getAlphaMean()); varianceVec.push_back(0); //setting parameter 1 will change the resulting scatter plot dramatically - is it meant to? //this is unexpected since the value was taken from the original calculation //paramIdVec.push_back(1); meanVec.push_back(tcplq.getAlphaVariance()); varianceVec.push_back(0); paramIdVec.push_back(2); meanVec.push_back(tcplq.getAlphaBeta()); varianceVec.push_back(0); paramIdVec.push_back(3); meanVec.push_back(tcplq.getRho()); varianceVec.push_back(0); CHECK_THROW_EXPLICIT(models::getScatterPlotVaryParameters(tcplq, paramIdVec, meanVecTest, varianceVec, normalizationDose, 50), core::InvalidParameterException); ScatterPlotType scatterVary = models::getScatterPlotVaryParameters(tcplq, paramIdVec, meanVec, varianceVec, normalizationDose, 50); CHECK_EQUAL(50, scatterVary.size()); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); models::BioModelParamType variance = 0.00015; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq, 0, alpha, variance, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 0, alpha, variance, normalizationDose, 100); scatterCompare.setVariance(variance); scatterCompare.setCompareScatter(tcpScatter); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); //test also with other parameter tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 3, roh, variance, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(tcplq.getAlphaMean()); varianceVec.push_back(variance); //paramIdVec.push_back(1); meanVec.push_back(tcplq.getAlphaVariance()); varianceVec.push_back(variance); paramIdVec.push_back(2); meanVec.push_back(tcplq.getAlphaBeta()); varianceVec.push_back(variance); paramIdVec.push_back(3); meanVec.push_back(tcplq.getRho()); varianceVec.push_back(variance); scatterVary = models::getScatterPlotVaryParameters(tcplq, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); models::BioModelParamType alphaBeta = 10; tcplq.setParameters(alpha, alphaBeta, roh, 0.08); tcplq.init(); normalizationDose = 40; curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq, 0, alpha, 0, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 0, alpha, 0, normalizationDose, 100); scatterCompare.setReferenceCurve(curve); scatterCompare.setVariance(0); //do not allow larger deviations scatterCompare.setAllowExceptions(false); scatterCompare.setCompareScatter(tcpScatter); CHECK_TESTER(scatterCompare); variance = 0.25; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq, 0, alpha, variance, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 0, alpha, variance, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); scatterCompare.setVariance(variance); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); /*TCP LQ Test*/ alpha = 0.3; beta = 0.03; roh = 10000000; numFractions = 20; rttb::models::TCPLQModel tcplq_test = rttb::models::TCPLQModel(dvh_test_tv, alpha, beta, roh, numFractions); CHECK_NO_THROW(tcplq_test.init()); normalizationDose = 60; curve = models::getCurveDoseVSBioModel(tcplq_test, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq_test, 0, alpha, 0, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq_test, 0, alpha, 0, normalizationDose, 100); scatterCompare.setReferenceCurve(curve); scatterCompare.setVariance(0); scatterCompare.setCompareScatter(tcpScatter); CHECK_TESTER(scatterCompare); //test also with other parameter tcpScatter = models::getScatterPlotVary1Parameter(tcplq_test, 3, roh, 0, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(tcplq_test.getAlphaMean()); varianceVec.push_back(0); //paramIdVec.push_back(1); meanVec.push_back(tcplq_test.getAlphaVariance()); varianceVec.push_back(0); paramIdVec.push_back(2); meanVec.push_back(tcplq_test.getAlphaBeta()); varianceVec.push_back(0); paramIdVec.push_back(3); meanVec.push_back(tcplq_test.getRho()); varianceVec.push_back(0); scatterVary = models::getScatterPlotVaryParameters(tcplq_test, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); variance = 0.00025; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq_test, 0, alpha, variance, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq_test, 0, alpha, variance, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); scatterCompare.setVariance(variance); CHECK_TESTER(scatterCompare); //test also with other parameter tcpScatter = models::getScatterPlotVary1Parameter(tcplq_test, 3, roh, variance, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); scatterCompare.setAllowExceptions(false); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(tcplq_test.getAlphaMean()); varianceVec.push_back(variance); //paramIdVec.push_back(1); meanVec.push_back(tcplq_test.getAlphaVariance()); varianceVec.push_back(variance); paramIdVec.push_back(2); meanVec.push_back(tcplq_test.getAlphaBeta()); varianceVec.push_back(variance); paramIdVec.push_back(3); meanVec.push_back(tcplq_test.getRho()); varianceVec.push_back(variance); scatterVary = models::getScatterPlotVaryParameters(tcplq_test, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); //DVH HT 1 rttb::io::other::DVHTxtFileReader dvhReader2 = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_NT1); DVHPointer dvhPtr2 = dvhReader2.generateDVH(); CHECK_CLOSE(1.07920836034015810000e+001, models::getEUD(dvhPtr2, 10), toleranceEUD); //test RTNTCPLKBModel rttb::models::NTCPLKBModel lkb = rttb::models::NTCPLKBModel(); models::BioModelParamType aVal = 10; models::BioModelParamType mVal = 0.16; models::BioModelParamType d50Val = 55; lkb.setDVH(dvhPtr2); lkb.setA(aVal); lkb.setM(mVal); lkb.setD50(d50Val); CHECK_NO_THROW(lkb.init()); normalizationDose = 60; curve = models::getCurveDoseVSBioModel(lkb, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(lkb, 2, aVal, 0, normalizationDose, 100)); ScatterPlotType scatter = models::getScatterPlotVary1Parameter(lkb, 2, aVal, 0, normalizationDose, 100); scatterCompare.setReferenceCurve(curve); scatterCompare.setVariance(0); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); //"d50":0,"m":1,"a":2 //test also with other parameter scatter = models::getScatterPlotVary1Parameter(lkb, 0, d50Val, 0, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(lkb.getD50()); varianceVec.push_back(0); paramIdVec.push_back(1); meanVec.push_back(lkb.getM()); varianceVec.push_back(0); paramIdVec.push_back(2); meanVec.push_back(lkb.getA()); varianceVec.push_back(0); scatterVary = models::getScatterPlotVaryParameters(lkb, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); variance = 0.00025; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(lkb, 2, aVal, variance, normalizationDose, 100)); scatter = models::getScatterPlotVary1Parameter(lkb, 2, aVal, variance, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); scatterCompare.setVariance(variance); CHECK_TESTER(scatterCompare); //test also with other parameter scatter = models::getScatterPlotVary1Parameter(lkb, 0, d50Val, variance, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(lkb.getD50()); varianceVec.push_back(variance); paramIdVec.push_back(1); meanVec.push_back(lkb.getM()); varianceVec.push_back(variance); paramIdVec.push_back(2); meanVec.push_back(lkb.getA()); varianceVec.push_back(variance); scatterVary = models::getScatterPlotVaryParameters(lkb, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); //test RTNTCPRSModel rttb::models::NTCPRSModel rs = rttb::models::NTCPRSModel(); models::BioModelParamType gammaVal = 1.7; models::BioModelParamType sVal = 1; rs.setDVH(dvhPtr2); rs.setD50(d50Val); rs.setGamma(gammaVal); rs.setS(sVal); CHECK_NO_THROW(rs.init()); normalizationDose = 60; curve = models::getCurveDoseVSBioModel(rs, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(rs, 0, d50Val, 0, normalizationDose, 100)); scatter = models::getScatterPlotVary1Parameter(rs, 0, d50Val, 0, normalizationDose, 100); scatterCompare.setReferenceCurve(curve); scatterCompare.setVariance(0); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); //"d50":0,"gamma":1,"s":2 //test also with other parameter scatter = models::getScatterPlotVary1Parameter(rs, 1, gammaVal, 0, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(rs.getD50()); varianceVec.push_back(0); paramIdVec.push_back(1); meanVec.push_back(rs.getGamma()); varianceVec.push_back(0); paramIdVec.push_back(2); meanVec.push_back(rs.getS()); varianceVec.push_back(0); scatterVary = models::getScatterPlotVaryParameters(rs, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); variance = 0.0075; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(rs, 0, d50Val, variance, normalizationDose, 100)); scatter = models::getScatterPlotVary1Parameter(rs, 0, d50Val, variance, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); scatterCompare.setVariance(variance); CHECK_TESTER(scatterCompare); //test also with other parameter scatter = models::getScatterPlotVary1Parameter(rs, 2, sVal, variance, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(rs.getD50()); varianceVec.push_back(variance); paramIdVec.push_back(1); meanVec.push_back(rs.getGamma()); varianceVec.push_back(variance); paramIdVec.push_back(2); meanVec.push_back(rs.getS()); varianceVec.push_back(variance); scatterVary = models::getScatterPlotVaryParameters(rs, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing -}//rttb \ No newline at end of file +}//rttb diff --git a/testing/examples/RTDVHTest.cpp b/testing/examples/RTDVHTest.cpp index 233d9b4..2d9cbb8 100644 --- a/testing/examples/RTDVHTest.cpp +++ b/testing/examples/RTDVHTest.cpp @@ -1,115 +1,109 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include "litCheckMacros.h" #include "rttbDVH.h" #include "rttbDVHTxtFileReader.h" #include "rttbBaseType.h" #include "rttbDvhBasedModels.h" namespace rttb { namespace testing { /*! @brief RTDVHTest. Max, min, mean, modal, median, Vx, Dx, EUD, BED, LQED2 are tested. Test if calculation in new architecture returns similar results to the original implementation. Comparison of actual DVH values is performed in DVHCalculatorComparisonTest.cpp. WARNING: The values for comparison need to be adjusted if the input files are changed! */ int RTDVHTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: dvh file name std::string RTDVH_FILENAME_PTV; if (argc > 1) { RTDVH_FILENAME_PTV = argv[1]; } typedef core::DVH::DVHPointer DVHPointer; /*test RT dvh*/ rttb::io::other::DVHTxtFileReader dvhReader = rttb::io::other::DVHTxtFileReader(RTDVH_FILENAME_PTV); const DoseCalcType expectedValue = 0.01305; //dvhReader DVHPointer dvh = dvhReader.generateDVH(); CHECK_CLOSE(expectedValue, models::getEUD(dvh, 10), errorConstant); std::cout << models::getEUD(dvh, 10) << std::endl; std::map bedMap = models::calcBEDDVH(dvh, 15, 15); std::map LqedMap = models::calcLQED2DVH(dvh, 15, 10); CHECK_NO_THROW(models::calcBEDDVH(dvh, 15, 15)); CHECK_NO_THROW(models::calcLQED2DVH(dvh, 15, 10)); CHECK_NO_THROW(dvh->getDataDifferential(true)); CHECK_EQUAL(1, dvh->calcCumulativeDVH(true).at(0)); CHECK_NO_THROW(models::calcBEDDVH(dvh, 15, 15, true)); CHECK_NO_THROW(models::calcLQED2DVH(dvh, 15, 10, true)); - //test statistics (relative comulative data) + //test statistics (relative cumulative data) CHECK_CLOSE(expectedValue, dvh->getMaximum(), errorConstant); CHECK_CLOSE(expectedValue, dvh->getMinimum(), errorConstant); CHECK_CLOSE(expectedValue, dvh->getMean(), errorConstant); CHECK_CLOSE(expectedValue, dvh->getMedian(), errorConstant); CHECK_CLOSE(expectedValue, dvh->getModal(), errorConstant); CHECK_EQUAL(0, dvh->getVx(0.014)); CHECK_EQUAL(0.125, dvh->getVx(0.01)); - rttb::DoseTypeGy dTest = dvh->getDx(100); CHECK_CLOSE(0.0131, dvh->getDx(100), errorConstant + errorConstant * 10); CHECK_CLOSE(0.013, dvh->getDx(249), errorConstant); - dTest = dvh->getDx(249); CHECK_EQUAL(0, dvh->getDx(251)); - dTest = dvh->getDx(251); //test statistics (absolute comulative data) CHECK_EQUAL(2000, dvh->calcCumulativeDVH(false).at(0)); CHECK_EQUAL(0, dvh->getVx(0.014)); CHECK_EQUAL(250, dvh->getVx(0.01)); - dTest = dvh->getDx(100); CHECK_CLOSE(0.0131, dvh->getDx(100), errorConstant + errorConstant * 10); CHECK_CLOSE(0.013, dvh->getDx(249), errorConstant); - dTest = dvh->getDx(249); CHECK_EQUAL(0, dvh->getDx(251)); - dTest = dvh->getDx(251); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/examples/VoxelizationValidationTest.cpp b/testing/examples/VoxelizationValidationTest.cpp index 8b33bb8..d26ebc7 100644 --- a/testing/examples/VoxelizationValidationTest.cpp +++ b/testing/examples/VoxelizationValidationTest.cpp @@ -1,201 +1,201 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVHCalculator.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbOTBMaskAccessor.h" #include "rttbDVHTxtFileReader.h" #include "rttbBoostMaskAccessor.h" #include "rttbITKImageMaskAccessorConverter.h" #include "rttbImageWriter.h" namespace rttb { namespace testing { /*! @brief VoxelizationValidationTest. Compare two differnt voxelizations: OTB and Boost. Check dvh maximum and minimum for each structure. Check write mask to itk file for further validation. */ int VoxelizationValidationTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: structure file name // 2: dose1 file name // 3: directory name to write boost mask of all structures // 4: directory name to write OTB mask of all structures std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string BoostMask_DIRNAME; std::string OTBMask_DIRNAME; if (argc > 4) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; BoostMask_DIRNAME = argv[3]; OTBMask_DIRNAME = argv[4]; } OFCondition status; DcmFileFormat fileformat; /* read dicom-rt dose */ io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create a vector of MaskAccessors (one for each structure) StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); std::vector rtStructSetMaskAccessorVec; if (rtStructureSet->getNumberOfStructures() > 0) { - for (int j = 0; j < rtStructureSet->getNumberOfStructures(); j++) + for (size_t j = 0; j < rtStructureSet->getNumberOfStructures(); j++) { std::cout << j << ": " << rtStructureSet->getStructure(j)->getLabel() << std::endl; clock_t start(clock()); //create OTB MaskAccessor ::boost::shared_ptr spOTBMaskAccessor = ::boost::make_shared(rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); spOTBMaskAccessor->updateMask(); MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor); ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(spMaskAccessor, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvh = *(calc.generateDVH()); clock_t finish(clock()); std::cout << "OTB Mask Calculation time: " << finish - start << " ms" << std::endl; //Write the mask image to a file. /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! */ /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter(spOTBMaskAccessor); CHECK(itkConverter.process()); std::stringstream fileNameSstr; fileNameSstr << OTBMask_DIRNAME << j << ".mhd"; rttb::io::itk::ImageWriter writer(fileNameSstr.str(), itkConverter.getITKImage()); CHECK(writer.writeFile());*/ clock_t start2(clock()); //create Boost MaskAccessor MaskAccessorPointer boostMaskAccessorPtr = ::boost::make_shared (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); CHECK_NO_THROW(boostMaskAccessorPtr->updateMask()); ::boost::shared_ptr spMaskedDoseIteratorTmp2 = ::boost::make_shared(boostMaskAccessorPtr, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator2(spMaskedDoseIteratorTmp2); rttb::core::DVHCalculator calc2(spMaskedDoseIterator2, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvh2 = *(calc2.generateDVH()); clock_t finish2(clock()); std::cout << "Boost Mask Calculation and write file time: " << finish2 - start2 << " ms" << std::endl; //Write the mask image to a file. /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! */ /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter2(boostMaskAccessorPtr); CHECK(itkConverter2.process()); std::stringstream fileNameSstr2; fileNameSstr2 << BoostMask_DIRNAME << j << ".mhd"; rttb::io::itk::ImageWriter writer2(fileNameSstr2.str(), itkConverter2.getITKImage()); CHECK(writer2.writeFile());*/ //check close of 2 voxelizatin: OTB and Boost CHECK_CLOSE(dvh.getMaximum(), dvh2.getMaximum(), 0.1); CHECK_CLOSE(dvh.getMinimum(), dvh2.getMinimum(), 0.1); if (j != 7) //7: Ref.Pkt, mean = -1.#IND { CHECK_CLOSE(dvh.getMean(), dvh2.getMean(), 0.1); } CHECK_CLOSE(dvh.getMedian(), dvh2.getMedian(), 0.1); CHECK_CLOSE(dvh.getModal(), dvh2.getModal(), 0.1); //0: Aussenkontur and 3: Niere li. failed. if (j != 0 && j != 3) { CHECK_CLOSE(dvh.getVx(0), dvh2.getVx(0), dvh.getVx(0) * 0.05); //check volume difference < 5% } } } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/interpolation/InterpolationITKTransformation/SimpleMappableDoseAccessorWithITKTest.cpp b/testing/interpolation/InterpolationITKTransformation/SimpleMappableDoseAccessorWithITKTest.cpp index 1de49a1..5606149 100644 --- a/testing/interpolation/InterpolationITKTransformation/SimpleMappableDoseAccessorWithITKTest.cpp +++ b/testing/interpolation/InterpolationITKTransformation/SimpleMappableDoseAccessorWithITKTest.cpp @@ -1,192 +1,192 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbGenericDoseIterator.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbSimpleMappableDoseAccessor.h" #include "rttbNearestNeighborInterpolation.h" #include "rttbLinearInterpolation.h" #include "rttbTransformationInterface.h" #include "rttbITKTransformation.h" #include "rttbNullPointerException.h" #include "itkTranslationTransform.h" namespace rttb { namespace testing { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef rttb::interpolation::SimpleMappableDoseAccessor SimpleMappableDoseAccessor; typedef rttb::interpolation::NearestNeighborInterpolation NearestNeighborInterpolation; typedef rttb::interpolation::LinearInterpolation LinearInterpolation; typedef rttb::interpolation::TransformationInterface TransformationInterface; typedef rttb::interpolation::ITKTransformation ITKTransformation; typedef itk::TranslationTransform TranslationTransformType; /*! @brief SimpleMappableDoseAccessorWithITKTest - test the API of SimpleMappableDoseAccessor with ITK transformation 1) Test constructor 2) test getDoseAt() a) with Identity transform b) with translation transform */ int SimpleMappableDoseAccessorWithITKTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string RTDOSE_FILENAME_CONSTANT_TWO; std::string RTDOSE_FILENAME_INCREASE_X; if (argc > 2) { RTDOSE_FILENAME_CONSTANT_TWO = argv[1]; RTDOSE_FILENAME_INCREASE_X = argv[2]; } else { std::cout << "at least two parameters for SimpleMappableDoseAccessorWithITKTest are expected" << std::endl; return -1; } const double doseGridScaling = 2.822386e-005; rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1( RTDOSE_FILENAME_CONSTANT_TWO.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator2( RTDOSE_FILENAME_INCREASE_X.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); core::GeometricInfo doseAccessor1GeometricInfo = doseAccessor1->getGeometricInfo(); NearestNeighborInterpolation::Pointer interpolationNN = NearestNeighborInterpolation::Pointer(new NearestNeighborInterpolation()); LinearInterpolation::Pointer interpolationLinear = LinearInterpolation::Pointer (new LinearInterpolation()); //auto transformITKIdentity = TransformationInterface::Pointer(new ITKTransformation( // IdentityTransformType::New())); TranslationTransformType::Pointer transformITKIdentityTemporary = TranslationTransformType::New(); TranslationTransformType::OutputVectorType translationIdentity; translationIdentity[0] = 0.0; translationIdentity[1] = 0.0; translationIdentity[2] = 0.0; transformITKIdentityTemporary->Translate(translationIdentity); TransformationInterface::Pointer transformITKIdentity = TransformationInterface::Pointer( new ITKTransformation( transformITKIdentityTemporary)); TranslationTransformType::Pointer transformITKTranslationTemporary = TranslationTransformType::New(); TranslationTransformType::OutputVectorType translation; translation[0] = 5.0; translation[1] = 5.0; translation[2] = 5.0; transformITKTranslationTemporary->Translate(translation); TransformationInterface::Pointer transformITKTranslation = TransformationInterface::Pointer( new ITKTransformation( transformITKTranslationTemporary)); SimpleMappableDoseAccessor::Pointer aSimpleMappableDoseAccessorITKIdentity = SimpleMappableDoseAccessor::Pointer( new SimpleMappableDoseAccessor( doseAccessor1GeometricInfo, doseAccessor2, transformITKIdentity, interpolationLinear)); SimpleMappableDoseAccessor::Pointer aSimpleMappableDoseAccessorITKTranslation = SimpleMappableDoseAccessor::Pointer( new SimpleMappableDoseAccessor( doseAccessor1GeometricInfo, doseAccessor2, transformITKTranslation, interpolationLinear)); //1) Test constructor CHECK_NO_THROW(SimpleMappableDoseAccessor( doseAccessor1GeometricInfo, doseAccessor2, transformITKIdentity, interpolationLinear)); CHECK_NO_THROW(SimpleMappableDoseAccessor( doseAccessor1GeometricInfo, doseAccessor2, transformITKTranslation, interpolationLinear)); CHECK_NO_THROW(ITKTransformation(transformITKTranslationTemporary.GetPointer())); CHECK_THROW_EXPLICIT(ITKTransformation(NULL), core::NullPointerException); //2) test getDoseAt() // a) with Identity transform //test valid voxels std::vector voxelsAsIndexToTest; std::vector expectedValues; voxelsAsIndexToTest.push_back(VoxelGridIndex3D(5, 9, 8)); voxelsAsIndexToTest.push_back(VoxelGridIndex3D(22, 15, 10)); voxelsAsIndexToTest.push_back(VoxelGridIndex3D(30, 20, 7)); expectedValues.push_back(5.0 * doseGridScaling); expectedValues.push_back(22.0 * doseGridScaling); expectedValues.push_back(30.0 * doseGridScaling); //convert VoxelGridIndex3D to VoxelGridID - for (int i = 0; i < voxelsAsIndexToTest.size(); i++) + for (size_t i = 0; i < voxelsAsIndexToTest.size(); i++) { VoxelGridID currentId; doseAccessor1GeometricInfo.convert(voxelsAsIndexToTest.at(i), currentId); //test if the expected interpolation values are computed CHECK_EQUAL(aSimpleMappableDoseAccessorITKIdentity->getValueAt(voxelsAsIndexToTest.at(i)), expectedValues.at(i)); //test if getDoseAt(VoxelGridIndex3D) and getDoseAt(VoxelGridD) lead to the same results CHECK_EQUAL(aSimpleMappableDoseAccessorITKIdentity->getValueAt(voxelsAsIndexToTest.at(i)), aSimpleMappableDoseAccessorITKIdentity->getValueAt(currentId)); } //no tests with invalid IDs, this has been done already tested in SimpleMappableDoseAccessorTest //b) with translation transform //translation of one voxel in each direction expectedValues.clear(); expectedValues.push_back(6.0 * doseGridScaling); expectedValues.push_back(23.0 * doseGridScaling); expectedValues.push_back(31.0 * doseGridScaling); - for (int i = 0; i < voxelsAsIndexToTest.size(); i++) + for (size_t i = 0; i < voxelsAsIndexToTest.size(); i++) { VoxelGridID currentId; doseAccessor1GeometricInfo.convert(voxelsAsIndexToTest.at(i), currentId); //test if the expected interpolation values are computed CHECK_EQUAL(aSimpleMappableDoseAccessorITKTranslation->getValueAt(voxelsAsIndexToTest.at(i)), expectedValues.at(i)); //test if getDoseAt(VoxelGridIndex3D) and getDoseAt(VoxelGridD) lead to the same results CHECK_EQUAL(aSimpleMappableDoseAccessorITKTranslation->getValueAt(voxelsAsIndexToTest.at(i)), aSimpleMappableDoseAccessorITKTranslation->getValueAt(currentId)); } RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/io/dicom/DicomFileReaderHelperTest.cpp b/testing/io/dicom/DicomFileReaderHelperTest.cpp index 3ad5651..208edbc 100644 --- a/testing/io/dicom/DicomFileReaderHelperTest.cpp +++ b/testing/io/dicom/DicomFileReaderHelperTest.cpp @@ -1,98 +1,98 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDicomFileReaderHelper.h" namespace rttb { namespace testing { /*!@brief DicomDoseAccessorGeneratorTest - test the generators for dicom data 1) test getFileNamesWithSameUID() with a directory name 2) test getFileNames() with a RTDOSE file name and check equal with getFileNamesWithSameUID() 3) test getFileNames() with a RTSTRUCT file name */ int DicomFileReaderHelperTest(int argc, char* argv[]) { typedef boost::shared_ptr DRTDoseIODPtr; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: // 1: helax directory name // 2: dose file name // 3: structure file name std::string RT_DIRNAME; std::string RTDOSE_FILENAME; std::string RTStr_FILENAME; if (argc > 3) { RT_DIRNAME = argv[1]; RTDOSE_FILENAME = argv[2]; RTStr_FILENAME = argv[3]; } rttb::io::dicom::Modality doseModality = {rttb::io::dicom::Modality::RTDOSE}; rttb::io::dicom::Modality strModality = {rttb::io::dicom::Modality::RTSTRUCT}; //1) test getFileNamesWithSameUID() with a directory name std::vector fileVector = rttb::io::dicom::getFileNamesWithSameUID(RT_DIRNAME, doseModality); CHECK_EQUAL(fileVector.size(), 52); //2) test getFileNames() with a RTDOSE file name and check equal with getFileNamesWithSameUID() std::vector fileVector2 = rttb::io::dicom::getFileNames(RTDOSE_FILENAME); - for (int i = 0; i < fileVector.size(); i++) + for (size_t i = 0; i < fileVector.size(); i++) { CHECK_EQUAL(fileVector.at(i), fileVector2.at(i)); } //3) test getFileNames() with a RTSTRUCT file name fileVector2 = rttb::io::dicom::getFileNames(RTStr_FILENAME); CHECK_EQUAL(fileVector2.size(), 1); fileVector = rttb::io::dicom::getFileNamesWithSameUID(RT_DIRNAME, strModality); CHECK_EQUAL(fileVector.size(), 1); CHECK_EQUAL(fileVector.at(0), fileVector2.at(0)); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/io/itk/ITKMaskAccessorConverterTest.cpp b/testing/io/itk/ITKMaskAccessorConverterTest.cpp index dca4301..49da530 100644 --- a/testing/io/itk/ITKMaskAccessorConverterTest.cpp +++ b/testing/io/itk/ITKMaskAccessorConverterTest.cpp @@ -1,180 +1,180 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include "litCheckMacros.h" #include "litImageTester.h" #include "litTestImageIO.h" #include "itkImage.h" #include "itkImageFileReader.h" #include "rttbBaseType.h" #include "rttbDoseIteratorInterface.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomDoseAccessor.h" #include "rttbInvalidDoseException.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbITKImageMaskAccessorConverter.h" #include "rttbITKImageFileMaskAccessorGenerator.h" #include "rttbOTBMaskAccessor.h" #include "rttbBoostMaskAccessor.h" namespace rttb { namespace testing { /*!@brief MaskAccessorConverterTest - test the conversion rttb dose accessor ->itk 1) test with dicom file (DicomDoseAccessorGenerator) 2) test with mhd file (ITKImageFileDoseAccessorGenerator) */ int ITKMaskAccessorConverterTest(int argc, char* argv[]) { typedef core::DoseIteratorInterface::DoseAccessorPointer DoseAccessorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; typedef core::MaskAccessorInterface::MaskAccessorPointer MaskAccessorPointer; typedef io::itk::ITKImageMaskAccessor::ITKMaskImageType::Pointer ITKImageTypePointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: //ARGUMENTS: 1: structure file name // 2: dose1 file name std::string RTStr_FILENAME; std::string RTDose_FILENAME; std::string Mask_FILENAME; if (argc > 3) { RTStr_FILENAME = argv[1]; RTDose_FILENAME = argv[2]; Mask_FILENAME = argv[3]; } //1) Read dicomFile and test getITKImage() io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDose_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTStr_FILENAME.c_str()).generateStructureSet(); MaskAccessorPointer maskAccessorPtr = boost::make_shared (rtStructureSet->getStructure(9), doseAccessor1->getGeometricInfo()); maskAccessorPtr->updateMask();//!Important: Update the mask before conversion. io::itk::ITKImageMaskAccessorConverter maskAccessorConverter(maskAccessorPtr); CHECK_NO_THROW(maskAccessorConverter.process()); CHECK_NO_THROW(maskAccessorConverter.getITKImage()); //2) Read itk image, generate mask and convert it back to itk image, check equal MaskAccessorPointer maskAccessorPtr2 = io::itk::ITKImageFileMaskAccessorGenerator( Mask_FILENAME.c_str()).generateMaskAccessor(); maskAccessorPtr2->updateMask();//!Important: Update the mask before conversion. io::itk::ITKImageMaskAccessorConverter maskAccessorConverter2(maskAccessorPtr2); maskAccessorConverter2.process(); typedef itk::Image< DoseTypeGy, 3 > MaskImageType; typedef itk::ImageFileReader ReaderType; ITKImageTypePointer convertedImagePtr = maskAccessorConverter2.getITKImage(); io::itk::ITKImageMaskAccessor::ITKMaskImageType::Pointer expectedImage = lit::TestImageIO::readImage( Mask_FILENAME); CHECK_EQUAL(convertedImagePtr->GetOrigin()[0], expectedImage->GetOrigin()[0]); CHECK_EQUAL(convertedImagePtr->GetOrigin()[1], expectedImage->GetOrigin()[1]); CHECK_EQUAL(convertedImagePtr->GetOrigin()[2], expectedImage->GetOrigin()[2]); CHECK_EQUAL(convertedImagePtr->GetSpacing()[0], expectedImage->GetSpacing()[0]); CHECK_EQUAL(convertedImagePtr->GetSpacing()[1], expectedImage->GetSpacing()[1]); CHECK_EQUAL(convertedImagePtr->GetSpacing()[2], expectedImage->GetSpacing()[2]); CHECK_EQUAL(convertedImagePtr->GetLargestPossibleRegion().GetSize()[0], expectedImage->GetLargestPossibleRegion().GetSize()[0]); CHECK_EQUAL(convertedImagePtr->GetLargestPossibleRegion().GetSize()[1], expectedImage->GetLargestPossibleRegion().GetSize()[1]); CHECK_EQUAL(convertedImagePtr->GetLargestPossibleRegion().GetSize()[2], expectedImage->GetLargestPossibleRegion().GetSize()[2]); - int sizeX = convertedImagePtr->GetLargestPossibleRegion().GetSize()[0]; - int sizeY = convertedImagePtr->GetLargestPossibleRegion().GetSize()[1]; - int sizeZ = convertedImagePtr->GetLargestPossibleRegion().GetSize()[2]; + unsigned int sizeX = convertedImagePtr->GetLargestPossibleRegion().GetSize()[0]; + unsigned int sizeY = convertedImagePtr->GetLargestPossibleRegion().GetSize()[1]; + unsigned int sizeZ = convertedImagePtr->GetLargestPossibleRegion().GetSize()[2]; io::itk::ITKImageMaskAccessor::ITKMaskImageType::IndexType index; for (unsigned int i = 0; i < 20 && i < sizeX && i < sizeY && i < sizeZ; i++) { index[0] = i; index[1] = i; index[2] = i; if (expectedImage->GetPixel(index) >= 0 && expectedImage->GetPixel(index) <= 1) { CHECK_EQUAL(convertedImagePtr->GetPixel(index), expectedImage->GetPixel(index)); } } for (unsigned int i = 0; i < 20; i++) { index[0] = sizeX - 1 - i; index[1] = sizeY - 1 - i; index[2] = sizeZ - 1 - i; if (expectedImage->GetPixel(index) >= 0 && expectedImage->GetPixel(index) <= 1) { CHECK_EQUAL(convertedImagePtr->GetPixel(index), expectedImage->GetPixel(index)); } } for (unsigned int i = 0; i < 20 && (sizeX / 2 - i) < sizeX && (sizeY / 2 - i) < sizeY && (sizeZ / 2 - i) < sizeZ; i++) { index[0] = sizeX / 2 - i; index[1] = sizeY / 2 - i; index[2] = sizeZ / 2 - i; if (expectedImage->GetPixel(index) >= 0 && expectedImage->GetPixel(index) <= 1) { CHECK_EQUAL(convertedImagePtr->GetPixel(index), expectedImage->GetPixel(index)); } } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/models/BioModelScatterPlotTest.cpp b/testing/models/BioModelScatterPlotTest.cpp index d5ade24..6e068d8 100644 --- a/testing/models/BioModelScatterPlotTest.cpp +++ b/testing/models/BioModelScatterPlotTest.cpp @@ -1,243 +1,239 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbBaseTypeModels.h" #include "rttbBioModel.h" #include "rttbBioModelScatterPlots.h" #include "../core/DummyDVHGenerator.h" #include "DummyModel.h" #include "rttbIntegration.h" #include #include #include namespace rttb { namespace testing { typedef models::ScatterPlotType ScatterPlotType; /*! @brief BioModelScatterPlotTest. Test the scatterplot methods. 1) test if parameters are set correctly 2) test if scatterData contain each model value exactly once */ int BioModelScatterPlotTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; //generate artificial DVH and corresponding statistical values - DoseTypeGy binSize = DoseTypeGy(0.1); - DoseVoxelVolumeType voxelVolume = 8; - - DoseCalcType value = 1000; DummyDVHGenerator dummyDVH; const IDType structureID = "myStructure"; const IDType doseID = "myDose"; core::DVH::DVHPointer dvhPtr = boost::make_shared(dummyDVH.generateDVH(structureID, doseID, 0, 2000)); //test Dummy Model models::BioModelParamType param1 = 0.35; models::BioModelParamType param2 = 0.023333333333333; models::BioModelParamType param3 = 10000000; //model values will be dose/normalisationDose DoseTypeGy normalisationDose = 10; rttb::models::DummyModel model(dvhPtr); model.resetCounters(); //default number of points is 100 CHECK_NO_THROW(models::getScatterPlotVary1Parameter(model, 0, param1, 0.5, normalisationDose)); model.resetCounters(); ScatterPlotType scatter = models::getScatterPlotVary1Parameter(model, 0, param1, 0.5, normalisationDose); //only 1st parameter should gave been changed CHECK_EQUAL(100, model.getSetParam1Count()); CHECK_EQUAL(0, model.getSetParam2Count()); CHECK_EQUAL(0, model.getSetParam3Count()); CHECK_EQUAL(100, model.getCalcCount()); // all model values should match the corresponding dose/normalizationDose ScatterPlotType::iterator iter; for (iter = scatter.begin(); iter != scatter.end(); ++iter) { double currentKey = iter->first; if ((currentKey / normalisationDose) != iter->second.first) { CHECK_EQUAL((currentKey / normalisationDose), iter->second.first); } } model.resetCounters(); //number of points is 50 CHECK_NO_THROW(models::getScatterPlotVary1Parameter(model, 1, param2, 0.25, normalisationDose, 50)); model.resetCounters(); scatter = models::getScatterPlotVary1Parameter(model, 1, param2, 0.25, normalisationDose, 50); //only 1st parameter should gave been changed CHECK_EQUAL(0, model.getSetParam1Count()); CHECK_EQUAL(50, model.getSetParam2Count()); CHECK_EQUAL(0, model.getSetParam3Count()); CHECK_EQUAL(50, model.getCalcCount()); // all model values should match the corresponding dose/normalizationDose for (iter = scatter.begin(); iter != scatter.end(); ++iter) { double currentKey = iter->first; if ((currentKey / normalisationDose) != iter->second.first) { CHECK_EQUAL((currentKey / normalisationDose), iter->second.first); } } model.resetCounters(); //number of points is 150 CHECK_NO_THROW(models::getScatterPlotVary1Parameter(model, 2, param3, 0.75, normalisationDose, 150)); model.resetCounters(); scatter = models::getScatterPlotVary1Parameter(model, 2, param3, 0.75, normalisationDose, 150); //only 1st parameter should gave been changed CHECK_EQUAL(0, model.getSetParam1Count()); CHECK_EQUAL(0, model.getSetParam2Count()); CHECK_EQUAL(150, model.getSetParam3Count()); CHECK_EQUAL(150, model.getCalcCount()); // all model values should match the corresponding dose/normalizationDose for (iter = scatter.begin(); iter != scatter.end(); ++iter) { double currentKey = iter->first; if ((currentKey / normalisationDose) != iter->second.first) { CHECK_EQUAL((currentKey / normalisationDose), iter->second.first); } } model.resetCounters(); std::vector paramIDVec; models::BioModel::ParamVectorType meanVec; models::BioModel::ParamVectorType varianceVec; paramIDVec.push_back(0); meanVec.push_back(1); varianceVec.push_back(0.6); paramIDVec.push_back(1); meanVec.push_back(10); varianceVec.push_back(0.5); paramIDVec.push_back(2); meanVec.push_back(100); varianceVec.push_back(0.4); //number of points is 50 CHECK_NO_THROW(models::getScatterPlotVaryParameters(model, paramIDVec, meanVec, varianceVec, normalisationDose, 50)); model.resetCounters(); scatter = models::getScatterPlotVaryParameters(model, paramIDVec, meanVec, varianceVec, normalisationDose, 50); //only 1st parameter should gave been changed CHECK_EQUAL(50, model.getSetParam1Count()); CHECK_EQUAL(50, model.getSetParam2Count()); CHECK_EQUAL(50, model.getSetParam3Count()); CHECK_EQUAL(50, model.getCalcCount()); // all model values should match the corresponding dose/normalizationDose for (iter = scatter.begin(); iter != scatter.end(); ++iter) { double currentKey = iter->first; if ((currentKey / normalisationDose) != iter->second.first) { CHECK_EQUAL((currentKey / normalisationDose), iter->second.first); } } model.resetCounters(); //test 2000 points CHECK_NO_THROW(models::getScatterPlotVaryParameters(model, paramIDVec, meanVec, varianceVec, normalisationDose, 10000)); //test iterativeIntegration /*std::string bedFileName="d:/Temp/bed.txt"; std::ifstream dvh_ifstr(bedFileName,std::ios::in); std::vector volumeV2; std::vector bedV2; //std::cout << "iterative integration"<< std::endl; if ( !dvh_ifstr.is_open() ) { std::cerr<< "DVH file name invalid: could not open the dvh file!"<> volume; volumeV2.push_back(volume); //std::cout << "("<< volume << ","; std::string lineSub2 = line.substr(found2+1,found3-found2-1); //std::cout << lineSub2 << std::endl; std::stringstream ss2(lineSub2); double bed; ss2 >> bed; bedV2.push_back(bed); //std::cout << bed << ");"; } } } struct rttb::models::TcpParams params={0.302101,0.08,10000000,volumeV2,bedV2}; double result = rttb::models::integrateTCP(0, params);*/ RETURN_AND_REPORT_TEST_SUCCESS; } }//testing -}//rttb \ No newline at end of file +}//rttb