diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp index 84bed4e..7c56b94 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,554 +1,556 @@ // ----------------------------------------------------------------------- // 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(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) + const std::vector& precomputeVolumeValues, unsigned int maxNumberMinimaPositions, + unsigned int maxNumberMaximaPositions) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } //"simple" dose statistics are mandatory - calculateSimpleDoseStatistics(); + calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); if (computeComplexMeasures) { //more complex dose statistics are optional calculateComplexDoseStatistics(precomputeDoseValues, precomputeVolumeValues); } return _statistics; } - void DoseStatisticsCalculator::calculateSimpleDoseStatistics() + 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(100); - ResultListPointer maximumVoxelPositions = computeMaximumPositions(100); + ResultListPointer minimumVoxelPositions = computeMinimumPositions(maxNumberMinimaPositions); + ResultListPointer maximumVoxelPositions = computeMaximumPositions(maxNumberMaximaPositions); _statistics->setMinimumVoxelPositions(minimumVoxelPositions); _statistics->setMaximumVoxelPositions(maximumVoxelPositions); } void DoseStatisticsCalculator::calculateComplexDoseStatistics(const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues) { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call calculateComplexDoseStatistics()"); } std::vector precomputeDoseValuesNonConst = precomputeDoseValues; std::vector precomputeVolumeValuesNonConst = precomputeVolumeValues; //set default values if (precomputeDoseValues.empty()) { DoseTypeGy maxDose = _statistics->getMaximum(); std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02 * maxDose)(0.05 * maxDose)(0.1 * maxDose)( 0.9 * maxDose)( 0.95 * maxDose)(0.98 * maxDose); precomputeDoseValuesNonConst = defaultPrecomputeDoseValues; } if (precomputeVolumeValues.empty()) { VolumeType volume = _statistics->getVolume(); std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02 * volume)( 0.05 * volume)(0.1 * volume)(0.9 * volume) (0.95 * volume)(0.98 * volume); precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues; } - DoseToVolumeFunctionType Vx = computeDoseToVolumeMulti(precomputeDoseValuesNonConst, DoseStatistics::Vx); + DoseToVolumeFunctionType Vx = computeDoseToVolumeFunctionMulti(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) 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::computeDoseToVolumeMulti( + DoseStatisticsCalculator::DoseToVolumeFunctionType DoseStatisticsCalculator::computeDoseToVolumeFunctionMulti( const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const { DoseToVolumeFunctionType VxMulti; for (int i = 0; i < precomputeDoseValues.size(); ++i) { if (name == DoseStatistics::Vx) { VxMulti.insert(std::pair(precomputeDoseValues.at(i), computeVx(precomputeDoseValues.at(i)))); } } return VxMulti; } DoseStatisticsCalculator::VolumeToDoseFunctionType DoseStatisticsCalculator::computeVolumeToDoseFunctionMulti( const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const { VolumeToDoseFunctionType multiValues; for (int i = 0; i < precomputeVolumeValues.size(); ++i) { switch (name) { case DoseStatistics::Dx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeDx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MOHx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeMOHx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MOCx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeMOCx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MaxOHx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeMaxOHx(precomputeVolumeValues.at(i)))); break; case DoseStatistics::MinOCx: multiValues.insert(std::pair(precomputeVolumeValues.at(i), computeMinOCx(precomputeVolumeValues.at(i)))); break; default: throw core::InvalidParameterException("unknown DoseStatistics name!"); } } return multiValues; } }//end namespace algorithms }//end namespace rttb diff --git a/code/algorithms/rttbDoseStatisticsCalculator.h b/code/algorithms/rttbDoseStatisticsCalculator.h index 91d64da..ede836b 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.h +++ b/code/algorithms/rttbDoseStatisticsCalculator.h @@ -1,141 +1,154 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 707 $ (last changed revision) // @date $Date: 2014-09-04 16:37:24 +0200 (Do, 04 Sep 2014) $ (last change date) // @author $Author: floca $ (last changed by) */ #ifndef __DOSE_STATISTICS_CALCULATOR_H #define __DOSE_STATISTICS_CALCULATOR_H #include #include "rttbBaseType.h" #include "rttbDoseIteratorInterface.h" #include "rttbDoseStatistics.h" namespace rttb { namespace algorithms { /*! @class DoseStatisticsCalculator @brief 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; - ResultListPointer computeMaximumPositions(unsigned int maxNumberMaxima) const; - ResultListPointer computeMinimumPositions(unsigned int maxNumberMinima) const; + /*! @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 computeDoseToVolumeMulti(const std::vector& precomputeDoseValues, + DoseToVolumeFunctionType computeDoseToVolumeFunctionMulti(const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const; VolumeToDoseFunctionType computeVolumeToDoseFunctionMulti(const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const; - /*! @brief Calculatues simple dose statistics (min, mean, max, stdDev, minDosePositions, maxDosePositions) + /*! @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(); - /*! @brief Calculatues complex dose statistics (Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx) + 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, 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! */ DoseStatisticsPointer calculateDoseStatistics(bool computeComplexMeasures = false, const std::vector& precomputeDoseValues = std::vector(), - const std::vector& precomputeVolumeValues = std::vector()); + const std::vector& precomputeVolumeValues = std::vector(), unsigned int maxNumberMinimaPositions = 100, + unsigned int maxNumberMaximaPositions = 100); }; } } #endif diff --git a/code/io/other/rttbDoseStatisticsXMLWriter.cpp b/code/io/other/rttbDoseStatisticsXMLWriter.cpp index 39bf48e..a3da094 100644 --- a/code/io/other/rttbDoseStatisticsXMLWriter.cpp +++ b/code/io/other/rttbDoseStatisticsXMLWriter.cpp @@ -1,464 +1,401 @@ // ----------------------------------------------------------------------- // 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 "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbDataNotAvailableException.h" -#include "rttbException.h" - -#include "boost/lexical_cast.hpp" namespace rttb { namespace io { namespace other { static const std::string xmlattrTag = ".x"; static const std::string statisticsTag = "statistics"; static const std::string columnSeparator = "@"; boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose) { using boost::property_tree::ptree; ptree pt; - /*rttb::core::DoseIteratorInterface::DoseIteratorPointer spDoseIterator = - aDoseStatistics->getDoseIterator(); - rttb::core::DoseIteratorInterface::DoseAccessorPointer spDoseAccessor = - spDoseIterator->getDoseAccessor(); - - boost::shared_ptr< std::vector > > myResultPairs = - boost::make_shared< std::vector > >(); - ResultListPointer spMyResultPairs(myResultPairs); - boost::shared_ptr< std::vector > > myResultPairs2 = - boost::make_shared< std::vector > >(); - ResultListPointer spMyResultPairs2(myResultPairs2);*/ - pt.put(statisticsTag + ".numberOfVoxels", aDoseStatistics->getNumberOfVoxels()); pt.put(statisticsTag + ".volume", aDoseStatistics->getVolume()); pt.put(statisticsTag + ".minimum", aDoseStatistics->getMinimum()); 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; pt.add(statisticsTag + ".minimum.voxelGridID", voxelID); - //VoxelGridIndex3D voxelIndex3D; - //spDoseAccessor->getGeometricInfo().convert(voxelID, voxelIndex3D); - //std::string voxelIndex3DStr; - //voxelIndex3DStr = boost::lexical_cast(voxelIndex3D.x()) + "," - // + boost::lexical_cast(voxelIndex3D.y()) + "," - // + boost::lexical_cast(voxelIndex3D.z()); - //pt.add(statisticsTag + ".minimum.voxelIndex3D", voxelIndex3DStr); - - //WorldCoordinate3D worldCoor; - //spDoseAccessor->getGeometricInfo().indexToWorldCoordinate(voxelIndex3D, worldCoor); - //std::string worldCoorStr; - //worldCoorStr = boost::lexical_cast(worldCoor.x()) + "," - // + boost::lexical_cast(worldCoor.y()) + "," - // + boost::lexical_cast(worldCoor.z()); - //pt.add(statisticsTag + ".minimum.worldCoordinate", worldCoorStr); - count++; } pt.put(statisticsTag + ".maximum", aDoseStatistics->getMaximum()); 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; pt.add(statisticsTag + ".maximum.voxelGridID", voxelID); - //VoxelGridIndex3D voxelIndex3D; - //spDoseAccessor->getGeometricInfo().convert(voxelID, voxelIndex3D); - //std::string voxelIndex3DStr; - //voxelIndex3DStr = boost::lexical_cast(voxelIndex3D.x()) + "," - // + boost::lexical_cast(voxelIndex3D.y()) + "," - // + boost::lexical_cast(voxelIndex3D.z()); - //pt.add(statisticsTag + ".maximum.voxelIndex3D", voxelIndex3DStr); - - - //WorldCoordinate3D worldCoor; - //spDoseAccessor->getGeometricInfo().indexToWorldCoordinate(voxelIndex3D, worldCoor); - //std::string worldCoorStr; - //worldCoorStr = boost::lexical_cast(worldCoor.x()) + "," - // + boost::lexical_cast(worldCoor.y()) + "," - // + boost::lexical_cast(worldCoor.z()); - //pt.add(statisticsTag + ".maximum.worldCoordinate", worldCoorStr); - count ++; } pt.put(statisticsTag + ".mean", aDoseStatistics->getMean()); pt.put(statisticsTag + ".standardDeviation", aDoseStatistics->getStdDeviation()); pt.put(statisticsTag + ".variance", aDoseStatistics->getVariance()); double absoluteVolume = aDoseStatistics->getVolume(); try { //Dx boost::property_tree::ptree dxNode1; boost::property_tree::ptree dxNode2; boost::property_tree::ptree dxNode3; boost::property_tree::ptree dxNode4; boost::property_tree::ptree dxNode5; boost::property_tree::ptree dxNode6; dxNode1.put("", aDoseStatistics->getDx(absoluteVolume * 0.02)); dxNode2.put("", aDoseStatistics->getDx(absoluteVolume * 0.05)); dxNode3.put("", aDoseStatistics->getDx(absoluteVolume * 0.10)); dxNode4.put("", aDoseStatistics->getDx(absoluteVolume * 0.90)); dxNode5.put("", aDoseStatistics->getDx(absoluteVolume * 0.95)); dxNode6.put("", aDoseStatistics->getDx(absoluteVolume * 0.98)); dxNode1.put(xmlattrTag, 2); dxNode2.put(xmlattrTag, 5); dxNode3.put(xmlattrTag, 10); dxNode4.put(xmlattrTag, 90); dxNode5.put(xmlattrTag, 95); dxNode6.put(xmlattrTag, 98); pt.add_child(statisticsTag + ".Dx", dxNode1); pt.add_child(statisticsTag + ".Dx", dxNode2); pt.add_child(statisticsTag + ".Dx", dxNode3); pt.add_child(statisticsTag + ".Dx", dxNode4); pt.add_child(statisticsTag + ".Dx", dxNode5); pt.add_child(statisticsTag + ".Dx", dxNode6); } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } if (aReferenceDose <= 0) { throw core::InvalidParameterException("aReferenceDose should be >0!"); } try { //Vx boost::property_tree::ptree vxNode1; boost::property_tree::ptree vxNode2; boost::property_tree::ptree vxNode3; boost::property_tree::ptree vxNode4; boost::property_tree::ptree vxNode5; boost::property_tree::ptree vxNode6; vxNode1.put("", aDoseStatistics->getVx(aReferenceDose * 0.02)); vxNode2.put("", aDoseStatistics->getVx(aReferenceDose * 0.05)); vxNode3.put("", aDoseStatistics->getVx(aReferenceDose * 0.10)); vxNode4.put("", aDoseStatistics->getVx(aReferenceDose * 0.90)); vxNode5.put("", aDoseStatistics->getVx(aReferenceDose * 0.95)); vxNode6.put("", aDoseStatistics->getVx(aReferenceDose * 0.98)); vxNode1.put(xmlattrTag, 2); vxNode2.put(xmlattrTag, 5); vxNode3.put(xmlattrTag, 10); vxNode4.put(xmlattrTag, 90); vxNode5.put(xmlattrTag, 95); vxNode6.put(xmlattrTag, 98); pt.add_child(statisticsTag + ".Vx", vxNode1); pt.add_child(statisticsTag + ".Vx", vxNode2); pt.add_child(statisticsTag + ".Vx", vxNode3); pt.add_child(statisticsTag + ".Vx", vxNode4); pt.add_child(statisticsTag + ".Vx", vxNode5); pt.add_child(statisticsTag + ".Vx", vxNode6); } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } try { //MOHx boost::property_tree::ptree mohxNode1; boost::property_tree::ptree mohxNode2; boost::property_tree::ptree mohxNode3; mohxNode1.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.02)); mohxNode2.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.05)); mohxNode3.put("", aDoseStatistics->getMOHx(absoluteVolume * 0.10)); mohxNode1.put(xmlattrTag, 2); mohxNode2.put(xmlattrTag, 5); mohxNode3.put(xmlattrTag, 10); pt.add_child(statisticsTag + ".MOHx", mohxNode1); pt.add_child(statisticsTag + ".MOHx", mohxNode2); pt.add_child(statisticsTag + ".MOHx", mohxNode3); } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } try { //MOCx boost::property_tree::ptree mocxNode1; boost::property_tree::ptree mocxNode2; boost::property_tree::ptree mocxNode3; mocxNode1.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.02)); mocxNode2.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.05)); mocxNode3.put("", aDoseStatistics->getMOCx(absoluteVolume * 0.10)); mocxNode1.put(xmlattrTag, 2); mocxNode2.put(xmlattrTag, 5); mocxNode3.put(xmlattrTag, 10); pt.add_child(statisticsTag + ".MOCx", mocxNode1); pt.add_child(statisticsTag + ".MOCx", mocxNode2); pt.add_child(statisticsTag + ".MOCx", mocxNode3); } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } try { //MaxOHx boost::property_tree::ptree maxohxNode1; boost::property_tree::ptree maxohxNode2; boost::property_tree::ptree maxohxNode3; maxohxNode1.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.02)); maxohxNode2.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.05)); maxohxNode3.put("", aDoseStatistics->getMaxOHx(absoluteVolume * 0.10)); maxohxNode1.put(xmlattrTag, 2); maxohxNode2.put(xmlattrTag, 5); maxohxNode3.put(xmlattrTag, 10); pt.add_child(statisticsTag + ".MaxOHx", maxohxNode1); pt.add_child(statisticsTag + ".MaxOHx", maxohxNode2); pt.add_child(statisticsTag + ".MaxOHx", maxohxNode3); } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } try { //MinOCx boost::property_tree::ptree minocxNode1; boost::property_tree::ptree minocxNode2; boost::property_tree::ptree minocxNode3; minocxNode1.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.02)); minocxNode2.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.05)); minocxNode3.put("", aDoseStatistics->getMinOCx(absoluteVolume * 0.10)); minocxNode1.put(xmlattrTag, 2); minocxNode2.put(xmlattrTag, 5); minocxNode3.put(xmlattrTag, 10); pt.add_child(statisticsTag + ".MinOCx", minocxNode1); pt.add_child(statisticsTag + ".MinOCx", minocxNode2); pt.add_child(statisticsTag + ".MinOCx", minocxNode3); } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } return pt; } void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName, DoseTypeGy aReferenceDose) { boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics, aReferenceDose); try { boost::property_tree::xml_parser::write_xml(aFileName, pt); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml failed: xml_parser_error!"); } } XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics, DoseTypeGy aReferenceDose) { boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics, aReferenceDose); std::stringstream sstr; try { boost::property_tree::xml_parser::write_xml(sstr, pt); } 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) { std::stringstream sstr; - - //rttb::core::DoseIteratorInterface::DoseIteratorPointer spDoseIterator = - // aDoseStatistics->getDoseIterator(); - //rttb::core::DoseIteratorInterface::DoseAccessorPointer spDoseAccessor = - // spDoseIterator->getDoseAccessor(); - - //boost::shared_ptr< std::vector > > myResultPairs = - // boost::make_shared< std::vector > >(); - //ResultListPointer spMyResultPairs(myResultPairs); - //boost::shared_ptr< std::vector > > myResultPairs2 = - // boost::make_shared< std::vector > >(); - //ResultListPointer spMyResultPairs2(myResultPairs2); - 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; /*to do: x should be defined based on the user's feedback*/ if (aReferenceDose <= 0) { throw core::InvalidParameterException("aReferenceDose should be >0!"); } //Vx try { sstr << aDoseStatistics->getVx(aReferenceDose * 0.02) * 1000 << columnSeparator; // cm3 to mm3 sstr << aDoseStatistics->getVx(aReferenceDose * 0.05) * 1000 << columnSeparator; // cm3 to mm3 sstr << aDoseStatistics->getVx(aReferenceDose * 0.10) * 1000 << columnSeparator; // cm3 to mm3 sstr << aDoseStatistics->getVx(aReferenceDose * 0.90) * 1000 << columnSeparator; // cm3 to mm3 sstr << aDoseStatistics->getVx(aReferenceDose * 0.95) * 1000 << columnSeparator; // cm3 to mm3 sstr << aDoseStatistics->getVx(aReferenceDose * 0.98) * 1000 << columnSeparator; // cm3 to mm3 } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } //Dx double absoluteVolume = aDoseStatistics->getVolume(); try { sstr << aDoseStatistics->getDx(absoluteVolume * 0.02) << columnSeparator; sstr << aDoseStatistics->getDx(absoluteVolume * 0.05) << columnSeparator; sstr << aDoseStatistics->getDx(absoluteVolume * 0.10) << columnSeparator; sstr << aDoseStatistics->getDx(absoluteVolume * 0.90) << columnSeparator; sstr << aDoseStatistics->getDx(absoluteVolume * 0.95) << columnSeparator; sstr << aDoseStatistics->getDx(absoluteVolume * 0.98) << columnSeparator; } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } try { //MOHx sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.02) << columnSeparator; sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.05) << columnSeparator; sstr << aDoseStatistics->getMOHx(absoluteVolume * 0.10) << columnSeparator; } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } try { //MOCx sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.02) << columnSeparator; sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.05) << columnSeparator; sstr << aDoseStatistics->getMOCx(absoluteVolume * 0.10) << columnSeparator; } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } try { //MaxOHx sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.02) << columnSeparator; sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.05) << columnSeparator; sstr << aDoseStatistics->getMaxOHx(absoluteVolume * 0.10) << columnSeparator; } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } try { //MinOCx sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.02) << columnSeparator; sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.05) << columnSeparator; sstr << aDoseStatistics->getMinOCx(absoluteVolume * 0.10) << columnSeparator; } catch (core::DataNotAvailableException) { //as data is not available (was not computed by doseStatistics), it cannot be written } return sstr.str(); } } } } \ No newline at end of file