diff --git a/code/core/rttbBaseType.h b/code/core/rttbBaseType.h index 7a6ca49..8219b4e 100644 --- a/code/core/rttbBaseType.h +++ b/code/core/rttbBaseType.h @@ -1,717 +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. */ class UnsignedIndex2D: public boost::numeric::ublas::vector { 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. */ class UnsignedIndex3D: public boost::numeric::ublas::vector { 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; } 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)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const UnsignedIndex3D& aVector) { 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) {} WorldCoordinate2D(const WorldCoordinate value) : boost::numeric::ublas::vector(2, value) {} 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(); } friend bool operator==(const WorldCoordinate2D& wc1, const WorldCoordinate2D& wc2) { if (wc1.size() != wc2.size()) { return false; } for (int i = 0; i < wc1.size(); i++) { if (wc1(i) != wc2(i)) { return false; } } 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) {} 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; } WorldCoordinate3D(const WorldCoordinate3D& w): boost::numeric::ublas::vector(3) { (*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) WorldCoordinate3D cross(WorldCoordinate3D aVector) const { 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); return result; } WorldCoordinate2D getXY() const { WorldCoordinate2D result; result(0) = (*this)(0); result(1) = (*this)(1); return result; } 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); } WorldCoordinate3D& operator=(const boost::numeric::ublas::vector wc) { (*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)); } WorldCoordinate3D operator+(const boost::numeric::ublas::vector wc) { 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)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const WorldCoordinate3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; + typedef WorldCoordinate3D DoubleVoxelGridIndex3D; + 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. */ class SpacingVectorType3D: public boost::numeric::ublas::vector { 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; } SpacingVectorType3D(const SpacingVectorType3D& w): boost::numeric::ublas::vector(3) { (*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(); } SpacingVectorType3D& operator=(const SpacingVectorType3D& wc) { (*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); } SpacingVectorType3D& operator=(const boost::numeric::ublas::vector wc) { (*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)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const SpacingVectorType3D& aVector) { s << "[ " << aVector(0) << ", " << aVector(1) << ", " << aVector(2) << " ]"; return s; } }; /*! @class OrientationMatrix @brief Used to store image orientation information */ class OrientationMatrix : public boost::numeric::ublas::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; } } 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++) { for (std::size_t n = 0; n < anOrientationMatrix.size2(); n++) { if (!((*this)(m, n) == anOrientationMatrix(m, n))) { return false; } } }// end element comparison } else { return false; } } else { return false; } 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) { s << anOrientationMatrix(m, n); } else { s << ", " << anOrientationMatrix(m, n); } } s << " ]"; } s << " ]"; return s; } }; /*! base for 2D and 3D VoxelIndex; Therefore required beside VoxelGridID */ 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) {} VoxelGridIndex3D(const GridIndexType value) : boost::numeric::ublas::vector(3, value) {} VoxelGridIndex3D(const GridIndexType xValue, const GridIndexType yValue, const GridIndexType zValue) : boost::numeric::ublas::vector(3, xValue) { (*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(); } VoxelGridIndex3D& operator=(const UnsignedIndex3D& ui) { (*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; } 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 VoxelGridIndex3D& aVector) { 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) {} VoxelGridIndex2D(const GridIndexType value) : boost::numeric::ublas::vector(2, value) {} 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(); } 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)) { return false; } } return true; } friend std::ostream& operator<<(std::ostream& s, const VoxelGridIndex2D& aVector) { 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 boost::numeric::ublas::vector VoxelGridDimensionVectorType; typedef double FractionType, VoxelSizeType, DVHVoxelNumber; typedef double DoseCalcType, DoseTypeGy, GenericValueType, DoseVoxelVolumeType, VolumeType, GridVolumeType, VoxelNumberType, BEDType, LQEDType; 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 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; void Init() { x_begin = -1000000; x_end = -1000000; y_begin = -1000000; y_end = -1000000; index_begin(0) = 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 #endif diff --git a/code/core/rttbGeometricInfo.cpp b/code/core/rttbGeometricInfo.cpp index efb2a14..7d31cde 100644 --- a/code/core/rttbGeometricInfo.cpp +++ b/code/core/rttbGeometricInfo.cpp @@ -1,303 +1,330 @@ // ----------------------------------------------------------------------- // 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 "rttbGeometricInfo.h" namespace rttb { namespace core { void GeometricInfo::setSpacing(const SpacingVectorType3D& aSpacingVector) { _spacing = aSpacingVector; } const SpacingVectorType3D& GeometricInfo::getSpacing() const { return _spacing; } void GeometricInfo::setImagePositionPatient(const WorldCoordinate3D& aImagePositionPatient) { _imagePositionPatient = aImagePositionPatient; } const WorldCoordinate3D& GeometricInfo::getImagePositionPatient() const { return _imagePositionPatient; } void GeometricInfo::setOrientationMatrix(const OrientationMatrix& anOrientationMatrix) { _orientationMatrix = anOrientationMatrix; computeInvertOrientation(); } bool GeometricInfo::computeInvertOrientation() { typedef boost::numeric::ublas::permutation_matrix pmatrix; boost::numeric::ublas::matrix A(_orientationMatrix); // create a permutation matrix for the LU-factorization pmatrix pm(A.size1()); size_t res = boost::numeric::ublas::lu_factorize(A, pm); if (res != 0) { return false; } _invertedOrientationMatrix.assign(boost::numeric::ublas::identity_matrix (A.size1())); // backsubstitute to get the inverse boost::numeric::ublas::lu_substitute(A, pm, _invertedOrientationMatrix); return true; } const ImageOrientation GeometricInfo::getImageOrientationRow() const { ImageOrientation _imageOrientationRow(0); _imageOrientationRow(0) = _orientationMatrix(0, 0); _imageOrientationRow(1) = _orientationMatrix(1, 0); _imageOrientationRow(2) = _orientationMatrix(2, 0); return _imageOrientationRow; } const ImageOrientation GeometricInfo::getImageOrientationColumn() const { ImageOrientation _imageOrientationColumn(0); _imageOrientationColumn(0) = _orientationMatrix(0, 1); _imageOrientationColumn(1) = _orientationMatrix(1, 1); _imageOrientationColumn(2) = _orientationMatrix(2, 1); return _imageOrientationColumn; } void GeometricInfo::setPixelSpacingRow(const GridVolumeType aValue) { _spacing(0) = aValue; } const GridVolumeType GeometricInfo::getPixelSpacingRow() const { return _spacing(0); } void GeometricInfo::setPixelSpacingColumn(const GridVolumeType aValue) { _spacing(1) = aValue; } const GridVolumeType GeometricInfo::getPixelSpacingColumn() const { return _spacing(1); } void GeometricInfo::setSliceThickness(const GridVolumeType aValue) { _spacing(2) = aValue; } const GridVolumeType GeometricInfo::getSliceThickness() const { return _spacing(2); } void GeometricInfo::setImageSize(const ImageSize& aSize) { setNumColumns(aSize(0)); setNumRows(aSize(1)); setNumSlices(aSize(2)); } const ImageSize GeometricInfo::getImageSize() const { return ImageSize(getNumColumns(), getNumRows(), getNumSlices()); } void GeometricInfo::setNumColumns(const VoxelGridDimensionType aValue) { _numberOfColumns = aValue; } const VoxelGridDimensionType GeometricInfo::getNumColumns() const { return _numberOfColumns; } void GeometricInfo::setNumRows(const VoxelGridDimensionType aValue) { _numberOfRows = aValue; } const VoxelGridDimensionType GeometricInfo::getNumRows() const { return _numberOfRows; } void GeometricInfo::setNumSlices(const VoxelGridDimensionType aValue) { _numberOfFrames = aValue; } const VoxelGridDimensionType GeometricInfo::getNumSlices() const { return _numberOfFrames; } bool operator==(const GeometricInfo& gInfo, const GeometricInfo& gInfo1) { return (gInfo.getImagePositionPatient() == gInfo1.getImagePositionPatient() && gInfo.getOrientationMatrix() == gInfo1.getOrientationMatrix() && gInfo.getSpacing() == gInfo1.getSpacing() && gInfo.getNumColumns() == gInfo1.getNumColumns() && gInfo.getNumRows() == gInfo1.getNumRows() && gInfo.getNumSlices() == gInfo1.getNumSlices()); } - bool GeometricInfo::worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, - VoxelGridIndex3D& aIndex) - const + bool GeometricInfo::worldCoordinateToDoubleGridIndex(const WorldCoordinate3D& aWorldCoordinate, + DoubleVoxelGridIndex3D& aIndex) + const { WorldCoordinate3D temp; temp = aWorldCoordinate - _imagePositionPatient; boost::numeric::ublas::vector temp2 = boost::numeric::ublas::prod( - _invertedOrientationMatrix, - temp); + _invertedOrientationMatrix, + temp); boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_div(temp2, - _spacing); + _spacing); - aIndex = VoxelGridIndex3D(GridIndexType(resultS(0) + 0.5), GridIndexType(resultS(1) + 0.5), - GridIndexType(resultS(2) + 0.5)); + aIndex = DoubleVoxelGridIndex3D(resultS(0), resultS(1), resultS(2)); - return isInside(aIndex); + //check if it is inside + VoxelGridIndex3D indexInt = VoxelGridIndex3D(GridIndexType(aIndex(0)+0.5), GridIndexType(aIndex(1)+0.5), + GridIndexType(aIndex(2)+0.5)); + + return isInside(indexInt); } - bool GeometricInfo::indexToWorldCoordinate(const VoxelGridIndex3D& aIndex, - WorldCoordinate3D& aWorldCoordinate) + bool GeometricInfo::worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, + VoxelGridIndex3D& aIndex) const { - WorldCoordinate3D centerVoxel(aIndex(0) - 0.5, aIndex(1) - 0.5, aIndex(2) - 0.5); + DoubleVoxelGridIndex3D doubleIndex; + bool inside = worldCoordinateToDoubleGridIndex(aWorldCoordinate, doubleIndex); + + aIndex = VoxelGridIndex3D(GridIndexType(doubleIndex(0)+0.5), GridIndexType(doubleIndex(1)+0.5), + GridIndexType(doubleIndex(2)+0.5)); + + return inside; + } + + + bool GeometricInfo::DoubleGridIndexToWorldCoordinate(const DoubleVoxelGridIndex3D& aIndex, + WorldCoordinate3D& aWorldCoordinate) + const + { boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_prod( - centerVoxel, - _spacing); + aIndex, + _spacing); boost::numeric::ublas::vector result = boost::numeric::ublas::prod( - _orientationMatrix, - resultS); + _orientationMatrix, + resultS); aWorldCoordinate = result + _imagePositionPatient; - return isInside(aIndex); + VoxelGridIndex3D indexInt = VoxelGridIndex3D(GridIndexType(aIndex(0)+0.5), GridIndexType(aIndex(1)+0.5), + GridIndexType(aIndex(2)+0.5)); + return isInside(indexInt); + } + + + bool GeometricInfo::indexToWorldCoordinate(const VoxelGridIndex3D& aIndex, + WorldCoordinate3D& aWorldCoordinate) + const + { + DoubleVoxelGridIndex3D indexDouble = DoubleVoxelGridIndex3D(aIndex(0) - 0.5, aIndex(1) - 0.5, aIndex(2) - 0.5); + return DoubleGridIndexToWorldCoordinate(indexDouble, aWorldCoordinate); } bool GeometricInfo::isInside(const VoxelGridIndex3D& aIndex) const { return (aIndex(0) >= 0 && aIndex(1) >= 0 && aIndex(2) >= 0 && aIndex(0) < static_cast(_numberOfColumns) && aIndex(1) < static_cast(_numberOfRows) && aIndex(2) < static_cast(_numberOfFrames)); } bool GeometricInfo::isInside(const WorldCoordinate3D& aWorldCoordinate) { VoxelGridIndex3D currentIndex; return (worldCoordinateToIndex(aWorldCoordinate, currentIndex)); } const GridSizeType GeometricInfo::getNumberOfVoxels() const { return _numberOfRows * _numberOfColumns * _numberOfFrames; } bool GeometricInfo::convert(const VoxelGridID& gridID, VoxelGridIndex3D& gridIndex) const { if (validID(gridID)) { gridIndex(0) = gridID % getNumColumns(); VoxelGridID tempID = (gridID - gridIndex.x()) / getNumColumns(); gridIndex(1) = tempID % getNumRows(); gridIndex(2) = (tempID - gridIndex.y()) / getNumRows(); return true; } return false; } bool GeometricInfo::convert(const VoxelGridIndex3D& gridIndex, VoxelGridID& gridID) const { if ((gridIndex.x() >= static_cast(getNumColumns())) || (gridIndex.y() >= static_cast(getNumRows())) || (gridIndex.z() >= static_cast(getNumSlices()))) { return false; } else { gridID = gridIndex.z() * getNumColumns() * getNumRows() + gridIndex.y() * getNumColumns() + gridIndex.x(); return validID(gridID); } } bool GeometricInfo::validID(const VoxelGridID aID) const { return (aID >= 0 && aID < getNumberOfVoxels()); } bool GeometricInfo::validIndex(const VoxelGridIndex3D& aIndex) const { VoxelGridID aID; if (!convert(aIndex, aID)) { return false; } else { return validID(aID); } } std::ostream& operator<<(std::ostream& s, const GeometricInfo& anGeometricInfo) { s << "[ " << anGeometricInfo.getImagePositionPatient() << "; " << anGeometricInfo.getOrientationMatrix() << "; " << anGeometricInfo.getSpacing() << "; " << "; " << anGeometricInfo.getNumColumns() << "; " << anGeometricInfo.getNumRows() << "; " << anGeometricInfo.getNumSlices() << " ]"; return s; } }//end namespace core }//end namespace rttb diff --git a/code/core/rttbGeometricInfo.h b/code/core/rttbGeometricInfo.h index d6119cc..2b53109 100644 --- a/code/core/rttbGeometricInfo.h +++ b/code/core/rttbGeometricInfo.h @@ -1,172 +1,187 @@ // ----------------------------------------------------------------------- // 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. */ 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), _numberOfFrames(0), _numberOfRows(0), _numberOfColumns(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. */ 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; - /*! @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. + + /*! @brief convert world coordinates to double voxel grid index. The world coordinate of the image position patient (center of the first voxel) will be convert to the double voxel (0.0, 0.0, 0.0) + 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 worldCoordinateToDoubleGridIndex(const WorldCoordinate3D& aWorldCoordinate, + DoubleVoxelGridIndex3D& aIndex) const; + + /*! @brief convert double voxel grid index to world coordinates. The double voxel index (0.0, 0.0, 0.0) will be convert to the world coordinate of the image postion patient (center of the first voxel) + 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 DoubleGridIndexToWorldCoordinate(const DoubleVoxelGridIndex3D& aIndex, + WorldCoordinate3D& aWorldCoordinate) const; + + /*! @brief convert int voxel grid index to 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; + bool indexToWorldCoordinate(const VoxelGridIndex3D& aIndex, + 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/testing/core/GeometricInfoTest.cpp b/testing/core/GeometricInfoTest.cpp index 3041657..079b999 100644 --- a/testing/core/GeometricInfoTest.cpp +++ b/testing/core/GeometricInfoTest.cpp @@ -1,362 +1,406 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbGeometricInfo.h" namespace rttb { namespace testing { /*!@brief GeometricInfoTest - test the API of GeometricInfo @note ITK pixel indexing: Index[0] = col, Index[1] = row, Index[2] = slice. 1) test default constructor (values as expected?) 2) test set/getImagePositionPatient 3) test set/getPixelSpacingColumn/Row/SliceThickness 4) test set/getSpacing 5) test set/getNumColumns/Rows/Slices 6) test get/setOrientationMatrix 7) test getImageOrientationRow/Column 8) test operators "==" 9) test world to index coordinate conversion 10) test isInside and index to world coordinate conversion 11) test getNumberOfVoxels */ int GeometricInfoTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; //1) test default constructor (values as expected?) CHECK_NO_THROW(core::GeometricInfo()); core::GeometricInfo geoInfo; SpacingVectorType3D testNullSV(0); WorldCoordinate3D testNullWC(0); OrientationMatrix testNullOM(0); CHECK_EQUAL(testNullSV, geoInfo.getSpacing()); CHECK_EQUAL(testNullWC, geoInfo.getImagePositionPatient()); CHECK_EQUAL(testNullOM, geoInfo.getOrientationMatrix()); //2) test set/getImagePositionPatient WorldCoordinate3D testIPP(1.2, 3.4, 5.6); CHECK_NO_THROW(geoInfo.setImagePositionPatient(testIPP)); geoInfo.setImagePositionPatient(testIPP); CHECK_EQUAL(testIPP, geoInfo.getImagePositionPatient()); //3) test set/getPixelSpacingColumn/Row/SliceThickness SpacingVectorType3D expectedSpacing(1, 2.2, 3.3); CHECK_NO_THROW(geoInfo.setPixelSpacingRow(expectedSpacing(0))); CHECK_EQUAL(expectedSpacing(0), geoInfo.getPixelSpacingRow()); CHECK_EQUAL(SpacingVectorType3D(expectedSpacing(0), 0, 0), geoInfo.getSpacing()); CHECK_EQUAL(geoInfo.getPixelSpacingRow(), geoInfo.getSpacing()(0)); CHECK_NO_THROW(geoInfo.setPixelSpacingColumn(expectedSpacing(1))); CHECK_EQUAL(expectedSpacing(1), geoInfo.getPixelSpacingColumn()); CHECK_EQUAL(SpacingVectorType3D(expectedSpacing(0), expectedSpacing(1), 0), geoInfo.getSpacing()); CHECK_EQUAL(geoInfo.getPixelSpacingColumn(), geoInfo.getSpacing()(1)); CHECK_NO_THROW(geoInfo.setSliceThickness(expectedSpacing(2))); CHECK_EQUAL(expectedSpacing(2), geoInfo.getSliceThickness()); CHECK_EQUAL(expectedSpacing, geoInfo.getSpacing()); CHECK_EQUAL(geoInfo.getSliceThickness(), geoInfo.getSpacing()(2)); //4) test set/getSpacing //negative spacing does not make sense! /*!@is related to #2028 Should SpacingTypeVector/GridVolumeType/OrientationMatrix be forced to be non-negative?*/ expectedSpacing(0) = 4.15; expectedSpacing(1) = 2.35; expectedSpacing(2) = 100; CHECK_NO_THROW(geoInfo.setSpacing(expectedSpacing)); geoInfo.setSpacing(expectedSpacing); CHECK_EQUAL(expectedSpacing, geoInfo.getSpacing()); //5) test set/getNumColumns/Rows/Slices const VoxelGridIndex3D expectedVoxelDims(10, 5, 3); //CHECK_THROW(geoInfo.setNumColumns(1.2)); -> implicit conversion will prevent exception CHECK_NO_THROW(geoInfo.setNumColumns(expectedVoxelDims(0))); geoInfo.setNumColumns(expectedVoxelDims(0)); CHECK_NO_THROW(geoInfo.setNumRows(expectedVoxelDims(1))); geoInfo.setNumRows(expectedVoxelDims(1)); //CHECK_THROW(geoInfo.setNumSlices(4.2)); -> implicit conversion will prevent exception CHECK_NO_THROW(geoInfo.setNumSlices(expectedVoxelDims(2))); geoInfo.setNumSlices(expectedVoxelDims(2)); ImageSize rttbSize = geoInfo.getImageSize(); CHECK_EQUAL(rttbSize(0), geoInfo.getNumColumns()); CHECK_EQUAL(rttbSize(1), geoInfo.getNumRows()); CHECK_EQUAL(rttbSize(2), geoInfo.getNumSlices()); rttbSize = ImageSize(11, 99, 6); core::GeometricInfo geoInfo3; geoInfo3.setImageSize(rttbSize); CHECK_EQUAL(rttbSize(0), geoInfo3.getNumColumns()); CHECK_EQUAL(rttbSize(1), geoInfo3.getNumRows()); CHECK_EQUAL(rttbSize(2), geoInfo3.getNumSlices()); //6) test get/setOrientationMatrix CHECK_EQUAL(testNullOM, geoInfo.getOrientationMatrix()); OrientationMatrix testOM(0); const WorldCoordinate3D testIORow(5.5, 4.7, 3.2); const WorldCoordinate3D testIOColumn(2.5, 1.8, 9.1); WorldCoordinate3D ortho = testIORow.cross(testIOColumn); testOM(0, 0) = testIORow(0); testOM(1, 0) = testIORow(1); testOM(2, 0) = testIORow(2); CHECK_NO_THROW(geoInfo.setOrientationMatrix(testOM)); geoInfo.setOrientationMatrix(testOM); CHECK_EQUAL(testOM, geoInfo.getOrientationMatrix()); testOM(0, 1) = testIOColumn(0); testOM(1, 1) = testIOColumn(1); testOM(2, 1) = testIOColumn(2); CHECK_NO_THROW(geoInfo.setOrientationMatrix(testOM)); geoInfo.setOrientationMatrix(testOM); CHECK_EQUAL(testOM, geoInfo.getOrientationMatrix()); testOM(0, 2) = ortho(0); testOM(1, 2) = ortho(1); testOM(2, 2) = ortho(2); CHECK_NO_THROW(geoInfo.setOrientationMatrix(testOM)); geoInfo.setOrientationMatrix(testOM); CHECK_EQUAL(testOM, geoInfo.getOrientationMatrix()); //7) test getImageOrientationRow/Column CHECK_EQUAL(testIORow, geoInfo.getImageOrientationRow()); CHECK_EQUAL(testIOColumn, geoInfo.getImageOrientationColumn()); //8) test operators "==" core::GeometricInfo geoInfo2; CHECK_EQUAL(geoInfo, geoInfo); CHECK(!(geoInfo == geoInfo2)); CHECK_EQUAL(geoInfo.getOrientationMatrix(), testOM); CHECK(!(geoInfo.getOrientationMatrix() == testNullOM)); //9) test world to index coordinate converion //use unit matrix as orientation matrix CHECK_NO_THROW(geoInfo.setOrientationMatrix(OrientationMatrix())); //origin (inside) WorldCoordinate3D insideTestWC1 = geoInfo.getImagePositionPatient(); //inside const VoxelGridIndex3D expectedIndex(8, 3, 2); WorldCoordinate3D insideTestWC2(expectedIndex(0)*expectedSpacing(0) + testIPP(0), expectedIndex(1)*expectedSpacing(1) + testIPP(1), expectedIndex(2)*expectedSpacing(2) + testIPP(2)); //outside WorldCoordinate3D insideTestWC3(-33.12, 0, 14); // outside (dimension of grid) WorldCoordinate3D insideTestWC4(expectedVoxelDims(0)*expectedSpacing(0) + testIPP(0), expectedVoxelDims(1)*expectedSpacing(1) + testIPP(1), expectedVoxelDims(2)*expectedSpacing(2) + testIPP(2)); CHECK(geoInfo.isInside(insideTestWC1)); CHECK(geoInfo.isInside(insideTestWC2)); CHECK(!(geoInfo.isInside(insideTestWC3))); CHECK(!(geoInfo.isInside(insideTestWC4))); VoxelGridIndex3D testConvert(0); CHECK(geoInfo.worldCoordinateToIndex(insideTestWC1, testConvert)); CHECK(geoInfo.isInside(testConvert)); CHECK_EQUAL(VoxelGridIndex3D(0), testConvert); CHECK(geoInfo.worldCoordinateToIndex(insideTestWC2, testConvert)); CHECK(geoInfo.isInside(testConvert)); CHECK_EQUAL(expectedIndex, testConvert); CHECK(!(geoInfo.worldCoordinateToIndex(insideTestWC3, testConvert))); //CHECK_EQUAL(VoxelGridIndex3D(0),testConvert); //if value is in a negative grid position it will be converted //to a very large unrelated number. CHECK(!(geoInfo.isInside(testConvert))); CHECK(!(geoInfo.worldCoordinateToIndex(insideTestWC4, testConvert))); CHECK_EQUAL(expectedVoxelDims, testConvert); CHECK(!(geoInfo.isInside(testConvert))); //use a more complicated orientation matrix OrientationMatrix newOrientation(0); newOrientation(0, 0) = 0.5; newOrientation(1, 2) = -3; newOrientation(2, 1) = 1; CHECK_NO_THROW(geoInfo.setOrientationMatrix(newOrientation)); testIPP = WorldCoordinate3D(20, 100, -1000); CHECK_NO_THROW(geoInfo.setImagePositionPatient(testIPP)); CHECK_NO_THROW(geoInfo.setSpacing(SpacingVectorType3D(1))); //values for testing were generated with a dedicated MeVisLab routine insideTestWC1 = geoInfo.getImagePositionPatient(); //origin (inside) const VoxelGridIndex3D expectedIndexWC1(0, 0, 0); insideTestWC2 = WorldCoordinate3D(22.5, 97, -998); //inside const VoxelGridIndex3D expectedIndexWC2(5, 2, 1); insideTestWC3 = WorldCoordinate3D(26, 88, -996); //outside const VoxelGridIndex3D expectedIndexWC3(12, 4, 4); insideTestWC4 = WorldCoordinate3D(25, 91, -995); // outside: Grid dimension = [10,5,3] const VoxelGridIndex3D expectedIndexWC4(10, 5, 3); CHECK(geoInfo.isInside(insideTestWC1)); CHECK_EQUAL(geoInfo.isInside(insideTestWC1), geoInfo.isInside(expectedIndexWC1)); CHECK(geoInfo.isInside(insideTestWC2)); CHECK_EQUAL(geoInfo.isInside(insideTestWC2), geoInfo.isInside(expectedIndexWC2)); CHECK(!(geoInfo.isInside(insideTestWC3))); CHECK_EQUAL(geoInfo.isInside(insideTestWC3), geoInfo.isInside(expectedIndexWC3)); CHECK(!(geoInfo.isInside(insideTestWC4))); CHECK_EQUAL(geoInfo.isInside(insideTestWC4), geoInfo.isInside(expectedIndexWC4)); CHECK(geoInfo.worldCoordinateToIndex(insideTestWC1, testConvert)); CHECK(geoInfo.isInside(testConvert)); CHECK_EQUAL(expectedIndexWC1, testConvert); CHECK(geoInfo.worldCoordinateToIndex(insideTestWC2, testConvert)); CHECK(geoInfo.isInside(testConvert)); CHECK_EQUAL(expectedIndexWC2, testConvert); CHECK(!(geoInfo.worldCoordinateToIndex(insideTestWC3, testConvert))); CHECK(!(geoInfo.isInside(testConvert))); CHECK_EQUAL(expectedIndexWC3, testConvert); CHECK(!(geoInfo.worldCoordinateToIndex(insideTestWC4, testConvert))); CHECK(!(geoInfo.isInside(testConvert))); CHECK_EQUAL(expectedIndexWC4, testConvert); //10) test isInside and index to world coordinate converion //use unit matrix as orientation matrix CHECK_NO_THROW(geoInfo.setOrientationMatrix(OrientationMatrix())); VoxelGridIndex3D insideTest1(0, 0, 0); //origin (inside) VoxelGridIndex3D insideTest2(2, 3, 1); //inside VoxelGridIndex3D insideTest3(0, 6, 14); //outside VoxelGridIndex3D insideTest4 = expectedVoxelDims; // outside CHECK(geoInfo.isInside(insideTest1)); CHECK(geoInfo.isInside(insideTest2)); CHECK(!(geoInfo.isInside(insideTest3))); CHECK(!(geoInfo.isInside(insideTest4))); WorldCoordinate3D testInside(0); CHECK(geoInfo.indexToWorldCoordinate(insideTest1, testInside)); CHECK(geoInfo.isInside(testInside)); //CHECK_EQUAL(geoInfo.getImagePositionPatient(),testInside); //half voxel shift prevents equality! CHECK(geoInfo.indexToWorldCoordinate(insideTest2, testInside)); CHECK(geoInfo.isInside(testInside)); CHECK(!(geoInfo.indexToWorldCoordinate(insideTest3, testInside))); CHECK(!(geoInfo.isInside(testInside))); CHECK(!(geoInfo.indexToWorldCoordinate(insideTest4, testInside))); CHECK(!(geoInfo.isInside(testInside))); + + DoubleVoxelGridIndex3D doubleIndex1 = DoubleVoxelGridIndex3D(0.1, 0, -0.3); + const WorldCoordinate3D expectedDoubleIndex1(20.1, 100, -1000.3); + DoubleVoxelGridIndex3D doubleIndex2 = DoubleVoxelGridIndex3D(11, 6, 15); //outside + const WorldCoordinate3D expectedDoubleIndex2(31, 106, -985); + DoubleVoxelGridIndex3D doubleIndex3 = DoubleVoxelGridIndex3D(9.6, 5.0, 3.0); // outside: Grid dimension = [10,5,3] + const WorldCoordinate3D expectedDoubleIndex3(29.6, 105, -997); + DoubleVoxelGridIndex3D doubleIndex4 = DoubleVoxelGridIndex3D(0.0, 0.0, 0.0); + const WorldCoordinate3D expectedDoubleIndex4 = geoInfo.getImagePositionPatient(); + + //test double index to world coordinate + WorldCoordinate3D test; + geoInfo.DoubleGridIndexToWorldCoordinate(doubleIndex1, test); + CHECK_EQUAL(test, expectedDoubleIndex1); + geoInfo.DoubleGridIndexToWorldCoordinate(doubleIndex2, test); + CHECK_EQUAL(test, expectedDoubleIndex2); + geoInfo.DoubleGridIndexToWorldCoordinate(doubleIndex3, test); + CHECK_EQUAL(test, expectedDoubleIndex3); + geoInfo.DoubleGridIndexToWorldCoordinate(doubleIndex4, test); + CHECK_EQUAL(test, expectedDoubleIndex4); + + DoubleVoxelGridIndex3D testDoubleIndex; + geoInfo.worldCoordinateToDoubleGridIndex(expectedDoubleIndex4, testDoubleIndex); + CHECK_EQUAL(testDoubleIndex, doubleIndex4); + geoInfo.worldCoordinateToDoubleGridIndex(expectedDoubleIndex3, testDoubleIndex); + CHECK_CLOSE(testDoubleIndex(0), doubleIndex3(0), errorConstant); + CHECK_CLOSE(testDoubleIndex(1), doubleIndex3(1), errorConstant); + CHECK_CLOSE(testDoubleIndex(2), doubleIndex3(2), errorConstant); + geoInfo.worldCoordinateToDoubleGridIndex(expectedDoubleIndex2, testDoubleIndex); + CHECK_EQUAL(testDoubleIndex, doubleIndex2); + geoInfo.worldCoordinateToDoubleGridIndex(expectedDoubleIndex1, testDoubleIndex); + CHECK_CLOSE(testDoubleIndex(0), doubleIndex1(0), errorConstant); + CHECK_CLOSE(testDoubleIndex(1), doubleIndex1(1), errorConstant); + CHECK_CLOSE(testDoubleIndex(2), doubleIndex1(2), errorConstant); + + VoxelGridIndex3D testIntIndex; + geoInfo.worldCoordinateToIndex(expectedDoubleIndex4, testIntIndex); + CHECK_EQUAL(testIntIndex, insideTest1); + geoInfo.worldCoordinateToIndex(expectedDoubleIndex1, testIntIndex); + CHECK_EQUAL(testIntIndex, insideTest1); + geoInfo.worldCoordinateToIndex(expectedDoubleIndex3, testIntIndex); + CHECK_EQUAL(testIntIndex, expectedVoxelDims); + //use a more complicated orientation matrix newOrientation = OrientationMatrix(0); newOrientation(0, 0) = 0.5; newOrientation(1, 2) = -3; newOrientation(2, 1) = 1; CHECK_NO_THROW(geoInfo.setOrientationMatrix(newOrientation)); testIPP = WorldCoordinate3D(20, 100, -1000); CHECK_NO_THROW(geoInfo.setImagePositionPatient(testIPP)); CHECK_NO_THROW(geoInfo.setSpacing(SpacingVectorType3D(1))); //values for testing were generated with a dedicated MeVisLab routine //expected world coordinates need to accommodate the half voxel shift //caused by the use of the central position in the voxel for conversion. insideTest1 = VoxelGridIndex3D(0, 0, 0); //origin (inside) const WorldCoordinate3D expectedIndex1(19.75, 101.5, -1000.5); insideTest2 = VoxelGridIndex3D(6, 0, 2); //inside const WorldCoordinate3D expectedIndex2(22.75, 95.5, -1000.5); insideTest3 = VoxelGridIndex3D(11, 6, 15); //outside const WorldCoordinate3D expectedIndex3(25.25, 56.5, -994.5); insideTest4 = VoxelGridIndex3D(10, 5, 3); // outside: Grid dimension = [10,5,3] const WorldCoordinate3D expectedIndex4(24.75, 92.5, -995.5); CHECK(geoInfo.isInside(insideTest1)); CHECK_EQUAL(geoInfo.isInside(insideTest1), geoInfo.isInside(expectedIndex1)); CHECK(geoInfo.isInside(insideTest2)); CHECK_EQUAL(geoInfo.isInside(insideTest2), geoInfo.isInside(expectedIndex2)); CHECK(!(geoInfo.isInside(insideTest3))); CHECK_EQUAL(geoInfo.isInside(insideTest3), geoInfo.isInside(expectedIndex3)); CHECK(!(geoInfo.isInside(insideTest4))); CHECK_EQUAL(geoInfo.isInside(insideTest4), geoInfo.isInside(expectedIndex4)); CHECK(geoInfo.indexToWorldCoordinate(insideTest1, testInside)); CHECK(geoInfo.isInside(testInside)); CHECK_EQUAL(expectedIndex1, testInside); CHECK(geoInfo.indexToWorldCoordinate(insideTest2, testInside)); CHECK(geoInfo.isInside(testInside)); CHECK_EQUAL(expectedIndex2, testInside); CHECK(!(geoInfo.indexToWorldCoordinate(insideTest3, testInside))); CHECK(!(geoInfo.isInside(testInside))); CHECK_EQUAL(expectedIndex3, testInside); CHECK(!(geoInfo.indexToWorldCoordinate(insideTest4, testInside))); CHECK(!(geoInfo.isInside(testInside))); CHECK_EQUAL(expectedIndex4, testInside); + //11) test getNumberOfVoxels CHECK_EQUAL(expectedVoxelDims(0)*expectedVoxelDims(1)*expectedVoxelDims(2), geoInfo.getNumberOfVoxels()); //12) Test convert, validID and validIndex geoInfo.setNumColumns(50); geoInfo.setNumRows(30); geoInfo.setNumSlices(40); VoxelGridIndex3D startIndex(0, 0, 0); VoxelGridID startId(0); VoxelGridIndex3D endIndex(geoInfo.getNumColumns() - 1, geoInfo.getNumRows() - 1, geoInfo.getNumSlices() - 1); VoxelGridID endId((geoInfo.getNumColumns()*geoInfo.getNumRows()*geoInfo.getNumSlices()) - 1); VoxelGridIndex3D indexInvalid(geoInfo.getNumColumns(), geoInfo.getNumRows(), geoInfo.getNumSlices()); VoxelGridID idInvalid(geoInfo.getNumColumns()*geoInfo.getNumRows()*geoInfo.getNumSlices()); CHECK(geoInfo.validID(startId)); CHECK(geoInfo.validIndex(startIndex)); VoxelGridIndex3D aIndex; VoxelGridID aId; CHECK(geoInfo.convert(startIndex, aId)); CHECK(geoInfo.convert(startId, aIndex)); CHECK_EQUAL(aId, startId); CHECK_EQUAL(aIndex, startIndex); CHECK(geoInfo.validID(endId)); CHECK(geoInfo.validIndex(endIndex)); CHECK(geoInfo.convert(endIndex, aId)); CHECK(geoInfo.convert(endId, aIndex)); CHECK_EQUAL(aId, endId); CHECK_EQUAL(aIndex, endIndex); CHECK(!geoInfo.validID(idInvalid)); CHECK(!geoInfo.validIndex(indexInvalid)); CHECK(!geoInfo.convert(idInvalid, aIndex)); CHECK(!geoInfo.convert(indexInvalid, aId)); RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb