diff --git a/code/core/rttbBaseType.h b/code/core/rttbBaseType.h index 618efa8..041bd02 100644 --- a/code/core/rttbBaseType.h +++ b/code/core/rttbBaseType.h @@ -1,719 +1,719 @@ // ----------------------------------------------------------------------- // 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. @todo Both UnsignedIndex2D and VoxelGridIndex2D required? */ 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. @todo Both UnsignedIndex3D and VoxelGridIndex3D required? */ 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++) { 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++) { 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++) { if (wc1(i)!=wc2(i)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const WorldCoordinate3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; 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. @todo should SpacingVectorType/GridVolumeType be forced to be non-negative? If yes add corresponding tests */ 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++) { if (wc1(i)!=wc2(i)) { return false; } } return true; } 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 @todo should OrientationMatrix be forced to be non-negative? If yes add corresponding tests */ 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 equal(const OrientationMatrix& anOrientationMatrix) 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 (!((*this)(m,n) == anOrientationMatrix(m,n))) { return false; } } }// end element comparison } else { return false; } } else { return false; } return true; } friend bool operator==(const OrientationMatrix& om1, const OrientationMatrix& om2) { return om1.equal(om2); } 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++) { 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++) { 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, DoseVoxelVolumeType, VolumeType, GridVolumeType; + typedef double DoseCalcType, DoseTypeGy, 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 a12b7c3..f43adb8 100644 --- a/code/core/rttbDVH.cpp +++ b/code/core/rttbDVH.cpp @@ -1,372 +1,371 @@ // ----------------------------------------------------------------------- // 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 ©){ _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(); } DVH &DVH::operator=(const DVH ©) { 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() << ", " < 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; DataDifferentialType::iterator it; size_t size=_dataDifferential.size(); for(int i=0;igetNumberOfVoxels()); } } if(!relativeVolume) { return _dataCumulative; } else { return _dataCumulativeRelative; } } DoseStatisticType DVH::getMedian() const{ double median_voxel=0; int median_i=0; for(GridIndexType i=0;i_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_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 +}//end namespace rttb diff --git a/code/core/rttbDVH.h b/code/core/rttbDVH.h index 214e50b..ae74a40 100644 --- a/code/core/rttbDVH.h +++ b/code/core/rttbDVH.h @@ -1,193 +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) */ #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; 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 ©); /*! @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ DVH &operator=(const DVH ©); /*! 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; /*! @todo Should it be possible to set the DoseID externally? Should this be only possible in the DVHCalculator and corresponding IO classes? */ void setDoseID(IDType aDoseID); /*! @todo Should it be possible to set the StructureID externally? */ 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/models/rttbDvhBasedModels.cpp b/code/models/rttbDvhBasedModels.cpp index 57cf5b5..932d2df 100644 --- a/code/models/rttbDvhBasedModels.cpp +++ b/code/models/rttbDvhBasedModels.cpp @@ -1,153 +1,153 @@ // ----------------------------------------------------------------------- // 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 "rttbDvhBasedModels.h" #include "rttbInvalidParameterException.h" namespace rttb{ - namespace models{ - DoseStatisticType getEUD(const DVHPointer dvh, const DoseCalcType aA) { - if(aA==0) - { - throw core::InvalidParameterException("Parameter invalid: aA should not be zero!"); - } - DataDifferentialType dataDifferential = dvh->getDataDifferential(); - if(dataDifferential.empty()) - { - throw core::InvalidParameterException("Parameter invalid: DVH data differential should not be empty!"); - } - - double eud=0; - - DoseTypeGy deltaD = dvh->getDeltaD(); - DVHVoxelNumber numberOfVoxels = dvh->getNumberOfVoxels(); - for(GridIndexType i=0; i calcBEDDVH(const DVHPointer dvh, const int numberOfFractions, const DoseCalcType alpha_beta, const bool relativeVolume){ - std::map dataBED; - std::map dataBEDRelative; - - DataDifferentialType dataDifferential = dvh->getDataDifferential(); - if(dataDifferential.empty()) - { - throw core::InvalidParameterException("Parameter invalid: DVH data differential should not be empty!"); - } - if(alpha_beta<=0) - { - throw core::InvalidParameterException("Parameter invalid: alpha_beta must be >0!"); - } - if(numberOfFractions<=0) - { - throw core::InvalidParameterException("Parameter invalid: numberOfFractions must be >0!"); - } - - DoseTypeGy deltaD = dvh->getDeltaD(); - DVHVoxelNumber numberOfVoxels = dvh->getNumberOfVoxels(); - - std::deque::iterator it; - int i=0; - - for(it=dataDifferential.begin();it!=dataDifferential.end();it++ ) - { - DoseTypeGy doseGyi=((i+0.5)*deltaD); - DoseTypeGy bedi=0; - bedi=(doseGyi*(1+doseGyi/(numberOfFractions*alpha_beta))); - if(!relativeVolume) - { - dataBED.insert(std::pair(bedi,(*it))); - } - else - { - dataBEDRelative.insert(std::pair(bedi,(*it)/numberOfVoxels)); - } - i++; - } - - if(!relativeVolume) - { - return dataBED; - } - else - { - return dataBEDRelative; - } - } - - std::map calcLQED2DVH(const DVHPointer dvh, const int numberOfFractions, const DoseCalcType alpha_beta, const bool relativeVolume) { - std::map dataLQED2; - std::map dataLQED2Relative; - - DataDifferentialType dataDifferential = dvh->getDataDifferential(); - if(dataDifferential.empty()) - { - throw core::InvalidParameterException("Parameter invalid: DVH data differential should not be empty!"); - } - if(alpha_beta<=0) - { - throw core::InvalidParameterException("Parameter invalid: alpha_beta must be >0!"); - } - if(numberOfFractions<=0) - { - throw core::InvalidParameterException("Parameter invalid: numberOfFractions must be >0!"); - } - - DoseTypeGy deltaD = dvh->getDeltaD(); - DVHVoxelNumber numberOfVoxels = dvh->getNumberOfVoxels(); - - std::deque::iterator it; - int i=0; - - for(it=dataDifferential.begin();it!=dataDifferential.end();it++ ) - { - DoseTypeGy doseGyi=((i+0.5)*deltaD); - DoseTypeGy lqed2i=0; - lqed2i=(doseGyi*((alpha_beta+doseGyi/numberOfFractions)/(alpha_beta+2))); - - if(!relativeVolume) - { - dataLQED2.insert(std::pair(lqed2i,*it)); - } - else - { - dataLQED2Relative.insert(std::pair(lqed2i,(*it)/numberOfVoxels)); - } - i++; - } - - if(!relativeVolume) - { - return dataLQED2; - } - else - { - return dataLQED2Relative; - } - } - } + namespace models{ + DoseStatisticType getEUD(const DVHPointer dvh, const DoseCalcType aA) { + if(aA==0) + { + throw core::InvalidParameterException("Parameter invalid: aA should not be zero!"); + } + DataDifferentialType dataDifferential = dvh->getDataDifferential(); + if(dataDifferential.empty()) + { + throw core::InvalidParameterException("Parameter invalid: DVH data differential should not be empty!"); + } + + double eud=0; + + DoseTypeGy deltaD = dvh->getDeltaD(); + DVHVoxelNumber numberOfVoxels = dvh->getNumberOfVoxels(); + for(GridIndexType i=0; i dataBED; + std::map dataBEDRelative; + + DataDifferentialType dataDifferential = dvh->getDataDifferential(); + if(dataDifferential.empty()) + { + throw core::InvalidParameterException("Parameter invalid: DVH data differential should not be empty!"); + } + if(alpha_beta<=0) + { + throw core::InvalidParameterException("Parameter invalid: alpha_beta must be >0!"); + } + if(numberOfFractions<=0) + { + throw core::InvalidParameterException("Parameter invalid: numberOfFractions must be >1! The dvh should be an accumulated-dvh of all fractions, not a single fraction-dvh!"); + } + + DoseTypeGy deltaD = dvh->getDeltaD(); + DVHVoxelNumber numberOfVoxels = dvh->getNumberOfVoxels(); + + std::deque::iterator it; + int i=0; + + for(it=dataDifferential.begin();it!=dataDifferential.end();it++ ) + { + DoseTypeGy doseGyi=((i+0.5)*deltaD); + DoseTypeGy bedi=0; + bedi=(doseGyi*(1+doseGyi/(numberOfFractions*alpha_beta))); + if(!relativeVolume) + { + dataBED.insert(std::pair(bedi,(*it))); + } + else + { + dataBEDRelative.insert(std::pair(bedi,(*it)/numberOfVoxels)); + } + i++; + } + + if(!relativeVolume) + { + return dataBED; + } + else + { + return dataBEDRelative; + } + } + + LQEDDVHType calcLQED2DVH(const DVHPointer dvh, const int numberOfFractions, const DoseCalcType alpha_beta, const bool relativeVolume) { + std::map dataLQED2; + std::map dataLQED2Relative; + + DataDifferentialType dataDifferential = dvh->getDataDifferential(); + if(dataDifferential.empty()) + { + throw core::InvalidParameterException("Parameter invalid: DVH data differential should not be empty!"); + } + if(alpha_beta<=0) + { + throw core::InvalidParameterException("Parameter invalid: alpha_beta must be >0!"); + } + if(numberOfFractions<=1) + { + throw core::InvalidParameterException("Parameter invalid: numberOfFractions must be >1! The dvh should be an accumulated-dvh of all fractions, not a single fraction-dvh!"); + } + + DoseTypeGy deltaD = dvh->getDeltaD(); + DVHVoxelNumber numberOfVoxels = dvh->getNumberOfVoxels(); + + std::deque::iterator it; + int i=0; + + for(it=dataDifferential.begin();it!=dataDifferential.end();it++ ) + { + DoseTypeGy doseGyi=((i+0.5)*deltaD); + DoseTypeGy lqed2i=0; + lqed2i=(doseGyi*((alpha_beta+doseGyi/numberOfFractions)/(alpha_beta+2))); + + if(!relativeVolume) + { + dataLQED2.insert(std::pair(lqed2i,*it)); + } + else + { + dataLQED2Relative.insert(std::pair(lqed2i,(*it)/numberOfVoxels)); + } + i++; + } + + if(!relativeVolume) + { + return dataLQED2; + } + else + { + return dataLQED2Relative; + } + } + } } \ No newline at end of file diff --git a/code/models/rttbDvhBasedModels.h b/code/models/rttbDvhBasedModels.h index 856bb7d..c10b073 100644 --- a/code/models/rttbDvhBasedModels.h +++ b/code/models/rttbDvhBasedModels.h @@ -1,45 +1,50 @@ #include #include "rttbDVH.h" #include "rttbBaseType.h" namespace rttb{ namespace models{ typedef core::DVH::DataDifferentialType DataDifferentialType; typedef core::DVH::DVHPointer DVHPointer; + typedef std::map BEDDVHType; + typedef std::map LQEDDVHType; + typedef double BioModelParameterType; /*! @brief Get Equivalent Uniform Dose (EUD) - @pre dvh data differential is not empty, - @pre aA is not zero, - @return Return calculated EUD value, - @exception InvalidParameterException Thrown if parameters were not set correctly. + @pre dvh data differential is not empty, + @pre aA is not zero, + @return Return calculated EUD value, + @exception InvalidParameterException Thrown if parameters were not set correctly. */ DoseStatisticType getEUD(const DVHPointer dvh, const DoseCalcType aA); /*! @brief Calculate Biological Effective/Equivalent Dose (BED) of dvh - @param relativeVolume default false-> the corresponding bed value is the voxel number of the dose bin; - if true-> the corresponding bed value is the relative volume % between 0 and 1, - (the voxel number of this dose bin)/(number of voxels) - @pre dvh data differential is not empty - @pre alpha_beta > 0 - @pre numberOfFractions > 0 - @return map of BED values - @exception InvalidParameterException Thrown if parameters were not set correctly. + @param relativeVolume default false-> the corresponding volume value is the voxel number of the dose bin; + if true-> the corresponding volume value is the relative volume % between 0 and 1, + (the voxel number of this dose bin)/(number of voxels) + @pre dvh should be an accumulated dvh of all fractions, not a single fraction dvh + @pre dvh data differential is not empty + @pre alpha_beta > 0 + @pre numberOfFractions > 1 + @return Return map: keys are BEDi in Gy, values are the volume of the dose bin + @exception InvalidParameterException Thrown if parameters were not set correctly. */ - std::map calcBEDDVH(const DVHPointer dvh, const int numberOfFractions, + BEDDVHType calcBEDDVH(const DVHPointer dvh, const int numberOfFractions, const DoseCalcType alpha_beta, const bool relativeVolume=false); /*! @brief Calculate Linear-quadratic equivalent dose for 2-Gy (LQED2) of dvh - @param relativeVolume default false-> the corresponding LQED2 value is the voxel number of the dose bin; - if true-> the corresponding LQED2 value is the relative volume % between 0 and 1, - (the voxel number of this dose bin)/(number of voxels) - @pre dvh data differential is not empty - @pre alpha_beta > 0 - @pre numberOfFractions > 0 - @return map of LQED2 values - @exception InvalidParameterException Thrown if parameters were not set correctly. + @param relativeVolume default false-> the corresponding volume value is the voxel number of the dose bin; + if true-> the corresponding volume value is the relative volume % between 0 and 1, + (the voxel number of this dose bin)/(number of voxels) + @pre dvh should be an accumulated dvh of all fractions, not a single fraction dvh + @pre dvh data differential is not empty + @pre alpha_beta > 0 + @pre numberOfFractions > 1 + @return Return map: keys are LQED2 in Gy, values are the volume of the dose bin; return empty map if not initialized + @exception InvalidParameterException Thrown if parameters were not set correctly. */ - std::map calcLQED2DVH(const DVHPointer dvh, const int numberOfFractions, + LQEDDVHType calcLQED2DVH(const DVHPointer dvh, const int numberOfFractions, const DoseCalcType alpha_beta, const bool relativeVolume=false); - } - } \ No newline at end of file + } +} \ No newline at end of file diff --git a/code/models/rttbTCPLQModel.cpp b/code/models/rttbTCPLQModel.cpp index 50f99da..9bddc66 100644 --- a/code/models/rttbTCPLQModel.cpp +++ b/code/models/rttbTCPLQModel.cpp @@ -1,229 +1,235 @@ // ----------------------------------------------------------------------- // 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 #define _USE_MATH_DEFINES #include #include #include #include #include "rttbTCPLQModel.h" #include "rttbDvhBasedModels.h" #include "rttbBaseTypeModels.h" #include "rttbIntegration.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb{ namespace models{ TCPLQModel::TCPLQModel(): TCPModel(),_alphaMean(0), _alphaVariance(0), _alpha_beta(0), _rho(0){} TCPLQModel::TCPLQModel(DVHPointer aDVH, BioModelParamType aAlphaMean, BioModelParamType aBeta, BioModelParamType aRho, int aNumberOfFractions): TCPModel(aDVH,aNumberOfFractions),_alphaMean(aAlphaMean), _alphaVariance(0), _alpha_beta(aAlphaMean/aBeta), _rho(aRho){} TCPLQModel::TCPLQModel(DVHPointer aDVH, BioModelParamType aRho, int aNumberOfFractions, BioModelParamType aAlpha_Beta, BioModelParamType aAlphaMean, BioModelParamType aAlphaVariance): TCPModel(aDVH,aNumberOfFractions),_alphaMean(aAlphaMean), _alphaVariance(aAlphaVariance), _alpha_beta(aAlpha_Beta), _rho(aRho){} void TCPLQModel::setParameters(const BioModelParamType aAlphaMean, const BioModelParamType aAlpha_Beta, const BioModelParamType aRho, const BioModelParamType aAlphaVariance){ _alphaMean=aAlphaMean; _alphaVariance=aAlphaVariance; _alpha_beta=aAlpha_Beta; _rho=aRho; //reset _value, because parameters have changed. _value=0; } void TCPLQModel::setAlpha(const BioModelParamType aAlphaMean, const BioModelParamType aAlphaVariance){ _alphaVariance=aAlphaVariance; _alphaMean=aAlphaMean; } void TCPLQModel::setAlphaBeta(const BioModelParamType aAlpha_Beta){ _alpha_beta=aAlpha_Beta; } void TCPLQModel::setRho(const BioModelParamType aRho){_rho=aRho;} const BioModelParamType TCPLQModel::getAlpahBeta(){return _alpha_beta;} const BioModelParamType TCPLQModel::getAlphaMean(){return _alphaMean;} const BioModelParamType TCPLQModel::getAlphaVariance(){return _alphaVariance;} const BioModelParamType TCPLQModel::getRho(){return _rho;} long double TCPLQModel::calcTCPi(BioModelParamType aRho, BioModelParamType aAlphaMean, double vj, double bedj){ return exp(-aRho*vj*exp(-aAlphaMean*bedj)); } long double TCPLQModel::calcTCP(std::map aBEDDVH, BioModelParamType aRho, BioModelParamType aAlphaMean, double aDeltaV){ std::map::iterator it; long double tcp=1; for(it=aBEDDVH.begin(); it!=aBEDDVH.end();it++) { long double tcpi=this->calcTCPi(aRho,aAlphaMean,(*it).second*aDeltaV,(*it).first); tcp=tcp*tcpi; } return tcp; } long double TCPLQModel::calcTCPAlphaNormalDistribution(std::map aBEDDVH, BioModelParamType aRho, BioModelParamType aAlphaMean, BioModelParamType aAlphaVariance, double aDeltaV){ std::map::iterator it; int bedSize=aBEDDVH.size(); std::vector volumeV2; std::vector bedV2; int i=0; for(it=aBEDDVH.begin();it!=aBEDDVH.end();it++) { volumeV2.push_back((*it).second*aDeltaV); bedV2.push_back((*it).first); i++; } struct TcpParams params={aAlphaMean,aAlphaVariance,aRho,volumeV2,bedV2}; double result = integrateTCP(0, params); if(result == -100) { std::cerr << "Integration error!\n"; return -1; } else{ long double tcp=1/(pow(2*M_PI, 0.5)*_alphaVariance)*result; return tcp; } } BioModelValueType TCPLQModel::calcModel(const double doseFactor){ core::DVH variantDVH=core::DVH(_dvh->getDataDifferential(),(DoseTypeGy)(_dvh->getDeltaD()*doseFactor), _dvh->getDeltaV(),"temporary","temporary"); boost::shared_ptr spDVH = boost::make_shared(variantDVH); BioModelValueType value=0; if(_alphaVariance==0){ - if(_alphaMean<=0 && _alpha_beta<=0 && _rho<=0 && _numberOfFractions<=0) + if(_alphaMean<=0 && _alpha_beta<=0 && _rho<=0 ) { throw core::InvalidParameterException("Parameter invalid: alpha, alpha/beta, rho and number of fractions must >0!"); } + if(_numberOfFractions<=1){ + throw core::InvalidParameterException("Parameter invalid: numberOfFractions must be >1! The dvh should be an accumulated-dvh of all fractions, not a single fraction-dvh!"); + } long double tcp=1; std::map dataBED=calcBEDDVH(spDVH,_numberOfFractions, _alpha_beta); value=(BioModelValueType)this->calcTCP(dataBED,_rho,_alphaMean,variantDVH.getDeltaV()); return value; } //if alpha normal distribution else { if(this->_alpha_beta<=0 && this->_alphaMean<=0 && this->_alphaVariance<0 && _rho<=0 && _numberOfFractions<=0) { throw core::InvalidParameterException("Parameter invalid: alpha/beta, alphaMean, rho and number of fractions must >0!"); } + if(_numberOfFractions<=1){ + throw core::InvalidParameterException("Parameter invalid: numberOfFractions must be >1! The dvh should be an accumulated-dvh of all fractions, not a single fraction-dvh!"); + } std::map dataBED=calcBEDDVH(spDVH,_numberOfFractions, _alpha_beta); value=(BioModelValueType)(this->calcTCPAlphaNormalDistribution(dataBED,_rho,_alphaMean,_alphaVariance, variantDVH.getDeltaV())); return value; } } void TCPLQModel::setParameterVector(const ParamVectorType aParameterVector){ if(aParameterVector.size()!=4) { throw core::InvalidParameterException("Parameter invalid: aParameterVector.size must be 4! "); } else { _alphaMean=aParameterVector.at(0); _alphaVariance=aParameterVector.at(1); _alpha_beta=aParameterVector.at(2); _rho=aParameterVector.at(3); } } void TCPLQModel::setParameterByID(const int aParamId, const BioModelParamType aValue){ if(aParamId==0) { _alphaMean=aValue; } else if(aParamId==1) { _alphaVariance=aValue; } else if(aParamId==2) { _alpha_beta=aValue; } else if(aParamId==3) { _rho=aValue; } else { throw core::InvalidParameterException("Parameter invalid: aParamID must be 0(alphaMean) or 1(alphaVariance) or 2(alpha_beta) or 3(rho)! "); } } const int TCPLQModel::getParameterID(const std::string aParamName) const{ if(aParamName=="alphaMean") { return 0; } else if(aParamName=="alphaVariance") { return 1; } else if(aParamName=="alpha_beta") { return 2; } else if(aParamName=="rho") { return 3; } else { rttbExceptionMacro(core::InvalidParameterException, << "Parameter name "< #include #include #include "rttbBaseType.h" #include "rttbTCPModel.h" #include "rttbBaseTypeModels.h" namespace rttb{ namespace models{ /*! @class TCPLQModel @brief This class represents a TCP(Tumor Control Probability) LQ model (Nahum and Sanchez-Nieto 2001, Hall and Giaccia 2006) @see TCPModel */ class TCPLQModel: public TCPModel { public: typedef TCPModel::ParamVectorType ParamVectorType; typedef TCPModel::DVHPointer DVHPointer; private: /*! @brief Calculate intermediate tcp using alpha constant. This is a helper function for calcTCP() @see calcTCP */ long double calcTCPi(BioModelParamType aRho, BioModelParamType aAlphaMean, double vj, double bedj); /*! @brief Calculate tcp using alpha constant. */ long double calcTCP(std::map aBEDDVH, BioModelParamType aRho, BioModelParamType aAlphaMean, double aDeltaV); /*! @brief Calculate tcp using a normal distribution for alpha. */ long double calcTCPAlphaNormalDistribution(std::map aBEDDVH, BioModelParamType aRho, BioModelParamType aAlphaMean, BioModelParamType aAlphaVariance, double aDeltaV); protected: BioModelParamType _alphaMean; BioModelParamType _alphaVariance; BioModelParamType _alpha_beta; /*! Roh is the initial clonogenic cell density */ BioModelParamType _rho; /*! @brief Calculate the model value @param doseFactor scaling factor for prescribed dose. The model calculation will use the dvh with each di=old di*doseFactor. @pre _alphaMean >0 @pre _alphaVariance >= 0 @pre _alpha_beta > 0 @pre _rho > 0 + @pre _numberOfFractions > 1 @exception InvalidParameterException Thrown if parameters were not set correctly. */ BioModelValueType calcModel(const double doseFactor=1); public: TCPLQModel(); /*! @brief Constructor initializes member variables with given parameters. + @pre aAlphaMean >0 + @pre aBeta > 0 + @pre aRho > 0 + @pre aNumberOfFractions > 1 */ TCPLQModel(DVHPointer aDVH, BioModelParamType aAlphaMean, BioModelParamType aBeta, BioModelParamType aRho, int aNumberOfFractions); /*! @brief Constructor for alpha distribution initializes member variables with given parameters. + @pre aAlphaMean >0 + @pre aAlphaVariance >0 + @pre aAlpha_Beta > 0 + @pre aRho > 0 + @pre aNumberOfFractions > 1 */ TCPLQModel(DVHPointer aDVH, BioModelParamType aRho, int aNumberOfFractions, BioModelParamType aAlpha_Beta, BioModelParamType aAlphaMean, BioModelParamType aAlphaVariance); const BioModelParamType getRho(); void setRho(const BioModelParamType aRho); - + const BioModelParamType getAlphaMean(); const BioModelParamType getAlphaVariance(); /*! @brief The distribution of the parameter alpha, which is characteristic for a population of cells, is described by the its mean and variance. If alpha is constant the variance is 0. @param aAlphaVariance The variance of alpha can be given, the default value is 0 resulting in constant alpha. */ void setAlpha(const BioModelParamType aAlphaMean, const BioModelParamType aAlphaVariance=0); const BioModelParamType getAlpahBeta(); void setAlphaBeta(const BioModelParamType aAlpha_Beta); /*! @brief Set parameters for the TCP model. _value will be reset to 0. @param aAlpha_Beta alpha/beta constant . @param aAlphaMean mean of alpha distribution. @param aAlphaVariance variance of alpha distribution. */ void setParameters(const BioModelParamType aAlphaMean, const BioModelParamType aAlpha_Beta, const BioModelParamType aRho, const BioModelParamType aAlphaVariance=0); /*! @brief Set parameter with ID. "alphaMean":0,"alphaVariance":1,"alpha_beta":2, "rho":3 @exception InvalidParameterException Thrown if aParamId is not 0 or 1 or 2 or 3. */ virtual void setParameterByID(const int aParamId, const BioModelParamType aValue); /*! @brief Set parameter vector, where index of vector is the parameter id. "alphaMean":0,"alphaVariance":1,"alpha_beta":2, "rho":3 @exception InvalidParameterException Thrown if aParamterVector.size()!=4. */ virtual void setParameterVector(const ParamVectorType aParameterVector); /*! @brief Get parameter id. "alphaMean":0,"alphaVariance":1,"alpha_beta":2, "rho":3 @return 0 for "alphaMean", 1 for "alphaVariance", 2 for "alpha_beta", 3 for "rho" @exception InvalidParameterException Thrown if aParamName is not alphaMean or alphaVariance or alpha_beta or rho. */ virtual const int getParameterID(const std::string aParamName) const; }; }//end algorithms - }//end rttb +}//end rttb #endif diff --git a/code/models/rttbTCPModel.h b/code/models/rttbTCPModel.h index 2cafa04..8eef010 100644 --- a/code/models/rttbTCPModel.h +++ b/code/models/rttbTCPModel.h @@ -1,60 +1,63 @@ // ----------------------------------------------------------------------- // 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 __TCP_MODEL_H #define __TCP_MODEL_H #include "rttbBaseType.h" #include "rttbBioModel.h" #include "rttbBaseTypeModels.h" namespace rttb{ namespace models{ /*! @class TCPModel @brief This is the interface class for TCP(Tumor Control Probability) models */ class TCPModel: public BioModel { public: typedef BioModel::ParamVectorType ParamVectorType; typedef BioModel::DVHPointer DVHPointer; protected: int _numberOfFractions; public: + TCPModel(): BioModel(), _numberOfFractions(0){}; + TCPModel(int aNum): BioModel(), _numberOfFractions(aNum){}; + TCPModel(DVHPointer aDvh, int aNum): BioModel(aDvh), _numberOfFractions(aNum){}; void setNumberOfFractions(const int aNumberOfFractions); const int getNumberOfFractions(); }; }//end namespace models }//end namespace rttb #endif diff --git a/testing/core/DVHTest.cpp b/testing/core/DVHTest.cpp index 0ed832f..78b0712 100644 --- a/testing/core/DVHTest.cpp +++ b/testing/core/DVHTest.cpp @@ -1,179 +1,178 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVH.h" #include "DummyDoseAccessor.h" #include "DummyMaskAccessor.h" namespace rttb{ namespace testing{ typedef core::DVH::DataDifferentialType DataDifferentialType; /*! @brief DVHTest - test the API of DVH - 1) test constructors (values as expected?) - 2) test asignement - 3) test set/getLabel - 4) test set/get - 5) test equality + 1) test constructors (values as expected?) + 2) test asignement + 3) test set/getLabel + 4) test set/get + 5) test equality */ int DVHTest(int argc, char* argv[] ) - { + { PREPARE_DEFAULT_TEST_REPORTING; //generate artificial DVH and corresponding statistical values DoseTypeGy binSize = DoseTypeGy(0.1); DoseVoxelVolumeType voxelVolume = 8; DataDifferentialType anEmptyDataDifferential; DataDifferentialType aDataDifferential; DataDifferentialType aDataDifferential2; DataDifferentialType aDataDifferentialRelative; DoseStatisticType maximum = 0; DoseStatisticType minimum = 0; double sum = 0; double squareSum = 0; DoseCalcType value = 0; DVHVoxelNumber numberOfVoxels = 0; // creat default values [0,100) for(int i = 0; i < 100; i++){ value = DoseCalcType((double(rand())/RAND_MAX)*1000); numberOfVoxels+=value; aDataDifferential.push_back( value ); aDataDifferential2.push_back( value*10 ); if (value > 0){ maximum = (i+0.5)*binSize; if (minimum==0){ minimum = (i+0.5)*binSize; - } } + } sum+=value*(i+0.5)*binSize; squareSum+=value*pow((i+0.5)*binSize,2); - } + } DoseStatisticType mean = sum/numberOfVoxels; DoseStatisticType variance=(squareSum/numberOfVoxels-mean*mean); DoseStatisticType stdDeviation=pow(variance,0.5); std::deque::iterator it; for(it=aDataDifferential.begin();it!=aDataDifferential.end();it++) - { + { aDataDifferentialRelative.push_back((*it)/numberOfVoxels); - } + } const IDType structureID = "myStructure"; const IDType doseID = "myDose"; const IDType voxelizationID = "myVoxelization"; //1) test default constructor (values as expected?) CHECK_THROW(core::DVH(anEmptyDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(anEmptyDataDifferential, binSize, voxelVolume, structureID, doseID)); - CHECK_NO_THROW(core::DVH(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID)); CHECK_NO_THROW(core::DVH(aDataDifferential, binSize, voxelVolume, structureID, doseID)); CHECK_THROW(core::DVH(aDataDifferential, 0, voxelVolume, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(aDataDifferential, 0, voxelVolume, structureID, doseID)); CHECK_THROW(core::DVH(aDataDifferential, binSize, 0, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(aDataDifferential, binSize, 0, structureID, doseID)); //2) test asignement core::DVH myTestDVH(aDataDifferential, binSize, voxelVolume, structureID, doseID); CHECK_NO_THROW(core::DVH myOtherDVH = myTestDVH); const core::DVH myOtherDVH = myTestDVH; CHECK_NO_THROW(core::DVH aDVH(myOtherDVH)); //3) test set/getLabel core::DVH myDVH(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID); StructureLabel label = ""; CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); label = "myLabel"; CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); label = "myLabel2"; CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); //4) test set/get //IDs CHECK_EQUAL(myDVH.getStructureID(), structureID); CHECK_EQUAL(myDVH.getDoseID(), doseID); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); /*! @todo Should you be able to set the DoseID?*/ myDVH.setDoseID("somethingElse"); CHECK_EQUAL(myDVH.getDoseID(), "somethingElse"); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); CHECK_EQUAL(myDVH.getStructureID(), structureID); /*! @todo Should you be able to set the DoseID?*/ myDVH.setStructureID("somethingOther"); CHECK_EQUAL(myDVH.getDoseID(), "somethingElse"); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); CHECK_EQUAL(myDVH.getStructureID(), "somethingOther"); //dataDifferential CHECK(myDVH.getDataDifferential()==aDataDifferential); CHECK(myDVH.getDataDifferential(false) == aDataDifferential); CHECK(myDVH.getDataDifferential(true)==aDataDifferentialRelative); - + CHECK_EQUAL(myDVH.getNumberOfVoxels(), numberOfVoxels); CHECK_EQUAL(myDVH.getDeltaV(), voxelVolume); CHECK_EQUAL(myDVH.getDeltaD(), binSize); CHECK_EQUAL(myDVH.getMaximum(), maximum); CHECK_EQUAL(myDVH.getMinimum(), minimum); CHECK_EQUAL(myDVH.getMean(), mean); CHECK_EQUAL(myDVH.getVariance(), variance); CHECK_EQUAL(myDVH.getStdDeviation(), stdDeviation); int percentage = 20; VolumeType absVol = VolumeType(percentage*numberOfVoxels*voxelVolume/100.0); CHECK_EQUAL(myDVH.getAbsoluteVolume(percentage), absVol); //5) test equality core::DVH myDVH2(aDataDifferential2, binSize, voxelVolume, structureID, doseID); CHECK(!(myDVH==myDVH2)); CHECK_EQUAL(myDVH,myDVH); core::DVH aDVH(myOtherDVH); CHECK_EQUAL(aDVH,myOtherDVH); RETURN_AND_REPORT_TEST_SUCCESS; - } + } - }//end namespace testing - }//end namespace rttb + }//end namespace testing +}//end namespace rttb diff --git a/testing/examples/RTBioModelExampleTest.cpp b/testing/examples/RTBioModelExampleTest.cpp index bd2180b..38acdf3 100644 --- a/testing/examples/RTBioModelExampleTest.cpp +++ b/testing/examples/RTBioModelExampleTest.cpp @@ -1,415 +1,374 @@ // ----------------------------------------------------------------------- // 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 "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 - { +namespace rttb{ + namespace testing{ - /*! @brief RTBioModelTest. - TCP calculated using a DVH PTV and LQ Model. + /*! @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 + 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[]) + int RTBioModelExampleTest(int argc, char* argv[] ) { PREPARE_DEFAULT_TEST_REPORTING; typedef rttb::models::CurveDataType CurveDataType; - typedef std::multimap > ScatterPlotType; + 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) + if (argc>1) { DVH_FILENAME_PTV = argv[1]; } - - if (argc > 2) + if (argc>2) { DVH_FILENAME_NT1 = argv[2]; } - - if (argc > 3) + if (argc>3) { DVH_FILENAME_NT2 = argv[3]; } - - if (argc > 4) + if (argc>4) { DVH_FILENAME_NT3 = argv[4]; } - - if (argc > 5) + if (argc>5) { DVH_FILENAME_TV_TEST = argv[5]; } - - if (argc > 6) + if(argc>6) { - DVH_Virtuos_Lung = argv[6]; + DVH_Virtuos_Lung=argv[6]; } - - if (argc > 7) + if(argc>7) { - DVH_Virtuos_Target = argv[7]; + DVH_Virtuos_Target=argv[7]; } - //DVH PTV - rttb::io::other::DVHTxtFileReader dvhReader = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_PTV); - DVHPointer dvhPtr = dvhReader.generateDVH(); + //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); + 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(); + 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 = 1; + 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.getAlpahBeta()); - CHECK_EQUAL(roh, tcplq.getRho()); + rttb::models::TCPLQModel tcplq=rttb::models::TCPLQModel(dvhPtr,alpha, beta, roh, numFractions); + CHECK_EQUAL(alpha,tcplq.getAlphaMean()); + CHECK_EQUAL(alpha/beta,tcplq.getAlpahBeta()); + CHECK_EQUAL(roh,tcplq.getRho()); CHECK_NO_THROW(tcplq.init()); - - if (tcplq.init()) - { - CHECK_CLOSE(1.00497232941856940000e-127, tcplq.getValue(), tolerance); + if(tcplq.init()){ + CHECK_CLOSE(1.00497232941856940000e-127,tcplq.getValue(),tolerance); } - CurveDataType curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); + CurveDataType curve=models::getCurveDoseVSBioModel(tcplq,normalizationDose); CurveDataType::iterator it; - - for (it = curve.begin(); it != curve.end(); ++it) - { - if ((*it).first < 62) - { - CHECK_EQUAL(0, (*it).second); + for(it=curve.begin();it!=curve.end();it++){ + if ((*it).first < 72){ + CHECK_EQUAL(0,(*it).second); } - else if ((*it).first > 114) - { - CHECK((*it).second > 0.9); + 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.getAlpahBeta()); - CHECK_EQUAL(roh, tcplq.getRho()); - - if (tcplq.init()) - { - CHECK_CLOSE(1.58632626255825960000e-002, tcplq.getValue(), tolerance); + tcplq.setParameters(alpha,alphaBeta,roh,0.08); + CHECK_EQUAL(alpha,tcplq.getAlphaMean()); + CHECK_EQUAL(alphaBeta,tcplq.getAlpahBeta()); + CHECK_EQUAL(roh,tcplq.getRho()); + if(tcplq.init()){ + CHECK_CLOSE(1.84e-005,tcplq.getValue(),tolerance); } normalizationDose = 40; - curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); + 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.getAlpahBeta()); - CHECK_EQUAL(roh, tcplq.getRho()); - - if (tcplq.init()) - { - CHECK_CLOSE(9.99338781766346610000e-001, tcplq.getValue(), tolerance); + CHECK_EQUAL(alpha,tcplq.getAlphaMean()); + CHECK_EQUAL(alphaBeta,tcplq.getAlpahBeta()); + 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.getAlpahBeta()); - CHECK_EQUAL(roh, tcplq.getRho()); - - if (tcplq.init()) - { - CHECK_EQUAL(1, tcplq.getValue()); + CHECK_EQUAL(alpha,tcplq.getAlphaMean()); + CHECK_EQUAL(alphaBeta,tcplq.getAlpahBeta()); + 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.getAlpahBeta()); - CHECK_EQUAL(roh, tcplq_test.getRho()); + 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.getAlpahBeta()); + 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); + if(tcplq_test.init()){ + CHECK_CLOSE(9.79050278878883180000e-001,tcplq_test.getValue(),tolerance); } - normalizationDose = 60; - curve = models::getCurveDoseVSBioModel(tcplq_test, normalizationDose); + 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(); + 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); + CHECK_CLOSE(1.07920836034015810000e+001,models::getEUD(dvhPtr2,10),toleranceEUD); //test RTNTCPLKBModel - rttb::models::NTCPLKBModel lkb = rttb::models::NTCPLKBModel(); + 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()); + CHECK_EQUAL(0,lkb.getA()); + CHECK_EQUAL(0,lkb.getM()); + CHECK_EQUAL(0,lkb.getD50()); lkb.setDVH(dvhPtr2); - CHECK_EQUAL(dvhPtr2, lkb.getDVH()); + CHECK_EQUAL(dvhPtr2,lkb.getDVH()); lkb.setA(aVal); - CHECK_EQUAL(aVal, lkb.getA()); + CHECK_EQUAL(aVal,lkb.getA()); lkb.setM(mVal); - CHECK_EQUAL(mVal, lkb.getM()); + CHECK_EQUAL(mVal,lkb.getM()); lkb.setD50(d50Val); - CHECK_EQUAL(d50Val, lkb.getD50()); + CHECK_EQUAL(d50Val,lkb.getD50()); CHECK_NO_THROW(lkb.init()); - - if (lkb.init()) - { - CHECK_CLOSE(2.53523522831366570000e-007, lkb.getValue(), tolerance); + if(lkb.init()){ + CHECK_CLOSE(2.53523522831366570000e-007,lkb.getValue(),tolerance); } //test RTNTCPRSModel - rttb::models::NTCPRSModel rs = rttb::models::NTCPRSModel(); + 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()); + CHECK_EQUAL(0,rs.getGamma()); + CHECK_EQUAL(0,rs.getS()); + CHECK_EQUAL(0,rs.getD50()); rs.setDVH(dvhPtr2); - CHECK_EQUAL(dvhPtr2, rs.getDVH()); + CHECK_EQUAL(dvhPtr2,rs.getDVH()); rs.setD50(d50Val); - CHECK_EQUAL(d50Val, rs.getD50()); + CHECK_EQUAL(d50Val,rs.getD50()); rs.setGamma(gammaVal); - CHECK_EQUAL(gammaVal, rs.getGamma()); + CHECK_EQUAL(gammaVal,rs.getGamma()); rs.setS(sVal); - CHECK_EQUAL(sVal, rs.getS()); + CHECK_EQUAL(sVal,rs.getS()); CHECK_NO_THROW(rs.init()); - - if (rs.init()) - { - CHECK_CLOSE(3.70385888626145740000e-009, rs.getValue(), tolerance); + 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); + 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()); + CHECK_EQUAL(dvhPtr3,lkb.getDVH()); lkb.setA(aVal); - CHECK_EQUAL(aVal, lkb.getA()); + CHECK_EQUAL(aVal,lkb.getA()); lkb.setM(mVal); - CHECK_EQUAL(mVal, lkb.getM()); + CHECK_EQUAL(mVal,lkb.getM()); lkb.setD50(d50Val); - CHECK_EQUAL(d50Val, lkb.getD50()); - - if (lkb.init()) - { - CHECK_CLOSE(7.36294657754956700000e-007, lkb.getValue(), tolerance); + CHECK_EQUAL(d50Val,lkb.getD50()); + if(lkb.init()){ + CHECK_CLOSE(7.36294657754956700000e-007,lkb.getValue(),tolerance); } //test RTNTCPRSModel - rs = rttb::models::NTCPRSModel(); + rs=rttb::models::NTCPRSModel(); gammaVal = 1.7; sVal = 1; - CHECK_EQUAL(0, rs.getGamma()); - CHECK_EQUAL(0, rs.getS()); - CHECK_EQUAL(0, rs.getD50()); + CHECK_EQUAL(0,rs.getGamma()); + CHECK_EQUAL(0,rs.getS()); + CHECK_EQUAL(0,rs.getD50()); rs.setDVH(dvhPtr3); - CHECK_EQUAL(dvhPtr3, rs.getDVH()); + CHECK_EQUAL(dvhPtr3,rs.getDVH()); rs.setD50(d50Val); - CHECK_EQUAL(d50Val, rs.getD50()); + CHECK_EQUAL(d50Val,rs.getD50()); rs.setGamma(gammaVal); - CHECK_EQUAL(gammaVal, rs.getGamma()); + CHECK_EQUAL(gammaVal,rs.getGamma()); rs.setS(sVal); - CHECK_EQUAL(sVal, rs.getS()); - - if (rs.init()) - { - CHECK_CLOSE(1.76778795490939440000e-007, rs.getValue(), tolerance); + 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); + 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()); + CHECK_EQUAL(dvhPtr4,lkb.getDVH()); lkb.setA(aVal); - CHECK_EQUAL(aVal, lkb.getA()); + CHECK_EQUAL(aVal,lkb.getA()); lkb.setM(mVal); - CHECK_EQUAL(mVal, lkb.getM()); + CHECK_EQUAL(mVal,lkb.getM()); lkb.setD50(d50Val); - CHECK_EQUAL(d50Val, lkb.getD50()); - - if (lkb.init()) - { - CHECK_CLOSE(8.15234192641929420000e-005, lkb.getValue(), tolerance); + CHECK_EQUAL(d50Val,lkb.getD50()); + if(lkb.init()){ + CHECK_CLOSE(8.15234192641929420000e-005,lkb.getValue(),tolerance); } //test RTNTCPRSModel - rs = rttb::models::NTCPRSModel(); + rs=rttb::models::NTCPRSModel(); gammaVal = 1.7; sVal = 1; - CHECK_EQUAL(0, rs.getGamma()); - CHECK_EQUAL(0, rs.getS()); - CHECK_EQUAL(0, rs.getD50()); + CHECK_EQUAL(0,rs.getGamma()); + CHECK_EQUAL(0,rs.getS()); + CHECK_EQUAL(0,rs.getD50()); rs.setDVH(dvhPtr4); - CHECK_EQUAL(dvhPtr4, rs.getDVH()); + CHECK_EQUAL(dvhPtr4,rs.getDVH()); rs.setD50(d50Val); - CHECK_EQUAL(d50Val, rs.getD50()); + CHECK_EQUAL(d50Val,rs.getD50()); rs.setGamma(gammaVal); - CHECK_EQUAL(gammaVal, rs.getGamma()); + CHECK_EQUAL(gammaVal,rs.getGamma()); rs.setS(sVal); - CHECK_EQUAL(sVal, rs.getS()); - - if (rs.init()) - { - CHECK_CLOSE(2.02607985020919480000e-004, rs.getValue(), tolerance); + 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 + //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_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(); + 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 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()) + 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); + 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; + 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); + 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 diff --git a/testing/examples/RTBioModelScatterPlotExampleTest.cpp b/testing/examples/RTBioModelScatterPlotExampleTest.cpp index a1b8c79..ac78bd8 100644 --- a/testing/examples/RTBioModelScatterPlotExampleTest.cpp +++ b/testing/examples/RTBioModelScatterPlotExampleTest.cpp @@ -1,480 +1,480 @@ // ----------------------------------------------------------------------- // 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 = 1; + 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.getAlpahBeta()); 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.getAlpahBeta()); 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.getAlpahBeta()); 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.getAlpahBeta()); 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 diff --git a/testing/models/BioModelTest.cpp b/testing/models/BioModelTest.cpp index 94993f8..f79ef4f 100644 --- a/testing/models/BioModelTest.cpp +++ b/testing/models/BioModelTest.cpp @@ -1,288 +1,288 @@ // ----------------------------------------------------------------------- // 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 "rttbDVH.h" #include "rttbTCPLQModel.h" #include "rttbNTCPLKBModel.h" #include "rttbNTCPRSModel.h" #include "rttbBaseTypeModels.h" #include "rttbBioModelCurve.h" #include "rttbInvalidParameterException.h" #include "rttbBioModelScatterPlots.h" namespace rttb{ namespace testing{ typedef core::DVH::DataDifferentialType DataDifferentialType; /*! @brief RTBioModelTest. TCP calculated using a DVH PTV and LQ Model. NTCP tested using 3 Normal Tissue DVHs and LKB/RS Model. TCPLQ: 1) test constructors (values as expected?) 2) test init (calcTCPxxx) 3) test set/get NTCP (LKB): 1) test constructors (values as expected?) 2) test init (calcxxx) 3) test set/get NTCP (RS): 1) test constructors (values as expected?) 2) test init (calcxxx) 3) test set/get */ int BioModelTest(int argc, char* argv[] ) { PREPARE_DEFAULT_TEST_REPORTING; //generate artificial DVH and corresponding statistical values DoseTypeGy binSize = DoseTypeGy(0.1); DoseVoxelVolumeType voxelVolume = 8; DataDifferentialType aDataDifferential; DoseCalcType value = 0; DVHVoxelNumber numberOfVoxels = 0; // creat default values for(int i = 0; i < 98; i++){ value = 0; numberOfVoxels+=value; aDataDifferential.push_back( value ); } aDataDifferential.push_back( 10 ); aDataDifferential.push_back( 20 ); const IDType structureID = "myStructure"; const IDType doseID = "myDose"; const IDType voxelizationID = "myVoxelization"; core::DVH::DVHPointer dvhPtr = boost::make_shared(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID); //test TCP LQ Model models::BioModelParamType alpha = 0.35; models::BioModelParamType beta = 0.023333333333333; models::BioModelParamType roh = 10000000; - int numFractions = 1; + int numFractions = 8; DoseTypeGy normalizationDose = 50; //1) test constructors (values as expected?) rttb::models::TCPLQModel tcplq=rttb::models::TCPLQModel(); CHECK_EQUAL(0,tcplq.getAlphaMean()); CHECK_EQUAL(0,tcplq.getAlpahBeta()); CHECK_EQUAL(0,tcplq.getRho()); CHECK_EQUAL(0,tcplq.getValue()); tcplq=rttb::models::TCPLQModel(dvhPtr,roh, numFractions,alpha/beta, alpha,0.08); CHECK_EQUAL(alpha,tcplq.getAlphaMean()); CHECK_EQUAL(alpha/beta,tcplq.getAlpahBeta()); CHECK_EQUAL(roh,tcplq.getRho()); CHECK_EQUAL(0,tcplq.getValue()); tcplq=rttb::models::TCPLQModel(); CHECK_EQUAL(0,tcplq.getAlphaMean()); CHECK_EQUAL(0,tcplq.getAlpahBeta()); CHECK_EQUAL(0,tcplq.getRho()); CHECK_EQUAL(0,tcplq.getValue()); tcplq=rttb::models::TCPLQModel(dvhPtr,alpha, beta, roh, numFractions); CHECK_EQUAL(alpha,tcplq.getAlphaMean()); CHECK_EQUAL(alpha/beta,tcplq.getAlpahBeta()); CHECK_EQUAL(roh,tcplq.getRho()); CHECK_EQUAL(0,tcplq.getValue()); //2) test init (calcTCPxxx) CHECK_NO_THROW(tcplq.init(1)); //3) test set/get CHECK_EQUAL(0,tcplq.getValue()); CHECK_NO_THROW(tcplq.setParameters(alpha,10,roh,0.08)); CHECK_EQUAL(10,tcplq.getAlpahBeta()); CHECK_EQUAL(0.08,tcplq.getAlphaVariance()); CHECK_EQUAL(alpha,tcplq.getAlphaMean()); CHECK_EQUAL(roh,tcplq.getRho()); CHECK_NO_THROW(models::getCurveDoseVSBioModel(tcplq,normalizationDose)); std::vector aParameterVector; aParameterVector.push_back(alpha+0.02); CHECK_THROW_EXPLICIT(tcplq.setParameterVector(aParameterVector),core::InvalidParameterException); aParameterVector.push_back(0.06); aParameterVector.push_back(8); aParameterVector.push_back(roh/10); CHECK_NO_THROW(tcplq.setParameterVector(aParameterVector)); CHECK_EQUAL(8,tcplq.getAlpahBeta()); CHECK_EQUAL(0.06,tcplq.getAlphaVariance()); CHECK_EQUAL(alpha+0.02,tcplq.getAlphaMean()); CHECK_EQUAL(roh/10,tcplq.getRho()); for (int i = 0; i < 4; i++){ CHECK_NO_THROW(tcplq.setParameterByID(i,models::BioModelParamType(i))); } CHECK_THROW_EXPLICIT(tcplq.setParameterByID(4,4.0),core::InvalidParameterException); CHECK_EQUAL(0,tcplq.getParameterID("alphaMean")); CHECK_EQUAL(0,tcplq.getAlphaMean()); CHECK_EQUAL(1,tcplq.getParameterID("alphaVariance")); CHECK_EQUAL(1,tcplq.getAlphaVariance()); CHECK_EQUAL(2,tcplq.getParameterID("alpha_beta")); CHECK_EQUAL(2,tcplq.getAlpahBeta()); CHECK_EQUAL(3,tcplq.getParameterID("rho")); CHECK_EQUAL(3,tcplq.getRho()); //test NTCPLKBModel //1) test constructors (values as expected?) models::BioModelParamType aVal = 10; models::BioModelParamType mVal = 0.16; models::BioModelParamType d50Val = 35; CHECK_NO_THROW(rttb::models::NTCPLKBModel()); rttb::models::NTCPLKBModel lkb=rttb::models::NTCPLKBModel(); CHECK_EQUAL(0,lkb.getA()); CHECK_EQUAL(0,lkb.getM()); CHECK_EQUAL(0,lkb.getD50()); CHECK_EQUAL(0,lkb.getValue()); CHECK_NO_THROW(rttb::models::NTCPLKBModel(dvhPtr, d50Val, mVal, aVal)); lkb=rttb::models::NTCPLKBModel(dvhPtr, d50Val, mVal, aVal); CHECK_EQUAL(0,lkb.getValue()); CHECK_EQUAL(dvhPtr,lkb.getDVH()); CHECK_EQUAL(aVal,lkb.getA()); CHECK_EQUAL(mVal,lkb.getM()); CHECK_EQUAL(d50Val,lkb.getD50()); //2) test init (calcxxx) CHECK_NO_THROW(lkb.init(1)); lkb.getValue(); //3) test set/get lkb=rttb::models::NTCPLKBModel(); CHECK_EQUAL(0,lkb.getA()); CHECK_EQUAL(0,lkb.getM()); CHECK_EQUAL(0,lkb.getD50()); lkb.setDVH(dvhPtr); CHECK_EQUAL(dvhPtr,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(models::getCurveEUDVSBioModel(lkb)); aParameterVector.clear(); aParameterVector.push_back(d50Val+5); CHECK_THROW_EXPLICIT(lkb.setParameterVector(aParameterVector),core::InvalidParameterException); aParameterVector.push_back(mVal+0.2); aParameterVector.push_back(aVal+0.5); CHECK_NO_THROW(lkb.setParameterVector(aParameterVector)); CHECK_EQUAL(aVal+0.5,lkb.getA()); CHECK_EQUAL(mVal+0.2,lkb.getM()); CHECK_EQUAL(d50Val+5,lkb.getD50()); for (int i = 0; i < 3; i++){ CHECK_NO_THROW(lkb.setParameterByID(i,models::BioModelParamType(i))); } CHECK_THROW_EXPLICIT(lkb.setParameterByID(3,4.0),core::InvalidParameterException); CHECK_EQUAL(0,lkb.getParameterID("d50")); CHECK_EQUAL(0,lkb.getD50()); CHECK_EQUAL(1,lkb.getParameterID("m")); CHECK_EQUAL(1,lkb.getM()); CHECK_EQUAL(2,lkb.getParameterID("a")); CHECK_EQUAL(2,lkb.getA()); //test NTCPRSModel //1) test constructors (values as expected?) CHECK_NO_THROW(rttb::models::NTCPRSModel()); models::BioModelParamType gammaVal = 1.7; models::BioModelParamType sVal = 1; CHECK_NO_THROW(rttb::models::NTCPRSModel(dvhPtr,d50Val,gammaVal,sVal)); rttb::models::NTCPRSModel rs=rttb::models::NTCPRSModel(dvhPtr,d50Val,gammaVal,sVal); CHECK_EQUAL(dvhPtr,rs.getDVH()); CHECK_EQUAL(d50Val,rs.getD50()); CHECK_EQUAL(gammaVal,rs.getGamma()); CHECK_EQUAL(sVal,rs.getS()); rs=rttb::models::NTCPRSModel(); CHECK_EQUAL(0,rs.getGamma()); CHECK_EQUAL(0,rs.getS()); CHECK_EQUAL(0,rs.getD50()); //3) test set/get rs.setDVH(dvhPtr); CHECK_EQUAL(dvhPtr,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()); //2) test init (calcxxx) CHECK_NO_THROW(rs.init(1)); //3) test set/get continued aParameterVector.clear(); aParameterVector.push_back(d50Val+5); CHECK_THROW_EXPLICIT(rs.setParameterVector(aParameterVector),core::InvalidParameterException); aParameterVector.push_back(gammaVal+0.2); aParameterVector.push_back(sVal+0.5); CHECK_NO_THROW(rs.setParameterVector(aParameterVector)); CHECK_EQUAL(gammaVal+0.2,rs.getGamma()); CHECK_EQUAL(sVal+0.5,rs.getS()); CHECK_EQUAL(d50Val+5,rs.getD50()); for (int i = 0; i < 3; i++){ CHECK_NO_THROW(rs.setParameterByID(i,models::BioModelParamType(i))); } CHECK_THROW_EXPLICIT(rs.setParameterByID(3,4.0),core::InvalidParameterException); CHECK_EQUAL(0,rs.getParameterID("d50")); CHECK_EQUAL(0,rs.getD50()); CHECK_EQUAL(1,rs.getParameterID("gamma")); CHECK_EQUAL(1,rs.getGamma()); CHECK_EQUAL(2,rs.getParameterID("s")); CHECK_EQUAL(2,rs.getS()); //Scatter plot tests CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq,0,alpha,0,normalizationDose));//variance=0, will be set to 1e-30 CHECK_THROW_EXPLICIT(models::getScatterPlotVary1Parameter(tcplq,0,alpha,alpha*0.1,0),core::InvalidParameterException);//normalisationdose=0 CHECK_THROW_EXPLICIT(models::getScatterPlotVary1Parameter(tcplq,0,alpha,alpha*0.1,normalizationDose,10000,0,0),core::InvalidParameterException);//maxDose-minDose=0 RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb \ No newline at end of file diff --git a/testing/models/CMakeLists.txt b/testing/models/CMakeLists.txt index af42f50..504fd5f 100644 --- a/testing/models/CMakeLists.txt +++ b/testing/models/CMakeLists.txt @@ -1,20 +1,21 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(MODELS_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbModelTests) SET(MODELS_HEADER_TEST ${EXECUTABLE_OUTPUT_PATH}/rttbModelsHeaderTest) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(BioModelTest ${MODELS_TESTS} BioModelTest) ADD_TEST(BioModelScatterPlotTest ${MODELS_TESTS} BioModelScatterPlotTest) +ADD_TEST(DvhBasedModelsTest ${MODELS_TESTS} DvhBasedModelsTest) RTTB_CREATE_TEST_MODULE(rttbModel DEPENDS RTTBModels RTTBOtherIO PACKAGE_DEPENDS Boost Litmus DCMTK) diff --git a/testing/models/DvhBasedModelsTest.cpp b/testing/models/DvhBasedModelsTest.cpp new file mode 100644 index 0000000..ce96902 --- /dev/null +++ b/testing/models/DvhBasedModelsTest.cpp @@ -0,0 +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: 4047 $ (last changed revision) +// @date $Date: 2012-10-29 16:19:15 +0100 (Mo, 29 Okt 2012) $ (last change date) +// @author $Author: mang $ (last changed by) +*/ +#include +#include + +#include "litCheckMacros.h" + +#include "rttbDVH.h" +#include "rttbDvhBasedModels.h" +#include "rttbInvalidParameterException.h" + +namespace rttb{ + namespace testing{ + + typedef core::DVH::DataDifferentialType DataDifferentialType; + + /*! @brief DvhBasedModelsTest. + 1) Test bed und lqed2 + */ + int DvhBasedModelsTest(int argc, char* argv[] ) + { + PREPARE_DEFAULT_TEST_REPORTING; + + //1) test calcBEDDVH and calcLQED2DVH + //generate artificial DVH and corresponding statistical values + DoseTypeGy binSize = DoseTypeGy(0.1); + DoseVoxelVolumeType voxelVolume = 8; + const IDType structureID = "myStructure"; + const IDType doseID = "myDose"; + const IDType voxelizationID = "myVoxelization"; + DataDifferentialType aDataDifferential; + std::vector bedVector; + std::vector lqed2Vector; + int numberOfFractions = 2; + double alpha_beta = 10; + for(int i=0; i<100; i++){ + double volume = DoseCalcType((double(rand())/RAND_MAX)*1000); + double dose = (i+0.5)* binSize; + aDataDifferential.push_back( volume); + bedVector.push_back(dose * (1 + dose / (numberOfFractions * alpha_beta))); + lqed2Vector.push_back(dose * ((alpha_beta + (dose / numberOfFractions)) / (alpha_beta+2))); + } + core::DVH myDVH(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID); + core::DVH::DVHPointer dvhPtr = boost::make_shared(myDVH); + + CHECK_THROW_EXPLICIT(rttb::models::calcBEDDVH(dvhPtr, 0, 10), core::InvalidParameterException); + CHECK_THROW_EXPLICIT(rttb::models::calcBEDDVH(dvhPtr, 10,-1), core::InvalidParameterException); + CHECK_NO_THROW(rttb::models::calcBEDDVH(dvhPtr,10,10)); + CHECK_EQUAL(rttb::models::calcBEDDVH(dvhPtr,2,10).size(), myDVH.getDataDifferential().size()); + rttb::models::BEDDVHType bedDVH = rttb::models::calcBEDDVH(dvhPtr,numberOfFractions, alpha_beta); + + CHECK_THROW_EXPLICIT(rttb::models::calcLQED2DVH(dvhPtr,1,10), core::InvalidParameterException); + CHECK_THROW_EXPLICIT(rttb::models::calcLQED2DVH(dvhPtr,10,-1), core::InvalidParameterException); + CHECK_NO_THROW(rttb::models::calcLQED2DVH(dvhPtr,10,10,true)); + CHECK_EQUAL(rttb::models::calcLQED2DVH(dvhPtr,2,10).size(), myDVH.getDataDifferential().size()); + rttb::models::BEDDVHType lqed2DVH = rttb::models::calcLQED2DVH(dvhPtr,numberOfFractions, alpha_beta); + + //check the calculation + rttb::models::BEDDVHType::iterator itBED, itLQED2; + std::vector::iterator itBEDVec, itLQED2Vec; + DataDifferentialType::iterator itDiff; + for(itBED = bedDVH.begin(), itLQED2 = lqed2DVH.begin(), itBEDVec = bedVector.begin(), itLQED2Vec = lqed2Vector.begin(), itDiff = aDataDifferential.begin(); + itBED != bedDVH.end(), itLQED2 != lqed2DVH.end(), itBEDVec != bedVector.end(), itLQED2Vec != lqed2Vector.end(), itDiff != aDataDifferential.end(); + ++itBED, ++itLQED2, ++itBEDVec, ++itLQED2Vec, ++itDiff){ + + //check volume + CHECK_EQUAL(*itDiff, (*itBED).second); + CHECK_EQUAL((*itBED).second, (*itLQED2).second); + + //check bed + CHECK_EQUAL(*itBEDVec, (*itBED).first); + + //check lqed2 + CHECK_EQUAL(*itLQED2Vec, (*itLQED2).first); + } + + RETURN_AND_REPORT_TEST_SUCCESS; + + } + + }//testing +}//rttb \ No newline at end of file diff --git a/testing/models/files.cmake b/testing/models/files.cmake index 18b857c..1e5c4fa 100644 --- a/testing/models/files.cmake +++ b/testing/models/files.cmake @@ -1,14 +1,15 @@ SET(CPP_FILES rttbModelsTest.cpp BioModelTest.cpp rttbScatterTester.cpp BioModelScatterPlotTest.cpp DummyModel.cpp ../core/DummyDVHGenerator.cpp + DvhBasedModelsTest.cpp ) SET(H_FILES rttbScatterTester.h DummyModel.h ../core/DummyDVHGenerator.h ) diff --git a/testing/models/rttbModelsTest.cpp b/testing/models/rttbModelsTest.cpp index f581daa..a285790 100644 --- a/testing/models/rttbModelsTest.cpp +++ b/testing/models/rttbModelsTest.cpp @@ -1,62 +1,63 @@ // ----------------------------------------------------------------------- // 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 #if defined(_MSC_VER) #pragma warning ( disable : 4786 ) #endif #include "litMultiTestsMain.h" namespace rttb{ namespace testing{ void registerTests() { LIT_REGISTER_TEST(BioModelTest); LIT_REGISTER_TEST(BioModelScatterPlotTest); + LIT_REGISTER_TEST(DvhBasedModelsTest); } } } int main(int argc, char* argv[]) { int result = 0; rttb::testing::registerTests(); try { result = lit::multiTestsMain(argc,argv); } catch(const std::exception& e) { result = -1; } catch(...) { result = -1; } return result; }