diff --git a/code/core/rttbGeometricInfo.h b/code/core/rttbGeometricInfo.h index e6aa52e..2f1f6e8 100644 --- a/code/core/rttbGeometricInfo.h +++ b/code/core/rttbGeometricInfo.h @@ -1,189 +1,190 @@ // ----------------------------------------------------------------------- // 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 __GEOMETRIC_INFO_NEW_H #define __GEOMETRIC_INFO_NEW_H #include #include #include #include #include "rttbBaseType.h" namespace rttb { namespace core { /*! @brief GeometricInfo objects contain all the information required for transformations between voxel grid coordinates and world coordinates. Corresponding converter functions are also available. @note ITK Pixel Indexing used (http://www.itk.org/Doxygen45/html/classitk_1_1Image.html): The Index type reverses the order so that with Index[0] = col, Index[1] = row, Index[2] = slice. */ class GeometricInfo { private: WorldCoordinate3D _imagePositionPatient; OrientationMatrix _orientationMatrix; OrientationMatrix _invertedOrientationMatrix; SpacingVectorType3D _spacing; VoxelGridDimensionType _numberOfColumns; VoxelGridDimensionType _numberOfRows; VoxelGridDimensionType _numberOfFrames; /* @brief Matrix inversion routine. Uses lu_factorize and lu_substitute in uBLAS to invert a matrix http://savingyoutime.wordpress.com/2009/09/21/c-matrix-inversion-boostublas/ */ bool computeInvertOrientation(); public: /*! @brief Constructor, initializes orientation matrix, spacing vector and patient position with zeros. */ - GeometricInfo() : _orientationMatrix(0), _spacing(0), _imagePositionPatient(0), _numberOfFrames(0), _numberOfRows(0), + GeometricInfo() : _orientationMatrix(0), _spacing(0), _imagePositionPatient(0), _numberOfFrames(0), + _numberOfRows(0), _numberOfColumns(0) {} void setSpacing(const SpacingVectorType3D& aSpacingVector); const SpacingVectorType3D& getSpacing() const; void setImagePositionPatient(const WorldCoordinate3D& aImagePositionPatient); const WorldCoordinate3D& getImagePositionPatient() const; void setOrientationMatrix(const OrientationMatrix& anOrientationMatrix); const OrientationMatrix getOrientationMatrix() const { return _orientationMatrix; }; /*! @brief Returns the 1st row of the OrientationMatrix. @deprecated please use getOrientationMatrix() (x,0) instead*/ const ImageOrientation getImageOrientationRow() const; /*! @brief Returns the 2nd row of the OrientationMatrix. @deprecated please use getOrientationMatrix() (x,1) instead*/ const ImageOrientation getImageOrientationColumn() const; void setPixelSpacingRow(const GridVolumeType aValue); /*! @brief Returns the spacing in the x-dimension (rows) of the data grid. @deprecated please use getSpacing() (0) instead*/ const GridVolumeType getPixelSpacingRow() const; void setPixelSpacingColumn(const GridVolumeType aValue); /*! @brief Returns the spacing in the y-dimension (columns) of the data grid. @deprecated please use getSpacing() (1) instead*/ const GridVolumeType getPixelSpacingColumn() const; void setSliceThickness(const GridVolumeType aValue); /*! @brief Returns the spacing in the z-dimension (slices) of the data grid. @deprecated please use getSpacing() (2) instead*/ const GridVolumeType getSliceThickness() const; void setImageSize(const ImageSize& aSize); const ImageSize getImageSize() const; void setNumColumns(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumColumns() const; void setNumRows(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumRows() const; void setNumSlices(const VoxelGridDimensionType aValue); const VoxelGridDimensionType getNumSlices() const; /*! @brief determines equality of two GeometricInfo objects. */ friend bool operator==(const GeometricInfo& gInfo, const GeometricInfo& gInfo1); bool equalsAlmost(const GeometricInfo& another, double errorConstant = 1e-5) const; /*! @brief convert world coordinates to voxel grid index. The conversion of values is done even if the target index is not inside the given voxel grid (return false). If the target is inside the grid return true. */ bool worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, VoxelGridIndex3D& aIndex) const; /*! @brief convert world coordinates to double geometry coordinate. The world coordinate of the image position patient (center of the first voxel) will be convert to the double voxel (0.0, 0.0, 0.0) The conversion of values is done even if the target index is not inside the given voxel grid (return false). If the target is inside the grid return true. */ bool worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, - DoubleVoxelGridIndex3D& aIndex) const; + DoubleVoxelGridIndex3D& aIndex) const; - /*! @brief convert double geometry coordinate to world coordinates. The double voxel index (0.0, 0.0, 0.0) will be convert to the world coordinate of the image postion patient (center of the first voxel) + /*! @brief convert double geometry coordinate to world coordinates. The double voxel index (0.0, 0.0, 0.0) will be convert to the world coordinate of the image postion patient (center of the first voxel) The conversion of values is done even if the target is not inside the given voxel grid (return false). If the target is inside the voxel grid return true. */ bool geometryCoordinateToWorldCoordinate(const DoubleVoxelGridIndex3D& aIndex, - WorldCoordinate3D& aWorldCoordinate) const; + WorldCoordinate3D& aWorldCoordinate) const; - /*! @brief convert int voxel grid index to world coordinates. + /*! @brief convert int voxel grid index to world coordinates. It is the upper left corner of the voxel. The conversion of values is done even if the target is not inside the given voxel grid (return false). If the target is inside the voxel grid return true. */ bool indexToWorldCoordinate(const VoxelGridIndex3D& aIndex, - WorldCoordinate3D& aWorldCoordinate) const; + WorldCoordinate3D& aWorldCoordinate) const; /*! @brief check if a given voxel grid index is inside the given voxel grid.*/ bool isInside(const VoxelGridIndex3D& aIndex) const; /*! @brief check if a given world coordinate is inside the given voxel grid.*/ bool isInside(const WorldCoordinate3D& aWorldCoordinate); const GridSizeType getNumberOfVoxels() const; bool convert(const VoxelGridID& gridID, VoxelGridIndex3D& gridIndex) const; bool convert(const VoxelGridIndex3D& gridIndex, VoxelGridID& gridID) const; /*! @brief test if given ID is inside current dose grid */ bool validID(const VoxelGridID aID) const; /*! @brief test if given index is inside current dose grid */ bool validIndex(const VoxelGridIndex3D& aIndex) const; /*!@ brief generates string stream representation of the GeometricInfo object. */ friend std::ostream& operator<<(std::ostream& s, const GeometricInfo& anGeometricInfo); }; } } #endif \ No newline at end of file diff --git a/code/masks/boost/files.cmake b/code/masks/boost/files.cmake index 161bcf0..2e2b40c 100644 --- a/code/masks/boost/files.cmake +++ b/code/masks/boost/files.cmake @@ -1,11 +1,13 @@ SET(CPP_FILES rttbBoostMask.cpp rttbBoostMaskAccessor.cpp rttbBoostMaskRedesign.cpp + rttbBoostMaskRedesignAccessor.cpp ) SET(H_FILES rttbBoostMask.h rttbBoostMaskAccessor.h rttbBoostMaskRedesign.h + rttbBoostMaskRedesignAccessor.h ) diff --git a/code/masks/boost/rttbBoostMaskAccessor.cpp b/code/masks/boost/rttbBoostMaskAccessor.cpp index c02b583..9050fbb 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.cpp +++ b/code/masks/boost/rttbBoostMaskAccessor.cpp @@ -1,161 +1,164 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #include "rttbBoostMaskAccessor.h" #include "rttbBoostMask.h" #include #include #include #include namespace rttb { namespace masks { namespace boost { - BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, bool strict) + BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, + const core::GeometricInfo& aGeometricInfo, bool strict) : _spStructure(aStructurePointer), _geoInfo(aGeometricInfo), _strict(strict) { _spRelevantVoxelVector = MaskVoxelListPointer(); //generate new structure set uid ::boost::uuids::uuid id; ::boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _maskUID = "BoostMask_" + ss.str(); } BoostMaskAccessor::~BoostMaskAccessor() { }; void BoostMaskAccessor::updateMask() { MaskVoxelList newRelevantVoxelVector; if (_spRelevantVoxelVector) { return; // already calculated } - BoostMask mask(::boost::make_shared(_geoInfo), _spStructure, _strict); + BoostMask mask(::boost::make_shared(_geoInfo), + _spStructure, _strict); _spRelevantVoxelVector = mask.getRelevantVoxelVector(); return; } BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector() { // if not already generated start voxelization here updateMask(); return _spRelevantVoxelVector; } - BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector(float lowerThreshold) + BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector( + float lowerThreshold) { MaskVoxelListPointer filteredVoxelVectorPointer(new MaskVoxelList); updateMask(); // filter relevant voxels BoostMaskAccessor::MaskVoxelList::iterator it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getRelevantVolumeFraction() > lowerThreshold) { filteredVoxelVectorPointer->push_back(*it); } ++it; } // if mask calculation was not successful this is empty! return filteredVoxelVectorPointer; } bool BoostMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const { //initialize return voxel voxel.setRelevantVolumeFraction(0); //check if ID is valid if (!_geoInfo.validID(aID)) { return false; } //determine how a given voxel on the dose grid is masked if (_spRelevantVoxelVector) { BoostMaskAccessor::MaskVoxelList::iterator it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getVoxelGridID() == aID) { voxel = (*it); return true; } ++it; } } // returns false if mask was not calculated without triggering calculation (otherwise not const!) else { return false; } return false; } bool BoostMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const { //convert VoxelGridIndex3D to VoxelGridID VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getMaskAt(aVoxelGridID, voxel); } else { return false; } } const core::GeometricInfo& BoostMaskAccessor::getGeometricInfo() const { return _geoInfo; }; } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskAccessor.h b/code/masks/boost/rttbBoostMaskAccessor.h index 0fdf506..96f1f42 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/boost/rttbBoostMaskAccessor.h @@ -1,127 +1,128 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __BOOST_MASK_ACCESSOR__H #define __BOOST_MASK_ACCESSOR__H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include "rttbGenericDoseIterator.h" #include "rttbStructure.h" #include namespace rttb { namespace masks { namespace boost { /*! @class BoostMaskAccessor * @brief Using the voxelization based on boost::geometry and generate the mask accessor. * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. * (A structure without self interseciton means all contours of the structure have no self intersection, and - * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) - * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored - * and the mask will be calculated, however, it may cause errors in the mask results. + * 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 BoostMaskAccessor: public core::MaskAccessorInterface { public: typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::Structure::StructTypePointer StructTypePointer; private: core::GeometricInfo _geoInfo; /*! vector containing list of mask voxels*/ MaskVoxelListPointer _spRelevantVoxelVector; StructTypePointer _spStructure; IDType _maskUID; bool _strict; public: /*! @brief Constructor with a structure pointer and a geometric info pointer * @param aStructurePointer smart pointer of the structure * @param aGeometricInfoPtr smart pointer of the geometricinfo of the dose * @param 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); + BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, + bool strict = true); /*! @brief destructor*/ ~BoostMaskAccessor(); /*! @brief voxelization of the given structures using boost algorithms*/ void updateMask(); /*! @brief get vector containing all relevant voxels that are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(); /*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold); /*!@brief determine how a given voxel on the dose grid is masked * @param aID ID of the voxel in grid. * @param voxel Reference to the voxel. * @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined. * The relevant volume fraction will be set to zero. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const; /*!@brief determine how a given voxel on the dose grid is masked * @param aIndex 3d index of the voxel in grid. * @param voxel Reference to the voxel. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const; /*! @brief give access to GeometricInfo*/ const core::GeometricInfo& getGeometricInfo() const; /* @ brief is true if dose is on a homogeneous grid * @remark Inhomogeneous grids are not supported at the moment, but if they will be supported in the future the interface does not need to change.*/ bool isGridHomogeneous() const { return true; }; IDType getMaskUID() const { return _maskUID; }; }; } } } #endif diff --git a/code/masks/boost/rttbBoostMaskRedesign.cpp b/code/masks/boost/rttbBoostMaskRedesign.cpp index fd71e53..c1f1854 100644 --- a/code/masks/boost/rttbBoostMaskRedesign.cpp +++ b/code/masks/boost/rttbBoostMaskRedesign.cpp @@ -1,511 +1,695 @@ #include #include #include #include #include #include #include #include #include #include "rttbBoostMaskRedesign.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace masks { namespace boostRedesign { BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, BoostMask::StructPointer aStructure, bool strict) : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _voxelInStructure(::boost::make_shared()) { _isUpToDate = false; if (!_geometricInfo) { throw rttb::core::NullPointerException("Error: Geometric info is NULL!"); } else if (!_structure) { throw rttb::core::NullPointerException("Error: Structure is NULL!"); } } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { if (!_isUpToDate) { calcMask(); } return _voxelInStructure; } void BoostMask::calcMask() { preprocessing(); + std::cout << "preprocessing done" << std::endl; voxelization(); - - if (calcVoxelizationThickness(_voxelizationThickness)) - { - - } - + std::cout << "voxelization done" << std::endl; + generateMaskVoxelList(); + std::cout << "generateMaskVoxelList done" << std::endl; _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::lowest(), std::numeric_limits::lowest(), std::numeric_limits::lowest()); rttb::DoubleVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(), _geometricInfo->getNumRows(), 0); for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) { PolygonType rttbPolygon = *it; PolygonType geometryCoordinatePolygon; //1. convert polygon to geometry coordinate polygons //2. calculate global min/max //3. check if polygon is planar if (!preprocessingPolygon(rttbPolygon, geometryCoordinatePolygon, globalMinGridIndex, globalMaxGridIndex, errorConstant)) { throw rttb::core::Exception("TiltedMaskPlaneException"); } geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon); } rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) + 0.5), GridIndexType(globalMinGridIndex(1) + 0.5), GridIndexType(globalMinGridIndex(2) + 0.5)); rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) + 0.5), GridIndexType(globalMaxGridIndex(1) + 0.5), GridIndexType(globalMaxGridIndex(2) + 0.5)); - std::cout << "global min: " << minIndex.toString() << ", globa max: " << maxIndex.toString() << - std::endl; + //std::cout << "global min: " << minIndex.toString() << ", globa max: " << maxIndex.toString() << std::endl; _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) BoostRingMap ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); //check donut and convert to a map of z index and a vector of boost polygon 2d (with or without holes) _geometryCoordinateBoostPolygonMap.clear(); BoostRingMap::iterator itMap; for (itMap = ringMap.begin(); itMap != ringMap.end(); ++itMap) { BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); _geometryCoordinateBoostPolygonMap.insert(std::pair((*itMap).first, polygonVector)); } } void BoostMask::voxelization() { BoostPolygonMap::iterator it; if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; for (it = _geometryCoordinateBoostPolygonMap.begin(); it != _geometryCoordinateBoostPolygonMap.end(); ++it) { 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); if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } maskArray[x][y] = volumeFraction; - std::cout << "(" << x << "," << y << "): " << volumeFraction << std::endl; + //std::cout << "(" << x << "," << y << "): " << volumeFraction << std::endl; } } //insert into voxelization map _voxelizationMap.insert(std::pair((*it).first, maskArray)); } + } + + void BoostMask::generateMaskVoxelList() + { + if (_globalBoundingBox.size() < 2) + { + throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); + } + + //check homogeneus of the voxelization plane (the contours plane) + if (!calcVoxelizationThickness(_voxelizationThickness)) + { + throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneus!"); + } + + rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); + rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); + int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; + int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; + + + for (unsigned int indexZ = 0; indexZ < _geometricInfo->getNumSlices(); ++indexZ) + { + //calculate weight vector + std::map weightVectorForZ; + calcWeightVector(indexZ, weightVectorForZ); + /*std::cout << "weight vector: " << std::endl; + std::map::iterator itTest; + + for (itTest = weightVectorForZ.begin(); itTest != weightVectorForZ.end(); ++itTest) + { + std::cout << "(" << (*itTest).first << ", " << (*itTest).second << ") "; + } + + std::cout << std::endl;*/ + + //For each x,y, calc sum of all voxelization plane, use weight vector + for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) + { + for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) + { + rttb::VoxelGridIndex3D currentIndex; + currentIndex[0] = x + minIndex[0]; + currentIndex[1] = y + minIndex[1]; + currentIndex[2] = indexZ; + rttb::VoxelGridID gridID; + _geometricInfo->convert(currentIndex, gridID); + double volumeFraction = 0; + + std::map::iterator it; + + for (it = weightVectorForZ.begin(); it != weightVectorForZ.end(); ++it) + { + BoostArrayMap::iterator findVoxelizationIt = _voxelizationMap.find((*it).first); + double weight = (*it).second; + + if (findVoxelizationIt == _voxelizationMap.end()) + { + throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneus!"); + } + + BoostArray2D voxelizationArray = (*findVoxelizationIt).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 maskVoxel = core::MaskVoxel(gridID, volumeFraction); + _voxelInStructure->push_back(maskVoxel);//push back the mask voxel in structure + } + } + + } + + } + + + + } bool BoostMask::preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const { double minZ = _geometricInfo->getNumSlices(); double maxZ = 0.0; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D worldCoordinatePoint = aRTTBPolygon.at(i); //convert to geometry coordinate polygon rttb::DoubleVoxelGridIndex3D geometryCoordinatePoint; _geometricInfo->worldCoordinateToGeometryCoordinate(worldCoordinatePoint, geometryCoordinatePoint); geometryCoordinatePolygon.push_back(geometryCoordinatePoint); //calculate the current global min/max //min and max for x if (geometryCoordinatePoint(0) < minimum(0)) { minimum(0) = geometryCoordinatePoint(0); } if (geometryCoordinatePoint(0) > maximum(0)) { maximum(0) = geometryCoordinatePoint(0); } //min and max for y if (geometryCoordinatePoint(1) < minimum(1)) { minimum(1) = geometryCoordinatePoint(1); } if (geometryCoordinatePoint(1) > maximum(1)) { maximum(1) = geometryCoordinatePoint(1); } //min and max for z if (geometryCoordinatePoint(2) < minimum(2)) { minimum(2) = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maximum(2)) { maximum(2) = geometryCoordinatePoint(2); } //check planar if (geometryCoordinatePoint(2) < minZ) { minZ = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maxZ) { maxZ = geometryCoordinatePoint(2); } } return (abs(maxZ - minZ) <= aErrorConstant); } BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const { BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); if (i == 0) { firstPoint = boostPoint; } ::boost::geometry::append(polygon2D, boostPoint); } ::boost::geometry::append(polygon2D, firstPoint); return polygon2D; } BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( - const rttb::PolygonSequenceType& aRTTBPolygonVector) const + const rttb::PolygonSequenceType& aRTTBPolygonVector) { rttb::PolygonSequenceType::const_iterator it; BoostMask::BoostRingMap aRingMap; for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it) { rttb::PolygonType rttbPolygon = *it; double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon + bool isFirstZ = true; if (!aRingMap.empty()) { - BoostMask::BoostRingMap::const_iterator findIt = findNearestKey(aRingMap, zIndex, errorConstant); + BoostMask::BoostRingMap::iterator findIt = findNearestKey(aRingMap, zIndex, errorConstant); //if the z index is found (same slice), add the polygon to vector if (findIt != aRingMap.end()) { - BoostRingVector ringVector = (*findIt).second; - ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); + //BoostRingVector ringVector = ; + (*findIt).second.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); + isFirstZ = false; } } //if it is the first z index in the map, insert vector with the polygon - BoostRingVector ringVector; - ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); - aRingMap.insert(std::pair(zIndex, ringVector)); + if (isFirstZ) + { + BoostRingVector ringVector; + ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); + aRingMap.insert(std::pair(zIndex, ringVector)); + } } return aRingMap; } - BoostMask::BoostRingMap::const_iterator BoostMask::findNearestKey(const BoostMask::BoostRingMap& - aBoostRingMap, double aIndex, double aErrorConstant) const + BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap& + aBoostRingMap, double aIndex, double aErrorConstant) { - BoostMask::BoostRingMap::const_iterator find = aBoostRingMap.find(aIndex); + BoostMask::BoostRingMap::iterator find = aBoostRingMap.find(aIndex); //if find a key equivalent to aIndex, found if (find != aBoostRingMap.end()) { return find; } else { - BoostMask::BoostRingMap::const_iterator lowerBound = aBoostRingMap.lower_bound(aIndex); + BoostMask::BoostRingMap::iterator lowerBound = aBoostRingMap.lower_bound(aIndex); + + //if all keys go before aIndex, check the last key + if (lowerBound == aBoostRingMap.end()) + { + lowerBound = --aBoostRingMap.end(); + } //if the lower bound very close to aIndex, found if (abs((*lowerBound).first - aIndex) <= aErrorConstant) { return lowerBound; } else { //if the lower bound is the beginning, not found if (lowerBound == aBoostRingMap.begin()) { return aBoostRingMap.end(); } else { - BoostMask::BoostRingMap::const_iterator lowerBound1 = --lowerBound;//the key before the lower bound + BoostMask::BoostRingMap::iterator lowerBound1 = --lowerBound;//the key before the lower bound //if the key before the lower bound very close to a Index, found if (abs((*lowerBound1).first - aIndex) <= aErrorConstant) { return lowerBound1; } //else, not found else { return aBoostRingMap.end(); } } } } } BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector) const { //check donut BoostMask::BoostRingVector::const_iterator it1; BoostMask::BoostRingVector::const_iterator it2; BoostMask::BoostPolygonVector boostPolygonVector; std::vector donutIndexVector;//store the outer and inner ring index BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon //Get donut index and donut polygon unsigned int index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } //if not jet, check now if (!it1IsDonut) { bool it2IsDonut = false; unsigned int index2 = 0; for (it2 = aRingVector.begin(); it2 != aRingVector.end(); it2++, index2++) { if (it2 != it1) { BoostMask::BoostPolygon2D polygon2D; if (::boost::geometry::within(*it1, *it2)) { ::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring it2IsDonut = true; } //if donut else if (::boost::geometry::within(*it2, *it1)) { ::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring it2IsDonut = true; } if (it2IsDonut) { donutIndexVector.push_back(index1); donutIndexVector.push_back(index2); donutVector.push_back(polygon2D);//store donut polygon break;//Only store the first donut! } } } } } //Store no donut polygon to boostPolygonVector index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) { bool it1IsDonut = false; //check if the ring is the outer or inner of a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } if (!it1IsDonut) { BoostMask::BoostPolygon2D polygon2D; ::boost::geometry::append(polygon2D, *it1); boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut } } //Append donut polygon to boostPolygonVector BoostMask::BoostPolygonVector::iterator itDonut; for (itDonut = donutVector.begin(); itDonut != donutVector.end(); itDonut++) { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } BoostMask::BoostRing2D BoostMask::get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) const { BoostMask::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; } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) const { BoostMask::BoostPolygonDeque polygonDeque; BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); ::boost::geometry::correct(voxelPolygon); BoostPolygonVector::const_iterator it; for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) { BoostPolygon2D contour = *it; ::boost::geometry::correct(contour); BoostPolygonDeque intersection; ::boost::geometry::intersection(voxelPolygon, contour, intersection); polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end()); } return polygonDeque; } /*Calculate the intersection area*/ double BoostMask::calcArea(const BoostPolygonDeque& aPolygonDeque) const { double area = 0; BoostPolygonDeque::const_iterator it; for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) { area += ::boost::geometry::area(*it); } return area; } + bool BoostMask::calcVoxelizationThickness(double& aThickness) const + { + BoostArrayMap::const_iterator it = _voxelizationMap.begin(); + BoostArrayMap::const_iterator it2 = ++_voxelizationMap.begin(); + double lastContourIndexZ; + + if (_voxelizationMap.size() <= 1) + { + aThickness = 1; + return true; + } + + double thickness = 0; + + for (; + it != _voxelizationMap.end(), it2 != _voxelizationMap.end(); ++it, ++it2) + { + if (thickness == 0) + { + thickness = (*it2).first - (*it).first; + } + else + { + //if no homogeneous, return false + if (thickness != ((*it2).first - (*it).first)) + { + return false; + } + } + + } + + if (thickness != 0) + { + aThickness = thickness; + } + else + { + aThickness = 1; + } + + return true; + } + + void BoostMask::calcWeightVector(const rttb::VoxelGridID& aIndexZ, + std::map& weightVector) const + { + BoostArrayMap::const_iterator it; + + double indexZMin = aIndexZ - 0.5; + double indexZMax = aIndexZ + 0.5; + + for (it = _voxelizationMap.begin(); it != _voxelizationMap.end(); ++it) + { + double voxelizationPlaneIndexMin = (*it).first - 0.5 * _voxelizationThickness; + double voxelizationPlaneIndexMax = (*it).first + 0.5 * _voxelizationThickness; + double weight = 0; + + if ((voxelizationPlaneIndexMin < indexZMin) && (voxelizationPlaneIndexMax > indexZMin)) + { + if (voxelizationPlaneIndexMax < indexZMax) + { + weight = voxelizationPlaneIndexMax - indexZMin; + } + else + { + weight = 1; + } + } + else if ((voxelizationPlaneIndexMin >= indexZMin) && (voxelizationPlaneIndexMin < indexZMax)) + { + if (voxelizationPlaneIndexMax < indexZMax) + { + weight = _voxelizationThickness; + } + else + { + weight = indexZMax - voxelizationPlaneIndexMin; + } + } + + weightVector.insert(std::pair((*it).first, weight)); + } + } + } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskRedesign.h b/code/masks/boost/rttbBoostMaskRedesign.h index 07004b4..633451c 100644 --- a/code/masks/boost/rttbBoostMaskRedesign.h +++ b/code/masks/boost/rttbBoostMaskRedesign.h @@ -1,210 +1,224 @@ // ----------------------------------------------------------------------- // 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_H -#define __BOOST_MASK_H +#ifndef __BOOST_MASK_R_H +#define __BOOST_MASK_R_H #include "rttbBaseType.h" #include "rttbStructure.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include #include #include #include #include #include namespace rttb { namespace masks { namespace boostRedesign { /*! @class BoostMask * @brief Implementation of voxelization using boost::geometry. * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. * (A structure without self interseciton means all contours of the structure have no self intersection, and * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored * and the mask will be calculated, however, it may cause errors in the mask results. */ class BoostMask { public: typedef ::boost::shared_ptr GeometricInfoPointer; typedef core::Structure::StructTypePointer StructPointer; typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; /*! @brief Constructor * @exception rttb::core::NullPointerException thrown if aDoseGeoInfo or aStructure is NULL * @param strict indicates whether to allow self intersection in the structure. If it is set to true, an exception will be thrown when the given structure has self intersection. * @exception InvalidParameterException thrown if strict is true and the structure has self intersections */ BoostMask(GeometricInfoPointer aDoseGeoInfo, StructPointer aStructure, bool strict = true); /*! @brief Generate mask and return the voxels in the mask * @exception rttb::core::InvalidParameterException thrown if the structure has self intersections */ MaskVoxelListPointer getRelevantVoxelVector(); private: typedef ::boost::geometry::model::d2::point_xy BoostPoint2D; typedef ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy > BoostPolygon2D; typedef ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy > BoostRing2D; typedef std::deque BoostPolygonDeque; typedef std::vector BoostRingVector;//polygon without holes typedef std::vector BoostPolygonVector;//polygon with or without holes typedef std::vector VoxelIndexVector; typedef std::map BoostPolygonMap;//map of the z index with the vector of boost 2d polygon typedef std::map BoostRingMap;//map of the z index with the vector of boost 2d ring - typedef boost::multi_array BoostArray2D; + typedef ::boost::multi_array BoostArray2D; typedef std::map BoostArrayMap; - typedef boost::array BoostArray1D; + typedef ::boost::array BoostArray1D; + typedef std::vector> WeightVector; GeometricInfoPointer _geometricInfo; StructPointer _structure; /*! @brief The map of z index and a vector of all boost polygons with the same z index. This will be used to calculate the mask. * Key: the double z grid index * Value: the vector of all boost 2d polygons with the same z grid index (donut polygon is accepted). */ BoostPolygonMap _geometryCoordinateBoostPolygonMap; /*! @brief The min and max index of the global bounding box. * The first index has the minimum for x/y/z of the global bounding box. * The second index has the maximum for x/y/z of the global bounding index. */ VoxelIndexVector _globalBoundingBox; /*! @brief The voxelization map - * key: the double z grid index + * key: the converted double z grid index of a contour plane * value: the 2d mask, array[i][j] = the mask value of the position (i,j) in the global bounding box, * i: 0 - (_globalBoundingBoxSize0-1), j: 0 - (_globalBoundingBoxSize1-1) */ BoostArrayMap _voxelizationMap; - //The thickness of the voxelization plane (the contour plane) + //The thickness of the voxelization plane (the contour plane), in double dose grid index + //(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2) double _voxelizationThickness; bool _strict; //vector of the MaskVoxel inside the structure MaskVoxelListPointer _voxelInStructure; /*! @brief If the mask is up to date */ bool _isUpToDate; /*! @brief Voxelization and generate mask */ void calcMask(); /*! @brief The preprocessing step, wich consists of the following logic and Sub setps: * For all contours in a struct: * 1) Transfer the contour polygons into boost::geometry structures * 1a) Convert the contur points from world coordinates into geometry coordinates. * 1b) get min and max for x/y/z of a contour * 2) Tilt check: if difference of z_min and z_max is larger then a tolerance value -> there is a tilt. Throw rttb::TiltedMaskPlaneException. * 3) Get struct-bounding-box: get x_min_struct, y_min_struct, x_max_struct, y_max_struct to define the bounding box that containes all contours of a struct in x-y-dimensions. */ void preprocessing(); /*! @brief The voxelization step, wich computes the voxelization planes (in x/y) for all contours of an struct. - *For each contour (that is in the z-Range of the reference geometry) of the struct: - *1) Allocate result array (voxelization plane) based on the bounding box (see Preprocessing Step 3) - *2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry). - *3) Add result Array (key is the z-Value of the contour) + * For each contour (that is in the z-Range of the reference geometry) of the struct: + * 1) Allocate result array (voxelization plane) based on the bounding box (see Preprocessing Step 3) + * 2) Generate voxelization plane for the contour (based on the x-y-raster of the reference geometry). + * 3) Add result Array (key is the z-Value of the contour) */ void voxelization(); + /*! @final mask voxel Generation step which transfers the voxilization planes into the (z-)geometry of the reference geometry. + * It consists of following Sub steps : + * For all "slices" in the reference geometry : + * 1) generate weight vector for all voxelization planes for a given z - value of a slice + * Itterate over the bounding box of a struct.For each voxel : + * 2) Compute weighted sum of all voxelization planes(use weight vector, step 1) + * 2a) If sum > 0 : Add mask voxel for the current x / y(inner Loop) and z value(outer Loop). + * 3) return mask voxel list. + */ + void generateMaskVoxelList(); + /*! @brief Convert the rttb polygon with world corrdinate to the rttb polygon with double geometry coordinate, calculate the current min/max * and check if the polygon is planar * @param minimum the current global minimum * @param maximum the current global maximum * @return Return true if the polygon is planar, which means that the minminal and maximal z-coordinate of the polygon is not larger than a error constant */ bool preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, rttb::PolygonType& geometryCoordinatePolygon, rttb::DoubleVoxelGridIndex3D& minimum, rttb::DoubleVoxelGridIndex3D& maximum, double aErrorConstant) const; /*! @brief Convert a rttb 3d polygon to a 2d boost ring*/ BoostRing2D convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const; /*! @brief Convert a rttb 3d polygon to a map of z index with a vector of boost 2d ring, because of tilt check use the first z index of the polygon as the map key*/ BoostRingMap convertRTTBPolygonSequenceToBoostRingMap(const rttb::PolygonSequenceType& - aRTTBPolygonVector) const; + aRTTBPolygonVector); /*! @brief Find the key with error constant to aIndex * @pre aBoostRingMap should not be empty * @return Return aBoostRingMap.end() if the key is not found */ - BoostMask::BoostRingMap::const_iterator findNearestKey(const BoostMask::BoostRingMap& aBoostRingMap, - double aIndex, double aErrorConstant) const; + BoostMask::BoostRingMap::iterator findNearestKey(BoostMask::BoostRingMap& aBoostRingMap, + double aIndex, double aErrorConstant); /*! @brief If 2 rings in the vector build a donut, convert the 2 rings to a donut polygon, other rings unchanged*/ BoostPolygonVector checkDonutAndConvert(const BoostRingVector& aRingVector) const; /*! @brief Get the voxel 2d contour polygon in geometry coordinate*/ BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) const; /*! @brief Get intersection polygons of the contour and a voxel polygon * @param aVoxelIndex3D The 3d grid index of the voxel * @param intersectionSlicePolygons The polygons of the slice intersecting the voxel * @return Return all intersetion polygons of the structure and the voxel */ BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) const; /*! @brief Calculate the area of all polygons * @param aPolygonDeque The deque of polygons * @return Return the area of all polygons */ double calcArea(const BoostPolygonDeque& aPolygonDeque) const; /*! @brief Calculate the voxelization thickness. Return false, if the voxelization plane is not homogeneous */ - bool calcVoxelizationThickness(double aThickness); + bool calcVoxelizationThickness(double& aThickness) const; /*! @brief For each dose grid index z, calculate the weight vector for each structure contour */ - void calcWeightVector(const rttb::VoxelGridIndex3D& aVoxelIndex3D, BoostArray1D& weightArray) const; + void calcWeightVector(const rttb::VoxelGridID& aIndexZ, + std::map& weightVector) const; }; } } } #endif \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskAccessor.cpp b/code/masks/boost/rttbBoostMaskRedesignAccessor.cpp similarity index 84% copy from code/masks/boost/rttbBoostMaskAccessor.cpp copy to code/masks/boost/rttbBoostMaskRedesignAccessor.cpp index c02b583..bf8fd80 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.cpp +++ b/code/masks/boost/rttbBoostMaskRedesignAccessor.cpp @@ -1,161 +1,164 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) +// @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 "rttbBoostMaskAccessor.h" -#include "rttbBoostMask.h" +#include "rttbBoostMaskRedesignAccessor.h" +#include "rttbBoostMaskRedesign.h" #include #include #include #include namespace rttb { namespace masks { - namespace boost + namespace boostRedesign { - BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, bool strict) + BoostMaskAccessor::BoostMaskAccessor(StructTypePointer aStructurePointer, + const core::GeometricInfo& aGeometricInfo, bool strict) : _spStructure(aStructurePointer), _geoInfo(aGeometricInfo), _strict(strict) { _spRelevantVoxelVector = MaskVoxelListPointer(); //generate new structure set uid ::boost::uuids::uuid id; ::boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _maskUID = "BoostMask_" + ss.str(); } BoostMaskAccessor::~BoostMaskAccessor() { }; void BoostMaskAccessor::updateMask() { MaskVoxelList newRelevantVoxelVector; if (_spRelevantVoxelVector) { return; // already calculated } - BoostMask mask(::boost::make_shared(_geoInfo), _spStructure, _strict); + BoostMask mask(::boost::make_shared(_geoInfo), + _spStructure, _strict); _spRelevantVoxelVector = mask.getRelevantVoxelVector(); return; } BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector() { // if not already generated start voxelization here updateMask(); return _spRelevantVoxelVector; } - BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector(float lowerThreshold) + BoostMaskAccessor::MaskVoxelListPointer BoostMaskAccessor::getRelevantVoxelVector( + float lowerThreshold) { MaskVoxelListPointer filteredVoxelVectorPointer(new MaskVoxelList); updateMask(); // filter relevant voxels BoostMaskAccessor::MaskVoxelList::iterator it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getRelevantVolumeFraction() > lowerThreshold) { filteredVoxelVectorPointer->push_back(*it); } ++it; } // if mask calculation was not successful this is empty! return filteredVoxelVectorPointer; } bool BoostMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const { //initialize return voxel voxel.setRelevantVolumeFraction(0); //check if ID is valid if (!_geoInfo.validID(aID)) { return false; } //determine how a given voxel on the dose grid is masked if (_spRelevantVoxelVector) { BoostMaskAccessor::MaskVoxelList::iterator it = _spRelevantVoxelVector->begin(); while (it != _spRelevantVoxelVector->end()) { if ((*it).getVoxelGridID() == aID) { voxel = (*it); return true; } ++it; } } // returns false if mask was not calculated without triggering calculation (otherwise not const!) else { return false; } return false; } bool BoostMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const { //convert VoxelGridIndex3D to VoxelGridID VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getMaskAt(aVoxelGridID, voxel); } else { return false; } } const core::GeometricInfo& BoostMaskAccessor::getGeometricInfo() const { return _geoInfo; }; } } } \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMaskAccessor.h b/code/masks/boost/rttbBoostMaskRedesignAccessor.h similarity index 89% copy from code/masks/boost/rttbBoostMaskAccessor.h copy to code/masks/boost/rttbBoostMaskRedesignAccessor.h index 0fdf506..850c36a 100644 --- a/code/masks/boost/rttbBoostMaskAccessor.h +++ b/code/masks/boost/rttbBoostMaskRedesignAccessor.h @@ -1,127 +1,128 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author $Author$ (last changed by) +// @version $Revision: 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_ACCESSOR__H -#define __BOOST_MASK_ACCESSOR__H +#ifndef __BOOST_MASK_R_ACCESSOR__H +#define __BOOST_MASK_R_ACCESSOR__H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskVoxel.h" #include "rttbMaskAccessorInterface.h" #include "rttbGenericDoseIterator.h" #include "rttbStructure.h" #include namespace rttb { namespace masks { - namespace boost + namespace boostRedesign { /*! @class BoostMaskAccessor * @brief Using the voxelization based on boost::geometry and generate the mask accessor. * @attention If "strict" is set to true, an exception will be thrown when the given structure has self intersection. * (A structure without self interseciton means all contours of the structure have no self intersection, and - * the polygons on the same slice have no intersection between each other, unless the case of a donut. A donut is accepted.) - * If "strict" is set to false, debug information will be displayed when the given structure has self intersection. Self intersections will be ignored - * and the mask will be calculated, however, it may cause errors in the mask results. + * 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 BoostMaskAccessor: public core::MaskAccessorInterface { public: typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::Structure::StructTypePointer StructTypePointer; private: core::GeometricInfo _geoInfo; /*! vector containing list of mask voxels*/ MaskVoxelListPointer _spRelevantVoxelVector; StructTypePointer _spStructure; IDType _maskUID; bool _strict; public: /*! @brief Constructor with a structure pointer and a geometric info pointer * @param aStructurePointer smart pointer of the structure * @param aGeometricInfoPtr smart pointer of the geometricinfo of the dose * @param 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); + BoostMaskAccessor(StructTypePointer aStructurePointer, const core::GeometricInfo& aGeometricInfo, + bool strict = true); /*! @brief destructor*/ ~BoostMaskAccessor(); /*! @brief voxelization of the given structures using boost algorithms*/ void updateMask(); /*! @brief get vector containing all relevant voxels that are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(); /*! @brief get vector containing all relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold); /*!@brief determine how a given voxel on the dose grid is masked * @param aID ID of the voxel in grid. * @param voxel Reference to the voxel. * @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined. * The relevant volume fraction will be set to zero. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const; /*!@brief determine how a given voxel on the dose grid is masked * @param aIndex 3d index of the voxel in grid. * @param voxel Reference to the voxel. * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const; /*! @brief give access to GeometricInfo*/ const core::GeometricInfo& getGeometricInfo() const; /* @ brief is true if dose is on a homogeneous grid * @remark Inhomogeneous grids are not supported at the moment, but if they will be supported in the future the interface does not need to change.*/ bool isGridHomogeneous() const { return true; }; IDType getMaskUID() const { return _maskUID; }; }; } } } #endif diff --git a/testing/examples/CMakeLists.txt b/testing/examples/CMakeLists.txt index d36fa75..3249a16 100644 --- a/testing/examples/CMakeLists.txt +++ b/testing/examples/CMakeLists.txt @@ -1,42 +1,42 @@ #----------------------------------------------------------------------------- # Setup the system information test. Write out some basic failsafe # information in case the test doesn't run. #----------------------------------------------------------------------------- SET(CORE_TEST_EXAMPLES ${EXECUTABLE_OUTPUT_PATH}/rttbExamplesTests) SET(TEST_DATA_ROOT ${RTTBTesting_SOURCE_DIR}/data) SET(TEMP ${RTTBTesting_BINARY_DIR}/temporary) #----------------------------------------------------------------------------- ADD_TEST(RTBioModelExampleTest ${CORE_TEST_EXAMPLES} RTBioModelExampleTest "${TEST_DATA_ROOT}/TestDVH/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT3.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_virtuos_diff_trunk6.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_virtuos_diff_trunk8.txt") ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwoGridScaling.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwoGridScaling05.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantFiftyGridScaling01.dcm") ADD_TEST(RTDoseIndexTest ${CORE_TEST_EXAMPLES} RTDoseIndexTest "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT2.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT3.txt") ADD_TEST(RTDoseStatisticsTest ${CORE_TEST_EXAMPLES} RTDoseStatisticsTest "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo_withDoseGridScaling.dcm" "${TEST_DATA_ROOT}/Virtuos/prostate_ac/prostate_ac101.dos.gz" "${TEST_DATA_ROOT}/Virtuos/prostate_ac/prostate_ac000.vdx" "${TEST_DATA_ROOT}/Virtuos/prostate_ac/prostate_ac101.pln") ADD_TEST(RTDVHTest ${CORE_TEST_EXAMPLES} RTDVHTest "${TEST_DATA_ROOT}/TestDVH/dvh_test.txt") ADD_TEST(DVHCalculatorExampleTest ${CORE_TEST_EXAMPLES} DVHCalculatorExampleTest "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/ConstantTwo_withDoseGridScaling.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncreaseX.dcm") ADD_TEST(RTBioModelScatterPlotExampleTest ${CORE_TEST_EXAMPLES} RTBioModelScatterPlotExampleTest "${TEST_DATA_ROOT}/TestDVH/dvh_PTV_HIT.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_HT1.txt" "${TEST_DATA_ROOT}/TestDVH/dvh_test_TV.txt") ADD_TEST(VoxelizationValidationTest ${CORE_TEST_EXAMPLES} VoxelizationValidationTest "${TEST_DATA_ROOT}/DICOM/StructureSet/RS1.3.6.1.4.1.2452.6.841242143.1311652612.1170940299.4217870819.dcm" -"${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/boostMask/" "${TEST_DATA_ROOT}/OTBMask/") +"${TEST_DATA_ROOT}/DICOM/TestDose/LinearIncrease3D.dcm" "${TEST_DATA_ROOT}/BoostMask/" "${TEST_DATA_ROOT}/OTBMask/" "${TEST_DATA_ROOT}/BoostMaskRedesign/") RTTB_CREATE_TEST_MODULE(rttbExamples DEPENDS RTTBCore RTTBAlgorithms RTTBMasks RTTBOTBMask RTTBBoostMask RTTBIndices RTTBDicomIO RTTBVirtuosIO RTTBITKIO RTTBOtherIO RTTBModels PACKAGE_DEPENDS Litmus) diff --git a/testing/examples/VoxelizationValidationTest.cpp b/testing/examples/VoxelizationValidationTest.cpp index 95cdd92..fa53eca 100644 --- a/testing/examples/VoxelizationValidationTest.cpp +++ b/testing/examples/VoxelizationValidationTest.cpp @@ -1,199 +1,250 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVHCalculator.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileStructureSetGenerator.h" #include "rttbOTBMaskAccessor.h" #include "rttbDVHTxtFileReader.h" #include "rttbBoostMaskAccessor.h" #include "rttbITKImageMaskAccessorConverter.h" #include "rttbImageWriter.h" +#include "rttbBoostMaskRedesign.h" +#include "rttbBoostMaskRedesignAccessor.h" namespace rttb { namespace testing { /*! @brief VoxelizationValidationTest. Compare two differnt voxelizations: OTB and Boost. Check dvh maximum and minimum for each structure. Check write mask to itk file for further validation. */ int VoxelizationValidationTest(int argc, char* argv[]) { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: structure file name // 2: dose1 file name // 3: directory name to write boost mask of all structures // 4: directory name to write OTB mask of all structures std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string BoostMask_DIRNAME; std::string OTBMask_DIRNAME; - + std::string BoostMaskRedesign_DIRNAME; if (argc > 4) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; BoostMask_DIRNAME = argv[3]; OTBMask_DIRNAME = argv[4]; + BoostMaskRedesign_DIRNAME = argv[5]; } OFCondition status; DcmFileFormat fileformat; /* read dicom-rt dose */ io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator1(RTDOSE_FILENAME.c_str()); DoseAccessorPointer doseAccessor1(doseAccessorGenerator1.generateDoseAccessor()); //create a vector of MaskAccessors (one for each structure) StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); std::vector rtStructSetMaskAccessorVec; if (rtStructureSet->getNumberOfStructures() > 0) { for (int j = 0; j < rtStructureSet->getNumberOfStructures(); j++) { - std::cout << j << ": " << rtStructureSet->getStructure(j)->getLabel() << std::endl; - clock_t start(clock()); - //create OTB MaskAccessor - ::boost::shared_ptr spOTBMaskAccessor = - ::boost::make_shared(rtStructureSet->getStructure(j), - doseAccessor1->getGeometricInfo()); - spOTBMaskAccessor->updateMask(); - MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor); - - ::boost::shared_ptr spMaskedDoseIteratorTmp = - ::boost::make_shared(spMaskAccessor, doseAccessor1); - DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); - rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), - doseAccessor1->getUID()); - rttb::core::DVH dvh = *(calc.generateDVH()); - - clock_t finish(clock()); - std::cout << "OTB Mask Calculation time: " << finish - start << " ms" << std::endl; - - //Write the mask image to a file. - /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. + if (j != 2 && j != 3) + { + std::cout << j << ": " << rtStructureSet->getStructure(j)->getLabel() << std::endl; + clock_t start(clock()); + //create OTB MaskAccessor + ::boost::shared_ptr spOTBMaskAccessor = + ::boost::make_shared(rtStructureSet->getStructure(j), + doseAccessor1->getGeometricInfo()); + spOTBMaskAccessor->updateMask(); + MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor); + + ::boost::shared_ptr spMaskedDoseIteratorTmp = + ::boost::make_shared(spMaskAccessor, doseAccessor1); + DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); + rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), + doseAccessor1->getUID()); + rttb::core::DVH dvh = *(calc.generateDVH()); + + clock_t finish(clock()); + std::cout << "OTB Mask Calculation time: " << finish - start << " ms" << std::endl; + + //Write the mask image to a file. + /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. + To write all mask files, please use the outcommented code and call the .exe directly! + */ + /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter(spOTBMaskAccessor); + CHECK(itkConverter.process()); + std::stringstream fileNameSstr; + fileNameSstr << OTBMask_DIRNAME << j << ".mhd"; + rttb::io::itk::ImageWriter writer(fileNameSstr.str(), itkConverter.getITKImage()); + CHECK(writer.writeFile());*/ + + + + clock_t start2(clock()); + //create Boost MaskAccessor + MaskAccessorPointer boostMaskAccessorPtr + = ::boost::make_shared + (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); + + CHECK_NO_THROW(boostMaskAccessorPtr->updateMask()); + + ::boost::shared_ptr spMaskedDoseIteratorTmp2 = + ::boost::make_shared(boostMaskAccessorPtr, doseAccessor1); + DoseIteratorPointer spMaskedDoseIterator2(spMaskedDoseIteratorTmp2); + rttb::core::DVHCalculator calc2(spMaskedDoseIterator2, (rtStructureSet->getStructure(j))->getUID(), + doseAccessor1->getUID()); + rttb::core::DVH dvh2 = *(calc2.generateDVH()); + + clock_t finish2(clock()); + std::cout << "Boost Mask Calculation and write file time: " << finish2 - start2 << " ms" << + std::endl; + + //Write the mask image to a file. + /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. + To write all mask files, please use the outcommented code and call the .exe directly! + */ + + /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter2(boostMaskAccessorPtr); + CHECK(itkConverter2.process()); + std::stringstream fileNameSstr2; + fileNameSstr2 << BoostMask_DIRNAME << j << ".mhd"; + rttb::io::itk::ImageWriter writer2(fileNameSstr2.str(), itkConverter2.getITKImage()); + CHECK(writer2.writeFile());*/ + + //create Boost MaskAccessor redesign + clock_t startR(clock()); + + MaskAccessorPointer boostMaskRPtr + = ::boost::make_shared + (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); + CHECK_NO_THROW(boostMaskRPtr->updateMask()); + ::boost::shared_ptr spMaskedDoseIteratorTmpR = + ::boost::make_shared(boostMaskRPtr, doseAccessor1); + DoseIteratorPointer spMaskedDoseIteratorR(spMaskedDoseIteratorTmpR); + rttb::core::DVHCalculator calcR(spMaskedDoseIteratorR, (rtStructureSet->getStructure(j))->getUID(), + doseAccessor1->getUID()); + rttb::core::DVH dvhR = *(calcR.generateDVH()); + clock_t finishR(clock()); + std::cout << "Boost Mask Redesign Calculation and write file time: " << finishR - startR << " ms" << + std::endl; + + //Write the mask image to a file. + /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. To write all mask files, please use the outcommented code and call the .exe directly! - */ - /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter(spOTBMaskAccessor); - CHECK(itkConverter.process()); - std::stringstream fileNameSstr; - fileNameSstr << OTBMask_DIRNAME << j << ".mhd"; - rttb::io::itk::ImageWriter writer(fileNameSstr.str(), itkConverter.getITKImage()); - CHECK(writer.writeFile());*/ - - - - clock_t start2(clock()); - //create Boost MaskAccessor - MaskAccessorPointer boostMaskAccessorPtr = ::boost::make_shared - (rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); + */ + rttb::io::itk::ITKImageMaskAccessorConverter itkConverterR(boostMaskRPtr); + CHECK(itkConverterR.process()); - CHECK_NO_THROW(boostMaskAccessorPtr->updateMask()); - ::boost::shared_ptr spMaskedDoseIteratorTmp2 = - ::boost::make_shared(boostMaskAccessorPtr, doseAccessor1); - DoseIteratorPointer spMaskedDoseIterator2(spMaskedDoseIteratorTmp2); - rttb::core::DVHCalculator calc2(spMaskedDoseIterator2, (rtStructureSet->getStructure(j))->getUID(), - doseAccessor1->getUID()); - rttb::core::DVH dvh2 = *(calc2.generateDVH()); + std::stringstream fileNameSstrR; + fileNameSstrR << BoostMaskRedesign_DIRNAME << j << ".mhd"; + rttb::io::itk::ImageWriter writerR(fileNameSstrR.str(), itkConverterR.getITKImage()); + CHECK(writerR.writeFile()); - clock_t finish2(clock()); - std::cout << "Boost Mask Calculation and write file time: " << finish2 - start2 << " ms" << std::endl; + //check close of 2 voxelizatin: OTB and Boost + CHECK_CLOSE(dvh.getMaximum(), dvh2.getMaximum(), 0.1); + CHECK_CLOSE(dvh.getMinimum(), dvh2.getMinimum(), 0.1); - //Write the mask image to a file. - /*! It takes a long time to write all mask files so that RUN_TESTS causes a timeout error. - To write all mask files, please use the outcommented code and call the .exe directly! - */ + if (j != 7) + { + CHECK_CLOSE(dvh.getMean(), dvh2.getMean(), 0.1); + } - /*rttb::io::itk::ITKImageMaskAccessorConverter itkConverter2(boostMaskAccessorPtr); - CHECK(itkConverter2.process()); - std::stringstream fileNameSstr2; - fileNameSstr2 << BoostMask_DIRNAME << j << ".mhd"; - rttb::io::itk::ImageWriter writer2(fileNameSstr2.str(), itkConverter2.getITKImage()); - CHECK(writer2.writeFile());*/ + CHECK_CLOSE(dvh.getMedian(), dvh2.getMedian(), 0.1); + CHECK_CLOSE(dvh.getModal(), dvh2.getModal(), 0.1); + //0: Aussenkontur and 3: Niere li. failed. + if (j != 0 && j != 3) + { + CHECK_CLOSE(dvh.getVx(0), dvh2.getVx(0), dvh.getVx(0) * 0.05); //check volume difference < 5% + } + //check close of 2 voxelization: Boost and BoostRedesign + CHECK_CLOSE(dvhR.getMaximum(), dvh2.getMaximum(), 0.1); + CHECK_CLOSE(dvhR.getMinimum(), dvh2.getMinimum(), 0.1); - //check close of 2 voxelizatin: OTB and Boost - CHECK_CLOSE(dvh.getMaximum(), dvh2.getMaximum(), 0.1); - CHECK_CLOSE(dvh.getMinimum(), dvh2.getMinimum(), 0.1); + if (j != 7) + { + CHECK_CLOSE(dvhR.getMean(), dvh2.getMean(), 0.1); + } - if (j != 7) //7: Ref.Pkt, mean = -1.#IND - { - CHECK_CLOSE(dvh.getMean(), dvh2.getMean(), 0.1); - } + CHECK_CLOSE(dvhR.getMedian(), dvh2.getMedian(), 0.1); + CHECK_CLOSE(dvhR.getModal(), dvh2.getModal(), 0.1); - CHECK_CLOSE(dvh.getMedian(), dvh2.getMedian(), 0.1); - CHECK_CLOSE(dvh.getModal(), dvh2.getModal(), 0.1); + //0: Aussenkontur and 3: Niere li. failed. + CHECK_CLOSE(dvhR.getVx(0), dvh2.getVx(0), dvhR.getVx(0) * 0.05); //check volume difference < 5% - //0: Aussenkontur and 3: Niere li. failed. - if (j != 0 && j != 3) - { - CHECK_CLOSE(dvh.getVx(0), dvh2.getVx(0), dvh.getVx(0) * 0.05); //check volume difference < 5% } - } } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/masks/boost/BoostMaskRedesignTest.cpp b/testing/masks/boost/BoostMaskRedesignTest.cpp index e174b20..6e94c0c 100644 --- a/testing/masks/boost/BoostMaskRedesignTest.cpp +++ b/testing/masks/boost/BoostMaskRedesignTest.cpp @@ -1,83 +1,98 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 1127 $ (last changed revision) // @date $Date: 2015-10-01 13:33:33 +0200 (Do, 01 Okt 2015) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "../../core/DummyStructure.h" #include "../../core/DummyDoseAccessor.h" +#include "rttbDicomDoseAccessor.h" +#include "rttbDicomFileDoseAccessorGenerator.h" +#include "rttbDicomFileStructureSetGenerator.h" +#include "rttbDicomFileStructureSetGenerator.h" +#include "rttbGenericDoseIterator.h" +#include "rttbDVHCalculator.h" +#include "rttbGenericMaskedDoseIterator.h" #include "rttbBoostMaskRedesign.h" +#include "rttbBoostMask.h" +#include "rttbBoostMaskAccessor.h" +#include "rttbBoostMaskRedesignAccessor.h" namespace rttb { namespace testing { /*! @brief BoostMaskRedesignTest. 1) test constructors 2) test getRelevantVoxelVector 3) test getMaskAt */ int BoostMaskRedesignTest(int argc, char* argv[]) { + typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; + typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; + typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; + typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; + PREPARE_DEFAULT_TEST_REPORTING; typedef core::Structure::StructTypePointer StructTypePointer; // generate test structure set boost::shared_ptr spTestDoseAccessor = boost::make_shared(); DummyStructure myStructGenerator(spTestDoseAccessor->getGeometricInfo()); GridIndexType zPlane = 4; core::Structure myTestStruct = myStructGenerator.CreateRectangularStructureCentered(zPlane); StructTypePointer spMyStruct = boost::make_shared(myTestStruct); boost::shared_ptr geometricPtr = boost::make_shared (spTestDoseAccessor->getGeometricInfo()); //1) test BoostMask constructor CHECK_NO_THROW(rttb::masks::boostRedesign::BoostMask(geometricPtr, spMyStruct)); rttb::masks::boostRedesign::BoostMask boostMask = rttb::masks::boostRedesign::BoostMask( geometricPtr, spMyStruct); //2) test getRelevantVoxelVector CHECK_NO_THROW(boostMask.getRelevantVoxelVector()); rttb::masks::boostRedesign::BoostMask::MaskVoxelListPointer list = boostMask.getRelevantVoxelVector(); for (int i = 0; i < list->size(); ++i) { std::cout << "id: " << list->at(i).getVoxelGridID() << ", " << list->at( i).getRelevantVolumeFraction() << std::endl; } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb