diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp index 9d6786d..3c2bbed 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/boost/rttbBoostMask.cpp @@ -1,561 +1,564 @@ // ----------------------------------------------------------------------- // 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) */ #include #include #include #include #include #include "rttbBoostMask.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbBoostMaskVoxelizationThread.h" namespace rttb { namespace masks { namespace boost { BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, BoostMask::StructPointer aStructure, bool strict, unsigned int numberOfThreads) : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _numberOfThreads(numberOfThreads), _voxelizationThickness(0.0), _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!"); } if (_numberOfThreads == 0) { _numberOfThreads = ::boost::thread::hardware_concurrency(); if (_numberOfThreads == 0) { throw rttb::core::InvalidParameterException("Error: detection of the number of hardware threads is not possible. Please specify number of threads for voxelization explicitly as parameter in BoostMask."); } } } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { if (!_isUpToDate) { calcMask(); } return _voxelInStructure; } void BoostMask::calcMask() { 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::min(), std::numeric_limits::min(), std::numeric_limits::min()); 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) _ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); } void BoostMask::voxelization() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } BoostRingMap::iterator itMap; size_t mapSizeInAThread = _ringMap.size() / _numberOfThreads; 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) { //the vector of all boost 2d polygons with the same z grid index(donut polygon is accepted). BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); if (count == mapSizeInAThread && countThread < (_numberOfThreads - 1)) { polygonMapVector.push_back(polygonMap); polygonMap.clear(); count = 0; countThread++; } polygonMap.insert(std::pair((*itMap).first, polygonVector)); count++; } _voxelizationMap = ::boost::make_shared >(); polygonMapVector.push_back(polygonMap); //insert the last one //generate voxelization map, multi-threading ::boost::thread_group threads; auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); for (int i = 0; i < polygonMapVector.size(); ++i) { BoostMaskVoxelizationThread t(polygonMapVector.at(i), _globalBoundingBox, _voxelizationMap, aMutex); threads.create_thread(t); } threads.join_all(); } void BoostMask::generateMaskVoxelList() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } //check homogeneous of the voxelization plane (the contours plane) if (!calcVoxelizationThickness(_voxelizationThickness)) { throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneous!"); } + + ::boost::thread_group threads; auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); unsigned int sliceNumberInAThread = _geometricInfo->getNumSlices() / _numberOfThreads; //generate mask voxel list, multi-threading for (unsigned int i = 0; i < _numberOfThreads; ++i) { unsigned int beginSlice = i * sliceNumberInAThread; unsigned int endSlice; if (i < _numberOfThreads - 1) { endSlice = (i + 1) * sliceNumberInAThread; } else { endSlice = _geometricInfo->getNumSlices(); } BoostMaskGenerateMaskVoxelListThread t(_globalBoundingBox, _geometricInfo, _voxelizationMap, _voxelizationThickness, beginSlice, endSlice, _voxelInStructure, aMutex); threads.create_thread(t); } threads.join_all(); } 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); + return (std::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) const { 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) const { 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) + if (std::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) + if (std::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; } bool BoostMask::calcVoxelizationThickness(double& aThickness) const { if (_voxelizationMap->size() <= 1) { aThickness = 1; return true; } double thickness = 0; auto it = _voxelizationMap->cbegin(); auto it2 = ++_voxelizationMap->cbegin(); for (; it != _voxelizationMap->cend(), it2 != _voxelizationMap->cend(); ++it, ++it2) { if (thickness == 0) { thickness = it2->first - it->first; } else { double curThickness = it2->first - it->first; //if not homogeneous (leave out double imprecisions), return false - if (abs(thickness-curThickness)>errorConstant) + if (std::abs(thickness-curThickness)>errorConstant) { return false; } } } if (thickness != 0) { aThickness = thickness; } else { aThickness = 1; } return true; } } } } diff --git a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp index fa25303..98bfd95 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp +++ b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp @@ -1,155 +1,156 @@ // ----------------------------------------------------------------------- // 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 "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbInvalidParameterException.h" #include namespace rttb { namespace masks { namespace boost { BoostMaskGenerateMaskVoxelListThread::BoostMaskGenerateMaskVoxelListThread( const VoxelIndexVector& aGlobalBoundingBox, GeometricInfoPointer aGeometricInfo, BoostArrayMapPointer aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, MaskVoxelListPointer aMaskVoxelList, ::boost::shared_ptr<::boost::shared_mutex> aMutex) : _globalBoundingBox(aGlobalBoundingBox), _geometricInfo(aGeometricInfo), _voxelizationMap(aVoxelizationMap), _voxelizationThickness(aVoxelizationThickness), _beginSlice(aBeginSlice), _endSlice(aEndSlice), _resultMaskVoxelList(aMaskVoxelList), _mutex(aMutex) {} void BoostMaskGenerateMaskVoxelListThread::operator()() { rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); unsigned int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; unsigned int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; std::vector maskVoxelsInThread; for (unsigned int indexZ = _beginSlice; indexZ < _endSlice; ++indexZ) { //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 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; auto it = weightVectorForZ.cbegin(); auto itMap = _voxelizationMap->cbegin(); for (; it != weightVectorForZ.cend() && itMap != _voxelizationMap->cend(); ++it, ++itMap) { double weight = it->second; if (weight > 0){ BoostArray2DPointer 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 = core::MaskVoxel(gridID, volumeFraction); maskVoxelsInThread.push_back(maskVoxelPtr); } } } } ::boost::unique_lock<::boost::shared_mutex> lock(*_mutex); _resultMaskVoxelList->insert(_resultMaskVoxelList->end(), maskVoxelsInThread.begin(), maskVoxelsInThread.end()); } void BoostMaskGenerateMaskVoxelListThread::calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const { double indexZMin = aIndexZ - 0.5; double indexZMax = aIndexZ + 0.5; for (auto 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)); } } } } } diff --git a/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp b/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp index 8aa32fd..4ff1d8f 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp +++ b/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp @@ -1,159 +1,159 @@ // ----------------------------------------------------------------------- // 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 "rttbBoostMaskVoxelizationThread.h" #include "rttbInvalidParameterException.h" #include #include #include namespace rttb { namespace masks { namespace boost { BoostMaskVoxelizationThread::BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex) : _geometryCoordinateBoostPolygonMap(APolygonMap), _globalBoundingBox(aGlobalBoundingBox), _resultVoxelization(anArrayMap), _mutex(aMutex) { } void BoostMaskVoxelizationThread::operator()() { rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); const unsigned int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; const unsigned int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; std::map > voxelizationMapInThread; for (BoostPolygonMap::iterator 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; } } voxelizationMapInThread.insert(std::pair(it->first, ::boost::make_shared(maskArray))); } //insert gathered values into voxelization map ::boost::unique_lock<::boost::shared_mutex> lock(*_mutex); _resultVoxelization->insert(voxelizationMapInThread.begin(), voxelizationMapInThread.end()); } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMaskVoxelizationThread::BoostPolygonDeque BoostMaskVoxelizationThread::getIntersections( const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) { 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) { 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) { 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/testing/core/CreateTestStructure.cpp b/testing/core/CreateTestStructure.cpp index 98e82a5..43db0a4 100644 --- a/testing/core/CreateTestStructure.cpp +++ b/testing/core/CreateTestStructure.cpp @@ -1,691 +1,692 @@ // ----------------------------------------------------------------------- // 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 rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include "CreateTestStructure.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace testing { CreateTestStructure::~CreateTestStructure() {} CreateTestStructure::CreateTestStructure(const core::GeometricInfo& aGeoInfo) { _geoInfo = aGeoInfo; } PolygonType CreateTestStructure::createPolygonLeftUpper(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); polygon.push_back(p1); std::cout << "(" << p1.x() << "," << p1.y() << "," << p1.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonCenter(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() / 2; p(1) = p1.y() + delta.y() / 2; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonCenterOnPlaneCenter(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; DoubleVoxelGridIndex3D indexDouble = DoubleVoxelGridIndex3D((aVoxelVector.at(i)).x(), (aVoxelVector.at(i)).y(), sliceNumber); WorldCoordinate3D p1; _geoInfo.geometryCoordinateToWorldCoordinate(indexDouble, p1); polygon.push_back(p1); + std::cout << "(" << p1.x() << "," << p1.y() << "," << p1.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonBetweenUpperLeftAndCenter( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() / 4; p(1) = p1.y() + delta.y() / 4; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonBetweenLowerRightAndCenter( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() * 0.75; p(1) = p1.y() + delta.y() * 0.75; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonLeftEdgeMiddle(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x(); p(1) = p1.y() + delta.y() * 0.5; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonUpperCenter(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() * 0.5; p(1) = p1.y(); p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonIntermediatePoints(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; VoxelGridIndex3D voxelIndex2; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; if (i < (aVoxelVector.size() - 1)) { voxelIndex2(0) = (aVoxelVector.at(i + 1)).x(); voxelIndex2(1) = (aVoxelVector.at(i + 1)).y(); voxelIndex2(2) = sliceNumber; } else { voxelIndex2(0) = (aVoxelVector.at(0)).x(); voxelIndex2(1) = (aVoxelVector.at(0)).y(); voxelIndex2(2) = sliceNumber; } WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p2; _geoInfo.indexToWorldCoordinate(voxelIndex2, p2); SpacingVectorType3D delta2 = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x() + delta.x() * 0.75; wp1(1) = p1.y() + delta.y() * 0.75; wp1(2) = p1.z(); WorldCoordinate3D wp2; wp2(0) = p2.x() + delta.x() * 0.75; wp2(1) = p2.y() + delta.y() * 0.75; wp2(2) = p2.z(); polygon.push_back(wp1); double diffX = (wp2.x() - wp1.x()) / 1000.0; double diffY = (wp2.y() - wp1.y()) / 1000.0; WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; for (int i = 0 ; i < 1000 ; i++) { wp_intermediate(0) = wp1.x() + i * diffX; wp_intermediate(1) = wp1.y() + i * diffY; polygon.push_back(wp_intermediate); } } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonCircle(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); double radius = 2 * delta.x(); double frac_radius = (radius * 0.001); double correct_y = (delta.x() / delta.y()); for (int i = 0 ; i <= 1000 ; i++) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() + frac_radius * i; wp_intermediate(1) = wp1.y() - (y_offset * correct_y) ; polygon.push_back(wp_intermediate); } for (int i = 1000 ; i >= 0 ; i--) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))) ; wp_intermediate(0) = wp1.x() + frac_radius * i; wp_intermediate(1) = wp1.y() + (y_offset * correct_y); polygon.push_back(wp_intermediate); } for (int i = 0 ; i <= 1000 ; i++) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() - frac_radius * i; wp_intermediate(1) = wp1.y() + y_offset * correct_y ; polygon.push_back(wp_intermediate); } for (int i = 1000 ; i >= 0 ; i--) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() - frac_radius * i ; wp_intermediate(1) = wp1.y() - (y_offset * correct_y); polygon.push_back(wp_intermediate); } } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSeveralSectionsInsideOneVoxelA( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x() + (delta.x() * 0.25); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.25); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.75); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.75); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.25); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.25); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.5); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.75); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.75); wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSelfTouchingA(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y(); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSelfTouchingB(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y(); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } }//testing }//rttb