diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp index d51bfbc..80e32be 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,625 +1,620 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "rttbDoseStatisticsCalculator.h" #include #include #include #include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace algorithms { DoseStatisticsCalculator::DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator) { if (aDoseIterator == NULL) { throw core::NullPointerException("DoseIterator must not be NULL"); } else { _doseIterator = aDoseIterator; } _simpleDoseStatisticsCalculated = false; } DoseStatisticsCalculator::~DoseStatisticsCalculator() { } DoseStatisticsCalculator::DoseIteratorPointer DoseStatisticsCalculator::getDoseIterator() const { return _doseIterator; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( bool computeComplexMeasures, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } //"simple" dose statistics are mandatory calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); if (computeComplexMeasures) { //more complex dose statistics are optional with default maximum dose and default relative x values calculateComplexDoseStatistics(_statistics->getMaximum(), std::vector(), std::vector()); } - return _statistics; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } if (referenceDose <= 0) { throw rttb::core::InvalidParameterException("Reference dose must be > 0 !"); } //simple dose statistics calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); //more complex dose statistics with given reference dose and default x values calculateComplexDoseStatistics(referenceDose, std::vector(), std::vector()); return _statistics; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues, DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (!_doseIterator) { throw core::NullPointerException("_doseIterator must not be NULL!"); } //"simple" dose statistics calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); - if (referenceDose <= 0) { //more complex dose statistics with default maximum dose and relative x values calculateComplexDoseStatistics(_statistics->getMaximum(), precomputeDoseValues, precomputeVolumeValues); } else { //more complex dose statistics with given reference dose and relative x values calculateComplexDoseStatistics(referenceDose, precomputeDoseValues, precomputeVolumeValues); } - return _statistics; } void DoseStatisticsCalculator::calculateSimpleDoseStatistics(unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { _doseVector.clear(); _voxelProportionVector.clear(); std::multimap doseValueVSIndexMap; std::vector voxelProportionVectorTemp; DoseStatisticType maximumDose = 0; DoseStatisticType minimumDose = std::numeric_limits::max(); DoseStatisticType meanDose; DoseStatisticType stdDeviationDose; DoseTypeGy sum = 0; VolumeType numVoxels = 0.0; DoseTypeGy squareSum = 0; VolumeType volume = 0; _doseIterator->reset(); int i = 0; DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid()) { doseValue = _doseIterator->getCurrentDoseValue(); if (i == 0) { minimumDose = doseValue; volume = _doseIterator->getCurrentVoxelVolume(); } rttb::FractionType voxelProportion = _doseIterator->getCurrentRelevantVolumeFraction(); sum += doseValue * voxelProportion; numVoxels += voxelProportion; squareSum += doseValue * doseValue * voxelProportion; if (doseValue > maximumDose) { maximumDose = doseValue; } else if (doseValue < minimumDose) { minimumDose = doseValue; } voxelProportionVectorTemp.push_back(voxelProportion); doseValueVSIndexMap.insert(std::pair(doseValue, i)); i++; _doseIterator->next(); } if (numVoxels != 0) { meanDose = sum / numVoxels; //standard deviation is defined only for n>=2 if (numVoxels >= 2) { //uncorrected variance is calculated DoseStatisticType varianceDose = (squareSum / numVoxels - meanDose * meanDose); if (varianceDose < errorConstant) { stdDeviationDose = 0; } else { stdDeviationDose = pow(varianceDose, 0.5); } } else { stdDeviationDose = 0; } } //sort dose values and corresponding volume fractions in member variables for (auto it = doseValueVSIndexMap.begin(); it != doseValueVSIndexMap.end(); ++it) { _doseVector.push_back((float)(*it).first); _voxelProportionVector.push_back(voxelProportionVectorTemp.at((*it).second)); } volume *= numVoxels; _statistics = boost::make_shared(minimumDose, maximumDose, meanDose, stdDeviationDose, numVoxels, volume); _simpleDoseStatisticsCalculated = true; ResultListPointer minimumVoxelPositions = computeMinimumPositions(maxNumberMinimaPositions); ResultListPointer maximumVoxelPositions = computeMaximumPositions(maxNumberMaximaPositions); _statistics->setMinimumVoxelPositions(minimumVoxelPositions); _statistics->setMaximumVoxelPositions(maximumVoxelPositions); } void DoseStatisticsCalculator::calculateComplexDoseStatistics(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues) { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call calculateComplexDoseStatistics()"); } std::vector precomputeDoseValuesNonConst = precomputeDoseValues; std::vector precomputeVolumeValuesNonConst = precomputeVolumeValues; //set default values if (precomputeDoseValues.empty()) { std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)( 0.95)(0.98); precomputeDoseValuesNonConst = defaultPrecomputeDoseValues; } if (precomputeVolumeValues.empty()) { std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)( 0.95)(0.98); precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues; } DoseToVolumeFunctionType Vx = computeDoseToVolumeFunctionMulti(referenceDose, precomputeDoseValuesNonConst, DoseStatistics::Vx); VolumeToDoseFunctionType Dx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::Dx); VolumeToDoseFunctionType MOHx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MOHx); VolumeToDoseFunctionType MOCx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MOCx); VolumeToDoseFunctionType MaxOHx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MaxOHx); VolumeToDoseFunctionType MinOCx = computeVolumeToDoseFunctionMulti(precomputeVolumeValuesNonConst, DoseStatistics::MinOCx); _statistics->setVx(Vx); _statistics->setDx(Dx); _statistics->setMOHx(MOHx); _statistics->setMOCx(MOCx); _statistics->setMaxOHx(MaxOHx); _statistics->setMinOCx(MinOCx); _statistics->setReferenceDose(referenceDose); } DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMaximumPositions( unsigned int maxNumberMaxima) const { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMaximumPositions()"); } ResultListPointer maxVoxelVector = boost::make_shared > >(); unsigned int count = 0; this->_doseIterator->reset(); DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid() && count < maxNumberMaxima) { doseValue = _doseIterator->getCurrentDoseValue(); if (doseValue == _statistics->getMaximum()) { VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID(); std::pair voxel(doseValue, currentID); maxVoxelVector->push_back(voxel); count++; } _doseIterator->next(); } return maxVoxelVector; } DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMinimumPositions( unsigned int maxNumberMinima) const { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMinimumPositions()"); } ResultListPointer minVoxelVector = boost::make_shared > >(); /*! @todo: Architecture Annotation: Finding the positions for the minimum only once reduces computation time, but will require sensible use by the programmers. To be save the output vector minVoxelVector will be always cleared here to garantee that no false values are presented. This change may be revoced to increase computation speed later on (only compute if(minVoxelVector->size()==0)). */ unsigned int count = 0; this->_doseIterator->reset(); DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid() && count < maxNumberMinima) { doseValue = _doseIterator->getCurrentDoseValue(); if (doseValue == _statistics->getMinimum()) { VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID(); std::pair voxel(doseValue, currentID); minVoxelVector->push_back(voxel); count++; } _doseIterator->next(); } return minVoxelVector; } VolumeType DoseStatisticsCalculator::computeVx(DoseTypeGy xDoseAbsolute) const { rttb::FractionType count = 0; _doseIterator->reset(); DoseTypeGy currentDose = 0; while (_doseIterator->isPositionValid()) { currentDose = _doseIterator->getCurrentDoseValue(); if (currentDose >= xDoseAbsolute) { count += _doseIterator->getCurrentRelevantVolumeFraction(); } _doseIterator->next(); } return count * this->_doseIterator->getCurrentVoxelVolume(); } DoseTypeGy DoseStatisticsCalculator::computeDx(VolumeType xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; - int i = static_cast(_doseVector.size() - 1); - + int i = _doseVector.size() - 1; for (; i >= 0; i--) { countVoxels += _voxelProportionVector.at(i); if (countVoxels >= noOfVoxel) { break; } } if (i >= 0) { resultDose = _doseVector.at(i); } else { resultDose = _statistics->getMinimum(); } return resultDose; } DoseTypeGy DoseStatisticsCalculator::computeMOHx(VolumeType xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); if (noOfVoxel == 0) { return 0; } else { double countVoxels = 0; double sum = 0; - for (int i = static_cast(_doseVector.size() - 1); i >= 0; i--) + for (size_t i = _doseVector.size() - 1; i >= 0; i--) { double voxelProportion = _voxelProportionVector.at(i); countVoxels += voxelProportion; sum += _doseVector.at(i) * voxelProportion; if (countVoxels >= noOfVoxel) { break; } } return (DoseTypeGy)(sum / noOfVoxel); } } DoseTypeGy DoseStatisticsCalculator::computeMOCx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); if (noOfVoxel == 0) { return 0; } else { double countVoxels = 0; double sum = 0; std::vector::const_iterator it = _doseVector.begin(); std::vector::const_iterator itD = _voxelProportionVector.begin(); for (; it != _doseVector.end(); ++it, ++itD) { double voxelProportion = *itD; countVoxels += voxelProportion; sum += (*it) * voxelProportion; if (countVoxels >= noOfVoxel) { break; } } return (DoseTypeGy)(sum / noOfVoxel); } } DoseTypeGy DoseStatisticsCalculator::computeMaxOHx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; - int i = static_cast(_doseVector.size() - 1); - + int i = _doseVector.size() - 1; for (; i >= 0; i--) { countVoxels += _voxelProportionVector.at(i); if (countVoxels >= noOfVoxel) { break; } } if (i - 1 >= 0) { resultDose = _doseVector.at(i - 1); } return resultDose; } DoseTypeGy DoseStatisticsCalculator::computeMinOCx(DoseTypeGy xVolumeAbsolute) const { double noOfVoxel = xVolumeAbsolute / _doseIterator->getCurrentVoxelVolume(); DoseTypeGy resultDose = 0; double countVoxels = 0; std::vector::const_iterator it = _doseVector.begin(); std::vector::const_iterator itD = _voxelProportionVector.begin(); for (; itD != _voxelProportionVector.end(); ++itD, ++it) { countVoxels += *itD; if (countVoxels >= noOfVoxel) { break; } } if (it != _doseVector.end()) { ++it; if (it != _doseVector.end()) { resultDose = *it; } else { resultDose = (DoseTypeGy)_statistics->getMaximum(); } } else { resultDose = (DoseTypeGy)_statistics->getMinimum(); } return resultDose; } DoseStatisticsCalculator::DoseToVolumeFunctionType DoseStatisticsCalculator::computeDoseToVolumeFunctionMulti(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, DoseStatistics::complexStatistics name) const { DoseToVolumeFunctionType VxMulti; - for (int i = 0; i < precomputeDoseValues.size(); ++i) + for (size_t i = 0; i < precomputeDoseValues.size(); ++i) { if (name == DoseStatistics::Vx) { double xAbsolue = precomputeDoseValues.at(i) * referenceDose; VxMulti.insert(std::pair(xAbsolue, computeVx(xAbsolue))); } else { throw core::InvalidParameterException("unknown DoseStatistics name!"); } } return VxMulti; } DoseStatisticsCalculator::VolumeToDoseFunctionType DoseStatisticsCalculator::computeVolumeToDoseFunctionMulti( const std::vector& precomputeVolumeValues, DoseStatistics::complexStatistics name) const { VolumeToDoseFunctionType multiValues; VolumeType volume = _statistics->getVolume(); - for (int i = 0; i < precomputeVolumeValues.size(); ++i) + for (size_t i = 0; i < precomputeVolumeValues.size(); ++i) { double xAbsolute = precomputeVolumeValues.at(i) * volume; switch (name) { case DoseStatistics::Dx: multiValues.insert(std::pair(xAbsolute, computeDx(xAbsolute))); break; case DoseStatistics::MOHx: multiValues.insert(std::pair(xAbsolute, computeMOHx(xAbsolute))); break; case DoseStatistics::MOCx: multiValues.insert(std::pair(xAbsolute, computeMOCx(xAbsolute))); break; case DoseStatistics::MaxOHx: multiValues.insert(std::pair(xAbsolute, computeMaxOHx(xAbsolute))); break; case DoseStatistics::MinOCx: multiValues.insert(std::pair(xAbsolute, computeMinOCx(xAbsolute))); break; default: throw core::InvalidParameterException("unknown DoseStatistics name!"); } } return multiValues; } }//end namespace algorithms }//end namespace rttb diff --git a/code/core/rttbBaseType.h b/code/core/rttbBaseType.h index b338a3e..5096805 100644 --- a/code/core/rttbBaseType.h +++ b/code/core/rttbBaseType.h @@ -1,742 +1,742 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __BASE_TYPE_NEW_H #define __BASE_TYPE_NEW_H #include #include #include #include #include #include #include #include #include #include namespace rttb { const double errorConstant = 1e-5; typedef unsigned short UnsignedIndex1D; /*! @class UnsignedIndex2D @brief 2D index. */ class UnsignedIndex2D: public boost::numeric::ublas::vector { public: UnsignedIndex2D() : boost::numeric::ublas::vector(2) {} UnsignedIndex2D(const UnsignedIndex1D value) : boost::numeric::ublas::vector(2, value) {} const UnsignedIndex1D x() const { return (*this)(0); } const UnsignedIndex1D y() const { return (*this)(1); } }; /*! @class UnsignedIndex3D @brief 3D index. */ class UnsignedIndex3D: public boost::numeric::ublas::vector { public: UnsignedIndex3D() : boost::numeric::ublas::vector(3) {} UnsignedIndex3D(const UnsignedIndex1D value) : boost::numeric::ublas::vector(3, value) {} UnsignedIndex3D(const UnsignedIndex1D xValue, const UnsignedIndex1D yValue, const UnsignedIndex1D zValue) : boost::numeric::ublas::vector(3, xValue) { (*this)(1) = yValue; (*this)(2) = zValue; } const UnsignedIndex1D x() const { return (*this)(0); } const UnsignedIndex1D y() const { return (*this)(1); } const UnsignedIndex1D z() const { return (*this)(2); } friend bool operator==(const UnsignedIndex3D& gi1, const UnsignedIndex3D& gi2) { if (gi1.size() != gi2.size()) { return false; } - for (int i = 0; i < gi1.size(); i++) + for (size_t i = 0; i < gi1.size(); i++) { if (gi1(i) != gi2(i)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const UnsignedIndex3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; typedef std::list UnsignedIndexList; typedef std::string FileNameString; typedef std::string ContourGeometricTypeString; typedef double WorldCoordinate; /*! @class WorldCoordinate2D @brief 2D coordinate in real world coordinates. */ class WorldCoordinate2D: public boost::numeric::ublas::vector { public: WorldCoordinate2D() : boost::numeric::ublas::vector (2) {} WorldCoordinate2D(const WorldCoordinate value) : boost::numeric::ublas::vector(2, value) {} WorldCoordinate2D(const WorldCoordinate xValue, const WorldCoordinate yValue) : boost::numeric::ublas::vector(2, xValue) { (*this)(1) = yValue; } const WorldCoordinate x() const { return (*this)(0); } const WorldCoordinate y() const { return (*this)(1); } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y(); return ss.str(); } friend bool operator==(const WorldCoordinate2D& wc1, const WorldCoordinate2D& wc2) { if (wc1.size() != wc2.size()) { return false; } - for (int i = 0; i < wc1.size(); i++) + for (size_t i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } }; /*! @class WorldCoordinate3D @brief 3D coordinate in real world coordinates. */ class WorldCoordinate3D: public boost::numeric::ublas::vector { public: WorldCoordinate3D() : boost::numeric::ublas::vector(3) {} WorldCoordinate3D(const WorldCoordinate value) : boost::numeric::ublas::vector(3, value) {} WorldCoordinate3D(const WorldCoordinate xValue, const WorldCoordinate yValue, const WorldCoordinate zValue) : boost::numeric::ublas::vector(3, xValue) { (*this)(1) = yValue; (*this)(2) = zValue; } WorldCoordinate3D(const WorldCoordinate3D& w): boost::numeric::ublas::vector(3) { (*this)(0) = w.x(); (*this)(1) = w.y(); (*this)(2) = w.z(); } const WorldCoordinate x() const { return (*this)(0); } const WorldCoordinate y() const { return (*this)(1); } const WorldCoordinate z() const { return (*this)(2); } //vector cross product (not included in boost.ublas) WorldCoordinate3D cross(WorldCoordinate3D aVector) const { WorldCoordinate3D result; WorldCoordinate x = (*this)(0); WorldCoordinate y = (*this)(1); WorldCoordinate z = (*this)(2); result(0) = y * aVector(2) - z * aVector(1); result(1) = z * aVector(0) - x * aVector(2); result(2) = x * aVector(1) - y * aVector(0); return result; } WorldCoordinate2D getXY() const { WorldCoordinate2D result; result(0) = (*this)(0); result(1) = (*this)(1); return result; } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y() << ' ' << z(); return ss.str(); } WorldCoordinate3D& operator=(const WorldCoordinate3D& wc) { (*this)(0) = wc.x(); (*this)(1) = wc.y(); (*this)(2) = wc.z(); return (*this); } WorldCoordinate3D& operator=(const boost::numeric::ublas::vector wc) { (*this)(0) = wc(0); (*this)(1) = wc(1); (*this)(2) = wc(2); return (*this); } WorldCoordinate3D operator-(const boost::numeric::ublas::vector wc) { return WorldCoordinate3D((*this)(0) - wc(0), (*this)(1) - wc(1), (*this)(2) - wc(2)); } WorldCoordinate3D operator+(const boost::numeric::ublas::vector wc) { return WorldCoordinate3D((*this)(0) + wc(0), (*this)(1) + wc(1), (*this)(2) + wc(2)); } friend bool operator==(const WorldCoordinate3D& wc1, const WorldCoordinate3D& wc2) { if (wc1.size() != wc2.size()) { return false; } - for (int i = 0; i < wc1.size(); i++) + for (size_t i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } bool equalsAlmost(const WorldCoordinate3D& another, double errorConstant = 1e-5) const { if (size() != another.size()) { return false; } double dist = norm_2(*this - another); return dist < errorConstant; } friend std::ostream& operator<<(std::ostream& s, const WorldCoordinate3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; typedef WorldCoordinate3D DoubleVoxelGridIndex3D; typedef UnsignedIndex3D ImageSize; /*! @deprecated Use OrientationMatrix instead. */ typedef WorldCoordinate3D ImageOrientation; typedef double GridVolumeType; /*! @class SpacingVectorType3D @brief 3D spacing vector. @pre values of this vector may not be negative. */ class SpacingVectorType3D: public boost::numeric::ublas::vector { public: SpacingVectorType3D() : boost::numeric::ublas::vector(3) {} SpacingVectorType3D(const GridVolumeType value) : boost::numeric::ublas::vector(3, value) {} SpacingVectorType3D(const GridVolumeType xValue, const GridVolumeType yValue, const GridVolumeType zValue) : boost::numeric::ublas::vector(3, xValue) { (*this)(1) = yValue; (*this)(2) = zValue; } SpacingVectorType3D(const SpacingVectorType3D& w): boost::numeric::ublas::vector(3) { (*this)(0) = w.x(); (*this)(1) = w.y(); (*this)(2) = w.z(); } const GridVolumeType x() const { return (*this)(0); } const GridVolumeType y() const { return (*this)(1); } const GridVolumeType z() const { return (*this)(2); } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y() << ' ' << z(); return ss.str(); } SpacingVectorType3D& operator=(const SpacingVectorType3D& wc) { (*this)(0) = wc.x(); (*this)(1) = wc.y(); (*this)(2) = wc.z(); return (*this); } SpacingVectorType3D& operator=(const WorldCoordinate3D& wc) { (*this)(0) = GridVolumeType(wc.x()); (*this)(1) = GridVolumeType(wc.y()); (*this)(2) = GridVolumeType(wc.z()); return (*this); } SpacingVectorType3D& operator=(const boost::numeric::ublas::vector wc) { (*this)(0) = wc(0); (*this)(1) = wc(1); (*this)(2) = wc(2); return (*this); } friend bool operator==(const SpacingVectorType3D& wc1, const SpacingVectorType3D& wc2) { if (wc1.size() != wc2.size()) { return false; } - for (int i = 0; i < wc1.size(); i++) + for (size_t i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } return true; } bool equalsAlmost(const SpacingVectorType3D& another, double errorConstant) const { if ((*this).size() != another.size()) { return false; } double dist = norm_2(*this - another); return dist < errorConstant; } friend std::ostream& operator<<(std::ostream& s, const SpacingVectorType3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; /*! @class OrientationMatrix @brief Used to store image orientation information */ class OrientationMatrix : public boost::numeric::ublas::matrix { public: /*! The default constructor generates a 3x3 unit matrix */ OrientationMatrix() : boost::numeric::ublas::matrix(3, 3, 0) { for (std::size_t m = 0; m < (*this).size1(); m++) { (*this)(m, m) = 1; } } OrientationMatrix(const WorldCoordinate value) : boost::numeric::ublas::matrix(3, 3, value) {} bool equalsAlmost(const OrientationMatrix& anOrientationMatrix, double errorConstant) const { if (anOrientationMatrix.size1() == (*this).size1()) { if (anOrientationMatrix.size2() == (*this).size2()) { for (std::size_t m = 0; m < anOrientationMatrix.size1(); m++) { for (std::size_t n = 0; n < anOrientationMatrix.size2(); n++) { if ((abs((*this)(m, n) - anOrientationMatrix(m, n)) > errorConstant)) { return false; } } }// end element comparison } else { return false; } } else { return false; } return true; } friend bool operator==(const OrientationMatrix& om1, const OrientationMatrix& om2) { return om1.equalsAlmost(om2, 0.0); } friend std::ostream& operator<<(std::ostream& s, const OrientationMatrix& anOrientationMatrix) { s << "[ "; for (std::size_t m = 0; m < anOrientationMatrix.size1(); m++) { s << "[ "; for (std::size_t n = 0; n < anOrientationMatrix.size2(); n++) { if (n == 0) { s << anOrientationMatrix(m, n); } else { s << ", " << anOrientationMatrix(m, n); } } s << " ]"; } s << " ]"; return s; } }; /*! base for 2D and 3D VoxelIndex; Therefore required beside VoxelGridID */ typedef unsigned int GridIndexType; /*! @class VoxelGridIndex3D @brief 3D voxel grid index. */ class VoxelGridIndex3D: public boost::numeric::ublas::vector { public: VoxelGridIndex3D() : boost::numeric::ublas::vector(3) {} VoxelGridIndex3D(const GridIndexType value) : boost::numeric::ublas::vector(3, value) {} VoxelGridIndex3D(const GridIndexType xValue, const GridIndexType yValue, const GridIndexType zValue) : boost::numeric::ublas::vector(3, xValue) { (*this)(1) = yValue; (*this)(2) = zValue; } const GridIndexType x() const { return (*this)(0); } const GridIndexType y() const { return (*this)(1); } const GridIndexType z() const { return (*this)(2); } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y() << ' ' << z(); return ss.str(); } VoxelGridIndex3D& operator=(const UnsignedIndex3D& ui) { (*this)(0) = ui(0); (*this)(1) = ui(1); (*this)(2) = ui(2); return (*this); } friend bool operator==(const VoxelGridIndex3D& gi1, const VoxelGridIndex3D& gi2) { if (gi1.size() != gi2.size()) { return false; } - for (int i = 0; i < gi1.size(); i++) + for (size_t i = 0; i < gi1.size(); i++) { if (gi1(i) != gi2(i)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const VoxelGridIndex3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; /*! @class VoxelGridIndex3D @brief 2D voxel grid index. */ class VoxelGridIndex2D: public boost::numeric::ublas::vector { public: VoxelGridIndex2D() : boost::numeric::ublas::vector(2) {} VoxelGridIndex2D(const GridIndexType value) : boost::numeric::ublas::vector(2, value) {} VoxelGridIndex2D(const GridIndexType xValue, const GridIndexType yValue) : boost::numeric::ublas::vector(2, xValue) { (*this)(1) = yValue; } const GridIndexType x() const { return (*this)(0); } const GridIndexType y() const { return (*this)(1); } const std::string toString() const { std::stringstream ss; ss << x() << ' ' << y(); return ss.str(); } friend bool operator==(const VoxelGridIndex2D& gi1, const VoxelGridIndex2D& gi2) { if (gi1.size() != gi2.size()) { return false; } - for (int i = 0; i < gi1.size(); i++) + for (size_t i = 0; i < gi1.size(); i++) { if (gi1(i) != gi2(i)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const VoxelGridIndex2D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << " ]"; return s; } }; typedef long GridSizeType; typedef int VoxelGridID; //starts from 0 and is continuously counting all positions on the grid typedef int VoxelGridDimensionType; typedef boost::numeric::ublas::vector VoxelGridDimensionVectorType; typedef double FractionType, VoxelSizeType, DVHVoxelNumber; typedef double DoseCalcType, DoseTypeGy, GenericValueType, DoseVoxelVolumeType, VolumeType, GridVolumeType, VoxelNumberType, BEDType, LQEDType; typedef std::string IDType; typedef std::string StructureLabel; struct DVHRole { enum Type { TargetVolume = 1, HealthyTissue = 2, WholeVolume = 4, UserDefined = 128 } Type; }; struct DVHType { enum Type { Differential = 1, Cumulative = 2 } Type; }; typedef std::string PatientInfoString; typedef std::string TimeString; typedef std::string FileNameType; typedef std::vector PolygonType; typedef std::vector PolygonSequenceType; typedef double IndexValueType; typedef double DoseStatisticType; typedef std::string DBStringType; typedef std::string VirtuosFileNameString, DICOMRTFileNameString; typedef unsigned short Uint16; struct Area2D { WorldCoordinate x_begin; WorldCoordinate x_end; WorldCoordinate y_begin; WorldCoordinate y_end; VoxelGridIndex2D index_begin; VoxelGridIndex2D index_end; void Init() { x_begin = -1000000; x_end = -1000000; y_begin = -1000000; y_end = -1000000; index_begin(0) = 0; index_begin(1) = 0; index_end(0) = 0; index_end(1) = 0; } }; struct Segment2D { WorldCoordinate2D begin; WorldCoordinate2D end; }; struct Segment3D { WorldCoordinate3D begin; WorldCoordinate3D end; }; typedef int DistributionType; typedef std::string PathType; typedef std::string XMLString, StatisticsString; }//end: namespace rttb #endif diff --git a/code/core/rttbDVH.cpp b/code/core/rttbDVH.cpp index 74c0ac0..f2d8157 100644 --- a/code/core/rttbDVH.cpp +++ b/code/core/rttbDVH.cpp @@ -1,418 +1,418 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include "rttbDVH.h" #include "rttbException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace core { DVH::~DVH() {} DVH::DVH(const DataDifferentialType& aDataDifferential, const DoseTypeGy& aDeltaD, const DoseVoxelVolumeType& aDeltaV, IDType aStructureID, IDType aDoseID): _deltaD(aDeltaD), _deltaV(aDeltaV), _structureID(aStructureID), _doseID(aDoseID), _voxelizationID("NONE") { _dataDifferential.clear(); _dataDifferential = aDataDifferential; this->init(); } DVH::DVH(const DataDifferentialType& aDataDifferential, DoseTypeGy aDeltaD, DoseVoxelVolumeType aDeltaV, IDType aStructureID, IDType aDoseID, IDType aVoxelizationID): _deltaD(aDeltaD), _deltaV(aDeltaV), _structureID(aStructureID), _doseID(aDoseID), _voxelizationID(aVoxelizationID) { _dataDifferential.clear(); _dataDifferential = aDataDifferential; this->init(); } DVH::DVH(const DVH& copy) : _structureID(copy._structureID), _doseID(copy._doseID), _voxelizationID(copy._voxelizationID), _label(copy._label) { _deltaD = copy._deltaD; _deltaV = copy._deltaV; _dataDifferential.clear(); _dataDifferential = copy._dataDifferential; this->init(); } 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) { if (aDVH.getStructureID() != otherDVH.getStructureID()) { return false; } if (aDVH.getDoseID() != otherDVH.getDoseID()) { return false; } if (aDVH.getVoxelizationID() != otherDVH.getVoxelizationID()) { return false; } if (aDVH.getNumberOfVoxels() != otherDVH.getNumberOfVoxels()) { return false; } return true; } std::ostream& operator<<(std::ostream& s, const DVH& aDVH) { s << "[ " << aDVH.getStructureID() << ", " << aDVH.getDoseID() << ", " << aDVH.getVoxelizationID() << ", " << aDVH.getNumberOfVoxels() << " ]"; return s; } std::deque DVH::getDataDifferential(bool relativeVolume) const { if (!relativeVolume) { return _dataDifferential; } else { return _dataDifferentialRelative; } } 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(); this->_dataCumulativeRelative.clear(); this->_dataDifferentialRelative.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; for (it = _dataDifferential.begin(); it != _dataDifferential.end(); ++it) { DoseCalcType datai = ((*it) * 1.0 / _numberOfVoxels); _dataDifferentialRelative.push_back(datai); } _variance = (squareSum / _numberOfVoxels - _mean * _mean); _stdDeviation = pow(_variance, 0.5); _dataCumulative = this->calcCumulativeDVH(); } std::deque DVH::calcCumulativeDVH(bool relativeVolume) { _dataCumulative.clear(); _dataCumulativeRelative.clear(); DoseCalcType cumulativeDVHi = 0; size_t size = _dataDifferential.size(); - for (int i = 0; i < size; i++) + for (size_t i = 0; i < size; i++) { cumulativeDVHi += _dataDifferential.at(size - i - 1); if (!relativeVolume) { _dataCumulative.push_front(cumulativeDVHi); } else { _dataCumulativeRelative.push_front(cumulativeDVHi / this->getNumberOfVoxels()); } } if (!relativeVolume) { return _dataCumulative; } else { return _dataCumulativeRelative; } } 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) { GridIndexType i = static_cast(xDoseAbsolute / _deltaD); if (i < _dataCumulative.size()) { VolumeType vx = (_dataCumulative.at(i)); vx = (vx * this->_deltaV); return vx; } else if (i < _dataCumulativeRelative.size()) { VolumeType vx = (_dataCumulativeRelative.at(i)); vx = (vx * this->_deltaV); return vx; } else { return 0; } } DoseTypeGy DVH::getDx(VolumeType xVolumeAbsolute) { GridIndexType i = 0; if (!_dataCumulative.empty()) { for (; i < _dataCumulative.size(); i++) { double volumeAbsoluteI = _dataCumulative[i] * this->_deltaV; if (xVolumeAbsolute > volumeAbsoluteI) { break; } } } else { for (; i < _dataCumulativeRelative.size(); i++) { double volumeAbsoluteI = _dataCumulativeRelative[i] * this->_deltaV; if (xVolumeAbsolute / this->getNumberOfVoxels() > volumeAbsoluteI) { break; } } } if (i <= _dataCumulative.size() && i > 0) { DoseTypeGy dx = (i - 1) * this->_deltaD; return dx; } else if (i < _dataCumulativeRelative.size() && i > 0) { DoseTypeGy dx = (i - 1) * this->_deltaD; return dx; } else { return 0; } } VolumeType DVH::getAbsoluteVolume(int relativePercent) { return (relativePercent * getNumberOfVoxels() * getDeltaV() / 100.0); } void DVH::setLabel(StructureLabel aLabel) { _label = aLabel; } StructureLabel DVH::getLabel() const { return _label; } }//end namespace core }//end namespace rttb diff --git a/code/core/rttbDVH.h b/code/core/rttbDVH.h index 5c49c66..f9d9a82 100644 --- a/code/core/rttbDVH.h +++ b/code/core/rttbDVH.h @@ -1,191 +1,194 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __DVH_H #define __DVH_H #include #include #include "rttbBaseType.h" #include "rttbStructure.h" namespace rttb { namespace core { class Structure; /*! @class DVH @brief This is a class representing a dose volume histogram (DVH) */ class DVH { public: typedef std::deque DataDifferentialType; typedef boost::shared_ptr DVHPointer; private: - IDType _structureID; - IDType _doseID; - IDType _voxelizationID; /*! @brief Differential dvh data index is the dose bin, value is the voxel number (sub voxel accuracy) of the dose bin */ DataDifferentialType _dataDifferential; /*! @brief Differential dvh data relative to the total number of voxels */ DataDifferentialType _dataDifferentialRelative; - /*! @brief Volume of a voxel in cm3 - */ - DoseVoxelVolumeType _deltaV; /*! @brief Absolute dose value of a dose-bin in Gy */ DoseTypeGy _deltaD; + /*! @brief Volume of a voxel in cm3 + */ + DoseVoxelVolumeType _deltaV; + + IDType _structureID; + IDType _doseID; + IDType _voxelizationID; + + StructureLabel _label; DoseStatisticType _maximum; DoseStatisticType _minimum; DoseStatisticType _mean; DoseStatisticType _modal; DVHVoxelNumber _numberOfVoxels; DoseStatisticType _median; DoseStatisticType _stdDeviation; DoseStatisticType _variance; DataDifferentialType _dataCumulative; DataDifferentialType _dataCumulativeRelative; /*! @brief DVH initialisation The DVH is initialized and all statistical values are calculated. @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ void init(); public: ~DVH(); /*! @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ DVH(const DataDifferentialType& aDataDifferential, const DoseTypeGy& aDeltaD, const DoseVoxelVolumeType& aDeltaV, IDType aStructureID, IDType aDoseID); /*! @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ DVH(const DataDifferentialType& aDataDifferential, DoseTypeGy aDeltaD, DoseVoxelVolumeType aDeltaV, IDType aStructureID, IDType aDoseID, IDType aVoxelizationID); DVH(const DVH& copy); /*! @throw if _deltaV or _deltaD are zero @throw is _data differential is empty */ DVH& operator=(const DVH& copy); /*! equality operator DVHs are considered equal if they have the same structure, dose and voxelization ID and the same number of voxels. */ bool friend operator==(const DVH& aDVH, const DVH& otherDVH); friend std::ostream& operator<<(std::ostream& s, const DVH& aDVH); void setLabel(StructureLabel aLabel); StructureLabel getLabel() const; /*! @param relativeVolume default false-> Value is the voxel number of the dose bin; if true-> value is the relative volume % between 0 and 1, (the voxel number of this dose bin)/(number of voxels) @return Return differential data of the dvh (relative or abolute depending on the input parameter). */ DataDifferentialType getDataDifferential(bool relativeVolume = false) const; DoseVoxelVolumeType getDeltaV() const; DoseTypeGy getDeltaD() const; IDType getStructureID() const; IDType getDoseID() const; IDType getVoxelizationID() const; void setDoseID(IDType aDoseID); void setStructureID(IDType aStrID); /*! @brief Calculate number of the voxels (with sub voxel accuracy) @return Return -1 if not initialized */ DVHVoxelNumber getNumberOfVoxels() const; /*! @brief Get the maximum dose in Gy from dvh @return Return the maximum dose in Gy (i+0.5)*deltaD, i-the maximal dose-bin with volume>0 Return -1 if not initialized */ DoseStatisticType getMaximum() const; /*! @brief Get the minimum dose in Gy from dvh @return Return the minimum dose (i+0.5)*deltaD, i-the minimal dose-bin with volume>0 Return -1 if not initialized */ DoseStatisticType getMinimum() const; DoseStatisticType getMean() const; DoseStatisticType getMedian() const; DoseStatisticType getModal() const; DoseStatisticType getStdDeviation() const; DoseStatisticType getVariance() const; /*! @brief Calculate the cumulative data of dvh @param relativeVolume default false-> Value is the voxel number of the dose bin; if true-> value is the relative volume % between 0 and 1, (the voxel number of this dose bin)/(number of voxels) @return Return cumulative dvh value i-voxel number in dose-bin i */ DataDifferentialType calcCumulativeDVH(bool relativeVolume = false); /*! @brief Get Vx the volume irradiated to >= x @return Return absolute Volume in absolute cm3 Return -1 if not initialized */ VolumeType getVx(DoseTypeGy xDoseAbsolute); /*! @brief Get Dx the minimal dose delivered to x @return Return absolute dose value in Gy Return -1 if not initialized */ DoseTypeGy getDx(VolumeType xVolumeAbsolute); /*! @brief Calculate the absolute volume in cm3 @param relativePercent 0~100, the percent of the whole volume */ VolumeType getAbsoluteVolume(int relativePercent); }; } } #endif diff --git a/code/core/rttbDVHSet.cpp b/code/core/rttbDVHSet.cpp index 2410755..6eac447 100644 --- a/code/core/rttbDVHSet.cpp +++ b/code/core/rttbDVHSet.cpp @@ -1,223 +1,223 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "rttbDVHSet.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace core { DVHSet::DVHSet(IDType aStructureSetID, IDType aDoseID) : _structureSetID(aStructureSetID), _doseID(aDoseID) {} DVHSet::DVHSet(DVHSetType aDVHTVSet, DVHSetType aDVHHTSet, IDType aStructureSetID, IDType aDoseID) { _dvhHTSet = aDVHHTSet; _dvhTVSet = aDVHTVSet; this->_structureSetID = aStructureSetID; this->_doseID = aDoseID; } DVHSet::DVHSet(DVHSetType aDVHTVSet, DVHSetType aDVHHTSet, DVHSetType aDVHWVSet, IDType aStructureSetID, IDType aDoseID) { _dvhHTSet = aDVHHTSet; _dvhTVSet = aDVHTVSet; _dvhWVSet = aDVHWVSet; this->_structureSetID = aStructureSetID; this->_doseID = aDoseID; } std::size_t DVHSet::size() const { return _dvhHTSet.size() + _dvhTVSet.size() + _dvhWVSet.size(); } void DVHSet::setStrSetID(IDType aStrSetID) { _structureSetID = aStrSetID; } void DVHSet::setDoseID(IDType aDoseID) { _doseID = aDoseID; } IDType DVHSet::getStrSetID() const { return _structureSetID; } IDType DVHSet::getDoseID() const { return _doseID; } DVH* DVHSet::getDVH(IDType structureID) { DVHSetType::iterator itTV = _dvhTVSet.begin(); for (; itTV != _dvhTVSet.end(); ++itTV) { if ((*(itTV)).getStructureID() == structureID) { return &(*itTV); } } DVHSetType::iterator itHT = _dvhHTSet.begin(); for (; itHT != _dvhHTSet.end(); ++itHT) { if ((*(itHT)).getStructureID() == structureID) { return &(*itHT); } } DVHSetType::iterator itWV = _dvhWVSet.begin(); for (; itWV != _dvhWVSet.end(); ++itWV) { if ((*(itWV)).getStructureID() == structureID) { return &(*itWV); } } std::cout << "No DVH with the structure id: " << structureID << " was found!" << std::endl; return NULL; } void DVHSet::insert(DVH& aDvh, DVHRole aDVHRole) { if (aDVHRole.Type == DVHRole::TargetVolume) { _dvhTVSet.push_back(aDvh); } else if (aDVHRole.Type == DVHRole::HealthyTissue) { _dvhHTSet.push_back(aDvh); } else if (aDVHRole.Type == DVHRole::WholeVolume) { _dvhWVSet.push_back(aDvh); } else { throw core::InvalidParameterException("aDVHType must be TV or HT or WV!"); } } VolumeType DVHSet::getWholeVolume(DoseTypeGy aDoseAbsolute) const { VolumeType volume = this->getHealthyTissueVolume(aDoseAbsolute) + this->getTargetVolume( aDoseAbsolute); return volume; } VolumeType DVHSet::getHealthyTissueVolume(DoseTypeGy aDoseAbsolute) const { DVHSetType::const_iterator itHT = _dvhHTSet.begin(); VolumeType volume = 0; while (itHT != _dvhHTSet.end()) { VolumeType testVol = 0; DVH dvh = *(itHT); testVol = dvh.getVx(aDoseAbsolute); if (testVol >= 0) { volume += testVol; } ++itHT; } return volume; } VolumeType DVHSet::getTargetVolume(DoseTypeGy aDoseAbsolute) const { DVHSetType::const_iterator itTV = _dvhTVSet.begin(); VolumeType volume = 0; while (itTV != _dvhTVSet.end()) { VolumeType testVol = 0; DVH dvh = *(itTV); testVol = dvh.getVx(aDoseAbsolute); if (testVol >= 0) { volume += testVol; } ++itTV; } return volume; } bool operator==(const DVHSet& aDVHSet, const DVHSet& otherDVHSet) { if (aDVHSet.getStrSetID() != otherDVHSet.getStrSetID()) { return false; } if (aDVHSet.getDoseID() != otherDVHSet.getDoseID()) { return false; } if (aDVHSet.size() != otherDVHSet.size()) { return false; } return true; } std::ostream& operator<<(std::ostream& s, const DVHSet& aDVHSet) { s << "[ " << aDVHSet.getStrSetID() << ", " << aDVHSet.getDoseID() << " ]"; return s; } std::ostream& operator<<(std::ostream& s, const DVHSet::DVHSetType& aDVHSet) { s << "[ "; - for (int i = 0; i < aDVHSet.size(); i++) + for (size_t i = 0; i < aDVHSet.size(); i++) { s << aDVHSet.at(i); } s << " ]"; return s; } } } diff --git a/code/core/rttbGeometricInfo.h b/code/core/rttbGeometricInfo.h index 2f1f6e8..3e457c5 100644 --- a/code/core/rttbGeometricInfo.h +++ b/code/core/rttbGeometricInfo.h @@ -1,190 +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. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __GEOMETRIC_INFO_NEW_H #define __GEOMETRIC_INFO_NEW_H #include #include #include #include #include "rttbBaseType.h" namespace rttb { namespace core { /*! @brief GeometricInfo objects contain all the information required for transformations between voxel grid coordinates and world coordinates. Corresponding converter functions are also available. @note ITK Pixel Indexing used (http://www.itk.org/Doxygen45/html/classitk_1_1Image.html): The Index type reverses the order so that with Index[0] = col, Index[1] = row, Index[2] = slice. */ class GeometricInfo { private: WorldCoordinate3D _imagePositionPatient; OrientationMatrix _orientationMatrix; OrientationMatrix _invertedOrientationMatrix; SpacingVectorType3D _spacing; VoxelGridDimensionType _numberOfColumns; VoxelGridDimensionType _numberOfRows; VoxelGridDimensionType _numberOfFrames; /* @brief Matrix inversion routine. Uses lu_factorize and lu_substitute in uBLAS to invert a matrix http://savingyoutime.wordpress.com/2009/09/21/c-matrix-inversion-boostublas/ */ bool computeInvertOrientation(); public: /*! @brief Constructor, initializes orientation matrix, spacing vector and patient position with zeros. */ - GeometricInfo() : _orientationMatrix(0), _spacing(0), _imagePositionPatient(0), _numberOfFrames(0), - _numberOfRows(0), - _numberOfColumns(0) {} + GeometricInfo() : _imagePositionPatient(0),_orientationMatrix(0), _spacing(0), _numberOfColumns(0), + _numberOfRows(0), _numberOfFrames(0) + {} void setSpacing(const SpacingVectorType3D& aSpacingVector); const SpacingVectorType3D& getSpacing() const; void setImagePositionPatient(const WorldCoordinate3D& aImagePositionPatient); const WorldCoordinate3D& getImagePositionPatient() const; void setOrientationMatrix(const OrientationMatrix& anOrientationMatrix); const OrientationMatrix getOrientationMatrix() const { return _orientationMatrix; }; /*! @brief Returns the 1st row of the OrientationMatrix. @deprecated please use getOrientationMatrix() (x,0) instead*/ const ImageOrientation getImageOrientationRow() const; /*! @brief Returns the 2nd row of the OrientationMatrix. @deprecated please use getOrientationMatrix() (x,1) instead*/ const ImageOrientation getImageOrientationColumn() const; void setPixelSpacingRow(const GridVolumeType aValue); /*! @brief Returns the spacing in the x-dimension (rows) of the data grid. @deprecated please use getSpacing() (0) instead*/ const GridVolumeType getPixelSpacingRow() const; void setPixelSpacingColumn(const GridVolumeType aValue); /*! @brief Returns the spacing in the y-dimension (columns) of the data grid. @deprecated please use getSpacing() (1) instead*/ const GridVolumeType getPixelSpacingColumn() const; void setSliceThickness(const GridVolumeType aValue); /*! @brief Returns the spacing in the z-dimension (slices) of the data grid. @deprecated please use getSpacing() (2) instead*/ const GridVolumeType getSliceThickness() const; void setImageSize(const ImageSize& aSize); const ImageSize getImageSize() const; void setNumColumns(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumColumns() const; void setNumRows(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumRows() const; void setNumSlices(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumSlices() const; /*! @brief determines equality of two GeometricInfo objects. */ friend bool operator==(const GeometricInfo& gInfo, const GeometricInfo& gInfo1); bool equalsAlmost(const GeometricInfo& another, double errorConstant = 1e-5) const; /*! @brief convert world coordinates to voxel grid index. The conversion of values is done even if the target index is not inside the given voxel grid (return false). If the target is inside the grid return true. */ bool worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, VoxelGridIndex3D& aIndex) const; /*! @brief convert world coordinates to double geometry coordinate. The world coordinate of the image position patient (center of the first voxel) will be convert to the double voxel (0.0, 0.0, 0.0) The conversion of values is done even if the target index is not inside the given voxel grid (return false). If the target is inside the grid return true. */ bool worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, DoubleVoxelGridIndex3D& aIndex) const; /*! @brief convert double geometry coordinate to world coordinates. The double voxel index (0.0, 0.0, 0.0) will be convert to the world coordinate of the image postion patient (center of the first voxel) The conversion of values is done even if the target is not inside the given voxel grid (return false). If the target is inside the voxel grid return true. */ bool geometryCoordinateToWorldCoordinate(const DoubleVoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const; /*! @brief convert int voxel grid index to world coordinates. It is the upper left corner of the voxel. The conversion of values is done even if the target is not inside the given voxel grid (return false). If the target is inside the voxel grid return true. */ bool indexToWorldCoordinate(const VoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const; /*! @brief check if a given voxel grid index is inside the given voxel grid.*/ bool isInside(const VoxelGridIndex3D& aIndex) const; /*! @brief check if a given world coordinate is inside the given voxel grid.*/ bool isInside(const WorldCoordinate3D& aWorldCoordinate); const GridSizeType getNumberOfVoxels() const; bool convert(const VoxelGridID& gridID, VoxelGridIndex3D& gridIndex) const; bool convert(const VoxelGridIndex3D& gridIndex, VoxelGridID& gridID) const; /*! @brief test if given ID is inside current dose grid */ bool validID(const VoxelGridID aID) const; /*! @brief test if given index is inside current dose grid */ bool validIndex(const VoxelGridIndex3D& aIndex) const; /*!@ brief generates string stream representation of the GeometricInfo object. */ friend std::ostream& operator<<(std::ostream& s, const GeometricInfo& anGeometricInfo); }; } } -#endif \ No newline at end of file +#endif diff --git a/code/indices/rttbCoverageIndex.cpp b/code/indices/rttbCoverageIndex.cpp index 426086f..f79ec24 100644 --- a/code/indices/rttbCoverageIndex.cpp +++ b/code/indices/rttbCoverageIndex.cpp @@ -1,80 +1,79 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "rttbCoverageIndex.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace indices { CoverageIndex::CoverageIndex(DVHSetPtr 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(); - VolumeType Vref = _dvhSet->getWholeVolume(_doseReference); 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/rttbHomogeneityIndex.cpp b/code/indices/rttbHomogeneityIndex.cpp index 556c24d..fd6ec79 100644 --- a/code/indices/rttbHomogeneityIndex.cpp +++ b/code/indices/rttbHomogeneityIndex.cpp @@ -1,102 +1,103 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ +#include #include "rttbHomogeneityIndex.h" #include "rttbInvalidParameterException.h" #include "rttbExceptionMacros.h" namespace rttb { namespace indices { HomogeneityIndex::HomogeneityIndex(DVHSetPtr dvhSet, DoseTypeGy aDoseReference) : DvhBasedDoseIndex(dvhSet, aDoseReference) { init(); } bool HomogeneityIndex::calcIndex() { double max = 0; - double min; + double min = std::numeric_limits::max(); std::vector dvhTVSet = this->_dvhSet->getTargetVolumeSet(); std::vector::iterator it; for (it = dvhTVSet.begin(); it != dvhTVSet.end(); ++it) { core::DVH dvh = *(it); if (it == dvhTVSet.begin()) { min = dvh.getMinimum(); } if (dvh.getMaximum() > max) { max = dvh.getMaximum(); } if (dvh.getMinimum() < min) { min = dvh.getMinimum(); } } if (this->getDoseReference() != 0) { _value = (max - min) / this->getDoseReference(); } else { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } return true; } IndexValueType HomogeneityIndex::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); if (this->getDoseReference() <= 0) { rttbExceptionMacro(core::InvalidParameterException, << "Reference dose " << this->getDoseReference() << " invalid: Volume of reference dose should not be 0!"); } return (dvh.getMaximum() - dvh.getMinimum()) / this->getDoseReference(); } }//end namespace indices -}//end namespace rttb \ No newline at end of file +}//end namespace rttb diff --git a/code/interpolation/rttbInterpolationBase.cpp b/code/interpolation/rttbInterpolationBase.cpp index f255271..d847323 100644 --- a/code/interpolation/rttbInterpolationBase.cpp +++ b/code/interpolation/rttbInterpolationBase.cpp @@ -1,196 +1,196 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include "rttbInterpolationBase.h" #include "rttbInvalidParameterException.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace interpolation { void InterpolationBase::setAccessorPointer(const AccessorPointer originalData) { if (originalData != NULL) { _spOriginalData = originalData; } else { throw core::NullPointerException("originalDose is NULL!"); } }; void InterpolationBase::getNeighborhoodVoxelValues( const WorldCoordinate3D& aWorldCoordinate, unsigned int neighborhood, boost::array& target, boost::shared_ptr values) const { if (_spOriginalData == NULL) { throw core::NullPointerException("originalDose is NULL!"); } //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 (int i = 0; i < 3; i++) + 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: this is a workaround, not a good solution 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 (int zIncr = 0; zIncr < 2; zIncr++) + for (unsigned int zIncr = 0; zIncr < 2; zIncr++) { - for (int yIncr = 0; yIncr < 2; yIncr++) + for (unsigned int yIncr = 0; yIncr < 2; yIncr++) { - for (int xIncr = 0; xIncr < 2; xIncr++) + for (unsigned int xIncr = 0; xIncr < 2; xIncr++) { cornerPoints.push_back(VoxelGridIndex3D(leftTopFrontCoordinate[0] + xIncr, leftTopFrontCoordinate[1] + yIncr, leftTopFrontCoordinate[2] + zIncr)); } } } //target range has to be always [0,1] - for (int i = 0; i < 3; i++) + 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}; - int runningIndex; + 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; - int replacementVoxelIndex = 0; + 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/io/dicom/rttbDicomIODStructureSetGenerator.cpp b/code/io/dicom/rttbDicomIODStructureSetGenerator.cpp index a1c028c..5172949 100644 --- a/code/io/dicom/rttbDicomIODStructureSetGenerator.cpp +++ b/code/io/dicom/rttbDicomIODStructureSetGenerator.cpp @@ -1,198 +1,198 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include #include #include #include #include #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbStructure.h" #include "rttbDicomIODStructureSetGenerator.h" #include "rttbDcmrtException.h" namespace rttb { - namespace io - { - namespace dicom - { - DicomIODStructureSetGenerator::DicomIODStructureSetGenerator(DRTStrSetIODPtr aDRTStructureSetIOD) - { - _drtStrSetIOD = aDRTStructureSetIOD; - } + namespace io + { + namespace dicom + { + DicomIODStructureSetGenerator::DicomIODStructureSetGenerator(DRTStrSetIODPtr aDRTStructureSetIOD) + { + _drtStrSetIOD = aDRTStructureSetIOD; + } - void DicomIODStructureSetGenerator::readStrSet() - { - OFString uid; - _drtStrSetIOD->getSeriesInstanceUID(uid); - _UID = uid.c_str(); + void DicomIODStructureSetGenerator::readStrSet() + { + OFString uid; + _drtStrSetIOD->getSeriesInstanceUID(uid); + _UID = uid.c_str(); - OFString uid2; - _drtStrSetIOD->getPatientID(uid2); - _patientUID = uid2.c_str(); + OFString uid2; + _drtStrSetIOD->getPatientID(uid2); + _patientUID = uid2.c_str(); - DRTStructureSetROISequence* rois = &_drtStrSetIOD->getStructureSetROISequence(); + DRTStructureSetROISequence* rois = &_drtStrSetIOD->getStructureSetROISequence(); //generate map of relevant ROIs std::map filteredROIs; ::boost::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()); if (!this->getStructureLabelFilterActive() || boost::regex_match(roiName, e)) { filteredROIs.insert(std::make_pair(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; + /*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(); + long numberOfStructures = rcs->getNumberOfItems(); + bool isEmpty = rcs->isEmpty(); - if (numberOfStructures == 0 || isEmpty) - { - throw core::InvalidParameterException("Empty Structure Set!") ; - } + if (numberOfStructures == 0 || isEmpty) + { + throw core::InvalidParameterException("Empty Structure Set!") ; + } - int structureNo = 0; + int structureNo = 0; - for (rcs->gotoFirstItem(); (rcs->getCurrentItem(rcsItem)).good(); rcs->gotoNextItem()) - { - OFString refROINumber; - rcsItem->getReferencedROINumber(refROINumber); + 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(); - long no2 = cs->getNumberOfItems(); - - PolygonSequenceType structureVector; - - for (int 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); - - int numberOfContourPointsInt; - std::stringstream is(numberOfContourPoints.c_str()); - is >> numberOfContourPointsInt; - OFString countourNumber; - csItem->getContourNumber(countourNumber); - - PolygonType contourVector; - char* pEnd; - - for (int k = 0; k < numberOfContourPointsInt; k++) - { - WorldCoordinate3D point; - - for (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); + 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() << std::endl; - std::stringstream sstr; - sstr << structureNo; - str->setUID(sstr.str()); + std::stringstream sstr; + sstr << structureNo; + str->setUID(sstr.str()); - _strVector.push_back(str); - } + _strVector.push_back(str); + } ++structureNo; - } - } - - DicomIODStructureSetGenerator::~DicomIODStructureSetGenerator() - { + } } - DicomIODStructureSetGenerator::StructureSetPointer - DicomIODStructureSetGenerator::generateStructureSet() - { - this->readStrSet(); - return boost::make_shared(_strVector, _patientUID, _UID); - } - }//end namespace dicom - }//end namespace io + DicomIODStructureSetGenerator::~DicomIODStructureSetGenerator() + { + } + + DicomIODStructureSetGenerator::StructureSetPointer + DicomIODStructureSetGenerator::generateStructureSet() + { + this->readStrSet(); + return boost::make_shared(_strVector, _patientUID, _UID); + } + }//end namespace dicom + }//end namespace io }//end namespace rttb diff --git a/code/io/helax/rttbDicomHelaxDoseAccessor.cpp b/code/io/helax/rttbDicomHelaxDoseAccessor.cpp index 13c1015..0157587 100644 --- a/code/io/helax/rttbDicomHelaxDoseAccessor.cpp +++ b/code/io/helax/rttbDicomHelaxDoseAccessor.cpp @@ -1,329 +1,329 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include "boost/filesystem/operations.hpp" #include "boost/filesystem/path.hpp" #include "boost/progress.hpp" #include #include "rttbDicomHelaxDoseAccessor.h" #include "rttbInvalidDoseException.h" #include "rttbDcmrtException.h" #include "rttbIndexOutOfBoundsException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace io { namespace helax { DicomHelaxDoseAccessor::~DicomHelaxDoseAccessor() { } DicomHelaxDoseAccessor::DicomHelaxDoseAccessor(std::vector aDICOMRTDoseVector) { - for (int i = 0; i < aDICOMRTDoseVector.size(); i++) + for (size_t i = 0; i < aDICOMRTDoseVector.size(); i++) { _doseVector.push_back(aDICOMRTDoseVector.at(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 (int i = 0; i < _doseVector.size(); i++) + for (size_t i = 0; i < _doseVector.size(); i++) { DRTDoseIODPtr dose = _doseVector.at(i); 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 i = 0; - i < static_cast(_geoInfo.getNumColumns()*_geoInfo.getNumRows()); i++) + for (unsigned int j = 0; + j < static_cast(_geoInfo.getNumColumns()*_geoInfo.getNumRows()); j++) { - Uint16 data = static_cast(pixelData[i] * currentDoseGridScaling / + Uint16 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; } bool 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 ImageOrientation 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.getPixelSpacingRow() == 0 || _geoInfo.getPixelSpacingColumn() == 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.getPixelSpacingRow() << "!" << std::endl; spacingVector(2) = spacingVector(0); } } _geoInfo.setSpacing(spacingVector); return true; } 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 228dc8d..1882187 100644 --- a/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp +++ b/code/io/helax/rttbDicomHelaxFileDoseAccessorGenerator.cpp @@ -1,94 +1,94 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include #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() {} 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 (int i = 0; i < fileVector.size(); i++) + for (size_t i = 0; i < fileVector.size(); i++) { DRTDoseIODPtr dose = boost::make_shared(); status = fileformat.loadFile(fileVector.at(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/rttbITKImageAccessor.cpp b/code/io/itk/rttbITKImageAccessor.cpp index 92c6fa1..f455ae0 100644 --- a/code/io/itk/rttbITKImageAccessor.cpp +++ b/code/io/itk/rttbITKImageAccessor.cpp @@ -1,118 +1,118 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include "rttbITKImageAccessor.h" #include "rttbException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace io { namespace itk { ITKImageAccessor::~ITKImageAccessor() { } 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; } } bool 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 (int col = 0; col < 3; ++col) + for (unsigned int col = 0; col < 3; ++col) { - for (int row = 0; row < 3; ++row) + 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!") ; } return true; } } } } diff --git a/code/io/itk/rttbITKImageAccessorConverter.cpp b/code/io/itk/rttbITKImageAccessorConverter.cpp index c6de577..0559cc2 100644 --- a/code/io/itk/rttbITKImageAccessorConverter.cpp +++ b/code/io/itk/rttbITKImageAccessorConverter.cpp @@ -1,111 +1,111 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include "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 (int col = 0; col < 3; ++col) + for (unsigned int col = 0; col < 3; ++col) { - for (int row = 0; row < 3; ++row) + 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/rttbITKImageMaskAccessor.cpp b/code/io/itk/rttbITKImageMaskAccessor.cpp index 0a7f4c2..129616c 100644 --- a/code/io/itk/rttbITKImageMaskAccessor.cpp +++ b/code/io/itk/rttbITKImageMaskAccessor.cpp @@ -1,176 +1,176 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include "rttbITKImageMaskAccessor.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() { }; 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 (int col = 0; col < 3; ++col) + for (unsigned int col = 0; col < 3; ++col) { - for (int row = 0; row < 3; ++row) + 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) { MaskVoxelListPointer filteredVoxelVectorPointer(new MaskVoxelList); updateMask(); - int size = _geoInfo->getNumColumns() * _geoInfo->getNumRows() * _geoInfo->getNumSlices(); + unsigned int size = _geoInfo->getNumColumns() * _geoInfo->getNumRows() * _geoInfo->getNumSlices(); filteredVoxelVectorPointer->reserve(size); - for (int gridIndex = 0 ; gridIndex < size; gridIndex++) + 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 1400def..420c37c 100644 --- a/code/io/itk/rttbITKImageMaskAccessorConverter.cpp +++ b/code/io/itk/rttbITKImageMaskAccessorConverter.cpp @@ -1,129 +1,129 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include "rttbITKImageMaskAccessorConverter.h" #include "rttbInvalidDoseException.h" #include "rttbGeometricInfo.h" #include "itkImageRegionIterator.h" #include "rttbITKImageMaskAccessor.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::RegionType region; 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 (int col = 0; col < 3; ++col) + for (unsigned int col = 0; col < 3; ++col) { - for (int row = 0; row < 3; ++row) + 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 = ITKImageMaskAccessor::ITKMaskImageType::New(); _itkImage->SetRegions(region); _itkImage->SetSpacing(spacing); _itkImage->SetDirection(direction); _itkImage->SetOrigin(origin); _itkImage->Allocate(); ::itk::ImageRegionIterator imageIterator(_itkImage, region); VoxelGridID id = 0; //Transfer Mask values to itk image //Large advantage: rttbVoxelGridId ordering is the same as itk ordering while (!imageIterator.IsAtEnd()) { VoxelGridIndex3D aIndex; if (_maskAccessor->getGeometricInfo().validID(id)) { rttb::core::MaskVoxel maskVoxel = core::MaskVoxel(id);; _maskAccessor->getMaskAt(id, maskVoxel); // Set the current pixel imageIterator.Set(maskVoxel.getRelevantVolumeFraction()); } else { if (failsOnInvalidIDs()) { throw core::InvalidDoseException("invalid Mask id!"); return false; } else { imageIterator.Set(_invalidDoseValue); } } ++imageIterator; ++id; } return true; } }//end namespace mask }//end namespace io }//end namespace rttb diff --git a/code/io/other/rttbDVHTxtFileWriter.cpp b/code/io/other/rttbDVHTxtFileWriter.cpp index 87651bd..3bf0254 100644 --- a/code/io/other/rttbDVHTxtFileWriter.cpp +++ b/code/io/other/rttbDVHTxtFileWriter.cpp @@ -1,128 +1,126 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include #include "rttbDVHTxtFileWriter.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbException.h" namespace rttb { namespace io { namespace other { DVHTxtFileWriter::DVHTxtFileWriter(FileNameString aFileName, DVHType aDVHType) { this->setFileName(aFileName); this->setDVHType(aDVHType); } void DVHTxtFileWriter::setDVHType(DVHType aDVHType) { _dvhType = aDVHType; } FileNameString DVHTxtFileWriter::getFileName() const { return _fileName; } void DVHTxtFileWriter::setFileName(FileNameString aFileName) { _fileName = aFileName; } DVHType DVHTxtFileWriter::getDVHType() const { return _dvhType; } void DVHTxtFileWriter::writeDVH(DVHPointer aDvh) { if (!aDvh) { throw core::NullPointerException("aDvh must not be NULL! "); } std::ofstream out_dvh_ofstream(this->_fileName.c_str(), std::ios::out); if (!out_dvh_ofstream.is_open() || !out_dvh_ofstream.good()) { throw core::InvalidParameterException("Invalid dvh file name: could not open write file"); } else { //setting string stream precission explicitly is mandatory to guarantee that tests //run sucessfully on different systems! out_dvh_ofstream.precision(10); if (_dvhType.Type != DVHType::Differential && _dvhType.Type != DVHType::Cumulative) { throw core::InvalidParameterException("DVH Type not acceptable: Only: DIFFERENTIAL/CUMULATIVE!"); } if (_dvhType.Type == DVHType::Differential) { out_dvh_ofstream << "DVH Type: DIFFERENTIAL\n"; } else if (_dvhType.Type == DVHType::Cumulative) { out_dvh_ofstream << "DVH Type: CUMULATIVE\n"; } out_dvh_ofstream << "DeltaD: " << aDvh->getDeltaD() << "\n"; out_dvh_ofstream << "DeltaV: " << aDvh->getDeltaV() << "\n"; out_dvh_ofstream << "StructureID: " << aDvh->getStructureID() << "\n"; out_dvh_ofstream << "DoseID: " << aDvh->getDoseID() << "\n"; out_dvh_ofstream << "DVH Data: " << "\n"; if (_dvhType.Type == DVHType::Differential) { DataDifferentialType dataDifferential = aDvh->getDataDifferential(); - unsigned int numberOfBins = static_cast(dataDifferential.size()); - for (unsigned int i = 0; i < numberOfBins; i++) + for (size_t i = 0; i < dataDifferential.size(); i++) { out_dvh_ofstream << i << "," << dataDifferential[i] << "\n"; } } else if (_dvhType.Type == DVHType::Cumulative) { DataDifferentialType dataCumulative = aDvh->calcCumulativeDVH(); - unsigned int numberOfBins = static_cast(dataCumulative.size()); - for (unsigned int i = 0; i < numberOfBins; i++) + for (size_t i = 0; i < dataCumulative.size(); i++) { out_dvh_ofstream << i << "," << dataCumulative[i] << "\n"; } } } } } } -} \ No newline at end of file +} diff --git a/code/io/other/rttbDVHXMLFileWriter.cpp b/code/io/other/rttbDVHXMLFileWriter.cpp index 4fbe3eb..f8f416d 100644 --- a/code/io/other/rttbDVHXMLFileWriter.cpp +++ b/code/io/other/rttbDVHXMLFileWriter.cpp @@ -1,140 +1,140 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "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(DVHPointer aDvh) { if (!aDvh) { throw core::NullPointerException("aDvh must not be NULL! "); } using boost::property_tree::ptree; ptree pt; if (_dvhType.Type == DVHType::Differential) { pt.put("dvh.type", "DIFFERENTIAL"); } else if (_dvhType.Type == DVHType::Cumulative) { pt.put("dvh.type", "CUMULATIVE"); } else { throw core::InvalidParameterException("DVH Type not acceptable: Only: DIFFERENTIAL/CUMULATIVE!"); } DataDifferentialType dataDifferential = aDvh->getDataDifferential(); size_t numberOfBins = dataDifferential.size(); pt.put("dvh.deltaD", aDvh->getDeltaD()); pt.put("dvh.deltaV", aDvh->getDeltaV()); pt.put("dvh.structureID", aDvh->getStructureID()); pt.put("dvh.doseID", aDvh->getDoseID()); if (_dvhType.Type == DVHType::Differential) { - for (int i = 0; i < numberOfBins; i++) + for (size_t i = 0; i < numberOfBins; i++) { boost::property_tree::ptree node; node.put("", dataDifferential[i]); node.put(".dosebin", i); pt.add_child("dvh.data.volume", node); } } else if (_dvhType.Type == DVHType::Cumulative) { DataDifferentialType dataCumulative = aDvh->calcCumulativeDVH(); - for (int i = 0; i < numberOfBins; i++) + for (size_t i = 0; i < numberOfBins; i++) { boost::property_tree::ptree node; node.put("", dataCumulative[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!"); } } } } -} \ No newline at end of file +} diff --git a/code/masks/boost/rttbBoostMaskAccessor.h b/code/masks/boost/rttbBoostMaskAccessor.h index 96f1f42..96cc8d6 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/boost/rttbBoostMaskAccessor.h @@ -1,128 +1,128 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __BOOST_MASK_ACCESSOR__H #define __BOOST_MASK_ACCESSOR__H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include "rttbGenericDoseIterator.h" #include "rttbStructure.h" #include namespace rttb { namespace masks { namespace boost { /*! @class BoostMaskAccessor * @brief Using the voxelization based on boost::geometry and generate the mask accessor. * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. * (A structure without self interseciton means all contours of the structure have no self intersection, and * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored * and the mask will be calculated, however, it may cause errors in the mask results. */ class BoostMaskAccessor: public core::MaskAccessorInterface { public: typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::Structure::StructTypePointer StructTypePointer; private: - + StructTypePointer _spStructure; core::GeometricInfo _geoInfo; /*! vector containing list of mask voxels*/ MaskVoxelListPointer _spRelevantVoxelVector; - StructTypePointer _spStructure; + IDType _maskUID; bool _strict; public: /*! @brief Constructor with a structure pointer and a geometric info pointer * @param aStructurePointer smart pointer of the structure * @param aGeometricInfoPtr smart pointer of the geometricinfo of the dose * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. * @exception InvalidParameterException thrown if strict is true and the structure has self intersections */ BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, bool strict = true); /*! @brief destructor*/ ~BoostMaskAccessor(); /*! @brief voxelization of the given structures using boost algorithms*/ void updateMask(); /*! @brief get vector containing all relevant voxels that are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(); /*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold); /*!@brief determine how a given voxel on the dose grid is masked * @param aID ID of the voxel in grid. * @param voxel Reference to the voxel. * @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined. * The relevant volume fraction will be set to zero. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const; /*!@brief determine how a given voxel on the dose grid is masked * @param aIndex 3d index of the voxel in grid. * @param voxel Reference to the voxel. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const; /*! @brief give access to GeometricInfo*/ const core::GeometricInfo& getGeometricInfo() const; /* @ brief is true if dose is on a homogeneous grid * @remark Inhomogeneous grids are not supported at the moment, but if they will be supported in the future the interface does not need to change.*/ bool isGridHomogeneous() const { return true; }; IDType getMaskUID() const { return _maskUID; }; }; } } } #endif diff --git a/code/masks/rttbVOIindexIdentifier.cpp b/code/masks/rttbVOIindexIdentifier.cpp index 5b9f874..be15378 100644 --- a/code/masks/rttbVOIindexIdentifier.cpp +++ b/code/masks/rttbVOIindexIdentifier.cpp @@ -1,103 +1,103 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: // @date $Date: // @author $Author: */ #include "rttbVOIindexIdentifier.h" #include #include "rttbExceptionMacros.h" namespace rttb { namespace masks { const std::vector VOIindexIdentifier::getIndicesByVoiRegex( StructSetTypePointer spStructSet, const std::string& nameAsRegEx) { if (!spStructSet) { rttbDefaultExceptionStaticMacro("spStructSet is NULL"); } std::vector voiLabelList; std::vector resultVOiIndices; boost::regex e(nameAsRegEx); - for (int i = 0; i < spStructSet->getNumberOfStructures(); i++) + for (size_t i = 0; i < spStructSet->getNumberOfStructures(); i++) { voiLabelList.push_back(spStructSet->getStructure(i)->getLabel()); std::string s = spStructSet->getStructure(i)->getLabel(); if (boost::regex_match(s, e)) { resultVOiIndices.push_back(i); } } return resultVOiIndices; } const unsigned int VOIindexIdentifier::getIndexByVoiName(StructSetTypePointer spStructSet, const std::string& name) { if (!spStructSet) { rttbDefaultExceptionStaticMacro("spStructSet is NULL"); } for (unsigned int i = 0; i < spStructSet->getNumberOfStructures(); i++) { if (spStructSet->getStructure(i)->getLabel() == name) { return i; } } rttbDefaultExceptionStaticMacro("no VOI was found with the given name"); } const std::string VOIindexIdentifier::getVoiNameByIndex(StructSetTypePointer spStructSet, unsigned int index) { if (!spStructSet) { rttbDefaultExceptionStaticMacro("spStructSet is NULL!"); } if (index >= spStructSet->getNumberOfStructures()) { rttbDefaultExceptionStaticMacro("invalid index, voiLabelList out of range"); } return spStructSet->getStructure(index)->getLabel(); } } } diff --git a/code/models/rttbIntegration.cpp b/code/models/rttbIntegration.cpp index 52e780f..f7cd95f 100644 --- a/code/models/rttbIntegration.cpp +++ b/code/models/rttbIntegration.cpp @@ -1,192 +1,192 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision $ (last changed revision) // @date $Date $ (last change date) // @author $Author $ (last changed by) */ #include #include #include "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 (int i = 0; i < tcp_params.volumeVector.size(); ++i) + 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; - int maxSteps = 50; + unsigned int maxSteps = 50; double eps = 1e-6; - int i = 1; + 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! "); } } -} \ No newline at end of file +} diff --git a/demoapps/BioModelCalc/BioModelCalc.cpp b/demoapps/BioModelCalc/BioModelCalc.cpp index b05f4ae..a51e823 100644 --- a/demoapps/BioModelCalc/BioModelCalc.cpp +++ b/demoapps/BioModelCalc/BioModelCalc.cpp @@ -1,133 +1,133 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "BioModelCalcApplicationData.h" #include "BioModelCalcHelper.h" #include "BioModelCmdLineParser.h" #include "boost/shared_ptr.hpp" #include "boost/make_shared.hpp" #include "RTToolboxConfigure.h" #include "rttbException.h" rttb::apps::bioModelCalc::ApplicationData appData; int main(int argc, const char** argv) { int result = 0; boost::shared_ptr argParser; try { std::string appName = "BioModelCalc"; std::string appVersion = RTTB_FULL_VERSION_STRING; argParser = boost::make_shared(argc, argv, appName, appVersion); } catch (const std::exception& e) { std::cerr << e.what() << std::endl; return -1; } // This is vital. The application needs to exit if the "help" or "version" parameter is set // because this means the other parameters won't be parsed. if (argParser->isSet(argParser->OPTION_HELP) || argParser->isSet(argParser->OPTION_VERSION)) { return 0; } rttb::apps::bioModelCalc::populateAppData(argParser, appData); std::cout << std::endl << "*******************************************" << std::endl; std::cout << "Dose file: " << appData._doseFileName << std::endl; std::cout << "Bio model output file: " << appData._outputFileName << std::endl; std::cout << "Model: " << appData._model << std::endl; std::cout << "Model parameters: "; - for (unsigned int i = 0; i < appData._modelParameters.size(); i++) + for (size_t i = 0; i < appData._modelParameters.size(); i++) { if (i != 0) { std::cout << ", "; } std::cout << appData._modelParameters.at(i); } std::cout << std::endl; std::cout << "Dose scaling: " << appData._doseScaling << std::endl; std::cout << std::endl; try { appData._dose = rttb::apps::bioModelCalc::loadDose(appData._doseFileName, appData._doseLoadStyle); } catch (rttb::core::Exception& e) { std::cerr << "RTTB Error!!!" << std::endl; std::cerr << e.what() << std::endl; return 1; } catch (const std::exception& e) { std::cerr << "Error!!!" << std::endl; std::cerr << e.what() << std::endl; return 1; } catch (...) { std::cerr << "Error!!! unknown error while reading input image." << std::endl; return 1; } try { rttb::apps::bioModelCalc::processData(appData); } catch (rttb::core::Exception& e) { std::cerr << "RTTB Error!!!" << std::endl; std::cerr << e.what() << std::endl; return 2; } catch (std::exception& e) { std::cerr << "Error!!!" << std::endl; std::cerr << e.what() << std::endl; return 2; } catch (...) { std::cerr << "Error!!! unknown error while calculating the bioModel or writing the image." << std::endl; return 2; } return result; } diff --git a/demoapps/DoseTool/DoseToolHelper.cpp b/demoapps/DoseTool/DoseToolHelper.cpp index 680d8a4..8bda333 100644 --- a/demoapps/DoseTool/DoseToolHelper.cpp +++ b/demoapps/DoseTool/DoseToolHelper.cpp @@ -1,369 +1,369 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "DoseToolHelper.h" #include "boost/make_shared.hpp" #include "boost/shared_ptr.hpp" #include "boost/property_tree/ptree.hpp" #include "boost/filesystem.hpp" #include "rttbExceptionMacros.h" #include "rttbVirtuosPlanFileDoseAccessorGenerator.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomHelaxFileDoseAccessorGenerator.h" #include "rttbITKImageFileAccessorGenerator.h" #include "rttbStructureSetGeneratorInterface.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbVirtuosFileStructureSetGenerator.h" #include "rttbITKImageFileMaskAccessorGenerator.h" #include "rttbDVHCalculator.h" #include "rttbDVHXMLFileWriter.h" #include "rttbDoseStatisticsCalculator.h" #include "rttbBoostMaskAccessor.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbDoseStatisticsXMLWriter.h" #include "rttbVOIindexIdentifier.h" rttb::core::DoseAccessorInterface::DoseAccessorPointer rttb::apps::doseTool::loadDose(const std::string& fileName, const rttb::apps::doseTool::ApplicationData::LoadingStyleArgType& args) { rttb::core::DoseAccessorInterface::DoseAccessorPointer result; std::cout << std::endl << "read dose file... "; if (args.empty() || args[0] == "dicom") { std::cout << "use RTTB dicom IO... "; result = loadDicomDose(fileName); } else if (args[0] == "helax") { std::cout << "use RTTB Helax IO... "; result = loadHelaxDose(fileName); } else if (args[0] == "itk") { std::cout << "use RTTB itk IO... "; result = loadITKDose(fileName); } else if (args[0] == "virtuos") { if (args.size() < 2) { rttbDefaultExceptionStaticMacro( << "Cannot load virtuos dose. Plan file is missing. Specify plan file as 2nd io stlye argument."); } std::cout << "use RTTB virtuos IO... " << std::endl; std::cout << " virtuos plan file: " << args[1] << " ... "; result = loadVirtuosDose(fileName, args[1]); } else { rttbDefaultExceptionStaticMacro( << "Unknown io style selected. Cannot load data. Selected style: " << args[0]); } std::cout << "done." << std::endl; return result; }; rttb::core::DoseAccessorInterface::DoseAccessorPointer rttb::apps::doseTool::loadDicomDose(const std::string& fileName) { rttb::io::dicom::DicomFileDoseAccessorGenerator generator(fileName); return generator.generateDoseAccessor(); }; rttb::core::DoseAccessorInterface::DoseAccessorPointer rttb::apps::doseTool::loadHelaxDose(const std::string& path) { rttb::io::helax::DicomHelaxFileDoseAccessorGenerator generator(path); return generator.generateDoseAccessor(); }; rttb::core::DoseAccessorInterface::DoseAccessorPointer rttb::apps::doseTool::loadITKDose(const std::string& fileName) { rttb::io::itk::ITKImageFileAccessorGenerator generator(fileName); return generator.generateDoseAccessor(); }; rttb::core::DoseAccessorInterface::DoseAccessorPointer rttb::apps::doseTool::loadVirtuosDose(const std::string& fileName, const std::string& planFileName) { rttb::io::virtuos::VirtuosPlanFileDoseAccessorGenerator generator(fileName, planFileName); return generator.generateDoseAccessor(); }; rttb::core::StructureSetGeneratorInterface::StructureSetPointer rttb::apps::doseTool::loadStruct( const std::string& fileName, const ApplicationData::LoadingStyleArgType& args, const std::string& structNameRegex) { rttb::core::StructureSetGeneratorInterface::StructureSetPointer result; std::cout << std::endl << "read struct file... "; if (args.empty() || args[0] == "dicom") { std::cout << "use RTTB dicom IO... "; result = loadDicomStruct(fileName, structNameRegex); } else if (args[0] == "virtuos") { std::cout << "use RTTB virtuos IO... " << std::endl; std::cout << " virtuos CTX file: " << args[1] << " ... "; result = loadVirtuosStruct(fileName, args[1]); } else { rttbDefaultExceptionStaticMacro( << "Unknown io style selected. Cannot load data. Selected style: " << args[0]); } std::cout << "done." << std::endl; return result; } rttb::core::StructureSetGeneratorInterface::StructureSetPointer rttb::apps::doseTool::loadDicomStruct( const std::string& fileName, const std::string& structNameRegex) { rttb::io::dicom::DicomFileStructureSetGenerator generator(fileName); if (!structNameRegex.empty()) { generator.setStructureLableFilterActive(true); generator.setFilterRegEx(structNameRegex); } return generator.generateStructureSet(); } rttb::core::StructureSetGeneratorInterface::StructureSetPointer rttb::apps::doseTool::loadVirtuosStruct( const std::string& fileName, const std::string& ctxFileName) { rttb::io::virtuos::VirtuosFileStructureSetGenerator generator(fileName, ctxFileName.c_str()); return generator.generateStructureSet(); } void rttb::apps::doseTool::processData(rttb::apps::doseTool::ApplicationData& appData) { std::cout << std::endl << "generating masks... "; std::vector maskAccessorPtrVector = generateMasks( appData); std::cout << "done." << std::endl; - for (int i = 0; i < maskAccessorPtrVector.size(); i++) + for (size_t i = 0; i < maskAccessorPtrVector.size(); i++) { core::DoseIteratorInterface::DoseIteratorPointer spDoseIterator(generateMaskedDoseIterator( maskAccessorPtrVector.at(i), appData._dose)); if (appData._computeDoseStatistics) { std::cout << std::endl << "computing dose statistics... "; algorithms::DoseStatistics::DoseStatisticsPointer statistics = calculateDoseStatistics( spDoseIterator, appData._computeComplexDoseStatistics, appData._prescribedDose); std::cout << "done." << std::endl; std::cout << std::endl << "writing dose statistics to file... "; std::string outputFilename; if (appData._multipleStructsMode) { outputFilename = assembleFilenameWithStruct(appData._doseStatisticOutputFileName, appData._structNames.at(i)); } else { outputFilename = appData._doseStatisticOutputFileName; } writeDoseStatisticsFile(statistics, outputFilename, appData._structNames.at(i), appData); std::cout << "done." << std::endl; } if (appData._computeDVH) { std::cout << std::endl << "computing DVH... "; core::DVH::DVHPointer dvh = calculateDVH(spDoseIterator, appData._struct->getUID(), appData._dose->getUID()); std::cout << "done." << std::endl; std::cout << std::endl << "writing DVH to file... "; std::string outputFilename; if (appData._multipleStructsMode) { outputFilename = assembleFilenameWithStruct(appData._dvhOutputFilename, appData._structNames.at(i)); } else { outputFilename = appData._dvhOutputFilename; } writeDVHFile(dvh, outputFilename); std::cout << "done." << std::endl; } } } std::vector rttb::apps::doseTool::generateMasks( rttb::apps::doseTool::ApplicationData& appData) { std::vector maskAccessorPtrVector; if (appData._structLoadStyle.front() == "itk") { maskAccessorPtrVector.push_back(rttb::io::itk::ITKImageFileMaskAccessorGenerator( appData._structFileName).generateMaskAccessor()); } else { std::vector foundIndices = rttb::masks::VOIindexIdentifier::getIndicesByVoiRegex( appData._struct, appData._structNameRegex); std::vector relevantIndices; if (appData._multipleStructsMode) { relevantIndices = foundIndices; } else { if (!foundIndices.empty()) { relevantIndices.push_back(foundIndices.front()); } } appData._structIndices = relevantIndices; bool strict = !appData._allowSelfIntersection; - for (int i = 0; i < appData._structIndices.size(); i++) + for (size_t i = 0; i < appData._structIndices.size(); i++) { maskAccessorPtrVector.push_back( boost::make_shared (appData._struct->getStructure(appData._structIndices.at(i)), appData._dose->getGeometricInfo(), strict)); maskAccessorPtrVector.at(i)->updateMask(); appData._structNames.push_back(appData._struct->getStructure(appData._structIndices.at( i))->getLabel()); } } return maskAccessorPtrVector; } rttb::core::DoseIteratorInterface::DoseIteratorPointer rttb::apps::doseTool::generateMaskedDoseIterator( rttb::core::MaskAccessorInterface::MaskAccessorPointer maskAccessorPtr, rttb::core::DoseAccessorInterface::DoseAccessorPointer doseAccessorPtr) { boost::shared_ptr maskedDoseIterator = boost::make_shared(maskAccessorPtr, doseAccessorPtr); rttb::core::DoseIteratorInterface::DoseIteratorPointer doseIterator(maskedDoseIterator); return doseIterator; } rttb::algorithms::DoseStatistics::DoseStatisticsPointer rttb::apps::doseTool::calculateDoseStatistics( core::DoseIteratorInterface::DoseIteratorPointer doseIterator, bool calculateComplexDoseStatistics, DoseTypeGy prescribedDose) { rttb::algorithms::DoseStatisticsCalculator doseStatsCalculator(doseIterator); if (calculateComplexDoseStatistics) { return doseStatsCalculator.calculateDoseStatistics(prescribedDose); } else { return doseStatsCalculator.calculateDoseStatistics(); } } rttb::core::DVH::DVHPointer rttb::apps::doseTool::calculateDVH( core::DoseIteratorInterface::DoseIteratorPointer doseIterator, IDType structUID, IDType doseUID) { rttb::core::DVHCalculator calc(doseIterator, structUID, doseUID); rttb::core::DVH::DVHPointer dvh = calc.generateDVH(); return dvh; } void rttb::apps::doseTool::writeDoseStatisticsFile( rttb::algorithms::DoseStatistics::DoseStatisticsPointer statistics, const std::string& filename, const std::string& structName, rttb::apps::doseTool::ApplicationData& appData) { boost::property_tree::ptree originalTree = rttb::io::other::writeDoseStatistics(statistics); //add config part to xml originalTree.add("statistics.config.requestedStructRegex", appData._structNameRegex); originalTree.add("statistics.config.structName", structName); originalTree.add("statistics.config.doseUID", appData._dose->getUID()); originalTree.add("statistics.config.doseFile", appData._doseFileName); originalTree.add("statistics.config.structFile", appData._structFileName); boost::property_tree::ptree reorderedTree, configTree, resultsTree; configTree = originalTree.get_child("statistics.config"); resultsTree = originalTree.get_child("statistics.results"); reorderedTree.add_child("statistics.config", configTree); reorderedTree.add_child("statistics.results", resultsTree); boost::property_tree::write_xml(filename, reorderedTree, std::locale(), boost::property_tree::xml_writer_make_settings('\t', 1)); } std::string rttb::apps::doseTool::assembleFilenameWithStruct(const std::string& originalFilename, const std::string& structName) { boost::filesystem::path originalFile(originalFilename); std::string newFilename = originalFile.stem().string() + "_" + structName + originalFile.extension().string(); boost::filesystem::path newFile(originalFile.parent_path() / newFilename); return newFile.string(); } void rttb::apps::doseTool::writeDVHFile(core::DVH::DVHPointer dvh, const std::string& filename) { DVHType typeCum = { DVHType::Cumulative }; rttb::io::other::DVHXMLFileWriter dvhWriter(filename, typeCum); dvhWriter.writeDVH(dvh); } diff --git a/demoapps/VoxelizerTool/rttbMaskWriter.cpp b/demoapps/VoxelizerTool/rttbMaskWriter.cpp index 951f411..3be65a4 100644 --- a/demoapps/VoxelizerTool/rttbMaskWriter.cpp +++ b/demoapps/VoxelizerTool/rttbMaskWriter.cpp @@ -1,133 +1,133 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "rttbMaskWriter.h" #include "itkImage.h" #include "itkImageFileWriter.h" #include "itkBinaryThresholdImageFilter.h" #include "itkAddImageFilter.h" namespace rttb { namespace apps { namespace voxelizer { MaskWriter::MaskWriter(std::vector maskPointer, bool voxelization) : _maskPointerVector(maskPointer), _booleanvoxelization(voxelization) { } void MaskWriter::writeMaskToFile(const std::string& outputFileName) const { if (!_maskPointerVector.empty()) { ITKImageTypeConstPointer itkImage; if (_maskPointerVector.size() > 1) { itkImage = addMultipleStructsToImage(); } else { io::itk::ITKImageMaskAccessorConverter maskAccessorConverter(_maskPointerVector.at(0)); maskAccessorConverter.process(); itkImage = maskAccessorConverter.getITKImage(); } if (_booleanvoxelization) { itkImage = applyThresholdFilter(itkImage); } writeITKImageToFile(itkImage, outputFileName); } } MaskWriter::ITKImageTypeConstPointer MaskWriter::addMultipleStructsToImage() const { std::vector listOfITKImages; - for (int i = 0; i < _maskPointerVector.size(); i++) + for (size_t i = 0; i < _maskPointerVector.size(); i++) { io::itk::ITKImageMaskAccessorConverter maskAccessorConverter(_maskPointerVector.at(i)); maskAccessorConverter.process(); listOfITKImages.push_back(maskAccessorConverter.getITKImage()); } typedef itk::AddImageFilter AddImageFilterType; AddImageFilterType::Pointer addFilter = AddImageFilterType::New(); ITKImageTypePointer filterResult; - for (int k = 1; k < listOfITKImages.size(); k++) + for (size_t k = 1; k < listOfITKImages.size(); k++) { if (k == 1) { addFilter->SetInput1(listOfITKImages.at(0)); } else { addFilter->SetInput1(filterResult); } addFilter->SetInput2(listOfITKImages.at(k)); addFilter->Update(); filterResult = addFilter->GetOutput(); } ITKImageTypeConstPointer filterResultConst(filterResult); return filterResultConst; } MaskWriter::ITKImageTypeConstPointer MaskWriter::applyThresholdFilter( ITKImageTypeConstPointer itkImage) const { typedef itk::BinaryThresholdImageFilter< ITKMaskImageType, ITKMaskImageType > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput(itkImage); filter->SetLowerThreshold(0.5); filter->SetUpperThreshold(1.0); filter->SetInsideValue(1.0); filter->Update(); return filter->GetOutput(); } void MaskWriter::writeITKImageToFile(ITKImageTypeConstPointer itkImage, const std::string& outputfilename, bool useCompression) const { typedef itk::ImageFileWriter< ITKMaskImageType > WriterType; WriterType::Pointer writer = WriterType::New(); writer->SetFileName(outputfilename); writer->SetUseCompression(useCompression); writer->SetInput(itkImage); writer->Update(); } } } } diff --git a/demoapps/VoxelizerTool/rttbStructDataReader.cpp b/demoapps/VoxelizerTool/rttbStructDataReader.cpp index 0433868..d1c29e7 100644 --- a/demoapps/VoxelizerTool/rttbStructDataReader.cpp +++ b/demoapps/VoxelizerTool/rttbStructDataReader.cpp @@ -1,117 +1,117 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "rttbStructDataReader.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbITKImageFileAccessorGenerator.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace apps { namespace voxelizer { StructDataReader::StructDataReader(const std::string& structFileName, const std::string& referenceFileName, const std::vector& referenceFileLoadingStyle) : _referenceFilename(referenceFileName), _structFilename(structFileName), _referenceFileLoadingStyle(referenceFileLoadingStyle) { } void StructDataReader::read() { _doseAccessor = readReferenceFile(_referenceFilename, _referenceFileLoadingStyle); _rtStructureSet = readStructFile(_structFilename); } std::vector StructDataReader::getAllLabels() const { std::vector allLabels; if (_rtStructureSet != NULL) { - for (int j = 0; j < _rtStructureSet->getNumberOfStructures(); j++) + for (size_t j = 0; j < _rtStructureSet->getNumberOfStructures(); j++) { allLabels.push_back(_rtStructureSet->getStructure(j)->getLabel()); } } return allLabels; } StructDataReader::StructureSetPointer StructDataReader::getStructureSetPointer() const { return _rtStructureSet; } StructDataReader::DoseAccessorPointer StructDataReader::getDoseAccessorPointer() const { return _doseAccessor; } StructDataReader::DoseAccessorPointer StructDataReader::readReferenceFile( const std::string& filename, const std::vector& fileLoadingStyle) const { if (fileLoadingStyle.at(0) == "dicom") { return readDicomFile(filename); } else if (fileLoadingStyle.at(0) == "itk") { return readITKFile(filename); } else { return NULL; } } StructDataReader::DoseAccessorPointer StructDataReader::readDicomFile( const std::string& filename) const { rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator(filename.c_str()); return doseAccessorGenerator.generateDoseAccessor(); } StructDataReader::DoseAccessorPointer StructDataReader::readITKFile(const std::string& filename) const { rttb::io::itk::ITKImageFileAccessorGenerator generator(filename); return generator.generateDoseAccessor(); } StructDataReader::StructureSetPointer StructDataReader::readStructFile( const std::string& filename) const { StructureSetPointer rtStructureSet = rttb::io::dicom::DicomFileStructureSetGenerator( filename.c_str()).generateStructureSet(); return rtStructureSet; } } } -} \ No newline at end of file +} diff --git a/demoapps/VoxelizerTool/rttbVoxelizerTool.cpp b/demoapps/VoxelizerTool/rttbVoxelizerTool.cpp index d4efc32..b6fbc66 100644 --- a/demoapps/VoxelizerTool/rttbVoxelizerTool.cpp +++ b/demoapps/VoxelizerTool/rttbVoxelizerTool.cpp @@ -1,181 +1,181 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include "itkMacro.h" #include "rttbVoxelizerHelper.h" #include "rttbMaskProcess.h" #include "rttbMaskWriter.h" #include "rttbStructDataReader.h" #include "rttbCommandOptions.h" int main(int argc, char* argv[]) { typedef rttb::core::MaskAccessorInterface::MaskAccessorPointer MaskAccessorPointer; boost::shared_ptr co = boost::make_shared(); if (!co->command(argc, argv)) { return 1; } if (co->isReturnAfterHelp()) { return 0; } rttb::apps::voxelizer::Parameters params = co->getParameters(); boost::shared_ptr reader = boost::make_shared(params.structFile, params.referenceFile, params.referenceFileLoadStyle); std::cout << "reading reference and structure file..."; try { reader->read(); } catch (rttb::core::Exception& e) { std::cerr << "RTTB Error!!!" << std::endl; std::cerr << e.what() << std::endl; return 2; } catch (const std::exception& e) { std::cerr << "Error!!!" << std::endl; std::cerr << e.what() << std::endl; return 2; } catch (...) { std::cerr << "Error!!! unknown error while reading input image." << std::endl; return 2; } std::cout << "done." << std::endl; std::cout << "searching for structs..."; std::vector listOfCorrectElements; - for (int i = 0; i < params.regEx.size(); i++) + for (size_t i = 0; i < params.regEx.size(); i++) { std::vector indexOfCorrectElements; indexOfCorrectElements = rttb::apps::voxelizer::filterForExpression(reader->getAllLabels(), params.regEx.at(i)); std::copy(indexOfCorrectElements.begin(), indexOfCorrectElements.end(), std::back_inserter(listOfCorrectElements)); } std::cout << "done." << std::endl; boost::shared_ptr maskProcessor = boost::make_shared(reader->getStructureSetPointer(), reader->getDoseAccessorPointer(), params.legacyVoxelization, params.allowSelfIntersections); if (!listOfCorrectElements.empty()) { std::string fileName = rttb::apps::voxelizer::getFilenameWithoutEnding( params.outputFilename); std::string fileEnding = rttb::apps::voxelizer::getFileEnding(params.outputFilename); std::vector maskVector; if (params.addStructures) { std::string labelName; - for (int i = 0; i < listOfCorrectElements.size(); i++) + for (size_t i = 0; i < listOfCorrectElements.size(); i++) { int labelIndex = listOfCorrectElements.at(i); maskVector.push_back(maskProcessor->createMask(labelIndex)); std::vector labelVector = reader->getAllLabels(); std::string labelOfInterest = labelVector.at(labelIndex); rttb::apps::voxelizer::removeSpecialCharacters(labelOfInterest); labelName += "_" + labelOfInterest; } boost::shared_ptr writer = boost::make_shared(maskVector, params.booleanVoxelization); writer->writeMaskToFile(fileName + labelName + fileEnding); } else { unsigned int maxIterationCount = 1; if (params.multipleStructs) { maxIterationCount = listOfCorrectElements.size(); } for (unsigned int i = 0; i < maxIterationCount; i++) { std::cout << "creating mask..."; maskVector.push_back(maskProcessor->createMask(listOfCorrectElements.at(i))); std::cout << "done" << std::endl; int labelIndex = listOfCorrectElements.at(i); std::vector labelVector = reader->getAllLabels(); std::string labelOfInterest = labelVector.at(labelIndex); rttb::apps::voxelizer::removeSpecialCharacters(labelOfInterest); boost::shared_ptr MW = boost::make_shared(maskVector, params.booleanVoxelization); try { MW->writeMaskToFile(fileName + "_" + labelOfInterest + fileEnding); } catch (itk::ExceptionObject& err) { std::cerr << "ExceptionObject caught !" << std::endl; std::cerr << err << std::endl; return 3; } catch (const std::exception& e) { std::cerr << "Error!!!" << std::endl; std::cerr << e.what() << std::endl; return 3; } catch (...) { std::cerr << "Error!!! unknown error while reading input image." << std::endl; return 3; } } } } else { std::cout << "No struct found" << std::endl; } return 0; }