diff --git a/code/masks/boost/files.cmake b/code/masks/boost/files.cmake index 616ac4e..161bcf0 100644 --- a/code/masks/boost/files.cmake +++ b/code/masks/boost/files.cmake @@ -1,9 +1,11 @@ SET(CPP_FILES rttbBoostMask.cpp rttbBoostMaskAccessor.cpp + rttbBoostMaskRedesign.cpp ) SET(H_FILES rttbBoostMask.h rttbBoostMaskAccessor.h + rttbBoostMaskRedesign.h ) diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp index daa529b..3791beb 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/boost/rttbBoostMask.cpp @@ -1,390 +1,392 @@ #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!"); } } void BoostMask::calcMask(){ 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++){ 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 */ 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::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++){ rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = indexX; currentIndex[1] = indexY; currentIndex[2] = indexZ; 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){ 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()); + + _isUpToDate = true; } bool BoostMask::hasSelfIntersections(){ bool hasSelfIntersection = false; for (unsigned int indexZ = 0; indexZ < _geometricInfo->getNumSlices(); 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){ ::boost::geometry::correct(*it); //test if polygon has self intersection if (::boost::geometry::detail::overlay::has_self_intersections(*it, false)){ hasSelfIntersection = true; std::cerr << _structure->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){ ::boost::geometry::correct(*it2); BoostPolygonDeque intersection; ::boost::geometry::intersection(*it, *it2, intersection); 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; 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()){ throw rttb::core::InvalidParameterException(std::string("Error: slice number is invalid!")); } rttb::VoxelGridIndex3D currentIndex(0, 0, sliceNumber); double xMin = 0; double yMin = 0; double xMax = 0; double yMax = 0; 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()){ 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); 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){ 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; } VoxelGridIndex3D BoostMask::getGridIndex3D(const core::MaskVoxel& aMaskVoxel){ rttb::VoxelGridIndex3D gridIndex3D; _geometricInfo->convert(aMaskVoxel.getVoxelGridID(), gridIndex3D); return gridIndex3D; } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector(){ if (!_isUpToDate){ calcMask(); } return _voxelInStructure; } BoostMask::BoostRing2D BoostMask::convert(const rttb::PolygonType& aRTTBPolygon){ BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++){ rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); if (i == 0){ firstPoint = boostPoint; } ::boost::geometry::append(polygon2D, boostPoint); } ::boost::geometry::append(polygon2D, firstPoint); return polygon2D; } BoostMask::BoostPolygonVector BoostMask::getIntersectionSlicePolygons(const rttb::GridIndexType aVoxelGridIndexZ){ BoostMask::BoostRingVector boostRingVector; rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); //get the polygons in the slice and convert to boost rings rttb::PolygonSequenceType::iterator it; for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it){ PolygonType rttbPolygon = *it; if (!rttbPolygon.empty()){ double polygonZCoor = rttbPolygon.at(0)[2]; rttb::WorldCoordinate3D polygonPoint(0, 0, polygonZCoor); rttb::VoxelGridIndex3D polygonPointIndex3D; _geometricInfo->worldCoordinateToIndex(polygonPoint, polygonPointIndex3D); //if the z if (aVoxelGridIndexZ == polygonPointIndex3D[2]){ boostRingVector.push_back(convert(rttbPolygon)); } } } 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 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/rttbBoostMaskRedesign.cpp b/code/masks/boost/rttbBoostMaskRedesign.cpp new file mode 100644 index 0000000..39c9ee8 --- /dev/null +++ b/code/masks/boost/rttbBoostMaskRedesign.cpp @@ -0,0 +1,287 @@ +#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()){ + + _isUpToDate = false; + if (!_geometricInfo){ + throw rttb::core::NullPointerException("Error: Geometric info is NULL!"); + } + else if (!_structure){ + throw rttb::core::NullPointerException("Error: Structure is NULL!"); + } + + } + + BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector(){ + if (!_isUpToDate){ + calcMask(); + } + return _voxelInStructure; + } + + void BoostMask::calcMask(){ + preprocessing(); + _isUpToDate = true; + } + + void BoostMask::preprocessing(){ + + rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); + + //Convert world coordinate polygons to the polygons with geometry coordinate + rttb::PolygonSequenceType geometryCoordinatePolygonVector; + rttb::PolygonSequenceType::iterator it; + 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){ + 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)){ + 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)){ + globalMinGridIndex(0) = minGridIndex(0); + } + if (maxGridIndex(0) > globalMaxGridIndex(0)){ + globalMaxGridIndex(0) = maxGridIndex(0); + } + //min and max for y + if (minGridIndex(1) < globalMinGridIndex(1)){ + globalMinGridIndex(1) = minGridIndex(1); + } + if (maxGridIndex(1) > globalMaxGridIndex(1)){ + globalMaxGridIndex(1) = maxGridIndex(1); + } + //min and max for z + if (minGridIndex(2) < globalMinGridIndex(2)){ + globalMinGridIndex(2) = minGridIndex(2); + } + if (maxGridIndex(2) > globalMaxGridIndex(2)){ + globalMaxGridIndex(2) = maxGridIndex(2); + } + + } + + _globalBoundingBox.push_back(globalMinGridIndex); + _globalBoundingBox.push_back(globalMaxGridIndex); + + //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){ + BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); + _geometryCoordinateBoostPolygonMap.insert(std::pair((*itMap).first, polygonVector)); + } + } + + rttb::PolygonType BoostMask::worldCoordinateToGeometryCoordinatePolygon(const rttb::PolygonType& aRTTBPolygon){ + rttb::PolygonType geometryCoordinatePolygon; + 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){ + return (abs(maximum(2) - minimum(2)) > aErrorConstant); + } + + 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){ + rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint = *polygonIt; + //min and max for x + if (geometryCoordinatePoint(0) < minimum(0)){ + minimum(0) = geometryCoordinatePoint(0); + } + if (geometryCoordinatePoint(0) > maximum(0)){ + maximum(0) = geometryCoordinatePoint(0); + } + //min and max for y + if (geometryCoordinatePoint(1) < minimum(1)){ + minimum(1) = geometryCoordinatePoint(1); + } + if (geometryCoordinatePoint(1) > maximum(1)){ + maximum(1) = geometryCoordinatePoint(1); + } + //min and max for z + if (geometryCoordinatePoint(2) < minimum(2)){ + minimum(2) = geometryCoordinatePoint(2); + } + if (geometryCoordinatePoint(2) > maximum(2)){ + maximum(2) = geometryCoordinatePoint(2); + } + } + } + + + BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon){ + BoostMask::BoostRing2D polygon2D; + BoostPoint2D firstPoint; + for (unsigned int i = 0; i < aRTTBPolygon.size(); i++){ + rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); + BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); + if (i == 0){ + firstPoint = boostPoint; + } + ::boost::geometry::append(polygon2D, boostPoint); + } + ::boost::geometry::append(polygon2D, firstPoint); + return polygon2D; + } + + BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap(const rttb::PolygonSequenceType& aRTTBPolygonVector){ + rttb::PolygonSequenceType::const_iterator it; + BoostMask::BoostRingMap aRingMap; + for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it){ + rttb::PolygonType rttbPolygon = *it; + double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon + BoostMask::BoostRingMap::iterator findIt = aRingMap.find(zIndex); + //if the z index is found (same slice), add the polygon to vector + 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{ + BoostRingVector ringVector; + ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); + aRingMap.insert(std::pair(zIndex, ringVector)); + } + } + return aRingMap; + } + + 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/rttbBoostMaskRedesign.h b/code/masks/boost/rttbBoostMaskRedesign.h new file mode 100644 index 0000000..967d5db --- /dev/null +++ b/code/masks/boost/rttbBoostMaskRedesign.h @@ -0,0 +1,155 @@ +// ----------------------------------------------------------------------- +// 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 + +namespace rttb +{ + namespace masks + { + namespace boostRedesign + { + /*! @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, 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; + 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 + + 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 + * 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. + * 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; + + 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. + *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 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); + + /*! @brief Calculate minimum and maximum for x/y/z of the polygon + */ + 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); + + /*! @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/testing/masks/boost/BoostMaskRedesignTest.cpp b/testing/masks/boost/BoostMaskRedesignTest.cpp new file mode 100644 index 0000000..4319f06 --- /dev/null +++ b/testing/masks/boost/BoostMaskRedesignTest.cpp @@ -0,0 +1,74 @@ +// ----------------------------------------------------------------------- +// 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) +*/ + +#include +#include + +#include "litCheckMacros.h" + +#include "rttbBaseType.h" + +#include "../../core/DummyStructure.h" +#include "../../core/DummyDoseAccessor.h" +#include "rttbBoostMaskRedesign.h" + + +namespace rttb +{ + namespace testing + { + + /*! @brief BoostMaskRedesignTest. + 1) test constructors + 2) test getRelevantVoxelVector + 3) test getMaskAt + */ + int BoostMaskRedesignTest(int argc, char* argv[]) + { + PREPARE_DEFAULT_TEST_REPORTING; + + typedef core::Structure::StructTypePointer StructTypePointer; + + // generate test structure set + boost::shared_ptr spTestDoseAccessor = + boost::make_shared(); + + DummyStructure myStructGenerator(spTestDoseAccessor->getGeometricInfo()); + + GridIndexType zPlane = 4; + core::Structure myTestStruct = myStructGenerator.CreateRectangularStructureCentered(zPlane); + StructTypePointer spMyStruct = boost::make_shared(myTestStruct); + boost::shared_ptr geometricPtr = boost::make_shared + (spTestDoseAccessor->getGeometricInfo()); + + //1) test BoostMask constructor + CHECK_NO_THROW(rttb::masks::boostRedesign::BoostMask(geometricPtr, spMyStruct)); + rttb::masks::boostRedesign::BoostMask boostMask = rttb::masks::boostRedesign::BoostMask(geometricPtr, spMyStruct); + + //2) test getRelevantVoxelVector + CHECK_NO_THROW(boostMask.getRelevantVoxelVector()); + + + RETURN_AND_REPORT_TEST_SUCCESS; + } + }//testing +}//rttb + diff --git a/testing/masks/boost/CMakeLists.txt b/testing/masks/boost/CMakeLists.txt index 8bc11ea..6dd75d4 100644 --- a/testing/masks/boost/CMakeLists.txt +++ b/testing/masks/boost/CMakeLists.txt @@ -1,21 +1,22 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(Boost_Mask_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbBoostMaskTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(BoostMaskAccessorTest ${Boost_Mask_TESTS} BoostMaskAccessorTest ) - +ADD_TEST(BoostMaskRedesignTest ${Boost_Mask_TESTS} BoostMaskRedesignTest +) RTTB_CREATE_TEST_MODULE(rttbBoostMask DEPENDS RTTBDicomIO RTTBMasks RTTBBoostMask PACKAGE_DEPENDS Boost BoostBinaries Litmus DCMTK) diff --git a/testing/masks/boost/files.cmake b/testing/masks/boost/files.cmake index c401dd6..ea3c821 100644 --- a/testing/masks/boost/files.cmake +++ b/testing/masks/boost/files.cmake @@ -1,13 +1,14 @@ SET(CPP_FILES ../../core/DummyStructure.cpp ../../core/CreateTestStructure.cpp ../../core/DummyDoseAccessor.cpp - BoostMaskAccessorTest.cpp + BoostMaskAccessorTest.cpp + BoostMaskRedesignTest.cpp rttbBoostMaskTests.cpp ) SET(H_FILES ../../core/DummyStructure.h ../../core/CreateTestStructure.h ../../core/DummyDoseAccessor.h ) diff --git a/testing/masks/boost/rttbBoostMaskTests.cpp b/testing/masks/boost/rttbBoostMaskTests.cpp index dcf6e16..1f552c0 100644 --- a/testing/masks/boost/rttbBoostMaskTests.cpp +++ b/testing/masks/boost/rttbBoostMaskTests.cpp @@ -1,63 +1,64 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ // this file defines the rttbAlgorithmsTests for the test driver // and all it expects is that you have a function called RegisterTests #if defined(_MSC_VER) #pragma warning ( disable : 4786 ) #endif #include "litMultiTestsMain.h" namespace rttb { namespace testing { void registerTests() { LIT_REGISTER_TEST(BoostMaskAccessorTest); + LIT_REGISTER_TEST(BoostMaskRedesignTest); } } } int main(int argc, char* argv[]) { int result = 0; rttb::testing::registerTests(); try { result = lit::multiTestsMain(argc, argv); } catch (const std::exception& /*e*/) { result = -1; } catch (...) { result = -1; } return result; }