diff --git a/code/models/rttbBaseTypeModels.h b/code/models/rttbBaseTypeModels.h index 07285cd..d5d6370 100644 --- a/code/models/rttbBaseTypeModels.h +++ b/code/models/rttbBaseTypeModels.h @@ -1,52 +1,46 @@ // ----------------------------------------------------------------------- // 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_MODEL_H #define __BASE_TYPE_MODEL_H #include #include "rttbBaseType.h" namespace rttb { namespace models { using BioModelParamType = double; using BioModelValueType = double; const double infinity = 1e30; struct TcpParams { double alphaMean; double alphaVariance; double rho; std::vector volumeVector; std::vector bedVector; }; } } #endif \ No newline at end of file diff --git a/code/models/rttbBioModel.cpp b/code/models/rttbBioModel.cpp index 4a4b9f7..3417517 100644 --- a/code/models/rttbBioModel.cpp +++ b/code/models/rttbBioModel.cpp @@ -1,62 +1,55 @@ // ----------------------------------------------------------------------- // 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 "rttbBioModel.h" - namespace rttb { namespace models { bool BioModel::init(const double doseFactor) { if (_dvh != nullptr) { _value = this->calcModel(doseFactor); return true; } return false; } void BioModel::setDVH(const DVHPointer aDVH) { _dvh = aDVH; } const BioModel::DVHPointer BioModel::getDVH() const { return _dvh; } const BioModelValueType BioModel::getValue() const { return _value; } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbBioModel.h b/code/models/rttbBioModel.h index 3d0d091..6fef77d 100644 --- a/code/models/rttbBioModel.h +++ b/code/models/rttbBioModel.h @@ -1,110 +1,104 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) -*/ #ifndef __BIO_MODEL_H #define __BIO_MODEL_H #include "rttbDVH.h" #include "rttbBaseTypeModels.h" #include #include "RTTBModelsExports.h" #include namespace rttb { namespace models { /*! @class BioModel @brief This is the interface class for biological models */ class RTTBModels_EXPORT BioModel { public: rttbClassMacroNoParent(BioModel) using ParamVectorType = std::vector; using DVHPointer = core::DVH::Pointer; protected: DVHPointer _dvh; BioModelValueType _value{0}; /*! @brief Calculate the model value @param doseFactor scaling factor for the dose. The model calculation will use the dvh with each di=old di*doseFactor. */ virtual BioModelValueType calcModel(const double doseFactor = 1) = 0; /* Map of all parameters */ std::map parameterMap; std::string _name; public: BioModel() = default; BioModel(DVHPointer aDvh): _dvh(aDvh), _value(0) {}; virtual ~BioModel()= default; /*! @brief Start the calculation. If any parameter changed, init() should be called again and return =true before getValue() is called! @return Return true if successful */ bool init(const double doseFactor = 1); /*! @param aDVH must be a DVH calculated by a cumulative dose distribution, not a fraction DVH! */ void setDVH(const DVHPointer aDVH); const DVHPointer getDVH() const; /*! @brief Set parameter vector, where index of vector is the parameter ID. */ virtual void setParameterVector(const ParamVectorType& aParameterVector) = 0; virtual void setParameterByID(const int aParamId, const BioModelParamType aValue) = 0; /*! @brief Get parameter by ID. @return Return -1 if ID is not found. */ virtual const int getParameterID(const std::string& aParamName) const = 0; virtual std::map getParameterMap() const = 0; virtual void fillParameterMap() = 0 ; virtual std::string getModelType() const = 0; /*! @brief Get the value of biological model @pre init() must be called and =true! */ const BioModelValueType getValue() const; }; }//end namespace models }//end namespace rttb #endif diff --git a/code/models/rttbBioModelCurve.cpp b/code/models/rttbBioModelCurve.cpp index a4565d2..24ca50c 100644 --- a/code/models/rttbBioModelCurve.cpp +++ b/code/models/rttbBioModelCurve.cpp @@ -1,93 +1,87 @@ // ----------------------------------------------------------------------- // 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 "rttbBioModel.h" #include "rttbNTCPLKBModel.h" #include "rttbBioModelCurve.h" #include "rttbDvhBasedModels.h" namespace rttb { namespace models { CurveDataType getCurveDoseVSBioModel(BioModel& aModel, double normalisationDose, int aBin, double minDose, double maxDose) { CurveDataType _curveData; double minFactor = 0, maxFactor = 0; minFactor = minDose / normalisationDose; maxFactor = maxDose / normalisationDose; double lastValue = 0; for (int i = 0; i < aBin; i++) { double factor = minFactor + i * (maxFactor - minFactor) / (aBin - 1); if (lastValue < 1) { aModel.init(factor); BioModelValueType value = aModel.getValue(); lastValue = value; _curveData.insert(CurvePointType(factor * normalisationDose, value)); } else { _curveData.insert(CurvePointType(factor * normalisationDose, 1)); } } return _curveData; } CurveDataType getCurveEUDVSBioModel(NTCPLKBModel& aModel, DoseCalcType maxFactor, DoseCalcType minFactor, int aBin) { CurveDataType _curveData; BioModel::DVHPointer _dvh = aModel.getDVH(); for (int i = 0; i < aBin; i++) { DoseCalcType factor = minFactor + i * (maxFactor - minFactor) / (aBin - 1); core::DVH variantDVH = core::DVH(_dvh->getDataDifferential(), _dvh->getDeltaD() * factor, _dvh->getDeltaV(), "temporary", "temporary"); auto spDVH = boost::make_shared(variantDVH); double eud = getEUD(spDVH, aModel.getA()); aModel.init(factor); BioModelValueType value = aModel.getValue(); _curveData.insert(CurvePointType(eud, value)); } return _curveData; } } } \ No newline at end of file diff --git a/code/models/rttbBioModelCurve.h b/code/models/rttbBioModelCurve.h index 8bd6f5a..7e4e44d 100644 --- a/code/models/rttbBioModelCurve.h +++ b/code/models/rttbBioModelCurve.h @@ -1,55 +1,49 @@ // ----------------------------------------------------------------------- // 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 __MODEL_CURVE_H #define __MODEL_CURVE_H namespace rttb { namespace models { class BioModel; class NTCPLKBModel; //map of dose value and model value typedef std::map CurveDataType; //pair of dose value and model value typedef std::pair CurvePointType; /*! @brief Get the curve TCP/NTCP Value vs normalisationDose, normalisationDose variant between minDose and maxDose. @param aBin the size of the map @param minDose min value for x axis @param maxDose max value for x axis @param normalisationDose prescribed dose of the current _dvh or mean/maximum. */ CurveDataType getCurveDoseVSBioModel(BioModel& aModel, double normalisationDose, int aBin = 201, double minDose = 0.1, double maxDose = 150); /*! @brief Get the curve NTCP Value vs EUD, dvh variant between minFactor*deltaD and maxFactor*deltaD. @param aBin the size of the map @param minFactor min factor for dvh deltaD @param maxFactor max factor for dvh deltaD */ CurveDataType getCurveEUDVSBioModel(NTCPLKBModel& aModel, DoseCalcType maxFactor = 10, DoseCalcType minFactor = 0.1, int aBin = 201); } } #endif \ No newline at end of file diff --git a/code/models/rttbBioModelScatterPlots.cpp b/code/models/rttbBioModelScatterPlots.cpp index b33a133..67565ef 100644 --- a/code/models/rttbBioModelScatterPlots.cpp +++ b/code/models/rttbBioModelScatterPlots.cpp @@ -1,220 +1,214 @@ // ----------------------------------------------------------------------- // 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 "rttbBioModel.h" #include "rttbBioModelScatterPlots.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace models { /* Initiate Random Number generator with current time */ boost::random::mt19937 rng(static_cast(time(nullptr))); /* Generate random number between 0 and 1 */ boost::random::uniform_01 uniDist(rng); ScatterPlotType getScatterPlotVary1Parameter(BioModel& aModel, int aParamId, BioModelParamType aMean, BioModelParamType aVariance, DoseTypeGy aNormalisationDose, int numberOfPoints, DoseTypeGy aMinDose, DoseTypeGy aMaxDose) { ScatterPlotType scatterPlotData; if (aVariance == 0) { //set to small positive value to avoid negative infinity! aVariance = 1e-30; } if (aMaxDose <= aMinDose) { throw core::InvalidParameterException("Parameter invalid: aMaxDose must be > aMinDose!"); } if (aNormalisationDose <= 0) { throw core::InvalidParameterException("Parameter invalid: aNormalisationDose must be > 0!"); } /* Choose Normal Distribution */ boost::random::normal_distribution gaussian_dist(0, aVariance); /* Create a Gaussian Random Number generator * by binding with previously defined * normal distribution object */ boost::random::variate_generator > generator( rng, gaussian_dist); int i = 0; while (i < numberOfPoints) { double paramValue, probability; double randomValue = generator(); paramValue = randomValue + aMean; probability = normal_pdf(randomValue, aVariance); if (probability > 1) { probability = 1; } //randomly select a dose between aMinDose and aMaxDose double dose = uniDist() * (aMaxDose - aMinDose) + aMinDose; if (probability > 0) { try { aModel.setParameterByID(aParamId, paramValue); aModel.init(dose / aNormalisationDose); double value = aModel.getValue(); std::pair modelProbPair = std::make_pair(value, probability); scatterPlotData.insert(std::pair >(dose, modelProbPair)); i++; } catch (core::InvalidParameterException& /*e*/) { //repeat evaluation to guarantee the correct number of scatter values continue; } } } return scatterPlotData; } double normal_pdf(double aValue, double aVariance) { static const double inv_sqrt_2pi = 0.3989422804014327; double a = (aValue) / aVariance; return inv_sqrt_2pi / aVariance * std::exp(-0.5f * a * a); } ScatterPlotType getScatterPlotVaryParameters(BioModel& aModel, std::vector aParamIdVec, BioModel::ParamVectorType aMeanVec, BioModel::ParamVectorType aVarianceVec, DoseTypeGy aNormalisationDose, int numberOfPoints, DoseTypeGy aMinDose, DoseTypeGy aMaxDose) { ScatterPlotType scatterPlotData; if (aMaxDose <= aMinDose) { throw core::InvalidParameterException("Parameter invalid: aMaxDose must be > aMinDose!"); } if (aNormalisationDose <= 0) { throw core::InvalidParameterException("Parameter invalid: aNormalisationDose must be > 0!"); } //all input vectors need to have the same size if (((aVarianceVec.size() != aMeanVec.size()) || (aVarianceVec.size() != aParamIdVec.size()))) { throw core::InvalidParameterException("Parameter vectors have different sizes!"); } for (double & v : aVarianceVec) { //set to small positive value to avoid negative infinity! if (v == 0) { v = 1e-30; } } double paramValue; // vary all parameters for each scattered point int i = 0; while (i < numberOfPoints) { double probability = 1; for (GridIndexType j = 0; j < aParamIdVec.size(); j++) { /* Choose Normal Distribution */ boost::random::normal_distribution gaussian_dist(0, aVarianceVec.at(j)); /* Create a Gaussian Random Number generator * by binding with previously defined * normal distribution object */ boost::random::variate_generator > generator( rng, gaussian_dist); double randomValue = generator(); paramValue = randomValue + aMeanVec.at(j); if (aVarianceVec.at(j) != 0) { /* calculate combined probability */ probability = probability * normal_pdf(randomValue, aVarianceVec.at(j)); } else { throw core::InvalidParameterException("Parameter invalid: Variance should not be 0!"); } aModel.setParameterByID(aParamIdVec.at(j), paramValue); } //randomly select a dose between aMinDose and aMaxDose double dose = uniDist() * (aMaxDose - aMinDose) + aMinDose; if (probability > 0) { try { aModel.init(dose / aNormalisationDose); double value = aModel.getValue(); std::pair modelProbPair = std::make_pair(value, probability); scatterPlotData.insert(std::pair >(dose, modelProbPair)); i++; } catch (core::InvalidParameterException& /*e*/) { //repeat evaluation to guarantee the correct number of scatter values continue; } } } return scatterPlotData; } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbBioModelScatterPlots.h b/code/models/rttbBioModelScatterPlots.h index ea3d31b..9268b63 100644 --- a/code/models/rttbBioModelScatterPlots.h +++ b/code/models/rttbBioModelScatterPlots.h @@ -1,96 +1,90 @@ // ----------------------------------------------------------------------- // 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 __MODEL_SCATTER_H #define __MODEL_SCATTER_H #include "rttbBaseType.h" #include "rttbBaseTypeModels.h" namespace rttb { namespace models { class BioModel; // maps dose to a pair of model value and probability typedef std::multimap > ScatterPlotType; using ParamVectorType = std::vector; /*! @brief Get the points (TCP/NTCP Value, probability of the value) if 1 parameter vary from a normal- distribution with mean=aMean, variance=aVariance. @param aModel biological model for which the scatter plot will be generated @param aParamId ID of the parameter to be varied to generate the scatter plot @param aMean mean value for the distribution of the varied parameter @param aVariance variance of the varied parameter. The variance may not be exactly zero. If so, it is set to 1e-30 to avoid numerical instability. @param aNormalisationDose prescribed dose of the current _dvh @param numberOfPoints the size of the map, number of points to be calculated @param aMinDose dose will be randomly selected from [aMinDose] (uniform distribution). They will define the minvalue for x axis @param aMaxDose dose will be randomly selected from [aMaxDose] (uniform distribution). They will define the max value for x axis @return Map of scattered values. If all parameters are valid, this map contains numberOfPoints valid scatter values. If aMaxDose<=aMinDose, the scatter plot cannot be generated. The map will therefore be empty. @warning This method is slow, do not use with too many points. Because the scatter plot map must contain numberOfPoints the scatter plot generation may run more often (producing invalid values). In tests the generation process runs on average approximately 20% more often. @exception InvalidParameterException Thrown if aNormalisationDose<=0 or aMinDose<=aMaxiDose */ ScatterPlotType getScatterPlotVary1Parameter(BioModel& aModel, int aParamId, BioModelParamType aMean, BioModelParamType aVariance, DoseTypeGy aNormalisationDose, int numberOfPoints = 100, DoseTypeGy aMinDose = 0, DoseTypeGy aMaxDose = 150); /*! @brief Get the points (TCP/NTCP Value, probability of the value) if >=1 parameter vary from a normal- distribution with mean of parameter aParamIdVec.at(i)=aMeanVec.at(i), variance of parameter aParamIdVec.at(i)= aVarianceVec.at(i). @param aModel biological model for which the scatter plot will be generated @param aParamIdVec a vector containing the IDs of the parameters to be varied to generate the scatter plot @param aMeanVec a vector of mean values for the distribution of individually the varied parameters @param aVarianceVec a vector of variance values of the individually varied parameter. The variance may not be exactly zero for any parameter. If so, it is set to 1e-30 to avoid numerical instability. @param aNormalisationDose prescribed dose of the current _dvh @param numberOfPoints the size of the map, number of points to be calculated @param aMinDose dose will be randomly selected from [aMinDose] (uniform distribution). They will define the min value for x axis @param aMaxDose dose will be randomly selected from [aMaxDose] (uniform distribution). They will define the max value for x axis @throw InvalidParameterException is thrown if the parameter vectors do not have the same size. @return Map of scattered values. If all parameters are valid, this map contains numberOfPoints valid scatter values. If aMaxDose<=aMinDose, the scatter plot cannot be generated. The map will therefore be empty. @warning This method is very slow do not use with too many points. Because the scatter plot map must contain numberOfPoints the scatter plot generation may run more often (producing invalid values). In tests the generation process runs on average approximately 20% more often. @exception InvalidParameterException Thrown if aNormalisationDose<=0 or aMinDose<=aMaxiDose */ ScatterPlotType getScatterPlotVaryParameters(BioModel& aModel, std::vector aParamIdVec, ParamVectorType aMeanVec, ParamVectorType aVarianceVec, DoseTypeGy aNormalisationDose, int numberOfPoints = 50, DoseTypeGy aMinDose = 0, DoseTypeGy aMaxDose = 150); /*! Compute normal probability density function for zero mean at aValue with aVariance. */ double normal_pdf(double aValue, double aVariance); } } #endif \ No newline at end of file diff --git a/code/models/rttbDoseBasedModels.cpp b/code/models/rttbDoseBasedModels.cpp index b473c4a..ad62861 100644 --- a/code/models/rttbDoseBasedModels.cpp +++ b/code/models/rttbDoseBasedModels.cpp @@ -1,41 +1,36 @@ // ----------------------------------------------------------------------- // 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 "rttbDoseBasedModels.h" #include "rttbInvalidParameterException.h" #include namespace rttb { namespace models { rttb::models::BioModelValueType calcLQ(DoseTypeGy dose, DoseCalcType alpha, DoseCalcType beta, unsigned int nFractions) { if (dose < 0 || alpha < 0 || beta < 0) { throw core::InvalidParameterException("Parameter invalid: dose, alpha, beta must be >=0!"); } return exp(-((alpha * dose) + (beta * dose * dose / DoseCalcType(nFractions)))); } } } \ No newline at end of file diff --git a/code/models/rttbDoseBasedModels.h b/code/models/rttbDoseBasedModels.h index fa9a7d6..0101de3 100644 --- a/code/models/rttbDoseBasedModels.h +++ b/code/models/rttbDoseBasedModels.h @@ -1,43 +1,37 @@ // ----------------------------------------------------------------------- // 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 "rttbBaseTypeModels.h" +#include "rttbBaseTypeModels.h" namespace rttb { namespace models { /*! @brief Calculate biological LinearQuadratic Model of a dose @details \f$LQ = exp(-\frac{alpha*d+beta*d^2}{nFractions})\f$ @param dose @param alpha @param beta @param nFractions the number of fractions @pre dose>=0 @pre alpha>=0 @pre beta>=0 @return The LQ value @exception rttb::core::InvalidParameterException Thrown if parameters were not set correctly. */ BioModelValueType calcLQ(DoseTypeGy dose, DoseCalcType alpha, DoseCalcType beta, unsigned int nFractions=1); } } \ No newline at end of file diff --git a/code/models/rttbDvhBasedModels.cpp b/code/models/rttbDvhBasedModels.cpp index f72f97b..0eb56c2 100644 --- a/code/models/rttbDvhBasedModels.cpp +++ b/code/models/rttbDvhBasedModels.cpp @@ -1,177 +1,172 @@ // ----------------------------------------------------------------------- // 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(core::DVH::ConstPointer 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 < dataDifferential.size(); i++) { double doseGyi = (i + 0.5) * deltaD; double relativeVolumei = dataDifferential[i] / numberOfVoxels; eud += pow((double)doseGyi, (double)aA) * relativeVolumei; } eud = pow((double)eud, (double)(1 / aA)); return eud; } BEDDVHType calcBEDDVH(core::DVH::ConstPointer 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 >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(core::DVH::ConstPointer 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 875f927..7d5334e 100644 --- a/code/models/rttbDvhBasedModels.h +++ b/code/models/rttbDvhBasedModels.h @@ -1,50 +1,65 @@ +// ----------------------------------------------------------------------- +// 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. +// +//------------------------------------------------------------------------ + #include #include "rttbDVH.h" #include "rttbBaseType.h" namespace rttb { namespace models { using DataDifferentialType = core::DVH::DataDifferentialType; typedef std::map BEDDVHType; typedef std::map LQEDDVHType; /*! @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. */ DoseStatisticType getEUD(core::DVH::ConstPointer dvh, const DoseCalcType aA); /*! @brief Calculate Biological Effective/Equivalent Dose (BED) of dvh @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. */ BEDDVHType calcBEDDVH(core::DVH::ConstPointer 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 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. */ LQEDDVHType calcLQED2DVH(core::DVH::ConstPointer dvh, const int numberOfFractions, const DoseCalcType alpha_beta, const bool relativeVolume = false); } } \ No newline at end of file diff --git a/code/models/rttbIntegration.cpp b/code/models/rttbIntegration.cpp index 2ec03c5..115467c 100644 --- a/code/models/rttbIntegration.cpp +++ b/code/models/rttbIntegration.cpp @@ -1,192 +1,186 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ -/*! -// @file -// @version $Revision $ (last changed revision) -// @date $Date $ (last change date) -// @author $Author $ (last changed by) -*/ + #include #include #include "rttbIntegration.h" #include "rttbInvalidParameterException.h" - namespace rttb { namespace models { double tcpModelFunctor::calculate(double x) const { if (x == 0) { x = 1e-30; } return tcpFunction(a + (1 - x) / x, params) / (x * x); } double LkbModelFunctor::calculate(double x) const { if (x == 0) { x = 1e-30; } return lkbFunction(b - (1 - x) / x) / (x * x); } double tcpFunction(double x, const TcpParams& tcp_params) { double alphaVariance = tcp_params.alphaVariance; if (alphaVariance == 0) { alphaVariance = 1e-30; } double f = exp(-pow((x - tcp_params.alphaMean) / alphaVariance, 2) / 2); for (size_t i = 0; i < tcp_params.volumeVector.size(); ++i) { double tmp, tmp1, tmp2, tmp3; tmp1 = exp(-x * tcp_params.bedVector.at(i)); tmp2 = -(tcp_params.rho) * tcp_params.volumeVector.at(i); tmp3 = tmp2 * tmp1; tmp = exp(tmp3); if (tmp != 1) { f = f * tmp; } } return f; } double integrateTCP(double a, const TcpParams& params) { double aNew = 1e-30; double bNew = 1.0; tcpModelFunctor BMFunction; BMFunction.params = params; BMFunction.a = a; return iterativeIntegration(BMFunction, aNew, bNew); } double lkbFunction(double x) { double tmp = -pow(x, 2) / 2; double step = exp(tmp); return step; } double integrateLKB(double b) { double aNew = 1e-30; double bNew = 1.0; LkbModelFunctor BMFunction; BMFunction.b = b; return iterativeIntegration(BMFunction, aNew, bNew); } //returns the nth stage of refinement of the extended trapezoidal rule template integrationType trapzd(const FunctorType& BMfunction, integrationType a, integrationType b, int stepNum) { static integrationType result; if (stepNum == 1) { result = 0.5 * (b - a) * (BMfunction.calculate(a) + BMfunction.calculate(b)); } else { integrationType x, tnm, sum, del; int it, j; for (it = 1, j = 1; j < stepNum - 1; j++) { it <<= 1; } tnm = it; del = (b - a) / tnm; x = a + 0.5 * del; for (sum = 0.0, j = 0; j < it; j++, x += del) { sum += BMfunction.calculate(x); } result = 0.5 * (result + (b - a) * sum / tnm); } return result; } template integrationType iterativeIntegration(const FunctorType& BMfunction, integrationType a, integrationType b) { integrationType ost = 0.0; integrationType os = 0.0; unsigned int maxSteps = 50; double eps = 1e-6; unsigned int i = 1; for (; i <= maxSteps; ++i) { integrationType st = trapzd(BMfunction, a, b, i); integrationType s = (4.0 * st - ost) / 3.0; if (i > 5) { if (fabs(s - os) < eps * fabs(os) || (s == 0.0 && os == 0.0)) { return s; } } os = s; ost = st; } //too many iterations, this should never be reachable! throw rttb::core::InvalidParameterException("Integral calculation failed: too many iterations! "); } } } diff --git a/code/models/rttbIntegration.h b/code/models/rttbIntegration.h index 5878aa0..23b7ece 100644 --- a/code/models/rttbIntegration.h +++ b/code/models/rttbIntegration.h @@ -1,119 +1,114 @@ // ----------------------------------------------------------------------- // 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 __INTEGRATION_H #define __INTEGRATION_H #include "rttbBaseTypeModels.h" namespace rttb { namespace models { using integrationType = double; /*! @class LkbModelFunctor @brief A FunctorType: calculate the transformed LKB-Model @details LBK Model is calculated using the transformation \f$x = b - \frac{1-t}{t}\f$. \f$\int_{-\infty}^{b} dx\f$, \f$f(x) = \int_0^1 dt \frac{f(b - \frac{(1-t)}{t})}{t^2}\f$ */ class LkbModelFunctor { public: /*!b: upper bound of the lkb integration*/ double b; /*! calculate the transformed LKB-Model using the transformation \f$x = b - \frac{1-t}{t}.\f$ */ double calculate(double x) const; }; /*! @class tcpModelFunctor @brief A FunctorType: calculate the transformed TCP-LQ-Model @details LBK Model is calculated using the transformation \f$x = a + \frac{1-t}{t}\f$. \f$\int_{a}^{+\infty} dx f(x) =\int_0^1 dt \frac{f(a + \frac{(1-t)}{t}}{t^2}\f$ */ class tcpModelFunctor { public: //TCP parameters including alpha mean, alpha variance, rho, volume vector and bed vector TcpParams params; /*! a: lower bound of the tcp integration*/ double a; /*! calculate the transformed TCP-LQ-Model using the transformation \f$x = a + \frac{1-t}{t}\f$. */ double calculate(double x) const; }; /*! @brief Function to be integrated for TCP LQ model. @param x: variable of the TCP LQ model @param tcp_params: TCP parameters including alpha mean, alpha variance, rho, volume vector and bed vector @return Return the function value */ double tcpFunction(double x, const TcpParams& tcp_params); /*! @brief Compute integration step for \f$f(x) = exp(-\frac{t^2}{2})\f$. @param x: variable of the lkb function @return Return the function value */ double lkbFunction(double x); /*! @brief Calculate LKB Integration over \f$(-\infty,b)\f$. The integral is mapped onto the semi-open interval \f$(0,1]\f$ using the transformation \f$x = b - \frac{1-t}{t}\f$ @param b: upper bound of the lkb integration */ double integrateLKB(double b); /*! @brief Calculate TCP integration over \f$(a, \infty)\f$. The integral is mapped onto the semi-open interval \f$(0,1]\f$ using the transformation \f$x = a + \frac{1-t}{t}\f$ @param a: lower bound of the tcp integration @param params: TCP parameters including alpha mean, alpha variance, rho, volume vector and bed vector */ double integrateTCP(double a, const TcpParams& params); /* @brief This function returns the nth stage of refinement of the extended trapezoidal rule. @param BMfunction: function to be integrated, for example a LkbModelFunctor or a tcpModelFunctor @param a: lower bound of the integral @param b: upper bound of the integral @param stepNum: the nth stage @param result: the current result */ template integrationType trapzd(const FunctorType& BMfunction, integrationType a, integrationType b, int stepNum); /*! @brief Iterative integration routine @param BMfunction: function to be integrated, for example a LkbModelFunctor or a tcpModelFunctor @param a: lower bound of the integral @param b: upper bound of the integral @exception throw InvalidParameterException if integral calculation failed. */ template integrationType iterativeIntegration(const FunctorType& BMfunction, integrationType a, integrationType b); } } #endif \ No newline at end of file diff --git a/code/models/rttbLQModelAccessor.cpp b/code/models/rttbLQModelAccessor.cpp index edc564d..850e246 100644 --- a/code/models/rttbLQModelAccessor.cpp +++ b/code/models/rttbLQModelAccessor.cpp @@ -1,100 +1,94 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) -*/ #include "rttbLQModelAccessor.h" #include "rttbDoseBasedModels.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace models { LQModelAccessor::~LQModelAccessor() = default; LQModelAccessor::LQModelAccessor(DoseAccessorPointer dose, BioModelParamType alpha, BioModelParamType beta, unsigned int nFractions, double doseScaling) : _dose(dose), _alpha(alpha), _beta(beta), _nFractions(nFractions), _alphaMap(nullptr), _betaMap(nullptr), _doseScaling(doseScaling), _withAlphaBetaMaps(false) { if (_dose == nullptr) { throw core::InvalidDoseException("Dose is nullptr"); } if (_doseScaling < 0) { throw core::InvalidParameterException("Dose Scaling must be >0"); } assembleGeometricInfo(); } LQModelAccessor::LQModelAccessor(DoseAccessorPointer dose, DoseAccessorPointer alphaMap, DoseAccessorPointer betaMap, unsigned int nFractions, double doseScaling) :_dose(dose), _alpha(-1.), _beta(-1.), _nFractions(nFractions), _alphaMap(betaMap), _betaMap(alphaMap), _doseScaling(doseScaling), _withAlphaBetaMaps(true) { if (_dose == nullptr || _alphaMap == nullptr || _betaMap == nullptr) { throw core::InvalidDoseException("Dose or alphaMap or betaMap is nullptr"); } if (_doseScaling < 0) { throw core::InvalidParameterException("Dose Scaling must be >0"); } assembleGeometricInfo(); } GenericValueType LQModelAccessor::getValueAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex3D; if (_geoInfo.convert(aID, aVoxelGridIndex3D)) { return getValueAt(aVoxelGridIndex3D); } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); } } GenericValueType LQModelAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { if (_withAlphaBetaMaps){ return calcLQ(_dose->getValueAt(aIndex) * _doseScaling, _alphaMap->getValueAt(aIndex), _betaMap->getValueAt(aIndex), _nFractions); } else { return calcLQ(_dose->getValueAt(aIndex) * _doseScaling, _alpha, _beta, _nFractions); } } void LQModelAccessor::assembleGeometricInfo() { _geoInfo = _dose->getGeometricInfo(); } } } diff --git a/code/models/rttbLQModelAccessor.h b/code/models/rttbLQModelAccessor.h index 3c6ab46..c97ba08 100644 --- a/code/models/rttbLQModelAccessor.h +++ b/code/models/rttbLQModelAccessor.h @@ -1,92 +1,87 @@ // ----------------------------------------------------------------------- // 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 __LQ_MODEL_ACCESSOR_H #define __LQ_MODEL_ACCESSOR_H #include "rttbAccessorWithGeoInfoBase.h" #include "rttbDoseAccessorInterface.h" #include "rttbBaseTypeModels.h" namespace rttb { namespace models { /*! @class LQModelAccessor @brief This class gives access to the LQ Model information in an image */ class LQModelAccessor: public core::AccessorWithGeoInfoBase { public: using DoseAccessorPointer = core::DoseAccessorInterface::Pointer; private: DoseAccessorPointer _dose; BioModelParamType _alpha; BioModelParamType _beta; unsigned int _nFractions; DoseAccessorPointer _alphaMap; DoseAccessorPointer _betaMap; double _doseScaling; bool _withAlphaBetaMaps; IDType _bioModelUID; LQModelAccessor() = delete; /*! @brief get all required data from the dose geometric info */ void assembleGeometricInfo() override; public: ~LQModelAccessor() override; /*! @brief Constructor. @pre dose must be a valid instance (and != nullptr) @exception InvalidDoseException if _dose is nullptr */ LQModelAccessor(DoseAccessorPointer dose, BioModelParamType alpha, BioModelParamType beta, unsigned int nFractions=1, double doseScaling = 1.0); /*! @brief Constructor. @pre dose must be a valid instance (and != nullptr) @exception InvalidDoseException if dose is nullptr, if alphaMap is nullptr or if betaMap is nullptr */ LQModelAccessor(DoseAccessorPointer dose, DoseAccessorPointer alphaMap, DoseAccessorPointer betaMap, unsigned int nFractions = 1, double doseScaling = 1.0); /*! @brief returns the LQ Model value for an id */ GenericValueType getValueAt(const VoxelGridID aID) const override; /*! @brief returns the LQ Model value for an index */ GenericValueType getValueAt(const VoxelGridIndex3D& aIndex) const override; const IDType getUID() const override { return _bioModelUID; }; }; } } #endif diff --git a/code/models/rttbNTCPLKBModel.cpp b/code/models/rttbNTCPLKBModel.cpp index 34fd4ad..b7cc495 100644 --- a/code/models/rttbNTCPLKBModel.cpp +++ b/code/models/rttbNTCPLKBModel.cpp @@ -1,181 +1,175 @@ // ----------------------------------------------------------------------- // 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) -*/ + #define _USE_MATH_DEFINES #include #include #include #include #include #include "rttbIntegration.h" #include "rttbNTCPLKBModel.h" #include "rttbDvhBasedModels.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" - namespace rttb { namespace models { NTCPLKBModel::NTCPLKBModel() : NTCPModel(){ _name = "NTCPLKBModel"; fillParameterMap(); } NTCPLKBModel::NTCPLKBModel(DVHPointer aDvh, BioModelParamType aD50, BioModelParamType aM, BioModelParamType aA): NTCPModel(aDvh, aD50), _m(aM), _a(aA) { _name = "NTCPLKBModel"; fillParameterMap(); } void NTCPLKBModel::setA(const BioModelParamType aA) { _a = aA; } const BioModelParamType NTCPLKBModel::getA() { return _a; } void NTCPLKBModel::setM(const BioModelParamType aM) { _m = aM; } const BioModelParamType NTCPLKBModel::getM() { return _m; } void NTCPLKBModel::setParameterVector(const ParamVectorType& aParameterVector) { if (aParameterVector.size() != 3) { throw core::InvalidParameterException("Parameter invalid: aParameterVector.size must be 3! "); } else { _d50 = aParameterVector.at(0); _m = aParameterVector.at(1); _a = aParameterVector.at(2); } } void NTCPLKBModel::setParameterByID(const int aParamId, const BioModelParamType aValue) { if (aParamId == 0) { _d50 = aValue; } else if (aParamId == 1) { _m = aValue; } else if (aParamId == 2) { _a = aValue; } else { throw core::InvalidParameterException("Parameter invalid: aParamID must be 0(for d50) or 1(for m) or 2(for a)! "); } } const int NTCPLKBModel::getParameterID(const std::string& aParamName) const { if (aParamName == "d50") { return 0; } else if (aParamName == "m") { return 1; } else if (aParamName == "a") { return 2; } else { rttbExceptionMacro(core::InvalidParameterException, << "Parameter name " << aParamName << " invalid: it should be d50 or m or a!"); } } std::map NTCPLKBModel::getParameterMap() const{ return parameterMap; } void NTCPLKBModel::fillParameterMap(){ parameterMap["d50"] = getD50(); parameterMap["m"] = getM(); parameterMap["a"] = getA(); } std::string NTCPLKBModel::getModelType() const{ return _name; } BioModelValueType NTCPLKBModel::calcModel(const double doseFactor) { if (_a == 0) { throw core::InvalidParameterException("_a must not be zero"); } if (_m == 0) { throw core::InvalidParameterException("_m must not be zero"); } core::DVH variantDVH = core::DVH(_dvh->getDataDifferential(), (DoseTypeGy)(_dvh->getDeltaD() * doseFactor), _dvh->getDeltaV(), "temporary", "temporary"); auto spDVH = boost::make_shared(variantDVH); double eud = getEUD(spDVH, this->_a); //_m must not be zero double t = (eud - this->_d50) / (this->_m * this->_d50); double value = 1 / pow(2 * M_PI, 0.5); double result = integrateLKB(t); if (result != -100) { value *= result; return value; } else { return false; } } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbNTCPLKBModel.h b/code/models/rttbNTCPLKBModel.h index 23c6c71..56b8690 100644 --- a/code/models/rttbNTCPLKBModel.h +++ b/code/models/rttbNTCPLKBModel.h @@ -1,104 +1,97 @@ // ----------------------------------------------------------------------- // 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 __NTCP_LKB_MODEL_H #define __NTCP_LKB_MODEL_H #include #include "rttbNTCPModel.h" #include "rttbBaseTypeModels.h" - - namespace rttb { namespace models { /*! @class NTCPLKBModel @brief This class represents a NTCP(Normal Tissue Complication Probability) LKB model (Lyman 1985, Kutcher and Burman 1989) @see NTCPModel */ class NTCPLKBModel: public NTCPModel { public: using ParamVectorType = NTCPModel::ParamVectorType; private: /*! The steepness of the dose-response curve. Must not be zero on model evaluation. */ BioModelParamType _m{0}; /*! Tumor or normal tissue-specific parameter that describes the dose-volume effect, e.g. -10 for prostate (Wu 2002). Must not be zero on model evaluation, because EUD calculation will fail. */ BioModelParamType _a{0}; protected: /*! @brief Calculate the model value * @param doseFactor: scaling factor for the dose. The model calculation will use the dvh with each di=old di*doseFactor. * @throw if either _a or _m is zero for the model calculation */ BioModelValueType calcModel(const double doseFactor = 1) override; public: NTCPLKBModel(); NTCPLKBModel(core::DVH::Pointer aDvh, BioModelParamType aD50, BioModelParamType aM, BioModelParamType aA); void setM(const BioModelParamType aM); const BioModelParamType getM(); void setA(const BioModelParamType aA); const BioModelParamType getA(); /*! @brief Set parameter with ID. "d50":0,"m":1,"a":2 @exception InvalidParameterException Thrown if aParamId is not 0 or 1 or 2. */ void setParameterByID(const int aParamId, const BioModelParamType aValue) override; /*! @brief Set parameter vector, where index of vector is the parameter ID. "d50":0,"m":1,"a":2 @exception InvalidParameterException Thrown if aParamterVector.size()!=3. */ void setParameterVector(const ParamVectorType& aParameterVector) override; /*! @brief Get parameter ID. "d50":0,"m":1,"a":2 @return 0 for "d50", 1 for "m", 2 for "a" @exception InvalidParameterException Thrown if aParamName is not d50 or m or a. */ const int getParameterID(const std::string& aParamName) const override; std::map getParameterMap() const override; void fillParameterMap() override; std::string getModelType() const override; }; } } #endif diff --git a/code/models/rttbNTCPModel.h b/code/models/rttbNTCPModel.h index 6e509d8..8648bfa 100644 --- a/code/models/rttbNTCPModel.h +++ b/code/models/rttbNTCPModel.h @@ -1,66 +1,57 @@ // ----------------------------------------------------------------------- // 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 __NTCP_MODEL_H #define __NTCP_MODEL_H - #include "rttbBioModel.h" #include "rttbBaseTypeModels.h" - - namespace rttb { namespace models { /*! @class NTCPModel @brief This is the interface class for NTCP(Normal Tissue Complication Probability) models */ class NTCPModel: public BioModel { public: using ParamVectorType = BioModel::ParamVectorType; protected: BioModelParamType _d50{0}; public: NTCPModel(): BioModel() {} NTCPModel(const BioModelParamType aD50): BioModel(), _d50(aD50) {} NTCPModel(core::DVH::Pointer aDvh, const BioModelParamType aD50): BioModel(aDvh), _d50(aD50) {} const BioModelParamType getD50() { return _d50; } void setD50(const BioModelParamType aD50) { _d50 = aD50; } }; }//end namespace models }//end namespace rttb #endif diff --git a/code/models/rttbNTCPRSModel.cpp b/code/models/rttbNTCPRSModel.cpp index 372ed39..e2b549e 100644 --- a/code/models/rttbNTCPRSModel.cpp +++ b/code/models/rttbNTCPRSModel.cpp @@ -1,173 +1,168 @@ // ----------------------------------------------------------------------- // 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) -*/ + #define _USE_MATH_DEFINES #include #include #include "rttbNTCPRSModel.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace models { NTCPRSModel::NTCPRSModel() : NTCPModel() { _name = "NTCPRSModel"; fillParameterMap(); } NTCPRSModel::NTCPRSModel(DVHPointer aDvh, BioModelParamType aD50, BioModelParamType aGamma, BioModelParamType aS): NTCPModel(aDvh, aD50), _gamma(aGamma), _s(aS) { _name = "NTCPRSModel"; fillParameterMap(); } void NTCPRSModel::setGamma(const BioModelParamType aGamma) { _gamma = aGamma; } const BioModelParamType NTCPRSModel::getGamma() { return _gamma; } void NTCPRSModel::setS(const BioModelParamType aS) { _s = aS; } const BioModelParamType NTCPRSModel::getS() { return _s; } void NTCPRSModel::setParameterVector(const ParamVectorType& aParameterVector) { if (aParameterVector.size() != 3) { throw core::InvalidParameterException("Parameter invalid: aParameterVector.size must be 3! "); } else { _d50 = aParameterVector.at(0); _gamma = aParameterVector.at(1); _s = aParameterVector.at(2); } } void NTCPRSModel::setParameterByID(const int aParamId, const BioModelParamType aValue) { if (aParamId == 0) { _d50 = aValue; } else if (aParamId == 1) { _gamma = aValue; } else if (aParamId == 2) { _s = aValue; } else { throw core::InvalidParameterException("Parameter invalid: aParamID must be 0(for d50) or 1(for gamma) or 2(for s)! "); } } const int NTCPRSModel::getParameterID(const std::string& aParamName) const { if (aParamName == "d50") { return 0; } else if (aParamName == "gamma") { return 1; } else if (aParamName == "s") { return 2; } else { rttbExceptionMacro(core::InvalidParameterException, << "Parameter name " << aParamName << " invalid: it should be d50 or gamma or s!"); } } std::map NTCPRSModel::getParameterMap() const{ return parameterMap; } void NTCPRSModel::fillParameterMap() { parameterMap["d50"] = getD50(); parameterMap["gamma"] = getGamma(); parameterMap["s"] = getS(); } std::string NTCPRSModel::getModelType() const{ return _name; } const double NTCPRSModel::poissonModel(const double dose) { //_d50 must not be zero return pow(2, -exp(M_E * this->_gamma * (1 - dose / this->_d50))); } BioModelValueType NTCPRSModel::calcModel(double doseFactor) { if (_d50 == 0) { throw core::InvalidParameterException("d50 must not be zero"); } if (_s == 0) { throw core::InvalidParameterException("s must not be zero"); } std::deque dataDifferential = this->_dvh->getDataDifferential(); double ntcp = 1; for (GridIndexType i = 0; i < dataDifferential.size(); i++) { double pd = pow(this->poissonModel(i * this->_dvh->getDeltaD() * doseFactor), this->_s); double vi = dataDifferential[i] / this->_dvh->getNumberOfVoxels(); ntcp *= pow((1 - pd), vi); } //_s must not be zero return (BioModelValueType)(pow((1 - ntcp), 1 / this->_s)); } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbNTCPRSModel.h b/code/models/rttbNTCPRSModel.h index 9aba15b..387871b 100644 --- a/code/models/rttbNTCPRSModel.h +++ b/code/models/rttbNTCPRSModel.h @@ -1,111 +1,104 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) -*/ #ifndef __NTCP_RS_MODEL_H #define __NTCP_RS_MODEL_H #include #include "rttbNTCPModel.h" #include "rttbBaseTypeModels.h" - namespace rttb { namespace models { /*! @class NTCPRSModel @brief This class represents a NTCP(Normal Tissue Complication Probability) relative seriality model (Kaellman 1992) @see NTCPModel */ class NTCPRSModel: public NTCPModel { public: using ParamVectorType = NTCPModel::ParamVectorType; using DVHPointer = NTCPModel::DVHPointer; private: /*! _gamma The normalised dose-response gradient, values between 1.7 and 2.0 are typical for human tumours. (Kaellman 1992) */ BioModelParamType _gamma{0}; /*! _s The relative seriality factor, e.g. s=3.4 for the esophagus (highly serial structure) and s=0.0061 for the lung(highly parallel structure). Must not be zero on model evaluation. */ BioModelParamType _s{0}; const double poissonModel(const double dose); protected: /*! @brief Calculate the model value @param doseFactor scaling factor for the dose. The model calculation will use the dvh with each di=old di*doseFactor. @throw if either _s or _d50 is zero for the model calculation. */ BioModelValueType calcModel(const double doseFactor = 1) override; public: NTCPRSModel(); /*!@brief Constructor initializing all member variables with given parameters. */ NTCPRSModel(DVHPointer aDvh, BioModelParamType aD50, BioModelParamType aGamma, BioModelParamType aS); void setGamma(const BioModelParamType aGamma); const BioModelParamType getGamma(); void setS(const BioModelParamType aS); const BioModelParamType getS(); /*! @brief Set parameter with ID. "d50":0,"gamma":1,"s":2 @exception InvalidParameterException Thrown if aParamId is not 0 or 1 or 2. */ void setParameterByID(const int aParamId, const BioModelParamType aValue) override; /*! @brief Set parameter vector, where index of vector is the parameter Id. "d50":0,"gamma":1,"s":2 @exception InvalidParameterException Thrown if aParamterVector.size()!=3. */ void setParameterVector(const ParamVectorType& aParameterVector) override; /*! @brief Get parameter ID. "d50":0,"gamma":1,"s":2 @return 0 for "d50", 1 for "gamma", 2 for "s" @exception InvalidParameterException Thrown if aParamName is not d50 or gamma or s. */ const int getParameterID(const std::string& aParamName) const override; std::map getParameterMap() const override; void fillParameterMap() override; std::string getModelType() const override; }; } } #endif diff --git a/code/models/rttbTCPLQModel.cpp b/code/models/rttbTCPLQModel.cpp index ac82e94..9651baa 100644 --- a/code/models/rttbTCPLQModel.cpp +++ b/code/models/rttbTCPLQModel.cpp @@ -1,311 +1,304 @@ // ----------------------------------------------------------------------- // 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) -*/ #define _USE_MATH_DEFINES #include #include #include #include #include #include "rttbTCPLQModel.h" #include "rttbDvhBasedModels.h" #include "rttbIntegration.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" - namespace rttb { namespace models { TCPLQModel::TCPLQModel(): TCPModel() { _name = "TCPLQModel"; fillParameterMap(); } 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) { _name = "TCPLQModel"; fillParameterMap(); } 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) { fillParameterMap(); _name = "TCPLQModel"; } 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::getAlphaBeta() { 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; 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"); auto spDVH = boost::make_shared(variantDVH); BioModelValueType value = 0; if (_alphaVariance == 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!"); } 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) { 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 " << aParamName << " invalid: it should be alphaMean or alphaVariance or alpha_beta or rho!"); } } std::map TCPLQModel::getParameterMap() const{ return parameterMap; } void TCPLQModel::fillParameterMap(){ parameterMap["numberOfFraction"] = getNumberOfFractions(); parameterMap["alphaMean"] = getAlphaMean(); parameterMap["alphaVariance"] = getAlphaVariance(); parameterMap["alpha_beta"] = getAlphaBeta(); parameterMap["rho"] = getRho(); } std::string TCPLQModel::getModelType() const{ return _name; } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbTCPLQModel.h b/code/models/rttbTCPLQModel.h index 3a96812..9b7eba2 100644 --- a/code/models/rttbTCPLQModel.h +++ b/code/models/rttbTCPLQModel.h @@ -1,169 +1,163 @@ // ----------------------------------------------------------------------- // 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_LQ_MODEL_H #define __TCP_LQ_MODEL_H #include #include #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: using ParamVectorType = TCPModel::ParamVectorType; using DVHPointer = core::DVH::Pointer; 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{0}; BioModelParamType _alphaVariance{0}; BioModelParamType _alpha_beta{0}; /*! Roh is the initial clonogenic cell density */ BioModelParamType _rho{0}; /*! @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) override; 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 getAlphaBeta(); 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. */ void setParameterByID(const int aParamId, const BioModelParamType aValue) override; /*! @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. */ void setParameterVector(const ParamVectorType& aParameterVector) override; /*! @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. */ const int getParameterID(const std::string& aParamName) const override; std::map getParameterMap() const override; void fillParameterMap() override; std::string getModelType() const override; }; }//end algorithms }//end rttb #endif diff --git a/code/models/rttbTCPModel.cpp b/code/models/rttbTCPModel.cpp index 6eb54d3..ba528ac 100644 --- a/code/models/rttbTCPModel.cpp +++ b/code/models/rttbTCPModel.cpp @@ -1,42 +1,36 @@ // ----------------------------------------------------------------------- // 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 "rttbTCPModel.h" namespace rttb { namespace models { void TCPModel::setNumberOfFractions(const int aNumberOfFractions) { _numberOfFractions = aNumberOfFractions; } const int TCPModel::getNumberOfFractions() { return _numberOfFractions; } } } diff --git a/code/models/rttbTCPModel.h b/code/models/rttbTCPModel.h index b82577d..652903e 100644 --- a/code/models/rttbTCPModel.h +++ b/code/models/rttbTCPModel.h @@ -1,62 +1,56 @@ // ----------------------------------------------------------------------- // 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 "rttbBioModel.h" namespace rttb { namespace models { /*! @class TCPModel @brief This is the interface class for TCP(Tumor Control Probability) models */ class TCPModel: public BioModel { public: using ParamVectorType = BioModel::ParamVectorType; protected: int _numberOfFractions{0}; public: TCPModel(): BioModel() {}; TCPModel(int aNum): BioModel(), _numberOfFractions(aNum) {}; TCPModel(core::DVH::Pointer 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/models/BioModelScatterPlotTest.cpp b/testing/models/BioModelScatterPlotTest.cpp index f620327..2b4c4fd 100644 --- a/testing/models/BioModelScatterPlotTest.cpp +++ b/testing/models/BioModelScatterPlotTest.cpp @@ -1,239 +1,233 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) -*/ #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbBaseTypeModels.h" #include "rttbBioModel.h" #include "rttbBioModelScatterPlots.h" #include "../core/DummyDVHGenerator.h" #include "DummyModel.h" #include "rttbIntegration.h" #include #include #include namespace rttb { namespace testing { typedef models::ScatterPlotType ScatterPlotType; /*! @brief BioModelScatterPlotTest. Test the scatterplot methods. 1) test if parameters are set correctly 2) test if scatterData contain each model value exactly once */ int BioModelScatterPlotTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; //generate artificial DVH and corresponding statistical values DummyDVHGenerator dummyDVH; const IDType structureID = "myStructure"; const IDType doseID = "myDose"; core::DVH::Pointer dvhPtr = boost::make_shared(dummyDVH.generateDVH(structureID, doseID, 0, 2000)); //test Dummy Model models::BioModelParamType param1 = 0.35; models::BioModelParamType param2 = 0.023333333333333; models::BioModelParamType param3 = 10000000; //model values will be dose/normalisationDose DoseTypeGy normalisationDose = 10; rttb::models::DummyModel model(dvhPtr); model.resetCounters(); //default number of points is 100 CHECK_NO_THROW(models::getScatterPlotVary1Parameter(model, 0, param1, 0.5, normalisationDose)); model.resetCounters(); ScatterPlotType scatter = models::getScatterPlotVary1Parameter(model, 0, param1, 0.5, normalisationDose); //only 1st parameter should gave been changed CHECK_EQUAL(100, model.getSetParam1Count()); CHECK_EQUAL(0, model.getSetParam2Count()); CHECK_EQUAL(0, model.getSetParam3Count()); CHECK_EQUAL(100, model.getCalcCount()); // all model values should match the corresponding dose/normalizationDose ScatterPlotType::iterator iter; for (iter = scatter.begin(); iter != scatter.end(); ++iter) { double currentKey = iter->first; if ((currentKey / normalisationDose) != iter->second.first) { CHECK_EQUAL((currentKey / normalisationDose), iter->second.first); } } model.resetCounters(); //number of points is 50 CHECK_NO_THROW(models::getScatterPlotVary1Parameter(model, 1, param2, 0.25, normalisationDose, 50)); model.resetCounters(); scatter = models::getScatterPlotVary1Parameter(model, 1, param2, 0.25, normalisationDose, 50); //only 1st parameter should gave been changed CHECK_EQUAL(0, model.getSetParam1Count()); CHECK_EQUAL(50, model.getSetParam2Count()); CHECK_EQUAL(0, model.getSetParam3Count()); CHECK_EQUAL(50, model.getCalcCount()); // all model values should match the corresponding dose/normalizationDose for (iter = scatter.begin(); iter != scatter.end(); ++iter) { double currentKey = iter->first; if ((currentKey / normalisationDose) != iter->second.first) { CHECK_EQUAL((currentKey / normalisationDose), iter->second.first); } } model.resetCounters(); //number of points is 150 CHECK_NO_THROW(models::getScatterPlotVary1Parameter(model, 2, param3, 0.75, normalisationDose, 150)); model.resetCounters(); scatter = models::getScatterPlotVary1Parameter(model, 2, param3, 0.75, normalisationDose, 150); //only 1st parameter should gave been changed CHECK_EQUAL(0, model.getSetParam1Count()); CHECK_EQUAL(0, model.getSetParam2Count()); CHECK_EQUAL(150, model.getSetParam3Count()); CHECK_EQUAL(150, model.getCalcCount()); // all model values should match the corresponding dose/normalizationDose for (iter = scatter.begin(); iter != scatter.end(); ++iter) { double currentKey = iter->first; if ((currentKey / normalisationDose) != iter->second.first) { CHECK_EQUAL((currentKey / normalisationDose), iter->second.first); } } model.resetCounters(); std::vector paramIDVec; models::BioModel::ParamVectorType meanVec; models::BioModel::ParamVectorType varianceVec; paramIDVec.push_back(0); meanVec.push_back(1); varianceVec.push_back(0.6); paramIDVec.push_back(1); meanVec.push_back(10); varianceVec.push_back(0.5); paramIDVec.push_back(2); meanVec.push_back(100); varianceVec.push_back(0.4); //number of points is 50 CHECK_NO_THROW(models::getScatterPlotVaryParameters(model, paramIDVec, meanVec, varianceVec, normalisationDose, 50)); model.resetCounters(); scatter = models::getScatterPlotVaryParameters(model, paramIDVec, meanVec, varianceVec, normalisationDose, 50); //only 1st parameter should gave been changed CHECK_EQUAL(50, model.getSetParam1Count()); CHECK_EQUAL(50, model.getSetParam2Count()); CHECK_EQUAL(50, model.getSetParam3Count()); CHECK_EQUAL(50, model.getCalcCount()); // all model values should match the corresponding dose/normalizationDose for (iter = scatter.begin(); iter != scatter.end(); ++iter) { double currentKey = iter->first; if ((currentKey / normalisationDose) != iter->second.first) { CHECK_EQUAL((currentKey / normalisationDose), iter->second.first); } } model.resetCounters(); //test 2000 points CHECK_NO_THROW(models::getScatterPlotVaryParameters(model, paramIDVec, meanVec, varianceVec, normalisationDose, 10000)); //test iterativeIntegration /*std::string bedFileName="d:/Temp/bed.txt"; std::ifstream dvh_ifstr(bedFileName,std::ios::in); std::vector volumeV2; std::vector bedV2; //std::cout << "iterative integration"<< std::endl; if ( !dvh_ifstr.is_open() ) { std::cerr<< "DVH file name invalid: could not open the dvh file!"<> volume; volumeV2.push_back(volume); //std::cout << "("<< volume << ","; std::string lineSub2 = line.substr(found2+1,found3-found2-1); //std::cout << lineSub2 << std::endl; std::stringstream ss2(lineSub2); double bed; ss2 >> bed; bedV2.push_back(bed); //std::cout << bed << ");"; } } } struct rttb::models::TcpParams params={0.302101,0.08,10000000,volumeV2,bedV2}; double result = rttb::models::integrateTCP(0, params);*/ RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/models/BioModelTest.cpp b/testing/models/BioModelTest.cpp index a4a6a36..6f4e9be 100644 --- a/testing/models/BioModelTest.cpp +++ b/testing/models/BioModelTest.cpp @@ -1,306 +1,297 @@ // ----------------------------------------------------------------------- // 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::Pointer 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 = 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.getAlphaBeta()); 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.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); CHECK_EQUAL(0, tcplq.getValue()); tcplq = rttb::models::TCPLQModel(); CHECK_EQUAL(0, tcplq.getAlphaMean()); CHECK_EQUAL(0, tcplq.getAlphaBeta()); 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.getAlphaBeta()); 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.getAlphaBeta()); 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.getAlphaBeta()); 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.getAlphaBeta()); 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/DummyModel.cpp b/testing/models/DummyModel.cpp index d3848f3..a5a025b 100644 --- a/testing/models/DummyModel.cpp +++ b/testing/models/DummyModel.cpp @@ -1,126 +1,121 @@ // ----------------------------------------------------------------------- // 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 "DummyModel.h" namespace rttb { namespace models { DummyModel::DummyModel(DVHPointer aDvh): BioModel(aDvh) { _calcCount = 0; _setParam1Count = 0; _setParam2Count = 0; _setParam3Count = 0; _name = "DummyModel"; } BioModelValueType DummyModel::calcModel(const double doseFactor) { _calcCount++; _value = doseFactor; return _value; } void DummyModel::setParameterVector(const ParamVectorType& aParameterVector) { if (aParameterVector.size() != 3) { std::cerr << "aParameterVector.size must be 3! Nothing will be done." << std::endl; } else { _param1 = aParameterVector.at(0); _setParam1Count++; _param2 = aParameterVector.at(1); _setParam2Count++; _param3 = aParameterVector.at(2); _setParam3Count++; } } void DummyModel::setParameterByID(const int aParamId, const BioModelParamType aValue) { if (aParamId == 0) { _param1 = aValue; _setParam1Count++; } else if (aParamId == 1) { _param2 = aValue; _setParam2Count++; } else if (aParamId == 2) { _param3 = aValue; _setParam3Count++; } else { std::cerr << "aParamID must be 0(DummyParam1) or 1(DummyParam2) or 2(DummyParam3)! Nothing will be done." << std::endl; } } const int DummyModel::getParameterID(const std::string& aParamName) const { if (aParamName == "DummyParam1") { return 0; } else if (aParamName == "DummyParam2") { return 1; } else if (aParamName == "DummyParam3") { return 2; } else { std::cerr << "There is no parameter with the name " << aParamName << "!" << std::endl; return -1; } } void DummyModel::resetCounters() { _calcCount = 0; _setParam1Count = 0; _setParam2Count = 0; _setParam3Count = 0; } std::map DummyModel::getParameterMap() const { return parameterMap; } void DummyModel::fillParameterMap(){ parameterMap["Dummy"] = 2; } std::string DummyModel::getModelType() const{ return _name; } } } \ No newline at end of file diff --git a/testing/models/DummyModel.h b/testing/models/DummyModel.h index 742ba1c..2499831 100644 --- a/testing/models/DummyModel.h +++ b/testing/models/DummyModel.h @@ -1,97 +1,92 @@ // ----------------------------------------------------------------------- // 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 __DUMMY_MODEL_H #define __DUMMY_MODEL_H #include "rttbBaseType.h" #include "rttbBaseTypeModels.h" #include "rttbBioModel.h" namespace rttb { namespace models { /*! @class DummyModel @brief generates a dummy model object used for unit tests */ class DummyModel: public BioModel { private: BioModelParamType _param1; BioModelParamType _param2; BioModelParamType _param3; int _calcCount; int _setParam1Count; int _setParam2Count; int _setParam3Count; protected: /*! @brief Calculate the model value @param doseFactor scaling factor for the dose. The model calculation will use the dvh with each di=old di*doseFactor. */ BioModelValueType calcModel(const double doseFactor = 1); std::map getParameterMap() const; void fillParameterMap() override; std::string getModelType() const; public: /*!@constructor initializes DVH pointer */ DummyModel(DVHPointer aDvh); /*! @brief Set parameter vector, where index of vector is the parameter ID. The length of the vector has to be 3. */ void setParameterVector(const ParamVectorType& aParameterVector); void setParameterByID(const int aParamId, const BioModelParamType aValue); /*! @brief Get parameter by ID. @return Return -1 if ID is not found. */ const int getParameterID(const std::string& aParamName) const; /*! return counting values */ const int getSetParam1Count() const { return _setParam1Count; }; const int getSetParam2Count() const { return _setParam2Count; }; const int getSetParam3Count() const { return _setParam3Count; }; const int getCalcCount() const { return _calcCount; }; void resetCounters(); }; } } #endif \ No newline at end of file diff --git a/testing/models/DvhBasedModelsTest.cpp b/testing/models/DvhBasedModelsTest.cpp index ee4ea26..318bc4d 100644 --- a/testing/models/DvhBasedModelsTest.cpp +++ b/testing/models/DvhBasedModelsTest.cpp @@ -1,110 +1,105 @@ // ----------------------------------------------------------------------- // 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 "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::Pointer 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/LQModelAccessorTest.cpp b/testing/models/LQModelAccessorTest.cpp index 58627e7..c581f08 100644 --- a/testing/models/LQModelAccessorTest.cpp +++ b/testing/models/LQModelAccessorTest.cpp @@ -1,155 +1,150 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) -*/ + #include "litCheckMacros.h" #include #include "rttbLQModelAccessor.h" #include "rttbGenericDoseIterator.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace testing { /*! @brief LQModelAccessorTest. 1) Test constructors 2) Test getGeometricInfo() 3) Test getValueAt() */ int LQModelAccessorTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string RTDOSE_FILENAME_CONSTANT_TWO; std::string RTDOSE_FILENAME_INCREASE_X; if (argc > 2) { RTDOSE_FILENAME_CONSTANT_TWO = argv[1]; RTDOSE_FILENAME_INCREASE_X = argv[2]; } else { std::cout << "at least two parameters for LQModelAccessorTest are expected" << std::endl; return -1; } typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1( RTDOSE_FILENAME_CONSTANT_TWO.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator2( RTDOSE_FILENAME_INCREASE_X.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); core::GeometricInfo doseAccessor1GeometricInfo = doseAccessor1->getGeometricInfo(); core::GeometricInfo doseAccessor2GeometricInfo = doseAccessor2->getGeometricInfo(); DoseAccessorPointer doseAccessorNull; core::AccessorInterface::Pointer LQWithConstantDose, LQWithConstantDoseDoseScalingTwo, LQWithConstantNFractionsTwo, LQWithIncreaseXDose, LQWithConstantDoseAndMap; //1) test constructor CHECK_THROW_EXPLICIT(models::LQModelAccessor(doseAccessorNull, 0, 0), core::InvalidDoseException); CHECK_THROW_EXPLICIT(models::LQModelAccessor(doseAccessor1, 0.2, 0.02, 1, -1), core::InvalidParameterException); CHECK_THROW_EXPLICIT(models::LQModelAccessor(doseAccessorNull, doseAccessorNull, doseAccessorNull), core::InvalidDoseException); CHECK_THROW_EXPLICIT(models::LQModelAccessor(doseAccessor1, doseAccessor1, doseAccessorNull), core::InvalidDoseException); CHECK_THROW_EXPLICIT(models::LQModelAccessor(doseAccessor1, doseAccessorNull, doseAccessor1), core::InvalidDoseException); CHECK_THROW_EXPLICIT(models::LQModelAccessor(doseAccessorNull, doseAccessor1, doseAccessor1), core::InvalidDoseException); CHECK_NO_THROW(LQWithConstantDose = boost::make_shared(doseAccessor1, 0.2, 0.02)); CHECK_NO_THROW(LQWithConstantNFractionsTwo = boost::make_shared(doseAccessor1, 0.2, 0.02, 2)); CHECK_NO_THROW(LQWithConstantDoseDoseScalingTwo = boost::make_shared (doseAccessor1, 0.2, 0.02, 1, 2.0)); CHECK_NO_THROW(LQWithIncreaseXDose = boost::make_shared(doseAccessor2, 0.3, 0.01)); CHECK_NO_THROW(LQWithConstantDoseAndMap = boost::make_shared(doseAccessor1, doseAccessor2, doseAccessor2)); //2) Test getGeometricInfo() CHECK_EQUAL(LQWithConstantDose->getGeometricInfo(), doseAccessor1GeometricInfo); CHECK_EQUAL(LQWithIncreaseXDose->getGeometricInfo(), doseAccessor2GeometricInfo); //3) Test getValueAt() models::BioModelParamType expectedLQWithDoseTwo = exp(-(0.2 * 2 + (0.02 * 2 * 2))); CHECK_EQUAL(LQWithConstantDose->getValueAt(0), expectedLQWithDoseTwo); CHECK_EQUAL(LQWithConstantDose->getValueAt(LQWithConstantDose->getGridSize() - 1), expectedLQWithDoseTwo); CHECK_EQUAL(LQWithConstantDose->getValueAt(VoxelGridIndex3D(1, 2, 6)), expectedLQWithDoseTwo); CHECK_EQUAL(LQWithConstantDose->getValueAt(VoxelGridIndex3D(65, 40, 60)), expectedLQWithDoseTwo); models::BioModelParamType expectedLQWithDoseTwoNFractionsTwo = exp(-(0.2 * 2 + (0.02 * 2 * 2/2))); CHECK_EQUAL(LQWithConstantNFractionsTwo->getValueAt(0), expectedLQWithDoseTwoNFractionsTwo); CHECK_EQUAL(LQWithConstantNFractionsTwo->getValueAt(LQWithConstantNFractionsTwo->getGridSize() - 1), expectedLQWithDoseTwoNFractionsTwo); CHECK_EQUAL(LQWithConstantNFractionsTwo->getValueAt(VoxelGridIndex3D(1, 2, 6)), expectedLQWithDoseTwoNFractionsTwo); CHECK_EQUAL(LQWithConstantNFractionsTwo->getValueAt(VoxelGridIndex3D(65, 40, 60)), expectedLQWithDoseTwoNFractionsTwo); models::BioModelParamType expectedLQWithDoseFour = exp(-(0.2 * 4 + (0.02 * 4 * 4))); CHECK_EQUAL(LQWithConstantDoseDoseScalingTwo->getValueAt(0), expectedLQWithDoseFour); CHECK_EQUAL(LQWithConstantDoseDoseScalingTwo->getValueAt(LQWithConstantDose->getGridSize() - 1), expectedLQWithDoseFour); CHECK_EQUAL(LQWithConstantDoseDoseScalingTwo->getValueAt(VoxelGridIndex3D(1, 2, 6)), expectedLQWithDoseFour); CHECK_EQUAL(LQWithConstantDoseDoseScalingTwo->getValueAt(VoxelGridIndex3D(65, 40, 60)), expectedLQWithDoseFour); models::BioModelParamType expectedLQWithIncreaseX = exp(-(0.3 * 66 * 2.822386e-5 + (0.01 * 66 * 2.822386e-5 * 66 * 2.822386e-5))); CHECK_EQUAL(LQWithIncreaseXDose->getValueAt(0), 1); CHECK_CLOSE(LQWithIncreaseXDose->getValueAt(LQWithIncreaseXDose->getGridSize() - 1), expectedLQWithIncreaseX, errorConstant); expectedLQWithIncreaseX = exp(-(0.3 * 1 * 2.822386e-5 + (0.01 * 1 * 2.822386e-5 * 1 * 2.822386e-5))); CHECK_CLOSE(LQWithIncreaseXDose->getValueAt(VoxelGridIndex3D(1, 2, 6)), expectedLQWithIncreaseX, errorConstant); expectedLQWithIncreaseX = exp(-(0.3 * 45 * 2.822386e-5 + (0.01 * 45 * 2.822386e-5 * 45 * 2.822386e-5))); CHECK_CLOSE(LQWithIncreaseXDose->getValueAt(VoxelGridIndex3D(45, 40, 60)), expectedLQWithIncreaseX, errorConstant); models::BioModelParamType expectedLQWithDoseAndMap = exp(-(66 * 2.822386e-5 * 2 + (66 * 2.822386e-5 * 2 * 2))); CHECK_EQUAL(LQWithConstantDoseAndMap->getValueAt(0), 1); CHECK_CLOSE(LQWithConstantDoseAndMap->getValueAt(LQWithIncreaseXDose->getGridSize() - 1), expectedLQWithDoseAndMap, errorConstant); CHECK_EQUAL(LQWithConstantDoseAndMap->getValueAt(0), LQWithConstantDoseAndMap->getValueAt(VoxelGridIndex3D(0, 0, 0))); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb \ No newline at end of file diff --git a/testing/models/rttbModelsTest.cpp b/testing/models/rttbModelsTest.cpp index ee706c1..f2245e3 100644 --- a/testing/models/rttbModelsTest.cpp +++ b/testing/models/rttbModelsTest.cpp @@ -1,66 +1,59 @@ // ----------------------------------------------------------------------- // 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); LIT_REGISTER_TEST(LQModelAccessorTest); } } } 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; }