diff --git a/code/algorithms/rttbDoseStatistics.cpp b/code/algorithms/rttbDoseStatistics.cpp index 14fe157..ac8f926 100644 --- a/code/algorithms/rttbDoseStatistics.cpp +++ b/code/algorithms/rttbDoseStatistics.cpp @@ -1,314 +1,315 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "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()*/, DoseToVolumeFunctionType Dx /*= std::map()*/, VolumeToDoseFunctionType Vx /*= std::map()*/, VolumeToDoseFunctionType MOHx /*= std::map()*/, VolumeToDoseFunctionType MOCx /*= std::map()*/, VolumeToDoseFunctionType MaxOHx /*= std::map()*/, VolumeToDoseFunctionType MinOCx /*= std::map()*/) : _minimum(minimum), _maximum(maximum), _mean(mean), _stdDeviation(stdDeviation), _numVoxels(numVoxels), _volume(volume), _maximumVoxelPositions(maximumVoxelPositions), _minimumVoxelPositions(minimumVoxelPositions), _Dx(Dx), _Vx(Vx), _MOHx(MOHx), _MOCx(MOCx), _MaxOHx(MaxOHx), _MinOCx(MinOCx) { } DoseStatistics::~DoseStatistics() { } void DoseStatistics::setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions) { _minimumVoxelPositions = minimumVoxelPositions; } void DoseStatistics::setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions) { _maximumVoxelPositions = maximumVoxelPositions; } void DoseStatistics::setDx(const DoseToVolumeFunctionType& DxValues) { _Dx = DxValues; } void DoseStatistics::setVx(const VolumeToDoseFunctionType& VxValues) { _Vx = VxValues; } void DoseStatistics::setMOHx(const VolumeToDoseFunctionType& MOHxValues) { _MOHx = MOHxValues; } void DoseStatistics::setMOCx(const VolumeToDoseFunctionType& MOCxValues) { _MOCx = MOCxValues; } void DoseStatistics::setMaxOHx(const VolumeToDoseFunctionType& MaxOHValues) { _MaxOHx = MaxOHValues; } void DoseStatistics::setMinOCx(const VolumeToDoseFunctionType& MinOCValues) { _MinOCx = MinOCValues; } unsigned int DoseStatistics::getNumberOfVoxels() const { return _numVoxels; } VolumeType DoseStatistics::getVolume() const { return _volume; } DoseStatisticType DoseStatistics::getMaximum() const { return _maximum; } DoseStatisticType DoseStatistics::getMinimum() const { return _minimum; } DoseStatisticType DoseStatistics::getMean() const { return _mean; } DoseStatisticType DoseStatistics::getStdDeviation() const { return _stdDeviation; } DoseStatisticType DoseStatistics::getVariance() const { return _stdDeviation * _stdDeviation; } VolumeType DoseStatistics::getVx(DoseTypeGy xDoseAbsolute, bool findNearestValue, DoseTypeGy& nearestXDose) const { return getValue(_Vx, xDoseAbsolute, findNearestValue, nearestXDose); } VolumeType DoseStatistics::getVx(DoseTypeGy xDoseAbsolute) const { DoseTypeGy dummy; return getValue(_Vx, xDoseAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getDx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const { + return getValue(_Dx, xVolumeAbsolute, findNearestValue, nearestXVolume); } DoseTypeGy DoseStatistics::getDx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_Dx, xVolumeAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getMOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& nearestXVolume) const { return getValue(_MOHx, xVolumeAbsolute, findNearestValue, nearestXVolume); } DoseTypeGy DoseStatistics::getMOHx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_MOHx, xVolumeAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getMOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { return getValue(_MOCx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } DoseTypeGy DoseStatistics::getMOCx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_MOCx, xVolumeAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getMaxOHx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { return getValue(_MaxOHx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } DoseTypeGy DoseStatistics::getMaxOHx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_MaxOHx, xVolumeAbsolute, false, dummy); } DoseTypeGy DoseStatistics::getMinOCx(VolumeType xVolumeAbsolute, bool findNearestValue, VolumeType& maximumDistanceFromOriginalVolume) const { return getValue(_MinOCx, xVolumeAbsolute, findNearestValue, maximumDistanceFromOriginalVolume); } DoseTypeGy DoseStatistics::getMinOCx(VolumeType xVolumeAbsolute) const { VolumeType dummy; return getValue(_MinOCx, xVolumeAbsolute, false, dummy); } double DoseStatistics::getValue(std::map aMap, double key, bool findNearestValueInstead, double& storedKey) const { if (aMap.find(key) != aMap.end()) { return aMap.find(key)->second; } else { //value not in map. We have to find the nearest value if (aMap.empty()) { throw core::DataNotAvailableException("No Vx values are defined"); } else { if (findNearestValueInstead) { auto iterator = findNearestKeyInMap(aMap, key); storedKey = iterator->first; return iterator->second; } else { throw core::DataNotAvailableException("No Vx value with required dose is defined"); } } } } std::map::const_iterator DoseStatistics::findNearestKeyInMap(std::map aMap, 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 5a497cc..6003a1e 100644 --- a/code/algorithms/rttbDoseStatistics.h +++ b/code/algorithms/rttbDoseStatistics.h @@ -1,205 +1,205 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __DOSE_STATISTICS_H #define __DOSE_STATISTICS_H #include #include #include #include "rttbBaseType.h" #include "rttbDoseIteratorInterface.h" namespace rttb { namespace algorithms { /*! @class DoseStatistics @brief This is a data class storing different statistical values from a rt dose distribution @sa DoseStatisticsCalculator */ class DoseStatistics { public: enum complexStatistics { Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx}; typedef boost::shared_ptr > > ResultListPointer; typedef boost::shared_ptr DoseStatisticsPointer; typedef std::map DoseToVolumeFunctionType; typedef std::map VolumeToDoseFunctionType; private: double getValue(std::map aMap, double key, bool findNearestValueInstead, double& storedKey) const; std::map::const_iterator findNearestKeyInMap(std::map aMap, double key) const; DoseStatisticType _maximum; DoseStatisticType _minimum; ResultListPointer _maximumVoxelPositions; ResultListPointer _minimumVoxelPositions; DoseStatisticType _mean; DoseStatisticType _stdDeviation; unsigned int _numVoxels; VolumeType _volume; - DoseToVolumeFunctionType _Dx; - VolumeToDoseFunctionType _Vx; + 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 = boost::make_shared > > (std::vector >()), ResultListPointer maximumVoxelPositions = boost::make_shared > > (std::vector >()), DoseToVolumeFunctionType Dx = DoseToVolumeFunctionType(), VolumeToDoseFunctionType Vx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MOHx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MOCx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MaxOHx = VolumeToDoseFunctionType(), VolumeToDoseFunctionType MinOCx = VolumeToDoseFunctionType()); ~DoseStatistics(); void setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions); void setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions); void setDx(const DoseToVolumeFunctionType& DxValues); void setVx(const VolumeToDoseFunctionType& VxValues); void setMOHx(const VolumeToDoseFunctionType& MOHxValues); void setMOCx(const VolumeToDoseFunctionType& MOCxValues); void setMaxOHx(const VolumeToDoseFunctionType& MaxOHxValues); void setMinOCx(const VolumeToDoseFunctionType& MinOCxValues); /*! @brief Get number of voxels in doseIterator, with sub-voxel accuracy. */ unsigned int getNumberOfVoxels() const; VolumeType getVolume() 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 2039fcf..b2de8f2 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,554 +1,564 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 877 $ (last changed revision) // @date $Date: 2015-01-09 10:51:10 +0100 (Fr, 09 Jan 2015) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #include "rttbDoseStatisticsCalculator.h" #include #include #include #include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace algorithms { DoseStatisticsCalculator::DoseStatisticsCalculator() { simpleDoseStatisticsCalculated = false; _doseIterator = NULL; } DoseStatisticsCalculator::DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator) { setDoseIterator(aDoseIterator); simpleDoseStatisticsCalculated = false; } DoseStatisticsCalculator::~DoseStatisticsCalculator() { } void DoseStatisticsCalculator::setDoseIterator(DoseIteratorPointer aDoseIterator) { if (aDoseIterator == NULL) { throw core::NullPointerException("DoseIterator must not be NULL"); } else { _doseIterator = aDoseIterator; } } DoseStatisticsCalculator::DoseIteratorPointer DoseStatisticsCalculator::getDoseIterator() const { return _doseIterator; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( bool computeComplexMeasures, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } //"simple" dose statistics are obligatory calculateSimpleDoseStatistics(); if (computeComplexMeasures) { //more complex dose statistics are optional calculateComplexDoseStatistics(precomputeDoseValues, precomputeVolumeValues); } return _statistics; } void DoseStatisticsCalculator::calculateSimpleDoseStatistics() { _doseVector.clear(); _voxelProportionVector.clear(); std::multimap doseValueVSIndexMap; std::vector voxelProportionVectorTemp; DoseStatisticType maximumDose = 0; DoseStatisticType minimumDose = 0; DoseStatisticType meanDose = 0; DoseStatisticType stdDeviationDose = 0; float sum = 0; - unsigned int numVoxels = 0; + rttb::DoseStatisticType numVoxels = 0.0; float squareSum = 0; - VolumeType volume; + 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; - DoseStatisticType varianceDose = (squareSum / numVoxels - meanDose * meanDose); //standard deviation is defined only for n>=2 - if (varianceDose < errorConstant || numVoxels <= 2) + if (numVoxels >= 2) { - stdDeviationDose = 0; + //uncorrected variance is calculated + DoseStatisticType varianceDose = (squareSum / numVoxels - meanDose * meanDose); + + if (varianceDose < errorConstant) + { + stdDeviationDose = 0; + } + else + { + stdDeviationDose = pow(varianceDose, 0.5); + } } else { - stdDeviationDose = pow(varianceDose, 0.5); + stdDeviationDose = 0; } } //sort dose values and corresponding volume fractions in member variables std::multimap::iterator it; for (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(100); ResultListPointer maximumVoxelPositions = computeMaximumPositions(100); _statistics->setMinimumVoxelPositions(minimumVoxelPositions); _statistics->setMaximumVoxelPositions(maximumVoxelPositions); } void DoseStatisticsCalculator::calculateComplexDoseStatistics(const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues) { std::vector precomputeDoseValuesNonConst = precomputeDoseValues; std::vector precomputeVolumeValuesNonConst = precomputeVolumeValues; //set default values if (precomputeDoseValues.empty()) { - std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02 * _statistics->getMaximum())( - 0.05 * _statistics->getMaximum())(0.1 * _statistics->getMaximum())(0.9 * _statistics->getMaximum())( - 0.95 * _statistics->getMaximum())(0.98 * _statistics->getMaximum()); + DoseTypeGy maxDose = _statistics->getMaximum(); + std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02 * maxDose)(0.05 * maxDose)(0.1 * maxDose)( + 0.9 * maxDose)( + 0.95 * maxDose)(0.98 * maxDose); precomputeDoseValuesNonConst = defaultPrecomputeDoseValues; } if (precomputeVolumeValues.empty()) { - std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02 * _statistics->getNumberOfVoxels())( - 0.05 * _statistics->getNumberOfVoxels())(0.1 * _statistics->getNumberOfVoxels())(0.9 * _statistics->getNumberOfVoxels()) - (0.95 * _statistics->getNumberOfVoxels())(0.98 * _statistics->getNumberOfVoxels()); + VolumeType volume = _statistics->getVolume(); + std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02 * volume)( + 0.05 * volume)(0.1 * volume)(0.9 * volume) + (0.95 * volume)(0.98 * volume); precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues; } DoseToVolumeFunctionType Vx = computeDoseToVolumeMulti(precomputeDoseValuesNonConst, DoseStatistics::Vx); VolumeToDoseFunctionType Dx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::Dx); VolumeToDoseFunctionType MOHx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MOHx); VolumeToDoseFunctionType MOCx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MOCx); VolumeToDoseFunctionType MaxOHx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MaxOHx); VolumeToDoseFunctionType MinOCx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MinOCx); _statistics->setVx(Vx); _statistics->setDx(Dx); _statistics->setMOHx(MOHx); _statistics->setMOCx(MOCx); _statistics->setMaxOHx(MaxOHx); _statistics->setMinOCx(MinOCx); } DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMaximumPositions( unsigned int maxNumberMaxima) { if (!simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMaximumAndPosition()"); } 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) { if (!simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMinimumAndPosition()"); } ResultListPointer minVoxelVector = boost::make_shared > >(); /*! @todo: Architecture Annotation: Finding the positions for the minimum only once reduces computation time, but will require sensible use by the programmers. To be save the output vector minVoxelVector will be always cleared here to garantee that no false values are presented. This change may be revoced to increase computation speed later on (only compute if(minVoxelVector->size()==0)). */ unsigned int count = 0; this->_doseIterator->reset(); DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid() && count < maxNumberMinima) { doseValue = _doseIterator->getCurrentDoseValue(); if (doseValue == _statistics->getMinimum()) { VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID(); std::pair voxel(doseValue, currentID); minVoxelVector->push_back(voxel); count++; } _doseIterator->next(); } return minVoxelVector; } VolumeType DoseStatisticsCalculator::computeVx(DoseTypeGy xDoseAbsolute) const { rttb::FractionType count = 0; _doseIterator->reset(); DoseTypeGy currentDose = 0; while (_doseIterator->isPositionValid()) { currentDose = _doseIterator->getCurrentDoseValue(); if (currentDose >= xDoseAbsolute) { count += _doseIterator->getCurrentRelevantVolumeFraction(); } _doseIterator->next(); } return count * this->_doseIterator->getCurrentVoxelVolume(); } DoseTypeGy DoseStatisticsCalculator::computeDx(VolumeType xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; int i = _doseVector.size() - 1; for (; i >= 0; i--) { countVoxels += _voxelProportionVector.at(i); if (countVoxels >= noOfVoxel) { break; } } if (i >= 0) { resultDose = _doseVector.at(i); } else { resultDose = _statistics->getMinimum(); } return resultDose; } DoseTypeGy DoseStatisticsCalculator::computeMOHx(VolumeType xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); if (noOfVoxel == 0) { return 0; } else { double countVoxels = 0; double sum = 0; for (int i = _doseVector.size() - 1; i >= 0; i--) { double voxelProportion = _voxelProportionVector.at(i); countVoxels += voxelProportion; sum += _doseVector.at(i) * voxelProportion; if (countVoxels >= noOfVoxel) { break; } } return (DoseTypeGy)(sum / noOfVoxel); } } DoseTypeGy DoseStatisticsCalculator::computeMOCx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); if (noOfVoxel == 0) { return 0; } else { double countVoxels = 0; double sum = 0; std::vector::const_iterator it = _doseVector.begin(); std::vector::const_iterator itD = _voxelProportionVector.begin(); for (; it != _doseVector.end(); ++it, ++itD) { double voxelProportion = *itD; countVoxels += voxelProportion; sum += (*it) * voxelProportion; if (countVoxels >= noOfVoxel) { break; } } return (DoseTypeGy)(sum / noOfVoxel); } } DoseTypeGy DoseStatisticsCalculator::computeMaxOHx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; int i = _doseVector.size() - 1; for (; i >= 0; i--) { countVoxels += _voxelProportionVector.at(i); if (countVoxels >= noOfVoxel) { break; } } if (i - 1 >= 0) { resultDose = _doseVector.at(i - 1); } return resultDose; } DoseTypeGy DoseStatisticsCalculator::computeMinOCx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; std::vector::const_iterator it = _doseVector.begin(); std::vector::const_iterator itD = _voxelProportionVector.begin(); for (; itD != _voxelProportionVector.end(); ++itD, ++it) { countVoxels += *itD; if (countVoxels >= noOfVoxel) { break; } } if (it != _doseVector.end()) { ++it; if (it != _doseVector.end()) { resultDose = *it; } else { resultDose = (DoseTypeGy)_statistics->getMaximum(); } } else { resultDose = (DoseTypeGy)_statistics->getMinimum(); } return resultDose; } DoseStatisticsCalculator::DoseToVolumeFunctionType DoseStatisticsCalculator::computeDoseToVolumeMulti( const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const { DoseToVolumeFunctionType VxMulti; - DoseTypeGy maxDose = _statistics->getMaximum(); for (int i = 0; i < precomputeDoseValues.size(); ++i) { if (name == DoseStatistics::Vx) { VxMulti.insert(std::pair(precomputeDoseValues.at(i), - computeVx(maxDose * (precomputeDoseValues.at(i) / 100.0)))); + computeVx(precomputeDoseValues.at(i)))); } } return VxMulti; } DoseStatisticsCalculator::VolumeToDoseFunctionType DoseStatisticsCalculator::computeVolumeToDoseFunctionMulti( const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const { VolumeToDoseFunctionType multiValues; - VolumeType maxVolume = _statistics->getNumberOfVoxels(); for (int i = 0; i < precomputeVolumeValues.size(); ++i) { switch (name) { case DoseStatistics::Dx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeDx(maxVolume * (precomputeVolumeValues.at(i) / 100.0)))); + computeDx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MOHx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeMOHx(maxVolume * (precomputeVolumeValues.at(i) / 100.0)))); + computeMOHx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MOCx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeMOCx(maxVolume * (precomputeVolumeValues.at(i) / 100.0)))); + computeMOCx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MaxOHx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeMaxOHx(maxVolume * (precomputeVolumeValues.at(i) / 100.0)))); + computeMaxOHx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MinOCx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), - computeMinOCx(maxVolume * (precomputeVolumeValues.at(i) / 100.0)))); + computeMinOCx(precomputeVolumeValues.at(i)))); break; default: throw core::InvalidParameterException("unknown DoseStatistics name!"); } } return multiValues; } }//end namespace algorithms }//end namespace rttb diff --git a/testing/algorithms/DoseStatisticsCalculatorTest.cpp b/testing/algorithms/DoseStatisticsCalculatorTest.cpp index 507ce7c..1a083d8 100644 --- a/testing/algorithms/DoseStatisticsCalculatorTest.cpp +++ b/testing/algorithms/DoseStatisticsCalculatorTest.cpp @@ -1,314 +1,313 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 707 $ (last changed revision) // @date $Date: 2014-09-04 16:37:24 +0200 (Do, 04 Sep 2014) $ (last change date) // @author $Author: floca $ (last changed by) */ #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbGenericDoseIterator.h" #include "rttbDoseIteratorInterface.h" #include "rttbNullPointerException.h" #include "rttbDoseStatisticsCalculator.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" #include "rttbDataNotAvailableException.h" #include "../core/DummyDoseAccessor.h" namespace rttb { namespace testing { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; typedef rttb::algorithms::DoseStatistics::ResultListPointer ResultListPointer; typedef rttb::algorithms::DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer; /*! @brief DoseStatisticsCalculatorTest - test the API of DoseStatisticsCalculator 1) test constructors 2) test setDoseIterator 3) test calculateDoseSatistics 4) get statistical values */ int DoseStatisticsCalculatorTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; boost::shared_ptr spTestDoseAccessor = boost::make_shared(); DoseAccessorPointer spDoseAccessor(spTestDoseAccessor); const std::vector* doseVals = spTestDoseAccessor->getDoseVector(); boost::shared_ptr spTestDoseIterator = boost::make_shared(spDoseAccessor); DoseIteratorPointer spDoseIterator(spTestDoseIterator); DoseIteratorPointer spDoseIteratorNull; //1) test constructors // the values cannot be accessed from outside, therefore correct default values are not tested CHECK_NO_THROW(rttb::algorithms::DoseStatisticsCalculator myEmptyDoseStatCalculatur); rttb::algorithms::DoseStatisticsCalculator myEmptyDoseStatCalculatur; rttb::algorithms::DoseStatisticsCalculator myDoseStatCalculaturToBeFilled; CHECK_THROW_EXPLICIT(rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIteratorNull), core::NullPointerException); CHECK_NO_THROW(rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator)); rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator); //2) test setDoseIterator CHECK_THROW_EXPLICIT(myDoseStatCalculaturToBeFilled.setDoseIterator(spDoseIteratorNull), core::NullPointerException); CHECK_NO_THROW(myDoseStatCalculaturToBeFilled.setDoseIterator(spDoseIterator)); //3) test calculateDoseStatistics DoseStatisticsPointer theStatistics, theStatistics2; CHECK_THROW_EXPLICIT(myEmptyDoseStatCalculatur.calculateDoseStatistics(), core::NullPointerException); CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics()); CHECK_NO_THROW(theStatistics2 = myDoseStatCalculaturToBeFilled.calculateDoseStatistics()); - CHECK_EQUAL(theStatistics->getMinimumPositions()->empty(), true); - CHECK_EQUAL(theStatistics->getMaximumPositions()->empty(), true); + 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.2 * theStatisticsDefault->getMaximum())); - CHECK_NO_THROW(theStatisticsDefault->getVx(0.8 * theStatisticsDefault->getMaximum())); CHECK_NO_THROW(theStatisticsDefault->getVx(0.9 * theStatisticsDefault->getMaximum())); CHECK_NO_THROW(theStatisticsDefault->getVx(0.95 * theStatisticsDefault->getMaximum())); CHECK_NO_THROW(theStatisticsDefault->getVx(0.98 * theStatisticsDefault->getMaximum())); CHECK_NO_THROW(theStatisticsDefault->getDx(0.02 * theStatisticsDefault->getNumberOfVoxels())); CHECK_NO_THROW(theStatisticsDefault->getDx(0.05 * theStatisticsDefault->getNumberOfVoxels())); CHECK_NO_THROW(theStatisticsDefault->getDx(0.1 * theStatisticsDefault->getNumberOfVoxels())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.2 * theStatisticsDefault->getNumberOfVoxels())); - CHECK_NO_THROW(theStatisticsDefault->getDx(0.8 * theStatisticsDefault->getNumberOfVoxels())); CHECK_NO_THROW(theStatisticsDefault->getDx(0.9 * theStatisticsDefault->getNumberOfVoxels())); CHECK_NO_THROW(theStatisticsDefault->getDx(0.95 * theStatisticsDefault->getNumberOfVoxels())); CHECK_NO_THROW(theStatisticsDefault->getDx(0.98 * theStatisticsDefault->getNumberOfVoxels())); //check manually set values for omputeComplexMeasures=true std::vector precomputeDoseValues, precomputeVolumeValues; precomputeDoseValues.push_back(0.01 * theStatistics->getMaximum()); precomputeDoseValues.push_back(0.02 * theStatistics->getMaximum()); precomputeDoseValues.push_back(0.05 * theStatistics->getMaximum()); precomputeVolumeValues.push_back(0.9 * theStatistics->getNumberOfVoxels()); precomputeVolumeValues.push_back(0.95 * theStatistics->getNumberOfVoxels()); precomputeVolumeValues.push_back(0.99 * theStatistics->getNumberOfVoxels()); CHECK_NO_THROW(theStatistics = myDoseStatsCalculator.calculateDoseStatistics(true, precomputeDoseValues, precomputeVolumeValues)); CHECK_NO_THROW(theStatistics->getVx(0.01 * theStatistics->getMaximum())); CHECK_NO_THROW(theStatistics->getVx(0.02 * theStatistics->getMaximum())); CHECK_NO_THROW(theStatistics->getVx(0.05 * theStatistics->getMaximum())); CHECK_THROW_EXPLICIT(theStatistics->getVx(0.03 * theStatistics->getMaximum()), core::DataNotAvailableException); CHECK_NO_THROW(theStatistics->getDx(0.9 * theStatistics->getNumberOfVoxels())); CHECK_NO_THROW(theStatistics->getDx(0.95 * theStatistics->getNumberOfVoxels())); CHECK_NO_THROW(theStatistics->getDx(0.99 * theStatistics->getNumberOfVoxels())); CHECK_THROW_EXPLICIT(theStatistics->getDx(0.03 * theStatistics->getNumberOfVoxels()), core::DataNotAvailableException); CHECK_EQUAL(theStatistics->getVx(0.02 * theStatistics->getMaximum()), theStatisticsDefault->getVx(0.02 * theStatistics->getMaximum())); CHECK_EQUAL(theStatistics->getVx(0.05 * theStatistics->getMaximum()), theStatisticsDefault->getVx(0.05 * theStatistics->getMaximum())); CHECK_EQUAL(theStatistics->getDx(0.9 * theStatistics->getNumberOfVoxels()), theStatisticsDefault->getDx(0.9 * theStatistics->getNumberOfVoxels())); CHECK_EQUAL(theStatistics->getDx(0.95 * theStatistics->getNumberOfVoxels()), theStatisticsDefault->getDx(0.95 * theStatistics->getNumberOfVoxels())); //MOHx, MOCx, MaxOHx and MinOCx are computed analogous to Dx, they will not be checked. //4) get statistical values CHECK_EQUAL(theStatistics->getNumberOfVoxels(), doseVals->size()); CHECK_EQUAL(theStatistics2->getNumberOfVoxels(), doseVals->size()); //compute simple statistical values (min, mean, max, stddev) for comparison DoseStatisticType maximum = 0; DoseStatisticType minimum = 1000000; DoseStatisticType mean = 0; DoseStatisticType variance = 0; std::vector::const_iterator doseIt = doseVals->begin(); while (doseIt != doseVals->end()) { if (maximum < *doseIt) { maximum = *doseIt; } if (minimum > *doseIt) { minimum = *doseIt; } mean += *doseIt; ++doseIt; } mean /= doseVals->size(); float compMean = (int(mean * 100)) / 100; doseIt = doseVals->begin(); while (doseIt != doseVals->end()) { variance += pow(*doseIt - mean, 2); ++doseIt; } - variance /= doseVals->size() - 1; + variance /= doseVals->size(); DoseStatisticType stdDev = pow(variance, 0.5); + //we have some precision problems here... + double errorConstantLarger = 1e-2; + CHECK_EQUAL(theStatistics->getMaximum(), maximum); CHECK_EQUAL(theStatistics2->getMaximum(), maximum); CHECK_EQUAL(theStatistics->getMinimum(), minimum); CHECK_EQUAL(theStatistics2->getMinimum(), minimum); - CHECK_EQUAL(theStatistics->getMean(), mean); - CHECK_EQUAL(theStatistics2->getMean(), mean); + CHECK_CLOSE(theStatistics->getMean(), mean, errorConstantLarger); + CHECK_CLOSE(theStatistics2->getMean(), mean, errorConstantLarger); - CHECK_EQUAL(theStatistics->getStdDeviation(), stdDev); - CHECK_EQUAL(theStatistics2->getStdDeviation(), stdDev); + CHECK_CLOSE(theStatistics->getStdDeviation(), stdDev, errorConstantLarger); + CHECK_CLOSE(theStatistics2->getStdDeviation(), stdDev, errorConstantLarger); - CHECK_EQUAL(theStatistics->getVariance(), variance); - CHECK_EQUAL(theStatistics2->getVariance(), variance); + CHECK_CLOSE(theStatistics->getVariance(), variance, errorConstantLarger); + CHECK_CLOSE(theStatistics2->getVariance(), variance, errorConstantLarger); //check for complex doseStatistics (maximumPositions, minimumPositions, Vx, Dx, MOHx, MOCx, MAXOHx, MinOCx) unsigned int nMax = 0, nMin = 0; doseIt = doseVals->begin(); while (doseIt != doseVals->end()) { if (*doseIt == theStatistics->getMaximum()) { nMax++; } if (*doseIt == theStatistics->getMinimum()) { nMin++; } ++doseIt; } //only 100 positions are stored if (nMax > 100) { nMax = 100; } if (nMin > 100) { nMin = 100; } auto maximaPositions = theStatistics->getMaximumPositions(); auto minimaPositions = theStatistics->getMinimumPositions(); CHECK_EQUAL(maximaPositions->size(), nMax); CHECK_EQUAL(minimaPositions->size(), nMin); for (auto maximaPositionsIterator = std::begin(*maximaPositions); maximaPositionsIterator != std::end(*maximaPositions); ++maximaPositionsIterator) { CHECK_EQUAL(maximaPositionsIterator->first, theStatistics->getMaximum()); } for (auto minimaPositionsIterator = std::begin(*minimaPositions); minimaPositionsIterator != std::end(*minimaPositions); ++minimaPositionsIterator) { CHECK_EQUAL(minimaPositionsIterator->first, theStatistics->getMinimum()); } //generate specific example dose maximum = 9.5; minimum = 2.5; mean = 6; int sizeTemplate = 500; std::vector aDoseVector; for (int i = 0; i < sizeTemplate; i++) { aDoseVector.push_back(maximum); aDoseVector.push_back(minimum); } core::GeometricInfo geoInfo = spTestDoseAccessor->getGeometricInfo(); geoInfo.setNumRows(20); geoInfo.setNumColumns(10); geoInfo.setNumSlices(5); boost::shared_ptr spTestDoseAccessor2 = boost::make_shared(aDoseVector, geoInfo); DoseAccessorPointer spDoseAccessor2(spTestDoseAccessor2); boost::shared_ptr spTestDoseIterator2 = boost::make_shared(spDoseAccessor2); DoseIteratorPointer spDoseIterator2(spTestDoseIterator2); rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator2(spDoseIterator2); DoseStatisticsPointer theStatistics3 = myDoseStatsCalculator2.calculateDoseStatistics(); CHECK_EQUAL(theStatistics3->getMaximum(), maximum); CHECK_EQUAL(theStatistics3->getMinimum(), minimum); CHECK_EQUAL(theStatistics3->getMean(), mean); maximaPositions = theStatistics3->getMaximumPositions(); minimaPositions = theStatistics3->getMinimumPositions(); CHECK_EQUAL(maximaPositions->empty(), false); CHECK_EQUAL(minimaPositions->empty(), false); for (auto maximaPositionsIterator = std::begin(*maximaPositions); maximaPositionsIterator != std::end(*maximaPositions); ++maximaPositionsIterator) { CHECK_EQUAL(maximaPositionsIterator->first, theStatistics3->getMaximum()); } for (auto minimaPositionsIterator = std::begin(*minimaPositions); minimaPositionsIterator != std::end(*minimaPositions); ++minimaPositionsIterator) { CHECK_EQUAL(minimaPositionsIterator->first, theStatistics3->getMinimum()); } RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/examples/CMakeLists.txt b/testing/examples/CMakeLists.txt index 1b3c46a..32d2d2b 100644 --- a/testing/examples/CMakeLists.txt +++ b/testing/examples/CMakeLists.txt @@ -1,37 +1,39 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(CORE_TEST_EXAMPLES ${EXECUTABLE_OUTPUT_PATH}/rttbExamplesTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(RTBioModelExampleTest ${CORE_TEST_EXAMPLES} RTBioModelExampleTest "${TEST_DATA_ROOT}/TestDVH/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT3.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_ah/dvh_diff_trunk6.txt" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_ah/dvh_diff_trunk8.txt") ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwoGridScaling.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwoGridScaling05.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantFiftyGridScaling01.dcm") ADD_TEST(RTDoseIndexTest ${CORE_TEST_EXAMPLES} RTDoseIndexTest "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT3.txt") ADD_TEST(RTDoseStatisticsTest ${CORE_TEST_EXAMPLES} RTDoseStatisticsTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo_withDoseGridScaling.dcm" -"${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm") +"${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa101.dos.gz" +"${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa000.vdx" +"${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa101.pln") ADD_TEST(RTDVHTest ${CORE_TEST_EXAMPLES} RTDVHTest "${TEST_DATA_ROOT}/TestDVH/dvh_test.txt") ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo_withDoseGridScaling.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm") ADD_TEST(RTBioModelScatterPlotExampleTest ${CORE_TEST_EXAMPLES} RTBioModelScatterPlotExampleTest "${TEST_DATA_ROOT}/TestDVH/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt") -RTTB_CREATE_TEST_MODULE(rttbExamples DEPENDS RTTBCore RTTBAlgorithms RTTBMasks RTTBIndices RTTBDicomIO RTTBOtherIO RTTBModels PACKAGE_DEPENDS BoostBinaries Litmus DCMTK) +RTTB_CREATE_TEST_MODULE(rttbExamples DEPENDS RTTBCore RTTBAlgorithms RTTBMasks RTTBIndices RTTBDicomIO RTTBVirtuosIO RTTBOtherIO RTTBModels PACKAGE_DEPENDS BoostBinaries Litmus DCMTK) diff --git a/testing/examples/RTDoseStatisticsTest.cpp b/testing/examples/RTDoseStatisticsTest.cpp index 9b363ae..5282d4f 100644 --- a/testing/examples/RTDoseStatisticsTest.cpp +++ b/testing/examples/RTDoseStatisticsTest.cpp @@ -1,132 +1,200 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include "dcmtk/config/osconfig.h" /* make sure OS specific configuration is included first */ #include "dcmtk/dcmrt/drtdose.h" #include "dcmtk/dcmdata/dcfilefo.h" #include "litCheckMacros.h" #include "rttbDoseStatistics.h" +#include "rttbDoseStatisticsCalculator.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" +#include "rttbVirtuosPlanFileDoseAccessorGenerator.h" +#include "rttbVirtuosFileStructureSetGenerator.h" +#include "rttbOTBMaskAccessor.h" +#include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbNullPointerException.h" #include "rttbBaseType.h" namespace rttb { namespace testing { /*! @brief RTDoseStatisticsTest. Max, min, mean, standard deviation, variance, Vx, Dx, MOHx, MOCx, MaxOHx, MinOCx are tested. Test if calculation in new architecture returns similar results to the original implementation. WARNING: The values for comparison need to be adjusted if the input files are changed!*/ int RTDoseStatisticsTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericDoseIterator::DoseIteratorPointer DoseIteratorPointer; typedef algorithms::DoseStatistics::ResultListPointer ResultListPointer; //ARGUMENTS: 1: dose1 file name // 2: dose2 file name + // 3: RTStruct filename std::string RTDOSE_FILENAME; std::string RTDOSE2_FILENAME; + std::string RTSTRUCT_FILENAME; + std::string RTPLAN_FILENAME; - if (argc > 1) + if (argc > 4) { RTDOSE_FILENAME = argv[1]; + RTDOSE2_FILENAME = argv[2]; + RTSTRUCT_FILENAME = argv[3]; + RTPLAN_FILENAME = argv[4]; } - - if (argc > 2) + else { - RTDOSE2_FILENAME = argv[2]; + std::cout << "at least four arguments required for RTDoseStatisticsTest" << std::endl; + return -1; } OFCondition status; DcmFileFormat fileformat; - double max, min; const double expectedVal = 5.64477e-005; const double reducedErrorConstant = 0.0001; - /*Test NullPointerException*/ - - rttb::algorithms::DoseStatistics doseStatistics; + rttb::algorithms::DoseStatisticsCalculator doseStatisticsCalculator; typedef boost::shared_ptr > > ResultsVectorPointer; - ResultsVectorPointer spResults = - boost::make_shared > >(); - - ResultListPointer minListPtr(spResults); - ResultListPointer maxListPtr(spResults); - ::DRTDoseIOD rtdoseDKFZ; io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create corresponding DoseIterator boost::shared_ptr spDoseIteratorTmp = boost::make_shared(doseAccessor1); DoseIteratorPointer spDoseIterator(spDoseIteratorTmp); - doseStatistics.setDoseIterator(spDoseIterator); + doseStatisticsCalculator.setDoseIterator(spDoseIterator); - CHECK_CLOSE(doseStatistics.getMean(), expectedVal, errorConstant); - CHECK_CLOSE(doseStatistics.getStdDeviation(), 0, errorConstant); - CHECK_CLOSE(doseStatistics.getVariance(), 0, errorConstant); + std::vector precomputedDoseValues; + precomputedDoseValues.push_back(0); + precomputedDoseValues.push_back(expectedVal); + precomputedDoseValues.push_back(expectedVal + reducedErrorConstant); + std::vector precomputedVolumeValues; + precomputedVolumeValues.push_back(20000); + precomputedVolumeValues.push_back(24120); - DoseTypeGy vx = doseStatistics.getVx(expectedVal); - CHECK_EQUAL(vx, doseStatistics.getVx(0)); - CHECK_CLOSE(expectedVal, doseStatistics.getDx(vx), reducedErrorConstant); + rttb::algorithms::DoseStatistics::DoseStatisticsPointer doseStatistics = + doseStatisticsCalculator.calculateDoseStatistics(true, precomputedDoseValues, precomputedVolumeValues); - max = doseStatistics.getMaximum(maxListPtr); - CHECK_CLOSE(doseStatistics.getMaximum(maxListPtr), expectedVal, errorConstant); - CHECK_EQUAL(doseStatistics.getVx(max + reducedErrorConstant), 0); + CHECK_CLOSE(doseStatistics->getMean(), expectedVal, errorConstant); + CHECK_CLOSE(doseStatistics->getStdDeviation(), 0, errorConstant); + CHECK_CLOSE(doseStatistics->getVariance(), 0, errorConstant); - min = doseStatistics.getMinimum(minListPtr); - CHECK_CLOSE(min, expectedVal, errorConstant); + DoseTypeGy vx = doseStatistics->getVx(expectedVal); + CHECK_EQUAL(vx, doseStatistics->getVx(0)); + CHECK_CLOSE(expectedVal, doseStatistics->getDx(vx), reducedErrorConstant); + + CHECK_CLOSE(doseStatistics->getMaximum(), expectedVal, errorConstant); + + CHECK_CLOSE(doseStatistics->getMinimum(), expectedVal, errorConstant); + auto minListPtr = doseStatistics->getMinimumPositions(); + auto maxListPtr = doseStatistics->getMaximumPositions(); + CHECK_EQUAL(maxListPtr->size(), 100); CHECK_EQUAL(minListPtr->size(), 100); - min = doseStatistics.getMinimum(minListPtr, 50); - CHECK_EQUAL(minListPtr->size(), 50); - - DoseTypeGy wholeVolume = doseStatistics.getVx(min); - CHECK_CLOSE(doseStatistics.getDx(wholeVolume), min, 0.001); - CHECK_CLOSE(doseStatistics.getMOHx(vx), doseStatistics.getMean(), reducedErrorConstant); - CHECK_CLOSE(doseStatistics.getMOCx(20000), doseStatistics.getMean(), reducedErrorConstant); - CHECK_CLOSE(doseStatistics.getMinOCx(20000), doseStatistics.getMean(), reducedErrorConstant); + + 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); + + //test with real data set (virtuos) + typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; + typedef rttb::algorithms::DoseStatistics::DoseStatisticsPointer DoseStatisticsPointer; + + auto virtuosDoseAccessor = io::virtuos::VirtuosPlanFileDoseAccessorGenerator(RTDOSE2_FILENAME.c_str(), + RTPLAN_FILENAME.c_str()).generateDoseAccessor(); + auto virtuosStructureSet = io::virtuos::VirtuosFileStructureSetGenerator( + RTSTRUCT_FILENAME.c_str(), RTDOSE2_FILENAME.c_str()).generateStructureSet(); + + //Structure 2 is RUECKENMARK + boost::shared_ptr spOTBMaskAccessorVirtuos = + boost::make_shared(virtuosStructureSet->getStructure(2), + virtuosDoseAccessor->getGeometricInfo()); + + spOTBMaskAccessorVirtuos->updateMask(); + + MaskAccessorPointer spMaskAccessor(spOTBMaskAccessorVirtuos); + + //create corresponding MaskedDoseIterator + + boost::shared_ptr spMaskedDoseIteratorTmp = + boost::make_shared(spMaskAccessor, virtuosDoseAccessor); + + DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); + + rttb::algorithms::DoseStatisticsCalculator doseStatisticsCalculatorVirtuos(spMaskedDoseIterator); + + + DoseStatisticsPointer doseStatisticsVirtuos = doseStatisticsCalculatorVirtuos.calculateDoseStatistics(true); + + //comparison values computed with "old" DoseStatistics implementation + CHECK_EQUAL(doseStatisticsVirtuos->getMinimum(), 0); + CHECK_CLOSE(doseStatisticsVirtuos->getMaximum(), 35.3075, reducedErrorConstant); + CHECK_CLOSE(doseStatisticsVirtuos->getMean(), 16.1803, reducedErrorConstant); + CHECK_CLOSE(doseStatisticsVirtuos->getStdDeviation(), 9.84436, reducedErrorConstant); + + auto maxPositions = doseStatisticsVirtuos->getMaximumPositions(); + auto minPositions = doseStatisticsVirtuos->getMinimumPositions(); + CHECK_EQUAL(maxPositions->size(), 1); + CHECK_EQUAL(minPositions->size(), 100); + CHECK_EQUAL(maxPositions->begin()->first, doseStatisticsVirtuos->getMaximum()); + CHECK_EQUAL(maxPositions->begin()->second, 2138227); + CHECK_EQUAL(minPositions->begin()->first, doseStatisticsVirtuos->getMinimum()); + CHECK_EQUAL(minPositions->begin()->second, 39026); + + CHECK_CLOSE(doseStatisticsVirtuos->getDx(0.02 * doseStatisticsVirtuos->getVolume()), 30.1528, reducedErrorConstant); + CHECK_CLOSE(doseStatisticsVirtuos->getVx(0.9 * doseStatisticsVirtuos->getMaximum()), 0.483529, + reducedErrorConstant); + CHECK_CLOSE(doseStatisticsVirtuos->getMOHx(0.1 * doseStatisticsVirtuos->getVolume()), 28.7019, + reducedErrorConstant); + CHECK_CLOSE(doseStatisticsVirtuos->getMOCx(0.05 * doseStatisticsVirtuos->getVolume()), 0, + reducedErrorConstant); + CHECK_CLOSE(doseStatisticsVirtuos->getMaxOHx(0.95 * doseStatisticsVirtuos->getVolume()), 0, + reducedErrorConstant); + CHECK_CLOSE(doseStatisticsVirtuos->getMinOCx(0.98 * doseStatisticsVirtuos->getVolume()), 30.1756, + reducedErrorConstant); + std::cout << "-->" << 0.1 * doseStatisticsVirtuos->getVolume() << std::endl; RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb \ No newline at end of file diff --git a/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.cpp b/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.cpp index 9ea6c4a..b8810f4 100644 --- a/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.cpp +++ b/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.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: 741 $ (last changed revision) // @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #undef MAP_SEAL_ALGORITHMS #include "simpleRegistrationWorkflow.h" #include "registrationHelper.h" simpleRegistrationWorkflow::simpleRegistrationWorkflow(std::string targetFilename, std::string movingFilename, bool isDirectory) { _targetFilename = targetFilename; _movingFilename = movingFilename; setImageFileNames(_targetFilename, _movingFilename, isDirectory, globals); loadData(globals); _spAlgorithmEuler = NULL; } vnl_vector simpleRegistrationWorkflow::getRegistrationParameters( Registration3D3DPointer reg) { typedef map::core::ModelBasedRegistrationKernel<3, 3> ModelBasedRegistrationKernel3D3D; const ModelBasedRegistrationKernel3D3D* pModelBasedDirectKernel3D3D = dynamic_cast(&(reg->getDirectMapping())); if (pModelBasedDirectKernel3D3D) { ModelBasedRegistrationKernel3D3D::ParametersType params = - pModelBasedDirectKernel3D3D->getTransformModel()->getTransform()->GetParameters(); + pModelBasedDirectKernel3D3D->getTransformModel()->GetParameters(); return params; } else { return vnl_vector(); } } void simpleRegistrationWorkflow::initializeAndPerformRegistration() { _spAlgorithmEuler = AlgorithmTypeEuler::New(); _spAlgorithmEuler->setProperty("PreinitTransform", map::core::MetaProperty::New(true)); _spAlgorithmEuler->setMovingImage(globals.spMovingImage); _spAlgorithmEuler->setTargetImage(globals.spTargetImage); AlgorithmTypeEuler::RegistrationType::Pointer spRegistration; try { spRegistration = _spAlgorithmEuler->getRegistration(); } catch (const map::core::ExceptionObject& e) { std::cerr << "caught an MatchPoint exception:\n"; e.Print(std::cerr); std::cerr << "\n"; } catch (const itk::ExceptionObject& e) { std::cerr << "caught an ITK exception:\n"; std::cerr << e.GetFile() << ":" << e.GetLine() << ":\n" << e.GetDescription() << "\n"; } catch (const std::exception& e) { std::cerr << "caught an exception:\n"; std::cerr << e.what() << "\n"; } catch (...) { std::cerr << "caught an unknown exception!!!\n"; } } map::core::Registration<3, 3>::Pointer simpleRegistrationWorkflow::getRegistration() { if (_spAlgorithmEuler.IsNull()) { initializeAndPerformRegistration(); } return _spAlgorithmEuler->getRegistration(); }; const itk::Image* simpleRegistrationWorkflow::getTargetImage() { if (_spAlgorithmEuler.IsNull()) { initializeAndPerformRegistration(); } return _spAlgorithmEuler->getTargetImage(); };