diff --git a/code/core/rttbGeometricInfo.cpp b/code/core/rttbGeometricInfo.cpp index f161f52..08c1dfa 100644 --- a/code/core/rttbGeometricInfo.cpp +++ b/code/core/rttbGeometricInfo.cpp @@ -1,304 +1,304 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbGeometricInfo.h" #include #include #include 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() { using pmatrix = boost::numeric::ublas::permutation_matrix; 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; } void GeometricInfo::setImageSize(const ImageSize& aSize) { setNumColumns(aSize(0)); setNumRows(aSize(1)); setNumSlices(aSize(2)); } const ImageSize GeometricInfo::getImageSize() const { return ImageSize(static_cast(getNumColumns()), static_cast(getNumRows()), static_cast(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::equalsAlmost(const GeometricInfo& another, double errorConstantGI /*= 1e-5*/) const { return (getImagePositionPatient().equalsAlmost(another.getImagePositionPatient(), errorConstantGI) && getOrientationMatrix().equalsAlmost(another.getOrientationMatrix(), errorConstantGI) && getSpacing().equalsAlmost(another.getSpacing(), errorConstantGI) && getNumColumns() == another.getNumColumns() && getNumRows() == another.getNumRows() && getNumSlices() == another.getNumSlices()); } - bool GeometricInfo::worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, + bool GeometricInfo::worldCoordinateToContinousIndex(const WorldCoordinate3D& aWorldCoordinate, ContinousVoxelGridIndex3D& aIndex) const { WorldCoordinate3D distanceToIP; distanceToIP = aWorldCoordinate - _imagePositionPatient; boost::numeric::ublas::vector result = boost::numeric::ublas::prod( _invertedOrientationMatrix, distanceToIP); boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_div(result, _spacing); aIndex = ContinousVoxelGridIndex3D(resultS(0), resultS(1), resultS(2)); //if we convert ContinousVoxelGridIndex3D (double) to VoxelGridIndex3D (unsigned int), we can't find out if it's negative. //So we have to check before. if (aIndex(0) < -0.5 || aIndex(1) < -0.5 || aIndex(2) < -0.5){ return false; } else { //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::worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, VoxelGridIndex3D& aIndex) const { ContinousVoxelGridIndex3D doubleIndex; - bool inside = worldCoordinateToGeometryCoordinate(aWorldCoordinate, doubleIndex); + bool inside = worldCoordinateToContinousIndex(aWorldCoordinate, doubleIndex); aIndex = VoxelGridIndex3D(GridIndexType(doubleIndex(0)+0.5), GridIndexType(doubleIndex(1)+0.5), GridIndexType(doubleIndex(2)+0.5)); return inside; } bool GeometricInfo::continousIndexToWorldCoordinate(const ContinousVoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const { boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_prod( aIndex, _spacing); boost::numeric::ublas::vector result = boost::numeric::ublas::prod( _orientationMatrix, resultS); aWorldCoordinate = result + _imagePositionPatient; //if we convert ContinousVoxelGridIndex3D (double) to VoxelGridIndex3D (unsigned int), we can't find out if it's negative. //So we have to check before. if (aIndex(0) < -0.5 || aIndex(1) < -0.5 || aIndex(2) < -0.5){ return false; } else { 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 { ContinousVoxelGridIndex3D indexDouble = ContinousVoxelGridIndex3D(aIndex(0), aIndex(1), aIndex(2)); return continousIndexToWorldCoordinate(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) const { VoxelGridIndex3D currentIndex; return (worldCoordinateToIndex(aWorldCoordinate, currentIndex)); } const GridSizeType GeometricInfo::getNumberOfVoxels() const { auto nVoxels = static_cast(_numberOfRows * _numberOfColumns * _numberOfFrames); return nVoxels; } 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 e694bd9..b275347 100644 --- a/code/core/rttbGeometricInfo.h +++ b/code/core/rttbGeometricInfo.h @@ -1,187 +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. // //------------------------------------------------------------------------ #ifndef __GEOMETRIC_INFO_NEW_H #define __GEOMETRIC_INFO_NEW_H #include #include #include "rttbBaseType.h" #include #include "rttbCommon.h" #include "RTTBCoreExports.h" #ifdef _MSC_VER #pragma warning(push) #pragma warning(disable: 4251) #endif 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 RTTBCore_EXPORT GeometricInfo { public: rttbClassMacroNoParent(GeometricInfo) private: WorldCoordinate3D _imagePositionPatient{ 0 }; OrientationMatrix _orientationMatrix{ 0 }; OrientationMatrix _invertedOrientationMatrix{ 0 }; SpacingVectorType3D _spacing{ 0 }; VoxelGridDimensionType _numberOfColumns{0}; VoxelGridDimensionType _numberOfRows{0}; VoxelGridDimensionType _numberOfFrames{0}; /* @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() = default; 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; }; 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 RTTBCore_EXPORT operator == (const GeometricInfo& gInfo, const GeometricInfo& gInfo1); bool equalsAlmost(const GeometricInfo& another, double errorConstantGI = 1e-5) const; /*! @brief converts world coordinates to voxel grid index. @details the voxels world coordinates are defined by spacing, orientation and imagePositionPatient. (-0.5/-0.5/-0.5) --> (0/0/0) and (0.4999/0.4999/0.4999) --> (0/0/0) define the outer coordinates of a voxel with spacing=1, orientation= x y z (identity matrix) and imagePositionPatient=(0/0/0). @sa WorldCoordinate3D VoxelGridIndex3D @note The conversion of values is done even if the target index is not inside the given voxel grid. @returns false if aWorldCoordinate is outside the voxel grid, true otherwise. */ bool worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, VoxelGridIndex3D& aIndex) const; /*! @brief converts world coordinates to double geometry coordinate. @details This is needed because of a double precision voxel coordinate system for voxelization. The world coordinate of the image position patient is the center of the first voxel (0.0/0.0/0.0). (-0.5/-0.5/-0.5) --> (-0.5/-0.5/-0.5) and (0.4999/0.4999/0.4999) --> (0.4999/0.4999/0.4999) with spacing=1, orientation= x y z (identity matrix) and imagePositionPatient=(0/0/0). @sa WorldCoordinate3D, ContinousVoxelGridIndex3D @note The conversion of values is done even if the target index is not inside the given voxel grid. @returns false if aWorldCoordinate is outside the voxel grid, true otherwise. */ - bool worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, + bool worldCoordinateToContinousIndex(const WorldCoordinate3D& aWorldCoordinate, ContinousVoxelGridIndex3D& aIndex) const; /*! @brief converts double geometry coordinate to world coordinates. @details This is needed because of a double precision voxel coordinate system for voxelization. The world coordinate of the image position patient is the center of the first voxel (0.0/0.0/0.0). (-0.5/-0.5/-0.5) --> (-0.5/-0.5/-0.5) and (5.5/3.2/1.0) --> (5.5/3.2/1.0) with spacing=1, orientation= x y z (identity matrix) and imagePositionPatient=(0/0/0). @sa ContinousVoxelGridIndex3D, WorldCoordinate3D @note The conversion of values is done even if the target index is not inside the given voxel grid. @returns false if aWorldCoordinate is outside the voxel grid, true otherwise. */ bool continousIndexToWorldCoordinate(const ContinousVoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const; /*! @brief convert voxel grid index to world coordinates @details The world coordinate of the image position patient (center of the first voxel) is the center of the first voxel (0.0/0.0/0.0) (0/0/0) --> (0.0/0.0/0.0) and (1/1/2) --> (1.0/1.0/2.0) with spacing=1, orientation= x y z (identity matrix) and imagePositionPatient=(0/0/0). Thus, the center of the voxel is taken and converted. @sa VoxelGridIndex3D, WorldCoordinate3D @note The conversion of values is done even if the target index is not inside the given voxel grid. @returns false if aWorldCoordinate is outside the voxel grid, true otherwise. */ 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; 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); }; } } #ifdef _MSC_VER #pragma warning(pop) #endif #endif diff --git a/code/masks/rttbBoostMask.cpp b/code/masks/rttbBoostMask.cpp index 10e8fde..6d72f2d 100644 --- a/code/masks/rttbBoostMask.cpp +++ b/code/masks/rttbBoostMask.cpp @@ -1,558 +1,558 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include #include #include #include #include #include "rttbBoostMask.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbBoostMaskVoxelizationThread.h" namespace rttb { namespace masks { namespace boost { BoostMask::BoostMask(core::GeometricInfo::Pointer aDoseGeoInfo, core::Structure::Pointer aStructure, bool strict, unsigned int numberOfThreads) : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _numberOfThreads(numberOfThreads), _voxelizationThickness(0.0), _voxelInStructure(::boost::make_shared()) { _isUpToDate = false; if (_geometricInfo == nullptr) { throw rttb::core::NullPointerException("Error: Geometric info is nullptr!"); } else if (_structure == nullptr) { throw rttb::core::NullPointerException("Error: Structure is nullptr!"); } if (_numberOfThreads == 0) { _numberOfThreads = ::boost::thread::hardware_concurrency(); if (_numberOfThreads == 0) { throw rttb::core::InvalidParameterException("Error: detection of the number of hardware threads is not possible. Please specify number of threads for voxelization explicitly as parameter in BoostMask."); } } } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { if (!_isUpToDate) { calcMask(); } return _voxelInStructure; } void BoostMask::calcMask() { preprocessing(); voxelization(); generateMaskVoxelList(); _isUpToDate = true; } void BoostMask::preprocessing() { rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); //Convert world coordinate polygons to the polygons with geometry coordinate rttb::PolygonSequenceType geometryCoordinatePolygonVector; rttb::PolygonSequenceType::iterator it; rttb::ContinousVoxelGridIndex3D globalMaxGridIndex(std::numeric_limits::min(), std::numeric_limits::min(), std::numeric_limits::min()); rttb::ContinousVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(), _geometricInfo->getNumRows(), 0); for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) { PolygonType rttbPolygon = *it; PolygonType geometryCoordinatePolygon; //1. convert polygon to geometry coordinate polygons //2. calculate global min/max //3. check if polygon is planar if (!preprocessingPolygon(rttbPolygon, geometryCoordinatePolygon, globalMinGridIndex, globalMaxGridIndex, errorConstant)) { throw rttb::core::Exception("TiltedMaskPlaneException"); } geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon); } rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) ), GridIndexType(globalMinGridIndex(1) ), GridIndexType(globalMinGridIndex(2) )); rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) ), GridIndexType(globalMaxGridIndex(1) ), GridIndexType(globalMaxGridIndex(2) )); _globalBoundingBox.push_back(minIndex); _globalBoundingBox.push_back(maxIndex); //convert rttb polygon sequence to a map of z index and a vector of boost ring 2d (without holes) _ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); } void BoostMask::voxelization() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } BoostRingMap::iterator itMap; size_t mapSizeInAThread = _ringMap.size() / _numberOfThreads; unsigned int count = 0; unsigned int countThread = 0; BoostPolygonMap polygonMap; std::vector polygonMapVector; //check donut and convert to a map of z index and a vector of boost polygon 2d (with or without holes) for (itMap = _ringMap.begin(); itMap != _ringMap.end(); ++itMap) { //the vector of all boost 2d polygons with the same z grid index(donut polygon is accepted). BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); if (count == mapSizeInAThread && countThread < (_numberOfThreads - 1)) { polygonMapVector.push_back(polygonMap); polygonMap.clear(); count = 0; countThread++; } polygonMap.insert(std::pair((*itMap).first, polygonVector)); count++; } _voxelizationMap = ::boost::make_shared >(); polygonMapVector.push_back(polygonMap); //insert the last one //generate voxelization map, multi-threading ::boost::thread_group threads; auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); for (const auto & i : polygonMapVector) { BoostMaskVoxelizationThread t(i, _globalBoundingBox, _voxelizationMap, aMutex, _strict); threads.create_thread(t); } threads.join_all(); } void BoostMask::generateMaskVoxelList() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } //check homogeneous of the voxelization plane (the contours plane) if (!calcVoxelizationThickness(_voxelizationThickness)) { throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneous!"); } ::boost::thread_group threads; auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); unsigned int sliceNumberInAThread = _geometricInfo->getNumSlices() / _numberOfThreads; //generate mask voxel list, multi-threading for (unsigned int i = 0; i < _numberOfThreads; ++i) { unsigned int beginSlice = i * sliceNumberInAThread; unsigned int endSlice; if (i < _numberOfThreads - 1) { endSlice = (i + 1) * sliceNumberInAThread; } else { endSlice = _geometricInfo->getNumSlices(); } BoostMaskGenerateMaskVoxelListThread t(_globalBoundingBox, _geometricInfo, _voxelizationMap, _voxelizationThickness, beginSlice, endSlice, _voxelInStructure, _strict, aMutex); threads.create_thread(t); } threads.join_all(); } bool BoostMask::preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, rttb::PolygonType& geometryCoordinatePolygon, rttb::ContinousVoxelGridIndex3D& minimum, rttb::ContinousVoxelGridIndex3D& maximum, double aErrorConstant) const { double minZ = _geometricInfo->getNumSlices(); double maxZ = 0.0; for (auto worldCoordinatePoint : aRTTBPolygon) { //convert to geometry coordinate polygon rttb::ContinousVoxelGridIndex3D geometryCoordinatePoint; - _geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, geometryCoordinatePoint); + _geometricInfo->worldCoordinateToContinousIndex(worldCoordinatePoint, geometryCoordinatePoint); geometryCoordinatePolygon.push_back(geometryCoordinatePoint); //calculate the current global min/max //min and max for x if (geometryCoordinatePoint(0) < minimum(0)) { minimum(0) = geometryCoordinatePoint(0); } if (geometryCoordinatePoint(0) > maximum(0)) { maximum(0) = geometryCoordinatePoint(0); } //min and max for y if (geometryCoordinatePoint(1) < minimum(1)) { minimum(1) = geometryCoordinatePoint(1); } if (geometryCoordinatePoint(1) > maximum(1)) { maximum(1) = geometryCoordinatePoint(1); } //min and max for z if (geometryCoordinatePoint(2) < minimum(2)) { minimum(2) = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maximum(2)) { maximum(2) = geometryCoordinatePoint(2); } //check planar if (geometryCoordinatePoint(2) < minZ) { minZ = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maxZ) { maxZ = geometryCoordinatePoint(2); } } return (std::abs(maxZ - minZ) <= aErrorConstant); } BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const { BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); if (i == 0) { firstPoint = boostPoint; } ::boost::geometry::append(polygon2D, boostPoint); } ::boost::geometry::append(polygon2D, firstPoint); return polygon2D; } BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( const rttb::PolygonSequenceType& aRTTBPolygonVector) const { rttb::PolygonSequenceType::const_iterator it; BoostMask::BoostRingMap aRingMap; for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it) { rttb::PolygonType rttbPolygon = *it; double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon bool isFirstZ = true; if (!aRingMap.empty()) { auto findIt = findNearestKey(aRingMap, zIndex, errorConstant); //if the z index is found (same slice), add the polygon to vector if (findIt != aRingMap.end()) { //BoostRingVector ringVector = ; (*findIt).second.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); isFirstZ = false; } } //if it is the first z index in the map, insert vector with the polygon if (isFirstZ) { BoostRingVector ringVector; ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); aRingMap.insert(std::pair(zIndex, ringVector)); } } return aRingMap; } BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap& aBoostRingMap, double aIndex, double aErrorConstant) const { auto find = aBoostRingMap.find(aIndex); //if find a key equivalent to aIndex, found if (find != aBoostRingMap.end()) { return find; } else { auto lowerBound = aBoostRingMap.lower_bound(aIndex); //if all keys go before aIndex, check the last key if (lowerBound == aBoostRingMap.end()) { lowerBound = --aBoostRingMap.end(); } //if the lower bound very close to aIndex, found if (std::abs((*lowerBound).first - aIndex) <= aErrorConstant) { return lowerBound; } else { //if the lower bound is the beginning, not found if (lowerBound == aBoostRingMap.begin()) { return aBoostRingMap.end(); } else { auto lowerBound1 = --lowerBound;//the key before the lower bound //if the key before the lower bound very close to a Index, found if (std::abs((*lowerBound1).first - aIndex) <= aErrorConstant) { return lowerBound1; } //else, not found else { return aBoostRingMap.end(); } } } } } BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector) const { //check donut BoostMask::BoostRingVector::const_iterator it1; BoostMask::BoostRingVector::const_iterator it2; BoostMask::BoostPolygonVector boostPolygonVector; std::vector donutIndexVector;//store the outer and inner ring index BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon //Get donut index and donut polygon unsigned int index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i : donutIndexVector) { if (i == index1) { it1IsDonut = true; break; } } //if not jet, check now if (!it1IsDonut) { bool it2IsDonut = false; unsigned int index2 = 0; for (it2 = aRingVector.begin(); it2 != aRingVector.end(); ++it2, index2++) { if (it2 != it1) { BoostMask::BoostPolygon2D polygon2D; if (::boost::geometry::within(*it1, *it2)) { ::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring it2IsDonut = true; } //if donut else if (::boost::geometry::within(*it2, *it1)) { ::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring it2IsDonut = true; } if (it2IsDonut) { donutIndexVector.push_back(index1); donutIndexVector.push_back(index2); donutVector.push_back(polygon2D);//store donut polygon break;//Only store the first donut! } } } } } //Store no donut polygon to boostPolygonVector index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++) { bool it1IsDonut = false; //check if the ring is the outer or inner of a donut for (unsigned int i : donutIndexVector) { if (i == index1) { it1IsDonut = true; break; } } if (!it1IsDonut) { BoostMask::BoostPolygon2D polygon2D; ::boost::geometry::append(polygon2D, *it1); boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut } } //Append donut polygon to boostPolygonVector BoostMask::BoostPolygonVector::iterator itDonut; for (itDonut = donutVector.begin(); itDonut != donutVector.end(); ++itDonut) { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } bool BoostMask::calcVoxelizationThickness(double& aThickness) const { if (_voxelizationMap->size() <= 1) { aThickness = 1; return true; } double thickness = 0; auto it = _voxelizationMap->cbegin(); auto it2 = ++_voxelizationMap->cbegin(); for (; it != _voxelizationMap->cend() && it2 != _voxelizationMap->cend(); ++it, ++it2) { if (thickness == 0) { thickness = it2->first - it->first; } else { double curThickness = it2->first - it->first; //if not homogeneous (leave out double imprecisions), return false if (std::abs(thickness-curThickness)>errorConstant) { //return false; std::cout << "Two polygons are far from each other?" << std::endl; } } } if (thickness != 0) { aThickness = thickness; } else { aThickness = 1; } return true; } } } } diff --git a/testing/core/GeometricInfoTest.cpp b/testing/core/GeometricInfoTest.cpp index df99308..735ae9b 100644 --- a/testing/core/GeometricInfoTest.cpp +++ b/testing/core/GeometricInfoTest.cpp @@ -1,555 +1,555 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #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 4) test set/getSpacing 5) test set/getNumColumns/Rows/Slices 6) test get/setOrientationMatrix 8) test operators "==" 9) test equalsAlmost 10) test world to index coordinate conversion 11) test isInside and index to world coordinate conversion - 12) test with simple Geometry: isInside, continousIndexToWorldCoordinate(), worldCoordinateToGeometryCoordinate(), indexToWorldCoordinate() + 12) test with simple Geometry: isInside, continousIndexToWorldCoordinate(), worldCoordinateToContinousIndex(), indexToWorldCoordinate() 13) test getNumberOfVoxels 14) 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()); //4) test set/getSpacing //negative spacing does not make sense! /*!@is related to #2028 Should SpacingTypeVector/GridVolumeType/OrientationMatrix be forced to be non-negative?*/ SpacingVectorType3D expectedSpacing(4.15, 2.35, 100); 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()); //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))); WorldCoordinate3D testWorldCoordinate; ContinousVoxelGridIndex3D testDoubleIndex; ContinousVoxelGridIndex3D doubleIndex1 = ContinousVoxelGridIndex3D(0.1, 0, -0.3); const WorldCoordinate3D expectedDoubleIndex1(20.1, 100, -1000.3); ContinousVoxelGridIndex3D doubleIndex2 = ContinousVoxelGridIndex3D(11, 6, 15); //outside const WorldCoordinate3D expectedDoubleIndex2(31, 106, -985); ContinousVoxelGridIndex3D doubleIndex3 = ContinousVoxelGridIndex3D(10.1, 5.0, 3.0); // outside: Grid dimension = [10,5,3] const WorldCoordinate3D expectedDoubleIndex3(30.1, 105, -997); ContinousVoxelGridIndex3D doubleIndex4 = ContinousVoxelGridIndex3D(0.0, 0.0, 0.0); const WorldCoordinate3D expectedDoubleIndex4 = geoInfo.getImagePositionPatient(); //test double index to world coordinate geoInfo.continousIndexToWorldCoordinate(doubleIndex1, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, expectedDoubleIndex1); geoInfo.continousIndexToWorldCoordinate(doubleIndex2, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, expectedDoubleIndex2); geoInfo.continousIndexToWorldCoordinate(doubleIndex3, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, expectedDoubleIndex3); geoInfo.continousIndexToWorldCoordinate(doubleIndex4, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, expectedDoubleIndex4); - geoInfo.worldCoordinateToGeometryCoordinate(expectedDoubleIndex4, testDoubleIndex); + geoInfo.worldCoordinateToContinousIndex(expectedDoubleIndex4, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndex4); - geoInfo.worldCoordinateToGeometryCoordinate(expectedDoubleIndex3, testDoubleIndex); + geoInfo.worldCoordinateToContinousIndex(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); + geoInfo.worldCoordinateToContinousIndex(expectedDoubleIndex2, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndex2); - geoInfo.worldCoordinateToGeometryCoordinate(expectedDoubleIndex1, testDoubleIndex); + geoInfo.worldCoordinateToContinousIndex(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 //no half voxel shift anymore because we changed indexToWorldCoordinate/worldCoordinateToIndex insideTest1 = VoxelGridIndex3D(0, 0, 0); //origin (inside) const WorldCoordinate3D expectedIndex1(20, 100, -1000); insideTest2 = VoxelGridIndex3D(6, 0, 2); //inside const WorldCoordinate3D expectedIndex2(23, 94, -1000); insideTest3 = VoxelGridIndex3D(11, 6, 15); //outside const WorldCoordinate3D expectedIndex3(25.5, 55, -994); insideTest4 = VoxelGridIndex3D(10, 5, 3); // outside: Grid dimension = [10,5,3] const WorldCoordinate3D expectedIndex4(25, 91, -995); 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 with simple Geometry: isInside, continousIndexToWorldCoordinate(), worldCoordinateToGeometryCoordinate(), indexToWorldCoordinate() + //12) test with simple Geometry: isInside, continousIndexToWorldCoordinate(), worldCoordinateToContinousIndex(), indexToWorldCoordinate() core::GeometricInfo geoInfoSimple; ImageSize rttbSimpleSize = ImageSize(10, 10, 10); geoInfoSimple.setImageSize(rttbSimpleSize); SpacingVectorType3D spacingSimple(1, 1, 1); geoInfoSimple.setSpacing(spacingSimple); OrientationMatrix OMOnes; geoInfoSimple.setOrientationMatrix(OMOnes); const ContinousVoxelGridIndex3D doubleIndexPixelOutside1 = ContinousVoxelGridIndex3D(-0.501, 0.0, 0.0); const ContinousVoxelGridIndex3D doubleIndexPixelOutside2 = ContinousVoxelGridIndex3D(0.0, 9.501, 0.0); const ContinousVoxelGridIndex3D doubleIndexPixelZero1 = ContinousVoxelGridIndex3D(0.0, 0.0, 0.0); const ContinousVoxelGridIndex3D doubleIndexPixelZero2 = ContinousVoxelGridIndex3D(-0.5, -0.5, -0.5); const ContinousVoxelGridIndex3D doubleIndexPixelZero3 = ContinousVoxelGridIndex3D(0.499999, 0.499999, 0.499999); const ContinousVoxelGridIndex3D doubleIndexPixelOne1 = ContinousVoxelGridIndex3D(1.0, 0.0, 0.0); const ContinousVoxelGridIndex3D doubleIndexPixelOne2 = ContinousVoxelGridIndex3D(0.5, 0.499999, 0.499999); const ContinousVoxelGridIndex3D doubleIndexPixelOne3 = ContinousVoxelGridIndex3D(1.49, -0.5, -0.5); const ContinousVoxelGridIndex3D doubleIndexPixelLast1 = ContinousVoxelGridIndex3D(9.0, 9.0, 9.0); const ContinousVoxelGridIndex3D doubleIndexPixelLast2 = ContinousVoxelGridIndex3D(9.4999, 9.4999, 9.4999); const ContinousVoxelGridIndex3D doubleIndexPixelLast3 = ContinousVoxelGridIndex3D(8.501, 8.501, 8.501); const VoxelGridIndex3D indexPixelOutside = VoxelGridIndex3D(11, 0, 0); const VoxelGridIndex3D indexPixelZero = VoxelGridIndex3D(0, 0, 0); const VoxelGridIndex3D indexPixelOne = VoxelGridIndex3D(1, 0, 0); const VoxelGridIndex3D indexPixelLast = VoxelGridIndex3D(9, 9, 9); const WorldCoordinate3D worldCoordinateOutside1(doubleIndexPixelOutside1(0), doubleIndexPixelOutside1(1), doubleIndexPixelOutside1(2)); const WorldCoordinate3D worldCoordinateOutside2(doubleIndexPixelOutside2(0), doubleIndexPixelOutside2(1), doubleIndexPixelOutside2(2)); const WorldCoordinate3D worldCoordinatePixelZero1(doubleIndexPixelZero1(0), doubleIndexPixelZero1(1), doubleIndexPixelZero1(2)); const WorldCoordinate3D worldCoordinatePixelZero2(doubleIndexPixelZero2(0), doubleIndexPixelZero2(1), doubleIndexPixelZero2(2)); const WorldCoordinate3D worldCoordinatePixelZero3(doubleIndexPixelZero3(0), doubleIndexPixelZero3(1), doubleIndexPixelZero3(2)); const WorldCoordinate3D worldCoordinatePixelOne1(doubleIndexPixelOne1(0), doubleIndexPixelOne1(1), doubleIndexPixelOne1(2)); const WorldCoordinate3D worldCoordinatePixelOne2(doubleIndexPixelOne2(0), doubleIndexPixelOne2(1), doubleIndexPixelOne2(2)); const WorldCoordinate3D worldCoordinatePixelOne3(doubleIndexPixelOne3(0), doubleIndexPixelOne3(1), doubleIndexPixelOne3(2)); const WorldCoordinate3D worldCoordinatePixelLast1(doubleIndexPixelLast1(0), doubleIndexPixelLast1(1), doubleIndexPixelLast1(2)); const WorldCoordinate3D worldCoordinatePixelLast2(doubleIndexPixelLast2(0), doubleIndexPixelLast2(1), doubleIndexPixelLast2(2)); const WorldCoordinate3D worldCoordinatePixelLast3(doubleIndexPixelLast3(0), doubleIndexPixelLast3(1), doubleIndexPixelLast3(2)); bool isInside; isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelOutside1, testWorldCoordinate); CHECK(!isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelOutside2, testWorldCoordinate); CHECK(!isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelZero1, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelZero1); CHECK(isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelZero2, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelZero2); CHECK(isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelZero3, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelZero3); CHECK(isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelOne1, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelOne1); CHECK(isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelOne2, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelOne2); CHECK(isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelOne3, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelOne3); CHECK(isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelLast1, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelLast1); CHECK(isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelLast2, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelLast2); CHECK(isInside); isInside = geoInfoSimple.continousIndexToWorldCoordinate(doubleIndexPixelLast3, testWorldCoordinate); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelLast3); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinateOutside1, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinateOutside1, testDoubleIndex); CHECK(!isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinateOutside2, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinateOutside2, testDoubleIndex); CHECK(!isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelZero1, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelZero1, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelZero1); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelZero2, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelZero2, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelZero2); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelZero3, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelZero3, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelZero3); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelOne1, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelOne1, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelOne1); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelOne2, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelOne2, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelOne2); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelOne3, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelOne3, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelOne3); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelLast1, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelLast1, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelLast1); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelLast2, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelLast2, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelLast2); CHECK(isInside); - isInside = geoInfoSimple.worldCoordinateToGeometryCoordinate(worldCoordinatePixelLast3, testDoubleIndex); + isInside = geoInfoSimple.worldCoordinateToContinousIndex(worldCoordinatePixelLast3, testDoubleIndex); CHECK_EQUAL(testDoubleIndex, doubleIndexPixelLast3); CHECK(isInside); isInside = geoInfoSimple.indexToWorldCoordinate(indexPixelOutside, testWorldCoordinate); CHECK(!isInside); isInside = geoInfoSimple.indexToWorldCoordinate(indexPixelZero, testWorldCoordinate); CHECK(isInside); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelZero1); isInside = geoInfoSimple.indexToWorldCoordinate(indexPixelOne, testWorldCoordinate); CHECK(isInside); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelOne1); isInside = geoInfoSimple.indexToWorldCoordinate(indexPixelLast, testWorldCoordinate); CHECK(isInside); CHECK_EQUAL(testWorldCoordinate, worldCoordinatePixelLast1); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinateOutside1, testIntIndex); CHECK(!isInside); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinateOutside2, testIntIndex); CHECK(!isInside); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelZero1, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelZero); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelZero2, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelZero); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelZero3, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelZero); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelOne1, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelOne); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelOne2, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelOne); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelOne3, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelOne); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelLast1, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelLast); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelLast2, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelLast); isInside = geoInfoSimple.worldCoordinateToIndex(worldCoordinatePixelLast3, testIntIndex); CHECK(isInside); CHECK_EQUAL(testIntIndex, indexPixelLast); //13) test getNumberOfVoxels CHECK_EQUAL(expectedVoxelDims(0)*expectedVoxelDims(1)*expectedVoxelDims(2), geoInfo.getNumberOfVoxels()); //14) 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