diff --git a/code/algorithms/rttbDoseStatistics.h b/code/algorithms/rttbDoseStatistics.h index 60a223c..713f9c4 100644 --- a/code/algorithms/rttbDoseStatistics.h +++ b/code/algorithms/rttbDoseStatistics.h @@ -1,229 +1,229 @@ // ----------------------------------------------------------------------- // 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 __DOSE_STATISTICS_H #define __DOSE_STATISTICS_H #include #include #include "boost/shared_ptr.hpp" #include "rttbBaseType.h" #include "RTTBAlgorithmsExports.h" #include "rttbVolumeToDoseMeasure.h" namespace rttb { namespace algorithms { /*! @class DoseStatistics @brief This is a data class storing different statistical values from a rt dose distribution @sa DoseStatisticsCalculator */ class RTTBAlgorithms_EXPORT DoseStatistics { public: enum complexStatistics { Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx }; typedef boost::shared_ptr > > ResultListPointer; typedef boost::shared_ptr DoseStatisticsPointer; typedef std::map DoseToVolumeFunctionType; typedef std::map VolumeToDoseFunctionType; private: double getValue(const std::map& aMap, double key, bool findNearestValueInstead, double& storedKey) const; std::map::const_iterator findNearestKeyInMap(const std::map& aMap, double key) const; DoseStatisticType _maximum; DoseStatisticType _minimum; ResultListPointer _maximumVoxelPositions; ResultListPointer _minimumVoxelPositions; DoseStatisticType _mean; DoseStatisticType _stdDeviation; VoxelNumberType _numVoxels; VolumeType _volume; DoseTypeGy _referenceDose; //for Vx computation VolumeToDoseMeasure _Dx; DoseToVolumeFunctionType _Vx; VolumeToDoseFunctionType _MOHx; VolumeToDoseFunctionType _MOCx; VolumeToDoseFunctionType _MaxOHx; VolumeToDoseFunctionType _MinOCx; public: /*! @brief Standard Constructor */ //DoseStatistics(); /*! @brief Constructor @details the dose statistic values are set. Complex values maximumVoxelLocation, maximumVoxelLocation, Dx, Vx, MOHx, MOCx, MaxOHx and MinOCx are optional */ DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean, DoseStatisticType stdDeviation, VoxelNumberType numVoxels, VolumeType volume, ResultListPointer minimumVoxelPositions = NULL, ResultListPointer maximumVoxelPositions = NULL, - VolumeToDoseMeasure _Dx = VolumeToDoseMeasure("Dx", std::map(), 0.0), + VolumeToDoseMeasure _Dx = VolumeToDoseMeasure("Dx"), DoseToVolumeFunctionType Vx = DoseToVolumeFunctionType(), VolumeToDoseFunctionType MOHx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MOCx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MaxOHx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MinOCx = VolumeToDoseFunctionType(), DoseTypeGy referenceDose = -1); ~DoseStatistics(); void setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions); void setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions); void setDx(const VolumeToDoseMeasure& DxValues); void setVx(const VolumeToDoseFunctionType& VxValues); void setMOHx(const VolumeToDoseFunctionType& MOHxValues); void setMOCx(const VolumeToDoseFunctionType& MOCxValues); void setMaxOHx(const VolumeToDoseFunctionType& MaxOHxValues); void setMinOCx(const VolumeToDoseFunctionType& MinOCxValues); void setReferenceDose(DoseTypeGy referenceDose); /*! @brief Get number of voxels in doseIterator, with sub-voxel accuracy. */ VoxelNumberType getNumberOfVoxels() const; /*! @brief Get the volume of the voxels in doseIterator (in cm3), with sub-voxel accuracy */ VolumeType getVolume() const; /*! @brief Get the reference dose for Vx computation */ DoseTypeGy getReferenceDose() const; /*! @brief Get the maximum of the current dose distribution. @return Return the maximum dose in Gy */ DoseStatisticType getMaximum() const; /*! @brief Get a vector of the the maximum dose VoxelGridIDs together with their dose value in Gy @exception InvalidDoseException if the vector has not been set (i.e. is empty) */ ResultListPointer getMaximumVoxelPositions() const; /*! @brief Get the minimum of the current dose distribution. @return Return the minimum dose in Gy */ DoseStatisticType getMinimum() const; /*! @brief Get a vector of the the minimum dose VoxelGridIDs together with their dose value in Gy @exception InvalidDoseException if the vector has not been set (i.e. is empty) */ ResultListPointer getMinimumVoxelPositions() const; /*! @brief Get the mean of the current dose distribution. @return Return the mean dose in Gy */ DoseStatisticType getMean() const; /*! @brief Get the standard deviation of the current dose distribution. @return Return the standard deviation in Gy */ DoseStatisticType getStdDeviation() const; /*! @brief Get the variance of of the current dose distribution. @return Return the variance in Gy */ DoseStatisticType getVariance() const; /*! @brief Get Vx: the volume irradiated with a dose >= x. @return Return absolute volume in absolute cm^3. @exception NoDataException if the Vx values have not been set (i.e. the vector is empty) @exception NoDataException if the requested Dose is not in the vector */ VolumeType getVx(DoseTypeGy xDoseAbsolute) const; VolumeType getVx(DoseTypeGy xDoseAbsolute, bool findNearestValue, DoseTypeGy& nearestXDose) const; VolumeType getVxRelative(DoseTypeGy xDoseRelative) const; VolumeType getVxRelative(DoseTypeGy xDoseRelative, bool findNearestValue, DoseTypeGy& nearestXDose) const; DoseToVolumeFunctionType getAllVx() const; VolumeToDoseMeasure getDx() const; /*! @brief Get MOHx: mean dose of the hottest x voxels. @return Return dose value in Gy. @exception InvalidDoseException if the values have not been set (i.e. the vector is empty) */ DoseTypeGy getMOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOHx(VolumeType xVolumeAbsolute) const; DoseTypeGy getMOHxRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOHxRelative(VolumeType xDoseRelative) const; VolumeToDoseFunctionType getAllMOHx() const; /*! @brief Get MOCx: mean dose of the coldest x voxels. @return Return dose value in Gy. @exception InvalidDoseException if the values have not been set (i.e. the vector is empty) */ DoseTypeGy getMOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOCx(VolumeType xVolumeAbsolute) const; DoseTypeGy getMOCxRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOCxRelative(VolumeType xDoseRelative) const; VolumeToDoseFunctionType getAllMOCx() const; /*! @brief Get MaxOHx: Maximum outside of the hottest x voxels. @return Return dose value in Gy. @exception InvalidDoseException if the values have not been set (i.e. the vector is empty) */ DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute) const; DoseTypeGy getMaxOHxRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXVolum) const; DoseTypeGy getMaxOHxRelative(VolumeType xDoseRelative) const; VolumeToDoseFunctionType getAllMaxOHx() const; /*! @brief Get MinOCx: Minimum outside of the coldest x voxels. @return Return dose value in Gy. @exception InvalidDoseException if the values have not been set (i.e. the vector is empty) */ DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute) const; DoseTypeGy getMinOCxRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMinOCxRelative(VolumeType xDoseRelative) const; VolumeToDoseFunctionType getAllMinOCx() const; }; } } #endif diff --git a/code/algorithms/rttbVolumeToDoseMeasure.cpp b/code/algorithms/rttbVolumeToDoseMeasure.cpp index 46579ef..12aff39 100644 --- a/code/algorithms/rttbVolumeToDoseMeasure.cpp +++ b/code/algorithms/rttbVolumeToDoseMeasure.cpp @@ -1,133 +1,138 @@ #include "rttbVolumeToDoseMeasure.h" #include "rttbDataNotAvailableException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace algorithms { VolumeToDoseMeasure::VolumeToDoseMeasure(std::string name, VolumeToDoseFunctionType values, VolumeType volume) : name(name), values(values), _volume(volume) {} std::string VolumeToDoseMeasure::getName() const { return this->name; } void VolumeToDoseMeasure::insertValue(std::pair value) { this->values.insert(value); } DoseTypeGy VolumeToDoseMeasure::getValue(VolumeType xVolumeAbsolute) { VolumeType dummy; return getSpecificValue(xVolumeAbsolute, false, dummy); } DoseTypeGy VolumeToDoseMeasure::getValue(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType & nearestXVolume) { return getSpecificValue(xVolumeAbsolute, findNearestValue, nearestXVolume); } DoseTypeGy VolumeToDoseMeasure::getValueRelative(VolumeType xVolumeRelative) { if (xVolumeRelative >= 0 && xVolumeRelative <= 1) { DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; VolumeType dummy; return getSpecificValue(xVolumeAbsolute, false, dummy); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } DoseTypeGy VolumeToDoseMeasure::getValueRelative(VolumeType xVolumeRelative, bool findNearestValue, VolumeType & nearestXVolume) { if (xVolumeRelative >= 0 && xVolumeRelative <= 1) { DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; return getSpecificValue(xVolumeAbsolute, findNearestValue, nearestXVolume); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } VolumeToDoseMeasure::VolumeToDoseFunctionType VolumeToDoseMeasure::getAllValues() const { return this->values; } + void VolumeToDoseMeasure::setVolume(VolumeType volume) + { + this->_volume = volume; + } + double VolumeToDoseMeasure::getSpecificValue(double key, bool findNearestValueInstead, double& storedKey) const { if (values.find(key) != std::end(values)) { return values.find(key)->second; } else { //value not in map. We have to find the nearest value if (values.empty()) { throw core::DataNotAvailableException("No Vx values are defined"); } else { if (findNearestValueInstead) { auto iterator = findNearestKeyInMap(values, key); storedKey = iterator->first; return iterator->second; } else { throw core::DataNotAvailableException("No Vx value with required dose is defined"); } } } } std::map::const_iterator VolumeToDoseMeasure::findNearestKeyInMap( const std::map& aMap, double key) const { double minDistance = 1e19; double minDistanceLast = 1e20; auto iterator = std::begin(aMap); while (iterator != std::end(aMap)) { minDistanceLast = minDistance; minDistance = fabs(iterator->first - key); if (minDistanceLast > minDistance) { ++iterator; } else { if (iterator != std::begin(aMap)) { --iterator; return iterator; } else { return std::begin(aMap); } } } --iterator; return iterator; } - + bool operator==(const VolumeToDoseMeasure& volumeToDoseMesure,const VolumeToDoseMeasure& otherVolumeToDoseMesure) { if (volumeToDoseMesure.getName() == otherVolumeToDoseMesure.getName() && volumeToDoseMesure.getAllValues() == otherVolumeToDoseMesure.getAllValues()) { return true; } return false; } -} + } } diff --git a/code/algorithms/rttbVolumeToDoseMeasure.h b/code/algorithms/rttbVolumeToDoseMeasure.h index c4bf284..0e7c9f0 100644 --- a/code/algorithms/rttbVolumeToDoseMeasure.h +++ b/code/algorithms/rttbVolumeToDoseMeasure.h @@ -1,68 +1,67 @@ // ----------------------------------------------------------------------- // 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: 1674 $ (last changed revision) // @date $Date: 2017-01-27 10:34:46 +0100 (Fr, 27 Jan 2017) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __VOLUME_TO_DOSE_MEASURE_H #define __VOLUME_TO_DOSE_MEASURE_H #include #include #include "rttbBaseType.h" #include "RTTBAlgorithmsExports.h" -#include "rttbDoseStatistics.h" - namespace rttb { namespace algorithms { class RTTBAlgorithms_EXPORT VolumeToDoseMeasure { public: typedef std::map VolumeToDoseFunctionType; private: std::string name; VolumeToDoseFunctionType values; VolumeType _volume; public: - VolumeToDoseMeasure(std::string name, VolumeToDoseFunctionType values, VolumeType volume); + VolumeToDoseMeasure(std::string name, VolumeToDoseFunctionType values = std::map(), VolumeType volume = -1); std::string getName() const; void insertValue(std::pair value); DoseTypeGy getValue(VolumeType xVolumeAbsolute); DoseTypeGy getValue(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXDose); DoseTypeGy getValueRelative(VolumeType xDoseRelative); DoseTypeGy getValueRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXDose); VolumeToDoseFunctionType getAllValues() const; friend bool operator==(const VolumeToDoseMeasure& volumeToDoseMesure, const VolumeToDoseMeasure& otherVolumeToDoseMesure); + void setVolume(VolumeType volume); private: double getSpecificValue(double key, bool findNearestValueInstead, double& storedKey) const; std::map::const_iterator findNearestKeyInMap(const std::map& aMap, double key) const; }; } } #endif diff --git a/code/algorithms/rttbVolumeToDoseMeasureCalculator.cpp b/code/algorithms/rttbVolumeToDoseMeasureCalculator.cpp index 96f6c8d..73ab35c 100644 --- a/code/algorithms/rttbVolumeToDoseMeasureCalculator.cpp +++ b/code/algorithms/rttbVolumeToDoseMeasureCalculator.cpp @@ -1,46 +1,46 @@ #include "rttbVolumeToDoseMeasureCalculator.h" #include //#include namespace rttb { namespace algorithms { VolumeToDoseMeasureCalculator::VolumeToDoseMeasureCalculator(const std::vector& precomputeVolumeValues, const VolumeType& volume, const std::vector& doseVector, const std::vector& voxelProportionVector, const DoseVoxelVolumeType& currentVoxelVolume) : - measure(VolumeToDoseMeasure("dx", std::map(), volume)), _precomputeVolumeValues(precomputeVolumeValues), _volume(volume), + measure(VolumeToDoseMeasure("dx")), _precomputeVolumeValues(precomputeVolumeValues), _volume(volume), _doseVector(doseVector), _voxelProportionVector(voxelProportionVector), _currentVoxelVolume(currentVoxelVolume) {} void VolumeToDoseMeasureCalculator::compute() { std::vector threads; for (size_t i = 0; i < _precomputeVolumeValues.size(); ++i) { if (false)//_multiThreading) { threads.push_back(boost::thread(&VolumeToDoseMeasureCalculator::computeSpecificValue, this, _precomputeVolumeValues.at(i) * _volume)); } else { this->computeSpecificValue(_precomputeVolumeValues.at(i) * _volume); } } for (unsigned int i = 0; i(xAbsolute, resultDose)); } } } diff --git a/code/algorithms/rttbVolumeToDoseMeasureCalculator.h b/code/algorithms/rttbVolumeToDoseMeasureCalculator.h index 5c06192..e53d614 100644 --- a/code/algorithms/rttbVolumeToDoseMeasureCalculator.h +++ b/code/algorithms/rttbVolumeToDoseMeasureCalculator.h @@ -1,69 +1,69 @@ // ----------------------------------------------------------------------- // 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: 1674 $ (last changed revision) // @date $Date: 2017-01-27 10:34:46 +0100 (Fr, 27 Jan 2017) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __VOLUME_TO_DOSE_MEASURE_CALCULATOR_H #define __VOLUME_TO_DOSE_MEASURE_CALCULATOR_H #include #include #include "rttbBaseType.h" #include "RTTBAlgorithmsExports.h" #include "rttbVolumeToDoseMeasure.h" namespace rttb { namespace algorithms { class RTTBAlgorithms_EXPORT VolumeToDoseMeasureCalculator { public: - typedef DoseStatistics::VolumeToDoseFunctionType VolumeToDoseFunctionType; + typedef std::map VolumeToDoseFunctionType; protected: std::vector _doseVector; DoseVoxelVolumeType _currentVoxelVolume; std::vector _voxelProportionVector; private: VolumeType _volume; VolumeToDoseMeasure measure; std::vector _precomputeVolumeValues; public: VolumeToDoseMeasureCalculator(const std::vector& precomputeVolumeValues, const VolumeType& volume, const std::vector& doseVector, const std::vector& voxelProportionVector, const DoseVoxelVolumeType& currentVoxelVolume); void compute(); VolumeToDoseMeasure getMeasure(); virtual void computeSpecificValue(double xAbsolute) = 0; protected: void insertIntoMeasure(VolumeType xAbsolute, DoseTypeGy resultDose); }; } } #endif diff --git a/code/io/other/rttbDoseStatisticsXMLReader.cpp b/code/io/other/rttbDoseStatisticsXMLReader.cpp index eb1e369..6990de5 100644 --- a/code/io/other/rttbDoseStatisticsXMLReader.cpp +++ b/code/io/other/rttbDoseStatisticsXMLReader.cpp @@ -1,215 +1,218 @@ // ----------------------------------------------------------------------- // 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 #include #include #include #include #include #include "rttbDoseStatisticsXMLReader.h" #include "rttbInvalidParameterException.h" +#include "rttbVolumeToDoseMeasure.h" + namespace rttb { namespace io { namespace other { typedef boost::shared_ptr DoseStatisticsPtr; DoseStatisticsXMLReader::DoseStatisticsXMLReader(const std::string& filename) : _filename(filename), _newFile(true) { } DoseStatisticsXMLReader::~DoseStatisticsXMLReader() { } void DoseStatisticsXMLReader::setFilename(const std::string& filename) { _filename = filename; _newFile = true; } DoseStatisticsPtr DoseStatisticsXMLReader::generateDoseStatistic() { if (_newFile) { this->createDoseStatistic(); } return _doseStatistic; } void DoseStatisticsXMLReader::createDoseStatistic() { boost::property_tree::ptree pt; // Load the XML file into the property tree. If reading fails // (cannot open file, parse error), an exception is thrown. try { read_xml(_filename, pt); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw rttb::core::InvalidParameterException("DoseStatistics file name invalid: could not open the xml file!"); } // data to fill the parameter for the DoseStatistics std::string name; std::string datum; unsigned int x; std::pair voxelid; std::vector < std::pair> vec; //initialize all parameters for the DoseStatistics double minimum=-1; double maximum=-1; double numVoxels=-1; double volume=-1; double referenceDose = -1; double mean=-1; double stdDeviation=-1; boost::shared_ptr > > minimumVoxelPositions = nullptr; boost::shared_ptr > > maximumVoxelPositions = nullptr; - std::map Dx; + rttb::algorithms::VolumeToDoseMeasure Dx = rttb::algorithms::VolumeToDoseMeasure("Dx"); std::map Vx; std::map MOHx; std::map MOCx; std::map MaxOHx; std::map MinOCx; BOOST_FOREACH(boost::property_tree::ptree::value_type & data, pt.get_child("statistics.results")) { datum = data.second.data(); BOOST_FOREACH(boost::property_tree::ptree::value_type & middle, data.second) { BOOST_FOREACH(boost::property_tree::ptree::value_type & innernode, middle.second) { std::string mia = innernode.first; if (innernode.first == "name") { name = innernode.second.data(); } else if (innernode.first == "voxelGridID") { boost::replace_all(datum, "\r\n", ""); boost::replace_all(datum, "\n", ""); boost::trim(datum); voxelid.first = boost::lexical_cast(datum); voxelid.second = boost::lexical_cast(innernode.second.data()); vec.push_back(voxelid); } else if (innernode.first == "x") { x = boost::lexical_cast(innernode.second.data()); } } } // fill with the extracted data if (name == "numberOfVoxels") { numVoxels = boost::lexical_cast(datum); } else if (name == "volume") { volume = boost::lexical_cast(datum); + Dx.setVolume(volume); } else if (name == "referenceDose") { referenceDose = boost::lexical_cast(datum); } else if (name == "mean") { mean = boost::lexical_cast(datum); } else if (name == "standardDeviation") { stdDeviation = boost::lexical_cast(datum); } else if (name == "minimum") { minimum = boost::lexical_cast(datum); if (!vec.empty()) { minimumVoxelPositions = boost::make_shared>>(vec); vec.clear(); } } else if (name == "maximum") { maximum = boost::lexical_cast(datum); if (!vec.empty()) { maximumVoxelPositions = boost::make_shared>>(vec); vec.clear(); } } else if (name == "Dx") { - Dx[boost::lexical_cast(x)*volume / 100] = boost::lexical_cast(datum); + Dx.insertValue(std::pair(boost::lexical_cast(x)*volume / 100, boost::lexical_cast(datum))); } else if (name == "Vx") { Vx[boost::lexical_cast(x)*referenceDose / 100] = boost::lexical_cast(datum); } else if (name == "MOHx") { MOHx[boost::lexical_cast(x)*volume / 100] = boost::lexical_cast(datum); } else if (name == "MOCx") { MOCx[boost::lexical_cast(x)*volume / 100] = boost::lexical_cast(datum); } else if (name == "MaxOHx") { MaxOHx[boost::lexical_cast(x)*volume / 100] = boost::lexical_cast(datum); } else if (name == "MinOCx") { MinOCx[boost::lexical_cast(x)*volume / 100] = boost::lexical_cast(datum); } } // make DoseStatistcs _doseStatistic = boost::make_shared( minimum, maximum, mean, stdDeviation, numVoxels, volume, minimumVoxelPositions, maximumVoxelPositions , Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx, referenceDose); } }//end namespace other }//end namespace io }//end namespace rttb diff --git a/code/io/other/rttbDoseStatisticsXMLWriter.cpp b/code/io/other/rttbDoseStatisticsXMLWriter.cpp index e18c864..4b968a6 100644 --- a/code/io/other/rttbDoseStatisticsXMLWriter.cpp +++ b/code/io/other/rttbDoseStatisticsXMLWriter.cpp @@ -1,302 +1,302 @@ // ----------------------------------------------------------------------- // 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 "rttbDoseStatisticsXMLWriter.h" #include #include #include #include "rttbInvalidParameterException.h" #include "rttbNullPointerException.h" namespace rttb { namespace io { namespace other { static const std::string xmlattrXTag = ".x"; static const std::string xmlattrNameTag = ".name"; static const std::string statisticsTag = "statistics.results"; static const std::string propertyTag = "property"; static const std::string columnSeparator = "@"; boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics) { using boost::property_tree::ptree; ptree pt; if (aDoseStatistics == nullptr){ throw core::NullPointerException("dose statistics is nullptr!"); } ptree numberOfVoxelsNode = createNodeWithNameAttribute(aDoseStatistics->getNumberOfVoxels(), "numberOfVoxels"); pt.add_child(statisticsTag + "." + propertyTag, numberOfVoxelsNode); ptree volumeNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getVolume()), "volume"); pt.add_child(statisticsTag + "." + propertyTag, volumeNode); ptree referenceNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getReferenceDose()), "referenceDose"); pt.add_child(statisticsTag + "." + propertyTag, referenceNode); ptree minimumNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getMinimum()), "minimum"); auto minimumPositions = aDoseStatistics->getMinimumVoxelPositions(); std::vector >::iterator pairItMin = minimumPositions->begin(); int count = 0; for (; pairItMin != minimumPositions->end() && count < 100; ++pairItMin) //output max. 100 minimum { VoxelGridID voxelID = pairItMin->second; ptree voxelMinPositions; voxelMinPositions.add("voxelGridID", voxelID); minimumNode.add_child("voxel", voxelMinPositions); count++; } pt.add_child(statisticsTag + "." + propertyTag, minimumNode); ptree maximumNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getMaximum()), "maximum"); auto maximumPositions = aDoseStatistics->getMaximumVoxelPositions(); std::vector >::iterator pairItMax = maximumPositions->begin(); count = 0; for (; pairItMax != maximumPositions->end() && count < 100; ++pairItMax) //output max. 100 maximum { VoxelGridID voxelID = pairItMax->second; ptree voxelMaxPositions; voxelMaxPositions.add("voxelGridID", voxelID); maximumNode.add_child("voxel", voxelMaxPositions); count++; } pt.add_child(statisticsTag + "." + propertyTag, maximumNode); ptree meanNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getMean()), "mean"); pt.add_child(statisticsTag + "." + propertyTag, meanNode); ptree sdNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getStdDeviation()), "standardDeviation"); pt.add_child(statisticsTag + "." + propertyTag, sdNode); ptree varianceNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getVariance()), "variance"); pt.add_child(statisticsTag + "." + propertyTag, varianceNode); double absoluteVolume = aDoseStatistics->getVolume(); double referenceDose = aDoseStatistics->getReferenceDose(); rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType AllVx = aDoseStatistics->getAllVx(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllDx = aDoseStatistics->getAllDx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllDx = aDoseStatistics->getDx().getAllValues(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getAllMOHx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getAllMOCx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getAllMaxOHx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getAllMinOCx(); rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType::iterator vxIt; rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType::iterator it; for (it = AllDx.begin(); it != AllDx.end(); ++it) { ptree DxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "Dx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, DxNode); } for (vxIt = AllVx.begin(); vxIt != AllVx.end(); ++vxIt) { ptree VxNode = createNodeWithNameAndXAttribute(static_cast(vxIt->second), "Vx", std::lround(convertToPercent(vxIt->first, referenceDose))); pt.add_child(statisticsTag + "." + propertyTag, VxNode); } for (it = AllMOHx.begin(); it != AllMOHx.end(); ++it) { ptree mohxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "MOHx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, mohxNode); } for (it = AllMOCx.begin(); it != AllMOCx.end(); ++it) { ptree mocxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "MOCx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, mocxNode); } for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); ++it) { ptree maxOhxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "MaxOHx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, maxOhxNode); } for (it = AllMinOCx.begin(); it != AllMinOCx.end(); ++it) { ptree minOCxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "MinOCx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, minOCxNode); } return pt; } void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName) { boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics); try { boost::property_tree::xml_parser::write_xml(aFileName, pt, std::locale(), boost::property_tree::xml_writer_make_settings('\t', 1)); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml failed: xml_parser_error!"); } } XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics) { boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics); std::stringstream sstr; try { boost::property_tree::xml_parser::write_xml(sstr, pt, boost::property_tree::xml_writer_make_settings('\t', 1)); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml to string failed: xml_parser_error!"); } return sstr.str(); } StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics) { if (aDoseStatistics == nullptr){ throw core::NullPointerException("dose statistics is nullptr!"); } std::stringstream sstr; sstr << static_cast(aDoseStatistics->getVolume() * 1000) << columnSeparator; // cm3 to mm3 sstr << static_cast(aDoseStatistics->getMaximum()) << columnSeparator; sstr << static_cast(aDoseStatistics->getMinimum()) << columnSeparator; sstr << static_cast(aDoseStatistics->getMean()) << columnSeparator; sstr << static_cast(aDoseStatistics->getStdDeviation()) << columnSeparator; sstr << static_cast(aDoseStatistics->getVariance()) << columnSeparator; rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType AllVx = aDoseStatistics->getAllVx(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllDx = aDoseStatistics->getAllDx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllDx = aDoseStatistics->getDx().getAllValues(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getAllMOHx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getAllMOCx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getAllMaxOHx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getAllMinOCx(); rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType::iterator vxIt; rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType::iterator it; for (it = AllDx.begin(); it != AllDx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } for (vxIt = AllVx.begin(); vxIt != AllVx.end(); ++vxIt) { // *1000 because of conversion cm3 to mm3 sstr << static_cast((*vxIt).second * 1000) << columnSeparator; } for (it = AllMOHx.begin(); it != AllMOHx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } for (it = AllMOCx.begin(); it != AllMOCx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } for (it = AllMinOCx.begin(); it != AllMinOCx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } return sstr.str(); } boost::property_tree::ptree createNodeWithNameAttribute(DoseTypeGy doseValue, const std::string& attributeName) { boost::property_tree::ptree node; node.put("", doseValue); node.put(xmlattrNameTag, attributeName); return node; } boost::property_tree::ptree createNodeWithNameAndXAttribute(DoseTypeGy doseValue, const std::string& attributeName, int xValue) { boost::property_tree::ptree node; node.put("", doseValue); node.put(xmlattrNameTag, attributeName); node.put(xmlattrXTag, xValue); return node; } double convertToPercent(double value, double maximum) { return (value / maximum) * 100; } }//end namespace other }//end namespace io }//end namespace rttb \ No newline at end of file diff --git a/testing/examples/RTDoseStatisticsDicomTest.cpp b/testing/examples/RTDoseStatisticsDicomTest.cpp index 540c25f..cd5d68e 100644 --- a/testing/examples/RTDoseStatisticsDicomTest.cpp +++ b/testing/examples/RTDoseStatisticsDicomTest.cpp @@ -1,118 +1,118 @@ // ----------------------------------------------------------------------- // 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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include "litCheckMacros.h" #include "rttbDoseStatistics.h" #include "rttbDoseStatisticsCalculator.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbBoostMaskAccessor.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbBaseType.h" namespace rttb { namespace testing { /*! @brief RTDoseStatisticsTest. Max, min, mean, standard deviation, variance, Vx, Dx, MOHx, MOCx, MaxOHx, MinOCx are tested. Test if calculation in new architecture returns similar results to the original implementation. WARNING: The values for comparison need to be adjusted if the input files are changed!*/ int RTDoseStatisticsDicomTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string RTDOSE_FILENAME; if (argc > 1) { RTDOSE_FILENAME = argv[1]; } else { std::cout << "at least one arguments required for RTDoseStatisticsDicomTest" << std::endl; return -1; } double expectedVal = 5.64477e-005; double volume = 24120.; typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericDoseIterator::DoseIteratorPointer DoseIteratorPointer; typedef algorithms::DoseStatisticsCalculator::ResultListPointer ResultListPointer; io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create corresponding DoseIterator ::boost::shared_ptr spDoseIteratorTmp = ::boost::make_shared(doseAccessor1); DoseIteratorPointer spDoseIterator(spDoseIteratorTmp); rttb::algorithms::DoseStatisticsCalculator doseStatisticsCalculator(spDoseIterator); std::vector precomputedDoseValues; precomputedDoseValues.push_back(0); precomputedDoseValues.push_back(0.5); precomputedDoseValues.push_back(1); std::vector precomputedVolumeValues; precomputedVolumeValues.push_back(20000 / volume); precomputedVolumeValues.push_back(1); rttb::algorithms::DoseStatistics::DoseStatisticsPointer doseStatistics = doseStatisticsCalculator.calculateDoseStatistics(precomputedDoseValues, precomputedVolumeValues); CHECK_CLOSE(doseStatistics->getMean(), expectedVal, errorConstant); CHECK_CLOSE(doseStatistics->getStdDeviation(), 0, errorConstant); CHECK_CLOSE(doseStatistics->getVariance(), 0, errorConstant); double dummy; DoseTypeGy vx = doseStatistics->getVx(expectedVal, true, dummy); CHECK_EQUAL(vx, doseStatistics->getVx(0)); - CHECK_CLOSE(expectedVal, doseStatistics->getDx(vx), reducedErrorConstant); + CHECK_CLOSE(expectedVal, doseStatistics->getDx().getValue(vx), reducedErrorConstant); CHECK_CLOSE(doseStatistics->getMaximum(), expectedVal, errorConstant); CHECK_CLOSE(doseStatistics->getMinimum(), expectedVal, errorConstant); ResultListPointer minListPtr = doseStatistics->getMinimumVoxelPositions(); ResultListPointer maxListPtr = doseStatistics->getMaximumVoxelPositions(); CHECK_EQUAL(maxListPtr->size(), 10); CHECK_EQUAL(minListPtr->size(), 10); - CHECK_CLOSE(doseStatistics->getDx(24120), doseStatistics->getMinimum(), 0.001); + CHECK_CLOSE(doseStatistics->getDx().getValue(24120), doseStatistics->getMinimum(), 0.001); CHECK_CLOSE(doseStatistics->getMOHx(24120), doseStatistics->getMean(), reducedErrorConstant); CHECK_CLOSE(doseStatistics->getMOCx(20000), doseStatistics->getMean(), reducedErrorConstant); CHECK_CLOSE(doseStatistics->getMinOCx(20000), doseStatistics->getMean(), reducedErrorConstant); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/io/other/CompareDoseStatistic.cpp b/testing/io/other/CompareDoseStatistic.cpp index 6c556b5..761a31f 100644 --- a/testing/io/other/CompareDoseStatistic.cpp +++ b/testing/io/other/CompareDoseStatistic.cpp @@ -1,116 +1,116 @@ // ----------------------------------------------------------------------- // 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: 1221 $ (last changed revision) // @date $Date: 2015-12-01 13:43:31 +0100 (Di, 01 Dez 2015) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #include "CompareDoseStatistic.h" namespace rttb { namespace testing { const double errorConstant = 1e-3;// errorConstant is so small because the double are castet to float when they are written bool checkEqualDoseStatistic(DoseStatisticsPtr aDoseStatistc1, DoseStatisticsPtr aDoseStatistic2){ bool result = lit::AreClose(aDoseStatistc1->getNumberOfVoxels(), aDoseStatistic2->getNumberOfVoxels(), 0.01); result = result && lit::AreClose(aDoseStatistc1->getVolume(), aDoseStatistic2->getVolume(), errorConstant); result = result && lit::AreClose(aDoseStatistc1->getMinimum(), aDoseStatistic2->getMinimum(), errorConstant); result = result && lit::AreClose(aDoseStatistc1->getMaximum(), aDoseStatistic2->getMaximum(), errorConstant); result = result && lit::AreClose(aDoseStatistc1->getMean(), aDoseStatistic2->getMean(), errorConstant); result = result && lit::AreClose(aDoseStatistc1->getStdDeviation(), aDoseStatistic2->getStdDeviation(), errorConstant); - result = result && mapCompare(aDoseStatistc1->getAllDx(), aDoseStatistic2->getAllDx()); + result = result && mapCompare(aDoseStatistc1->getDx().getAllValues(), aDoseStatistic2->getDx().getAllValues()); result = result && mapCompare(aDoseStatistc1->getAllVx(), aDoseStatistic2->getAllVx()); result = result && mapCompare(aDoseStatistc1->getAllMaxOHx(), aDoseStatistic2->getAllMaxOHx()); result = result && mapCompare(aDoseStatistc1->getAllMinOCx(), aDoseStatistic2->getAllMinOCx()); result = result && mapCompare(aDoseStatistc1->getAllMOCx(), aDoseStatistic2->getAllMOCx()); result = result && mapCompare(aDoseStatistc1->getAllMOHx(), aDoseStatistic2->getAllMOHx()); return result; } std::map::const_iterator findNearestKeyInMap(const std::map& aMap, double key) { double minDistance = 1e19; double minDistanceLast = 1e20; auto iterator = std::begin(aMap); while (iterator != std::end(aMap)) { minDistanceLast = minDistance; minDistance = fabs(iterator->first - key); if (minDistanceLast > minDistance) { ++iterator; } else { if (iterator != std::begin(aMap)) { --iterator; return iterator; } else { return std::begin(aMap); } } } --iterator; return iterator; } double getValue(const std::map& aMap, double key) { if (aMap.find(key) != std::end(aMap)) { return aMap.find(key)->second; } else { auto iterator = findNearestKeyInMap(aMap, key); return iterator->second; } } bool mapCompare(const std::map& lhs, const std::map& rhs) { if (lhs.size() != rhs.size()){ return false; } for (std::map::const_iterator i = rhs.cbegin(); i != rhs.cend(); ++i) { double a = i->second; double b = getValue(lhs, i->first) ; if (std::abs(a - b ) > errorConstant){// errorConstant is 1e-3 because the double-->float cast when they are written return false; } } return true; } }//testing }//rttb