diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp index 76bc0ca..9d6786d 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/boost/rttbBoostMask.cpp @@ -1,580 +1,561 @@ -#include -#include -#include +// ----------------------------------------------------------------------- +// 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 #include #include "rttbBoostMask.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbBoostMaskVoxelizationThread.h" - namespace rttb { namespace masks { namespace boost { BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, BoostMask::StructPointer aStructure, bool strict, unsigned int numberOfThreads) : _geometricInfo(aDoseGeoInfo), _structure(aStructure), - _strict(strict), _numberOfThreads(numberOfThreads), + _strict(strict), _numberOfThreads(numberOfThreads), _voxelizationThickness(0.0), _voxelInStructure(::boost::make_shared()) { _isUpToDate = false; if (!_geometricInfo) { throw rttb::core::NullPointerException("Error: Geometric info is NULL!"); } else if (!_structure) { throw rttb::core::NullPointerException("Error: Structure is NULL!"); } if (_numberOfThreads == 0) { _numberOfThreads = ::boost::thread::hardware_concurrency(); + if (_numberOfThreads == 0) + { + throw rttb::core::InvalidParameterException("Error: detection of the number of hardware threads is not possible. Please specify number of threads for voxelization explicitly as parameter in BoostMask."); + } } - if (_numberOfThreads <= 0) - { - throw rttb::core::InvalidParameterException("Error: the given number of threads must > 0, or detection of the number of hardwore thread is not possible."); - } - - std::cout << "number of threads: " << _numberOfThreads << std::endl; - } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { if (!_isUpToDate) { calcMask(); } return _voxelInStructure; } void BoostMask::calcMask() { preprocessing(); voxelization(); generateMaskVoxelList(); _isUpToDate = true; } void BoostMask::preprocessing() { rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); //Convert world coordinate polygons to the polygons with geometry coordinate rttb::PolygonSequenceType geometryCoordinatePolygonVector; rttb::PolygonSequenceType::iterator it; rttb::DoubleVoxelGridIndex3D globalMaxGridIndex(std::numeric_limits::min(), std::numeric_limits::min(), std::numeric_limits::min()); rttb::DoubleVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(), _geometricInfo->getNumRows(), 0); for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) { PolygonType rttbPolygon = *it; PolygonType geometryCoordinatePolygon; //1. convert polygon to geometry coordinate polygons //2. calculate global min/max //3. check if polygon is planar if (!preprocessingPolygon(rttbPolygon, geometryCoordinatePolygon, globalMinGridIndex, globalMaxGridIndex, errorConstant)) { throw rttb::core::Exception("TiltedMaskPlaneException"); } geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon); } rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) + 0.5), GridIndexType(globalMinGridIndex(1) + 0.5), GridIndexType(globalMinGridIndex(2) + 0.5)); rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) + 0.5), GridIndexType(globalMaxGridIndex(1) + 0.5), GridIndexType(globalMaxGridIndex(2) + 0.5)); _globalBoundingBox.push_back(minIndex); _globalBoundingBox.push_back(maxIndex); //convert rttb polygon sequence to a map of z index and a vector of boost ring 2d (without holes) _ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); } void BoostMask::voxelization() { - BoostPolygonMap::iterator it; if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } - unsigned int ringMapSize = _ringMap.size(); - BoostRingMap::iterator itMap; - unsigned int mapSizeInAThread = _ringMap.size() / _numberOfThreads; + 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; - BoostMaskVoxelizationThread::VoxelizationQueuePointer resultQueue - = ::boost::make_shared<::boost::lockfree::queue> - (ringMapSize); - BoostMaskVoxelizationThread::VoxelizationIndexQueuePointer resultIndexQueue - = ::boost::make_shared<::boost::lockfree::queue> - (ringMapSize); + auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); for (int i = 0; i < polygonMapVector.size(); ++i) { - BoostMaskVoxelizationThread t(polygonMapVector.at(i), _globalBoundingBox, resultIndexQueue, - resultQueue); + BoostMaskVoxelizationThread t(polygonMapVector.at(i), _globalBoundingBox, + _voxelizationMap, aMutex); threads.create_thread(t); } threads.join_all(); - - - BoostArray2D* array; - double index; - - //generate voxelization map using resultQueue and resultIndexQueue - while (resultIndexQueue->pop(index) && resultQueue->pop(array)) - { - _voxelizationMap.insert(std::pair(index, *array)); - } - } void BoostMask::generateMaskVoxelList() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } - //check homogeneus of the voxelization plane (the contours plane) + //check homogeneous of the voxelization plane (the contours plane) if (!calcVoxelizationThickness(_voxelizationThickness)) { - throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneus!"); + throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneous!"); } - - ::boost::shared_ptr<::boost::lockfree::queue> resultQueue - = ::boost::make_shared<::boost::lockfree::queue> - (_geometricInfo->getNumberOfVoxels()); - - ::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, - resultQueue); + _voxelInStructure, aMutex); threads.create_thread(t); } threads.join_all(); - core::MaskVoxel* voxel; - - //generate _voxelInStructure using resultQueue - while (resultQueue->pop(voxel)) - { - _voxelInStructure->push_back(*voxel);//push back the mask voxel in structure - } - } bool BoostMask::preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const { double minZ = _geometricInfo->getNumSlices(); double maxZ = 0.0; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D worldCoordinatePoint = aRTTBPolygon.at(i); //convert to geometry coordinate polygon rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint; - _geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, - geometryCoordinatePoint); + _geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, geometryCoordinatePoint); + geometryCoordinatePolygon.push_back(geometryCoordinatePoint); //calculate the current global min/max //min and max for x if (geometryCoordinatePoint(0) < minimum(0)) { minimum(0) = geometryCoordinatePoint(0); } if (geometryCoordinatePoint(0) > maximum(0)) { maximum(0) = geometryCoordinatePoint(0); } //min and max for y if (geometryCoordinatePoint(1) < minimum(1)) { minimum(1) = geometryCoordinatePoint(1); } if (geometryCoordinatePoint(1) > maximum(1)) { maximum(1) = geometryCoordinatePoint(1); } //min and max for z if (geometryCoordinatePoint(2) < minimum(2)) { minimum(2) = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maximum(2)) { maximum(2) = geometryCoordinatePoint(2); } //check planar if (geometryCoordinatePoint(2) < minZ) { minZ = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maxZ) { maxZ = geometryCoordinatePoint(2); } } return (abs(maxZ - minZ) <= aErrorConstant); } BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const { BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); if (i == 0) { firstPoint = boostPoint; } ::boost::geometry::append(polygon2D, boostPoint); } ::boost::geometry::append(polygon2D, firstPoint); return polygon2D; } BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( - const rttb::PolygonSequenceType& aRTTBPolygonVector) + const rttb::PolygonSequenceType& aRTTBPolygonVector) const { rttb::PolygonSequenceType::const_iterator it; BoostMask::BoostRingMap aRingMap; for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it) { rttb::PolygonType rttbPolygon = *it; double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon bool isFirstZ = true; if (!aRingMap.empty()) { BoostMask::BoostRingMap::iterator findIt = findNearestKey(aRingMap, zIndex, errorConstant); //if the z index is found (same slice), add the polygon to vector if (findIt != aRingMap.end()) { //BoostRingVector ringVector = ; (*findIt).second.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); isFirstZ = false; } } //if it is the first z index in the map, insert vector with the polygon if (isFirstZ) { BoostRingVector ringVector; ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); aRingMap.insert(std::pair(zIndex, ringVector)); } } return aRingMap; } BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap& - aBoostRingMap, double aIndex, double aErrorConstant) + aBoostRingMap, double aIndex, double aErrorConstant) const { BoostMask::BoostRingMap::iterator find = aBoostRingMap.find(aIndex); //if find a key equivalent to aIndex, found if (find != aBoostRingMap.end()) { return find; } else { BoostMask::BoostRingMap::iterator lowerBound = aBoostRingMap.lower_bound(aIndex); //if all keys go before aIndex, check the last key if (lowerBound == aBoostRingMap.end()) { lowerBound = --aBoostRingMap.end(); } //if the lower bound very close to aIndex, found if (abs((*lowerBound).first - aIndex) <= aErrorConstant) { return lowerBound; } else { //if the lower bound is the beginning, not found if (lowerBound == aBoostRingMap.begin()) { return aBoostRingMap.end(); } else { BoostMask::BoostRingMap::iterator lowerBound1 = --lowerBound;//the key before the lower bound //if the key before the lower bound very close to a Index, found if (abs((*lowerBound1).first - aIndex) <= aErrorConstant) { return lowerBound1; } //else, not found else { return aBoostRingMap.end(); } } } } } BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector) const { //check donut BoostMask::BoostRingVector::const_iterator it1; BoostMask::BoostRingVector::const_iterator it2; BoostMask::BoostPolygonVector boostPolygonVector; std::vector donutIndexVector;//store the outer and inner ring index BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon //Get donut index and donut polygon unsigned int index1 = 0; - for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) + for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } //if not jet, check now if (!it1IsDonut) { bool it2IsDonut = false; unsigned int index2 = 0; - for (it2 = aRingVector.begin(); it2 != aRingVector.end(); it2++, index2++) + 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++) + for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++) { bool it1IsDonut = false; //check if the ring is the outer or inner of a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } if (!it1IsDonut) { BoostMask::BoostPolygon2D polygon2D; ::boost::geometry::append(polygon2D, *it1); boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut } } //Append donut polygon to boostPolygonVector BoostMask::BoostPolygonVector::iterator itDonut; - for (itDonut = donutVector.begin(); itDonut != donutVector.end(); itDonut++) + for (itDonut = donutVector.begin(); itDonut != donutVector.end(); ++itDonut) { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } bool BoostMask::calcVoxelizationThickness(double& aThickness) const { - BoostArrayMap::const_iterator it = _voxelizationMap.begin(); - BoostArrayMap::const_iterator it2 = ++_voxelizationMap.begin(); - if (_voxelizationMap.size() <= 1) + if (_voxelizationMap->size() <= 1) { aThickness = 1; return true; } double thickness = 0; + auto it = _voxelizationMap->cbegin(); + auto it2 = ++_voxelizationMap->cbegin(); for (; - it != _voxelizationMap.end(), it2 != _voxelizationMap.end(); ++it, ++it2) + it != _voxelizationMap->cend(), it2 != _voxelizationMap->cend(); ++it, ++it2) { if (thickness == 0) { - thickness = (*it2).first - (*it).first; + thickness = it2->first - it->first; } else { - double curThickness = (*it2).first - (*it).first; - //if no homogeneous (leave out double imprecisions), return false + double curThickness = it2->first - it->first; + //if not homogeneous (leave out double imprecisions), return false if (abs(thickness-curThickness)>errorConstant) { return false; } } } if (thickness != 0) { aThickness = thickness; } else { aThickness = 1; } return true; } } } } diff --git a/code/masks/boost/rttbBoostMask.h b/code/masks/boost/rttbBoostMask.h index b484539..41ceed3 100644 --- a/code/masks/boost/rttbBoostMask.h +++ b/code/masks/boost/rttbBoostMask.h @@ -1,210 +1,205 @@ // ----------------------------------------------------------------------- // 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_H -#define __BOOST_MASK_H +#ifndef __BOOST_MASK_R_H +#define __BOOST_MASK_R_H #include "rttbBaseType.h" #include "rttbStructure.h" #include "rttbGeometricInfo.h" -#include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include #include -#include -#include #include #include namespace rttb { namespace masks { namespace 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 interseciton means all contours of the structure have no self intersection, and + * (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: typedef ::boost::shared_ptr GeometricInfoPointer; typedef core::Structure::StructTypePointer StructPointer; typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; /*! @brief Constructor * @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is NULL * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. * @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(GeometricInfoPointer aDoseGeoInfo, StructPointer 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: typedef ::boost::geometry::model::d2::point_xy BoostPoint2D; typedef ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy > BoostPolygon2D; typedef ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy > BoostRing2D; - typedef std::deque BoostPolygonDeque; typedef std::vector BoostRingVector;//polygon without holes typedef std::vector BoostPolygonVector;//polygon with or without holes typedef std::vector VoxelIndexVector; typedef std::map BoostPolygonMap;//map of the z index with the vector of boost 2d polygon typedef std::map BoostRingMap;//map of the z index with the vector of boost 2d ring typedef ::boost::multi_array BoostArray2D; - typedef std::map BoostArrayMap; - typedef ::boost::array BoostArray1D; - typedef std::vector > WeightVector; + typedef ::boost::shared_ptr BoostArray2DPointer; + typedef ::boost::shared_ptr > BoostArrayMapPointer; GeometricInfoPointer _geometricInfo; StructPointer _structure; /*! @brief The map of z index and a vector of boost ring 2d (without holes) * 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. * The first index has the minimum for x/y/z of the global bounding box. * The second index has the maximum for x/y/z of the global bounding index. */ VoxelIndexVector _globalBoundingBox; /*! @brief The voxelization map * key: the converted double z grid index of a contour plane * value: the 2d mask, array[i][j] = the mask value of the position (i,j) in the global bounding box, * i: 0 - (_globalBoundingBoxSize0-1), j: 0 - (_globalBoundingBoxSize1-1) */ - BoostArrayMap _voxelizationMap; + BoostArrayMapPointer _voxelizationMap; - //The thickness of the voxelization plane (the contour plane), in double dose grid index + //@brief The thickness of the voxelization plane (the contour plane), in double dose grid index //(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2) double _voxelizationThickness; bool _strict; - //vector of the MaskVoxel inside the structure + //@brief vector of the MaskVoxel inside the structure MaskVoxelListPointer _voxelInStructure; /*! @brief If the mask is up to date */ bool _isUpToDate; /*! @brief The number of threads */ unsigned int _numberOfThreads; /*! @brief Voxelization and generate mask */ void calcMask(); /*! @brief The preprocessing step, wich consists of the following logic and Sub setps: * For all contours in a struct: * 1) Transfer the contour polygons into boost::geometry structures * 1a) Convert the contur points from world coordinates into geometry coordinates. * 1b) get min and max for x/y/z of a contour * 2) Tilt check: if difference of z_min and z_max is larger then a tolerance value -> there is a tilt. Throw rttb::TiltedMaskPlaneException. * 3) Get struct-bounding-box: get x_min_struct, y_min_struct, x_max_struct, y_max_struct to define the bounding box that containes all contours of a struct in x-y-dimensions. */ void preprocessing(); - /*! @brief The voxelization step, wich computes the voxelization planes (in x/y) for all contours of an struct. + /*! @brief The voxelization step, which computes the voxelization planes (in x/y) for all contours of an struct. * For each contour (that is in the z-Range of the reference geometry) of the struct: * 1) Allocate result array (voxelization plane) based on the bounding box (see Preprocessing Step 3) * 2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry). * 3) Add result Array (key is the z-Value of the contour) */ void voxelization(); - /*! @final mask voxel Generation step which transfers the voxilization planes into the (z-)geometry of the reference geometry. + /*! @final mask voxel Generation step which transfers the voxelization planes into the (z-)geometry of the reference geometry. * It consists of following Sub steps : * For all "slices" in the reference geometry : * 1) generate weight vector for all voxelization planes for a given z - value of a slice - * Itterate over the bounding box of a struct.For each voxel : + * 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 corrdinate to the rttb polygon with double geometry coordinate, calculate the current min/max + /*! @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 minminal and maximal z-coordinate of the polygon is not larger than a error constant + * @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); + 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); + 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.h b/code/masks/boost/rttbBoostMaskAccessor.h index 360891b..22990d2 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/boost/rttbBoostMaskAccessor.h @@ -1,129 +1,127 @@ // ----------------------------------------------------------------------- // 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_ACCESSOR__H -#define __BOOST_MASK_ACCESSOR__H +#ifndef __BOOST_MASK_R_ACCESSOR__H +#define __BOOST_MASK_R_ACCESSOR__H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include "rttbGenericDoseIterator.h" #include "rttbStructure.h" -#include - #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 interseciton means all contours of the structure have no self intersection, and + * (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: typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::Structure::StructTypePointer StructTypePointer; private: core::GeometricInfo _geoInfo; /*! vector containing list of mask voxels*/ MaskVoxelListPointer _spRelevantVoxelVector; StructTypePointer _spStructure; IDType _maskUID; bool _strict; public: /*! @brief Constructor with a structure pointer and a geometric info pointer * @param aStructurePointer smart pointer of the structure - * @param aGeometricInfoPtr smart pointer of the geometricinfo of the dose + * @param aGeometricInfoPtr 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(); /*! @brief voxelization of the given structures using boost algorithms*/ void updateMask(); /*! @brief get vector containing all relevant voxels that are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(); /*! @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); /*!@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; /*!@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; /*! @brief give access to GeometricInfo*/ const core::GeometricInfo& getGeometricInfo() const; /* @ 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 { return true; }; IDType getMaskUID() const { return _maskUID; }; }; } } } #endif diff --git a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp index 4edbf29..9ecccdc 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp +++ b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp @@ -1,133 +1,155 @@ +// ----------------------------------------------------------------------- +// RTToolbox - DKFZ radiotherapy quantitative evaluation library +// +// Copyright (c) German Cancer Research Center (DKFZ), +// Software development for Integrated Diagnostics and Therapy (SIDT). +// ALL RIGHTS RESERVED. +// See rttbCopyright.txt or +// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notices for more information. +// +//------------------------------------------------------------------------ +/*! +// @file +// @version $Revision: 1127 $ (last changed revision) +// @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ #include "rttbBoostMaskGenerateMaskVoxelListThread.h" - #include "rttbInvalidParameterException.h" -#include +#include namespace rttb { namespace masks { namespace boost { BoostMaskGenerateMaskVoxelListThread::BoostMaskGenerateMaskVoxelListThread( const VoxelIndexVector& aGlobalBoundingBox, GeometricInfoPointer aGeometricInfo, - const BoostArrayMap& aVoxelizationMap, + BoostArrayMapPointer aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, - MaskVoxelQueuePointer aResultMaskVoxelQueue) : + MaskVoxelListPointer aMaskVoxelList, + ::boost::shared_ptr<::boost::shared_mutex> aMutex) : _globalBoundingBox(aGlobalBoundingBox), _geometricInfo(aGeometricInfo), _voxelizationMap(aVoxelizationMap), _voxelizationThickness(aVoxelizationThickness), _beginSlice(aBeginSlice), _endSlice(aEndSlice), - _resultMaskVoxelQueue(aResultMaskVoxelQueue) + _resultMaskVoxelList(aMaskVoxelList), _mutex(aMutex) {} void BoostMaskGenerateMaskVoxelListThread::operator()() { rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); - int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; - int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; + 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; - std::map::iterator it; - BoostArrayMap::iterator itMap; + auto it = weightVectorForZ.cbegin(); + auto itMap = _voxelizationMap->cbegin(); - for (it = weightVectorForZ.begin(), itMap = _voxelizationMap.begin(); it != weightVectorForZ.end() - && itMap != _voxelizationMap.end(); ++it, ++itMap) + for (; it != weightVectorForZ.cend() + && itMap != _voxelizationMap->cend(); ++it, ++itMap) { - double weight = (*it).second; - - BoostArray2D voxelizationArray = (*itMap).second; - //calc sum of all voxelization plane, use weight - volumeFraction += voxelizationArray[x][y] * weight; + 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 = new core::MaskVoxel(gridID, volumeFraction); - _resultMaskVoxelQueue->push(maskVoxelPtr);//push back the mask voxel in structure + 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 { - BoostArrayMap::const_iterator it; - double indexZMin = aIndexZ - 0.5; double indexZMax = aIndexZ + 0.5; - for (it = _voxelizationMap.begin(); it != _voxelizationMap.end(); ++it) + for (auto it = _voxelizationMap->begin(); it != _voxelizationMap->end(); ++it) { double voxelizationPlaneIndexMin = (*it).first - 0.5 * _voxelizationThickness; double voxelizationPlaneIndexMax = (*it).first + 0.5 * _voxelizationThickness; double weight = 0; if ((voxelizationPlaneIndexMin < indexZMin) && (voxelizationPlaneIndexMax > indexZMin)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = voxelizationPlaneIndexMax - indexZMin; } else { weight = 1; } } else if ((voxelizationPlaneIndexMin >= indexZMin) && (voxelizationPlaneIndexMin < indexZMax)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = _voxelizationThickness; } else { weight = indexZMax - voxelizationPlaneIndexMin; } } - weightVector.insert(std::pair((*it).first, weight)); + weightVector.insert(std::pair(it->first, weight)); } } } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h index a311a4e..e157429 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h +++ b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h @@ -1,88 +1,90 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 1127 $ (last changed revision) // @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H #define __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" -#include "rttbMaskVoxel.h" +#include "rttbMaskAccessorInterface.h" #include #include -#include +#include namespace rttb { namespace masks { namespace boost { /*! @class BoostMaskGenerateMaskVoxelListThread * */ class BoostMaskGenerateMaskVoxelListThread { public: typedef ::boost::shared_ptr GeometricInfoPointer; typedef ::boost::multi_array BoostArray2D; - typedef std::map BoostArrayMap; - typedef ::boost::shared_ptr<::boost::lockfree::queue> MaskVoxelQueuePointer; + typedef ::boost::shared_ptr BoostArray2DPointer; + typedef ::boost::shared_ptr > BoostArrayMapPointer; typedef std::vector VoxelIndexVector; + typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; BoostMaskGenerateMaskVoxelListThread(const VoxelIndexVector& aGlobalBoundingBox, GeometricInfoPointer aGeometricInfo, - const BoostArrayMap& aVoxelizationMap, + BoostArrayMapPointer aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, - MaskVoxelQueuePointer aResultMaskVoxelQueue); + MaskVoxelListPointer aMaskVoxelList, ::boost::shared_ptr<::boost::shared_mutex> aMutex); void operator()(); private: VoxelIndexVector _globalBoundingBox; GeometricInfoPointer _geometricInfo; - BoostArrayMap _voxelizationMap; + 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; - MaskVoxelQueuePointer _resultMaskVoxelQueue; + 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/boost/rttbBoostMaskVoxelizationThread.cpp index b402fa8..8aa32fd 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp +++ b/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp @@ -1,135 +1,159 @@ +// ----------------------------------------------------------------------- +// RTToolbox - DKFZ radiotherapy quantitative evaluation library +// +// Copyright (c) German Cancer Research Center (DKFZ), +// Software development for Integrated Diagnostics and Therapy (SIDT). +// ALL RIGHTS RESERVED. +// See rttbCopyright.txt or +// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html +// +// This software is distributed WITHOUT ANY WARRANTY; without even +// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +// PURPOSE. See the above copyright notices for more information. +// +//------------------------------------------------------------------------ +/*! +// @file +// @version $Revision: 1127 $ (last changed revision) +// @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date) +// @author $Author: hentsch $ (last changed by) +*/ + #include "rttbBoostMaskVoxelizationThread.h" +#include "rttbInvalidParameterException.h" + +#include +#include +#include + namespace rttb { namespace masks { namespace boost { BoostMaskVoxelizationThread::BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, - const VoxelIndexVector& aGlobalBoundingBox, VoxelizationIndexQueuePointer aResultIndexQueue, - VoxelizationQueuePointer aVoxelizationQueue) : _geometryCoordinateBoostPolygonMap(APolygonMap), - _globalBoundingBox(aGlobalBoundingBox), _resultIndexQueue(aResultIndexQueue), - _resultVoxelizationQueue(aVoxelizationQueue) + const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex) : _geometryCoordinateBoostPolygonMap(APolygonMap), + _globalBoundingBox(aGlobalBoundingBox), _resultVoxelization(anArrayMap), _mutex(aMutex) { } - void BoostMaskVoxelizationThread::operator()() { rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); - int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; - int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; + const unsigned int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; + const unsigned int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; - BoostPolygonMap::iterator it; + std::map > voxelizationMapInThread; - for (it = _geometryCoordinateBoostPolygonMap.begin(); + for (BoostPolygonMap::iterator it = _geometryCoordinateBoostPolygonMap.begin(); it != _geometryCoordinateBoostPolygonMap.end(); ++it) { - BoostArray2D* maskArray = new BoostArray2D( - ::boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); + BoostArray2D maskArray(::boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); - BoostPolygonVector boostPolygonVec = (*it).second; + BoostPolygonVector boostPolygonVec = it->second; for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) { for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = x + minIndex[0]; currentIndex[1] = y + minIndex[1]; currentIndex[2] = 0; //Get intersection polygons of the dose voxel and the structure BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec); //Calc areas of all intersection polygons double volumeFraction = calcArea(polygons); if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } - (*maskArray)[x][y] = volumeFraction; + maskArray[x][y] = volumeFraction; } } - //insert into voxelization map - _resultIndexQueue->push((*it).first); - _resultVoxelizationQueue->push(maskArray); + voxelizationMapInThread.insert(std::pair(it->first, ::boost::make_shared(maskArray))); } + //insert gathered values into voxelization map + ::boost::unique_lock<::boost::shared_mutex> lock(*_mutex); + _resultVoxelization->insert(voxelizationMapInThread.begin(), voxelizationMapInThread.end()); } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMaskVoxelizationThread::BoostPolygonDeque BoostMaskVoxelizationThread::getIntersections( const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) { BoostMaskVoxelizationThread::BoostPolygonDeque polygonDeque; BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); ::boost::geometry::correct(voxelPolygon); BoostPolygonVector::const_iterator it; for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) { BoostPolygon2D contour = *it; ::boost::geometry::correct(contour); BoostPolygonDeque intersection; ::boost::geometry::intersection(voxelPolygon, contour, intersection); polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end()); } return polygonDeque; } BoostMaskVoxelizationThread::BoostRing2D BoostMaskVoxelizationThread::get2DContour( const rttb::VoxelGridIndex3D& aVoxelGrid3D) { BoostRing2D polygon; BoostPoint2D point1(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] - 0.5); ::boost::geometry::append(polygon, point1); BoostPoint2D point2(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] - 0.5); ::boost::geometry::append(polygon, point2); BoostPoint2D point3(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] + 0.5); ::boost::geometry::append(polygon, point3); BoostPoint2D point4(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] + 0.5); ::boost::geometry::append(polygon, point4); ::boost::geometry::append(polygon, point1); return polygon; } /*Calculate the intersection area*/ double BoostMaskVoxelizationThread::calcArea(const BoostPolygonDeque& aPolygonDeque) { double area = 0; BoostPolygonDeque::const_iterator it; for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) { area += ::boost::geometry::area(*it); } return area; } } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskVoxelizationThread.h b/code/masks/boost/rttbBoostMaskVoxelizationThread.h index fda6ff4..3ae8bcf 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.h +++ b/code/masks/boost/rttbBoostMaskVoxelizationThread.h @@ -1,114 +1,103 @@ // ----------------------------------------------------------------------- // 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 "rttbGeometricInfo.h" -#include "rttbMaskVoxel.h" -#include "rttbInvalidParameterException.h" #include #include -#include -#include +#include #include #include -#include -#include -#include -#include -#include namespace rttb { namespace masks { namespace boost { /*! @class BoostMaskGenerateMaskVoxelListThread * */ class BoostMaskVoxelizationThread { public: typedef ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy > BoostPolygon2D; typedef std::vector BoostPolygonVector;//polygon with or without holes typedef std::map BoostPolygonMap;//map of the z index with the vector of boost 2d polygon typedef std::vector VoxelIndexVector; typedef ::boost::multi_array BoostArray2D; - typedef ::boost::shared_ptr<::boost::lockfree::queue> - VoxelizationQueuePointer; - typedef ::boost::shared_ptr<::boost::lockfree::queue> - VoxelizationIndexQueuePointer; + typedef ::boost::shared_ptr BoostArray2DPointer; + typedef ::boost::shared_ptr > BoostArrayMapPointer; BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, - const VoxelIndexVector& aGlobalBoundingBox, VoxelizationIndexQueuePointer aResultIndexQueue, - VoxelizationQueuePointer aVoxelizationQueue); - + const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex); void operator()(); private: typedef std::deque BoostPolygonDeque; typedef ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy > BoostRing2D; typedef ::boost::geometry::model::d2::point_xy BoostPoint2D; BoostPolygonMap _geometryCoordinateBoostPolygonMap; VoxelIndexVector _globalBoundingBox; - VoxelizationQueuePointer _resultVoxelizationQueue; - VoxelizationIndexQueuePointer _resultIndexQueue; + BoostArrayMapPointer _resultVoxelization; + ::boost::shared_ptr<::boost::shared_mutex> _mutex; /*! @brief Get intersection polygons of the contour and a voxel polygon * @param aVoxelIndex3D The 3d grid index of the voxel * @param intersectionSlicePolygons The polygons of the slice intersecting the voxel - * @return Return all intersetion polygons of the structure and the voxel + * @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); }; } } } #endif \ No newline at end of file diff --git a/testing/examples/VoxelizationValidationTest.cpp b/testing/examples/VoxelizationValidationTest.cpp index e27ecc4..9365973 100644 --- a/testing/examples/VoxelizationValidationTest.cpp +++ b/testing/examples/VoxelizationValidationTest.cpp @@ -1,298 +1,301 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVHCalculator.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbNullPointerException.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbOTBMaskAccessor.h" #include "rttbDVHTxtFileReader.h" #include "rttbBoostMaskAccessor.h" #include "rttbITKImageMaskAccessorConverter.h" #include "rttbImageWriter.h" #include "rttbITKImageAccessorGenerator.h" #include "rttbITKImageFileAccessorGenerator.h" #include "rttbInvalidParameterException.h" #include "rttbBoostMask_LEGACY.h" #include "rttbBoostMaskAccessor_LEGACY.h" namespace rttb { namespace testing { /*! @brief VoxelizationValidationTest. Compare three differnt voxelizations: OTB, BoostLegacy and Boost. Check dvh maximum and minimum for each structure. Check write mask to itk file for further validation. */ int VoxelizationValidationTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: file name of test structure 1 // 2: file name of test dose 1 // 3: directory name to write boost mask of all structures // 4: directory name to write OTB mask of all structures // 5: file name of test dose 2: data with different z spacing of dose and structure, BooseLegacy throws exceptions by handling this data // 6: file name of test structure 2: data with different z spacing of dose and structure, BooseLegacy throws exceptions by handling this data std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string BoostMaskLegacy_DIRNAME; std::string OTBMask_DIRNAME; std::string BoostMask_DIRNAME; std::string RTDOSE2_FILENAME; std::string RTSTRUCT2_FILENAME; if (argc > 4) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; BoostMaskLegacy_DIRNAME = argv[3]; OTBMask_DIRNAME = argv[4]; BoostMask_DIRNAME = argv[5]; RTDOSE2_FILENAME = argv[6]; RTSTRUCT2_FILENAME = argv[7]; } ::boost::filesystem::create_directories(BoostMask_DIRNAME); ::boost::filesystem::create_directories(OTBMask_DIRNAME); ::boost::filesystem::create_directories(BoostMaskLegacy_DIRNAME); OFCondition status; DcmFileFormat fileformat; /* read dicom-rt dose */ io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create a vector of MaskAccessors (one for each structure) StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); + ::boost::filesystem::create_directories(BoostMask_DIRNAME); + ::boost::filesystem::create_directories(OTBMask_DIRNAME); + ::boost::filesystem::create_directories(BoostMaskRedesign_DIRNAME); if (rtStructureSet->getNumberOfStructures() > 0) { for (size_t j = 0; j < rtStructureSet->getNumberOfStructures(); j++) { std::cout << j << ": " << rtStructureSet->getStructure(j)->getLabel() << std::endl; clock_t start(clock()); //create OTB MaskAccessor ::boost::shared_ptr spOTBMaskAccessor = ::boost::make_shared(rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); spOTBMaskAccessor->updateMask(); MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor); ::boost::shared_ptr spMaskedDoseIteratorTmp = ::boost::make_shared(spMaskAccessor, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvh = *(calc.generateDVH()); clock_t finish(clock()); std::cout << "OTB Mask Calculation time: " << finish - start << " ms" << std::endl; //Write the mask image to a file. /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! */ rttb::io::itk::ITKImageMaskAccessorConverter itkConverter(spOTBMaskAccessor); CHECK(itkConverter.process()); std::stringstream fileNameSstr; fileNameSstr << OTBMask_DIRNAME << j << ".mhd"; rttb::io::itk::ImageWriter writer(fileNameSstr.str(), itkConverter.getITKImage()); CHECK(writer.writeFile()); clock_t start2(clock()); //create Boost MaskAccessor_LEGACY MaskAccessorPointer boostMaskAccessorPtr = ::boost::make_shared (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); CHECK_NO_THROW(boostMaskAccessorPtr->updateMask()); ::boost::shared_ptr spMaskedDoseIteratorTmp2 = ::boost::make_shared(boostMaskAccessorPtr, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator2(spMaskedDoseIteratorTmp2); rttb::core::DVHCalculator calc2(spMaskedDoseIterator2, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvh2 = *(calc2.generateDVH()); clock_t finish2(clock()); std::cout << "Boost Mask Calculation and write file time: " << finish2 - start2 << " ms" << std::endl; //Write the mask image to a file. /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! */ rttb::io::itk::ITKImageMaskAccessorConverter itkConverter2(boostMaskAccessorPtr); CHECK(itkConverter2.process()); std::stringstream fileNameSstr2; fileNameSstr2 << BoostMaskLegacy_DIRNAME << j << ".mhd"; rttb::io::itk::ImageWriter writer2(fileNameSstr2.str(), itkConverter2.getITKImage()); CHECK(writer2.writeFile()); //create Boost MaskAccessor redesign clock_t startR(clock()); MaskAccessorPointer boostMaskRPtr = ::boost::make_shared (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); CHECK_NO_THROW(boostMaskRPtr->updateMask()); ::boost::shared_ptr spMaskedDoseIteratorTmpR = ::boost::make_shared(boostMaskRPtr, doseAccessor1); DoseIteratorPointer spMaskedDoseIteratorR(spMaskedDoseIteratorTmpR); rttb::core::DVHCalculator calcR(spMaskedDoseIteratorR, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvhR = *(calcR.generateDVH()); clock_t finishR(clock()); std::cout << "Boost Mask Redesign Calculation and write file time: " << finishR - startR << " ms" << std::endl; //Write the mask image to a file. /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! */ rttb::io::itk::ITKImageMaskAccessorConverter itkConverterR(boostMaskRPtr); CHECK(itkConverterR.process()); std::stringstream fileNameSstrR; fileNameSstrR << BoostMask_DIRNAME << j << ".mhd"; rttb::io::itk::ImageWriter writerR(fileNameSstrR.str(), itkConverterR.getITKImage()); CHECK(writerR.writeFile()); //check close of 2 voxelizatin: OTB and Boost CHECK_CLOSE(dvh.getMaximum(), dvh2.getMaximum(), 0.1); CHECK_CLOSE(dvh.getMinimum(), dvh2.getMinimum(), 0.1); if (j != 7) { CHECK_CLOSE(dvh.getMean(), dvh2.getMean(), 0.1); } CHECK_CLOSE(dvh.getMedian(), dvh2.getMedian(), 0.1); CHECK_CLOSE(dvh.getModal(), dvh2.getModal(), 0.1); //0: Aussenkontur and 3: Niere li. failed. if (j != 0 && j != 3) { CHECK_CLOSE(dvh.getVx(0), dvh2.getVx(0), dvh.getVx(0) * 0.05); //check volume difference < 5% } //check close of 2 voxelization: Boost and BoostRedesign CHECK_CLOSE(dvhR.getMaximum(), dvh2.getMaximum(), 0.1); CHECK_CLOSE(dvhR.getMinimum(), dvh2.getMinimum(), 0.1); if (j != 7) { CHECK_CLOSE(dvhR.getMean(), dvh2.getMean(), 0.1); } CHECK_CLOSE(dvhR.getMedian(), dvh2.getMedian(), 0.1); CHECK_CLOSE(dvhR.getModal(), dvh2.getModal(), 0.1); //0: Aussenkontur and 3: Niere li. failed. CHECK_CLOSE(dvhR.getVx(0), dvh2.getVx(0), dvhR.getVx(0) * 0.05); //check volume difference < 5% } } /* Exception tests using data with different z spacing of dose and structure */ io::itk::ITKImageFileAccessorGenerator doseAccessorGenerator2(RTDOSE2_FILENAME.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); StructureSetPointer rtStructureSet2 = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT2_FILENAME.c_str()).generateStructureSet(); if (rtStructureSet2->getNumberOfStructures() > 0) { for (size_t j = 1; j < 20; j++) { std::cout << j << ": " << rtStructureSet2->getStructure(j)->getLabel() << std::endl; clock_t start(clock()); //create Boost MaskAccessor Legacy MaskAccessorPointer boostMaskAccessorPtr = ::boost::make_shared (rtStructureSet2->getStructure(j), doseAccessor2->getGeometricInfo()); //Two polygons in the same slice exception using boost mask, because of the different z spacing of dose and structure if (j != 9) { CHECK_THROW_EXPLICIT(boostMaskAccessorPtr->updateMask(), rttb::core::InvalidParameterException); } //create Boost MaskAccessor redesign MaskAccessorPointer boostMaskRPtr = ::boost::make_shared (rtStructureSet2->getStructure(j), doseAccessor2->getGeometricInfo()); //No exception using redesigned boost mask CHECK_NO_THROW(boostMaskRPtr->updateMask()); } } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/validation/CMakeLists.txt b/testing/validation/CMakeLists.txt index 400e75e..a117647 100644 --- a/testing/validation/CMakeLists.txt +++ b/testing/validation/CMakeLists.txt @@ -1,28 +1,28 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(CORE_TEST_VALIDATION ${EXECUTABLE_OUTPUT_PATH}/rttbValidationTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/Temporary) #----------------------------------------------------------------------------- ADD_TEST(VoxelizationDVHComparisonTest ${CORE_TEST_VALIDATION} VoxelizationDVHComparisonTest "${TEST_DATA_ROOT}/StructureSet/DICOM/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/Dose/DICOM/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/DVH/XML/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/OTBMask/" "${TEST_DATA_ROOT}/DVH/XML/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/BoostMask/" -"${TEST_DATA_ROOT}/DVH/XML/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/BoostMaskRedesign/") +"${RTTBTesting_BINARY_DIR}/validation/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/BoostMaskRedesign/") ADD_TEST(VoxelizationValidationTest ${CORE_TEST_VALIDATION} VoxelizationValidationTest "${TEST_DATA_ROOT}/StructureSet/DICOM/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/Dose/DICOM/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/Mask/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/BoostMask/" -"${TEST_DATA_ROOT}/Mask/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/OTBMask/" "${TEST_DATA_ROOT}/Mask/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/BoostMaskRedesign/" +"${TEST_DATA_ROOT}/Mask/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/OTBMask/" "${RTTBTesting_BINARY_DIR}/validation/1.3.6.1.4.1.2452.6.2037938358.1234393433.864109958.30410275/BoostMaskRedesign/" "${TEST_DATA_ROOT}/Dose/DICOM/PatBM116__RTDOSE_Main_corr.mhd" "${TEST_DATA_ROOT}/StructureSet/DICOM/__1__mm__1.2.276.0.28.19.640188827224102750711261838091512821903075578763.4195312189") RTTB_CREATE_TEST_MODULE(rttbValidation DEPENDS RTTBCore RTTBBoostMask RTTBDicomIO RTTBITKIO RTTBOtherIO PACKAGE_DEPENDS Litmus ITK) diff --git a/testing/validation/VoxelizationDVHComparisonTest.cpp b/testing/validation/VoxelizationDVHComparisonTest.cpp index c2b7c38..e8e5e28 100644 --- a/testing/validation/VoxelizationDVHComparisonTest.cpp +++ b/testing/validation/VoxelizationDVHComparisonTest.cpp @@ -1,178 +1,181 @@ // ----------------------------------------------------------------------- // 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: 1333 $ (last changed revision) // @date $Date: 2016-04-22 11:12:14 +0200 (Fr, 22 Apr 2016) $ (last change date) // @author $Author: hentsch $ (last changed by) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include #include "litCheckMacros.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbDVHCalculator.h" #include "rttbBoostMaskAccessor.h" #include "rttbDVHXMLFileWriter.h" #include "rttbDVHXMLFileReader.h" #include "../io/other/CompareDVH.h" namespace rttb { namespace testing { /*! @brief VoxelizationDVHComparisonTests. Computes the DVH difference of different voxelizations (OTB/Boost legacy/Boost) Writes the difference in a DVH and saves in a file */ core::DoseIteratorInterface::DoseIteratorPointer createMaskDoseIterator(masks::boost::BoostMaskAccessor::StructTypePointer rtstruct, core::GenericDoseIterator::DoseAccessorPointer doseAccessor, const std::string& voxelizationType) { core::GenericMaskedDoseIterator::MaskAccessorPointer spMaskAccessor; if (voxelizationType == "BoostRedesign"){ auto spBoostRedesignMaskAccessor = ::boost::make_shared(rtstruct, doseAccessor->getGeometricInfo()); spBoostRedesignMaskAccessor->updateMask(); spMaskAccessor = spBoostRedesignMaskAccessor; } auto spMaskedDoseIteratorTmp = ::boost::make_shared(spMaskAccessor, doseAccessor); core::DoseIteratorInterface::DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); return spMaskedDoseIterator; } rttb::core::DVH calcDVH(core::DVHCalculator::DoseIteratorPointer doseIterator, const IDType& structUID, const IDType& doseUID) { rttb::core::DVHCalculator calc(doseIterator, structUID, doseUID); rttb::core::DVH dvh = *(calc.generateDVH()); return dvh; } void writeCumulativeDVH(const std::string& filename, rttb::core::DVH dvh) { DVHType typeCum = { DVHType::Cumulative }; io::other::DVHXMLFileWriter dvhWriter(filename, typeCum); dvhWriter.writeDVH(boost::make_shared(dvh)); } int VoxelizationDVHComparisonTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::DVHCalculator::MaskedDoseIteratorPointer MaskedDoseIteratorPointer; typedef masks::boost::BoostMaskAccessor::StructTypePointer StructTypePointer; typedef core::DVH::DVHPointer DVHPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; PREPARE_DEFAULT_TEST_REPORTING; std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string RTDVH_XML_OTB_DIRECTORY; std::string RTDVH_XML_BOOST_DIRECTORY; std::string RTDVH_XML_BOOSTREDESIGN_DIRECTORY; if (argc > 5) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; RTDVH_XML_OTB_DIRECTORY = argv[3]; RTDVH_XML_BOOST_DIRECTORY = argv[4]; RTDVH_XML_BOOSTREDESIGN_DIRECTORY = argv[5]; } // read dicom-rt dose io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create a vector of MaskAccessors (one for each structure) StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); + //create directory + boost::filesystem::create_directories(RTDVH_XML_BOOSTREDESIGN_DIRECTORY); + //start evaluation clock_t start(clock()); if (rtStructureSet->getNumberOfStructures() > 0) { for (int j = 0; j < static_cast(rtStructureSet->getNumberOfStructures()); j++) { std::cout << rtStructureSet->getStructure(j)->getLabel() << std::endl; auto spMaskedDoseIteratorBoostRedesign = createMaskDoseIterator(rtStructureSet->getStructure(j), doseAccessor1, "BoostRedesign"); auto label = rtStructureSet->getStructure(j)->getLabel(); ::boost::replace_all(label, "/", "_"); boost::filesystem::path dvhOTBFilename(RTDVH_XML_OTB_DIRECTORY); dvhOTBFilename /= "DVH_" + label + ".xml"; boost::filesystem::path dvhBoostFilename(RTDVH_XML_BOOST_DIRECTORY); dvhBoostFilename /= "DVH_" + label + ".xml"; io::other::DVHXMLFileReader dvhReaderOTB(dvhOTBFilename.string()); auto dvhOTB = dvhReaderOTB.generateDVH(); io::other::DVHXMLFileReader dvhReaderBoost(dvhBoostFilename.string()); auto dvhBoost = dvhReaderBoost.generateDVH(); auto dvhBoostRedesign = calcDVH(spMaskedDoseIteratorBoostRedesign, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); boost::filesystem::path dvhBoostRedesignFilename(RTDVH_XML_BOOSTREDESIGN_DIRECTORY); dvhBoostRedesignFilename /= "DVH_" + label + ".xml"; writeCumulativeDVH(dvhBoostRedesignFilename.string(), dvhBoostRedesign); std::cout << "=== Dose 1 Structure " << j << "===" << std::endl; std::cout << "with OTB voxelization: " << std::endl; std::cout << dvhOTB << std::endl; std::cout << "with Boost voxelization: " << std::endl; std::cout << dvhBoost << std::endl; std::cout << "with BoostRedesign voxelization: " << std::endl; std::cout << dvhBoostRedesign << std::endl; //compare DVH for different voxelizations auto diffDVH = computeDiffDVH(dvhOTB, boost::make_shared(dvhBoostRedesign)); boost::filesystem::path dvhBoostRedesignDiffFilename(RTDVH_XML_BOOSTREDESIGN_DIRECTORY); dvhBoostRedesignDiffFilename /= "DVHDiff_" + label + ".xml"; writeCumulativeDVH(dvhBoostRedesignDiffFilename.string(), *diffDVH); } } clock_t finish(clock()); std::cout << "DVH Calculation time: " << finish - start << " ms" << std::endl; RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/validation/VoxelizationValidationTest.cpp b/testing/validation/VoxelizationValidationTest.cpp index f93aacd..f27e268 100644 --- a/testing/validation/VoxelizationValidationTest.cpp +++ b/testing/validation/VoxelizationValidationTest.cpp @@ -1,191 +1,194 @@ // ----------------------------------------------------------------------- // 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: 1495 $ (last changed revision) // @date $Date: 2016-09-29 16:17:47 +0200 (Do, 29 Sep 2016) $ (last change date) // @author $Author: hentsch $ (last changed by) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbBoostMaskAccessor.h" #include "rttbITKImageMaskAccessorConverter.h" #include "rttbImageWriter.h" #include "rttbBoostMask.h" #include "rttbBoostMaskAccessor.h" #include "rttbITKImageAccessorGenerator.h" #include "rttbITKImageFileAccessorGenerator.h" #include "rttbInvalidParameterException.h" #include "itkSubtractImageFilter.h" #include "itkImageFileReader.h" namespace rttb { namespace testing { io::itk::ITKImageMaskAccessor::ITKMaskImageType::Pointer subtractImages(const io::itk::ITKImageMaskAccessor::ITKMaskImageType::Pointer image1, const io::itk::ITKImageMaskAccessor::ITKMaskImageType::Pointer image2) { typedef itk::SubtractImageFilter SubtractImageFilterType; SubtractImageFilterType::Pointer subtractFilter = SubtractImageFilterType::New(); subtractFilter->SetInput1(image1); subtractFilter->SetInput2(image2); //since the image origin may be modified through the writing process a bit subtractFilter->SetCoordinateTolerance(5e-2); subtractFilter->Update(); return subtractFilter->GetOutput(); } /*! @brief VoxelizationValidationTest. Compare the new boost voxelization to the OTB voxelization Check the creating of new boost masks for files where the old boost voxelization failed. */ int VoxelizationValidationTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string BoostMask_DIRNAME; std::string OTBMask_DIRNAME; std::string BoostMaskRedesign_DIRNAME; std::string RTDose_BoostRedesign; std::string RTStr_BoostRedesign; if (argc > 7) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; BoostMask_DIRNAME = argv[3]; OTBMask_DIRNAME = argv[4]; BoostMaskRedesign_DIRNAME = argv[5]; RTDose_BoostRedesign = argv[6]; RTStr_BoostRedesign = argv[7]; } - OFCondition status; - DcmFileFormat fileformat; + //create directory + boost::filesystem::create_directories(BoostMaskRedesign_DIRNAME); /* read dicom-rt dose */ io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create a vector of MaskAccessors (one for each structure) StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); if (rtStructureSet->getNumberOfStructures() > 0) { - //do not compute struct "Aussenkontur" since it is very large (15000 cm³) + //do not compute structure "Aussenkontur" since it is very large (15000 cm³) for (size_t j = 1; j < rtStructureSet->getNumberOfStructures(); j++) { std::cout << j << ": " << rtStructureSet->getStructure(j)->getLabel() << std::endl; //read OTB mask image boost::filesystem::path otbMaskFilename(OTBMask_DIRNAME); otbMaskFilename /= boost::lexical_cast(j)+".mhd"; typedef itk::ImageFileReader ReaderType; ReaderType::Pointer readerOTB = ReaderType::New(); readerOTB->SetFileName(otbMaskFilename.string()); readerOTB->Update(); const io::itk::ITKImageMaskAccessor::ITKMaskImageType::Pointer otbMaskImage = readerOTB->GetOutput(); //create Boost MaskAccessor redesign clock_t startR(clock()); MaskAccessorPointer boostMaskRPtr = ::boost::make_shared (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); CHECK_NO_THROW(boostMaskRPtr->updateMask()); clock_t finishR(clock()); std::cout << "Boost Mask Redesign Calculation: " << finishR - startR << " ms" << std::endl; rttb::io::itk::ITKImageMaskAccessorConverter itkConverterR(boostMaskRPtr); CHECK(itkConverterR.process()); boost::filesystem::path redesignFilename(BoostMaskRedesign_DIRNAME); redesignFilename /= boost::lexical_cast(j)+".nrrd"; rttb::io::itk::ImageWriter writerR(redesignFilename.string(), itkConverterR.getITKImage()); CHECK(writerR.writeFile()); auto subtractedRedesignImage = subtractImages(otbMaskImage, itkConverterR.getITKImage()); boost::filesystem::path subtractRedesignFilename(BoostMaskRedesign_DIRNAME); subtractRedesignFilename /= boost::lexical_cast(j) + "_subtracted.nrrd"; rttb::io::itk::ImageWriter writerRSubtracted(subtractRedesignFilename.string(), subtractedRedesignImage); CHECK(writerRSubtracted.writeFile()); } } - + /* Exception tests using data with different z spacing of dose and structure */ io::itk::ITKImageFileAccessorGenerator doseAccessorGenerator2(RTDose_BoostRedesign.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); StructureSetPointer rtStructureSet2 = io::dicom::DicomFileStructureSetGenerator( RTStr_BoostRedesign.c_str()).generateStructureSet(); if (rtStructureSet2->getNumberOfStructures() > 0) { - for (size_t j = 12; j < 26; j++) - { - std::cout << j << ": " << rtStructureSet2->getStructure(j)->getLabel() << std::endl; - - //create Boost MaskAccessor redesign - MaskAccessorPointer boostMaskAccessorRedesignPtr - = ::boost::make_shared - (rtStructureSet2->getStructure(j), doseAccessor2->getGeometricInfo()); - - //No exception using redesigned boost mask - CHECK_NO_THROW(boostMaskAccessorRedesignPtr->updateMask()); + for (size_t j = 12; j < 26; j++) + { + //do not compute the largest structure + if (j != 18){ + std::cout << j << ": " << rtStructureSet2->getStructure(j)->getLabel() << std::endl; + + //create Boost MaskAccessor redesign + MaskAccessorPointer boostMaskAccessorRedesignPtr + = ::boost::make_shared + (rtStructureSet2->getStructure(j), doseAccessor2->getGeometricInfo()); + + //No exception using redesigned boost mask + CHECK_NO_THROW(boostMaskAccessorRedesignPtr->updateMask()); + } } } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb