diff --git a/code/masks/boost/rttbBoostMaskRedesign.cpp b/code/masks/boost/rttbBoostMaskRedesign.cpp index e942830..b21db9d 100644 --- a/code/masks/boost/rttbBoostMaskRedesign.cpp +++ b/code/masks/boost/rttbBoostMaskRedesign.cpp @@ -1,306 +1,504 @@ #include #include #include #include #include #include #include #include #include "rttbBoostMaskRedesign.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace masks { namespace boostRedesign { - BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, BoostMask::StructPointer aStructure, bool strict) - :_geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _voxelInStructure(::boost::make_shared()){ + BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, + BoostMask::StructPointer aStructure, bool strict) + : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), + _voxelInStructure(::boost::make_shared()) + { _isUpToDate = false; - if (!_geometricInfo){ + + if (!_geometricInfo) + { throw rttb::core::NullPointerException("Error: Geometric info is NULL!"); } - else if (!_structure){ + else if (!_structure) + { throw rttb::core::NullPointerException("Error: Structure is NULL!"); } } - BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector(){ - if (!_isUpToDate){ + BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() + { + if (!_isUpToDate) + { calcMask(); } + return _voxelInStructure; } - void BoostMask::calcMask(){ + void BoostMask::calcMask() + { preprocessing(); _isUpToDate = true; } - void BoostMask::preprocessing(){ + 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; - for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it){ + + for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) + { PolygonType rttbPolygon = *it; PolygonType geometryCoordinatePolygon = worldCoordinateToGeometryCoordinatePolygon(rttbPolygon); geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon); } //Get global bounding box rttb::DoubleVoxelGridIndex3D globalMaxGridIndex(0.0, 0.0, 0); - rttb::DoubleVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(), _geometricInfo->getNumRows(), 0); - for (it = geometryCoordinatePolygonVector.begin(); it != geometryCoordinatePolygonVector.end(); ++it){ + rttb::DoubleVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(), + _geometricInfo->getNumRows(), 0); + + for (it = geometryCoordinatePolygonVector.begin(); it != geometryCoordinatePolygonVector.end(); + ++it) + { PolygonType geometryCoordinatePolygon = *it; rttb::DoubleVoxelGridIndex3D maxGridIndex; rttb::DoubleVoxelGridIndex3D minGridIndex; //get min/max for x/y/z of the contour calcMinMax(geometryCoordinatePolygon, minGridIndex, maxGridIndex); //check tilt, if more than the error constant throw exception - if (checkTilt(minGridIndex, maxGridIndex, errorConstant)){ + if (checkTilt(minGridIndex, maxGridIndex, errorConstant)) + { throw rttb::core::Exception("TiltedMaskPlaneException"); } //get global min/max for x/y/z of all contours //min and max for x - if (minGridIndex(0) < globalMinGridIndex(0)){ + if (minGridIndex(0) < globalMinGridIndex(0)) + { globalMinGridIndex(0) = minGridIndex(0); } - if (maxGridIndex(0) > globalMaxGridIndex(0)){ + + if (maxGridIndex(0) > globalMaxGridIndex(0)) + { globalMaxGridIndex(0) = maxGridIndex(0); } + //min and max for y - if (minGridIndex(1) < globalMinGridIndex(1)){ + if (minGridIndex(1) < globalMinGridIndex(1)) + { globalMinGridIndex(1) = minGridIndex(1); } - if (maxGridIndex(1) > globalMaxGridIndex(1)){ + + if (maxGridIndex(1) > globalMaxGridIndex(1)) + { globalMaxGridIndex(1) = maxGridIndex(1); } + //min and max for z - if (minGridIndex(2) < globalMinGridIndex(2)){ + if (minGridIndex(2) < globalMinGridIndex(2)) + { globalMinGridIndex(2) = minGridIndex(2); } - if (maxGridIndex(2) > globalMaxGridIndex(2)){ + + if (maxGridIndex(2) > globalMaxGridIndex(2)) + { globalMaxGridIndex(2) = maxGridIndex(2); } - + } _globalBoundingBox.push_back(globalMinGridIndex); _globalBoundingBox.push_back(globalMaxGridIndex); - rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) + 0.5), GridIndexType(globalMinGridIndex(1) + 0.5), GridIndexType(globalMinGridIndex(1) + 0.5)); - rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) + 0.5), GridIndexType(globalMaxGridIndex(1) + 0.5), GridIndexType(globalMaxGridIndex(1) + 0.5)); + rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) + 0.5), + GridIndexType(globalMinGridIndex(1) + 0.5), GridIndexType(globalMinGridIndex(2) + 0.5)); + rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) + 0.5), + GridIndexType(globalMaxGridIndex(1) + 0.5), GridIndexType(globalMaxGridIndex(2) + 0.5)); _globalBoundingBoxSize0 = maxIndex(0) - minIndex(0) + 1; _globalBoundingBoxSize1 = maxIndex(1) - minIndex(1) + 1; + _globalBoundingBoxInt.push_back(minIndex); + _globalBoundingBoxInt.push_back(maxIndex); //convert rttb polygon sequence to a map of z index and a vector of boost ring 2d (without holes) BoostRingMap ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); //check donut and convert to a map of z index and a vector of boost polygon 2d (with or without holes) _geometryCoordinateBoostPolygonMap.clear(); BoostRingMap::iterator itMap; - for (itMap = ringMap.begin(); itMap != ringMap.end(); ++itMap){ + + for (itMap = ringMap.begin(); itMap != ringMap.end(); ++itMap) + { BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); - _geometryCoordinateBoostPolygonMap.insert(std::pair((*itMap).first, polygonVector)); + _geometryCoordinateBoostPolygonMap.insert(std::pair((*itMap).first, + polygonVector)); } } - void BoostMask::voxelization(){ + void BoostMask::voxelization() + { BoostPolygonMap::iterator it; - for (it = _geometryCoordinateBoostPolygonMap.begin(); it != _geometryCoordinateBoostPolygonMap.end(); ++it){ + + for (it = _geometryCoordinateBoostPolygonMap.begin(); + it != _geometryCoordinateBoostPolygonMap.end(); ++it) + { BoostArray2D maskArray(boost::extents[_globalBoundingBoxSize0][_globalBoundingBoxSize1]); - - - - + BoostPolygonVector boostPolygonVec = (*it).second; + + for (unsigned int x = 0; x < _globalBoundingBoxSize0; ++x) + { + for (unsigned y = 0; y < _globalBoundingBoxSize1; ++y) + { + rttb::VoxelGridIndex3D currentIndex; + currentIndex[0] = x + _globalBoundingBoxInt.at(0)[0]; + currentIndex[1] = y + _globalBoundingBoxInt.at(0)[1]; + currentIndex[2] = 0; + + //Get intersection polygons of the dose voxel and the structure + BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec); + + //Calc areas of all intersection polygons + double intersectionArea = calcArea(polygons); + + double voxelSize = _geometricInfo->getSpacing()[0] * _geometricInfo->getSpacing()[1]; + + + double volumeFraction = intersectionArea / voxelSize; + + if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) + { + volumeFraction = 1; + } + + maskArray[x][y] = volumeFraction; + } + } + //insert into voxelization map _voxelizationMap.insert(std::pair((*it).first, maskArray)); } - + } - rttb::PolygonType BoostMask::worldCoordinateToGeometryCoordinatePolygon(const rttb::PolygonType& aRTTBPolygon){ + rttb::PolygonType BoostMask::worldCoordinateToGeometryCoordinatePolygon( + const rttb::PolygonType& aRTTBPolygon) + { rttb::PolygonType geometryCoordinatePolygon; - for (unsigned int i = 0; i < aRTTBPolygon.size(); i++){ + + for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) + { rttb::WorldCoordinate3D worldCoordinatePoint = aRTTBPolygon.at(i); rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint; _geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, geometryCoordinatePoint); geometryCoordinatePolygon.push_back(geometryCoordinatePoint); } + return geometryCoordinatePolygon; } - bool BoostMask::checkTilt(const rttb::DoubleVoxelGridIndex3D& minimum, const rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant){ + bool BoostMask::checkTilt(const rttb::DoubleVoxelGridIndex3D& minimum, + const rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) + { return (abs(maximum(2) - minimum(2)) > aErrorConstant); } - void BoostMask::calcMinMax(const rttb::PolygonType& aRTTBPolygon, rttb::DoubleVoxelGridIndex3D minimum, rttb::DoubleVoxelGridIndex3D maximum){ + void BoostMask::calcMinMax(const rttb::PolygonType& aRTTBPolygon, + rttb::DoubleVoxelGridIndex3D minimum, rttb::DoubleVoxelGridIndex3D maximum) + { maximum(0) = 0.0; maximum(1) = 0.0; maximum(2) = 0.0; minimum(0) = _geometricInfo->getNumColumns(); minimum(1) = _geometricInfo->getNumRows(); minimum(2) = _geometricInfo->getNumSlices(); //get min and max for x/y/z of a contour PolygonType::const_iterator polygonIt; - for (polygonIt = aRTTBPolygon.begin(); polygonIt < aRTTBPolygon.end(); ++polygonIt){ + + for (polygonIt = aRTTBPolygon.begin(); polygonIt < aRTTBPolygon.end(); ++polygonIt) + { rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint = *polygonIt; + //min and max for x - if (geometryCoordinatePoint(0) < minimum(0)){ + if (geometryCoordinatePoint(0) < minimum(0)) + { minimum(0) = geometryCoordinatePoint(0); } - if (geometryCoordinatePoint(0) > maximum(0)){ + + if (geometryCoordinatePoint(0) > maximum(0)) + { maximum(0) = geometryCoordinatePoint(0); } + //min and max for y - if (geometryCoordinatePoint(1) < minimum(1)){ + if (geometryCoordinatePoint(1) < minimum(1)) + { minimum(1) = geometryCoordinatePoint(1); } - if (geometryCoordinatePoint(1) > maximum(1)){ + + if (geometryCoordinatePoint(1) > maximum(1)) + { maximum(1) = geometryCoordinatePoint(1); } + //min and max for z - if (geometryCoordinatePoint(2) < minimum(2)){ + if (geometryCoordinatePoint(2) < minimum(2)) + { minimum(2) = geometryCoordinatePoint(2); } - if (geometryCoordinatePoint(2) > maximum(2)){ + + if (geometryCoordinatePoint(2) > maximum(2)) + { maximum(2) = geometryCoordinatePoint(2); } } } - BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon){ + BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& + aRTTBPolygon) + { BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; - for (unsigned int i = 0; i < aRTTBPolygon.size(); i++){ + + for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) + { rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); - if (i == 0){ + + 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){ + BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( + const rttb::PolygonSequenceType& aRTTBPolygonVector) + { rttb::PolygonSequenceType::const_iterator it; BoostMask::BoostRingMap aRingMap; - for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it){ + + 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 BoostMask::BoostRingMap::iterator findIt = aRingMap.find(zIndex); + //if the z index is found (same slice), add the polygon to vector - if (findIt != aRingMap.end()){ + if (findIt != aRingMap.end()) + { BoostRingVector ringVector = (*findIt).second; ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); } //if it is the first z index in the map, insert vector with the polygon - else{ + else + { BoostRingVector ringVector; ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); aRingMap.insert(std::pair(zIndex, ringVector)); } } + return aRingMap; } - BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector){ + BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& + aRingVector) + { //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++){ + + 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 = 0; i < donutIndexVector.size(); i++){ - if (donutIndexVector.at(i) == index1){ + for (unsigned int i = 0; i < donutIndexVector.size(); i++) + { + if (donutIndexVector.at(i) == index1) + { it1IsDonut = true; break; } } + //if not jet, check now - if (!it1IsDonut){ + if (!it1IsDonut) + { bool it2IsDonut = false; unsigned int index2 = 0; - for (it2 = aRingVector.begin(); it2 != aRingVector.end(); it2++, index2++){ - if (it2 != it1){ + + for (it2 = aRingVector.begin(); it2 != aRingVector.end(); it2++, index2++) + { + if (it2 != it1) + { BoostMask::BoostPolygon2D polygon2D; - if (::boost::geometry::within(*it1, *it2)){ + + 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)){ + 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){ + + 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++){ + + 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 = 0; i < donutIndexVector.size(); i++){ - if (donutIndexVector.at(i) == index1){ + for (unsigned int i = 0; i < donutIndexVector.size(); i++) + { + if (donutIndexVector.at(i) == index1) + { it1IsDonut = true; break; } } - if (!it1IsDonut){ + + 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++){ + + for (itDonut = donutVector.begin(); itDonut != donutVector.end(); itDonut++) + { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } + BoostMask::BoostRing2D BoostMask::get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) + { + BoostMask::BoostRing2D polygon; + rttb::WorldCoordinate3D rttbPoint; + _geometricInfo->indexToWorldCoordinate(aVoxelGrid3D, rttbPoint); + + BoostPoint2D point1(rttbPoint[0], rttbPoint[1]); + ::boost::geometry::append(polygon, point1); + + double xSpacing = _geometricInfo->getSpacing()[0]; + double ySpacing = _geometricInfo->getSpacing()[1]; + + BoostPoint2D point2(rttbPoint[0] + xSpacing, rttbPoint[1]); + ::boost::geometry::append(polygon, point2); + + BoostPoint2D point3(rttbPoint[0] + xSpacing, rttbPoint[1] + ySpacing); + ::boost::geometry::append(polygon, point3); + + BoostPoint2D point4(rttbPoint[0], rttbPoint[1] + ySpacing); + ::boost::geometry::append(polygon, point4); + + ::boost::geometry::append(polygon, point1); + + return polygon; + + } + + /*Get intersection polygons of the contour and a voxel polygon*/ + BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& + aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) + { + BoostMask::BoostPolygonDeque polygonDeque; + + BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); + ::boost::geometry::correct(voxelPolygon); + + BoostPolygonVector::const_iterator it; + + for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) + { + BoostPolygon2D contour = *it; + ::boost::geometry::correct(contour); + BoostPolygonDeque intersection; + ::boost::geometry::intersection(voxelPolygon, contour, intersection); + + polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end()); + } + + return polygonDeque; + } + + /*Calculate the intersection area*/ + double BoostMask::calcArea(const BoostPolygonDeque& aPolygonDeque) + { + double area = 0; + + BoostPolygonDeque::const_iterator it; + + for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) + { + area += ::boost::geometry::area(*it); + } + + return area; + } + } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskRedesign.h b/code/masks/boost/rttbBoostMaskRedesign.h index 2712e04..75c7465 100644 --- a/code/masks/boost/rttbBoostMaskRedesign.h +++ b/code/masks/boost/rttbBoostMaskRedesign.h @@ -1,178 +1,204 @@ // ----------------------------------------------------------------------- // 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: 1127 $ (last changed revision) // @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __BOOST_MASK_H #define __BOOST_MASK_H #include "rttbBaseType.h" #include "rttbStructure.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include #include #include #include #include #include namespace rttb { namespace masks { namespace boostRedesign { /*! @class BoostMask - * @brief Implementation of voxelization using boost::geometry. + * @brief Implementation of voxelization using boost::geometry. * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. * (A structure without self interseciton means all contours of the structure have no self intersection, and - * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) - * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored - * and the mask will be calculated, however, it may cause errors in the mask results. + * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) + * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored + * and the mask will be calculated, however, it may cause errors in the mask results. */ class BoostMask { public: typedef ::boost::shared_ptr GeometricInfoPointer; typedef core::Structure::StructTypePointer StructPointer; typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; /*! @brief Constructor * @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is NULL * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. * @exception InvalidParameterException thrown if strict is true and the structure has self intersections */ BoostMask(GeometricInfoPointer aDoseGeoInfo, StructPointer aStructure, bool strict = true); /*! @brief Generate mask and return the voxels in the mask * @exception rttb::core::InvalidParameterException thrown if the structure has self intersections */ MaskVoxelListPointer getRelevantVoxelVector(); private: typedef ::boost::geometry::model::d2::point_xy BoostPoint2D; - typedef ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy > BoostPolygon2D; - typedef ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy > BoostRing2D; + typedef ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy > + BoostPolygon2D; + typedef ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy > + BoostRing2D; typedef std::deque BoostPolygonDeque; typedef std::vector BoostRingVector;//polygon without holes typedef std::vector BoostPolygonVector;//polygon with or without holes typedef std::vector VoxelIndexVector; - typedef std::map BoostPolygonMap;//map of the z index with the vector of boost 2d polygon - typedef std::map BoostRingMap;//map of the z index with the vector of boost 2d ring + typedef std::map + BoostPolygonMap;//map of the z index with the vector of boost 2d polygon + typedef std::map + BoostRingMap;//map of the z index with the vector of boost 2d ring typedef boost::multi_array BoostArray2D; typedef std::map BoostArrayMap; GeometricInfoPointer _geometricInfo; StructPointer _structure; /*! @brief The map of z index and a vector of all boost polygons with the same z index. This will be used to calculate the mask. - * Key: the double z grid index + * Key: the double z grid index * Value: the vector of all boost 2d polygons with the same z grid index (donut polygon is accepted). */ BoostPolygonMap _geometryCoordinateBoostPolygonMap; - /*! @brief The min and max index of the global bounding box. + /*! @brief The min and max index of the global bounding box. * The first index has the minimum for x/y/z of the global bounding box. * The second index has the maximum for x/y/z of the global bounding index. */ std::vector _globalBoundingBox; /*! @brief The bounding box x/y size*/ int _globalBoundingBoxSize0; int _globalBoundingBoxSize1; + VoxelIndexVector _globalBoundingBoxInt; + /*! @brief The voxelization map * key: the double z grid index - * value: the 2d mask, array[i][j] = the mask value of the position (i,j) in the global bounding box, + * value: the 2d mask, array[i][j] = the mask value of the position (i,j) in the global bounding box, * i: 0 - (_globalBoundingBoxSize0-1), j: 0 - (_globalBoundingBoxSize1-1) */ BoostArrayMap _voxelizationMap; bool _strict; //vector of the MaskVoxel inside the structure MaskVoxelListPointer _voxelInStructure; /*! @brief If the mask is up to date */ bool _isUpToDate; /*! @brief Voxelization and generate mask */ void calcMask(); /*! @brief The preprocessing step, wich consists of the following logic and Sub setps: *For all contours in a struct: *1) Transfer the contour polygons into boost::geometry structures - *1a) Convert the contur points from world coordinates into geometry coordinates. + *1a) Convert the contur points from world coordinates into geometry coordinates. *1b) get min and max for x/y/z of a contour *2) Tilt check: if difference of z_min and z_max is larger then a tolerance value -> there is a tilt. Throw rttb::TiltedMaskPlaneException. *3) Get struct-bounding-box: get x_min_struct, y_min_struct, x_max_struct, y_max_struct to define the bounding box that containes all contours of a struct in x-y-dimensions. */ void preprocessing(); /*! @brief The voxelization step, wich computes the voxelization planes (in x/y) for all contours of an struct. *For each contour (that is in the z-Range of the reference geometry) of the struct: *1) Allocate result array (voxelization plane) based on the bounding box (see Proprocessing Step 3) - *2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry). + *2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry). *3) Add result Array (key is the z-Value of the contoure) */ void voxelization(); /*! @brief Convert the rttb polygon with world corrdinate to the rttb polygon with double geometry coordinate */ rttb::PolygonType worldCoordinateToGeometryCoordinatePolygon(const rttb::PolygonType& aRTTBPolygon); /*! @brief Check if the polygon with the minimum and maximum is tilted more than a error constant @return Return true if tilted > aErrorConstant */ - bool checkTilt(const rttb::DoubleVoxelGridIndex3D& minimum, const rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant); + bool checkTilt(const rttb::DoubleVoxelGridIndex3D& minimum, + const rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant); /*! @brief Calculate minimum and maximum for x/y/z of the polygon */ - void calcMinMax(const rttb::PolygonType& aRTTBPolygon, rttb::DoubleVoxelGridIndex3D minimum, rttb::DoubleVoxelGridIndex3D maximum); + void calcMinMax(const rttb::PolygonType& aRTTBPolygon, rttb::DoubleVoxelGridIndex3D minimum, + rttb::DoubleVoxelGridIndex3D maximum); /*! @brief If 2 rttb polygons in the vector build a donut, convert the 2 rttb polygons to a donut boost polygon, other rttb polygons unchanged convert to boost ring*/ - + /*! @brief Convert a rttb 3d polygon to a 2d boost ring*/ BoostRing2D convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon); /*! @brief Convert a rttb 3d polygon to a map of z index with a vector of boost 2d ring, because of tilt check use the first z index of the polygon as the map key*/ - BoostRingMap convertRTTBPolygonSequenceToBoostRingMap(const rttb::PolygonSequenceType& aRTTBPolygonVector); - + BoostRingMap convertRTTBPolygonSequenceToBoostRingMap(const rttb::PolygonSequenceType& + aRTTBPolygonVector); + /*! @brief If 2 rings in the vector build a donut, convert the 2 rings to a donut polygon, other rings unchanged*/ BoostPolygonVector checkDonutAndConvert(const BoostRingVector& aRingVector); + + /*! @brief Get the voxel 2d contour polygon*/ + BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D); + + /*! @brief Get intersection polygons of the contour and a voxel polygon + * @param aVoxelIndex3D The 3d grid index of the voxel + * @param intersectionSlicePolygons The polygons of the slice intersecting the voxel + * @return Return all intersetion polygons of the structure and the voxel + */ + BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, + const BoostPolygonVector& intersectionSlicePolygons); + + /*! @brief Calculate the area of all polygons + * @param aPolygonDeque The deque of polygons + * @return Return the area of all polygons + */ + double calcArea(const BoostPolygonDeque& aPolygonDeque); }; } } } #endif \ No newline at end of file