diff --git a/code/masks/rttbBoostMask.cpp b/code/masks/rttbBoostMask.cpp
index af20a7d..05eeba3 100644
--- a/code/masks/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 <limits>
-
-#include <boost/geometry/geometries/register/point.hpp>
-#include <boost/geometry/geometries/register/ring.hpp>
-#include <boost/make_shared.hpp>
-#include <boost/thread/thread.hpp>
-
-#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<MaskVoxelList>())
- {
-
- _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<double>::min(),
- std::numeric_limits<double>::min(), std::numeric_limits<double>::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<BoostPolygonMap> 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<double, BoostPolygonVector>((*itMap).first,
- polygonVector));
- count++;
-
- }
-
- _voxelizationMap = ::boost::make_shared<std::map<double, BoostArray2DPointer> >();
-
- 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<double, BoostRingVector>(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<unsigned int> 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 <limits>
+
+#include <boost/geometry/geometries/register/point.hpp>
+#include <boost/geometry/geometries/register/ring.hpp>
+#include <boost/make_shared.hpp>
+#include <boost/thread/thread.hpp>
+
+#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<MaskVoxelList>())
+ {
+
+ _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<double>::min(),
+ std::numeric_limits<double>::min(), std::numeric_limits<double>::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<BoostPolygonMap> 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<double, BoostPolygonVector>((*itMap).first,
+ polygonVector));
+ count++;
+
+ }
+
+ _voxelizationMap = ::boost::make_shared<std::map<double, BoostArray2DPointer> >();
+
+ 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, _strict, 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<double, BoostRingVector>(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<unsigned int> 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/rttbBoostMaskGenerateMaskVoxelListThread.cpp b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp
index 02c39ce..db6db0b 100644
--- a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp
+++ b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp
@@ -1,156 +1,157 @@
-// -----------------------------------------------------------------------
-// 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 <boost/thread.hpp>
-
-
-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<core::MaskVoxel> maskVoxelsInThread;
-
- for (unsigned int indexZ = _beginSlice; indexZ < _endSlice; ++indexZ)
- {
- //calculate weight vector
- std::map<double, double> 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<double, double>& 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<double, double>(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 <boost/thread.hpp>
+
+
+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,
+ bool strictVoxelization,
+ ::boost::shared_ptr<::boost::shared_mutex> aMutex) :
+ _globalBoundingBox(aGlobalBoundingBox), _geometricInfo(aGeometricInfo),
+ _voxelizationMap(aVoxelizationMap), _voxelizationThickness(aVoxelizationThickness),
+ _beginSlice(aBeginSlice), _endSlice(aEndSlice),
+ _resultMaskVoxelList(aMaskVoxelList), _strictVoxelization(strictVoxelization), _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<core::MaskVoxel> maskVoxelsInThread;
+
+ for (unsigned int indexZ = _beginSlice; indexZ < _endSlice; ++indexZ)
+ {
+ //calculate weight vector
+ std::map<double, double> 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 || !_strictVoxelization))
+ {
+ 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<double, double>& 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<double, double>(it.first, weight));
+ }
+ }
+ }
+ }
+}
diff --git a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h
index 76dc429..7be9b2f 100644
--- a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h
+++ b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h
@@ -1,89 +1,91 @@
-// -----------------------------------------------------------------------
-// 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 <boost/multi_array.hpp>
-#include <boost/shared_ptr.hpp>
-#include <boost/thread/shared_mutex.hpp>
-
-namespace rttb
-{
- namespace masks
- {
- namespace boost
- {
- /*! @class BoostMaskGenerateMaskVoxelListThread
- *
- */
- class BoostMaskGenerateMaskVoxelListThread
- {
-
- public:
- typedef ::boost::multi_array<double, 2> BoostArray2D;
- using BoostArray2DPointer = ::boost::shared_ptr<BoostArray2D>;
- typedef ::boost::shared_ptr<std::map<double, BoostArray2DPointer> > BoostArrayMapPointer;
- using VoxelIndexVector = std::vector<rttb::VoxelGridIndex3D>;
- 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<double, double>& 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 <boost/multi_array.hpp>
+#include <boost/shared_ptr.hpp>
+#include <boost/thread/shared_mutex.hpp>
+
+namespace rttb
+{
+ namespace masks
+ {
+ namespace boost
+ {
+ /*! @class BoostMaskGenerateMaskVoxelListThread
+ *
+ */
+ class BoostMaskGenerateMaskVoxelListThread
+ {
+
+ public:
+ typedef ::boost::multi_array<double, 2> BoostArray2D;
+ using BoostArray2DPointer = ::boost::shared_ptr<BoostArray2D>;
+ typedef ::boost::shared_ptr<std::map<double, BoostArray2DPointer> > BoostArrayMapPointer;
+ using VoxelIndexVector = std::vector<rttb::VoxelGridIndex3D>;
+ 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, bool strictVoxelization, ::boost::shared_ptr<::boost::shared_mutex> aMutex);
+ void operator()();
+
+ private:
+ VoxelIndexVector _globalBoundingBox;
+ core::GeometricInfo::Pointer _geometricInfo;
+ BoostArrayMapPointer _voxelizationMap;
+ bool _strictVoxelization=true;
+ //(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<double, double>& weightVector) const;
+
+ };
+
+ }
+
+
+ }
+}
+
#endif
\ No newline at end of file