diff --git a/code/interpolation/rttbRosuMappableDoseAccessor.cpp b/code/interpolation/rttbRosuMappableDoseAccessor.cpp index 466c403..bcb9a2c 100644 --- a/code/interpolation/rttbRosuMappableDoseAccessor.cpp +++ b/code/interpolation/rttbRosuMappableDoseAccessor.cpp @@ -1,176 +1,167 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) -*/ #include "rttbRosuMappableDoseAccessor.h" #include #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" #include "rttbLinearInterpolation.h" namespace rttb { namespace interpolation { RosuMappableDoseAccessor::RosuMappableDoseAccessor(const core::GeometricInfo& geoInfoTargetImage, const DoseAccessorPointer doseMovingImage, const TransformationPointer aTransformation, bool acceptPadding, DoseTypeGy defaultOutsideValue): MappableDoseAccessorInterface(geoInfoTargetImage, doseMovingImage, aTransformation, acceptPadding, defaultOutsideValue) { //define linear interpolation auto interpolationLinear = ::boost::make_shared(); _spInterpolation = interpolationLinear; _spInterpolation->setAccessorPointer(_spOriginalDoseDataMovingImage); } GenericValueType RosuMappableDoseAccessor::getValueAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex3D; if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) { return getValueAt(aVoxelGridIndex3D); } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); - 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)) { 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 15403a9..a2d47c5 100644 --- a/code/interpolation/rttbSimpleMappableDoseAccessor.cpp +++ b/code/interpolation/rttbSimpleMappableDoseAccessor.cpp @@ -1,122 +1,113 @@ // ----------------------------------------------------------------------- // 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 == nullptr) { throw core::NullPointerException("Pointer to interpolation cannot be nullptr."); } else { _spInterpolation->setAccessorPointer(_spOriginalDoseDataMovingImage); } } GenericValueType SimpleMappableDoseAccessor::getValueAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex3D; if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) { return getValueAt(aVoxelGridIndex3D); } else { if (_acceptPadding) { return _defaultOutsideValue; } else { throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); - 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)) { //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 diff --git a/testing/interpolation/InterpolationTest.cpp b/testing/interpolation/InterpolationTest.cpp index 75ba483..6be3f13 100644 --- a/testing/interpolation/InterpolationTest.cpp +++ b/testing/interpolation/InterpolationTest.cpp @@ -1,214 +1,221 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbNearestNeighborInterpolation.h" #include "rttbLinearInterpolation.h" #include "rttbDicomDoseAccessor.h" #include "rttbGenericDoseIterator.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace testing { typedef rttb::interpolation::NearestNeighborInterpolation NearestNeighborInterpolation; typedef rttb::interpolation::LinearInterpolation LinearInterpolation; typedef rttb::core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; /*! @brief InterpolationTest - tests only interpolation 1) test both interpolation types with simple image (Dose = 2) 2) test both interpolation types with increasing x image values image (Dose = y value) - 3) test exception handling + 3) test right corner interpolation + 4) test exception handling */ int InterpolationTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string RTDOSE_FILENAME_CONSTANT_TWO; std::string RTDOSE_FILENAME_INCREASE_X; if (argc > 2) { RTDOSE_FILENAME_CONSTANT_TWO = argv[1]; RTDOSE_FILENAME_INCREASE_X = argv[2]; } else { std::cout << "at least two parameters for InterpolationTest expected" << std::endl; return -1; } rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1( RTDOSE_FILENAME_CONSTANT_TWO.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); DoseAccessorPointer doseAccessorNull; rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator2( RTDOSE_FILENAME_INCREASE_X.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); //doseAccessor1 is used as dose image - auto interpolationNN = boost::make_shared(); + auto interpolationNN = boost::make_shared(); interpolationNN->setAccessorPointer(doseAccessor1); auto interpolationLinear = boost::make_shared(); interpolationLinear->setAccessorPointer(doseAccessor1); //doseAccessor2 is used as dose image. //RTDOSE_FILENAME_INCREASE_X and RTDOSE_FILENAME_CONSTANT_TWO have the same GeometricInfo auto interpolationNN2 = boost::make_shared(); interpolationNN2->setAccessorPointer(doseAccessor2); auto interpolationLinear2 = boost::make_shared(); interpolationLinear2->setAccessorPointer(doseAccessor2); auto interpolationNullNN = boost::make_shared(); auto interpolationNullLinear = boost::make_shared(); rttb::WorldCoordinate3D imagePositionPatient = doseAccessor1->getGeometricInfo().getImagePositionPatient(); rttb::SpacingVectorType3D pixelSpacing = doseAccessor1->getGeometricInfo().getSpacing(); unsigned int size[] = {doseAccessor1->getGeometricInfo().getNumColumns(), doseAccessor1->getGeometricInfo().getNumRows(), doseAccessor1->getGeometricInfo().getNumSlices()}; //Which voxels to check is irrelevant. The following three situations are checked: // - exactly in-between two voxels: v_i=(0.5*v1 + 0.5*v2) --> target 0.5 // - the middle of a voxel: v_i = 1*v1 --> target 0 (semantically equivalent with target 1 and left/right flipped) // - 10% shifted from the middle: v_i=0.1*v1 + 0.9*v2 --> target 0.9 std::vector coordinatesX; coordinatesX.push_back(imagePositionPatient.x() - (pixelSpacing.x() * 0.5)); coordinatesX.push_back(imagePositionPatient.x() + 5 * pixelSpacing.x()); coordinatesX.push_back(imagePositionPatient.x() + 9.9 * pixelSpacing.x()); std::vector coordinatesY; coordinatesY.push_back(imagePositionPatient.y() - (pixelSpacing.y() * 0.5)); coordinatesY.push_back(imagePositionPatient.y() + 5 * pixelSpacing.y()); coordinatesY.push_back(imagePositionPatient.x() + 9.9 * pixelSpacing.y()); std::vector coordinatesZ; coordinatesZ.push_back(imagePositionPatient.z() - (pixelSpacing.z() * 0.5)); coordinatesZ.push_back(imagePositionPatient.z() + 5 * pixelSpacing.z()); coordinatesZ.push_back(imagePositionPatient.z() + 9.9 * pixelSpacing.z()); std::vector coordinatesToCheck; for (unsigned int x = 0; x < coordinatesX.size(); x++) { for (unsigned int y = 0; y < coordinatesY.size(); y++) { for (unsigned int z = 0; z < coordinatesZ.size(); z++) { coordinatesToCheck.push_back(rttb::WorldCoordinate3D(coordinatesX.at(x), coordinatesY.at(y), coordinatesZ.at(z))); } } } rttb::WorldCoordinate3D positionOutsideOfImageLeft = imagePositionPatient - rttb::WorldCoordinate3D( pixelSpacing.x() * 2.0, pixelSpacing.y() * 2.0, pixelSpacing.z() * 2.0); rttb::WorldCoordinate3D positionOutsideOfImageRight = imagePositionPatient + rttb::WorldCoordinate3D(size[0] * pixelSpacing.x(), size[1] * pixelSpacing.y(), size[2] * pixelSpacing.z()); + rttb::WorldCoordinate3D positionLastInsightImageRight = rttb::WorldCoordinate3D( + positionOutsideOfImageRight.x() - 0.5 * pixelSpacing.x() - 0.000001, + positionOutsideOfImageRight.y() - 0.5 * pixelSpacing.y() - 0.000001, + positionOutsideOfImageRight.z() - 0.5 * pixelSpacing.z() - 0.000001 + ); + + //precomputed values for Nearest neighbor + Linear interpolator double expectedDoseIncreaseXNearest[27]; double expectedDoseIncreaseXLinear[27]; for (int i = 0; i < 27; i++) { if (i < 9) { expectedDoseIncreaseXNearest[i] = 0.0; expectedDoseIncreaseXLinear[i] = 1.41119e-005; } else if (i < 18) { expectedDoseIncreaseXNearest[i] = 0.000141119; expectedDoseIncreaseXLinear[i] = 0.000141119; } else { expectedDoseIncreaseXNearest[i] = 0.000282239; expectedDoseIncreaseXLinear[i] = 0.000279416; } } //TEST 1) RTDOSE_FILENAME_CONSTANT_TWO contains Dose 2.0 Gy for each pixel, so for every interpolation, 2.0 has to be the result std::vector::const_iterator iterCoordinates = coordinatesToCheck.cbegin(); while (iterCoordinates != coordinatesToCheck.cend()) { CHECK_EQUAL(interpolationNN->getValue(*iterCoordinates), 2.0); CHECK_EQUAL(interpolationLinear->getValue(*iterCoordinates), 2.0); ++iterCoordinates; } //TEST 2) RTDOSE_FILENAME_INCREASE_X contains a Dose gradient file, correct interpolation values have been computed manually. To avoid floating point inaccuracy, CHECK_CLOSE is used iterCoordinates = coordinatesToCheck.cbegin(); unsigned int index = 0; while (iterCoordinates != coordinatesToCheck.cend() && index < 27) { CHECK_CLOSE(interpolationNN2->getValue(*iterCoordinates), expectedDoseIncreaseXNearest[index], errorConstant); CHECK_CLOSE(interpolationLinear2->getValue(*iterCoordinates), expectedDoseIncreaseXLinear[index], errorConstant); ++iterCoordinates; ++index; } - //TEST 3) Exception handling + //TEST 3) Right corner interpolation + //Checks if the interpolation works at the corner without an error + CHECK_NO_THROW(interpolationLinear->getValue(positionLastInsightImageRight)); + + + //TEST 4) Exception handling //Check that core::MappingOutOfImageException is thrown if requested position is outside image CHECK_THROW_EXPLICIT(interpolationNN->getValue(positionOutsideOfImageLeft), core::MappingOutsideOfImageException); CHECK_THROW_EXPLICIT(interpolationNN->getValue(positionOutsideOfImageRight), core::MappingOutsideOfImageException); CHECK_THROW_EXPLICIT(interpolationLinear->getValue(positionOutsideOfImageLeft), core::MappingOutsideOfImageException); CHECK_THROW_EXPLICIT(interpolationLinear->getValue(positionOutsideOfImageRight), core::MappingOutsideOfImageException); //Check that core::NullPointerException is thrown if Null Pointers are given to interpolator CHECK_THROW_EXPLICIT(interpolationNullLinear->setAccessorPointer(doseAccessorNull), core::NullPointerException); CHECK_THROW_EXPLICIT(interpolationNullNN->setAccessorPointer(doseAccessorNull), core::NullPointerException); CHECK_THROW_EXPLICIT(interpolationNullLinear->getValue(coordinatesToCheck.front()), core::NullPointerException); CHECK_THROW_EXPLICIT(interpolationNullNN->getValue(coordinatesToCheck.front()), core::NullPointerException); //Check that no exception is thrown otherwise CHECK_NO_THROW(boost::make_shared()); CHECK_NO_THROW(boost::make_shared()); RETURN_AND_REPORT_TEST_SUCCESS; } } } diff --git a/testing/interpolation/rttbInterpolationTests.cpp b/testing/interpolation/rttbInterpolationTests.cpp index c100175..a4180bc 100644 --- a/testing/interpolation/rttbInterpolationTests.cpp +++ b/testing/interpolation/rttbInterpolationTests.cpp @@ -1,65 +1,55 @@ // ----------------------------------------------------------------------- // 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) -*/ // this file defines the rttbAlgorithmsTests for the test driver // and all it expects is that you have a function called RegisterTests #if defined(_MSC_VER) #pragma warning ( disable : 4786 ) #endif #include "litMultiTestsMain.h" namespace rttb { namespace testing { void registerTests() { LIT_REGISTER_TEST(SimpleMappableDoseAccessorTest); LIT_REGISTER_TEST(RosuMappableDoseAccessorTest); LIT_REGISTER_TEST(InterpolationTest); } } } int main(int argc, char* argv[]) { int result = 0; rttb::testing::registerTests(); try { result = lit::multiTestsMain(argc, argv); } - catch (const std::exception& /*e*/) - { - result = -1; - } catch (...) { result = -1; } return result; }