diff --git a/code/core/rttbBaseType.h b/code/core/rttbBaseType.h index 4763c46..618efa8 100644 --- a/code/core/rttbBaseType.h +++ b/code/core/rttbBaseType.h @@ -1,719 +1,719 @@ // ----------------------------------------------------------------------- // 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. @todo Both UnsignedIndex2D and VoxelGridIndex2D required? */ class UnsignedIndex2D: public boost::numeric::ublas::vector - { - public: - UnsignedIndex2D() : boost::numeric::ublas::vector(2) {} + { + 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. @todo Both UnsignedIndex3D and VoxelGridIndex3D required? */ class UnsignedIndex3D: public boost::numeric::ublas::vector - { - public: - UnsignedIndex3D() : boost::numeric::ublas::vector(3) {} + { + 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; - } + :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++) - { - if (gi1(i) != gi2(i)) - { + if (gi1.size()!=gi2.size()) + { return false; - } - } + } - return true; - } + for (int 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; - } - }; + 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) {} + { + public: + WorldCoordinate2D() : boost::numeric::ublas::vector (2) {} WorldCoordinate2D(const WorldCoordinate value) : boost::numeric::ublas::vector(2, value) {} - WorldCoordinate2D(const WorldCoordinate xValue, const WorldCoordinate yValue) + 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(); - } + 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; - } + if (wc1.size()!=wc2.size()) + { + return false; + } - for (int i = 0; i < wc1.size(); i++) - { + for (int i = 0; i < wc1.size(); i++) + { if (wc1(i) != wc2(i)) { - return false; - } - } + return false; + } + } - return true; - } - }; + 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) {} + { + 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; - } + :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(); - } + (*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) + //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); + 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); + 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; - } + return result; + } WorldCoordinate2D getXY() const { - WorldCoordinate2D result; + WorldCoordinate2D result; - result(0) = (*this)(0); - result(1) = (*this)(1); + result(0) = (*this)(0); + result(1) = (*this)(1); - return result; - } + return result; + } - const std::string toString() const - { - std::stringstream ss; - ss << x() << ' ' << y() << ' ' << z(); - return ss.str(); - } + 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); - } + (*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); - } + (*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)); - } + 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)); - } + 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++) - { - if (wc1(i) != wc2(i)) - { + if (wc1.size()!=wc2.size()) + { return false; - } - } + } - return true; - } + for (int i = 0; i < wc1.size(); i++) + { + if (wc1(i)!=wc2(i)) + { + return false; + } + } + + return true; + } friend std::ostream& operator<<(std::ostream& s, const WorldCoordinate3D& aVector) { - s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; - return s; - } + s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; + return s; + } - }; + }; 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. @todo should SpacingVectorType/GridVolumeType be forced to be non-negative? If yes add corresponding tests */ class SpacingVectorType3D: public boost::numeric::ublas::vector - { - public: - SpacingVectorType3D() : boost::numeric::ublas::vector(3) {} + { + 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; - } + (*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(); - } + (*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(); - } + 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); - } + (*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); - } + (*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); - } + (*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++) - { - if (wc1(i) != wc2(i)) - { + if (wc1.size()!=wc2.size()) + { return false; - } - } + } - return true; - } + for (int i = 0; i < wc1.size(); i++) + { + if (wc1(i)!=wc2(i)) + { + return false; + } + } + + return true; + } friend std::ostream& operator<<(std::ostream& s, const SpacingVectorType3D& aVector) { - s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; - return s; - } - }; + s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; + return s; + } + }; /*! @class OrientationMatrix @brief Used to store image orientation information @todo should OrientationMatrix be forced to be non-negative? If yes add corresponding tests */ class OrientationMatrix : public boost::numeric::ublas::matrix - { - public: - /*! The default constructor generates a 3x3 unit 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; - } - } + (*this)(m,m) = 1; + } + } OrientationMatrix(const WorldCoordinate value) : boost::numeric::ublas::matrix(3, 3, value) {} bool equal(const OrientationMatrix& anOrientationMatrix) const { - if (anOrientationMatrix.size1() == (*this).size1()) - { - if (anOrientationMatrix.size2() == (*this).size2()) - { - for (std::size_t m = 0; m < anOrientationMatrix.size1(); m++) + if (anOrientationMatrix.size1()==(*this).size1()) { - for (std::size_t n = 0; n < anOrientationMatrix.size2(); n++) + if (anOrientationMatrix.size2()==(*this).size2()) { - if (!((*this)(m, n) == anOrientationMatrix(m, n))) + for (std::size_t m = 0; m < anOrientationMatrix.size1(); m++) { - return false; - } + for (std::size_t n = 0; n < anOrientationMatrix.size2(); n++) + { + if (!((*this)(m,n) == anOrientationMatrix(m,n))) + { + return false; + } + } + }// end element comparison } - }// end element comparison - } - else - { + else + { + return false; + } + } + else + { return false; - } - } - else - { - return false; - } + } - return true; - } + return true; + } friend bool operator==(const OrientationMatrix& om1, const OrientationMatrix& om2) { return om1.equal(om2); } 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) + for (std::size_t m = 0; m < anOrientationMatrix.size1(); m++) { - s << anOrientationMatrix(m, n); - } - else - { - s << ", " << anOrientationMatrix(m, n); + 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 << " ]"; - } - - s << " ]"; - return s; - } - }; + return s; + } + }; - /*! base for 2D and 3D VoxelIndex; Therefore required beside VoxelGridID + /*! base for 2D and 3D VoxelIndex; Therefore required beside VoxelGridID */ - typedef unsigned int GridIndexType; + 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) {} + { + 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) + VoxelGridIndex3D( const GridIndexType xValue, const GridIndexType yValue, const GridIndexType zValue) : boost::numeric::ublas::vector(3, xValue) { - (*this)(1) = yValue; - (*this)(2) = zValue; - } + (*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(); - } + 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); - } + (*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; - } + if (gi1.size()!=gi2.size()) + { + return false; + } - for (int i = 0; i < gi1.size(); i++) - { + for (int i = 0; i < gi1.size(); i++) + { if (gi1(i) != gi2(i)) { - return false; - } - } + return false; + } + } - return true; - } + return true; + } friend std::ostream& operator<<(std::ostream& s, const VoxelGridIndex3D& aVector) { - s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; - return s; - } - }; + 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) {} + { + public: + VoxelGridIndex2D() : boost::numeric::ublas::vector(2) {} VoxelGridIndex2D(const GridIndexType value) : boost::numeric::ublas::vector(2, value) {} - VoxelGridIndex2D(const GridIndexType xValue, const GridIndexType yValue) + 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(); - } + 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++) - { - if (gi1(i) != gi2(i)) - { + if (gi1.size()!=gi2.size()) + { return false; - } - } + } - return true; - } + for (int 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; - } - }; + 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 int VoxelGridDimensionType; typedef boost::numeric::ublas::vector VoxelGridDimensionVectorType; typedef double FractionType, VoxelSizeType, DVHVoxelNumber; typedef double DoseCalcType, DoseTypeGy, DoseVoxelVolumeType, VolumeType, GridVolumeType; 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 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; + VoxelGridIndex2D index_end; void Init() { x_begin = -1000000; x_end = -1000000; y_begin = -1000000; - y_end = -1000000; + y_end = -1000000; index_begin(0) = 0; - index_begin(1) = 0; - index_end(0) = 0; - index_end(1) = 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 + }//end: namespace rttb #endif diff --git a/code/core/rttbGeometricInfo.h b/code/core/rttbGeometricInfo.h index 3470ea5..140d1c5 100644 --- a/code/core/rttbGeometricInfo.h +++ b/code/core/rttbGeometricInfo.h @@ -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) */ #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. @todo Similarity for core::GeometricInfo needs to be defined. Orientation of the world matrix needs to be the same, but the extend of the grid may differ. This is for example the case if the slice thickness is different in Dose and Mask. This may cause a difference in the grid extend. */ 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) {} 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. @todo Should equality be less restrictive and allow similar Orientation (different values) or equal orientation and different extend in world coordinates/different grid size */ friend bool operator==(const GeometricInfo& gInfo, const GeometricInfo& gInfo1); /*! @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; + VoxelGridIndex3D& aIndex) const; /*! @brief voxel grid index to convert world coordinates. 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; + 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 diff --git a/code/interpolation/CMakeLists.txt b/code/interpolation/CMakeLists.txt index 3aa2474..d511e1c 100644 --- a/code/interpolation/CMakeLists.txt +++ b/code/interpolation/CMakeLists.txt @@ -1,7 +1,12 @@ OPTION(BUILD_MatchPointBinding "Determine if the MatchPoint binding for dose interpolation classes will be generated." OFF) +OPTION(BUILD_ITKBinding + "Determine if the ITK binding for dose interpolation classes will be generated." OFF) IF(BUILD_MatchPointBinding) -ADD_SUBDIRECTORY(MatchPointBinding) + ADD_SUBDIRECTORY(MatchPointBinding) ENDIF(BUILD_MatchPointBinding) +IF(BUILD_ITKBinding) + ADD_SUBDIRECTORY(ITKBinding) +ENDIF(BUILD_ITKBinding) RTTB_CREATE_MODULE(RTTBInterpolation DEPENDS RTTBCore PACKAGE_DEPENDS Boost) \ No newline at end of file diff --git a/code/interpolation/ITKBinding/CMakeLists.txt b/code/interpolation/ITKBinding/CMakeLists.txt new file mode 100644 index 0000000..11c45f1 --- /dev/null +++ b/code/interpolation/ITKBinding/CMakeLists.txt @@ -0,0 +1 @@ +RTTB_CREATE_MODULE(RTTBITKBinding DEPENDS RTTBCore RTTBInterpolation RTTBITKIO PACKAGE_DEPENDS Boost ITK) \ No newline at end of file diff --git a/code/interpolation/ITKBinding/files.cmake b/code/interpolation/ITKBinding/files.cmake new file mode 100644 index 0000000..2121a29 --- /dev/null +++ b/code/interpolation/ITKBinding/files.cmake @@ -0,0 +1,7 @@ +SET(CPP_FILES + rttbITKTransformation.cpp + ) + +SET(H_FILES + rttbITKTransformation.h + ) diff --git a/code/interpolation/ITKBinding/rttbITKTransformation.cpp b/code/interpolation/ITKBinding/rttbITKTransformation.cpp new file mode 100644 index 0000000..4691acc --- /dev/null +++ b/code/interpolation/ITKBinding/rttbITKTransformation.cpp @@ -0,0 +1,91 @@ +// ----------------------------------------------------------------------- +// 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: 484 $ (last changed revision) +// @date $Date: 2014-03-26 16:16:16 +0100 (Mi, 26 Mrz 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ + +#include "rttbITKTransformation.h" +#include "rttbNullPointerException.h" + +namespace rttb +{ + namespace interpolation + { + ITKTransformation::ITKTransformation(const Transform3D3DType* aTransformation): + _pTransformation(aTransformation) + { + //handle null pointer + if (aTransformation == NULL) + { + throw core::NullPointerException("Pointer to registration is NULL."); + } + + } + + void ITKTransformation::convert(const WorldCoordinate3D& aWorldCoordinate, + InputPointType& aInputPoint) const + { + assert(aWorldCoordinate.size() == 3); + assert(aInputPoint.Length == 3); + + for (unsigned int i = 0; i < aInputPoint.Length; ++i) + { + aInputPoint[i] = aWorldCoordinate[i]; + } + } + + void ITKTransformation::convert(const OutputPointType& aOutputPoint, + WorldCoordinate3D& aWorldCoordinate) const + { + assert(aWorldCoordinate.size() == 3); + assert(aOutputPoint.Length == 3); + + for (unsigned int i = 0; i < aOutputPoint.Length; ++i) + { + aWorldCoordinate[i] = aOutputPoint[i]; + } + } + + bool ITKTransformation::transformInverse(const WorldCoordinate3D& + worldCoordinateTarget, WorldCoordinate3D& worldCoordinateMoving) const + { + InputPointType aTargetPoint; + OutputPointType aMovingPoint; + + convert(worldCoordinateTarget, aTargetPoint); + aMovingPoint = _pTransformation->TransformPoint(aTargetPoint); + convert(aMovingPoint, worldCoordinateMoving); + //TransformPoint has no return value... + return true; + } + + bool ITKTransformation::transform(const WorldCoordinate3D& + worldCoordinateMoving, WorldCoordinate3D& worldCoordinateTarget) const + { + OutputPointType aTargetPoint; + InputPointType aMovingPoint; + + convert(worldCoordinateMoving, aMovingPoint); + aTargetPoint = _pTransformation->TransformPoint(aMovingPoint); + convert(aTargetPoint, worldCoordinateTarget); + //TransformPoint has no return value... + return true; + } + + } +} diff --git a/code/interpolation/ITKBinding/rttbITKTransformation.h b/code/interpolation/ITKBinding/rttbITKTransformation.h new file mode 100644 index 0000000..5d682e5 --- /dev/null +++ b/code/interpolation/ITKBinding/rttbITKTransformation.h @@ -0,0 +1,81 @@ +// ----------------------------------------------------------------------- +// 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: 484 $ (last changed revision) +// @date $Date: 2014-03-26 16:16:16 +0100 (Mi, 26 Mrz 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ +#ifndef __ITK_MAPPABLE_DOSE_ACCESSOR_H +#define __ITK_MAPPABLE_DOSE_ACCESSOR_H + +#include + +#include "itkTransform.h" + +#include "rttbDoseAccessorInterface.h" +#include "rttbTransformationInterface.h" + +namespace rttb +{ + namespace interpolation + { + /*! @class ITKMappableDoseAccessor + @brief This class can deal with dose information that has to be transformed into another geometry than the original dose image (transformation specified by ITK transformation object). + */ + class ITKTransformation: public TransformationInterface + { + public: + static const unsigned int InputDimension3D = 3; + static const unsigned int OutputDimension3D = 3; + typedef double TransformScalarType; + typedef itk::Transform Transform3D3DType; + typedef Transform3D3DType::InputPointType InputPointType; + typedef Transform3D3DType::OutputPointType OutputPointType; + typedef boost::shared_ptr Pointer; + + private: + //! Has to be a Pointer type because of inheritance issues with itkSmartPointer (that doesn't recognize the inheritance) + const Transform3D3DType* _pTransformation; + + protected: + void convert(const WorldCoordinate3D& aWorldCoordinate, InputPointType& aInputPoint) const; + void convert(const OutputPointType& aOutputPoint, WorldCoordinate3D& aWorldCoordinate) const; + + public: + /*! @brief Constructor. + @param aRegistration registration given in MatchPoint format (note the pointer format since itkSmartPointer does not support inheritance) + @sa MappableDoseAccessorBase + @pre all input parameters have to be valid + @exception core::NullPointerException if one input parameter is NULL + @exception core::PaddingException if the transformation is undefined and if _acceptPadding==false + */ + ITKTransformation(const Transform3D3DType* aTransformation); + + ~ITKTransformation() {}; + + /*! @brief performs a transformation targetImage --> movingImage + */ + bool transformInverse(const WorldCoordinate3D& worldCoordinateTarget, + WorldCoordinate3D& worldCoordinateMoving) const; + /*! @brief performs a transformation movingImage --> targetImage + */ + bool transform(const WorldCoordinate3D& worldCoordinateMoving, + WorldCoordinate3D& worldCoordinateTarget) const; + }; + } +} + +#endif diff --git a/code/interpolation/MatchPointBinding/rttbMatchPointTransformation.cpp b/code/interpolation/MatchPointBinding/rttbMatchPointTransformation.cpp index a663564..513ca05 100644 --- a/code/interpolation/MatchPointBinding/rttbMatchPointTransformation.cpp +++ b/code/interpolation/MatchPointBinding/rttbMatchPointTransformation.cpp @@ -1,89 +1,88 @@ // ----------------------------------------------------------------------- // 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 "rttbMatchPointTransformation.h" -#include "rttbPaddingException.h" #include "rttbNullPointerException.h" namespace rttb { namespace interpolation { MatchPointTransformation::MatchPointTransformation( const Registration3D3DType* aRegistration): _pRegistration(aRegistration) { //handle null pointers if (aRegistration == NULL) { throw core::NullPointerException("Pointer to registration is NULL."); } } void MatchPointTransformation::convert(const WorldCoordinate3D& aWorldCoordinate, TargetPointType& aTargetPoint) const { assert(aWorldCoordinate.size() == 3); assert(aTargetPoint.Length == 3); for (unsigned int i = 0; i < aTargetPoint.Length; ++i) { aTargetPoint[i] = aWorldCoordinate[i]; } } void MatchPointTransformation::convert(const MovingPointType& aMovingPoint, WorldCoordinate3D& aWorldCoordinate) const { assert(aWorldCoordinate.size() == 3); assert(aMovingPoint.Length == 3); for (unsigned int i = 0; i < aMovingPoint.Length; ++i) { aWorldCoordinate[i] = aMovingPoint[i]; } } bool MatchPointTransformation::transformInverse(const WorldCoordinate3D& worldCoordinateTarget, WorldCoordinate3D& worldCoordinateMoving) const { TargetPointType aTargetPoint; MovingPointType aMovingPoint; convert(worldCoordinateTarget, aTargetPoint); bool ok = _pRegistration->mapPointInverse(aTargetPoint, aMovingPoint); convert(aMovingPoint, worldCoordinateMoving); return ok; } bool MatchPointTransformation::transform(const WorldCoordinate3D& worldCoordinateMoving, WorldCoordinate3D& worldCoordinateTarget) const { TargetPointType aTargetPoint; MovingPointType aMovingPoint; convert(worldCoordinateMoving, aMovingPoint); bool ok = _pRegistration->mapPoint(aMovingPoint, aTargetPoint); convert(aTargetPoint, worldCoordinateTarget); return ok; } } } diff --git a/code/interpolation/files.cmake b/code/interpolation/files.cmake index f4120aa..0851292 100644 --- a/code/interpolation/files.cmake +++ b/code/interpolation/files.cmake @@ -1,18 +1,17 @@ SET(CPP_FILES - rttbMappableDoseAccessorBase.cpp rttbRosuMappableDoseAccessor.cpp rttbSimpleMappableDoseAccessor.cpp rttbInterpolationBase.cpp rttbNearestNeighborInterpolation.cpp rttbLinearInterpolation.cpp ) SET(H_FILES - rttbMappableDoseAccessorBase.h + rttbMappableDoseAccessorInterface.h rttbRosuMappableDoseAccessor.h rttbSimpleMappableDoseAccessor.h rttbInterpolationBase.h rttbNearestNeighborInterpolation.h rttbLinearInterpolation.h rttbTransformationInterface.h ) diff --git a/code/interpolation/rttbInterpolationBase.cpp b/code/interpolation/rttbInterpolationBase.cpp index 0bce389..65f134a 100644 --- a/code/interpolation/rttbInterpolationBase.cpp +++ b/code/interpolation/rttbInterpolationBase.cpp @@ -1,190 +1,191 @@ // ----------------------------------------------------------------------- // 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 "rttbInterpolationBase.h" #include "rttbInvalidParameterException.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace interpolation { void InterpolationBase::setDoseAccessorPointer(const DoseAccessorPointer originalDose) { if (originalDose != NULL) { _spOriginalDose = originalDose; } else { throw core::NullPointerException("originalDose is NULL!"); } }; void InterpolationBase::getNeighborhoodVoxelValues( const WorldCoordinate3D& aWorldCoordinate, unsigned int neighborhood, boost::array& target, boost::shared_array values) const { if (_spOriginalDose == 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 (_spOriginalDose->getGeometricInfo().worldCoordinateToIndex(aWorldCoordinate, aIndex)) { //determine just the nearest voxel to the world coordinate if (neighborhood == 0) { values[0] = _spOriginalDose->getDoseAt(aIndex); } //determine the 8 voxels around the world coordinate else if (neighborhood == 8) { std::vector cornerPoints; WorldCoordinate3D theNextVoxel; _spOriginalDose->getGeometricInfo().indexToWorldCoordinate(aIndex, theNextVoxel); SpacingVectorType3D pixelSpacing = (_spOriginalDose->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); } //now just get the values of all (dose) voxels and store them in values for (int i = 0; i < cornerPoints.size(); ++i) { if (_spOriginalDose->getGeometricInfo().isInside(cornerPoints.at(i))) { values[i] = _spOriginalDose->getDoseAt(cornerPoints.at(i)); } else { //outside value! boundary treatment values[i] = getNearestInsideVoxelValue(cornerPoints.at(i)); } assert(values[i] != -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 (_spOriginalDose->getGeometricInfo().validIndex(voxelChangedXYZ[replacementVoxelIndex])) { return _spOriginalDose->getDoseAt(voxelChangedXYZ[replacementVoxelIndex]); } ++replacementVoxelIndex; } return -1; } }//end namespace core }//end namespace rttb diff --git a/code/interpolation/rttbMappableDoseAccessorBase.cpp b/code/interpolation/rttbMappableDoseAccessorBase.cpp deleted file mode 100644 index afed435..0000000 --- a/code/interpolation/rttbMappableDoseAccessorBase.cpp +++ /dev/null @@ -1,67 +0,0 @@ -// ----------------------------------------------------------------------- -// 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 "rttbMappableDoseAccessorBase.h" -#include "rttbNullPointerException.h" -#include "rttbMappingOutsideOfImageException.h" - -namespace rttb -{ - namespace interpolation - { - MappableDoseAccessorBase::MappableDoseAccessorBase(const core::GeometricInfo& geoInfoTargetImage, - const DoseAccessorPointer doseMovingImage, - const TransformationInterface::Pointer aTransformation, bool acceptPadding, - double defaultOutsideValue): _spOriginalDoseDataMovingImage(doseMovingImage), - _spTransformation(aTransformation), _geoInfoTargetImage(geoInfoTargetImage), - _acceptPadding(acceptPadding), _defaultOutsideValue(defaultOutsideValue) - { - //handle null pointers - if (doseMovingImage == NULL || aTransformation == NULL) - { - throw core::NullPointerException("Pointers to input accessors/transformation cannot be NULL."); - } - } - - DoseTypeGy MappableDoseAccessorBase::getDoseAt(const VoxelGridID aID) const - { - VoxelGridIndex3D aVoxelGridIndex3D; - - if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) - { - return getDoseAt(aVoxelGridIndex3D); - } - else - { - if (_acceptPadding) - { - return _defaultOutsideValue; - } - else - { - throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); - return -1; - } - } - } - - } -} diff --git a/code/interpolation/rttbMappableDoseAccessorBase.h b/code/interpolation/rttbMappableDoseAccessorInterface.h similarity index 67% rename from code/interpolation/rttbMappableDoseAccessorBase.h rename to code/interpolation/rttbMappableDoseAccessorInterface.h index b03a2fd..0d18e33 100644 --- a/code/interpolation/rttbMappableDoseAccessorBase.h +++ b/code/interpolation/rttbMappableDoseAccessorInterface.h @@ -1,99 +1,109 @@ // ----------------------------------------------------------------------- // 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 __MAPPABLE_DOSE_ACCESSOR_BASE_H #define __MAPPABLE_DOSE_ACCESSOR_BASE_H #include #include "rttbDoseAccessorInterface.h" #include "rttbGeometricInfo.h" #include "rttbBaseType.h" #include "rttbTransformationInterface.h" +#include "rttbNullPointerException.h" namespace rttb { namespace interpolation { - /*! @class MappableDoseAccessorBase - @brief Base class for dealing with dose information that has to be transformed into another geometry than the original dose image + /*! @class MappableDoseAccessorInterface + @brief Interface for dealing with dose information that has to be transformed into another geometry than the original dose image @details implementation of strategy is done by derived class (e.g. SimpleMappableDoseAccessor or RosuMappableDoseAccessor. Transformation is defined in TransformationInterface @ingroup interpolation */ - class MappableDoseAccessorBase: public core::DoseAccessorInterface + class MappableDoseAccessorInterface: public core::DoseAccessorInterface { public: - typedef boost::shared_ptr Pointer; + typedef boost::shared_ptr Pointer; protected: DoseAccessorPointer _spOriginalDoseDataMovingImage; TransformationInterface::Pointer _spTransformation; core::GeometricInfo _geoInfoTargetImage; bool _acceptPadding; DoseTypeGy _defaultOutsideValue; public: /*! @brief Constructor. @param geoInfoTargetImage target image geometry @param doseMovingImage dose of moving image @param aTransformation the transformation @param acceptPadding is mapping outside the image allowed @param defaultOutsideValue the default outside voxel value if accepptPadding=true @pre all input parameters have to be valid @exception core::NullPointerException if one input parameter is NULL */ - MappableDoseAccessorBase(const core::GeometricInfo& geoInfoTargetImage, - const DoseAccessorPointer doseMovingImage, const TransformationInterface::Pointer aTransformation, - bool acceptPadding = true, - DoseTypeGy defaultOutsideValue = 0.0); + MappableDoseAccessorInterface(const core::GeometricInfo& geoInfoTargetImage, + const DoseAccessorPointer doseMovingImage, const TransformationInterface::Pointer aTransformation, + bool acceptPadding = true, + DoseTypeGy defaultOutsideValue = 0.0): _spOriginalDoseDataMovingImage(doseMovingImage), + _spTransformation(aTransformation), _geoInfoTargetImage(geoInfoTargetImage), + _acceptPadding(acceptPadding), _defaultOutsideValue(defaultOutsideValue) + { + //handle null pointers + if (doseMovingImage == NULL || aTransformation == NULL) + { + throw core::NullPointerException("Pointers to input accessors/transformation cannot be NULL."); + } + } /*! @brief Virtual destructor of base class */ - virtual ~MappableDoseAccessorBase() {}; + virtual ~MappableDoseAccessorInterface() {}; inline const core::GeometricInfo& getGeometricInfo() const { return _geoInfoTargetImage; }; inline GridSizeType getGridSize() const { return _geoInfoTargetImage.getNumberOfVoxels(); }; /*! @brief Returns the dose for a given VoxelGridID (convenience function that handles conversion VoxelGridID->VoxelGridIndex3D) @sa getDoseAt(const VoxelGridIndex3D& aIndex) */ - DoseTypeGy getDoseAt(const VoxelGridID aID) const; + DoseTypeGy getDoseAt(const VoxelGridID aID) const = 0; /*! @brief Returns the dose for a given VoxelGridIndex3D */ virtual DoseTypeGy getDoseAt(const VoxelGridIndex3D& aIndex) const = 0 ; const IDType getDoseUID() const { return _spOriginalDoseDataMovingImage->getDoseUID(); }; }; } } #endif diff --git a/code/interpolation/rttbRosuMappableDoseAccessor.cpp b/code/interpolation/rttbRosuMappableDoseAccessor.cpp index 22cddab..2ad8876 100644 --- a/code/interpolation/rttbRosuMappableDoseAccessor.cpp +++ b/code/interpolation/rttbRosuMappableDoseAccessor.cpp @@ -1,152 +1,174 @@ // ----------------------------------------------------------------------- // 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): MappableDoseAccessorBase(geoInfoTargetImage, doseMovingImage, + DoseTypeGy defaultOutsideValue): MappableDoseAccessorInterface(geoInfoTargetImage, doseMovingImage, aTransformation, acceptPadding, defaultOutsideValue) { //define linear interpolation InterpolationBase::Pointer interpolationLinear = InterpolationBase::Pointer( new LinearInterpolation()); _spInterpolation = interpolationLinear; _spInterpolation->setDoseAccessorPointer(_spOriginalDoseDataMovingImage); } + DoseTypeGy RosuMappableDoseAccessor::getDoseAt(const VoxelGridID aID) const + { + VoxelGridIndex3D aVoxelGridIndex3D; + + if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) + { + return getDoseAt(aVoxelGridIndex3D); + } + else + { + if (_acceptPadding) + { + return _defaultOutsideValue; + } + else + { + throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); + return -1; + } + } + } + DoseTypeGy RosuMappableDoseAccessor::getDoseAt(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 = std::begin(octants); octantIterator != std::end(octants); ++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.75; xOct += 0.5) + for (double xOct = -0.25; xOct <= 0.25; xOct += 0.5) { - for (double yOct = 0.25; yOct <= 0.75; yOct += 0.5) + for (double yOct = -0.25; yOct <= 0.25; yOct += 0.5) { - for (double zOct = 0.25; zOct <= 0.75; zOct += 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/rttbRosuMappableDoseAccessor.h b/code/interpolation/rttbRosuMappableDoseAccessor.h index 76cf50c..60e3786 100644 --- a/code/interpolation/rttbRosuMappableDoseAccessor.h +++ b/code/interpolation/rttbRosuMappableDoseAccessor.h @@ -1,80 +1,82 @@ // ----------------------------------------------------------------------- // 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 __ROSU_MAPPABLE_DOSE_ACCESSOR_H #define __ROSU_MAPPABLE_DOSE_ACCESSOR_H #include #include "rttbBaseType.h" #include "rttbInterpolationBase.h" #include "rttbTransformationInterface.h" -#include "rttbMappableDoseAccessorBase.h" +#include "rttbMappableDoseAccessorInterface.h" namespace rttb { namespace interpolation { /*! @class RosuMappableDoseAccessor @brief Class for dose mapping based on interpolation described in the Rosu2005 paper @details implementation of the following paper: Rosu, M., Chetty, I. J., Balter, J. M., Kessler, M. L., McShan, D. L., & Ten Haken, R. K. (2005). Dose reconstruction in deforming lung anatomy: Dose grid size effects and clinical implications. Medical Physics, 32(8), 2487. @ingroup interpolation */ - class RosuMappableDoseAccessor: public MappableDoseAccessorBase + class RosuMappableDoseAccessor: public MappableDoseAccessorInterface { private: InterpolationBase::Pointer _spInterpolation; public: typedef boost::shared_ptr Pointer; /*! @brief Constructor. Just hands values over to base class constructor. @note no interpolation as parameter since linear interpolation is fixed. @sa MappableDoseAccessorBase */ RosuMappableDoseAccessor(const core::GeometricInfo& geoInfoTargetImage, const DoseAccessorPointer doseMovingImage, const TransformationInterface::Pointer aTransformation, bool acceptPadding = true, DoseTypeGy defaultOutsideValue = 0.0); /*! @brief Virtual destructor. */ virtual ~RosuMappableDoseAccessor() {}; + DoseTypeGy getDoseAt(const VoxelGridID aID) const; + /*! @brief Returns the dose for a given voxel grid index. The computation of the octant around the voxel is done and the interpolation is performed. @details Boundary treatment: if more than 6 subvoxels are outside: return _defaultOutsideValue. Otherwise: ignore the outside values. @return the dose or if (isOutside==true && _acceptPadding==true) then _defaultValue @exception core::MappingOutsideOfImageException if the point is mapped outside and if _acceptPadding==false, possibly returning _defaultValue) */ DoseTypeGy getDoseAt(const VoxelGridIndex3D& aIndex) const; private: /*! @brief returns the octant coordinates around a coordinate. @details i.e. coordinate is the center of a virtual voxel. Then, each side is divided into equal parts. The centers of the new subvoxels are then returned. @return a vector of the octant coordinates. */ std::vector getOctants(const WorldCoordinate3D& aCoordinate) const; }; } } #endif diff --git a/code/interpolation/rttbSimpleMappableDoseAccessor.cpp b/code/interpolation/rttbSimpleMappableDoseAccessor.cpp index 776cb2c..4a56a74 100644 --- a/code/interpolation/rttbSimpleMappableDoseAccessor.cpp +++ b/code/interpolation/rttbSimpleMappableDoseAccessor.cpp @@ -1,100 +1,122 @@ // ----------------------------------------------------------------------- // 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): MappableDoseAccessorBase(geoInfoTargetImage, doseMovingImage, + 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->setDoseAccessorPointer(_spOriginalDoseDataMovingImage); } } + DoseTypeGy SimpleMappableDoseAccessor::getDoseAt(const VoxelGridID aID) const + { + VoxelGridIndex3D aVoxelGridIndex3D; + + if (_geoInfoTargetImage.convert(aID, aVoxelGridIndex3D)) + { + return getDoseAt(aVoxelGridIndex3D); + } + else + { + if (_acceptPadding) + { + return _defaultOutsideValue; + } + else + { + throw core::MappingOutsideOfImageException("Error in conversion from index to world coordinates"); + return -1; + } + } + } + DoseTypeGy SimpleMappableDoseAccessor::getDoseAt(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/code/interpolation/rttbSimpleMappableDoseAccessor.h b/code/interpolation/rttbSimpleMappableDoseAccessor.h index a0f2959..e9ae94a 100644 --- a/code/interpolation/rttbSimpleMappableDoseAccessor.h +++ b/code/interpolation/rttbSimpleMappableDoseAccessor.h @@ -1,71 +1,75 @@ // ----------------------------------------------------------------------- // 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 __SIMPLE_MAPPABLE_DOSE_ACCESSOR_H #define __SIMPLE_MAPPABLE_DOSE_ACCESSOR_H #include #include "rttbBaseType.h" #include "rttbInterpolationBase.h" #include "rttbTransformationInterface.h" -#include "rttbMappableDoseAccessorBase.h" +#include "rttbMappableDoseAccessorInterface.h" namespace rttb { namespace interpolation { /*! @class SimpleMappableDoseAccessor @brief Class for dose mapping based on simple trilinear interpolation @ingroup interpolation */ - class SimpleMappableDoseAccessor: public MappableDoseAccessorBase + class SimpleMappableDoseAccessor: public MappableDoseAccessorInterface { private: InterpolationBase::Pointer _spInterpolation; public: typedef boost::shared_ptr Pointer; /*! @brief Constructor. Just hands values over to base class constructor. @param aInterpolation the used interpolation. @sa MappableDoseAccessorBase */ SimpleMappableDoseAccessor(const core::GeometricInfo& geoInfoTargetImage, const DoseAccessorPointer doseMovingImage, const TransformationInterface::Pointer aTransformation, const InterpolationBase::Pointer aInterpolation, bool acceptPadding = true, DoseTypeGy defaultOutsideValue = 0.0); /*! @brief Virtual destructor of class */ virtual ~SimpleMappableDoseAccessor() {}; + /*! @brief Returns the dose for a given voxel grid id. Plain trilinear interpolation is performed. + @sa getDoseAt(const VoxelGridIndex3D& aIndex) + */ + DoseTypeGy getDoseAt(const VoxelGridID aID) const; /*! @brief Returns the dose for a given voxel grid index. Plain trilinear interpolation is performed. @return the dose or if (isOutside==true && _acceptPadding==true) then _defaultValue @exception core::MappingOutsideOfImageException if the point is mapped outside and if _acceptPadding==false, possibly returning _defaultValue) */ DoseTypeGy getDoseAt(const VoxelGridIndex3D& aIndex) const; }; } } #endif diff --git a/code/io/itk/rttbDoseAccessorProcessorBase.h b/code/io/itk/rttbDoseAccessorProcessorBase.h index 9bbd057..1f07a32 100644 --- a/code/io/itk/rttbDoseAccessorProcessorBase.h +++ b/code/io/itk/rttbDoseAccessorProcessorBase.h @@ -1,58 +1,61 @@ // ----------------------------------------------------------------------- // 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 __DOSE_ACCESSOR_PROCESSOR_BASE_H #define __DOSE_ACCESSOR_PROCESSOR_BASE_H #include #include "rttbDoseAccessorProcessorInterface.h" -namespace rttb{ +namespace rttb +{ namespace core - { + { /*! @class DoseAccessorProcessorBase @brief Abstract class for all DoseAccessor generating classes */ class DoseAccessorProcessorBase: public DoseAccessorProcessorInterface + { + public: + typedef core::DoseAccessorInterface::DoseAccessorPointer DoseAccessorPointer; + + virtual void setDoseAccessor(DoseAccessorPointer accessor) { - public: - typedef core::DoseAccessorInterface::DoseAccessorPointer DoseAccessorPointer; - - virtual void setDoseAccessor(DoseAccessorPointer accessor) - { - _doseAccessor = accessor; - }; - - private: - DoseAccessorProcessorBase(const DoseAccessorProcessorBase&); //not implemented on purpose -> non-copyable - DoseAccessorProcessorBase& operator=(const DoseAccessorProcessorBase&);//not implemented on purpose -> non-copyable - - protected: - DoseAccessorProcessorBase(){}; - virtual ~DoseAccessorProcessorBase(){}; - - /*! @brief Dose accessor which should be generated */ - DoseAccessorPointer _doseAccessor; + _doseAccessor = accessor; }; - } + + private: + DoseAccessorProcessorBase(const + DoseAccessorProcessorBase&); //not implemented on purpose -> non-copyable + DoseAccessorProcessorBase& operator=(const + DoseAccessorProcessorBase&);//not implemented on purpose -> non-copyable + + protected: + DoseAccessorProcessorBase() {}; + virtual ~DoseAccessorProcessorBase() {}; + + /*! @brief Dose accessor which should be generated */ + DoseAccessorPointer _doseAccessor; + }; } +} #endif diff --git a/code/io/itk/rttbITKImageDoseAccessor.cpp b/code/io/itk/rttbITKImageDoseAccessor.cpp index 182adfe..51870de 100644 --- a/code/io/itk/rttbITKImageDoseAccessor.cpp +++ b/code/io/itk/rttbITKImageDoseAccessor.cpp @@ -1,133 +1,119 @@ // ----------------------------------------------------------------------- // 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 "rttbITKImageDoseAccessor.h" #include "rttbException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace io { namespace itk { ITKImageDoseAccessor::~ITKImageDoseAccessor() { } ITKImageDoseAccessor::ITKImageDoseAccessor(ITKDoseImageType::ConstPointer doseImage): _doseGridScaling(1.0) { _dose = doseImage; if (_dose.IsNull()) { throw core::InvalidDoseException("Dose image = 0!") ; } - _geoInfo = convertToGeometricInfo(_dose); + assembleGeometricInfo(); } DoseTypeGy ITKImageDoseAccessor::getDoseAt(const VoxelGridID aID) const { VoxelGridIndex3D aVoxelGridIndex; if (_geoInfo.convert(aID, aVoxelGridIndex)) { return getDoseAt(aVoxelGridIndex); } else { return -1; } } DoseTypeGy ITKImageDoseAccessor::getDoseAt(const VoxelGridIndex3D& aIndex) const { if (_geoInfo.validIndex(aIndex)) { const ITKDoseImageType::IndexType pixelIndex = {{aIndex[0], aIndex[1], aIndex[2]}}; return _dose->GetPixel(pixelIndex) * _doseGridScaling; } else { return -1; } } - bool isValid(const core::GeometricInfo& geoInfo) + bool ITKImageDoseAccessor::assembleGeometricInfo() { - if (geoInfo.getSpacing()[0] == 0 || geoInfo.getSpacing()[1] == 0 || geoInfo.getSpacing()[2] == 0) - { - throw core::InvalidDoseException("Dose spacing = 0!") ; - return false; - } + _geoInfo.setSpacing(SpacingVectorType3D(_dose->GetSpacing()[0], _dose->GetSpacing()[1], + _dose->GetSpacing()[2])); - if (geoInfo.getNumRows() == 0 || geoInfo.getNumColumns() == 0 || geoInfo.getNumSlices() == 0) + if (_geoInfo.getSpacing()[0] == 0 || _geoInfo.getSpacing()[1] == 0 || _geoInfo.getSpacing()[2] == 0) { - throw core::InvalidDoseException("Dose row/col/slices = 0!") ; - return false; + throw core::InvalidDoseException("Pixel spacing = 0!"); } - return true; - } - - core::GeometricInfo convertToGeometricInfo(const ITKImageBaseType* image) - { - core::GeometricInfo geoInfo; - geoInfo.setSpacing(SpacingVectorType3D(image->GetSpacing()[0], image->GetSpacing()[1], - image->GetSpacing()[2])); - geoInfo.setImagePositionPatient(WorldCoordinate3D(image->GetOrigin()[0], image->GetOrigin()[1], - image->GetOrigin()[2])); + _geoInfo.setImagePositionPatient(WorldCoordinate3D(_dose->GetOrigin()[0], _dose->GetOrigin()[1], + _dose->GetOrigin()[2])); OrientationMatrix OM(0); for (int col = 0; col < 3; ++col) { for (int row = 0; row < 3; ++row) { - OM(col, row) = image->GetDirection()(col, row); + OM(col, row) = _dose->GetDirection()(col, row); } } - geoInfo.setOrientationMatrix(OM); - geoInfo.setNumColumns(image->GetLargestPossibleRegion().GetSize()[0]); - geoInfo.setNumRows(image->GetLargestPossibleRegion().GetSize()[1]); - geoInfo.setNumSlices(image->GetLargestPossibleRegion().GetSize()[2]); + _geoInfo.setOrientationMatrix(OM); + _geoInfo.setNumColumns(_dose->GetLargestPossibleRegion().GetSize()[0]); + _geoInfo.setNumRows(_dose->GetLargestPossibleRegion().GetSize()[1]); + _geoInfo.setNumSlices(_dose->GetLargestPossibleRegion().GetSize()[2]); - if (isValid(geoInfo)) + if (_geoInfo.getNumColumns() == 0 || _geoInfo.getNumRows() == 0 || _geoInfo.getNumSlices() == 0) { - return geoInfo; - } - else - { - throw core::InvalidDoseException("no valid GeometricInfo after conversion from itk"); - return core::GeometricInfo(); + throw core::InvalidDoseException("Empty dicom dose!") ; } + + return true; + } } } } diff --git a/code/io/itk/rttbITKImageDoseAccessor.h b/code/io/itk/rttbITKImageDoseAccessor.h index ded44fe..827795a 100644 --- a/code/io/itk/rttbITKImageDoseAccessor.h +++ b/code/io/itk/rttbITKImageDoseAccessor.h @@ -1,103 +1,97 @@ // ----------------------------------------------------------------------- // 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 __ITK_DOSE_ACCESSOR_H #define __ITK_DOSE_ACCESSOR_H #include #include #include "rttbDoseAccessorInterface.h" #include "rttbGeometricInfo.h" #include "rttbBaseType.h" #include "itkImage.h" namespace rttb { namespace io { namespace itk { typedef ::itk::Image ITKDoseImageType; typedef ::itk::ImageBase<3> ITKImageBaseType; /*! @class ITKImageDoseAccessor @brief This class gives access to dose information stored in an itk image @details _doseGridScaling is always 1.0. Thus, it is assumed that the values in the itkImage are absolute. */ class ITKImageDoseAccessor: public core::DoseAccessorInterface { private: /** @brief The dose as itkImage */ ITKDoseImageType::ConstPointer _dose; IDType _doseUID; /** @brief The dosegridscaling * @note is always 1.0 */ double _doseGridScaling; /*! @brief constructor @exception InvalidDoseException if _dose is NULL */ ITKImageDoseAccessor(); + /*! @brief get all required data from the itk image contained in _dose + @exception InvalidDoseException if PixelSpacing is 0 or size in any dimension is 0. + */ + bool assembleGeometricInfo(); + public: ~ITKImageDoseAccessor(); /*! @brief Constructor. Initialization with a itk image containing the dose @pre doseImage must be a valid instance (and not null) @note the doseImage pixels are dose (i.e. _doseGridScaling=1.0 is assumed always) */ ITKImageDoseAccessor(ITKDoseImageType::ConstPointer doseImage); /*! @brief returns the dose for an id */ DoseTypeGy getDoseAt(const VoxelGridID aID) const; /*! @brief returns the dose for an index */ DoseTypeGy getDoseAt(const VoxelGridIndex3D& aIndex) const; const IDType getDoseUID() const { return _doseUID; }; }; - - /*! @brief get all required data from the itk image contained in _dose - */ - core::GeometricInfo convertToGeometricInfo(const ITKImageBaseType* image); - - /*! @brief check if GeometricInfo is valid - @details e.g. if spacing != 0 and size[i] != 0 - @exception InvalidDoseException - @todo add as member function in GeometricInfo? - */ - bool isValid(const core::GeometricInfo& geoInfo); } } } #endif diff --git a/code/io/other/rttbDVHTxtFileReader.cpp b/code/io/other/rttbDVHTxtFileReader.cpp index 76e478b..a126a4a 100644 --- a/code/io/other/rttbDVHTxtFileReader.cpp +++ b/code/io/other/rttbDVHTxtFileReader.cpp @@ -1,196 +1,225 @@ // ----------------------------------------------------------------------- // 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 "rttbDVHTxtFileReader.h" #include "rttbInvalidParameterException.h" #include "rttbBaseType.h" -namespace rttb{ - namespace io{ - namespace other{ - - DVHTxtFileReader::DVHTxtFileReader(FileNameString aFileName){ - _fileName=aFileName; - _resetFile=true; +namespace rttb +{ + namespace io + { + namespace other + { + + DVHTxtFileReader::DVHTxtFileReader(FileNameString aFileName) + { + _fileName = aFileName; + _resetFile = true; } - DVHTxtFileReader::~DVHTxtFileReader(){} + DVHTxtFileReader::~DVHTxtFileReader() {} - void DVHTxtFileReader::resetFileName(FileNameString aFileName){ - _fileName=aFileName; - _resetFile=true; + void DVHTxtFileReader::resetFileName(FileNameString aFileName) + { + _fileName = aFileName; + _resetFile = true; } - void DVHTxtFileReader::createDVH(){ - std::ifstream dvh_ifstr(this->_fileName.c_str(),std::ios::in); + void DVHTxtFileReader::createDVH() + { + std::ifstream dvh_ifstr(this->_fileName.c_str(), std::ios::in); std::string structureLabel; std::string dvhType; int numberOfBins; DoseTypeGy prescribedDose; double _estimated_max_dose_prescribed_dose_ratio; int voxelsInStructure; std::deque dataDifferential; std::deque dataCumulative; - DoseTypeGy deltaD=0; - DoseVoxelVolumeType deltaV=0; + DoseTypeGy deltaD = 0; + DoseVoxelVolumeType deltaV = 0; IDType strID; IDType doseID; - if ( !dvh_ifstr.is_open() ) + if (!dvh_ifstr.is_open()) { throw core::InvalidParameterException("DVH file name invalid: could not open the dvh file!"); } - else { - bool data_begin=false; + else + { + bool data_begin = false; - while(!dvh_ifstr.eof()) + while (!dvh_ifstr.eof()) { std::string line; - std::getline(dvh_ifstr,line); + std::getline(dvh_ifstr, line); + + if (line.find("DVH Data:") != std::string::npos) + { + data_begin = true; + } - if(line.find("DVH Data:")!= std::string::npos) - data_begin=true; - if(data_begin && (line.find(",")!= std::string::npos)) - { - std::stringstream ss(line.substr(line.find(",")+1)); + if (data_begin && (line.find(",") != std::string::npos)) + { + std::stringstream ss(line.substr(line.find(",") + 1)); DoseCalcType dvh_i; ss >> dvh_i; - if(dvhType=="DIFFERENTIAL") + + if (dvhType == "DIFFERENTIAL") { dataDifferential.push_back(dvh_i); } - else if(dvhType=="CUMULATIVE") + else if (dvhType == "CUMULATIVE") { dataCumulative.push_back(dvh_i); } } - if(line.find("DeltaD: ")!= std::string::npos) + + if (line.find("DeltaD: ") != std::string::npos) { std::stringstream ss(line.substr(8)); ss >> deltaD; } - if(line.find("DeltaV: ")!= std::string::npos) + + if (line.find("DeltaV: ") != std::string::npos) { std::stringstream ss(line.substr(8)); ss >> deltaV; } - if(line.find("StructureID: ")!= std::string::npos) + + if (line.find("StructureID: ") != std::string::npos) { std::stringstream ss(line.substr(13)); ss >> strID; } - if(line.find("DoseID: ")!= std::string::npos) + + if (line.find("DoseID: ") != std::string::npos) { std::stringstream ss(line.substr(8)); ss >> doseID; } - if(line.find("Number of bins: ")!= std::string::npos) + + if (line.find("Number of bins: ") != std::string::npos) { std::stringstream ss(line.substr(16)); ss >> numberOfBins; } - if(line.find("Structure Label: ")!= std::string::npos) + + if (line.find("Structure Label: ") != std::string::npos) { std::stringstream ss(line.substr(17)); ss >> structureLabel; } - if(line.find("DVH Type: ")!= std::string::npos) + + if (line.find("DVH Type: ") != std::string::npos) { std::stringstream ss(line.substr(10)); ss >> dvhType; } - if(line.find("Prescribed Dose: ")!=std::string::npos) + + if (line.find("Prescribed Dose: ") != std::string::npos) { std::stringstream ss(line.substr(17)); ss >> prescribedDose; } - if(line.find("Estimated_max_dose_prescribed_dose_ratio: ")!=std::string::npos) + + if (line.find("Estimated_max_dose_prescribed_dose_ratio: ") != std::string::npos) { std::stringstream ss(line.substr(42)); ss >> _estimated_max_dose_prescribed_dose_ratio; } - if(line.find("Voxels In Structure: ")!=std::string::npos) + + if (line.find("Voxels In Structure: ") != std::string::npos) { std::stringstream ss(line.substr(21)); ss >> voxelsInStructure; } } } - numberOfBins=std::max(dataDifferential.size(),dataCumulative.size()); + numberOfBins = std::max(dataDifferential.size(), dataCumulative.size()); - if(dvhType=="CUMULATIVE") + if (dvhType == "CUMULATIVE") { DoseCalcType differentialDVHi = 0; std::deque::iterator it; - for(it=dataCumulative.begin();it!=dataCumulative.end();it++ ) + for (it = dataCumulative.begin(); it != dataCumulative.end(); it++) { - if(dataDifferential.size()==numberOfBins-1) + if (dataDifferential.size() == numberOfBins - 1) { - differentialDVHi=*it; + differentialDVHi = *it; } - else{ - differentialDVHi = *it-*(it+1); + else + { + differentialDVHi = *it - *(it + 1); } dataDifferential.push_back(differentialDVHi); } } - if(numberOfBins==0) + + if (numberOfBins == 0) { throw core::InvalidParameterException("Invalid dvh file: empty dvh data!"); } - if(deltaD==0) + + if (deltaD == 0) { - deltaD=prescribedDose*_estimated_max_dose_prescribed_dose_ratio/numberOfBins; + deltaD = prescribedDose * _estimated_max_dose_prescribed_dose_ratio / numberOfBins; } - if(deltaV==0) + + if (deltaV == 0) { - deltaV=0.027; + deltaV = 0.027; } - if(deltaD==0 || deltaV==0) + + if (deltaD == 0 || deltaV == 0) { throw core::InvalidParameterException("Invalid dvh file: deltaD or deltaV must not be zero!"); } - _dvh=boost::make_shared(dataDifferential,deltaD,deltaV,strID,doseID); - _resetFile=false; + _dvh = boost::make_shared(dataDifferential, deltaD, deltaV, strID, doseID); + _resetFile = false; } - DVHTxtFileReader::DVHPointer DVHTxtFileReader::generateDVH(){ - if(_resetFile){ + DVHTxtFileReader::DVHPointer DVHTxtFileReader::generateDVH() + { + if (_resetFile) + { this->createDVH(); } + return _dvh; } }//end namepsace other }//end namespace io }//end namespace rttb diff --git a/testing/interpolation/CMakeLists.txt b/testing/interpolation/CMakeLists.txt index d74fd91..816a47b 100644 --- a/testing/interpolation/CMakeLists.txt +++ b/testing/interpolation/CMakeLists.txt @@ -1,19 +1,26 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(INTERPOLATION_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbInterpolationTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(SimpleMappableDoseAccessorTest ${INTERPOLATION_TESTS} SimpleMappableDoseAccessorTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm") ADD_TEST(RosuMappableDoseAccessorTest ${INTERPOLATION_TESTS} RosuMappableDoseAccessorTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm") ADD_TEST(InterpolationTest ${INTERPOLATION_TESTS} InterpolationTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm") -RTTB_CREATE_TEST_MODULE(rttbInterpolation DEPENDS RTTBInterpolation RTTBMatchPointBinding RTTBDicomIO RTTBITKIO RTTBHelaxIO PACKAGE_DEPENDS Boost Litmus MatchPoint DCMTK) +IF(BUILD_ITKBinding) + ADD_SUBDIRECTORY(ITKBinding) +ENDIF(BUILD_ITKBinding) +IF(BUILD_MatchPointBinding) + ADD_SUBDIRECTORY(MatchPointBinding) +ENDIF(BUILD_MatchPointBinding) + +RTTB_CREATE_TEST_MODULE(rttbInterpolation DEPENDS RTTBInterpolation RTTBDicomIO PACKAGE_DEPENDS Litmus) diff --git a/testing/interpolation/DummyTransformation.cpp b/testing/interpolation/DummyTransformation.cpp new file mode 100644 index 0000000..5bc54c4 --- /dev/null +++ b/testing/interpolation/DummyTransformation.cpp @@ -0,0 +1,42 @@ +// ----------------------------------------------------------------------- +// 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: 3987 $ (last changed revision) +// @date $Date: 2012-10-08 09:50:17 +0200 (Mo, 08 Okt 2012) $ (last change date) +// @author $Author: mang $ (last changed by) +*/ + +#include "DummyTransformation.h" + +namespace rttb +{ + namespace testing + { + bool DummyTransformation::transformInverse(const WorldCoordinate3D& worldCoordinateTarget, + WorldCoordinate3D& worldCoordinateMoving) const + { + worldCoordinateMoving = worldCoordinateTarget; + return true; + } + + bool DummyTransformation::transform(const WorldCoordinate3D& worldCoordinateMoving, + WorldCoordinate3D& worldCoordinateTarget) const + { + worldCoordinateTarget = worldCoordinateMoving; + return true; + } + } +} diff --git a/testing/interpolation/DummyTransformation.h b/testing/interpolation/DummyTransformation.h new file mode 100644 index 0000000..41c9ac2 --- /dev/null +++ b/testing/interpolation/DummyTransformation.h @@ -0,0 +1,54 @@ +// ----------------------------------------------------------------------- +// 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: 3987 $ (last changed revision) +// @date $Date: 2012-10-08 09:50:17 +0200 (Mo, 08 Okt 2012) $ (last change date) +// @author $Author: mang $ (last changed by) +*/ +#ifndef __DUMMY_TRANSFORMATION_H +#define __DUMMY_TRANSFORMATION_H + +#include "rttbTransformationInterface.h" + +namespace rttb +{ + namespace testing + { + /*! @class DummyTransformation + @brief implements a dummy transformation (return the points as they were given) + */ + class DummyTransformation : public interpolation::TransformationInterface + { + public: + /*! @brief Constructor + */ + DummyTransformation() {}; + + /*! @brief target equals to moving + */ + bool transformInverse(const WorldCoordinate3D& worldCoordinateTarget, + WorldCoordinate3D& worldCoordinateMoving) const; + /*! @brief moving equals to target + */ + bool transform(const WorldCoordinate3D& worldCoordinateMoving, + WorldCoordinate3D& worldCoordinateTarget) const; + + }; + + } +} + +#endif \ No newline at end of file diff --git a/testing/interpolation/ITKBinding/CMakeLists.txt b/testing/interpolation/ITKBinding/CMakeLists.txt new file mode 100644 index 0000000..40d69e0 --- /dev/null +++ b/testing/interpolation/ITKBinding/CMakeLists.txt @@ -0,0 +1,17 @@ +#----------------------------------------------------------------------------- +# Setup the system information test. Write out some basic failsafe +# information in case the test doesn't run. +#----------------------------------------------------------------------------- + + +SET(INTERPOLATION_ITK_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbInterpolationITKTests) + +SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) + +SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) + + +#----------------------------------------------------------------------------- +ADD_TEST(SimpleMappableDoseAccessorWithITKTest ${INTERPOLATION_ITK_TESTS} SimpleMappableDoseAccessorWithITKTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm") + +RTTB_CREATE_TEST_MODULE(rttbInterpolationITK DEPENDS RTTBInterpolation RTTBITKBinding RTTBDicomIO PACKAGE_DEPENDS Litmus ITK) diff --git a/testing/interpolation/ITKBinding/SimpleMappableDoseAccessorWithITKTest.cpp b/testing/interpolation/ITKBinding/SimpleMappableDoseAccessorWithITKTest.cpp new file mode 100644 index 0000000..eeb7a39 --- /dev/null +++ b/testing/interpolation/ITKBinding/SimpleMappableDoseAccessorWithITKTest.cpp @@ -0,0 +1,186 @@ +// ----------------------------------------------------------------------- +// 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: 741 $ (last changed revision) +// @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ + +#include "litCheckMacros.h" + +#include "rttbBaseType.h" +#include "rttbGenericDoseIterator.h" +#include "rttbDicomFileDoseAccessorGenerator.h" +#include "rttbSimpleMappableDoseAccessor.h" +#include "rttbNearestNeighborInterpolation.h" +#include "rttbLinearInterpolation.h" +#include "rttbTransformationInterface.h" +#include "rttbITKTransformation.h" + +#include "rttbNullPointerException.h" + +#include "itkTranslationTransform.h" + +namespace rttb +{ + namespace testing + { + typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; + typedef rttb::interpolation::SimpleMappableDoseAccessor SimpleMappableDoseAccessor; + typedef rttb::interpolation::NearestNeighborInterpolation NearestNeighborInterpolation; + typedef rttb::interpolation::LinearInterpolation LinearInterpolation; + typedef rttb::interpolation::TransformationInterface TransformationInterface; + typedef rttb::interpolation::ITKTransformation ITKTransformation; + typedef itk::TranslationTransform TranslationTransformType; + + /*! @brief SimpleMappableDoseAccessorWithITKTest - test the API of SimpleMappableDoseAccessor with ITK transformation + 1) Test constructor + 2) test getDoseAt() + a) with Identity transform + b) with translation transform + */ + + int SimpleMappableDoseAccessorWithITKTest(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 SimpleMappableDoseAccessorWithITKTest are expected" << + std::endl; + return -1; + } + + const double doseGridScaling = 2.822386e-005; + + rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1( + RTDOSE_FILENAME_CONSTANT_TWO.c_str()); + DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); + + rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator2( + RTDOSE_FILENAME_INCREASE_X.c_str()); + DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); + + core::GeometricInfo doseAccessor1GeometricInfo = doseAccessor1->getGeometricInfo(); + + auto interpolationNN = + NearestNeighborInterpolation::Pointer(new NearestNeighborInterpolation()); + auto interpolationLinear = boost::shared_ptr + (new LinearInterpolation()); + + //auto transformITKIdentity = TransformationInterface::Pointer(new ITKTransformation( + // IdentityTransformType::New())); + + TranslationTransformType::Pointer transformITKIdentityTemporary = + TranslationTransformType::New(); + TranslationTransformType::OutputVectorType translationIdentity; + translationIdentity[0] = 0.0; + translationIdentity[1] = 0.0; + translationIdentity[2] = 0.0; + transformITKIdentityTemporary->Translate(translationIdentity); + + auto transformITKIdentity = TransformationInterface::Pointer(new ITKTransformation( + transformITKIdentityTemporary)); + + TranslationTransformType::Pointer transformITKTranslationTemporary = + TranslationTransformType::New(); + TranslationTransformType::OutputVectorType translation; + translation[0] = 5.0; + translation[1] = 5.0; + translation[2] = 5.0; + transformITKTranslationTemporary->Translate(translation); + + auto transformITKTranslation = TransformationInterface::Pointer(new ITKTransformation( + transformITKTranslationTemporary)); + + auto aSimpleMappableDoseAccessorITKIdentity = new SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformITKIdentity, interpolationLinear); + auto aSimpleMappableDoseAccessorITKTranslation = new SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformITKTranslation, interpolationLinear); + + //1) Test constructor + CHECK_NO_THROW(SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformITKIdentity, interpolationLinear)); + CHECK_NO_THROW(SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformITKTranslation, interpolationLinear)); + CHECK_NO_THROW(ITKTransformation(transformITKTranslationTemporary.GetPointer())); + CHECK_THROW_EXPLICIT(ITKTransformation(NULL), core::NullPointerException); + + //2) test getDoseAt() + // a) with Identity transform + + //test valid voxels + std::vector voxelsAsIndexToTest; + std::vector expectedValues; + + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(5, 9, 8)); + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(22, 15, 10)); + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(30, 20, 7)); + + expectedValues.push_back(5.0 * doseGridScaling); + expectedValues.push_back(22.0 * doseGridScaling); + expectedValues.push_back(30.0 * doseGridScaling); + + //convert VoxelGridIndex3D to VoxelGridID + for (int i = 0; i < voxelsAsIndexToTest.size(); i++) + { + VoxelGridID currentId; + doseAccessor1GeometricInfo.convert(voxelsAsIndexToTest.at(i), currentId); + //test if the expected interpolation values are computed + CHECK_EQUAL(aSimpleMappableDoseAccessorITKIdentity->getDoseAt(voxelsAsIndexToTest.at(i)), + expectedValues.at(i)); + //test if getDoseAt(VoxelGridIndex3D) and getDoseAt(VoxelGridD) lead to the same results + CHECK_EQUAL(aSimpleMappableDoseAccessorITKIdentity->getDoseAt(voxelsAsIndexToTest.at(i)), + aSimpleMappableDoseAccessorITKIdentity->getDoseAt(currentId)); + } + + + //no tests with invalid IDs, this has been done already tested in SimpleMappableDoseAccessorTest + + //b) with translation transform + + //translation of one voxel in each direction + expectedValues.clear(); + expectedValues.push_back(6.0 * doseGridScaling); + expectedValues.push_back(23.0 * doseGridScaling); + expectedValues.push_back(31.0 * doseGridScaling); + + for (int i = 0; i < voxelsAsIndexToTest.size(); i++) + { + VoxelGridID currentId; + doseAccessor1GeometricInfo.convert(voxelsAsIndexToTest.at(i), currentId); + //test if the expected interpolation values are computed + CHECK_EQUAL(aSimpleMappableDoseAccessorITKTranslation->getDoseAt(voxelsAsIndexToTest.at(i)), + expectedValues.at(i)); + //test if getDoseAt(VoxelGridIndex3D) and getDoseAt(VoxelGridD) lead to the same results + CHECK_EQUAL(aSimpleMappableDoseAccessorITKTranslation->getDoseAt(voxelsAsIndexToTest.at(i)), + aSimpleMappableDoseAccessorITKTranslation->getDoseAt(currentId)); + } + + + RETURN_AND_REPORT_TEST_SUCCESS; + } + + }//end namespace testing +}//end namespace rttb diff --git a/testing/interpolation/ITKBinding/files.cmake b/testing/interpolation/ITKBinding/files.cmake new file mode 100644 index 0000000..91a4abc --- /dev/null +++ b/testing/interpolation/ITKBinding/files.cmake @@ -0,0 +1,7 @@ +SET(CPP_FILES + SimpleMappableDoseAccessorWithITKTest.cpp + rttbITKBindingTests.cpp + ) + +SET(H_FILES + ) diff --git a/testing/interpolation/ITKBinding/rttbITKBindingTests.cpp b/testing/interpolation/ITKBinding/rttbITKBindingTests.cpp new file mode 100644 index 0000000..9695ef5 --- /dev/null +++ b/testing/interpolation/ITKBinding/rttbITKBindingTests.cpp @@ -0,0 +1,63 @@ +// ----------------------------------------------------------------------- +// 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(SimpleMappableDoseAccessorWithITKTest); + } + } +} + +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; +} diff --git a/testing/interpolation/MatchPointBinding/CMakeLists.txt b/testing/interpolation/MatchPointBinding/CMakeLists.txt new file mode 100644 index 0000000..6b6ff81 --- /dev/null +++ b/testing/interpolation/MatchPointBinding/CMakeLists.txt @@ -0,0 +1,17 @@ +#----------------------------------------------------------------------------- +# Setup the system information test. Write out some basic failsafe +# information in case the test doesn't run. +#----------------------------------------------------------------------------- + + +SET(INTERPOLATION_MP_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbInterpolationMPTests) + +SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) + +SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) + + +#----------------------------------------------------------------------------- +ADD_TEST(SimpleMappableDoseAccessorWithMatchPointTest ${INTERPOLATION_MP_TESTS} SimpleMappableDoseAccessorWithMatchPointTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm") + +RTTB_CREATE_TEST_MODULE(rttbInterpolationMP DEPENDS RTTBInterpolation RTTBMatchPointBinding RTTBDicomIO PACKAGE_DEPENDS Litmus MatchPoint) diff --git a/testing/interpolation/MatchPointBinding/SimpleMappableDoseAccessorWithMatchPointTest.cpp b/testing/interpolation/MatchPointBinding/SimpleMappableDoseAccessorWithMatchPointTest.cpp new file mode 100644 index 0000000..ac38381 --- /dev/null +++ b/testing/interpolation/MatchPointBinding/SimpleMappableDoseAccessorWithMatchPointTest.cpp @@ -0,0 +1,248 @@ +// ----------------------------------------------------------------------- +// 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: 741 $ (last changed revision) +// @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ + +#include "litCheckMacros.h" + +#include "rttbBaseType.h" +#include "rttbGenericDoseIterator.h" +#include "rttbDicomFileDoseAccessorGenerator.h" +#include "rttbSimpleMappableDoseAccessor.h" +#include "rttbNearestNeighborInterpolation.h" +#include "rttbLinearInterpolation.h" +#include "rttbGeometricInfo.h" +#include "rttbTransformationInterface.h" +#include "rttbMatchPointTransformation.h" + +#include "rttbNullPointerException.h" + +#include "registrationTest.h" +#include "simpleRegistrationWorkflow.h" + +namespace rttb +{ + namespace testing + { + static const unsigned int TargetDimension3D = 3; + static const unsigned int MovingDimension3D = 3; + typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; + typedef rttb::interpolation::SimpleMappableDoseAccessor SimpleMappableDoseAccessor; + typedef map::core::RegistrationTest Registration3D3DTypeTest; + typedef Registration3D3DTypeTest::Pointer Registration3D3DTypeTestPointer; + typedef map::core::Registration Registration3D3DType; + typedef Registration3D3DType::Pointer Registration3D3DPointer; + typedef rttb::interpolation::NearestNeighborInterpolation NearestNeighborInterpolation; + typedef rttb::interpolation::LinearInterpolation LinearInterpolation; + typedef rttb::interpolation::TransformationInterface TransformationInterface; + typedef rttb::interpolation::MatchPointTransformation MatchPointTransformation; + + /*! @brief SimpleMappableDoseAccessorWithMatchPointTest - test the API of SimpleMappableDoseAccessor with MatchPoint transform + 1) Test constructor + 2) test getDoseAt() + a) with Identity transform + b) with translation transform + [3) test with rigid registration optional (if filenames are given as argument)] + */ + + int SimpleMappableDoseAccessorWithMatchPointTest(int argc, char* argv[]) + { + PREPARE_DEFAULT_TEST_REPORTING; + + std::string RTDOSE_FILENAME_CONSTANT_TWO; + std::string RTDOSE_FILENAME_INCREASE_X; + std::string RTDOSE_FILENAME_REALISTIC = ""; + std::string CT_PLANNING = ""; + std::string CT_FRACTION = ""; + + if (argc > 2) + { + RTDOSE_FILENAME_CONSTANT_TWO = argv[1]; + RTDOSE_FILENAME_INCREASE_X = argv[2]; + } + else + { + std::cout << "at least two parameters for SimpleMappableDoseAccessorWithMatchPointTest are expected" + << std::endl; + return -1; + } + + if (argc > 5) + { + RTDOSE_FILENAME_REALISTIC = argv[3]; + CT_PLANNING = argv[4]; + CT_FRACTION = argv[5]; + } + + 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()); + + core::GeometricInfo doseAccessor1GeometricInfo = doseAccessor1->getGeometricInfo(); + + Registration3D3DTypeTestPointer registration = Registration3D3DTypeTest::New(); + double translation[] = {0.0, 0.0, 0.0}; + registration->_translation = translation; + registration->_limitedTarget = false; + + auto interpolationNN = + NearestNeighborInterpolation::Pointer(new NearestNeighborInterpolation()); + auto interpolationLinear = boost::shared_ptr + (new LinearInterpolation()); + boost::shared_ptr interpolationNull; + + auto transformMP = TransformationInterface::Pointer(new MatchPointTransformation( + registration.GetPointer())); + + auto aSimpleMappableDoseAccessorMPIdentityLinear = new SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformMP, interpolationLinear); + auto aSimpleMappableDoseAccessorMPIdentityNN = new SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformMP, interpolationNN); + + //1) Test constructor + CHECK_NO_THROW(SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformMP, interpolationLinear)); + CHECK_NO_THROW(SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformMP, interpolationNN)); + CHECK_NO_THROW(MatchPointTransformation(registration.GetPointer())); + CHECK_THROW_EXPLICIT(MatchPointTransformation(NULL), core::NullPointerException); + + //2) test getDoseAt() + // a) with Identity transform + + double vectorDoseAccessorStartEnd = 0.0; + + while (vectorDoseAccessorStartEnd <= 1.0) + { + VoxelGridID runningID = (VoxelGridID)(vectorDoseAccessorStartEnd * + (double)aSimpleMappableDoseAccessorMPIdentityLinear->getGridSize()); + + CHECK_EQUAL(aSimpleMappableDoseAccessorMPIdentityLinear->getDoseAt(runningID), + doseAccessor2->getDoseAt(runningID)); + CHECK_EQUAL(aSimpleMappableDoseAccessorMPIdentityNN->getDoseAt(runningID), + doseAccessor2->getDoseAt(runningID)); + vectorDoseAccessorStartEnd += 0.1; + } + + + // b) with translation transform + + + //Second: Translation (5mm/5mm/5mm) --> in voxel: (1/1/1) as pixelspacing = 5 mm + translation[0] = translation[1] = translation[2] = 5.0; + registration->_translation = translation; + auto aSimpleMappableDoseAccessorMPTranslationLinear = new SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformMP, interpolationLinear); + auto aSimpleMappableDoseAccessorMPTranslationNN = new SimpleMappableDoseAccessor( + doseAccessor1GeometricInfo, doseAccessor2, transformMP, interpolationNN); + + rttb::VoxelGridIndex3D aIndexBeforeTransformation(0, 0, 0); + rttb::VoxelGridIndex3D aIndexAfterTransformation(1, 1, 1); + CHECK_EQUAL(aSimpleMappableDoseAccessorMPTranslationLinear->getDoseAt(aIndexBeforeTransformation), + doseAccessor2->getDoseAt(aIndexAfterTransformation)); + CHECK_EQUAL(aSimpleMappableDoseAccessorMPTranslationNN->getDoseAt(aIndexBeforeTransformation), + doseAccessor2->getDoseAt(aIndexAfterTransformation)); + + rttb::VoxelGridIndex3D aIndexBeforeTransformation2(20, 10, 10); + rttb::VoxelGridIndex3D aIndexAfterTransformation2(21, 11, 11); + CHECK_EQUAL(aSimpleMappableDoseAccessorMPTranslationLinear->getDoseAt(aIndexBeforeTransformation2), + doseAccessor2->getDoseAt(aIndexAfterTransformation2)); + CHECK_EQUAL(aSimpleMappableDoseAccessorMPTranslationNN->getDoseAt(aIndexBeforeTransformation2), + doseAccessor2->getDoseAt(aIndexAfterTransformation2)); + + rttb::VoxelGridIndex3D aIndexBeforeTransformation3( + aSimpleMappableDoseAccessorMPTranslationLinear->getGeometricInfo().getNumColumns() - 2, + aSimpleMappableDoseAccessorMPTranslationLinear->getGeometricInfo().getNumRows() - 2, + aSimpleMappableDoseAccessorMPTranslationLinear->getGeometricInfo().getNumSlices() - 2); + rttb::VoxelGridIndex3D aIndexAfterTransformation3( + aSimpleMappableDoseAccessorMPTranslationLinear->getGeometricInfo().getNumColumns() - 1, + aSimpleMappableDoseAccessorMPTranslationLinear->getGeometricInfo().getNumRows() - 1, + aSimpleMappableDoseAccessorMPTranslationLinear->getGeometricInfo().getNumSlices() - 1); + CHECK_EQUAL(aSimpleMappableDoseAccessorMPTranslationLinear->getDoseAt(aIndexBeforeTransformation3), + doseAccessor2->getDoseAt(aIndexAfterTransformation3)); + CHECK_EQUAL(aSimpleMappableDoseAccessorMPTranslationNN->getDoseAt(aIndexBeforeTransformation3), + doseAccessor2->getDoseAt(aIndexAfterTransformation3)); + + if (RTDOSE_FILENAME_REALISTIC != "" && CT_FRACTION != "" && CT_PLANNING != "") + { + //3) test with rigid registration + //realistic background: registration from BP-CT to fraction CT, apply on planning dose that is based on BP-CT (proof of concept) + + //Target image: fraction CT, Moving image: planning CT + simpleRegistrationWorkflow prepareRegistrationRealisticScenario(CT_FRACTION, CT_PLANNING, true); + Registration3D3DPointer registrationRealisticScenario = + prepareRegistrationRealisticScenario.getRegistration(); + + auto transformRealistic = TransformationInterface::Pointer(new MatchPointTransformation( + registrationRealisticScenario)); + + io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator3(RTDOSE_FILENAME_REALISTIC.c_str()); + DoseAccessorPointer doseAccessor3(doseAccessorGenerator3.generateDoseAccessor()); + + core::GeometricInfo geoInfoRealistic; + geoInfoRealistic.setNumColumns( + prepareRegistrationRealisticScenario.getTargetImage()->GetLargestPossibleRegion().GetSize()[0]); + geoInfoRealistic.setNumRows( + prepareRegistrationRealisticScenario.getTargetImage()->GetLargestPossibleRegion().GetSize()[1]); + geoInfoRealistic.setNumSlices( + prepareRegistrationRealisticScenario.getTargetImage()->GetLargestPossibleRegion().GetSize()[2]); + + //Dose is related to BP-CT, map dose to fraction CT geometry + auto aSimpleMappableDoseAccessorRealisticScenarioLinear = + new SimpleMappableDoseAccessor(geoInfoRealistic, + doseAccessor3, transformRealistic, interpolationLinear); + + //combination of 0, size/2 and size to check as coordinates + std::vector coordinatesToCheckX, coordinatesToCheckY, coordinatesToCheckZ; + coordinatesToCheckX.push_back(0); + coordinatesToCheckX.push_back(geoInfoRealistic.getNumColumns() / 2); + coordinatesToCheckX.push_back(geoInfoRealistic.getNumColumns() - 1); + coordinatesToCheckY.push_back(0); + coordinatesToCheckY.push_back(geoInfoRealistic.getNumRows() / 2); + coordinatesToCheckY.push_back(geoInfoRealistic.getNumRows() - 1); + coordinatesToCheckZ.push_back(0); + coordinatesToCheckZ.push_back(geoInfoRealistic.getNumSlices() / 2); + coordinatesToCheckZ.push_back(geoInfoRealistic.getNumSlices() - 1); + + //Pixels are inside the fraction CT image and mapping should work (even if they map outside of doseAccessor3) + for (unsigned int i = 0; i < coordinatesToCheckX.size(); ++i) + { + for (unsigned int j = 0; j < coordinatesToCheckY.size(); ++j) + { + for (unsigned int k = 0; k < coordinatesToCheckZ.size(); ++k) + { + CHECK_NO_THROW(aSimpleMappableDoseAccessorRealisticScenarioLinear->getDoseAt(VoxelGridIndex3D( + coordinatesToCheckX.at(i), coordinatesToCheckY.at(j), coordinatesToCheckZ.at(k)))); + } + } + } + + } + + RETURN_AND_REPORT_TEST_SUCCESS; + } + + }//end namespace testing +}//end namespace rttb diff --git a/testing/interpolation/MatchPointBinding/files.cmake b/testing/interpolation/MatchPointBinding/files.cmake new file mode 100644 index 0000000..dba2d96 --- /dev/null +++ b/testing/interpolation/MatchPointBinding/files.cmake @@ -0,0 +1,12 @@ +SET(CPP_FILES + SimpleMappableDoseAccessorWithMatchPointTest.cpp + rttbMatchPointBindingTests.cpp + registrationHelper.cpp + simpleRegistrationWorkflow.cpp + ) + +SET(H_FILES + simpleRegistrationWorkflow.h + registrationTest.h + registrationHelper.h + ) diff --git a/testing/interpolation/MatchPointBinding/registrationHelper.cpp b/testing/interpolation/MatchPointBinding/registrationHelper.cpp new file mode 100644 index 0000000..4aa1dd9 --- /dev/null +++ b/testing/interpolation/MatchPointBinding/registrationHelper.cpp @@ -0,0 +1,88 @@ +// ----------------------------------------------------------------------- +// 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: 741 $ (last changed revision) +// @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ + +#if defined(_MSC_VER) +#pragma warning ( disable : 4786 ) +#endif + +#include "registrationHelper.h" + +#include "litTestImageIO.h" +#include "litCheckMacros.h" +#include "litImageTester.h" +#include "litPointSetTester.h" + +#include + +int setImageFileNames(std::string targetImage, std::string movingImage, bool isDirectory, + AppGlobals& globals) +{ + globals.targetImageFileName = targetImage; + globals.movingImageFileName = movingImage; + globals.isDirectory = isDirectory; + + return EXIT_SUCCESS; +} + +int loadData(AppGlobals& globals) +{ + if (!globals.isDirectory) + { + globals.spTargetImage = + lit::TestImageIO::InternalImageType>::readImage( + globals.targetImageFileName); + } + else + { + globals.spTargetImage = map::io::readImage + (globals.targetImageFileName, map::io::ImageSeriesReadStyle::Dicom); + } + + if (globals.spTargetImage.IsNull()) + { + std::cerr << "Error. Cannot load target image: " << globals.targetImageFileName << std::endl; + return EXIT_FAILURE; + } + + if (!globals.isDirectory) + { + globals.spMovingImage = + lit::TestImageIO::InternalImageType>::readImage( + globals.movingImageFileName); + } + else + { + globals.spMovingImage = map::io::readImage + (globals.movingImageFileName, map::io::ImageSeriesReadStyle::Dicom); + } + + if (globals.spMovingImage.IsNull()) + { + std::cerr << "Error. Cannot load moving image: " << globals.movingImageFileName << std::endl; + return EXIT_FAILURE; + } + + return EXIT_SUCCESS; +} + +AppGlobals::AppGlobals() +{ +}; \ No newline at end of file diff --git a/testing/interpolation/mapDemoHelloWorldRegistration1Helper.h b/testing/interpolation/MatchPointBinding/registrationHelper.h similarity index 58% rename from testing/interpolation/mapDemoHelloWorldRegistration1Helper.h rename to testing/interpolation/MatchPointBinding/registrationHelper.h index 8444a4d..dc788b3 100644 --- a/testing/interpolation/mapDemoHelloWorldRegistration1Helper.h +++ b/testing/interpolation/MatchPointBinding/registrationHelper.h @@ -1,59 +1,56 @@ // ----------------------------------------------------------------------- -// MatchPoint - DKFZ translational registration framework +// 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 mapCopyright.txt or -// http://www.dkfz.de/en/sidt/projects/MatchPoint/copyright.html +// 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. +// 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) -// Subversion HeadURL: $HeadURL: https://svn/sbr/Sources/SBR-Projects/MatchPoint/trunk/Examples/Algorithms/mapDemoHelloWorldRegistration1Helper.h $ +// @version $Revision: 741 $ (last changed revision) +// @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) */ #if defined(_MSC_VER) #pragma warning ( disable : 4786 ) #endif -#ifndef __MAP_DEMO_HELLO_WORLD_REGISTRATION1_HELPER_H -#define __MAP_DEMO_HELLO_WORLD_REGISTRATION1_HELPER_H +#ifndef __REGISTRATION_HELPER_H +#define __REGISTRATION_HELPER_H #include "mapContinuousElements.h" #include "mapDiscreteElements.h" #include "mapImageReader.h" typedef map::core::discrete::Elements<3>::InternalImageType ImageType; typedef map::core::continuous::Elements<3>::InternalPointSetType LandmarksType; struct AppGlobals { std::string targetImageFileName; std::string movingImageFileName; - std::string resultImageFileName; bool isDirectory; ImageType::Pointer spTargetImage; ImageType::Pointer spMovingImage; ImageType::Pointer spResultImage; AppGlobals(); }; -int setImageFileNames(std::string targetImage, std::string movingImage, std::string resultImage, bool isDirectory, AppGlobals& globals); +int setImageFileNames(std::string targetImage, std::string movingImage, bool isDirectory, + AppGlobals& globals); int loadData(AppGlobals& globals); -int saveResults(AppGlobals& globals); - #endif \ No newline at end of file diff --git a/testing/interpolation/registrationTest.h b/testing/interpolation/MatchPointBinding/registrationTest.h similarity index 50% rename from testing/interpolation/registrationTest.h rename to testing/interpolation/MatchPointBinding/registrationTest.h index d551e99..ab3a4a9 100644 --- a/testing/interpolation/registrationTest.h +++ b/testing/interpolation/MatchPointBinding/registrationTest.h @@ -1,59 +1,79 @@ -#ifndef __TESTREGISTRATION_H -#define __TESTREGISTRATION_H +// ----------------------------------------------------------------------- +// 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: 741 $ (last changed revision) +// @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ +#ifndef __REGISTRATION_TEST_H +#define __REGISTRATION_TEST_H #include "mapRegistration.h" namespace map { namespace core { - /*! @class TestRegistration + /*! @class TestRegistration @brief Simple implementation of MatchPoint Registration class with direct access to mapping. */ template class RegistrationTest: public Registration { public: typedef RegistrationTest Self; typedef RegistrationBase Superclass; typedef itk::SmartPointer Pointer; typedef itk::SmartPointer ConstPointer; itkTypeMacro(RegistrationTest, Registration); itkNewMacro(Self); bool _limitedTarget; - double* translation; + double* _translation; RegistrationTest() {}; ~RegistrationTest() { }; virtual bool mapPointInverse(const TargetPointType& inPoint, MovingPointType& outPoint) const { for (unsigned int i = 0; i < VTargetDimensions; i++) { - outPoint[i] = inPoint[i] + translation[i]; + outPoint[i] = inPoint[i] + _translation[i]; } return true; }; virtual bool hasLimitedTargetRepresentation() const { return _limitedTarget; } private: RegistrationTest(const Self& source); //purposely not implemented void operator=(const Self&); //purposely not implemented }; } } #endif \ No newline at end of file diff --git a/testing/interpolation/MatchPointBinding/rttbMatchPointBindingTests.cpp b/testing/interpolation/MatchPointBinding/rttbMatchPointBindingTests.cpp new file mode 100644 index 0000000..473ee84 --- /dev/null +++ b/testing/interpolation/MatchPointBinding/rttbMatchPointBindingTests.cpp @@ -0,0 +1,63 @@ +// ----------------------------------------------------------------------- +// 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(SimpleMappableDoseAccessorWithMatchPointTest); + } + } +} + +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; +} diff --git a/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.cpp b/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.cpp new file mode 100644 index 0000000..9ea6c4a --- /dev/null +++ b/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.cpp @@ -0,0 +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: 741 $ (last changed revision) +// @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ + +#undef MAP_SEAL_ALGORITHMS + +#include "simpleRegistrationWorkflow.h" +#include "registrationHelper.h" + +simpleRegistrationWorkflow::simpleRegistrationWorkflow(std::string targetFilename, + std::string movingFilename, + bool isDirectory) +{ + _targetFilename = targetFilename; + _movingFilename = movingFilename; + + setImageFileNames(_targetFilename, _movingFilename, isDirectory, globals); + loadData(globals); + _spAlgorithmEuler = NULL; +} + +vnl_vector simpleRegistrationWorkflow::getRegistrationParameters( + Registration3D3DPointer reg) +{ + typedef map::core::ModelBasedRegistrationKernel<3, 3> ModelBasedRegistrationKernel3D3D; + + const ModelBasedRegistrationKernel3D3D* pModelBasedDirectKernel3D3D = + dynamic_cast(&(reg->getDirectMapping())); + + if (pModelBasedDirectKernel3D3D) + { + ModelBasedRegistrationKernel3D3D::ParametersType params = + pModelBasedDirectKernel3D3D->getTransformModel()->getTransform()->GetParameters(); + + return params; + } + else + { + return vnl_vector(); + } +} + +void simpleRegistrationWorkflow::initializeAndPerformRegistration() +{ + _spAlgorithmEuler = AlgorithmTypeEuler::New(); + _spAlgorithmEuler->setProperty("PreinitTransform", map::core::MetaProperty::New(true)); + _spAlgorithmEuler->setMovingImage(globals.spMovingImage); + _spAlgorithmEuler->setTargetImage(globals.spTargetImage); + AlgorithmTypeEuler::RegistrationType::Pointer spRegistration; + + try + { + spRegistration = _spAlgorithmEuler->getRegistration(); + } + catch (const map::core::ExceptionObject& e) + { + std::cerr << "caught an MatchPoint exception:\n"; + e.Print(std::cerr); + std::cerr << "\n"; + } + catch (const itk::ExceptionObject& e) + { + std::cerr << "caught an ITK exception:\n"; + std::cerr << e.GetFile() << ":" << e.GetLine() << ":\n" + << e.GetDescription() << "\n"; + } + catch (const std::exception& e) + { + std::cerr << "caught an exception:\n"; + std::cerr << e.what() << "\n"; + } + catch (...) + { + std::cerr << "caught an unknown exception!!!\n"; + } +} + +map::core::Registration<3, 3>::Pointer simpleRegistrationWorkflow::getRegistration() +{ + if (_spAlgorithmEuler.IsNull()) + { + initializeAndPerformRegistration(); + } + + return _spAlgorithmEuler->getRegistration(); +}; + +const itk::Image* simpleRegistrationWorkflow::getTargetImage() +{ + if (_spAlgorithmEuler.IsNull()) + { + initializeAndPerformRegistration(); + } + + return _spAlgorithmEuler->getTargetImage(); +}; diff --git a/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.h b/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.h new file mode 100644 index 0000000..73db052 --- /dev/null +++ b/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.h @@ -0,0 +1,64 @@ +// ----------------------------------------------------------------------- +// 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: 741 $ (last changed revision) +// @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ +#ifndef __SIMPLE_REGISTRATION_WORKFLOW_H +#define __SIMPLE_REGISTRATION_WORKFLOW_H + +#include +#include + +#include "registrationHelper.h" + +#include "mapImageMappingTask.h" +#include "mapITKEuler3DMattesMIRegistrationAlgorithmTemplate.h" +#include "mapExceptionObject.h" + +/*! @class simpleRegistrationWorkflow + @brief implements a concrete registration algorithm of MatchPoint +*/ +class simpleRegistrationWorkflow +{ +public: + typedef map::core::Registration<3, 3> Registration3D3DType; + typedef Registration3D3DType::Pointer Registration3D3DPointer; + typedef map::algorithm::boxed::ITKEuler3DMattesMIRegistrationAlgorithm + AlgorithmTypeEuler; +private: + std::string _targetFilename; + std::string _movingFilename; + std::string _targetDirectory; + std::string _movingDirectory; + AppGlobals globals; + AlgorithmTypeEuler::Pointer _spAlgorithmEuler; + +public: + /*! @brief Constructor + */ + simpleRegistrationWorkflow(std::string targetFilename, std::string movingFilename, + bool isDirectory = false); + map::core::Registration<3, 3>::Pointer getRegistration(); + const itk::Image* getTargetImage(); + vnl_vector getRegistrationParameters(Registration3D3DPointer reg); + +protected: + void initializeAndPerformRegistration(); +}; + +#endif \ No newline at end of file diff --git a/testing/interpolation/RosuMappableDoseAccessorTest.cpp b/testing/interpolation/RosuMappableDoseAccessorTest.cpp index 79a32eb..fd9dc90 100644 --- a/testing/interpolation/RosuMappableDoseAccessorTest.cpp +++ b/testing/interpolation/RosuMappableDoseAccessorTest.cpp @@ -1,162 +1,153 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbGenericDoseIterator.h" -#include "rttbDoseIteratorInterface.h" #include "rttbDicomFileDoseAccessorGenerator.h" -#include "rttbMatchPointTransformation.h" #include "rttbRosuMappableDoseAccessor.h" -#include "rttbNearestNeighborInterpolation.h" -#include "rttbLinearInterpolation.h" -#include "rttbITKImageDoseAccessor.h" + +#include "DummyTransformation.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" -#include "itkImage.h" - -#include "registrationTest.h" -#include "registrationAlgorithm.h" - namespace rttb { namespace testing { - static const unsigned int TargetDimension3D = 3; - static const unsigned int MovingDimension3D = 3; typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; - typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; typedef rttb::interpolation::RosuMappableDoseAccessor RosuMappableDoseAccessor; - typedef map::core::RegistrationTest Registration3D3DTypeTest; - typedef Registration3D3DTypeTest::Pointer Registration3D3DTypeTestPointer; - typedef map::core::Registration Registration3D3DType; - typedef Registration3D3DType::Pointer Registration3D3DPointer; typedef rttb::interpolation::TransformationInterface TransformationInterface; - typedef rttb::interpolation::MatchPointTransformation MatchPointTransformation; - typedef RosuMappableDoseAccessor::Pointer RosuMappableDoseAccessorPointer; - typedef rttb::interpolation::MappableDoseAccessorBase MappableDoseAccessorBase; /*! @brief RosuMappableDoseAccessorTest - test the API of RosuMappableDoseAccessor - 1) test constructor - 2) test with registration stub + 1) Constructor + 2) test getDoseAt() */ int RosuMappableDoseAccessorTest(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 MappableDoseAccessorTest are expected" << std::endl; + std::cout << "at least two parameters for RosuMappableDoseAccessorTest are expected" << std::endl; return -1; } + const double doseGridScaling = 2.822386e-005; + 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()); - const double doseGridScalingDose2 = 2.822386e-005; - - Registration3D3DTypeTestPointer registration = Registration3D3DTypeTest::New(); - double translation[] = {0.0, 0.0, 0.0}; - registration->translation = translation; - registration->_limitedTarget = false; - - TransformationInterface::Pointer transformMP = TransformationInterface::Pointer( - new MatchPointTransformation(registration)); - boost::shared_ptr transformNull; + DoseAccessorPointer doseAccessorNull; - MappableDoseAccessorBase::Pointer aMappableDoseAccessorDefault = - MappableDoseAccessorBase::Pointer(new RosuMappableDoseAccessor( - doseAccessor1->getGeometricInfo(), doseAccessor2, transformMP)); + auto transformDummy = TransformationInterface::Pointer(new DummyTransformation()); + TransformationInterface::Pointer transformNull; - //1) test constructors - CHECK_NO_THROW(RosuMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, - transformMP)); + auto aRosuMappableDoseAccessorDefault = + RosuMappableDoseAccessor::Pointer(new RosuMappableDoseAccessor( + doseAccessor1->getGeometricInfo(), doseAccessor2, transformDummy)); + auto aRosuMappableDoseAccessorNoPadding = + RosuMappableDoseAccessor::Pointer(new RosuMappableDoseAccessor( + doseAccessor1->getGeometricInfo(), doseAccessor2, transformDummy, false)); - //2) test with registration stub - //First: Identity Transformation (nothing changes). Should lead to the same value although the subvoxel values are different (but symmetric with respect to the voxel center) + //1) Constructor - //outside: should have default value 0 - CHECK_EQUAL(aMappableDoseAccessorDefault->getDoseAt( - aMappableDoseAccessorDefault->getGridSize() + 1), 0.0); + CHECK_NO_THROW(RosuMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, + transformDummy)); + CHECK_NO_THROW(RosuMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, + transformDummy, false)); + CHECK_NO_THROW(RosuMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, + transformDummy, true, 5.0)); - const core::GeometricInfo geoInfoTarget = doseAccessor1->getGeometricInfo(); + CHECK_THROW_EXPLICIT(RosuMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), + doseAccessorNull, transformDummy), core::NullPointerException); + CHECK_THROW_EXPLICIT(RosuMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), + doseAccessor2, transformNull), core::NullPointerException); + CHECK_THROW_EXPLICIT(RosuMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), + doseAccessorNull, transformNull), core::NullPointerException); - const VoxelGridIndex3D indexValid(geoInfoTarget.getNumColumns() - 1, geoInfoTarget.getNumRows() - 1, - geoInfoTarget.getNumSlices() - 1); + //2) test getDoseAt() - const VoxelGridID idValid(indexValid[2] * geoInfoTarget.getNumColumns() * geoInfoTarget.getNumRows() - + indexValid[1] * geoInfoTarget.getNumColumns() - + indexValid[0]); + std::vector voxelsAsIndexToTest; + std::vector voxelsAsIdToTest; + std::vector expectedValues; - //check if getDoseAt(VoxelGridID) and getDoseAt(VoxelGridIndex3D) lead to same result - CHECK_NO_THROW(aMappableDoseAccessorDefault->getDoseAt(indexValid)); - CHECK_EQUAL(aMappableDoseAccessorDefault->getDoseAt(indexValid), - aMappableDoseAccessorDefault->getDoseAt(idValid)); - CHECK_CLOSE(aMappableDoseAccessorDefault->getDoseAt(indexValid), - doseGridScalingDose2 * indexValid.x(), - errorConstant); - //dose value equals x-position. 0.5*(0.25+0.75)=0.5 (compute interpolation value at x-pos 0.25 and 0.75) - CHECK_CLOSE(aMappableDoseAccessorDefault->getDoseAt(VoxelGridID(0)), 0.5 * doseGridScalingDose2, - errorConstant); + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(5, 9, 8)); + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(22, 15, 10)); + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(30, 20, 7)); - //Second: Translation (2.5mm/2.5mm/2.5mm) --> in voxel: (0.5/0.5/0.5) as pixelspacing = 5 mm - translation[0] = translation[1] = translation[2] = 2.5; - registration->translation = translation; - rttb::VoxelGridIndex3D aIndex(23, 11, 9); + expectedValues.push_back(5.0 * doseGridScaling); + expectedValues.push_back(22.0 * doseGridScaling); + expectedValues.push_back(30.0 * doseGridScaling); - //value equals x-position. 0.5*((indexValid.x()+0.25+translation_x)+(indexValid.x()+0.75+translation_x))=1 - CHECK_CLOSE(aMappableDoseAccessorDefault->getDoseAt(aIndex), - (aIndex.x() + 1.0)*doseGridScalingDose2, errorConstant); + //convert VoxelGridIndex3D to VoxelGridID + for (int i = 0; i < voxelsAsIndexToTest.size(); i++) + { + VoxelGridID currentId; + doseAccessor1->getGeometricInfo().convert(voxelsAsIndexToTest.at(i), currentId); + voxelsAsIdToTest.push_back(currentId); + } - rttb::VoxelGridIndex3D aIndex2(49, 32, 29); + for (int i = 0; i < voxelsAsIndexToTest.size(); i++) + { + //test if the expected interpolation values are computed + CHECK_CLOSE(aRosuMappableDoseAccessorDefault->getDoseAt(voxelsAsIndexToTest.at(i)), + expectedValues.at(i), errorConstant); + //test if getDoseAt(VoxelGridIndex3D) and getDoseAt(VoxelGridD) lead to the same results + CHECK_EQUAL(aRosuMappableDoseAccessorDefault->getDoseAt(voxelsAsIndexToTest.at(i)), + aRosuMappableDoseAccessorDefault->getDoseAt(voxelsAsIdToTest.at(i))); + } - CHECK_CLOSE(aMappableDoseAccessorDefault->getDoseAt(aIndex2), - (aIndex2.x() + 1.0)*doseGridScalingDose2, errorConstant); + //test invalid voxels + VoxelGridID invalidID(doseAccessor1->getGeometricInfo().getNumberOfVoxels() + 1); + VoxelGridIndex3D invalidIndex(doseAccessor1->getGeometricInfo().getNumColumns() + 1, + doseAccessor1->getGeometricInfo().getNumRows() + 1, + doseAccessor1->getGeometricInfo().getNumSlices() + 1); + CHECK_EQUAL(aRosuMappableDoseAccessorDefault->getDoseAt(invalidID), 0.0); + CHECK_EQUAL(aRosuMappableDoseAccessorDefault->getDoseAt(invalidIndex), 0.0); + CHECK_THROW_EXPLICIT(aRosuMappableDoseAccessorNoPadding->getDoseAt(invalidID), + core::MappingOutsideOfImageException); + CHECK_THROW_EXPLICIT(aRosuMappableDoseAccessorNoPadding->getDoseAt(invalidIndex), + core::MappingOutsideOfImageException); - //@todo: more tests needed? RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/interpolation/SimpleMappableDoseAccessorTest.cpp b/testing/interpolation/SimpleMappableDoseAccessorTest.cpp index 81ed6e7..84ad44d 100644 --- a/testing/interpolation/SimpleMappableDoseAccessorTest.cpp +++ b/testing/interpolation/SimpleMappableDoseAccessorTest.cpp @@ -1,341 +1,165 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbGenericDoseIterator.h" #include "rttbDicomFileDoseAccessorGenerator.h" -#include "rttbDicomHelaxDoseAccessor.h" -#include "rttbDicomHelaxFileDoseAccessorGenerator.h" #include "rttbSimpleMappableDoseAccessor.h" #include "rttbNearestNeighborInterpolation.h" #include "rttbLinearInterpolation.h" -#include "rttbITKImageDoseAccessor.h" #include "rttbTransformationInterface.h" -#include "rttbMatchPointTransformation.h" +#include "DummyTransformation.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" -#include "itkImage.h" - -#include "registrationTest.h" -#include "registrationAlgorithm.h" - namespace rttb { namespace testing { - static const unsigned int TargetDimension3D = 3; - static const unsigned int MovingDimension3D = 3; typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef rttb::interpolation::SimpleMappableDoseAccessor SimpleMappableDoseAccessor; - typedef map::core::RegistrationTest Registration3D3DTypeTest; - typedef Registration3D3DTypeTest::Pointer Registration3D3DTypeTestPointer; - typedef map::core::Registration Registration3D3DType; - typedef Registration3D3DType::Pointer Registration3D3DPointer; - typedef rttb::interpolation::NearestNeighborInterpolation NearestNeighborInterpolation; - typedef rttb::interpolation::LinearInterpolation LinearInterpolation; typedef rttb::interpolation::TransformationInterface TransformationInterface; - typedef rttb::interpolation::MatchPointTransformation MatchPointTransformation; - typedef rttb::interpolation::MappableDoseAccessorBase MappableDoseAccessorBase; + typedef rttb::interpolation::LinearInterpolation LinearInterpolation; + typedef rttb::interpolation::NearestNeighborInterpolation NearestNeighborInterpolation; - /*! @brief SimpleMappableDoseAccessorTest - test the API of SimpleMappableDoseAccessor + MappableDoseAccessorBase - 1) test constructor - 2) test getGeometricInfo(), validId(), validIndex() function - 3) test with registration stub - optional (depends on input parameters): - 4) test with rigid registration + /*! @brief SimpleMappableDoseAccessorTest - test the API of SimpleMappableDoseAccessor + 1) Test constructor + 2) test getDoseAt() */ int SimpleMappableDoseAccessorTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string RTDOSE_FILENAME_CONSTANT_TWO; std::string RTDOSE_FILENAME_INCREASE_X; - std::string RTDOSE_FILENAME_REALISTIC; - std::string CT_PLANNING; - std::string CT_FRACTION; if (argc > 2) { RTDOSE_FILENAME_CONSTANT_TWO = argv[1]; RTDOSE_FILENAME_INCREASE_X = argv[2]; - RTDOSE_FILENAME_REALISTIC = ""; - CT_PLANNING = ""; - CT_FRACTION = ""; } else { - std::cout << "at least two parameters for MappableDoseAccessorTest are expected" << std::endl; + std::cout << "at least two parameters for SimpleMappableDoseAccessorTest are expected" << std::endl; return -1; } - if (argc > 5) - { - RTDOSE_FILENAME_CONSTANT_TWO = argv[1]; - RTDOSE_FILENAME_INCREASE_X = argv[2]; - RTDOSE_FILENAME_REALISTIC = argv[3]; - CT_PLANNING = argv[4]; - CT_FRACTION = argv[5]; - } + const double doseGridScaling = 2.822386e-005; rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1( RTDOSE_FILENAME_CONSTANT_TWO.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); - - DoseAccessorPointer doseAccessorNull; + core::GeometricInfo doseAccessor1GeometricInfo = doseAccessor1->getGeometricInfo(); rttb::io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator2( RTDOSE_FILENAME_INCREASE_X.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); - Registration3D3DTypeTestPointer registration = Registration3D3DTypeTest::New(); - double translation[] = {0.0, 0.0, 0.0}; - registration->translation = translation; - registration->_limitedTarget = false; + DoseAccessorPointer doseAccessorNull; + + auto transformDummy = TransformationInterface::Pointer(new DummyTransformation()); + TransformationInterface::Pointer transformNull; auto interpolationNN = NearestNeighborInterpolation::Pointer(new NearestNeighborInterpolation()); auto interpolationLinear = boost::shared_ptr (new LinearInterpolation()); boost::shared_ptr interpolationNull; - auto transformMP = TransformationInterface::Pointer(new MatchPointTransformation(registration)); - boost::shared_ptr transformNull; - - auto aSimpleMappableDoseAccessorMPDefaultNN = MappableDoseAccessorBase::Pointer( - new SimpleMappableDoseAccessor(doseAccessor1->getGeometricInfo(), doseAccessor2, transformMP, - interpolationNN)); - auto aSimpleMappableDoseAccessorMPNoPaddingNN = MappableDoseAccessorBase::Pointer( - new SimpleMappableDoseAccessor(doseAccessor1->getGeometricInfo(), doseAccessor2, transformMP, - interpolationNN, false)); - auto aSimpleMappableDoseAccessorMPPaddingNN = MappableDoseAccessorBase::Pointer( - new SimpleMappableDoseAccessor(doseAccessor1->getGeometricInfo(), doseAccessor2, transformMP, - interpolationNN, true, 42.0)); - auto aSimpleMappableDoseAccessorMPDefaultLinear = MappableDoseAccessorBase::Pointer( - new SimpleMappableDoseAccessor(doseAccessor1->getGeometricInfo(), doseAccessor2, transformMP, - interpolationLinear)); - auto aSimpleMappableDoseAccessorMPNoPaddingLinear = MappableDoseAccessorBase::Pointer( - new SimpleMappableDoseAccessor(doseAccessor1->getGeometricInfo(), doseAccessor2, transformMP, - interpolationLinear, false)); - auto aSimpleMappableDoseAccessorMPPaddingLinear = MappableDoseAccessorBase::Pointer( - new SimpleMappableDoseAccessor(doseAccessor1->getGeometricInfo(), doseAccessor2, transformMP, - interpolationLinear, true, 42.0)); + auto aSimpleMappableDoseAccessorDefault = new SimpleMappableDoseAccessor( + doseAccessor1->getGeometricInfo(), doseAccessor2, transformDummy, interpolationLinear); + auto aSimpleMappableDoseAccessorNoPadding = new SimpleMappableDoseAccessor( + doseAccessor1->getGeometricInfo(), doseAccessor2, transformDummy, interpolationLinear, false); - const core::GeometricInfo geoInfoMoving = doseAccessor1->getGeometricInfo(); - core::GeometricInfo geoInfoTarget = geoInfoMoving; - geoInfoTarget.setNumColumns(geoInfoTarget.getNumColumns() / 2); - geoInfoTarget.setNumRows(geoInfoTarget.getNumRows() / 2); - geoInfoTarget.setNumSlices(geoInfoTarget.getNumSlices() / 2); - geoInfoTarget.setPixelSpacingRow(geoInfoTarget.getPixelSpacingRow() * 2.0); - geoInfoTarget.setPixelSpacingColumn(geoInfoTarget.getPixelSpacingColumn() * 2.0); - geoInfoTarget.setSliceThickness(geoInfoTarget.getSliceThickness() * 2.0); + //1) Test constructor - const VoxelGridIndex3D indexValid(geoInfoTarget.getNumColumns() - 1, geoInfoTarget.getNumRows() - 1, - geoInfoTarget.getNumSlices() - 1); - const VoxelGridIndex3D indexInvalid(geoInfoMoving.getNumColumns() - 1, - geoInfoMoving.getNumRows() - 1, geoInfoMoving.getNumSlices() - 1); - - const VoxelGridID idValid(indexValid[2] * geoInfoTarget.getNumColumns() * geoInfoTarget.getNumRows() - + - indexValid[1] * geoInfoTarget.getNumColumns() - + indexValid[0]); - const VoxelGridID idInvalid(indexInvalid[2] * geoInfoTarget.getNumColumns() * - geoInfoTarget.getNumRows() + - indexInvalid[1] * geoInfoTarget.getNumColumns() - + indexInvalid[0]); - - auto aMappableDoseAccessorMPChangedGeometricInfo = MappableDoseAccessorBase::Pointer( - new SimpleMappableDoseAccessor(geoInfoTarget, doseAccessor1, transformMP, - interpolationNN)); - - //1) test constructors - CHECK_NO_THROW(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, - transformMP, interpolationNN)); - CHECK_NO_THROW(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, - transformMP, interpolationNN, false)); CHECK_NO_THROW(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, - transformMP, interpolationNN, true, 5.0)); + transformDummy, interpolationLinear)); CHECK_NO_THROW(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, - transformMP, interpolationLinear)); + transformDummy, interpolationLinear, false)); CHECK_NO_THROW(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, - transformMP, interpolationLinear, false)); - CHECK_NO_THROW(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), doseAccessor2, - transformMP, interpolationLinear, true, 5.0)); + transformDummy, interpolationLinear, true, 5.0)); + CHECK_NO_THROW(SimpleMappableDoseAccessor(doseAccessor1GeometricInfo, doseAccessor2, + transformDummy, interpolationNN)); CHECK_THROW_EXPLICIT(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), - doseAccessorNull, transformMP, interpolationNN), core::NullPointerException); + doseAccessorNull, transformDummy, interpolationLinear), core::NullPointerException); CHECK_THROW_EXPLICIT(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), - doseAccessor2, transformNull, interpolationNN), core::NullPointerException); + doseAccessor2, transformNull, interpolationLinear), core::NullPointerException); CHECK_THROW_EXPLICIT(SimpleMappableDoseAccessor test(doseAccessor1->getGeometricInfo(), - doseAccessor2, - transformMP, interpolationNull), core::NullPointerException); - - //2) test getGeometricInfo(), validId(), validIndex() function - - CHECK_EQUAL(aMappableDoseAccessorMPChangedGeometricInfo->getGeometricInfo(), geoInfoTarget); - - - //3) test with registration stub - //First: Identity Transformation (nothing changes) - double vectorDoseAccessorStartEnd = 0.0; + doseAccessorNull, transformNull, interpolationLinear), core::NullPointerException); + CHECK_THROW_EXPLICIT(SimpleMappableDoseAccessor(doseAccessor1GeometricInfo, doseAccessor2, + transformDummy, interpolationNull), core::NullPointerException); - while (vectorDoseAccessorStartEnd <= 1.0) - { - VoxelGridID id = (VoxelGridID)(vectorDoseAccessorStartEnd * - (double)aSimpleMappableDoseAccessorMPDefaultNN->getGridSize()); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultNN->getDoseAt(id), doseAccessor2->getDoseAt(id)); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultLinear->getDoseAt(id), - doseAccessor2->getDoseAt(id)); - vectorDoseAccessorStartEnd += 0.1; - } + //2) test getGeometricInfo(), getGridSize(), getDoseUID() function - CHECK_THROW_EXPLICIT(aSimpleMappableDoseAccessorMPNoPaddingNN->getDoseAt( - aSimpleMappableDoseAccessorMPNoPaddingNN->getGridSize() + 1), core::MappingOutsideOfImageException); - CHECK_NO_THROW(aSimpleMappableDoseAccessorMPDefaultNN->getDoseAt( - aSimpleMappableDoseAccessorMPDefaultNN->getGridSize() + 1)); - CHECK_NO_THROW(aSimpleMappableDoseAccessorMPPaddingNN->getDoseAt( - aSimpleMappableDoseAccessorMPPaddingNN->getGridSize() + 1)); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultNN->getDoseAt( - aSimpleMappableDoseAccessorMPDefaultNN->getGridSize() + 1), 0.0); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPPaddingNN->getDoseAt( - aSimpleMappableDoseAccessorMPPaddingNN->getGridSize() + 1), 42.0); + CHECK_EQUAL(aSimpleMappableDoseAccessorDefault->getGeometricInfo(), + doseAccessor1->getGeometricInfo()); + CHECK_EQUAL(aSimpleMappableDoseAccessorDefault->getGridSize(), + doseAccessor1->getGeometricInfo().getNumberOfVoxels()); + CHECK_EQUAL(aSimpleMappableDoseAccessorDefault->getDoseUID(), doseAccessor2->getDoseUID()); - CHECK_THROW_EXPLICIT(aSimpleMappableDoseAccessorMPNoPaddingLinear->getDoseAt( - aSimpleMappableDoseAccessorMPNoPaddingLinear->getGridSize() + 1), - core::MappingOutsideOfImageException); - CHECK_NO_THROW(aSimpleMappableDoseAccessorMPDefaultLinear->getDoseAt( - aSimpleMappableDoseAccessorMPDefaultLinear->getGridSize() + 1)); - CHECK_NO_THROW(aSimpleMappableDoseAccessorMPPaddingLinear->getDoseAt( - aSimpleMappableDoseAccessorMPPaddingLinear->getGridSize() + 1)); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultLinear->getDoseAt( - aSimpleMappableDoseAccessorMPDefaultLinear->getGridSize() + 1), 0.0); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPPaddingLinear->getDoseAt( - aSimpleMappableDoseAccessorMPPaddingLinear->getGridSize() + 1), 42.0); + //3) test getDoseAt() - //check if getDoseAt(VoxelGridID) and getDoseAt(VoxelGridIndex3D) lead to same result - CHECK_NO_THROW(aMappableDoseAccessorMPChangedGeometricInfo->getDoseAt(indexValid)); - CHECK_NO_THROW(aMappableDoseAccessorMPChangedGeometricInfo->getDoseAt(indexInvalid)); - CHECK_EQUAL(aMappableDoseAccessorMPChangedGeometricInfo->getDoseAt(indexValid), - aMappableDoseAccessorMPChangedGeometricInfo->getDoseAt(idValid)); - CHECK_EQUAL(aMappableDoseAccessorMPChangedGeometricInfo->getDoseAt(indexInvalid), 0.0); + //test valid voxels + std::vector voxelsAsIndexToTest; + std::vector expectedValues; - //Second: Translation (5mm/5mm/5mm) --> in voxel: (1/1/1) as pixelspacing = 5 mm - translation[0] = translation[1] = translation[2] = 5.0; - registration->translation = translation; - rttb::VoxelGridIndex3D aIndexBeforeTransformation(0, 0, 0); - rttb::VoxelGridIndex3D aIndexAfterTransformation(1, 1, 1); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultNN->getDoseAt(aIndexBeforeTransformation), - doseAccessor2->getDoseAt(aIndexAfterTransformation)); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultLinear->getDoseAt(aIndexBeforeTransformation), - doseAccessor2->getDoseAt(aIndexAfterTransformation)); + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(5, 9, 8)); + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(22, 15, 10)); + voxelsAsIndexToTest.push_back(VoxelGridIndex3D(30, 20, 7)); - rttb::VoxelGridIndex3D aIndexBeforeTransformation2(20, 10, 10); - rttb::VoxelGridIndex3D aIndexAfterTransformation2(21, 11, 11); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultNN->getDoseAt(aIndexBeforeTransformation2), - doseAccessor2->getDoseAt(aIndexAfterTransformation2)); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultLinear->getDoseAt(aIndexBeforeTransformation2), - doseAccessor2->getDoseAt(aIndexAfterTransformation2)); + expectedValues.push_back(5.0 * doseGridScaling); + expectedValues.push_back(22.0 * doseGridScaling); + expectedValues.push_back(30.0 * doseGridScaling); - rttb::VoxelGridIndex3D aIndexBeforeTransformation3( - aSimpleMappableDoseAccessorMPDefaultNN->getGeometricInfo().getNumColumns() - 2, - aSimpleMappableDoseAccessorMPDefaultNN->getGeometricInfo().getNumRows() - 2, - aSimpleMappableDoseAccessorMPDefaultNN->getGeometricInfo().getNumSlices() - 2); - rttb::VoxelGridIndex3D aIndexAfterTransformation3( - aSimpleMappableDoseAccessorMPDefaultNN->getGeometricInfo().getNumColumns() - 1, - aSimpleMappableDoseAccessorMPDefaultNN->getGeometricInfo().getNumRows() - 1, - aSimpleMappableDoseAccessorMPDefaultNN->getGeometricInfo().getNumSlices() - 1); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultNN->getDoseAt(aIndexBeforeTransformation3), - doseAccessor2->getDoseAt(aIndexAfterTransformation3)); - CHECK_EQUAL(aSimpleMappableDoseAccessorMPDefaultLinear->getDoseAt(aIndexBeforeTransformation3), - doseAccessor2->getDoseAt(aIndexAfterTransformation3)); - - if (RTDOSE_FILENAME_REALISTIC != "" && CT_FRACTION != "" && CT_PLANNING != "") + for (int i = 0; i < voxelsAsIndexToTest.size(); i++) { - //4) test with rigid registration - //realistic background: registration from BP-CT to fraction CT, apply on planning dose that is based on BP-CT (proof of concept) - - //to prevent error if files are not in repository - if (boost::filesystem::exists(CT_FRACTION) && boost::filesystem::exists(CT_PLANNING)) - { - - //Target image: fraction CT, Moving image: planning CT - registrationAlgorithm prepareRegistrationRealisticScenario(CT_FRACTION, CT_PLANNING, true); - Registration3D3DPointer registrationRealisticScenario = - prepareRegistrationRealisticScenario.getRegistration(); - - auto transformRealistic = TransformationInterface::Pointer(new MatchPointTransformation( - registrationRealisticScenario)); - - io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator3(RTDOSE_FILENAME_REALISTIC.c_str()); - DoseAccessorPointer doseAccessor3(doseAccessorGenerator3.generateDoseAccessor()); - - core::GeometricInfo geoInfoRealistic = io::itk::convertToGeometricInfo( - prepareRegistrationRealisticScenario.getTargetImage()); - - //Dose is related to BP-CT, map dose to fraction CT geometry - auto aSimpleMappableDoseAccessorRealisticScenarioLinear = MappableDoseAccessorBase::Pointer( - new SimpleMappableDoseAccessor(geoInfoRealistic, - doseAccessor3, transformRealistic, interpolationLinear)); - - //combination of 0, size/2 and size to check as coordinates - std::vector coordinatesToCheckX, coordinatesToCheckY, coordinatesToCheckZ; - coordinatesToCheckX.push_back(0); - coordinatesToCheckX.push_back(geoInfoRealistic.getNumColumns() / 2); - coordinatesToCheckX.push_back(geoInfoRealistic.getNumColumns() - 1); - coordinatesToCheckY.push_back(0); - coordinatesToCheckY.push_back(geoInfoRealistic.getNumRows() / 2); - coordinatesToCheckY.push_back(geoInfoRealistic.getNumRows() - 1); - coordinatesToCheckZ.push_back(0); - coordinatesToCheckZ.push_back(geoInfoRealistic.getNumSlices() / 2); - coordinatesToCheckZ.push_back(geoInfoRealistic.getNumSlices() - 1); - - //Pixels are inside the fraction CT image and mapping should work (even if they map outside of doseAccessor3) - for (unsigned int i = 0; i < coordinatesToCheckX.size(); ++i) - { - for (unsigned int j = 0; j < coordinatesToCheckY.size(); ++j) - { - for (unsigned int k = 0; k < coordinatesToCheckZ.size(); ++k) - { - CHECK_NO_THROW(aSimpleMappableDoseAccessorRealisticScenarioLinear->getDoseAt(VoxelGridIndex3D( - coordinatesToCheckX.at(i), coordinatesToCheckY.at(j), coordinatesToCheckZ.at(k)))); - } - } - } - } - else - { - std::cout << CT_PLANNING << " and " << CT_FRACTION << - " files are missing. Rigid registration tests are not performed." << std::endl; - } + VoxelGridID currentId; + doseAccessor1GeometricInfo.convert(voxelsAsIndexToTest.at(i), currentId); + //test if the expected interpolation values are computed + CHECK_EQUAL(aSimpleMappableDoseAccessorDefault->getDoseAt(currentId), + expectedValues.at(i)); + //test if getDoseAt(VoxelGridIndex3D) and getDoseAt(VoxelGridD) lead to the same results + CHECK_EQUAL(aSimpleMappableDoseAccessorDefault->getDoseAt(currentId), + aSimpleMappableDoseAccessorDefault->getDoseAt(voxelsAsIndexToTest.at(i))); } + //test invalid voxels + VoxelGridID invalidID(doseAccessor1GeometricInfo.getNumberOfVoxels() + 1); + VoxelGridIndex3D invalidIndex(doseAccessor1GeometricInfo.getNumColumns() + 1, + doseAccessor1GeometricInfo.getNumRows() + 1, doseAccessor1GeometricInfo.getNumSlices() + 1); + CHECK_EQUAL(aSimpleMappableDoseAccessorDefault->getDoseAt(invalidID), 0.0); + CHECK_EQUAL(aSimpleMappableDoseAccessorDefault->getDoseAt(invalidIndex), 0.0); + CHECK_THROW_EXPLICIT(aSimpleMappableDoseAccessorNoPadding->getDoseAt(invalidID), + core::MappingOutsideOfImageException); + CHECK_THROW_EXPLICIT(aSimpleMappableDoseAccessorNoPadding->getDoseAt(invalidIndex), + core::MappingOutsideOfImageException); + RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/interpolation/files.cmake b/testing/interpolation/files.cmake index c006bd5..5b079ca 100644 --- a/testing/interpolation/files.cmake +++ b/testing/interpolation/files.cmake @@ -1,14 +1,11 @@ SET(CPP_FILES SimpleMappableDoseAccessorTest.cpp RosuMappableDoseAccessorTest.cpp InterpolationTest.cpp + DummyTransformation.cpp rttbInterpolationTests.cpp - registrationAlgorithm.cpp - mapDemoHelloWorldRegistration1Helper.cpp ) SET(H_FILES - registrationAlgorithm.h - mapDemoHelloWorldRegistration1Helper.h - registrationTest.h + DummyTransformation.h ) diff --git a/testing/interpolation/mapDemoHelloWorldRegistration1Helper.cpp b/testing/interpolation/mapDemoHelloWorldRegistration1Helper.cpp deleted file mode 100644 index 8317ab9..0000000 --- a/testing/interpolation/mapDemoHelloWorldRegistration1Helper.cpp +++ /dev/null @@ -1,93 +0,0 @@ -// ----------------------------------------------------------------------- -// MatchPoint - DKFZ translational registration framework -// -// Copyright (c) German Cancer Research Center (DKFZ), -// Software development for Integrated Diagnostics and Therapy (SIDT). -// ALL RIGHTS RESERVED. -// See mapCopyright.txt or -// http://www.dkfz.de/en/sidt/projects/MatchPoint/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) -// Subversion HeadURL: $HeadURL: https://svn/sbr/Sources/SBR-Projects/MatchPoint/trunk/Examples/Algorithms/mapDemoHelloWorldRegistration1Helper.cpp $ -*/ - -#if defined(_MSC_VER) -#pragma warning ( disable : 4786 ) -#endif - -#include "mapDemoHelloWorldRegistration1Helper.h" - -#include "litTestImageIO.h" -#include "litCheckMacros.h" -#include "litImageTester.h" -#include "litPointSetTester.h" - -#include - -int setImageFileNames(std::string targetImage, std::string movingImage, std::string resultImage, bool isDirectory, AppGlobals& globals) -{ - globals.targetImageFileName = targetImage; - globals.movingImageFileName = movingImage; - globals.resultImageFileName = resultImage; - globals.isDirectory = isDirectory; - - return EXIT_SUCCESS; -} - -int loadData(AppGlobals& globals) -{ - if (!globals.isDirectory){ - globals.spTargetImage = - lit::TestImageIO::InternalImageType>::readImage( - globals.targetImageFileName); - } - else { - globals.spTargetImage = map::io::readImage(globals.targetImageFileName, map::io::ImageSeriesReadStyle::Dicom); - } - - if (globals.spTargetImage.IsNull()) - { - std::cerr << "Error. Cannot load target image: " << globals.targetImageFileName << std::endl; - return EXIT_FAILURE; - } - - if (!globals.isDirectory){ - globals.spMovingImage = - lit::TestImageIO::InternalImageType>::readImage( - globals.movingImageFileName); - } - else { - globals.spMovingImage = map::io::readImage(globals.movingImageFileName, map::io::ImageSeriesReadStyle::Dicom); - } - - if (globals.spMovingImage.IsNull()) - { - std::cerr << "Error. Cannot load moving image: " << globals.movingImageFileName << std::endl; - return EXIT_FAILURE; - } - - return EXIT_SUCCESS; -} - -int saveResults(AppGlobals& globals) -{ - std::cout << std::endl << "Save result data..." << std::endl; - - lit::TestImageIO::InternalImageType>::writeImage( - globals.spResultImage, globals.resultImageFileName); - - return EXIT_SUCCESS; -} - -AppGlobals::AppGlobals() -{ -}; \ No newline at end of file diff --git a/testing/interpolation/registrationAlgorithm.cpp b/testing/interpolation/registrationAlgorithm.cpp deleted file mode 100644 index ad0eaa1..0000000 --- a/testing/interpolation/registrationAlgorithm.cpp +++ /dev/null @@ -1,128 +0,0 @@ -// ----------------------------------------------------------------------- -// MatchPoint - DKFZ translational registration framework -// -// Copyright (c) German Cancer Research Center (DKFZ), -// Software development for Integrated Diagnostics and Therapy (SIDT). -// ALL RIGHTS RESERVED. -// See mapCopyright.txt or -// http://www.dkfz.de/en/sidt/projects/MatchPoint/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) -// Subversion HeadURL: $HeadURL: https://svn/sbr/Sources/SBR-Projects/MatchPoint/trunk/Examples/Algorithms/mapDemoHelloWorldRegistration1.cpp $ -*/ - -#undef MAP_SEAL_ALGORITHMS - -#include "registrationAlgorithm.h" -#include "mapDemoHelloWorldRegistration1Helper.h" - -registrationAlgorithm::registrationAlgorithm(std::string targetFilename, std::string movingFilename, bool isDirectory, bool writeResultsFile) -{ - _targetFilename = targetFilename; - _movingFilename = movingFilename; - AppGlobals globals; - int result = EXIT_FAILURE; - - _saveResults = writeResultsFile; - - setImageFileNames(_targetFilename, _movingFilename, "result.dcm", isDirectory, globals); - - //load image and landmark files - std::cout << "Load images and landmarks..." << std::endl; - - loadData(globals); - - std::cout << "Establish registration algorithm..." << std::endl; - - _spAlgorithmEuler = AlgorithmTypeEuler::New(); - //to ensure that the Geometric centers of both images are matched - //all other parameters are default - _spAlgorithmEuler->setProperty("PreinitTransform",map::core::MetaProperty::New(true)); - - //Finally we set moving and target image for that should be - //used by the image based registration algorithm - _spAlgorithmEuler->setMovingImage(globals.spMovingImage); - _spAlgorithmEuler->setTargetImage(globals.spTargetImage); - - //The algorithm is set up and ready to run... - std::cout << "Starting to determine the registration..." << std::endl; - AlgorithmTypeEuler::RegistrationType::Pointer spRegistration; - try { - spRegistration = _spAlgorithmEuler->getRegistration(); - } - catch (const map::core::ExceptionObject& e) - { - std::cerr << "caught an MatchPoint exception:\n"; - e.Print(std::cerr); - std::cerr << "\n"; - result = -1; - } - catch (const itk::ExceptionObject& e) - { - std::cerr << "caught an ITK exception:\n"; - std::cerr << e.GetFile() << ":" << e.GetLine() << ":\n" - << e.GetDescription() << "\n"; - result = -1; - } - catch (const std::exception& e) - { - std::cerr << "caught an exception:\n"; - std::cerr << e.what() << "\n"; - result = -1; - } - catch (...) - { - std::cerr << "caught an unknown exception!!!\n"; - result = -1; - } - - std::cout << std::endl << "Registration determined..." << std::endl; - //Thats all. Now we have a registration that map the moving image into the target image - //space. But the algorithm only delivers the registration as a mapping function between - //moving and target image space. So the moving image must be mapped... - - if (_saveResults){ - //define registration task - std::cout << "Define registration task..." << std::endl; - - ImageMappingTaskType::Pointer spImageTask = ImageMappingTaskType::New(); - spImageTask->setInputImage(globals.spMovingImage); - - spImageTask->setRegistration(spRegistration); - spImageTask->setResultImageDescriptor(map::core::createFieldRepresentation(*(_spAlgorithmEuler->getTargetImage().GetPointer()))); - - //Process the task - std::cout << "Process registration task..." << std::endl; - spImageTask->execute(); - - globals.spResultImage = spImageTask->getResultImage(); - - saveResults(globals); - } - -} - -vnl_vector registrationAlgorithm::getRegistrationParameters(Registration3D3DPointer reg){ - typedef map::core::ModelBasedRegistrationKernel<3, 3> ModelBasedRegistrationKernel3D3D; - - const ModelBasedRegistrationKernel3D3D* pModelBasedDirectKernel3D3D = dynamic_cast(&(reg->getDirectMapping())); - - if (pModelBasedDirectKernel3D3D) - { - ModelBasedRegistrationKernel3D3D::ParametersType params = - pModelBasedDirectKernel3D3D->getTransformModel()->getTransform()->GetParameters(); - - return params; - } - else - return vnl_vector(); -} diff --git a/testing/interpolation/registrationAlgorithm.h b/testing/interpolation/registrationAlgorithm.h deleted file mode 100644 index c54bf28..0000000 --- a/testing/interpolation/registrationAlgorithm.h +++ /dev/null @@ -1,51 +0,0 @@ -#ifndef __MYREGISTRATION_H -#define __MYREGISTRATION_H - -#include -#include - -#include "mapDemoHelloWorldRegistration1Helper.h" -#include "mapImageMappingTask.h" -#include "mapITKEuler3DMattesMIRegistrationAlgorithmTemplate.h" -#include "mapExceptionObject.h" - -/*! @class registrationAlgorithm - @brief implements a concrete registration algorithm of MatchPoint -*/ -class registrationAlgorithm -{ -private: - std::string _targetFilename; - std::string _movingFilename; - std::string _targetDirectory; - std::string _movingDirectory; - - bool _saveResults; - -public: - typedef map::core::Registration<3, 3> Registration3D3DType; - typedef Registration3D3DType::Pointer Registration3D3DPointer; - typedef map::algorithm::boxed::ITKEuler3DMattesMIRegistrationAlgorithm - AlgorithmTypeEuler; - typedef map::core::ImageMappingTask - ImageMappingTaskType; - - /*! @brief Constructor - */ - registrationAlgorithm(std::string targetFilename, std::string movingFilename, - bool isDirectory = false, bool saveResult = false); - map::core::Registration<3, 3>::Pointer getRegistration() - { - return _spAlgorithmEuler->getRegistration(); - }; - const itk::Image* getTargetImage() - { - return _spAlgorithmEuler->getTargetImage(); - }; - vnl_vector getRegistrationParameters(Registration3D3DPointer reg); - -private: - AlgorithmTypeEuler::Pointer _spAlgorithmEuler; -}; - -#endif \ No newline at end of file diff --git a/testing/interpolation/testRegistration.h b/testing/interpolation/testRegistration.h deleted file mode 100644 index c23424c..0000000 --- a/testing/interpolation/testRegistration.h +++ /dev/null @@ -1,56 +0,0 @@ -#ifndef __TESTREGISTRATION_H -#define __TESTREGISTRATION_H - -#include "mapRegistration.h" - -namespace map -{ - namespace core - { - template - class TestRegistration: public - Registration - { - public: - typedef TestRegistration Self; - typedef RegistrationBase Superclass; - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer ConstPointer; - - itkTypeMacro(TestRegistration, Registration); - itkNewMacro(Self); - - bool _limitedTarget; - double* translation; - - TestRegistration() {}; - - ~TestRegistration() - { - }; - - - virtual bool mapPointInverse(const TargetPointType& inPoint, MovingPointType& outPoint) const - { - for (unsigned int i = 0; i < VTargetDimensions; i++) - { - outPoint[i] = inPoint[i] + translation[i]; - } - - return true; - }; - - virtual bool hasLimitedTargetRepresentation() const - { - return _limitedTarget; - } - - - private: - TestRegistration(const Self& source); //purposely not implemented - void operator=(const Self&); //purposely not implemented - }; - } -} - -#endif \ No newline at end of file