diff --git a/code/core/rttbBaseType.h b/code/core/rttbBaseType.h index 8d48a16..ef3b651 100644 --- a/code/core/rttbBaseType.h +++ b/code/core/rttbBaseType.h @@ -1,739 +1,741 @@ // ----------------------------------------------------------------------- // 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; } bool equalsAlmost(const WorldCoordinate3D& another, double errorConstant = 1e-5) const { if (size() != another.size()) { return false; } double dist = norm_2(*this - another); return dist < errorConstant; } 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; } bool equalsAlmost(const SpacingVectorType3D& another, double errorConstant) const { if ((*this).size() != another.size()) { return false; } double dist = norm_2(*this - another); return dist < errorConstant; } 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 equalsAlmost(const OrientationMatrix& anOrientationMatrix, double errorConstant) 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 ((abs((*this)(m, n) - anOrientationMatrix(m, n)) > errorConstant)) { return false; } } }// end element comparison } else { return false; } } else { return false; } return true; } friend bool operator==(const OrientationMatrix& om1, const OrientationMatrix& om2) { return om1.equalsAlmost(om2, 0.0); } 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 c33b882..aae96fa 100644 --- a/code/core/rttbGeometricInfo.cpp +++ b/code/core/rttbGeometricInfo.cpp @@ -1,313 +1,339 @@ // ----------------------------------------------------------------------- // 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())); + (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()); + && gInfo.getOrientationMatrix() == gInfo1.getOrientationMatrix() + && gInfo.getSpacing() == gInfo1.getSpacing() && gInfo.getNumColumns() == gInfo1.getNumColumns() + && gInfo.getNumRows() == gInfo1.getNumRows() && gInfo.getNumSlices() == gInfo1.getNumSlices()); } bool GeometricInfo::equalsAlmost(const GeometricInfo& another, double errorConstant /*= 1e-5*/) const { return (getImagePositionPatient().equalsAlmost(another.getImagePositionPatient(), errorConstant) - && getOrientationMatrix().equalsAlmost(another.getOrientationMatrix(), errorConstant) - && getSpacing().equalsAlmost(another.getSpacing(), errorConstant) && getNumColumns() == another.getNumColumns() - && getNumRows() == another.getNumRows() && getNumSlices() == another.getNumSlices()); + && getOrientationMatrix().equalsAlmost(another.getOrientationMatrix(), errorConstant) + && getSpacing().equalsAlmost(another.getSpacing(), errorConstant) && getNumColumns() == another.getNumColumns() + && getNumRows() == another.getNumRows() && getNumSlices() == another.getNumSlices()); } - - bool GeometricInfo::worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, - VoxelGridIndex3D& aIndex) - const + bool GeometricInfo::worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, + DoubleVoxelGridIndex3D& aIndex) + const { - WorldCoordinate3D temp; - temp = aWorldCoordinate - _imagePositionPatient; + WorldCoordinate3D distanceToIP; + distanceToIP = aWorldCoordinate - _imagePositionPatient; - boost::numeric::ublas::vector temp2 = boost::numeric::ublas::prod( - _invertedOrientationMatrix, - temp); + boost::numeric::ublas::vector result = boost::numeric::ublas::prod( + _invertedOrientationMatrix, + distanceToIP); - boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_div(temp2, - _spacing); + boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_div(result, + _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) - const + bool GeometricInfo::worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, + VoxelGridIndex3D& aIndex) + const + { + DoubleVoxelGridIndex3D doubleIndex; + bool inside = worldCoordinateToGeometryCoordinate(aWorldCoordinate, doubleIndex); + + aIndex = VoxelGridIndex3D(GridIndexType(doubleIndex(0) + 0.5), GridIndexType(doubleIndex(1) + 0.5), + GridIndexType(doubleIndex(2) + 0.5)); + + return inside; + } + + + bool GeometricInfo::geometryCoordinateToWorldCoordinate(const DoubleVoxelGridIndex3D& aIndex, + WorldCoordinate3D& aWorldCoordinate) + const { - WorldCoordinate3D centerVoxel(aIndex(0) - 0.5, aIndex(1) - 0.5, aIndex(2) - 0.5); 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 geometryCoordinateToWorldCoordinate(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)); + && 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()))) + || (gridIndex.y() >= static_cast(getNumRows())) + || (gridIndex.z() >= static_cast(getNumSlices()))) { return false; } else { gridID = gridIndex.z() * getNumColumns() * getNumRows() + gridIndex.y() * - getNumColumns() - + gridIndex.x(); + 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() << " ]"; + 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 b2a95d9..e6aa52e 100644 --- a/code/core/rttbGeometricInfo.h +++ b/code/core/rttbGeometricInfo.h @@ -1,174 +1,189 @@ // ----------------------------------------------------------------------- // 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); bool equalsAlmost(const GeometricInfo& another, double errorConstant = 1e-5) const; /*! @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. + + /*! @brief convert world coordinates to double geometry coordinate. 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 worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, + DoubleVoxelGridIndex3D& aIndex) const; + + /*! @brief convert double geometry coordinate 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 indexToWorldCoordinate(const VoxelGridIndex3D& aIndex , + bool geometryCoordinateToWorldCoordinate(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; + /*! @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 d207c15..84d1019 100644 --- a/testing/core/GeometricInfoTest.cpp +++ b/testing/core/GeometricInfoTest.cpp @@ -1,390 +1,433 @@ // ----------------------------------------------------------------------- // 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 equalsAlmost 10) test world to index coordinate conversion 11) test isInside and index to world coordinate conversion 12) test getNumberOfVoxels 13) Test convert, validID and validIndex */ 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 equalsALmost OrientationMatrix testOM2 = testOM; SpacingVectorType3D testSPV2 = expectedSpacing; WorldCoordinate3D testIPP2 = testIPP; core::GeometricInfo testGI2, testGIEmpty; testGI2.setImagePositionPatient(testIPP2); testGI2.setOrientationMatrix(testOM2); testGI2.setSpacing(testSPV2); double smallValue = 0.000000001; testOM(0, 0) += smallValue; testSPV2(2) += smallValue; testIPP2(1) += smallValue; core::GeometricInfo testGI2similar; testGI2similar.setImagePositionPatient(testIPP2); testGI2similar.setOrientationMatrix(testOM2); testGI2similar.setSpacing(testSPV2); CHECK_EQUAL(testGI2.equalsAlmost(testGI2similar), true); CHECK_EQUAL(testGI2similar.equalsAlmost(testGI2), true); CHECK_EQUAL(testGI2.equalsAlmost(testGI2similar, smallValue * 0.001), false); CHECK_EQUAL(testGIEmpty.equalsAlmost(testGI2), false); CHECK_EQUAL(testGI2.equalsAlmost(testGIEmpty), false); //10) test world to index coordinate conversion //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); //11) test isInside and index to world coordinate conversion //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.geometryCoordinateToWorldCoordinate(doubleIndex1, test); + CHECK_EQUAL(test, expectedDoubleIndex1); + geoInfo.geometryCoordinateToWorldCoordinate(doubleIndex2, test); + CHECK_EQUAL(test, expectedDoubleIndex2); + geoInfo.geometryCoordinateToWorldCoordinate(doubleIndex3, test); + CHECK_EQUAL(test, expectedDoubleIndex3); + geoInfo.geometryCoordinateToWorldCoordinate(doubleIndex4, test); + CHECK_EQUAL(test, expectedDoubleIndex4); + + DoubleVoxelGridIndex3D testDoubleIndex; + geoInfo.worldCoordinateToGeometryCoordinate(expectedDoubleIndex4, testDoubleIndex); + CHECK_EQUAL(testDoubleIndex, doubleIndex4); + geoInfo.worldCoordinateToGeometryCoordinate(expectedDoubleIndex3, testDoubleIndex); + CHECK_CLOSE(testDoubleIndex(0), doubleIndex3(0), errorConstant); + CHECK_CLOSE(testDoubleIndex(1), doubleIndex3(1), errorConstant); + CHECK_CLOSE(testDoubleIndex(2), doubleIndex3(2), errorConstant); + geoInfo.worldCoordinateToGeometryCoordinate(expectedDoubleIndex2, testDoubleIndex); + CHECK_EQUAL(testDoubleIndex, doubleIndex2); + geoInfo.worldCoordinateToGeometryCoordinate(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); - + //12) test getNumberOfVoxels CHECK_EQUAL(expectedVoxelDims(0)*expectedVoxelDims(1)*expectedVoxelDims(2), geoInfo.getNumberOfVoxels()); //13) 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