diff --git a/code/masks/boost/rttbBoostMaskRedesign.cpp b/code/masks/boost/rttbBoostMaskRedesign.cpp index c75a2a4..88b60ec 100644 --- a/code/masks/boost/rttbBoostMaskRedesign.cpp +++ b/code/masks/boost/rttbBoostMaskRedesign.cpp @@ -1,580 +1,580 @@ #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); - + _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/testing/core/CreateTestStructure.cpp b/testing/core/CreateTestStructure.cpp index 79df5dc..98e82a5 100644 --- a/testing/core/CreateTestStructure.cpp +++ b/testing/core/CreateTestStructure.cpp @@ -1,670 +1,691 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include "CreateTestStructure.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace testing { CreateTestStructure::~CreateTestStructure() {} CreateTestStructure::CreateTestStructure(const core::GeometricInfo& aGeoInfo) { _geoInfo = aGeoInfo; } PolygonType CreateTestStructure::createPolygonLeftUpper(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); polygon.push_back(p1); std::cout << "(" << p1.x() << "," << p1.y() << "," << p1.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonCenter(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() / 2; p(1) = p1.y() + delta.y() / 2; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } + PolygonType CreateTestStructure::createPolygonCenterOnPlaneCenter(std::vector aVoxelVector, + GridIndexType sliceNumber) + { + + PolygonType polygon; + + for (size_t i = 0; i < aVoxelVector.size(); i++) + { + VoxelGridIndex3D voxelIndex; + DoubleVoxelGridIndex3D indexDouble = DoubleVoxelGridIndex3D((aVoxelVector.at(i)).x(), (aVoxelVector.at(i)).y(), + sliceNumber); + + WorldCoordinate3D p1; + _geoInfo.geometryCoordinateToWorldCoordinate(indexDouble, p1); + + polygon.push_back(p1); + } + + std::cout << std::endl; + return polygon; + } PolygonType CreateTestStructure::createPolygonBetweenUpperLeftAndCenter( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() / 4; p(1) = p1.y() + delta.y() / 4; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonBetweenLowerRightAndCenter( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() * 0.75; p(1) = p1.y() + delta.y() * 0.75; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonLeftEdgeMiddle(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x(); p(1) = p1.y() + delta.y() * 0.5; p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonUpperCenter(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p; p(0) = p1.x() + delta.x() * 0.5; p(1) = p1.y(); p(2) = p1.z(); polygon.push_back(p); std::cout << "(" << p.x() << "," << p.y() << "," << p.z() << ")" << "; "; } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonIntermediatePoints(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; for (size_t i = 0; i < aVoxelVector.size(); i++) { VoxelGridIndex3D voxelIndex; VoxelGridIndex3D voxelIndex2; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; if (i < (aVoxelVector.size() - 1)) { voxelIndex2(0) = (aVoxelVector.at(i + 1)).x(); voxelIndex2(1) = (aVoxelVector.at(i + 1)).y(); voxelIndex2(2) = sliceNumber; } else { voxelIndex2(0) = (aVoxelVector.at(0)).x(); voxelIndex2(1) = (aVoxelVector.at(0)).y(); voxelIndex2(2) = sliceNumber; } WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D p2; _geoInfo.indexToWorldCoordinate(voxelIndex2, p2); SpacingVectorType3D delta2 = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x() + delta.x() * 0.75; wp1(1) = p1.y() + delta.y() * 0.75; wp1(2) = p1.z(); WorldCoordinate3D wp2; wp2(0) = p2.x() + delta.x() * 0.75; wp2(1) = p2.y() + delta.y() * 0.75; wp2(2) = p2.z(); polygon.push_back(wp1); double diffX = (wp2.x() - wp1.x()) / 1000.0; double diffY = (wp2.y() - wp1.y()) / 1000.0; WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; for (int i = 0 ; i < 1000 ; i++) { wp_intermediate(0) = wp1.x() + i * diffX; wp_intermediate(1) = wp1.y() + i * diffY; polygon.push_back(wp_intermediate); } } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createPolygonCircle(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); double radius = 2 * delta.x(); double frac_radius = (radius * 0.001); double correct_y = (delta.x() / delta.y()); for (int i = 0 ; i <= 1000 ; i++) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() + frac_radius * i; wp_intermediate(1) = wp1.y() - (y_offset * correct_y) ; polygon.push_back(wp_intermediate); } for (int i = 1000 ; i >= 0 ; i--) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))) ; wp_intermediate(0) = wp1.x() + frac_radius * i; wp_intermediate(1) = wp1.y() + (y_offset * correct_y); polygon.push_back(wp_intermediate); } for (int i = 0 ; i <= 1000 ; i++) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() - frac_radius * i; wp_intermediate(1) = wp1.y() + y_offset * correct_y ; polygon.push_back(wp_intermediate); } for (int i = 1000 ; i >= 0 ; i--) { double y_offset = sqrt((radius * radius) - ((frac_radius * i) * (frac_radius * i))); wp_intermediate(0) = wp1.x() - frac_radius * i ; wp_intermediate(1) = wp1.y() - (y_offset * correct_y); polygon.push_back(wp_intermediate); } } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSeveralSectionsInsideOneVoxelA( std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x() + (delta.x() * 0.25); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.25); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.75); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.75); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.25); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.25); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.5); wp_intermediate(1) = wp1.y() + (delta.y() * 2.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.75); wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.75); wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 3.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.75); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSelfTouchingA(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y(); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } PolygonType CreateTestStructure::createStructureSelfTouchingB(std::vector aVoxelVector, GridIndexType sliceNumber) { PolygonType polygon; if (aVoxelVector.size() > 0) { int i = 0; VoxelGridIndex3D voxelIndex; voxelIndex(0) = (aVoxelVector.at(i)).x(); voxelIndex(1) = (aVoxelVector.at(i)).y(); voxelIndex(2) = sliceNumber; WorldCoordinate3D p1; _geoInfo.indexToWorldCoordinate(voxelIndex, p1); SpacingVectorType3D delta = _geoInfo.getSpacing(); WorldCoordinate3D wp1; wp1(0) = p1.x(); wp1(1) = p1.y(); wp1(2) = p1.z(); WorldCoordinate3D wp_intermediate; wp_intermediate(0) = 0; wp_intermediate(1) = 0; wp_intermediate(2) = p1.z(); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y(); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 1.0); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x(); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() + (delta.x() * 0.5); wp_intermediate(1) = wp1.y() + (delta.y() * 1.0); polygon.push_back(wp_intermediate); wp_intermediate(0) = wp1.x() ; wp_intermediate(1) = wp1.y() + (delta.y() * 0.5); polygon.push_back(wp_intermediate); } std::cout << std::endl; return polygon; } }//testing }//rttb diff --git a/testing/core/CreateTestStructure.h b/testing/core/CreateTestStructure.h index e8deecf..627c5cf 100644 --- a/testing/core/CreateTestStructure.h +++ b/testing/core/CreateTestStructure.h @@ -1,74 +1,81 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "rttbStructure.h" #include "rttbBaseType.h" #include "rttbGeometricInfo.h" namespace rttb { namespace testing { /*!@class CreateTestStructure @brief Create dummy structures for testing. */ class CreateTestStructure { private: core::GeometricInfo _geoInfo; CreateTestStructure() {}; public: ~CreateTestStructure(); CreateTestStructure(const core::GeometricInfo& aGeoInfo); PolygonType createPolygonLeftUpper(std::vector aVoxelVector, GridIndexType sliceNumber); + + /*The z coordinate of the polygon is the upper side of the voxel on the slice*/ PolygonType createPolygonCenter(std::vector aVoxelVector, GridIndexType sliceNumber); + + /*The z coordinate of the polygon is the center of the voxel on the slice*/ + PolygonType createPolygonCenterOnPlaneCenter(std::vector aVoxelVector, + GridIndexType sliceNumber); + PolygonType createPolygonBetweenUpperLeftAndCenter(std::vector aVoxelVector, GridIndexType sliceNumber); PolygonType createPolygonBetweenLowerRightAndCenter(std::vector aVoxelVector, GridIndexType sliceNumber); PolygonType createPolygonIntermediatePoints(std::vector aVoxelVector, GridIndexType sliceNumber); PolygonType createPolygonUpperCenter(std::vector aVoxelVector, GridIndexType sliceNumber); PolygonType createPolygonLeftEdgeMiddle(std::vector aVoxelVector, GridIndexType sliceNumber); PolygonType createPolygonCircle(std::vector aVoxelVector, GridIndexType sliceNumber); PolygonType createStructureSeveralSectionsInsideOneVoxelA(std::vector aVoxelVector, GridIndexType sliceNumber); PolygonType createStructureSelfTouchingA(std::vector aVoxelVector, GridIndexType sliceNumber); PolygonType createStructureSelfTouchingB(std::vector aVoxelVector, GridIndexType sliceNumber); }; }//testing }//rttb diff --git a/testing/core/DummyStructure.cpp b/testing/core/DummyStructure.cpp index f8dd8dc..a8980a9 100644 --- a/testing/core/DummyStructure.cpp +++ b/testing/core/DummyStructure.cpp @@ -1,655 +1,710 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include #include "DummyStructure.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace testing { DummyStructure::~DummyStructure() {} DummyStructure::DummyStructure(const core::GeometricInfo& aGeoInfo) { _geoInfo = aGeoInfo; } core::Structure DummyStructure::CreateRectangularStructureCentered(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 1); VoxelGridIndex2D another_i2(5, 1); VoxelGridIndex2D another_i3(5, 5); VoxelGridIndex2D another_i4(2, 5); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure_rectangular_centered = core::Structure(another_polySeq); return test_structure_rectangular_centered; } + core::Structure DummyStructure::CreateRectangularStructureCenteredContourPlaneThicknessNotEqualDosePlaneThickness(GridIndexType zPlane) + { + CreateTestStructure another_cts = CreateTestStructure(_geoInfo); + + std::vector another_voxelVector; + VoxelGridIndex2D another_i1(2, 1); + VoxelGridIndex2D another_i2(5, 1); + VoxelGridIndex2D another_i3(5, 5); + VoxelGridIndex2D another_i4(2, 5); + + another_voxelVector.push_back(another_i1); + another_voxelVector.push_back(another_i2); + another_voxelVector.push_back(another_i3); + another_voxelVector.push_back(another_i4); + PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector, zPlane); + + PolygonType another_polygon2 = another_cts.createPolygonCenterOnPlaneCenter(another_voxelVector, zPlane); + + PolygonSequenceType another_polySeq; + another_polySeq.push_back(another_polygon1); + + another_polySeq.push_back(another_polygon2); + + core::Structure test_structure_rectangular_centered = core::Structure(another_polySeq); + + return test_structure_rectangular_centered; + } + + core::Structure DummyStructure::CreateRectangularStructureCentered(GridIndexType fromZPlane, GridIndexType toZPlane) + { + CreateTestStructure another_cts = CreateTestStructure(_geoInfo); + + std::vector another_voxelVector; + VoxelGridIndex2D another_i1(2, 1); + VoxelGridIndex2D another_i2(5, 1); + VoxelGridIndex2D another_i3(5, 5); + VoxelGridIndex2D another_i4(2, 5); + PolygonSequenceType another_polySeq; + + for (int i = fromZPlane; i <= toZPlane; ++i){ + another_voxelVector.clear(); + another_voxelVector.push_back(another_i1); + another_voxelVector.push_back(another_i2); + another_voxelVector.push_back(another_i3); + another_voxelVector.push_back(another_i4); + PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector, i); + + another_polySeq.push_back(another_polygon1); + } + + core::Structure test_structure_rectangular_centered = core::Structure(another_polySeq); + + return test_structure_rectangular_centered; + } + core::Structure DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacement( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(5, 1); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(5, 7); VoxelGridIndex2D another_i4(2, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonBetweenUpperLeftAndCenter( another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRight() { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(5, 1); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(5, 7); VoxelGridIndex2D another_i4(2, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonSequenceType another_polySeq; core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRightCounterClock( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(5, 7); VoxelGridIndex2D another_i3(8, 4); VoxelGridIndex2D another_i4(5, 1); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonBetweenLowerRightAndCenter( another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRightCounterClockIntermediatePoints( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(5, 7); VoxelGridIndex2D another_i3(8, 4); VoxelGridIndex2D another_i4(5, 1); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonIntermediatePoints(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureSeveralSeperateSectionsInsideOneVoxel( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 2); another_voxelVector.push_back(another_i1); PolygonType another_polygon1 = another_cts.createStructureSeveralSectionsInsideOneVoxelA( another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureSelfTouchingA(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 2); another_voxelVector.push_back(another_i1); PolygonType another_polygon1 = another_cts.createStructureSelfTouchingA(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureIntersectingTwoPolygonsInDifferentSlices( GridIndexType zPlane1, GridIndexType zPlane2) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); CreateTestStructure one_more_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(8, 6); VoxelGridIndex2D another_i4(2, 6); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector , zPlane1); std::vector one_more_voxelVector; VoxelGridIndex2D one_more_i1(3, 5); VoxelGridIndex2D one_more_i2(9, 5); VoxelGridIndex2D one_more_i3(9, 7); VoxelGridIndex2D one_more_i4(3, 7); one_more_voxelVector.push_back(one_more_i1); one_more_voxelVector.push_back(one_more_i2); one_more_voxelVector.push_back(one_more_i3); one_more_voxelVector.push_back(one_more_i4); PolygonType another_polygon2 = one_more_cts.createPolygonCenter(one_more_voxelVector , zPlane2); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); another_polySeq.push_back(another_polygon2); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureIntersectingTwoPolygons(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); CreateTestStructure one_more_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(8, 6); VoxelGridIndex2D another_i4(2, 6); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector , zPlane); std::vector one_more_voxelVector; VoxelGridIndex2D one_more_i1(3, 5); VoxelGridIndex2D one_more_i2(9, 5); VoxelGridIndex2D one_more_i3(9, 7); VoxelGridIndex2D one_more_i4(3, 7); one_more_voxelVector.push_back(one_more_i1); one_more_voxelVector.push_back(one_more_i2); one_more_voxelVector.push_back(one_more_i3); one_more_voxelVector.push_back(one_more_i4); PolygonType another_polygon2 = one_more_cts.createPolygonCenter(one_more_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); another_polySeq.push_back(another_polygon2); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureIntersecting(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(2, 6); VoxelGridIndex2D another_i4(8, 6); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonCenter(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouches(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 4); VoxelGridIndex2D another_i2(2, 8); VoxelGridIndex2D another_i3(4, 8); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); PolygonType another_polygon1 = another_cts.createPolygonUpperCenter(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesRotatedPointDoubeled( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 4); VoxelGridIndex2D another_i4(2, 4); VoxelGridIndex2D another_i2(4, 4); VoxelGridIndex2D another_i3(3, 8); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); PolygonType another_polygon1 = another_cts.createPolygonUpperCenter(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesCounterClockRotatedOnePointFivePi( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 5); VoxelGridIndex2D another_i2(3, 5); VoxelGridIndex2D another_i3(7, 6); VoxelGridIndex2D another_i4(7, 6); VoxelGridIndex2D another_i5(7, 4); VoxelGridIndex2D another_i6(7, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftEdgeMiddle(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesCounterClockRotatedQuaterPi( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(7, 5); VoxelGridIndex2D another_i2(7, 5); VoxelGridIndex2D another_i3(3, 6); VoxelGridIndex2D another_i4(3, 6); VoxelGridIndex2D another_i5(3, 4); VoxelGridIndex2D another_i6(3, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftEdgeMiddle(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesRotatedQuaterPi( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 4); VoxelGridIndex2D another_i2(3, 4); VoxelGridIndex2D another_i3(3, 6); VoxelGridIndex2D another_i4(3, 6); VoxelGridIndex2D another_i5(7, 5); VoxelGridIndex2D another_i6(7, 5); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftEdgeMiddle(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureCircle(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(4, 4); another_voxelVector.push_back(another_i1); PolygonType another_polygon1 = another_cts.createPolygonCircle(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesUpperRight( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 4); VoxelGridIndex2D another_i3(4, 5); VoxelGridIndex2D another_i5(5, 3); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i5); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesLowerRight( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(2, 2); VoxelGridIndex2D another_i2(2, 2); VoxelGridIndex2D another_i3(3, 1); VoxelGridIndex2D another_i4(3, 1); VoxelGridIndex2D another_i5(4, 3); VoxelGridIndex2D another_i6(4, 3); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesLowerLeft( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(3, 3); VoxelGridIndex2D another_i2(3, 3); VoxelGridIndex2D another_i3(4, 4); VoxelGridIndex2D another_i4(4, 4); VoxelGridIndex2D another_i5(2, 5); VoxelGridIndex2D another_i6(2, 5); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } core::Structure DummyStructure::CreateTestStructureInsideInsideTouchesUpperLeft( GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(5, 5); VoxelGridIndex2D another_i2(5, 5); VoxelGridIndex2D another_i3(4, 6); VoxelGridIndex2D another_i4(4, 6); VoxelGridIndex2D another_i5(3, 4); VoxelGridIndex2D another_i6(3, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); another_voxelVector.push_back(another_i5); another_voxelVector.push_back(another_i6); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } void DummyStructure::ShowTestStructure(core::Structure aStructure) { WorldCoordinate3D aPoint(0); PolygonSequenceType strVector = aStructure.getStructureVector(); for (size_t struct_index = 0 ; struct_index < strVector.size() ; struct_index++) { for (size_t point_index = 0 ; point_index < strVector.at(struct_index).size() ; point_index++) { aPoint = strVector.at(struct_index).at(point_index); std::cout << " aPoint.x " << aPoint.x() << std::endl; std::cout << " aPoint.y " << aPoint.y() << std::endl; std::cout << " aPoint.z " << aPoint.z() << std::endl; } } } core::Structure DummyStructure::CreateRectangularStructureUpperLeftRotated(GridIndexType zPlane) { CreateTestStructure another_cts = CreateTestStructure(_geoInfo); std::vector another_voxelVector; VoxelGridIndex2D another_i1(5, 1); VoxelGridIndex2D another_i2(8, 4); VoxelGridIndex2D another_i3(5, 7); VoxelGridIndex2D another_i4(2, 4); another_voxelVector.push_back(another_i1); another_voxelVector.push_back(another_i2); another_voxelVector.push_back(another_i3); another_voxelVector.push_back(another_i4); PolygonType another_polygon1 = another_cts.createPolygonLeftUpper(another_voxelVector , zPlane); PolygonSequenceType another_polySeq; another_polySeq.push_back(another_polygon1); core::Structure test_structure = core::Structure(another_polySeq); return test_structure; } }//testing }//rttb diff --git a/testing/core/DummyStructure.h b/testing/core/DummyStructure.h index c59ea11..53ee40b 100644 --- a/testing/core/DummyStructure.h +++ b/testing/core/DummyStructure.h @@ -1,100 +1,108 @@ // ----------------------------------------------------------------------- // 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 "rttbStructure.h" #include "CreateTestStructure.h" #include "rttbBaseType.h" #include "rttbGeometricInfo.h" namespace rttb { namespace testing { /*! @class DummyStructure @brief generate simple geometric testing structures. The maximal x coordinate used is 9 and the maximal y coordinate is 8. Make sure the geometricInfo corresponds to a sufficiently large data grid. @see CreateTestStructures */ class DummyStructure { private: core::GeometricInfo _geoInfo; public: ~DummyStructure(); DummyStructure(const core::GeometricInfo& aGeoInfo); const core::GeometricInfo& getGeometricInfo() { return _geoInfo; }; core::Structure CreateRectangularStructureCentered(GridIndexType zPlane); + /* Generate rectangular structure for the z slice and another slice between z and z+1. So the structure has a smaller z spacing than the dose + */ + core::Structure CreateRectangularStructureCenteredContourPlaneThicknessNotEqualDosePlaneThickness(GridIndexType zPlane); + + /* Generate rectangular structure for the z slices fromm fromZPlane(included) to toZPlane(included) + */ + core::Structure CreateRectangularStructureCentered(GridIndexType fromZPlane, GridIndexType toZPlane); + core::Structure CreateTestStructureCircle(GridIndexType zPlane); core::Structure CreateRectangularStructureUpperLeftRotated(GridIndexType zPlane); core::Structure CreateTestStructureSeveralSeperateSectionsInsideOneVoxel(GridIndexType zPlane); core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacement( GridIndexType zPlane); core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRight(); core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRightCounterClock( GridIndexType zPlane); core::Structure CreateRectangularStructureCenteredRotatedIntermediatePlacementLowerRightCounterClockIntermediatePoints( GridIndexType zPlane); core::Structure CreateTestStructureSelfTouchingA(GridIndexType zPlane); core::Structure CreateTestStructureIntersecting(GridIndexType zPlane); core::Structure CreateTestStructureIntersectingTwoPolygons(GridIndexType zPlane); core::Structure CreateTestStructureIntersectingTwoPolygonsInDifferentSlices(GridIndexType zPlane1, GridIndexType zPlane2); core::Structure CreateTestStructureInsideInsideTouches(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesRotatedQuaterPi(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesCounterClockRotatedQuaterPi( GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesCounterClockRotatedOnePointFivePi( GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesRotatedPointDoubeled(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesUpperLeft(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesLowerLeft(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesLowerRight(GridIndexType zPlane); core::Structure CreateTestStructureInsideInsideTouchesUpperRight(GridIndexType zPlane); void ShowTestStructure(core::Structure aStructure); }; }//testing }//rttb diff --git a/testing/examples/CMakeLists.txt b/testing/examples/CMakeLists.txt index 1111855..f039934 100644 --- a/testing/examples/CMakeLists.txt +++ b/testing/examples/CMakeLists.txt @@ -1,48 +1,49 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(CORE_TEST_EXAMPLES ${EXECUTABLE_OUTPUT_PATH}/rttbExamplesTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(RTBioModelExampleTest ${CORE_TEST_EXAMPLES} RTBioModelExampleTest "${TEST_DATA_ROOT}/DVH/Text/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_HT3.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_TV.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_virtuos_diff_trunk6.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_virtuos_diff_trunk8.txt") ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest "${TEST_DATA_ROOT}/StructureSet/DICOM/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/Dose/DICOM/ConstantTwoGridScaling.dcm" "${TEST_DATA_ROOT}/Dose/DICOM/ConstantTwoGridScaling05.dcm" "${TEST_DATA_ROOT}/Dose/DICOM/ConstantFiftyGridScaling01.dcm") ADD_TEST(RTDoseIndexTest ${CORE_TEST_EXAMPLES} RTDoseIndexTest "${TEST_DATA_ROOT}/DVH/Text/dvh_test_TV.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_HT3.txt") ADD_TEST(RTDoseStatisticsDicomTest ${CORE_TEST_EXAMPLES} RTDoseStatisticsDicomTest "${TEST_DATA_ROOT}/Dose/DICOM/ConstantTwo_withDoseGridScaling.dcm") IF(RTTB_VIRTUOS_SUPPORT AND BUILD_IO_Virtuos) ADD_TEST(RTDoseStatisticsVirtuosTest ${CORE_TEST_EXAMPLES} RTDoseStatisticsVirtuosTest "${TEST_DATA_ROOT}/Virtuos/prostate_ac/prostate_ac101.dos.gz" "${TEST_DATA_ROOT}/Virtuos/prostate_ac/prostate_ac000.vdx" "${TEST_DATA_ROOT}/Virtuos/prostate_ac/prostate_ac101.pln") ENDIF() ADD_TEST(RTDVHTest ${CORE_TEST_EXAMPLES} RTDVHTest "${TEST_DATA_ROOT}/DVH/Text/dvh_test.txt") ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest "${TEST_DATA_ROOT}/StructureSet/DICOM/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/Dose/DICOM/ConstantTwo_withDoseGridScaling.dcm" "${TEST_DATA_ROOT}/Dose/DICOM/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/Dose/DICOM/LinearIncreaseX.dcm") ADD_TEST(RTBioModelScatterPlotExampleTest ${CORE_TEST_EXAMPLES} RTBioModelScatterPlotExampleTest "${TEST_DATA_ROOT}/DVH/Text/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/DVH/Text/dvh_test_TV.txt") ADD_TEST(VoxelizationValidationTest ${CORE_TEST_EXAMPLES} VoxelizationValidationTest "${TEST_DATA_ROOT}/StructureSet/DICOM/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" -"${TEST_DATA_ROOT}/Dose/DICOM/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/BoostMask/" "${TEST_DATA_ROOT}/OTBMask/" "${TEST_DATA_ROOT}/BoostMaskRedesign/") +"${TEST_DATA_ROOT}/Dose/DICOM/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/Mask/BoostMask/" "${TEST_DATA_ROOT}/Mask/OTBMask/" "${TEST_DATA_ROOT}/Mask/BoostMaskRedesign/" +"${TEST_DATA_ROOT}/Dose/DICOM/PatBM116__RTDOSE_Main_corr.mhd" "${TEST_DATA_ROOT}/StructureSet/DICOM/__1__mm__1.2.276.0.28.19.640188827224102750711261838091512821903075578763.4195312189") IF(RTTB_VIRTUOS_SUPPORT AND BUILD_IO_Virtuos) RTTB_CREATE_TEST_MODULE(rttbExamples DEPENDS RTTBCore RTTBAlgorithms RTTBMasks RTTBOTBMask RTTBBoostMask RTTBIndices RTTBDicomIO RTTBVirtuosIO RTTBITKIO RTTBOtherIO RTTBModels PACKAGE_DEPENDS Litmus) ELSE() RTTB_CREATE_TEST_MODULE(rttbExamples DEPENDS RTTBCore RTTBAlgorithms RTTBMasks RTTBOTBMask RTTBBoostMask RTTBIndices RTTBDicomIO RTTBITKIO RTTBOtherIO RTTBModels PACKAGE_DEPENDS Litmus) ENDIF() diff --git a/testing/examples/VoxelizationValidationTest.cpp b/testing/examples/VoxelizationValidationTest.cpp index 4b1347a..0648cff 100644 --- a/testing/examples/VoxelizationValidationTest.cpp +++ b/testing/examples/VoxelizationValidationTest.cpp @@ -1,247 +1,291 @@ // ----------------------------------------------------------------------- // 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" namespace rttb { namespace testing { /*! @brief VoxelizationValidationTest. Compare two differnt voxelizations: OTB 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 // 3: directory name to write boost mask of all structures // 4: directory name to write OTB mask of all structures std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string BoostMask_DIRNAME; std::string OTBMask_DIRNAME; std::string BoostMaskRedesign_DIRNAME; + std::string RTDose_BoostRedesign; + std::string RTStr_BoostRedesign; + if (argc > 4) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; BoostMask_DIRNAME = argv[3]; OTBMask_DIRNAME = argv[4]; BoostMaskRedesign_DIRNAME = argv[5]; + + RTDose_BoostRedesign = argv[6]; + RTStr_BoostRedesign = argv[7]; } OFCondition status; DcmFileFormat fileformat; /* 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++) { - if (j != 2 && j != 3) - { 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); + */ + 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());*/ + CHECK(writer.writeFile()); clock_t start2(clock()); //create Boost MaskAccessor MaskAccessorPointer boostMaskAccessorPtr = ::boost::make_shared (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); CHECK_NO_THROW(boostMaskAccessorPtr->updateMask()); ::boost::shared_ptr spMaskedDoseIteratorTmp2 = ::boost::make_shared(boostMaskAccessorPtr, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator2(spMaskedDoseIteratorTmp2); rttb::core::DVHCalculator calc2(spMaskedDoseIterator2, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getUID()); rttb::core::DVH dvh2 = *(calc2.generateDVH()); clock_t finish2(clock()); std::cout << "Boost Mask Calculation and write file time: " << finish2 - start2 << " ms" << std::endl; //Write the mask image to a file. /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! - */ + */ - /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter2(boostMaskAccessorPtr); + rttb::io::itk::ITKImageMaskAccessorConverter itkConverter2(boostMaskAccessorPtr); CHECK(itkConverter2.process()); std::stringstream fileNameSstr2; fileNameSstr2 << BoostMask_DIRNAME << j << ".mhd"; rttb::io::itk::ImageWriter writer2(fileNameSstr2.str(), itkConverter2.getITKImage()); - CHECK(writer2.writeFile());*/ - - //create Boost MaskAccessor redesign - clock_t startR(clock()); - - MaskAccessorPointer boostMaskRPtr - = ::boost::make_shared - (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); - CHECK_NO_THROW(boostMaskRPtr->updateMask()); - ::boost::shared_ptr spMaskedDoseIteratorTmpR = - ::boost::make_shared(boostMaskRPtr, doseAccessor1); - DoseIteratorPointer spMaskedDoseIteratorR(spMaskedDoseIteratorTmpR); - rttb::core::DVHCalculator calcR(spMaskedDoseIteratorR, (rtStructureSet->getStructure(j))->getUID(), - doseAccessor1->getUID()); - rttb::core::DVH dvhR = *(calcR.generateDVH()); - clock_t finishR(clock()); - std::cout << "Boost Mask Redesign Calculation and write file time: " << finishR - startR << " ms" << - std::endl; - - //Write the mask image to a file. - /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. - To write all mask files, please use the outcommented code and call the .exe directly! - */ - /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverterR(boostMaskRPtr); - CHECK(itkConverterR.process()); + CHECK(writer2.writeFile()); + + + //create Boost MaskAccessor redesign + clock_t startR(clock()); + MaskAccessorPointer boostMaskRPtr + = ::boost::make_shared + (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); + CHECK_NO_THROW(boostMaskRPtr->updateMask()); + ::boost::shared_ptr spMaskedDoseIteratorTmpR = + ::boost::make_shared(boostMaskRPtr, doseAccessor1); + DoseIteratorPointer spMaskedDoseIteratorR(spMaskedDoseIteratorTmpR); + rttb::core::DVHCalculator calcR(spMaskedDoseIteratorR, (rtStructureSet->getStructure(j))->getUID(), + doseAccessor1->getUID()); + rttb::core::DVH dvhR = *(calcR.generateDVH()); + clock_t finishR(clock()); + std::cout << "Boost Mask Redesign Calculation and write file time: " << finishR - startR << " ms" << + std::endl; + + //Write the mask image to a file. + /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. + To write all mask files, please use the outcommented code and call the .exe directly! + */ + rttb::io::itk::ITKImageMaskAccessorConverter itkConverterR(boostMaskRPtr); + CHECK(itkConverterR.process()); + + + std::stringstream fileNameSstrR; + fileNameSstrR << BoostMaskRedesign_DIRNAME << j << ".mhd"; + rttb::io::itk::ImageWriter writerR(fileNameSstrR.str(), itkConverterR.getITKImage()); + CHECK(writerR.writeFile()); - std::stringstream fileNameSstrR; - fileNameSstrR << BoostMaskRedesign_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) + 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); + //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); - } + 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); + 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% + //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()); + DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); + StructureSetPointer rtStructureSet2 = io::dicom::DicomFileStructureSetGenerator( + RTStr_BoostRedesign.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 + MaskAccessorPointer boostMaskAccessorPtr + = ::boost::make_shared + (rtStructureSet2->getStructure(j), doseAccessor2->getGeometricInfo()); + + //Two polygons in the same slice exception using boost mask, because of the different z spacing of dose and structure + if (j != 9) + { + CHECK_THROW_EXPLICIT(boostMaskAccessorPtr->updateMask(), rttb::core::InvalidParameterException); + } + + + //create Boost MaskAccessor redesign + MaskAccessorPointer boostMaskRPtr + = ::boost::make_shared + (rtStructureSet2->getStructure(j), doseAccessor2->getGeometricInfo()); + + //No exception using redesigned boost mask + CHECK_NO_THROW(boostMaskRPtr->updateMask()); + + } + } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/masks/boost/BoostMaskRedesignTest.cpp b/testing/masks/boost/BoostMaskRedesignTest.cpp index 6e94c0c..0784097 100644 --- a/testing/masks/boost/BoostMaskRedesignTest.cpp +++ b/testing/masks/boost/BoostMaskRedesignTest.cpp @@ -1,98 +1,252 @@ // ----------------------------------------------------------------------- // 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[]) { 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; - typedef core::Structure::StructTypePointer StructTypePointer; - - // generate test structure set + // 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()); - GridIndexType zPlane = 4; - core::Structure myTestStruct = myStructGenerator.CreateRectangularStructureCentered(zPlane); + //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); - boost::shared_ptr geometricPtr = boost::make_shared - (spTestDoseAccessor->getGeometricInfo()); - - //1) test BoostMask constructor + + //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( geometricPtr, spMyStruct); - //2) test getRelevantVoxelVector - CHECK_NO_THROW(boostMask.getRelevantVoxelVector()); - rttb::masks::boostRedesign::BoostMask::MaskVoxelListPointer list = - boostMask.getRelevantVoxelVector(); - - for (int i = 0; i < list->size(); ++i) - { - std::cout << "id: " << list->at(i).getVoxelGridID() << ", " << list->at( - i).getRelevantVolumeFraction() << std::endl; - } - - - RETURN_AND_REPORT_TEST_SUCCESS; + CHECK_NO_THROW(rttb::masks::boostRedesign::BoostMaskAccessor(spMyStruct, + spTestDoseAccessor->getGeometricInfo(), true)); + rttb::masks::boostRedesign::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( + geometricPtr, spMyStruct2); + CHECK_NO_THROW(boostMask2.getRelevantVoxelVector()); + rttb::masks::boostRedesign::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