diff --git a/code/algorithms/rttbDoseStatistics.cpp b/code/algorithms/rttbDoseStatistics.cpp index caaa119..6669f52 100644 --- a/code/algorithms/rttbDoseStatistics.cpp +++ b/code/algorithms/rttbDoseStatistics.cpp @@ -1,330 +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()*/) : _minimum(minimum), - _maximum(maximum), _mean(mean), _stdDeviation(stdDeviation), _numVoxels(numVoxels), _volume(volume), + 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 2f6e9e2..8049864 100644 --- a/code/algorithms/rttbDoseStatistics.h +++ b/code/algorithms/rttbDoseStatistics.h @@ -1,204 +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()); + 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 e633a1f..cfbcb33 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,560 +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, const std::vector& precomputeDoseValues, - const std::vector& precomputeVolumeValues, unsigned int maxNumberMinimaPositions, - unsigned int maxNumberMaximaPositions) - { - + 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 - calculateComplexDoseStatistics(precomputeDoseValues, precomputeVolumeValues); + 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(const std::vector& precomputeDoseValues, + 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()) { - 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); + std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)(0.95)(0.98); 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); + std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)(0.95)(0.98); precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues; } - DoseToVolumeFunctionType Vx = computeDoseToVolumeFunctionMulti(precomputeDoseValuesNonConst, DoseStatistics::Vx); + 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( + 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) { - VxMulti.insert(std::pair(precomputeDoseValues.at(i), - computeVx(precomputeDoseValues.at(i)))); + 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(precomputeVolumeValues.at(i), - computeDx(precomputeVolumeValues.at(i)))); + multiValues.insert(std::pair(xAbsolute, + computeDx(xAbsolute))); break; case DoseStatistics::MOHx: - multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeMOHx(precomputeVolumeValues.at(i)))); + multiValues.insert(std::pair(xAbsolute, + computeMOHx(xAbsolute))); break; case DoseStatistics::MOCx: - multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeMOCx(precomputeVolumeValues.at(i)))); + multiValues.insert(std::pair(xAbsolute, + computeMOCx(xAbsolute))); break; case DoseStatistics::MaxOHx: - multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeMaxOHx(precomputeVolumeValues.at(i)))); + multiValues.insert(std::pair(xAbsolute, + computeMaxOHx(xAbsolute))); break; case DoseStatistics::MinOCx: - multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeMinOCx(precomputeVolumeValues.at(i)))); + 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/code/algorithms/rttbDoseStatisticsCalculator.h b/code/algorithms/rttbDoseStatisticsCalculator.h index 2ff92b0..99f1af2 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.h +++ b/code/algorithms/rttbDoseStatisticsCalculator.h @@ -1,155 +1,175 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __DOSE_STATISTICS_CALCULATOR_H #define __DOSE_STATISTICS_CALCULATOR_H #include #include "rttbBaseType.h" #include "rttbDoseIteratorInterface.h" #include "rttbDoseStatistics.h" namespace rttb { namespace algorithms { /*! @class DoseStatisticsCalculator @brief Class for calculating different statistical values from a RT dose distribution @details These values range from standard statistical values such as minimum, maximum and mean to more complex dose specific measures such as Vx (volume irradiated with a dose >=x), Dx (minimal dose delivered to x% of the VOI) or MOHx (mean in the hottest volume). For a complete list, see calculateDoseStatistics(). @note the complex dose statistics are precomputed and cannot be computed "on the fly" lateron! The doses/volumes that should be used for precomputation have to be set in calculateDoseStatistics() */ class DoseStatisticsCalculator { public: typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; typedef DoseStatistics::ResultListPointer ResultListPointer; typedef DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer; typedef DoseStatistics::DoseToVolumeFunctionType DoseToVolumeFunctionType; typedef DoseStatistics::VolumeToDoseFunctionType VolumeToDoseFunctionType; private: DoseIteratorPointer _doseIterator; /*! @brief Contains relevant dose values sorted in descending order. */ std::vector _doseVector; /*! @brief Contains the corresponding voxel proportions to the values in doseVector. */ std::vector _voxelProportionVector; /*! @brief The doseStatistics are stored here. */ DoseStatisticsPointer _statistics; bool _simpleDoseStatisticsCalculated; /*! @brief Calculates the positions where the dose has its maximum @param maxNumberMaximaPositions the maximal amount of computed positions @pre maximumDose must be defined in _statistics with the correct value */ ResultListPointer computeMaximumPositions(unsigned int maxNumberMaximaPositions) const; /*! @brief Calculates the positions where the dose has its minimum @param maxNumberMinimaPositions the maximal amount of computed positions (they are read sequentially using the iterator until maxNumberMinimaPositions have been read, other positions are not considered) @pre minimumDose must be defined in _statistics with the correct value */ ResultListPointer computeMinimumPositions(unsigned int maxNumberMinimaPositions) const; VolumeType computeVx(DoseTypeGy xDoseAbsolute) const; 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 computeDoseToVolumeFunctionMulti(const std::vector& precomputeDoseValues, + DoseToVolumeFunctionType computeDoseToVolumeFunctionMulti(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const; VolumeToDoseFunctionType computeVolumeToDoseFunctionMulti(const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const; /*! @brief Calculates simple dose statistics (min, mean, max, stdDev, minDosePositions, maxDosePositions) @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed */ void calculateSimpleDoseStatistics(unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions); /*! @brief Calculates complex dose statistics (Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx) @warning computations can take quite long (>1 min) for large structures as many statistics are precomputed */ - void calculateComplexDoseStatistics(const std::vector& precomputeDoseValues, + void calculateComplexDoseStatistics(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues); public: ~DoseStatisticsCalculator(); /*! @brief Constructor @param aDoseIterator the dose to be analyzed */ DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator); DoseIteratorPointer getDoseIterator() const; - /*! @brief 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 - @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed - @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed - @warning if computeComplexMeasures==true, computations can take quite long (>1 min) for large structures as many statistics are precomputed - @note the complex dose statistics are precomputed and cannot be computed "on the fly" lateron! 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! + /*! @brief Compute simple or complex dose statistics with default relative x values and the maximum dose as default reference dose (for Vx computation) + @details The following statistics are calculated always (i.e. also if computeComplexMeasures=false): +
    +
  • minimum dose +
  • mean dose +
  • maximum dose +
  • standard deviation dose +
  • voxel positions of minimum dose +
  • voxel positions of maximum dose +
