diff --git a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp index 114c484..ee08a0f 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp +++ b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp @@ -1,143 +1,134 @@ #include "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbInvalidParameterException.h" #include namespace rttb { namespace masks { namespace boostRedesign { BoostMaskGenerateMaskVoxelListThread::BoostMaskGenerateMaskVoxelListThread( - int aGlobalBoundingBoxSize0, - int aGlobalBoundingBoxSize1, - const rttb::VoxelGridIndex3D& aMinIndex, - const rttb::VoxelGridIndex3D& aMaxIndex, + const VoxelIndexVector& aGlobalBoundingBox, GeometricInfoPointer aGeometricInfo, const BoostArrayMap& aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, MaskVoxelQueuePointer aResultMaskVoxelQueue) : - _globalBoundingBoxSize0(aGlobalBoundingBoxSize0), - _globalBoundingBoxSize1(aGlobalBoundingBoxSize1), - _minIndex(aMinIndex), _maxIndex(aMaxIndex), _geometricInfo(aGeometricInfo), + _globalBoundingBox(aGlobalBoundingBox), _geometricInfo(aGeometricInfo), _voxelizationMap(aVoxelizationMap), _voxelizationThickness(aVoxelizationThickness), _beginSlice(aBeginSlice), _endSlice(aEndSlice), _resultMaskVoxelQueue(aResultMaskVoxelQueue) {} void BoostMaskGenerateMaskVoxelListThread::operator()() { - std::cout << "BoostMaskGenerateMaskVoxelListThread " << _beginSlice << " - " << _endSlice << - ": running" << std::endl; for (unsigned int indexZ = _beginSlice; indexZ < _endSlice; ++indexZ) { + rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); + rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); + int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; + int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; //calculate weight vector std::map weightVectorForZ; calcWeightVector(indexZ, weightVectorForZ); //For each x,y, calc sum of all voxelization plane, use weight vector - for (unsigned int x = 0; x < _globalBoundingBoxSize0; ++x) + for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) { - for (unsigned int y = 0; y < _globalBoundingBoxSize1; ++y) + for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) { rttb::VoxelGridIndex3D currentIndex; - currentIndex[0] = x + _minIndex[0]; - currentIndex[1] = y + _minIndex[1]; + currentIndex[0] = x + minIndex[0]; + currentIndex[1] = y + minIndex[1]; currentIndex[2] = indexZ; rttb::VoxelGridID gridID; _geometricInfo->convert(currentIndex, gridID); double volumeFraction = 0; std::map::iterator it; + BoostArrayMap::iterator itMap; - for (it = weightVectorForZ.begin(); it != weightVectorForZ.end(); ++it) + for (it = weightVectorForZ.begin(), itMap = _voxelizationMap.begin(); it != weightVectorForZ.end() + && itMap != _voxelizationMap.end(); ++it, ++itMap) { - BoostArrayMap::iterator findVoxelizationIt = _voxelizationMap.find((*it).first); double weight = (*it).second; - if (findVoxelizationIt == _voxelizationMap.end()) - { - throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneus!"); - } - - BoostArray2D voxelizationArray = (*findVoxelizationIt).second; + BoostArray2D voxelizationArray = (*itMap).second; //calc sum of all voxelization plane, use weight volumeFraction += voxelizationArray[x][y] * weight; } if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } //insert mask voxel if volumeFraction > 0 if (volumeFraction > 0) { core::MaskVoxel* maskVoxelPtr = new core::MaskVoxel(gridID, volumeFraction); _resultMaskVoxelQueue->push(maskVoxelPtr);//push back the mask voxel in structure } } } } - std::cout << "BoostMaskGenerateMaskVoxelListThread" << _beginSlice << " - " << _endSlice << - ": finished" << std::endl; } void BoostMaskGenerateMaskVoxelListThread::calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const { BoostArrayMap::const_iterator it; double indexZMin = aIndexZ - 0.5; double indexZMax = aIndexZ + 0.5; for (it = _voxelizationMap.begin(); it != _voxelizationMap.end(); ++it) { double voxelizationPlaneIndexMin = (*it).first - 0.5 * _voxelizationThickness; double voxelizationPlaneIndexMax = (*it).first + 0.5 * _voxelizationThickness; double weight = 0; if ((voxelizationPlaneIndexMin < indexZMin) && (voxelizationPlaneIndexMax > indexZMin)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = voxelizationPlaneIndexMax - indexZMin; } else { weight = 1; } } else if ((voxelizationPlaneIndexMin >= indexZMin) && (voxelizationPlaneIndexMin < indexZMax)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = _voxelizationThickness; } else { weight = indexZMax - voxelizationPlaneIndexMin; } } weightVector.insert(std::pair((*it).first, weight)); } } } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h index 2c69e66..9172498 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h +++ b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h @@ -1,92 +1,88 @@ // ----------------------------------------------------------------------- // 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_GENERATE_MASK_VOXEL_LIST_H #define __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include #include #include namespace rttb { namespace masks { namespace boostRedesign { /*! @class BoostMaskGenerateMaskVoxelListThread * */ class BoostMaskGenerateMaskVoxelListThread { public: typedef ::boost::shared_ptr GeometricInfoPointer; typedef ::boost::multi_array BoostArray2D; typedef std::map BoostArrayMap; - typedef ::boost::shared_ptr MaskVoxelPointer; typedef ::boost::shared_ptr<::boost::lockfree::queue> MaskVoxelQueuePointer; + typedef std::vector VoxelIndexVector; - BoostMaskGenerateMaskVoxelListThread(int aGlobalBoundingBoxSize0, - int aGlobalBoundingBoxSize1, - const rttb::VoxelGridIndex3D& aMinIndex, - const rttb::VoxelGridIndex3D& aMaxIndex, + BoostMaskGenerateMaskVoxelListThread(const VoxelIndexVector& aGlobalBoundingBox, GeometricInfoPointer aGeometricInfo, const BoostArrayMap& aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, MaskVoxelQueuePointer aResultMaskVoxelQueue); void operator()(); private: - int _globalBoundingBoxSize0; - int _globalBoundingBoxSize1; - rttb::VoxelGridIndex3D _minIndex; - rttb::VoxelGridIndex3D _maxIndex; + VoxelIndexVector _globalBoundingBox; GeometricInfoPointer _geometricInfo; BoostArrayMap _voxelizationMap; //(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2) double _voxelizationThickness; unsigned int _beginSlice; + /*! @brief _endSlice is the first index not to be processed (like a end iterator) + */ unsigned int _endSlice; MaskVoxelQueuePointer _resultMaskVoxelQueue; /*! @brief For each dose grid index z, calculate the weight vector for each structure contour */ void calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const; }; } } } #endif \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskRedesign.cpp b/code/masks/boost/rttbBoostMaskRedesign.cpp index c4a1477..068d6d4 100644 --- a/code/masks/boost/rttbBoostMaskRedesign.cpp +++ b/code/masks/boost/rttbBoostMaskRedesign.cpp @@ -1,738 +1,782 @@ #include #include #include #include #include #include #include #include #include #include #include "rttbBoostMaskRedesign.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbBoostMaskGenerateMaskVoxelListThread.h" +#include "rttbBoostMaskVoxelizationThread.h" namespace rttb { namespace masks { namespace boostRedesign { + BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, - BoostMask::StructPointer aStructure, bool strict) - : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), + BoostMask::StructPointer aStructure, unsigned int threadSize, bool strict) + : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _threadSize(threadSize), _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() { + _threadSize = 8; preprocessing(); voxelization(); generateMaskVoxelList(); _isUpToDate = true; } void BoostMask::preprocessing() { rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); //Convert world coordinate polygons to the polygons with geometry coordinate rttb::PolygonSequenceType geometryCoordinatePolygonVector; rttb::PolygonSequenceType::iterator it; rttb::DoubleVoxelGridIndex3D globalMaxGridIndex(std::numeric_limits::lowest(), std::numeric_limits::lowest(), std::numeric_limits::lowest()); rttb::DoubleVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(), _geometricInfo->getNumRows(), 0); for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) { PolygonType rttbPolygon = *it; PolygonType geometryCoordinatePolygon; //1. convert polygon to geometry coordinate polygons //2. calculate global min/max //3. check if polygon is planar if (!preprocessingPolygon(rttbPolygon, geometryCoordinatePolygon, globalMinGridIndex, globalMaxGridIndex, errorConstant)) { throw rttb::core::Exception("TiltedMaskPlaneException"); } geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon); } rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) + 0.5), GridIndexType(globalMinGridIndex(1) + 0.5), GridIndexType(globalMinGridIndex(2) + 0.5)); rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) + 0.5), GridIndexType(globalMaxGridIndex(1) + 0.5), GridIndexType(globalMaxGridIndex(2) + 0.5)); _globalBoundingBox.push_back(minIndex); _globalBoundingBox.push_back(maxIndex); //convert rttb polygon sequence to a map of z index and a vector of boost ring 2d (without holes) - 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; + _ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); - for (itMap = ringMap.begin(); itMap != ringMap.end(); ++itMap) - { - BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); - _geometryCoordinateBoostPolygonMap.insert(std::pair((*itMap).first, - polygonVector)); - } } void BoostMask::voxelization() { BoostPolygonMap::iterator it; if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } - rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); + unsigned int ringMapSize = _ringMap.size(); + + BoostRingMap::iterator itMap; + + unsigned int mapSizeInAThread = _ringMap.size() / _threadSize; + unsigned int count = 0; + unsigned int countThread = 0; + BoostPolygonMap polygonMap; + std::vector polygonMapVector; + + //check donut and convert to a map of z index and a vector of boost polygon 2d (with or without holes) + for (itMap = _ringMap.begin(); itMap != _ringMap.end(); ++itMap) + { + BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); + + if (count == mapSizeInAThread && countThread < (_threadSize - 1)) + { + polygonMapVector.push_back(polygonMap); + polygonMap.clear(); + count = 0; + countThread++; + } + + polygonMap.insert(std::pair((*itMap).first, + polygonVector)); + count++; + + } + + polygonMapVector.push_back(polygonMap); //insert the last one + + //generate voxelization map, multi-threading + ::boost::thread_group threads; + BoostMaskVoxelizationThread::VoxelizationQueuePointer resultQueue + = ::boost::make_shared<::boost::lockfree::queue> + (ringMapSize); + BoostMaskVoxelizationThread::VoxelizationIndexQueuePointer resultIndexQueue + = ::boost::make_shared<::boost::lockfree::queue> + (ringMapSize); + + for (int i = 0; i < polygonMapVector.size(); ++i) + { + BoostMaskVoxelizationThread t(polygonMapVector.at(i), _globalBoundingBox, resultIndexQueue, + resultQueue); + threads.create_thread(t); + } + + threads.join_all(); + + + BoostArray2D* array; + double index; + + //generate voxelization map using resultQueue and resultIndexQueue + while (resultIndexQueue->pop(index) && resultQueue->pop(array)) + { + _voxelizationMap.insert(std::pair(index, *array)); + } + + + + /*rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; for (it = _geometryCoordinateBoostPolygonMap.begin(); it != _geometryCoordinateBoostPolygonMap.end(); ++it) { BoostArray2D maskArray(boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); BoostPolygonVector boostPolygonVec = (*it).second; for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) { for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = x + minIndex[0]; currentIndex[1] = y + minIndex[1]; currentIndex[2] = 0; //Get intersection polygons of the dose voxel and the structure BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec); //Calc areas of all intersection polygons double volumeFraction = calcArea(polygons); if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } maskArray[x][y] = volumeFraction; } } //insert into voxelization map _voxelizationMap.insert(std::pair((*it).first, maskArray)); - } + }*/ } void BoostMask::generateMaskVoxelList() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } //check homogeneus of the voxelization plane (the contours plane) if (!calcVoxelizationThickness(_voxelizationThickness)) { throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneus!"); } - rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); - rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); - int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; - int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; ::boost::shared_ptr<::boost::lockfree::queue> resultQueue = ::boost::make_shared<::boost::lockfree::queue> (_geometricInfo->getNumberOfVoxels()); - //BoostMaskGenerateMaskVoxelListThread::MaskVoxelQueuePointer resultQueue; - unsigned int thread_size = 8; + ::boost::thread_group threads; - unsigned int sliceNumberInAThread = _geometricInfo->getNumSlices() / thread_size ; - std::cout << "number of slices: " << _geometricInfo->getNumSlices() << - ", slice number in a thread: " << - sliceNumberInAThread << std::endl; + unsigned int sliceNumberInAThread = _geometricInfo->getNumSlices() / _threadSize; - for (unsigned int i = 0; i < thread_size; ++i) + for (unsigned int i = 0; i < _threadSize; ++i) { unsigned int beginSlice = i * sliceNumberInAThread; unsigned int endSlice; - if (i < thread_size - 1) + if (i < _threadSize - 1) { endSlice = (i + 1) * sliceNumberInAThread; } else { endSlice = _geometricInfo->getNumSlices(); } - std::cout << "BoostMask create slice" << std::endl; - BoostMaskGenerateMaskVoxelListThread t(globalBoundingBoxSize0, globalBoundingBoxSize1, - minIndex, maxIndex, _geometricInfo, _voxelizationMap, _voxelizationThickness, beginSlice, endSlice, + BoostMaskGenerateMaskVoxelListThread t(_globalBoundingBox, _geometricInfo, _voxelizationMap, + _voxelizationThickness, beginSlice, endSlice, resultQueue); threads.create_thread(t); } threads.join_all(); - std::cout << "BoostMask: all threads finished." << std::endl; - core::MaskVoxel* voxel; + //generate _voxelInStructure using resultQueue while (resultQueue->pop(voxel)) { _voxelInStructure->push_back(*voxel);//push back the mask voxel in structure } //calculate weight vector /* for (unsigned int indexZ = 0; indexZ < thread_size; ++indexZ) { std::map weightVectorForZ; calcWeightVector(indexZ, weightVectorForZ); //For each x,y, calc sum of all voxelization plane, use weight vector for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) { for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = x + minIndex[0]; currentIndex[1] = y + minIndex[1]; currentIndex[2] = indexZ; rttb::VoxelGridID gridID; _geometricInfo->convert(currentIndex, gridID); double volumeFraction = 0; std::map::iterator it; for (it = weightVectorForZ.begin(); it != weightVectorForZ.end(); ++it) { BoostArrayMap::iterator findVoxelizationIt = _voxelizationMap.find((*it).first); double weight = (*it).second; if (findVoxelizationIt == _voxelizationMap.end()) { throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneus!"); } BoostArray2D voxelizationArray = (*findVoxelizationIt).second; //calc sum of all voxelization plane, use weight volumeFraction += voxelizationArray[x][y] * weight; } if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } //insert mask voxel if volumeFraction > 0 if (volumeFraction > 0) { core::MaskVoxel maskVoxel = core::MaskVoxel(gridID, volumeFraction); _voxelInStructure->push_back(maskVoxel);//push back the mask voxel in structure } } } } */ } bool BoostMask::preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const { double minZ = _geometricInfo->getNumSlices(); double maxZ = 0.0; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D worldCoordinatePoint = aRTTBPolygon.at(i); //convert to geometry coordinate polygon rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint; _geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, geometryCoordinatePoint); geometryCoordinatePolygon.push_back(geometryCoordinatePoint); //calculate the current global min/max //min and max for x if (geometryCoordinatePoint(0) < minimum(0)) { minimum(0) = geometryCoordinatePoint(0); } if (geometryCoordinatePoint(0) > maximum(0)) { maximum(0) = geometryCoordinatePoint(0); } //min and max for y if (geometryCoordinatePoint(1) < minimum(1)) { minimum(1) = geometryCoordinatePoint(1); } if (geometryCoordinatePoint(1) > maximum(1)) { maximum(1) = geometryCoordinatePoint(1); } //min and max for z if (geometryCoordinatePoint(2) < minimum(2)) { minimum(2) = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maximum(2)) { maximum(2) = geometryCoordinatePoint(2); } //check planar if (geometryCoordinatePoint(2) < minZ) { minZ = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maxZ) { maxZ = geometryCoordinatePoint(2); } } return (abs(maxZ - minZ) <= aErrorConstant); } BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const { BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); if (i == 0) { firstPoint = boostPoint; } ::boost::geometry::append(polygon2D, boostPoint); } ::boost::geometry::append(polygon2D, firstPoint); return polygon2D; } BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( const rttb::PolygonSequenceType& aRTTBPolygonVector) { rttb::PolygonSequenceType::const_iterator it; BoostMask::BoostRingMap aRingMap; for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it) { rttb::PolygonType rttbPolygon = *it; double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon bool isFirstZ = true; if (!aRingMap.empty()) { BoostMask::BoostRingMap::iterator findIt = findNearestKey(aRingMap, zIndex, errorConstant); //if the z index is found (same slice), add the polygon to vector if (findIt != aRingMap.end()) { //BoostRingVector ringVector = ; (*findIt).second.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); isFirstZ = false; } } //if it is the first z index in the map, insert vector with the polygon if (isFirstZ) { BoostRingVector ringVector; ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); aRingMap.insert(std::pair(zIndex, ringVector)); } } return aRingMap; } BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap& aBoostRingMap, double aIndex, double aErrorConstant) { BoostMask::BoostRingMap::iterator find = aBoostRingMap.find(aIndex); //if find a key equivalent to aIndex, found if (find != aBoostRingMap.end()) { return find; } else { BoostMask::BoostRingMap::iterator lowerBound = aBoostRingMap.lower_bound(aIndex); //if all keys go before aIndex, check the last key if (lowerBound == aBoostRingMap.end()) { lowerBound = --aBoostRingMap.end(); } //if the lower bound very close to aIndex, found if (abs((*lowerBound).first - aIndex) <= aErrorConstant) { return lowerBound; } else { //if the lower bound is the beginning, not found if (lowerBound == aBoostRingMap.begin()) { return aBoostRingMap.end(); } else { BoostMask::BoostRingMap::iterator lowerBound1 = --lowerBound;//the key before the lower bound //if the key before the lower bound very close to a Index, found if (abs((*lowerBound1).first - aIndex) <= aErrorConstant) { return lowerBound1; } //else, not found else { return aBoostRingMap.end(); } } } } } BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector) const { //check donut BoostMask::BoostRingVector::const_iterator it1; BoostMask::BoostRingVector::const_iterator it2; BoostMask::BoostPolygonVector boostPolygonVector; std::vector donutIndexVector;//store the outer and inner ring index BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon //Get donut index and donut polygon unsigned int index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i = 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; } BoostMask::BoostRing2D BoostMask::get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) const { BoostMask::BoostRing2D polygon; BoostPoint2D point1(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] - 0.5); ::boost::geometry::append(polygon, point1); BoostPoint2D point2(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] - 0.5); ::boost::geometry::append(polygon, point2); BoostPoint2D point3(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] + 0.5); ::boost::geometry::append(polygon, point3); BoostPoint2D point4(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] + 0.5); ::boost::geometry::append(polygon, point4); ::boost::geometry::append(polygon, point1); return polygon; } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) const { 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) const { double area = 0; BoostPolygonDeque::const_iterator it; for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) { area += ::boost::geometry::area(*it); } return area; } bool BoostMask::calcVoxelizationThickness(double& aThickness) const { BoostArrayMap::const_iterator it = _voxelizationMap.begin(); BoostArrayMap::const_iterator it2 = ++_voxelizationMap.begin(); double lastContourIndexZ; if (_voxelizationMap.size() <= 1) { aThickness = 1; return true; } double thickness = 0; for (; it != _voxelizationMap.end(), it2 != _voxelizationMap.end(); ++it, ++it2) { if (thickness == 0) { thickness = (*it2).first - (*it).first; } else { //if no homogeneous, return false if (thickness != ((*it2).first - (*it).first)) { return false; } } } if (thickness != 0) { aThickness = thickness; } else { aThickness = 1; } return true; } void BoostMask::calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const { BoostArrayMap::const_iterator it; double indexZMin = aIndexZ - 0.5; double indexZMax = aIndexZ + 0.5; for (it = _voxelizationMap.begin(); it != _voxelizationMap.end(); ++it) { double voxelizationPlaneIndexMin = (*it).first - 0.5 * _voxelizationThickness; double voxelizationPlaneIndexMax = (*it).first + 0.5 * _voxelizationThickness; double weight = 0; if ((voxelizationPlaneIndexMin < indexZMin) && (voxelizationPlaneIndexMax > indexZMin)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = voxelizationPlaneIndexMax - indexZMin; } else { weight = 1; } } else if ((voxelizationPlaneIndexMin >= indexZMin) && (voxelizationPlaneIndexMin < indexZMax)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = _voxelizationThickness; } else { weight = indexZMax - voxelizationPlaneIndexMin; } } weightVector.insert(std::pair((*it).first, weight)); } } } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskRedesign.h b/code/masks/boost/rttbBoostMaskRedesign.h index 633451c..cd670f7 100644 --- a/code/masks/boost/rttbBoostMaskRedesign.h +++ b/code/masks/boost/rttbBoostMaskRedesign.h @@ -1,224 +1,234 @@ // ----------------------------------------------------------------------- // 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_R_H #define __BOOST_MASK_R_H #include "rttbBaseType.h" #include "rttbStructure.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include #include #include #include #include #include namespace rttb { namespace masks { namespace boostRedesign { /*! @class BoostMask * @brief Implementation of voxelization using boost::geometry. * @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); + BoostMask(GeometricInfoPointer aDoseGeoInfo, StructPointer aStructure, unsigned int threadSize = 8, + 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 typedef ::boost::multi_array BoostArray2D; typedef std::map BoostArrayMap; typedef ::boost::array BoostArray1D; typedef std::vector> WeightVector; GeometricInfoPointer _geometricInfo; StructPointer _structure; + /*! @brief A map of z index and a vector of boost ring 2d (without holes) + */ + BoostRingMap _ringMap; + /*! @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; + //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. */ VoxelIndexVector _globalBoundingBox; /*! @brief The voxelization map * key: the converted double z grid index of a contour plane * value: the 2d mask, array[i][j] = the mask value of the position (i,j) in the global bounding box, * i: 0 - (_globalBoundingBoxSize0-1), j: 0 - (_globalBoundingBoxSize1-1) */ BoostArrayMap _voxelizationMap; //The thickness of the voxelization plane (the contour plane), in double dose grid index //(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2) double _voxelizationThickness; bool _strict; //vector of the MaskVoxel inside the structure MaskVoxelListPointer _voxelInStructure; /*! @brief If the mask is up to date */ bool _isUpToDate; + + /*! @brief The thread size + */ + unsigned int _threadSize; + /*! @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 The voxelization step, wich computes the voxelization planes (in x/y) for all contours of an struct. * For each contour (that is in the z-Range of the reference geometry) of the struct: * 1) Allocate result array (voxelization plane) based on the bounding box (see Preprocessing Step 3) * 2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry). * 3) Add result Array (key is the z-Value of the contour) */ void voxelization(); /*! @final mask voxel Generation step which transfers the voxilization planes into the (z-)geometry of the reference geometry. * It consists of following Sub steps : * For all "slices" in the reference geometry : * 1) generate weight vector for all voxelization planes for a given z - value of a slice * Itterate over the bounding box of a struct.For each voxel : * 2) Compute weighted sum of all voxelization planes(use weight vector, step 1) * 2a) If sum > 0 : Add mask voxel for the current x / y(inner Loop) and z value(outer Loop). * 3) return mask voxel list. */ void generateMaskVoxelList(); /*! @brief Convert the rttb polygon with world corrdinate to the rttb polygon with double geometry coordinate, calculate the current min/max * and check if the polygon is planar * @param minimum the current global minimum * @param maximum the current global maximum * @return Return true if the polygon is planar, which means that the minminal and maximal z-coordinate of the polygon is not larger than a error constant */ bool preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const; /*! @brief Convert a rttb 3d polygon to a 2d boost ring*/ BoostRing2D convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const; /*! @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 Find the key with error constant to aIndex * @pre aBoostRingMap should not be empty * @return Return aBoostRingMap.end() if the key is not found */ BoostMask::BoostRingMap::iterator findNearestKey(BoostMask::BoostRingMap& aBoostRingMap, double aIndex, double aErrorConstant); /*! @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) const; /*! @brief Get the voxel 2d contour polygon in geometry coordinate*/ BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) const; /*! @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) const; /*! @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) const; /*! @brief Calculate the voxelization thickness. Return false, if the voxelization plane is not homogeneous */ bool calcVoxelizationThickness(double& aThickness) const; /*! @brief For each dose grid index z, calculate the weight vector for each structure contour */ void calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const; }; } } } #endif \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp b/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp index e69de29..592ba31 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp +++ b/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp @@ -0,0 +1,138 @@ +#include "rttbBoostMaskVoxelizationThread.h" + +namespace rttb +{ + namespace masks + { + namespace boostRedesign + { + BoostMaskVoxelizationThread::BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, + const VoxelIndexVector& aGlobalBoundingBox, VoxelizationIndexQueuePointer aResultIndexQueue, + VoxelizationQueuePointer aVoxelizationQueue) : _geometryCoordinateBoostPolygonMap(APolygonMap), + _globalBoundingBox(aGlobalBoundingBox), _resultIndexQueue(aResultIndexQueue), + _resultVoxelizationQueue(aVoxelizationQueue) + { + } + + + void BoostMaskVoxelizationThread::operator()() + { + + rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); + rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); + int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; + int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; + + BoostPolygonMap::iterator it; + + for (it = _geometryCoordinateBoostPolygonMap.begin(); + it != _geometryCoordinateBoostPolygonMap.end(); ++it) + { + BoostArray2D* maskArray = new BoostArray2D( + boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); + + BoostPolygonVector boostPolygonVec = (*it).second; + + for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) + { + for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) + { + rttb::VoxelGridIndex3D currentIndex; + currentIndex[0] = x + minIndex[0]; + currentIndex[1] = y + minIndex[1]; + currentIndex[2] = 0; + + //Get intersection polygons of the dose voxel and the structure + BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec); + + //Calc areas of all intersection polygons + double volumeFraction = calcArea(polygons); + + if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) + { + volumeFraction = 1; + } + else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) + { + throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); + } + + (*maskArray)[x][y] = volumeFraction; + } + } + + //insert into voxelization map + _resultIndexQueue->push((*it).first); + _resultVoxelizationQueue->push(maskArray); + + } + + + } + + /*Get intersection polygons of the contour and a voxel polygon*/ + BoostMaskVoxelizationThread::BoostPolygonDeque BoostMaskVoxelizationThread::getIntersections( + const rttb::VoxelGridIndex3D& + aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) const + { + BoostMaskVoxelizationThread::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; + } + + BoostMaskVoxelizationThread::BoostRing2D BoostMaskVoxelizationThread::get2DContour( + const rttb::VoxelGridIndex3D& aVoxelGrid3D) const + { + BoostRing2D polygon; + + + BoostPoint2D point1(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] - 0.5); + ::boost::geometry::append(polygon, point1); + + BoostPoint2D point2(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] - 0.5); + ::boost::geometry::append(polygon, point2); + + BoostPoint2D point3(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] + 0.5); + ::boost::geometry::append(polygon, point3); + + BoostPoint2D point4(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] + 0.5); + ::boost::geometry::append(polygon, point4); + + ::boost::geometry::append(polygon, point1); + + return polygon; + + } + + /*Calculate the intersection area*/ + double BoostMaskVoxelizationThread::calcArea(const BoostPolygonDeque& aPolygonDeque) const + { + double area = 0; + + BoostPolygonDeque::const_iterator it; + + for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) + { + area += ::boost::geometry::area(*it); + } + + return area; + } + } + } +} \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskVoxelizationThread.h b/code/masks/boost/rttbBoostMaskVoxelizationThread.h index a79ac18..a9c20e5 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.h +++ b/code/masks/boost/rttbBoostMaskVoxelizationThread.h @@ -1,111 +1,114 @@ // ----------------------------------------------------------------------- // 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_VOXELIZATION_THREAD_H #define __BOOST_MASK_VOXELIZATION_THREAD_H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" +#include "rttbInvalidParameterException.h" #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include namespace rttb { namespace masks { namespace boostRedesign { /*! @class BoostMaskGenerateMaskVoxelListThread * */ class BoostMaskVoxelizationThread { public: + typedef ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy > + BoostPolygon2D; + typedef std::vector BoostPolygonVector;//polygon with or without holes + typedef std::map + BoostPolygonMap;//map of the z index with the vector of boost 2d polygon + typedef std::vector VoxelIndexVector; + typedef ::boost::multi_array BoostArray2D; + typedef ::boost::shared_ptr<::boost::lockfree::queue> + VoxelizationQueuePointer; + typedef ::boost::shared_ptr<::boost::lockfree::queue> + VoxelizationIndexQueuePointer; + BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, + const VoxelIndexVector& aGlobalBoundingBox, VoxelizationIndexQueuePointer aResultIndexQueue, + VoxelizationQueuePointer aVoxelizationQueue); - BoostMaskVoxelizationThread(); - - - void operator()() - { - - /*for (it = _geometryCoordinateBoostPolygonMap.begin(); - it != _geometryCoordinateBoostPolygonMap.end(); ++it) - { - BoostArray2D maskArray(boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); - - BoostPolygonVector boostPolygonVec = (*it).second; - - for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) - { - for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) - { - rttb::VoxelGridIndex3D currentIndex; - currentIndex[0] = x + minIndex[0]; - currentIndex[1] = y + minIndex[1]; - currentIndex[2] = 0; - //Get intersection polygons of the dose voxel and the structure - BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec); + void operator()(); - //Calc areas of all intersection polygons - double volumeFraction = calcArea(polygons); - if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) - { - volumeFraction = 1; - } - else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) - { - throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); - } + private: - maskArray[x][y] = volumeFraction; - } - } - //insert into voxelization map - _voxelizationMap.insert(std::pair((*it).first, maskArray)); - }*/ - } + typedef std::deque BoostPolygonDeque; + typedef ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy > + BoostRing2D; + typedef ::boost::geometry::model::d2::point_xy BoostPoint2D; - private: - typedef ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy > - BoostPolygon2D; - typedef std::vector BoostPolygonVector;//polygon with or without holes - typedef std::map - BoostPolygonMap;//map of the z index with the vector of boost 2d polygon BoostPolygonMap _geometryCoordinateBoostPolygonMap; + VoxelIndexVector _globalBoundingBox; + VoxelizationQueuePointer _resultVoxelizationQueue; + VoxelizationIndexQueuePointer _resultIndexQueue; + + /*! @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) const; + + /*! @brief Get the voxel 2d contour polygon in geometry coordinate*/ + BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) const; + + /*! @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) const; }; } } } #endif \ No newline at end of file