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 <boost/make_shared.hpp>
#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<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))
{
std::vector<WorldCoordinate3D> octants = getOctants(worldCoordinateTarget);
if (octants.size() > 2)
{
DoseTypeGy interpolatedDoseValue = 0.0;
//get trilinear interpolation value of every octant point
for (std::vector<WorldCoordinate3D>::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<WorldCoordinate3D> RosuMappableDoseAccessor::getOctants(
const WorldCoordinate3D& aCoordinate) const
{
std::vector<WorldCoordinate3D> 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 <vector>
#include <boost/shared_ptr.hpp>
#include <boost/make_shared.hpp>
#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<rttb::interpolation::NearestNeighborInterpolation>();
+ auto interpolationNN = boost::make_shared<rttb::interpolation::NearestNeighborInterpolation>();
interpolationNN->setAccessorPointer(doseAccessor1);
auto interpolationLinear = boost::make_shared<rttb::interpolation::LinearInterpolation>();
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<rttb::interpolation::NearestNeighborInterpolation>();
interpolationNN2->setAccessorPointer(doseAccessor2);
auto interpolationLinear2 = boost::make_shared<rttb::interpolation::LinearInterpolation>();
interpolationLinear2->setAccessorPointer(doseAccessor2);
auto interpolationNullNN = boost::make_shared<rttb::interpolation::NearestNeighborInterpolation>();
auto interpolationNullLinear = boost::make_shared<rttb::interpolation::LinearInterpolation>();
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<double> 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<double> 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<double> 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<rttb::WorldCoordinate3D> 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<WorldCoordinate3D>::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<NearestNeighborInterpolation>());
CHECK_NO_THROW(boost::make_shared<LinearInterpolation>());
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;
}