diff --git a/code/interpolation/rttbInterpolationBase.cpp b/code/interpolation/rttbInterpolationBase.cpp index 65f134a..3d95507 100644 --- a/code/interpolation/rttbInterpolationBase.cpp +++ b/code/interpolation/rttbInterpolationBase.cpp @@ -1,191 +1,196 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ +#include + #include "rttbInterpolationBase.h" #include "rttbInvalidParameterException.h" #include "rttbNullPointerException.h" #include "rttbMappingOutsideOfImageException.h" namespace rttb { namespace interpolation { void InterpolationBase::setDoseAccessorPointer(const DoseAccessorPointer originalDose) { if (originalDose != NULL) { _spOriginalDose = originalDose; } else { throw core::NullPointerException("originalDose is NULL!"); } }; void InterpolationBase::getNeighborhoodVoxelValues( const WorldCoordinate3D& aWorldCoordinate, unsigned int neighborhood, boost::array& target, - boost::shared_array values) const + boost::shared_ptr values) const { if (_spOriginalDose == NULL) { throw core::NullPointerException("originalDose is NULL!"); } //Determine target (abs(desired worldCoordinate- corner pixel world coordinate/pixel spacing) and values of corner pixels (from originalDose) VoxelGridIndex3D aIndex; if (_spOriginalDose->getGeometricInfo().worldCoordinateToIndex(aWorldCoordinate, aIndex)) { //determine just the nearest voxel to the world coordinate if (neighborhood == 0) { values[0] = _spOriginalDose->getDoseAt(aIndex); } //determine the 8 voxels around the world coordinate else if (neighborhood == 8) { - std::vector cornerPoints; + std::list cornerPoints; WorldCoordinate3D theNextVoxel; _spOriginalDose->getGeometricInfo().indexToWorldCoordinate(aIndex, theNextVoxel); SpacingVectorType3D pixelSpacing = (_spOriginalDose->getGeometricInfo()).getSpacing(); VoxelGridIndex3D leftTopFrontCoordinate; //find the voxel with the smallest coordinate values in each dimension. This defines the standard cube for (int i = 0; i < 3; i++) { if (aWorldCoordinate[i] < theNextVoxel[i]) { if (aIndex[i] > 0) { leftTopFrontCoordinate[i] = aIndex[i] - 1; target[i] = (aWorldCoordinate[i] - (theNextVoxel[i] - pixelSpacing[i])) / pixelSpacing[i]; } //@todo: this is a workaround, not a good solution else { leftTopFrontCoordinate[i] = aIndex[i]; target[i] = (aWorldCoordinate[i] - (theNextVoxel[i] - pixelSpacing[i])) / pixelSpacing[i]; } } else { leftTopFrontCoordinate[i] = aIndex[i]; target[i] = (aWorldCoordinate[i] - theNextVoxel[i]) / pixelSpacing[i]; } } for (int zIncr = 0; zIncr < 2; zIncr++) { for (int yIncr = 0; yIncr < 2; yIncr++) { for (int xIncr = 0; xIncr < 2; xIncr++) { cornerPoints.push_back(VoxelGridIndex3D(leftTopFrontCoordinate[0] + xIncr, leftTopFrontCoordinate[1] + yIncr, leftTopFrontCoordinate[2] + zIncr)); } } } //target range has to be always [0,1] for (int i = 0; i < 3; i++) { assert(target[i] >= 0.0 && target[i] <= 1.0); } + unsigned int count = 0; + //now just get the values of all (dose) voxels and store them in values - for (int i = 0; i < cornerPoints.size(); ++i) + for (auto cornerPointsIterator = std::begin(cornerPoints); cornerPointsIterator != std::end(cornerPoints); + ++cornerPointsIterator, ++count) { - if (_spOriginalDose->getGeometricInfo().isInside(cornerPoints.at(i))) + if (_spOriginalDose->getGeometricInfo().isInside(*cornerPointsIterator)) { - values[i] = _spOriginalDose->getDoseAt(cornerPoints.at(i)); + values[count] = _spOriginalDose->getDoseAt(*cornerPointsIterator); } else { //outside value! boundary treatment - values[i] = getNearestInsideVoxelValue(cornerPoints.at(i)); + values[count] = getNearestInsideVoxelValue(*cornerPointsIterator); } - assert(values[i] != -1); + assert(values[count] != -1); } } else { throw core::InvalidParameterException("neighborhoods other than 0 and 8 not yet supported in Interpolation"); } } else { throw core::MappingOutsideOfImageException("Error in conversion from world coordinates to index"); } } DoseTypeGy InterpolationBase::getNearestInsideVoxelValue(const VoxelGridIndex3D& currentVoxelIndex) const { VoxelGridIndex3D voxelChangedXYZ[] = {currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex, currentVoxelIndex}; int runningIndex; //x,y,z for (runningIndex = 0; runningIndex < 3; ++runningIndex) { voxelChangedXYZ[runningIndex][runningIndex] -= 1; } //xy voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][1] -= 1; ++runningIndex; //xz voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; //yz voxelChangedXYZ[runningIndex][1] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; //xyz voxelChangedXYZ[runningIndex][0] -= 1; voxelChangedXYZ[runningIndex][1] -= 1; voxelChangedXYZ[runningIndex][2] -= 1; ++runningIndex; int replacementVoxelIndex = 0; while (replacementVoxelIndex < runningIndex) { if (_spOriginalDose->getGeometricInfo().validIndex(voxelChangedXYZ[replacementVoxelIndex])) { return _spOriginalDose->getDoseAt(voxelChangedXYZ[replacementVoxelIndex]); } ++replacementVoxelIndex; } return -1; } }//end namespace core }//end namespace rttb diff --git a/code/interpolation/rttbInterpolationBase.h b/code/interpolation/rttbInterpolationBase.h index b0a28a8..2f0d46f 100644 --- a/code/interpolation/rttbInterpolationBase.h +++ b/code/interpolation/rttbInterpolationBase.h @@ -1,95 +1,92 @@ // ----------------------------------------------------------------------- // 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 __INTERPOLATION_BASE_H #define __INTERPOLATION_BASE_H -#include #include #include -#include #include #include "rttbDoseAccessorInterface.h" #include "rttbBaseType.h" namespace rttb { namespace interpolation { /*! @class InterpolationBase @brief Base class for interpolation. @ingroup interpolation */ class InterpolationBase { public: typedef boost::shared_ptr Pointer; typedef rttb::core::DoseAccessorInterface::DoseAccessorPointer DoseAccessorPointer; /*! @brief Constructor */ InterpolationBase() {}; /*! @brief Virtual destructor of base class */ virtual ~InterpolationBase() {}; /*! @brief Sets the DoseAccessPointer @pre originalDose initialized @exception core::NullPointerException if originalDose==NULL */ void setDoseAccessorPointer(const DoseAccessorPointer originalDose); /*! @brief Returns the interpolated value for the given world coordinate */ virtual DoseTypeGy getValue(const WorldCoordinate3D& aWorldCoordinate) const = 0; protected: DoseAccessorPointer _spOriginalDose; /*! @brief determines voxels in a certain neighborhood of a physical based coordinate and converts in a standard cube with corner points [0 0 0], [1 0 0], [0 1 0], [1 1 0], [0 0 1], [1 0 1], [0 1 1], [1 1 1]. @param aWorldCoordinate the coordinate where to start @param neighborhood voxel around coordinate (currently only 0 and 8 implemented) @param target coordinates inside the standard cube with values [0 1] in each dimension. - @param values dose values at all corner points of the standard cube. Is of type boost:shared_array[neighborhood] - @pre target and values have to be correctly initialized (e.g. boost::array target = {0.0, 0.0, 0.0}; boost::shared_array values(new DoseTypeGy[8]()); where 8 is neighborhood) + @param values dose values at all corner points of the standard cube. Is of type boost:shared_ptr[neighborhood] + @pre target and values have to be correctly initialized (e.g. boost::array target = {0.0, 0.0, 0.0}; boost::shared_ptr values(new DoseTypeGy[8]()); where 8 is neighborhood) @exception core::InvalidParameterException if neighborhood =! 0 && !=8 @exception core::MappingOutsideOfImageException if initial mapping of aWorldCoordinate is outside image @exception core::NullPointerException if dose is NULL - @todo use shared_ptr[] with boost >1.52 instead of shared_array */ void getNeighborhoodVoxelValues(const WorldCoordinate3D& aWorldCoordinate, - unsigned int neighborhood, boost::array& target, - boost::shared_array values) const; + unsigned int neighborhood, boost::array& target, + boost::shared_ptr values) const; /*! @brief returns the nearest inside voxel value @pre the voxelGridIndex is outside the image and voxelGridIndex>image.size() for all dimensions. Also voxelGridIndex[]>=0 for all dimensions @note used for virtually expanding the image by one voxel as edge handling */ DoseTypeGy getNearestInsideVoxelValue(const VoxelGridIndex3D& currentVoxelIndex) const; }; } } #endif diff --git a/code/interpolation/rttbLinearInterpolation.cpp b/code/interpolation/rttbLinearInterpolation.cpp index e54ab53..60a1ae1 100644 --- a/code/interpolation/rttbLinearInterpolation.cpp +++ b/code/interpolation/rttbLinearInterpolation.cpp @@ -1,57 +1,57 @@ // ----------------------------------------------------------------------- // 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 "rttbLinearInterpolation.h" namespace rttb { namespace interpolation { DoseTypeGy LinearInterpolation::trilinear(boost::array target, - boost::shared_array values) const + boost::shared_ptr values) const { //4 linear interpolation in x direction DoseTypeGy c_00 = values[0] * (1.0 - target[0]) + values[1] * target[0]; DoseTypeGy c_10 = values[2] * (1.0 - target[0]) + values[3] * target[0]; DoseTypeGy c_01 = values[4] * (1.0 - target[0]) + values[5] * target[0]; DoseTypeGy c_11 = values[6] * (1.0 - target[0]) + values[7] * target[0]; //combine result in y direction DoseTypeGy c_0 = c_00 * (1.0 - target[1]) + c_10 * target[1]; DoseTypeGy c_1 = c_01 * (1.0 - target[1]) + c_11 * target[1]; //finally incorporate z direction return (c_0 * (1.0 - target[2]) + c_1 * target[2]); } DoseTypeGy LinearInterpolation::getValue(const WorldCoordinate3D& aWorldCoordinate) const { //proper initialization of target and values boost::array target = {0.0, 0.0, 0.0}; - boost::shared_array values(new DoseTypeGy[8]()); + boost::shared_ptr values(new DoseTypeGy[8]()); getNeighborhoodVoxelValues(aWorldCoordinate, 8, target, values); return trilinear(target, values); } } } diff --git a/code/interpolation/rttbLinearInterpolation.h b/code/interpolation/rttbLinearInterpolation.h index 0556017..ea3f3c3 100644 --- a/code/interpolation/rttbLinearInterpolation.h +++ b/code/interpolation/rttbLinearInterpolation.h @@ -1,62 +1,61 @@ // ----------------------------------------------------------------------- // 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 __LINEAR_INTERPOLATION_H #define __LINEAR_INTERPOLATION_H -#include #include #include "rttbInterpolationBase.h" namespace rttb { namespace interpolation { /*! @class LinearInterpolation @brief Linear interpolation. @ingroup interpolation */ class LinearInterpolation : public InterpolationBase { public: typedef rttb::core::DoseAccessorInterface::DoseAccessorPointer DoseAccessorPointer; /*! @brief Constructor */ LinearInterpolation() {}; /*! @brief Returns the interpolated value */ DoseTypeGy getValue(const WorldCoordinate3D& aWorldCoordinate) const; private: /*! @brief Trilinar interpolation @sa InterpolationBase for details about target and values @note Source: http://en.wikipedia.org/wiki/Trilinear_interpolation */ - DoseTypeGy trilinear(boost::array target, boost::shared_array values) const; + DoseTypeGy trilinear(boost::array target, boost::shared_ptr values) const; }; } } #endif diff --git a/code/interpolation/rttbNearestNeighborInterpolation.cpp b/code/interpolation/rttbNearestNeighborInterpolation.cpp index e3e0484..69edbf0 100644 --- a/code/interpolation/rttbNearestNeighborInterpolation.cpp +++ b/code/interpolation/rttbNearestNeighborInterpolation.cpp @@ -1,38 +1,40 @@ // ----------------------------------------------------------------------- // 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 "rttbNearestNeighborInterpolation.h" +#include + namespace rttb { namespace interpolation { DoseTypeGy NearestNeighborInterpolation::getValue(const WorldCoordinate3D& aWorldCoordinate) const { //proper initialization of target and values (although target is irrelevant in nearest neighbor case) boost::array target = {0.0, 0.0, 0.0}; - boost::shared_array values(new DoseTypeGy[8]()); + boost::shared_ptr values(new DoseTypeGy[8]()); getNeighborhoodVoxelValues(aWorldCoordinate, 0, target, values); return values[0]; } } } diff --git a/code/interpolation/rttbNearestNeighborInterpolation.h b/code/interpolation/rttbNearestNeighborInterpolation.h index 1baf3fa..ffc8079 100644 --- a/code/interpolation/rttbNearestNeighborInterpolation.h +++ b/code/interpolation/rttbNearestNeighborInterpolation.h @@ -1,54 +1,51 @@ // ----------------------------------------------------------------------- // 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 __NEAREST_NEIGHBOR_INTERPOLATION_H #define __NEAREST_NEIGHBOR_INTERPOLATION_H -#include -#include - #include "rttbInterpolationBase.h" namespace rttb { namespace interpolation { /*! @class NearestNeighborInterpolation @brief Nearest Neighbor interpolation @ingroup interpolation */ class NearestNeighborInterpolation : public InterpolationBase { public: /*! @brief Constructor */ NearestNeighborInterpolation() {}; /*! @brief Returns the interpolated value (the nearest voxel value given by _spOriginalDose->getGeometricInfo().worldCoordinateToIndex()) */ DoseTypeGy getValue(const WorldCoordinate3D& aWorldCoordinate) const; }; } } #endif diff --git a/testing/examples/DVHCalculatorComparisonTest.cpp b/testing/examples/DVHCalculatorComparisonTest.cpp index 4ecf496..e562487 100644 --- a/testing/examples/DVHCalculatorComparisonTest.cpp +++ b/testing/examples/DVHCalculatorComparisonTest.cpp @@ -1,373 +1,373 @@ // ----------------------------------------------------------------------- // 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 "dcmtk/config/osconfig.h" /* make sure OS specific configuration is included first */ #include "dcmtk/dcmrt/drtdose.h" #include "dcmtk/dcmrt/drtstrct.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" namespace rttb { namespace testing { /*! @brief DVHCalculatorTest. Test if calculation in new architecture returns similar results to the original implementation. WARNING: The values for comparison need to be adjusted if the input files are changed! This test can be used to get more detailed information, but it will always fail, because differences in voxelization accuracy, especially the ones caused by the change from double to float precission will not cause considerable deviations in the structure sizes, which correspond to considerable differences in the calculated DVHs. Even in double precission differences up to 0.005 between values from old and new implementation can occure. */ int DVHCalculatorComparisonTest(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: dose2 file name // 4: dose3 file name std::string RTSTRUCT_FILENAME; std::string RTDOSE_FILENAME; std::string RTDOSE2_FILENAME; std::string RTDOSE3_FILENAME; std::string COMPARISON_DVH_FOLDER; if (argc > 1) { RTSTRUCT_FILENAME = argv[1]; } if (argc > 2) { RTDOSE_FILENAME = argv[2]; } if (argc > 3) { RTDOSE2_FILENAME = argv[3]; } if (argc > 4) { RTDOSE3_FILENAME = argv[4]; } if (argc > 5) { COMPARISON_DVH_FOLDER = argv[5]; } OFCondition status; DcmFileFormat fileformat; /* read dicom-rt dose */ ::DRTDoseIOD rtdose1; 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; ::DRTDoseIOD rtdose2; io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator2(RTDOSE2_FILENAME.c_str()); DoseAccessorPointer doseAccessor2(doseAccessorGenerator2.generateDoseAccessor()); ::DRTDoseIOD rtdose3; io::dicom::DicomFileDoseAccessorGenerator doseAccessorGenerator3(RTDOSE3_FILENAME.c_str()); DoseAccessorPointer doseAccessor3(doseAccessorGenerator3.generateDoseAccessor()); double maxDifference = 0; double difference = 0; double minDifference = 1000; clock_t start(clock()); if (rtStructureSet->getNumberOfStructures() > 0) { for (int j = 0; j < rtStructureSet->getNumberOfStructures(); j++) { //create MaskAccessor boost::shared_ptr spOTBMaskAccessor = boost::make_shared(rtStructureSet->getStructure(j), doseAccessor1->getGeometricInfo()); spOTBMaskAccessor->updateMask(); MaskAccessorPointer spMaskAccessor(spOTBMaskAccessor); //create corresponding MaskedDoseIterator boost::shared_ptr spMaskedDoseIteratorTmp = boost::make_shared(spMaskAccessor, doseAccessor1); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); //store MaskAccessor rtStructSetMaskAccessorVec.push_back(spMaskAccessor); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor1->getDoseUID()); std::string dvhFileName = "dvh1"; std::string label = (rtStructureSet->getStructure(j))->getLabel(); dvhFileName.append(label); if (dvhFileName.find("/") != std::string::npos) { dvhFileName.replace(dvhFileName.find("/"), 1, ""); } std::cout << "=== Dose 1: " << label << "===" << std::endl; rttb::io::other::DVHTxtFileReader dvhOriginalReader = rttb::io::other::DVHTxtFileReader( COMPARISON_DVH_FOLDER + dvhFileName + ".txt"); rttb::core::DVH dvhOrig = *(dvhOriginalReader.generateDVH()); rttb::core::DVH::DataDifferentialType dvhOrigData = dvhOrig.getDataDifferential(); rttb::core::DVH dvh = *(calc.generateDVH()); rttb::core::DVH::DataDifferentialType dvhData = dvh.getDataDifferential(); CHECK_EQUAL(dvhOrig, dvh); CHECK_CLOSE(dvhOrig.getMaximum(), dvh.getMaximum(), errorConstant); CHECK_CLOSE(dvhOrig.getMinimum(), dvh.getMinimum(), errorConstant); CHECK_CLOSE(dvhOrig.getMean(), dvh.getMean(), errorConstant); rttb::core::DVH::DataDifferentialType::iterator it; rttb::core::DVH::DataDifferentialType::iterator itOrig; itOrig = dvhOrigData.begin(); - for (it = dvhData.begin(); it != dvhData.end() && itOrig != dvhOrigData.end(); ++it, ++itOrig) + for (it = dvhData.begin(); it != dvhData.end() || itOrig != dvhOrigData.end(); ++it, ++itOrig) { CHECK_CLOSE(*(itOrig), *(it), errorConstant); std::cout << std::setprecision(20) << "difference: " << abs(*(itOrig) - * (it)) << std::endl; difference = abs(*(itOrig) - * (it)); if (difference < minDifference) { minDifference = difference; } if (difference > maxDifference) { maxDifference = difference; } } } } clock_t finish(clock()); std::cout << std::setprecision(20) << "max(difference): " << maxDifference << std::endl; std::cout << std::setprecision(20) << "min(difference): " << minDifference << std::endl; std::cout << "DVH Calculation time: " << finish - start << " ms" << std::endl; maxDifference = 0; minDifference = 1000; clock_t start2(clock()); for (int j = 0; j < rtStructSetMaskAccessorVec.size(); j++) { //create corresponding MaskedDoseIterator boost::shared_ptr spMaskedDoseIteratorTmp = boost::make_shared(rtStructSetMaskAccessorVec.at(j), doseAccessor2); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor2->getDoseUID()); std::string dvhFileName = "dvh2"; std::string label = (rtStructureSet->getStructure(j))->getLabel(); dvhFileName.append(label); if (dvhFileName.find("/") != std::string::npos) { dvhFileName.replace(dvhFileName.find("/"), 1, ""); } std::cout << "=== Dose 2: " << label << "===" << std::endl; rttb::io::other::DVHTxtFileReader dvhOriginalReader = rttb::io::other::DVHTxtFileReader( COMPARISON_DVH_FOLDER + dvhFileName + ".txt"); rttb::core::DVH dvhOrig = *(dvhOriginalReader.generateDVH()); rttb::core::DVH::DataDifferentialType dvhOrigData = dvhOrig.getDataDifferential(); rttb::core::DVH dvh = *(calc.generateDVH()); rttb::core::DVH::DataDifferentialType dvhData = dvh.getDataDifferential(); CHECK_EQUAL(dvhOrig, dvh); CHECK_CLOSE(dvhOrig.getMaximum(), dvh.getMaximum(), errorConstant); CHECK_CLOSE(dvhOrig.getMinimum(), dvh.getMinimum(), errorConstant); CHECK_CLOSE(dvhOrig.getMean(), dvh.getMean(), errorConstant); rttb::core::DVH::DataDifferentialType::iterator it; rttb::core::DVH::DataDifferentialType::iterator itOrig; itOrig = dvhOrigData.begin(); - for (it = dvhData.begin(); it != dvhData.end() && itOrig != dvhOrigData.end(); ++it, ++itOrig) + for (it = dvhData.begin(); it != dvhData.end() || itOrig != dvhOrigData.end(); ++it, ++itOrig) { CHECK_CLOSE(*(itOrig), *(it), errorConstant); std::cout << std::setprecision(20) << "difference: " << abs(*(itOrig) - * (it)) << std::endl; difference = abs(*(itOrig) - * (it)); if (difference > 10) { std::cout << "large difference: " << abs(*(itOrig) - * (it)) << " = " << *(itOrig) << " - " << * (it) << std::endl; } if (difference < minDifference) { minDifference = difference; } if (difference > maxDifference) { maxDifference = difference; } } } clock_t finish2(clock()); std::cout << std::setprecision(20) << "max(difference): " << maxDifference << std::endl; std::cout << std::setprecision(20) << "min(difference): " << minDifference << std::endl; std::cout << "Reset dose 2, DVH Calculation time: " << finish2 - start2 << " ms" << std::endl; maxDifference = 0; minDifference = 1000; clock_t start3(clock()); for (int j = 0; j < rtStructSetMaskAccessorVec.size(); j++) { //create corresponding MaskedDoseIterator boost::shared_ptr spMaskedDoseIteratorTmp = boost::make_shared(rtStructSetMaskAccessorVec.at(j), doseAccessor3); DoseIteratorPointer spMaskedDoseIterator(spMaskedDoseIteratorTmp); rttb::core::DVHCalculator calc(spMaskedDoseIterator, (rtStructureSet->getStructure(j))->getUID(), doseAccessor3->getDoseUID()); std::string dvhFileName = "dvh3"; std::string label = (rtStructureSet->getStructure(j))->getLabel(); dvhFileName.append(label); if (dvhFileName.find("/") != std::string::npos) { dvhFileName.replace(dvhFileName.find("/"), 1, ""); } std::cout << "=== Dose 3: " << label << "===" << std::endl; rttb::io::other::DVHTxtFileReader dvhOriginalReader = rttb::io::other::DVHTxtFileReader( COMPARISON_DVH_FOLDER + dvhFileName + ".txt"); rttb::core::DVH dvhOrig = *(dvhOriginalReader.generateDVH()); rttb::core::DVH::DataDifferentialType dvhOrigData = dvhOrig.getDataDifferential(); rttb::core::DVH dvh = *(calc.generateDVH()); rttb::core::DVH::DataDifferentialType dvhData = dvh.getDataDifferential(); CHECK_EQUAL(dvhOrig, dvh); CHECK_CLOSE(dvhOrig.getMaximum(), dvh.getMaximum(), errorConstant); CHECK_CLOSE(dvhOrig.getMinimum(), dvh.getMinimum(), errorConstant); CHECK_CLOSE(dvhOrig.getMean(), dvh.getMean(), errorConstant); rttb::core::DVH::DataDifferentialType::iterator it; rttb::core::DVH::DataDifferentialType::iterator itOrig; itOrig = dvhOrigData.begin(); - for (it = dvhData.begin(); it != dvhData.end() && itOrig != dvhOrigData.end(); ++it, ++itOrig) + for (it = dvhData.begin(); it != dvhData.end() || itOrig != dvhOrigData.end(); ++it, ++itOrig) { CHECK_CLOSE(*(itOrig), *(it), errorConstant); std::cout << std::setprecision(20) << "difference: " << abs(*(itOrig) - * (it)) << std::endl; if (difference > 10) { std::cout << "large difference: " << abs(*(itOrig) - * (it)) << " = " << *(itOrig) << " - " << * (it) << std::endl; } difference = abs(*(itOrig) - * (it)); if (difference < minDifference) { minDifference = difference; } if (difference > maxDifference) { maxDifference = difference; } } } clock_t finish3(clock()); std::cout << std::setprecision(20) << "max(difference): " << maxDifference << std::endl; std::cout << std::setprecision(20) << "min(difference): " << minDifference << std::endl; std::cout << "Reset dose 3, DVH Calculation time: " << finish3 - start3 << " ms" << std::endl; RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.h b/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.h index 73db052..71912da 100644 --- a/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.h +++ b/testing/interpolation/MatchPointBinding/simpleRegistrationWorkflow.h @@ -1,64 +1,63 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision: 741 $ (last changed revision) // @date $Date: 2014-09-16 16:34:22 +0200 (Di, 16 Sep 2014) $ (last change date) // @author $Author: hentsch $ (last changed by) */ #ifndef __SIMPLE_REGISTRATION_WORKFLOW_H #define __SIMPLE_REGISTRATION_WORKFLOW_H #include -#include #include "registrationHelper.h" #include "mapImageMappingTask.h" #include "mapITKEuler3DMattesMIRegistrationAlgorithmTemplate.h" #include "mapExceptionObject.h" /*! @class simpleRegistrationWorkflow @brief implements a concrete registration algorithm of MatchPoint */ class simpleRegistrationWorkflow { public: typedef map::core::Registration<3, 3> Registration3D3DType; typedef Registration3D3DType::Pointer Registration3D3DPointer; typedef map::algorithm::boxed::ITKEuler3DMattesMIRegistrationAlgorithm AlgorithmTypeEuler; private: std::string _targetFilename; std::string _movingFilename; std::string _targetDirectory; std::string _movingDirectory; AppGlobals globals; AlgorithmTypeEuler::Pointer _spAlgorithmEuler; public: /*! @brief Constructor */ simpleRegistrationWorkflow(std::string targetFilename, std::string movingFilename, bool isDirectory = false); map::core::Registration<3, 3>::Pointer getRegistration(); const itk::Image* getTargetImage(); vnl_vector getRegistrationParameters(Registration3D3DPointer reg); protected: void initializeAndPerformRegistration(); }; #endif \ No newline at end of file