diff --git a/code/algorithms/rttbDoseStatistics.cpp b/code/algorithms/rttbDoseStatistics.cpp index b3140e1..b289693 100644 --- a/code/algorithms/rttbDoseStatistics.cpp +++ b/code/algorithms/rttbDoseStatistics.cpp @@ -1,327 +1,292 @@ // ----------------------------------------------------------------------- // 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 "rttbDoseStatistics.h" #include "boost/make_shared.hpp" #include "rttbDataNotAvailableException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace algorithms { DoseStatistics::DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean, DoseStatisticType stdDeviation, VoxelNumberType numVoxels, VolumeType volume, ResultListPointer maximumVoxelPositions /*= ResultListPointer()*/, ResultListPointer minimumVoxelPositions /*= ResultListPointer()*/, VolumeToDoseMeasure Dx, - DoseToVolumeFunctionType Vx, + DoseToVolumeMeasure Vx, VolumeToDoseMeasure MOHx, VolumeToDoseMeasure MOCx, VolumeToDoseMeasure MaxOHx, VolumeToDoseMeasure MinOCx, DoseTypeGy referenceDose /*=-1*/): _minimum(minimum), _maximum(maximum), _mean(mean), _stdDeviation(stdDeviation), _numVoxels(numVoxels), _volume(volume), _Dx(Dx), _Vx(Vx), _MOHx(MOHx), _MOCx(MOCx), _MaxOHx(MaxOHx), _MinOCx(MinOCx) { if (maximumVoxelPositions == NULL) { _maximumVoxelPositions = boost::make_shared > > (std::vector >()); } else { _maximumVoxelPositions = maximumVoxelPositions; } if (minimumVoxelPositions == NULL) { _minimumVoxelPositions = boost::make_shared > > (std::vector >()); } else { _minimumVoxelPositions = minimumVoxelPositions; } if (referenceDose <= 0) { _referenceDose = _maximum; } else { _referenceDose = referenceDose; } } DoseStatistics::~DoseStatistics() { } void DoseStatistics::setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions) { _minimumVoxelPositions = minimumVoxelPositions; } void DoseStatistics::setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions) { _maximumVoxelPositions = maximumVoxelPositions; } void DoseStatistics::setDx(const VolumeToDoseMeasure& DxValues) { _Dx = DxValues; } - void DoseStatistics::setVx(const VolumeToDoseFunctionType& VxValues) + void DoseStatistics::setVx(const DoseToVolumeMeasure& VxValues) { _Vx = VxValues; } void DoseStatistics::setMOHx(const VolumeToDoseMeasure& MOHxValues) { _MOHx = MOHxValues; } void DoseStatistics::setMOCx(const VolumeToDoseMeasure& MOCxValues) { _MOCx = MOCxValues; } void DoseStatistics::setMaxOHx(const VolumeToDoseMeasure& MaxOHValues) { _MaxOHx = MaxOHValues; } void DoseStatistics::setMinOCx(const VolumeToDoseMeasure& MinOCValues) { _MinOCx = MinOCValues; } void DoseStatistics::setReferenceDose(DoseTypeGy referenceDose) { if (referenceDose <= 0) { _referenceDose = _maximum; } else { _referenceDose = referenceDose; } } VoxelNumberType DoseStatistics::getNumberOfVoxels() const { return _numVoxels; } VolumeType DoseStatistics::getVolume() const { return _volume; } DoseTypeGy DoseStatistics::getReferenceDose() const { return _referenceDose; } DoseStatisticType DoseStatistics::getMaximum() const { return _maximum; } DoseStatisticType DoseStatistics::getMinimum() const { return _minimum; } DoseStatisticType DoseStatistics::getMean() const { return _mean; } DoseStatisticType DoseStatistics::getStdDeviation() const { return _stdDeviation; } DoseStatisticType DoseStatistics::getVariance() const { return _stdDeviation * _stdDeviation; } - - VolumeType DoseStatistics::getVx(DoseTypeGy xDoseAbsolute, bool findNearestValue, - DoseTypeGy& nearestXDose) const - { - return getValue(_Vx, xDoseAbsolute, findNearestValue, nearestXDose); - } - - VolumeType DoseStatistics::getVx(DoseTypeGy xDoseAbsolute) const - { - DoseTypeGy dummy; - return getValue(_Vx, xDoseAbsolute, false, dummy); - } - VolumeType DoseStatistics::getVxRelative(DoseTypeGy xDoseRelative) const - { - if (_referenceDose != -1 && xDoseRelative >=0 && xDoseRelative <=1){ - DoseTypeGy xDoseAbsolute = xDoseRelative * _referenceDose; - DoseTypeGy dummy; - return getValue(_Vx, xDoseAbsolute, false, dummy); - } - else { - throw rttb::core::InvalidParameterException("Reference dose must be > 0 and 0 <= relative Dose <= 1"); - } - } - VolumeType DoseStatistics::getVxRelative(DoseTypeGy xDoseRelative, bool findNearestValue, - DoseTypeGy& nearestXDose) const - { - if (_referenceDose != -1 && xDoseRelative >= 0 && xDoseRelative <= 1){ - DoseTypeGy xDoseAbsolute = xDoseRelative * _referenceDose; - return getValue(_Vx, xDoseAbsolute, findNearestValue, nearestXDose); - } - else { - throw rttb::core::InvalidParameterException("Reference dose must be > 0 and 0 <= relative Dose <= 1"); - } - } - double DoseStatistics::getValue(const std::map& aMap, double key, bool findNearestValueInstead, double& storedKey) const { if (aMap.find(key) != std::end(aMap)) { return aMap.find(key)->second; } else { //value not in map. We have to find the nearest value if (aMap.empty()) { throw core::DataNotAvailableException("No Vx values are defined"); } else { if (findNearestValueInstead) { auto iterator = findNearestKeyInMap(aMap, key); storedKey = iterator->first; return iterator->second; } else { throw core::DataNotAvailableException("No Vx value with required dose is defined"); } } } } std::map::const_iterator DoseStatistics::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; } DoseStatistics::ResultListPointer DoseStatistics::getMaximumVoxelPositions() const { return _maximumVoxelPositions; } DoseStatistics::ResultListPointer DoseStatistics::getMinimumVoxelPositions() const { return _minimumVoxelPositions; } - DoseStatistics::DoseToVolumeFunctionType DoseStatistics::getAllVx() const + DoseToVolumeMeasure DoseStatistics::getVx() const { return _Vx; } VolumeToDoseMeasure DoseStatistics::getDx() const { return _Dx; } VolumeToDoseMeasure DoseStatistics::getMOHx() const { return _MOHx; } VolumeToDoseMeasure DoseStatistics::getMOCx() const { return _MOCx; } VolumeToDoseMeasure DoseStatistics::getMaxOHx() const { return _MaxOHx; } VolumeToDoseMeasure DoseStatistics::getMinOCx() const { return _MinOCx; } }//end namespace algorithms }//end namespace rttb diff --git a/code/algorithms/rttbDoseStatistics.h b/code/algorithms/rttbDoseStatistics.h index 84d6a8f..7138e0f 100644 --- a/code/algorithms/rttbDoseStatistics.h +++ b/code/algorithms/rttbDoseStatistics.h @@ -1,183 +1,170 @@ // ----------------------------------------------------------------------- // 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" +#include "rttbDoseToVolumeMeasure.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; + DoseToVolumeMeasure _Vx; VolumeToDoseMeasure _MOHx; VolumeToDoseMeasure _MOCx; VolumeToDoseMeasure _MaxOHx; VolumeToDoseMeasure _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(VolumeToDoseMeasure::Dx), - DoseToVolumeFunctionType Vx = DoseToVolumeFunctionType(), + DoseToVolumeMeasure Vx = DoseToVolumeMeasure(DoseToVolumeMeasure::Vx), VolumeToDoseMeasure MOHx = VolumeToDoseMeasure(VolumeToDoseMeasure::MOHx), VolumeToDoseMeasure MOCx = VolumeToDoseMeasure(VolumeToDoseMeasure::MOCx), VolumeToDoseMeasure MaxOHx = VolumeToDoseMeasure(VolumeToDoseMeasure::MaxOHx), VolumeToDoseMeasure MinOCx = VolumeToDoseMeasure(VolumeToDoseMeasure::MinOCx), DoseTypeGy referenceDose = -1); ~DoseStatistics(); void setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions); void setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions); void setDx(const VolumeToDoseMeasure& DxValues); - void setVx(const VolumeToDoseFunctionType& VxValues); + void setVx(const DoseToVolumeMeasure& VxValues); void setMOHx(const VolumeToDoseMeasure& MOHxValues); void setMOCx(const VolumeToDoseMeasure& MOCxValues); void setMaxOHx(const VolumeToDoseMeasure& MaxOHxValues); void setMinOCx(const VolumeToDoseMeasure& 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; + DoseToVolumeMeasure getVx() const; VolumeToDoseMeasure getMOHx() const; VolumeToDoseMeasure getMOCx() const; VolumeToDoseMeasure getMaxOHx() const; VolumeToDoseMeasure getMinOCx() const; }; } } #endif diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp index 4aed7b3..b5a5fb4 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,473 +1,400 @@ // ----------------------------------------------------------------------- // 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 "rttbDoseStatisticsCalculator.h" #include #include #include #include #include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" #include "rttbDxVolumeToDoseMeasureCalculator.h" +#include "rttbVxDoseToVolumeMeasureCalculator.h" #include "rttbMOHxVolumeToDoseMeasureCalculator.h" #include "rttbMOCxVolumeToDoseMeasureCalculator.h" #include "rttbMaxOHxVolumeToDoseMeasureCalculator.h" #include "rttbMinOCxVolumeToDoseMeasureCalculator.h" #include namespace rttb { namespace algorithms { DoseStatisticsCalculator::DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator) { if (aDoseIterator == NULL) { throw core::NullPointerException("DoseIterator must not be NULL"); } else { _doseIterator = aDoseIterator; } _simpleDoseStatisticsCalculated = false; _multiThreading = false; _mutex = ::boost::make_shared<::boost::shared_mutex>(); } DoseStatisticsCalculator::~DoseStatisticsCalculator() { } DoseStatisticsCalculator::DoseIteratorPointer DoseStatisticsCalculator::getDoseIterator() const { return _doseIterator; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( bool computeComplexMeasures, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } //"simple" dose statistics are mandatory calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); if (computeComplexMeasures) { //more complex dose statistics are optional with default maximum dose and default relative x values calculateComplexDoseStatistics(_statistics->getMaximum(), std::vector(), std::vector()); } return _statistics; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } if (referenceDose <= 0) { throw rttb::core::InvalidParameterException("Reference dose must be > 0 !"); } //simple dose statistics calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); //more complex dose statistics with given reference dose and default x values calculateComplexDoseStatistics(referenceDose, std::vector(), std::vector()); return _statistics; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues, DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } //"simple" dose statistics calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); if (referenceDose <= 0) { //more complex dose statistics with default maximum dose and relative x values calculateComplexDoseStatistics(_statistics->getMaximum(), precomputeDoseValues, precomputeVolumeValues); } else { //more complex dose statistics with given reference dose and relative x values calculateComplexDoseStatistics(referenceDose, precomputeDoseValues, precomputeVolumeValues); } return _statistics; } void DoseStatisticsCalculator::calculateSimpleDoseStatistics(unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { _doseVector.clear(); _voxelProportionVector.clear(); std::multimap doseValueVSIndexMap; std::vector voxelProportionVectorTemp; DoseStatisticType maximumDose = 0; DoseStatisticType minimumDose = std::numeric_limits::max(); DoseStatisticType meanDose; DoseStatisticType stdDeviationDose; DoseTypeGy sum = 0; VolumeType numVoxels = 0.0; DoseTypeGy squareSum = 0; VolumeType volume = 0; _doseIterator->reset(); int i = 0; DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid()) { doseValue = _doseIterator->getCurrentDoseValue(); if (i == 0) { minimumDose = doseValue; volume = _doseIterator->getCurrentVoxelVolume(); } rttb::FractionType voxelProportion = _doseIterator->getCurrentRelevantVolumeFraction(); sum += doseValue * voxelProportion; numVoxels += voxelProportion; squareSum += doseValue * doseValue * voxelProportion; if (doseValue > maximumDose) { maximumDose = doseValue; } else if (doseValue < minimumDose) { minimumDose = doseValue; } voxelProportionVectorTemp.push_back(voxelProportion); doseValueVSIndexMap.insert(std::pair(doseValue, i)); i++; _doseIterator->next(); } if (numVoxels != 0) { meanDose = sum / numVoxels; //standard deviation is defined only for n>=2 if (numVoxels >= 2) { //uncorrected variance is calculated DoseStatisticType varianceDose = (squareSum / numVoxels - meanDose * meanDose); if (varianceDose < errorConstant) { stdDeviationDose = 0; } else { stdDeviationDose = pow(varianceDose, 0.5); } } else { stdDeviationDose = 0; } } //sort dose values and corresponding volume fractions in member variables for (auto it = doseValueVSIndexMap.begin(); it != doseValueVSIndexMap.end(); ++it) { _doseVector.push_back((float)(*it).first); _voxelProportionVector.push_back(voxelProportionVectorTemp.at((*it).second)); } volume *= numVoxels; _statistics = boost::make_shared(minimumDose, maximumDose, meanDose, stdDeviationDose, numVoxels, volume); _simpleDoseStatisticsCalculated = true; ResultListPointer minimumVoxelPositions = computeMinimumPositions(maxNumberMinimaPositions); ResultListPointer maximumVoxelPositions = computeMaximumPositions(maxNumberMaximaPositions); _statistics->setMinimumVoxelPositions(minimumVoxelPositions); _statistics->setMaximumVoxelPositions(maximumVoxelPositions); } void DoseStatisticsCalculator::calculateComplexDoseStatistics(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues) { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call calculateComplexDoseStatistics()"); } std::vector precomputeDoseValuesNonConst = precomputeDoseValues; std::vector precomputeVolumeValuesNonConst = precomputeVolumeValues; //set default values if (precomputeDoseValues.empty()) { std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)( 0.95)(0.98); precomputeDoseValuesNonConst = defaultPrecomputeDoseValues; } if (precomputeVolumeValues.empty()) { std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)( 0.95)(0.98); precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues; } DxVolumeToDoseMeasureCalculator Dx = DxVolumeToDoseMeasureCalculator(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), _statistics->getMinimum(), VolumeToDoseMeasure::Dx); Dx.compute(); - DoseToVolumeFunctionType Vx = computeDoseToVolumeFunctionMulti(referenceDose, - precomputeDoseValuesNonConst, DoseStatistics::Vx); + VxDoseToVolumeMeasureCalculator Vx = VxDoseToVolumeMeasureCalculator(precomputeDoseValuesNonConst, referenceDose, _doseIterator, DoseToVolumeMeasure::Vx); + Vx.compute(); MOHxVolumeToDoseMeasureCalculator MOHx = MOHxVolumeToDoseMeasureCalculator(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), VolumeToDoseMeasure::MOHx); MOHx.compute(); MOCxVolumeToDoseMeasureCalculator MOCx = MOCxVolumeToDoseMeasureCalculator(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), VolumeToDoseMeasure::MOCx); MOCx.compute(); MaxOHxVolumeToDoseMeasureCalculator MaxOHx = MaxOHxVolumeToDoseMeasureCalculator(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), VolumeToDoseMeasure::MaxOHx); MaxOHx.compute(); MinOCxVolumeToDoseMeasureCalculator MinOCx = MinOCxVolumeToDoseMeasureCalculator(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), _statistics->getMinimum(), _statistics->getMaximum(), VolumeToDoseMeasure::MinOCx); MinOCx.compute(); - _statistics->setVx(Vx); + _statistics->setVx(Vx.getMeasure()); _statistics->setDx(Dx.getMeasure()); _statistics->setMOHx(MOHx.getMeasure()); _statistics->setMOCx(MOCx.getMeasure()); _statistics->setMaxOHx(MaxOHx.getMeasure()); _statistics->setMinOCx(MinOCx.getMeasure()); _statistics->setReferenceDose(referenceDose); } DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMaximumPositions( unsigned int maxNumberMaxima) const { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMaximumPositions()"); } ResultListPointer maxVoxelVector = boost::make_shared > >(); unsigned int count = 0; this->_doseIterator->reset(); DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid() && count < maxNumberMaxima) { doseValue = _doseIterator->getCurrentDoseValue(); if (doseValue == _statistics->getMaximum()) { VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID(); std::pair voxel(doseValue, currentID); maxVoxelVector->push_back(voxel); count++; } _doseIterator->next(); } return maxVoxelVector; } DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMinimumPositions( unsigned int maxNumberMinima) const { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMinimumPositions()"); } ResultListPointer minVoxelVector = boost::make_shared > >(); /*! @todo: Architecture Annotation: Finding the positions for the minimum only once reduces computation time, but will require sensible use by the programmers. To be save the output vector minVoxelVector will be always cleared here to garantee that no false values are presented. This change may be revoced to increase computation speed later on (only compute if(minVoxelVector->size()==0)). */ unsigned int count = 0; this->_doseIterator->reset(); DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid() && count < maxNumberMinima) { doseValue = _doseIterator->getCurrentDoseValue(); if (doseValue == _statistics->getMinimum()) { VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID(); std::pair voxel(doseValue, currentID); minVoxelVector->push_back(voxel); count++; } _doseIterator->next(); } return minVoxelVector; } - VolumeType DoseStatisticsCalculator::computeVx(DoseTypeGy xDoseAbsolute) const - { - rttb::FractionType count = 0; - _doseIterator->reset(); - - DoseTypeGy currentDose = 0; - - while (_doseIterator->isPositionValid()) - { - currentDose = _doseIterator->getCurrentDoseValue(); - - if (currentDose >= xDoseAbsolute) - { - count += _doseIterator->getCurrentRelevantVolumeFraction(); - } - - _doseIterator->next(); - } - - return count * this->_doseIterator->getCurrentVoxelVolume(); - } - - DoseStatisticsCalculator::DoseToVolumeFunctionType - DoseStatisticsCalculator::computeDoseToVolumeFunctionMulti(DoseTypeGy referenceDose, - const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const - { - std::vector threads; - - DoseToVolumeFunctionType VxMulti; - - for (size_t i = 0; i < precomputeDoseValues.size(); ++i) - { - if (_multiThreading) - { - threads.push_back(boost::thread(&DoseStatisticsCalculator::computeDoseToVolumeSingle, this, - referenceDose, precomputeDoseValues.at(i), name, std::ref(VxMulti))); - } - else - { - DoseStatisticsCalculator::computeDoseToVolumeSingle(referenceDose, precomputeDoseValues.at(i), - name, std::ref(VxMulti)); - } - } - - for (unsigned int i=0; i lock(*_mutex); - VxMulti.insert(std::pair(xAbsolue, - computeVx(xAbsolue))); - } - else { - VxMulti.insert(std::pair(xAbsolue, - computeVx(xAbsolue))); - } - } - else - { - throw core::InvalidParameterException("unknown DoseStatistics name!"); - } - } - void DoseStatisticsCalculator::setMultiThreading(const bool choice) { _multiThreading = choice; } }//end namespace algorithms }//end namespace rttb diff --git a/code/algorithms/rttbDoseStatisticsCalculator.h b/code/algorithms/rttbDoseStatisticsCalculator.h index 16edf49..ceea464 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.h +++ b/code/algorithms/rttbDoseStatisticsCalculator.h @@ -1,185 +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) */ #ifndef __DOSE_STATISTICS_CALCULATOR_H #define __DOSE_STATISTICS_CALCULATOR_H #include #include #include #include "rttbDoseIteratorInterface.h" #include "rttbDoseStatistics.h" #include "RTTBAlgorithmsExports.h" namespace rttb { namespace algorithms { /*! @class DoseStatisticsCalculator @brief Class for calculating different statistical values from a RT dose distribution @details These values range from standard statistical values such as minimum, maximum and mean to more complex dose specific measures such as Vx (volume irradiated with a dose >=x), Dx (minimal dose delivered to x% of the VOI) or MOHx (mean in the hottest volume). For a complete list, see calculateDoseStatistics(). @note the complex dose statistics are precomputed and cannot be computed "on the fly" lateron! The doses/volumes that should be used for precomputation have to be set in calculateDoseStatistics() */ class RTTBAlgorithms_EXPORT DoseStatisticsCalculator { public: typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; typedef DoseStatistics::ResultListPointer ResultListPointer; typedef DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer; - typedef DoseStatistics::DoseToVolumeFunctionType DoseToVolumeFunctionType; - typedef DoseStatistics::VolumeToDoseFunctionType VolumeToDoseFunctionType; private: DoseIteratorPointer _doseIterator; /*! @brief Contains relevant dose values sorted in descending order. */ std::vector _doseVector; /*! @brief Contains the corresponding voxel proportions to the values in doseVector. */ std::vector _voxelProportionVector; /*! @brief The doseStatistics are stored here. */ DoseStatisticsPointer _statistics; bool _simpleDoseStatisticsCalculated; bool _multiThreading; ::boost::shared_ptr _mutex; /*! @brief Calculates the positions where the dose has its maximum @param maxNumberMaximaPositions the maximal amount of computed positions @pre maximumDose must be defined in _statistics with the correct value */ ResultListPointer computeMaximumPositions(unsigned int maxNumberMaximaPositions) const; /*! @brief Calculates the positions where the dose has its minimum @param maxNumberMinimaPositions the maximal amount of computed positions (they are read sequentially using the iterator until maxNumberMinimaPositions have been read, other positions are not considered) @pre minimumDose must be defined in _statistics with the correct value */ ResultListPointer computeMinimumPositions(unsigned int maxNumberMinimaPositions) const; - - VolumeType computeVx(DoseTypeGy xDoseAbsolute) const; - - DoseToVolumeFunctionType computeDoseToVolumeFunctionMulti(DoseTypeGy referenceDose, - const std::vector& precomputeDoseValues, - DoseStatistics::complexStatistics name) const; - - void computeDoseToVolumeSingle(DoseTypeGy referenceDose, - double precomputeDoseValue, DoseStatistics::complexStatistics name, DoseToVolumeFunctionType& VxMulti) const; - + /*! @brief Calculates simple dose statistics (min, mean, max, stdDev, minDosePositions, maxDosePositions) @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed */ void calculateSimpleDoseStatistics(unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions); /*! @brief Calculates complex dose statistics (Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx) @warning computations can take quite long (>1 min) for large structures as many statistics are precomputed */ void calculateComplexDoseStatistics(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues); public: ~DoseStatisticsCalculator(); /*! @brief Constructor @param aDoseIterator the dose to be analyzed */ DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator); DoseIteratorPointer getDoseIterator() const; /*! @brief Compute simple or complex dose statistics with default relative x values and the maximum dose as default reference dose (for Vx computation) @details The following statistics are calculated always (i.e. also if computeComplexMeasures=false):
  • minimum dose
  • mean dose
  • maximum dose
  • standard deviation dose
  • voxel positions of minimum dose
  • voxel positions of maximum dose
