diff --git a/code/masks/rttbBoostMask.cpp b/code/masks/rttbBoostMask.cpp
index af20a7d..cbdf59d 100644
--- a/code/masks/rttbBoostMask.cpp
+++ b/code/masks/rttbBoostMask.cpp
@@ -1,565 +1,558 @@
// -----------------------------------------------------------------------
// RTToolbox - DKFZ radiotherapy quantitative evaluation library
//
// Copyright (c) German Cancer Research Center (DKFZ),
// Software development for Integrated Diagnostics and Therapy (SIDT).
// ALL RIGHTS RESERVED.
// See rttbCopyright.txt or
// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notices for more information.
//
//------------------------------------------------------------------------
-/*!
-// @file
-// @version $Revision$ (last changed revision)
-// @date $Date$ (last change date)
-// @author $Author$ (last changed by)
-*/
#include <limits>
#include <boost/geometry/geometries/register/point.hpp>
#include <boost/geometry/geometries/register/ring.hpp>
#include <boost/make_shared.hpp>
#include <boost/thread/thread.hpp>
#include "rttbBoostMask.h"
#include "rttbNullPointerException.h"
#include "rttbInvalidParameterException.h"
#include "rttbBoostMaskGenerateMaskVoxelListThread.h"
#include "rttbBoostMaskVoxelizationThread.h"
-
namespace rttb
{
namespace masks
{
namespace boost
{
BoostMask::BoostMask(core::GeometricInfo::Pointer aDoseGeoInfo,
core::Structure::Pointer aStructure, bool strict, unsigned int numberOfThreads)
: _geometricInfo(aDoseGeoInfo), _structure(aStructure),
_strict(strict), _numberOfThreads(numberOfThreads), _voxelizationThickness(0.0),
_voxelInStructure(::boost::make_shared<MaskVoxelList>())
{
_isUpToDate = false;
if (_geometricInfo == nullptr)
{
throw rttb::core::NullPointerException("Error: Geometric info is nullptr!");
}
else if (_structure == nullptr)
{
throw rttb::core::NullPointerException("Error: Structure is nullptr!");
}
if (_numberOfThreads == 0)
{
_numberOfThreads = ::boost::thread::hardware_concurrency();
if (_numberOfThreads == 0)
{
throw rttb::core::InvalidParameterException("Error: detection of the number of hardware threads is not possible. Please specify number of threads for voxelization explicitly as parameter in BoostMask.");
}
}
}
BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector()
{
if (!_isUpToDate)
{
calcMask();
}
return _voxelInStructure;
}
void BoostMask::calcMask()
{
preprocessing();
voxelization();
generateMaskVoxelList();
_isUpToDate = true;
}
void BoostMask::preprocessing()
{
rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector();
//Convert world coordinate polygons to the polygons with geometry coordinate
rttb::PolygonSequenceType geometryCoordinatePolygonVector;
rttb::PolygonSequenceType::iterator it;
rttb::DoubleVoxelGridIndex3D globalMaxGridIndex(std::numeric_limits<double>::min(),
std::numeric_limits<double>::min(), std::numeric_limits<double>::min());
rttb::DoubleVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(),
_geometricInfo->getNumRows(), 0);
for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it)
{
PolygonType rttbPolygon = *it;
PolygonType geometryCoordinatePolygon;
//1. convert polygon to geometry coordinate polygons
//2. calculate global min/max
//3. check if polygon is planar
if (!preprocessingPolygon(rttbPolygon, geometryCoordinatePolygon, globalMinGridIndex,
globalMaxGridIndex, errorConstant))
{
throw rttb::core::Exception("TiltedMaskPlaneException");
}
geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon);
}
rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) ),
GridIndexType(globalMinGridIndex(1) ), GridIndexType(globalMinGridIndex(2) ));
rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) ),
GridIndexType(globalMaxGridIndex(1) ), GridIndexType(globalMaxGridIndex(2) ));
_globalBoundingBox.push_back(minIndex);
_globalBoundingBox.push_back(maxIndex);
//convert rttb polygon sequence to a map of z index and a vector of boost ring 2d (without holes)
_ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector);
}
void BoostMask::voxelization()
{
if (_globalBoundingBox.size() < 2)
{
throw rttb::core::InvalidParameterException("Bounding box calculation failed! ");
}
BoostRingMap::iterator itMap;
size_t mapSizeInAThread = _ringMap.size() / _numberOfThreads;
unsigned int count = 0;
unsigned int countThread = 0;
BoostPolygonMap polygonMap;
std::vector<BoostPolygonMap> polygonMapVector;
//check donut and convert to a map of z index and a vector of boost polygon 2d (with or without holes)
for (itMap = _ringMap.begin(); itMap != _ringMap.end(); ++itMap)
{
//the vector of all boost 2d polygons with the same z grid index(donut polygon is accepted).
BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second);
if (count == mapSizeInAThread && countThread < (_numberOfThreads - 1))
{
polygonMapVector.push_back(polygonMap);
polygonMap.clear();
count = 0;
countThread++;
}
polygonMap.insert(std::pair<double, BoostPolygonVector>((*itMap).first,
polygonVector));
count++;
}
_voxelizationMap = ::boost::make_shared<std::map<double, BoostArray2DPointer> >();
polygonMapVector.push_back(polygonMap); //insert the last one
//generate voxelization map, multi-threading
::boost::thread_group threads;
auto aMutex = ::boost::make_shared<::boost::shared_mutex>();
for (const auto & i : polygonMapVector)
{
BoostMaskVoxelizationThread t(i, _globalBoundingBox,
_voxelizationMap, aMutex, _strict);
threads.create_thread(t);
}
threads.join_all();
}
void BoostMask::generateMaskVoxelList()
{
if (_globalBoundingBox.size() < 2)
{
throw rttb::core::InvalidParameterException("Bounding box calculation failed! ");
}
//check homogeneous of the voxelization plane (the contours plane)
if (!calcVoxelizationThickness(_voxelizationThickness))
{
throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneous!");
}
::boost::thread_group threads;
auto aMutex = ::boost::make_shared<::boost::shared_mutex>();
unsigned int sliceNumberInAThread = _geometricInfo->getNumSlices() / _numberOfThreads;
//generate mask voxel list, multi-threading
for (unsigned int i = 0; i < _numberOfThreads; ++i)
{
unsigned int beginSlice = i * sliceNumberInAThread;
unsigned int endSlice;
if (i < _numberOfThreads - 1)
{
endSlice = (i + 1) * sliceNumberInAThread;
}
else
{
endSlice = _geometricInfo->getNumSlices();
}
BoostMaskGenerateMaskVoxelListThread t(_globalBoundingBox, _geometricInfo, _voxelizationMap,
_voxelizationThickness, beginSlice, endSlice,
_voxelInStructure, aMutex);
threads.create_thread(t);
}
threads.join_all();
}
bool BoostMask::preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon,
rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum,
rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const
{
double minZ = _geometricInfo->getNumSlices();
double maxZ = 0.0;
for (auto worldCoordinatePoint : aRTTBPolygon)
{
//convert to geometry coordinate polygon
rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint;
_geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, geometryCoordinatePoint);
geometryCoordinatePolygon.push_back(geometryCoordinatePoint);
//calculate the current global min/max
//min and max for x
if (geometryCoordinatePoint(0) < minimum(0))
{
minimum(0) = geometryCoordinatePoint(0);
}
if (geometryCoordinatePoint(0) > maximum(0))
{
maximum(0) = geometryCoordinatePoint(0);
}
//min and max for y
if (geometryCoordinatePoint(1) < minimum(1))
{
minimum(1) = geometryCoordinatePoint(1);
}
if (geometryCoordinatePoint(1) > maximum(1))
{
maximum(1) = geometryCoordinatePoint(1);
}
//min and max for z
if (geometryCoordinatePoint(2) < minimum(2))
{
minimum(2) = geometryCoordinatePoint(2);
}
if (geometryCoordinatePoint(2) > maximum(2))
{
maximum(2) = geometryCoordinatePoint(2);
}
//check planar
if (geometryCoordinatePoint(2) < minZ)
{
minZ = geometryCoordinatePoint(2);
}
if (geometryCoordinatePoint(2) > maxZ)
{
maxZ = geometryCoordinatePoint(2);
}
}
return (std::abs(maxZ - minZ) <= aErrorConstant);
}
BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType&
aRTTBPolygon) const
{
BoostMask::BoostRing2D polygon2D;
BoostPoint2D firstPoint;
for (unsigned int i = 0; i < aRTTBPolygon.size(); i++)
{
rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i);
BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]);
if (i == 0)
{
firstPoint = boostPoint;
}
::boost::geometry::append(polygon2D, boostPoint);
}
::boost::geometry::append(polygon2D, firstPoint);
return polygon2D;
}
BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap(
const rttb::PolygonSequenceType& aRTTBPolygonVector) const
{
rttb::PolygonSequenceType::const_iterator it;
BoostMask::BoostRingMap aRingMap;
for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it)
{
rttb::PolygonType rttbPolygon = *it;
double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon
bool isFirstZ = true;
if (!aRingMap.empty())
{
auto findIt = findNearestKey(aRingMap, zIndex, errorConstant);
//if the z index is found (same slice), add the polygon to vector
if (findIt != aRingMap.end())
{
//BoostRingVector ringVector = ;
(*findIt).second.push_back(convertRTTBPolygonToBoostRing(rttbPolygon));
isFirstZ = false;
}
}
//if it is the first z index in the map, insert vector with the polygon
if (isFirstZ)
{
BoostRingVector ringVector;
ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon));
aRingMap.insert(std::pair<double, BoostRingVector>(zIndex, ringVector));
}
}
return aRingMap;
}
BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap&
aBoostRingMap, double aIndex, double aErrorConstant) const
{
auto find = aBoostRingMap.find(aIndex);
//if find a key equivalent to aIndex, found
if (find != aBoostRingMap.end())
{
return find;
}
else
{
auto lowerBound = aBoostRingMap.lower_bound(aIndex);
//if all keys go before aIndex, check the last key
if (lowerBound == aBoostRingMap.end())
{
lowerBound = --aBoostRingMap.end();
}
//if the lower bound very close to aIndex, found
if (std::abs((*lowerBound).first - aIndex) <= aErrorConstant)
{
return lowerBound;
}
else
{
//if the lower bound is the beginning, not found
if (lowerBound == aBoostRingMap.begin())
{
return aBoostRingMap.end();
}
else
{
auto lowerBound1 = --lowerBound;//the key before the lower bound
//if the key before the lower bound very close to a Index, found
if (std::abs((*lowerBound1).first - aIndex) <= aErrorConstant)
{
return lowerBound1;
}
//else, not found
else
{
return aBoostRingMap.end();
}
}
}
}
}
BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector&
aRingVector) const
{
//check donut
BoostMask::BoostRingVector::const_iterator it1;
BoostMask::BoostRingVector::const_iterator it2;
BoostMask::BoostPolygonVector boostPolygonVector;
std::vector<unsigned int> donutIndexVector;//store the outer and inner ring index
BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon
//Get donut index and donut polygon
unsigned int index1 = 0;
for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++)
{
bool it1IsDonut = false;
//check if the ring is already determined as a donut
for (unsigned int i : donutIndexVector)
{
if (i == index1)
{
it1IsDonut = true;
break;
}
}
//if not jet, check now
if (!it1IsDonut)
{
bool it2IsDonut = false;
unsigned int index2 = 0;
for (it2 = aRingVector.begin(); it2 != aRingVector.end(); ++it2, index2++)
{
if (it2 != it1)
{
BoostMask::BoostPolygon2D polygon2D;
if (::boost::geometry::within(*it1, *it2))
{
::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon
::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring
::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring
it2IsDonut = true;
}
//if donut
else if (::boost::geometry::within(*it2, *it1))
{
::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon
::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring
::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring
it2IsDonut = true;
}
if (it2IsDonut)
{
donutIndexVector.push_back(index1);
donutIndexVector.push_back(index2);
donutVector.push_back(polygon2D);//store donut polygon
break;//Only store the first donut!
}
}
}
}
}
//Store no donut polygon to boostPolygonVector
index1 = 0;
for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++)
{
bool it1IsDonut = false;
//check if the ring is the outer or inner of a donut
for (unsigned int i : donutIndexVector)
{
if (i == index1)
{
it1IsDonut = true;
break;
}
}
if (!it1IsDonut)
{
BoostMask::BoostPolygon2D polygon2D;
::boost::geometry::append(polygon2D, *it1);
boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut
}
}
//Append donut polygon to boostPolygonVector
BoostMask::BoostPolygonVector::iterator itDonut;
for (itDonut = donutVector.begin(); itDonut != donutVector.end(); ++itDonut)
{
boostPolygonVector.push_back(*itDonut);//append donuts
}
return boostPolygonVector;
}
bool BoostMask::calcVoxelizationThickness(double& aThickness) const
{
if (_voxelizationMap->size() <= 1)
{
aThickness = 1;
return true;
}
double thickness = 0;
auto it = _voxelizationMap->cbegin();
auto it2 = ++_voxelizationMap->cbegin();
for (;
it != _voxelizationMap->cend() && it2 != _voxelizationMap->cend(); ++it, ++it2)
{
if (thickness == 0)
{
thickness = it2->first - it->first;
}
else
{
double curThickness = it2->first - it->first;
//if not homogeneous (leave out double imprecisions), return false
if (std::abs(thickness-curThickness)>errorConstant)
{
//return false;
std::cout << "Two polygons are far from each other?" << std::endl;
}
}
}
if (thickness != 0)
{
aThickness = thickness;
}
else
{
aThickness = 1;
}
return true;
}
}
}
}
diff --git a/code/masks/rttbBoostMask.h b/code/masks/rttbBoostMask.h
index cf994bb..1db206a 100644
--- a/code/masks/rttbBoostMask.h
+++ b/code/masks/rttbBoostMask.h
@@ -1,202 +1,195 @@
// -----------------------------------------------------------------------
// RTToolbox - DKFZ radiotherapy quantitative evaluation library
//
// Copyright (c) German Cancer Research Center (DKFZ),
// Software development for Integrated Diagnostics and Therapy (SIDT).
// ALL RIGHTS RESERVED.
// See rttbCopyright.txt or
// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notices for more information.
//
//------------------------------------------------------------------------
-/*!
-// @file
-// @version $Revision$ (last changed revision)
-// @date $Date$ (last change date)
-// @author $Author$ (last changed by)
-*/
-
#ifndef __BOOST_MASK_R_H
#define __BOOST_MASK_R_H
#include "rttbBaseType.h"
#include "rttbStructure.h"
#include "rttbGeometricInfo.h"
#include "rttbMaskAccessorInterface.h"
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/multi_array.hpp>
namespace rttb
{
namespace masks
{
namespace boost
{
/*! @class BoostMask
* @brief Implementation of voxelization using boost::geometry.
* @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection.
* (A structure without self intersection means all contours of the structure have no self intersection, and
* the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.)
* If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored
* and the mask will be calculated, however, it may cause errors in the mask results.
*/
class BoostMask
{
public:
using MaskVoxelList = core::MaskAccessorInterface::MaskVoxelList;
using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer;
/*! @brief Constructor
* @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is nullptr
* @param aDoseGeoInfo the GeometricInfo
* @param aStructure the structure set
* @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection.
* @param numberOfThreads number of threads used for voxelization. default value 0 means automatic detection, using the number of Hardware thread/cores
* @exception InvalidParameterException thrown if strict is true and the structure has self intersections
*/
BoostMask(core::GeometricInfo::Pointer aDoseGeoInfo, core::Structure::Pointer aStructure,
bool strict = true, unsigned int numberOfThreads = 0);
/*! @brief Generate mask and return the voxels in the mask
* @exception rttb::core::InvalidParameterException thrown if the structure has self intersections
*/
MaskVoxelListPointer getRelevantVoxelVector();
private:
using BoostPoint2D = ::boost::geometry::model::d2::point_xy<double>;
using BoostPolygon2D = ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy<double> >;
using BoostRing2D = ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy<double> >;
using BoostRingVector = std::vector<BoostRing2D>;//polygon without holes
using BoostPolygonVector = std::vector<BoostPolygon2D>;//polygon with or without holes
using VoxelIndexVector = std::vector<rttb::VoxelGridIndex3D>;
typedef std::map<double, BoostPolygonVector>
BoostPolygonMap;//map of the z index with the vector of boost 2d polygon
typedef std::map<double, BoostRingVector>
BoostRingMap;//map of the z index with the vector of boost 2d ring
typedef ::boost::multi_array<double, 2> BoostArray2D;
using BoostArray2DPointer = ::boost::shared_ptr<BoostArray2D>;
typedef ::boost::shared_ptr<std::map<double, BoostArray2DPointer> > BoostArrayMapPointer;
core::GeometricInfo::Pointer _geometricInfo;
core::Structure::Pointer _structure;
bool _strict;
/*! @brief The number of threads
*/
unsigned int _numberOfThreads;
//@brief The thickness of the voxelization plane (the contour plane), in double dose grid index
//@details for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2
double _voxelizationThickness;
//@brief vector of the MaskVoxel inside the structure
MaskVoxelListPointer _voxelInStructure;
/*! @brief The map of z index and a vector of boost ring 2d (without holes)
* @details Key: the double z grid index
* Value: the vector of boost ring 2d (without holes)
*/
BoostRingMap _ringMap;
/*! @brief The min and max index of the global bounding box.
* @details The first index has the minimum for x/y/z of the global bounding box.
* The second index has the maximum for x/y/z of the global bounding index.
*/
VoxelIndexVector _globalBoundingBox;
/*! @brief The voxelization map
* @details key: the converted double z grid index of a contour plane
* value: the 2d mask, array[i][j] = the mask value of the position (i,j) in the global bounding box,
* i: 0 - (_globalBoundingBoxSize0-1), j: 0 - (_globalBoundingBoxSize1-1)
*/
BoostArrayMapPointer _voxelizationMap;
/*! @brief If the mask is up to date
*/
bool _isUpToDate;
/*! @brief Voxelization and generate mask
*/
void calcMask();
/*! @brief The preprocessing step, wich consists of the following logic and Sub setps:
* @details For all contours in a struct:
* 1) Transfer the contour polygons into boost::geometry structures
* 1a) Convert the contur points from world coordinates into geometry coordinates.
* 1b) get min and max for x/y/z of a contour
* 2) Tilt check: if difference of z_min and z_max is larger then a tolerance value -> there is a tilt. Throw rttb::TiltedMaskPlaneException.
* 3) Get struct-bounding-box: get x_min_struct, y_min_struct, x_max_struct, y_max_struct to define the bounding box that containes all contours of a struct in x-y-dimensions.
*/
void preprocessing();
/*! @brief The voxelization step, which computes the voxelization planes (in x/y) for all contours of an struct.
* @details For each contour (that is in the z-Range of the reference geometry) of the struct:
* 1) Allocate result array (voxelization plane) based on the bounding box (see Preprocessing Step 3)
* 2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry).
* 3) Add result Array (key is the z-Value of the contour)
*/
void voxelization();
/*! @brief mask voxel Generation step which transfers the voxelization planes into the (z-)geometry of the reference geometry.
* @details It consists of following Sub steps :
* For all "slices" in the reference geometry :
* 1) generate weight vector for all voxelization planes for a given z - value of a slice
* Iterate over the bounding box of a struct.For each voxel :
* 2) Compute weighted sum of all voxelization planes(use weight vector, step 1)
* 2a) If sum > 0 : Add mask voxel for the current x / y(inner Loop) and z value(outer Loop).
* 3) return mask voxel list.
*/
void generateMaskVoxelList();
/*! @brief Convert the rttb polygon with world coordinate to the rttb polygon with double geometry coordinate, calculate the current min/max
* and check if the polygon is planar
* @param minimum the current global minimum
* @param maximum the current global maximum
* @return Return true if the polygon is planar, which means that the minimal and maximal z-coordinate of the polygon is not larger than a error constant
*/
bool preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon,
rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum,
rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const;
/*! @brief Convert a rttb 3d polygon to a 2d boost ring*/
BoostRing2D convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const;
/*! @brief Convert a rttb 3d polygon to a map of z index with a vector of boost 2d ring, because of tilt check use the first z index of the polygon as the map key*/
BoostRingMap convertRTTBPolygonSequenceToBoostRingMap(const rttb::PolygonSequenceType&
aRTTBPolygonVector) const;
/*! @brief Find the key with error constant to aIndex
* @pre aBoostRingMap should not be empty
* @return Return aBoostRingMap.end() if the key is not found
*/
BoostMask::BoostRingMap::iterator findNearestKey(BoostMask::BoostRingMap& aBoostRingMap,
double aIndex, double aErrorConstant) const;
/*! @brief If 2 rings in the vector build a donut, convert the 2 rings to a donut polygon, other rings unchanged*/
BoostPolygonVector checkDonutAndConvert(const BoostRingVector& aRingVector) const;
/*! @brief Calculate the voxelization thickness.
Return false, if the voxelization plane is not homogeneous
*/
bool calcVoxelizationThickness(double& aThickness) const;
};
}
}
}
#endif
diff --git a/code/masks/rttbBoostMaskAccessor.cpp b/code/masks/rttbBoostMaskAccessor.cpp
index 0985637..231f1a7 100644
--- a/code/masks/rttbBoostMaskAccessor.cpp
+++ b/code/masks/rttbBoostMaskAccessor.cpp
@@ -1,161 +1,155 @@
// -----------------------------------------------------------------------
// RTToolbox - DKFZ radiotherapy quantitative evaluation library
//
// Copyright (c) German Cancer Research Center (DKFZ),
// Software development for Integrated Diagnostics and Therapy (SIDT).
// ALL RIGHTS RESERVED.
// See rttbCopyright.txt or
// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notices for more information.
//
//------------------------------------------------------------------------
-/*!
-// @file
-// @version $Revision$ (last changed revision)
-// @date $Date$ (last change date)
-// @author $Author$ (last changed by)
-*/
#include "rttbBoostMaskAccessor.h"
#include "rttbBoostMask.h"
#include <boost/make_shared.hpp>
#include <boost/uuid/uuid.hpp>
#include <boost/uuid/uuid_generators.hpp>
#include <boost/uuid/uuid_io.hpp>
namespace rttb
{
namespace masks
{
namespace boost
{
BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer,
const core::GeometricInfo& aGeometricInfo, bool strict)
: _spStructure(aStructurePointer), _geoInfo(aGeometricInfo), _strict(strict)
{
_spRelevantVoxelVector = MaskVoxelListPointer();
//generate new structure set uid
::boost::uuids::uuid id;
::boost::uuids::random_generator generator;
id = generator();
std::stringstream ss;
ss << id;
_maskUID = "BoostMask_" + ss.str();
}
BoostMaskAccessor::~BoostMaskAccessor()
= default;
void BoostMaskAccessor::updateMask()
{
MaskVoxelList newRelevantVoxelVector;
if (_spRelevantVoxelVector)
{
return; // already calculated
}
BoostMask mask(::boost::make_shared<core::GeometricInfo>(_geoInfo),
_spStructure, _strict);
_spRelevantVoxelVector = mask.getRelevantVoxelVector();
}
BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector()
{
// if not already generated start voxelization here
updateMask();
return _spRelevantVoxelVector;
}
BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector(
float lowerThreshold)
{
auto filteredVoxelVectorPointer = ::boost::make_shared<MaskVoxelList>();
updateMask();
// filter relevant voxels
auto it = _spRelevantVoxelVector->begin();
while (it != _spRelevantVoxelVector->end())
{
if ((*it).getRelevantVolumeFraction() > lowerThreshold)
{
filteredVoxelVectorPointer->push_back(*it);
}
++it;
}
// if mask calculation was not successful this is empty!
return filteredVoxelVectorPointer;
}
bool BoostMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const
{
//initialize return voxel
voxel.setRelevantVolumeFraction(0);
//check if ID is valid
if (!_geoInfo.validID(aID))
{
return false;
}
//determine how a given voxel on the dose grid is masked
if (_spRelevantVoxelVector)
{
auto it = _spRelevantVoxelVector->begin();
while (it != _spRelevantVoxelVector->end())
{
if ((*it).getVoxelGridID() == aID)
{
voxel = (*it);
return true;
}
++it;
}
}
// returns false if mask was not calculated without triggering calculation (otherwise not const!)
else
{
return false;
}
return false;
}
bool BoostMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const
{
//convert VoxelGridIndex3D to VoxelGridID
VoxelGridID aVoxelGridID;
if (_geoInfo.convert(aIndex, aVoxelGridID))
{
return getMaskAt(aVoxelGridID, voxel);
}
else
{
return false;
}
}
const core::GeometricInfo& BoostMaskAccessor::getGeometricInfo() const
{
return _geoInfo;
};
}
}
}
\ No newline at end of file
diff --git a/code/masks/rttbBoostMaskAccessor.h b/code/masks/rttbBoostMaskAccessor.h
index 83bf37c..5a02dd4 100644
--- a/code/masks/rttbBoostMaskAccessor.h
+++ b/code/masks/rttbBoostMaskAccessor.h
@@ -1,121 +1,116 @@
// -----------------------------------------------------------------------
// RTToolbox - DKFZ radiotherapy quantitative evaluation library
//
// Copyright (c) German Cancer Research Center (DKFZ),
// Software development for Integrated Diagnostics and Therapy (SIDT).
// ALL RIGHTS RESERVED.
// See rttbCopyright.txt or
// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notices for more information.
//
//------------------------------------------------------------------------
-/*!
-// @file
-// @version $Revision$ (last changed revision)
-// @date $Date$ (last change date)
-// @author $Author$ (last changed by)
-*/
+
#ifndef __BOOST_MASK_R_ACCESSOR__H
#define __BOOST_MASK_R_ACCESSOR__H
#include "rttbBaseType.h"
#include "rttbGeometricInfo.h"
#include "rttbMaskAccessorInterface.h"
#include "rttbStructure.h"
#include "RTTBBoostMaskExports.h"
namespace rttb
{
namespace masks
{
namespace boost
{
/*! @class BoostMaskAccessor
* @brief Using the voxelization based on boost::geometry and generate the mask accessor.
* @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection.
* (A structure without self intersection means all contours of the structure have no self intersection, and
* the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.)
* If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored
* and the mask will be calculated, however, it may cause errors in the mask results.
*/
class RTTBBoostMask_EXPORT BoostMaskAccessor : public core::MaskAccessorInterface
{
public:
using MaskVoxelList = core::MaskAccessorInterface::MaskVoxelList;
using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer;
using StructTypePointer = core::Structure::Pointer;
private:
StructTypePointer _spStructure;
core::GeometricInfo _geoInfo;
bool _strict;
/*! vector containing list of mask voxels*/
MaskVoxelListPointer _spRelevantVoxelVector;
IDType _maskUID;
public:
/*! @brief Constructor with a structure pointer and a geometric info pointer
* @param aStructurePointer smart pointer of the structure
* @param aGeometricInfo smart pointer of the geometricInfo of the dose
* @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection.
* @exception InvalidParameterException thrown if strict is true and the structure has self intersections
*/
BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo,
bool strict = true);
/*! @brief destructor*/
~BoostMaskAccessor() override;
/*! @brief voxelization of the given structures using boost algorithms*/
void updateMask() override;
/*! @brief get vector containing all relevant voxels that are inside the given structure*/
MaskVoxelListPointer getRelevantVoxelVector() override;
/*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/
MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold) override;
/*!@brief determine how a given voxel on the dose grid is masked
* @param aID ID of the voxel in grid.
* @param voxel Reference to the voxel.
* @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined.
* The relevant volume fraction will be set to zero.
* @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/
bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const override;
/*!@brief determine how a given voxel on the dose grid is masked
* @param aIndex 3d index of the voxel in grid.
* @param voxel Reference to the voxel.
* @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/
bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const override;
/*! @brief give access to GeometricInfo*/
const core::GeometricInfo& getGeometricInfo() const override;
/* @ brief is true if dose is on a homogeneous grid
* @remark Inhomogeneous grids are not supported at the moment, but if they will be supported in the future the interface does not need to change.*/
bool isGridHomogeneous() const override
{
return true;
};
IDType getMaskUID() const override
{
return _maskUID;
};
};
}
}
}
#endif
diff --git a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp
index 02c39ce..479c086 100644
--- a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp
+++ b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp
@@ -1,156 +1,149 @@
// -----------------------------------------------------------------------
// RTToolbox - DKFZ radiotherapy quantitative evaluation library
//
// Copyright (c) German Cancer Research Center (DKFZ),
// Software development for Integrated Diagnostics and Therapy (SIDT).
// ALL RIGHTS RESERVED.
// See rttbCopyright.txt or
// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notices for more information.
//
//------------------------------------------------------------------------
-/*!
-// @file
-// @version $Revision: 1127 $ (last changed revision)
-// @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date)
-// @author $Author: hentsch $ (last changed by)
-*/
#include "rttbBoostMaskGenerateMaskVoxelListThread.h"
#include "rttbInvalidParameterException.h"
#include <boost/thread.hpp>
-
namespace rttb
{
namespace masks
{
namespace boost
{
BoostMaskGenerateMaskVoxelListThread::BoostMaskGenerateMaskVoxelListThread(
const VoxelIndexVector& aGlobalBoundingBox,
core::GeometricInfo::Pointer aGeometricInfo,
BoostArrayMapPointer aVoxelizationMap,
double aVoxelizationThickness,
unsigned int aBeginSlice,
unsigned int aEndSlice,
MaskVoxelListPointer aMaskVoxelList,
::boost::shared_ptr<::boost::shared_mutex> aMutex) :
_globalBoundingBox(aGlobalBoundingBox), _geometricInfo(aGeometricInfo),
_voxelizationMap(aVoxelizationMap), _voxelizationThickness(aVoxelizationThickness),
_beginSlice(aBeginSlice), _endSlice(aEndSlice),
_resultMaskVoxelList(aMaskVoxelList), _mutex(aMutex)
{}
void BoostMaskGenerateMaskVoxelListThread::operator()()
{
rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0);
rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1);
unsigned int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1;
unsigned int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1;
std::vector<core::MaskVoxel> maskVoxelsInThread;
for (unsigned int indexZ = _beginSlice; indexZ < _endSlice; ++indexZ)
{
//calculate weight vector
std::map<double, double> weightVectorForZ;
calcWeightVector(indexZ, weightVectorForZ);
//For each x,y, calc sum of all voxelization plane, use weight vector
for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x)
{
for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y)
{
rttb::VoxelGridIndex3D currentIndex;
currentIndex[0] = x + minIndex[0];
currentIndex[1] = y + minIndex[1];
currentIndex[2] = indexZ;
rttb::VoxelGridID gridID;
_geometricInfo->convert(currentIndex, gridID);
double volumeFraction = 0;
auto it = weightVectorForZ.cbegin();
auto itMap = _voxelizationMap->cbegin();
for (; it != weightVectorForZ.cend()
&& itMap != _voxelizationMap->cend(); ++it, ++itMap)
{
double weight = it->second;
if (weight > 0){
BoostArray2DPointer voxelizationArray = itMap->second;
//calc sum of all voxelization plane, use weight
volumeFraction += (*voxelizationArray)[x][y] * weight;
}
}
if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant)
{
volumeFraction = 1;
}
else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant)
{
throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!");
}
//insert mask voxel if volumeFraction > 0
if (volumeFraction > 0)
{
core::MaskVoxel maskVoxelPtr = core::MaskVoxel(gridID, volumeFraction);
maskVoxelsInThread.push_back(maskVoxelPtr);
}
}
}
}
::boost::unique_lock<::boost::shared_mutex> lock(*_mutex);
_resultMaskVoxelList->insert(_resultMaskVoxelList->end(), maskVoxelsInThread.begin(), maskVoxelsInThread.end());
}
void BoostMaskGenerateMaskVoxelListThread::calcWeightVector(const rttb::VoxelGridID& aIndexZ,
std::map<double, double>& weightVector) const
{
double indexZMin = aIndexZ - 0.5;
double indexZMax = aIndexZ + 0.5;
for (auto & it : *_voxelizationMap)
{
double voxelizationPlaneIndexMin = it.first - 0.5 * _voxelizationThickness;
double voxelizationPlaneIndexMax = it.first + 0.5 * _voxelizationThickness;
double weight = 0;
if ((voxelizationPlaneIndexMin < indexZMin) && (voxelizationPlaneIndexMax > indexZMin))
{
if (voxelizationPlaneIndexMax < indexZMax)
{
weight = voxelizationPlaneIndexMax - indexZMin;
}
else
{
weight = 1;
}
}
else if ((voxelizationPlaneIndexMin >= indexZMin) && (voxelizationPlaneIndexMin < indexZMax))
{
if (voxelizationPlaneIndexMax < indexZMax)
{
weight = _voxelizationThickness;
}
else
{
weight = indexZMax - voxelizationPlaneIndexMin;
}
}
weightVector.insert(std::pair<double, double>(it.first, weight));
}
}
}
}
}
diff --git a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h
index 76dc429..822a6cf 100644
--- a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h
+++ b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h
@@ -1,89 +1,82 @@
// -----------------------------------------------------------------------
// RTToolbox - DKFZ radiotherapy quantitative evaluation library
//
// Copyright (c) German Cancer Research Center (DKFZ),
// Software development for Integrated Diagnostics and Therapy (SIDT).
// ALL RIGHTS RESERVED.
// See rttbCopyright.txt or
// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notices for more information.
//
//------------------------------------------------------------------------
-/*!
-// @file
-// @version $Revision: 1127 $ (last changed revision)
-// @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date)
-// @author $Author: hentsch $ (last changed by)
-*/
-
#ifndef __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H
#define __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H
#include "rttbBaseType.h"
#include "rttbGeometricInfo.h"
#include "rttbMaskAccessorInterface.h"
#include <boost/multi_array.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/thread/shared_mutex.hpp>
namespace rttb
{
namespace masks
{
namespace boost
{
/*! @class BoostMaskGenerateMaskVoxelListThread
*
*/
class BoostMaskGenerateMaskVoxelListThread
{
public:
typedef ::boost::multi_array<double, 2> BoostArray2D;
using BoostArray2DPointer = ::boost::shared_ptr<BoostArray2D>;
typedef ::boost::shared_ptr<std::map<double, BoostArray2DPointer> > BoostArrayMapPointer;
using VoxelIndexVector = std::vector<rttb::VoxelGridIndex3D>;
using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer;
BoostMaskGenerateMaskVoxelListThread(const VoxelIndexVector& aGlobalBoundingBox,
core::GeometricInfo::Pointer aGeometricInfo,
BoostArrayMapPointer aVoxelizationMap,
double aVoxelizationThickness,
unsigned int aBeginSlice,
unsigned int aEndSlice,
MaskVoxelListPointer aMaskVoxelList, ::boost::shared_ptr<::boost::shared_mutex> aMutex);
void operator()();
private:
VoxelIndexVector _globalBoundingBox;
core::GeometricInfo::Pointer _geometricInfo;
BoostArrayMapPointer _voxelizationMap;
//(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2)
double _voxelizationThickness;
unsigned int _beginSlice;
/*! @brief _endSlice is the first index not to be processed (like a end iterator)
*/
unsigned int _endSlice;
MaskVoxelListPointer _resultMaskVoxelList;
::boost::shared_ptr<::boost::shared_mutex> _mutex;
/*! @brief For each dose grid index z, calculate the weight vector for each structure contour
*/
void calcWeightVector(const rttb::VoxelGridID& aIndexZ,
std::map<double, double>& weightVector) const;
};
}
}
}
#endif
\ No newline at end of file
diff --git a/code/masks/rttbBoostMaskVoxelizationThread.cpp b/code/masks/rttbBoostMaskVoxelizationThread.cpp
index d184883..693898b 100644
--- a/code/masks/rttbBoostMaskVoxelizationThread.cpp
+++ b/code/masks/rttbBoostMaskVoxelizationThread.cpp
@@ -1,174 +1,168 @@
// -----------------------------------------------------------------------
// RTToolbox - DKFZ radiotherapy quantitative evaluation library
//
// Copyright (c) German Cancer Research Center (DKFZ),
// Software development for Integrated Diagnostics and Therapy (SIDT).
// ALL RIGHTS RESERVED.
// See rttbCopyright.txt or
// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html
//
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notices for more information.
//
//------------------------------------------------------------------------
-/*!
-// @file
-// @version $Revision: 1127 $ (last changed revision)
-// @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date)
-// @author $Author: hentsch $ (last changed by)
-*/
#include "rttbBoostMaskVoxelizationThread.h"
#include "rttbInvalidParameterException.h"
#include <boost/geometry.hpp>
#include <boost/thread.hpp>
#include <boost/make_shared.hpp>
namespace rttb
{
namespace masks
{
namespace boost
{
BoostMaskVoxelizationThread::BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap,
const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex, bool strict) : _geometryCoordinateBoostPolygonMap(APolygonMap),
_globalBoundingBox(aGlobalBoundingBox), _resultVoxelization(anArrayMap), _mutex(aMutex), _strict(strict)
{
}
void BoostMaskVoxelizationThread::operator()()
{
rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0);
rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1);
const unsigned int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1;
const unsigned int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1;
std::map<double, ::boost::shared_ptr<BoostArray2D> > voxelizationMapInThread;
for (auto & it : _geometryCoordinateBoostPolygonMap)
{
BoostArray2D maskArray(::boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]);
BoostPolygonVector boostPolygonVec = it.second;
for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x)
{
for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y)
{
rttb::VoxelGridIndex3D currentIndex;
currentIndex[0] = x + minIndex[0];
currentIndex[1] = y + minIndex[1];
currentIndex[2] = 0;
//Get intersection polygons of the dose voxel and the structure
BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec);
//Calc areas of all intersection polygons
double volumeFraction = calcArea(polygons);
volumeFraction = correctForErrorAndStrictness(volumeFraction, _strict);
if (volumeFraction < 0 || volumeFraction > 1 )
{
throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!");
}
maskArray[x][y] = volumeFraction;
}
}
voxelizationMapInThread.insert(std::pair<double, BoostArray2DPointer>(it.first, ::boost::make_shared<BoostArray2D>(maskArray)));
}
//insert gathered values into voxelization map
::boost::unique_lock<::boost::shared_mutex> lock(*_mutex);
_resultVoxelization->insert(voxelizationMapInThread.begin(), voxelizationMapInThread.end());
}
/*Get intersection polygons of the contour and a voxel polygon*/
BoostMaskVoxelizationThread::BoostPolygonDeque BoostMaskVoxelizationThread::getIntersections(
const rttb::VoxelGridIndex3D&
aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons)
{
BoostMaskVoxelizationThread::BoostPolygonDeque polygonDeque;
BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D);
::boost::geometry::correct(voxelPolygon);
BoostPolygonVector::const_iterator it;
for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it)
{
BoostPolygon2D contour = *it;
::boost::geometry::correct(contour);
BoostPolygonDeque intersection;
::boost::geometry::intersection(voxelPolygon, contour, intersection);
polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end());
}
return polygonDeque;
}
BoostMaskVoxelizationThread::BoostRing2D BoostMaskVoxelizationThread::get2DContour(
const rttb::VoxelGridIndex3D& aVoxelGrid3D)
{
BoostRing2D polygon;
BoostPoint2D point1(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] - 0.5);
::boost::geometry::append(polygon, point1);
BoostPoint2D point2(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] - 0.5);
::boost::geometry::append(polygon, point2);
BoostPoint2D point3(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] + 0.5);
::boost::geometry::append(polygon, point3);
BoostPoint2D point4(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] + 0.5);
::boost::geometry::append(polygon, point4);
::boost::geometry::append(polygon, point1);
return polygon;
}
/*Calculate the intersection area*/
double BoostMaskVoxelizationThread::calcArea(const BoostPolygonDeque& aPolygonDeque)
{
double area = 0;
BoostPolygonDeque::const_iterator it;
for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it)
{
area += ::boost::geometry::area(*it);
}
return area;
}
double BoostMaskVoxelizationThread::correctForErrorAndStrictness(double volumeFraction, bool strict) const
{
if (strict){
if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant)
{
volumeFraction = 1;
}
}
else {
if (volumeFraction > 1){
volumeFraction = 1;
}
else if (volumeFraction < 0){
volumeFraction = 0;
}
}
return volumeFraction;
}
}
}
}
diff --git a/code/masks/rttbBoostMaskVoxelizationThread.h b/code/masks/rttbBoostMaskVoxelizationThread.h
index 3ffd14c..b31079f 100644
--- a/code/masks/rttbBoostMaskVoxelizationThread.h
+++ b/code/masks/rttbBoostMaskVoxelizationThread.h
@@ -1,109 +1,102 @@
// -----------------------------------------------------------------------
// 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 <deque>
#include "rttbBaseType.h"
#include <boost/multi_array.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/thread/shared_mutex.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
namespace rttb
{
namespace masks
{
namespace boost
{
/*! @class BoostMaskGenerateMaskVoxelListThread
*
*/
class BoostMaskVoxelizationThread
{
public:
using BoostPolygon2D = ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy<double> >;
using BoostPolygonVector = std::vector<BoostPolygon2D>;//polygon with or without holes
typedef std::map<double, BoostPolygonVector>
BoostPolygonMap;//map of the z index with the vector of boost 2d polygon
using VoxelIndexVector = std::vector<rttb::VoxelGridIndex3D>;
typedef ::boost::multi_array<double, 2> BoostArray2D;
using BoostArray2DPointer = ::boost::shared_ptr<BoostArray2D>;
typedef ::boost::shared_ptr<std::map<double, BoostArray2DPointer > > BoostArrayMapPointer;
/*! @brief Constructor
* @param aMutex a mutex for thread-safe handling of the _resultVoxelization
* @param strict true means that volumeFractions of <0 and >1 are NOT corrected. Otherwise, they are automatically corrected to 0 or 1, respectively.
*/
BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap,
const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex, bool strict);
void operator()();
private:
using BoostPolygonDeque = std::deque<BoostPolygon2D>;
using BoostRing2D = ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy<double> >;
using BoostPoint2D = ::boost::geometry::model::d2::point_xy<double>;
BoostPolygonMap _geometryCoordinateBoostPolygonMap;
VoxelIndexVector _globalBoundingBox;
BoostArrayMapPointer _resultVoxelization;
::boost::shared_ptr<::boost::shared_mutex> _mutex;
bool _strict;
/*! @brief Get intersection polygons of the contour and a voxel polygon
* @param aVoxelIndex3D The 3d grid index of the voxel
* @param intersectionSlicePolygons The polygons of the slice intersecting the voxel
* @return Return all intersection polygons of the structure and the voxel
*/
static BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D,
const BoostPolygonVector& intersectionSlicePolygons);
/*! @brief Get the voxel 2d contour polygon in geometry coordinate*/
static BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D);
/*! @brief Calculate the area of all polygons
* @param aPolygonDeque The deque of polygons
* @return Return the area of all polygons
*/
static double calcArea(const BoostPolygonDeque& aPolygonDeque);
/*! @brief Corrects the volumeFraction
* @details the volume fraction is corrected in case of strict=true. Otherwise, it's only corrected for double imprecision
* @return The corrected volumeFraction
*/
double correctForErrorAndStrictness(double volumeFraction, bool strict) const;
};
}
}
}
#endif
\ No newline at end of file
diff --git a/code/masks/rttbGenericMutableMaskAccessor.cpp b/code/masks/rttbGenericMutableMaskAccessor.cpp
index 3e0fd97..16f2eac 100644
--- a/code/masks/rttbGenericMutableMaskAccessor.cpp
+++ b/code/masks/rttbGenericMutableMaskAccessor.cpp
@@ -1,179 +1,173 @@
// -----------------------------------------------------------------------
// 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 "rttbGenericMutableMaskAccessor.h"
#include "rttbMaskVoxel.h"
#include "rttbNullPointerException.h"
#include <boost/make_shared.hpp>
#include <boost/uuid/uuid.hpp>
#include <boost/uuid/uuid_generators.hpp>
#include <boost/uuid/uuid_io.hpp>
namespace rttb
{
namespace masks
{
GenericMutableMaskAccessor::GenericMutableMaskAccessor(const core::GeometricInfo& aGeometricInfo) :
_geoInfo(aGeometricInfo), _spRelevantVoxelVector(MaskVoxelListPointer())
{
//generate new structure set uid
boost::uuids::uuid id;
boost::uuids::random_generator generator;
id = generator();
std::stringstream ss;
ss << id;
_maskUID = "GenericMutableMask_" + ss.str();
}
GenericMutableMaskAccessor::~GenericMutableMaskAccessor() = default;
void GenericMutableMaskAccessor::updateMask() {}
GenericMutableMaskAccessor::MaskVoxelListPointer
GenericMutableMaskAccessor::getRelevantVoxelVector()
{
return _spRelevantVoxelVector;
}
GenericMutableMaskAccessor::MaskVoxelListPointer GenericMutableMaskAccessor::getRelevantVoxelVector(
float lowerThreshold)
{
auto filteredVoxelVectorPointer = boost::make_shared<MaskVoxelList>();
// filter relevant voxels
auto it = _spRelevantVoxelVector->begin();
while (it != _spRelevantVoxelVector->end())
{
if ((*it).getRelevantVolumeFraction() > lowerThreshold)
{
filteredVoxelVectorPointer->push_back(*it);
}
++it;
}
// if mask calculation was not successful this is empty!
return filteredVoxelVectorPointer;
}
bool GenericMutableMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const
{
//initialize return voxel
voxel.setRelevantVolumeFraction(0);
//check if ID is valid
if (!_geoInfo.validID(aID))
{
return false;
}
//determine how a given voxel on the dose grid is masked
if (_spRelevantVoxelVector)
{
auto it = _spRelevantVoxelVector->begin();
while (it != _spRelevantVoxelVector->end())
{
if ((*it).getVoxelGridID() == aID)
{
voxel = (*it);
return true;
}
++it;
}
//aID is not in mask
voxel.setRelevantVolumeFraction(0);
}
return false;
}
bool GenericMutableMaskAccessor::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;
}
}
void GenericMutableMaskAccessor::setMaskAt(VoxelGridID aID, const core::MaskVoxel& voxel)
{
//check if ID is valid
if (!_geoInfo.validID(aID))
{
return;
}
//determine how a given voxel on the dose grid is masked
if (_spRelevantVoxelVector)
{
auto it = _spRelevantVoxelVector->begin();
while (it != _spRelevantVoxelVector->end())
{
if ((*it).getVoxelGridID() == aID)
{
(*it) = voxel;
return;
}
++it;
}
//not sID is not found in existing voxels
_spRelevantVoxelVector->push_back(voxel);
}
}
void GenericMutableMaskAccessor::setMaskAt(const VoxelGridIndex3D& aIndex,
const core::MaskVoxel& voxel)
{
//convert VoxelGridIndex3D to VoxelGridID
VoxelGridID aVoxelGridID;
if (_geoInfo.convert(aIndex, aVoxelGridID))
{
setMaskAt(aVoxelGridID, voxel);
}
}
void GenericMutableMaskAccessor::setRelevantVoxelVector(MaskVoxelListPointer aVoxelListPointer)
{
_spRelevantVoxelVector = MaskVoxelListPointer(aVoxelListPointer);
}
}
}
\ No newline at end of file
diff --git a/code/masks/rttbGenericMutableMaskAccessor.h b/code/masks/rttbGenericMutableMaskAccessor.h
index ed5e2cf..f85cb7d 100644
--- a/code/masks/rttbGenericMutableMaskAccessor.h
+++ b/code/masks/rttbGenericMutableMaskAccessor.h
@@ -1,107 +1,102 @@
// -----------------------------------------------------------------------
// 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 __GENERIC_MUTABLE_MASK_ACCESSOR_H
#define __GENERIC_MUTABLE_MASK_ACCESSOR_H
#include "rttbMutableMaskAccessorInterface.h"
#include "rttbBaseType.h"
#include "rttbGeometricInfo.h"
namespace rttb
{
namespace masks
{
/*! @class GenericMutableMaskAccessor
@brief Default implementation of MutableMaskAccessorInterface.
@see MutableMaskAccessorInterface
*/
class GenericMutableMaskAccessor: public core::MutableMaskAccessorInterface
{
public:
using MaskVoxelList = core::MutableMaskAccessorInterface::MaskVoxelList;
using MaskVoxelListPointer = core::MutableMaskAccessorInterface::MaskVoxelListPointer;
private:
core::GeometricInfo _geoInfo;
/*! vector containing list of mask voxels*/
MaskVoxelListPointer _spRelevantVoxelVector;
IDType _maskUID;
GenericMutableMaskAccessor(const
GenericMutableMaskAccessor&) = delete; //not implemented on purpose -> non-copyable
GenericMutableMaskAccessor& operator=(const
GenericMutableMaskAccessor&) = delete;//not implemented on purpose -> non-copyable
public:
~GenericMutableMaskAccessor() override;
GenericMutableMaskAccessor(const core::GeometricInfo& aGeometricInfo);
/*! @brief initialize mask structure if _spRelevantVoxelVector was not previously initialized*/
void updateMask() override;
/*! @brief get vector containing all relevant voxels that are inside the given structure*/
MaskVoxelListPointer getRelevantVoxelVector() override;
/*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/
MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold) override;
/*!@brief determine how a given voxel on the dose grid is masked
* @param aID ID of the voxel in grid.
* @param voxel Reference to the voxel.
* @post after a valid call voxel contains 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 if the voxel exists and therefore if parameter voxel contains valid values.*/
bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const override;
bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const override;
/* @ brief is true if dose is on a homogeneous grid */
// Inhomogeneous grids are not supported at the moment, but if they will
// be supported in the future the interface does not need to change.
bool isGridHomogeneous() const override
{
return true;
};
/*! @brief give access to GeometricInfo*/
inline const core::GeometricInfo& getGeometricInfo() const override
{
return _geoInfo;
};
IDType getMaskUID() const override
{
return _maskUID;
};
void setMaskAt(VoxelGridID aID, const core::MaskVoxel& voxel) override;
void setMaskAt(const VoxelGridIndex3D& gridIndex, const core::MaskVoxel& voxel) override;
void setRelevantVoxelVector(MaskVoxelListPointer aVoxelListPointer) override;
};
}
}
#endif
\ No newline at end of file
diff --git a/testing/core/CreateTestStructure.cpp b/testing/core/CreateTestStructure.cpp
index 6c81a22..d310ac4 100644
--- a/testing/core/CreateTestStructure.cpp
+++ b/testing/core/CreateTestStructure.cpp
@@ -1,688 +1,681 @@
// -----------------------------------------------------------------------
// 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 <math.h>
#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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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);
WorldCoordinate3D p;
p(0) = p1.x() ;
p(1) = p1.y() ;
//This can be done directly for x/y coordinates, but not for z. Thus, we do it in this function.
p(2) = p1.z() - _geoInfo.getSpacing().z() / 2;
polygon.push_back(p);
}
return polygon;
}
PolygonType CreateTestStructure::createPolygonCenterOnPlaneCenter(std::vector<VoxelGridIndex2D> 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);
}
return polygon;
}
PolygonType CreateTestStructure::createPolygonBetweenUpperLeftAndCenter(
std::vector<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D>
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<VoxelGridIndex2D>
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<VoxelGridIndex2D>
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 j = 0; j < 1000; j++)
{
wp_intermediate(0) = wp1.x() + j * diffX;
wp_intermediate(1) = wp1.y() + j * diffY;
polygon.push_back(wp_intermediate);
}
}
std::cout << std::endl;
return polygon;
}
PolygonType CreateTestStructure::createPolygonCircle(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber)
{
PolygonType polygon;
if (aVoxelVector.size() > 0)
{
unsigned 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 (unsigned int j = 0; j <= 1000; j++)
{
double y_offset = sqrt((radius * radius) - ((frac_radius * j) * (frac_radius * j)));
wp_intermediate(0) = wp1.x() + frac_radius * j;
wp_intermediate(1) = wp1.y() - (y_offset * correct_y);
polygon.push_back(wp_intermediate);
}
for (unsigned int j = 1000; j <= 1000; j--)
{
double y_offset = sqrt((radius * radius) - ((frac_radius * j) * (frac_radius * j)));
wp_intermediate(0) = wp1.x() + frac_radius * j;
wp_intermediate(1) = wp1.y() + (y_offset * correct_y);
polygon.push_back(wp_intermediate);
}
for (unsigned int j = 0; j <= 1000; j++)
{
double y_offset = sqrt((radius * radius) - ((frac_radius * j) * (frac_radius * j)));
wp_intermediate(0) = wp1.x() - frac_radius * j;
wp_intermediate(1) = wp1.y() + y_offset * correct_y;
polygon.push_back(wp_intermediate);
}
for (unsigned int j = 1000; j <= 1000; j--)
{
double y_offset = sqrt((radius * radius) - ((frac_radius * j) * (frac_radius * j)));
wp_intermediate(0) = wp1.x() + frac_radius * j;
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<VoxelGridIndex2D> 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<VoxelGridIndex2D>
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<VoxelGridIndex2D>
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 627c5cf..cb0631f 100644
--- a/testing/core/CreateTestStructure.h
+++ b/testing/core/CreateTestStructure.h
@@ -1,81 +1,75 @@
// -----------------------------------------------------------------------
// 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<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
/*The z coordinate of the polygon is the upper side of the voxel on the slice*/
PolygonType createPolygonCenter(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
/*The z coordinate of the polygon is the center of the voxel on the slice*/
PolygonType createPolygonCenterOnPlaneCenter(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
PolygonType createPolygonBetweenUpperLeftAndCenter(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
PolygonType createPolygonBetweenLowerRightAndCenter(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
PolygonType createPolygonIntermediatePoints(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
PolygonType createPolygonUpperCenter(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
PolygonType createPolygonLeftEdgeMiddle(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
PolygonType createPolygonCircle(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
PolygonType createStructureSeveralSectionsInsideOneVoxelA(std::vector<VoxelGridIndex2D>
aVoxelVector, GridIndexType sliceNumber);
PolygonType createStructureSelfTouchingA(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
PolygonType createStructureSelfTouchingB(std::vector<VoxelGridIndex2D> aVoxelVector,
GridIndexType sliceNumber);
};
}//testing
}//rttb
diff --git a/testing/core/DummyStructure.cpp b/testing/core/DummyStructure.cpp
index 4a496a9..3345ff0 100644
--- a/testing/core/DummyStructure.cpp
+++ b/testing/core/DummyStructure.cpp
@@ -1,748 +1,742 @@
// -----------------------------------------------------------------------
// 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 <math.h>
#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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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 (unsigned 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::CreateRectangularStructureCentered(GridIndexType fromZPlane, GridIndexType toZPlane, GridIndexType fromZPlane2, GridIndexType toZPlane2)
{
CreateTestStructure another_cts = CreateTestStructure(_geoInfo);
std::vector<VoxelGridIndex2D> 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 (unsigned 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);
}
for (unsigned int i = fromZPlane2; i <= toZPlane2; ++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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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<VoxelGridIndex2D> 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 dea6e31..052283f 100644
--- a/testing/core/DummyStructure.h
+++ b/testing/core/DummyStructure.h
@@ -1,112 +1,105 @@
// -----------------------------------------------------------------------
// 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 from fromZPlane(included) to toZPlane(included)
*/
core::Structure CreateRectangularStructureCentered(GridIndexType fromZPlane, GridIndexType toZPlane);
/* Generate rectangular structure for the z slices from fromZPlane(included) to toZPlane(included) and from fromZPlane2(included) to toZPlane2(included)
*/
core::Structure CreateRectangularStructureCentered(GridIndexType fromZPlane, GridIndexType toZPlane, GridIndexType fromZPlane2, GridIndexType toZPlane2);
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/masks/boost/BoostMaskTest.cpp b/testing/masks/boost/BoostMaskTest.cpp
index 9582a2d..d72d06e 100644
--- a/testing/masks/boost/BoostMaskTest.cpp
+++ b/testing/masks/boost/BoostMaskTest.cpp
@@ -1,255 +1,248 @@
// -----------------------------------------------------------------------
// 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 <boost/make_shared.hpp>
#include <boost/shared_ptr.hpp>
#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 "rttbBoostMask.h"
#include "rttbBoostMaskAccessor.h"
#include "rttbInvalidParameterException.h"
-
namespace rttb
{
namespace testing
{
/*! @brief BoostMaskRedesignTest.
1) test constructors
2) test getRelevantVoxelVector
3) test getMaskAt
*/
int BoostMaskTest(int argc, char* argv[])
{
PREPARE_DEFAULT_TEST_REPORTING;
typedef core::Structure::Pointer StructTypePointer;
// generate test dose. geometric info: patient position (-25, -2, 35), center of the 1st.voxel
boost::shared_ptr<DummyDoseAccessor> spTestDoseAccessor =
boost::make_shared<DummyDoseAccessor>();
boost::shared_ptr<core::GeometricInfo> geometricPtr = boost::make_shared<core::GeometricInfo>
(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<core::Structure>(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<core::Structure>(myTestStruct2);
//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);
//(-20, 0.5, 48.75); (-12.5, 0.5, 48.75); (-12.5, 10.5, 48.75); (-20, 10.5, 48.75);
//(-20, 0.5, 51.25); (-12.5, 0.5, 51.25); (-12.5, 10.5, 51.25); (-20, 10.5, 51.25);
//to test calcVoxelizationThickness bug
core::Structure myTestStruct3 = myStructGenerator.CreateRectangularStructureCentered(2, 3, 6, 7);
StructTypePointer spMyStruct3 = boost::make_shared<core::Structure>(myTestStruct3);
//1) test BoostMask & 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,
spTestDoseAccessor->getGeometricInfo(), true));
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 errorConstantBoostMask = 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(), errorConstantBoostMask);
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(), errorConstantBoostMask);
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(), errorConstantBoostMask);
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(), errorConstantBoostMask);
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(), errorConstantBoostMask);
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());
rttb::masks::boost::BoostMask boostMask2 = rttb::masks::boost::BoostMask(
geometricPtr, spMyStruct2);
CHECK_NO_THROW(boostMask2.getRelevantVoxelVector());
rttb::masks::boost::BoostMaskAccessor boostMaskAccessor2(spMyStruct2,
spTestDoseAccessor->getGeometricInfo(), true);
CHECK_NO_THROW(boostMaskAccessor2.getRelevantVoxelVector());
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(), errorConstantBoostMask);
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(), errorConstantBoostMask);
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(), errorConstantBoostMask);
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(), errorConstantBoostMask);
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
//3) test calcVoxelizationThickness bug
rttb::masks::boost::BoostMask boostMask3 = rttb::masks::boost::BoostMask(
geometricPtr, spMyStruct3);
CHECK_NO_THROW(boostMask3.getRelevantVoxelVector());
rttb::masks::boost::BoostMaskAccessor boostMaskAccessor3(spMyStruct3,
spTestDoseAccessor->getGeometricInfo(), true);
CHECK_NO_THROW(boostMaskAccessor3.getRelevantVoxelVector());
RETURN_AND_REPORT_TEST_SUCCESS;
}
}//testing
}//rttb
diff --git a/testing/masks/boost/rttbBoostMaskTests.cpp b/testing/masks/boost/rttbBoostMaskTests.cpp
index 4dfa247..887d994 100644
--- a/testing/masks/boost/rttbBoostMaskTests.cpp
+++ b/testing/masks/boost/rttbBoostMaskTests.cpp
@@ -1,63 +1,56 @@
// -----------------------------------------------------------------------
// 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(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;
}