diff --git a/code/interpolation/rttbInterpolationBase.cpp b/code/interpolation/rttbInterpolationBase.cpp index 282bb3f..0e5cc56 100644 --- a/code/interpolation/rttbInterpolationBase.cpp +++ b/code/interpolation/rttbInterpolationBase.cpp @@ -1,191 +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. // //------------------------------------------------------------------------ #include #include #include "rttbInterpolationBase.h" #include "rttbInvalidParameterException.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace interpolation { void InterpolationBase::setAccessorPointer(core::AccessorInterface::ConstPointer originalData) { if (originalData != nullptr) { _spOriginalData = originalData; } else { throw core::NullPointerException("originalDose is nullptr!"); } }; + core::AccessorInterface::ConstPointer InterpolationBase::getAccessorPointer() const + { + return _spOriginalData; + } + void InterpolationBase::getNeighborhoodVoxelValues( const WorldCoordinate3D& aWorldCoordinate, unsigned int neighborhood, std::array& target, boost::shared_ptr values) const { if (_spOriginalData == nullptr) { throw core::NullPointerException("originalDose is nullptr!"); } //Determine target (abs(desired worldCoordinate- corner pixel world coordinate/pixel spacing) and values of corner pixels (from originalDose) VoxelGridIndex3D aIndex; if (_spOriginalData->getGeometricInfo().worldCoordinateToIndex(aWorldCoordinate, aIndex)) { //determine just the nearest voxel to the world coordinate if (neighborhood == 0) { values[0] = _spOriginalData->getValueAt(aIndex); } //determine the 8 voxels around the world coordinate else if (neighborhood == 8) { std::list cornerPoints; WorldCoordinate3D theNextVoxel; _spOriginalData->getGeometricInfo().indexToWorldCoordinate(aIndex, theNextVoxel); SpacingVectorType3D pixelSpacing = (_spOriginalData->getGeometricInfo()).getSpacing(); VoxelGridIndex3D leftTopFrontCoordinate; //find the voxel with the smallest coordinate values in each dimension. This defines the standard cube for (unsigned int i = 0; i < 3; i++) { if (aWorldCoordinate[i] < theNextVoxel[i]) { if (aIndex[i] > 0) { leftTopFrontCoordinate[i] = aIndex[i] - 1; target[i] = (aWorldCoordinate[i] - (theNextVoxel[i] - pixelSpacing[i])) / pixelSpacing[i]; } //@todo: see T22315 else { leftTopFrontCoordinate[i] = aIndex[i]; target[i] = (aWorldCoordinate[i] - (theNextVoxel[i] - pixelSpacing[i])) / pixelSpacing[i]; } } else { leftTopFrontCoordinate[i] = aIndex[i]; target[i] = (aWorldCoordinate[i] - theNextVoxel[i]) / pixelSpacing[i]; } } for (unsigned int zIncr = 0; zIncr < 2; zIncr++) { for (unsigned int yIncr = 0; yIncr < 2; yIncr++) { for (unsigned int xIncr = 0; xIncr < 2; xIncr++) { cornerPoints.emplace_back(leftTopFrontCoordinate[0] + xIncr, leftTopFrontCoordinate[1] + yIncr, leftTopFrontCoordinate[2] + zIncr); } } } //target range has to be always [0,1] for (unsigned int i = 0; i < 3; i++) { assert(target[i] >= 0.0 && target[i] <= 1.0); } unsigned int count = 0; //now just get the values of all (dose) voxels and store them in values for (auto cornerPointsIterator = cornerPoints.begin(); cornerPointsIterator != cornerPoints.end(); ++cornerPointsIterator, ++count) { if (_spOriginalData->getGeometricInfo().isInside(*cornerPointsIterator)) { values[count] = _spOriginalData->getValueAt(*cornerPointsIterator); } else { //outside value! boundary treatment values[count] = getNearestInsideVoxelValue(*cornerPointsIterator); } assert(values[count] != -1); } } else { throw core::InvalidParameterException("neighborhoods other than 0 and 8 not yet supported in Interpolation"); } } else { throw core::MappingOutsideOfImageException("Error in conversion from world coordinates to index"); } } DoseTypeGy InterpolationBase::getNearestInsideVoxelValue(const VoxelGridIndex3D& currentVoxelIndex) const { VoxelGridIndex3D voxelChangedXYZ[] = {currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex}; unsigned int runningIndex; //x,y,z for (runningIndex = 0; runningIndex < 3; ++runningIndex) { voxelChangedXYZ[runningIndex][runningIndex] -= 1; } //xy voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][1] -= 1; ++runningIndex; //xz voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; //yz voxelChangedXYZ[runningIndex][1] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; //xyz voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][1] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; unsigned int replacementVoxelIndex = 0; while (replacementVoxelIndex < runningIndex) { if (_spOriginalData->getGeometricInfo().validIndex(voxelChangedXYZ[replacementVoxelIndex])) { return _spOriginalData->getValueAt(voxelChangedXYZ[replacementVoxelIndex]); } ++replacementVoxelIndex; } return -1; } }//end namespace core }//end namespace rttb diff --git a/code/interpolation/rttbInterpolationBase.h b/code/interpolation/rttbInterpolationBase.h index 990f775..f77c8ed 100644 --- a/code/interpolation/rttbInterpolationBase.h +++ b/code/interpolation/rttbInterpolationBase.h @@ -1,94 +1,96 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #ifndef __INTERPOLATION_BASE_H #define __INTERPOLATION_BASE_H #include #include #include #include "rttbAccessorInterface.h" #include "RTTBInterpolationExports.h" #ifdef _MSC_VER #pragma warning(push) #pragma warning(disable: 4251) #endif namespace rttb { namespace interpolation { /*! @class InterpolationBase @brief Base class for interpolation. @ingroup interpolation */ class RTTBInterpolation_EXPORT InterpolationBase { public: rttbClassMacroNoParent(InterpolationBase) /*! @brief Constructor */ InterpolationBase() = default; /*! @brief Virtual destructor of base class */ virtual ~InterpolationBase() = default; /*! @brief Sets the AccessorPointer @pre originalData initialized @exception core::NullPointerException if originalData==nullptr */ void setAccessorPointer(core::AccessorInterface::ConstPointer originalData); + core::AccessorInterface::ConstPointer getAccessorPointer() const; + /*! @brief Returns the interpolated value for the given world coordinate */ virtual DoseTypeGy getValue(const WorldCoordinate3D& aWorldCoordinate) const = 0; protected: rttb::core::AccessorInterface::ConstPointer _spOriginalData; /*! @brief determines voxels in a certain neighborhood of a physical based coordinate and converts in a standard cube with corner points [0 0 0], [1 0 0], [0 1 0], [1 1 0], [0 0 1], [1 0 1], [0 1 1], [1 1 1]. @param aWorldCoordinate the coordinate where to start @param neighborhood voxel around coordinate (currently only 0 and 8 implemented) @param target coordinates inside the standard cube with values [0 1] in each dimension. @param values dose values at all corner points of the standard cube. Is of type boost:shared_ptr[neighborhood] @pre target and values have to be correctly initialized (e.g. std::array target = {0.0, 0.0, 0.0}; boost::shared_ptr values(new DoseTypeGy[8]()); where 8 is neighborhood) @exception core::InvalidParameterException if neighborhood =! 0 && !=8 @exception core::MappingOutsideOfImageException if initial mapping of aWorldCoordinate is outside image @exception core::NullPointerException if dose is nullptr */ void getNeighborhoodVoxelValues(const WorldCoordinate3D& aWorldCoordinate, unsigned int neighborhood, std::array& target, boost::shared_ptr values) const; /*! @brief returns the nearest inside voxel value @pre the voxelGridIndex is outside the image and voxelGridIndex>image.size() for all dimensions. Also voxelGridIndex[]>=0 for all dimensions @note used for virtually expanding the image by one voxel as edge handling */ DoseTypeGy getNearestInsideVoxelValue(const VoxelGridIndex3D& currentVoxelIndex) const; }; } } #ifdef _MSC_VER #pragma warning(pop) #endif #endif