diff --git a/code/algorithms/rttbDoseStatistics.cpp b/code/algorithms/rttbDoseStatistics.cpp index a1147a0..8d02a9b 100644 --- a/code/algorithms/rttbDoseStatistics.cpp +++ b/code/algorithms/rttbDoseStatistics.cpp @@ -1,287 +1,283 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseStatistics.h" - -#include "boost/make_shared.hpp" - #include "rttbDataNotAvailableException.h" -#include "rttbInvalidParameterException.h" namespace rttb { namespace algorithms { DoseStatistics::DoseStatistics(DoseStatisticType minimum, DoseStatisticType maximum, DoseStatisticType mean, DoseStatisticType stdDeviation, VoxelNumberType numVoxels, VolumeType volume, ResultListPointer maximumVoxelPositions /*= ResultListPointer()*/, ResultListPointer minimumVoxelPositions /*= ResultListPointer()*/, VolumeToDoseMeasureCollection Dx, DoseToVolumeMeasureCollection Vx, VolumeToDoseMeasureCollection MOHx, VolumeToDoseMeasureCollection MOCx, VolumeToDoseMeasureCollection MaxOHx, VolumeToDoseMeasureCollection MinOCx, DoseTypeGy referenceDose /*=-1*/): _minimum(minimum), _maximum(maximum), _mean(mean), _stdDeviation(stdDeviation), _numVoxels(numVoxels), _volume(volume), _Dx(::boost::make_shared(Dx)), _Vx(::boost::make_shared(Vx)), _MOHx(::boost::make_shared(MOHx)), _MOCx(::boost::make_shared(MOCx)), _MaxOHx(::boost::make_shared(MaxOHx)), _MinOCx(::boost::make_shared(MinOCx)) { if (maximumVoxelPositions == nullptr) { _maximumVoxelPositions = boost::make_shared > > (std::vector >()); } else { _maximumVoxelPositions = maximumVoxelPositions; } if (minimumVoxelPositions == nullptr) { _minimumVoxelPositions = boost::make_shared > > (std::vector >()); } else { _minimumVoxelPositions = minimumVoxelPositions; } if (referenceDose <= 0) { _referenceDose = _maximum; } else { _referenceDose = referenceDose; } } DoseStatistics::~DoseStatistics() = default; void DoseStatistics::setMinimumVoxelPositions(ResultListPointer minimumVoxelPositions) { _minimumVoxelPositions = minimumVoxelPositions; } void DoseStatistics::setMaximumVoxelPositions(ResultListPointer maximumVoxelPositions) { _maximumVoxelPositions = maximumVoxelPositions; } void DoseStatistics::setDx(VolumeToDoseMeasureCollection::Pointer DxValues) { _Dx = DxValues; } void DoseStatistics::setVx(DoseToVolumeMeasureCollection::Pointer VxValues) { _Vx = VxValues; } void DoseStatistics::setMOHx(VolumeToDoseMeasureCollection::Pointer MOHxValues) { _MOHx = MOHxValues; } void DoseStatistics::setMOCx(VolumeToDoseMeasureCollection::Pointer MOCxValues) { _MOCx = MOCxValues; } void DoseStatistics::setMaxOHx(VolumeToDoseMeasureCollection::Pointer MaxOHValues) { _MaxOHx = MaxOHValues; } void DoseStatistics::setMinOCx(VolumeToDoseMeasureCollection::Pointer MinOCValues) { _MinOCx = MinOCValues; } void DoseStatistics::setReferenceDose(DoseTypeGy referenceDose) { if (referenceDose <= 0) { _referenceDose = _maximum; } else { _referenceDose = referenceDose; } } VoxelNumberType DoseStatistics::getNumberOfVoxels() const { return _numVoxels; } VolumeType DoseStatistics::getVolume() const { return _volume; } DoseTypeGy DoseStatistics::getReferenceDose() const { return _referenceDose; } DoseStatisticType DoseStatistics::getMaximum() const { return _maximum; } DoseStatisticType DoseStatistics::getMinimum() const { return _minimum; } DoseStatisticType DoseStatistics::getMean() const { return _mean; } DoseStatisticType DoseStatistics::getStdDeviation() const { return _stdDeviation; } DoseStatisticType DoseStatistics::getVariance() const { return _stdDeviation * _stdDeviation; } double DoseStatistics::getValue(const std::map& aMap, double key, bool findNearestValueInstead, double& storedKey) const { if (aMap.find(key) != std::end(aMap)) { return aMap.find(key)->second; } else { //value not in map. We have to find the nearest value if (aMap.empty()) { throw core::DataNotAvailableException("No Vx values are defined"); } else { if (findNearestValueInstead) { auto iterator = findNearestKeyInMap(aMap, key); storedKey = iterator->first; return iterator->second; } else { throw core::DataNotAvailableException("No Vx value with required dose is defined"); } } } } std::map::const_iterator DoseStatistics::findNearestKeyInMap( const std::map& aMap, double key) const { double minDistance = 1e19; double minDistanceLast = 1e20; auto iterator = std::begin(aMap); while (iterator != std::end(aMap)) { minDistanceLast = minDistance; minDistance = fabs(iterator->first - key); if (minDistanceLast > minDistance) { ++iterator; } else { if (iterator != std::begin(aMap)) { --iterator; return iterator; } else { return std::begin(aMap); } } } --iterator; return iterator; } DoseStatistics::ResultListPointer DoseStatistics::getMaximumVoxelPositions() const { return _maximumVoxelPositions; } DoseStatistics::ResultListPointer DoseStatistics::getMinimumVoxelPositions() const { return _minimumVoxelPositions; } DoseToVolumeMeasureCollection DoseStatistics::getVx() const { return *_Vx; } VolumeToDoseMeasureCollection DoseStatistics::getDx() const { return *_Dx; } VolumeToDoseMeasureCollection DoseStatistics::getMOHx() const { return *_MOHx; } VolumeToDoseMeasureCollection DoseStatistics::getMOCx() const { return *_MOCx; } VolumeToDoseMeasureCollection DoseStatistics::getMaxOHx() const { return *_MaxOHx; } VolumeToDoseMeasureCollection DoseStatistics::getMinOCx() const { return *_MinOCx; } }//end namespace algorithms }//end namespace rttb diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp index 6969acb..85ca1c7 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,413 +1,407 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseStatisticsCalculator.h" - -#include -#include -#include -#include - -#include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" +#include "rttbNullPointerException.h" -#include +#include namespace rttb { namespace algorithms { DoseStatisticsCalculator::DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator) { if (aDoseIterator == nullptr) { throw core::NullPointerException("DoseIterator must not be nullptr"); } else { _doseIterator = aDoseIterator; } _simpleDoseStatisticsCalculated = false; _complexDoseStatisticsCalculated = false; _multiThreading = false; _mutex = ::boost::make_shared<::boost::shared_mutex>(); } DoseStatisticsCalculator::~DoseStatisticsCalculator() = default; DoseStatisticsCalculator::DoseIteratorPointer DoseStatisticsCalculator::getDoseIterator() const { return _doseIterator; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( bool computeComplexMeasures, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (_doseIterator == nullptr) { throw core::NullPointerException("_doseIterator must not be nullptr!"); } //"simple" dose statistics are mandatory calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); if (computeComplexMeasures) { //more complex dose statistics are optional with default maximum dose and default relative x values calculateComplexDoseStatistics(_statistics->getMaximum(), std::vector(), std::vector()); } return _statistics; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (_doseIterator == nullptr) { throw core::NullPointerException("_doseIterator must not be nullptr!"); } if (referenceDose <= 0) { throw rttb::core::InvalidParameterException("Reference dose must be > 0 !"); } //simple dose statistics calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); //more complex dose statistics with given reference dose and default x values calculateComplexDoseStatistics(referenceDose, std::vector(), std::vector()); return _statistics; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues, DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (_doseIterator == nullptr) { throw core::NullPointerException("_doseIterator must not be nullptr!"); } //"simple" dose statistics calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); if (referenceDose <= 0) { //more complex dose statistics with default maximum dose and relative x values calculateComplexDoseStatistics(_statistics->getMaximum(), precomputeDoseValues, precomputeVolumeValues); } else { //more complex dose statistics with given reference dose and relative x values calculateComplexDoseStatistics(referenceDose, precomputeDoseValues, precomputeVolumeValues); } return _statistics; } void DoseStatisticsCalculator::calculateSimpleDoseStatistics(unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { _doseVector.clear(); _voxelProportionVector.clear(); std::multimap doseValueVSIndexMap; std::vector voxelProportionVectorTemp; DoseStatisticType maximumDose = 0; DoseStatisticType minimumDose = std::numeric_limits::max(); DoseStatisticType meanDose; DoseStatisticType stdDeviationDose; DoseTypeGy sum = 0; VolumeType numVoxels = 0.0; DoseTypeGy squareSum = 0; VolumeType volume = 0; _doseIterator->reset(); int i = 0; DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid()) { doseValue = _doseIterator->getCurrentDoseValue(); if (i == 0) { minimumDose = doseValue; volume = _doseIterator->getCurrentVoxelVolume(); } rttb::FractionType voxelProportion = _doseIterator->getCurrentRelevantVolumeFraction(); sum += doseValue * voxelProportion; numVoxels += voxelProportion; squareSum += doseValue * doseValue * voxelProportion; if (doseValue > maximumDose) { maximumDose = doseValue; } else if (doseValue < minimumDose) { minimumDose = doseValue; } voxelProportionVectorTemp.push_back(voxelProportion); doseValueVSIndexMap.insert(std::pair(doseValue, i)); i++; _doseIterator->next(); } if (numVoxels != 0) { meanDose = sum / numVoxels; //standard deviation is defined only for n>=2 if (numVoxels >= 2) { //uncorrected variance is calculated DoseStatisticType varianceDose = (squareSum / numVoxels - meanDose * meanDose); if (varianceDose < errorConstant) { stdDeviationDose = 0; } else { stdDeviationDose = pow(varianceDose, 0.5); } } else { stdDeviationDose = 0; } } //sort dose values and corresponding volume fractions in member variables for (auto & it : doseValueVSIndexMap) { _doseVector.push_back((float)it.first); _voxelProportionVector.push_back(voxelProportionVectorTemp.at(it.second)); } volume *= numVoxels; _statistics = boost::make_shared(minimumDose, maximumDose, meanDose, stdDeviationDose, numVoxels, volume); _simpleDoseStatisticsCalculated = true; ResultListPointer minimumVoxelPositions = computeMinimumPositions(maxNumberMinimaPositions); ResultListPointer maximumVoxelPositions = computeMaximumPositions(maxNumberMaximaPositions); _statistics->setMinimumVoxelPositions(minimumVoxelPositions); _statistics->setMaximumVoxelPositions(maximumVoxelPositions); } void DoseStatisticsCalculator::calculateComplexDoseStatistics(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues) { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call calculateComplexDoseStatistics()"); } std::vector precomputeDoseValuesNonConst = precomputeDoseValues; std::vector precomputeVolumeValuesNonConst = precomputeVolumeValues; //set default values if (precomputeDoseValues.empty()) { std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)( 0.95)(0.98); precomputeDoseValuesNonConst = defaultPrecomputeDoseValues; } if (precomputeVolumeValues.empty()) { std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)( 0.95)(0.98); precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues; } _Vx = ::boost::make_shared(precomputeDoseValuesNonConst, referenceDose, _doseIterator); _Vx->compute(); _Dx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), _statistics->getMinimum()); _Dx->compute(); _MOHx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume()); _MOHx->compute(); _MOCx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume()); _MOCx->compute(); _MaxOHx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume()); _MaxOHx->compute(); _MinOCx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), _statistics->getMinimum(), _statistics->getMaximum()); _MinOCx->compute(); _statistics->setVx(_Vx->getMeasureCollection()); _statistics->setDx(_Dx->getMeasureCollection()); _statistics->setMOHx(_MOHx->getMeasureCollection()); _statistics->setMOCx(_MOCx->getMeasureCollection()); _statistics->setMaxOHx(_MaxOHx->getMeasureCollection()); _statistics->setMinOCx(_MinOCx->getMeasureCollection()); _statistics->setReferenceDose(referenceDose); _complexDoseStatisticsCalculated = true; } void DoseStatisticsCalculator::addPrecomputeValues(const std::vector& values) { if (!_complexDoseStatisticsCalculated) { throw core::InvalidDoseException("Complex DoseStatistics have to be computed in order to call addPrecomputeDoseValues()"); } _Vx->addPrecomputeDoseValues(values); _Dx->addPrecomputeVolumeValues(values); _MOHx->addPrecomputeVolumeValues(values); _MOCx->addPrecomputeVolumeValues(values); _MaxOHx->addPrecomputeVolumeValues(values); _MinOCx->addPrecomputeVolumeValues(values); } void DoseStatisticsCalculator::recalculateDoseStatistics() { if (!_complexDoseStatisticsCalculated) { throw core::InvalidDoseException("Complex DoseStatistics have to be computed in order to call recalculateDoseStatistics()"); } _Vx->compute(); _Dx->compute(); _MOHx->compute(); _MOCx->compute(); _MaxOHx->compute(); _MinOCx->compute(); } 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; } void DoseStatisticsCalculator::setMultiThreading(const bool choice) { _multiThreading = choice; } }//end namespace algorithms }//end namespace rttb diff --git a/code/algorithms/rttbDoseToVolumeMeasureCollection.cpp b/code/algorithms/rttbDoseToVolumeMeasureCollection.cpp index bc97a13..068c7d3 100644 --- a/code/algorithms/rttbDoseToVolumeMeasureCollection.cpp +++ b/code/algorithms/rttbDoseToVolumeMeasureCollection.cpp @@ -1,85 +1,84 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseToVolumeMeasureCollection.h" #include "rttbInvalidParameterException.h" -#include "rttbDataNotAvailableException.h" namespace rttb { namespace algorithms { DoseToVolumeMeasureCollection::DoseToVolumeMeasureCollection(complexStatistics name, DoseTypeGy referenceDose) : _name(name), _referenceDose(referenceDose), _values(std::map()) {} void DoseToVolumeMeasureCollection::setReferenceDose(DoseTypeGy referenceDose) { this->_referenceDose = referenceDose; } void DoseToVolumeMeasureCollection::insertValue(DoseTypeGy dose, VolumeType volume) { this->_values.insert(std::pair(dose, volume)); } VolumeType DoseToVolumeMeasureCollection::getValue(DoseTypeGy xVolumeAbsolute) const { VolumeType dummy; return getSpecificValue(_values, xVolumeAbsolute, false, dummy); } VolumeType DoseToVolumeMeasureCollection::getValue(DoseTypeGy xDoseAbsolute, bool findNearestValue, DoseTypeGy & nearestXDose) const { return getSpecificValue(_values, xDoseAbsolute, findNearestValue, nearestXDose); } VolumeType DoseToVolumeMeasureCollection::getValueRelative(DoseTypeGy xDoseRelative) const { if (_referenceDose != -1 && xDoseRelative >= 0 && xDoseRelative <= 1) { DoseTypeGy xDoseAbsolute = xDoseRelative * _referenceDose; DoseTypeGy dummy; return getSpecificValue(_values, xDoseAbsolute, false, dummy); } else { throw rttb::core::InvalidParameterException("Reference dose must be > 0 and 0 <= relative Dose <= 1"); } } VolumeType DoseToVolumeMeasureCollection::getValueRelative(DoseTypeGy xDoseRelative, bool findNearestValue, DoseTypeGy & nearestXDose) const { if (_referenceDose != -1 && xDoseRelative >= 0 && xDoseRelative <= 1) { DoseTypeGy xDoseAbsolute = xDoseRelative * _referenceDose; return getSpecificValue(_values, xDoseAbsolute, findNearestValue, nearestXDose); } else { throw rttb::core::InvalidParameterException("Reference dose must be > 0 and 0 <= relative Dose <= 1"); } } DoseToVolumeMeasureCollection::DoseToVolumeFunctionType DoseToVolumeMeasureCollection::getAllValues() const { return this->_values; } bool operator==(const DoseToVolumeMeasureCollection& volumeToDoseMesureCollection, const DoseToVolumeMeasureCollection& otherVolumeToDoseMesureCollection) { if (volumeToDoseMesureCollection.getName() == otherVolumeToDoseMesureCollection.getName() && volumeToDoseMesureCollection.getAllValues() == otherVolumeToDoseMesureCollection.getAllValues()) { return true; } return false; } } } diff --git a/code/algorithms/rttbDoseToVolumeMeasureCollectionCalculator.cpp b/code/algorithms/rttbDoseToVolumeMeasureCollectionCalculator.cpp index 419c861..d84fd39 100644 --- a/code/algorithms/rttbDoseToVolumeMeasureCollectionCalculator.cpp +++ b/code/algorithms/rttbDoseToVolumeMeasureCollectionCalculator.cpp @@ -1,86 +1,83 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseToVolumeMeasureCollectionCalculator.h" -#include #include "rttbInvalidParameterException.h" #include "rttbUtils.h" -#include - -//#include +#include namespace rttb { namespace algorithms { DoseToVolumeMeasureCollectionCalculator::DoseToVolumeMeasureCollectionCalculator(const std::vector& precomputeDoseValues, const DoseTypeGy referenceDose, const core::DoseIteratorInterface::Pointer doseIterator, DoseToVolumeMeasureCollection::complexStatistics name, bool multiThreading) : _doseIterator(doseIterator), _referenceDose(referenceDose), _measureCollection(::boost::make_shared(name)), _multiThreading(multiThreading) { addPrecomputeDoseValues(precomputeDoseValues); } void DoseToVolumeMeasureCollectionCalculator::compute() { std::vector threads; for (double _precomputeDoseValue : _precomputeDoseValues) { double xAbsolute = _precomputeDoseValue * _referenceDose; if (!rttb::core::isKey(_measureCollection->getAllValues(), xAbsolute)) { if (_multiThreading) { throw rttb::core::InvalidParameterException("MultiThreading is not implemented yet."); //threads.push_back(boost::thread(&DoseToVolumeMeasureCollectionCalculator::insertIntoMeasureCollection, this , xAbsolute, computeSpecificValue(xAbsolute))); } else { insertIntoMeasureCollection(xAbsolute, this->computeSpecificValue(xAbsolute)); } } } for (auto & thread : threads) { thread.join(); } } void DoseToVolumeMeasureCollectionCalculator::addPrecomputeDoseValues(const std::vector& values) { for (double value : values) { if (value > 1 || value < 0) { throw rttb::core::InvalidParameterException("Values must be between 1 and 0!"); } if (!rttb::core::isKey(_precomputeDoseValues, value)) { _precomputeDoseValues.push_back(value); } } } DoseToVolumeMeasureCollection::Pointer DoseToVolumeMeasureCollectionCalculator::getMeasureCollection() { return _measureCollection; } void DoseToVolumeMeasureCollectionCalculator::insertIntoMeasureCollection(DoseTypeGy xAbsolute, VolumeType resultVolume) { _measureCollection->insertValue(xAbsolute, resultVolume); } } } diff --git a/code/algorithms/rttbVolumeToDoseMeasureCollectionCalculator.cpp b/code/algorithms/rttbVolumeToDoseMeasureCollectionCalculator.cpp index abc876b..ad4ed48 100644 --- a/code/algorithms/rttbVolumeToDoseMeasureCollectionCalculator.cpp +++ b/code/algorithms/rttbVolumeToDoseMeasureCollectionCalculator.cpp @@ -1,84 +1,82 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbVolumeToDoseMeasureCollectionCalculator.h" -#include #include "rttbInvalidParameterException.h" #include "rttbUtils.h" -#include -//#include +#include namespace rttb { namespace algorithms { VolumeToDoseMeasureCollectionCalculator::VolumeToDoseMeasureCollectionCalculator(const std::vector& precomputeVolumeValues, const VolumeType volume, const std::vector& doseVector, const std::vector& voxelProportionVector, const DoseVoxelVolumeType currentVoxelVolume, VolumeToDoseMeasureCollection::complexStatistics name, bool multiThreading) : _doseVector(doseVector), _currentVoxelVolume(currentVoxelVolume), _voxelProportionVector(voxelProportionVector), _volume(volume), _measureCollection(::boost::make_shared(name)), _multiThreading(multiThreading) { addPrecomputeVolumeValues(precomputeVolumeValues); } void VolumeToDoseMeasureCollectionCalculator::compute() { std::vector threads; for (double _precomputeVolumeValue : _precomputeVolumeValues) { double xAbsolute = _precomputeVolumeValue * _volume; if (!rttb::core::isKey(_measureCollection->getAllValues(), xAbsolute)) { if (_multiThreading) { throw rttb::core::InvalidParameterException("MultiThreading is not implemented yet."); //threads.push_back(boost::thread(&VolumeToDoseMeasureCollectionCalculator::insertIntoMeasureCollection, this, xAbsolute, computeSpecificValue(xAbsolute))); } else { insertIntoMeasureCollection(xAbsolute, this->computeSpecificValue(xAbsolute)); } } } for (auto & thread : threads) { thread.join(); } } void VolumeToDoseMeasureCollectionCalculator::addPrecomputeVolumeValues(const std::vector& values) { for (double value : values) { if (value > 1 || value < 0) { throw rttb::core::InvalidParameterException("Values must be between 1 and 0!"); } if (!rttb::core::isKey(_precomputeVolumeValues, value)) { _precomputeVolumeValues.push_back(value); } } } VolumeToDoseMeasureCollection::Pointer VolumeToDoseMeasureCollectionCalculator::getMeasureCollection() { return _measureCollection; } void VolumeToDoseMeasureCollectionCalculator::insertIntoMeasureCollection(VolumeType xAbsolute, DoseTypeGy resultDose) { _measureCollection->insertValue(xAbsolute, resultDose); } } } diff --git a/code/core/rttbDVH.cpp b/code/core/rttbDVH.cpp index 880fe12..5a93e8e 100644 --- a/code/core/rttbDVH.cpp +++ b/code/core/rttbDVH.cpp @@ -1,388 +1,386 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbDVH.h" -#include "rttbException.h" #include "rttbInvalidParameterException.h" #include "rttbUtils.h" namespace rttb { namespace core { DVH::~DVH() = default; DVH::DVH(const DataDifferentialType& aDataDifferential, const DoseTypeGy& aDeltaD, const DoseVoxelVolumeType& aDeltaV, const IDType& aStructureID, const IDType& aDoseID): DVH(aDataDifferential, aDeltaD, aDeltaV, aStructureID, aDoseID, "") { } DVH::DVH(const DataDifferentialType& aDataDifferential, DoseTypeGy aDeltaD, DoseVoxelVolumeType aDeltaV, const IDType& aStructureID, const IDType& aDoseID, const IDType& aVoxelizationID): _deltaD(aDeltaD), _deltaV(aDeltaV), _structureID(aStructureID), _doseID(aDoseID), _voxelizationID(aVoxelizationID) { _dataDifferential.clear(); _dataDifferential = aDataDifferential; this->init(); } DVH::DVH(const DVH& copy) : DVH(copy._dataDifferential, copy._deltaD, copy._deltaV, copy._structureID, copy._doseID, copy._voxelizationID) { _label = copy._label; } DVH& DVH::operator=(const DVH& copy) { if (this != ©) { _deltaD = copy._deltaD; _deltaV = copy._deltaV; _structureID = copy._structureID; _doseID = copy._doseID; _voxelizationID = copy._voxelizationID; _label = copy._label; _dataDifferential.clear(); _dataDifferential = copy._dataDifferential; } this->init(); return *this; } bool operator==(const DVH& aDVH, const DVH& otherDVH) { bool result; //larger error constant because especially numberOfVoxels differ quite a bit after serialization const double errorConstantDVH = 1e-4; result = valueIsClose(aDVH.getDeltaD(), otherDVH.getDeltaD(), errorConstantDVH); result = result && valueIsClose(aDVH.getDeltaV(), otherDVH.getDeltaV(), errorConstantDVH); result = result && (aDVH.getDoseID() == otherDVH.getDoseID()); result = result && (aDVH.getStructureID() == otherDVH.getStructureID()); result = result && (aDVH.getVoxelizationID() == otherDVH.getVoxelizationID()); result = result && valueIsClose(aDVH.getNumberOfVoxels(),otherDVH.getNumberOfVoxels(), errorConstantDVH); result = result && valueIsClose(aDVH.getMaximum(), otherDVH.getMaximum(), errorConstantDVH); result = result && valueIsClose(aDVH.getMinimum(), otherDVH.getMinimum(), errorConstantDVH); result = result && valueIsClose(aDVH.getMean(), otherDVH.getMean(), errorConstantDVH); result = result && (aDVH.getDataDifferential().size() == otherDVH.getDataDifferential().size()); if (!result) { return result; } for (size_t i = 0; i < aDVH.getDataDifferential().size(); i++) { result = result && valueIsClose(aDVH.getDataDifferential().at(i), otherDVH.getDataDifferential().at(i), errorConstantDVH); } return result; } std::ostream& operator<<(std::ostream& s, const DVH& aDVH) { s << "[ " << aDVH.getStructureID() << ", " << aDVH.getDoseID() << ", " << aDVH.getVoxelizationID() << "\n " << "Number of Voxels: " << aDVH.getNumberOfVoxels() << " " << "Minimum/Maximum/Mean/Standard deviation: " << aDVH.getMinimum() << ", " << aDVH.getMaximum() << ", " << aDVH.getMean() << ", " << aDVH.getStdDeviation() << " ]"; return s; } std::deque DVH::getDataDifferential() const { return _dataDifferential; } std::deque DVH::getDataCumulative() const { return _dataCumulative; } DoseVoxelVolumeType DVH::getDeltaV() const { return _deltaV; } DoseTypeGy DVH::getDeltaD() const { return _deltaD; } IDType DVH::getDoseID() const { return this->_doseID; } IDType DVH::getStructureID() const { return this->_structureID; } IDType DVH::getVoxelizationID() const { return this->_voxelizationID; } void DVH::setDoseID(IDType aDoseID) { _doseID = aDoseID; } void DVH::setStructureID(IDType aStrID) { _structureID = aStrID; } DoseStatisticType DVH::getMaximum() const { return _maximum; } DoseStatisticType DVH::getMinimum() const { return _minimum; } DoseStatisticType DVH::getMean() const { return _mean; } DVHVoxelNumber DVH::getNumberOfVoxels() const { return _numberOfVoxels; } DoseStatisticType DVH::getStdDeviation() const { return _stdDeviation; } DoseStatisticType DVH::getVariance() const { return _variance; } void DVH::init() { if (_deltaD == 0 || _deltaV == 0) { throw InvalidParameterException("DVH init error: neither _deltaD nor _deltaV must be zero!"); } if (this->_dataDifferential.empty()) { throw InvalidParameterException("DVH init error: data differential is empty!"); } double sum = 0; double squareSum = 0; _numberOfVoxels = 0; _maximum = 0; _minimum = 0; _dataCumulative.clear(); DataDifferentialType::iterator it; int i = 0; for (it = _dataDifferential.begin(); it != _dataDifferential.end(); ++it) { _numberOfVoxels += (*it); if ((*it) > 0) { _maximum = (i + 0.5) * this->_deltaD; } if ((_minimum == 0.0f) && ((*it) > 0)) { _minimum = (i + 0.5) * this->_deltaD; } sum += (*it) * (i + 0.5) * this->_deltaD; squareSum += (*it) * pow((i + 0.5) * this->_deltaD, 2); i++; } _mean = sum / _numberOfVoxels; _variance = (squareSum / _numberOfVoxels - _mean * _mean); _stdDeviation = pow(_variance, 0.5); this->calcCumulativeDVH(); } void DVH::calcCumulativeDVH() { _dataCumulative.clear(); DoseCalcType cumulativeDVHi = 0; for (auto valueItr = _dataDifferential.rbegin(); valueItr != _dataDifferential.rend(); ++valueItr) { cumulativeDVHi += *valueItr; _dataCumulative.push_front(cumulativeDVHi); } } DoseStatisticType DVH::getMedian() const { double median_voxel = 0; int median_i = 0; for (GridIndexType i = 0; i < this->_dataDifferential.size(); i++) { if (median_voxel < (_numberOfVoxels - median_voxel)) { median_voxel += _dataDifferential[i]; median_i = i; } } double median = (median_i + 0.5) * this->_deltaD; return median; } DoseStatisticType DVH::getModal() const { double modal_voxel = 0; int modal_i = 0; for (GridIndexType i = 0; i < this->_dataDifferential.size(); i++) { if (modal_voxel < _dataDifferential[i]) { modal_voxel = _dataDifferential[i]; modal_i = i; } } double modal = (modal_i + 0.5) * this->_deltaD; return modal; } VolumeType DVH::getVx(DoseTypeGy xDoseAbsolute) const { auto i = static_cast(xDoseAbsolute / _deltaD); if (i < _dataCumulative.size()) { VolumeType vx = (_dataCumulative.at(i)); vx = (vx * this->_deltaV); return vx; } else { return 0; } } DoseTypeGy DVH::getDx(VolumeType xVolumeAbsolute) const { GridIndexType i = 0; for (; i < _dataCumulative.size(); i++) { double volumeAbsoluteI = _dataCumulative[i] * this->_deltaV; if (xVolumeAbsolute > volumeAbsoluteI) { break; } } if (i <= _dataCumulative.size() && i > 0) { DoseTypeGy dx = (i - 1) * this->_deltaD; return dx; } else { return 0; } } VolumeType DVH::getAbsoluteVolume(int relativePercent) const { return (relativePercent * getNumberOfVoxels() * getDeltaV() / 100.0); } std::deque DVH::convertAbsoluteToRelative(bool isCumulative) const { DataDifferentialType relativeData, absoluteData; if (isCumulative) { absoluteData = getDataCumulative(); } else { absoluteData = getDataDifferential(); } for (std::deque::iterator it = absoluteData.begin(); it != absoluteData.end(); ++it) { relativeData.push_back((*it) / getNumberOfVoxels()); } return relativeData; } void DVH::setLabel(StructureLabel aLabel) { _label = aLabel; } StructureLabel DVH::getLabel() const { return _label; } std::map DVH::getNormalizedDVH(DVHType dvhType) const { std::map normalizedDVH; DataDifferentialType data; if (dvhType.Type == DVHType::Cumulative) { data = getDataCumulative(); } else { data = getDataDifferential(); } if (data.empty()) { throw InvalidParameterException("DVH data is empty. Can't retrieve normalized DVH"); } for (size_t i = 0; i < data.size(); i++) { normalizedDVH.insert(std::pair(i * getDeltaD(), data[i] * getDeltaV())); } return normalizedDVH; } }//end namespace core }//end namespace rttb diff --git a/code/core/rttbGeometricInfo.cpp b/code/core/rttbGeometricInfo.cpp index 5d78d25..0d1af7f 100644 --- a/code/core/rttbGeometricInfo.cpp +++ b/code/core/rttbGeometricInfo.cpp @@ -1,304 +1,302 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbGeometricInfo.h" -#include #include -#include namespace rttb { namespace core { void GeometricInfo::setSpacing(const SpacingVectorType3D& aSpacingVector) { _spacing = aSpacingVector; } const SpacingVectorType3D& GeometricInfo::getSpacing() const { return _spacing; } void GeometricInfo::setImagePositionPatient(const WorldCoordinate3D& aImagePositionPatient) { _imagePositionPatient = aImagePositionPatient; } const WorldCoordinate3D& GeometricInfo::getImagePositionPatient() const { return _imagePositionPatient; } void GeometricInfo::setOrientationMatrix(const OrientationMatrix& anOrientationMatrix) { _orientationMatrix = anOrientationMatrix; computeInvertOrientation(); } bool GeometricInfo::computeInvertOrientation() { using pmatrix = boost::numeric::ublas::permutation_matrix; boost::numeric::ublas::matrix A(_orientationMatrix); // create a permutation matrix for the LU-factorization pmatrix pm(A.size1()); size_t res = boost::numeric::ublas::lu_factorize(A, pm); if (res != 0) { return false; } _invertedOrientationMatrix.assign(boost::numeric::ublas::identity_matrix (A.size1())); // backsubstitute to get the inverse boost::numeric::ublas::lu_substitute(A, pm, _invertedOrientationMatrix); return true; } void GeometricInfo::setImageSize(const ImageSize& aSize) { setNumColumns(aSize(0)); setNumRows(aSize(1)); setNumSlices(aSize(2)); } const ImageSize GeometricInfo::getImageSize() const { return ImageSize(static_cast(getNumColumns()), static_cast(getNumRows()), static_cast(getNumSlices())); } void GeometricInfo::setNumColumns(const VoxelGridDimensionType aValue) { _numberOfColumns = aValue; } const VoxelGridDimensionType GeometricInfo::getNumColumns() const { return _numberOfColumns; } void GeometricInfo::setNumRows(const VoxelGridDimensionType aValue) { _numberOfRows = aValue; } const VoxelGridDimensionType GeometricInfo::getNumRows() const { return _numberOfRows; } void GeometricInfo::setNumSlices(const VoxelGridDimensionType aValue) { _numberOfFrames = aValue; } const VoxelGridDimensionType GeometricInfo::getNumSlices() const { return _numberOfFrames; } bool operator==(const GeometricInfo& gInfo, const GeometricInfo& gInfo1) { return (gInfo.getImagePositionPatient() == gInfo1.getImagePositionPatient() && gInfo.getOrientationMatrix() == gInfo1.getOrientationMatrix() && gInfo.getSpacing() == gInfo1.getSpacing() && gInfo.getNumColumns() == gInfo1.getNumColumns() && gInfo.getNumRows() == gInfo1.getNumRows() && gInfo.getNumSlices() == gInfo1.getNumSlices()); } bool GeometricInfo::equalsAlmost(const GeometricInfo& another, double errorConstantGI /*= 1e-5*/) const { return (getImagePositionPatient().equalsAlmost(another.getImagePositionPatient(), errorConstantGI) && getOrientationMatrix().equalsAlmost(another.getOrientationMatrix(), errorConstantGI) && getSpacing().equalsAlmost(another.getSpacing(), errorConstantGI) && getNumColumns() == another.getNumColumns() && getNumRows() == another.getNumRows() && getNumSlices() == another.getNumSlices()); } bool GeometricInfo::worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, DoubleVoxelGridIndex3D& aIndex) const { WorldCoordinate3D distanceToIP; distanceToIP = aWorldCoordinate - _imagePositionPatient; boost::numeric::ublas::vector result = boost::numeric::ublas::prod( _invertedOrientationMatrix, distanceToIP); boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_div(result, _spacing); aIndex = DoubleVoxelGridIndex3D(resultS(0), resultS(1), resultS(2)); //if we convert DoubleVoxelGridIndex3D (double) to VoxelGridIndex3D (unsigned int), we can't find out if it's negative. //So we have to check before. if (aIndex(0) < -0.5 || aIndex(1) < -0.5 || aIndex(2) < -0.5){ return false; } else { //check if it is inside VoxelGridIndex3D indexInt = VoxelGridIndex3D(GridIndexType(aIndex(0) + 0.5), GridIndexType(aIndex(1) + 0.5), GridIndexType(aIndex(2) + 0.5)); return isInside(indexInt); } } bool GeometricInfo::worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, VoxelGridIndex3D& aIndex) const { DoubleVoxelGridIndex3D doubleIndex; bool inside = worldCoordinateToGeometryCoordinate(aWorldCoordinate, doubleIndex); aIndex = VoxelGridIndex3D(GridIndexType(doubleIndex(0)+0.5), GridIndexType(doubleIndex(1)+0.5), GridIndexType(doubleIndex(2)+0.5)); return inside; } bool GeometricInfo::geometryCoordinateToWorldCoordinate(const DoubleVoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const { boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_prod( aIndex, _spacing); boost::numeric::ublas::vector result = boost::numeric::ublas::prod( _orientationMatrix, resultS); aWorldCoordinate = result + _imagePositionPatient; //if we convert DoubleVoxelGridIndex3D (double) to VoxelGridIndex3D (unsigned int), we can't find out if it's negative. //So we have to check before. if (aIndex(0) < -0.5 || aIndex(1) < -0.5 || aIndex(2) < -0.5){ return false; } else { VoxelGridIndex3D indexInt = VoxelGridIndex3D(GridIndexType(aIndex(0) + 0.5), GridIndexType(aIndex(1) + 0.5), GridIndexType(aIndex(2) + 0.5)); return isInside(indexInt); } } bool GeometricInfo::indexToWorldCoordinate(const VoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const { DoubleVoxelGridIndex3D indexDouble = DoubleVoxelGridIndex3D(aIndex(0), aIndex(1), aIndex(2)); return geometryCoordinateToWorldCoordinate(indexDouble, aWorldCoordinate); } bool GeometricInfo::isInside(const VoxelGridIndex3D& aIndex) const { return (aIndex(0) >= 0 && aIndex(1) >= 0 && aIndex(2) >= 0 && aIndex(0) < static_cast(_numberOfColumns) && aIndex(1) < static_cast(_numberOfRows) && aIndex(2) < static_cast(_numberOfFrames)); } bool GeometricInfo::isInside(const WorldCoordinate3D& aWorldCoordinate) const { VoxelGridIndex3D currentIndex; return (worldCoordinateToIndex(aWorldCoordinate, currentIndex)); } const GridSizeType GeometricInfo::getNumberOfVoxels() const { auto nVoxels = static_cast(_numberOfRows * _numberOfColumns * _numberOfFrames); return nVoxels; } bool GeometricInfo::convert(const VoxelGridID& gridID, VoxelGridIndex3D& gridIndex) const { if (validID(gridID)) { gridIndex(0) = gridID % getNumColumns(); VoxelGridID tempID = (gridID - gridIndex.x()) / getNumColumns(); gridIndex(1) = tempID % getNumRows(); gridIndex(2) = (tempID - gridIndex.y()) / getNumRows(); return true; } return false; } bool GeometricInfo::convert(const VoxelGridIndex3D& gridIndex, VoxelGridID& gridID) const { if ((gridIndex.x() >= static_cast(getNumColumns())) || (gridIndex.y() >= static_cast(getNumRows())) || (gridIndex.z() >= static_cast(getNumSlices()))) { return false; } else { gridID = gridIndex.z() * getNumColumns() * getNumRows() + gridIndex.y() * getNumColumns() + gridIndex.x(); return validID(gridID); } } bool GeometricInfo::validID(const VoxelGridID aID) const { return (aID >= 0 && aID < getNumberOfVoxels()); } bool GeometricInfo::validIndex(const VoxelGridIndex3D& aIndex) const { VoxelGridID aID; if (!convert(aIndex, aID)) { return false; } else { return validID(aID); } } std::ostream& operator<<(std::ostream& s, const GeometricInfo& anGeometricInfo) { s << "[ " << anGeometricInfo.getImagePositionPatient() << "; " << anGeometricInfo.getOrientationMatrix() << "; " << anGeometricInfo.getSpacing() << "; " << "; " << anGeometricInfo.getNumColumns() << "; " << anGeometricInfo.getNumRows() << "; " << anGeometricInfo.getNumSlices() << " ]"; return s; } }//end namespace core }//end namespace rttb diff --git a/code/core/rttbMaskedDoseIteratorInterface.cpp b/code/core/rttbMaskedDoseIteratorInterface.cpp index cc89049..0195659 100644 --- a/code/core/rttbMaskedDoseIteratorInterface.cpp +++ b/code/core/rttbMaskedDoseIteratorInterface.cpp @@ -1,43 +1,41 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbMaskedDoseIteratorInterface.h" -#include "rttbMaskAccessorInterface.h" #include "rttbNullPointerException.h" -#include "rttbException.h" namespace rttb { namespace core { MaskedDoseIteratorInterface::MaskedDoseIteratorInterface(MaskAccessorPointer aMaskAccessor, DoseAccessorPointer aDoseAccessor): DoseIteratorInterface(aDoseAccessor) { if (aMaskAccessor == nullptr) { throw NullPointerException(" mask pointer must not be nullptr!"); } else { _spMask = aMaskAccessor; } if (!(_spMask->getGeometricInfo() == _spDoseAccessor->getGeometricInfo())) { throw Exception("Mask and Dose need to be defined on the same grid"); } } } } \ No newline at end of file diff --git a/code/core/rttbStrVectorStructureSetGenerator.cpp b/code/core/rttbStrVectorStructureSetGenerator.cpp index e0651f7..59cd183 100644 --- a/code/core/rttbStrVectorStructureSetGenerator.cpp +++ b/code/core/rttbStrVectorStructureSetGenerator.cpp @@ -1,61 +1,61 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbStrVectorStructureSetGenerator.h" -#include #include "boost/make_shared.hpp" +#include namespace rttb { namespace core { StrVectorStructureSetGenerator::StrVectorStructureSetGenerator(std::vector& aStructureVector, IDType aPatientUID) { _patientUID = aPatientUID; _strVector = aStructureVector; } StrVectorStructureSetGenerator::StructureSetPointer StrVectorStructureSetGenerator::generateStructureSet() { std::vector _filteredStructs = _strVector; if (this->getStructureLabelFilterActive()) { _filteredStructs.clear(); std::regex e(this->getFilterRegEx()); std::vector::iterator it; for(it= _strVector.begin();it!=_strVector.end();++it) { if (std::regex_match((*it)->getLabel(), e)) { _filteredStructs.push_back((*it)); } } } return boost::make_shared(_filteredStructs, _patientUID); } } }//end namespace rttb diff --git a/code/core/rttbStructure.cpp b/code/core/rttbStructure.cpp index 22dfe89..a0bc0d7 100644 --- a/code/core/rttbStructure.cpp +++ b/code/core/rttbStructure.cpp @@ -1,139 +1,134 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include -#include -#include #include #include #include "rttbStructure.h" -#include "rttbNullPointerException.h" -#include "rttbInvalidParameterException.h" namespace rttb { namespace core { /*! Compares two polygons in the same plane. Helper function for sorting of polygons. */ bool comparePolygon(PolygonType A, PolygonType B) { PolygonType::iterator it; for (it = A.begin(); it != A.end(); ++it) { if ((*it)(2) != A.at(0)(2)) { throw std::range_error("Error: A must in the same _z plane!"); } } PolygonType::iterator it2; for (it2 = B.begin(); it2 != B.end(); ++it2) { if ((*it2)(2) != B.at(0)(2)) { throw std::range_error("Error: B must in the same _z plane!"); } } if (A.size() == 0 || B.size() == 0) { throw std::range_error("Error: A and B must not be empty!"); } return (A.at(0)(2) < B.at(0)(2)); } Structure::Structure() : _structureVector(0), _label("") { if (_strUID.empty()) { boost::uuids::uuid id; boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _strUID = ss.str(); } } Structure::Structure(PolygonSequenceType strVector) : Structure() { _structureVector = strVector; sort(_structureVector.begin(), _structureVector.end(), comparePolygon); } Structure::Structure(const Structure& copy) : _structureVector(copy.getStructureVector()), _strUID(copy.getUID()), _label(copy.getLabel()) { } Structure::~Structure() = default; const PolygonSequenceType& Structure::getStructureVector() const { return _structureVector; } int Structure::getNumberOfEndpoints() const { int count = 0; PolygonSequenceType::const_iterator itVV; for (itVV = _structureVector.begin(); itVV != _structureVector.end(); ++itVV) { count += (int)(*itVV).size(); } return count; } IDType Structure::getUID() const { return _strUID; } void Structure::setUID(const IDType& aUID) { _strUID = aUID; } void Structure::setLabel(const StructureLabel& aLabel) { _label = aLabel; } StructureLabel Structure::getLabel() const { return _label; } }//end namespace core }//end namespace rttb diff --git a/code/core/rttbStructureSet.cpp b/code/core/rttbStructureSet.cpp index 10366cd..565e06c 100644 --- a/code/core/rttbStructureSet.cpp +++ b/code/core/rttbStructureSet.cpp @@ -1,78 +1,76 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include -#include #include #include #include "rttbStructureSet.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace core { StructureSet::StructureSet(const std::vector& aStructureVector, IDType aPatientUID, IDType aUID) { _structureSetVector = aStructureVector; _patientUID = aPatientUID; _UID = aUID; if (_UID == "") { boost::uuids::uuid id; boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _UID = ss.str(); } } Structure::Pointer StructureSet::getStructure(size_t aStructureNo) const { auto size = this->getNumberOfStructures(); if (aStructureNo >= size) { std::stringstream sstr; sstr << "aStructureNo must be between 0 and " << size; throw InvalidParameterException(sstr.str()); } return _structureSetVector.at(aStructureNo); } StructureSet::NumberOfStructuresType StructureSet::getNumberOfStructures() const { return _structureSetVector.size(); } IDType StructureSet::getUID() const { return _UID; } IDType StructureSet::getPatientUID() const { return _patientUID; } }//end namespace core }//end namespace rttb diff --git a/code/indices/rttbConformalIndex.cpp b/code/indices/rttbConformalIndex.cpp index a474ef2..464a5ef 100644 --- a/code/indices/rttbConformalIndex.cpp +++ b/code/indices/rttbConformalIndex.cpp @@ -1,120 +1,119 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbConformalIndex.h" -#include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace indices { ConformalIndex::ConformalIndex(core::DVHSet::Pointer dvhSet, DoseTypeGy aDoseReference) : DvhBasedDoseIndex(dvhSet, aDoseReference) { init(); } bool ConformalIndex::calcIndex() { VolumeType TV = _dvhSet->getTargetVolume(0); VolumeType Vref = _dvhSet->getWholeVolume(_doseReference); if (TV != 0 && Vref != 0) { _value = (_dvhSet->getTargetVolume(_doseReference) / TV) * (_dvhSet->getTargetVolume(_doseReference) / Vref); std::vector dvhHTSet = this->_dvhSet->getHealthyTissueSet(); std::vector::iterator it; for (it = dvhHTSet.begin(); it != dvhHTSet.end(); ++it) { core::DVH dvh = *(it); VolumeType HT = dvh.getVx(0); if (HT != 0) { _value *= (1 - dvh.getVx(this->_doseReference) / HT); } } } else if (TV == 0) { throw core::InvalidParameterException("DVH Set invalid: Target volume should not be 0!"); } else { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } return true; } IndexValueType ConformalIndex::getValueAt(core::DVHSet::IndexType tvIndex) { std::vector dvhTVSet = this->_dvhSet->getTargetVolumeSet(); VolumeType Vref = _dvhSet->getWholeVolume(_doseReference); if (tvIndex >= dvhTVSet.size()) { rttbExceptionMacro(core::InvalidParameterException, << "tvIndex invalid: it should be <" << dvhTVSet.size() << "!"); } else { core::DVH dvh = dvhTVSet.at(tvIndex); VolumeType TV = dvh.getVx(0); if (TV == 0) { throw core::InvalidParameterException("DVH invalid: Volume of tvIndex should not be 0!"); } else if (Vref == 0) { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } double value = dvh.getVx(_doseReference) / TV; //the irradiation factor of i-th target volume value = value * dvh.getVx(_doseReference) / Vref; //conformation number std::vector dvhHTSet = this->_dvhSet->getHealthyTissueSet(); std::vector::iterator it; for (it = dvhHTSet.begin(); it != dvhHTSet.end(); ++it) { dvh = *(it); VolumeType HT = dvh.getVx(0); if (HT != 0) { value *= (1 - dvh.getVx(this->_doseReference) / HT); } } return value; } } }//end namespace indices }//end namespace rttb diff --git a/code/indices/rttbConformationNumber.cpp b/code/indices/rttbConformationNumber.cpp index d034a23..340a019 100644 --- a/code/indices/rttbConformationNumber.cpp +++ b/code/indices/rttbConformationNumber.cpp @@ -1,91 +1,90 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbConformationNumber.h" -#include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace indices { ConformationNumber::ConformationNumber(core::DVHSet::Pointer dvhSet, DoseTypeGy aDoseReference) : DvhBasedDoseIndex(dvhSet, aDoseReference) { init(); } bool ConformationNumber::calcIndex() { VolumeType TV = _dvhSet->getTargetVolume(0); VolumeType Vref = _dvhSet->getWholeVolume(_doseReference); if (TV != 0 && Vref != 0) { _value = (_dvhSet->getTargetVolume(_doseReference) / TV) * (_dvhSet->getTargetVolume(_doseReference) / Vref); } else if (TV == 0) { throw core::InvalidParameterException("DVH Set invalid: Target volume should not be 0!"); } else { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } return true; } IndexValueType ConformationNumber::getValueAt(core::DVHSet::IndexType tvIndex) { std::vector dvhTVSet = this->_dvhSet->getTargetVolumeSet(); VolumeType Vref = _dvhSet->getWholeVolume(_doseReference); if (tvIndex >= dvhTVSet.size()) { rttbExceptionMacro(core::InvalidParameterException, << "tvIndex invalid: it should be <" << dvhTVSet.size() << "!"); } else { core::DVH dvh = dvhTVSet.at(tvIndex); VolumeType TV = dvh.getVx(0); if (TV == 0) { throw core::InvalidParameterException("DVH invalid: Volume of tvIndex should not be 0!"); } else if (Vref == 0) { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } IndexValueType value = dvh.getVx(_doseReference) / TV; //the irradiation factor of i-th target volume value = value * dvh.getVx(_doseReference) / Vref; return value; } } }//end namespace indices }//end namespace rttb \ No newline at end of file diff --git a/code/indices/rttbConformityIndex.cpp b/code/indices/rttbConformityIndex.cpp index 5353a18..c1dc102 100644 --- a/code/indices/rttbConformityIndex.cpp +++ b/code/indices/rttbConformityIndex.cpp @@ -1,90 +1,89 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbConformityIndex.h" -#include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace indices { ConformityIndex::ConformityIndex(core::DVHSet::Pointer dvhSet, DoseTypeGy aDoseReference) : DvhBasedDoseIndex(dvhSet, aDoseReference) { init(); } bool ConformityIndex::calcIndex() { VolumeType TV = _dvhSet->getTargetVolume(0); VolumeType Vref = _dvhSet->getWholeVolume(_doseReference); if (TV != 0 && Vref != 0) { _value = (_dvhSet->getTargetVolume(this->_doseReference) / TV) * (1 - _dvhSet->getHealthyTissueVolume(_doseReference) / Vref); } else if (TV == 0) { throw core::InvalidParameterException("DVH Set invalid: Target volume should not be 0!"); } else { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } return true; } IndexValueType ConformityIndex::getValueAt(core::DVHSet::IndexType tvIndex) { std::vector dvhTVSet = this->_dvhSet->getTargetVolumeSet(); VolumeType Vref = _dvhSet->getWholeVolume(_doseReference); if (tvIndex >= dvhTVSet.size()) { rttbExceptionMacro(core::InvalidParameterException, << "tvIndex invalid: it should be <" << dvhTVSet.size() << "!"); } else { core::DVH dvh = dvhTVSet.at(tvIndex); VolumeType TV = dvh.getVx(0); if (TV == 0) { throw core::InvalidParameterException("DVH invalid: Volume of tvIndex should not be 0!"); } else if (Vref == 0) { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } double value = dvh.getVx(_doseReference) / TV; //the irradiation factor of i-th treated volume value = value * (1 - _dvhSet->getHealthyTissueVolume(_doseReference) / Vref); return value; } } }//end namespace indices }//end namespace rttb diff --git a/code/indices/rttbCoverageIndex.cpp b/code/indices/rttbCoverageIndex.cpp index 93bec42..23d0029 100644 --- a/code/indices/rttbCoverageIndex.cpp +++ b/code/indices/rttbCoverageIndex.cpp @@ -1,73 +1,72 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbCoverageIndex.h" -#include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace indices { CoverageIndex::CoverageIndex(core::DVHSet::Pointer dvhSet, DoseTypeGy aDoseReference) : DvhBasedDoseIndex(dvhSet, aDoseReference) { init(); } bool CoverageIndex::calcIndex() { VolumeType TV = _dvhSet->getTargetVolume(0); if (TV != 0) { _value = _dvhSet->getTargetVolume(this->_doseReference) / TV; } else { throw core::InvalidParameterException("DVH Set invalid: Target volume should not be 0!"); } return true; } IndexValueType CoverageIndex::getValueAt(core::DVHSet::IndexType tvIndex) { std::vector dvhTVSet = this->_dvhSet->getTargetVolumeSet(); if (tvIndex >= dvhTVSet.size()) { rttbExceptionMacro(core::InvalidParameterException, << "tvIndex invalid: it should be <" << dvhTVSet.size() << "!"); } core::DVH dvh = dvhTVSet.at(tvIndex); VolumeType TV = dvh.getVx(0); if (TV == 0) { throw core::InvalidParameterException("DVH invalid: Volume of tvIndex should not be 0!"); } IndexValueType value = dvh.getVx(_doseReference) / TV; //the irradiation factor of i-th treated volume return value; } }//end namespace indices }//end namespace rttb diff --git a/code/indices/rttbDoseIndex.cpp b/code/indices/rttbDoseIndex.cpp index 309e8da..47117c7 100644 --- a/code/indices/rttbDoseIndex.cpp +++ b/code/indices/rttbDoseIndex.cpp @@ -1,75 +1,74 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseIndex.h" -#include "rttbException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace indices { DoseIndex::DoseIndex(DoseTypeGy aDoseReference) : _doseReference(aDoseReference), _initSuccess(false) {} DoseIndex::~DoseIndex()= default; bool DoseIndex::init() { if (!(this->checkInputs())) { throw core::InvalidParameterException("Check inputs failed: invalid parameters! "); } if (this->calcIndex()) { _initSuccess = true; } else { throw core::InvalidParameterException("Index calculation failed! "); } return _initSuccess; } void DoseIndex::setDoseReference(DoseTypeGy aDoseReference) { _doseReference = aDoseReference; _initSuccess = false; init(); } DoseTypeGy DoseIndex::getDoseReference() const { return _doseReference; } IndexValueType DoseIndex::getValue() const { if (_initSuccess) { return _value; } else { throw core::Exception("DoseIndex init error: init() must be called first!"); } } } } diff --git a/code/interpolation/rttbInterpolationBase.cpp b/code/interpolation/rttbInterpolationBase.cpp index 282bb3f..6fd316f 100644 --- a/code/interpolation/rttbInterpolationBase.cpp +++ b/code/interpolation/rttbInterpolationBase.cpp @@ -1,191 +1,190 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbInterpolationBase.h" #include "rttbInvalidParameterException.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace interpolation { void InterpolationBase::setAccessorPointer(core::AccessorInterface::ConstPointer originalData) { if (originalData != nullptr) { _spOriginalData = originalData; } else { throw core::NullPointerException("originalDose is nullptr!"); } }; void InterpolationBase::getNeighborhoodVoxelValues( const WorldCoordinate3D& aWorldCoordinate, unsigned int neighborhood, std::array& target, boost::shared_ptr values) const { if (_spOriginalData == nullptr) { throw core::NullPointerException("originalDose is nullptr!"); } //Determine target (abs(desired worldCoordinate- corner pixel world coordinate/pixel spacing) and values of corner pixels (from originalDose) VoxelGridIndex3D aIndex; if (_spOriginalData->getGeometricInfo().worldCoordinateToIndex(aWorldCoordinate, aIndex)) { //determine just the nearest voxel to the world coordinate if (neighborhood == 0) { values[0] = _spOriginalData->getValueAt(aIndex); } //determine the 8 voxels around the world coordinate else if (neighborhood == 8) { std::list cornerPoints; WorldCoordinate3D theNextVoxel; _spOriginalData->getGeometricInfo().indexToWorldCoordinate(aIndex, theNextVoxel); SpacingVectorType3D pixelSpacing = (_spOriginalData->getGeometricInfo()).getSpacing(); VoxelGridIndex3D leftTopFrontCoordinate; //find the voxel with the smallest coordinate values in each dimension. This defines the standard cube for (unsigned int i = 0; i < 3; i++) { if (aWorldCoordinate[i] < theNextVoxel[i]) { if (aIndex[i] > 0) { leftTopFrontCoordinate[i] = aIndex[i] - 1; target[i] = (aWorldCoordinate[i] - (theNextVoxel[i] - pixelSpacing[i])) / pixelSpacing[i]; } //@todo: see T22315 else { leftTopFrontCoordinate[i] = aIndex[i]; target[i] = (aWorldCoordinate[i] - (theNextVoxel[i] - pixelSpacing[i])) / pixelSpacing[i]; } } else { leftTopFrontCoordinate[i] = aIndex[i]; target[i] = (aWorldCoordinate[i] - theNextVoxel[i]) / pixelSpacing[i]; } } for (unsigned int zIncr = 0; zIncr < 2; zIncr++) { for (unsigned int yIncr = 0; yIncr < 2; yIncr++) { for (unsigned int xIncr = 0; xIncr < 2; xIncr++) { cornerPoints.emplace_back(leftTopFrontCoordinate[0] + xIncr, leftTopFrontCoordinate[1] + yIncr, leftTopFrontCoordinate[2] + zIncr); } } } //target range has to be always [0,1] for (unsigned int i = 0; i < 3; i++) { assert(target[i] >= 0.0 && target[i] <= 1.0); } unsigned int count = 0; //now just get the values of all (dose) voxels and store them in values for (auto cornerPointsIterator = cornerPoints.begin(); cornerPointsIterator != cornerPoints.end(); ++cornerPointsIterator, ++count) { if (_spOriginalData->getGeometricInfo().isInside(*cornerPointsIterator)) { values[count] = _spOriginalData->getValueAt(*cornerPointsIterator); } else { //outside value! boundary treatment values[count] = getNearestInsideVoxelValue(*cornerPointsIterator); } assert(values[count] != -1); } } else { throw core::InvalidParameterException("neighborhoods other than 0 and 8 not yet supported in Interpolation"); } } else { throw core::MappingOutsideOfImageException("Error in conversion from world coordinates to index"); } } DoseTypeGy InterpolationBase::getNearestInsideVoxelValue(const VoxelGridIndex3D& currentVoxelIndex) const { VoxelGridIndex3D voxelChangedXYZ[] = {currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex}; unsigned int runningIndex; //x,y,z for (runningIndex = 0; runningIndex < 3; ++runningIndex) { voxelChangedXYZ[runningIndex][runningIndex] -= 1; } //xy voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][1] -= 1; ++runningIndex; //xz voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; //yz voxelChangedXYZ[runningIndex][1] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; //xyz voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][1] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; unsigned int replacementVoxelIndex = 0; while (replacementVoxelIndex < runningIndex) { if (_spOriginalData->getGeometricInfo().validIndex(voxelChangedXYZ[replacementVoxelIndex])) { return _spOriginalData->getValueAt(voxelChangedXYZ[replacementVoxelIndex]); } ++replacementVoxelIndex; } return -1; } }//end namespace core }//end namespace rttb diff --git a/code/interpolation/rttbNearestNeighborInterpolation.cpp b/code/interpolation/rttbNearestNeighborInterpolation.cpp index 6493de5..871f172 100644 --- a/code/interpolation/rttbNearestNeighborInterpolation.cpp +++ b/code/interpolation/rttbNearestNeighborInterpolation.cpp @@ -1,35 +1,34 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbNearestNeighborInterpolation.h" -#include #include namespace rttb { namespace interpolation { DoseTypeGy NearestNeighborInterpolation::getValue(const WorldCoordinate3D& aWorldCoordinate) const { //proper initialization of target and values (although target is irrelevant in nearest neighbor case) std::array target = {{0.0, 0.0, 0.0}}; auto values = boost::make_shared(8); getNeighborhoodVoxelValues(aWorldCoordinate, 0, target, values); return values[0]; } } } diff --git a/code/interpolation/rttbRosuMappableDoseAccessor.cpp b/code/interpolation/rttbRosuMappableDoseAccessor.cpp index 5e75151..aca020a 100644 --- a/code/interpolation/rttbRosuMappableDoseAccessor.cpp +++ b/code/interpolation/rttbRosuMappableDoseAccessor.cpp @@ -1,167 +1,165 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ -#include "rttbRosuMappableDoseAccessor.h" - #include -#include "rttbNullPointerException.h" +#include "rttbRosuMappableDoseAccessor.h" #include "rttbMappingOutsideOfImageException.h" #include "rttbLinearInterpolation.h" namespace rttb { namespace interpolation { RosuMappableDoseAccessor::RosuMappableDoseAccessor(const core::GeometricInfo& geoInfoTargetImage, core::DoseAccessorInterface::ConstPointer doseMovingImage, const TransformationInterface::Pointer aTransformation, bool acceptPadding, DoseTypeGy defaultOutsideValue): MappableDoseAccessorInterface(geoInfoTargetImage, doseMovingImage, aTransformation, acceptPadding, defaultOutsideValue) { //define linear interpolation auto interpolationLinear = ::boost::make_shared(); _spInterpolation = interpolationLinear; _spInterpolation->setAccessorPointer(_spOriginalDoseDataMovingImage); } GenericValueType RosuMappableDoseAccessor::getValueAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex3D; if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) { return getValueAt(aVoxelGridIndex3D); } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); } } } GenericValueType RosuMappableDoseAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { //Transform requested voxel coordinates of original image into world coordinates with RTTB WorldCoordinate3D worldCoordinateTarget; if (_geoInfoTargetImage.indexToWorldCoordinate(aIndex, worldCoordinateTarget)) { std::vector octants = getOctants(worldCoordinateTarget); if (octants.size() > 2) { DoseTypeGy interpolatedDoseValue = 0.0; //get trilinear interpolation value of every octant point for (std::vector::const_iterator octantIterator = octants.begin(); octantIterator != octants.end(); ++octantIterator) { //transform coordinates WorldCoordinate3D worldCoordinateMoving; WorldCoordinate3D worldCoordinateTargetOctant = *octantIterator; _spTransformation->transformInverse(worldCoordinateTargetOctant, worldCoordinateMoving); try { interpolatedDoseValue += _spInterpolation->getValue(worldCoordinateMoving); } //Mapped outside of image? Check if padding is allowed catch (core::MappingOutsideOfImageException& /*e*/) { if (_acceptPadding) { interpolatedDoseValue += _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Mapping outside of image"); } } catch (core::Exception& e) { std::cout << e.what() << std::endl; return -1; } } return interpolatedDoseValue / (DoseTypeGy)octants.size(); } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Too many samples are mapped outside the image!"); } } } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); } } } std::vector RosuMappableDoseAccessor::getOctants( const WorldCoordinate3D& aCoordinate) const { std::vector octants; SpacingVectorType3D spacingTargetImage = _geoInfoTargetImage.getSpacing(); core::GeometricInfo geometricInfoDoseData = _spOriginalDoseDataMovingImage->getGeometricInfo(); //as the corner point is the coordinate of the voxel (grid), 0.25 and 0.75 are the center of the subvoxels for (double xOct = -0.25; xOct <= 0.25; xOct += 0.5) { for (double yOct = -0.25; yOct <= 0.25; yOct += 0.5) { for (double zOct = -0.25; zOct <= 0.25; zOct += 0.5) { WorldCoordinate3D aWorldCoordinate(aCoordinate.x() + (xOct * spacingTargetImage.x()), aCoordinate.y() + (yOct * spacingTargetImage.y()), aCoordinate.z() + (zOct * spacingTargetImage.z())); if (geometricInfoDoseData.isInside(aWorldCoordinate)) { octants.push_back(aWorldCoordinate); } } } } return octants; } }//end namespace interpolation }//end namespace rttb diff --git a/code/interpolation/rttbSimpleMappableDoseAccessor.cpp b/code/interpolation/rttbSimpleMappableDoseAccessor.cpp index c0c3a76..273bf8f 100644 --- a/code/interpolation/rttbSimpleMappableDoseAccessor.cpp +++ b/code/interpolation/rttbSimpleMappableDoseAccessor.cpp @@ -1,113 +1,112 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbSimpleMappableDoseAccessor.h" -#include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace interpolation { SimpleMappableDoseAccessor::SimpleMappableDoseAccessor(const core::GeometricInfo& geoInfoTargetImage, core::DoseAccessorInterface::ConstPointer doseMovingImage, const TransformationInterface::Pointer aTransformation, const InterpolationBase::Pointer aInterpolation, bool acceptPadding, double defaultOutsideValue): MappableDoseAccessorInterface(geoInfoTargetImage, doseMovingImage, aTransformation, acceptPadding, defaultOutsideValue), _spInterpolation(aInterpolation) { //handle null pointers if (aInterpolation == nullptr) { throw core::NullPointerException("Pointer to interpolation cannot be nullptr."); } else { _spInterpolation->setAccessorPointer(_spOriginalDoseDataMovingImage); } } GenericValueType SimpleMappableDoseAccessor::getValueAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex3D; if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) { return getValueAt(aVoxelGridIndex3D); } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); } } } GenericValueType SimpleMappableDoseAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { //Transform requested voxel coordinates of original image into world coordinates with RTTB WorldCoordinate3D worldCoordinateTarget; if (_geoInfoTargetImage.indexToWorldCoordinate(aIndex, worldCoordinateTarget)) { //transform coordinates WorldCoordinate3D worldCoordinateMoving; _spTransformation->transformInverse(worldCoordinateTarget, worldCoordinateMoving); //Use Interpolation to compute dose at mappedImage try { return _spInterpolation->getValue(worldCoordinateMoving); } //Mapped outside of image? Check if padding is allowed catch (core::MappingOutsideOfImageException& /*e*/) { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); } } catch (core::Exception& e) { std::cout << e.what() << std::endl; return -1; } } //ok, if that fails, throw exception. Makes no sense to go further else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); } } } }//end namespace interpolation }//end namespace rttb diff --git a/code/io/dicom/rttbDicomDoseAccessor.cpp b/code/io/dicom/rttbDicomDoseAccessor.cpp index a018e77..dee085f 100644 --- a/code/io/dicom/rttbDicomDoseAccessor.cpp +++ b/code/io/dicom/rttbDicomDoseAccessor.cpp @@ -1,285 +1,287 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ +#include "rttbDicomDoseAccessor.h" #include "drtdose.h" +#include "rttbDcmrtException.h" +#include "rttbInvalidDoseException.h" #include -#include #include "rttbDicomDoseAccessor.h" #include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbDcmrtException.h" #include "rttbIndexOutOfBoundsException.h" namespace rttb { namespace io { namespace dicom { DicomDoseAccessor::~DicomDoseAccessor() = default; DicomDoseAccessor::DicomDoseAccessor(DRTDoseIODPtr aDRTDoseIODP, DcmItemPtr aDcmDataset) { _dose = aDRTDoseIODP; _dataSet = aDcmDataset; OFString uid; _dose->getSeriesInstanceUID(uid); _doseUID = uid.c_str(); this->begin(); } bool DicomDoseAccessor::begin() { assembleGeometricInfo(); doseData.clear(); OFString doseGridScalingStr; this->_dose->getDoseGridScaling(doseGridScalingStr); try { _doseGridScaling = boost::lexical_cast(doseGridScalingStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Dose grid scaling not readable or = 0!") ; } OFCondition status; unsigned long count; const Uint16* pixelData; status = _dataSet->findAndGetUint16Array(DcmTagKey(0x7fe0, 0x0010), pixelData, &count); if (status.good()) { for (unsigned int i = 0; i < static_cast(this->_geoInfo.getNumberOfVoxels()); i++) { this->doseData.push_back(pixelData[i]); } return true; } else { throw io::dicom::DcmrtException("Read Pixel Data (7FE0,0010) failed!"); } } void DicomDoseAccessor::assembleGeometricInfo() { Uint16 temp = 0; this->_dose->getColumns(temp); _geoInfo.setNumColumns(temp); temp = 0; this->_dose->getRows(temp); _geoInfo.setNumRows(temp); if (_geoInfo.getNumColumns() == 0 || _geoInfo.getNumRows() == 0) { throw core::InvalidDoseException("Empty dicom dose!") ; } OFString numberOfFramesStr; OFString imageOrientationRowX, imageOrientationRowY, imageOrientationRowZ; OFString imageOrientationColumnX, imageOrientationColumnY, imageOrientationColumnZ; WorldCoordinate3D imageOrientationRow; WorldCoordinate3D imageOrientationColumn; try { this->_dose->getNumberOfFrames(numberOfFramesStr); _geoInfo.setNumSlices(boost::lexical_cast(numberOfFramesStr.c_str())); _dose->getImageOrientationPatient(imageOrientationRowX, 0); _dose->getImageOrientationPatient(imageOrientationRowY, 1); _dose->getImageOrientationPatient(imageOrientationRowZ, 2); _dose->getImageOrientationPatient(imageOrientationColumnX, 3); _dose->getImageOrientationPatient(imageOrientationColumnY, 4); _dose->getImageOrientationPatient(imageOrientationColumnZ, 5); imageOrientationRow(0) = boost::lexical_cast(imageOrientationRowX.c_str()); imageOrientationRow(1) = boost::lexical_cast(imageOrientationRowY.c_str()); imageOrientationRow(2) = boost::lexical_cast(imageOrientationRowZ.c_str()); imageOrientationColumn(0) = boost::lexical_cast(imageOrientationColumnX.c_str()); imageOrientationColumn(1) = boost::lexical_cast(imageOrientationColumnY.c_str()); imageOrientationColumn(2) = boost::lexical_cast(imageOrientationColumnZ.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("boost::lexical_cast failed! Empty dicom dose!") ; } /*Get orientation*/ OrientationMatrix orientation; orientation(0, 0) = imageOrientationRow.x(); orientation(1, 0) = imageOrientationRow.y(); orientation(2, 0) = imageOrientationRow.z(); orientation(0, 1) = imageOrientationColumn.x(); orientation(1, 1) = imageOrientationColumn.y(); orientation(2, 1) = imageOrientationColumn.z(); WorldCoordinate3D perpendicular = imageOrientationRow.cross(imageOrientationColumn); orientation(0, 2) = perpendicular.x(); orientation(1, 2) = perpendicular.y(); orientation(2, 2) = perpendicular.z(); _geoInfo.setOrientationMatrix(orientation); OFString imagePositionX, imagePositionY, imagePositionZ; _dose->getImagePositionPatient(imagePositionX, 0); _dose->getImagePositionPatient(imagePositionY, 1); _dose->getImagePositionPatient(imagePositionZ, 2); WorldCoordinate3D imagePositionPatient; try { imagePositionPatient(0) = boost::lexical_cast(imagePositionX.c_str()); imagePositionPatient(1) = boost::lexical_cast(imagePositionY.c_str()); imagePositionPatient(2) = boost::lexical_cast(imagePositionZ.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Can not read image position X/Y/Z!") ; } _geoInfo.setImagePositionPatient(imagePositionPatient); /*Get spacing*/ SpacingVectorType3D spacingVector; OFString pixelSpacingRowStr, pixelSpacingColumnStr, sliceThicknessStr; _dose->getPixelSpacing(pixelSpacingRowStr, 0); _dose->getPixelSpacing(pixelSpacingColumnStr, 1); try { spacingVector(1) = boost::lexical_cast(pixelSpacingRowStr.c_str()); spacingVector(0) = boost::lexical_cast(pixelSpacingColumnStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Can not read Pixel Spacing Row/Column!") ; } _geoInfo.setSpacing(spacingVector); if (_geoInfo.getSpacing()(0) == 0 || _geoInfo.getSpacing()(1) == 0) { throw core::InvalidDoseException("Pixel spacing is 0!"); } _dose->getSliceThickness(sliceThicknessStr); try { spacingVector(2) = boost::lexical_cast(sliceThicknessStr.c_str()); } catch (boost::bad_lexical_cast&) { spacingVector(2) = 0 ; } if (spacingVector(2) == 0) { OFVector gridFrameOffsetVector; _dose->getGridFrameOffsetVector(gridFrameOffsetVector); if (gridFrameOffsetVector.size() >= 2) { spacingVector(2) = gridFrameOffsetVector.at(1) - gridFrameOffsetVector.at( 0); //read slice thickness from GridFrameOffsetVector (3004,000c) } if (spacingVector(2) == 0) { OFCondition status; DcmItem doseitem; OFString pixelSpacingBetweenSlices; status = _dose->write(doseitem); if (status.good()) { status = doseitem.findAndGetOFString(DcmTagKey(0x0018, 0x0088), pixelSpacingBetweenSlices); try { spacingVector(2) = boost::lexical_cast (pixelSpacingBetweenSlices.c_str());//read slice thickness from PixelSpacingBetweenSlices (0018,0088) } catch (boost::bad_lexical_cast&) { spacingVector(2) = 0 ; } } //if no useful tags to compute slicing -> set slice thickness to spacingVector(0) if (spacingVector(2) == 0) { std::cerr << "sliceThickness == 0! It wird be replaced with pixelSpacingRow=" << _geoInfo.getSpacing()(0) << "!" << std::endl; spacingVector(2) = spacingVector(0); } } } _geoInfo.setSpacing(spacingVector); } GenericValueType DicomDoseAccessor::getValueAt(const VoxelGridID aID) const { return doseData.at(aID) * _doseGridScaling; } GenericValueType DicomDoseAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getValueAt(aVoxelGridID); } else { return -1; } } } } } diff --git a/code/io/dicom/rttbDicomFileDoseAccessorGenerator.cpp b/code/io/dicom/rttbDicomFileDoseAccessorGenerator.cpp index 422d923..8746fa0 100644 --- a/code/io/dicom/rttbDicomFileDoseAccessorGenerator.cpp +++ b/code/io/dicom/rttbDicomFileDoseAccessorGenerator.cpp @@ -1,103 +1,99 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomDoseAccessor.h" -#include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" -#include "rttbDcmrtException.h" -#include "rttbIndexOutOfBoundsException.h" #include "rttbDicomFileReaderHelper.h" #include "rttbInvalidParameterException.h" #include "rttbUtils.h" namespace rttb { namespace io { namespace dicom { DicomFileDoseAccessorGenerator::~DicomFileDoseAccessorGenerator() = default; DicomFileDoseAccessorGenerator::DicomFileDoseAccessorGenerator(FileNameType aDICOMRTDoseFileName) { _dicomDoseFileName = aDICOMRTDoseFileName; } core::DoseAccessorGeneratorInterface::DoseAccessorPointer DicomFileDoseAccessorGenerator::generateDoseAccessor() { std::vector fileVector; //if a file if (core::isFile(_dicomDoseFileName)) { fileVector.push_back(_dicomDoseFileName); } //if a directory else if (core::isDirectory(_dicomDoseFileName)) { rttb::io::dicom::Modality doseModality = {rttb::io::dicom::Modality::RTDOSE}; fileVector = getFileNamesWithSameUID(_dicomDoseFileName, doseModality); } else { throw rttb::core::InvalidParameterException("Invalid file/directory name!"); } if (fileVector.size() < 1) { throw rttb::core::InvalidParameterException("There is no structure set files in the directory!"); } OFCondition status; DcmFileFormat fileformat; DRTDoseIODPtr dose = boost::make_shared(); status = fileformat.loadFile(fileVector.at(0).c_str()); if (!status.good()) { std::cerr << "Error: load rtdose loadFile(" << fileVector.at(0) << ") failed!" << std::endl; throw core::InvalidDoseException("Invalid dicom dose!"); } DcmItemPtr dataSet = boost::make_shared(*fileformat.getDataset()); status = dose->read(*dataSet); if (!status.good()) { std::cerr << "Error: read DRTDoseIOD failed!" << std::endl; throw core::InvalidDoseException("Invalid dicom dose!"); } _doseAccessor = boost::make_shared(dose, dataSet); return _doseAccessor; } } } } diff --git a/code/io/dicom/rttbDicomFileDoseAccessorWriter.cpp b/code/io/dicom/rttbDicomFileDoseAccessorWriter.cpp index 9543f80..ef0aeb9 100644 --- a/code/io/dicom/rttbDicomFileDoseAccessorWriter.cpp +++ b/code/io/dicom/rttbDicomFileDoseAccessorWriter.cpp @@ -1,250 +1,247 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include -#include #include "rttbDicomFileDoseAccessorWriter.h" #include "rttbInvalidDoseException.h" -#include "rttbGeometricInfo.h" #include "rttbGenericDoseIterator.h" #include "rttbDoseStatisticsCalculator.h" namespace rttb { namespace io { namespace dicom { DicomFileDoseAccessorWriter::DicomFileDoseAccessorWriter() { _doseIOD = boost::make_shared(); _dataset = _fileformat.getDataset(); } void DicomFileDoseAccessorWriter::setFileName(DICOMRTFileNameString aFileName) { _fileName = aFileName; } bool DicomFileDoseAccessorWriter::process() { OFCondition status; /* Prepare dcmtk */ DcmItem* dcm_item = nullptr; //get geometric info rttb::core::GeometricInfo geometricInfo = _doseAccessor->getGeometricInfo(); /* ----------------------------------------------------------------- */ /* Part 1 -- General header */ /* ----------------------------------------------------------------- */ OFString CreationUID(_doseAccessor->getUID().c_str()); _dataset->putAndInsertString(DCM_ImageType, R"(DERIVED\SECONDARY\REFORMATTED)"); _dataset->putAndInsertOFStringArray(DCM_InstanceCreationDate, "");//Creation Date _dataset->putAndInsertOFStringArray(DCM_InstanceCreationTime, "");//Creation Time _dataset->putAndInsertOFStringArray(DCM_InstanceCreatorUID, CreationUID); _dataset->putAndInsertString(DCM_SOPClassUID, UID_RTDoseStorage); _dataset->putAndInsertString(DCM_SOPInstanceUID, _doseAccessor->getUID().c_str()); _dataset->putAndInsertOFStringArray(DCM_StudyDate, ""); _dataset->putAndInsertOFStringArray(DCM_StudyTime, ""); _dataset->putAndInsertOFStringArray(DCM_AccessionNumber, ""); _dataset->putAndInsertOFStringArray(DCM_Modality, "RTDOSE"); _dataset->putAndInsertString(DCM_Manufacturer, "RTToolbox"); _dataset->putAndInsertString(DCM_InstitutionName, ""); _dataset->putAndInsertString(DCM_ReferringPhysicianName, ""); _dataset->putAndInsertString(DCM_StationName, ""); _dataset->putAndInsertString(DCM_ManufacturerModelName, "RTToolbox"); /* (0008,1140) DCM_ReferencedImageSequence -- MIM likes this */ dcm_item = nullptr; _dataset->findOrCreateSequenceItem( DCM_ReferencedImageSequence, dcm_item, -2); dcm_item->putAndInsertString(DCM_ReferencedSOPClassUID, UID_CTImageStorage); dcm_item->putAndInsertString(DCM_ReferencedSOPInstanceUID, ""); _dataset->putAndInsertString(DCM_PatientName, ""); _dataset->putAndInsertString(DCM_PatientID, ""); _dataset->putAndInsertString(DCM_PatientBirthDate, ""); _dataset->putAndInsertString(DCM_PatientSex, "O"); _dataset->putAndInsertString(DCM_SliceThickness, boost::lexical_cast(geometricInfo.getSpacing()(2)).c_str()); _dataset->putAndInsertString(DCM_SoftwareVersions, ""); _dataset->putAndInsertString(DCM_StudyInstanceUID, ""); _dataset->putAndInsertString(DCM_SeriesInstanceUID, ""); _dataset->putAndInsertString(DCM_StudyID, "10001"); _dataset->putAndInsertString(DCM_SeriesNumber, ""); _dataset->putAndInsertString(DCM_InstanceNumber, "1"); /* GCS FIX: PatientOrientation */ std::ostringstream sstr; sstr << geometricInfo.getImagePositionPatient().x() << R"(\)" << geometricInfo.getImagePositionPatient().y() << R"(\)" << geometricInfo.getImagePositionPatient().z(); _dataset->putAndInsertString(DCM_PatientOrientation, "L/P"); _dataset->putAndInsertString(DCM_ImagePositionPatient, sstr.str().c_str()); auto orientationMatrix = geometricInfo.getOrientationMatrix(); sstr.str(""); sstr << orientationMatrix(0,0) << R"(\)" << orientationMatrix(1,0) << R"(\)" << orientationMatrix(2,0) << R"(\)" << orientationMatrix(0,1) << R"(\)" << orientationMatrix(1,1) << R"(\)" << orientationMatrix(2,1); _dataset->putAndInsertString(DCM_ImageOrientationPatient, sstr.str().c_str()); _dataset->putAndInsertString(DCM_FrameOfReferenceUID, ""); _dataset->putAndInsertString(DCM_SamplesPerPixel, "1"); _dataset->putAndInsertString(DCM_PhotometricInterpretation, "MONOCHROME2"); sstr.str(""); sstr << geometricInfo.getNumSlices(); _dataset->putAndInsertString(DCM_NumberOfFrames, sstr.str().c_str()); /* GCS FIX: Add FrameIncrementPointer */ _dataset->putAndInsertString(DCM_FrameIncrementPointer, "(3004,000c)"); _dataset->putAndInsertUint16(DCM_Rows, static_cast(geometricInfo.getNumRows())); _dataset->putAndInsertUint16(DCM_Columns, static_cast(geometricInfo.getNumColumns())); sstr.str(""); sstr << geometricInfo.getSpacing()(1) << R"(\)" << geometricInfo.getSpacing()(0); _dataset->putAndInsertString(DCM_PixelSpacing, sstr.str().c_str()); _dataset->putAndInsertString(DCM_BitsAllocated, "32"); _dataset->putAndInsertString(DCM_BitsStored, "32"); _dataset->putAndInsertString(DCM_HighBit, "31"); _dataset->putAndInsertString(DCM_DoseUnits, "GY"); _dataset->putAndInsertString(DCM_DoseSummationType, "PLAN"); sstr.str("0"); for (unsigned int i = 1; i < geometricInfo.getNumSlices(); i++) { sstr << R"(\)" << i* geometricInfo.getSpacing()(2); } _dataset->putAndInsertString(DCM_GridFrameOffsetVector, sstr.str().c_str()); /* We need to convert image to uint16_t, but first we need to scale it so that the maximum dose fits in a 16-bit unsigned integer. Compute an appropriate scaling factor based on the maximum dose. */ /* Find the maximum value in the image */ boost::shared_ptr spTestDoseIterator = boost::make_shared(_doseAccessor); rttb::core::GenericDoseIterator::DoseIteratorPointer spDoseIterator(spTestDoseIterator); rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator); auto doseStatistics = myDoseStatsCalculator.calculateDoseStatistics(); double maxDose = doseStatistics->getMaximum(); /* Find scale factor */ double dose_scale; dose_scale = maxDose / PixelDataMaxValue; /* Scale the image and add scale factor to _dataset */ sstr.str(""); sstr << dose_scale; _dataset->putAndInsertString(DCM_DoseGridScaling, sstr.str().c_str()); /* (300c,0002) ReferencedRTPlanSequence -- for future expansion */ dcm_item = nullptr; _dataset->findOrCreateSequenceItem( DCM_ReferencedRTPlanSequence, dcm_item, -2); dcm_item->putAndInsertString(DCM_ReferencedSOPClassUID, UID_RTPlanStorage); dcm_item->putAndInsertString(DCM_ReferencedSOPInstanceUID, ""); /* (300c,0060) DCM_ReferencedStructureSetSequence -- MIM likes this */ dcm_item = nullptr; _dataset->findOrCreateSequenceItem( DCM_ReferencedStructureSetSequence, dcm_item, -2); dcm_item->putAndInsertString(DCM_ReferencedSOPClassUID, UID_RTStructureSetStorage); dcm_item->putAndInsertString(DCM_ReferencedSOPInstanceUID, ""); /* Convert image bytes to integer, then add to _dataset */ Uint16* pixelData; auto pixelCount = static_cast(geometricInfo.getNumRows() * geometricInfo.getNumColumns() * geometricInfo.getNumSlices()); pixelData = new Uint16[pixelCount]; for (unsigned int i = 0; i < pixelCount; ++i) { double doseValue = _doseAccessor->getValueAt(i); double pixelValue = doseValue / dose_scale; if (pixelValue > PixelDataMaxValue) { pixelValue = PixelDataMaxValue; } pixelData[i] = boost::numeric_cast(pixelValue); } status = _dataset->putAndInsertUint16Array(DCM_PixelData, pixelData, pixelCount); if (!status.good()) { throw core::InvalidDoseException("Error: put and insert pixel data failed!"); } //Write dose to file status = _fileformat.saveFile(_fileName.c_str(), EXS_LittleEndianExplicit); if (status.bad()) { std::cerr << "Error: cannot write DICOM RTDOSE!" << std::endl; } return true; } }//end namespace itk }//end namespace io }//end namespace rttb diff --git a/code/io/dicom/rttbDicomFileReaderHelper.cpp b/code/io/dicom/rttbDicomFileReaderHelper.cpp index cace20d..3a084f0 100644 --- a/code/io/dicom/rttbDicomFileReaderHelper.cpp +++ b/code/io/dicom/rttbDicomFileReaderHelper.cpp @@ -1,219 +1,216 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDicomFileReaderHelper.h" #include -#include -#include #include "boost/filesystem/operations.hpp" -#include "boost/filesystem/path.hpp" #include "rttbDcmrtException.h" #include "rttbInvalidParameterException.h" #include "rttbUtils.h" namespace rttb { namespace io { namespace dicom { OFString getModality(DcmDataSetPtr aDcmDataSet) { OFString modality; OFCondition status; status = aDcmDataSet->findAndGetOFString(DcmTagKey(0x0008, 0x0060), modality); if (!status.good()) { throw DcmrtException("Error: get modality failed!"); } return modality; } OFString getUID(DcmDataSetPtr aDcmDataSet) { OFString uid; OFCondition status; status = aDcmDataSet->findAndGetOFString(DcmTagKey(0x0020, 0x000e), uid); if (!status.good()) { throw DcmrtException("Error: get uid failed!"); } return uid; } std::vector getFileNamesWithSameUID(FileNameType aDirName, Modality aModality) { std::vector fileNameVector; std::string modalityStrArray[] = {"RTDOSE", "RTSTRUCT", "RTPLAN"}; boost::filesystem::path path = boost::filesystem::path(aDirName); OFCondition status; DcmFileFormat fileformat; IDType uid = ""; DcmDataSetPtr datasetPtr; if (aModality.Type != 1 && aModality.Type != 2 && aModality.Type != 3) { throw core::InvalidParameterException("Error: invalid modality! The modality should be RTDOSE(1)/RTSTRUCT(2)/RTPLAN(3)."); } //if a directory if (core::isDirectory(aDirName)) { boost::filesystem::directory_iterator end_iter; bool isFirst = true; for (boost::filesystem::directory_iterator dir_itr(path); dir_itr != end_iter; ++dir_itr) { if (boost::filesystem::is_regular_file(dir_itr->status())) { boost::filesystem::path filePath(dir_itr->path().filename().string()); filePath = boost::filesystem::system_complete(dir_itr->path()); status = fileformat.loadFile(filePath.string().c_str()); if (!status.good()) { throw DcmrtException("Error: load dose file " + filePath.string() + " failed!"); } datasetPtr = boost::make_shared(*fileformat.getDataset()); OFString modalityOFS = getModality(datasetPtr); for (unsigned int i = 0; i < 3; i++) { if (aModality.Type == (i + 1) && modalityOFS == modalityStrArray[i].c_str()) { OFString currentUID = getUID(datasetPtr); //get the first uid of the given modality if (isFirst) { uid = currentUID.c_str(); isFirst = false; } if (uid == currentUID.c_str()) { fileNameVector.emplace_back(filePath.string().c_str()); } break; } } } } } else if (core::isFile(aDirName)) { std::cout << "Important: the given name " + aDirName + " is a file name, not a directory name. Given modality will be ignored, use the modality of the file." << std::endl; fileNameVector = getFileNames(aDirName); } else { throw core::InvalidParameterException("Error: file/directory does not exist!"); } return fileNameVector; } std::vector getFileNames(FileNameType aFileName) { if (!core::isFile(aFileName)) { throw core::InvalidParameterException("Error: file does not exist!"); } std::vector fileNameVector; boost::filesystem::path path = boost::filesystem::path(aFileName); OFCondition status; DcmFileFormat fileformat; DcmDataSetPtr datasetPtr; OFString modality;//modality of the given file OFString uid;//uid of the given file status = fileformat.loadFile(aFileName.c_str()); if (!status.good()) { throw DcmrtException("Error: fileformat.loadFile failed!"); } datasetPtr = boost::make_shared(*fileformat.getDataset()); modality = getModality(datasetPtr); uid = getUID(datasetPtr); //get parent directory boost::filesystem::path parentDir = path.parent_path(); if (boost::filesystem::is_directory(parentDir)) { boost::filesystem::directory_iterator end_iter; for (boost::filesystem::directory_iterator dir_itr(parentDir); dir_itr != end_iter; ++dir_itr) { if (boost::filesystem::is_regular_file(dir_itr->status())) { boost::filesystem::path currentFilePath(dir_itr->path().filename().string()); currentFilePath = boost::filesystem::system_complete(dir_itr->path()); status = fileformat.loadFile(currentFilePath.string().c_str()); if (!status.good()) { throw DcmrtException("Error: load dose fileformat.loadFile failed!"); } datasetPtr = boost::make_shared(*fileformat.getDataset()); OFString currentModality = getModality(datasetPtr); OFString currentUID = getUID(datasetPtr); //if the same modality if (modality == currentModality && uid == currentUID) { fileNameVector.emplace_back(currentFilePath.string().c_str()); } } } } return fileNameVector; } } } } diff --git a/code/io/dicom/rttbDicomFileStructureSetGenerator.cpp b/code/io/dicom/rttbDicomFileStructureSetGenerator.cpp index 05252d8..b876a1a 100644 --- a/code/io/dicom/rttbDicomFileStructureSetGenerator.cpp +++ b/code/io/dicom/rttbDicomFileStructureSetGenerator.cpp @@ -1,104 +1,100 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include -#include -#include #include "rttbInvalidParameterException.h" -#include "rttbStructure.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbDicomIODStructureSetGenerator.h" #include "rttbDcmrtException.h" #include "rttbDicomFileReaderHelper.h" #include "rttbUtils.h" namespace rttb { namespace io { namespace dicom { DicomFileStructureSetGenerator::DicomFileStructureSetGenerator(DICOMRTFileNameString aDICOMRTStrSetFileName) { _fileName = aDICOMRTStrSetFileName; } DicomFileStructureSetGenerator::~DicomFileStructureSetGenerator() = default; DicomFileStructureSetGenerator::StructureSetPointer DicomFileStructureSetGenerator::generateStructureSet() { std::vector fileVector; //if a file if (core::isFile(_fileName)) { fileVector.push_back(_fileName); } //if a directory else if (core::isDirectory(_fileName)) { rttb::io::dicom::Modality strModality = {rttb::io::dicom::Modality::RTSTRUCT}; fileVector = getFileNamesWithSameUID(_fileName, strModality); } else { throw rttb::core::InvalidParameterException("Invalid file/directory name!"); } if (fileVector.size() < 1) { throw rttb::core::InvalidParameterException("There is no structure set files in the directory!"); } OFCondition status; DcmFileFormat fileformat; DRTStrSetIODPtr drtStrSetIODPtr = boost::make_shared(); //get the first structure set file status = fileformat.loadFile(fileVector.at(0).c_str()); if (!status.good()) { throw DcmrtException("Load rt structure set loadFile() failed!"); } status = drtStrSetIODPtr->read(*fileformat.getDataset()); if (!status.good()) { throw DcmrtException("Read DRTStructureSetIOD DRTStructureSetIOD.read() failed!"); } io::dicom::DicomIODStructureSetGenerator iodGenerator(drtStrSetIODPtr); iodGenerator.setStructureLabelFilterActive(this->getStructureLabelFilterActive()); iodGenerator.setFilterRegEx(this->getFilterRegEx()); return iodGenerator.generateStructureSet(); } }//end namespace dicom }//end namespace io }//end namespace rttb diff --git a/code/io/dicom/rttbDicomIODDoseAccessorGenerator.cpp b/code/io/dicom/rttbDicomIODDoseAccessorGenerator.cpp index 700514c..bda7d93 100644 --- a/code/io/dicom/rttbDicomIODDoseAccessorGenerator.cpp +++ b/code/io/dicom/rttbDicomIODDoseAccessorGenerator.cpp @@ -1,65 +1,61 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbDicomIODDoseAccessorGenerator.h" #include "rttbDicomDoseAccessor.h" -#include "rttbNullPointerException.h" -#include "rttbInvalidDoseException.h" #include "rttbDcmrtException.h" -#include "rttbIndexOutOfBoundsException.h" namespace rttb { namespace io { namespace dicom { DicomIODDoseAccessorGenerator::~DicomIODDoseAccessorGenerator() = default; DicomIODDoseAccessorGenerator::DicomIODDoseAccessorGenerator(DRTDoseIODPtr aDRTDoseIODP) { _doseIODPtr = aDRTDoseIODP; } core::DoseAccessorGeneratorInterface::DoseAccessorPointer DicomIODDoseAccessorGenerator::generateDoseAccessor() { DcmItem doseitem; OFCondition status; status = _doseIODPtr->write(doseitem);//write DoseIOD to DcmItem to get pixel data if (status.good()) { DcmItemPtr dataSetPtr = boost::make_shared(doseitem); _doseAccessor = boost::make_shared(_doseIODPtr, dataSetPtr); return _doseAccessor; } else { throw io::dicom::DcmrtException("Write DICOM RT Dose to DcmItem failed!"); } } } } } diff --git a/code/io/dicom/rttbDicomIODStructureSetGenerator.cpp b/code/io/dicom/rttbDicomIODStructureSetGenerator.cpp index 70ff824..7f395d6 100644 --- a/code/io/dicom/rttbDicomIODStructureSetGenerator.cpp +++ b/code/io/dicom/rttbDicomIODStructureSetGenerator.cpp @@ -1,201 +1,195 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include -#include -#include #include #include - #include -#include "rttbStructure.h" #include "rttbDicomIODStructureSetGenerator.h" -#include "rttbDcmrtException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace io { namespace dicom { DicomIODStructureSetGenerator::DicomIODStructureSetGenerator(DRTStrSetIODPtr aDRTStructureSetIOD) { _drtStrSetIOD = aDRTStructureSetIOD; } void DicomIODStructureSetGenerator::readStrSet() { OFString uid; _drtStrSetIOD->getSeriesInstanceUID(uid); _UID = uid.c_str(); OFString uid2; _drtStrSetIOD->getPatientID(uid2); _patientUID = uid2.c_str(); DRTStructureSetROISequence* rois = &_drtStrSetIOD->getStructureSetROISequence(); //generate map of relevant ROIs std::map filteredROIs; std::regex e(this->getFilterRegEx()); for (unsigned long i = 0; i < rois->getNumberOfItems(); i++) { DRTStructureSetROISequence::Item* roisItem = &rois->getItem(i); OFString roiNumber; roisItem->getROINumber(roiNumber); OFString ofRoiName; roisItem->getROIName(ofRoiName); std::string roiName(ofRoiName.c_str()); //replace wrongly 'á' character by ' ' in ROI name correctSpacesInROIName(roiName); if (!this->getStructureLabelFilterActive() || std::regex_match(roiName, e)) { filteredROIs.emplace(roiNumber, roiName); } } /*A structure is a DRTROIContourSequence::Item. Each Item defines a roi. Each ROI contains a sequence of one or more contours, where a contour is either a single point (for a point ROI) or more than one point (representing an open or closed polygon). */ DRTROIContourSequence* rcs; rcs = &_drtStrSetIOD->getROIContourSequence(); DRTROIContourSequence::Item* rcsItem; long numberOfStructures = rcs->getNumberOfItems(); bool isEmpty = rcs->isEmpty(); if (numberOfStructures == 0 || isEmpty) { throw core::InvalidParameterException("Empty Structure Set!") ; } int structureNo = 0; for (rcs->gotoFirstItem(); (rcs->getCurrentItem(rcsItem)).good(); rcs->gotoNextItem()) { OFString refROINumber; rcsItem->getReferencedROINumber(refROINumber); //check if ROI number is in the filtered ROIS if (filteredROIs.find(refROINumber) != filteredROIs.end()) { DRTContourSequence* cs; cs = &rcsItem->getContourSequence(); unsigned long no2 = cs->getNumberOfItems(); PolygonSequenceType structureVector; for (unsigned long j = 0; j < no2; j++) { /*DRTContourSequence::Item represents a contour (either a single point (for a point ROI) or more than one point (representing an open or closed polygon))*/ DRTContourSequence::Item* csItem; csItem = &cs->getItem(j); OFString contourData; OFString numberOfContourPoints; csItem->getNumberOfContourPoints(numberOfContourPoints); unsigned int numberOfContourPointsInt; std::stringstream is(numberOfContourPoints.c_str()); is >> numberOfContourPointsInt; OFString countourNumber; csItem->getContourNumber(countourNumber); PolygonType contourVector; char* pEnd; for (unsigned int k = 0; k < numberOfContourPointsInt; k++) { WorldCoordinate3D point; for (unsigned int i = 0; i < 3; i++) { csItem->getContourData(contourData, k * 3 + i); WorldCoordinate value = strtod(contourData.c_str(), &pEnd); if (*pEnd != '\0') { throw core::InvalidParameterException("Contour data not readable!") ; } if (i == 0) { point(0) = value; } else if (i == 1) { point(1) = value; } else { point(2) = value; } } contourVector.push_back(point); } structureVector.push_back(contourVector); } boost::shared_ptr spStruct = boost::make_shared(structureVector); StructTypePointer str(spStruct); str->setLabel(filteredROIs[refROINumber]); std::cout << filteredROIs[refROINumber].c_str() << " read" << std::endl; std::stringstream sstr; sstr << structureNo; str->setUID(sstr.str()); _strVector.push_back(str); } ++structureNo; } } DicomIODStructureSetGenerator::~DicomIODStructureSetGenerator() = default; DicomIODStructureSetGenerator::StructureSetPointer DicomIODStructureSetGenerator::generateStructureSet() { this->readStrSet(); return boost::make_shared(_strVector, _patientUID, _UID); } void DicomIODStructureSetGenerator::correctSpacesInROIName(std::string& roiName) { for (auto& character : roiName) { if (character == -96) { character = ' '; } } } }//end namespace dicom }//end namespace io }//end namespace rttb diff --git a/code/io/helax/rttbDicomHelaxDoseAccessor.cpp b/code/io/helax/rttbDicomHelaxDoseAccessor.cpp index 733cc4f..4679414 100644 --- a/code/io/helax/rttbDicomHelaxDoseAccessor.cpp +++ b/code/io/helax/rttbDicomHelaxDoseAccessor.cpp @@ -1,319 +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. // //------------------------------------------------------------------------ -#include "rttbDicomHelaxDoseAccessor.h" - #include "boost/lexical_cast.hpp" -#include "boost/filesystem/operations.hpp" -#include "boost/filesystem/path.hpp" -#include "boost/numeric/ublas/matrix.hpp" #include -#include "rttbInvalidDoseException.h" +#include "rttbDicomHelaxDoseAccessor.h" #include "rttbInvalidDoseException.h" #include "rttbDcmrtException.h" +#include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace io { namespace helax { DicomHelaxDoseAccessor::~DicomHelaxDoseAccessor() = default; DicomHelaxDoseAccessor::DicomHelaxDoseAccessor(std::vector aDICOMRTDoseVector) { for (const auto & i : aDICOMRTDoseVector) { _doseVector.push_back(i); } this->begin(); } bool DicomHelaxDoseAccessor::begin() { if (_doseVector.size() == 0) { throw core::InvalidParameterException(" The size of aDICOMRTDoseVector is 0!"); } assembleGeometricInfo(); _doseData.clear(); OFString doseGridScalingStr; _doseVector.at(0)->getDoseGridScaling( doseGridScalingStr);//get the first dose grid scaling as _doseGridScaling try { _doseGridScaling = boost::lexical_cast(doseGridScalingStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Dose grid scaling not readable or = 0!") ; } for (auto dose : _doseVector) { OFString currentDoseGridScalingStr; dose->getDoseGridScaling(currentDoseGridScalingStr); double currentDoseGridScaling; try { currentDoseGridScaling = boost::lexical_cast(currentDoseGridScalingStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Dose grid scaling not readable or = 0!") ; } OFCondition status; DcmFileFormat fileformat; DcmItem doseitem; status = dose->write(doseitem); if (status.good()) { unsigned long count; const Uint16* pixelData; status = doseitem.findAndGetUint16Array(DcmTagKey(0x7fe0, 0x0010), pixelData, &count); if (status.good()) { for (unsigned int j = 0; j < static_cast(_geoInfo.getNumColumns()*_geoInfo.getNumRows()); j++) { auto data = static_cast(pixelData[j] * currentDoseGridScaling / _doseGridScaling); this->_doseData.push_back(data); //recalculate dose data } } else { throw dicom::DcmrtException("Read Pixel Data (7FE0,0010) failed!"); } } else { throw dicom::DcmrtException("Read DICOM-RT Dose file failed!"); } } return true; } void DicomHelaxDoseAccessor::assembleGeometricInfo() { DRTDoseIODPtr dose = _doseVector.at(0); Uint16 temp = 0; dose->getColumns(temp); _geoInfo.setNumColumns(temp); temp = 0; dose->getRows(temp); _geoInfo.setNumRows(temp); OFString numberOfFramesStr; dose->getNumberOfFrames(numberOfFramesStr); if (!numberOfFramesStr.empty()) { _geoInfo.setNumSlices(boost::lexical_cast(numberOfFramesStr.c_str())); } else { _geoInfo.setNumSlices((VoxelGridDimensionType)_doseVector.size()); } if (_geoInfo.getNumColumns() == 0 || _geoInfo.getNumRows() == 0 || _geoInfo.getNumSlices() == 0) { throw core::InvalidDoseException("Empty dicom dose!") ; } OFString imageOrientationRowX; dose->getImageOrientationPatient(imageOrientationRowX, 0); OFString imageOrientationRowY; dose->getImageOrientationPatient(imageOrientationRowY, 1); OFString imageOrientationRowZ; dose->getImageOrientationPatient(imageOrientationRowZ, 2); OFString imageOrientationColumnX; dose->getImageOrientationPatient(imageOrientationColumnX, 3); OFString imageOrientationColumnY; dose->getImageOrientationPatient(imageOrientationColumnY, 4); OFString imageOrientationColumnZ; dose->getImageOrientationPatient(imageOrientationColumnZ, 5); WorldCoordinate3D imageOrientationRow; WorldCoordinate3D imageOrientationColumn; try { imageOrientationRow(0) = boost::lexical_cast(imageOrientationRowX.c_str()); imageOrientationRow(1) = boost::lexical_cast(imageOrientationRowY.c_str()); imageOrientationRow(2) = boost::lexical_cast(imageOrientationRowZ.c_str()); imageOrientationColumn(0) = boost::lexical_cast(imageOrientationColumnX.c_str()); imageOrientationColumn(1) = boost::lexical_cast(imageOrientationColumnY.c_str()); imageOrientationColumn(2) = boost::lexical_cast(imageOrientationColumnZ.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("boost::lexical_cast WorldCoordinate failed! Can not read image orientation X/Y/Z!") ; } OrientationMatrix orientation; orientation(0, 0) = imageOrientationRow.x(); orientation(1, 0) = imageOrientationRow.y(); orientation(2, 0) = imageOrientationRow.z(); orientation(0, 1) = imageOrientationColumn.x(); orientation(1, 1) = imageOrientationColumn.y(); orientation(2, 1) = imageOrientationColumn.z(); WorldCoordinate3D perpendicular = imageOrientationRow.cross(imageOrientationColumn); orientation(0, 2) = perpendicular.x(); orientation(1, 2) = perpendicular.y(); orientation(2, 2) = perpendicular.z(); _geoInfo.setOrientationMatrix(orientation); OFString imagePositionX; dose->getImagePositionPatient(imagePositionX, 0); OFString imagePositionY; dose->getImagePositionPatient(imagePositionY, 1); OFString imagePositionZ; dose->getImagePositionPatient(imagePositionZ, 2); WorldCoordinate3D imagePositionPatient; try { imagePositionPatient(0) = boost::lexical_cast(imagePositionX.c_str()); imagePositionPatient(1) = boost::lexical_cast(imagePositionY.c_str()); imagePositionPatient(2) = boost::lexical_cast(imagePositionZ.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("boost::lexical_cast ImagePosition failed! Can not read image position X/Y/Z!") ; } _geoInfo.setImagePositionPatient(imagePositionPatient); SpacingVectorType3D spacingVector; OFString pixelSpacingRowStr; dose->getPixelSpacing(pixelSpacingRowStr, 0); OFString pixelSpacingColumnStr; dose->getPixelSpacing(pixelSpacingColumnStr, 1); try { spacingVector(1) = boost::lexical_cast(pixelSpacingRowStr.c_str()); spacingVector(0) = boost::lexical_cast(pixelSpacingColumnStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Can not read Pixel Spacing Row/Column!") ; } _geoInfo.setSpacing(spacingVector); if (_geoInfo.getSpacing()(0) == 0 || _geoInfo.getSpacing()(1) == 0) { throw core::InvalidDoseException("Pixel spacing not readable or = 0!"); } OFString sliceThicknessStr; dose->getSliceThickness(sliceThicknessStr); try { spacingVector(2) = boost::lexical_cast(sliceThicknessStr.c_str()); } catch (boost::bad_lexical_cast&) { spacingVector(2) = 0 ;//if no information about slice thickness, set to 0 and calculate it using z coordinate difference between 1. and 2. dose } if (spacingVector(2) == 0) { if (_doseVector.size() > 1) { DRTDoseIODPtr dose2 = _doseVector.at(1);//get the 2. dose OFString imagePositionZ2; dose2->getImagePositionPatient(imagePositionZ2, 2); try { spacingVector(2) = boost::lexical_cast(imagePositionZ2.c_str()) - imagePositionPatient( 2); //caculate slicethickness } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Can not read image position Z of the 2. dose!"); } } else { std::cerr << "sliceThickness == 0! It will be replaced with pixelSpacingRow=" << _geoInfo.getSpacing()(0) << "!" << std::endl; spacingVector(2) = spacingVector(0); } } _geoInfo.setSpacing(spacingVector); } GenericValueType DicomHelaxDoseAccessor::getValueAt(const VoxelGridID aID) const { return _doseData.at(aID) * _doseGridScaling; } GenericValueType DicomHelaxDoseAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getValueAt(aVoxelGridID); } else { return -1; } } } } } diff --git a/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp b/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp index 2c15a74..cc11904 100644 --- a/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp +++ b/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp @@ -1,88 +1,80 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include -#include -#include "boost/filesystem/operations.hpp" -#include "boost/filesystem/path.hpp" -#include "boost/progress.hpp" #include "rttbDicomHelaxFileDoseAccessorGenerator.h" #include "rttbDicomHelaxDoseAccessor.h" #include "rttbInvalidDoseException.h" -#include "rttbDcmrtException.h" -#include "rttbIndexOutOfBoundsException.h" -#include "rttbInvalidParameterException.h" #include "rttbDicomFileReaderHelper.h" namespace rttb { namespace io { namespace helax { DicomHelaxFileDoseAccessorGenerator::~DicomHelaxFileDoseAccessorGenerator() = default; DicomHelaxFileDoseAccessorGenerator::DicomHelaxFileDoseAccessorGenerator( FileNameType aDICOMRTDoseDirName) { _doseDirName = aDICOMRTDoseDirName; } core::DoseAccessorGeneratorInterface::DoseAccessorPointer DicomHelaxFileDoseAccessorGenerator::generateDoseAccessor() { rttb::io::dicom::Modality doseModality = {rttb::io::dicom::Modality::RTDOSE}; std::vector fileVector = rttb::io::dicom::getFileNamesWithSameUID(_doseDirName, doseModality); OFCondition status; DcmFileFormat fileformat; std::vector doseVector; for (auto & i : fileVector) { DRTDoseIODPtr dose = boost::make_shared(); status = fileformat.loadFile(i.c_str()); if (!status.good()) { throw core::InvalidDoseException("Error: load dose fileformat.loadFile failed!"); } status = dose->read(*fileformat.getDataset()); if (!status.good()) { throw core::InvalidDoseException("Error: read DRTDoseIOD failed!"); } doseVector.push_back(dose); } _doseAccessor = boost::make_shared(doseVector); return _doseAccessor; } } } } diff --git a/code/io/itk/itkMaskAccessorImageSource.cpp b/code/io/itk/itkMaskAccessorImageSource.cpp index 2a523f9..6cd5e43 100644 --- a/code/io/itk/itkMaskAccessorImageSource.cpp +++ b/code/io/itk/itkMaskAccessorImageSource.cpp @@ -1,101 +1,98 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "itkMaskAccessorImageSource.h" #include "itkImageRegionIterator.h" -#include "itkImageRegionConstIteratorWithIndex.h" #include "itkProgressReporter.h" -#include "itkExceptionObject.h" -#include "rttbGeometricInfo.h" namespace itk { /** * Constructor */ MaskAccessorImageSource ::MaskAccessorImageSource() = default; void MaskAccessorImageSource ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType threadId) { ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); using OutputImageRegionIteratorType = ImageRegionIterator; OutputImagePointer outputPtr = dynamic_cast< OutputImageType* >(ProcessObject::GetOutput(0)); OutputImageRegionIteratorType outputItr; if (outputPtr) { outputItr = OutputImageRegionIteratorType(outputPtr, outputRegionForThread); for (; !(outputItr.IsAtEnd()); ++outputItr) { OutputImageType::IndexType index = outputItr.GetIndex(); rttb::VoxelGridIndex3D maskIndex(index[0], index[1], index[2]); rttb::VoxelGridID id = 0; if (m_Accessor->getGeometricInfo().convert(maskIndex, id)) { auto finding = m_idMap.find(id); if (finding != m_idMap.end()) { // Set the current pixel outputItr.Set(m_MaskedVoxels->at(finding->second).getRelevantVolumeFraction()); } } else { if (m_FailsOnInvalidIDs) { itkExceptionMacro(<<"invalid Mask index. Index:"<(ProcessObject::GetOutput(0)); outputPtr->FillBuffer(0.0); //The id map approach and working with relevant voxel vector is a workarround till task T22067 is solved and and can be used again. m_MaskedVoxels = m_Accessor->getRelevantVoxelVector(); m_idMap.clear(); for (rttb::core::MaskAccessorInterface::MaskVoxelList::size_type pos = 0; pos < m_MaskedVoxels->size(); ++pos) { m_idMap.insert(std::make_pair((*m_MaskedVoxels)[pos].getVoxelGridID(), pos)); } } } // end namespace itk diff --git a/code/io/itk/rttbITKIOHelper.cpp b/code/io/itk/rttbITKIOHelper.cpp index be8ee1d..095d87f 100644 --- a/code/io/itk/rttbITKIOHelper.cpp +++ b/code/io/itk/rttbITKIOHelper.cpp @@ -1,165 +1,164 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbITKIOHelper.h" -#include "rttbException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace io { namespace itk { ITKImageType::Pointer readITKDoubleImage(FileNameType aITKImageFile, bool isDicom) { ITKImageType::Pointer itkDoubleImage; GenericImageReader::Pointer spReader = GenericImageReader::New(); if (isDicom) { spReader->setSeriesReadStyle(ImageSeriesReadStyle::Type::Dicom); } spReader->setFileName(aITKImageFile); GenericImageReader::GenericOutputImageType::Pointer itkGenericImage; ITKImageType::ConstPointer itkDoubleImageConst; try { unsigned int loadedDimensions; GenericImageReader::LoadedPixelType loadedPixelType; GenericImageReader::LoadedComponentType loadedComponentType; itkGenericImage = spReader->GetOutput(loadedDimensions, loadedPixelType, loadedComponentType); if (loadedDimensions != 3) { throw core::InvalidParameterException("image dimensions != 3. Only dim = 3 supported."); } if (loadedPixelType != ::itk::ImageIOBase::SCALAR) { throw core::InvalidParameterException("image component type != SCALAR. Only SCALAR supported."); } if (loadedComponentType == ::itk::ImageIOBase::DOUBLE) { itkDoubleImage = dynamic_cast(itkGenericImage.GetPointer()); } else { itkDoubleImage = handleGenericImage(itkGenericImage, loadedComponentType); } if (itkDoubleImage.IsNull()) { throw core::InvalidDoseException("Error!!! unable to load input image. File is not existing or has an unsupported format."); } } catch (::itk::ExceptionObject& e) { std::cerr << "Error!!!" << std::endl; std::cerr << e << std::endl; throw rttb::core::InvalidDoseException(e.GetDescription()); } return itkDoubleImage; } ITKImageType::Pointer handleGenericImage( GenericImageReader::GenericOutputImageType* itkGenericImage, ::itk::ImageIOBase::IOComponentType& loadedComponentType) { ITKImageType::Pointer itkDoubleImage; switch (loadedComponentType) { case ::itk::ImageIOBase::UCHAR: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::CHAR: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::USHORT: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::SHORT: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::UINT: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::INT: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::ULONG: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::LONG: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::FLOAT: { itkDoubleImage = doCasting(itkGenericImage); break; } case ::itk::ImageIOBase::DOUBLE: { itkDoubleImage = doCasting(itkGenericImage); break; } default: { throw core::InvalidParameterException("image type unknown"); } } return itkDoubleImage; } }//end namespace itk }//end namespace io }//end namespace rttb diff --git a/code/io/itk/rttbITKImageAccessor.cpp b/code/io/itk/rttbITKImageAccessor.cpp index 9ce23cd..0b75928 100644 --- a/code/io/itk/rttbITKImageAccessor.cpp +++ b/code/io/itk/rttbITKImageAccessor.cpp @@ -1,108 +1,107 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include #include "rttbITKImageAccessor.h" -#include "rttbException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace io { namespace itk { ITKImageAccessor::~ITKImageAccessor() = default; ITKImageAccessor::ITKImageAccessor(ITKImageType::ConstPointer image) { _data = image; if (_data.IsNull()) { throw core::InvalidDoseException("Dose image = 0!") ; } assembleGeometricInfo(); } GenericValueType ITKImageAccessor::getValueAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex; if (_geoInfo.convert(aID, aVoxelGridIndex)) { return getValueAt(aVoxelGridIndex); } else { return -1; } } GenericValueType ITKImageAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { if (_geoInfo.validIndex(aIndex)) { const ITKImageType::IndexType pixelIndex = {{aIndex[0], aIndex[1], aIndex[2]}}; return _data->GetPixel(pixelIndex); } else { return -1; } } void ITKImageAccessor::assembleGeometricInfo() { _geoInfo.setSpacing(SpacingVectorType3D(_data->GetSpacing()[0], _data->GetSpacing()[1], _data->GetSpacing()[2])); if (_geoInfo.getSpacing()[0] == 0 || _geoInfo.getSpacing()[1] == 0 || _geoInfo.getSpacing()[2] == 0) { throw core::InvalidDoseException("Pixel spacing = 0!"); } _geoInfo.setImagePositionPatient(WorldCoordinate3D(_data->GetOrigin()[0], _data->GetOrigin()[1], _data->GetOrigin()[2])); OrientationMatrix OM(0); for (unsigned int col = 0; col < 3; ++col) { for (unsigned int row = 0; row < 3; ++row) { OM(col, row) = _data->GetDirection()(col, row); } } _geoInfo.setOrientationMatrix(OM); _geoInfo.setNumColumns(_data->GetLargestPossibleRegion().GetSize()[0]); _geoInfo.setNumRows(_data->GetLargestPossibleRegion().GetSize()[1]); _geoInfo.setNumSlices(_data->GetLargestPossibleRegion().GetSize()[2]); if (_geoInfo.getNumColumns() == 0 || _geoInfo.getNumRows() == 0 || _geoInfo.getNumSlices() == 0) { throw core::InvalidDoseException("Empty dicom dose!") ; } } } } } diff --git a/code/io/itk/rttbITKImageAccessorConverter.cpp b/code/io/itk/rttbITKImageAccessorConverter.cpp index 64aead4..10190ad 100644 --- a/code/io/itk/rttbITKImageAccessorConverter.cpp +++ b/code/io/itk/rttbITKImageAccessorConverter.cpp @@ -1,105 +1,100 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include - #include "rttbITKImageAccessorConverter.h" -#include "rttbException.h" -#include "rttbInvalidDoseException.h" -#include "rttbGeometricInfo.h" #include "itkDoseAccessorImageFilter.h" namespace rttb { namespace io { namespace itk { ITKImageAccessorConverter::ITKImageAccessorConverter(DoseAccessorPointer accessor) { setDoseAccessor(accessor); } bool ITKImageAccessorConverter::process() { //Transfer GeometricInfo to ITK Properties core::GeometricInfo geoInfo = _doseAccessor->getGeometricInfo(); ITKImageType::RegionType region; ITKImageType::IndexType start; for (unsigned int i = 0; i < 3; ++i) { start[i] = 0; } ITKImageType::SizeType size; size[0] = geoInfo.getNumColumns(); size[1] = geoInfo.getNumRows(); size[2] = geoInfo.getNumSlices(); ITKImageType::SpacingType spacing; for (unsigned int i = 0; i < 3; ++i) { spacing[i] = geoInfo.getSpacing()[i]; } ITKImageType::PointType origin; for (unsigned int i = 0; i < 3; ++i) { origin[i] = geoInfo.getImagePositionPatient()[i]; } ITKImageType::DirectionType direction; OrientationMatrix OM = geoInfo.getOrientationMatrix(); for (unsigned int col = 0; col < 3; ++col) { for (unsigned int row = 0; row < 3; ++row) { direction(col, row) = OM(col, row); } } //Create image, assign properties region.SetSize(size); region.SetIndex(start); _itkImage = ITKImageType::New(); _itkImage->SetRegions(region); _itkImage->SetSpacing(spacing); _itkImage->SetDirection(direction); _itkImage->SetOrigin(origin); _itkImage->Allocate(); ::itk::DoseAccessorImageFilter::Pointer accessorFilter = ::itk::DoseAccessorImageFilter::New(); accessorFilter->SetInput(_itkImage); accessorFilter->SetAccessor(_doseAccessor); accessorFilter->Update(); _itkImage = accessorFilter->GetOutput(); return true; } }//end namespace itk }//end namespace io }//end namespace rttb diff --git a/code/io/itk/rttbITKImageAccessorGenerator.cpp b/code/io/itk/rttbITKImageAccessorGenerator.cpp index 3012752..d895ba2 100644 --- a/code/io/itk/rttbITKImageAccessorGenerator.cpp +++ b/code/io/itk/rttbITKImageAccessorGenerator.cpp @@ -1,53 +1,50 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include - #include -#include #include "rttbITKImageAccessor.h" #include "rttbITKImageAccessorGenerator.h" -#include "rttbException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace io { namespace itk { ITKImageAccessorGenerator::ITKImageAccessorGenerator(const ITKImageType* aDoseImage) { if (aDoseImage == nullptr) { throw core::InvalidDoseException("doseImage is nullptr"); } _dosePtr = aDoseImage; } core::DoseAccessorGeneratorBase::DoseAccessorPointer ITKImageAccessorGenerator::generateDoseAccessor() { _doseAccessor = boost::make_shared(_dosePtr); return _doseAccessor; } }//end namespace itk }//end namespace io }//end namespace rttb diff --git a/code/io/itk/rttbITKImageFileAccessorGenerator.cpp b/code/io/itk/rttbITKImageFileAccessorGenerator.cpp index 03f4395..a3978c7 100644 --- a/code/io/itk/rttbITKImageFileAccessorGenerator.cpp +++ b/code/io/itk/rttbITKImageFileAccessorGenerator.cpp @@ -1,48 +1,46 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbITKImageFileAccessorGenerator.h" -#include "rttbITKImageAccessorConverter.h" #include "rttbITKIOHelper.h" namespace rttb { namespace io { namespace itk { ITKImageFileAccessorGenerator::~ITKImageFileAccessorGenerator() = default; ITKImageFileAccessorGenerator::ITKImageFileAccessorGenerator(const FileNameType& fileName, const bool& useDicom) { _fileName = fileName; _useDicom = useDicom; } rttb::core::DoseAccessorGeneratorBase::DoseAccessorPointer ITKImageFileAccessorGenerator::generateDoseAccessor() { _itkDoubleImage = readITKDoubleImage(_fileName, _useDicom); _doseAccessor = boost::make_shared(_itkDoubleImage.GetPointer()); return _doseAccessor; } }//end namespace itk }//end namespace io }//end namespace rttb diff --git a/code/io/itk/rttbITKImageFileMaskAccessorGenerator.cpp b/code/io/itk/rttbITKImageFileMaskAccessorGenerator.cpp index b1dfd6f..bc80659 100644 --- a/code/io/itk/rttbITKImageFileMaskAccessorGenerator.cpp +++ b/code/io/itk/rttbITKImageFileMaskAccessorGenerator.cpp @@ -1,50 +1,47 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbITKImageFileMaskAccessorGenerator.h" -#include "rttbInvalidDoseException.h" -#include "rttbInvalidParameterException.h" -#include "../itk/rttbITKIOHelper.h" +#include "rttbITKIOHelper.h" namespace rttb { namespace io { namespace itk { ITKImageFileMaskAccessorGenerator::~ITKImageFileMaskAccessorGenerator() = default; ITKImageFileMaskAccessorGenerator::ITKImageFileMaskAccessorGenerator(const FileNameType& fileName) { _itkMaskFileName = fileName; } rttb::core::MaskAccessorGeneratorBase::MaskAccessorPointer ITKImageFileMaskAccessorGenerator::generateMaskAccessor() { _itkDoubleImage = rttb::io::itk::readITKDoubleImage(_itkMaskFileName); _maskAccessor = boost::make_shared(_itkDoubleImage.GetPointer()); return _maskAccessor; } }//end namespace itk }//end namespace io }//end namespace rttb diff --git a/code/io/itk/rttbITKImageMaskAccessor.cpp b/code/io/itk/rttbITKImageMaskAccessor.cpp index 830af89..9e0a9b8 100644 --- a/code/io/itk/rttbITKImageMaskAccessor.cpp +++ b/code/io/itk/rttbITKImageMaskAccessor.cpp @@ -1,170 +1,168 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbITKImageMaskAccessor.h" -#include "rttbGeometricInfo.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace io { namespace itk { ITKImageMaskAccessor::ITKImageMaskAccessor(ITKMaskImageType::ConstPointer aMaskImage) : _mask(aMaskImage) { if (_mask.IsNull()) { throw core::InvalidDoseException("Mask image = 0!") ; } assembleGeometricInfo(); } ITKImageMaskAccessor::~ITKImageMaskAccessor() = default; bool ITKImageMaskAccessor::assembleGeometricInfo() { _geoInfo = boost::make_shared(); _geoInfo->setSpacing(SpacingVectorType3D(_mask->GetSpacing()[0], _mask->GetSpacing()[1], _mask->GetSpacing()[2])); if (_geoInfo->getSpacing()[0] == 0 || _geoInfo->getSpacing()[1] == 0 || _geoInfo->getSpacing()[2] == 0) { throw core::InvalidDoseException("Pixel spacing = 0!"); } _geoInfo->setImagePositionPatient(WorldCoordinate3D(_mask->GetOrigin()[0], _mask->GetOrigin()[1], _mask->GetOrigin()[2])); OrientationMatrix OM(0); for (unsigned int col = 0; col < 3; ++col) { for (unsigned int row = 0; row < 3; ++row) { OM(col, row) = _mask->GetDirection()(col, row); } } _geoInfo->setOrientationMatrix(OM); _geoInfo->setNumColumns(_mask->GetLargestPossibleRegion().GetSize()[0]); _geoInfo->setNumRows(_mask->GetLargestPossibleRegion().GetSize()[1]); _geoInfo->setNumSlices(_mask->GetLargestPossibleRegion().GetSize()[2]); if (_geoInfo->getNumColumns() == 0 || _geoInfo->getNumRows() == 0 || _geoInfo->getNumSlices() == 0) { throw core::InvalidDoseException("Empty mask!") ; } return true; } void ITKImageMaskAccessor::updateMask() { return; } ITKImageMaskAccessor::MaskVoxelListPointer ITKImageMaskAccessor::getRelevantVoxelVector() { updateMask(); _relevantVoxelVector = getRelevantVoxelVector(0); return _relevantVoxelVector; } ITKImageMaskAccessor::MaskVoxelListPointer ITKImageMaskAccessor::getRelevantVoxelVector( float lowerThreshold) { auto filteredVoxelVectorPointer = boost::make_shared(); updateMask(); unsigned int size = _geoInfo->getNumColumns() * _geoInfo->getNumRows() * _geoInfo->getNumSlices(); filteredVoxelVectorPointer->reserve(size); for (unsigned int gridIndex = 0 ; gridIndex < size; gridIndex++) { core::MaskVoxel currentVoxel = core::MaskVoxel(gridIndex); if (getMaskAt(gridIndex, currentVoxel)) { if (currentVoxel.getRelevantVolumeFraction() > lowerThreshold) { filteredVoxelVectorPointer->push_back(currentVoxel); } } } return filteredVoxelVectorPointer; } bool ITKImageMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const { VoxelGridIndex3D aVoxelGridIndex; if (_geoInfo->convert(aID, aVoxelGridIndex)) { return getMaskAt(aVoxelGridIndex, voxel); } else { return false; } } bool ITKImageMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const { voxel.setRelevantVolumeFraction(0); if (_geoInfo->validIndex(aIndex)) { const ITKMaskImageType::IndexType pixelIndex = {{aIndex[0], aIndex[1], aIndex[2]}}; double value = _mask->GetPixel(pixelIndex); VoxelGridID gridId; _geoInfo->convert(aIndex, gridId); if (value >= 0 && value <= 1) { voxel.setRelevantVolumeFraction(value); return true; } else { std::cerr << "The pixel value of the mask should be >=0 and <=1!" << std::endl; return false; } } else { return false; } } const core::GeometricInfo& ITKImageMaskAccessor::getGeometricInfo() const { return *_geoInfo; }; } } } diff --git a/code/io/itk/rttbITKImageMaskAccessorConverter.cpp b/code/io/itk/rttbITKImageMaskAccessorConverter.cpp index 7982d28..c39d8b5 100644 --- a/code/io/itk/rttbITKImageMaskAccessorConverter.cpp +++ b/code/io/itk/rttbITKImageMaskAccessorConverter.cpp @@ -1,85 +1,84 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include #include "rttbITKImageMaskAccessorConverter.h" -#include "rttbGeometricInfo.h" #include "itkMaskAccessorImageSource.h" namespace rttb { namespace io { namespace itk { ITKImageMaskAccessorConverter::ITKImageMaskAccessorConverter(MaskAccessorPointer accessor) { setMaskAccessor(accessor); } bool ITKImageMaskAccessorConverter::process() { //Transfer GeometricInfo to ITK Properties core::GeometricInfo geoInfo = _maskAccessor->getGeometricInfo(); ITKImageMaskAccessor::ITKMaskImageType::IndexType start = {{0, 0, 0}}; ITKImageMaskAccessor::ITKMaskImageType::SizeType size = {{geoInfo.getNumColumns(), geoInfo.getNumRows(), geoInfo.getNumSlices()}}; ITKImageMaskAccessor::ITKMaskImageType::SpacingType spacing; for (unsigned int i = 0; i < 3; ++i) { spacing[i] = geoInfo.getSpacing()[i]; } ITKImageMaskAccessor::ITKMaskImageType::PointType origin; for (unsigned int i = 0; i < 3; ++i) { origin[i] = geoInfo.getImagePositionPatient()[i]; } ITKImageMaskAccessor::ITKMaskImageType::DirectionType direction; OrientationMatrix OM = geoInfo.getOrientationMatrix(); for (unsigned int col = 0; col < 3; ++col) { for (unsigned int row = 0; row < 3; ++row) { direction(col, row) = OM(col, row); } } ::itk::MaskAccessorImageSource::Pointer imagesource = ::itk::MaskAccessorImageSource::New(); imagesource->SetSize(size); imagesource->SetSpacing(spacing); imagesource->SetDirection(direction); imagesource->SetOrigin(origin); imagesource->SetAccessor(_maskAccessor); imagesource->SetFailsOnInvalidIDs(failsOnInvalidIDs()); imagesource->SetInvalidMaskValue(_invalidDoseValue); imagesource->Update(); _itkImage = imagesource->GetOutput(); return true; } }//end namespace mask }//end namespace io }//end namespace rttb diff --git a/code/io/itk/rttbITKImageMaskAccessorGenerator.cpp b/code/io/itk/rttbITKImageMaskAccessorGenerator.cpp index 102605d..5566ea8 100644 --- a/code/io/itk/rttbITKImageMaskAccessorGenerator.cpp +++ b/code/io/itk/rttbITKImageMaskAccessorGenerator.cpp @@ -1,51 +1,49 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbITKImageMaskAccessorGenerator.h" -#include "rttbException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace io { namespace itk { ITKImageMaskAccessorGenerator::ITKImageMaskAccessorGenerator(const ITKImageMaskAccessor::ITKMaskImageType* aMaskImage) { if (aMaskImage == nullptr) { throw core::InvalidDoseException("MaskImage is nullptr"); } _maskPtr = aMaskImage; } core::MaskAccessorGeneratorBase::MaskAccessorPointer ITKImageMaskAccessorGenerator::generateMaskAccessor() { _maskAccessor = boost::make_shared(_maskPtr); return _maskAccessor; } }//end namespace mask }//end namespace io }//end namespace rttb diff --git a/code/io/models/rttbModelXMLWriter.cpp b/code/io/models/rttbModelXMLWriter.cpp index b187d7c..be962c0 100644 --- a/code/io/models/rttbModelXMLWriter.cpp +++ b/code/io/models/rttbModelXMLWriter.cpp @@ -1,142 +1,135 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbModelXMLWriter.h" +#include "rttbDVHXMLFileWriter.h" +#include "rttbInvalidParameterException.h" -#include - -/*boost includes*/ -#include #include -#include - -#include "rttbInvalidParameterException.h" -#include "rttbDVHXMLFileWriter.h" -#include "rttbDoseStatisticsXMLWriter.h" namespace rttb { namespace io { namespace models { ModelXMLWriter::ModelXMLWriter(const std::string& filename, rttb::models::BioModel::Pointer model, bool printDVH) : _filename(filename), _model(model), _printDVH(printDVH) { } void ModelXMLWriter::setFilename(FileNameString filename) { _filename = filename; } FileNameString ModelXMLWriter::getFilename() const { return _filename; } void ModelXMLWriter::setModel(rttb::models::BioModel::Pointer model) { _model = model; } rttb::models::BioModel::Pointer ModelXMLWriter::getModel() const { return _model; } void ModelXMLWriter::setPrintDVH(bool printDVH) { _printDVH = printDVH; } bool ModelXMLWriter::getPrintDVH() const { return _printDVH; } void ModelXMLWriter::writeModel() { std::setprecision(2); boost::property_tree::ptree pt; std::string xmlattrNameTag = ".name"; std::string modelTag = "BioModel"; std::string propertyTag = "property"; std::string configTag = "config"; std::string resultsTag = "results"; std::string valueTag = "value"; auto dvh = _model->getDVH(); boost::property_tree::ptree propertynode; boost::property_tree::ptree confignode; confignode.put("BioModelType", _model->getModelType()); confignode.put("StructureID", dvh->getStructureID()); confignode.put("DoseID", dvh->getDoseID()); if (_printDVH) { FileNameString filename = "dvhfor" + _model->getModelType() + ".xml"; DVHType typeDiff = { DVHType::Differential }; io::other::DVHXMLFileWriter dvhWriter(filename, typeDiff); dvhWriter.writeDVH(dvh); confignode.put("DVHReference", filename); } pt.add_child(modelTag + "." + configTag, confignode); propertynode.put("", _model->getValue()); propertynode.put(xmlattrNameTag, valueTag); pt.add_child(modelTag + "." + resultsTag + "." + propertyTag, propertynode); std::map parameterMap = _model->getParameterMap(); for (std::map::const_iterator it = parameterMap.begin(); it != parameterMap.end(); ++it) { double value = it->second; if (value - static_cast(std::round(value)) < errorConstant) { propertynode.put("", static_cast(value)); } else { propertynode.put("", static_cast(value)); } propertynode.put(xmlattrNameTag, it->first); pt.add_child(modelTag + "." + resultsTag + "." + propertyTag, propertynode); } try { boost::property_tree::xml_parser::xml_writer_settings settings = boost::property_tree::xml_writer_make_settings('\t', 1); boost::property_tree::xml_parser::write_xml(_filename, pt, std::locale(), settings); } catch (const boost::property_tree::xml_parser_error& e) { std::cout << e.what(); throw core::InvalidParameterException("Write xml failed: xml_parser_error!"); } } } } } diff --git a/code/io/other/rttbDVHXMLFileReader.cpp b/code/io/other/rttbDVHXMLFileReader.cpp index 6616a50..42b9ac8 100644 --- a/code/io/other/rttbDVHXMLFileReader.cpp +++ b/code/io/other/rttbDVHXMLFileReader.cpp @@ -1,148 +1,147 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include #include -#include #include #include #include #include "rttbDVHXMLFileReader.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace io { namespace other { DVHXMLFileReader::DVHXMLFileReader(FileNameString aFileName) { _fileName = aFileName; _resetFile = true; } DVHXMLFileReader::~DVHXMLFileReader() = default; void DVHXMLFileReader::resetFileName(FileNameString aFileName) { _fileName = aFileName; _resetFile = true; } void DVHXMLFileReader::createDVH() { boost::property_tree::ptree pt; // Load the XML file into the property tree. If reading fails // (cannot open file, parse error), an exception is thrown. try { read_xml(_fileName, pt); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw rttb::core::InvalidParameterException("DVH file name invalid: could not open the xml file!"); } std::string dvhType; std::deque dataDifferential, dataCumulative; DoseTypeGy deltaD = 0; DoseVoxelVolumeType deltaV = 0; IDType strID; IDType doseID; bool normalized; dvhType = pt.get("dvh.type"); deltaD = pt.get("dvh.deltaD"); deltaV = pt.get("dvh.deltaV"); strID = pt.get("dvh.structureID"); doseID = pt.get("dvh.doseID"); boost::optional dvhNormalized = pt.get_optional("dvh.normalized"); if (dvhNormalized) { normalized = dvhNormalized.get(); } else { normalized = false; } if (dvhType != "DIFFERENTIAL" && dvhType != "CUMULATIVE") { throw core::InvalidParameterException("DVH Type invalid! Only: DIFFERENTIAL/CUMULATIVE!"); } BOOST_FOREACH(boost::property_tree::ptree::value_type & v, pt.get_child("dvh.data")) { if (dvhType == "DIFFERENTIAL") { dataDifferential.push_back(boost::lexical_cast(v.second.data()) / (normalized ? deltaV : 1)); } else if (dvhType == "CUMULATIVE") { dataCumulative.push_back(boost::lexical_cast(v.second.data()) / (normalized ? deltaV : 1)); } } unsigned int numberOfBins = static_cast(std::max(dataDifferential.size(), dataCumulative.size())); if (dvhType == "CUMULATIVE") //dataDifferential should be calculated { DoseCalcType differentialDVHi = 0; std::deque::iterator it; for (it = dataCumulative.begin(); it != dataCumulative.end(); ++it) { if (dataDifferential.size() == numberOfBins - 1) { differentialDVHi = *it; } else { differentialDVHi = *it - *(it + 1); } dataDifferential.push_back(differentialDVHi); } } _dvh = boost::make_shared(dataDifferential, deltaD, deltaV, strID, doseID); _resetFile = false; } core::DVH::Pointer DVHXMLFileReader::generateDVH() { if (_resetFile) { this->createDVH(); } return _dvh; } }//end namespace other }//end namespace io }//end namespace rttb diff --git a/code/io/other/rttbDVHXMLFileWriter.cpp b/code/io/other/rttbDVHXMLFileWriter.cpp index edede47..908741a 100644 --- a/code/io/other/rttbDVHXMLFileWriter.cpp +++ b/code/io/other/rttbDVHXMLFileWriter.cpp @@ -1,139 +1,137 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDVHXMLFileWriter.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" -#include "rttbException.h" -#include #include namespace rttb { namespace io { namespace other { DVHXMLFileWriter::DVHXMLFileWriter(FileNameString aFileName, DVHType aDVHType) { this->setFileName(aFileName); this->setDVHType(aDVHType); } void DVHXMLFileWriter::setDVHType(DVHType aDVHType) { _dvhType = aDVHType; } FileNameString DVHXMLFileWriter::getFileName() const { return _fileName; } void DVHXMLFileWriter::setFileName(FileNameString aFileName) { _fileName = aFileName; } DVHType DVHXMLFileWriter::getDVHType() const { return _dvhType; } void DVHXMLFileWriter::writeDVH(core::DVH::Pointer aDvh, bool normalized) { if (aDvh == nullptr) { throw core::NullPointerException("aDvh must not be nullptr! "); } using boost::property_tree::ptree; ptree pt; DataDifferentialType data; std::map normalizedData; if (_dvhType.Type == DVHType::Differential) { pt.put("dvh.type", "DIFFERENTIAL"); if (normalized) { normalizedData = aDvh->getNormalizedDVH({ DVHType::Cumulative }); } else { data = aDvh->getDataDifferential(); } } else if (_dvhType.Type == DVHType::Cumulative) { pt.put("dvh.type", "CUMULATIVE"); if (normalized) { normalizedData = aDvh->getNormalizedDVH(); } else { data = aDvh->getDataCumulative(); } } else { throw core::InvalidParameterException("DVH Type not acceptable: Only: DIFFERENTIAL/CUMULATIVE!"); } pt.put("dvh.deltaD", aDvh->getDeltaD()); pt.put("dvh.deltaV", aDvh->getDeltaV()); pt.put("dvh.structureID", aDvh->getStructureID()); pt.put("dvh.doseID", aDvh->getDoseID()); pt.put("dvh.normalized", normalized); if (normalized) { for (auto elem : normalizedData) { boost::property_tree::ptree node; node.put("", elem.second); node.put(".doseGy", elem.first); pt.add_child("dvh.data.volumeInCcm", node); } } else { for (size_t i = 0; i < data.size(); i++) { boost::property_tree::ptree node; node.put("", data[i]); node.put(".dosebin", i); pt.add_child("dvh.data.volume", node); } } try { boost::property_tree::xml_parser::xml_writer_settings settings = boost::property_tree::xml_writer_make_settings ('\t', 1); boost::property_tree::xml_parser::write_xml(_fileName, pt, std::locale(), settings); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml failed: xml_parser_error!"); } } } } } diff --git a/code/io/other/rttbDoseStatisticsXMLReader.cpp b/code/io/other/rttbDoseStatisticsXMLReader.cpp index b4530f5..db16dd4 100644 --- a/code/io/other/rttbDoseStatisticsXMLReader.cpp +++ b/code/io/other/rttbDoseStatisticsXMLReader.cpp @@ -1,214 +1,209 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ -#include -#include -#include -#include -#include -#include - #include "rttbDoseStatisticsXMLReader.h" - #include "rttbInvalidParameterException.h" -#include "rttbVolumeToDoseMeasureCollection.h" +#include +#include +#include +#include namespace rttb { namespace io { namespace other { DoseStatisticsXMLReader::DoseStatisticsXMLReader(const std::string& filename) : _filename(filename), _newFile(true) { } DoseStatisticsXMLReader::~DoseStatisticsXMLReader() = default; void DoseStatisticsXMLReader::setFilename(const std::string& filename) { _filename = filename; _newFile = true; } algorithms::DoseStatistics::Pointer DoseStatisticsXMLReader::generateDoseStatistic() { if (_newFile) { this->createDoseStatistic(); } return _doseStatistic; } void DoseStatisticsXMLReader::createDoseStatistic() { boost::property_tree::ptree pt; // Load the XML file into the property tree. If reading fails // (cannot open file, parse error), an exception is thrown. try { read_xml(_filename, pt); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw rttb::core::InvalidParameterException("DoseStatistics file name invalid: could not open the xml file!"); } // data to fill the parameter for the DoseStatistics std::string name; std::string datum; unsigned int x; std::pair voxelid; std::vector < std::pair> vec; //initialize all parameters for the DoseStatistics double minimum=-1; double maximum=-1; double numVoxels=-1; double volume=-1; double referenceDose = -1; double mean=-1; double stdDeviation=-1; boost::shared_ptr > > minimumVoxelPositions = nullptr; boost::shared_ptr > > maximumVoxelPositions = nullptr; rttb::algorithms::VolumeToDoseMeasureCollection Dx(rttb::algorithms::VolumeToDoseMeasureCollection::Dx); rttb::algorithms::DoseToVolumeMeasureCollection Vx(rttb::algorithms::DoseToVolumeMeasureCollection::Vx); rttb::algorithms::VolumeToDoseMeasureCollection MOHx(rttb::algorithms::VolumeToDoseMeasureCollection::MOHx); rttb::algorithms::VolumeToDoseMeasureCollection MOCx(rttb::algorithms::VolumeToDoseMeasureCollection::MOCx); rttb::algorithms::VolumeToDoseMeasureCollection MaxOHx(rttb::algorithms::VolumeToDoseMeasureCollection::MaxOHx); rttb::algorithms::VolumeToDoseMeasureCollection MinOCx(rttb::algorithms::VolumeToDoseMeasureCollection::MinOCx); BOOST_FOREACH(boost::property_tree::ptree::value_type & data, pt.get_child("statistics.results")) { datum = data.second.data(); BOOST_FOREACH(boost::property_tree::ptree::value_type & middle, data.second) { BOOST_FOREACH(boost::property_tree::ptree::value_type & innernode, middle.second) { std::string mia = innernode.first; if (innernode.first == "name") { name = innernode.second.data(); } else if (innernode.first == "voxelGridID") { boost::replace_all(datum, "\r\n", ""); boost::replace_all(datum, "\n", ""); boost::trim(datum); voxelid.first = boost::lexical_cast(datum); voxelid.second = boost::lexical_cast(innernode.second.data()); vec.push_back(voxelid); } else if (innernode.first == "x") { x = boost::lexical_cast(innernode.second.data()); } } } // fill with the extracted data if (name == "numberOfVoxels") { numVoxels = boost::lexical_cast(datum); } else if (name == "volume") { volume = boost::lexical_cast(datum); Dx.setVolume(volume); MOHx.setVolume(volume); MOCx.setVolume(volume); MaxOHx.setVolume(volume); } else if (name == "referenceDose") { referenceDose = boost::lexical_cast(datum); Vx.setReferenceDose(referenceDose); } else if (name == "mean") { mean = boost::lexical_cast(datum); } else if (name == "standardDeviation") { stdDeviation = boost::lexical_cast(datum); } else if (name == "minimum") { minimum = boost::lexical_cast(datum); if (!vec.empty()) { minimumVoxelPositions = boost::make_shared>>(vec); vec.clear(); } } else if (name == "maximum") { maximum = boost::lexical_cast(datum); if (!vec.empty()) { maximumVoxelPositions = boost::make_shared>>(vec); vec.clear(); } } else if (name == "Dx") { Dx.insertValue(static_cast(x)*volume / 100, boost::lexical_cast(datum)); } else if (name == "Vx") { Vx.insertValue(static_cast(x)*referenceDose / 100, boost::lexical_cast(datum)); } else if (name == "MOHx") { MOHx.insertValue(static_cast(x)*volume / 100, boost::lexical_cast(datum)); } else if (name == "MOCx") { MOCx.insertValue(static_cast(x)*volume / 100, boost::lexical_cast(datum)); } else if (name == "MaxOHx") { MaxOHx.insertValue(static_cast(x)*volume / 100, boost::lexical_cast(datum)); } else if (name == "MinOCx") { MinOCx.insertValue(static_cast(x)*volume / 100, boost::lexical_cast(datum)); } } // make DoseStatistcs _doseStatistic = boost::make_shared( minimum, maximum, mean, stdDeviation, numVoxels, volume, minimumVoxelPositions, maximumVoxelPositions , Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx, referenceDose); } }//end namespace other }//end namespace io }//end namespace rttb diff --git a/code/io/other/rttbDoseStatisticsXMLWriter.cpp b/code/io/other/rttbDoseStatisticsXMLWriter.cpp index 5ed7159..4ef59e7 100644 --- a/code/io/other/rttbDoseStatisticsXMLWriter.cpp +++ b/code/io/other/rttbDoseStatisticsXMLWriter.cpp @@ -1,296 +1,290 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseStatisticsXMLWriter.h" - -#include -#include - -#include - #include "rttbInvalidParameterException.h" #include "rttbNullPointerException.h" -#include "rttbVolumeToDoseMeasureCollection.h" -#include "rttbDoseToVolumeMeasureCollection.h" + +#include namespace rttb { namespace io { namespace other { static const std::string xmlattrXTag = ".x"; static const std::string xmlattrNameTag = ".name"; static const std::string statisticsTag = "statistics.results"; static const std::string propertyTag = "property"; static const std::string columnSeparator = "@"; boost::property_tree::ptree DoseStatisticsXMLWriter::writeDoseStatistics(DoseStatisticsPtr aDoseStatistics) { using boost::property_tree::ptree; ptree pt; if (aDoseStatistics == nullptr){ throw core::NullPointerException("dose statistics is nullptr!"); } ptree numberOfVoxelsNode = createNodeWithNameAttribute(aDoseStatistics->getNumberOfVoxels(), "numberOfVoxels"); pt.add_child(statisticsTag + "." + propertyTag, numberOfVoxelsNode); ptree volumeNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getVolume()), "volume"); pt.add_child(statisticsTag + "." + propertyTag, volumeNode); ptree referenceNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getReferenceDose()), "referenceDose"); pt.add_child(statisticsTag + "." + propertyTag, referenceNode); ptree minimumNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getMinimum()), "minimum"); auto minimumPositions = aDoseStatistics->getMinimumVoxelPositions(); auto pairItMin = minimumPositions->begin(); int count = 0; for (; pairItMin != minimumPositions->end() && count < 100; ++pairItMin) //output max. 100 minimum { VoxelGridID voxelID = pairItMin->second; ptree voxelMinPositions; voxelMinPositions.add("voxelGridID", voxelID); minimumNode.add_child("voxel", voxelMinPositions); count++; } pt.add_child(statisticsTag + "." + propertyTag, minimumNode); ptree maximumNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getMaximum()), "maximum"); auto maximumPositions = aDoseStatistics->getMaximumVoxelPositions(); auto pairItMax = maximumPositions->begin(); count = 0; for (; pairItMax != maximumPositions->end() && count < 100; ++pairItMax) //output max. 100 maximum { VoxelGridID voxelID = pairItMax->second; ptree voxelMaxPositions; voxelMaxPositions.add("voxelGridID", voxelID); maximumNode.add_child("voxel", voxelMaxPositions); count++; } pt.add_child(statisticsTag + "." + propertyTag, maximumNode); ptree meanNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getMean()), "mean"); pt.add_child(statisticsTag + "." + propertyTag, meanNode); ptree sdNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getStdDeviation()), "standardDeviation"); pt.add_child(statisticsTag + "." + propertyTag, sdNode); ptree varianceNode = createNodeWithNameAttribute(static_cast(aDoseStatistics->getVariance()), "variance"); pt.add_child(statisticsTag + "." + propertyTag, varianceNode); double absoluteVolume = aDoseStatistics->getVolume(); double referenceDose = aDoseStatistics->getReferenceDose(); rttb::algorithms::DoseToVolumeMeasureCollection::DoseToVolumeFunctionType AllVx = aDoseStatistics->getVx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllDx = aDoseStatistics->getDx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getMOHx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getMOCx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getMaxOHx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getMinOCx().getAllValues(); rttb::algorithms::DoseToVolumeMeasureCollection::DoseToVolumeFunctionType::iterator vxIt; rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType::iterator it; for (it = AllDx.begin(); it != AllDx.end(); ++it) { ptree DxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "Dx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, DxNode); } for (vxIt = AllVx.begin(); vxIt != AllVx.end(); ++vxIt) { ptree VxNode = createNodeWithNameAndXAttribute(static_cast(vxIt->second), "Vx", std::lround(convertToPercent(vxIt->first, referenceDose))); pt.add_child(statisticsTag + "." + propertyTag, VxNode); } for (it = AllMOHx.begin(); it != AllMOHx.end(); ++it) { ptree mohxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "MOHx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, mohxNode); } for (it = AllMOCx.begin(); it != AllMOCx.end(); ++it) { ptree mocxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "MOCx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, mocxNode); } for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); ++it) { ptree maxOhxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "MaxOHx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, maxOhxNode); } for (it = AllMinOCx.begin(); it != AllMinOCx.end(); ++it) { ptree minOCxNode = createNodeWithNameAndXAttribute(static_cast(it->second), "MinOCx", std::lround(convertToPercent(it->first, absoluteVolume))); pt.add_child(statisticsTag + "." + propertyTag, minOCxNode); } return pt; } void DoseStatisticsXMLWriter::writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName) { boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics); try { boost::property_tree::xml_parser::write_xml(aFileName, pt, std::locale(), boost::property_tree::xml_writer_make_settings('\t', 1)); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml failed: xml_parser_error!"); } } XMLString DoseStatisticsXMLWriter::writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics) { boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics); std::stringstream sstr; try { boost::property_tree::xml_parser::write_xml(sstr, pt, boost::property_tree::xml_writer_make_settings('\t', 1)); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml to string failed: xml_parser_error!"); } return sstr.str(); } StatisticsString DoseStatisticsXMLWriter::writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics) { if (aDoseStatistics == nullptr){ throw core::NullPointerException("dose statistics is nullptr!"); } std::stringstream sstr; sstr << static_cast(aDoseStatistics->getVolume() * 1000) << columnSeparator; // cm3 to mm3 sstr << static_cast(aDoseStatistics->getMaximum()) << columnSeparator; sstr << static_cast(aDoseStatistics->getMinimum()) << columnSeparator; sstr << static_cast(aDoseStatistics->getMean()) << columnSeparator; sstr << static_cast(aDoseStatistics->getStdDeviation()) << columnSeparator; sstr << static_cast(aDoseStatistics->getVariance()) << columnSeparator; rttb::algorithms::DoseToVolumeMeasureCollection::DoseToVolumeFunctionType AllVx = aDoseStatistics->getVx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllDx = aDoseStatistics->getDx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getMOHx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getMOCx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getMaxOHx().getAllValues(); rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getMinOCx().getAllValues(); rttb::algorithms::DoseToVolumeMeasureCollection::DoseToVolumeFunctionType::iterator vxIt; rttb::algorithms::VolumeToDoseMeasureCollection::VolumeToDoseFunctionType::iterator it; for (it = AllDx.begin(); it != AllDx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } for (vxIt = AllVx.begin(); vxIt != AllVx.end(); ++vxIt) { // *1000 because of conversion cm3 to mm3 sstr << static_cast((*vxIt).second * 1000) << columnSeparator; } for (it = AllMOHx.begin(); it != AllMOHx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } for (it = AllMOCx.begin(); it != AllMOCx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } for (it = AllMinOCx.begin(); it != AllMinOCx.end(); ++it) { sstr << static_cast((*it).second) << columnSeparator; } return sstr.str(); } boost::property_tree::ptree DoseStatisticsXMLWriter::createNodeWithNameAttribute(DoseTypeGy doseValue, const std::string& attributeName) { boost::property_tree::ptree node; node.put("", doseValue); node.put(xmlattrNameTag, attributeName); return node; } boost::property_tree::ptree DoseStatisticsXMLWriter::createNodeWithNameAndXAttribute(DoseTypeGy doseValue, const std::string& attributeName, int xValue) { boost::property_tree::ptree node; node.put("", doseValue); node.put(xmlattrNameTag, attributeName); node.put(xmlattrXTag, xValue); return node; } double DoseStatisticsXMLWriter::convertToPercent(double value, double maximum) { return (value / maximum) * 100; } }//end namespace other }//end namespace io }//end namespace rttb \ No newline at end of file diff --git a/code/io/utils/rttbDoseLoader.cpp b/code/io/utils/rttbDoseLoader.cpp index dfd892e..a4ffbd6 100644 --- a/code/io/utils/rttbDoseLoader.cpp +++ b/code/io/utils/rttbDoseLoader.cpp @@ -1,87 +1,85 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseLoader.h" - #include "rttbExceptionMacros.h" - #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomHelaxFileDoseAccessorGenerator.h" #include "rttbITKImageFileAccessorGenerator.h" namespace rttb { namespace io { namespace utils { rttb::core::DoseAccessorInterface::Pointer loadDicomDose(const std::string& fileName) { rttb::io::dicom::DicomFileDoseAccessorGenerator generator(fileName); return generator.generateDoseAccessor(); } rttb::core::DoseAccessorInterface::Pointer loadHelaxDose(const std::string& path) { rttb::io::helax::DicomHelaxFileDoseAccessorGenerator generator(path); return generator.generateDoseAccessor(); } rttb::core::DoseAccessorInterface::Pointer loadITKDose(const std::string& fileName) { rttb::io::itk::ITKImageFileAccessorGenerator generator(fileName); return generator.generateDoseAccessor(); } rttb::core::DoseAccessorInterface::Pointer loadITKDicomDose(const std::string& fileName) { rttb::io::itk::ITKImageFileAccessorGenerator generator(fileName, true); return generator.generateDoseAccessor(); } rttb::core::DoseAccessorInterface::Pointer loadDose(const std::string& fileName, const std::string& loadStyle) { rttb::core::DoseAccessorInterface::Pointer result; if (loadStyle == "" || loadStyle == "dicom") { result = loadDicomDose(fileName); } else if (loadStyle == "helax") { result = loadHelaxDose(fileName); } else if (loadStyle == "itk") { result = loadITKDose(fileName); } else if (loadStyle == "itkDicom") { result = loadITKDicomDose(fileName); } else { rttbDefaultExceptionStaticMacro(<< "Unknown io style selected. Cannot load data. Selected style: " << loadStyle); } return result; } } } } diff --git a/code/masks/rttbBoostMask.cpp b/code/masks/rttbBoostMask.cpp index 63c13ae..e507fd1 100644 --- a/code/masks/rttbBoostMask.cpp +++ b/code/masks/rttbBoostMask.cpp @@ -1,558 +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. // //------------------------------------------------------------------------ #include -#include -#include #include #include #include "rttbBoostMask.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbBoostMaskVoxelizationThread.h" namespace rttb { namespace masks { namespace boost { BoostMask::BoostMask(core::GeometricInfo::Pointer aDoseGeoInfo, core::Structure::Pointer aStructure, bool strict, unsigned int numberOfThreads) : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _numberOfThreads(numberOfThreads), _voxelizationThickness(0.0), _voxelInStructure(::boost::make_shared()) { _isUpToDate = false; if (_geometricInfo == nullptr) { throw rttb::core::NullPointerException("Error: Geometric info is nullptr!"); } else if (_structure == nullptr) { throw rttb::core::NullPointerException("Error: Structure is nullptr!"); } if (_numberOfThreads == 0) { _numberOfThreads = ::boost::thread::hardware_concurrency(); if (_numberOfThreads == 0) { throw rttb::core::InvalidParameterException("Error: detection of the number of hardware threads is not possible. Please specify number of threads for voxelization explicitly as parameter in BoostMask."); } } } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { if (!_isUpToDate) { calcMask(); } return _voxelInStructure; } void BoostMask::calcMask() { preprocessing(); voxelization(); generateMaskVoxelList(); _isUpToDate = true; } void BoostMask::preprocessing() { rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); //Convert world coordinate polygons to the polygons with geometry coordinate rttb::PolygonSequenceType geometryCoordinatePolygonVector; rttb::PolygonSequenceType::iterator it; rttb::DoubleVoxelGridIndex3D globalMaxGridIndex(std::numeric_limits::min(), std::numeric_limits::min(), std::numeric_limits::min()); rttb::DoubleVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(), _geometricInfo->getNumRows(), 0); for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) { PolygonType rttbPolygon = *it; PolygonType geometryCoordinatePolygon; //1. convert polygon to geometry coordinate polygons //2. calculate global min/max //3. check if polygon is planar if (!preprocessingPolygon(rttbPolygon, geometryCoordinatePolygon, globalMinGridIndex, globalMaxGridIndex, errorConstant)) { throw rttb::core::Exception("TiltedMaskPlaneException"); } geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon); } rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) ), GridIndexType(globalMinGridIndex(1) ), GridIndexType(globalMinGridIndex(2) )); rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) ), GridIndexType(globalMaxGridIndex(1) ), GridIndexType(globalMaxGridIndex(2) )); _globalBoundingBox.push_back(minIndex); _globalBoundingBox.push_back(maxIndex); //convert rttb polygon sequence to a map of z index and a vector of boost ring 2d (without holes) _ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); } void BoostMask::voxelization() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } BoostRingMap::iterator itMap; size_t mapSizeInAThread = _ringMap.size() / _numberOfThreads; unsigned int count = 0; unsigned int countThread = 0; BoostPolygonMap polygonMap; std::vector polygonMapVector; //check donut and convert to a map of z index and a vector of boost polygon 2d (with or without holes) for (itMap = _ringMap.begin(); itMap != _ringMap.end(); ++itMap) { //the vector of all boost 2d polygons with the same z grid index(donut polygon is accepted). BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); if (count == mapSizeInAThread && countThread < (_numberOfThreads - 1)) { polygonMapVector.push_back(polygonMap); polygonMap.clear(); count = 0; countThread++; } polygonMap.insert(std::pair((*itMap).first, polygonVector)); count++; } _voxelizationMap = ::boost::make_shared >(); polygonMapVector.push_back(polygonMap); //insert the last one //generate voxelization map, multi-threading ::boost::thread_group threads; auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); for (const auto & i : polygonMapVector) { BoostMaskVoxelizationThread t(i, _globalBoundingBox, _voxelizationMap, aMutex, _strict); threads.create_thread(t); } threads.join_all(); } void BoostMask::generateMaskVoxelList() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } //check homogeneous of the voxelization plane (the contours plane) if (!calcVoxelizationThickness(_voxelizationThickness)) { throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneous!"); } ::boost::thread_group threads; auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); unsigned int sliceNumberInAThread = _geometricInfo->getNumSlices() / _numberOfThreads; //generate mask voxel list, multi-threading for (unsigned int i = 0; i < _numberOfThreads; ++i) { unsigned int beginSlice = i * sliceNumberInAThread; unsigned int endSlice; if (i < _numberOfThreads - 1) { endSlice = (i + 1) * sliceNumberInAThread; } else { endSlice = _geometricInfo->getNumSlices(); } BoostMaskGenerateMaskVoxelListThread t(_globalBoundingBox, _geometricInfo, _voxelizationMap, _voxelizationThickness, beginSlice, endSlice, _voxelInStructure, _strict, aMutex); threads.create_thread(t); } threads.join_all(); } bool BoostMask::preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const { double minZ = _geometricInfo->getNumSlices(); double maxZ = 0.0; for (auto worldCoordinatePoint : aRTTBPolygon) { //convert to geometry coordinate polygon rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint; _geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, geometryCoordinatePoint); geometryCoordinatePolygon.push_back(geometryCoordinatePoint); //calculate the current global min/max //min and max for x if (geometryCoordinatePoint(0) < minimum(0)) { minimum(0) = geometryCoordinatePoint(0); } if (geometryCoordinatePoint(0) > maximum(0)) { maximum(0) = geometryCoordinatePoint(0); } //min and max for y if (geometryCoordinatePoint(1) < minimum(1)) { minimum(1) = geometryCoordinatePoint(1); } if (geometryCoordinatePoint(1) > maximum(1)) { maximum(1) = geometryCoordinatePoint(1); } //min and max for z if (geometryCoordinatePoint(2) < minimum(2)) { minimum(2) = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maximum(2)) { maximum(2) = geometryCoordinatePoint(2); } //check planar if (geometryCoordinatePoint(2) < minZ) { minZ = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maxZ) { maxZ = geometryCoordinatePoint(2); } } return (std::abs(maxZ - minZ) <= aErrorConstant); } BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const { BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); if (i == 0) { firstPoint = boostPoint; } ::boost::geometry::append(polygon2D, boostPoint); } ::boost::geometry::append(polygon2D, firstPoint); return polygon2D; } BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( const rttb::PolygonSequenceType& aRTTBPolygonVector) const { rttb::PolygonSequenceType::const_iterator it; BoostMask::BoostRingMap aRingMap; for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it) { rttb::PolygonType rttbPolygon = *it; double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon bool isFirstZ = true; if (!aRingMap.empty()) { auto findIt = findNearestKey(aRingMap, zIndex, errorConstant); //if the z index is found (same slice), add the polygon to vector if (findIt != aRingMap.end()) { //BoostRingVector ringVector = ; (*findIt).second.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); isFirstZ = false; } } //if it is the first z index in the map, insert vector with the polygon if (isFirstZ) { BoostRingVector ringVector; ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); aRingMap.insert(std::pair(zIndex, ringVector)); } } return aRingMap; } BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap& aBoostRingMap, double aIndex, double aErrorConstant) const { auto find = aBoostRingMap.find(aIndex); //if find a key equivalent to aIndex, found if (find != aBoostRingMap.end()) { return find; } else { auto lowerBound = aBoostRingMap.lower_bound(aIndex); //if all keys go before aIndex, check the last key if (lowerBound == aBoostRingMap.end()) { lowerBound = --aBoostRingMap.end(); } //if the lower bound very close to aIndex, found if (std::abs((*lowerBound).first - aIndex) <= aErrorConstant) { return lowerBound; } else { //if the lower bound is the beginning, not found if (lowerBound == aBoostRingMap.begin()) { return aBoostRingMap.end(); } else { auto lowerBound1 = --lowerBound;//the key before the lower bound //if the key before the lower bound very close to a Index, found if (std::abs((*lowerBound1).first - aIndex) <= aErrorConstant) { return lowerBound1; } //else, not found else { return aBoostRingMap.end(); } } } } } BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector) const { //check donut BoostMask::BoostRingVector::const_iterator it1; BoostMask::BoostRingVector::const_iterator it2; BoostMask::BoostPolygonVector boostPolygonVector; std::vector donutIndexVector;//store the outer and inner ring index BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon //Get donut index and donut polygon unsigned int index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i : donutIndexVector) { if (i == index1) { it1IsDonut = true; break; } } //if not jet, check now if (!it1IsDonut) { bool it2IsDonut = false; unsigned int index2 = 0; for (it2 = aRingVector.begin(); it2 != aRingVector.end(); ++it2, index2++) { if (it2 != it1) { BoostMask::BoostPolygon2D polygon2D; if (::boost::geometry::within(*it1, *it2)) { ::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring it2IsDonut = true; } //if donut else if (::boost::geometry::within(*it2, *it1)) { ::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring it2IsDonut = true; } if (it2IsDonut) { donutIndexVector.push_back(index1); donutIndexVector.push_back(index2); donutVector.push_back(polygon2D);//store donut polygon break;//Only store the first donut! } } } } } //Store no donut polygon to boostPolygonVector index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++) { bool it1IsDonut = false; //check if the ring is the outer or inner of a donut for (unsigned int i : donutIndexVector) { if (i == index1) { it1IsDonut = true; break; } } if (!it1IsDonut) { BoostMask::BoostPolygon2D polygon2D; ::boost::geometry::append(polygon2D, *it1); boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut } } //Append donut polygon to boostPolygonVector BoostMask::BoostPolygonVector::iterator itDonut; for (itDonut = donutVector.begin(); itDonut != donutVector.end(); ++itDonut) { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } bool BoostMask::calcVoxelizationThickness(double& aThickness) const { if (_voxelizationMap->size() <= 1) { aThickness = 1; return true; } double thickness = 0; auto it = _voxelizationMap->cbegin(); auto it2 = ++_voxelizationMap->cbegin(); for (; it != _voxelizationMap->cend() && it2 != _voxelizationMap->cend(); ++it, ++it2) { if (thickness == 0) { thickness = it2->first - it->first; } else { double curThickness = it2->first - it->first; //if not homogeneous (leave out double imprecisions), return false if (std::abs(thickness-curThickness)>errorConstant) { //return false; std::cout << "Two polygons are far from each other?" << std::endl; } } } if (thickness != 0) { aThickness = thickness; } else { aThickness = 1; } return true; } } } } diff --git a/code/masks/rttbBoostMaskAccessor.cpp b/code/masks/rttbBoostMaskAccessor.cpp index 231f1a7..2061fce 100644 --- a/code/masks/rttbBoostMaskAccessor.cpp +++ b/code/masks/rttbBoostMaskAccessor.cpp @@ -1,155 +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. // //------------------------------------------------------------------------ #include "rttbBoostMaskAccessor.h" #include "rttbBoostMask.h" #include -#include #include #include namespace rttb { namespace masks { namespace boost { BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, bool strict) : _spStructure(aStructurePointer), _geoInfo(aGeometricInfo), _strict(strict) { _spRelevantVoxelVector = MaskVoxelListPointer(); //generate new structure set uid ::boost::uuids::uuid id; ::boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _maskUID = "BoostMask_" + ss.str(); } BoostMaskAccessor::~BoostMaskAccessor() = default; void BoostMaskAccessor::updateMask() { MaskVoxelList newRelevantVoxelVector; if (_spRelevantVoxelVector) { return; // already calculated } BoostMask mask(::boost::make_shared(_geoInfo), _spStructure, _strict); _spRelevantVoxelVector = mask.getRelevantVoxelVector(); } BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector() { // if not already generated start voxelization here updateMask(); return _spRelevantVoxelVector; } BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector( float lowerThreshold) { auto filteredVoxelVectorPointer = ::boost::make_shared(); updateMask(); // filter relevant voxels auto it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getRelevantVolumeFraction() > lowerThreshold) { filteredVoxelVectorPointer->push_back(*it); } ++it; } // if mask calculation was not successful this is empty! return filteredVoxelVectorPointer; } bool BoostMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const { //initialize return voxel voxel.setRelevantVolumeFraction(0); //check if ID is valid if (!_geoInfo.validID(aID)) { return false; } //determine how a given voxel on the dose grid is masked if (_spRelevantVoxelVector) { auto it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getVoxelGridID() == aID) { voxel = (*it); return true; } ++it; } } // returns false if mask was not calculated without triggering calculation (otherwise not const!) else { return false; } return false; } bool BoostMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const { //convert VoxelGridIndex3D to VoxelGridID VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getMaskAt(aVoxelGridID, voxel); } else { return false; } } const core::GeometricInfo& BoostMaskAccessor::getGeometricInfo() const { return _geoInfo; }; } } } \ No newline at end of file diff --git a/code/masks/rttbBoostMaskVoxelizationThread.cpp b/code/masks/rttbBoostMaskVoxelizationThread.cpp index 693898b..801005f 100644 --- a/code/masks/rttbBoostMaskVoxelizationThread.cpp +++ b/code/masks/rttbBoostMaskVoxelizationThread.cpp @@ -1,168 +1,167 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbBoostMaskVoxelizationThread.h" #include "rttbInvalidParameterException.h" #include #include -#include namespace rttb { namespace masks { namespace boost { BoostMaskVoxelizationThread::BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex, bool strict) : _geometryCoordinateBoostPolygonMap(APolygonMap), _globalBoundingBox(aGlobalBoundingBox), _resultVoxelization(anArrayMap), _mutex(aMutex), _strict(strict) { } void BoostMaskVoxelizationThread::operator()() { rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); const unsigned int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; const unsigned int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; std::map > voxelizationMapInThread; for (auto & it : _geometryCoordinateBoostPolygonMap) { BoostArray2D maskArray(::boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); BoostPolygonVector boostPolygonVec = it.second; for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) { for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = x + minIndex[0]; currentIndex[1] = y + minIndex[1]; currentIndex[2] = 0; //Get intersection polygons of the dose voxel and the structure BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec); //Calc areas of all intersection polygons double volumeFraction = calcArea(polygons); volumeFraction = correctForErrorAndStrictness(volumeFraction, _strict); if (volumeFraction < 0 || volumeFraction > 1 ) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } maskArray[x][y] = volumeFraction; } } voxelizationMapInThread.insert(std::pair(it.first, ::boost::make_shared(maskArray))); } //insert gathered values into voxelization map ::boost::unique_lock<::boost::shared_mutex> lock(*_mutex); _resultVoxelization->insert(voxelizationMapInThread.begin(), voxelizationMapInThread.end()); } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMaskVoxelizationThread::BoostPolygonDeque BoostMaskVoxelizationThread::getIntersections( const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) { BoostMaskVoxelizationThread::BoostPolygonDeque polygonDeque; BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); ::boost::geometry::correct(voxelPolygon); BoostPolygonVector::const_iterator it; for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) { BoostPolygon2D contour = *it; ::boost::geometry::correct(contour); BoostPolygonDeque intersection; ::boost::geometry::intersection(voxelPolygon, contour, intersection); polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end()); } return polygonDeque; } BoostMaskVoxelizationThread::BoostRing2D BoostMaskVoxelizationThread::get2DContour( const rttb::VoxelGridIndex3D& aVoxelGrid3D) { BoostRing2D polygon; BoostPoint2D point1(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] - 0.5); ::boost::geometry::append(polygon, point1); BoostPoint2D point2(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] - 0.5); ::boost::geometry::append(polygon, point2); BoostPoint2D point3(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] + 0.5); ::boost::geometry::append(polygon, point3); BoostPoint2D point4(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] + 0.5); ::boost::geometry::append(polygon, point4); ::boost::geometry::append(polygon, point1); return polygon; } /*Calculate the intersection area*/ double BoostMaskVoxelizationThread::calcArea(const BoostPolygonDeque& aPolygonDeque) { double area = 0; BoostPolygonDeque::const_iterator it; for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) { area += ::boost::geometry::area(*it); } return area; } double BoostMaskVoxelizationThread::correctForErrorAndStrictness(double volumeFraction, bool strict) const { if (strict){ if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } } else { if (volumeFraction > 1){ volumeFraction = 1; } else if (volumeFraction < 0){ volumeFraction = 0; } } return volumeFraction; } } } } diff --git a/code/masks/rttbGenericMutableMaskAccessor.cpp b/code/masks/rttbGenericMutableMaskAccessor.cpp index 16f2eac..a3a8d68 100644 --- a/code/masks/rttbGenericMutableMaskAccessor.cpp +++ b/code/masks/rttbGenericMutableMaskAccessor.cpp @@ -1,173 +1,170 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbGenericMutableMaskAccessor.h" -#include "rttbMaskVoxel.h" -#include "rttbNullPointerException.h" #include -#include #include #include namespace rttb { namespace masks { GenericMutableMaskAccessor::GenericMutableMaskAccessor(const core::GeometricInfo& aGeometricInfo) : _geoInfo(aGeometricInfo), _spRelevantVoxelVector(MaskVoxelListPointer()) { //generate new structure set uid boost::uuids::uuid id; boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _maskUID = "GenericMutableMask_" + ss.str(); } GenericMutableMaskAccessor::~GenericMutableMaskAccessor() = default; void GenericMutableMaskAccessor::updateMask() {} GenericMutableMaskAccessor::MaskVoxelListPointer GenericMutableMaskAccessor::getRelevantVoxelVector() { return _spRelevantVoxelVector; } GenericMutableMaskAccessor::MaskVoxelListPointer GenericMutableMaskAccessor::getRelevantVoxelVector( float lowerThreshold) { auto filteredVoxelVectorPointer = boost::make_shared(); // filter relevant voxels auto it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getRelevantVolumeFraction() > lowerThreshold) { filteredVoxelVectorPointer->push_back(*it); } ++it; } // if mask calculation was not successful this is empty! return filteredVoxelVectorPointer; } bool GenericMutableMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const { //initialize return voxel voxel.setRelevantVolumeFraction(0); //check if ID is valid if (!_geoInfo.validID(aID)) { return false; } //determine how a given voxel on the dose grid is masked if (_spRelevantVoxelVector) { auto it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getVoxelGridID() == aID) { voxel = (*it); return true; } ++it; } //aID is not in mask voxel.setRelevantVolumeFraction(0); } return false; } bool GenericMutableMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const { //convert VoxelGridIndex3D to VoxelGridID VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getMaskAt(aVoxelGridID, voxel); } else { return false; } } void GenericMutableMaskAccessor::setMaskAt(VoxelGridID aID, const core::MaskVoxel& voxel) { //check if ID is valid if (!_geoInfo.validID(aID)) { return; } //determine how a given voxel on the dose grid is masked if (_spRelevantVoxelVector) { auto it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getVoxelGridID() == aID) { (*it) = voxel; return; } ++it; } //not sID is not found in existing voxels _spRelevantVoxelVector->push_back(voxel); } } void GenericMutableMaskAccessor::setMaskAt(const VoxelGridIndex3D& aIndex, const core::MaskVoxel& voxel) { //convert VoxelGridIndex3D to VoxelGridID VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { setMaskAt(aVoxelGridID, voxel); } } void GenericMutableMaskAccessor::setRelevantVoxelVector(MaskVoxelListPointer aVoxelListPointer) { _spRelevantVoxelVector = MaskVoxelListPointer(aVoxelListPointer); } } } \ No newline at end of file diff --git a/code/models/rttbBioModel.cpp b/code/models/rttbBioModel.cpp index 3417517..8eea30b 100644 --- a/code/models/rttbBioModel.cpp +++ b/code/models/rttbBioModel.cpp @@ -1,55 +1,54 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ -#include #define _USE_MATH_DEFINES -#include -#include + +#include #include "rttbBioModel.h" namespace rttb { namespace models { bool BioModel::init(const double doseFactor) { if (_dvh != nullptr) { _value = this->calcModel(doseFactor); return true; } return false; } void BioModel::setDVH(const DVHPointer aDVH) { _dvh = aDVH; } const BioModel::DVHPointer BioModel::getDVH() const { return _dvh; } const BioModelValueType BioModel::getValue() const { return _value; } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbBioModelCurve.cpp b/code/models/rttbBioModelCurve.cpp index 24ca50c..3fc1eb9 100644 --- a/code/models/rttbBioModelCurve.cpp +++ b/code/models/rttbBioModelCurve.cpp @@ -1,87 +1,85 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include - #include -#include "rttbBioModel.h" #include "rttbNTCPLKBModel.h" #include "rttbBioModelCurve.h" #include "rttbDvhBasedModels.h" namespace rttb { namespace models { CurveDataType getCurveDoseVSBioModel(BioModel& aModel, double normalisationDose, int aBin, double minDose, double maxDose) { CurveDataType _curveData; double minFactor = 0, maxFactor = 0; minFactor = minDose / normalisationDose; maxFactor = maxDose / normalisationDose; double lastValue = 0; for (int i = 0; i < aBin; i++) { double factor = minFactor + i * (maxFactor - minFactor) / (aBin - 1); if (lastValue < 1) { aModel.init(factor); BioModelValueType value = aModel.getValue(); lastValue = value; _curveData.insert(CurvePointType(factor * normalisationDose, value)); } else { _curveData.insert(CurvePointType(factor * normalisationDose, 1)); } } return _curveData; } CurveDataType getCurveEUDVSBioModel(NTCPLKBModel& aModel, DoseCalcType maxFactor, DoseCalcType minFactor, int aBin) { CurveDataType _curveData; BioModel::DVHPointer _dvh = aModel.getDVH(); for (int i = 0; i < aBin; i++) { DoseCalcType factor = minFactor + i * (maxFactor - minFactor) / (aBin - 1); core::DVH variantDVH = core::DVH(_dvh->getDataDifferential(), _dvh->getDeltaD() * factor, _dvh->getDeltaV(), "temporary", "temporary"); auto spDVH = boost::make_shared(variantDVH); double eud = getEUD(spDVH, aModel.getA()); aModel.init(factor); BioModelValueType value = aModel.getValue(); _curveData.insert(CurvePointType(eud, value)); } return _curveData; } } } \ No newline at end of file diff --git a/code/models/rttbBioModelScatterPlots.cpp b/code/models/rttbBioModelScatterPlots.cpp index 67565ef..360a01e 100644 --- a/code/models/rttbBioModelScatterPlots.cpp +++ b/code/models/rttbBioModelScatterPlots.cpp @@ -1,214 +1,212 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include - #include -#include #include "rttbBioModel.h" #include "rttbBioModelScatterPlots.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace models { /* Initiate Random Number generator with current time */ boost::random::mt19937 rng(static_cast(time(nullptr))); /* Generate random number between 0 and 1 */ boost::random::uniform_01 uniDist(rng); ScatterPlotType getScatterPlotVary1Parameter(BioModel& aModel, int aParamId, BioModelParamType aMean, BioModelParamType aVariance, DoseTypeGy aNormalisationDose, int numberOfPoints, DoseTypeGy aMinDose, DoseTypeGy aMaxDose) { ScatterPlotType scatterPlotData; if (aVariance == 0) { //set to small positive value to avoid negative infinity! aVariance = 1e-30; } if (aMaxDose <= aMinDose) { throw core::InvalidParameterException("Parameter invalid: aMaxDose must be > aMinDose!"); } if (aNormalisationDose <= 0) { throw core::InvalidParameterException("Parameter invalid: aNormalisationDose must be > 0!"); } /* Choose Normal Distribution */ boost::random::normal_distribution gaussian_dist(0, aVariance); /* Create a Gaussian Random Number generator * by binding with previously defined * normal distribution object */ boost::random::variate_generator > generator( rng, gaussian_dist); int i = 0; while (i < numberOfPoints) { double paramValue, probability; double randomValue = generator(); paramValue = randomValue + aMean; probability = normal_pdf(randomValue, aVariance); if (probability > 1) { probability = 1; } //randomly select a dose between aMinDose and aMaxDose double dose = uniDist() * (aMaxDose - aMinDose) + aMinDose; if (probability > 0) { try { aModel.setParameterByID(aParamId, paramValue); aModel.init(dose / aNormalisationDose); double value = aModel.getValue(); std::pair modelProbPair = std::make_pair(value, probability); scatterPlotData.insert(std::pair >(dose, modelProbPair)); i++; } catch (core::InvalidParameterException& /*e*/) { //repeat evaluation to guarantee the correct number of scatter values continue; } } } return scatterPlotData; } double normal_pdf(double aValue, double aVariance) { static const double inv_sqrt_2pi = 0.3989422804014327; double a = (aValue) / aVariance; return inv_sqrt_2pi / aVariance * std::exp(-0.5f * a * a); } ScatterPlotType getScatterPlotVaryParameters(BioModel& aModel, std::vector aParamIdVec, BioModel::ParamVectorType aMeanVec, BioModel::ParamVectorType aVarianceVec, DoseTypeGy aNormalisationDose, int numberOfPoints, DoseTypeGy aMinDose, DoseTypeGy aMaxDose) { ScatterPlotType scatterPlotData; if (aMaxDose <= aMinDose) { throw core::InvalidParameterException("Parameter invalid: aMaxDose must be > aMinDose!"); } if (aNormalisationDose <= 0) { throw core::InvalidParameterException("Parameter invalid: aNormalisationDose must be > 0!"); } //all input vectors need to have the same size if (((aVarianceVec.size() != aMeanVec.size()) || (aVarianceVec.size() != aParamIdVec.size()))) { throw core::InvalidParameterException("Parameter vectors have different sizes!"); } for (double & v : aVarianceVec) { //set to small positive value to avoid negative infinity! if (v == 0) { v = 1e-30; } } double paramValue; // vary all parameters for each scattered point int i = 0; while (i < numberOfPoints) { double probability = 1; for (GridIndexType j = 0; j < aParamIdVec.size(); j++) { /* Choose Normal Distribution */ boost::random::normal_distribution gaussian_dist(0, aVarianceVec.at(j)); /* Create a Gaussian Random Number generator * by binding with previously defined * normal distribution object */ boost::random::variate_generator > generator( rng, gaussian_dist); double randomValue = generator(); paramValue = randomValue + aMeanVec.at(j); if (aVarianceVec.at(j) != 0) { /* calculate combined probability */ probability = probability * normal_pdf(randomValue, aVarianceVec.at(j)); } else { throw core::InvalidParameterException("Parameter invalid: Variance should not be 0!"); } aModel.setParameterByID(aParamIdVec.at(j), paramValue); } //randomly select a dose between aMinDose and aMaxDose double dose = uniDist() * (aMaxDose - aMinDose) + aMinDose; if (probability > 0) { try { aModel.init(dose / aNormalisationDose); double value = aModel.getValue(); std::pair modelProbPair = std::make_pair(value, probability); scatterPlotData.insert(std::pair >(dose, modelProbPair)); i++; } catch (core::InvalidParameterException& /*e*/) { //repeat evaluation to guarantee the correct number of scatter values continue; } } } return scatterPlotData; } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbIntegration.cpp b/code/models/rttbIntegration.cpp index 115467c..9199fe6 100644 --- a/code/models/rttbIntegration.cpp +++ b/code/models/rttbIntegration.cpp @@ -1,186 +1,185 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include -#include #include "rttbIntegration.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace models { double tcpModelFunctor::calculate(double x) const { if (x == 0) { x = 1e-30; } return tcpFunction(a + (1 - x) / x, params) / (x * x); } double LkbModelFunctor::calculate(double x) const { if (x == 0) { x = 1e-30; } return lkbFunction(b - (1 - x) / x) / (x * x); } double tcpFunction(double x, const TcpParams& tcp_params) { double alphaVariance = tcp_params.alphaVariance; if (alphaVariance == 0) { alphaVariance = 1e-30; } double f = exp(-pow((x - tcp_params.alphaMean) / alphaVariance, 2) / 2); for (size_t i = 0; i < tcp_params.volumeVector.size(); ++i) { double tmp, tmp1, tmp2, tmp3; tmp1 = exp(-x * tcp_params.bedVector.at(i)); tmp2 = -(tcp_params.rho) * tcp_params.volumeVector.at(i); tmp3 = tmp2 * tmp1; tmp = exp(tmp3); if (tmp != 1) { f = f * tmp; } } return f; } double integrateTCP(double a, const TcpParams& params) { double aNew = 1e-30; double bNew = 1.0; tcpModelFunctor BMFunction; BMFunction.params = params; BMFunction.a = a; return iterativeIntegration(BMFunction, aNew, bNew); } double lkbFunction(double x) { double tmp = -pow(x, 2) / 2; double step = exp(tmp); return step; } double integrateLKB(double b) { double aNew = 1e-30; double bNew = 1.0; LkbModelFunctor BMFunction; BMFunction.b = b; return iterativeIntegration(BMFunction, aNew, bNew); } //returns the nth stage of refinement of the extended trapezoidal rule template integrationType trapzd(const FunctorType& BMfunction, integrationType a, integrationType b, int stepNum) { static integrationType result; if (stepNum == 1) { result = 0.5 * (b - a) * (BMfunction.calculate(a) + BMfunction.calculate(b)); } else { integrationType x, tnm, sum, del; int it, j; for (it = 1, j = 1; j < stepNum - 1; j++) { it <<= 1; } tnm = it; del = (b - a) / tnm; x = a + 0.5 * del; for (sum = 0.0, j = 0; j < it; j++, x += del) { sum += BMfunction.calculate(x); } result = 0.5 * (result + (b - a) * sum / tnm); } return result; } template integrationType iterativeIntegration(const FunctorType& BMfunction, integrationType a, integrationType b) { integrationType ost = 0.0; integrationType os = 0.0; unsigned int maxSteps = 50; double eps = 1e-6; unsigned int i = 1; for (; i <= maxSteps; ++i) { integrationType st = trapzd(BMfunction, a, b, i); integrationType s = (4.0 * st - ost) / 3.0; if (i > 5) { if (fabs(s - os) < eps * fabs(os) || (s == 0.0 && os == 0.0)) { return s; } } os = s; ost = st; } //too many iterations, this should never be reachable! throw rttb::core::InvalidParameterException("Integral calculation failed: too many iterations! "); } } } diff --git a/code/models/rttbNTCPLKBModel.cpp b/code/models/rttbNTCPLKBModel.cpp index b7cc495..ec4f65d 100644 --- a/code/models/rttbNTCPLKBModel.cpp +++ b/code/models/rttbNTCPLKBModel.cpp @@ -1,175 +1,170 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #define _USE_MATH_DEFINES #include -#include -#include - #include -#include - #include "rttbIntegration.h" #include "rttbNTCPLKBModel.h" #include "rttbDvhBasedModels.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace models { NTCPLKBModel::NTCPLKBModel() : NTCPModel(){ _name = "NTCPLKBModel"; fillParameterMap(); } NTCPLKBModel::NTCPLKBModel(DVHPointer aDvh, BioModelParamType aD50, BioModelParamType aM, BioModelParamType aA): NTCPModel(aDvh, aD50), _m(aM), _a(aA) { _name = "NTCPLKBModel"; fillParameterMap(); } void NTCPLKBModel::setA(const BioModelParamType aA) { _a = aA; } const BioModelParamType NTCPLKBModel::getA() { return _a; } void NTCPLKBModel::setM(const BioModelParamType aM) { _m = aM; } const BioModelParamType NTCPLKBModel::getM() { return _m; } void NTCPLKBModel::setParameterVector(const ParamVectorType& aParameterVector) { if (aParameterVector.size() != 3) { throw core::InvalidParameterException("Parameter invalid: aParameterVector.size must be 3! "); } else { _d50 = aParameterVector.at(0); _m = aParameterVector.at(1); _a = aParameterVector.at(2); } } void NTCPLKBModel::setParameterByID(const int aParamId, const BioModelParamType aValue) { if (aParamId == 0) { _d50 = aValue; } else if (aParamId == 1) { _m = aValue; } else if (aParamId == 2) { _a = aValue; } else { throw core::InvalidParameterException("Parameter invalid: aParamID must be 0(for d50) or 1(for m) or 2(for a)! "); } } const int NTCPLKBModel::getParameterID(const std::string& aParamName) const { if (aParamName == "d50") { return 0; } else if (aParamName == "m") { return 1; } else if (aParamName == "a") { return 2; } else { rttbExceptionMacro(core::InvalidParameterException, << "Parameter name " << aParamName << " invalid: it should be d50 or m or a!"); } } std::map NTCPLKBModel::getParameterMap() const{ return parameterMap; } void NTCPLKBModel::fillParameterMap(){ parameterMap["d50"] = getD50(); parameterMap["m"] = getM(); parameterMap["a"] = getA(); } std::string NTCPLKBModel::getModelType() const{ return _name; } BioModelValueType NTCPLKBModel::calcModel(const double doseFactor) { if (_a == 0) { throw core::InvalidParameterException("_a must not be zero"); } if (_m == 0) { throw core::InvalidParameterException("_m must not be zero"); } core::DVH variantDVH = core::DVH(_dvh->getDataDifferential(), (DoseTypeGy)(_dvh->getDeltaD() * doseFactor), _dvh->getDeltaV(), "temporary", "temporary"); auto spDVH = boost::make_shared(variantDVH); double eud = getEUD(spDVH, this->_a); //_m must not be zero double t = (eud - this->_d50) / (this->_m * this->_d50); double value = 1 / pow(2 * M_PI, 0.5); double result = integrateLKB(t); if (result != -100) { value *= result; return value; } else { return false; } } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbNTCPRSModel.cpp b/code/models/rttbNTCPRSModel.cpp index e2b549e..8779a33 100644 --- a/code/models/rttbNTCPRSModel.cpp +++ b/code/models/rttbNTCPRSModel.cpp @@ -1,168 +1,166 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #define _USE_MATH_DEFINES #include -#include - #include "rttbNTCPRSModel.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace models { NTCPRSModel::NTCPRSModel() : NTCPModel() { _name = "NTCPRSModel"; fillParameterMap(); } NTCPRSModel::NTCPRSModel(DVHPointer aDvh, BioModelParamType aD50, BioModelParamType aGamma, BioModelParamType aS): NTCPModel(aDvh, aD50), _gamma(aGamma), _s(aS) { _name = "NTCPRSModel"; fillParameterMap(); } void NTCPRSModel::setGamma(const BioModelParamType aGamma) { _gamma = aGamma; } const BioModelParamType NTCPRSModel::getGamma() { return _gamma; } void NTCPRSModel::setS(const BioModelParamType aS) { _s = aS; } const BioModelParamType NTCPRSModel::getS() { return _s; } void NTCPRSModel::setParameterVector(const ParamVectorType& aParameterVector) { if (aParameterVector.size() != 3) { throw core::InvalidParameterException("Parameter invalid: aParameterVector.size must be 3! "); } else { _d50 = aParameterVector.at(0); _gamma = aParameterVector.at(1); _s = aParameterVector.at(2); } } void NTCPRSModel::setParameterByID(const int aParamId, const BioModelParamType aValue) { if (aParamId == 0) { _d50 = aValue; } else if (aParamId == 1) { _gamma = aValue; } else if (aParamId == 2) { _s = aValue; } else { throw core::InvalidParameterException("Parameter invalid: aParamID must be 0(for d50) or 1(for gamma) or 2(for s)! "); } } const int NTCPRSModel::getParameterID(const std::string& aParamName) const { if (aParamName == "d50") { return 0; } else if (aParamName == "gamma") { return 1; } else if (aParamName == "s") { return 2; } else { rttbExceptionMacro(core::InvalidParameterException, << "Parameter name " << aParamName << " invalid: it should be d50 or gamma or s!"); } } std::map NTCPRSModel::getParameterMap() const{ return parameterMap; } void NTCPRSModel::fillParameterMap() { parameterMap["d50"] = getD50(); parameterMap["gamma"] = getGamma(); parameterMap["s"] = getS(); } std::string NTCPRSModel::getModelType() const{ return _name; } const double NTCPRSModel::poissonModel(const double dose) { //_d50 must not be zero return pow(2, -exp(M_E * this->_gamma * (1 - dose / this->_d50))); } BioModelValueType NTCPRSModel::calcModel(double doseFactor) { if (_d50 == 0) { throw core::InvalidParameterException("d50 must not be zero"); } if (_s == 0) { throw core::InvalidParameterException("s must not be zero"); } std::deque dataDifferential = this->_dvh->getDataDifferential(); double ntcp = 1; for (GridIndexType i = 0; i < dataDifferential.size(); i++) { double pd = pow(this->poissonModel(i * this->_dvh->getDeltaD() * doseFactor), this->_s); double vi = dataDifferential[i] / this->_dvh->getNumberOfVoxels(); ntcp *= pow((1 - pd), vi); } //_s must not be zero return (BioModelValueType)(pow((1 - ntcp), 1 / this->_s)); } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbTCPLQModel.cpp b/code/models/rttbTCPLQModel.cpp index 9651baa..fedebbb 100644 --- a/code/models/rttbTCPLQModel.cpp +++ b/code/models/rttbTCPLQModel.cpp @@ -1,304 +1,299 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #define _USE_MATH_DEFINES #include - -#include -#include - #include -#include #include "rttbTCPLQModel.h" #include "rttbDvhBasedModels.h" #include "rttbIntegration.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace models { TCPLQModel::TCPLQModel(): TCPModel() { _name = "TCPLQModel"; fillParameterMap(); } TCPLQModel::TCPLQModel(DVHPointer aDVH, BioModelParamType aAlphaMean, BioModelParamType aBeta, BioModelParamType aRho, int aNumberOfFractions): TCPModel(aDVH, aNumberOfFractions), _alphaMean(aAlphaMean), _alphaVariance(0), _alpha_beta(aAlphaMean / aBeta), _rho(aRho) { _name = "TCPLQModel"; fillParameterMap(); } TCPLQModel::TCPLQModel(DVHPointer aDVH, BioModelParamType aRho, int aNumberOfFractions, BioModelParamType aAlpha_Beta, BioModelParamType aAlphaMean, BioModelParamType aAlphaVariance): TCPModel(aDVH, aNumberOfFractions), _alphaMean(aAlphaMean), _alphaVariance(aAlphaVariance), _alpha_beta(aAlpha_Beta), _rho(aRho) { fillParameterMap(); _name = "TCPLQModel"; } void TCPLQModel::setParameters(const BioModelParamType aAlphaMean, const BioModelParamType aAlpha_Beta, const BioModelParamType aRho, const BioModelParamType aAlphaVariance) { _alphaMean = aAlphaMean; _alphaVariance = aAlphaVariance; _alpha_beta = aAlpha_Beta; _rho = aRho; //reset _value, because parameters have changed. _value = 0; } void TCPLQModel::setAlpha(const BioModelParamType aAlphaMean, const BioModelParamType aAlphaVariance) { _alphaVariance = aAlphaVariance; _alphaMean = aAlphaMean; } void TCPLQModel::setAlphaBeta(const BioModelParamType aAlpha_Beta) { _alpha_beta = aAlpha_Beta; } void TCPLQModel::setRho(const BioModelParamType aRho) { _rho = aRho; } const BioModelParamType TCPLQModel::getAlphaBeta() { return _alpha_beta; } const BioModelParamType TCPLQModel::getAlphaMean() { return _alphaMean; } const BioModelParamType TCPLQModel::getAlphaVariance() { return _alphaVariance; } const BioModelParamType TCPLQModel::getRho() { return _rho; } long double TCPLQModel::calcTCPi(BioModelParamType aRho, BioModelParamType aAlphaMean, double vj, double bedj) { return exp(-aRho * vj * exp(-aAlphaMean * bedj)); } long double TCPLQModel::calcTCP(std::map aBEDDVH, BioModelParamType aRho, BioModelParamType aAlphaMean, double aDeltaV) { std::map::iterator it; long double tcp = 1; for (it = aBEDDVH.begin(); it != aBEDDVH.end(); ++it) { long double tcpi = this->calcTCPi(aRho, aAlphaMean, (*it).second * aDeltaV, (*it).first); tcp = tcp * tcpi; } return tcp; } long double TCPLQModel::calcTCPAlphaNormalDistribution( std::map aBEDDVH, BioModelParamType aRho, BioModelParamType aAlphaMean, BioModelParamType aAlphaVariance, double aDeltaV) { std::map::iterator it; std::vector volumeV2; std::vector bedV2; int i = 0; for (it = aBEDDVH.begin(); it != aBEDDVH.end(); ++it) { volumeV2.push_back((*it).second * aDeltaV); bedV2.push_back((*it).first); i++; } struct TcpParams params = {aAlphaMean, aAlphaVariance, aRho, volumeV2, bedV2}; double result = integrateTCP(0, params); if (result == -100) { std::cerr << "Integration error!\n"; return -1; } else { long double tcp = 1 / (pow(2 * M_PI, 0.5) * _alphaVariance) * result; return tcp; } } BioModelValueType TCPLQModel::calcModel(const double doseFactor) { core::DVH variantDVH = core::DVH(_dvh->getDataDifferential(), (DoseTypeGy)(_dvh->getDeltaD() * doseFactor), _dvh->getDeltaV(), "temporary", "temporary"); auto spDVH = boost::make_shared(variantDVH); BioModelValueType value = 0; if (_alphaVariance == 0) { if (_alphaMean <= 0 || _alpha_beta <= 0 || _rho <= 0) { throw core::InvalidParameterException("Parameter invalid: alpha, alpha/beta, rho and number of fractions must >0!"); } if (_numberOfFractions <= 1) { throw core::InvalidParameterException("Parameter invalid: numberOfFractions must be >1! The dvh should be an accumulated-dvh of all fractions, not a single fraction-dvh!"); } std::map dataBED = calcBEDDVH(spDVH, _numberOfFractions, _alpha_beta); value = (BioModelValueType)this->calcTCP(dataBED, _rho, _alphaMean, variantDVH.getDeltaV()); return value; } //if alpha normal distribution else { if (this->_alpha_beta <= 0 || this->_alphaMean <= 0 || this->_alphaVariance < 0 || _rho <= 0) { throw core::InvalidParameterException("Parameter invalid: alpha/beta, alphaMean, rho and number of fractions must >0!"); } if (_numberOfFractions <= 1) { throw core::InvalidParameterException("Parameter invalid: numberOfFractions must be >1! The dvh should be an accumulated-dvh of all fractions, not a single fraction-dvh!"); } std::map dataBED = calcBEDDVH(spDVH, _numberOfFractions, _alpha_beta); value = (BioModelValueType)(this->calcTCPAlphaNormalDistribution(dataBED, _rho, _alphaMean, _alphaVariance, variantDVH.getDeltaV())); return value; } } void TCPLQModel::setParameterVector(const ParamVectorType& aParameterVector) { if (aParameterVector.size() != 4) { throw core::InvalidParameterException("Parameter invalid: aParameterVector.size must be 4! "); } else { _alphaMean = aParameterVector.at(0); _alphaVariance = aParameterVector.at(1); _alpha_beta = aParameterVector.at(2); _rho = aParameterVector.at(3); } } void TCPLQModel::setParameterByID(const int aParamId, const BioModelParamType aValue) { if (aParamId == 0) { _alphaMean = aValue; } else if (aParamId == 1) { _alphaVariance = aValue; } else if (aParamId == 2) { _alpha_beta = aValue; } else if (aParamId == 3) { _rho = aValue; } else { throw core::InvalidParameterException("Parameter invalid: aParamID must be 0(alphaMean) or 1(alphaVariance) or 2(alpha_beta) or 3(rho)! "); } } const int TCPLQModel::getParameterID(const std::string& aParamName) const { if (aParamName == "alphaMean") { return 0; } else if (aParamName == "alphaVariance") { return 1; } else if (aParamName == "alpha_beta") { return 2; } else if (aParamName == "rho") { return 3; } else { rttbExceptionMacro(core::InvalidParameterException, << "Parameter name " << aParamName << " invalid: it should be alphaMean or alphaVariance or alpha_beta or rho!"); } } std::map TCPLQModel::getParameterMap() const{ return parameterMap; } void TCPLQModel::fillParameterMap(){ parameterMap["numberOfFraction"] = getNumberOfFractions(); parameterMap["alphaMean"] = getAlphaMean(); parameterMap["alphaVariance"] = getAlphaVariance(); parameterMap["alpha_beta"] = getAlphaBeta(); parameterMap["rho"] = getRho(); } std::string TCPLQModel::getModelType() const{ return _name; } }//end namespace models }//end namespace rttb