diff --git a/code/io/models/rttbModelXMLWriter.cpp b/code/io/models/rttbModelXMLWriter.cpp index 6e966b4..cfe9c18 100644 --- a/code/io/models/rttbModelXMLWriter.cpp +++ b/code/io/models/rttbModelXMLWriter.cpp @@ -1,106 +1,106 @@ // ----------------------------------------------------------------------- // 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: 1328 $ (last changed revision) // @date $Date: 2016-04-22 09:50:01 +0200 (Fr, 22 Apr 2016) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #include "rttbModelXMLWriter.h" /*boost includes*/ #include #include #include #include #include "rttbInvalidParameterException.h" namespace rttb { namespace io { namespace models { - ModelXMLWriter::ModelXMLWriter(std::string& const filename, rttb::models::BioModel& const model) : _filename(filename), _model(model){ + ModelXMLWriter::ModelXMLWriter(const std::string& filename, boost::shared_ptr model) : _filename(filename), _model(model){ } void ModelXMLWriter::setFilename(FileNameString filename){ _filename = filename; } FileNameString ModelXMLWriter::getFilename() const{ return _filename; } - void ModelXMLWriter::setModel(rttb::models::BioModel& model){ + void ModelXMLWriter::setModel(boost::shared_ptr model){ _model = model; } - rttb::models::BioModel& ModelXMLWriter::getModel() const { + boost::shared_ptr ModelXMLWriter::getModel() const { return _model; } void ModelXMLWriter::writeModel(){ boost::property_tree::ptree pt; std::string xmlattrNameTag = ".name"; std::string modelTag = "BioModel"; std::string propertyTag = "property"; std::string configTag = "config"; std::string resultsTag = "results"; std::string valueTag = "value"; boost::property_tree::ptree propertynode; boost::property_tree::ptree confignode; - confignode.put("BioModelType", _model.getModelType()); - confignode.put("StructureID", _model.getDVH()->getStructureID()); - confignode.put("DoseID", _model.getDVH()->getDoseID()); + confignode.put("BioModelType", _model->getModelType()); + confignode.put("StructureID", _model->getDVH()->getStructureID()); + confignode.put("DoseID", _model->getDVH()->getDoseID()); pt.add_child(modelTag + "." + configTag, confignode); - propertynode.put("", _model.getValue()); + propertynode.put("", _model->getValue()); propertynode.put(xmlattrNameTag, valueTag); pt.add_child(modelTag + "."+ resultsTag + "." + propertyTag, propertynode); - std::map parameterMap = _model.getParameterMap(); + std::map parameterMap = _model->getParameterMap(); - for (std::map::const_iterator it = parameterMap.begin(); it != parameterMap.end(); it++){ + for (std::map::const_iterator it = parameterMap.begin(); it != parameterMap.end(); ++it){ propertynode.put("", it->second); propertynode.put(xmlattrNameTag, it->first); pt.add_child(modelTag + "." + resultsTag + "." + propertyTag, propertynode); } try { boost::property_tree::xml_parser::xml_writer_settings settings = boost::property_tree::xml_writer_make_settings('\t', 1); boost::property_tree::xml_parser::write_xml(_filename, pt, std::locale(), settings); } - catch (boost::property_tree::xml_parser_error& const e) + catch (const boost::property_tree::xml_parser_error& e) { std::cout << e.what(); throw core::InvalidParameterException("Write xml failed: xml_parser_error!"); } } } } } diff --git a/code/io/models/rttbModelXMLWriter.h b/code/io/models/rttbModelXMLWriter.h index 6758d23..4f8bc7a 100644 --- a/code/io/models/rttbModelXMLWriter.h +++ b/code/io/models/rttbModelXMLWriter.h @@ -1,56 +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: 1250 $ (last changed revision) // @date $Date: 2016-02-18 15:25:55 +0100 (Do, 18 Feb 2016) $ (last change date) // @author $Author: zhangl $ (last changed by) */ #ifndef __MODELS_XML_WRITER_H #define __MODELS_XML_WRITER_H #include "rttbBioModel.h" namespace rttb { namespace io { namespace models { class ModelXMLWriter{ private: std::string _filename; - rttb::models::BioModel& _model; + boost::shared_ptr _model; public: - ModelXMLWriter(std::string& const filename, rttb::models::BioModel& const model); + ModelXMLWriter(const std::string& filename, boost::shared_ptr model); void setFilename(std::string filename); - FileNameString getFilename() const; + std::string getFilename() const; - void setModel(rttb::models::BioModel& model); - rttb::models::BioModel& getModel() const; + void setModel(boost::shared_ptr model); + boost::shared_ptr getModel() const; void writeModel(); }; } } } #endif diff --git a/code/models/rttbBioModel.h b/code/models/rttbBioModel.h index d08b6ce..9b53b6d 100644 --- a/code/models/rttbBioModel.h +++ b/code/models/rttbBioModel.h @@ -1,107 +1,108 @@ // ----------------------------------------------------------------------- // 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 "rttbBaseType.h" #include "rttbDVH.h" #include "rttbBaseTypeModels.h" #include "RTTBModelsExports.h" #include namespace rttb { namespace models { /*! @class BioModel @brief This is the interface class for biological models */ class RTTBModels_EXPORT BioModel { public: typedef std::vector ParamVectorType; typedef core::DVH::DVHPointer DVHPointer; protected: DVHPointer _dvh; BioModelValueType _value; /*! @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(): _value(0) {}; BioModel(DVHPointer aDvh): _dvh(aDvh), _value(0) {}; /*! @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() = 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/rttbNTCPLKBModel.cpp b/code/models/rttbNTCPLKBModel.cpp index 539cb38..36df48f 100644 --- a/code/models/rttbNTCPLKBModel.cpp +++ b/code/models/rttbNTCPLKBModel.cpp @@ -1,178 +1,183 @@ // ----------------------------------------------------------------------- // 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 #include #define _USE_MATH_DEFINES #include #include "rttbIntegration.h" #include "rttbNTCPLKBModel.h" #include "rttbDvhBasedModels.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace models { NTCPLKBModel::NTCPLKBModel() : NTCPModel(), _m(0), _a(0){ _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(){ + std::map NTCPLKBModel::getParameterMap() const{ + return parameterMap; + } + + void NTCPLKBModel::fillParameterMap(){ parameterMap["d50"] = getD50(); parameterMap["m"] = getM(); parameterMap["a"] = getA(); - return parameterMap; } 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"); boost::shared_ptr 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 b9ba18a..6fcd713 100644 --- a/code/models/rttbNTCPLKBModel.h +++ b/code/models/rttbNTCPLKBModel.h @@ -1,105 +1,107 @@ // ----------------------------------------------------------------------- // 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 #include "rttbBaseType.h" #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: typedef NTCPModel::ParamVectorType ParamVectorType; typedef NTCPModel::DVHPointer DVHPointer; private: /*! The steepness of the dose-response curve. Must not be zero on model evaluation. */ BioModelParamType _m; /*! 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; 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); public: NTCPLKBModel(); NTCPLKBModel(DVHPointer 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. */ virtual void setParameterByID(const int aParamId, const BioModelParamType aValue); /*! @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. */ virtual void setParameterVector(const ParamVectorType& aParameterVector); /*! @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. */ virtual const int getParameterID(const std::string& aParamName) const; - virtual std::map getParameterMap(); + virtual std::map getParameterMap() const; - virtual std::string getModelType() const; + void fillParameterMap(); + + virtual std::string getModelType() const override; }; } } #endif diff --git a/code/models/rttbNTCPRSModel.cpp b/code/models/rttbNTCPRSModel.cpp index 93af17f..9da34d8 100644 --- a/code/models/rttbNTCPRSModel.cpp +++ b/code/models/rttbNTCPRSModel.cpp @@ -1,168 +1,174 @@ // ----------------------------------------------------------------------- // 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 "rttbNTCPRSModel.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace models { NTCPRSModel::NTCPRSModel() : NTCPModel(), _gamma(0), _s(0){ _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(){ + std::map NTCPRSModel::getParameterMap() const{ + return parameterMap; + } + + void NTCPRSModel::fillParameterMap() + { parameterMap["d50"] = getD50(); parameterMap["gamma"] = getGamma(); parameterMap["s"] = getS(); - return parameterMap; } 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 32068eb..491254f 100644 --- a/code/models/rttbNTCPRSModel.h +++ b/code/models/rttbNTCPRSModel.h @@ -1,111 +1,113 @@ // ----------------------------------------------------------------------- // 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 #include "rttbBaseType.h" #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 (Källman 1992) @see NTCPModel */ class NTCPRSModel: public NTCPModel { public: typedef NTCPModel::ParamVectorType ParamVectorType; typedef NTCPModel::DVHPointer DVHPointer; private: /*! _gamma The normalised dose-response gradient, values between 1.7 and 2.0 are typical for human tumours. (Källman 1992) */ BioModelParamType _gamma; /*! _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; 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); 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. */ virtual void setParameterByID(const int aParamId, const BioModelParamType aValue); /*! @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. */ virtual void setParameterVector(const ParamVectorType& aParameterVector); /*! @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. */ virtual const int getParameterID(const std::string& aParamName) const; - virtual std::map getParameterMap(); + virtual std::map getParameterMap() const; - virtual std::string getModelType() const; + void fillParameterMap(); + + virtual std::string getModelType() const override; }; } } #endif diff --git a/code/models/rttbTCPLQModel.cpp b/code/models/rttbTCPLQModel.cpp index eaceb5a..38d3b4f 100644 --- a/code/models/rttbTCPLQModel.cpp +++ b/code/models/rttbTCPLQModel.cpp @@ -1,303 +1,308 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #define _USE_MATH_DEFINES #include #include #include #include #include "rttbTCPLQModel.h" #include "rttbDvhBasedModels.h" #include "rttbBaseTypeModels.h" #include "rttbIntegration.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace models { TCPLQModel::TCPLQModel(): TCPModel(), _alphaMean(0), _alphaVariance(0), _alpha_beta(0), _rho(0) { _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) {} 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"); boost::shared_ptr 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(){ + 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(); - return parameterMap; } 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 5fc6b2d..a6946f2 100644 --- a/code/models/rttbTCPLQModel.h +++ b/code/models/rttbTCPLQModel.h @@ -1,169 +1,171 @@ // ----------------------------------------------------------------------- // 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 #include "rttbBaseType.h" #include "rttbTCPModel.h" #include "rttbBaseTypeModels.h" namespace rttb { namespace models { /*! @class TCPLQModel @brief This class represents a TCP(Tumor Control Probability) LQ model (Nahum and Sanchez-Nieto 2001, Hall and Giaccia 2006) @see TCPModel */ class TCPLQModel: public TCPModel { public: typedef TCPModel::ParamVectorType ParamVectorType; typedef TCPModel::DVHPointer DVHPointer; private: /*! @brief Calculate intermediate tcp using alpha constant. This is a helper function for calcTCP() @see calcTCP */ long double calcTCPi(BioModelParamType aRho, BioModelParamType aAlphaMean, double vj, double bedj); /*! @brief Calculate tcp using alpha constant. */ long double calcTCP(std::map aBEDDVH, BioModelParamType aRho, BioModelParamType aAlphaMean, double aDeltaV); /*! @brief Calculate tcp using a normal distribution for alpha. */ long double calcTCPAlphaNormalDistribution(std::map aBEDDVH, BioModelParamType aRho, BioModelParamType aAlphaMean, BioModelParamType aAlphaVariance, double aDeltaV); protected: BioModelParamType _alphaMean; BioModelParamType _alphaVariance; BioModelParamType _alpha_beta; /*! Roh is the initial clonogenic cell density */ BioModelParamType _rho; /*! @brief Calculate the model value @param doseFactor scaling factor for prescribed dose. The model calculation will use the dvh with each di=old di*doseFactor. @pre _alphaMean >0 @pre _alphaVariance >= 0 @pre _alpha_beta > 0 @pre _rho > 0 @pre _numberOfFractions > 1 @exception InvalidParameterException Thrown if parameters were not set correctly. */ BioModelValueType calcModel(const double doseFactor = 1); public: TCPLQModel(); /*! @brief Constructor initializes member variables with given parameters. @pre aAlphaMean >0 @pre aBeta > 0 @pre aRho > 0 @pre aNumberOfFractions > 1 */ TCPLQModel(DVHPointer aDVH, BioModelParamType aAlphaMean, BioModelParamType aBeta, BioModelParamType aRho, int aNumberOfFractions); /*! @brief Constructor for alpha distribution initializes member variables with given parameters. @pre aAlphaMean >0 @pre aAlphaVariance >0 @pre aAlpha_Beta > 0 @pre aRho > 0 @pre aNumberOfFractions > 1 */ TCPLQModel(DVHPointer aDVH, BioModelParamType aRho, int aNumberOfFractions, BioModelParamType aAlpha_Beta, BioModelParamType aAlphaMean, BioModelParamType aAlphaVariance); const BioModelParamType getRho(); void setRho(const BioModelParamType aRho); const BioModelParamType getAlphaMean(); const BioModelParamType getAlphaVariance(); /*! @brief The distribution of the parameter alpha, which is characteristic for a population of cells, is described by the its mean and variance. If alpha is constant the variance is 0. @param aAlphaVariance The variance of alpha can be given, the default value is 0 resulting in constant alpha. */ void setAlpha(const BioModelParamType aAlphaMean, const BioModelParamType aAlphaVariance = 0); const BioModelParamType 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. */ virtual void setParameterByID(const int aParamId, const BioModelParamType aValue); /*! @brief Set parameter vector, where index of vector is the parameter id. "alphaMean":0,"alphaVariance":1,"alpha_beta":2, "rho":3 @exception InvalidParameterException Thrown if aParamterVector.size()!=4. */ virtual void setParameterVector(const ParamVectorType& aParameterVector); /*! @brief Get parameter id. "alphaMean":0,"alphaVariance":1,"alpha_beta":2, "rho":3 @return 0 for "alphaMean", 1 for "alphaVariance", 2 for "alpha_beta", 3 for "rho" @exception InvalidParameterException Thrown if aParamName is not alphaMean or alphaVariance or alpha_beta or rho. */ virtual const int getParameterID(const std::string& aParamName) const; - virtual std::map getParameterMap(); + virtual std::map getParameterMap() const; - virtual std::string getModelType() const; + void fillParameterMap(); + + virtual std::string getModelType() const override; }; }//end algorithms }//end rttb #endif diff --git a/testing/io/models/ModelsIOTest.cpp b/testing/io/models/ModelsIOTest.cpp index ba7a82c..42a3640 100644 --- a/testing/io/models/ModelsIOTest.cpp +++ b/testing/io/models/ModelsIOTest.cpp @@ -1,128 +1,128 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // (c) Copyright 2007, DKFZ, Heidelberg, Germany // ALL RIGHTS RESERVED // // THIS FILE CONTAINS CONFIDENTIAL AND PROPRIETARY INFORMATION OF DKFZ. // ANY DUPLICATION, MODIFICATION, DISTRIBUTION, OR // DISCLOSURE IN ANY FORM, IN WHOLE, OR IN PART, IS STRICTLY PROHIBITED // WITHOUT THE PRIOR EXPRESS WRITTEN PERMISSION OF DKFZ. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 1221 $ (last changed revision) // @date $Date: 2015-12-01 13:43:31 +0100 (Di, 01 Dez 2015) $ (last change date) // @author zhangl (last changed by) // @author *none* (Reviewer) // @author zhangl (Programmer) // // Subversion HeadURL: $HeadURL: http://sidt-hpc1/dkfz_repository/NotMeVisLab/SIDT/RTToolbox/trunk/testing/core/DVHCalculatorTest.cpp $ */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include "litCheckMacros.h" #include "rttbBioModel.h" #include "rttbTCPLQModel.h" #include "rttbNTCPLKBModel.h" #include #include #include "boost/filesystem.hpp" #include "rttbModelXMLWriter.h" namespace rttb { namespace testing { static std::string readFile(const std::string& filename); int ModelsIOTest(int argc, char* argv[]) { typedef core::DVH::DataDifferentialType DataDifferentialType; PREPARE_DEFAULT_TEST_REPORTING; //generate artificial DVH and corresponding statistical values DoseTypeGy binSize = DoseTypeGy(0.1); DoseVoxelVolumeType voxelVolume = 8; DataDifferentialType aDataDifferential; DoseCalcType value = 0; DVHVoxelNumber numberOfVoxels = 0; // creat default values for (int i = 0; i < 98; i++) { value = 0; numberOfVoxels += value; aDataDifferential.push_back(value); } aDataDifferential.push_back(10); aDataDifferential.push_back(20); const IDType structureID = "myStructure"; const IDType doseID = "myDose"; const IDType voxelizationID = "myVoxelization"; core::DVH::DVHPointer dvhPtr = boost::make_shared(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID); //test TCP LQ Model models::BioModelParamType alpha = 0.35; models::BioModelParamType beta = 0.023333333333333; models::BioModelParamType roh = 10000000; int numFractions = 8; - rttb::models::TCPLQModel tcplq = rttb::models::TCPLQModel(dvhPtr, roh, numFractions, alpha / beta, alpha, 0.08); + boost::shared_ptr tcplq = boost::make_shared(dvhPtr, roh, numFractions, alpha / beta, alpha, 0.08); std::string filename = "BioModeltcpleqIOTest.xml"; rttb::io::models::ModelXMLWriter writer = rttb::io::models::ModelXMLWriter(filename, tcplq); CHECK_NO_THROW(writer.writeModel()); CHECK_EQUAL(boost::filesystem::exists(filename), true); CHECK_EQUAL(std::remove(filename.c_str()), 0); //test NTCPLKBModel models::BioModelParamType aVal = 10; models::BioModelParamType mVal = 0.16; models::BioModelParamType d50Val = 35; - rttb::models::NTCPLKBModel ntcplk = rttb::models::NTCPLKBModel(dvhPtr, d50Val, mVal, aVal); + boost::shared_ptr ntcplk= boost::make_shared(dvhPtr, d50Val, mVal, aVal); filename = "BioModelntcplkIOTest.xml"; rttb::io::models::ModelXMLWriter writer2 = rttb::io::models::ModelXMLWriter(filename, ntcplk); CHECK_NO_THROW(writer2.writeModel()); CHECK_EQUAL(boost::filesystem::exists(filename), true); std::string defaultAsIs = readFile(filename); std::string defaultExpected = readFile("referenceBioModelntcplkIOTest.xml"); CHECK_EQUAL(defaultAsIs, defaultExpected); CHECK_EQUAL(std::remove(filename.c_str()), 0); RETURN_AND_REPORT_TEST_SUCCESS; } std::string readFile(const std::string& filename) { std::ifstream fileStream(filename.c_str()); std::string content((std::istreambuf_iterator(fileStream)), (std::istreambuf_iterator())); return content; } } } \ No newline at end of file diff --git a/testing/models/DummyModel.cpp b/testing/models/DummyModel.cpp index 7537e1c..d3848f3 100644 --- a/testing/models/DummyModel.cpp +++ b/testing/models/DummyModel.cpp @@ -1,121 +1,126 @@ // ----------------------------------------------------------------------- // 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(){ + 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 fb02095..742ba1c 100644 --- a/testing/models/DummyModel.h +++ b/testing/models/DummyModel.h @@ -1,96 +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 __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(); + 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