diff --git a/code/algorithms/rttbDoseStatistics.cpp b/code/algorithms/rttbDoseStatistics.cpp index 86f21a1..6669f52 100644 --- a/code/algorithms/rttbDoseStatistics.cpp +++ b/code/algorithms/rttbDoseStatistics.cpp @@ -1,343 +1,352 @@ // ----------------------------------------------------------------------- // 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 "rttbDoseStatistics.h" #include "rttbDataNotAvailableException.h" namespace rttb { namespace algorithms { DoseStatistics::DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean, DoseStatisticType stdDeviation, unsigned int numVoxels, VolumeType volume, ResultListPointer maximumVoxelPositions /*= ResultListPointer()*/, ResultListPointer minimumVoxelPositions /*= ResultListPointer()*/, VolumeToDoseFunctionType 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), _maximumVoxelPositions(maximumVoxelPositions), _minimumVoxelPositions(minimumVoxelPositions), _Dx(Dx), _Vx(Vx), _MOHx(MOHx), _MOCx(MOCx), _MaxOHx(MaxOHx), _MinOCx(MinOCx) { if (maximumVoxelPositions == nullptr) { _maximumVoxelPositions = boost::make_shared > > (std::vector >()); } if (minimumVoxelPositions == nullptr) { _minimumVoxelPositions = boost::make_shared > > (std::vector >()); } 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) { _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; + } + } + unsigned int 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); } 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(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::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 4fe7e26..8049864 100644 --- a/code/algorithms/rttbDoseStatistics.h +++ b/code/algorithms/rttbDoseStatistics.h @@ -1,213 +1,214 @@ // ----------------------------------------------------------------------- // 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 "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 */ 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(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; unsigned int _numVoxels; VolumeType _volume; DoseTypeGy _referenceDose; //for Vx computation 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 */ DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean, DoseStatisticType stdDeviation, unsigned int numVoxels, VolumeType volume, ResultListPointer minimumVoxelPositions = nullptr, ResultListPointer maximumVoxelPositions = nullptr, VolumeToDoseFunctionType Dx = VolumeToDoseFunctionType(), 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 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. */ unsigned int 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 getMaximumPositions() 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 getMinimumPositions() 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; 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; 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) */ DoseTypeGy getMOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const; DoseTypeGy getMOHx(VolumeType xVolumeAbsolute) 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; 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; 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; VolumeToDoseFunctionType getAllMinOCx() const; }; } } #endif diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp index d70dc0b..cfbcb33 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,599 +1,600 @@ // ----------------------------------------------------------------------- // 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 "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace algorithms { DoseStatisticsCalculator::DoseStatisticsCalculator(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, 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; } 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->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; 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 = 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 = 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::computeDoseToVolumeFunctionMulti(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const { DoseToVolumeFunctionType VxMulti; for (int i = 0; i < precomputeDoseValues.size(); ++i) { if (name == DoseStatistics::Vx) { double xAbsolue = precomputeDoseValues.at(i) * referenceDose; VxMulti.insert(std::pair(xAbsolue, computeVx(xAbsolue))); } else { throw core::InvalidParameterException("unknown DoseStatistics name!"); } } return VxMulti; } DoseStatisticsCalculator::VolumeToDoseFunctionType DoseStatisticsCalculator::computeVolumeToDoseFunctionMulti( const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const { VolumeToDoseFunctionType multiValues; VolumeType volume = _statistics->getVolume(); for (int i = 0; i < precomputeVolumeValues.size(); ++i) { double xAbsolute = precomputeVolumeValues.at(i) * volume; 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!"); } } return multiValues; } }//end namespace algorithms }//end namespace rttb diff --git a/testing/algorithms/CMakeLists.txt b/testing/algorithms/CMakeLists.txt index 07cb2cc..03d03d7 100644 --- a/testing/algorithms/CMakeLists.txt +++ b/testing/algorithms/CMakeLists.txt @@ -1,24 +1,25 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(ALGORITHMS_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbAlgorithmsTests) SET(ALGORITHMS_HEADER_TEST ${EXECUTABLE_OUTPUT_PATH}/rttbAlgorithmsHeaderTest) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(DoseStatisticsTest ${ALGORITHMS_TESTS} DoseStatisticsTest) ADD_TEST(ArithmeticTest ${ALGORITHMS_TESTS} ArithmeticTest) +ADD_TEST(DoseStatisticsCalculatorTest ${ALGORITHMS_TESTS} DoseStatisticsCalculatorTest) ADD_TEST(BinaryFunctorAccessorTest ${ALGORITHMS_TESTS} BinaryFunctorAccessorTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/dicompylerTestDose.dcm") RTTB_CREATE_TEST_MODULE(rttbAlgorithms DEPENDS RTTBAlgorithms RTTBMasks RTTBDicomIO PACKAGE_DEPENDS Boost Litmus DCMTK) IF (NOT WIN32) #CMake 3.1 provides target_compile_features(RTTB_Interpolation cxx_auto_type cxx_nullptr cxx_override) to automatically add required compiler flags set(CMAKE_CXX_FLAGS "-std=c++11") ENDIF() diff --git a/testing/algorithms/DoseStatisticsCalculatorTest.cpp b/testing/algorithms/DoseStatisticsCalculatorTest.cpp index 296287b..379bce9 100644 --- a/testing/algorithms/DoseStatisticsCalculatorTest.cpp +++ b/testing/algorithms/DoseStatisticsCalculatorTest.cpp @@ -1,314 +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 #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_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->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->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 reference dose and the default x values CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(100.0)); - CHECK_NO_THROW(theStatistics->getVx(0.1 * theStatisticsDefault->getMaximum())); - CHECK_NO_THROW(theStatistics->getVx(0.9 * theStatisticsDefault->getMaximum())); - CHECK_NO_THROW(theStatistics->getDx(0.1 * theStatisticsDefault->getVolume())); - CHECK_NO_THROW(theStatistics->getDx(0.9 * theStatisticsDefault->getVolume())); - CHECK_NO_THROW(theStatistics->getMOHx(0.95 * theStatisticsDefault->getVolume())); - CHECK_NO_THROW(theStatistics->getMOCx(0.98 * theStatisticsDefault->getVolume())); + 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->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 * theStatistics->getMaximum()); - precomputeDoseValues.push_back(0.02 * theStatistics->getMaximum()); - precomputeDoseValues.push_back(0.05 * theStatistics->getMaximum()); + precomputeDoseValues.push_back(0.01); + precomputeDoseValues.push_back(0.02); + precomputeDoseValues.push_back(0.05); - precomputeVolumeValues.push_back(0.9 * theStatistics->getVolume()); - precomputeVolumeValues.push_back(0.95 * theStatistics->getVolume()); - precomputeVolumeValues.push_back(0.99 * theStatistics->getVolume()); + 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()), core::DataNotAvailableException); CHECK_EQUAL(theStatistics->getVx(0.02 * theStatistics->getMaximum()), - theStatisticsDefault->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())); + theStatisticsDefault->getVx(0.05 * theStatistics->getMaximum())); CHECK_EQUAL(theStatistics->getDx(0.9 * theStatistics->getVolume()), - theStatisticsDefault->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())); + theStatisticsDefault->getDx(0.95 * theStatistics->getVolume())); //check manually set reference dose and x values CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(precomputeDoseValues, precomputeVolumeValues, 100.0)); - CHECK_NO_THROW(theStatistics->getVx(0.01 * theStatistics->getMaximum())); + 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_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->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