diff --git a/code/algorithms/files.cmake b/code/algorithms/files.cmake index 3577b8e..6dc817d 100644 --- a/code/algorithms/files.cmake +++ b/code/algorithms/files.cmake @@ -1,18 +1,24 @@ SET(CPP_FILES rttbDoseStatistics.cpp rttbDoseStatisticsCalculator.cpp rttbArithmetic.cpp + rttbDxVolumeToDoseMeasureCalculator.cpp + rttbVolumeToDoseMeasureCalculator.cpp + rttbVolumeToDoseMeasure.cpp ) SET(H_FILES rttbDoseStatistics.h rttbDoseStatisticsCalculator.h rttbArithmetic.h rttbBinaryFunctorAccessor.h + rttbDxVolumeToDoseMeasureCalculator.h + rttbVolumeToDoseMeasureCalculator.h + rttbVolumeToDoseMeasure.h ) SET(TXX_FILES rttbArithmetic.tpp rttbBinaryFunctorAccessor.tpp ) diff --git a/code/algorithms/rttbDoseStatistics.cpp b/code/algorithms/rttbDoseStatistics.cpp index 40e5aed..ee67021 100644 --- a/code/algorithms/rttbDoseStatistics.cpp +++ b/code/algorithms/rttbDoseStatistics.cpp @@ -1,507 +1,471 @@ // ----------------------------------------------------------------------- // 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()*/, - VolumeToDoseFunctionType Dx /*= std::map()*/, + VolumeToDoseMeasure Dx /*= std::map()*/, DoseToVolumeFunctionType Vx /*= std::map()*/, VolumeToDoseFunctionType MOHx /*= std::map()*/, VolumeToDoseFunctionType MOCx /*= std::map()*/, VolumeToDoseFunctionType MaxOHx /*= std::map()*/, VolumeToDoseFunctionType MinOCx /*= std::map()*/, 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 DoseToVolumeFunctionType& DxValues) + void DoseStatistics::setDx(const VolumeToDoseMeasure& DxValues) { _Dx = DxValues; } void DoseStatistics::setVx(const VolumeToDoseFunctionType& VxValues) { _Vx = VxValues; } void DoseStatistics::setMOHx(const VolumeToDoseFunctionType& MOHxValues) { _MOHx = MOHxValues; } void DoseStatistics::setMOCx(const VolumeToDoseFunctionType& MOCxValues) { _MOCx = MOCxValues; } void DoseStatistics::setMaxOHx(const VolumeToDoseFunctionType& MaxOHValues) { _MaxOHx = MaxOHValues; } void DoseStatistics::setMinOCx(const VolumeToDoseFunctionType& 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"); } } - DoseTypeGy DoseStatistics::getDx(VolumeType xVolumeAbsolute, bool findNearestValue, - VolumeType& nearestXVolume) const - { - - return getValue(_Dx, xVolumeAbsolute, findNearestValue, nearestXVolume); - } - - DoseTypeGy DoseStatistics::getDx(VolumeType xVolumeAbsolute) const - { - VolumeType dummy; - return getValue(_Dx, xVolumeAbsolute, false, dummy); - } - DoseTypeGy DoseStatistics::getDxRelative(VolumeType xVolumeRelative, bool findNearestValue, - VolumeType& nearestXVolume) const - { - if ( xVolumeRelative>=0 && xVolumeRelative <= 1){ - DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; - return getValue(_Dx, xVolumeAbsolute, findNearestValue, nearestXVolume); - } - else { - throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); - } - } - - DoseTypeGy DoseStatistics::getDxRelative(VolumeType xVolumeRelative) const - { - if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ - DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; - VolumeType dummy; - return getValue(_Dx, xVolumeAbsolute, false, dummy); - } - else { - throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); - } - } - DoseTypeGy DoseStatistics::getMOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const { return getValue(_MOHx, xVolumeAbsolute, findNearestValue, nearestXVolume); } DoseTypeGy DoseStatistics::getMOHx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_MOHx, xVolumeAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getMOHxRelative(VolumeType xVolumeRelative, bool findNearestValue, VolumeType& nearestXVolume) const { if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; return getValue(_MOHx, xVolumeAbsolute, findNearestValue, nearestXVolume); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } DoseTypeGy DoseStatistics::getMOHxRelative(VolumeType xVolumeRelative) const { if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; VolumeType dummy; return getValue(_MOHx, xVolumeAbsolute, false, dummy); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } DoseTypeGy DoseStatistics::getMOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { return getValue(_MOCx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } DoseTypeGy DoseStatistics::getMOCx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_MOCx, xVolumeAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getMOCxRelative(VolumeType xVolumeRelative, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; return getValue(_MOCx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } DoseTypeGy DoseStatistics::getMOCxRelative(VolumeType xVolumeRelative) const { if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; VolumeType dummy; return getValue(_MOCx, xVolumeAbsolute, false, dummy); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } DoseTypeGy DoseStatistics::getMaxOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { return getValue(_MaxOHx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } DoseTypeGy DoseStatistics::getMaxOHx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_MaxOHx, xVolumeAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getMaxOHxRelative(VolumeType xVolumeRelative, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; return getValue(_MaxOHx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } DoseTypeGy DoseStatistics::getMaxOHxRelative(VolumeType xVolumeRelative) const { if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; VolumeType dummy; return getValue(_MaxOHx, xVolumeAbsolute, false, dummy); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } DoseTypeGy DoseStatistics::getMinOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { return getValue(_MinOCx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } DoseTypeGy DoseStatistics::getMinOCx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_MinOCx, xVolumeAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getMinOCxRelative(VolumeType xVolumeRelative, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; return getValue(_MinOCx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } } DoseTypeGy DoseStatistics::getMinOCxRelative(VolumeType xVolumeRelative) const { if (xVolumeRelative >= 0 && xVolumeRelative <= 1){ DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; VolumeType dummy; return getValue(_MinOCx, xVolumeAbsolute, false, dummy); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=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 { return _Vx; } - DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllDx() const + VolumeToDoseMeasure DoseStatistics::getDx() const { return _Dx; } DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllMOHx() const { return _MOHx; } DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllMOCx() const { return _MOCx; } DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllMaxOHx() const { return _MaxOHx; } DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllMinOCx() const { return _MinOCx; } }//end namespace algorithms }//end namespace rttb diff --git a/code/algorithms/rttbDoseStatistics.h b/code/algorithms/rttbDoseStatistics.h index aa3c8ce..60a223c 100644 --- a/code/algorithms/rttbDoseStatistics.h +++ b/code/algorithms/rttbDoseStatistics.h @@ -1,237 +1,229 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __DOSE_STATISTICS_H #define __DOSE_STATISTICS_H #include #include #include "boost/shared_ptr.hpp" #include "rttbBaseType.h" #include "RTTBAlgorithmsExports.h" +#include "rttbVolumeToDoseMeasure.h" + namespace rttb { namespace algorithms { /*! @class DoseStatistics @brief This is a data class storing different statistical values from a rt dose distribution @sa DoseStatisticsCalculator */ class RTTBAlgorithms_EXPORT DoseStatistics { public: enum complexStatistics { Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx }; typedef boost::shared_ptr > > ResultListPointer; typedef boost::shared_ptr DoseStatisticsPointer; typedef std::map DoseToVolumeFunctionType; typedef std::map VolumeToDoseFunctionType; private: double getValue(const std::map& aMap, double key, bool findNearestValueInstead, double& storedKey) const; std::map::const_iterator findNearestKeyInMap(const std::map& aMap, double key) const; DoseStatisticType _maximum; DoseStatisticType _minimum; ResultListPointer _maximumVoxelPositions; ResultListPointer _minimumVoxelPositions; DoseStatisticType _mean; DoseStatisticType _stdDeviation; VoxelNumberType _numVoxels; VolumeType _volume; DoseTypeGy _referenceDose; //for Vx computation - VolumeToDoseFunctionType _Dx; + VolumeToDoseMeasure _Dx; DoseToVolumeFunctionType _Vx; VolumeToDoseFunctionType _MOHx; VolumeToDoseFunctionType _MOCx; VolumeToDoseFunctionType _MaxOHx; VolumeToDoseFunctionType _MinOCx; public: /*! @brief Standard Constructor */ //DoseStatistics(); /*! @brief Constructor @details the dose statistic values are set. Complex values maximumVoxelLocation, maximumVoxelLocation, Dx, Vx, MOHx, MOCx, MaxOHx and MinOCx are optional */ DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean, DoseStatisticType stdDeviation, VoxelNumberType numVoxels, VolumeType volume, ResultListPointer minimumVoxelPositions = NULL, ResultListPointer maximumVoxelPositions = NULL, - VolumeToDoseFunctionType Dx = VolumeToDoseFunctionType(), + VolumeToDoseMeasure _Dx = VolumeToDoseMeasure("Dx", std::map(), 0.0), DoseToVolumeFunctionType Vx = DoseToVolumeFunctionType(), VolumeToDoseFunctionType MOHx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MOCx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MaxOHx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MinOCx = VolumeToDoseFunctionType(), DoseTypeGy referenceDose = -1); ~DoseStatistics(); void setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions); void setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions); - void setDx(const DoseToVolumeFunctionType& DxValues); + void setDx(const VolumeToDoseMeasure& DxValues); void setVx(const VolumeToDoseFunctionType& VxValues); void setMOHx(const VolumeToDoseFunctionType& MOHxValues); void setMOCx(const VolumeToDoseFunctionType& MOCxValues); void setMaxOHx(const VolumeToDoseFunctionType& MaxOHxValues); void setMinOCx(const VolumeToDoseFunctionType& MinOCxValues); void setReferenceDose(DoseTypeGy referenceDose); /*! @brief Get number of voxels in doseIterator, with sub-voxel accuracy. */ VoxelNumberType getNumberOfVoxels() const; /*! @brief Get the volume of the voxels in doseIterator (in cm3), with sub-voxel accuracy */ VolumeType getVolume() const; /*! @brief Get the reference dose for Vx computation */ DoseTypeGy getReferenceDose() const; /*! @brief Get the maximum of the current dose distribution. @return Return the maximum dose in Gy */ DoseStatisticType getMaximum() const; /*! @brief Get a vector of the the maximum dose VoxelGridIDs together with their dose value in Gy @exception InvalidDoseException if the vector has not been set (i.e. is empty) */ ResultListPointer getMaximumVoxelPositions() const; /*! @brief Get the minimum of the current dose distribution. @return Return the minimum dose in Gy */ DoseStatisticType getMinimum() const; /*! @brief Get a vector of the the minimum dose VoxelGridIDs together with their dose value in Gy @exception InvalidDoseException if the vector has not been set (i.e. is empty) */ ResultListPointer getMinimumVoxelPositions() const; /*! @brief Get the mean of the current dose distribution. @return Return the mean dose in Gy */ DoseStatisticType getMean() const; /*! @brief Get the standard deviation of the current dose distribution. @return Return the standard deviation in Gy */ DoseStatisticType getStdDeviation() const; /*! @brief Get the variance of of the current dose distribution. @return Return the variance in Gy */ DoseStatisticType getVariance() const; /*! @brief Get Vx: the volume irradiated with a dose >= x. @return Return absolute volume in absolute cm^3. @exception NoDataException if the Vx values have not been set (i.e. the vector is empty) @exception NoDataException if the requested Dose is not in the vector */ VolumeType getVx(DoseTypeGy xDoseAbsolute) const; VolumeType getVx(DoseTypeGy xDoseAbsolute, bool findNearestValue, DoseTypeGy& nearestXDose) const; VolumeType getVxRelative(DoseTypeGy xDoseRelative) const; VolumeType getVxRelative(DoseTypeGy xDoseRelative, bool findNearestValue, DoseTypeGy& nearestXDose) const; DoseToVolumeFunctionType getAllVx() const; - /*! @brief Get Dx: the minimal dose delivered to part x of the current volume. - @return Return dose value in Gy. - @exception InvalidDoseException if the Dx values have not been set (i.e. the vector is empty) - */ - DoseTypeGy getDx(VolumeType xVolumeAbsolute, bool findNearestValue, - VolumeType& nearestXVolume) const; - DoseTypeGy getDx(VolumeType xVolumeAbsolute) const; - DoseTypeGy getDxRelative(VolumeType xDoseRelative, bool findNearestValue, - VolumeType& nearestXVolume) const; - DoseTypeGy getDxRelative(VolumeType xDoseRelative) const; - VolumeToDoseFunctionType getAllDx() const; + VolumeToDoseMeasure getDx() const; /*! @brief Get MOHx: mean dose of the hottest x voxels. @return Return dose value in Gy. @exception InvalidDoseException if the values have not been set (i.e. the vector is empty) */ DoseTypeGy getMOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOHx(VolumeType xVolumeAbsolute) const; DoseTypeGy getMOHxRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOHxRelative(VolumeType xDoseRelative) const; VolumeToDoseFunctionType getAllMOHx() const; /*! @brief Get MOCx: mean dose of the coldest x voxels. @return Return dose value in Gy. @exception InvalidDoseException if the values have not been set (i.e. the vector is empty) */ DoseTypeGy getMOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOCx(VolumeType xVolumeAbsolute) const; DoseTypeGy getMOCxRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOCxRelative(VolumeType xDoseRelative) const; VolumeToDoseFunctionType getAllMOCx() const; /*! @brief Get MaxOHx: Maximum outside of the hottest x voxels. @return Return dose value in Gy. @exception InvalidDoseException if the values have not been set (i.e. the vector is empty) */ DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute) const; DoseTypeGy getMaxOHxRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXVolum) const; DoseTypeGy getMaxOHxRelative(VolumeType xDoseRelative) const; VolumeToDoseFunctionType getAllMaxOHx() const; /*! @brief Get MinOCx: Minimum outside of the coldest x voxels. @return Return dose value in Gy. @exception InvalidDoseException if the values have not been set (i.e. the vector is empty) */ DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute) const; DoseTypeGy getMinOCxRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMinOCxRelative(VolumeType xDoseRelative) const; VolumeToDoseFunctionType getAllMinOCx() const; }; } } #endif diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp index 5a90da5..0251167 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,683 +1,680 @@ // ----------------------------------------------------------------------- // 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 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 test = DxVolumeToDoseMeasureCalculator(precomputeVolumeValuesNonConst, _statistics->getVolume(), + this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), _statistics->getMinimum()); + test.compute(); DoseToVolumeFunctionType Vx = computeDoseToVolumeFunctionMulti(referenceDose, precomputeDoseValuesNonConst, DoseStatistics::Vx); - VolumeToDoseFunctionType Dx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, - DoseStatistics::Dx); VolumeToDoseFunctionType MOHx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MOHx); VolumeToDoseFunctionType MOCx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MOCx); VolumeToDoseFunctionType MaxOHx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MaxOHx); VolumeToDoseFunctionType MinOCx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MinOCx); _statistics->setVx(Vx); - _statistics->setDx(Dx); + _statistics->setDx(test.getMeasure()); _statistics->setMOHx(MOHx); _statistics->setMOCx(MOCx); _statistics->setMaxOHx(MaxOHx); _statistics->setMinOCx(MinOCx); _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(); } DoseTypeGy DoseStatisticsCalculator::computeDx(VolumeType xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; bool voxelOverflow = false; for (auto i = _doseVector.size() - 1; i != -1; i--) { countVoxels += _voxelProportionVector.at(i); if (countVoxels >= noOfVoxel) { voxelOverflow = true; resultDose = _doseVector.at(i); break; } } if (!voxelOverflow) { resultDose = _statistics->getMinimum(); } return resultDose; } DoseTypeGy DoseStatisticsCalculator::computeMOHx(VolumeType xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); if (noOfVoxel == 0) { return 0; } else { double countVoxels = 0; double sum = 0; for (auto i = _doseVector.size() - 1; i != -1; i--) { double voxelProportion = _voxelProportionVector.at(i); countVoxels += voxelProportion; sum += _doseVector.at(i) * voxelProportion; if (countVoxels >= noOfVoxel) { break; } } return (DoseTypeGy)(sum / noOfVoxel); } } DoseTypeGy DoseStatisticsCalculator::computeMOCx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); if (noOfVoxel == 0) { return 0; } else { double countVoxels = 0; double sum = 0; std::vector::const_iterator it = _doseVector.begin(); std::vector::const_iterator itD = _voxelProportionVector.begin(); for (; it != _doseVector.end(); ++it, ++itD) { double voxelProportion = *itD; countVoxels += voxelProportion; sum += (*it) * voxelProportion; if (countVoxels >= noOfVoxel) { break; } } return (DoseTypeGy)(sum / noOfVoxel); } } DoseTypeGy DoseStatisticsCalculator::computeMaxOHx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; for (auto i = _doseVector.size() - 1; i != -1; i--) { countVoxels += _voxelProportionVector.at(i); if (countVoxels >= noOfVoxel) { if (i > 0) { resultDose = _doseVector.at(i - 1); } break; } } return resultDose; } DoseTypeGy DoseStatisticsCalculator::computeMinOCx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; std::vector::const_iterator it = _doseVector.begin(); std::vector::const_iterator itD = _voxelProportionVector.begin(); for (; itD != _voxelProportionVector.end(); ++itD, ++it) { countVoxels += *itD; if (countVoxels >= noOfVoxel) { break; } } if (it != _doseVector.end()) { ++it; if (it != _doseVector.end()) { resultDose = *it; } else { resultDose = (DoseTypeGy)_statistics->getMaximum(); } } else { resultDose = (DoseTypeGy)_statistics->getMinimum(); } return resultDose; } 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!"); } } DoseStatisticsCalculator::VolumeToDoseFunctionType DoseStatisticsCalculator::computeVolumeToDoseFunctionMulti( const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const { std::vector threads; VolumeToDoseFunctionType multiValues; VolumeType volume = _statistics->getVolume(); for (size_t i = 0; i < precomputeVolumeValues.size(); ++i) { if (_multiThreading) { threads.push_back(boost::thread(&DoseStatisticsCalculator::computeVolumeToDoseSingle, this, precomputeVolumeValues.at(i), name, std::ref(multiValues), volume)); } else { DoseStatisticsCalculator::computeVolumeToDoseSingle(precomputeVolumeValues.at(i), name, std::ref(multiValues), volume); } } for (unsigned int i=0; i lock(*_mutex); switch (name) { - case DoseStatistics::Dx: - multiValues.insert(std::pair(xAbsolute, - computeDx(xAbsolute))); - break; - case DoseStatistics::MOHx: multiValues.insert(std::pair(xAbsolute, computeMOHx(xAbsolute))); break; case DoseStatistics::MOCx: multiValues.insert(std::pair(xAbsolute, computeMOCx(xAbsolute))); break; case DoseStatistics::MaxOHx: multiValues.insert(std::pair(xAbsolute, computeMaxOHx(xAbsolute))); break; case DoseStatistics::MinOCx: multiValues.insert(std::pair(xAbsolute, computeMinOCx(xAbsolute))); break; default: 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/rttbDxVolumeToDoseMeasureCalculator.h b/code/algorithms/rttbDxVolumeToDoseMeasureCalculator.h index a3dfbbf..d2780d9 100644 --- a/code/algorithms/rttbDxVolumeToDoseMeasureCalculator.h +++ b/code/algorithms/rttbDxVolumeToDoseMeasureCalculator.h @@ -1,58 +1,57 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 1674 $ (last changed revision) // @date $Date: 2017-01-27 10:34:46 +0100 (Fr, 27 Jan 2017) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __DX_VOLUME_TO_DOSE_MEASURE_CALCULATOR_H #define __DX_VOLUME_TO_DOSE_MEASURE_CALCULATOR_H #include #include #include "rttbBaseType.h" #include "RTTBAlgorithmsExports.h" #include "rttbvolumetodosemeasureCalculator.h" #include "rttbvolumetodosemeasure.h" namespace rttb { namespace algorithms { class RTTBAlgorithms_EXPORT DxVolumeToDoseMeasureCalculator: public VolumeToDoseMeasureCalculator { private: DoseStatisticType _minimum; public: DxVolumeToDoseMeasureCalculator(const std::vector& precomputeVolumeValues, const VolumeType& volume, const std::vector& doseVector, const std::vector& voxelProportionVector, const DoseVoxelVolumeType& currentVoxelVolume, const DoseStatisticType& minimum); private: void computeSpecificValue(double xAbsolute); }; } } #endif -#pragma once diff --git a/code/algorithms/rttbVolumeToDoseMeasure.cpp b/code/algorithms/rttbVolumeToDoseMeasure.cpp index 3091107..46579ef 100644 --- a/code/algorithms/rttbVolumeToDoseMeasure.cpp +++ b/code/algorithms/rttbVolumeToDoseMeasure.cpp @@ -1,127 +1,133 @@ #include "rttbVolumeToDoseMeasure.h" #include "rttbDataNotAvailableException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace algorithms { - VolumeToDoseMeasure::VolumeToDoseMeasure(DoseStatistics::complexStatistics name, VolumeToDoseFunctionType values) : - name(name), values(values) {} + VolumeToDoseMeasure::VolumeToDoseMeasure(std::string name, VolumeToDoseFunctionType values, VolumeType volume) : + name(name), values(values), _volume(volume) {} - DoseStatistics::complexStatistics VolumeToDoseMeasure::getName() + std::string VolumeToDoseMeasure::getName() const { return this->name; } void VolumeToDoseMeasure::insertValue(std::pair value) { this->values.insert(value); } DoseTypeGy VolumeToDoseMeasure::getValue(VolumeType xVolumeAbsolute) { VolumeType dummy; - return getValue(xVolumeAbsolute, false, dummy); - return DoseTypeGy(); + return getSpecificValue(xVolumeAbsolute, false, dummy); } + DoseTypeGy VolumeToDoseMeasure::getValue(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType & nearestXVolume) { - return getValue(xVolumeAbsolute, findNearestValue, nearestXVolume); - return DoseTypeGy(); + return getSpecificValue(xVolumeAbsolute, findNearestValue, nearestXVolume); } DoseTypeGy VolumeToDoseMeasure::getValueRelative(VolumeType xVolumeRelative) { if (xVolumeRelative >= 0 && xVolumeRelative <= 1) { DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; VolumeType dummy; - return getValue(xVolumeAbsolute, false, dummy); + return getSpecificValue(xVolumeAbsolute, false, dummy); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } - return DoseTypeGy(); } DoseTypeGy VolumeToDoseMeasure::getValueRelative(VolumeType xVolumeRelative, bool findNearestValue, VolumeType & nearestXVolume) { if (xVolumeRelative >= 0 && xVolumeRelative <= 1) { DoseTypeGy xVolumeAbsolute = xVolumeRelative*_volume; - return getValue(xVolumeAbsolute, findNearestValue, nearestXVolume); + return getSpecificValue(xVolumeAbsolute, findNearestValue, nearestXVolume); } else { throw rttb::core::InvalidParameterException("Relative Volume must be >= 0 and <=1"); } - return DoseTypeGy(); } - VolumeToDoseMeasure::VolumeToDoseFunctionType VolumeToDoseMeasure::getAllValues() + VolumeToDoseMeasure::VolumeToDoseFunctionType VolumeToDoseMeasure::getAllValues() const { return this->values; } - double VolumeToDoseMeasure::getValue(double key, + double VolumeToDoseMeasure::getSpecificValue(double key, bool findNearestValueInstead, double& storedKey) const { if (values.find(key) != std::end(values)) { return values.find(key)->second; } else { //value not in map. We have to find the nearest value if (values.empty()) { throw core::DataNotAvailableException("No Vx values are defined"); } else { if (findNearestValueInstead) { auto iterator = findNearestKeyInMap(values, key); storedKey = iterator->first; return iterator->second; } else { throw core::DataNotAvailableException("No Vx value with required dose is defined"); } } } } + std::map::const_iterator VolumeToDoseMeasure::findNearestKeyInMap( const std::map& aMap, double key) const { double minDistance = 1e19; double minDistanceLast = 1e20; auto iterator = std::begin(aMap); while (iterator != std::end(aMap)) { minDistanceLast = minDistance; minDistance = fabs(iterator->first - key); if (minDistanceLast > minDistance) { ++iterator; } else { if (iterator != std::begin(aMap)) { --iterator; return iterator; } else { return std::begin(aMap); } } } --iterator; return iterator; } - } + + bool operator==(const VolumeToDoseMeasure& volumeToDoseMesure,const VolumeToDoseMeasure& otherVolumeToDoseMesure) + { + if (volumeToDoseMesure.getName() == otherVolumeToDoseMesure.getName() && volumeToDoseMesure.getAllValues() == otherVolumeToDoseMesure.getAllValues()) { + return true; + } + return false; + } +} } diff --git a/code/algorithms/rttbVolumeToDoseMeasure.h b/code/algorithms/rttbVolumeToDoseMeasure.h index e14e7f3..c4bf284 100644 --- a/code/algorithms/rttbVolumeToDoseMeasure.h +++ b/code/algorithms/rttbVolumeToDoseMeasure.h @@ -1,70 +1,68 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 1674 $ (last changed revision) // @date $Date: 2017-01-27 10:34:46 +0100 (Fr, 27 Jan 2017) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __VOLUME_TO_DOSE_MEASURE_H #define __VOLUME_TO_DOSE_MEASURE_H #include #include -#include "boost/shared_ptr.hpp" - #include "rttbBaseType.h" #include "RTTBAlgorithmsExports.h" #include "rttbDoseStatistics.h" namespace rttb { namespace algorithms { class RTTBAlgorithms_EXPORT VolumeToDoseMeasure { public: typedef std::map VolumeToDoseFunctionType; private: - DoseStatistics::complexStatistics name; + std::string name; VolumeToDoseFunctionType values; + VolumeType _volume; public: - VolumeToDoseMeasure(DoseStatistics::complexStatistics name, VolumeToDoseFunctionType values); - DoseStatistics::complexStatistics getName(); + VolumeToDoseMeasure(std::string name, VolumeToDoseFunctionType values, VolumeType volume); + std::string getName() const; void insertValue(std::pair value); DoseTypeGy getValue(VolumeType xVolumeAbsolute); DoseTypeGy getValue(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXDose); DoseTypeGy getValueRelative(VolumeType xDoseRelative); DoseTypeGy getValueRelative(VolumeType xDoseRelative, bool findNearestValue, VolumeType& nearestXDose); - VolumeToDoseFunctionType getAllValues(); + VolumeToDoseFunctionType getAllValues() const; + friend bool operator==(const VolumeToDoseMeasure& volumeToDoseMesure, const VolumeToDoseMeasure& otherVolumeToDoseMesure); private: - double getValue(double key, bool findNearestValueInstead, double& storedKey) const; + double getSpecificValue(double key, bool findNearestValueInstead, double& storedKey) const; std::map::const_iterator findNearestKeyInMap(const std::map& aMap, double key) const; }; - } } #endif -#pragma once diff --git a/code/algorithms/rttbVolumeToDoseMeasureCalculator.cpp b/code/algorithms/rttbVolumeToDoseMeasureCalculator.cpp index c236107..96f6c8d 100644 --- a/code/algorithms/rttbVolumeToDoseMeasureCalculator.cpp +++ b/code/algorithms/rttbVolumeToDoseMeasureCalculator.cpp @@ -1,46 +1,46 @@ #include "rttbVolumeToDoseMeasureCalculator.h" #include //#include namespace rttb { namespace algorithms { VolumeToDoseMeasureCalculator::VolumeToDoseMeasureCalculator(const std::vector& precomputeVolumeValues, const VolumeType& volume, const std::vector& doseVector, const std::vector& voxelProportionVector, const DoseVoxelVolumeType& currentVoxelVolume) : - measure(VolumeToDoseMeasure(DoseStatistics::Dx, std::map())), _precomputeVolumeValues(precomputeVolumeValues), _volume(volume), + measure(VolumeToDoseMeasure("dx", std::map(), volume)), _precomputeVolumeValues(precomputeVolumeValues), _volume(volume), _doseVector(doseVector), _voxelProportionVector(voxelProportionVector), _currentVoxelVolume(currentVoxelVolume) {} void VolumeToDoseMeasureCalculator::compute() { std::vector threads; for (size_t i = 0; i < _precomputeVolumeValues.size(); ++i) { if (false)//_multiThreading) { threads.push_back(boost::thread(&VolumeToDoseMeasureCalculator::computeSpecificValue, this, _precomputeVolumeValues.at(i) * _volume)); } else { this->computeSpecificValue(_precomputeVolumeValues.at(i) * _volume); } } for (unsigned int i = 0; i(xAbsolute, resultDose)); } } } diff --git a/code/algorithms/rttbVolumeToDoseMeasureCalculator.h b/code/algorithms/rttbVolumeToDoseMeasureCalculator.h index 535a87d..5c06192 100644 --- a/code/algorithms/rttbVolumeToDoseMeasureCalculator.h +++ b/code/algorithms/rttbVolumeToDoseMeasureCalculator.h @@ -1,70 +1,69 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 1674 $ (last changed revision) // @date $Date: 2017-01-27 10:34:46 +0100 (Fr, 27 Jan 2017) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __VOLUME_TO_DOSE_MEASURE_CALCULATOR_H #define __VOLUME_TO_DOSE_MEASURE_CALCULATOR_H #include #include #include "rttbBaseType.h" #include "RTTBAlgorithmsExports.h" #include "rttbVolumeToDoseMeasure.h" namespace rttb { namespace algorithms { class RTTBAlgorithms_EXPORT VolumeToDoseMeasureCalculator { public: typedef DoseStatistics::VolumeToDoseFunctionType VolumeToDoseFunctionType; protected: std::vector _doseVector; DoseVoxelVolumeType _currentVoxelVolume; std::vector _voxelProportionVector; private: VolumeType _volume; VolumeToDoseMeasure measure; std::vector _precomputeVolumeValues; public: VolumeToDoseMeasureCalculator(const std::vector& precomputeVolumeValues, const VolumeType& volume, const std::vector& doseVector, const std::vector& voxelProportionVector, const DoseVoxelVolumeType& currentVoxelVolume); void compute(); VolumeToDoseMeasure getMeasure(); virtual void computeSpecificValue(double xAbsolute) = 0; protected: void insertIntoMeasure(VolumeType xAbsolute, DoseTypeGy resultDose); }; } } #endif -#pragma once diff --git a/testing/algorithms/DoseStatisticsCalculatorTest.cpp b/testing/algorithms/DoseStatisticsCalculatorTest.cpp index 22cbd7e..36b091b 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->getAllDx().empty(), true); + CHECK_EQUAL(theStatistics->getDx().getAllValues().empty(), true); CHECK_EQUAL(theStatistics->getAllVx().empty(), true); CHECK_EQUAL(theStatistics->getAllMaxOHx().empty(), true); CHECK_EQUAL(theStatistics->getAllMOHx().empty(), true); CHECK_EQUAL(theStatistics->getAllMOCx().empty(), true); CHECK_EQUAL(theStatistics->getAllMinOCx().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->getDx(0.02 * theStatisticsDefault->getVolume())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.05 * theStatisticsDefault->getVolume())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.1 * theStatisticsDefault->getVolume())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.9 * theStatisticsDefault->getVolume())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.95 * theStatisticsDefault->getVolume())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.98 * theStatisticsDefault->getVolume())); + 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()), core::DataNotAvailableException); CHECK_NO_THROW(theStatistics->getVx(0.1 * 100.0)); - CHECK_NO_THROW(theStatistics->getDx(0.1 * theStatistics->getVolume())); - CHECK_NO_THROW(theStatistics->getDx(0.9 * theStatistics->getVolume())); + 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(0.95 * theStatistics->getVolume())); CHECK_NO_THROW(theStatistics->getMOCx(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()), core::DataNotAvailableException); - CHECK_NO_THROW(theStatistics->getDx(0.9 * theStatistics->getVolume())); - CHECK_NO_THROW(theStatistics->getDx(0.95 * theStatistics->getVolume())); - CHECK_NO_THROW(theStatistics->getDx(0.99 * theStatistics->getVolume())); - CHECK_THROW_EXPLICIT(theStatistics->getDx(0.03 * theStatistics->getVolume()), + 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->getDx(0.9 * theStatistics->getVolume()), - theStatisticsDefault->getDx(0.9 * theStatistics->getVolume())); - CHECK_EQUAL(theStatistics->getDx(0.95 * theStatistics->getVolume()), - theStatisticsDefault->getDx(0.95 * theStatistics->getVolume())); + 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()), core::DataNotAvailableException); CHECK_NO_THROW(theStatistics->getVx(0.01 * 100.0)); - CHECK_NO_THROW(theStatistics->getDx(0.9 * theStatistics->getVolume())); + 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 6d6d1ec..5ed8d60 100644 --- a/testing/algorithms/DoseStatisticsTest.cpp +++ b/testing/algorithms/DoseStatisticsTest.cpp @@ -1,251 +1,249 @@ // ----------------------------------------------------------------------- // 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)); - - VolumeToDoseFunctionType Dx; - Dx.insert(std::make_pair(1000, 1.1)); - Dx.insert(std::make_pair(99000, 106.9)); + algorithms::VolumeToDoseMeasure Dx = algorithms::VolumeToDoseMeasure("Dx", std::map(), volume); + Dx.insertValue(std::make_pair(1000, 1.1)); + Dx.insertValue(std::make_pair(99000, 106.9)); VolumeToDoseFunctionType MOHx; MOHx.insert(std::make_pair(1000, 5)); MOHx.insert(std::make_pair(99000, 105.5)); VolumeToDoseFunctionType MOCx; MOCx.insert(std::make_pair(1000, 10)); MOCx.insert(std::make_pair(99000, 99)); VolumeToDoseFunctionType MaxOHx; MaxOHx.insert(std::make_pair(1000, 40)); MaxOHx.insert(std::make_pair(99000, 98.3)); VolumeToDoseFunctionType MinOCx; MinOCx.insert(std::make_pair(1000, 25.5)); MinOCx.insert(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.getAllDx().empty(), true); + CHECK_EQUAL(aDoseStatistic.getDx().getAllValues().empty(), true); CHECK_EQUAL(aDoseStatistic.getAllVx().empty(), true); CHECK_EQUAL(aDoseStatistic.getAllMOHx().empty(), true); CHECK_EQUAL(aDoseStatistic.getAllMOCx().empty(), true); CHECK_EQUAL(aDoseStatistic.getAllMaxOHx().empty(), true); CHECK_EQUAL(aDoseStatistic.getAllMinOCx().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.getAllDx() == Dx, true); + CHECK_EQUAL(aDoseStatisticComplex.getDx().getAllValues()== Dx.getAllValues(), true); CHECK_EQUAL(aDoseStatisticComplex.getAllVx() == Vx, true); CHECK_EQUAL(aDoseStatisticComplex.getAllMOHx() == MOHx, true); CHECK_EQUAL(aDoseStatisticComplex.getAllMOCx() == MOCx, true); CHECK_EQUAL(aDoseStatisticComplex.getAllMaxOHx() == MaxOHx, true); CHECK_EQUAL(aDoseStatisticComplex.getAllMinOCx() == 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.getAllDx() == Dx, true); + CHECK_EQUAL(aDoseStatistic.getDx().getAllValues() == Dx.getAllValues(), true); CHECK_EQUAL(aDoseStatistic.getAllVx() == Vx, true); CHECK_EQUAL(aDoseStatistic.getAllMOHx() == MOHx, true); CHECK_EQUAL(aDoseStatistic.getAllMOCx() == MOCx, true); CHECK_EQUAL(aDoseStatistic.getAllMaxOHx() == MaxOHx, true); CHECK_EQUAL(aDoseStatistic.getAllMinOCx() == 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)); - Dx.clear(); - Dx.insert(std::make_pair(1000, 1.1)); - Dx.insert(std::make_pair(2000, 2.0)); - Dx.insert(std::make_pair(5000, 10.8)); - Dx.insert(std::make_pair(90000, 89.5)); - Dx.insert(std::make_pair(98000, 104.4)); - Dx.insert(std::make_pair(99000, 106.9)); + 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)); 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.getDx(1000)); - CHECK_NO_THROW(aDoseStatisticNewValues.getDx(98000)); + 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.getDx(1000), Dx.find(1000)->second); - CHECK_EQUAL(aDoseStatisticNewValues.getDx(98000), Dx.find(98000)->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(1001), core::DataNotAvailableException); + CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getDx().getValue(1001), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getVx(10), core::DataNotAvailableException); double closestDxKey, closestVxKey; - CHECK_NO_THROW(aDoseStatisticNewValues.getDx(900, true, closestDxKey)); - CHECK_NO_THROW(aDoseStatisticNewValues.getDx(99001, true, closestDxKey)); + 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_EQUAL(aDoseStatisticNewValues.getDx(900, true, closestDxKey), Dx.find(1000)->second); - CHECK_EQUAL(aDoseStatisticNewValues.getDx(99001, true, closestDxKey), Dx.find(99000)->second); + 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(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.getDxRelative(1000 / aDoseStatistic.getVolume())); + CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValueRelative(1000 / aDoseStatistic.getVolume())); CHECK_THROW(aDoseStatisticNewValues.getVxRelative(-0.3)); CHECK_THROW(aDoseStatisticNewValues.getVxRelative(1.1)); - CHECK_THROW(aDoseStatisticNewValues.getDxRelative(0.5)); + CHECK_THROW(aDoseStatisticNewValues.getDx().getValueRelative(0.5)); - CHECK_NO_THROW(aDoseStatisticNewValues.getDxRelative(900 / aDoseStatistic.getVolume(), true, closestDxKey)); - CHECK_NO_THROW(aDoseStatisticNewValues.getDxRelative(0.5, true, closestDxKey)); + 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_EQUAL(aDoseStatisticNewValues.getDxRelative(900 / aDoseStatistic.getVolume(), true, closestDxKey), Dx.find(1000)->second); + 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(closestVxKey, 5); //equal distance to two values. First value is returned. - CHECK_NO_THROW(aDoseStatisticNewValues.getDx(1500, true, closestDxKey)); + CHECK_NO_THROW(aDoseStatisticNewValues.getDx().getValue(1500, true, closestDxKey)); CHECK_NO_THROW(aDoseStatisticNewValues.getVx(98.5, true, closestVxKey)); - CHECK_EQUAL(aDoseStatisticNewValues.getDx(1500, true, closestDxKey), Dx.find(1000)->second); + 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(closestDxKey, 1000); CHECK_EQUAL(closestVxKey, 90.0); double dummy; CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCx(25), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHx(9999), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCx(25, true, dummy), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHx(9999, true, dummy), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCxRelative(25 / aDoseStatistic.getVolume()), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHxRelative(9999 / aDoseStatistic.getVolume()), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMinOCxRelative(25 / aDoseStatistic.getVolume(), true, dummy), core::DataNotAvailableException); CHECK_THROW_EXPLICIT(aDoseStatisticNewValues.getMOHxRelative(9999 / aDoseStatistic.getVolume(), true, dummy), core::DataNotAvailableException); RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb