diff --git a/code/masks/boost/files.cmake b/code/masks/boost/files.cmake index 86bbcad..faa905d 100644 --- a/code/masks/boost/files.cmake +++ b/code/masks/boost/files.cmake @@ -1,17 +1,13 @@ SET(CPP_FILES rttbBoostMask.cpp rttbBoostMaskAccessor.cpp rttbBoostMaskGenerateMaskVoxelListThread.cpp - rttbBoostMaskRedesign.cpp - rttbBoostMaskRedesignAccessor.cpp rttbBoostMaskVoxelizationThread.cpp ) SET(H_FILES rttbBoostMask.h rttbBoostMaskAccessor.h rttbBoostMaskGenerateMaskVoxelListThread.h - rttbBoostMaskRedesign.h - rttbBoostMaskRedesignAccessor.h rttbBoostMaskVoxelizationThread.h ) diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp index 912c1fa..76bc0ca 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/boost/rttbBoostMask.cpp @@ -1,498 +1,580 @@ #include #include #include +#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) - : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), + BoostMask::StructPointer aStructure, bool strict, unsigned int numberOfThreads) + : _geometricInfo(aDoseGeoInfo), _structure(aStructure), + _strict(strict), _numberOfThreads(numberOfThreads), _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: 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; + } - void BoostMask::calcMask() + BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { - if (hasSelfIntersections()) + if (!_isUpToDate) { - if (_strict) - { - throw rttb::core::InvalidParameterException("Error: structure has self intersections!"); - } - else - { - std::cerr << _structure->getLabel() << - " has self intersections! It may cause errors in the voxelization results!" << std::endl; - } + calcMask(); } - std::vector - intersectionSlicePolygonsVector;//store the polygons of intersection slice of each index z + return _voxelInStructure; + } - unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); + void BoostMask::calcMask() + { + preprocessing(); + voxelization(); + generateMaskVoxelList(); + _isUpToDate = true; + } + + void BoostMask::preprocessing() + { + rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); - //For each dose slice, get intersection structure slice and the contours on the structure slice - for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) + //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) { - BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(indexZ); - BoostMask::BoostPolygonVector::iterator it; + PolygonType rttbPolygon = *it; + PolygonType geometryCoordinatePolygon; - for (it = boostPolygons.begin(); it != boostPolygons.end(); ++it) + //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)) { - ::boost::geometry::correct(*it); + throw rttb::core::Exception("TiltedMaskPlaneException"); } - intersectionSlicePolygonsVector.push_back(boostPolygons); + geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon); } - /* Use simple nearest neighbour interpolation (NNI) if dose and structure has different z spacing: - If dose slice (indexZ) has no intersection with structure slice, but the last and the next do, that means the dose - slice lies between two structure slices -> use the last slice intersection contours for this slice - */ - for (unsigned int indexZ = 1; indexZ < nSlices - 1; indexZ++) - { - if (intersectionSlicePolygonsVector.at(indexZ).size() == 0 - && intersectionSlicePolygonsVector.at(indexZ - 1).size() > 0 - && intersectionSlicePolygonsVector.at(indexZ + 1).size() > 0) - { - intersectionSlicePolygonsVector.at(indexZ) = intersectionSlicePolygonsVector.at(indexZ - 1); - } - } + 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); - for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) - { - BoostMask::BoostPolygonVector intersectionSlicePolygons = intersectionSlicePolygonsVector.at( - indexZ); + //convert rttb polygon sequence to a map of z index and a vector of boost ring 2d (without holes) + _ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); - //Get bounding box of this dose slice - VoxelIndexVector voxelList = getBoundingBox(indexZ, intersectionSlicePolygons); - rttb::VoxelGridIndex3D gridIndex3D0 = voxelList.at(0); - rttb::VoxelGridIndex3D gridIndex3D1 = voxelList.at(1); + } - for (unsigned int indexX = gridIndex3D0[0]; indexX <= gridIndex3D1[0]; indexX++) - { - for (unsigned int indexY = gridIndex3D0[1]; indexY <= gridIndex3D1[1]; indexY++) - { - rttb::VoxelGridIndex3D currentIndex; - currentIndex[0] = indexX; - currentIndex[1] = indexY; - currentIndex[2] = indexZ; + void BoostMask::voxelization() + { + BoostPolygonMap::iterator it; - rttb::VoxelGridID gridID; - _geometricInfo->convert(currentIndex, gridID); + if (_globalBoundingBox.size() < 2) + { + throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); + } - //Get intersection polygons of the dose voxel and the structure - BoostPolygonDeque polygons = getIntersections(currentIndex, intersectionSlicePolygons); + unsigned int ringMapSize = _ringMap.size(); - //Calc areas of all intersection polygons - double intersectionArea = calcArea(polygons); + BoostRingMap::iterator itMap; - double voxelSize = _geometricInfo->getSpacing()[0] * _geometricInfo->getSpacing()[1]; + unsigned int mapSizeInAThread = _ringMap.size() / _numberOfThreads; + unsigned int count = 0; + unsigned int countThread = 0; + BoostPolygonMap polygonMap; + std::vector polygonMapVector; - if (intersectionArea > 0) - { - double volumeFraction = intersectionArea / voxelSize; - - if (volumeFraction > 1 && (volumeFraction - 1) <= 0.0001) - { - volumeFraction = 1; - } + //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); - core::MaskVoxel maskVoxel = core::MaskVoxel(gridID, (volumeFraction)); - _voxelInStructure->push_back(maskVoxel);//push back the mask voxel in structure - } - } + if (count == mapSizeInAThread && countThread < (_numberOfThreads - 1)) + { + polygonMapVector.push_back(polygonMap); + polygonMap.clear(); + count = 0; + countThread++; } - } - sort(_voxelInStructure->begin(), _voxelInStructure->end()); + polygonMap.insert(std::pair((*itMap).first, + polygonVector)); + count++; - _isUpToDate = true; - } + } - bool BoostMask::hasSelfIntersections() - { - bool hasSelfIntersection = false; + polygonMapVector.push_back(polygonMap); //insert the last one - unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); + //generate voxelization map, multi-threading + ::boost::thread_group threads; + BoostMaskVoxelizationThread::VoxelizationQueuePointer resultQueue + = ::boost::make_shared<::boost::lockfree::queue> + (ringMapSize); + BoostMaskVoxelizationThread::VoxelizationIndexQueuePointer resultIndexQueue + = ::boost::make_shared<::boost::lockfree::queue> + (ringMapSize); - for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) + for (int i = 0; i < polygonMapVector.size(); ++i) { + BoostMaskVoxelizationThread t(polygonMapVector.at(i), _globalBoundingBox, resultIndexQueue, + resultQueue); + threads.create_thread(t); + } - rttb::VoxelGridIndex3D currentIndex(0, 0, indexZ); - - BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(currentIndex[2]); - - BoostMask::BoostPolygonVector::iterator it; - BoostMask::BoostPolygonVector::iterator it2; - - for (it = boostPolygons.begin(); it != boostPolygons.end(); ++it) - { - ::boost::geometry::correct(*it); - - //test if polygon has self intersection - if (::boost::geometry::detail::overlay::has_self_intersections(*it, false)) - { - hasSelfIntersection = true; - std::cerr << _structure->getLabel() << " has self intersection polygon! Slice: " << indexZ << ". " - << std::endl; - break; - } + threads.join_all(); - //test if the polygons on the same slice has intersection - for (it2 = boostPolygons.begin(); it2 != boostPolygons.end() && it2 != it; ++it2) - { - ::boost::geometry::correct(*it2); - BoostPolygonDeque intersection; - ::boost::geometry::intersection(*it, *it2, intersection); - if (intersection.size() > 0) - { - //if no donut - if (!(::boost::geometry::within(*it, *it2)) && !(::boost::geometry::within(*it2, *it))) - { - hasSelfIntersection = true; - std::cerr << _structure->getLabel() << ": Two polygons on the same slice has intersection! Slice: " - << indexZ << "." << std::endl; - break; - } - } - } - } + BoostArray2D* array; + double index; + //generate voxelization map using resultQueue and resultIndexQueue + while (resultIndexQueue->pop(index) && resultQueue->pop(array)) + { + _voxelizationMap.insert(std::pair(index, *array)); } - return hasSelfIntersection; } - - /*Get the 4 voxel index of the bounding box of the polygon in the slice with sliceNumber*/ - BoostMask::VoxelIndexVector BoostMask::getBoundingBox(unsigned int sliceNumber, - const BoostPolygonVector& intersectionSlicePolygons) + void BoostMask::generateMaskVoxelList() { - unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); + if (_globalBoundingBox.size() < 2) + { + throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); + } - if (sliceNumber < 0 || sliceNumber >= nSlices) + //check homogeneus of the voxelization plane (the contours plane) + if (!calcVoxelizationThickness(_voxelizationThickness)) { - throw rttb::core::InvalidParameterException(std::string("Error: slice number is invalid!")); + throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneus!"); } - rttb::VoxelGridIndex3D currentIndex(0, 0, sliceNumber); - double xMin = 0; - double yMin = 0; - double xMax = 0; - double yMax = 0; + ::boost::shared_ptr<::boost::lockfree::queue> resultQueue + = ::boost::make_shared<::boost::lockfree::queue> + (_geometricInfo->getNumberOfVoxels()); + - BoostPolygonVector::const_iterator it; + ::boost::thread_group threads; - for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) + unsigned int sliceNumberInAThread = _geometricInfo->getNumSlices() / _numberOfThreads; + + //generate mask voxel list, multi-threading + for (unsigned int i = 0; i < _numberOfThreads; ++i) { - ::boost::geometry::model::box box; - ::boost::geometry::envelope(*it, box); - BoostPoint2D minPoint = box.min_corner(); - BoostPoint2D maxPoint = box.max_corner(); + unsigned int beginSlice = i * sliceNumberInAThread; + unsigned int endSlice; - if (it == intersectionSlicePolygons.begin()) + if (i < _numberOfThreads - 1) + { + endSlice = (i + 1) * sliceNumberInAThread; + } + else { - xMin = minPoint.x(); - yMin = minPoint.y(); - xMax = maxPoint.x(); - yMax = maxPoint.y(); + endSlice = _geometricInfo->getNumSlices(); } - xMin = std::min(xMin, minPoint.x()); - yMin = std::min(yMin, minPoint.y()); - xMax = std::max(xMax, maxPoint.x()); - yMax = std::max(yMax, maxPoint.y()); + + BoostMaskGenerateMaskVoxelListThread t(_globalBoundingBox, _geometricInfo, _voxelizationMap, + _voxelizationThickness, beginSlice, endSlice, + resultQueue); + + threads.create_thread(t); + } - rttb::WorldCoordinate3D nullWorldCoord; - _geometricInfo->indexToWorldCoordinate(VoxelGridIndex3D(0, 0, 0), nullWorldCoord); - rttb::WorldCoordinate3D minWorldCoord(xMin, yMin, nullWorldCoord.z()); - rttb::WorldCoordinate3D maxWorldCoord(xMax, yMax, nullWorldCoord.z()); + threads.join_all(); - rttb::VoxelGridIndex3D index0; - rttb::VoxelGridIndex3D index1; - _geometricInfo->worldCoordinateToIndex(minWorldCoord, index0); - _geometricInfo->worldCoordinateToIndex(maxWorldCoord, index1); + core::MaskVoxel* voxel; - VoxelIndexVector voxelList; - voxelList.push_back(index0); - voxelList.push_back(index1); + //generate _voxelInStructure using resultQueue + while (resultQueue->pop(voxel)) + { + _voxelInStructure->push_back(*voxel);//push back the mask voxel in structure + } - return voxelList; } - - /*Get intersection polygons of the contour and a voxel polygon*/ - BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& - aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) + bool BoostMask::preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, + rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, + rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const { - BoostMask::BoostPolygonDeque polygonDeque; - - BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); - ::boost::geometry::correct(voxelPolygon); - BoostPolygonVector::const_iterator it; + double minZ = _geometricInfo->getNumSlices(); + double maxZ = 0.0; - for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) + for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { - BoostPolygon2D contour = *it; - ::boost::geometry::correct(contour); - BoostPolygonDeque intersection; - ::boost::geometry::intersection(voxelPolygon, contour, intersection); + rttb::WorldCoordinate3D worldCoordinatePoint = aRTTBPolygon.at(i); - polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end()); - } + //convert to geometry coordinate polygon + rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint; + _geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, + geometryCoordinatePoint); + geometryCoordinatePolygon.push_back(geometryCoordinatePoint); - return polygonDeque; - } + //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); + } - /*Calculate the intersection area*/ - double BoostMask::calcArea(const BoostPolygonDeque& aPolygonDeque) - { - double area = 0; + //min and max for y + if (geometryCoordinatePoint(1) < minimum(1)) + { + minimum(1) = geometryCoordinatePoint(1); + } - BoostPolygonDeque::const_iterator it; + if (geometryCoordinatePoint(1) > maximum(1)) + { + maximum(1) = geometryCoordinatePoint(1); + } - for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) - { - area += ::boost::geometry::area(*it); - } + //min and max for z + if (geometryCoordinatePoint(2) < minimum(2)) + { + minimum(2) = geometryCoordinatePoint(2); + } - return area; - } + if (geometryCoordinatePoint(2) > maximum(2)) + { + maximum(2) = geometryCoordinatePoint(2); + } - VoxelGridIndex3D BoostMask::getGridIndex3D(const core::MaskVoxel& aMaskVoxel) - { - rttb::VoxelGridIndex3D gridIndex3D; - _geometricInfo->convert(aMaskVoxel.getVoxelGridID(), gridIndex3D); - return gridIndex3D; - } + //check planar + if (geometryCoordinatePoint(2) < minZ) + { + minZ = geometryCoordinatePoint(2); + } + + if (geometryCoordinatePoint(2) > maxZ) + { + maxZ = geometryCoordinatePoint(2); + } - BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() - { - if (!_isUpToDate) - { - calcMask(); } - return _voxelInStructure; + return (abs(maxZ - minZ) <= aErrorConstant); } - BoostMask::BoostRing2D BoostMask::convert(const rttb::PolygonType& aRTTBPolygon) + + 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::BoostPolygonVector BoostMask::getIntersectionSlicePolygons( - const rttb::GridIndexType aVoxelGridIndexZ) + BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( + const rttb::PolygonSequenceType& aRTTBPolygonVector) { - BoostMask::BoostRingVector boostRingVector; + rttb::PolygonSequenceType::const_iterator it; + BoostMask::BoostRingMap aRingMap; - rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); - - //get the polygons in the slice and convert to boost rings - rttb::PolygonSequenceType::iterator it; - - for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) + for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it) { - PolygonType rttbPolygon = *it; + rttb::PolygonType rttbPolygon = *it; + double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon + bool isFirstZ = true; - if (!rttbPolygon.empty()) + if (!aRingMap.empty()) { - double polygonZCoor = rttbPolygon.at(0)[2]; - rttb::WorldCoordinate3D polygonPoint(0, 0, polygonZCoor); - rttb::VoxelGridIndex3D polygonPointIndex3D; - _geometricInfo->worldCoordinateToIndex(polygonPoint, polygonPointIndex3D); + BoostMask::BoostRingMap::iterator findIt = findNearestKey(aRingMap, zIndex, errorConstant); - //if the z - if (aVoxelGridIndexZ == polygonPointIndex3D[2]) + //if the z index is found (same slice), add the polygon to vector + if (findIt != aRingMap.end()) { - boostRingVector.push_back(convert(rttbPolygon)); + //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 checkDonutAndConvert(boostRingVector); + return aRingMap; } - BoostMask::BoostRing2D BoostMask::get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) + BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap& + aBoostRingMap, double aIndex, double aErrorConstant) { - BoostMask::BoostRing2D polygon; - rttb::WorldCoordinate3D rttbPoint; - _geometricInfo->indexToWorldCoordinate(aVoxelGrid3D, rttbPoint); - - BoostPoint2D point1(rttbPoint[0], rttbPoint[1]); - ::boost::geometry::append(polygon, point1); - - double xSpacing = _geometricInfo->getSpacing()[0]; - double ySpacing = _geometricInfo->getSpacing()[1]; + BoostMask::BoostRingMap::iterator find = aBoostRingMap.find(aIndex); - BoostPoint2D point2(rttbPoint[0] + xSpacing, rttbPoint[1]); - ::boost::geometry::append(polygon, point2); - - BoostPoint2D point3(rttbPoint[0] + xSpacing, rttbPoint[1] + ySpacing); - ::boost::geometry::append(polygon, point3); + //if find a key equivalent to aIndex, found + if (find != aBoostRingMap.end()) + { + return find; + } + else + { + BoostMask::BoostRingMap::iterator lowerBound = aBoostRingMap.lower_bound(aIndex); - BoostPoint2D point4(rttbPoint[0], rttbPoint[1] + ySpacing); - ::boost::geometry::append(polygon, point4); + //if all keys go before aIndex, check the last key + if (lowerBound == aBoostRingMap.end()) + { + lowerBound = --aBoostRingMap.end(); + } - ::boost::geometry::append(polygon, point1); + //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 - return polygon; + //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) + 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) + { + aThickness = 1; + return true; + } + + double thickness = 0; + + for (; + it != _voxelizationMap.end(), it2 != _voxelizationMap.end(); ++it, ++it2) + { + if (thickness == 0) + { + thickness = (*it2).first - (*it).first; + } + else + { + double curThickness = (*it2).first - (*it).first; + //if no homogeneous (leave out double imprecisions), return false + if (abs(thickness-curThickness)>errorConstant) + { + return false; + } + } + + } + + if (thickness != 0) + { + aThickness = thickness; + } + else + { + aThickness = 1; + } + + return true; + } } } -} \ No newline at end of file +} diff --git a/code/masks/boost/rttbBoostMask.h b/code/masks/boost/rttbBoostMask.h index af2ced7..b484539 100644 --- a/code/masks/boost/rttbBoostMask.h +++ b/code/masks/boost/rttbBoostMask.h @@ -1,156 +1,210 @@ // ----------------------------------------------------------------------- // 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 #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 * 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); + 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; 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; + + //The thickness of the voxelization plane (the contour plane), in double dose grid index + //(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2) + double _voxelizationThickness; + bool _strict; //vector of the MaskVoxel inside the structure MaskVoxelListPointer _voxelInStructure; /*! @brief If the mask is up to date */ bool _isUpToDate; + + /*! @brief The number of threads + */ + unsigned int _numberOfThreads; + /*! @brief Voxelization and generate mask */ void calcMask(); - /*! @brief Check if the structure has self intersections*/ - bool hasSelfIntersections(); - - /*! @brief Get the min/max voxel index of the bounding box of the polygon in the slice with sliceNumber - * @param sliceNumber slice number between 0 and _geometricInfo->getNumSlices() - * @param intersectionSlicePolygons Get the polygons intersecting the slice - * @exception InvalidParameterException thrown if sliceNumber < 0 or sliceNumber >= _geometricInfo->getNumSlices() - * @return Return the 4 voxel index of the bounding box + /*! @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. */ - VoxelIndexVector getBoundingBox(unsigned int sliceNumber, - const BoostPolygonVector& intersectionSlicePolygons); + void preprocessing(); - /*! @brief Get intersection polygons of the contour and a voxel polygon - * @param aVoxelIndex3D The 3d grid index of the voxel - * @param intersectionSlicePolygons The polygons of the slice intersecting the voxel - * @return Return all intersetion polygons of the structure and the voxel - */ - BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, - const BoostPolygonVector& intersectionSlicePolygons); + /*! @brief The voxelization step, wich computes the voxelization planes (in x/y) for all contours of an struct. - /*! @brief Calculate the area of all polygons - * @param aPolygonDeque The deque of polygons - * @return Return the area of all polygons + * For each contour (that is in the z-Range of the reference geometry) of the struct: + * 1) Allocate result array (voxelization plane) based on the bounding box (see Preprocessing Step 3) + * 2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry). + * 3) Add result Array (key is the z-Value of the contour) + */ + void voxelization(); + + /*! @final mask voxel Generation step which transfers the voxilization planes into the (z-)geometry of the reference geometry. + * It consists of following Sub steps : + * For all "slices" in the reference geometry : + * 1) generate weight vector for all voxelization planes for a given z - value of a slice + * Itterate over the bounding box of a struct.For each voxel : + * 2) Compute weighted sum of all voxelization planes(use weight vector, step 1) + * 2a) If sum > 0 : Add mask voxel for the current x / y(inner Loop) and z value(outer Loop). + * 3) return mask voxel list. */ - double calcArea(const BoostPolygonDeque& aPolygonDeque); + void generateMaskVoxelList(); - /*! @brief Get grid index of a mask voxel - * @param aMaskVoxel A mask voxel - * @return Return the 3d grid index of the mask voxel + /*! @brief Convert the rttb polygon with world corrdinate to the rttb polygon with double geometry coordinate, calculate the current min/max + * and check if the polygon is planar + * @param minimum the current global minimum + * @param maximum the current global maximum + * @return Return true if the polygon is planar, which means that the minminal and maximal z-coordinate of the polygon is not larger than a error constant */ - VoxelGridIndex3D getGridIndex3D(const core::MaskVoxel& aMaskVoxel); + bool preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, + rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, + rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const; - /*! @brief Convert RTTB polygon to boost polygon*/ - BoostRing2D convert(const rttb::PolygonType& aRTTBPolygon); + /*! @brief Convert a rttb 3d polygon to a 2d boost ring*/ + BoostRing2D convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const; - /*! @brief Get the intersection slice of the voxel, return the polygons (with (donut) or without holes) in the slice - * @param aVoxelGridIndexZ The z grid index (slice number) of the voxel - * @return Return the structure polygons intersecting the slice - */ - BoostPolygonVector getIntersectionSlicePolygons(const rttb::GridIndexType aVoxelGridIndexZ); + /*! @brief Convert a rttb 3d polygon to a map of z index with a vector of boost 2d ring, because of tilt check use the first z index of the polygon as the map key*/ + BoostRingMap convertRTTBPolygonSequenceToBoostRingMap(const rttb::PolygonSequenceType& + aRTTBPolygonVector); - /*! @brief Get the voxel 2d contour polygon*/ - BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D); + /*! @brief Find the key with error constant to aIndex + * @pre aBoostRingMap should not be empty + * @return Return aBoostRingMap.end() if the key is not found + */ + BoostMask::BoostRingMap::iterator findNearestKey(BoostMask::BoostRingMap& aBoostRingMap, + double aIndex, double aErrorConstant); /*! @brief If 2 rings in the vector build a donut, convert the 2 rings to a donut polygon, other rings unchanged*/ - BoostPolygonVector checkDonutAndConvert(const BoostRingVector& aRingVector); + 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 \ No newline at end of file +#endif diff --git a/code/masks/boost/rttbBoostMaskAccessor.h b/code/masks/boost/rttbBoostMaskAccessor.h index 36f462f..360891b 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/boost/rttbBoostMaskAccessor.h @@ -1,130 +1,129 @@ // ----------------------------------------------------------------------- // 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 #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 * 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 + class RTTBBoostMask_EXPORT BoostMaskAccessor : public core::MaskAccessorInterface { public: typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::Structure::StructTypePointer StructTypePointer; private: - StructTypePointer _spStructure; + 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 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 eca3ca5..4edbf29 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp +++ b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.cpp @@ -1,133 +1,133 @@ #include "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbInvalidParameterException.h" #include namespace rttb { namespace masks { - namespace boostRedesign + namespace boost { BoostMaskGenerateMaskVoxelListThread::BoostMaskGenerateMaskVoxelListThread( const VoxelIndexVector& aGlobalBoundingBox, GeometricInfoPointer aGeometricInfo, const BoostArrayMap& aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, MaskVoxelQueuePointer aResultMaskVoxelQueue) : _globalBoundingBox(aGlobalBoundingBox), _geometricInfo(aGeometricInfo), _voxelizationMap(aVoxelizationMap), _voxelizationThickness(aVoxelizationThickness), _beginSlice(aBeginSlice), _endSlice(aEndSlice), _resultMaskVoxelQueue(aResultMaskVoxelQueue) {} 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; 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; for (it = weightVectorForZ.begin(), itMap = _voxelizationMap.begin(); it != weightVectorForZ.end() && itMap != _voxelizationMap.end(); ++it, ++itMap) { double weight = (*it).second; BoostArray2D voxelizationArray = (*itMap).second; //calc sum of all voxelization plane, use weight volumeFraction += voxelizationArray[x][y] * weight; } if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } //insert mask voxel if volumeFraction > 0 if (volumeFraction > 0) { core::MaskVoxel* maskVoxelPtr = new core::MaskVoxel(gridID, volumeFraction); _resultMaskVoxelQueue->push(maskVoxelPtr);//push back the mask voxel in structure } } } } } void BoostMaskGenerateMaskVoxelListThread::calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const { BoostArrayMap::const_iterator it; double indexZMin = aIndexZ - 0.5; double indexZMax = aIndexZ + 0.5; for (it = _voxelizationMap.begin(); it != _voxelizationMap.end(); ++it) { double voxelizationPlaneIndexMin = (*it).first - 0.5 * _voxelizationThickness; double voxelizationPlaneIndexMax = (*it).first + 0.5 * _voxelizationThickness; double weight = 0; if ((voxelizationPlaneIndexMin < indexZMin) && (voxelizationPlaneIndexMax > indexZMin)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = voxelizationPlaneIndexMax - indexZMin; } else { weight = 1; } } else if ((voxelizationPlaneIndexMin >= indexZMin) && (voxelizationPlaneIndexMin < indexZMax)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = _voxelizationThickness; } else { weight = indexZMax - voxelizationPlaneIndexMin; } } weightVector.insert(std::pair((*it).first, weight)); } } } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h index 9172498..a311a4e 100644 --- a/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h +++ b/code/masks/boost/rttbBoostMaskGenerateMaskVoxelListThread.h @@ -1,88 +1,88 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 1127 $ (last changed revision) // @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H #define __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include #include #include namespace rttb { namespace masks { - namespace boostRedesign + 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 std::vector VoxelIndexVector; BoostMaskGenerateMaskVoxelListThread(const VoxelIndexVector& aGlobalBoundingBox, GeometricInfoPointer aGeometricInfo, const BoostArrayMap& aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, MaskVoxelQueuePointer aResultMaskVoxelQueue); void operator()(); private: VoxelIndexVector _globalBoundingBox; GeometricInfoPointer _geometricInfo; BoostArrayMap _voxelizationMap; //(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2) double _voxelizationThickness; unsigned int _beginSlice; /*! @brief _endSlice is the first index not to be processed (like a end iterator) */ unsigned int _endSlice; MaskVoxelQueuePointer _resultMaskVoxelQueue; /*! @brief For each dose grid index z, calculate the weight vector for each structure contour */ void calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const; }; } } } #endif \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskRedesign.cpp b/code/masks/boost/rttbBoostMaskRedesign.cpp deleted file mode 100644 index 88b60ec..0000000 --- a/code/masks/boost/rttbBoostMaskRedesign.cpp +++ /dev/null @@ -1,580 +0,0 @@ -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#include "rttbBoostMaskRedesign.h" -#include "rttbNullPointerException.h" -#include "rttbInvalidParameterException.h" -#include "rttbBoostMaskGenerateMaskVoxelListThread.h" -#include "rttbBoostMaskVoxelizationThread.h" - - - -namespace rttb -{ - namespace masks - { - namespace boostRedesign - { - - - BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, - BoostMask::StructPointer aStructure, bool strict, unsigned int numberOfThreads) - : _geometricInfo(aDoseGeoInfo), _structure(aStructure), - _strict(strict), _numberOfThreads(numberOfThreads), - _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: 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; - 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++; - - } - - polygonMapVector.push_back(polygonMap); //insert the last one - - //generate voxelization map, multi-threading - ::boost::thread_group threads; - BoostMaskVoxelizationThread::VoxelizationQueuePointer resultQueue - = ::boost::make_shared<::boost::lockfree::queue> - (ringMapSize); - BoostMaskVoxelizationThread::VoxelizationIndexQueuePointer resultIndexQueue - = ::boost::make_shared<::boost::lockfree::queue> - (ringMapSize); - - for (int i = 0; i < polygonMapVector.size(); ++i) - { - BoostMaskVoxelizationThread t(polygonMapVector.at(i), _globalBoundingBox, resultIndexQueue, - resultQueue); - threads.create_thread(t); - } - - threads.join_all(); - - - BoostArray2D* array; - double index; - - //generate voxelization map using resultQueue and resultIndexQueue - while (resultIndexQueue->pop(index) && resultQueue->pop(array)) - { - _voxelizationMap.insert(std::pair(index, *array)); - } - - } - - void BoostMask::generateMaskVoxelList() - { - if (_globalBoundingBox.size() < 2) - { - throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); - } - - //check homogeneus of the voxelization plane (the contours plane) - if (!calcVoxelizationThickness(_voxelizationThickness)) - { - throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneus!"); - } - - - ::boost::shared_ptr<::boost::lockfree::queue> resultQueue - = ::boost::make_shared<::boost::lockfree::queue> - (_geometricInfo->getNumberOfVoxels()); - - - ::boost::thread_group threads; - - 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); - - 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); - geometryCoordinatePolygon.push_back(geometryCoordinatePoint); - - //calculate the current global min/max - //min and max for x - if (geometryCoordinatePoint(0) < minimum(0)) - { - minimum(0) = geometryCoordinatePoint(0); - } - - if (geometryCoordinatePoint(0) > maximum(0)) - { - maximum(0) = geometryCoordinatePoint(0); - } - - //min and max for y - if (geometryCoordinatePoint(1) < minimum(1)) - { - minimum(1) = geometryCoordinatePoint(1); - } - - if (geometryCoordinatePoint(1) > maximum(1)) - { - maximum(1) = geometryCoordinatePoint(1); - } - - //min and max for z - if (geometryCoordinatePoint(2) < minimum(2)) - { - minimum(2) = geometryCoordinatePoint(2); - } - - if (geometryCoordinatePoint(2) > maximum(2)) - { - maximum(2) = geometryCoordinatePoint(2); - } - - //check planar - if (geometryCoordinatePoint(2) < minZ) - { - minZ = geometryCoordinatePoint(2); - } - - if (geometryCoordinatePoint(2) > maxZ) - { - maxZ = geometryCoordinatePoint(2); - } - - } - - return (abs(maxZ - minZ) <= aErrorConstant); - } - - - BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& - aRTTBPolygon) const - { - BoostMask::BoostRing2D polygon2D; - BoostPoint2D firstPoint; - - for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) - { - rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); - BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); - - if (i == 0) - { - firstPoint = boostPoint; - } - - ::boost::geometry::append(polygon2D, boostPoint); - } - - ::boost::geometry::append(polygon2D, firstPoint); - return polygon2D; - } - - BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( - const rttb::PolygonSequenceType& aRTTBPolygonVector) - { - rttb::PolygonSequenceType::const_iterator it; - BoostMask::BoostRingMap aRingMap; - - for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it) - { - rttb::PolygonType rttbPolygon = *it; - double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon - bool isFirstZ = true; - - if (!aRingMap.empty()) - { - BoostMask::BoostRingMap::iterator findIt = findNearestKey(aRingMap, zIndex, errorConstant); - - //if the z index is found (same slice), add the polygon to vector - if (findIt != aRingMap.end()) - { - //BoostRingVector ringVector = ; - (*findIt).second.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); - isFirstZ = false; - } - } - - //if it is the first z index in the map, insert vector with the polygon - if (isFirstZ) - { - BoostRingVector ringVector; - ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); - aRingMap.insert(std::pair(zIndex, ringVector)); - } - - } - - return aRingMap; - } - - BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap& - aBoostRingMap, double aIndex, double aErrorConstant) - { - BoostMask::BoostRingMap::iterator find = aBoostRingMap.find(aIndex); - - //if find a key equivalent to aIndex, found - if (find != aBoostRingMap.end()) - { - return find; - } - else - { - BoostMask::BoostRingMap::iterator lowerBound = aBoostRingMap.lower_bound(aIndex); - - //if all keys go before aIndex, check the last key - if (lowerBound == aBoostRingMap.end()) - { - lowerBound = --aBoostRingMap.end(); - } - - //if the lower bound very close to aIndex, found - if (abs((*lowerBound).first - aIndex) <= aErrorConstant) - { - return lowerBound; - } - else - { - //if the lower bound is the beginning, not found - if (lowerBound == aBoostRingMap.begin()) - { - return aBoostRingMap.end(); - } - else - { - BoostMask::BoostRingMap::iterator lowerBound1 = --lowerBound;//the key before the lower bound - - //if the key before the lower bound very close to a Index, found - if (abs((*lowerBound1).first - aIndex) <= aErrorConstant) - { - return lowerBound1; - } - //else, not found - else - { - return aBoostRingMap.end(); - } - } - } - - } - } - - BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& - aRingVector) const - { - //check donut - BoostMask::BoostRingVector::const_iterator it1; - BoostMask::BoostRingVector::const_iterator it2; - BoostMask::BoostPolygonVector boostPolygonVector; - std::vector donutIndexVector;//store the outer and inner ring index - BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon - - //Get donut index and donut polygon - unsigned int index1 = 0; - - for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) - { - bool it1IsDonut = false; - - //check if the ring is already determined as a donut - for (unsigned int i = 0; i < donutIndexVector.size(); i++) - { - if (donutIndexVector.at(i) == index1) - { - it1IsDonut = true; - break; - } - } - - //if not jet, check now - if (!it1IsDonut) - { - bool it2IsDonut = false; - unsigned int index2 = 0; - - for (it2 = aRingVector.begin(); it2 != aRingVector.end(); it2++, index2++) - { - if (it2 != it1) - { - BoostMask::BoostPolygon2D polygon2D; - - if (::boost::geometry::within(*it1, *it2)) - { - ::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon - ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring - ::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring - it2IsDonut = true; - } - //if donut - else if (::boost::geometry::within(*it2, *it1)) - { - ::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon - ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring - ::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring - it2IsDonut = true; - } - - if (it2IsDonut) - { - donutIndexVector.push_back(index1); - donutIndexVector.push_back(index2); - donutVector.push_back(polygon2D);//store donut polygon - break;//Only store the first donut! - } - } - } - } - } - - //Store no donut polygon to boostPolygonVector - index1 = 0; - - for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) - { - bool it1IsDonut = false; - - //check if the ring is the outer or inner of a donut - for (unsigned int i = 0; i < donutIndexVector.size(); i++) - { - if (donutIndexVector.at(i) == index1) - { - it1IsDonut = true; - break; - } - } - - if (!it1IsDonut) - { - BoostMask::BoostPolygon2D polygon2D; - ::boost::geometry::append(polygon2D, *it1); - boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut - } - } - - //Append donut polygon to boostPolygonVector - BoostMask::BoostPolygonVector::iterator itDonut; - - for (itDonut = donutVector.begin(); itDonut != donutVector.end(); itDonut++) - { - boostPolygonVector.push_back(*itDonut);//append donuts - } - - return boostPolygonVector; - } - - bool BoostMask::calcVoxelizationThickness(double& aThickness) const - { - BoostArrayMap::const_iterator it = _voxelizationMap.begin(); - BoostArrayMap::const_iterator it2 = ++_voxelizationMap.begin(); - - if (_voxelizationMap.size() <= 1) - { - aThickness = 1; - return true; - } - - double thickness = 0; - - for (; - it != _voxelizationMap.end(), it2 != _voxelizationMap.end(); ++it, ++it2) - { - if (thickness == 0) - { - thickness = (*it2).first - (*it).first; - } - else - { - double curThickness = (*it2).first - (*it).first; - //if no 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/rttbBoostMaskRedesign.h b/code/masks/boost/rttbBoostMaskRedesign.h deleted file mode 100644 index 4844a4c..0000000 --- a/code/masks/boost/rttbBoostMaskRedesign.h +++ /dev/null @@ -1,210 +0,0 @@ -// ----------------------------------------------------------------------- -// RTToolbox - DKFZ radiotherapy quantitative evaluation library -// -// Copyright (c) German Cancer Research Center (DKFZ), -// Software development for Integrated Diagnostics and Therapy (SIDT). -// ALL RIGHTS RESERVED. -// See rttbCopyright.txt or -// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html -// -// This software is distributed WITHOUT ANY WARRANTY; without even -// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -// PURPOSE. See the above copyright notices for more information. -// -//------------------------------------------------------------------------ -/*! -// @file -// @version $Revision: 1127 $ (last changed revision) -// @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date) -// @author $Author: hentsch $ (last changed by) -*/ - - -#ifndef __BOOST_MASK_R_H -#define __BOOST_MASK_R_H - -#include "rttbBaseType.h" -#include "rttbStructure.h" -#include "rttbGeometricInfo.h" -#include "rttbMaskVoxel.h" -#include "rttbMaskAccessorInterface.h" - -#include -#include -#include -#include -#include -#include - -namespace rttb -{ - namespace masks - { - namespace boostRedesign - { - /*! @class BoostMask - * @brief Implementation of voxelization using boost::geometry. - * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. - * (A structure without self interseciton means all contours of the structure have no self intersection, and - * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) - * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored - * and the mask will be calculated, however, it may cause errors in the mask results. - */ - class BoostMask - { - - public: - typedef ::boost::shared_ptr GeometricInfoPointer; - typedef core::Structure::StructTypePointer StructPointer; - typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; - typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; - - /*! @brief Constructor - * @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is NULL - * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. - * @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; - - 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; - - //The thickness of the voxelization plane (the contour plane), in double dose grid index - //(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2) - double _voxelizationThickness; - - bool _strict; - - //vector of the MaskVoxel inside the structure - MaskVoxelListPointer _voxelInStructure; - - /*! @brief If the mask is up to date - */ - bool _isUpToDate; - - - /*! @brief The 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. - - * For each contour (that is in the z-Range of the reference geometry) of the struct: - * 1) Allocate result array (voxelization plane) based on the bounding box (see Preprocessing Step 3) - * 2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry). - * 3) Add result Array (key is the z-Value of the contour) - */ - void voxelization(); - - /*! @final mask voxel Generation step which transfers the voxilization planes into the (z-)geometry of the reference geometry. - * It consists of following Sub steps : - * For all "slices" in the reference geometry : - * 1) generate weight vector for all voxelization planes for a given z - value of a slice - * Itterate over the bounding box of a struct.For each voxel : - * 2) Compute weighted sum of all voxelization planes(use weight vector, step 1) - * 2a) If sum > 0 : Add mask voxel for the current x / y(inner Loop) and z value(outer Loop). - * 3) return mask voxel list. - */ - void generateMaskVoxelList(); - - /*! @brief Convert the rttb polygon with world corrdinate to the rttb polygon with double geometry coordinate, calculate the current min/max - * and check if the polygon is planar - * @param minimum the current global minimum - * @param maximum the current global maximum - * @return Return true if the polygon is planar, which means that the minminal and maximal z-coordinate of the polygon is not larger than a error constant - */ - bool preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, - rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, - rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const; - - /*! @brief Convert a rttb 3d polygon to a 2d boost ring*/ - BoostRing2D convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const; - - /*! @brief Convert a rttb 3d polygon to a map of z index with a vector of boost 2d ring, because of tilt check use the first z index of the polygon as the map key*/ - BoostRingMap convertRTTBPolygonSequenceToBoostRingMap(const rttb::PolygonSequenceType& - aRTTBPolygonVector); - - /*! @brief Find the key with error constant to aIndex - * @pre aBoostRingMap should not be empty - * @return Return aBoostRingMap.end() if the key is not found - */ - BoostMask::BoostRingMap::iterator findNearestKey(BoostMask::BoostRingMap& aBoostRingMap, - double aIndex, double aErrorConstant); - - /*! @brief If 2 rings in the vector build a donut, convert the 2 rings to a donut polygon, other rings unchanged*/ - BoostPolygonVector checkDonutAndConvert(const BoostRingVector& aRingVector) const; - - /*! @brief 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/rttbBoostMaskVoxelizationThread.cpp b/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp index 2d218b1..b402fa8 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp +++ b/code/masks/boost/rttbBoostMaskVoxelizationThread.cpp @@ -1,135 +1,135 @@ #include "rttbBoostMaskVoxelizationThread.h" namespace rttb { namespace masks { - namespace boostRedesign + namespace boost { BoostMaskVoxelizationThread::BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, const VoxelIndexVector& aGlobalBoundingBox, VoxelizationIndexQueuePointer aResultIndexQueue, VoxelizationQueuePointer aVoxelizationQueue) : _geometryCoordinateBoostPolygonMap(APolygonMap), _globalBoundingBox(aGlobalBoundingBox), _resultIndexQueue(aResultIndexQueue), _resultVoxelizationQueue(aVoxelizationQueue) { } void BoostMaskVoxelizationThread::operator()() { rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; BoostPolygonMap::iterator it; for (it = _geometryCoordinateBoostPolygonMap.begin(); it != _geometryCoordinateBoostPolygonMap.end(); ++it) { BoostArray2D* maskArray = new BoostArray2D( - boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); + ::boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); BoostPolygonVector boostPolygonVec = (*it).second; for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) { for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = x + minIndex[0]; currentIndex[1] = y + minIndex[1]; currentIndex[2] = 0; //Get intersection polygons of the dose voxel and the structure BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec); //Calc areas of all intersection polygons double volumeFraction = calcArea(polygons); if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } (*maskArray)[x][y] = volumeFraction; } } //insert into voxelization map _resultIndexQueue->push((*it).first); _resultVoxelizationQueue->push(maskArray); } } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMaskVoxelizationThread::BoostPolygonDeque BoostMaskVoxelizationThread::getIntersections( const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) { 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 1bba438..fda6ff4 100644 --- a/code/masks/boost/rttbBoostMaskVoxelizationThread.h +++ b/code/masks/boost/rttbBoostMaskVoxelizationThread.h @@ -1,114 +1,114 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 1127 $ (last changed revision) // @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __BOOST_MASK_VOXELIZATION_THREAD_H #define __BOOST_MASK_VOXELIZATION_THREAD_H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbInvalidParameterException.h" #include #include #include #include #include #include #include #include #include #include #include namespace rttb { namespace masks { - namespace boostRedesign + 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; BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, const VoxelIndexVector& aGlobalBoundingBox, VoxelizationIndexQueuePointer aResultIndexQueue, VoxelizationQueuePointer aVoxelizationQueue); 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; /*! @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 */ 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/code/masks/legacy/files.cmake b/code/masks/legacy/files.cmake index 51b0fbf..e589bcb 100644 --- a/code/masks/legacy/files.cmake +++ b/code/masks/legacy/files.cmake @@ -1,28 +1,33 @@ SET(CPP_FILES rttbOTBMaskAccessor.cpp rttbMask.cpp DoseVoxel_LEGACY.cpp DoseIteratorInterface_LEGACY.cpp DoseIterator_LEGACY.cpp rttbContour_LEGACY.cpp rttbIterativePolygonInTermsOfIntersections_LEGACY.cpp rttbPolygonInfo_LEGACY.cpp rttbSelfIntersectingStructureException.cpp rttbStructure_LEGACY.cpp + rttbBoostMask_LEGACY.cpp + rttbBoostMaskAccessor_LEGACY.cpp ) SET(H_FILES rttbMaskType_LEGACY.h rttbMask.h rttbOTBMaskAccessor.h DoseVoxel_LEGACY.h DoseIteratorInterface_LEGACY.h DoseIterator_LEGACY.h rttbContour_LEGACY.h rttbIterativePolygonInTermsOfIntersections_LEGACY.h rttbField_LEGACY.h rttbPolygonInfo_LEGACY.h rttbBaseType_LEGACY.h rttbSelfIntersectingStructureException.h rttbStructure_LEGACY.h + rttbBoostMask_LEGACY.h + rttbBoostMaskAccessor_LEGACY.h + ) diff --git a/code/masks/boost/rttbBoostMaskRedesignAccessor.cpp b/code/masks/legacy/rttbBoostMaskAccessor_LEGACY.cpp similarity index 89% rename from code/masks/boost/rttbBoostMaskRedesignAccessor.cpp rename to code/masks/legacy/rttbBoostMaskAccessor_LEGACY.cpp index bf8fd80..95bf36f 100644 --- a/code/masks/boost/rttbBoostMaskRedesignAccessor.cpp +++ b/code/masks/legacy/rttbBoostMaskAccessor_LEGACY.cpp @@ -1,164 +1,164 @@ // ----------------------------------------------------------------------- // 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) +// @version $Revision: 1265 $ (last changed revision) +// @date $Date: 2016-03-08 15:11:43 +0100 (Di, 08 Mrz 2016) $ (last change date) +// @author $Author: zhangl $ (last changed by) */ -#include "rttbBoostMaskRedesignAccessor.h" -#include "rttbBoostMaskRedesign.h" +#include "rttbBoostMaskAccessor_LEGACY.h" +#include "rttbBoostMask_LEGACY.h" #include #include #include #include namespace rttb { namespace masks { - namespace boostRedesign + namespace boostLegacy { BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, bool strict) : _spStructure(aStructurePointer), _geoInfo(aGeometricInfo), _strict(strict) { _spRelevantVoxelVector = MaskVoxelListPointer(); //generate new structure set uid ::boost::uuids::uuid id; ::boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _maskUID = "BoostMask_" + ss.str(); } BoostMaskAccessor::~BoostMaskAccessor() { }; void BoostMaskAccessor::updateMask() { MaskVoxelList newRelevantVoxelVector; if (_spRelevantVoxelVector) { return; // already calculated } BoostMask mask(::boost::make_shared(_geoInfo), _spStructure, _strict); _spRelevantVoxelVector = mask.getRelevantVoxelVector(); return; } BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector() { // if not already generated start voxelization here updateMask(); return _spRelevantVoxelVector; } BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector( float lowerThreshold) { MaskVoxelListPointer filteredVoxelVectorPointer(new MaskVoxelList); updateMask(); // filter relevant voxels BoostMaskAccessor::MaskVoxelList::iterator it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getRelevantVolumeFraction() > lowerThreshold) { filteredVoxelVectorPointer->push_back(*it); } ++it; } // if mask calculation was not successful this is empty! return filteredVoxelVectorPointer; } bool BoostMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const { //initialize return voxel voxel.setRelevantVolumeFraction(0); //check if ID is valid if (!_geoInfo.validID(aID)) { return false; } //determine how a given voxel on the dose grid is masked if (_spRelevantVoxelVector) { BoostMaskAccessor::MaskVoxelList::iterator it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getVoxelGridID() == aID) { voxel = (*it); return true; } ++it; } } // returns false if mask was not calculated without triggering calculation (otherwise not const!) else { return false; } return false; } bool BoostMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const { //convert VoxelGridIndex3D to VoxelGridID VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getMaskAt(aVoxelGridID, voxel); } else { return false; } } const core::GeometricInfo& BoostMaskAccessor::getGeometricInfo() const { return _geoInfo; }; } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskRedesignAccessor.h b/code/masks/legacy/rttbBoostMaskAccessor_LEGACY.h similarity index 91% rename from code/masks/boost/rttbBoostMaskRedesignAccessor.h rename to code/masks/legacy/rttbBoostMaskAccessor_LEGACY.h index 8268e29..fac733e 100644 --- a/code/masks/boost/rttbBoostMaskRedesignAccessor.h +++ b/code/masks/legacy/rttbBoostMaskAccessor_LEGACY.h @@ -1,129 +1,130 @@ // ----------------------------------------------------------------------- // 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) +// @version $Revision: 1364 $ (last changed revision) +// @date $Date: 2016-05-23 10:49:20 +0200 (Mo, 23 Mai 2016) $ (last change date) // @author $Author: hentsch $ (last changed by) */ -#ifndef __BOOST_MASK_R_ACCESSOR__H -#define __BOOST_MASK_R_ACCESSOR__H +#ifndef __BOOST_MASK_ACCESSOR_L__H +#define __BOOST_MASK_ACCESSOR_L__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 boostRedesign + namespace boostLegacy { /*! @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 * 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 + class RTTBBoostMask_EXPORT BoostMaskAccessor : public core::MaskAccessorInterface { public: typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::Structure::StructTypePointer StructTypePointer; private: - + StructTypePointer _spStructure; 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 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/rttbBoostMask.cpp b/code/masks/legacy/rttbBoostMask_LEGACY.cpp similarity index 96% copy from code/masks/boost/rttbBoostMask.cpp copy to code/masks/legacy/rttbBoostMask_LEGACY.cpp index 912c1fa..9473e06 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/legacy/rttbBoostMask_LEGACY.cpp @@ -1,498 +1,498 @@ #include #include #include #include #include #include #include #include -#include "rttbBoostMask.h" +#include "rttbBoostMask_LEGACY.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace masks { - namespace boost + namespace boostLegacy { BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, BoostMask::StructPointer aStructure, bool strict) : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _voxelInStructure(::boost::make_shared()) { _isUpToDate = false; if (!_geometricInfo) { throw rttb::core::NullPointerException("Error: Geometric info is NULL!"); } else if (!_structure) { throw rttb::core::NullPointerException("Error: Structure is NULL!"); } } void BoostMask::calcMask() { if (hasSelfIntersections()) { if (_strict) { throw rttb::core::InvalidParameterException("Error: structure has self intersections!"); } else { std::cerr << _structure->getLabel() << " has self intersections! It may cause errors in the voxelization results!" << std::endl; } } std::vector intersectionSlicePolygonsVector;//store the polygons of intersection slice of each index z unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); //For each dose slice, get intersection structure slice and the contours on the structure slice for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) { BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(indexZ); BoostMask::BoostPolygonVector::iterator it; for (it = boostPolygons.begin(); it != boostPolygons.end(); ++it) { ::boost::geometry::correct(*it); } intersectionSlicePolygonsVector.push_back(boostPolygons); } /* Use simple nearest neighbour interpolation (NNI) if dose and structure has different z spacing: If dose slice (indexZ) has no intersection with structure slice, but the last and the next do, that means the dose slice lies between two structure slices -> use the last slice intersection contours for this slice */ for (unsigned int indexZ = 1; indexZ < nSlices - 1; indexZ++) { if (intersectionSlicePolygonsVector.at(indexZ).size() == 0 && intersectionSlicePolygonsVector.at(indexZ - 1).size() > 0 && intersectionSlicePolygonsVector.at(indexZ + 1).size() > 0) { intersectionSlicePolygonsVector.at(indexZ) = intersectionSlicePolygonsVector.at(indexZ - 1); } } for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) { BoostMask::BoostPolygonVector intersectionSlicePolygons = intersectionSlicePolygonsVector.at( indexZ); //Get bounding box of this dose slice VoxelIndexVector voxelList = getBoundingBox(indexZ, intersectionSlicePolygons); rttb::VoxelGridIndex3D gridIndex3D0 = voxelList.at(0); rttb::VoxelGridIndex3D gridIndex3D1 = voxelList.at(1); for (unsigned int indexX = gridIndex3D0[0]; indexX <= gridIndex3D1[0]; indexX++) { for (unsigned int indexY = gridIndex3D0[1]; indexY <= gridIndex3D1[1]; indexY++) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = indexX; currentIndex[1] = indexY; currentIndex[2] = indexZ; rttb::VoxelGridID gridID; _geometricInfo->convert(currentIndex, gridID); //Get intersection polygons of the dose voxel and the structure BoostPolygonDeque polygons = getIntersections(currentIndex, intersectionSlicePolygons); //Calc areas of all intersection polygons double intersectionArea = calcArea(polygons); double voxelSize = _geometricInfo->getSpacing()[0] * _geometricInfo->getSpacing()[1]; if (intersectionArea > 0) { double volumeFraction = intersectionArea / voxelSize; if (volumeFraction > 1 && (volumeFraction - 1) <= 0.0001) { volumeFraction = 1; } core::MaskVoxel maskVoxel = core::MaskVoxel(gridID, (volumeFraction)); _voxelInStructure->push_back(maskVoxel);//push back the mask voxel in structure } } } } sort(_voxelInStructure->begin(), _voxelInStructure->end()); _isUpToDate = true; } bool BoostMask::hasSelfIntersections() { bool hasSelfIntersection = false; unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) { rttb::VoxelGridIndex3D currentIndex(0, 0, indexZ); BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(currentIndex[2]); BoostMask::BoostPolygonVector::iterator it; BoostMask::BoostPolygonVector::iterator it2; for (it = boostPolygons.begin(); it != boostPolygons.end(); ++it) { ::boost::geometry::correct(*it); //test if polygon has self intersection if (::boost::geometry::detail::overlay::has_self_intersections(*it, false)) { hasSelfIntersection = true; std::cerr << _structure->getLabel() << " has self intersection polygon! Slice: " << indexZ << ". " << std::endl; break; } //test if the polygons on the same slice has intersection for (it2 = boostPolygons.begin(); it2 != boostPolygons.end() && it2 != it; ++it2) { ::boost::geometry::correct(*it2); BoostPolygonDeque intersection; ::boost::geometry::intersection(*it, *it2, intersection); if (intersection.size() > 0) { //if no donut if (!(::boost::geometry::within(*it, *it2)) && !(::boost::geometry::within(*it2, *it))) { hasSelfIntersection = true; std::cerr << _structure->getLabel() << ": Two polygons on the same slice has intersection! Slice: " << indexZ << "." << std::endl; break; } } } } } return hasSelfIntersection; } /*Get the 4 voxel index of the bounding box of the polygon in the slice with sliceNumber*/ BoostMask::VoxelIndexVector BoostMask::getBoundingBox(unsigned int sliceNumber, const BoostPolygonVector& intersectionSlicePolygons) { unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); if (sliceNumber < 0 || sliceNumber >= nSlices) { throw rttb::core::InvalidParameterException(std::string("Error: slice number is invalid!")); } rttb::VoxelGridIndex3D currentIndex(0, 0, sliceNumber); double xMin = 0; double yMin = 0; double xMax = 0; double yMax = 0; BoostPolygonVector::const_iterator it; for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) { ::boost::geometry::model::box box; ::boost::geometry::envelope(*it, box); BoostPoint2D minPoint = box.min_corner(); BoostPoint2D maxPoint = box.max_corner(); if (it == intersectionSlicePolygons.begin()) { xMin = minPoint.x(); yMin = minPoint.y(); xMax = maxPoint.x(); yMax = maxPoint.y(); } xMin = std::min(xMin, minPoint.x()); yMin = std::min(yMin, minPoint.y()); xMax = std::max(xMax, maxPoint.x()); yMax = std::max(yMax, maxPoint.y()); } rttb::WorldCoordinate3D nullWorldCoord; _geometricInfo->indexToWorldCoordinate(VoxelGridIndex3D(0, 0, 0), nullWorldCoord); rttb::WorldCoordinate3D minWorldCoord(xMin, yMin, nullWorldCoord.z()); rttb::WorldCoordinate3D maxWorldCoord(xMax, yMax, nullWorldCoord.z()); rttb::VoxelGridIndex3D index0; rttb::VoxelGridIndex3D index1; _geometricInfo->worldCoordinateToIndex(minWorldCoord, index0); _geometricInfo->worldCoordinateToIndex(maxWorldCoord, index1); VoxelIndexVector voxelList; voxelList.push_back(index0); voxelList.push_back(index1); return voxelList; } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) { BoostMask::BoostPolygonDeque polygonDeque; BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); ::boost::geometry::correct(voxelPolygon); BoostPolygonVector::const_iterator it; for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) { BoostPolygon2D contour = *it; ::boost::geometry::correct(contour); BoostPolygonDeque intersection; ::boost::geometry::intersection(voxelPolygon, contour, intersection); polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end()); } return polygonDeque; } /*Calculate the intersection area*/ double BoostMask::calcArea(const BoostPolygonDeque& aPolygonDeque) { double area = 0; BoostPolygonDeque::const_iterator it; for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) { area += ::boost::geometry::area(*it); } return area; } VoxelGridIndex3D BoostMask::getGridIndex3D(const core::MaskVoxel& aMaskVoxel) { rttb::VoxelGridIndex3D gridIndex3D; _geometricInfo->convert(aMaskVoxel.getVoxelGridID(), gridIndex3D); return gridIndex3D; } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { if (!_isUpToDate) { calcMask(); } return _voxelInStructure; } BoostMask::BoostRing2D BoostMask::convert(const rttb::PolygonType& aRTTBPolygon) { 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::BoostPolygonVector BoostMask::getIntersectionSlicePolygons( const rttb::GridIndexType aVoxelGridIndexZ) { BoostMask::BoostRingVector boostRingVector; rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); //get the polygons in the slice and convert to boost rings rttb::PolygonSequenceType::iterator it; for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) { PolygonType rttbPolygon = *it; if (!rttbPolygon.empty()) { double polygonZCoor = rttbPolygon.at(0)[2]; rttb::WorldCoordinate3D polygonPoint(0, 0, polygonZCoor); rttb::VoxelGridIndex3D polygonPointIndex3D; _geometricInfo->worldCoordinateToIndex(polygonPoint, polygonPointIndex3D); //if the z if (aVoxelGridIndexZ == polygonPointIndex3D[2]) { boostRingVector.push_back(convert(rttbPolygon)); } } } return checkDonutAndConvert(boostRingVector); } BoostMask::BoostRing2D BoostMask::get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) { BoostMask::BoostRing2D polygon; rttb::WorldCoordinate3D rttbPoint; _geometricInfo->indexToWorldCoordinate(aVoxelGrid3D, rttbPoint); BoostPoint2D point1(rttbPoint[0], rttbPoint[1]); ::boost::geometry::append(polygon, point1); double xSpacing = _geometricInfo->getSpacing()[0]; double ySpacing = _geometricInfo->getSpacing()[1]; BoostPoint2D point2(rttbPoint[0] + xSpacing, rttbPoint[1]); ::boost::geometry::append(polygon, point2); BoostPoint2D point3(rttbPoint[0] + xSpacing, rttbPoint[1] + ySpacing); ::boost::geometry::append(polygon, point3); BoostPoint2D point4(rttbPoint[0], rttbPoint[1] + ySpacing); ::boost::geometry::append(polygon, point4); ::boost::geometry::append(polygon, point1); return polygon; } BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector) { //check donut BoostMask::BoostRingVector::const_iterator it1; BoostMask::BoostRingVector::const_iterator it2; BoostMask::BoostPolygonVector boostPolygonVector; std::vector donutIndexVector;//store the outer and inner ring index BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon //Get donut index and donut polygon unsigned int index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, ++index1) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } //if not jet, check now if (!it1IsDonut) { bool it2IsDonut = false; unsigned int index2 = 0; for (it2 = aRingVector.begin(); it2 != aRingVector.end(); ++it2, ++index2) { if (it2 != it1) { BoostMask::BoostPolygon2D polygon2D; if (::boost::geometry::within(*it1, *it2)) { ::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring it2IsDonut = true; } //if donut else if (::boost::geometry::within(*it2, *it1)) { ::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring it2IsDonut = true; } if (it2IsDonut) { donutIndexVector.push_back(index1); donutIndexVector.push_back(index2); donutVector.push_back(polygon2D);//store donut polygon break;//Only store the first donut! } } } } } //Store no donut polygon to boostPolygonVector index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, ++index1) { bool it1IsDonut = false; //check if the ring is the outer or inner of a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } if (!it1IsDonut) { BoostMask::BoostPolygon2D polygon2D; ::boost::geometry::append(polygon2D, *it1); boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut } } //Append donut polygon to boostPolygonVector BoostMask::BoostPolygonVector::iterator itDonut; for (itDonut = donutVector.begin(); itDonut != donutVector.end(); ++itDonut) { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMask.h b/code/masks/legacy/rttbBoostMask_LEGACY.h similarity index 93% copy from code/masks/boost/rttbBoostMask.h copy to code/masks/legacy/rttbBoostMask_LEGACY.h index af2ced7..86bfaf4 100644 --- a/code/masks/boost/rttbBoostMask.h +++ b/code/masks/legacy/rttbBoostMask_LEGACY.h @@ -1,156 +1,156 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) +// @version $Revision: 1221 $ (last changed revision) +// @date $Date: 2015-12-01 13:43:31 +0100 (Di, 01 Dez 2015) $ (last change date) +// @author $Author: hentsch $ (last changed by) */ -#ifndef __BOOST_MASK_H -#define __BOOST_MASK_H +#ifndef __BOOST_MASK_L_H +#define __BOOST_MASK_L_H #include "rttbBaseType.h" #include "rttbStructure.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include #include #include #include #include namespace rttb { namespace masks { - namespace boost + namespace boostLegacy { /*! @class BoostMask * @brief Implementation of voxelization using boost::geometry. * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. * (A structure without self interseciton means all contours of the structure have no self intersection, and * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored * and the mask will be calculated, however, it may cause errors in the mask results. */ class BoostMask { public: typedef ::boost::shared_ptr GeometricInfoPointer; typedef core::Structure::StructTypePointer StructPointer; typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; /*! @brief Constructor * @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is NULL * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. * @exception InvalidParameterException thrown if strict is true and the structure has self intersections */ BoostMask(GeometricInfoPointer aDoseGeoInfo, StructPointer aStructure, bool strict = true); /*! @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; GeometricInfoPointer _geometricInfo; StructPointer _structure; bool _strict; //vector of the MaskVoxel inside the structure MaskVoxelListPointer _voxelInStructure; /*! @brief If the mask is up to date */ bool _isUpToDate; /*! @brief Voxelization and generate mask */ void calcMask(); /*! @brief Check if the structure has self intersections*/ bool hasSelfIntersections(); /*! @brief Get the min/max voxel index of the bounding box of the polygon in the slice with sliceNumber * @param sliceNumber slice number between 0 and _geometricInfo->getNumSlices() * @param intersectionSlicePolygons Get the polygons intersecting the slice * @exception InvalidParameterException thrown if sliceNumber < 0 or sliceNumber >= _geometricInfo->getNumSlices() * @return Return the 4 voxel index of the bounding box */ VoxelIndexVector getBoundingBox(unsigned int sliceNumber, const BoostPolygonVector& intersectionSlicePolygons); /*! @brief Get intersection polygons of the contour and a voxel polygon * @param aVoxelIndex3D The 3d grid index of the voxel * @param intersectionSlicePolygons The polygons of the slice intersecting the voxel * @return Return all intersetion polygons of the structure and the voxel */ BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons); /*! @brief Calculate the area of all polygons * @param aPolygonDeque The deque of polygons * @return Return the area of all polygons */ double calcArea(const BoostPolygonDeque& aPolygonDeque); /*! @brief Get grid index of a mask voxel * @param aMaskVoxel A mask voxel * @return Return the 3d grid index of the mask voxel */ VoxelGridIndex3D getGridIndex3D(const core::MaskVoxel& aMaskVoxel); /*! @brief Convert RTTB polygon to boost polygon*/ BoostRing2D convert(const rttb::PolygonType& aRTTBPolygon); /*! @brief Get the intersection slice of the voxel, return the polygons (with (donut) or without holes) in the slice * @param aVoxelGridIndexZ The z grid index (slice number) of the voxel * @return Return the structure polygons intersecting the slice */ BoostPolygonVector getIntersectionSlicePolygons(const rttb::GridIndexType aVoxelGridIndexZ); /*! @brief Get the voxel 2d contour polygon*/ BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D); /*! @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); }; } } } #endif \ No newline at end of file diff --git a/testing/examples/VoxelizationValidationTest.cpp b/testing/examples/VoxelizationValidationTest.cpp index 0648cff..3980650 100644 --- a/testing/examples/VoxelizationValidationTest.cpp +++ b/testing/examples/VoxelizationValidationTest.cpp @@ -1,291 +1,293 @@ // ----------------------------------------------------------------------- // 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 "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 "rttbBoostMaskRedesign.h" -#include "rttbBoostMaskRedesignAccessor.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 two differnt voxelizations: OTB and Boost. + 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: structure file name - // 2: dose1 file name + //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 BoostMask_DIRNAME; + std::string BoostMaskLegacy_DIRNAME; std::string OTBMask_DIRNAME; - std::string BoostMaskRedesign_DIRNAME; + std::string BoostMask_DIRNAME; - std::string RTDose_BoostRedesign; - std::string RTStr_BoostRedesign; + std::string RTDOSE2_FILENAME; + std::string RTSTRUCT2_FILENAME; if (argc > 4) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; - BoostMask_DIRNAME = argv[3]; + BoostMaskLegacy_DIRNAME = argv[3]; OTBMask_DIRNAME = argv[4]; - BoostMaskRedesign_DIRNAME = argv[5]; + BoostMask_DIRNAME = argv[5]; - RTDose_BoostRedesign = argv[6]; - RTStr_BoostRedesign = argv[7]; + RTDOSE2_FILENAME = argv[6]; + RTSTRUCT2_FILENAME = argv[7]; } 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(); 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 + //create Boost MaskAccessor_LEGACY MaskAccessorPointer boostMaskAccessorPtr - = ::boost::make_shared + = ::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 << BoostMask_DIRNAME << j << ".mhd"; + 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 + = ::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 << BoostMaskRedesign_DIRNAME << j << ".mhd"; + 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(RTDose_BoostRedesign.c_str()); + io::itk::ITKImageFileAccessorGenerator doseAccessorGenerator2(RTDOSE2_FILENAME.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); StructureSetPointer rtStructureSet2 = io::dicom::DicomFileStructureSetGenerator( - RTStr_BoostRedesign.c_str()).generateStructureSet(); + 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 + //create Boost MaskAccessor Legacy MaskAccessorPointer boostMaskAccessorPtr - = ::boost::make_shared + = ::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 + = ::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/masks/boost/BoostMaskRedesignTest.cpp b/testing/masks/boost/BoostMaskTest.cpp similarity index 89% rename from testing/masks/boost/BoostMaskRedesignTest.cpp rename to testing/masks/boost/BoostMaskTest.cpp index 0784097..5231a26 100644 --- a/testing/masks/boost/BoostMaskRedesignTest.cpp +++ b/testing/masks/boost/BoostMaskTest.cpp @@ -1,252 +1,243 @@ // ----------------------------------------------------------------------- // 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 #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "../../core/DummyStructure.h" #include "../../core/DummyDoseAccessor.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbGenericDoseIterator.h" #include "rttbDVHCalculator.h" #include "rttbGenericMaskedDoseIterator.h" -#include "rttbBoostMaskRedesign.h" #include "rttbBoostMask.h" #include "rttbBoostMaskAccessor.h" -#include "rttbBoostMaskRedesignAccessor.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace testing { /*! @brief BoostMaskRedesignTest. 1) test constructors 2) test getRelevantVoxelVector 3) test getMaskAt */ - int BoostMaskRedesignTest(int argc, char* argv[]) + int BoostMaskTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::Structure::StructTypePointer StructTypePointer; PREPARE_DEFAULT_TEST_REPORTING; // generate test dose. geometric info: patient position (-25, -2, 35), center of the 1st.voxel boost::shared_ptr spTestDoseAccessor = boost::make_shared(); boost::shared_ptr geometricPtr = boost::make_shared (spTestDoseAccessor->getGeometricInfo()); DummyStructure myStructGenerator(spTestDoseAccessor->getGeometricInfo()); //generate test structure. contours are (-20,0.5,38.75); (-12.5,0.5,38.75); (-12.5,10.5,38.75); (-20,10.5,38.75); //(-20, 0.5, 41.25); (-12.5, 0.5, 41.25); (-12.5, 10.5, 41.25); (-20, 10.5, 41.25); core::Structure myTestStruct = myStructGenerator.CreateRectangularStructureCentered(2,3); StructTypePointer spMyStruct = boost::make_shared(myTestStruct); //generate test structure 2. contours are (-20,0.5,38.75); (-12.5,0.5,38.75); (-12.5,10.5,38.75); (-20,10.5,38.75); //(-20, 0.5, 40); (-12.5, 0.5, 40); (-12.5, 10.5, 40); (-20, 10.5, 40); core::Structure myTestStruct2 = myStructGenerator.CreateRectangularStructureCenteredContourPlaneThicknessNotEqualDosePlaneThickness(2); StructTypePointer spMyStruct2 = boost::make_shared(myTestStruct2); //1) test BoostMask & BoostMaskAccessor constructor - CHECK_NO_THROW(rttb::masks::boostRedesign::BoostMask(geometricPtr, spMyStruct)); - rttb::masks::boostRedesign::BoostMask boostMask = rttb::masks::boostRedesign::BoostMask( + CHECK_NO_THROW(rttb::masks::boost::BoostMask(geometricPtr, spMyStruct)); + rttb::masks::boost::BoostMask boostMask = rttb::masks::boost::BoostMask( geometricPtr, spMyStruct); - CHECK_NO_THROW(rttb::masks::boostRedesign::BoostMaskAccessor(spMyStruct, + CHECK_NO_THROW(rttb::masks::boost::BoostMaskAccessor(spMyStruct, spTestDoseAccessor->getGeometricInfo(), true)); - rttb::masks::boostRedesign::BoostMaskAccessor boostMaskAccessor(spMyStruct, + rttb::masks::boost::BoostMaskAccessor boostMaskAccessor(spMyStruct, spTestDoseAccessor->getGeometricInfo(), true); //2) test getRelevantVoxelVector CHECK_NO_THROW(boostMask.getRelevantVoxelVector()); CHECK_NO_THROW(boostMaskAccessor.getRelevantVoxelVector()); //3) test getMaskAt const VoxelGridIndex3D inMask1(2, 1, 2); //corner between two contours slice -> volumeFraction = 0.25 const VoxelGridIndex3D inMask2(3, 4, 2); //inside between two contours slice ->volumeFraction = 1 const VoxelGridIndex3D inMask3(4, 5, 2); //side between two contours slice -> volumeFraction = 0.5 const VoxelGridIndex3D inMask4(2, 1, 1); //corner on the first contour slice -> volumeFraction = 0.25/2 = 0.125 const VoxelGridIndex3D inMask5(2, 1, 3); //corner on the last contour slice -> volumeFraction = 0.25/2 = 0.125 const VoxelGridIndex3D inMask6(3, 4, 1); //inside on the first contour slice ->volumeFraction = 1 /2 = 0.5 const VoxelGridIndex3D outMask1(7, 5, 4); const VoxelGridIndex3D outMask2(2, 1, 0); const VoxelGridIndex3D outMask3(2, 1, 4); VoxelGridID testId; double errorConstant = 1e-7; core::MaskVoxel tmpMV1(0), tmpMV2(0); CHECK(boostMaskAccessor.getMaskAt(inMask1, tmpMV1)); geometricPtr->convert(inMask1, testId); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_CLOSE(0.25, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor.getMaskAt(inMask2, tmpMV1)); CHECK(geometricPtr->convert(inMask2, testId)); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(1, tmpMV1.getRelevantVolumeFraction()); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor.getMaskAt(inMask3, tmpMV1)); CHECK(geometricPtr->convert(inMask3, testId)); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_CLOSE(0.5, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor.getMaskAt(inMask4, tmpMV1)); CHECK(geometricPtr->convert(inMask4, testId)); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_CLOSE(0.125, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor.getMaskAt(inMask5, tmpMV1)); CHECK(geometricPtr->convert(inMask5, testId)); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_CLOSE(0.125, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor.getMaskAt(inMask6, tmpMV1)); CHECK(geometricPtr->convert(inMask6, testId)); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_CLOSE(0.5, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(!boostMaskAccessor.getMaskAt(outMask1, tmpMV1)); CHECK(geometricPtr->convert(outMask1, testId)); CHECK(!boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(0, tmpMV1.getRelevantVolumeFraction()); CHECK(!boostMaskAccessor.getMaskAt(outMask2, tmpMV1)); CHECK(geometricPtr->convert(outMask2, testId)); CHECK(!boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(0, tmpMV1.getRelevantVolumeFraction()); CHECK(!boostMaskAccessor.getMaskAt(outMask3, tmpMV1)); CHECK(geometricPtr->convert(outMask3, testId)); CHECK(!boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(0, tmpMV1.getRelevantVolumeFraction()); - //4) Exception tests if using old boost mask - MaskAccessorPointer boostMaskAccessorPtr - = ::boost::make_shared - (spMyStruct2, spTestDoseAccessor->getGeometricInfo()); - CHECK_THROW_EXPLICIT(boostMaskAccessorPtr->updateMask(), rttb::core::InvalidParameterException); - - rttb::masks::boostRedesign::BoostMask boostMask2 = rttb::masks::boostRedesign::BoostMask( + rttb::masks::boost::BoostMask boostMask2 = rttb::masks::boost::BoostMask( geometricPtr, spMyStruct2); CHECK_NO_THROW(boostMask2.getRelevantVoxelVector()); - rttb::masks::boostRedesign::BoostMaskAccessor boostMaskAccessor2(spMyStruct2, + rttb::masks::boost::BoostMaskAccessor boostMaskAccessor2(spMyStruct2, spTestDoseAccessor->getGeometricInfo(), true); CHECK_NO_THROW(boostMaskAccessor2.getRelevantVoxelVector()); - //5) test getMaskAt CHECK(boostMaskAccessor2.getMaskAt(inMask1, tmpMV1)); geometricPtr->convert(inMask1, testId); CHECK(boostMaskAccessor2.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); //corner, the first contour weight 0.25, the second contour weights 0.5 -> volumeFraction = 0.25*0.25 + 1.25*0.5 = 0.1875 CHECK_CLOSE(0.1875, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor2.getMaskAt(inMask2, tmpMV1)); CHECK(geometricPtr->convert(inMask2, testId)); CHECK(boostMaskAccessor2.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); //inside, the first contour weight 0.25, the second contour weights 0.5 -> ->volumeFraction = 1*0.25 + 1*0.5 = 0.75 CHECK_EQUAL(0.75, tmpMV1.getRelevantVolumeFraction()); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor2.getMaskAt(inMask3, tmpMV1)); CHECK(geometricPtr->convert(inMask3, testId)); CHECK(boostMaskAccessor2.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); //side the first contour weight 0.25, the second contour weights 0.5 -> ->volumeFraction = 0.5*0.25 + 0.5*0.5 = 0.75 CHECK_CLOSE(0.375, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor2.getMaskAt(inMask4, tmpMV1)); CHECK(geometricPtr->convert(inMask4, testId)); CHECK(boostMaskAccessor2.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); //corner on the first contour slice, weight 0.25 -> volumeFraction = 0.25*0.25 = 0.0625 CHECK_CLOSE(0.0625, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor2.getMaskAt(inMask6, tmpMV1)); CHECK(geometricPtr->convert(inMask6, testId)); CHECK(boostMaskAccessor2.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); //inside on the first contour slice, weight 0.25 ->volumeFraction = 1 * 0.25 = 0.25 CHECK_CLOSE(0.25, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(!boostMaskAccessor2.getMaskAt(outMask1, tmpMV1)); CHECK(geometricPtr->convert(outMask1, testId)); CHECK(!boostMaskAccessor2.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(0, tmpMV1.getRelevantVolumeFraction()); //CHECK_EQUAL(testId,tmpMV1.getVoxelGridID()); -> return value will not be valid outside the mask CHECK(!boostMaskAccessor2.getMaskAt(outMask2, tmpMV1)); CHECK(geometricPtr->convert(outMask2, testId)); CHECK(!boostMaskAccessor2.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(0, tmpMV1.getRelevantVolumeFraction()); //CHECK_EQUAL(testId,tmpMV1.getVoxelGridID()); -> return value will not be valid outside the mask CHECK(!boostMaskAccessor2.getMaskAt(outMask3, tmpMV1)); CHECK(geometricPtr->convert(outMask3, testId)); CHECK(!boostMaskAccessor2.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(0, tmpMV1.getRelevantVolumeFraction()); //CHECK_EQUAL(testId,tmpMV1.getVoxelGridID()); -> return value will not be valid outside the mask RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/masks/boost/CMakeLists.txt b/testing/masks/boost/CMakeLists.txt index 6dd75d4..051b1e3 100644 --- a/testing/masks/boost/CMakeLists.txt +++ b/testing/masks/boost/CMakeLists.txt @@ -1,22 +1,20 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(Boost_Mask_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbBoostMaskTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- -ADD_TEST(BoostMaskAccessorTest ${Boost_Mask_TESTS} BoostMaskAccessorTest -) -ADD_TEST(BoostMaskRedesignTest ${Boost_Mask_TESTS} BoostMaskRedesignTest +ADD_TEST(BoostMaskTest ${Boost_Mask_TESTS} BoostMaskTest ) RTTB_CREATE_TEST_MODULE(rttbBoostMask DEPENDS RTTBDicomIO RTTBMasks RTTBBoostMask PACKAGE_DEPENDS Boost BoostBinaries Litmus DCMTK) diff --git a/testing/masks/boost/files.cmake b/testing/masks/boost/files.cmake index ea3c821..004661b 100644 --- a/testing/masks/boost/files.cmake +++ b/testing/masks/boost/files.cmake @@ -1,14 +1,13 @@ SET(CPP_FILES ../../core/DummyStructure.cpp ../../core/CreateTestStructure.cpp ../../core/DummyDoseAccessor.cpp - BoostMaskAccessorTest.cpp - BoostMaskRedesignTest.cpp + BoostMaskTest.cpp rttbBoostMaskTests.cpp ) SET(H_FILES ../../core/DummyStructure.h ../../core/CreateTestStructure.h ../../core/DummyDoseAccessor.h ) diff --git a/testing/masks/boost/rttbBoostMaskTests.cpp b/testing/masks/boost/rttbBoostMaskTests.cpp index 1f552c0..4dfa247 100644 --- a/testing/masks/boost/rttbBoostMaskTests.cpp +++ b/testing/masks/boost/rttbBoostMaskTests.cpp @@ -1,64 +1,63 @@ // ----------------------------------------------------------------------- // 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 rttbAlgorithmsTests for the test driver // and all it expects is that you have a function called RegisterTests #if defined(_MSC_VER) #pragma warning ( disable : 4786 ) #endif #include "litMultiTestsMain.h" namespace rttb { namespace testing { void registerTests() { - LIT_REGISTER_TEST(BoostMaskAccessorTest); - LIT_REGISTER_TEST(BoostMaskRedesignTest); + LIT_REGISTER_TEST(BoostMaskTest); } } } int main(int argc, char* argv[]) { int result = 0; rttb::testing::registerTests(); try { result = lit::multiTestsMain(argc, argv); } catch (const std::exception& /*e*/) { result = -1; } catch (...) { result = -1; } return result; } diff --git a/testing/masks/boost/BoostMaskAccessorTest.cpp b/testing/masks/legacy/BoostMaskLegacyTest.cpp similarity index 80% rename from testing/masks/boost/BoostMaskAccessorTest.cpp rename to testing/masks/legacy/BoostMaskLegacyTest.cpp index caefd8b..f3cf1e0 100644 --- a/testing/masks/boost/BoostMaskAccessorTest.cpp +++ b/testing/masks/legacy/BoostMaskLegacyTest.cpp @@ -1,130 +1,130 @@ // ----------------------------------------------------------------------- // 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) +// @version $Revision: 1221 $ (last changed revision) +// @date $Date: 2015-12-01 13:43:31 +0100 (Di, 01 Dez 2015) $ (last change date) +// @author $Author: hentsch $ (last changed by) */ #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbMaskVoxel.h" #include "../../core/DummyStructure.h" #include "../../core/DummyDoseAccessor.h" -#include "rttbBoostMask.h" -#include "rttbBoostMaskAccessor.h" +#include "rttbBoostMask_LEGACY.h" +#include "rttbBoostMaskAccessor_LEGACY.h" namespace rttb { namespace testing { /*! @brief BoostMaskAccessorTest. 1) test constructors 2) test getRelevantVoxelVector 3) test getMaskAt */ - int BoostMaskAccessorTest(int argc, char* argv[]) + int BoostMaskLegacyTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; typedef core::Structure::StructTypePointer StructTypePointer; - typedef masks::boost::BoostMaskAccessor::MaskVoxelListPointer MaskVoxelListPointer; - typedef masks::boost::BoostMaskAccessor::MaskVoxelList MaskVoxelList; + typedef masks::boostLegacy::BoostMaskAccessor::MaskVoxelListPointer MaskVoxelListPointer; + typedef masks::boostLegacy::BoostMaskAccessor::MaskVoxelList MaskVoxelList; // generate test structure set boost::shared_ptr spTestDoseAccessor = boost::make_shared(); DummyStructure myStructGenerator(spTestDoseAccessor->getGeometricInfo()); GridIndexType zPlane = 4; core::Structure myTestStruct = myStructGenerator.CreateRectangularStructureCentered(zPlane); StructTypePointer spMyStruct = boost::make_shared(myTestStruct); boost::shared_ptr geometricPtr = boost::make_shared (spTestDoseAccessor->getGeometricInfo()); //1) test BoostMask and BoostMaskAccessor constructor - CHECK_NO_THROW(rttb::masks::boost::BoostMask(geometricPtr, spMyStruct)); - rttb::masks::boost::BoostMask boostMask = rttb::masks::boost::BoostMask(geometricPtr, spMyStruct); - CHECK_NO_THROW(rttb::masks::boost::BoostMaskAccessor(spMyStruct, + CHECK_NO_THROW(rttb::masks::boostLegacy::BoostMask(geometricPtr, spMyStruct)); + rttb::masks::boostLegacy::BoostMask boostMask = rttb::masks::boostLegacy::BoostMask(geometricPtr, spMyStruct); + CHECK_NO_THROW(rttb::masks::boostLegacy::BoostMaskAccessor(spMyStruct, spTestDoseAccessor->getGeometricInfo())); - rttb::masks::boost::BoostMaskAccessor boostMaskAccessor(spMyStruct, + rttb::masks::boostLegacy::BoostMaskAccessor boostMaskAccessor(spMyStruct, spTestDoseAccessor->getGeometricInfo()); //2) test getRelevantVoxelVector CHECK_NO_THROW(boostMask.getRelevantVoxelVector()); CHECK_NO_THROW(boostMaskAccessor.getRelevantVoxelVector()); //3) test getMaskAt const VoxelGridIndex3D inMask1(2, 1, 4); //corner -> volumeFraction = 0.25 const VoxelGridIndex3D inMask2(3, 4, 4); //inside ->volumeFraction = 1 const VoxelGridIndex3D inMask3(4, 5, 4); //side -> volumeFraction = 0.5 const VoxelGridIndex3D outMask1(7, 5, 4); const VoxelGridIndex3D outMask2(2, 1, 2); VoxelGridID testId; double errorConstant = 1e-7; core::MaskVoxel tmpMV1(0), tmpMV2(0); CHECK(boostMaskAccessor.getMaskAt(inMask1, tmpMV1)); geometricPtr->convert(inMask1, testId); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_CLOSE(0.25, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor.getMaskAt(inMask2, tmpMV1)); CHECK(geometricPtr->convert(inMask2, testId)); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(1, tmpMV1.getRelevantVolumeFraction()); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(boostMaskAccessor.getMaskAt(inMask3, tmpMV1)); CHECK(geometricPtr->convert(inMask3, testId)); CHECK(boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_CLOSE(0.5, tmpMV1.getRelevantVolumeFraction(), errorConstant); CHECK_EQUAL(testId, tmpMV1.getVoxelGridID()); CHECK(!boostMaskAccessor.getMaskAt(outMask1, tmpMV1)); CHECK(geometricPtr->convert(outMask1, testId)); CHECK(!boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(0, tmpMV1.getRelevantVolumeFraction()); //CHECK_EQUAL(testId,tmpMV1.getVoxelGridID()); -> return value will not be valid outside the mask CHECK(!boostMaskAccessor.getMaskAt(outMask2, tmpMV1)); CHECK(geometricPtr->convert(outMask2, testId)); CHECK(!boostMaskAccessor.getMaskAt(testId, tmpMV2)); CHECK_EQUAL(tmpMV1, tmpMV2); CHECK_EQUAL(0, tmpMV1.getRelevantVolumeFraction()); //CHECK_EQUAL(testId,tmpMV1.getVoxelGridID()); -> return value will not be valid outside the mask RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/masks/legacy/CMakeLists.txt b/testing/masks/legacy/CMakeLists.txt index 443c87d..1ef9f98 100644 --- a/testing/masks/legacy/CMakeLists.txt +++ b/testing/masks/legacy/CMakeLists.txt @@ -1,20 +1,23 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(OTB_MASK_TESTS ${EXECUTABLE_OUTPUT_PATH}/rttbOTBMaskTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(OTBMaskAccessorTest ${OTB_MASK_TESTS} OTBMaskAccessorTest ) +ADD_TEST(BoostMaskLegacyTest ${OTB_MASK_TESTS} BoostMaskLegacyTest +) + RTTB_CREATE_TEST_MODULE(rttbOTBMask DEPENDS RTTBMasks RTTBOTBMask PACKAGE_DEPENDS Litmus DCMTK boost) diff --git a/testing/masks/legacy/files.cmake b/testing/masks/legacy/files.cmake index 5720825..144e8bb 100644 --- a/testing/masks/legacy/files.cmake +++ b/testing/masks/legacy/files.cmake @@ -1,17 +1,18 @@ SET(CPP_FILES ../../core/DummyStructure.cpp ../../core/CreateTestStructure.cpp ../../core/DummyDoseAccessor.cpp OTBMaskAccessorTest.cpp + BoostMaskLegacyTest.cpp ../rttbMaskVoxelListTester.cpp ../rttbMaskRectStructTester.cpp rttbOTBMaskTests.cpp ) SET(H_FILES ../../core/DummyStructure.h ../../core/CreateTestStructure.h ../../core/DummyDoseAccessor.h ../rttbMaskVoxelListTester.h ../rttbMaskRectStructTester.h ) diff --git a/testing/masks/legacy/rttbOTBMaskTests.cpp b/testing/masks/legacy/rttbOTBMaskTests.cpp index c865e78..23c8143 100644 --- a/testing/masks/legacy/rttbOTBMaskTests.cpp +++ b/testing/masks/legacy/rttbOTBMaskTests.cpp @@ -1,63 +1,64 @@ // ----------------------------------------------------------------------- // 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 rttbAlgorithmsTests for the test driver // and all it expects is that you have a function called RegisterTests #if defined(_MSC_VER) #pragma warning ( disable : 4786 ) #endif #include "litMultiTestsMain.h" namespace rttb { namespace testing { void registerTests() { LIT_REGISTER_TEST(OTBMaskAccessorTest); + LIT_REGISTER_TEST(BoostMaskLegacyTest); } } } int main(int argc, char* argv[]) { int result = 0; rttb::testing::registerTests(); try { result = lit::multiTestsMain(argc, argv); } catch (const std::exception& /*e*/) { result = -1; } catch (...) { result = -1; } return result; }