diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp index 58be82e..daa529b 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/boost/rttbBoostMask.cpp @@ -1,300 +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) - :_geometricInfo(aDoseGeoInfo), _structure(aStructure), _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 ){ - 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()){ - throw rttb::core::InvalidParameterException("Error: structure has self intersections!"); + 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 + 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++){ - BoostMask::BoostRingVector boostPolygons = getIntersectionSlicePolygons(indexZ); - BoostMask::BoostRingVector::iterator it; - for(it = boostPolygons.begin(); it != boostPolygons.end(); ++it){ + 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){ ::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++){ - BoostMask::BoostRingVector intersectionSlicePolygons = intersectionSlicePolygonsVector.at(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 - std::deque polygons = getIntersections(currentIndex, intersectionSlicePolygons); + 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++){ - rttb::VoxelGridIndex3D currentIndex(0,0,indexZ); + for (unsigned int indexZ = 0; indexZ < _geometricInfo->getNumSlices(); indexZ++){ + + rttb::VoxelGridIndex3D currentIndex(0, 0, indexZ); - BoostMask::BoostRingVector boostPolygons = getIntersectionSlicePolygons(currentIndex[2]); + BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(currentIndex[2]); - BoostMask::BoostRingVector::iterator it; - BoostMask::BoostRingVector::iterator it2; - for(it = boostPolygons.begin(); it != boostPolygons.end(); ++it){ + BoostMask::BoostPolygonVector::iterator it; + BoostMask::BoostPolygonVector::iterator it2; + 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)){ + if (::boost::geometry::detail::overlay::has_self_intersections(*it, false)){ hasSelfIntersection = true; - std::cout << "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){ - hasSelfIntersection = true; - std::cout << "Two polygons on the same slice has intersection! Slice: "<< indexZ << "." < 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; + 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 BoostRingVector& intersectionSlicePolygons){ - if(sliceNumber<0 || sliceNumber >= _geometricInfo->getNumSlices()){ + BoostMask::VoxelIndexVector BoostMask::getBoundingBox(unsigned int sliceNumber, const BoostPolygonVector& intersectionSlicePolygons){ + 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; - BoostRingVector::const_iterator it; - for(it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it){ + BoostPolygonVector::const_iterator 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 BoostRingVector& intersectionSlicePolygons){ + BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons){ BoostMask::BoostPolygonDeque polygonDeque; BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); ::boost::geometry::correct(voxelPolygon); - BoostRingVector::const_iterator it; - for(it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it){ - BoostRing2D contour = *it; + 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){ + 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]){ - boostPolygonVector.push_back(convert(rttbPolygon)); + if (aVoxelGridIndexZ == polygonPointIndex3D[2]){ + boostRingVector.push_back(convert(rttbPolygon)); } + } } - return boostPolygonVector; + 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]); + 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]); + BoostPoint2D point2(rttbPoint[0] + xSpacing, rttbPoint[1]); ::boost::geometry::append(polygon, point2); - BoostPoint2D point3(rttbPoint[0]+xSpacing, rttbPoint[1]+ySpacing); + BoostPoint2D point3(rttbPoint[0] + xSpacing, rttbPoint[1] + ySpacing); ::boost::geometry::append(polygon, point3); - BoostPoint2D point4(rttbPoint[0], rttbPoint[1]+ySpacing); + 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 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 = 0; i < donutIndexVector.size(); i++){ + if (donutIndexVector.at(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 = 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 + } + } + + //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 576b244..27be5fe 100644 --- a/code/masks/boost/rttbBoostMask.h +++ b/code/masks/boost/rttbBoostMask.h @@ -1,137 +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. + * @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); + 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; + 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 BoostRingVector& intersectionSlicePolygons); + 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 BoostRingVector& intersectionSlicePolygons); + 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 in the slice + /*! @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 */ - BoostRingVector getIntersectionSlicePolygons(const rttb::GridIndexType aVoxelGridIndexZ); + 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.cpp b/code/masks/boost/rttbBoostMaskAccessor.cpp index 764d645..87edd92 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.cpp +++ b/code/masks/boost/rttbBoostMaskAccessor.cpp @@ -1,162 +1,161 @@ // ----------------------------------------------------------------------- // 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) */ #include "rttbBoostMaskAccessor.h" #include "rttbBoostMask.h" #include #include #include #include namespace rttb { namespace masks { namespace boost { - BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo) - : _spStructure(aStructurePointer), _geoInfo(aGeometricInfo) + BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, bool strict) + : _spStructure(aStructurePointer), _geoInfo(aGeometricInfo), _strict(strict) { _spRelevantVoxelVector = MaskVoxelListPointer(); //generate new structure set uid ::boost::uuids::uuid id; ::boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _maskUID = "BoostMask_" + ss.str(); } BoostMaskAccessor::~BoostMaskAccessor() { }; void BoostMaskAccessor::updateMask() { MaskVoxelList newRelevantVoxelVector; if (_spRelevantVoxelVector) { return; // already calculated } - - BoostMask mask(::boost::make_shared(_geoInfo), _spStructure); + BoostMask mask(::boost::make_shared(_geoInfo), _spStructure, _strict); _spRelevantVoxelVector = mask.getRelevantVoxelVector(); return; } BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector() { // if not already generated start voxelization here updateMask(); return _spRelevantVoxelVector; } BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector(float lowerThreshold) { MaskVoxelListPointer filteredVoxelVectorPointer(new MaskVoxelList); updateMask(); // filter relevant voxels BoostMaskAccessor::MaskVoxelList::iterator it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getRelevantVolumeFraction() > lowerThreshold) { filteredVoxelVectorPointer->push_back(*it); } ++it; } // if mask calculation was not successful this is empty! return filteredVoxelVectorPointer; } bool BoostMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const { //initialize return voxel voxel.setRelevantVolumeFraction(0); //check if ID is valid if (!_geoInfo.validID(aID)) { return false; } //determine how a given voxel on the dose grid is masked if (_spRelevantVoxelVector) { BoostMaskAccessor::MaskVoxelList::iterator it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getVoxelGridID() == aID) { voxel = (*it); return true; } ++it; } } // returns false if mask was not calculated without triggering calculation (otherwise not const!) else { return false; } return false; } bool BoostMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const { //convert VoxelGridIndex3D to VoxelGridID VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getMaskAt(aVoxelGridID, voxel); } else { return false; } } const core::GeometricInfo& BoostMaskAccessor::getGeometricInfo() const { return _geoInfo; }; } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskAccessor.h b/code/masks/boost/rttbBoostMaskAccessor.h index f33d0dc..374eb26 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/boost/rttbBoostMaskAccessor.h @@ -1,118 +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 Implementation of the voxelization using boost + * @brief Using the voxelization based on boost::geometry and generate the mask accessor. + * @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 + /*! @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); + 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 diff --git a/testing/examples/CMakeLists.txt b/testing/examples/CMakeLists.txt index 83e2a6f..d1d5de7 100644 --- a/testing/examples/CMakeLists.txt +++ b/testing/examples/CMakeLists.txt @@ -1,47 +1,47 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(CORE_TEST_EXAMPLES ${EXECUTABLE_OUTPUT_PATH}/rttbExamplesTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(RTBioModelExampleTest ${CORE_TEST_EXAMPLES} RTBioModelExampleTest "${TEST_DATA_ROOT}/TestDVH/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT3.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_ah/dvh_diff_trunk6.txt" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_ah/dvh_diff_trunk8.txt") ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwoGridScaling.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwoGridScaling05.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantFiftyGridScaling01.dcm") ADD_TEST(RTDoseIndexTest ${CORE_TEST_EXAMPLES} RTDoseIndexTest "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT3.txt") ADD_TEST(RTDoseStatisticsTest ${CORE_TEST_EXAMPLES} RTDoseStatisticsTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo_withDoseGridScaling.dcm" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa101.dos.gz" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa000.vdx" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa101.pln") ADD_TEST(RTDVHTest ${CORE_TEST_EXAMPLES} RTDVHTest "${TEST_DATA_ROOT}/TestDVH/dvh_test.txt") ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo_withDoseGridScaling.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm") ADD_TEST(RTBioModelScatterPlotExampleTest ${CORE_TEST_EXAMPLES} RTBioModelScatterPlotExampleTest "${TEST_DATA_ROOT}/TestDVH/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt") ADD_TEST(VoxelizationValidationTest ${CORE_TEST_EXAMPLES} VoxelizationValidationTest "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/boostMask/" "${TEST_DATA_ROOT}/OTBMask/") ADD_TEST(DoseStatisticsIOVirtuosTest ${CORE_TEST_EXAMPLES} DoseStatisticsIOVirtuosTest "${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa000.vdx" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa000.ctx.gz" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa101.dos.gz" "${TEST_DATA_ROOT}/Virtuos/MPM_LR_aa/MPM_LR_aa101.pln" "CTV") -RTTB_CREATE_TEST_MODULE(rttbExamples DEPENDS RTTBCore RTTBAlgorithms RTTBMasks RTTBOTBMask RTTBBoostMask RTTBIndices RTTBDicomIO RTTBVirtuosIO RTTBOtherIO RTTBModels PACKAGE_DEPENDS Litmus) +RTTB_CREATE_TEST_MODULE(rttbExamples DEPENDS RTTBCore RTTBAlgorithms RTTBMasks RTTBOTBMask RTTBBoostMask RTTBIndices RTTBDicomIO RTTBVirtuosIO RTTBITKIO RTTBOtherIO RTTBModels PACKAGE_DEPENDS Litmus) diff --git a/testing/examples/VoxelizationValidationTest.cpp b/testing/examples/VoxelizationValidationTest.cpp index 77a6021..9596ee5 100644 --- a/testing/examples/VoxelizationValidationTest.cpp +++ b/testing/examples/VoxelizationValidationTest.cpp @@ -1,195 +1,199 @@ // ----------------------------------------------------------------------- // 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: 929 $ (last changed revision) // @date $Date: 2015-04-08 14:50:57 +0200 (Mi, 08 Apr 2015) $ (last change date) // @author $Author: zhangl $ (last changed by) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVHCalculator.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbOTBMaskAccessor.h" #include "rttbDVHTxtFileReader.h" #include "rttbBoostMaskAccessor.h" -//#include "rttbITKImageMaskAccessorConverter.h" -//#include "rttbImageWriter.h" +#include "rttbITKImageMaskAccessorConverter.h" +#include "rttbImageWriter.h" namespace rttb { namespace testing { /*! @brief VoxelizationValidationTest. Compare two differnt voxelizations: OTB and Boost. Check dvh maximum and minimum for each structure. Check write mask to itk file for further validation. */ int VoxelizationValidationTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: structure file name // 2: dose1 file name // 3: directory name to write boost mask of all structures // 4: directory name to write OTB mask of all structures - std::string RTSTRUCT_FILENAME ; + std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string BoostMask_DIRNAME; std::string OTBMask_DIRNAME; if (argc > 4) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; BoostMask_DIRNAME = argv[3]; OTBMask_DIRNAME = argv[4]; } OFCondition status; DcmFileFormat fileformat; /* read dicom-rt dose */ io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create a vector of MaskAccessors (one for each structure) StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); std::vector rtStructSetMaskAccessorVec; if (rtStructureSet->getNumberOfStructures() > 0) { for (int j = 0; j < rtStructureSet->getNumberOfStructures(); j++) { std::cout << j << ": " << rtStructureSet->getStructure(j)->getLabel() << std::endl; clock_t start(clock()); //create OTB MaskAccessor ::boost::shared_ptr spOTBMaskAccessor = ::boost::make_shared(rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); spOTBMaskAccessor->updateMask(); MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor); ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(spMaskAccessor, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvh = *(calc.generateDVH()); clock_t finish(clock()); std::cout << "OTB Mask Calculation time: " << finish - start << " ms" << std::endl; //Write the mask image to a file. /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! */ /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter(spOTBMaskAccessor); CHECK(itkConverter.process()); std::stringstream fileNameSstr; - fileNameSstr< (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); + CHECK_NO_THROW(boostMaskAccessorPtr->updateMask()); ::boost::shared_ptr spMaskedDoseIteratorTmp2 = ::boost::make_shared(boostMaskAccessorPtr, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator2(spMaskedDoseIteratorTmp2); rttb::core::DVHCalculator calc2(spMaskedDoseIterator2, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvh2 = *(calc2.generateDVH()); clock_t finish2(clock()); std::cout << "Boost Mask Calculation and write file time: " << finish2 - start2 << " ms" << std::endl; //Write the mask image to a file. /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! */ + /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter2(boostMaskAccessorPtr); CHECK(itkConverter2.process()); std::stringstream fileNameSstr2; - fileNameSstr2<