Additionally, these statistics are computed if computeComplexMeasures=true:
  • Dx (the minimal dose delivered to a volume >= x)
  • Vx (the volume irradiated with a dose >= x)
  • MOHx (mean dose of the hottest x volume)
  • MOCx (mean dose of the coldest x volume)
  • MaxOHx (Maximum outside of the hottest x volume)
  • MinOCx (Minimum outside of the coldest x volume)
Default x values for Vx are 0.02, 0.05, 0.1, 0.9, 0.95 and 0.98, with respect to maxDose. Default x values for Dx, MOHx, MOCx, MaxOHx and MinOCx are 0.02, 0.05, 0.1, 0.9, 0.95 and 0.98, with respect to volume. @param computeComplexMeasures should complex statistics be calculated? If it is true, the complex dose statistics will be calculated with default relative x values and the maximum dose as reference dose @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed @warning If computeComplexMeasures==true, computations can take quite long (>1 min) for large structures as many statistics are precomputed @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! Only the default x values can be requested in DoseStatistics! */ DoseStatisticsPointer calculateDoseStatistics(bool computeComplexMeasures = false, unsigned int maxNumberMinimaPositions = 10, unsigned int maxNumberMaximaPositions = 10); /*! @brief Compute complex dose statistics with given reference dose and default relative x values @param referenceDose the reference dose to compute Vx, normally it should be the prescribed dose @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed @exception InvalidParameterException thrown if referenceDose <= 0 @warning Computations can take quite long (>1 min) for large structures as many statistics are precomputed @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! Only the default x values can be requested in DoseStatistics! */ DoseStatisticsPointer calculateDoseStatistics(DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions = 10, unsigned int maxNumberMaximaPositions = 10); /*! @brief Compute complex dose statistics with given relative x values and reference dose @param precomputeDoseValues the relative dose values for Vx precomputation, e.g. 0.02, 0.05, 0.95... @param precomputeVolumeValues the relative volume values for Dx, MOHx, MOCx, MaxOHx and MinOCx precomputation, e.g. 0.02, 0.05, 0.95... @param referenceDose the reference dose to compute Vx, normally it should be the prescribed dose. Default value is the maximum dose. @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed @warning Computations can take quite long (>1 min) for large structures as many statistics are precomputed @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! The doses/volumes that should be used for precomputation have to be set by in precomputeDoseValues and precomputeVolumeValues. Only these values can be requested in DoseStatistics! */ DoseStatisticsPointer calculateDoseStatistics(const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues, DoseTypeGy referenceDose = -1, unsigned int maxNumberMinimaPositions = 10, unsigned int maxNumberMaximaPositions = 10); void setMultiThreading(bool choice); }; } } #endif diff --git a/code/io/other/rttbDoseStatisticsXMLReader.cpp b/code/io/other/rttbDoseStatisticsXMLReader.cpp index 1688303..3838d2c 100644 --- a/code/io/other/rttbDoseStatisticsXMLReader.cpp +++ b/code/io/other/rttbDoseStatisticsXMLReader.cpp @@ -1,218 +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; rttb::algorithms::VolumeToDoseMeasure Dx = rttb::algorithms::VolumeToDoseMeasure(rttb::algorithms::VolumeToDoseMeasure::Dx); + rttb::algorithms::DoseToVolumeMeasure Vx = rttb::algorithms::DoseToVolumeMeasure(rttb::algorithms::DoseToVolumeMeasure::Vx); rttb::algorithms::VolumeToDoseMeasure MOHx = rttb::algorithms::VolumeToDoseMeasure(rttb::algorithms::VolumeToDoseMeasure::MOHx); rttb::algorithms::VolumeToDoseMeasure MOCx = rttb::algorithms::VolumeToDoseMeasure(rttb::algorithms::VolumeToDoseMeasure::MOCx); rttb::algorithms::VolumeToDoseMeasure MaxOHx = rttb::algorithms::VolumeToDoseMeasure(rttb::algorithms::VolumeToDoseMeasure::MaxOHx); rttb::algorithms::VolumeToDoseMeasure MinOCx = rttb::algorithms::VolumeToDoseMeasure(rttb::algorithms::VolumeToDoseMeasure::MinOCx); - std::map Vx; 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.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); + Vx.insertValue(std::pair(boost::lexical_cast(x)*referenceDose / 100, boost::lexical_cast(datum))); } else if (name == "MOHx") { MOHx.insertValue(std::pair(boost::lexical_cast(x)*volume / 100, boost::lexical_cast(datum))); } else if (name == "MOCx") { MOCx.insertValue(std::pair(boost::lexical_cast(x)*volume / 100, boost::lexical_cast(datum))); } else if (name == "MaxOHx") { MaxOHx.insertValue(std::pair(boost::lexical_cast(x)*volume / 100, boost::lexical_cast(datum))); } else if (name == "MinOCx") { MinOCx.insertValue(std::pair(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 8e26b33..a12d786 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" - +#include "rttbVolumeToDoseMeasure.h" +#include "rttbDoseToVolumeMeasure.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->getDx().getAllValues(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getMOHx().getAllValues(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getMOCx().getAllValues(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMaxOHx = + rttb::algorithms::DoseToVolumeMeasure::DoseToVolumeFunctionType AllVx = aDoseStatistics->getVx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllDx = aDoseStatistics->getDx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getMOHx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getMOCx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getMaxOHx().getAllValues(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMinOCx = + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getMinOCx().getAllValues(); - rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType::iterator vxIt; - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType::iterator it; + rttb::algorithms::DoseToVolumeMeasure::DoseToVolumeFunctionType::iterator vxIt; + rttb::algorithms::VolumeToDoseMeasure::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->getDx().getAllValues(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getMOHx().getAllValues(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getMOCx().getAllValues(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMaxOHx = - aDoseStatistics->getMaxOHx().getAllValues(); - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMinOCx = - aDoseStatistics->getMinOCx().getAllValues(); - - - rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType::iterator vxIt; - rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType::iterator it; + rttb::algorithms::DoseToVolumeMeasure::DoseToVolumeFunctionType AllVx = aDoseStatistics->getVx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllDx = aDoseStatistics->getDx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getMOHx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getMOCx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllMaxOHx = + aDoseStatistics->getMaxOHx().getAllValues(); + rttb::algorithms::VolumeToDoseMeasure::VolumeToDoseFunctionType AllMinOCx = + aDoseStatistics->getMinOCx().getAllValues(); + + rttb::algorithms::DoseToVolumeMeasure::DoseToVolumeFunctionType::iterator vxIt; + rttb::algorithms::VolumeToDoseMeasure::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/algorithms/DoseStatisticsCalculatorTest.cpp b/testing/algorithms/DoseStatisticsCalculatorTest.cpp index 98e4568..1e40a2d 100644 --- a/testing/algorithms/DoseStatisticsCalculatorTest.cpp +++ b/testing/algorithms/DoseStatisticsCalculatorTest.cpp @@ -1,374 +1,374 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbGenericDoseIterator.h" #include "rttbDoseIteratorInterface.h" #include "rttbNullPointerException.h" #include "rttbDoseStatisticsCalculator.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" #include "rttbDataNotAvailableException.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbVOIindexIdentifier.h" #include "rttbBoostMaskAccessor.h" #include "rttbGenericMaskedDoseIterator.h" #include "../io/other/CompareDoseStatistic.h" #include "../../code/io/other/rttbDoseStatisticsXMLReader.h" #include "../core/DummyDoseAccessor.h" namespace rttb { namespace testing { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; typedef rttb::algorithms::DoseStatistics::ResultListPointer ResultListPointer; typedef rttb::algorithms::DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer; /*! @brief DoseStatisticsCalculatorTest - test the API of DoseStatisticsCalculator 1) test constructors 2) test setDoseIterator 3) test calculateDoseSatistics 4) get statistical values */ int DoseStatisticsCalculatorTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string referenceXMLFilename; std::string doseFilename, structFilename; boost::shared_ptr spTestDoseAccessor = boost::make_shared(); DoseAccessorPointer spDoseAccessor(spTestDoseAccessor); const std::vector* doseVals = spTestDoseAccessor->getDoseVector(); boost::shared_ptr spTestDoseIterator = boost::make_shared(spDoseAccessor); DoseIteratorPointer spDoseIterator(spTestDoseIterator); DoseIteratorPointer spDoseIteratorNull; if (argc > 3) { referenceXMLFilename = argv[1]; doseFilename = argv[2]; structFilename = argv[3]; } //1) test constructors // the values cannot be accessed from outside, therefore correct default values are not tested CHECK_THROW_EXPLICIT(rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator( spDoseIteratorNull), core::NullPointerException); CHECK_NO_THROW(rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator)); rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator); //2) test setDoseIterator //3) test calculateDoseStatistics DoseStatisticsPointer theStatistics; //simple dose statistics CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics()); CHECK_EQUAL(theStatistics->getMinimumVoxelPositions()->empty(), false); CHECK_EQUAL(theStatistics->getMaximumVoxelPositions()->empty(), false); - CHECK_EQUAL(theStatistics->getAllVx().empty(), true); + CHECK_EQUAL(theStatistics->getVx().getAllValues().empty(), true); CHECK_EQUAL(theStatistics->getDx().getAllValues().empty(), true); - CHECK_EQUAL(theStatistics->getAllVx().empty(), true); + CHECK_EQUAL(theStatistics->getVx().getAllValues().empty(), true); CHECK_EQUAL(theStatistics->getMaxOHx().getAllValues().empty(), true); CHECK_EQUAL(theStatistics->getMOHx().getAllValues().empty(), true); CHECK_EQUAL(theStatistics->getMOCx().getAllValues().empty(), true); CHECK_EQUAL(theStatistics->getMinOCx().getAllValues().empty(), true); //check default values for computeComplexMeasures=true DoseStatisticsPointer theStatisticsDefault; myDoseStatsCalculator.setMultiThreading(true); CHECK_NO_THROW(theStatisticsDefault = myDoseStatsCalculator.calculateDoseStatistics(true)); - CHECK_NO_THROW(theStatisticsDefault->getVx(0.02 * theStatisticsDefault->getMaximum())); - CHECK_NO_THROW(theStatisticsDefault->getVx(0.05 * theStatisticsDefault->getMaximum())); - CHECK_NO_THROW(theStatisticsDefault->getVx(0.1 * theStatisticsDefault->getMaximum())); - CHECK_NO_THROW(theStatisticsDefault->getVx(0.9 * theStatisticsDefault->getMaximum())); - CHECK_NO_THROW(theStatisticsDefault->getVx(0.95 * theStatisticsDefault->getMaximum())); - CHECK_NO_THROW(theStatisticsDefault->getVx(0.98 * theStatisticsDefault->getMaximum())); + CHECK_NO_THROW(theStatisticsDefault->getVx().getValue(0.02 * theStatisticsDefault->getMaximum())); + CHECK_NO_THROW(theStatisticsDefault->getVx().getValue(0.05 * theStatisticsDefault->getMaximum())); + CHECK_NO_THROW(theStatisticsDefault->getVx().getValue(0.1 * theStatisticsDefault->getMaximum())); + CHECK_NO_THROW(theStatisticsDefault->getVx().getValue(0.9 * theStatisticsDefault->getMaximum())); + CHECK_NO_THROW(theStatisticsDefault->getVx().getValue(0.95 * theStatisticsDefault->getMaximum())); + CHECK_NO_THROW(theStatisticsDefault->getVx().getValue(0.98 * theStatisticsDefault->getMaximum())); CHECK_NO_THROW(theStatisticsDefault->getDx().getValue(0.02 * theStatisticsDefault->getVolume())); CHECK_NO_THROW(theStatisticsDefault->getDx().getValue(0.05 * theStatisticsDefault->getVolume())); CHECK_NO_THROW(theStatisticsDefault->getDx().getValue(0.1 * theStatisticsDefault->getVolume())); CHECK_NO_THROW(theStatisticsDefault->getDx().getValue(0.9 * theStatisticsDefault->getVolume())); CHECK_NO_THROW(theStatisticsDefault->getDx().getValue(0.95 * theStatisticsDefault->getVolume())); CHECK_NO_THROW(theStatisticsDefault->getDx().getValue(0.98 * theStatisticsDefault->getVolume())); //check manually set reference dose and the default x values CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(100.0)); - CHECK_THROW_EXPLICIT(theStatistics->getVx(0.1 * theStatistics->getMaximum()), + CHECK_THROW_EXPLICIT(theStatistics->getVx().getValue(0.1 * theStatistics->getMaximum()), core::DataNotAvailableException); - CHECK_NO_THROW(theStatistics->getVx(0.1 * 100.0)); + CHECK_NO_THROW(theStatistics->getVx().getValue(0.1 * 100.0)); CHECK_NO_THROW(theStatistics->getDx().getValue(0.1 * theStatistics->getVolume())); CHECK_NO_THROW(theStatistics->getDx().getValue(0.9 * theStatistics->getVolume())); CHECK_NO_THROW(theStatistics->getMOHx().getValue(0.95 * theStatistics->getVolume())); CHECK_NO_THROW(theStatistics->getMOCx().getValue(0.98 * theStatistics->getVolume())); CHECK_EQUAL(theStatistics->getReferenceDose(), 100.0); //check manually set x values std::vector precomputeDoseValues, precomputeVolumeValues; precomputeDoseValues.push_back(0.01); precomputeDoseValues.push_back(0.02); precomputeDoseValues.push_back(0.05); precomputeVolumeValues.push_back(0.9); precomputeVolumeValues.push_back(0.95); precomputeVolumeValues.push_back(0.99); CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(precomputeDoseValues, precomputeVolumeValues)); - CHECK_NO_THROW(theStatistics->getVx(0.01 * theStatistics->getMaximum())); - CHECK_NO_THROW(theStatistics->getVx(0.02 * theStatistics->getMaximum())); - CHECK_NO_THROW(theStatistics->getVx(0.05 * theStatistics->getMaximum())); - CHECK_THROW_EXPLICIT(theStatistics->getVx(0.03 * theStatistics->getMaximum()), + CHECK_NO_THROW(theStatistics->getVx().getValue(0.01 * theStatistics->getMaximum())); + CHECK_NO_THROW(theStatistics->getVx().getValue(0.02 * theStatistics->getMaximum())); + CHECK_NO_THROW(theStatistics->getVx().getValue(0.05 * theStatistics->getMaximum())); + CHECK_THROW_EXPLICIT(theStatistics->getVx().getValue(0.03 * theStatistics->getMaximum()), core::DataNotAvailableException); CHECK_NO_THROW(theStatistics->getDx().getValue(0.9 * theStatistics->getVolume())); CHECK_NO_THROW(theStatistics->getDx().getValue(0.95 * theStatistics->getVolume())); CHECK_NO_THROW(theStatistics->getDx().getValue(0.99 * theStatistics->getVolume())); CHECK_THROW_EXPLICIT(theStatistics->getDx().getValue(0.03 * theStatistics->getVolume()), core::DataNotAvailableException); - CHECK_EQUAL(theStatistics->getVx(0.02 * theStatistics->getMaximum()), - theStatisticsDefault->getVx(0.02 * theStatistics->getMaximum())); - CHECK_EQUAL(theStatistics->getVx(0.05 * theStatistics->getMaximum()), - theStatisticsDefault->getVx(0.05 * theStatistics->getMaximum())); + CHECK_EQUAL(theStatistics->getVx().getValue(0.02 * theStatistics->getMaximum()), + theStatisticsDefault->getVx().getValue(0.02 * theStatistics->getMaximum())); + CHECK_EQUAL(theStatistics->getVx().getValue(0.05 * theStatistics->getMaximum()), + theStatisticsDefault->getVx().getValue(0.05 * theStatistics->getMaximum())); CHECK_EQUAL(theStatistics->getDx().getValue(0.9 * theStatistics->getVolume()), theStatisticsDefault->getDx().getValue(0.9 * theStatistics->getVolume())); CHECK_EQUAL(theStatistics->getDx().getValue(0.95 * theStatistics->getVolume()), theStatisticsDefault->getDx().getValue(0.95 * theStatistics->getVolume())); //check manually set reference dose and x values CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(precomputeDoseValues, precomputeVolumeValues, 100.0)); - CHECK_THROW_EXPLICIT(theStatistics->getVx(0.01 * theStatistics->getMaximum()), + CHECK_THROW_EXPLICIT(theStatistics->getVx().getValue(0.01 * theStatistics->getMaximum()), core::DataNotAvailableException); - CHECK_NO_THROW(theStatistics->getVx(0.01 * 100.0)); + CHECK_NO_THROW(theStatistics->getVx().getValue(0.01 * 100.0)); CHECK_NO_THROW(theStatistics->getDx().getValue(0.9 * theStatistics->getVolume())); CHECK_EQUAL(theStatistics->getReferenceDose(), 100.0); //MOHx, MOCx, MaxOHx and MinOCx are computed analogous to Dx, they will not be checked. //4) get statistical values CHECK_EQUAL(theStatistics->getNumberOfVoxels(), doseVals->size()); //compute simple statistical values (min, mean, max, stddev) for comparison DoseStatisticType maximum = 0; DoseStatisticType minimum = 1000000; DoseStatisticType mean = 0; DoseStatisticType variance = 0; std::vector::const_iterator doseIt = doseVals->begin(); while (doseIt != doseVals->end()) { if (maximum < *doseIt) { maximum = *doseIt; } if (minimum > *doseIt) { minimum = *doseIt; } mean += *doseIt; ++doseIt; } mean /= doseVals->size(); DoseTypeGy compMean = (int(mean * 100)) / 100; doseIt = doseVals->begin(); while (doseIt != doseVals->end()) { variance += pow(*doseIt - mean, 2); ++doseIt; } variance /= doseVals->size(); DoseStatisticType stdDev = pow(variance, 0.5); //we have some precision problems here... double errorConstantLarger = 1e-2; CHECK_EQUAL(theStatistics->getMaximum(), maximum); CHECK_EQUAL(theStatistics->getMinimum(), minimum); CHECK_CLOSE(theStatistics->getMean(), mean, errorConstantLarger); CHECK_CLOSE(theStatistics->getStdDeviation(), stdDev, errorConstantLarger); CHECK_CLOSE(theStatistics->getVariance(), variance, errorConstantLarger); //check for complex doseStatistics (maximumPositions, minimumPositions, Vx, Dx, MOHx, MOCx, MAXOHx, MinOCx) unsigned int nMax = 0, nMin = 0; doseIt = doseVals->begin(); while (doseIt != doseVals->end()) { if (*doseIt == theStatistics->getMaximum()) { nMax++; } if (*doseIt == theStatistics->getMinimum()) { nMin++; } ++doseIt; } //only 100 positions are stored if (nMax > 100) { nMax = 100; } if (nMin > 100) { nMin = 100; } auto maximaPositions = theStatistics->getMaximumVoxelPositions(); auto minimaPositions = theStatistics->getMinimumVoxelPositions(); CHECK_EQUAL(maximaPositions->size(), nMax); CHECK_EQUAL(minimaPositions->size(), nMin); for (auto maximaPositionsIterator = std::begin(*maximaPositions); maximaPositionsIterator != std::end(*maximaPositions); ++maximaPositionsIterator) { CHECK_EQUAL(maximaPositionsIterator->first, theStatistics->getMaximum()); } for (auto minimaPositionsIterator = std::begin(*minimaPositions); minimaPositionsIterator != std::end(*minimaPositions); ++minimaPositionsIterator) { CHECK_EQUAL(minimaPositionsIterator->first, theStatistics->getMinimum()); } //generate specific example dose maximum = 9.5; minimum = 2.5; mean = 6; int sizeTemplate = 500; std::vector aDoseVector; for (int i = 0; i < sizeTemplate; i++) { aDoseVector.push_back(maximum); aDoseVector.push_back(minimum); } core::GeometricInfo geoInfo = spTestDoseAccessor->getGeometricInfo(); geoInfo.setNumRows(20); geoInfo.setNumColumns(10); geoInfo.setNumSlices(5); boost::shared_ptr spTestDoseAccessor2 = boost::make_shared(aDoseVector, geoInfo); DoseAccessorPointer spDoseAccessor2(spTestDoseAccessor2); boost::shared_ptr spTestDoseIterator2 = boost::make_shared(spDoseAccessor2); DoseIteratorPointer spDoseIterator2(spTestDoseIterator2); rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator2(spDoseIterator2); DoseStatisticsPointer theStatistics3 = myDoseStatsCalculator2.calculateDoseStatistics(); CHECK_EQUAL(theStatistics3->getMaximum(), maximum); CHECK_EQUAL(theStatistics3->getMinimum(), minimum); CHECK_EQUAL(theStatistics3->getMean(), mean); maximaPositions = theStatistics3->getMaximumVoxelPositions(); minimaPositions = theStatistics3->getMinimumVoxelPositions(); CHECK_EQUAL(maximaPositions->empty(), false); CHECK_EQUAL(minimaPositions->empty(), false); for (auto maximaPositionsIterator = std::begin(*maximaPositions); maximaPositionsIterator != std::end(*maximaPositions); ++maximaPositionsIterator) { CHECK_EQUAL(maximaPositionsIterator->first, theStatistics3->getMaximum()); } for (auto minimaPositionsIterator = std::begin(*minimaPositions); minimaPositionsIterator != std::end(*minimaPositions); ++minimaPositionsIterator) { CHECK_EQUAL(minimaPositionsIterator->first, theStatistics3->getMinimum()); } // compare with actual XML io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator(doseFilename.c_str()); core::DoseAccessorInterface::DoseAccessorPointer doseAccessorPointer(doseAccessorGenerator.generateDoseAccessor()); rttb::io::dicom::DicomFileStructureSetGenerator structAccessorGenerator(structFilename.c_str()); core::StructureSetGeneratorInterface::StructureSetPointer structerSetGeneratorPointer = structAccessorGenerator.generateStructureSet(); std::vector foundIndices = rttb::masks::VOIindexIdentifier::getIndicesByVoiRegex( structerSetGeneratorPointer, "Heart"); CHECK_EQUAL(foundIndices.size(), 1); core::MaskAccessorInterface::MaskAccessorPointer maskAccessorPointer = boost::make_shared (structerSetGeneratorPointer->getStructure(foundIndices.at(0)), doseAccessorPointer->getGeometricInfo(), true); maskAccessorPointer->updateMask(); boost::shared_ptr maskedDoseIterator = boost::make_shared(maskAccessorPointer, doseAccessorPointer); rttb::core::DoseIteratorInterface::DoseIteratorPointer doseIteratorPointer(maskedDoseIterator); rttb::algorithms::DoseStatisticsCalculator doseStatisticsCalculator(doseIteratorPointer); DoseStatisticsPointer doseStatisticsActual = doseStatisticsCalculator.calculateDoseStatistics(14.0); io::other::DoseStatisticsXMLReader readerDefaultExpected(referenceXMLFilename); auto doseStatisticsExpected = readerDefaultExpected.generateDoseStatistic(); CHECK(checkEqualDoseStatistic(doseStatisticsExpected, doseStatisticsActual)); RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/algorithms/DoseStatisticsTest.cpp b/testing/algorithms/DoseStatisticsTest.cpp index 46ead60..52a9162 100644 --- a/testing/algorithms/DoseStatisticsTest.cpp +++ b/testing/algorithms/DoseStatisticsTest.cpp @@ -1,249 +1,251 @@ // ----------------------------------------------------------------------- // 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 "rttbDoseStatistics.h" #include "rttbDataNotAvailableException.h" #include "rttbVolumeToDoseMeasure.h" namespace rttb { namespace testing { typedef rttb::algorithms::DoseStatistics::ResultListPointer ResultListPointer; - typedef rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType DoseToVolumeFunctionType; - typedef rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType VolumeToDoseFunctionType; /*! @brief DoseStatisticsTest - test the API of DoseStatistics 1) test constructors 2) test setters 3) test getters of complex statistics (with stored key and without stored key) */ int DoseStatisticsTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; DoseStatisticType minimum = 1.0; DoseStatisticType mean = 5.5; DoseStatisticType maximum = 108.2; DoseStatisticType stdDeviation = 10.1; unsigned int numVoxels = 100000; VolumeType volume = numVoxels * (0.5 * 0.5 * 0.5); std::vector > minVoxels; std::vector > maxVoxels; minVoxels.push_back(std::make_pair(1.0, 11)); minVoxels.push_back(std::make_pair(1.0, 22)); minVoxels.push_back(std::make_pair(1.0, 33)); minVoxels.push_back(std::make_pair(1.0, 44)); maxVoxels.push_back(std::make_pair(108.2, 5)); maxVoxels.push_back(std::make_pair(108.2, 6)); maxVoxels.push_back(std::make_pair(108.2, 7)); maxVoxels.push_back(std::make_pair(108.2, 8)); ResultListPointer resultsMinVoxels = boost::make_shared > >(minVoxels); ResultListPointer resultsMaxVoxels = boost::make_shared > >(maxVoxels); - DoseToVolumeFunctionType Vx; - Vx.insert(std::make_pair(1.1, 1000)); - Vx.insert(std::make_pair(106.9, 99000)); + algorithms::DoseToVolumeMeasure Vx = algorithms::DoseToVolumeMeasure(algorithms::DoseToVolumeMeasure::Vx, std::map(), maximum); + Vx.insertValue(std::make_pair(1.1, 1000)); + Vx.insertValue(std::make_pair(106.9, 99000)); + algorithms::VolumeToDoseMeasure Dx = algorithms::VolumeToDoseMeasure(algorithms::VolumeToDoseMeasure::Dx, std::map(), volume); Dx.insertValue(std::make_pair(1000, 1.1)); Dx.insertValue(std::make_pair(99000, 106.9)); algorithms::VolumeToDoseMeasure MOHx = algorithms::VolumeToDoseMeasure(algorithms::VolumeToDoseMeasure::MOHx, std::map(), volume); MOHx.insertValue(std::make_pair(1000, 5)); MOHx.insertValue(std::make_pair(99000, 105.5)); algorithms::VolumeToDoseMeasure MOCx = algorithms::VolumeToDoseMeasure(algorithms::VolumeToDoseMeasure::MOCx, std::map(), volume); MOCx.insertValue(std::make_pair(1000, 10)); MOCx.insertValue(std::make_pair(99000, 99)); algorithms::VolumeToDoseMeasure MaxOHx = algorithms::VolumeToDoseMeasure(algorithms::VolumeToDoseMeasure::MaxOHx, std::map(), volume); MaxOHx.insertValue(std::make_pair(1000, 40)); MaxOHx.insertValue(std::make_pair(99000, 98.3)); algorithms::VolumeToDoseMeasure MinOCx = algorithms::VolumeToDoseMeasure(algorithms::VolumeToDoseMeasure::MinOCx, std::map(), volume); MinOCx.insertValue(std::make_pair(1000, 25.5)); MinOCx.insertValue(std::make_pair(99000, 102.7)); //1) test constructors CHECK_NO_THROW(rttb::algorithms::DoseStatistics aDoseStatistic(minimum, maximum, mean, stdDeviation, numVoxels, volume)); rttb::algorithms::DoseStatistics aDoseStatistic(minimum, maximum, mean, stdDeviation, numVoxels, volume); CHECK_EQUAL(aDoseStatistic.getMinimum(), minimum); CHECK_EQUAL(aDoseStatistic.getMaximum(), maximum); CHECK_EQUAL(aDoseStatistic.getMean(), mean); CHECK_EQUAL(aDoseStatistic.getStdDeviation(), stdDeviation); CHECK_EQUAL(aDoseStatistic.getVariance(), stdDeviation * stdDeviation); CHECK_EQUAL(aDoseStatistic.getNumberOfVoxels(), numVoxels); CHECK_EQUAL(aDoseStatistic.getVolume(), volume); //check default values for unset complex values CHECK_EQUAL(aDoseStatistic.getMaximumVoxelPositions()->empty(), true); CHECK_EQUAL(aDoseStatistic.getMinimumVoxelPositions()->empty(), true); CHECK_EQUAL(aDoseStatistic.getDx().getAllValues().empty(), true); - CHECK_EQUAL(aDoseStatistic.getAllVx().empty(), true); + CHECK_EQUAL(aDoseStatistic.getVx().getAllValues().empty(), true); CHECK_EQUAL(aDoseStatistic.getMOHx().getAllValues().empty(), true); CHECK_EQUAL(aDoseStatistic.getMOCx().getAllValues().empty(), true); CHECK_EQUAL(aDoseStatistic.getMaxOHx().getAllValues().empty(), true); CHECK_EQUAL(aDoseStatistic.getMinOCx().getAllValues().empty(), true); CHECK_NO_THROW(rttb::algorithms::DoseStatistics aDoseStatisticComplex(minimum, maximum, mean, stdDeviation, numVoxels, volume, resultsMaxVoxels, resultsMinVoxels, Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx)); rttb::algorithms::DoseStatistics aDoseStatisticComplex(minimum, maximum, mean, stdDeviation, numVoxels, volume, resultsMaxVoxels, resultsMinVoxels, Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx); CHECK_EQUAL(aDoseStatisticComplex.getMaximumVoxelPositions(), resultsMaxVoxels); CHECK_EQUAL(aDoseStatisticComplex.getMinimumVoxelPositions(), resultsMinVoxels); CHECK_EQUAL(aDoseStatisticComplex.getDx() == Dx, true); - CHECK_EQUAL(aDoseStatisticComplex.getAllVx() == Vx, true); + CHECK_EQUAL(aDoseStatisticComplex.getVx() == Vx, true); CHECK_EQUAL(aDoseStatisticComplex.getMOHx() == MOHx, true); CHECK_EQUAL(aDoseStatisticComplex.getMOCx() == MOCx, true); CHECK_EQUAL(aDoseStatisticComplex.getMaxOHx() == MaxOHx, true); CHECK_EQUAL(aDoseStatisticComplex.getMinOCx() == MinOCx, true); //2) test setters (only complex statistics have setters) CHECK_NO_THROW(aDoseStatistic.setMaximumVoxelPositions(resultsMaxVoxels)); CHECK_NO_THROW(aDoseStatistic.setMinimumVoxelPositions(resultsMinVoxels)); CHECK_NO_THROW(aDoseStatistic.setDx(Dx)); CHECK_NO_THROW(aDoseStatistic.setVx(Vx)); CHECK_NO_THROW(aDoseStatistic.setMOHx(MOHx)); CHECK_NO_THROW(aDoseStatistic.setMOCx(MOCx)); CHECK_NO_THROW(aDoseStatistic.setMaxOHx(MaxOHx)); CHECK_NO_THROW(aDoseStatistic.setMinOCx(MinOCx)); CHECK_EQUAL(aDoseStatistic.getMaximumVoxelPositions(), resultsMaxVoxels); CHECK_EQUAL(aDoseStatistic.getMinimumVoxelPositions(), resultsMinVoxels); CHECK_EQUAL(aDoseStatistic.getDx() == Dx, true); - CHECK_EQUAL(aDoseStatistic.getAllVx() == Vx, true); + CHECK_EQUAL(aDoseStatistic.getVx() == Vx, true); CHECK_EQUAL(aDoseStatistic.getMOHx() == MOHx, true); CHECK_EQUAL(aDoseStatistic.getMOCx() == MOCx, true); CHECK_EQUAL(aDoseStatistic.getMaxOHx() == MaxOHx, true); CHECK_EQUAL(aDoseStatistic.getMinOCx() == MinOCx, true); //3) test getters of complex statistics(with stored key and without stored key) //getAll*() already tested in (2) - Vx.clear(); - Vx.insert(std::make_pair(1.1, 1000)); - Vx.insert(std::make_pair(5.0, 2300)); - Vx.insert(std::make_pair(90, 90500)); - Vx.insert(std::make_pair(107, 99000)); + Vx = algorithms::DoseToVolumeMeasure(algorithms::DoseToVolumeMeasure::Vx, std::map(), maximum); + Vx.insertValue(std::make_pair(1.1, 1000)); + Vx.insertValue(std::make_pair(5.0, 2300)); + Vx.insertValue(std::make_pair(90, 90500)); + Vx.insertValue(std::make_pair(107, 99000)); + Dx = algorithms::VolumeToDoseMeasure(algorithms::VolumeToDoseMeasure::Dx, std::map(), volume); + Dx.insertValue(std::make_pair(1000, 1.1)); Dx.insertValue(std::make_pair(2000, 2.0)); Dx.insertValue(std::make_pair(5000, 10.8)); Dx.insertValue(std::make_pair(90000, 89.5)); Dx.insertValue(std::make_pair(98000, 104.4)); + Dx.insertValue(std::make_pair(99000, 106.9)); rttb::algorithms::DoseStatistics aDoseStatisticNewValues(minimum, maximum, mean, stdDeviation, numVoxels, volume); aDoseStatisticNewValues.setDx(Dx); aDoseStatisticNewValues.setVx(Vx); - CHECK_NO_THROW(aDoseStatisticNewValues.getVx(1.1)); - CHECK_NO_THROW(aDoseStatisticNewValues.getVx(90)); + CHECK_NO_THROW(aDoseStatisticNewValues.getVx().getValue(1.1)); + CHECK_NO_THROW(aDoseStatisticNewValues.getVx().getValue(90)); CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValue(1000)); CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValue(98000)); - CHECK_EQUAL(aDoseStatisticNewValues.getVx(1.1), Vx.find(1.1)->second); - CHECK_EQUAL(aDoseStatisticNewValues.getVx(90), Vx.find(90)->second); + CHECK_EQUAL(aDoseStatisticNewValues.getVx().getValue(1.1), Vx.getAllValues().find(1.1)->second); + CHECK_EQUAL(aDoseStatisticNewValues.getVx().getValue(90), Vx.getAllValues().find(90)->second); CHECK_EQUAL(aDoseStatisticNewValues.getDx().getValue(1000), Dx.getAllValues().find(1000)->second); CHECK_EQUAL(aDoseStatisticNewValues.getDx().getValue(98000), Dx.getAllValues().find(98000)->second); //test if key-value combination NOT in map CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getDx().getValue(1001), core::DataNotAvailableException); - CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getVx(10), core::DataNotAvailableException); + CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getVx().getValue(10), core::DataNotAvailableException); double closestDxKey, closestVxKey; CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValue(900, true, closestDxKey)); CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValue(99001, true, closestDxKey)); - CHECK_NO_THROW(aDoseStatisticNewValues.getVx(10, true, closestVxKey)); + CHECK_NO_THROW(aDoseStatisticNewValues.getVx().getValue(10, true, closestVxKey)); CHECK_EQUAL(aDoseStatisticNewValues.getDx().getValue(900, true, closestDxKey), Dx.getAllValues().find(1000)->second); CHECK_EQUAL(aDoseStatisticNewValues.getDx().getValue(99001, true, closestDxKey), Dx.getAllValues().find(99000)->second); - CHECK_EQUAL(aDoseStatisticNewValues.getVx(10, true, closestVxKey), Vx.find(5.0)->second); + CHECK_EQUAL(aDoseStatisticNewValues.getVx().getValue(10, true, closestVxKey), Vx.getAllValues().find(5.0)->second); CHECK_EQUAL(closestDxKey, 99000); CHECK_EQUAL(closestVxKey, 5); // relatives only between 0 and 1 - CHECK_NO_THROW(aDoseStatisticNewValues.getVxRelative(1.1 / aDoseStatistic.getReferenceDose())); + CHECK_NO_THROW(aDoseStatisticNewValues.getVx().getValueRelative(1.1 / aDoseStatistic.getReferenceDose())); CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValueRelative(1000 / aDoseStatistic.getVolume())); - CHECK_THROW(aDoseStatisticNewValues.getVxRelative(-0.3)); - CHECK_THROW(aDoseStatisticNewValues.getVxRelative(1.1)); + CHECK_THROW(aDoseStatisticNewValues.getVx().getValueRelative(-0.3)); + CHECK_THROW(aDoseStatisticNewValues.getVx().getValueRelative(1.1)); CHECK_THROW(aDoseStatisticNewValues.getDx().getValueRelative(0.5)); CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValueRelative(900 / aDoseStatistic.getVolume(), true, closestDxKey)); CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValueRelative(0.5, true, closestDxKey)); - CHECK_NO_THROW(aDoseStatisticNewValues.getVxRelative(10 / aDoseStatistic.getReferenceDose(), true, closestVxKey)); + CHECK_NO_THROW(aDoseStatisticNewValues.getVx().getValueRelative(10 / aDoseStatistic.getReferenceDose(), true, closestVxKey)); CHECK_EQUAL(aDoseStatisticNewValues.getDx().getValueRelative(900 / aDoseStatistic.getVolume(), true, closestDxKey), Dx.getAllValues().find(1000)->second); - CHECK_EQUAL(aDoseStatisticNewValues.getVxRelative(10 / aDoseStatistic.getReferenceDose(), true, closestVxKey), Vx.find(5.0)->second); + CHECK_EQUAL(aDoseStatisticNewValues.getVx().getValueRelative(10 / aDoseStatistic.getReferenceDose(), true, closestVxKey), Vx.getAllValues().find(5.0)->second); CHECK_EQUAL(closestVxKey, 5); //equal distance to two values. First value is returned. CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValue(1500, true, closestDxKey)); - CHECK_NO_THROW(aDoseStatisticNewValues.getVx(98.5, true, closestVxKey)); + CHECK_NO_THROW(aDoseStatisticNewValues.getVx().getValue(98.5, true, closestVxKey)); CHECK_EQUAL(aDoseStatisticNewValues.getDx().getValue(1500, true, closestDxKey), Dx.getAllValues().find(1000)->second); - CHECK_EQUAL(aDoseStatisticNewValues.getVx(98.5, true, closestVxKey), Vx.find(90.0)->second); + CHECK_EQUAL(aDoseStatisticNewValues.getVx().getValue(98.5, true, closestVxKey), Vx.getAllValues().find(90.0)->second); CHECK_EQUAL(closestDxKey, 1000); CHECK_EQUAL(closestVxKey, 90.0); double dummy; CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCx().getValue(25), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHx().getValue(9999), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCx().getValue(25, true, dummy), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHx().getValue(9999, true, dummy), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCx().getValueRelative(25 / aDoseStatistic.getVolume()), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHx().getValueRelative(9999 / aDoseStatistic.getVolume()), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCx().getValueRelative(25 / aDoseStatistic.getVolume(), true, dummy), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHx().getValueRelative(9999 / aDoseStatistic.getVolume(), true, dummy), core::DataNotAvailableException); RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/examples/RTDoseStatisticsDicomTest.cpp b/testing/examples/RTDoseStatisticsDicomTest.cpp index 7847b63..f052b0f 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)); + DoseTypeGy vx = doseStatistics->getVx().getValue(expectedVal, true, dummy); + CHECK_EQUAL(vx, doseStatistics->getVx().getValue(0)); 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().getValue(24120), doseStatistics->getMinimum(), 0.001); CHECK_CLOSE(doseStatistics->getMOHx().getValue(24120), doseStatistics->getMean(), reducedErrorConstant); CHECK_CLOSE(doseStatistics->getMOCx().getValue(20000), doseStatistics->getMean(), reducedErrorConstant); CHECK_CLOSE(doseStatistics->getMinOCx().getValue(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 32014f5..fc56bae 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->getDx().getAllValues(), aDoseStatistic2->getDx().getAllValues()); - result = result && mapCompare(aDoseStatistc1->getAllVx(), aDoseStatistic2->getAllVx()); + result = result && mapCompare(aDoseStatistc1->getVx().getAllValues(), aDoseStatistic2->getVx().getAllValues()); result = result && mapCompare(aDoseStatistc1->getMaxOHx().getAllValues(), aDoseStatistic2->getMaxOHx().getAllValues()); result = result && mapCompare(aDoseStatistc1->getMinOCx().getAllValues(), aDoseStatistic2->getMinOCx().getAllValues()); result = result && mapCompare(aDoseStatistc1->getMOCx().getAllValues(), aDoseStatistic2->getMOCx().getAllValues()); result = result && mapCompare(aDoseStatistc1->getMOHx().getAllValues(), aDoseStatistic2->getMOHx().getAllValues()); 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