diff --git a/code/algorithms/rttbDoseStatistics.cpp b/code/algorithms/rttbDoseStatistics.cpp index ac8f926..f6e357d 100644 --- a/code/algorithms/rttbDoseStatistics.cpp +++ b/code/algorithms/rttbDoseStatistics.cpp @@ -1,315 +1,315 @@ // ----------------------------------------------------------------------- // 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 "rttbDataNotAvailableException.h" namespace rttb { namespace algorithms { DoseStatistics::DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean, - DoseStatisticType stdDeviation, unsigned int numVoxels, VolumeType volume, + DoseStatisticType stdDeviation, VolumeType numVoxels, VolumeType volume, ResultListPointer maximumVoxelPositions /*= ResultListPointer()*/, ResultListPointer minimumVoxelPositions /*= ResultListPointer()*/, - DoseToVolumeFunctionType Dx /*= std::map()*/, - VolumeToDoseFunctionType Vx /*= std::map()*/, + VolumeToDoseFunctionType Dx /*= std::map()*/, + DoseToVolumeFunctionType Vx /*= std::map()*/, VolumeToDoseFunctionType MOHx /*= std::map()*/, VolumeToDoseFunctionType MOCx /*= std::map()*/, VolumeToDoseFunctionType MaxOHx /*= std::map()*/, VolumeToDoseFunctionType MinOCx /*= std::map()*/) : _minimum(minimum), _maximum(maximum), _mean(mean), _stdDeviation(stdDeviation), _numVoxels(numVoxels), _volume(volume), _maximumVoxelPositions(maximumVoxelPositions), _minimumVoxelPositions(minimumVoxelPositions), _Dx(Dx), _Vx(Vx), _MOHx(MOHx), _MOCx(MOCx), _MaxOHx(MaxOHx), _MinOCx(MinOCx) { } DoseStatistics::~DoseStatistics() { } void DoseStatistics::setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions) { _minimumVoxelPositions = minimumVoxelPositions; } void DoseStatistics::setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions) { _maximumVoxelPositions = maximumVoxelPositions; } void DoseStatistics::setDx(const DoseToVolumeFunctionType& 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; } - unsigned int DoseStatistics::getNumberOfVoxels() const + VolumeType DoseStatistics::getNumberOfVoxels() const { return _numVoxels; } VolumeType DoseStatistics::getVolume() const { return _volume; } 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); } 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::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::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::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::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); } - double DoseStatistics::getValue(std::map aMap, double key, + double DoseStatistics::getValue(const std::map& aMap, double key, bool findNearestValueInstead, double& storedKey) const { - if (aMap.find(key) != aMap.end()) + 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(std::map aMap, + 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::getMaximumPositions() const { return _maximumVoxelPositions; } DoseStatistics::ResultListPointer DoseStatistics::getMinimumPositions() const { return _minimumVoxelPositions; } DoseStatistics::DoseToVolumeFunctionType DoseStatistics::getAllVx() const { return _Vx; } DoseStatistics::VolumeToDoseFunctionType DoseStatistics::getAllDx() 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 6003a1e..4fe884f 100644 --- a/code/algorithms/rttbDoseStatistics.h +++ b/code/algorithms/rttbDoseStatistics.h @@ -1,205 +1,265 @@ // ----------------------------------------------------------------------- // 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 #include "rttbBaseType.h" #include "rttbDoseIteratorInterface.h" namespace rttb { namespace algorithms { /*! @class DoseStatistics - @brief This is a data class storing different statistical values from a rt dose distribution - @sa DoseStatisticsCalculator + @brief Data class for storing different statistical values from a RT dose distribution + @details DoseStatisticsCalculator is used to compute the DoseStatistics */ class 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(std::map aMap, double key, bool findNearestValueInstead, + /*! @brief Returns the value of a map belonging to a key + @details if findNearestValueInstead=true and the key is not in the map, the value belonging to the nearest key is returned. This key is then storedKey. + */ + double getValue(const std::map& aMap, double key, bool findNearestValueInstead, double& storedKey) const; - std::map::const_iterator findNearestKeyInMap(std::map aMap, double key) 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; - unsigned int _numVoxels; + VolumeType _numVoxels; VolumeType _volume; VolumeToDoseFunctionType _Dx; DoseToVolumeFunctionType _Vx; VolumeToDoseFunctionType _MOHx; VolumeToDoseFunctionType _MOCx; VolumeToDoseFunctionType _MaxOHx; VolumeToDoseFunctionType _MinOCx; public: - /*! @brief Standard Constructor - */ - //DoseStatistics(); - /*! @brief Constructor - @detail the dose statistic values are set. Complex values maximumVoxelLocation, maximumVoxelLocation, Dx, Vx, MOHx, MOCx, MaxOHx and MinOCx are optional + @details the dose statistic values are set. Minimum, maximum, mean, stdDeviation, numVoxels and volume are mandatory. minimumVoxelPositions, maximumVoxelPositions, Dx, Vx, MOHx, MOCx, MaxOHx and MinOCx are optional. + @param minimum the minimum dose in Gy + @param maximum the maximum dose in Gy + @param mean the mean dose in Gy + @param stdDeviation the stdDeviation of the dose in Gy + @param numVoxels the amount of voxels (sub-voxel accurate, i.e. not necessarily whole voxels) + @param volume the volume in cm^3 + @param minimumVoxelPositions a vector containing the function minDose --> VoxelGrid as a pair + @param maximumVoxelPositions a vector containing the function maxDose --> VoxelGrid as a pair + @param Dx a map containing the function Volume x --> Dose y with: the minimal dose y is delivered to the volume x. + @param Vx a map containing the function Dose x --> Volume y with: the volume y is irradiated with a dose >= x. + @param MOHx a map containing the function Volume x --> Dose y with: the mean dose y of the hottest voxels in volume x. + @param MOCx a map containing the function Volume x --> Dose y with: the mean dose y of the coldest voxels in volume x. + @param MaxOHx a map containing the function Volume x --> Dose y with: the maximum dose y outside of the hottest voxels in volume x. + @param MinOCx a map containing the function Volume x --> Dose y with: the minimum dose y outside of the coldest voxels in volume x. */ DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean, - DoseStatisticType stdDeviation, unsigned int numVoxels, VolumeType volume, + DoseStatisticType stdDeviation, VolumeType numVoxels, VolumeType volume, ResultListPointer minimumVoxelPositions = boost::make_shared > > (std::vector >()), ResultListPointer maximumVoxelPositions = boost::make_shared > > (std::vector >()), - DoseToVolumeFunctionType Dx = DoseToVolumeFunctionType(), - VolumeToDoseFunctionType Vx = VolumeToDoseFunctionType(), + VolumeToDoseFunctionType Dx = VolumeToDoseFunctionType(), + DoseToVolumeFunctionType Vx = DoseToVolumeFunctionType(), VolumeToDoseFunctionType MOHx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MOCx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MaxOHx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MinOCx = VolumeToDoseFunctionType()); ~DoseStatistics(); void setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions); void setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions); - void setDx(const DoseToVolumeFunctionType& DxValues); - void setVx(const VolumeToDoseFunctionType& VxValues); + void setDx(const VolumeToDoseFunctionType& DxValues); + void setVx(const DoseToVolumeFunctionType& VxValues); void setMOHx(const VolumeToDoseFunctionType& MOHxValues); void setMOCx(const VolumeToDoseFunctionType& MOCxValues); void setMaxOHx(const VolumeToDoseFunctionType& MaxOHxValues); void setMinOCx(const VolumeToDoseFunctionType& MinOCxValues); /*! @brief Get number of voxels in doseIterator, with sub-voxel accuracy. */ - unsigned int getNumberOfVoxels() const; + VolumeType getNumberOfVoxels() const; VolumeType getVolume() const; /*! @brief Get the maximum of the current dose distribution. - @return Return the maximum dose in Gy + @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) + /*! @brief Get a vector of pairs of the the maximum dose VoxelGridIDs together with their dose value in Gy */ ResultListPointer getMaximumPositions() const; /*! @brief Get the minimum of the current dose distribution. - @return Return the minimum dose in Gy + @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) + /*! @brief Get a vector of pairs of the the minimum dose VoxelGridIDs together with their dose value in Gy */ ResultListPointer getMinimumPositions() const; /*! @brief Get the mean of the current dose distribution. - @return Return the mean dose in Gy + @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 + @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 + @return the variance in Gy */ DoseStatisticType getVariance() const; - + /*! @overload + @details calls getVx(xDoseAbsolute, false, Dummy) + */ + VolumeType getVx(DoseTypeGy xDoseAbsolute) const; /*! @brief Get Vx: the volume irradiated with a dose >= x. - @return Return absolute volume in absolute cm^3. + @param xDoseAbsolute the dose + @param findNearestValue if false: only the value (=volume) is returned that corresponds to the exact key (=dose). + If true: the value corresponding to the nearest key given is returned if the exact value is not in the map + @param nearestXDose the nearest key (=dose) found + @return Absolute volume in 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 + @exception NoDataException if the requested Dose is not in the precomputed vector */ - VolumeType getVx(DoseTypeGy xDoseAbsolute) const; VolumeType getVx(DoseTypeGy xDoseAbsolute, bool findNearestValue, DoseTypeGy& nearestXDose) const; + /*! @brief Get all computed Dx values + */ 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) + /*! + @brief Get Dx: the minimal dose delivered to part x of the current volume. + @param xVolumeAbsolute the absolute volume in cm^3 + @param findNearestValue if false: only the value (=dose) is returned that corresponds to the exact key (=volume). + If true: the value corresponding to the nearest key given is returned if the exact value is not in the map + @param nearestXVolume the nearest key (=volume) found + @return dose value in Gy. + @exception NoDataException if the Dx values have not been set (i.e. the vector is empty) + @exception NoDataException if the requested volume is not in the precomputed vector */ DoseTypeGy getDx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; + /*! @overload + @details calls getDx(xVolumeAbsolute, false, Dummy) + */ DoseTypeGy getDx(VolumeType xVolumeAbsolute) const; + /*! @brief Get all computed Vx values + */ VolumeToDoseFunctionType getAllDx() 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) + @exception NoDataException if the MOHx values have not been set (i.e. the vector is empty) + @exception NoDataException if the requested volume is not in the precomputed vector + @sa getMOHx() */ DoseTypeGy getMOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; + /*! @overload + @details calls getMOHx(xVolumeAbsolute, false, Dummy) + */ DoseTypeGy getMOHx(VolumeType xVolumeAbsolute) const; + /*! @brief Get all computed MOHx values + @return Return dose value in Gy. + */ 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) + @exception NoDataException if the MOCx values have not been set (i.e. the vector is empty) + @exception NoDataException if the requested volume is not in the precomputed vector + @sa getDx() */ DoseTypeGy getMOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; + /*! @overload + @details calls getMOCx(xVolumeAbsolute, false, Dummy) + */ DoseTypeGy getMOCx(VolumeType xVolumeAbsolute) const; + /*! @brief Get all computed MOCx values + @return Return dose value in Gy. + */ 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) + @exception NoDataException if the MaxOHx values have not been set (i.e. the vector is empty) + @exception NoDataException if the requested volume is not in the precomputed vector + @sa getDx() */ DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; + /*! @overload + @details calls getMaxOHx(xVolumeAbsolute, false, Dummy) + */ DoseTypeGy getMaxOHx(VolumeType xVolumeAbsolute) const; + /*! @brief Get all computed MaxOHx values + */ 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) + @exception NoDataException if the MinOCx values have not been set (i.e. the vector is empty) + @exception NoDataException if the requested volume is not in the precomputed vector + @sa getDx() */ DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; + /*! @overload + @details calls getMinOCx(xVolumeAbsolute, false, Dummy) + */ DoseTypeGy getMinOCx(VolumeType xVolumeAbsolute) const; + /*! @brief Get all computed MinOCx values + */ VolumeToDoseFunctionType getAllMinOCx() const; }; } } #endif diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp index b2de8f2..84bed4e 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,564 +1,554 @@ // ----------------------------------------------------------------------- // 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: 877 $ (last changed revision) // @date $Date: 2015-01-09 10:51:10 +0100 (Fr, 09 Jan 2015) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #include "rttbDoseStatisticsCalculator.h" #include #include #include #include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace algorithms { - DoseStatisticsCalculator::DoseStatisticsCalculator() - { - simpleDoseStatisticsCalculated = false; - _doseIterator = NULL; - } DoseStatisticsCalculator::DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator) - { - setDoseIterator(aDoseIterator); - simpleDoseStatisticsCalculated = false; - } - - - DoseStatisticsCalculator::~DoseStatisticsCalculator() - { - - } - - - void DoseStatisticsCalculator::setDoseIterator(DoseIteratorPointer aDoseIterator) { if (aDoseIterator == NULL) { throw core::NullPointerException("DoseIterator must not be NULL"); } else { _doseIterator = aDoseIterator; } + + _simpleDoseStatisticsCalculated = false; + } + + + DoseStatisticsCalculator::~DoseStatisticsCalculator() + { + } DoseStatisticsCalculator::DoseIteratorPointer DoseStatisticsCalculator::getDoseIterator() const { return _doseIterator; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( bool computeComplexMeasures, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } - //"simple" dose statistics are obligatory + //"simple" dose statistics are mandatory calculateSimpleDoseStatistics(); if (computeComplexMeasures) { //more complex dose statistics are optional calculateComplexDoseStatistics(precomputeDoseValues, precomputeVolumeValues); } return _statistics; } void DoseStatisticsCalculator::calculateSimpleDoseStatistics() { _doseVector.clear(); _voxelProportionVector.clear(); std::multimap doseValueVSIndexMap; std::vector voxelProportionVectorTemp; DoseStatisticType maximumDose = 0; - DoseStatisticType minimumDose = 0; - DoseStatisticType meanDose = 0; - DoseStatisticType stdDeviationDose = 0; + DoseStatisticType minimumDose = std::numeric_limits::max(); + DoseStatisticType meanDose; + DoseStatisticType stdDeviationDose; - float sum = 0; - rttb::DoseStatisticType numVoxels = 0.0; - float squareSum = 0; + 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 - std::multimap::iterator it; - - for (it = doseValueVSIndexMap.begin(); it != doseValueVSIndexMap.end(); ++it) + 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; + _simpleDoseStatisticsCalculated = true; ResultListPointer minimumVoxelPositions = computeMinimumPositions(100); ResultListPointer maximumVoxelPositions = computeMaximumPositions(100); _statistics->setMinimumVoxelPositions(minimumVoxelPositions); _statistics->setMaximumVoxelPositions(maximumVoxelPositions); - - } void DoseStatisticsCalculator::calculateComplexDoseStatistics(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()) { DoseTypeGy maxDose = _statistics->getMaximum(); std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02 * maxDose)(0.05 * maxDose)(0.1 * maxDose)( 0.9 * maxDose)( 0.95 * maxDose)(0.98 * maxDose); precomputeDoseValuesNonConst = defaultPrecomputeDoseValues; } if (precomputeVolumeValues.empty()) { VolumeType volume = _statistics->getVolume(); std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02 * volume)( 0.05 * volume)(0.1 * volume)(0.9 * volume) (0.95 * volume)(0.98 * volume); precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues; } DoseToVolumeFunctionType Vx = computeDoseToVolumeMulti(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->setMOHx(MOHx); _statistics->setMOCx(MOCx); _statistics->setMaxOHx(MaxOHx); _statistics->setMinOCx(MinOCx); } DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMaximumPositions( - unsigned int maxNumberMaxima) + unsigned int maxNumberMaxima) const { - if (!simpleDoseStatisticsCalculated) + if (!_simpleDoseStatisticsCalculated) { - throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMaximumAndPosition()"); - + 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) + unsigned int maxNumberMinima) const { - if (!simpleDoseStatisticsCalculated) + if (!_simpleDoseStatisticsCalculated) { - throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMinimumAndPosition()"); + 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; - int i = _doseVector.size() - 1; + int i = static_cast(_doseVector.size() - 1); for (; i >= 0; i--) { countVoxels += _voxelProportionVector.at(i); if (countVoxels >= noOfVoxel) { break; } } if (i >= 0) { resultDose = _doseVector.at(i); } else { 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 (int i = _doseVector.size() - 1; i >= 0; i--) + for (int i = static_cast(_doseVector.size() - 1); i >= 0; 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; - int i = _doseVector.size() - 1; + int i = static_cast(_doseVector.size() - 1); for (; i >= 0; i--) { countVoxels += _voxelProportionVector.at(i); if (countVoxels >= noOfVoxel) { break; } } if (i - 1 >= 0) { resultDose = _doseVector.at(i - 1); } 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::computeDoseToVolumeMulti( const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const { DoseToVolumeFunctionType VxMulti; for (int i = 0; i < precomputeDoseValues.size(); ++i) { if (name == DoseStatistics::Vx) { VxMulti.insert(std::pair(precomputeDoseValues.at(i), computeVx(precomputeDoseValues.at(i)))); } } return VxMulti; } DoseStatisticsCalculator::VolumeToDoseFunctionType DoseStatisticsCalculator::computeVolumeToDoseFunctionMulti( const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const { VolumeToDoseFunctionType multiValues; for (int i = 0; i < precomputeVolumeValues.size(); ++i) { switch (name) { case DoseStatistics::Dx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeDx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MOHx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeMOHx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MOCx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeMOCx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MaxOHx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeMaxOHx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MinOCx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeMinOCx(precomputeVolumeValues.at(i)))); break; default: throw core::InvalidParameterException("unknown DoseStatistics name!"); } } return multiValues; } }//end namespace algorithms }//end namespace rttb diff --git a/code/algorithms/rttbDoseStatisticsCalculator.h b/code/algorithms/rttbDoseStatisticsCalculator.h index 8d01d25..91d64da 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.h +++ b/code/algorithms/rttbDoseStatisticsCalculator.h @@ -1,120 +1,141 @@ // ----------------------------------------------------------------------- // 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: 707 $ (last changed revision) // @date $Date: 2014-09-04 16:37:24 +0200 (Do, 04 Sep 2014) $ (last change date) // @author $Author: floca $ (last changed by) */ #ifndef __DOSE_STATISTICS_CALCULATOR_H #define __DOSE_STATISTICS_CALCULATOR_H #include #include "rttbBaseType.h" #include "rttbDoseIteratorInterface.h" #include "rttbDoseStatistics.h" namespace rttb { namespace algorithms { /*! @class DoseStatisticsCalculator - @brief This is a class calculating different statistical values from a rt dose distribution - These values range from standard statistical values such as minimum, maximum and mean to more - dose specific properties such as Vx (volume irradiated with a dose >=x), Dx (minimal dose delivered - to x% of the VOI, and MOHx (mean in the hottest volume). + @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 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; - /* Contains relevant dose values sorted in descending order. + /*! @brief Contains relevant dose values sorted in descending order. */ std::vector _doseVector; - /*! Contains the corresponding voxel proportions to the values in doseVector. + /*! @brief Contains the corresponding voxel proportions to the values in doseVector. */ std::vector _voxelProportionVector; - - - DoseStatisticsPointer _statistics; - - bool simpleDoseStatisticsCalculated; - - /*! @brief Calculation of basic dose statistics. It will be called in Constructor. - @exception NullPointerException Thrown if _doseIterator is NULL. + /*! @brief The doseStatistics are stored here. */ + DoseStatisticsPointer _statistics; + bool _simpleDoseStatisticsCalculated; - ResultListPointer computeMaximumPositions(unsigned int maxNumberMaxima); + ResultListPointer computeMaximumPositions(unsigned int maxNumberMaxima) const; + ResultListPointer computeMinimumPositions(unsigned int maxNumberMinima) const; - ResultListPointer computeMinimumPositions(unsigned int maxNumberMinima); VolumeType computeVx(DoseTypeGy xDoseAbsolute) const; DoseTypeGy computeDx(VolumeType xVolumeAbsolute) const; DoseTypeGy computeMOHx(VolumeType xVolumeAbsolute) const; DoseTypeGy computeMOCx(DoseTypeGy xVolumeAbsolute) const; DoseTypeGy computeMaxOHx(DoseTypeGy xVolumeAbsolute) const; DoseTypeGy computeMinOCx(DoseTypeGy xVolumeAbsolute) const; DoseToVolumeFunctionType computeDoseToVolumeMulti(const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const; VolumeToDoseFunctionType computeVolumeToDoseFunctionMulti(const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const; + /*! @brief Calculatues simple dose statistics (min, mean, max, stdDev, minDosePositions, maxDosePositions) + */ void calculateSimpleDoseStatistics(); + /*! @brief Calculatues 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(const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues); public: - /*! @brief Standard Constructor - */ - DoseStatisticsCalculator(); - ~DoseStatisticsCalculator(); /*! @brief Constructor + @param aDoseIterator the dose to be analyzed */ DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator); - /*! @brief Set new dose data for statistics. Statistics will be re-initialized. - */ - void setDoseIterator(DoseIteratorPointer aDoseIterator); - DoseIteratorPointer getDoseIterator() const; + /*! @brief Calculatues dose statistics + @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 part x of the current volume) +
  • Vx (the volume irradiated with a dose >= x) +
  • MOHx (mean dose of the hottest x voxels) +
  • MOCx (mean dose of the coldest x voxels) +
  • MaxOHx (Maximum outside of the hottest x voxels) +
  • MinOCx (Minimum outside of the coldest x voxels) +
