diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp index c38b963..daa529b 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/boost/rttbBoostMask.cpp @@ -1,366 +1,390 @@ #include #include #include #include #include #include #include #include #include "rttbBoostMask.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace masks { namespace boost { BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, BoostMask::StructPointer aStructure, bool strict) :_geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _voxelInStructure(::boost::make_shared()){ - _isUpToDate = false; - if(! _geometricInfo ){ - throw rttb::core::NullPointerException("Error: Geometric info is NULL!"); - } - else if(!_structure){ - throw rttb::core::NullPointerException("Error: Structure is NULL!"); - } + _isUpToDate = false; + if (!_geometricInfo){ + throw rttb::core::NullPointerException("Error: Geometric info is NULL!"); + } + else if (!_structure){ + throw rttb::core::NullPointerException("Error: Structure is NULL!"); + } } void BoostMask::calcMask(){ - if(hasSelfIntersections()){ + if (hasSelfIntersections()){ if (_strict){ throw rttb::core::InvalidParameterException("Error: structure has self intersections!"); } else{ std::cerr << _structure->getLabel() << " has self intersections! It may cause errors in the voxelization results!" << std::endl; } } std::vector intersectionSlicePolygonsVector;//store the polygons of intersection slice of each index z //For each dose slice, get intersection structure slice and the contours on the structure slice - for(unsigned int indexZ=0; indexZ< _geometricInfo->getNumSlices(); indexZ++){ + for (unsigned int indexZ = 0; indexZ < _geometricInfo->getNumSlices(); indexZ++){ BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(indexZ); BoostMask::BoostPolygonVector::iterator it; - for(it = boostPolygons.begin(); it != boostPolygons.end(); ++it){ + for (it = boostPolygons.begin(); it != boostPolygons.end(); ++it){ ::boost::geometry::correct(*it); } intersectionSlicePolygonsVector.push_back(boostPolygons); } /* Use simple nearest neighbour interpolation (NNI) if dose and structure has different z spacing: If dose slice (indexZ) has no intersection with structure slice, but the last and the next do, that means the dose - slice lies between two structure slices -> use the last slice intersection contours for this slice + slice lies between two structure slices -> use the last slice intersection contours for this slice */ - for(unsigned int indexZ=1; indexZ<_geometricInfo->getNumSlices()-1; indexZ++){ - if(intersectionSlicePolygonsVector.at(indexZ).size() == 0 && intersectionSlicePolygonsVector.at(indexZ-1).size() > 0 - && intersectionSlicePolygonsVector.at(indexZ+1).size() > 0) - intersectionSlicePolygonsVector.at(indexZ) = intersectionSlicePolygonsVector.at(indexZ-1); + for (unsigned int indexZ = 1; indexZ < _geometricInfo->getNumSlices() - 1; indexZ++){ + if (intersectionSlicePolygonsVector.at(indexZ).size() == 0 && intersectionSlicePolygonsVector.at(indexZ - 1).size() > 0 + && intersectionSlicePolygonsVector.at(indexZ + 1).size() > 0) + intersectionSlicePolygonsVector.at(indexZ) = intersectionSlicePolygonsVector.at(indexZ - 1); } - for(unsigned int indexZ=0; indexZ< _geometricInfo->getNumSlices(); indexZ++){ + for (unsigned int indexZ = 0; indexZ < _geometricInfo->getNumSlices(); indexZ++){ BoostMask::BoostPolygonVector intersectionSlicePolygons = intersectionSlicePolygonsVector.at(indexZ); //Get bounding box of this dose slice VoxelIndexVector voxelList = getBoundingBox(indexZ, intersectionSlicePolygons); rttb::VoxelGridIndex3D gridIndex3D0 = voxelList.at(0); rttb::VoxelGridIndex3D gridIndex3D1 = voxelList.at(1); - for(unsigned int indexX = gridIndex3D0[0] ; indexX <= gridIndex3D1[0]; indexX ++ ){ - for(unsigned int indexY = gridIndex3D0[1]; indexY <= gridIndex3D1[1]; indexY ++){ + for (unsigned int indexX = gridIndex3D0[0]; indexX <= gridIndex3D1[0]; indexX++){ + for (unsigned int indexY = gridIndex3D0[1]; indexY <= gridIndex3D1[1]; indexY++){ rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = indexX; currentIndex[1] = indexY; currentIndex[2] = indexZ; - rttb::VoxelGridID gridID; + rttb::VoxelGridID gridID; _geometricInfo->convert(currentIndex, gridID); //Get intersection polygons of the dose voxel and the structure BoostPolygonDeque polygons = getIntersections(currentIndex, intersectionSlicePolygons); //Calc areas of all intersection polygons double intersectionArea = calcArea(polygons); double voxelSize = _geometricInfo->getSpacing()[0] * _geometricInfo->getSpacing()[1]; - if(intersectionArea > 0) { - double volumeFraction = intersectionArea/voxelSize; - if(volumeFraction > 1 && (volumeFraction-1) <= 0.0001){ + if (intersectionArea > 0) { + double volumeFraction = intersectionArea / voxelSize; + if (volumeFraction > 1 && (volumeFraction - 1) <= 0.0001){ volumeFraction = 1; } core::MaskVoxel maskVoxel = core::MaskVoxel(gridID, (volumeFraction)); _voxelInStructure->push_back(maskVoxel);//push back the mask voxel in structure } } } } sort(_voxelInStructure->begin(), _voxelInStructure->end()); } bool BoostMask::hasSelfIntersections(){ bool hasSelfIntersection = false; - for(unsigned int indexZ=0; indexZ< _geometricInfo->getNumSlices(); indexZ++){ + for (unsigned int indexZ = 0; indexZ < _geometricInfo->getNumSlices(); indexZ++){ - rttb::VoxelGridIndex3D currentIndex(0,0,indexZ); + rttb::VoxelGridIndex3D currentIndex(0, 0, indexZ); BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(currentIndex[2]); BoostMask::BoostPolygonVector::iterator it; BoostMask::BoostPolygonVector::iterator it2; - for(it = boostPolygons.begin(); it != boostPolygons.end(); ++it){ + for (it = boostPolygons.begin(); it != boostPolygons.end(); ++it){ ::boost::geometry::correct(*it); //test if polygon has self intersection - if(::boost::geometry::detail::overlay::has_self_intersections(*it, false)){ - //if (!::boost::geometry::is_valid(*it)){ + if (::boost::geometry::detail::overlay::has_self_intersections(*it, false)){ hasSelfIntersection = true; - std::cerr << _structure->getLabel() << " has self intersection polygon! Slice: "<< indexZ << ". " <getLabel() << " has self intersection polygon! Slice: " << indexZ << ". " << std::endl; break; } //test if the polygons on the same slice has intersection - for(it2 = boostPolygons.begin(); it2 != boostPolygons.end() && it2 != it; ++it2){ + for (it2 = boostPolygons.begin(); it2 != boostPolygons.end() && it2 != it; ++it2){ ::boost::geometry::correct(*it2); BoostPolygonDeque intersection; ::boost::geometry::intersection(*it, *it2, intersection); - if(intersection.size()>0){ + if (intersection.size() > 0){ //if no donut if (!(::boost::geometry::within(*it, *it2)) && !(::boost::geometry::within(*it2, *it))) { hasSelfIntersection = true; - std::cerr << _structure->getLabel() << ": Two polygons on the same slice has intersection! Slice: " << indexZ << "." << std::endl; + std::cerr << _structure->getLabel() << ": Two polygons on the same slice has intersection! Slice: " << indexZ << "." << std::endl; break; } } } - } + } } return hasSelfIntersection; } /*Get the 4 voxel index of the bounding box of the polygon in the slice with sliceNumber*/ BoostMask::VoxelIndexVector BoostMask::getBoundingBox(unsigned int sliceNumber, const BoostPolygonVector& intersectionSlicePolygons){ - if(sliceNumber<0 || sliceNumber >= _geometricInfo->getNumSlices()){ + if (sliceNumber < 0 || sliceNumber >= _geometricInfo->getNumSlices()){ throw rttb::core::InvalidParameterException(std::string("Error: slice number is invalid!")); } - rttb::VoxelGridIndex3D currentIndex(0,0,sliceNumber); + rttb::VoxelGridIndex3D currentIndex(0, 0, sliceNumber); double xMin = 0; double yMin = 0; - double xMax = 0; + double xMax = 0; double yMax = 0; BoostPolygonVector::const_iterator it; - for(it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it){ + for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it){ ::boost::geometry::model::box box; ::boost::geometry::envelope(*it, box); BoostPoint2D minPoint = box.min_corner(); BoostPoint2D maxPoint = box.max_corner(); - if(it == intersectionSlicePolygons.begin()){ + if (it == intersectionSlicePolygons.begin()){ xMin = minPoint.x(); yMin = minPoint.y(); xMax = maxPoint.x(); yMax = maxPoint.y(); } xMin = std::min(xMin, minPoint.x()); yMin = std::min(yMin, minPoint.y()); xMax = std::max(xMax, maxPoint.x()); yMax = std::max(yMax, maxPoint.y()); } rttb::WorldCoordinate3D nullWorldCoord; - _geometricInfo->indexToWorldCoordinate(VoxelGridIndex3D(0,0,0),nullWorldCoord); + _geometricInfo->indexToWorldCoordinate(VoxelGridIndex3D(0, 0, 0), nullWorldCoord); rttb::WorldCoordinate3D minWorldCoord(xMin, yMin, nullWorldCoord.z()); rttb::WorldCoordinate3D maxWorldCoord(xMax, yMax, nullWorldCoord.z()); rttb::VoxelGridIndex3D index0; rttb::VoxelGridIndex3D index1; _geometricInfo->worldCoordinateToIndex(minWorldCoord, index0); _geometricInfo->worldCoordinateToIndex(maxWorldCoord, index1); VoxelIndexVector voxelList; voxelList.push_back(index0); voxelList.push_back(index1); return voxelList; - } + } /*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){ + 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; - BoostPolygonDeque::const_iterator it2; - for(it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it){ + for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it){ area += ::boost::geometry::area(*it); } return area; } VoxelGridIndex3D BoostMask::getGridIndex3D(const core::MaskVoxel& aMaskVoxel){ - rttb::VoxelGridIndex3D gridIndex3D; + rttb::VoxelGridIndex3D gridIndex3D; _geometricInfo->convert(aMaskVoxel.getVoxelGridID(), gridIndex3D); return gridIndex3D; } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector(){ - if(!_isUpToDate){ + if (!_isUpToDate){ calcMask(); } return _voxelInStructure; } BoostMask::BoostRing2D BoostMask::convert(const rttb::PolygonType& aRTTBPolygon){ BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; - for(unsigned int i=0; igetStructureVector(); + + //get the polygons in the slice and convert to boost rings rttb::PolygonSequenceType::iterator it; - for(it = polygonSequence.begin(); it != polygonSequence.end(); ++it){ + for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it){ PolygonType rttbPolygon = *it; - if(!rttbPolygon.empty()){ + if (!rttbPolygon.empty()){ double polygonZCoor = rttbPolygon.at(0)[2]; - rttb::WorldCoordinate3D polygonPoint(0,0, polygonZCoor); + rttb::WorldCoordinate3D polygonPoint(0, 0, polygonZCoor); rttb::VoxelGridIndex3D polygonPointIndex3D; _geometricInfo->worldCoordinateToIndex(polygonPoint, polygonPointIndex3D); //if the z - if(aVoxelGridIndexZ == polygonPointIndex3D[2]){ + if (aVoxelGridIndexZ == polygonPointIndex3D[2]){ boostRingVector.push_back(convert(rttbPolygon)); } - + } } - BoostMask::BoostRingVector::iterator it1; - BoostMask::BoostRingVector::iterator it2; + return checkDonutAndConvert(boostRingVector); + } + + 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; + + } + + 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 donutIndex;//store the outer and inner ring index + 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; - unsigned int index2; - for (it1 = boostRingVector.begin(); it1 != boostRingVector.end(); it1++){ + 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 < donutIndex.size(); i++){ - if (donutIndex.at(i) == index1){ + //check if the ring is already determined as a donut + for (unsigned int i = 0; i < donutIndexVector.size(); i++){ + if (donutIndexVector.at(i) == index1){ it1IsDonut = true; break; } } + //if not jet, check now if (!it1IsDonut){ - BoostMask::BoostPolygon2D polygon2D; - ::boost::geometry::append(polygon2D, *it1); - boostPolygonVector.push_back(polygon2D);//insert the ring - - index2 = 0; - for (it2 = boostRingVector.begin(); it2 != boostRingVector.end(); it2++){ + bool it2IsDonut = false; + unsigned int index2 = 0; + for (it2 = aRingVector.begin(); it2 != aRingVector.end(); it2++, index2++){ if (it2 != it1){ - //if donut + BoostMask::BoostPolygon2D polygon2D; if (::boost::geometry::within(*it1, *it2)){ - BoostMask::BoostPolygon2D polygon2D; ::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 - boostPolygonVector.pop_back();//delete the ring - boostPolygonVector.push_back(polygon2D);//insert the donut - donutIndex.push_back(index1); - donutIndex.push_back(index2); + it2IsDonut = true; } //if donut else if (::boost::geometry::within(*it2, *it1)){ - BoostMask::BoostPolygon2D polygon2D; ::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 - boostPolygonVector.pop_back();//delete the ring - boostPolygonVector.push_back(polygon2D);//insert the donut - donutIndex.push_back(index1); - donutIndex.push_back(index2); + 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! } } - index2++; } } - index1++; } - - 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); + //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 = 0; i < donutIndexVector.size(); i++){ + if (donutIndexVector.at(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 + } + } - return polygon; + //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; } } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMask.h b/code/masks/boost/rttbBoostMask.h index 471446b..27be5fe 100644 --- a/code/masks/boost/rttbBoostMask.h +++ b/code/masks/boost/rttbBoostMask.h @@ -1,150 +1,153 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 707 $ (last changed revision) // @date $Date: 2014-09-04 16:37:24 +0200 (Do, 04 Sep 2014) $ (last change date) // @author $Author: floca $ (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 namespace rttb { namespace masks { namespace boost { /*! @class BoostMask * @brief Implementation of voxelization using boost::geometry. - * !Important! If "strict" is set to true, an exception will be thrown when the given structure has self intersection. + * @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. */ 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 std::deque BoostPolygonDeque; typedef std::vector BoostRingVector;//polygon without holes typedef std::vector BoostPolygonVector;//polygon with or without holes typedef std::vector VoxelIndexVector; GeometricInfoPointer _geometricInfo; StructPointer _structure; 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 Check if the structure has self intersections*/ bool hasSelfIntersections(); /*! @brief Get the min/max voxel index of the bounding box of the polygon in the slice with sliceNumber * @param sliceNumber slice number between 0 and _geometricInfo->getNumSlices() * @param intersectionSlicePolygons Get the polygons intersecting the slice * @exception InvalidParameterException thrown if sliceNumber < 0 or sliceNumber >= _geometricInfo->getNumSlices() * @return Return the 4 voxel index of the bounding box */ VoxelIndexVector getBoundingBox(unsigned int sliceNumber, const BoostPolygonVector& intersectionSlicePolygons); /*! @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); /*! @brief Get grid index of a mask voxel * @param aMaskVoxel A mask voxel * @return Return the 3d grid index of the mask voxel */ VoxelGridIndex3D getGridIndex3D(const core::MaskVoxel& aMaskVoxel); /*! @brief Convert RTTB polygon to boost polygon*/ BoostRing2D convert(const rttb::PolygonType& aRTTBPolygon); /*! @brief Get the intersection slice of the voxel, return the polygons (with (donut) or without holes) in the slice * @param aVoxelGridIndexZ The z grid index (slice number) of the voxel * @return Return the structure polygons intersecting the slice */ BoostPolygonVector getIntersectionSlicePolygons(const rttb::GridIndexType aVoxelGridIndexZ); /*! @brief Get the voxel 2d contour polygon*/ BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D); + /*! @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); + }; } } } #endif \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskAccessor.h b/code/masks/boost/rttbBoostMaskAccessor.h index eeafd13..374eb26 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/boost/rttbBoostMaskAccessor.h @@ -1,127 +1,127 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 484 $ (last changed revision) // @date $Date: 2014-03-26 16:16:16 +0100 (Mi, 26 Mrz 2014) $ (last change date) // @author $Author: zhangl $ (last changed by) */ #ifndef __BOOST_MASK_ACCESSOR__H #define __BOOST_MASK_ACCESSOR__H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include "rttbGenericDoseIterator.h" #include "rttbStructure.h" #include namespace rttb { namespace masks { namespace boost { /*! @class BoostMaskAccessor * @brief Using the voxelization based on boost::geometry and generate the mask accessor. - * !Important! If "strict" is set to true, an exception will be thrown when the given structure has self intersection. + * @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. */ class BoostMaskAccessor: public core::MaskAccessorInterface { public: typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::Structure::StructTypePointer StructTypePointer; private: core::GeometricInfo _geoInfo; /*! vector containing list of mask voxels*/ MaskVoxelListPointer _spRelevantVoxelVector; StructTypePointer _spStructure; IDType _maskUID; bool _strict; public: /*! @brief Constructor with a structure pointer and a geometric info pointer * @param aStructurePointer smart pointer of the structure * @param aGeometricInfoPtr smart pointer of the geometricinfo of the dose * @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 */ BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, bool strict = true); /*! @brief destructor*/ ~BoostMaskAccessor(); /*! @brief voxelization of the given structures using boost algorithms*/ void updateMask(); /*! @brief get vector containing all relevant voxels that are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(); /*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold); /*!@brief determine how a given voxel on the dose grid is masked * @param aID ID of the voxel in grid. * @param voxel Reference to the voxel. * @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined. * The relevant volume fraction will be set to zero. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const; /*!@brief determine how a given voxel on the dose grid is masked * @param aIndex 3d index of the voxel in grid. * @param voxel Reference to the voxel. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const; /*! @brief give access to GeometricInfo*/ const core::GeometricInfo& getGeometricInfo() const; /* @ brief is true if dose is on a homogeneous grid * @remark Inhomogeneous grids are not supported at the moment, but if they will be supported in the future the interface does not need to change.*/ bool isGridHomogeneous() const { return true; }; IDType getMaskUID() const { return _maskUID; }; }; } } } #endif