+ Additionally, these statistics are computed if computeComplexMeasures=true: +
    +
  • Dx (the minimal dose delivered to a volume >= x) +
  • Vx (the volume irradiated with a dose >= x) +
  • MOHx (mean dose of the hottest x volume) +
  • MOCx (mean dose of the coldest x volume) +
  • MaxOHx (Maximum outside of the hottest x volume) +
  • MinOCx (Minimum outside of the coldest x volume) +
+ Default x values for Vx are 0.02, 0.05, 0.1, 0.9, 0.95 and 0.98, with respect to maxDose. + Default x values for Dx, MOHx, MOCx, MaxOHx and MinOCx are 0.02, 0.05, 0.1, 0.9, 0.95 and 0.98, with respect to volume. + @param computeComplexMeasures should complex statistics be calculated? If it is true, the complex dose statistics will be calculated with default relative x values and the maximum dose as reference dose + @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed + @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed + @warning If computeComplexMeasures==true, computations can take quite long (>1 min) for large structures as many statistics are precomputed + @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! Only the default x values can be requested in DoseStatistics! + */ + DoseStatisticsPointer calculateDoseStatistics(bool computeComplexMeasures = false, unsigned int maxNumberMinimaPositions = 10, + unsigned int maxNumberMaximaPositions = 10); + + /*! @brief Compute complex dose statistics with given reference dose and default relative x values + @param referenceDose the reference dose to compute Vx, normally it should be the prescribed dose + @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed + @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed + @exception InvalidParameterException thrown if referenceDose <= 0 + @warning Computations can take quite long (>1 min) for large structures as many statistics are precomputed + @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! Only the default x values can be requested in DoseStatistics! + */ + DoseStatisticsPointer calculateDoseStatistics(DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions = 10, + unsigned int maxNumberMaximaPositions = 10); + + /*! @brief Compute complex dose statistics with given relative x values and reference dose + @param precomputeDoseValues the relative dose values for Vx precomputation, e.g. 0.02, 0.05, 0.95... + @param precomputeVolumeValues the relative volume values for Dx, MOHx, MOCx, MaxOHx and MinOCx precomputation, e.g. 0.02, 0.05, 0.95... + @param referenceDose the reference dose to compute Vx, normally it should be the prescribed dose. Default value is the maximum dose. + @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed + @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed + @warning Computations can take quite long (>1 min) for large structures as many statistics are precomputed + @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! The doses/volumes that should be used for precomputation have to be set by in precomputeDoseValues and precomputeVolumeValues. Only these values can be requested in DoseStatistics! */ - DoseStatisticsPointer calculateDoseStatistics(bool computeComplexMeasures = false, - const std::vector& precomputeDoseValues = std::vector(), - const std::vector& precomputeVolumeValues = std::vector(), unsigned int maxNumberMinimaPositions = 100, - unsigned int maxNumberMaximaPositions = 100); + DoseStatisticsPointer calculateDoseStatistics(const std::vector& precomputeDoseValues, + const std::vector& precomputeVolumeValues, DoseTypeGy referenceDose = -1, unsigned int maxNumberMinimaPositions = 10, + unsigned int maxNumberMaximaPositions = 10); }; } } #endif diff --git a/code/io/other/rttbDoseStatisticsXMLWriter.cpp b/code/io/other/rttbDoseStatisticsXMLWriter.cpp index b25dbba..7fee2c7 100644 --- a/code/io/other/rttbDoseStatisticsXMLWriter.cpp +++ b/code/io/other/rttbDoseStatisticsXMLWriter.cpp @@ -1,360 +1,271 @@ // ----------------------------------------------------------------------- // 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 "rttbDoseStatisticsXMLWriter.h" #include "rttbInvalidParameterException.h" #include "rttbDataNotAvailableException.h" namespace rttb { namespace io { namespace other { static const std::string xmlattrXTag = ".x"; static const std::string xmlattrNameTag = ".name"; static const std::string statisticsTag = "statistics.results"; static const std::string propertyTag = "property"; static const std::string columnSeparator = "@"; - boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, - DoseTypeGy aReferenceDose) + boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics) { using boost::property_tree::ptree; ptree pt; ptree numberOfVoxelsNode = createNodeWithNameAttribute(aDoseStatistics->getNumberOfVoxels(), "numberOfVoxels"); pt.add_child(statisticsTag + "." + propertyTag, numberOfVoxelsNode); ptree volumeNode = createNodeWithNameAttribute(aDoseStatistics->getVolume(), "volume"); pt.add_child(statisticsTag + "." + propertyTag, volumeNode); ptree minimumNode = createNodeWithNameAttribute(aDoseStatistics->getMinimum(), "minimum"); auto minimumPositions = aDoseStatistics->getMinimumPositions(); std::vector >::iterator pairItMin = minimumPositions->begin(); int count = 0; for (; pairItMin != minimumPositions->end() && count < 100; ++pairItMin) //output max. 100 minimum { VoxelGridID voxelID = pairItMin->second; ptree voxelMinPositions; voxelMinPositions.add("voxelGridID", voxelID); minimumNode.add_child("voxel", voxelMinPositions); count++; } pt.add_child(statisticsTag + "." + propertyTag, minimumNode); ptree maximumNode = createNodeWithNameAttribute(aDoseStatistics->getMaximum(), "maximum"); auto maximumPositions = aDoseStatistics->getMaximumPositions(); std::vector >::iterator pairItMax = maximumPositions->begin(); count = 0; for (; pairItMax != maximumPositions->end() && count < 100; ++pairItMax) //output max. 100 maximum { VoxelGridID voxelID = pairItMax->second; ptree voxelMaxPositions; voxelMaxPositions.add("voxelGridID", voxelID); maximumNode.add_child("voxel", voxelMaxPositions); count++; } pt.add_child(statisticsTag + "." + propertyTag, maximumNode); ptree meanNode = createNodeWithNameAttribute(aDoseStatistics->getMean(), "mean"); pt.add_child(statisticsTag + "." + propertyTag, meanNode); ptree sdNode = createNodeWithNameAttribute(aDoseStatistics->getStdDeviation(), "standardDeviation"); pt.add_child(statisticsTag + "." + propertyTag, sdNode); ptree varianceNode = createNodeWithNameAttribute(aDoseStatistics->getVariance(), "variance"); pt.add_child(statisticsTag + "." + propertyTag, varianceNode); double absoluteVolume = aDoseStatistics->getVolume(); + double referenceDose = aDoseStatistics->getReferenceDose(); + rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType AllVx = aDoseStatistics->getAllVx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllDx = aDoseStatistics->getAllDx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getAllMOHx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getAllMOCx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getAllMaxOHx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getAllMinOCx(); - std::vector DxVxOutputX = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)(0.95)(0.98); - - try - { - for (int i = 0; i < DxVxOutputX.size(); i++) - { - ptree DxNode = createNodeWithNameAndXAttribute(aDoseStatistics->getDx(absoluteVolume * DxVxOutputX.at(i)), "Dx", - static_cast(DxVxOutputX.at(i) * 100)); - pt.add_child(statisticsTag + "." + propertyTag, DxNode); - } - } - catch (core::DataNotAvailableException) - { - //as data is not available (was not computed by doseStatistics), it cannot be written - } + rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType::iterator vxIt; + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType::iterator it; - if (aReferenceDose <= 0) + for (it = AllDx.begin(); it != AllDx.end(); it++) { - throw core::InvalidParameterException("aReferenceDose should be >0!"); + ptree DxNode = createNodeWithNameAndXAttribute((*it).second, "Dx", + static_cast((*it).first / absoluteVolume * 100)); + pt.add_child(statisticsTag + "." + propertyTag, DxNode); } - try - { - for (int i = 0; i < DxVxOutputX.size(); i++) - { - ptree VxNode = createNodeWithNameAndXAttribute(aDoseStatistics->getVx(aReferenceDose * DxVxOutputX.at(i)), "Vx", - static_cast(DxVxOutputX.at(i) * 100)); - pt.add_child(statisticsTag + "." + propertyTag, VxNode); - } - } - catch (core::DataNotAvailableException) + for (vxIt = AllVx.begin(); vxIt != AllVx.end(); vxIt++) { - //as data is not available (was not computed by doseStatistics), it cannot be written + ptree VxNode = createNodeWithNameAndXAttribute((*vxIt).second, "Vx", + static_cast((*vxIt).first / referenceDose * 100)); + pt.add_child(statisticsTag + "." + propertyTag, VxNode); } - std::vector OHOCOutputX = boost::assign::list_of(0.02)(0.05)(0.1); - - try + for (it = AllMOHx.begin(); it != AllMOHx.end(); it++) { - for (int i = 0; i < OHOCOutputX.size(); i++) - { - ptree mohxNode = createNodeWithNameAndXAttribute(aDoseStatistics->getMOHx(absoluteVolume * OHOCOutputX.at(i)), "MOHx", - static_cast(OHOCOutputX.at(i) * 100)); - pt.add_child(statisticsTag + "." + propertyTag, mohxNode); - } - } - catch (core::DataNotAvailableException) - { - //as data is not available (was not computed by doseStatistics), it cannot be written + ptree mohxNode = createNodeWithNameAndXAttribute((*it).second, "MOHx", + static_cast((*it).first / absoluteVolume * 100)); + pt.add_child(statisticsTag + "." + propertyTag, mohxNode); } - try - { - for (int i = 0; i < OHOCOutputX.size(); i++) - { - ptree mocxNode = createNodeWithNameAndXAttribute(aDoseStatistics->getMOCx(absoluteVolume * OHOCOutputX.at(i)), "MOCx", - static_cast(OHOCOutputX.at(i) * 100)); - pt.add_child(statisticsTag + "." + propertyTag, mocxNode); - } - } - catch (core::DataNotAvailableException) + for (it = AllMOCx.begin(); it != AllMOCx.end(); it++) { - //as data is not available (was not computed by doseStatistics), it cannot be written + ptree mocxNode = createNodeWithNameAndXAttribute((*it).second, "MOCx", + static_cast((*it).first / absoluteVolume * 100)); + pt.add_child(statisticsTag + "." + propertyTag, mocxNode); } - try + for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); it++) { - for (int i = 0; i < OHOCOutputX.size(); i++) - { - ptree maxOhxNode = createNodeWithNameAndXAttribute(aDoseStatistics->getMaxOHx(absoluteVolume * OHOCOutputX.at(i)), - "MaxOHx", - static_cast(OHOCOutputX.at(i) * 100)); - pt.add_child(statisticsTag + "." + propertyTag, maxOhxNode); - } - } - catch (core::DataNotAvailableException) - { - //as data is not available (was not computed by doseStatistics), it cannot be written + ptree maxOhxNode = createNodeWithNameAndXAttribute((*it).second, "MaxOHx", + static_cast((*it).first / absoluteVolume * 100)); + pt.add_child(statisticsTag + "." + propertyTag, maxOhxNode); } - try + for (it = AllMinOCx.begin(); it != AllMinOCx.end(); it++) { - for (int i = 0; i < OHOCOutputX.size(); i++) - { - ptree minOCxNode = createNodeWithNameAndXAttribute(aDoseStatistics->getMinOCx(absoluteVolume * OHOCOutputX.at(i)), - "MinOCx", - static_cast(OHOCOutputX.at(i) * 100)); - pt.add_child(statisticsTag + "." + propertyTag, minOCxNode); - } - } - catch (core::DataNotAvailableException) - { - //as data is not available (was not computed by doseStatistics), it cannot be written + ptree minOCxNode = createNodeWithNameAndXAttribute((*it).second, "MinOCx", + static_cast((*it).first / absoluteVolume * 100)); + pt.add_child(statisticsTag + "." + propertyTag, minOCxNode); } return pt; } - void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName, - DoseTypeGy aReferenceDose) + void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName) { - boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics, aReferenceDose); + boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics); try { boost::property_tree::xml_parser::write_xml(aFileName, pt, std::locale(), - boost::property_tree::xml_writer_make_settings('\t', 1)); + boost::property_tree::xml_writer_make_settings('\t', 1)); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml failed: xml_parser_error!"); } } - XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose) + XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics) { - boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics, aReferenceDose); + boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics); std::stringstream sstr; try { boost::property_tree::xml_parser::write_xml(sstr, pt, boost::property_tree::xml_writer_make_settings('\t', - 1)); + 1)); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml to string failed: xml_parser_error!"); } return sstr.str(); } - StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics, - DoseTypeGy aReferenceDose) + StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics) { std::stringstream sstr; sstr << aDoseStatistics->getVolume() * 1000 << columnSeparator; // cm3 to mm3 sstr << aDoseStatistics->getMaximum() << columnSeparator; sstr << aDoseStatistics->getMinimum() << columnSeparator; sstr << aDoseStatistics->getMean() << columnSeparator; sstr << aDoseStatistics->getStdDeviation() << columnSeparator; sstr << aDoseStatistics->getVariance() << columnSeparator; - /*todo: x should be defined based on the user's feedback*/ - if (aReferenceDose <= 0) - { - throw core::InvalidParameterException("aReferenceDose should be >0!"); - } + double absoluteVolume = aDoseStatistics->getVolume(); + rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType AllVx = aDoseStatistics->getAllVx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllDx = aDoseStatistics->getAllDx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getAllMOHx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getAllMOCx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getAllMaxOHx(); + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getAllMinOCx(); - std::vector DxVxOutputX = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)(0.95)(0.98); - double absoluteVolume = aDoseStatistics->getVolume(); + rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType::iterator vxIt; + rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType::iterator it; - try - { - for (int i = 0; i < DxVxOutputX.size(); i++) - { - sstr << aDoseStatistics->getDx(absoluteVolume * DxVxOutputX.at(i)) << columnSeparator; - } - } - catch (core::DataNotAvailableException) + for (it = AllDx.begin(); it != AllDx.end(); it++) { - //as data is not available (was not computed by doseStatistics), it cannot be written + sstr << (*it).second << columnSeparator; } - try - { - for (int i = 0; i < DxVxOutputX.size(); i++) - { - // *1000 because of conversion cm3 to mm3 - sstr << aDoseStatistics->getVx(aReferenceDose * DxVxOutputX.at(i)) * 1000 << columnSeparator; - } - } - catch (core::DataNotAvailableException) + for (vxIt = AllVx.begin(); vxIt != AllVx.end(); vxIt++) { - //as data is not available (was not computed by doseStatistics), it cannot be written + // *1000 because of conversion cm3 to mm3 + sstr << (*vxIt).second * 1000 << columnSeparator; } - std::vector OHOCOutputX = boost::assign::list_of(0.02)(0.05)(0.1); - - try + for (it = AllMOHx.begin(); it != AllMOHx.end(); it++) { - for (int i = 0; i < OHOCOutputX.size(); i++) - { - sstr << aDoseStatistics->getMOHx(absoluteVolume * OHOCOutputX.at(i)) << columnSeparator; - } - } - catch (core::DataNotAvailableException) - { - //as data is not available (was not computed by doseStatistics), it cannot be written + sstr << (*it).second << columnSeparator; } - try - { - for (int i = 0; i < OHOCOutputX.size(); i++) - { - sstr << aDoseStatistics->getMOCx(absoluteVolume * OHOCOutputX.at(i)) << columnSeparator; - } - } - catch (core::DataNotAvailableException) + for (it = AllMOCx.begin(); it != AllMOCx.end(); it++) { - //as data is not available (was not computed by doseStatistics), it cannot be written + sstr << (*it).second << columnSeparator; } - try + for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); it++) { - for (int i = 0; i < OHOCOutputX.size(); i++) - { - sstr << aDoseStatistics->getMaxOHx(absoluteVolume * OHOCOutputX.at(i)) << columnSeparator; - } - } - catch (core::DataNotAvailableException) - { - //as data is not available (was not computed by doseStatistics), it cannot be written + sstr << (*it).second << columnSeparator; } - try + for (it = AllMinOCx.begin(); it != AllMinOCx.end(); it++) { - for (int i = 0; i < OHOCOutputX.size(); i++) - { - sstr << aDoseStatistics->getMinOCx(absoluteVolume * OHOCOutputX.at(i)) << columnSeparator; - } + sstr << (*it).second << columnSeparator; } - catch (core::DataNotAvailableException) - { - //as data is not available (was not computed by doseStatistics), it cannot be written - } - return sstr.str(); } boost::property_tree::ptree createNodeWithNameAttribute(DoseTypeGy doseValue, const std::string& attributeName) { boost::property_tree::ptree node; node.put("", doseValue); node.put(xmlattrNameTag, attributeName); return node; } boost::property_tree::ptree createNodeWithNameAndXAttribute(DoseTypeGy doseValue, const std::string& attributeName, - int xValue) + int xValue) { boost::property_tree::ptree node; node.put("", doseValue); node.put(xmlattrNameTag, attributeName); node.put(xmlattrXTag, xValue); return node; } }//end namespace other }//end namespace io }//end namespace rttb \ No newline at end of file diff --git a/code/io/other/rttbDoseStatisticsXMLWriter.h b/code/io/other/rttbDoseStatisticsXMLWriter.h index 7fceb9d..02a0c99 100644 --- a/code/io/other/rttbDoseStatisticsXMLWriter.h +++ b/code/io/other/rttbDoseStatisticsXMLWriter.h @@ -1,97 +1,97 @@ // ----------------------------------------------------------------------- // 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 { 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); + boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics); /*! @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); + XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics); /*! @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); + void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName); boost::property_tree::ptree createNodeWithNameAttribute(DoseTypeGy doseValue, const std::string& attributeName); boost::property_tree::ptree createNodeWithNameAndXAttribute(DoseTypeGy doseValue, const std::string& attributeName, int xValue); /*! @brief Write statistics to String to generate a table: "Volume mm3@Max@Min@Mean@Std.Dev.@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 @pre aReferenceDose must >0 @exception InvalidParameterException Thrown if aReferenceDose<=0 or xml_parse_error @note is used for the Mevislab-Linking of RTTB */ - StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose); + StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics); } } } #endif 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 b65361f..379bce9 100644 --- a/testing/algorithms/DoseStatisticsCalculatorTest.cpp +++ b/testing/algorithms/DoseStatisticsCalculatorTest.cpp @@ -1,296 +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 values for omputeComplexMeasures=true + //check manually set reference dose and the default x values + CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(100.0)); + CHECK_THROW_EXPLICIT(theStatistics->getVx(0.1 * theStatistics->getMaximum()), core::DataNotAvailableException); + CHECK_NO_THROW(theStatistics->getVx(0.1 * 100.0)); + CHECK_NO_THROW(theStatistics->getDx(0.1 * theStatistics->getVolume())); + CHECK_NO_THROW(theStatistics->getDx(0.9 * theStatistics->getVolume())); + CHECK_NO_THROW(theStatistics->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(true, precomputeDoseValues, + 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_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 diff --git a/testing/examples/DoseStatisticsIOVirtuosTest.cpp b/testing/examples/DoseStatisticsIOVirtuosTest.cpp index ba92311..79f4ac0 100644 --- a/testing/examples/DoseStatisticsIOVirtuosTest.cpp +++ b/testing/examples/DoseStatisticsIOVirtuosTest.cpp @@ -1,113 +1,113 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDoseStatistics.h" #include "rttbDoseStatisticsCalculator.h" #include "rttbDoseStatisticsXMLWriter.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbDoseIteratorInterface.h" #include "rttbVirtuosPlanFileDoseAccessorGenerator.h" #include "rttbVirtuosDoseAccessor.h" #include "rttbVirtuosFileStructureSetGenerator.h" #include "rttbBoostMaskAccessor.h" namespace rttb { namespace testing { /*! @brief OtherIOTest - test the IO for dose statistics 1) test exception 2) test writing statistics to xml file */ int DoseStatisticsIOVirtuosTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; typedef boost::shared_ptr DoseStatisticsPtr; typedef core::MaskAccessorInterface::MaskAccessorPointer MaskAccessorPointer; typedef boost::shared_ptr GeometricInfoPointer; PREPARE_DEFAULT_TEST_REPORTING; std::string Struct_FILENAME; std::string BPLCT_FILENAME; std::string Dose_FILENAME; std::string Plan_FILENAME; std::string Struct_NAME; if (argc > 5) { Struct_FILENAME = argv[1]; BPLCT_FILENAME = argv[2]; Dose_FILENAME = argv[3]; Plan_FILENAME = argv[4]; Struct_NAME = argv[5]; } else { std::cout << "at least 5 arguments are expected in DoseStatisticsIOVirtuosTest." << std::endl; } /* generate dummy dose */ io::virtuos::VirtuosPlanFileDoseAccessorGenerator doseAccessorGenerator(Dose_FILENAME.c_str(), Plan_FILENAME.c_str()); DoseAccessorPointer doseAccessor(doseAccessorGenerator.generateDoseAccessor()); StructureSetPointer rtStructureSet = io::virtuos::VirtuosFileStructureSetGenerator( Struct_FILENAME.c_str(), BPLCT_FILENAME.c_str()).generateStructureSet(); GeometricInfoPointer geometricInfoPtr = boost::make_shared(doseAccessor->getGeometricInfo()); MaskAccessorPointer maskAccessorPtr = boost::make_shared (rtStructureSet->getStructure(3), doseAccessor->getGeometricInfo()); maskAccessorPtr->updateMask(); boost::shared_ptr spTestDoseIterator = boost::make_shared(maskAccessorPtr, doseAccessor); DoseIteratorPointer spDoseIterator(spTestDoseIterator); rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator); rttb::algorithms::DoseStatistics::DoseStatisticsPointer doseStatistics = myDoseStatsCalculator.calculateDoseStatistics(true); /* test writing statistcs to xml file */ FileNameString fN = "doseStatisticsVirtuosComplex.xml"; - CHECK_NO_THROW(io::other::writeDoseStatistics(doseStatistics, fN, 100)); + CHECK_NO_THROW(io::other::writeDoseStatistics(doseStatistics, fN)); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/examples/RTDoseStatisticsTest.cpp b/testing/examples/RTDoseStatisticsTest.cpp index 1545bd0..82d8ff6 100644 --- a/testing/examples/RTDoseStatisticsTest.cpp +++ b/testing/examples/RTDoseStatisticsTest.cpp @@ -1,202 +1,204 @@ // ----------------------------------------------------------------------- // 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 "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 "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; + const double volume = 24120.; 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[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericDoseIterator::DoseIteratorPointer DoseIteratorPointer; typedef algorithms::DoseStatistics::ResultListPointer ResultListPointer; 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); RETURN_AND_REPORT_TEST_SUCCESS; } void testWithDummyDoseData(const std::string& doseFilename) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericDoseIterator::DoseIteratorPointer DoseIteratorPointer; typedef ::boost::shared_ptr > > ResultsVectorPointer; 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); rttb::algorithms::DoseStatisticsCalculator doseStatisticsCalculator(spDoseIterator); std::vector precomputedDoseValues; precomputedDoseValues.push_back(0); - precomputedDoseValues.push_back(expectedVal); - precomputedDoseValues.push_back(expectedVal + reducedErrorConstant); + precomputedDoseValues.push_back(0.5); + precomputedDoseValues.push_back(1); std::vector precomputedVolumeValues; - precomputedVolumeValues.push_back(20000); - precomputedVolumeValues.push_back(24120); + precomputedVolumeValues.push_back(20000/volume); + precomputedVolumeValues.push_back(1); rttb::algorithms::DoseStatistics::DoseStatisticsPointer doseStatistics = - doseStatisticsCalculator.calculateDoseStatistics(true, precomputedDoseValues, precomputedVolumeValues); + doseStatisticsCalculator.calculateDoseStatistics(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); + double dummy; + DoseTypeGy vx = doseStatistics->getVx(expectedVal,true,dummy); 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_EQUAL(maxListPtr->size(), 10); + CHECK_EQUAL(minListPtr->size(), 10); 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; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef rttb::algorithms::DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer; auto virtuosDoseAccessor = io::virtuos::VirtuosPlanFileDoseAccessorGenerator(doseFilename.c_str(), planFilename.c_str()).generateDoseAccessor(); auto virtuosStructureSet = io::virtuos::VirtuosFileStructureSetGenerator( structFilename.c_str(), doseFilename.c_str()).generateStructureSet(); boost::shared_ptr spOTBMaskAccessorVirtuos = 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(minPositions->size(), 10); 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); } }//testing }//rttb \ No newline at end of file diff --git a/testing/io/other/DoseStatisticsIOTest.cpp b/testing/io/other/DoseStatisticsIOTest.cpp index fcc5bfe..3a68872 100644 --- a/testing/io/other/DoseStatisticsIOTest.cpp +++ b/testing/io/other/DoseStatisticsIOTest.cpp @@ -1,109 +1,104 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDoseStatisticsCalculator.h" #include "rttbDoseStatisticsXMLWriter.h" #include "rttbGenericDoseIterator.h" #include "rttbDoseIteratorInterface.h" #include "rttbInvalidParameterException.h" #include "rttbNullPointerException.h" #include "../../core/DummyDoseAccessor.h" namespace rttb { namespace testing { /*! @brief OtherIOTest - test the IO for dose statistics 1) test exception 2) test writing statistics to xml file */ int DoseStatisticsIOTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; typedef boost::shared_ptr DoseStatisticsCalculatorPtr; PREPARE_DEFAULT_TEST_REPORTING; /* generate dummy dose */ boost::shared_ptr spTestDoseAccessor = boost::make_shared(); DoseAccessorPointer spDoseAccessor(spTestDoseAccessor); boost::shared_ptr spTestDoseIterator = boost::make_shared(spDoseAccessor); DoseIteratorPointer spDoseIterator(spTestDoseIterator); rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator); auto myDoseStatsSimple = myDoseStatsCalculator.calculateDoseStatistics(); auto myDoseStatsComplex = myDoseStatsCalculator.calculateDoseStatistics(true); - /* test exception */ - CHECK_THROW_EXPLICIT(io::other::writeDoseStatistics(myDoseStatsSimple, "test.test", 0), - core::InvalidParameterException); - - /* test writing statistics to xml file */ FileNameString filenameSimple = "testStatisticsSimple.xml"; - CHECK_NO_THROW(io::other::writeDoseStatistics(myDoseStatsSimple, filenameSimple, 100)); + CHECK_NO_THROW(io::other::writeDoseStatistics(myDoseStatsSimple, filenameSimple)); FileNameString filenameComplex = "testStatisticsComplex.xml"; - CHECK_NO_THROW(io::other::writeDoseStatistics(myDoseStatsComplex, filenameComplex, 100)); + CHECK_NO_THROW(io::other::writeDoseStatistics(myDoseStatsComplex, filenameComplex)); /* test writing statistics to string */ - boost::property_tree::ptree ptSimple = io::other::writeDoseStatistics(myDoseStatsSimple, 100); - XMLString strSimple = io::other::writerDoseStatisticsToString(myDoseStatsSimple, 100); + boost::property_tree::ptree ptSimple = io::other::writeDoseStatistics(myDoseStatsSimple); + XMLString strSimple = io::other::writerDoseStatisticsToString(myDoseStatsSimple); - boost::property_tree::ptree ptComplex = io::other::writeDoseStatistics(myDoseStatsComplex, 100); - XMLString strComplex = io::other::writerDoseStatisticsToString(myDoseStatsComplex, 100); + boost::property_tree::ptree ptComplex = io::other::writeDoseStatistics(myDoseStatsComplex); + XMLString strComplex = io::other::writerDoseStatisticsToString(myDoseStatsComplex); std::stringstream sstrSimple; boost::property_tree::xml_parser::write_xml(sstrSimple, ptSimple, boost::property_tree::xml_writer_make_settings('\t', 1)); CHECK_EQUAL(strSimple, sstrSimple.str()); std::stringstream sstrComplex; boost::property_tree::xml_parser::write_xml(sstrComplex, ptComplex, boost::property_tree::xml_writer_make_settings('\t', 1)); CHECK_EQUAL(strComplex, sstrComplex.str()); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb