diff --git a/code/interpolation/rttbInterpolationBase.cpp b/code/interpolation/rttbInterpolationBase.cpp index 3024109..09def85 100644 --- a/code/interpolation/rttbInterpolationBase.cpp +++ b/code/interpolation/rttbInterpolationBase.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 "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); - /*_spOriginalData->getGeometricInfo().geometryCoordinateToWorldCoordinate(DoubleVoxelGridIndex3D( - aIndex[0], aIndex[1], aIndex[2]), theNextVoxel);*/ + //_spOriginalData->getGeometricInfo().indexToWorldCoordinate(aIndex, theNextVoxel); + _spOriginalData->getGeometricInfo().geometryCoordinateToWorldCoordinate(DoubleVoxelGridIndex3D( + aIndex[0], aIndex[1], aIndex[2]), 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++) { 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 (int yIncr = 0; yIncr < 2; yIncr++) { for (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++) { 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; //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; 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/rttbRosuMappableDoseAccessor.cpp b/code/interpolation/rttbRosuMappableDoseAccessor.cpp index 7e05512..2b7015e 100644 --- a/code/interpolation/rttbRosuMappableDoseAccessor.cpp +++ b/code/interpolation/rttbRosuMappableDoseAccessor.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 "rttbRosuMappableDoseAccessor.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" #include "rttbLinearInterpolation.h" namespace rttb { namespace interpolation { RosuMappableDoseAccessor::RosuMappableDoseAccessor(const core::GeometricInfo& geoInfoTargetImage, const DoseAccessorPointer doseMovingImage, const TransformationInterface::Pointer aTransformation, bool acceptPadding, DoseTypeGy defaultOutsideValue): MappableDoseAccessorInterface(geoInfoTargetImage, doseMovingImage, aTransformation, acceptPadding, defaultOutsideValue) { //define linear interpolation InterpolationBase::Pointer interpolationLinear = InterpolationBase::Pointer( new LinearInterpolation()); _spInterpolation = interpolationLinear; _spInterpolation->setAccessorPointer(_spOriginalDoseDataMovingImage); } GenericValueType RosuMappableDoseAccessor::getValueAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex3D; if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) { return getValueAt(aVoxelGridIndex3D); } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); return -1; } } } GenericValueType RosuMappableDoseAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { //Transform requested voxel coordinates of original image into world coordinates with RTTB WorldCoordinate3D worldCoordinateTarget; - if (_geoInfoTargetImage.indexToWorldCoordinate(aIndex, worldCoordinateTarget)) - /*if (_geoInfoTargetImage.geometryCoordinateToWorldCoordinate(DoubleVoxelGridIndex3D( - aIndex[0], aIndex[1], aIndex[2]), worldCoordinateTarget))*/ + //if (_geoInfoTargetImage.indexToWorldCoordinate(aIndex, worldCoordinateTarget)) + if (_geoInfoTargetImage.geometryCoordinateToWorldCoordinate(DoubleVoxelGridIndex3D( + aIndex[0], aIndex[1], aIndex[2]), worldCoordinateTarget)) { std::vector octants = getOctants(worldCoordinateTarget); if (octants.size() > 2) { DoseTypeGy interpolatedDoseValue = 0.0; //get trilinear interpolation value of every octant point for (std::vector::const_iterator octantIterator = octants.begin(); octantIterator != octants.end(); ++octantIterator) { //transform coordinates WorldCoordinate3D worldCoordinateMoving; WorldCoordinate3D worldCoordinateTargetOctant = *octantIterator; _spTransformation->transformInverse(worldCoordinateTargetOctant, worldCoordinateMoving); try { interpolatedDoseValue += _spInterpolation->getValue(worldCoordinateMoving); } //Mapped outside of image? Check if padding is allowed catch (core::MappingOutsideOfImageException& /*e*/) { if (_acceptPadding) { interpolatedDoseValue += _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Mapping outside of image"); } } catch (core::Exception& e) { std::cout << e.what() << std::endl; return -1; } } return interpolatedDoseValue / (DoseTypeGy)octants.size(); } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Too many samples are mapped outside the image!"); return -1; } } } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); return -1; } } } std::vector RosuMappableDoseAccessor::getOctants( const WorldCoordinate3D& aCoordinate) const { std::vector octants; SpacingVectorType3D spacingTargetImage = _geoInfoTargetImage.getSpacing(); core::GeometricInfo geometricInfoDoseData = _spOriginalDoseDataMovingImage->getGeometricInfo(); //as the corner point is the coordinate of the voxel (grid), 0.25 and 0.75 are the center of the subvoxels for (double xOct = -0.25; xOct <= 0.25; xOct += 0.5) { for (double yOct = -0.25; yOct <= 0.25; yOct += 0.5) { for (double zOct = -0.25; zOct <= 0.25; zOct += 0.5) { WorldCoordinate3D aWorldCoordinate(aCoordinate.x() + (xOct * spacingTargetImage.x()), aCoordinate.y() + (yOct * spacingTargetImage.y()), aCoordinate.z() + (zOct * spacingTargetImage.z())); if (geometricInfoDoseData.isInside(aWorldCoordinate)) { octants.push_back(aWorldCoordinate); } } } } return octants; } }//end namespace interpolation }//end namespace rttb diff --git a/code/interpolation/rttbSimpleMappableDoseAccessor.cpp b/code/interpolation/rttbSimpleMappableDoseAccessor.cpp index a0287de..e4b4856 100644 --- a/code/interpolation/rttbSimpleMappableDoseAccessor.cpp +++ b/code/interpolation/rttbSimpleMappableDoseAccessor.cpp @@ -1,124 +1,124 @@ // ----------------------------------------------------------------------- // 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 "rttbSimpleMappableDoseAccessor.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace interpolation { SimpleMappableDoseAccessor::SimpleMappableDoseAccessor(const core::GeometricInfo& geoInfoTargetImage, const DoseAccessorPointer doseMovingImage, TransformationInterface::Pointer aTransformation, const InterpolationBase::Pointer aInterpolation, bool acceptPadding, double defaultOutsideValue): MappableDoseAccessorInterface(geoInfoTargetImage, doseMovingImage, aTransformation, acceptPadding, defaultOutsideValue), _spInterpolation(aInterpolation) { //handle null pointers if (aInterpolation == NULL) { throw core::NullPointerException("Pointer to interpolation cannot be NULL."); } else { _spInterpolation->setAccessorPointer(_spOriginalDoseDataMovingImage); } } GenericValueType SimpleMappableDoseAccessor::getValueAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex3D; if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) { return getValueAt(aVoxelGridIndex3D); } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); return -1; } } } GenericValueType SimpleMappableDoseAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { //Transform requested voxel coordinates of original image into world coordinates with RTTB WorldCoordinate3D worldCoordinateTarget; - if (_geoInfoTargetImage.indexToWorldCoordinate(aIndex, worldCoordinateTarget)) - /*if (_geoInfoTargetImage.geometryCoordinateToWorldCoordinate(DoubleVoxelGridIndex3D( - aIndex[0], aIndex[1], aIndex[2]), worldCoordinateTarget))*/ + //if (_geoInfoTargetImage.indexToWorldCoordinate(aIndex, worldCoordinateTarget)) + if (_geoInfoTargetImage.geometryCoordinateToWorldCoordinate(DoubleVoxelGridIndex3D( + aIndex[0], aIndex[1], aIndex[2]), worldCoordinateTarget)) { //transform coordinates WorldCoordinate3D worldCoordinateMoving; _spTransformation->transformInverse(worldCoordinateTarget, worldCoordinateMoving); //Use Interpolation to compute dose at mappedImage try { return _spInterpolation->getValue(worldCoordinateMoving); } //Mapped outside of image? Check if padding is allowed catch (core::MappingOutsideOfImageException& /*e*/) { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); return -1; } } catch (core::Exception& e) { std::cout << e.what() << std::endl; return -1; } } //ok, if that fails, throw exception. Makes no sense to go further else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); return -1; } } } }//end namespace interpolation }//end namespace rttb