diff --git a/code/masks/CMakeLists.txt b/code/masks/CMakeLists.txt index 2ef7620..ea7fb93 100644 --- a/code/masks/CMakeLists.txt +++ b/code/masks/CMakeLists.txt @@ -1,5 +1,9 @@ -MESSAGE (STATUS "processing RTToolbox mask") - -RTTB_CREATE_MODULE(RTTBMasks DEPENDS RTTBCore PACKAGE_DEPENDS BoostBinaries) - -ADD_SUBDIRECTORY(boost) +MESSAGE (STATUS "processing RTToolbox mask") + +RTTB_CREATE_MODULE(RTTBMasks DEPENDS RTTBCore PACKAGE_DEPENDS BoostBinaries) + +SET(RTTB_Boost_ADDITIONAL_COMPONENT thread) +RTTB_CREATE_MODULE(RTTBBoostMask DEPENDS RTTBCore RTTBMasks PACKAGE_DEPENDS BoostBinaries) +IF (CMAKE_COMPILER_IS_GNUCC) + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fext-numeric-literals") +ENDIF() diff --git a/code/masks/boost/CMakeLists.txt b/code/masks/boost/CMakeLists.txt deleted file mode 100644 index 92bb848..0000000 --- a/code/masks/boost/CMakeLists.txt +++ /dev/null @@ -1,5 +0,0 @@ -SET(RTTB_Boost_ADDITIONAL_COMPONENT thread) -RTTB_CREATE_MODULE(RTTBBoostMask DEPENDS RTTBCore RTTBMasks PACKAGE_DEPENDS BoostBinaries) -IF (CMAKE_COMPILER_IS_GNUCC) - SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fext-numeric-literals") -ENDIF() diff --git a/code/masks/boost/files.cmake b/code/masks/boost/files.cmake deleted file mode 100644 index faa905d..0000000 --- a/code/masks/boost/files.cmake +++ /dev/null @@ -1,13 +0,0 @@ -SET(CPP_FILES - rttbBoostMask.cpp - rttbBoostMaskAccessor.cpp - rttbBoostMaskGenerateMaskVoxelListThread.cpp - rttbBoostMaskVoxelizationThread.cpp -) - -SET(H_FILES - rttbBoostMask.h - rttbBoostMaskAccessor.h - rttbBoostMaskGenerateMaskVoxelListThread.h - rttbBoostMaskVoxelizationThread.h -) diff --git a/code/masks/files.cmake b/code/masks/files.cmake index 0cc43df..a3fdf77 100644 --- a/code/masks/files.cmake +++ b/code/masks/files.cmake @@ -1,7 +1,15 @@ SET(CPP_FILES rttbGenericMutableMaskAccessor.cpp + rttbBoostMask.cpp + rttbBoostMaskAccessor.cpp + rttbBoostMaskGenerateMaskVoxelListThread.cpp + rttbBoostMaskVoxelizationThread.cpp ) SET(H_FILES rttbGenericMutableMaskAccessor.h + rttbBoostMask.h + rttbBoostMaskAccessor.h + rttbBoostMaskGenerateMaskVoxelListThread.h + rttbBoostMaskVoxelizationThread.h ) diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/rttbBoostMask.cpp similarity index 96% rename from code/masks/boost/rttbBoostMask.cpp rename to code/masks/rttbBoostMask.cpp index bc8cf17..af20a7d 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/rttbBoostMask.cpp @@ -1,565 +1,565 @@ -// ----------------------------------------------------------------------- -// 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(core::GeometricInfo::Pointer aDoseGeoInfo, - core::Structure::Pointer 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 == nullptr) - { - throw rttb::core::NullPointerException("Error: Geometric info is nullptr!"); - } - else if (_structure == nullptr) - { - throw rttb::core::NullPointerException("Error: Structure is nullptr!"); - } - - 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) ), - GridIndexType(globalMinGridIndex(1) ), GridIndexType(globalMinGridIndex(2) )); - rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) ), - GridIndexType(globalMaxGridIndex(1) ), GridIndexType(globalMaxGridIndex(2) )); - - _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 (const auto & i : polygonMapVector) - { - BoostMaskVoxelizationThread t(i, _globalBoundingBox, - _voxelizationMap, aMutex, _strict); - 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 (auto worldCoordinatePoint : aRTTBPolygon) - { - //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 (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()) - { - auto 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 - { - auto find = aBoostRingMap.find(aIndex); - - //if find a key equivalent to aIndex, found - if (find != aBoostRingMap.end()) - { - return find; - } - else - { - auto 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 (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 - { - auto lowerBound1 = --lowerBound;//the key before the lower bound - - //if the key before the lower bound very close to a Index, found - 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 : donutIndexVector) - { - if (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 : donutIndexVector) - { - if (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 (std::abs(thickness-curThickness)>errorConstant) - { - //return false; - std::cout << "Two polygons are far from each other?" << std::endl; - } - } - - } - - if (thickness != 0) - { - aThickness = thickness; - } - else - { - aThickness = 1; - } - - return true; - } - } - } -} +// ----------------------------------------------------------------------- +// 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(core::GeometricInfo::Pointer aDoseGeoInfo, + core::Structure::Pointer 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 == nullptr) + { + throw rttb::core::NullPointerException("Error: Geometric info is nullptr!"); + } + else if (_structure == nullptr) + { + throw rttb::core::NullPointerException("Error: Structure is nullptr!"); + } + + 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) ), + GridIndexType(globalMinGridIndex(1) ), GridIndexType(globalMinGridIndex(2) )); + rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) ), + GridIndexType(globalMaxGridIndex(1) ), GridIndexType(globalMaxGridIndex(2) )); + + _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 (const auto & i : polygonMapVector) + { + BoostMaskVoxelizationThread t(i, _globalBoundingBox, + _voxelizationMap, aMutex, _strict); + 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 (auto worldCoordinatePoint : aRTTBPolygon) + { + //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 (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()) + { + auto 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 + { + auto find = aBoostRingMap.find(aIndex); + + //if find a key equivalent to aIndex, found + if (find != aBoostRingMap.end()) + { + return find; + } + else + { + auto 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 (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 + { + auto lowerBound1 = --lowerBound;//the key before the lower bound + + //if the key before the lower bound very close to a Index, found + 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 : donutIndexVector) + { + if (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 : donutIndexVector) + { + if (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 (std::abs(thickness-curThickness)>errorConstant) + { + //return false; + std::cout << "Two polygons are far from each other?" << std::endl; + } + } + + } + + if (thickness != 0) + { + aThickness = thickness; + } + else + { + aThickness = 1; + } + + return true; + } + } + } +} diff --git a/code/masks/boost/rttbBoostMask.h b/code/masks/rttbBoostMask.h similarity index 97% rename from code/masks/boost/rttbBoostMask.h rename to code/masks/rttbBoostMask.h index 9481c17..cf994bb 100644 --- a/code/masks/boost/rttbBoostMask.h +++ b/code/masks/rttbBoostMask.h @@ -1,202 +1,202 @@ -// ----------------------------------------------------------------------- -// 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) -*/ - - -#ifndef __BOOST_MASK_R_H -#define __BOOST_MASK_R_H - -#include "rttbBaseType.h" -#include "rttbStructure.h" -#include "rttbGeometricInfo.h" -#include "rttbMaskAccessorInterface.h" - -#include -#include -#include -#include - -namespace rttb -{ - namespace masks - { - namespace boost - { - /*! @class BoostMask - * @brief Implementation of voxelization using boost::geometry. - * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. - * (A structure without self intersection 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: - using MaskVoxelList = core::MaskAccessorInterface::MaskVoxelList; - using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer; - - /*! @brief Constructor - * @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is nullptr - * @param aDoseGeoInfo the GeometricInfo - * @param aStructure the structure set - * @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. - * @param numberOfThreads number of threads used for voxelization. default value 0 means automatic detection, using the number of Hardware thread/cores - * @exception InvalidParameterException thrown if strict is true and the structure has self intersections - */ - BoostMask(core::GeometricInfo::Pointer aDoseGeoInfo, core::Structure::Pointer aStructure, - bool strict = true, unsigned int numberOfThreads = 0); - - /*! @brief Generate mask and return the voxels in the mask - * @exception rttb::core::InvalidParameterException thrown if the structure has self intersections - */ - MaskVoxelListPointer getRelevantVoxelVector(); - - private: - using BoostPoint2D = ::boost::geometry::model::d2::point_xy; - using BoostPolygon2D = ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy >; - using BoostRing2D = ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy >; - using BoostRingVector = std::vector;//polygon without holes - using BoostPolygonVector = std::vector;//polygon with or without holes - using VoxelIndexVector = std::vector; - 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; - using BoostArray2DPointer = ::boost::shared_ptr; - typedef ::boost::shared_ptr > BoostArrayMapPointer; - - core::GeometricInfo::Pointer _geometricInfo; - - core::Structure::Pointer _structure; - - bool _strict; - - /*! @brief The number of threads - */ - unsigned int _numberOfThreads; - - //@brief The thickness of the voxelization plane (the contour plane), in double dose grid index - //@details 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; - - //@brief vector of the MaskVoxel inside the structure - MaskVoxelListPointer _voxelInStructure; - - /*! @brief The map of z index and a vector of boost ring 2d (without holes) - * @details Key: the double z grid index - * Value: the vector of boost ring 2d (without holes) - */ - BoostRingMap _ringMap; - - /*! @brief The min and max index of the global bounding box. - * @details 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 - * @details 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) - */ - BoostArrayMapPointer _voxelizationMap; - - /*! @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: - * @details 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, which computes the voxelization planes (in x/y) for all contours of an struct. - - * @details 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(); - - /*! @brief mask voxel Generation step which transfers the voxelization planes into the (z-)geometry of the reference geometry. - * @details 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 - * Iterate 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 coordinate 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 minimal 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) const; - - /*! @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) const; - - /*! @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 Calculate the voxelization thickness. - Return false, if the voxelization plane is not homogeneous - */ - bool calcVoxelizationThickness(double& aThickness) const; - - }; - - } - - - } -} - -#endif +// ----------------------------------------------------------------------- +// 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) +*/ + + +#ifndef __BOOST_MASK_R_H +#define __BOOST_MASK_R_H + +#include "rttbBaseType.h" +#include "rttbStructure.h" +#include "rttbGeometricInfo.h" +#include "rttbMaskAccessorInterface.h" + +#include +#include +#include +#include + +namespace rttb +{ + namespace masks + { + namespace boost + { + /*! @class BoostMask + * @brief Implementation of voxelization using boost::geometry. + * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. + * (A structure without self intersection 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: + using MaskVoxelList = core::MaskAccessorInterface::MaskVoxelList; + using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer; + + /*! @brief Constructor + * @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is nullptr + * @param aDoseGeoInfo the GeometricInfo + * @param aStructure the structure set + * @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. + * @param numberOfThreads number of threads used for voxelization. default value 0 means automatic detection, using the number of Hardware thread/cores + * @exception InvalidParameterException thrown if strict is true and the structure has self intersections + */ + BoostMask(core::GeometricInfo::Pointer aDoseGeoInfo, core::Structure::Pointer aStructure, + bool strict = true, unsigned int numberOfThreads = 0); + + /*! @brief Generate mask and return the voxels in the mask + * @exception rttb::core::InvalidParameterException thrown if the structure has self intersections + */ + MaskVoxelListPointer getRelevantVoxelVector(); + + private: + using BoostPoint2D = ::boost::geometry::model::d2::point_xy; + using BoostPolygon2D = ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy >; + using BoostRing2D = ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy >; + using BoostRingVector = std::vector;//polygon without holes + using BoostPolygonVector = std::vector;//polygon with or without holes + using VoxelIndexVector = std::vector; + 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; + using BoostArray2DPointer = ::boost::shared_ptr; + typedef ::boost::shared_ptr > BoostArrayMapPointer; + + core::GeometricInfo::Pointer _geometricInfo; + + core::Structure::Pointer _structure; + + bool _strict; + + /*! @brief The number of threads + */ + unsigned int _numberOfThreads; + + //@brief The thickness of the voxelization plane (the contour plane), in double dose grid index + //@details 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; + + //@brief vector of the MaskVoxel inside the structure + MaskVoxelListPointer _voxelInStructure; + + /*! @brief The map of z index and a vector of boost ring 2d (without holes) + * @details Key: the double z grid index + * Value: the vector of boost ring 2d (without holes) + */ + BoostRingMap _ringMap; + + /*! @brief The min and max index of the global bounding box. + * @details 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 + * @details 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) + */ + BoostArrayMapPointer _voxelizationMap; + + /*! @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: + * @details 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, which computes the voxelization planes (in x/y) for all contours of an struct. + + * @details 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(); + + /*! @brief mask voxel Generation step which transfers the voxelization planes into the (z-)geometry of the reference geometry. + * @details 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 + * Iterate 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 coordinate 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 minimal 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) const; + + /*! @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) const; + + /*! @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 Calculate the voxelization thickness. + Return false, if the voxelization plane is not homogeneous + */ + bool calcVoxelizationThickness(double& aThickness) const; + + }; + + } + + + } +} + +#endif diff --git a/code/masks/boost/rttbBoostMaskAccessor.cpp b/code/masks/rttbBoostMaskAccessor.cpp similarity index 96% rename from code/masks/boost/rttbBoostMaskAccessor.cpp rename to code/masks/rttbBoostMaskAccessor.cpp index 61d2525..0985637 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.cpp +++ b/code/masks/rttbBoostMaskAccessor.cpp @@ -1,161 +1,161 @@ -// ----------------------------------------------------------------------- -// RTToolbox - DKFZ radiotherapy quantitative evaluation library -// -// Copyright (c) German Cancer Research Center (DKFZ), -// Software development for Integrated Diagnostics and Therapy (SIDT). -// ALL RIGHTS RESERVED. -// See rttbCopyright.txt or -// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html -// -// This software is distributed WITHOUT ANY WARRANTY; without even -// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -// PURPOSE. See the above copyright notices for more information. -// -//------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) -*/ - -#include "rttbBoostMaskAccessor.h" -#include "rttbBoostMask.h" - -#include - -#include -#include -#include - -namespace rttb -{ - namespace masks - { - namespace boost - { - - BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, - const core::GeometricInfo& aGeometricInfo, bool strict) - : _spStructure(aStructurePointer), _geoInfo(aGeometricInfo), _strict(strict) - { - _spRelevantVoxelVector = MaskVoxelListPointer(); - - //generate new structure set uid - ::boost::uuids::uuid id; - ::boost::uuids::random_generator generator; - id = generator(); - std::stringstream ss; - ss << id; - _maskUID = "BoostMask_" + ss.str(); - } - - - BoostMaskAccessor::~BoostMaskAccessor() - = default; - - void BoostMaskAccessor::updateMask() - { - MaskVoxelList newRelevantVoxelVector; - - if (_spRelevantVoxelVector) - { - return; // already calculated - } - - BoostMask mask(::boost::make_shared(_geoInfo), - _spStructure, _strict); - - _spRelevantVoxelVector = mask.getRelevantVoxelVector(); - } - - BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector() - { - // if not already generated start voxelization here - updateMask(); - return _spRelevantVoxelVector; - } - - BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector( - float lowerThreshold) - { - auto filteredVoxelVectorPointer = ::boost::make_shared(); - updateMask(); - // filter relevant voxels - auto it = _spRelevantVoxelVector->begin(); - - while (it != _spRelevantVoxelVector->end()) - { - if ((*it).getRelevantVolumeFraction() > lowerThreshold) - { - filteredVoxelVectorPointer->push_back(*it); - } - - ++it; - } - - // if mask calculation was not successful this is empty! - return filteredVoxelVectorPointer; - } - - bool BoostMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const - { - //initialize return voxel - voxel.setRelevantVolumeFraction(0); - - //check if ID is valid - if (!_geoInfo.validID(aID)) - { - return false; - } - - //determine how a given voxel on the dose grid is masked - if (_spRelevantVoxelVector) - { - auto it = _spRelevantVoxelVector->begin(); - - while (it != _spRelevantVoxelVector->end()) - { - if ((*it).getVoxelGridID() == aID) - { - voxel = (*it); - return true; - } - - ++it; - } - - } - // returns false if mask was not calculated without triggering calculation (otherwise not const!) - else - { - return false; - } - - return false; - - } - - bool BoostMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const - { - //convert VoxelGridIndex3D to VoxelGridID - VoxelGridID aVoxelGridID; - - if (_geoInfo.convert(aIndex, aVoxelGridID)) - { - return getMaskAt(aVoxelGridID, voxel); - } - else - { - return false; - } - } - - const core::GeometricInfo& BoostMaskAccessor::getGeometricInfo() const - { - return _geoInfo; - }; - } - - } +// ----------------------------------------------------------------------- +// 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 "rttbBoostMaskAccessor.h" +#include "rttbBoostMask.h" + +#include + +#include +#include +#include + +namespace rttb +{ + namespace masks + { + namespace boost + { + + BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, + const core::GeometricInfo& aGeometricInfo, bool strict) + : _spStructure(aStructurePointer), _geoInfo(aGeometricInfo), _strict(strict) + { + _spRelevantVoxelVector = MaskVoxelListPointer(); + + //generate new structure set uid + ::boost::uuids::uuid id; + ::boost::uuids::random_generator generator; + id = generator(); + std::stringstream ss; + ss << id; + _maskUID = "BoostMask_" + ss.str(); + } + + + BoostMaskAccessor::~BoostMaskAccessor() + = default; + + void BoostMaskAccessor::updateMask() + { + MaskVoxelList newRelevantVoxelVector; + + if (_spRelevantVoxelVector) + { + return; // already calculated + } + + BoostMask mask(::boost::make_shared(_geoInfo), + _spStructure, _strict); + + _spRelevantVoxelVector = mask.getRelevantVoxelVector(); + } + + BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector() + { + // if not already generated start voxelization here + updateMask(); + return _spRelevantVoxelVector; + } + + BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector( + float lowerThreshold) + { + auto filteredVoxelVectorPointer = ::boost::make_shared(); + updateMask(); + // filter relevant voxels + auto it = _spRelevantVoxelVector->begin(); + + while (it != _spRelevantVoxelVector->end()) + { + if ((*it).getRelevantVolumeFraction() > lowerThreshold) + { + filteredVoxelVectorPointer->push_back(*it); + } + + ++it; + } + + // if mask calculation was not successful this is empty! + return filteredVoxelVectorPointer; + } + + bool BoostMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const + { + //initialize return voxel + voxel.setRelevantVolumeFraction(0); + + //check if ID is valid + if (!_geoInfo.validID(aID)) + { + return false; + } + + //determine how a given voxel on the dose grid is masked + if (_spRelevantVoxelVector) + { + auto it = _spRelevantVoxelVector->begin(); + + while (it != _spRelevantVoxelVector->end()) + { + if ((*it).getVoxelGridID() == aID) + { + voxel = (*it); + return true; + } + + ++it; + } + + } + // returns false if mask was not calculated without triggering calculation (otherwise not const!) + else + { + return false; + } + + return false; + + } + + bool BoostMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const + { + //convert VoxelGridIndex3D to VoxelGridID + VoxelGridID aVoxelGridID; + + if (_geoInfo.convert(aIndex, aVoxelGridID)) + { + return getMaskAt(aVoxelGridID, voxel); + } + else + { + return false; + } + } + + const core::GeometricInfo& BoostMaskAccessor::getGeometricInfo() const + { + return _geoInfo; + }; + } + + } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskAccessor.h b/code/masks/rttbBoostMaskAccessor.h similarity index 97% rename from code/masks/boost/rttbBoostMaskAccessor.h rename to code/masks/rttbBoostMaskAccessor.h index 759b50c..83bf37c 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/rttbBoostMaskAccessor.h @@ -1,121 +1,121 @@ -// ----------------------------------------------------------------------- -// 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) -*/ -#ifndef __BOOST_MASK_R_ACCESSOR__H -#define __BOOST_MASK_R_ACCESSOR__H - -#include "rttbBaseType.h" -#include "rttbGeometricInfo.h" -#include "rttbMaskAccessorInterface.h" -#include "rttbStructure.h" - -#include "RTTBBoostMaskExports.h" - -namespace rttb -{ - namespace masks - { - namespace boost - { - /*! @class BoostMaskAccessor - * @brief Using the voxelization based on boost::geometry and generate the mask accessor. - * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. - * (A structure without self intersection 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 RTTBBoostMask_EXPORT BoostMaskAccessor : public core::MaskAccessorInterface - { - public: - using MaskVoxelList = core::MaskAccessorInterface::MaskVoxelList; - using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer; - - using StructTypePointer = core::Structure::Pointer; - - private: - StructTypePointer _spStructure; - core::GeometricInfo _geoInfo; - bool _strict; - - /*! vector containing list of mask voxels*/ - MaskVoxelListPointer _spRelevantVoxelVector; - - IDType _maskUID; - - public: - - - /*! @brief Constructor with a structure pointer and a geometric info pointer - * @param aStructurePointer smart pointer of the structure - * @param aGeometricInfo smart pointer of the geometricInfo of the dose - * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. - * @exception InvalidParameterException thrown if strict is true and the structure has self intersections - */ - BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, - bool strict = true); - - /*! @brief destructor*/ - ~BoostMaskAccessor() override; - - /*! @brief voxelization of the given structures using boost algorithms*/ - void updateMask() override; - - /*! @brief get vector containing all relevant voxels that are inside the given structure*/ - MaskVoxelListPointer getRelevantVoxelVector() override; - - /*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/ - MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold) override; - - /*!@brief determine how a given voxel on the dose grid is masked - * @param aID ID of the voxel in grid. - * @param voxel Reference to the voxel. - * @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined. - * The relevant volume fraction will be set to zero. - * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ - bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const override; - - /*!@brief determine how a given voxel on the dose grid is masked - * @param aIndex 3d index of the voxel in grid. - * @param voxel Reference to the voxel. - * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ - bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const override; - - /*! @brief give access to GeometricInfo*/ - const core::GeometricInfo& getGeometricInfo() const override; - - /* @ brief is true if dose is on a homogeneous grid - * @remark Inhomogeneous grids are not supported at the moment, but if they will be supported in the future the interface does not need to change.*/ - bool isGridHomogeneous() const override - { - return true; - }; - - IDType getMaskUID() const override - { - return _maskUID; - }; - - }; - } - } -} - -#endif +// ----------------------------------------------------------------------- +// 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) +*/ +#ifndef __BOOST_MASK_R_ACCESSOR__H +#define __BOOST_MASK_R_ACCESSOR__H + +#include "rttbBaseType.h" +#include "rttbGeometricInfo.h" +#include "rttbMaskAccessorInterface.h" +#include "rttbStructure.h" + +#include "RTTBBoostMaskExports.h" + +namespace rttb +{ + namespace masks + { + namespace boost + { + /*! @class BoostMaskAccessor + * @brief Using the voxelization based on boost::geometry and generate the mask accessor. + * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. + * (A structure without self intersection 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 RTTBBoostMask_EXPORT BoostMaskAccessor : public core::MaskAccessorInterface + { + public: + using MaskVoxelList = core::MaskAccessorInterface::MaskVoxelList; + using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer; + + using StructTypePointer = core::Structure::Pointer; + + private: + StructTypePointer _spStructure; + core::GeometricInfo _geoInfo; + bool _strict; + + /*! vector containing list of mask voxels*/ + MaskVoxelListPointer _spRelevantVoxelVector; + + IDType _maskUID; + + public: + + + /*! @brief Constructor with a structure pointer and a geometric info pointer + * @param aStructurePointer smart pointer of the structure + * @param aGeometricInfo smart pointer of the geometricInfo of the dose + * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. + * @exception InvalidParameterException thrown if strict is true and the structure has self intersections + */ + BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, + bool strict = true); + + /*! @brief destructor*/ + ~BoostMaskAccessor() override; + + /*! @brief voxelization of the given structures using boost algorithms*/ + void updateMask() override; + + /*! @brief get vector containing all relevant voxels that are inside the given structure*/ + MaskVoxelListPointer getRelevantVoxelVector() override; + + /*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/ + MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold) override; + + /*!@brief determine how a given voxel on the dose grid is masked + * @param aID ID of the voxel in grid. + * @param voxel Reference to the voxel. + * @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined. + * The relevant volume fraction will be set to zero. + * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ + bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const override; + + /*!@brief determine how a given voxel on the dose grid is masked + * @param aIndex 3d index of the voxel in grid. + * @param voxel Reference to the voxel. + * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ + bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const override; + + /*! @brief give access to GeometricInfo*/ + const core::GeometricInfo& getGeometricInfo() const override; + + /* @ brief is true if dose is on a homogeneous grid + * @remark Inhomogeneous grids are not supported at the moment, but if they will be supported in the future the interface does not need to change.*/ + bool isGridHomogeneous() const override + { + return true; + }; + + IDType getMaskUID() const override + { + return _maskUID; + }; + + }; + } + } +} + +#endif diff --git a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp similarity index 97% rename from code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp rename to code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp index cf0f169..02c39ce 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp +++ b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp @@ -1,156 +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, - core::GeometricInfo::Pointer 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) - { - 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)); - } - } - } - } -} +// ----------------------------------------------------------------------- +// 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, + core::GeometricInfo::Pointer 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) + { + 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/rttbBoostMaskGenerateMaskVoxelListThread.h b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h similarity index 97% rename from code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h rename to code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h index b77f70e..76dc429 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h +++ b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h @@ -1,89 +1,89 @@ -// ----------------------------------------------------------------------- -// 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 "rttbMaskAccessorInterface.h" - -#include -#include -#include - -namespace rttb -{ - namespace masks - { - namespace boost - { - /*! @class BoostMaskGenerateMaskVoxelListThread - * - */ - class BoostMaskGenerateMaskVoxelListThread - { - - public: - typedef ::boost::multi_array BoostArray2D; - using BoostArray2DPointer = ::boost::shared_ptr; - typedef ::boost::shared_ptr > BoostArrayMapPointer; - using VoxelIndexVector = std::vector; - using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer; - - BoostMaskGenerateMaskVoxelListThread(const VoxelIndexVector& aGlobalBoundingBox, - core::GeometricInfo::Pointer aGeometricInfo, - BoostArrayMapPointer aVoxelizationMap, - double aVoxelizationThickness, - unsigned int aBeginSlice, - unsigned int aEndSlice, - MaskVoxelListPointer aMaskVoxelList, ::boost::shared_ptr<::boost::shared_mutex> aMutex); - void operator()(); - - private: - VoxelIndexVector _globalBoundingBox; - core::GeometricInfo::Pointer _geometricInfo; - BoostArrayMapPointer _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; - - MaskVoxelListPointer _resultMaskVoxelList; - ::boost::shared_ptr<::boost::shared_mutex> _mutex; - - /*! @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; - }; - - } - - - } -} - +// ----------------------------------------------------------------------- +// 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 "rttbMaskAccessorInterface.h" + +#include +#include +#include + +namespace rttb +{ + namespace masks + { + namespace boost + { + /*! @class BoostMaskGenerateMaskVoxelListThread + * + */ + class BoostMaskGenerateMaskVoxelListThread + { + + public: + typedef ::boost::multi_array BoostArray2D; + using BoostArray2DPointer = ::boost::shared_ptr; + typedef ::boost::shared_ptr > BoostArrayMapPointer; + using VoxelIndexVector = std::vector; + using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer; + + BoostMaskGenerateMaskVoxelListThread(const VoxelIndexVector& aGlobalBoundingBox, + core::GeometricInfo::Pointer aGeometricInfo, + BoostArrayMapPointer aVoxelizationMap, + double aVoxelizationThickness, + unsigned int aBeginSlice, + unsigned int aEndSlice, + MaskVoxelListPointer aMaskVoxelList, ::boost::shared_ptr<::boost::shared_mutex> aMutex); + void operator()(); + + private: + VoxelIndexVector _globalBoundingBox; + core::GeometricInfo::Pointer _geometricInfo; + BoostArrayMapPointer _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; + + MaskVoxelListPointer _resultMaskVoxelList; + ::boost::shared_ptr<::boost::shared_mutex> _mutex; + + /*! @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/rttbBoostMaskVoxelizationThread.cpp similarity index 97% rename from code/masks/boost/rttbBoostMaskVoxelizationThread.cpp rename to code/masks/rttbBoostMaskVoxelizationThread.cpp index 9b4deec..d184883 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp +++ b/code/masks/rttbBoostMaskVoxelizationThread.cpp @@ -1,174 +1,174 @@ -// ----------------------------------------------------------------------- -// 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, bool strict) : _geometryCoordinateBoostPolygonMap(APolygonMap), - _globalBoundingBox(aGlobalBoundingBox), _resultVoxelization(anArrayMap), _mutex(aMutex), _strict(strict) - { - } - - 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 (auto & it : _geometryCoordinateBoostPolygonMap) - { - 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); - volumeFraction = correctForErrorAndStrictness(volumeFraction, _strict); - if (volumeFraction < 0 || volumeFraction > 1 ) - { - 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; - } - - double BoostMaskVoxelizationThread::correctForErrorAndStrictness(double volumeFraction, bool strict) const - { - if (strict){ - if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) - { - volumeFraction = 1; - } - } - else { - if (volumeFraction > 1){ - volumeFraction = 1; - } - else if (volumeFraction < 0){ - volumeFraction = 0; - } - } - return volumeFraction; - } - - } - } -} +// ----------------------------------------------------------------------- +// 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, bool strict) : _geometryCoordinateBoostPolygonMap(APolygonMap), + _globalBoundingBox(aGlobalBoundingBox), _resultVoxelization(anArrayMap), _mutex(aMutex), _strict(strict) + { + } + + 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 (auto & it : _geometryCoordinateBoostPolygonMap) + { + 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); + volumeFraction = correctForErrorAndStrictness(volumeFraction, _strict); + if (volumeFraction < 0 || volumeFraction > 1 ) + { + 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; + } + + double BoostMaskVoxelizationThread::correctForErrorAndStrictness(double volumeFraction, bool strict) const + { + if (strict){ + if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) + { + volumeFraction = 1; + } + } + else { + if (volumeFraction > 1){ + volumeFraction = 1; + } + else if (volumeFraction < 0){ + volumeFraction = 0; + } + } + return volumeFraction; + } + + } + } +} diff --git a/code/masks/boost/rttbBoostMaskVoxelizationThread.h b/code/masks/rttbBoostMaskVoxelizationThread.h similarity index 97% rename from code/masks/boost/rttbBoostMaskVoxelizationThread.h rename to code/masks/rttbBoostMaskVoxelizationThread.h index 9743c4e..3ffd14c 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.h +++ b/code/masks/rttbBoostMaskVoxelizationThread.h @@ -1,109 +1,109 @@ -// ----------------------------------------------------------------------- -// 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 - -#include "rttbBaseType.h" - -#include -#include -#include -#include -#include - -namespace rttb -{ - namespace masks - { - namespace boost - { - /*! @class BoostMaskGenerateMaskVoxelListThread - * - */ - class BoostMaskVoxelizationThread - { - - public: - using BoostPolygon2D = ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy >; - using BoostPolygonVector = std::vector;//polygon with or without holes - typedef std::map - BoostPolygonMap;//map of the z index with the vector of boost 2d polygon - using VoxelIndexVector = std::vector; - typedef ::boost::multi_array BoostArray2D; - using BoostArray2DPointer = ::boost::shared_ptr; - typedef ::boost::shared_ptr > BoostArrayMapPointer; - - /*! @brief Constructor - * @param aMutex a mutex for thread-safe handling of the _resultVoxelization - * @param strict true means that volumeFractions of <0 and >1 are NOT corrected. Otherwise, they are automatically corrected to 0 or 1, respectively. - */ - BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, - const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex, bool strict); - - void operator()(); - - - private: - using BoostPolygonDeque = std::deque; - using BoostRing2D = ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy >; - using BoostPoint2D = ::boost::geometry::model::d2::point_xy; - - - BoostPolygonMap _geometryCoordinateBoostPolygonMap; - VoxelIndexVector _globalBoundingBox; - BoostArrayMapPointer _resultVoxelization; - ::boost::shared_ptr<::boost::shared_mutex> _mutex; - bool _strict; - - /*! @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 intersection polygons of the structure and the voxel - */ - static BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, - const BoostPolygonVector& intersectionSlicePolygons); - - /*! @brief Get the voxel 2d contour polygon in geometry coordinate*/ - static BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D); - - /*! @brief Calculate the area of all polygons - * @param aPolygonDeque The deque of polygons - * @return Return the area of all polygons - */ - static double calcArea(const BoostPolygonDeque& aPolygonDeque); - /*! @brief Corrects the volumeFraction - * @details the volume fraction is corrected in case of strict=true. Otherwise, it's only corrected for double imprecision - * @return The corrected volumeFraction - */ - double correctForErrorAndStrictness(double volumeFraction, bool strict) const; - }; - - } - - - } -} - +// ----------------------------------------------------------------------- +// 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 + +#include "rttbBaseType.h" + +#include +#include +#include +#include +#include + +namespace rttb +{ + namespace masks + { + namespace boost + { + /*! @class BoostMaskGenerateMaskVoxelListThread + * + */ + class BoostMaskVoxelizationThread + { + + public: + using BoostPolygon2D = ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy >; + using BoostPolygonVector = std::vector;//polygon with or without holes + typedef std::map + BoostPolygonMap;//map of the z index with the vector of boost 2d polygon + using VoxelIndexVector = std::vector; + typedef ::boost::multi_array BoostArray2D; + using BoostArray2DPointer = ::boost::shared_ptr; + typedef ::boost::shared_ptr > BoostArrayMapPointer; + + /*! @brief Constructor + * @param aMutex a mutex for thread-safe handling of the _resultVoxelization + * @param strict true means that volumeFractions of <0 and >1 are NOT corrected. Otherwise, they are automatically corrected to 0 or 1, respectively. + */ + BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, + const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex, bool strict); + + void operator()(); + + + private: + using BoostPolygonDeque = std::deque; + using BoostRing2D = ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy >; + using BoostPoint2D = ::boost::geometry::model::d2::point_xy; + + + BoostPolygonMap _geometryCoordinateBoostPolygonMap; + VoxelIndexVector _globalBoundingBox; + BoostArrayMapPointer _resultVoxelization; + ::boost::shared_ptr<::boost::shared_mutex> _mutex; + bool _strict; + + /*! @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 intersection polygons of the structure and the voxel + */ + static BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, + const BoostPolygonVector& intersectionSlicePolygons); + + /*! @brief Get the voxel 2d contour polygon in geometry coordinate*/ + static BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D); + + /*! @brief Calculate the area of all polygons + * @param aPolygonDeque The deque of polygons + * @return Return the area of all polygons + */ + static double calcArea(const BoostPolygonDeque& aPolygonDeque); + /*! @brief Corrects the volumeFraction + * @details the volume fraction is corrected in case of strict=true. Otherwise, it's only corrected for double imprecision + * @return The corrected volumeFraction + */ + double correctForErrorAndStrictness(double volumeFraction, bool strict) const; + }; + + } + + + } +} + #endif \ No newline at end of file