+ Default values for precomputeDoseValues are 0.02, 0.05, 0.1, 0.9, 0.95 and 0.98 with respect to maxDose. + Default values for precomputeVolumeValues 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? + @param precomputeDoseValues the dose values for Vx precomputation + @param precomputeVolumeValues the volume values for Dx, MOHx, MOCx, MaxOHx and MinOCx precomputation + @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! 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(bool computeComplexMeasures = false, const std::vector& precomputeDoseValues = std::vector(), const std::vector& precomputeVolumeValues = std::vector()); }; } } #endif diff --git a/code/io/other/rttbDoseStatisticsXMLWriter.h b/code/io/other/rttbDoseStatisticsXMLWriter.h index 85ea353..6713b1b 100644 --- a/code/io/other/rttbDoseStatisticsXMLWriter.h +++ b/code/io/other/rttbDoseStatisticsXMLWriter.h @@ -1,89 +1,92 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __DOSE_STATISTICS_XML_WRITER_H #define __DOSE_STATISTICS_XML_WRITER_H #include "rttbDoseStatistics.h" #include "rttbBaseType.h" /*boost includes*/ #include #include #include #include -namespace rttb{ - namespace io{ - namespace other{ +namespace rttb +{ + namespace io + { + namespace other + { typedef boost::shared_ptr DoseStatisticsPtr; typedef rttb::algorithms::DoseStatistics::ResultListPointer ResultListPointer; /*! @brief Write statistics to boost::property_tree::ptree. @param aReferenceDose A reference dose for the calculation of Vx @param aDoseStatistics DoseStatistics to write @pre aReferenceDose must >0 @exception InvalidParameterException Thrown if aReferenceDose<=0 or xml_parse_error */ - boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose=100); + boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose = 100); /*! @brief Write statistics to String. @param aReferenceDose A reference dose for the calculation of Vx @param aDoseStatistics DoseStatistics to write @pre aReferenceDose must >0 @exception InvalidParameterException Thrown if aReferenceDose<=0 or xml_parse_error */ - XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose=100); + XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose = 100); - /*! @brief Write statistics to xml file, including - numberOfVoxels, - volume, - minimum, - maximum, - mean, - standard deviation, - variance, + /*! @brief Write statistics to xml file, including + numberOfVoxels, + volume, + minimum, + maximum, + mean, + standard deviation, + variance, D2,D5,D10,D90,D95,D98, V2,V5,V10,V90,V95,V98, MOH2,MOH5,MOH10, MOC2,MOC5,MOC10, MaxOH2,MaxOH5,MaxOH10, MinOC2,MinOC5,MinOC10. @param aReferenceDose A reference dose for the calculation of Vx @param aDoseStatistics DoseStatistics to write @param aFileName To write dose statistics @pre aReferenceDose must >0 @exception InvalidParameterException Thrown if aReferenceDose<=0 or xml_parse_error */ - void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName, DoseTypeGy aReferenceDose=100); + void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName, DoseTypeGy aReferenceDose = 100); /*! @brief Write statistics to String to generate a table: "Volume mm3@Max@Min@Mean@Std.Dev.@Variance@V2@V5@V10@V90@V95@V98@D2@D5@D10@D90@D95@D98@MOH2@MOH5@MOH10@MOC2@MOC5@MOC10@MaxOH2@MaxOH5@MaxOH10@MinOC2@MinOC5@MinOC10" @param aReferenceDose A reference dose for the calculation of Vx @param aDoseStatistics DoseStatistics to write @pre aReferenceDose must >0 @exception InvalidParameterException Thrown if aReferenceDose<=0 or xml_parse_error */ - StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose=100); + StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose = 100); } } } #endif diff --git a/testing/algorithms/DoseStatisticsCalculatorTest.cpp b/testing/algorithms/DoseStatisticsCalculatorTest.cpp index 1a083d8..958d443 100644 --- a/testing/algorithms/DoseStatisticsCalculatorTest.cpp +++ b/testing/algorithms/DoseStatisticsCalculatorTest.cpp @@ -1,313 +1,296 @@ // ----------------------------------------------------------------------- // 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: 707 $ (last changed revision) // @date $Date: 2014-09-04 16:37:24 +0200 (Do, 04 Sep 2014) $ (last change date) // @author $Author: floca $ (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 "../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; 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; //1) test constructors // the values cannot be accessed from outside, therefore correct default values are not tested - CHECK_NO_THROW(rttb::algorithms::DoseStatisticsCalculator myEmptyDoseStatCalculatur); - rttb::algorithms::DoseStatisticsCalculator myEmptyDoseStatCalculatur; - rttb::algorithms::DoseStatisticsCalculator myDoseStatCalculaturToBeFilled; 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 - CHECK_THROW_EXPLICIT(myDoseStatCalculaturToBeFilled.setDoseIterator(spDoseIteratorNull), core::NullPointerException); - CHECK_NO_THROW(myDoseStatCalculaturToBeFilled.setDoseIterator(spDoseIterator)); //3) test calculateDoseStatistics - DoseStatisticsPointer theStatistics, theStatistics2; - CHECK_THROW_EXPLICIT(myEmptyDoseStatCalculatur.calculateDoseStatistics(), core::NullPointerException); + DoseStatisticsPointer theStatistics; CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics()); - CHECK_NO_THROW(theStatistics2 = myDoseStatCalculaturToBeFilled.calculateDoseStatistics()); CHECK_EQUAL(theStatistics->getMinimumPositions()->empty(), false); CHECK_EQUAL(theStatistics->getMaximumPositions()->empty(), false); CHECK_EQUAL(theStatistics->getAllVx().empty(), true); CHECK_EQUAL(theStatistics->getAllDx().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; 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->getNumberOfVoxels())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.05 * theStatisticsDefault->getNumberOfVoxels())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.1 * theStatisticsDefault->getNumberOfVoxels())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.9 * theStatisticsDefault->getNumberOfVoxels())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.95 * theStatisticsDefault->getNumberOfVoxels())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.98 * theStatisticsDefault->getNumberOfVoxels())); + 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 manually set values for omputeComplexMeasures=true std::vector precomputeDoseValues, precomputeVolumeValues; precomputeDoseValues.push_back(0.01 * theStatistics->getMaximum()); precomputeDoseValues.push_back(0.02 * theStatistics->getMaximum()); precomputeDoseValues.push_back(0.05 * theStatistics->getMaximum()); - precomputeVolumeValues.push_back(0.9 * theStatistics->getNumberOfVoxels()); - precomputeVolumeValues.push_back(0.95 * theStatistics->getNumberOfVoxels()); - precomputeVolumeValues.push_back(0.99 * theStatistics->getNumberOfVoxels()); + precomputeVolumeValues.push_back(0.9 * theStatistics->getVolume()); + precomputeVolumeValues.push_back(0.95 * theStatistics->getVolume()); + precomputeVolumeValues.push_back(0.99 * theStatistics->getVolume()); CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(true, 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->getNumberOfVoxels())); - CHECK_NO_THROW(theStatistics->getDx(0.95 * theStatistics->getNumberOfVoxels())); - CHECK_NO_THROW(theStatistics->getDx(0.99 * theStatistics->getNumberOfVoxels())); - CHECK_THROW_EXPLICIT(theStatistics->getDx(0.03 * theStatistics->getNumberOfVoxels()), 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()), 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->getNumberOfVoxels()), - theStatisticsDefault->getDx(0.9 * theStatistics->getNumberOfVoxels())); - CHECK_EQUAL(theStatistics->getDx(0.95 * theStatistics->getNumberOfVoxels()), - theStatisticsDefault->getDx(0.95 * theStatistics->getNumberOfVoxels())); + 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())); //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()); - CHECK_EQUAL(theStatistics2->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(); - float compMean = (int(mean * 100)) / 100; + 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(theStatistics2->getMaximum(), maximum); - CHECK_EQUAL(theStatistics->getMinimum(), minimum); - CHECK_EQUAL(theStatistics2->getMinimum(), minimum); - CHECK_CLOSE(theStatistics->getMean(), mean, errorConstantLarger); - CHECK_CLOSE(theStatistics2->getMean(), mean, errorConstantLarger); - CHECK_CLOSE(theStatistics->getStdDeviation(), stdDev, errorConstantLarger); - CHECK_CLOSE(theStatistics2->getStdDeviation(), stdDev, errorConstantLarger); - CHECK_CLOSE(theStatistics->getVariance(), variance, errorConstantLarger); - CHECK_CLOSE(theStatistics2->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->getMaximumPositions(); auto minimaPositions = theStatistics->getMinimumPositions(); 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->getMaximumPositions(); minimaPositions = theStatistics3->getMinimumPositions(); 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()); } RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/examples/RTDoseStatisticsTest.cpp b/testing/examples/RTDoseStatisticsTest.cpp index 5282d4f..1a0745c 100644 --- a/testing/examples/RTDoseStatisticsTest.cpp +++ b/testing/examples/RTDoseStatisticsTest.cpp @@ -1,200 +1,203 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests -#include - #include #include -#include "dcmtk/config/osconfig.h" /* make sure OS specific configuration is included first */ -#include "dcmtk/dcmrt/drtdose.h" -#include "dcmtk/dcmdata/dcfilefo.h" - #include "litCheckMacros.h" + #include "rttbDoseStatistics.h" #include "rttbDoseStatisticsCalculator.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbVirtuosPlanFileDoseAccessorGenerator.h" #include "rttbVirtuosFileStructureSetGenerator.h" #include "rttbOTBMaskAccessor.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" -#include "rttbNullPointerException.h" #include "rttbBaseType.h" namespace rttb { namespace testing { + PREPARE_DEFAULT_TEST_REPORTING; /*! @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!*/ + + const double reducedErrorConstant = 0.0001; + const double expectedVal = 5.64477e-005; + + void testWithDummyDoseData(const std::string& doseFilename); + void testWithRealVirtuosDoseDataAndStructure(const std::string& doseFilename, const std::string& structFilename, + const std::string& planFilename, unsigned int structureNr); + int RTDoseStatisticsTest(int argc, char* argv[]) { - PREPARE_DEFAULT_TEST_REPORTING; typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericDoseIterator::DoseIteratorPointer DoseIteratorPointer; typedef algorithms::DoseStatistics::ResultListPointer ResultListPointer; - //ARGUMENTS: 1: dose1 file name - // 2: dose2 file name - // 3: RTStruct filename - std::string RTDOSE_FILENAME; std::string RTDOSE2_FILENAME; std::string RTSTRUCT_FILENAME; std::string RTPLAN_FILENAME; if (argc > 4) { RTDOSE_FILENAME = argv[1]; RTDOSE2_FILENAME = argv[2]; RTSTRUCT_FILENAME = argv[3]; RTPLAN_FILENAME = argv[4]; } else { std::cout << "at least four arguments required for RTDoseStatisticsTest" << std::endl; return -1; } + testWithDummyDoseData(RTDOSE_FILENAME); + //Structure 2 is RUECKENMARK + testWithRealVirtuosDoseDataAndStructure(RTDOSE2_FILENAME, RTSTRUCT_FILENAME, RTPLAN_FILENAME, 2); - OFCondition status; - DcmFileFormat fileformat; + RETURN_AND_REPORT_TEST_SUCCESS; + } - const double expectedVal = 5.64477e-005; - const double reducedErrorConstant = 0.0001; + void testWithDummyDoseData(const std::string& doseFilename) + { + typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; + typedef core::GenericDoseIterator::DoseIteratorPointer DoseIteratorPointer; - rttb::algorithms::DoseStatisticsCalculator doseStatisticsCalculator; typedef boost::shared_ptr > > ResultsVectorPointer; ::DRTDoseIOD rtdoseDKFZ; - io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); + io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(doseFilename.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create corresponding DoseIterator boost::shared_ptr spDoseIteratorTmp = boost::make_shared(doseAccessor1); DoseIteratorPointer spDoseIterator(spDoseIteratorTmp); - doseStatisticsCalculator.setDoseIterator(spDoseIterator); + rttb::algorithms::DoseStatisticsCalculator doseStatisticsCalculator(spDoseIterator); std::vector precomputedDoseValues; precomputedDoseValues.push_back(0); precomputedDoseValues.push_back(expectedVal); precomputedDoseValues.push_back(expectedVal + reducedErrorConstant); std::vector precomputedVolumeValues; precomputedVolumeValues.push_back(20000); precomputedVolumeValues.push_back(24120); rttb::algorithms::DoseStatistics::DoseStatisticsPointer doseStatistics = doseStatisticsCalculator.calculateDoseStatistics(true, precomputedDoseValues, precomputedVolumeValues); CHECK_CLOSE(doseStatistics->getMean(), expectedVal, errorConstant); CHECK_CLOSE(doseStatistics->getStdDeviation(), 0, errorConstant); CHECK_CLOSE(doseStatistics->getVariance(), 0, errorConstant); DoseTypeGy vx = doseStatistics->getVx(expectedVal); CHECK_EQUAL(vx, doseStatistics->getVx(0)); CHECK_CLOSE(expectedVal, doseStatistics->getDx(vx), reducedErrorConstant); CHECK_CLOSE(doseStatistics->getMaximum(), expectedVal, errorConstant); CHECK_CLOSE(doseStatistics->getMinimum(), expectedVal, errorConstant); auto minListPtr = doseStatistics->getMinimumPositions(); auto maxListPtr = doseStatistics->getMaximumPositions(); CHECK_EQUAL(maxListPtr->size(), 100); CHECK_EQUAL(minListPtr->size(), 100); CHECK_CLOSE(doseStatistics->getDx(24120), doseStatistics->getMinimum(), 0.001); CHECK_CLOSE(doseStatistics->getMOHx(24120), doseStatistics->getMean(), reducedErrorConstant); CHECK_CLOSE(doseStatistics->getMOCx(20000), doseStatistics->getMean(), reducedErrorConstant); CHECK_CLOSE(doseStatistics->getMinOCx(20000), doseStatistics->getMean(), reducedErrorConstant); + } + + void testWithRealVirtuosDoseDataAndStructure(const std::string& doseFilename, const std::string& structFilename, + const std::string& planFilename, unsigned int structureNr) + { + typedef core::GenericDoseIterator::DoseIteratorPointer DoseIteratorPointer; - //test with real data set (virtuos) typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef rttb::algorithms::DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer; - auto virtuosDoseAccessor = io::virtuos::VirtuosPlanFileDoseAccessorGenerator(RTDOSE2_FILENAME.c_str(), - RTPLAN_FILENAME.c_str()).generateDoseAccessor(); + auto virtuosDoseAccessor = io::virtuos::VirtuosPlanFileDoseAccessorGenerator(doseFilename.c_str(), + planFilename.c_str()).generateDoseAccessor(); auto virtuosStructureSet = io::virtuos::VirtuosFileStructureSetGenerator( - RTSTRUCT_FILENAME.c_str(), RTDOSE2_FILENAME.c_str()).generateStructureSet(); + structFilename.c_str(), doseFilename.c_str()).generateStructureSet(); - //Structure 2 is RUECKENMARK boost::shared_ptr spOTBMaskAccessorVirtuos = - boost::make_shared(virtuosStructureSet->getStructure(2), + boost::make_shared(virtuosStructureSet->getStructure(structureNr), virtuosDoseAccessor->getGeometricInfo()); spOTBMaskAccessorVirtuos->updateMask(); MaskAccessorPointer spMaskAccessor(spOTBMaskAccessorVirtuos); //create corresponding MaskedDoseIterator boost::shared_ptr spMaskedDoseIteratorTmp = boost::make_shared(spMaskAccessor, virtuosDoseAccessor); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); rttb::algorithms::DoseStatisticsCalculator doseStatisticsCalculatorVirtuos(spMaskedDoseIterator); DoseStatisticsPointer doseStatisticsVirtuos = doseStatisticsCalculatorVirtuos.calculateDoseStatistics(true); //comparison values computed with "old" DoseStatistics implementation CHECK_EQUAL(doseStatisticsVirtuos->getMinimum(), 0); CHECK_CLOSE(doseStatisticsVirtuos->getMaximum(), 35.3075, reducedErrorConstant); CHECK_CLOSE(doseStatisticsVirtuos->getMean(), 16.1803, reducedErrorConstant); CHECK_CLOSE(doseStatisticsVirtuos->getStdDeviation(), 9.84436, reducedErrorConstant); auto maxPositions = doseStatisticsVirtuos->getMaximumPositions(); auto minPositions = doseStatisticsVirtuos->getMinimumPositions(); CHECK_EQUAL(maxPositions->size(), 1); CHECK_EQUAL(minPositions->size(), 100); CHECK_EQUAL(maxPositions->begin()->first, doseStatisticsVirtuos->getMaximum()); CHECK_EQUAL(maxPositions->begin()->second, 2138227); CHECK_EQUAL(minPositions->begin()->first, doseStatisticsVirtuos->getMinimum()); CHECK_EQUAL(minPositions->begin()->second, 39026); CHECK_CLOSE(doseStatisticsVirtuos->getDx(0.02 * doseStatisticsVirtuos->getVolume()), 30.1528, reducedErrorConstant); CHECK_CLOSE(doseStatisticsVirtuos->getVx(0.9 * doseStatisticsVirtuos->getMaximum()), 0.483529, reducedErrorConstant); CHECK_CLOSE(doseStatisticsVirtuos->getMOHx(0.1 * doseStatisticsVirtuos->getVolume()), 28.7019, reducedErrorConstant); CHECK_CLOSE(doseStatisticsVirtuos->getMOCx(0.05 * doseStatisticsVirtuos->getVolume()), 0, reducedErrorConstant); CHECK_CLOSE(doseStatisticsVirtuos->getMaxOHx(0.95 * doseStatisticsVirtuos->getVolume()), 0, reducedErrorConstant); CHECK_CLOSE(doseStatisticsVirtuos->getMinOCx(0.98 * doseStatisticsVirtuos->getVolume()), 30.1756, reducedErrorConstant); - std::cout << "-->" << 0.1 * doseStatisticsVirtuos->getVolume() << std::endl; - - RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb \ No newline at end of file