diff --git a/code/core/rttbGeometricInfo.cpp b/code/core/rttbGeometricInfo.cpp index 8bfcfeb..336aeb0 100644 --- a/code/core/rttbGeometricInfo.cpp +++ b/code/core/rttbGeometricInfo.cpp @@ -1,344 +1,345 @@ // ----------------------------------------------------------------------- // 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 #include #include "rttbGeometricInfo.h" namespace rttb { namespace core { void GeometricInfo::setSpacing(const SpacingVectorType3D& aSpacingVector) { _spacing = aSpacingVector; } const SpacingVectorType3D& GeometricInfo::getSpacing() const { return _spacing; } void GeometricInfo::setImagePositionPatient(const WorldCoordinate3D& aImagePositionPatient) { _imagePositionPatient = aImagePositionPatient; } const WorldCoordinate3D& GeometricInfo::getImagePositionPatient() const { return _imagePositionPatient; } void GeometricInfo::setOrientationMatrix(const OrientationMatrix& anOrientationMatrix) { _orientationMatrix = anOrientationMatrix; computeInvertOrientation(); } bool GeometricInfo::computeInvertOrientation() { typedef boost::numeric::ublas::permutation_matrix pmatrix; boost::numeric::ublas::matrix A(_orientationMatrix); // create a permutation matrix for the LU-factorization pmatrix pm(A.size1()); size_t res = boost::numeric::ublas::lu_factorize(A, pm); if (res != 0) { return false; } _invertedOrientationMatrix.assign(boost::numeric::ublas::identity_matrix (A.size1())); // backsubstitute to get the inverse boost::numeric::ublas::lu_substitute(A, pm, _invertedOrientationMatrix); return true; } const ImageOrientation GeometricInfo::getImageOrientationRow() const { ImageOrientation _imageOrientationRow(0); _imageOrientationRow(0) = _orientationMatrix(0, 0); _imageOrientationRow(1) = _orientationMatrix(1, 0); _imageOrientationRow(2) = _orientationMatrix(2, 0); return _imageOrientationRow; } const ImageOrientation GeometricInfo::getImageOrientationColumn() const { ImageOrientation _imageOrientationColumn(0); _imageOrientationColumn(0) = _orientationMatrix(0, 1); _imageOrientationColumn(1) = _orientationMatrix(1, 1); _imageOrientationColumn(2) = _orientationMatrix(2, 1); return _imageOrientationColumn; } void GeometricInfo::setPixelSpacingRow(const GridVolumeType aValue) { _spacing(0) = aValue; } const GridVolumeType GeometricInfo::getPixelSpacingRow() const { return _spacing(0); } void GeometricInfo::setPixelSpacingColumn(const GridVolumeType aValue) { _spacing(1) = aValue; } const GridVolumeType GeometricInfo::getPixelSpacingColumn() const { return _spacing(1); } void GeometricInfo::setSliceThickness(const GridVolumeType aValue) { _spacing(2) = aValue; } const GridVolumeType GeometricInfo::getSliceThickness() const { return _spacing(2); } void GeometricInfo::setImageSize(const ImageSize& aSize) { setNumColumns(aSize(0)); setNumRows(aSize(1)); setNumSlices(aSize(2)); } const ImageSize GeometricInfo::getImageSize() const { return ImageSize(getNumColumns(), getNumRows(), getNumSlices()); } void GeometricInfo::setNumColumns(const VoxelGridDimensionType aValue) { _numberOfColumns = aValue; } const VoxelGridDimensionType GeometricInfo::getNumColumns() const { return _numberOfColumns; } void GeometricInfo::setNumRows(const VoxelGridDimensionType aValue) { _numberOfRows = aValue; } const VoxelGridDimensionType GeometricInfo::getNumRows() const { return _numberOfRows; } void GeometricInfo::setNumSlices(const VoxelGridDimensionType aValue) { _numberOfFrames = aValue; } const VoxelGridDimensionType GeometricInfo::getNumSlices() const { return _numberOfFrames; } bool operator==(const GeometricInfo& gInfo, const GeometricInfo& gInfo1) { return (gInfo.getImagePositionPatient() == gInfo1.getImagePositionPatient() && gInfo.getOrientationMatrix() == gInfo1.getOrientationMatrix() && gInfo.getSpacing() == gInfo1.getSpacing() && gInfo.getNumColumns() == gInfo1.getNumColumns() && gInfo.getNumRows() == gInfo1.getNumRows() && gInfo.getNumSlices() == gInfo1.getNumSlices()); } bool GeometricInfo::equalsAlmost(const GeometricInfo& another, double errorConstant /*= 1e-5*/) const { return (getImagePositionPatient().equalsAlmost(another.getImagePositionPatient(), errorConstant) && getOrientationMatrix().equalsAlmost(another.getOrientationMatrix(), errorConstant) && getSpacing().equalsAlmost(another.getSpacing(), errorConstant) && getNumColumns() == another.getNumColumns() && getNumRows() == another.getNumRows() && getNumSlices() == another.getNumSlices()); } bool GeometricInfo::worldCoordinateToGeometryCoordinate(const WorldCoordinate3D& aWorldCoordinate, DoubleVoxelGridIndex3D& aIndex) const { WorldCoordinate3D distanceToIP; distanceToIP = aWorldCoordinate - _imagePositionPatient; boost::numeric::ublas::vector result = boost::numeric::ublas::prod( _invertedOrientationMatrix, distanceToIP); boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_div(result, _spacing); aIndex = DoubleVoxelGridIndex3D(resultS(0), resultS(1), resultS(2)); //check if it is inside VoxelGridIndex3D indexInt = VoxelGridIndex3D(GridIndexType(aIndex(0) + 0.5), GridIndexType(aIndex(1) + 0.5), GridIndexType(aIndex(2) + 0.5)); return isInside(indexInt); } bool GeometricInfo::worldCoordinateToIndex(const WorldCoordinate3D& aWorldCoordinate, VoxelGridIndex3D& aIndex) const { DoubleVoxelGridIndex3D doubleIndex; bool inside = worldCoordinateToGeometryCoordinate(aWorldCoordinate, doubleIndex); aIndex = VoxelGridIndex3D(GridIndexType(doubleIndex(0) + 0.5), GridIndexType(doubleIndex(1) + 0.5), GridIndexType(doubleIndex(2) + 0.5)); return inside; } bool GeometricInfo::geometryCoordinateToWorldCoordinate(const DoubleVoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const { boost::numeric::ublas::vector resultS = boost::numeric::ublas::element_prod( aIndex, _spacing); boost::numeric::ublas::vector result = boost::numeric::ublas::prod( _orientationMatrix, resultS); aWorldCoordinate = result + _imagePositionPatient; VoxelGridIndex3D indexInt = VoxelGridIndex3D(GridIndexType(aIndex(0) + 0.5), GridIndexType(aIndex(1) + 0.5), GridIndexType(aIndex(2) + 0.5)); return isInside(indexInt); } bool GeometricInfo::indexToWorldCoordinate(const VoxelGridIndex3D& aIndex, WorldCoordinate3D& aWorldCoordinate) const { DoubleVoxelGridIndex3D indexDouble = DoubleVoxelGridIndex3D(aIndex(0) - 0.5, aIndex(1) - 0.5, aIndex(2) - 0.5); return geometryCoordinateToWorldCoordinate(indexDouble, aWorldCoordinate); } bool GeometricInfo::isInside(const VoxelGridIndex3D& aIndex) const { return (aIndex(0) >= 0 && aIndex(1) >= 0 && aIndex(2) >= 0 && aIndex(0) < static_cast(_numberOfColumns) && aIndex(1) < static_cast(_numberOfRows) && aIndex(2) < static_cast(_numberOfFrames)); } bool GeometricInfo::isInside(const WorldCoordinate3D& aWorldCoordinate) { VoxelGridIndex3D currentIndex; return (worldCoordinateToIndex(aWorldCoordinate, currentIndex)); } const GridSizeType GeometricInfo::getNumberOfVoxels() const { - return _numberOfRows * _numberOfColumns * _numberOfFrames; + GridSizeType nVoxels = static_cast(_numberOfRows * _numberOfColumns * _numberOfFrames); + return nVoxels; } bool GeometricInfo::convert(const VoxelGridID& gridID, VoxelGridIndex3D& gridIndex) const { if (validID(gridID)) { gridIndex(0) = gridID % getNumColumns(); VoxelGridID tempID = (gridID - gridIndex.x()) / getNumColumns(); gridIndex(1) = tempID % getNumRows(); gridIndex(2) = (tempID - gridIndex.y()) / getNumRows(); return true; } return false; } bool GeometricInfo::convert(const VoxelGridIndex3D& gridIndex, VoxelGridID& gridID) const { if ((gridIndex.x() >= static_cast(getNumColumns())) || (gridIndex.y() >= static_cast(getNumRows())) || (gridIndex.z() >= static_cast(getNumSlices()))) { return false; } else { gridID = gridIndex.z() * getNumColumns() * getNumRows() + gridIndex.y() * getNumColumns() + gridIndex.x(); return validID(gridID); } } bool GeometricInfo::validID(const VoxelGridID aID) const { return (aID >= 0 && aID < getNumberOfVoxels()); } bool GeometricInfo::validIndex(const VoxelGridIndex3D& aIndex) const { VoxelGridID aID; if (!convert(aIndex, aID)) { return false; } else { return validID(aID); } } std::ostream& operator<<(std::ostream& s, const GeometricInfo& anGeometricInfo) { s << "[ " << anGeometricInfo.getImagePositionPatient() << "; " << anGeometricInfo.getOrientationMatrix() << "; " << anGeometricInfo.getSpacing() << "; " << "; " << anGeometricInfo.getNumColumns() << "; " << anGeometricInfo.getNumRows() << "; " << anGeometricInfo.getNumSlices() << " ]"; return s; } }//end namespace core }//end namespace rttb diff --git a/code/io/itk/itkDoseAccessorImageFilter.cpp b/code/io/itk/itkDoseAccessorImageFilter.cpp index 183d2c8..748f776 100644 --- a/code/io/itk/itkDoseAccessorImageFilter.cpp +++ b/code/io/itk/itkDoseAccessorImageFilter.cpp @@ -1,89 +1,83 @@ // ----------------------------------------------------------------------- // 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 "itkDoseAccessorImageFilter.h" #include "itkImageRegionIterator.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkProgressReporter.h" namespace itk { /** * Constructor */ DoseAccessorImageFilter ::DoseAccessorImageFilter() { this->SetNumberOfRequiredInputs(1); } void DoseAccessorImageFilter ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType threadId) { ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); - const unsigned int numberOfInputImages = - static_cast< unsigned int >(this->GetNumberOfIndexedInputs()); - - const unsigned int numberOfOutputImages = - static_cast< unsigned int >(this->GetNumberOfIndexedOutputs()); - typedef ImageRegionConstIteratorWithIndex< InputImageType > ImageRegionConstIteratorType; typedef ImageRegionIterator< OutputImageType > OutputImageRegionIteratorType; InputImagePointer inputPtr = dynamic_cast< InputImageType* >(ProcessObject::GetInput(0)); ImageRegionConstIteratorType inputItr; if (inputPtr) { inputItr = ImageRegionConstIteratorType(inputPtr, outputRegionForThread); } OutputImagePointer outputPtr = dynamic_cast< OutputImageType* >(ProcessObject::GetOutput(0)); OutputImageRegionIteratorType outputItr; if (outputPtr) { outputItr = OutputImageRegionIteratorType(outputPtr, outputRegionForThread); } if (inputPtr && outputPtr) { while (!(outputItr.IsAtEnd())) { ImageRegionConstIteratorType::IndexType index = inputItr.GetIndex(); rttb::VoxelGridIndex3D doseIndex(index[0], index[1], index[2]); outputItr.Set(m_Accessor->getValueAt(doseIndex)); ++outputItr; ++inputItr; progress.CompletedPixel(); } } } } // end namespace itk diff --git a/code/io/other/rttbDVHTxtFileReader.cpp b/code/io/other/rttbDVHTxtFileReader.cpp index 366a631..79b90c2 100644 --- a/code/io/other/rttbDVHTxtFileReader.cpp +++ b/code/io/other/rttbDVHTxtFileReader.cpp @@ -1,226 +1,219 @@ // ----------------------------------------------------------------------- // 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 #include #include #include #include "rttbDVHTxtFileReader.h" #include "rttbInvalidParameterException.h" #include "rttbBaseType.h" namespace rttb { namespace io { namespace other { DVHTxtFileReader::DVHTxtFileReader(FileNameString aFileName) { _fileName = aFileName; _resetFile = true; } DVHTxtFileReader::~DVHTxtFileReader() {} void DVHTxtFileReader::resetFileName(FileNameString aFileName) { _fileName = aFileName; _resetFile = true; } void DVHTxtFileReader::createDVH() { std::ifstream dvh_ifstr(this->_fileName.c_str(), std::ios::in); std::string structureLabel; std::string dvhType; unsigned int numberOfBins; DoseTypeGy prescribedDose; double estimated_max_dose_prescribed_dose_ratio; - int voxelsInStructure; std::deque dataDifferential; std::deque dataCumulative; DoseTypeGy deltaD = 0; DoseVoxelVolumeType deltaV = 0; IDType strID; IDType doseID; if (!dvh_ifstr.is_open()) { throw core::InvalidParameterException("DVH file name invalid: could not open the dvh file!"); } else { bool data_begin = false; while (!dvh_ifstr.eof()) { std::string line; std::getline(dvh_ifstr, line); if (line.find("DVH Data:") != std::string::npos) { data_begin = true; } if (data_begin && (line.find(",") != std::string::npos)) { std::stringstream ss(line.substr(line.find(",") + 1)); DoseCalcType dvh_i; ss >> dvh_i; if (dvhType == "DIFFERENTIAL") { dataDifferential.push_back(dvh_i); } else if (dvhType == "CUMULATIVE") { dataCumulative.push_back(dvh_i); } } if (line.find("DeltaD: ") != std::string::npos) { std::stringstream ss(line.substr(8)); ss >> deltaD; } if (line.find("DeltaV: ") != std::string::npos) { std::stringstream ss(line.substr(8)); ss >> deltaV; } if (line.find("StructureID: ") != std::string::npos) { std::stringstream ss(line.substr(13)); ss >> strID; } if (line.find("DoseID: ") != std::string::npos) { std::stringstream ss(line.substr(8)); ss >> doseID; } if (line.find("Number of bins: ") != std::string::npos) { std::stringstream ss(line.substr(16)); ss >> numberOfBins; } if (line.find("Structure Label: ") != std::string::npos) { std::stringstream ss(line.substr(17)); ss >> structureLabel; } if (line.find("DVH Type: ") != std::string::npos) { std::stringstream ss(line.substr(10)); ss >> dvhType; } if (line.find("Prescribed Dose: ") != std::string::npos) { std::stringstream ss(line.substr(17)); ss >> prescribedDose; } if (line.find("Estimated_max_dose_prescribed_dose_ratio: ") != std::string::npos) { std::stringstream ss(line.substr(42)); ss >> estimated_max_dose_prescribed_dose_ratio; } - - if (line.find("Voxels In Structure: ") != std::string::npos) - { - std::stringstream ss(line.substr(21)); - ss >> voxelsInStructure; - } } } numberOfBins = static_cast(std::max(dataDifferential.size(), dataCumulative.size())); if (dvhType == "CUMULATIVE") { DoseCalcType differentialDVHi = 0; std::deque::iterator it; for (it = dataCumulative.begin(); it != dataCumulative.end(); ++it) { if (dataDifferential.size() == numberOfBins - 1) { differentialDVHi = *it; } else { differentialDVHi = *it - *(it + 1); } dataDifferential.push_back(differentialDVHi); } } if (numberOfBins == 0) { throw core::InvalidParameterException("Invalid dvh file: empty dvh data!"); } if (deltaD == 0) { deltaD = prescribedDose * estimated_max_dose_prescribed_dose_ratio / numberOfBins; } if (deltaV == 0) { deltaV = 0.027; } if (deltaD == 0 || deltaV == 0) { throw core::InvalidParameterException("Invalid dvh file: deltaD or deltaV must not be zero!"); } _dvh = boost::make_shared(dataDifferential, deltaD, deltaV, strID, doseID); _resetFile = false; } DVHTxtFileReader::DVHPointer DVHTxtFileReader::generateDVH() { if (_resetFile) { this->createDVH(); } return _dvh; } }//end namespace other }//end namespace io }//end namespace rttb diff --git a/code/io/other/rttbDoseStatisticsXMLWriter.cpp b/code/io/other/rttbDoseStatisticsXMLWriter.cpp index c3693b9..9529669 100644 --- a/code/io/other/rttbDoseStatisticsXMLWriter.cpp +++ b/code/io/other/rttbDoseStatisticsXMLWriter.cpp @@ -1,279 +1,278 @@ // ----------------------------------------------------------------------- // 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 #include "rttbDoseStatisticsXMLWriter.h" #include "rttbInvalidParameterException.h" #include "rttbDataNotAvailableException.h" namespace rttb { namespace io { namespace other { static const std::string xmlattrXTag = ".x"; static const std::string xmlattrNameTag = ".name"; static const std::string statisticsTag = "statistics.results"; static const std::string propertyTag = "property"; static const std::string columnSeparator = "@"; boost::property_tree::ptree writeDoseStatistics(DoseStatisticsPtr aDoseStatistics) { using boost::property_tree::ptree; ptree pt; ptree numberOfVoxelsNode = createNodeWithNameAttribute(aDoseStatistics->getNumberOfVoxels(), "numberOfVoxels"); pt.add_child(statisticsTag + "." + propertyTag, numberOfVoxelsNode); ptree volumeNode = createNodeWithNameAttribute(aDoseStatistics->getVolume(), "volume"); pt.add_child(statisticsTag + "." + propertyTag, volumeNode); ptree minimumNode = createNodeWithNameAttribute(aDoseStatistics->getMinimum(), "minimum"); auto minimumPositions = aDoseStatistics->getMinimumPositions(); std::vector >::iterator pairItMin = minimumPositions->begin(); int count = 0; for (; pairItMin != minimumPositions->end() && count < 100; ++pairItMin) //output max. 100 minimum { VoxelGridID voxelID = pairItMin->second; ptree voxelMinPositions; voxelMinPositions.add("voxelGridID", voxelID); minimumNode.add_child("voxel", voxelMinPositions); count++; } pt.add_child(statisticsTag + "." + propertyTag, minimumNode); ptree maximumNode = createNodeWithNameAttribute(aDoseStatistics->getMaximum(), "maximum"); auto maximumPositions = aDoseStatistics->getMaximumPositions(); std::vector >::iterator pairItMax = maximumPositions->begin(); count = 0; for (; pairItMax != maximumPositions->end() && count < 100; ++pairItMax) //output max. 100 maximum { VoxelGridID voxelID = pairItMax->second; ptree voxelMaxPositions; voxelMaxPositions.add("voxelGridID", voxelID); maximumNode.add_child("voxel", voxelMaxPositions); count++; } pt.add_child(statisticsTag + "." + propertyTag, maximumNode); ptree meanNode = createNodeWithNameAttribute(aDoseStatistics->getMean(), "mean"); pt.add_child(statisticsTag + "." + propertyTag, meanNode); ptree sdNode = createNodeWithNameAttribute(aDoseStatistics->getStdDeviation(), "standardDeviation"); pt.add_child(statisticsTag + "." + propertyTag, sdNode); ptree varianceNode = createNodeWithNameAttribute(aDoseStatistics->getVariance(), "variance"); pt.add_child(statisticsTag + "." + propertyTag, varianceNode); double absoluteVolume = aDoseStatistics->getVolume(); double referenceDose = aDoseStatistics->getReferenceDose(); rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType AllVx = aDoseStatistics->getAllVx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllDx = aDoseStatistics->getAllDx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getAllMOHx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getAllMOCx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getAllMaxOHx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getAllMinOCx(); rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType::iterator vxIt; rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType::iterator it; - for (it = AllDx.begin(); it != AllDx.end(); it++) + for (it = AllDx.begin(); it != AllDx.end(); ++it) { ptree DxNode = createNodeWithNameAndXAttribute((*it).second, "Dx", static_cast((*it).first / absoluteVolume * 100)); pt.add_child(statisticsTag + "." + propertyTag, DxNode); } - for (vxIt = AllVx.begin(); vxIt != AllVx.end(); vxIt++) + for (vxIt = AllVx.begin(); vxIt != AllVx.end(); ++vxIt) { ptree VxNode = createNodeWithNameAndXAttribute((*vxIt).second, "Vx", static_cast((*vxIt).first / referenceDose * 100)); pt.add_child(statisticsTag + "." + propertyTag, VxNode); } - for (it = AllMOHx.begin(); it != AllMOHx.end(); it++) + for (it = AllMOHx.begin(); it != AllMOHx.end(); ++it) { ptree mohxNode = createNodeWithNameAndXAttribute((*it).second, "MOHx", static_cast((*it).first / absoluteVolume * 100)); pt.add_child(statisticsTag + "." + propertyTag, mohxNode); } - for (it = AllMOCx.begin(); it != AllMOCx.end(); it++) + for (it = AllMOCx.begin(); it != AllMOCx.end(); ++it) { ptree mocxNode = createNodeWithNameAndXAttribute((*it).second, "MOCx", static_cast((*it).first / absoluteVolume * 100)); pt.add_child(statisticsTag + "." + propertyTag, mocxNode); } - for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); it++) + for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); ++it) { ptree maxOhxNode = createNodeWithNameAndXAttribute((*it).second, "MaxOHx", static_cast((*it).first / absoluteVolume * 100)); pt.add_child(statisticsTag + "." + propertyTag, maxOhxNode); } - for (it = AllMinOCx.begin(); it != AllMinOCx.end(); it++) + for (it = AllMinOCx.begin(); it != AllMinOCx.end(); ++it) { ptree minOCxNode = createNodeWithNameAndXAttribute((*it).second, "MinOCx", static_cast((*it).first / absoluteVolume * 100)); pt.add_child(statisticsTag + "." + propertyTag, minOCxNode); } return pt; } void writeDoseStatistics(DoseStatisticsPtr aDoseStatistics, FileNameString aFileName) { boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics); try { boost::property_tree::xml_parser::write_xml(aFileName, pt, std::locale(), boost::property_tree::xml_writer_make_settings('\t', 1)); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml failed: xml_parser_error!"); } } XMLString writerDoseStatisticsToString(DoseStatisticsPtr aDoseStatistics) { boost::property_tree::ptree pt = writeDoseStatistics(aDoseStatistics); std::stringstream sstr; try { boost::property_tree::xml_parser::write_xml(sstr, pt, boost::property_tree::xml_writer_make_settings('\t', 1)); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw core::InvalidParameterException("Write xml to string failed: xml_parser_error!"); } return sstr.str(); } StatisticsString writerDoseStatisticsToTableString(DoseStatisticsPtr aDoseStatistics) { std::stringstream sstr; sstr << aDoseStatistics->getVolume() * 1000 << columnSeparator; // cm3 to mm3 sstr << aDoseStatistics->getMaximum() << columnSeparator; sstr << aDoseStatistics->getMinimum() << columnSeparator; sstr << aDoseStatistics->getMean() << columnSeparator; sstr << aDoseStatistics->getStdDeviation() << columnSeparator; sstr << aDoseStatistics->getVariance() << columnSeparator; - double absoluteVolume = aDoseStatistics->getVolume(); rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType AllVx = aDoseStatistics->getAllVx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllDx = aDoseStatistics->getAllDx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOHx = aDoseStatistics->getAllMOHx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMOCx = aDoseStatistics->getAllMOCx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMaxOHx = aDoseStatistics->getAllMaxOHx(); rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType AllMinOCx = aDoseStatistics->getAllMinOCx(); rttb::algorithms::DoseStatistics::DoseToVolumeFunctionType::iterator vxIt; rttb::algorithms::DoseStatistics::VolumeToDoseFunctionType::iterator it; - for (it = AllDx.begin(); it != AllDx.end(); it++) + for (it = AllDx.begin(); it != AllDx.end(); ++it) { sstr << (*it).second << columnSeparator; } - for (vxIt = AllVx.begin(); vxIt != AllVx.end(); vxIt++) + for (vxIt = AllVx.begin(); vxIt != AllVx.end(); ++vxIt) { // *1000 because of conversion cm3 to mm3 sstr << (*vxIt).second * 1000 << columnSeparator; } - for (it = AllMOHx.begin(); it != AllMOHx.end(); it++) + for (it = AllMOHx.begin(); it != AllMOHx.end(); ++it) { sstr << (*it).second << columnSeparator; } - for (it = AllMOCx.begin(); it != AllMOCx.end(); it++) + for (it = AllMOCx.begin(); it != AllMOCx.end(); ++it) { sstr << (*it).second << columnSeparator; } - for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); it++) + for (it = AllMaxOHx.begin(); it != AllMaxOHx.end(); ++it) { sstr << (*it).second << columnSeparator; } - for (it = AllMinOCx.begin(); it != AllMinOCx.end(); it++) + for (it = AllMinOCx.begin(); it != AllMinOCx.end(); ++it) { sstr << (*it).second << columnSeparator; } return sstr.str(); } boost::property_tree::ptree createNodeWithNameAttribute(DoseTypeGy doseValue, const std::string& attributeName) { boost::property_tree::ptree node; node.put("", doseValue); node.put(xmlattrNameTag, attributeName); return node; } boost::property_tree::ptree createNodeWithNameAndXAttribute(DoseTypeGy doseValue, const std::string& attributeName, int xValue) { boost::property_tree::ptree node; node.put("", doseValue); node.put(xmlattrNameTag, attributeName); node.put(xmlattrXTag, xValue); return node; } }//end namespace other }//end namespace io }//end namespace rttb \ No newline at end of file diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp index 186b08f..faf697e 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/boost/rttbBoostMask.cpp @@ -1,496 +1,496 @@ #include #include #include #include #include #include #include #include #include "rttbBoostMask.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace masks { namespace boost { BoostMask::BoostMask(BoostMask::GeometricInfoPointer aDoseGeoInfo, BoostMask::StructPointer aStructure, bool strict) : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _voxelInStructure(::boost::make_shared()) { _isUpToDate = false; if (!_geometricInfo) { throw rttb::core::NullPointerException("Error: Geometric info is NULL!"); } else if (!_structure) { throw rttb::core::NullPointerException("Error: Structure is NULL!"); } } void BoostMask::calcMask() { if (hasSelfIntersections()) { if (_strict) { throw rttb::core::InvalidParameterException("Error: structure has self intersections!"); } else { std::cerr << _structure->getLabel() << " has self intersections! It may cause errors in the voxelization results!" << std::endl; } } std::vector intersectionSlicePolygonsVector;//store the polygons of intersection slice of each index z unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); //For each dose slice, get intersection structure slice and the contours on the structure slice for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) { BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(indexZ); BoostMask::BoostPolygonVector::iterator it; for (it = boostPolygons.begin(); it != boostPolygons.end(); ++it) { ::boost::geometry::correct(*it); } intersectionSlicePolygonsVector.push_back(boostPolygons); } /* Use simple nearest neighbour interpolation (NNI) if dose and structure has different z spacing: If dose slice (indexZ) has no intersection with structure slice, but the last and the next do, that means the dose slice lies between two structure slices -> use the last slice intersection contours for this slice */ for (unsigned int indexZ = 1; indexZ < nSlices - 1; indexZ++) { if (intersectionSlicePolygonsVector.at(indexZ).size() == 0 && intersectionSlicePolygonsVector.at(indexZ - 1).size() > 0 && intersectionSlicePolygonsVector.at(indexZ + 1).size() > 0) { intersectionSlicePolygonsVector.at(indexZ) = intersectionSlicePolygonsVector.at(indexZ - 1); } } for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) { BoostMask::BoostPolygonVector intersectionSlicePolygons = intersectionSlicePolygonsVector.at( indexZ); //Get bounding box of this dose slice VoxelIndexVector voxelList = getBoundingBox(indexZ, intersectionSlicePolygons); rttb::VoxelGridIndex3D gridIndex3D0 = voxelList.at(0); rttb::VoxelGridIndex3D gridIndex3D1 = voxelList.at(1); for (unsigned int indexX = gridIndex3D0[0]; indexX <= gridIndex3D1[0]; indexX++) { for (unsigned int indexY = gridIndex3D0[1]; indexY <= gridIndex3D1[1]; indexY++) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = indexX; currentIndex[1] = indexY; currentIndex[2] = indexZ; rttb::VoxelGridID gridID; _geometricInfo->convert(currentIndex, gridID); //Get intersection polygons of the dose voxel and the structure BoostPolygonDeque polygons = getIntersections(currentIndex, intersectionSlicePolygons); //Calc areas of all intersection polygons double intersectionArea = calcArea(polygons); double voxelSize = _geometricInfo->getSpacing()[0] * _geometricInfo->getSpacing()[1]; if (intersectionArea > 0) { double volumeFraction = intersectionArea / voxelSize; if (volumeFraction > 1 && (volumeFraction - 1) <= 0.0001) { volumeFraction = 1; } core::MaskVoxel maskVoxel = core::MaskVoxel(gridID, (volumeFraction)); _voxelInStructure->push_back(maskVoxel);//push back the mask voxel in structure } } } } sort(_voxelInStructure->begin(), _voxelInStructure->end()); } bool BoostMask::hasSelfIntersections() { bool hasSelfIntersection = false; unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); for (unsigned int indexZ = 0; indexZ < nSlices; indexZ++) { rttb::VoxelGridIndex3D currentIndex(0, 0, indexZ); BoostMask::BoostPolygonVector boostPolygons = getIntersectionSlicePolygons(currentIndex[2]); BoostMask::BoostPolygonVector::iterator it; BoostMask::BoostPolygonVector::iterator it2; for (it = boostPolygons.begin(); it != boostPolygons.end(); ++it) { ::boost::geometry::correct(*it); //test if polygon has self intersection if (::boost::geometry::detail::overlay::has_self_intersections(*it, false)) { hasSelfIntersection = true; std::cerr << _structure->getLabel() << " has self intersection polygon! Slice: " << indexZ << ". " << std::endl; break; } //test if the polygons on the same slice has intersection for (it2 = boostPolygons.begin(); it2 != boostPolygons.end() && it2 != it; ++it2) { ::boost::geometry::correct(*it2); BoostPolygonDeque intersection; ::boost::geometry::intersection(*it, *it2, intersection); if (intersection.size() > 0) { //if no donut if (!(::boost::geometry::within(*it, *it2)) && !(::boost::geometry::within(*it2, *it))) { hasSelfIntersection = true; std::cerr << _structure->getLabel() << ": Two polygons on the same slice has intersection! Slice: " << indexZ << "." << std::endl; break; } } } } } return hasSelfIntersection; } /*Get the 4 voxel index of the bounding box of the polygon in the slice with sliceNumber*/ BoostMask::VoxelIndexVector BoostMask::getBoundingBox(unsigned int sliceNumber, const BoostPolygonVector& intersectionSlicePolygons) { unsigned int nSlices = static_cast(_geometricInfo->getNumSlices()); if (sliceNumber < 0 || sliceNumber >= nSlices) { throw rttb::core::InvalidParameterException(std::string("Error: slice number is invalid!")); } rttb::VoxelGridIndex3D currentIndex(0, 0, sliceNumber); double xMin = 0; double yMin = 0; double xMax = 0; double yMax = 0; BoostPolygonVector::const_iterator it; for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) { ::boost::geometry::model::box box; ::boost::geometry::envelope(*it, box); BoostPoint2D minPoint = box.min_corner(); BoostPoint2D maxPoint = box.max_corner(); if (it == intersectionSlicePolygons.begin()) { xMin = minPoint.x(); yMin = minPoint.y(); xMax = maxPoint.x(); yMax = maxPoint.y(); } xMin = std::min(xMin, minPoint.x()); yMin = std::min(yMin, minPoint.y()); xMax = std::max(xMax, maxPoint.x()); yMax = std::max(yMax, maxPoint.y()); } rttb::WorldCoordinate3D nullWorldCoord; _geometricInfo->indexToWorldCoordinate(VoxelGridIndex3D(0, 0, 0), nullWorldCoord); rttb::WorldCoordinate3D minWorldCoord(xMin, yMin, nullWorldCoord.z()); rttb::WorldCoordinate3D maxWorldCoord(xMax, yMax, nullWorldCoord.z()); rttb::VoxelGridIndex3D index0; rttb::VoxelGridIndex3D index1; _geometricInfo->worldCoordinateToIndex(minWorldCoord, index0); _geometricInfo->worldCoordinateToIndex(maxWorldCoord, index1); VoxelIndexVector voxelList; voxelList.push_back(index0); voxelList.push_back(index1); return voxelList; } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMask::BoostPolygonDeque BoostMask::getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) { BoostMask::BoostPolygonDeque polygonDeque; BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); ::boost::geometry::correct(voxelPolygon); BoostPolygonVector::const_iterator it; for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) { BoostPolygon2D contour = *it; ::boost::geometry::correct(contour); BoostPolygonDeque intersection; ::boost::geometry::intersection(voxelPolygon, contour, intersection); polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end()); } return polygonDeque; } /*Calculate the intersection area*/ double BoostMask::calcArea(const BoostPolygonDeque& aPolygonDeque) { double area = 0; BoostPolygonDeque::const_iterator it; for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) { area += ::boost::geometry::area(*it); } return area; } VoxelGridIndex3D BoostMask::getGridIndex3D(const core::MaskVoxel& aMaskVoxel) { rttb::VoxelGridIndex3D gridIndex3D; _geometricInfo->convert(aMaskVoxel.getVoxelGridID(), gridIndex3D); return gridIndex3D; } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { if (!_isUpToDate) { calcMask(); } return _voxelInStructure; } BoostMask::BoostRing2D BoostMask::convert(const rttb::PolygonType& aRTTBPolygon) { BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); if (i == 0) { firstPoint = boostPoint; } ::boost::geometry::append(polygon2D, boostPoint); } ::boost::geometry::append(polygon2D, firstPoint); return polygon2D; } BoostMask::BoostPolygonVector BoostMask::getIntersectionSlicePolygons( const rttb::GridIndexType aVoxelGridIndexZ) { BoostMask::BoostRingVector boostRingVector; rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); //get the polygons in the slice and convert to boost rings rttb::PolygonSequenceType::iterator it; for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) { PolygonType rttbPolygon = *it; if (!rttbPolygon.empty()) { double polygonZCoor = rttbPolygon.at(0)[2]; rttb::WorldCoordinate3D polygonPoint(0, 0, polygonZCoor); rttb::VoxelGridIndex3D polygonPointIndex3D; _geometricInfo->worldCoordinateToIndex(polygonPoint, polygonPointIndex3D); //if the z if (aVoxelGridIndexZ == polygonPointIndex3D[2]) { boostRingVector.push_back(convert(rttbPolygon)); } } } return checkDonutAndConvert(boostRingVector); } BoostMask::BoostRing2D BoostMask::get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D) { BoostMask::BoostRing2D polygon; rttb::WorldCoordinate3D rttbPoint; _geometricInfo->indexToWorldCoordinate(aVoxelGrid3D, rttbPoint); BoostPoint2D point1(rttbPoint[0], rttbPoint[1]); ::boost::geometry::append(polygon, point1); double xSpacing = _geometricInfo->getSpacing()[0]; double ySpacing = _geometricInfo->getSpacing()[1]; BoostPoint2D point2(rttbPoint[0] + xSpacing, rttbPoint[1]); ::boost::geometry::append(polygon, point2); BoostPoint2D point3(rttbPoint[0] + xSpacing, rttbPoint[1] + ySpacing); ::boost::geometry::append(polygon, point3); BoostPoint2D point4(rttbPoint[0], rttbPoint[1] + ySpacing); ::boost::geometry::append(polygon, point4); ::boost::geometry::append(polygon, point1); return polygon; } BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector) { //check donut BoostMask::BoostRingVector::const_iterator it1; BoostMask::BoostRingVector::const_iterator it2; BoostMask::BoostPolygonVector boostPolygonVector; std::vector donutIndexVector;//store the outer and inner ring index BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon //Get donut index and donut polygon unsigned int index1 = 0; - for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) + for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, ++index1) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } //if not jet, check now if (!it1IsDonut) { bool it2IsDonut = false; unsigned int index2 = 0; - for (it2 = aRingVector.begin(); it2 != aRingVector.end(); it2++, index2++) + for (it2 = aRingVector.begin(); it2 != aRingVector.end(); ++it2, ++index2) { if (it2 != it1) { BoostMask::BoostPolygon2D polygon2D; if (::boost::geometry::within(*it1, *it2)) { ::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring it2IsDonut = true; } //if donut else if (::boost::geometry::within(*it2, *it1)) { ::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring it2IsDonut = true; } if (it2IsDonut) { donutIndexVector.push_back(index1); donutIndexVector.push_back(index2); donutVector.push_back(polygon2D);//store donut polygon break;//Only store the first donut! } } } } } //Store no donut polygon to boostPolygonVector index1 = 0; - for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) + for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, ++index1) { bool it1IsDonut = false; //check if the ring is the outer or inner of a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } if (!it1IsDonut) { BoostMask::BoostPolygon2D polygon2D; ::boost::geometry::append(polygon2D, *it1); boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut } } //Append donut polygon to boostPolygonVector BoostMask::BoostPolygonVector::iterator itDonut; - for (itDonut = donutVector.begin(); itDonut != donutVector.end(); itDonut++) + for (itDonut = donutVector.begin(); itDonut != donutVector.end(); ++itDonut) { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } } } } \ No newline at end of file diff --git a/demoapps/BioModelCalc/BioModelCalcApplicationData.cpp b/demoapps/BioModelCalc/BioModelCalcApplicationData.cpp index 545efc8..ea35a91 100644 --- a/demoapps/BioModelCalc/BioModelCalcApplicationData.cpp +++ b/demoapps/BioModelCalc/BioModelCalcApplicationData.cpp @@ -1,56 +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 "BioModelCalcApplicationData.h" namespace rttb { namespace apps { namespace bioModelCalc { ApplicationData:: ApplicationData() { this->Reset(); } void ApplicationData:: Reset() { _doseFileName = ""; _outputFileName = ""; + _doseScaling = 1.0; } void populateAppData(boost::shared_ptr argParser, ApplicationData& appData) { appData._doseFileName = argParser->get(argParser->OPTION_DOSE_FILE); appData._doseLoadStyle = argParser->get >(argParser->OPTION_LOAD_STYLE); appData._doseScaling = argParser->get(argParser->OPTION_DOSE_SCALING); appData._outputFileName = argParser->get(argParser->OPTION_OUTPUT_FILE); appData._model = argParser->get(argParser->OPTION_MODEL); appData._modelParameters = argParser->get >(argParser->OPTION_MODEL_PARAMETERS); } } } } diff --git a/testing/core/DVHCalculatorTest.cpp b/testing/core/DVHCalculatorTest.cpp index 8e8cc8d..23f51ef 100644 --- a/testing/core/DVHCalculatorTest.cpp +++ b/testing/core/DVHCalculatorTest.cpp @@ -1,128 +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) */ // this file defines the rttbCoreTests for the test driver // and all it expects is that you have a function called RegisterTests #include #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVHCalculator.h" #include "rttbGenericMaskedDoseIterator.h" #include "rttbGenericDoseIterator.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "DummyDoseAccessor.h" #include "DummyMaskAccessor.h" namespace rttb { namespace testing { typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; typedef core::GenericMaskedDoseIterator::MaskAccessorPointer MaskAccessorPointer; typedef core::DVHCalculator::DoseIteratorPointer DoseIteratorPointer; typedef core::DVHCalculator::MaskedDoseIteratorPointer MaskedDoseIteratorPointer; /*!@brief DVHTest - test the API of DVH 1) test constructors (values as expected?) */ int DVHCalculatorTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; //create null pointer to DoseIterator DoseIteratorPointer spDoseIteratorNull; const IDType structureID = "myStructure"; const IDType doseIDNull = "myDose"; //1) test constructors (values as expected?) CHECK_THROW_EXPLICIT(core::DVHCalculator myDVHCalc(spDoseIteratorNull, structureID, doseIDNull), core::NullPointerException); //create dummy DoseAccessor boost::shared_ptr spTestDoseAccessor = boost::make_shared(); DoseAccessorPointer spDoseAccessor(spTestDoseAccessor); const std::vector* doseVals = spTestDoseAccessor->getDoseVector(); //create corresponding DoseIterator boost::shared_ptr spTestDoseIterator = boost::make_shared(spDoseAccessor); DoseIteratorPointer spDoseIterator(spTestDoseIterator); const IDType doseID = spDoseAccessor->getUID(); CHECK_NO_THROW(core::DVHCalculator myDVHCalc(spDoseIterator, structureID, doseID)); CHECK_THROW_EXPLICIT(core::DVHCalculator myDVHCalc(spDoseIterator, structureID, doseID, -1), core::InvalidParameterException); int numBins = 21; DoseTypeGy binSize = 0; DoseStatisticType maximum = 0; std::vector::const_iterator doseIt = doseVals->begin(); while (doseIt != doseVals->end()) { if (maximum < *doseIt) { maximum = *doseIt; } - doseIt++; + ++doseIt; } binSize = maximum * 1.5 / numBins; CHECK_NO_THROW(core::DVHCalculator myDVHCalc(spDoseIterator, structureID, doseID, binSize, numBins)); CHECK_THROW_EXPLICIT(core::DVHCalculator myDVHCalc(spDoseIterator, structureID, doseID, 0, 0), core::InvalidParameterException);//aNumberOfBins must be >=0 core::DVHCalculator myDVHCalc(spDoseIterator, structureID, doseID, binSize, 1); CHECK_THROW_EXPLICIT(myDVHCalc.generateDVH(), core::InvalidParameterException);//_numberOfBins must be > max(aDoseIterator)/aDeltaD! //create dummy MaskAccessor boost::shared_ptr spTestMaskAccessor = boost::make_shared(spDoseAccessor->getGeometricInfo()); MaskAccessorPointer spMaskAccessor(spTestMaskAccessor); //create corresponding MaskedDoseIterator boost::shared_ptr spTestMaskedDoseIterator = boost::make_shared(spMaskAccessor, spDoseAccessor); DoseIteratorPointer spMaskedDoseIterator(spTestMaskedDoseIterator); CHECK_NO_THROW(core::DVHCalculator myDVHCalc(spMaskedDoseIterator, structureID, doseID)); //actual calculation is still missing RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/core/DVHTest.cpp b/testing/core/DVHTest.cpp index a338d42..1349013 100644 --- a/testing/core/DVHTest.cpp +++ b/testing/core/DVHTest.cpp @@ -1,191 +1,191 @@ // ----------------------------------------------------------------------- // 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 #include #include "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbDVH.h" #include "DummyDoseAccessor.h" #include "DummyMaskAccessor.h" namespace rttb { namespace testing { typedef core::DVH::DataDifferentialType DataDifferentialType; /*! @brief DVHTest - test the API of DVH 1) test constructors (values as expected?) 2) test asignement 3) test set/getLabel 4) test set/get 5) test equality */ int DVHTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; //generate artificial DVH and corresponding statistical values DoseTypeGy binSize = DoseTypeGy(0.1); DoseVoxelVolumeType voxelVolume = 8; DataDifferentialType anEmptyDataDifferential; DataDifferentialType aDataDifferential; DataDifferentialType aDataDifferential2; DataDifferentialType aDataDifferentialRelative; DoseStatisticType maximum = 0; DoseStatisticType minimum = 0; double sum = 0; double squareSum = 0; DoseCalcType value = 0; DVHVoxelNumber numberOfVoxels = 0; // creat default values [0,100) for (int i = 0; i < 100; i++) { value = DoseCalcType((double(rand()) / RAND_MAX) * 1000); numberOfVoxels += value; aDataDifferential.push_back(value); aDataDifferential2.push_back(value * 10); if (value > 0) { maximum = (i + 0.5) * binSize; if (minimum == 0) { minimum = (i + 0.5) * binSize; } } sum += value * (i + 0.5) * binSize; squareSum += value * pow((i + 0.5) * binSize, 2); } DoseStatisticType mean = sum / numberOfVoxels; DoseStatisticType variance = (squareSum / numberOfVoxels - mean * mean); DoseStatisticType stdDeviation = pow(variance, 0.5); std::deque::iterator it; - for (it = aDataDifferential.begin(); it != aDataDifferential.end(); it++) + for (it = aDataDifferential.begin(); it != aDataDifferential.end(); ++it) { aDataDifferentialRelative.push_back((*it) / numberOfVoxels); } const IDType structureID = "myStructure"; const IDType doseID = "myDose"; const IDType voxelizationID = "myVoxelization"; //1) test default constructor (values as expected?) CHECK_THROW(core::DVH(anEmptyDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(anEmptyDataDifferential, binSize, voxelVolume, structureID, doseID)); CHECK_NO_THROW(core::DVH(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID)); CHECK_NO_THROW(core::DVH(aDataDifferential, binSize, voxelVolume, structureID, doseID)); CHECK_THROW(core::DVH(aDataDifferential, 0, voxelVolume, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(aDataDifferential, 0, voxelVolume, structureID, doseID)); CHECK_THROW(core::DVH(aDataDifferential, binSize, 0, structureID, doseID, voxelizationID)); CHECK_THROW(core::DVH(aDataDifferential, binSize, 0, structureID, doseID)); //2) test asignement core::DVH myTestDVH(aDataDifferential, binSize, voxelVolume, structureID, doseID); CHECK_NO_THROW(core::DVH myOtherDVH = myTestDVH); const core::DVH myOtherDVH = myTestDVH; CHECK_NO_THROW(core::DVH aDVH(myOtherDVH)); //3) test set/getLabel core::DVH myDVH(aDataDifferential, binSize, voxelVolume, structureID, doseID, voxelizationID); StructureLabel label = ""; CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); label = "myLabel"; CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); label = "myLabel2"; CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_NO_THROW(myDVH.setLabel(label)); CHECK_EQUAL(myDVH.getLabel(), label); CHECK_EQUAL(myDVH.getLabel(), label); //4) test set/get //IDs CHECK_EQUAL(myDVH.getStructureID(), structureID); CHECK_EQUAL(myDVH.getDoseID(), doseID); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); /*! is related to #2029*/ myDVH.setDoseID("somethingElse"); CHECK_EQUAL(myDVH.getDoseID(), "somethingElse"); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); CHECK_EQUAL(myDVH.getStructureID(), structureID); /*! is related to #2029*/ myDVH.setStructureID("somethingOther"); CHECK_EQUAL(myDVH.getDoseID(), "somethingElse"); CHECK_EQUAL(myDVH.getVoxelizationID(), voxelizationID); CHECK_EQUAL(myDVH.getStructureID(), "somethingOther"); //dataDifferential CHECK(myDVH.getDataDifferential() == aDataDifferential); CHECK(myDVH.getDataDifferential(false) == aDataDifferential); CHECK(myDVH.getDataDifferential(true) == aDataDifferentialRelative); CHECK_EQUAL(myDVH.getNumberOfVoxels(), numberOfVoxels); CHECK_EQUAL(myDVH.getDeltaV(), voxelVolume); CHECK_EQUAL(myDVH.getDeltaD(), binSize); CHECK_EQUAL(myDVH.getMaximum(), maximum); CHECK_EQUAL(myDVH.getMinimum(), minimum); CHECK_EQUAL(myDVH.getMean(), mean); CHECK_EQUAL(myDVH.getVariance(), variance); CHECK_EQUAL(myDVH.getStdDeviation(), stdDeviation); int percentage = 20; VolumeType absVol = VolumeType(percentage * numberOfVoxels * voxelVolume / 100.0); CHECK_EQUAL(myDVH.getAbsoluteVolume(percentage), absVol); //5) test equality core::DVH myDVH2(aDataDifferential2, binSize, voxelVolume, structureID, doseID); CHECK(!(myDVH == myDVH2)); CHECK_EQUAL(myDVH, myDVH); core::DVH aDVH(myOtherDVH); CHECK_EQUAL(aDVH, myOtherDVH); RETURN_AND_REPORT_TEST_SUCCESS; } }//end namespace testing }//end namespace rttb diff --git a/testing/core/StructureTest.cpp b/testing/core/StructureTest.cpp index 42b550a..f6e7cf6 100644 --- a/testing/core/StructureTest.cpp +++ b/testing/core/StructureTest.cpp @@ -1,138 +1,138 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBaseType.h" #include "rttbStructure.h" #include "rttbInvalidParameterException.h" #include "DummyStructure.h" #include "DummyDoseAccessor.h" namespace rttb { namespace testing { /*! @brief StructureTest - tests the API for Structure 1) constructors 2) get/setXX */ int StructureTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; boost::shared_ptr spTestDoseAccessor = boost::make_shared(); DummyStructure myStructGenerator(spTestDoseAccessor->getGeometricInfo()); //1) constructors CHECK_NO_THROW(core::Structure()); core::Structure emptyTestStruct; CHECK_EQUAL("None", emptyTestStruct.getLabel()); CHECK_NO_THROW(emptyTestStruct.getUID()); GridIndexType zPlane = 4; core::Structure rect = myStructGenerator.CreateRectangularStructureCentered(zPlane); CHECK_NO_THROW(core::Structure(rect.getStructureVector())); core::Structure rect2 = core::Structure(rect.getStructureVector()); CHECK_EQUAL(rect.getLabel(), rect2.getLabel()); CHECK(rect.getUID() != rect2.getUID()); PolygonSequenceType rectVec = rect.getStructureVector(); PolygonSequenceType rect2Vec = rect2.getStructureVector(); CHECK_EQUAL(rectVec.size(), rect2Vec.size()); PolygonSequenceType::iterator it = rectVec.begin(); PolygonSequenceType::iterator it2 = rect2Vec.begin(); - for (; it != rectVec.end(); it++) + for (; it != rectVec.end(); ++it) { CHECK_EQUAL(it->size(), it2->size()); PolygonType::iterator pit = it->begin(); PolygonType::iterator pit2 = it2->begin(); - for (; pit != it->end(); pit++) + for (; pit != it->end(); ++pit) { CHECK_EQUAL(*(pit), *(pit2)); - pit2++; + ++pit2; } - it2++; + ++it2; } CHECK_NO_THROW(core::Structure rect3 = rect); core::Structure rect3 = rect; CHECK_EQUAL(rect.getLabel(), rect3.getLabel()); CHECK_EQUAL(rect.getUID(), rect3.getUID()); PolygonSequenceType rect3Vec = rect3.getStructureVector(); CHECK_EQUAL(rectVec.size(), rect3Vec.size()); it = rectVec.begin(); PolygonSequenceType::iterator it3 = rect3Vec.begin(); - for (; it != rectVec.end(); it++) + for (; it != rectVec.end(); ++it) { CHECK_EQUAL(it->size(), it3->size()); PolygonType::iterator pit = it->begin(); PolygonType::iterator pit3 = it3->begin(); - for (; pit != it->end(); pit++) + for (; pit != it->end(); ++pit) { CHECK_EQUAL(*(pit), *(pit3)); - pit3++; + ++pit3; } - it3++; + ++it3; } //2) get/setXX CHECK_EQUAL("None", emptyTestStruct.getLabel()); CHECK_NO_THROW(emptyTestStruct.setLabel("NEW Label")); CHECK_EQUAL("NEW Label", emptyTestStruct.getLabel()); CHECK_NO_THROW(emptyTestStruct.getUID()); CHECK_NO_THROW(emptyTestStruct.setUID("1.2.345.67.8.9")); CHECK_EQUAL("1.2.345.67.8.9", emptyTestStruct.getUID()); CHECK((emptyTestStruct.getStructureVector()).empty()); CHECK_EQUAL(0, emptyTestStruct.getNumberOfEndpoints()); CHECK_EQUAL(4, rect.getNumberOfEndpoints()); CHECK_EQUAL(rect.getNumberOfEndpoints(), rect2.getNumberOfEndpoints()); core::Structure circ = myStructGenerator.CreateTestStructureCircle(zPlane); CHECK_EQUAL(4004, circ.getNumberOfEndpoints()); core::Structure multiPoly = myStructGenerator.CreateTestStructureIntersectingTwoPolygonsInDifferentSlices(zPlane, zPlane + 1); CHECK_EQUAL(8, multiPoly.getNumberOfEndpoints()); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb \ No newline at end of file diff --git a/testing/demoapps/DoseTool/DoseToolDVHTest.cpp b/testing/demoapps/DoseTool/DoseToolDVHTest.cpp index b5e9425..6d71d25 100644 --- a/testing/demoapps/DoseTool/DoseToolDVHTest.cpp +++ b/testing/demoapps/DoseTool/DoseToolDVHTest.cpp @@ -1,107 +1,106 @@ // ----------------------------------------------------------------------- // 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 #include #include "litCheckMacros.h" #include "boost/filesystem.hpp" namespace rttb { namespace testing { //path to the current running directory. DoseTool is in the same directory (Debug/Release) extern const char* _callingAppPath; static std::string readFile(const std::string& filename); int DoseToolDVHTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string doseToolExecutable; std::string doseFilename; - std::string planFilename; std::string structFilename; std::string structName; std::string referenceXMLFilename; boost::filesystem::path callingPath(_callingAppPath); if (argc > 5) { doseToolExecutable = argv[1]; doseFilename = argv[2]; structFilename = argv[3]; structName = argv[4]; referenceXMLFilename = argv[5]; } std::string voxelizerToolExeWithPath = callingPath.parent_path().string() + "/" + doseToolExecutable; std::string defaultOutputFilename = "dicomDVHOutput.xml"; std::string baseCommand = voxelizerToolExeWithPath; baseCommand += " -d " + doseFilename; baseCommand += " -s " + structFilename; if (structName != "") { baseCommand += " -n " + structName; } else { baseCommand += " -u itk "; } std::string defaultDVHCommand = baseCommand + " -z " + defaultOutputFilename; std::cout << "Command line call: " + defaultDVHCommand << std::endl; CHECK_EQUAL(system(defaultDVHCommand.c_str()), 0); //check if file exists CHECK_EQUAL(boost::filesystem::exists(defaultOutputFilename), true); //check if file is the same than reference file std::string defaultAsIs = readFile(defaultOutputFilename); std::string defaultExpected = readFile(referenceXMLFilename); CHECK_EQUAL(defaultAsIs, defaultExpected); //delete file again CHECK_EQUAL(std::remove(defaultOutputFilename.c_str()), 0); RETURN_AND_REPORT_TEST_SUCCESS; } std::string readFile(const std::string& filename) { std::ifstream fileStream(filename.c_str()); std::string content((std::istreambuf_iterator(fileStream)), (std::istreambuf_iterator())); return content; } } //namespace testing } //namespace rttb diff --git a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerAllStructsTest.cpp b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerAllStructsTest.cpp index 6ce2031..324a947 100644 --- a/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerAllStructsTest.cpp +++ b/testing/demoapps/VoxelizerTool/rttbVoxelizerToolVoxelizerAllStructsTest.cpp @@ -1,111 +1,109 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include #include /*! @brief VoxelizerToolTest. Tests a selection of structs. */ namespace rttb { namespace testing { //path to the current running directory. VoxelizerTool is in the same directory (Debug/Release) extern const char* _callingAppPath; int VoxelizerToolVoxelizerAllStructsTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; std::string voxelizerToolExe; std::string tempDirectory; std::string structFile; std::string referenceFile; if (argc > 4) { voxelizerToolExe = argv[1]; tempDirectory = argv[2]; structFile = argv[3]; referenceFile = argv[4]; } boost::filesystem::path callingPath(_callingAppPath); std::string voxelizerToolExeWithPath = callingPath.parent_path().string() + "/" + voxelizerToolExe; std::vector structNames; structNames.push_back("Niere.*"); structNames.push_back("Magen/DD"); structNames.push_back("PTV"); //structure names will be used for file naming, BUT '.' in the end will be cropped and '/' will be replaced by '_'. Thus, the different filenames. std::vector filenames; filenames.push_back("Niere re"); filenames.push_back("Magen_DD"); filenames.push_back("PTV"); std::string baseCommand = voxelizerToolExeWithPath; baseCommand += " -s " + structFile; baseCommand += " -r " + referenceFile; baseCommand += " -e \""; for (int i = 0; i < structNames.size(); i++) { std::string command = baseCommand + structNames.at(i) + "\""; std::cout << "Command line call: " + command << std::endl; int returnValue = system(command.c_str()); CHECK_EQUAL(returnValue, 0); - const std::string HDRfileName = tempDirectory + "/out_" + filenames.at(i) + ".hdr"; boost::filesystem::path HDRFile(tempDirectory); HDRFile /= "out_" + filenames.at(i) + ".hdr"; - const std::string IMGfileName = tempDirectory + "/out_" + filenames.at(i) + ".img"; boost::filesystem::path IMGFile(tempDirectory); IMGFile /= "out_" + filenames.at(i) + ".img"; CHECK_EQUAL( boost::filesystem::exists(HDRFile), true); CHECK_EQUAL( boost::filesystem::exists(IMGFile), true); if (boost::filesystem::exists(IMGFile)) { boost::filesystem::remove(IMGFile); } if (boost::filesystem::exists(HDRFile)) { boost::filesystem::remove(HDRFile); } } RETURN_AND_REPORT_TEST_SUCCESS; } } } diff --git a/testing/examples/DoseStatisticsIOVirtuosTest.cpp b/testing/examples/DoseStatisticsIOVirtuosTest.cpp deleted file mode 100644 index fe450bf..0000000 --- a/testing/examples/DoseStatisticsIOVirtuosTest.cpp +++ /dev/null @@ -1,115 +0,0 @@ -// ----------------------------------------------------------------------- -// RTToolbox - DKFZ radiotherapy quantitative evaluation library -// -// Copyright (c) German Cancer Research Center (DKFZ), -// Software development for Integrated Diagnostics and Therapy (SIDT). -// ALL RIGHTS RESERVED. -// See rttbCopyright.txt or -// http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] -// -// This software is distributed WITHOUT ANY WARRANTY; without even -// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -// PURPOSE. See the above copyright notices for more information. -// -//------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (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 "litCheckMacros.h" - -#include "rttbBaseType.h" -#include "rttbDoseStatistics.h" -#include "rttbDoseStatisticsCalculator.h" -#include "rttbDoseStatisticsXMLWriter.h" -#include "rttbGenericMaskedDoseIterator.h" -#include "rttbDoseIteratorInterface.h" -#include "rttbVirtuosPlanFileDoseAccessorGenerator.h" -#include "rttbVirtuosDoseAccessor.h" -#include "rttbVirtuosFileStructureSetGenerator.h" -#include "rttbBoostMaskAccessor.h" - - -namespace rttb -{ - - namespace testing - { - - /*! @brief OtherIOTest - test the IO for dose statistics - 1) test exception - 2) test writing statistics to xml file - */ - - int DoseStatisticsIOVirtuosTest(int argc, char* argv[]) - { - typedef core::GenericDoseIterator::DoseAccessorPointer DoseAccessorPointer; - typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; - typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; - typedef boost::shared_ptr DoseStatisticsPtr; - typedef core::MaskAccessorInterface::MaskAccessorPointer MaskAccessorPointer; - typedef boost::shared_ptr GeometricInfoPointer; - - PREPARE_DEFAULT_TEST_REPORTING; - - std::string Struct_FILENAME; - std::string BPLCT_FILENAME; - std::string Dose_FILENAME; - std::string Plan_FILENAME; - std::string Struct_NAME; - - if (argc > 5) - { - Struct_FILENAME = argv[1]; - BPLCT_FILENAME = argv[2]; - Dose_FILENAME = argv[3]; - Plan_FILENAME = argv[4]; - Struct_NAME = argv[5]; - } - else - { - std::cout << "at least 5 arguments are expected in DoseStatisticsIOVirtuosTest." << std::endl; - } - - /* generate dummy dose */ - io::virtuos::VirtuosPlanFileDoseAccessorGenerator doseAccessorGenerator(Dose_FILENAME.c_str(), - Plan_FILENAME.c_str()); - - DoseAccessorPointer doseAccessor(doseAccessorGenerator.generateDoseAccessor()); - - StructureSetPointer rtStructureSet = io::virtuos::VirtuosFileStructureSetGenerator( - Struct_FILENAME.c_str(), BPLCT_FILENAME.c_str()).generateStructureSet(); - - GeometricInfoPointer geometricInfoPtr = boost::make_shared - (doseAccessor->getGeometricInfo()); - MaskAccessorPointer maskAccessorPtr = boost::make_shared - (rtStructureSet->getStructure(3), doseAccessor->getGeometricInfo()); - maskAccessorPtr->updateMask(); - - boost::shared_ptr spTestDoseIterator = - boost::make_shared(maskAccessorPtr, doseAccessor); - DoseIteratorPointer spDoseIterator(spTestDoseIterator); - - rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator); - rttb::algorithms::DoseStatistics::DoseStatisticsPointer doseStatistics = - myDoseStatsCalculator.calculateDoseStatistics(true); - - /* test writing statistcs to xml file */ - FileNameString fN = "doseStatisticsVirtuosComplex.xml"; - CHECK_NO_THROW(io::other::writeDoseStatistics(doseStatistics, fN)); - - RETURN_AND_REPORT_TEST_SUCCESS; - } - - }//testing -}//rttb - diff --git a/testing/examples/RTBioModelExampleTest.cpp b/testing/examples/RTBioModelExampleTest.cpp index c292641..cf86376 100644 --- a/testing/examples/RTBioModelExampleTest.cpp +++ b/testing/examples/RTBioModelExampleTest.cpp @@ -1,416 +1,416 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBioModel.h" #include "rttbDVHTxtFileReader.h" #include "rttbDVH.h" #include "rttbTCPLQModel.h" #include "rttbNTCPLKBModel.h" #include "rttbNTCPRSModel.h" #include "rttbBioModelScatterPlots.h" #include "rttbBioModelCurve.h" #include "rttbDvhBasedModels.h" #include "rttbDoseIteratorInterface.h" namespace rttb { namespace testing { /*! @brief RTBioModelTest. TCP calculated using a DVH PTV and LQ Model. NTCP tested using 3 Normal Tissue DVHs and LKB/RS Model. 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! */ int RTBioModelExampleTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; typedef rttb::models::CurveDataType CurveDataType; typedef std::multimap > ScatterPlotType; typedef core::DVH::DVHPointer DVHPointer; //increased accuracy requires double values in the calculation (rttbBaseType.h) double toleranceEUD = 1e-5; double tolerance = 1e-7; //ARGUMENTS: 1: ptv dvh file name // 2: normal tissue 1 dvh file name // 3: normal tissue 2 dvh file name // 4: normal tissue 3 dvh file name //...........5: Virtuos MPM_LR_ah dvh lung file name //...........6: Virtuos MPM_LR_ah dvh target file name std::string DVH_FILENAME_PTV; std::string DVH_FILENAME_NT1; std::string DVH_FILENAME_NT2; std::string DVH_FILENAME_NT3; std::string DVH_FILENAME_TV_TEST; std::string DVH_Virtuos_Target; std::string DVH_Virtuos_Lung; if (argc > 1) { DVH_FILENAME_PTV = argv[1]; } if (argc > 2) { DVH_FILENAME_NT1 = argv[2]; } if (argc > 3) { DVH_FILENAME_NT2 = argv[3]; } if (argc > 4) { DVH_FILENAME_NT3 = argv[4]; } if (argc > 5) { DVH_FILENAME_TV_TEST = argv[5]; } if (argc > 6) { DVH_Virtuos_Lung = argv[6]; } if (argc > 7) { DVH_Virtuos_Target = argv[7]; } //DVH PTV rttb::io::other::DVHTxtFileReader dvhReader = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_PTV); DVHPointer dvhPtr = dvhReader.generateDVH(); CHECK_CLOSE(6.04759613161786830000e+001, models::getEUD(dvhPtr, 10), toleranceEUD); rttb::io::other::DVHTxtFileReader dvhReader_test_tv = rttb::io::other::DVHTxtFileReader( DVH_FILENAME_TV_TEST); DVHPointer dvh_test_tv = dvhReader_test_tv.generateDVH(); //test TCP LQ Model models::BioModelParamType alpha = 0.35; models::BioModelParamType beta = 0.023333333333333; models::BioModelParamType roh = 10000000; int numFractions = 2; DoseTypeGy normalizationDose = 68; rttb::models::TCPLQModel tcplq = rttb::models::TCPLQModel(dvhPtr, alpha, beta, roh, numFractions); CHECK_EQUAL(alpha, tcplq.getAlphaMean()); CHECK_EQUAL(alpha / beta, tcplq.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); CHECK_NO_THROW(tcplq.init()); if (tcplq.init()) { CHECK_CLOSE(1.00497232941856940000e-127, tcplq.getValue(), tolerance); } CurveDataType curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); CurveDataType::iterator it; - for (it = curve.begin(); it != curve.end(); it++) + for (it = curve.begin(); it != curve.end(); ++it) { if ((*it).first < 72) { CHECK_EQUAL(0, (*it).second); } else if ((*it).first > 150) { CHECK((*it).second > 0.9); } } models::BioModelParamType alphaBeta = 10; tcplq.setParameters(alpha, alphaBeta, roh, 0.08); CHECK_EQUAL(alpha, tcplq.getAlphaMean()); CHECK_EQUAL(alphaBeta, tcplq.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); if (tcplq.init()) { CHECK_CLOSE(1.84e-005, tcplq.getValue(), tolerance); } normalizationDose = 40; curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); alpha = 1; alphaBeta = 14.5; tcplq.setAlpha(alpha); tcplq.setAlphaBeta(alphaBeta); tcplq.setRho(roh); CHECK_EQUAL(alpha, tcplq.getAlphaMean()); CHECK_EQUAL(alphaBeta, tcplq.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); if (tcplq.init()) { CHECK_CLOSE(0.954885, tcplq.getValue(), toleranceEUD); } alpha = 0.9; alphaBeta = 1; tcplq.setAlpha(alpha); tcplq.setAlphaBeta(alphaBeta); tcplq.setRho(roh); CHECK_EQUAL(alpha, tcplq.getAlphaMean()); CHECK_EQUAL(alphaBeta, tcplq.getAlphaBeta()); CHECK_EQUAL(roh, tcplq.getRho()); if (tcplq.init()) { CHECK_EQUAL(1, tcplq.getValue()); } //TCP LQ Test alpha = 0.3; beta = 0.03; roh = 10000000; numFractions = 20; rttb::models::TCPLQModel tcplq_test = rttb::models::TCPLQModel(dvh_test_tv, alpha, beta, roh, numFractions); CHECK_EQUAL(alpha, tcplq_test.getAlphaMean()); CHECK_EQUAL(alpha / beta, tcplq_test.getAlphaBeta()); CHECK_EQUAL(roh, tcplq_test.getRho()); CHECK_NO_THROW(tcplq_test.init()); if (tcplq_test.init()) { CHECK_CLOSE(9.79050278878883180000e-001, tcplq_test.getValue(), tolerance); } normalizationDose = 60; curve = models::getCurveDoseVSBioModel(tcplq_test, normalizationDose); //DVH HT 1 rttb::io::other::DVHTxtFileReader dvhReader2 = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_NT1); DVHPointer dvhPtr2 = dvhReader2.generateDVH(); CHECK_CLOSE(1.07920836034015810000e+001, models::getEUD(dvhPtr2, 10), toleranceEUD); //test RTNTCPLKBModel rttb::models::NTCPLKBModel lkb = rttb::models::NTCPLKBModel(); models::BioModelParamType aVal = 10; models::BioModelParamType mVal = 0.16; models::BioModelParamType d50Val = 55; CHECK_EQUAL(0, lkb.getA()); CHECK_EQUAL(0, lkb.getM()); CHECK_EQUAL(0, lkb.getD50()); lkb.setDVH(dvhPtr2); CHECK_EQUAL(dvhPtr2, lkb.getDVH()); lkb.setA(aVal); CHECK_EQUAL(aVal, lkb.getA()); lkb.setM(mVal); CHECK_EQUAL(mVal, lkb.getM()); lkb.setD50(d50Val); CHECK_EQUAL(d50Val, lkb.getD50()); CHECK_NO_THROW(lkb.init()); if (lkb.init()) { CHECK_CLOSE(2.53523522831366570000e-007, lkb.getValue(), tolerance); } //test RTNTCPRSModel rttb::models::NTCPRSModel rs = rttb::models::NTCPRSModel(); models::BioModelParamType gammaVal = 1.7; models::BioModelParamType sVal = 1; CHECK_EQUAL(0, rs.getGamma()); CHECK_EQUAL(0, rs.getS()); CHECK_EQUAL(0, rs.getD50()); rs.setDVH(dvhPtr2); CHECK_EQUAL(dvhPtr2, rs.getDVH()); rs.setD50(d50Val); CHECK_EQUAL(d50Val, rs.getD50()); rs.setGamma(gammaVal); CHECK_EQUAL(gammaVal, rs.getGamma()); rs.setS(sVal); CHECK_EQUAL(sVal, rs.getS()); CHECK_NO_THROW(rs.init()); if (rs.init()) { CHECK_CLOSE(3.70385888626145740000e-009, rs.getValue(), tolerance); } //DVH HT 2 rttb::io::other::DVHTxtFileReader dvhReader3 = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_NT2); DVHPointer dvhPtr3 = dvhReader3.generateDVH(); CHECK_CLOSE(1.26287047025885110000e+001, models::getEUD(dvhPtr3, 10), toleranceEUD); //test RTNTCPLKBModel aVal = 10; mVal = 0.16; d50Val = 55; lkb.setDVH(dvhPtr3); CHECK_EQUAL(dvhPtr3, lkb.getDVH()); lkb.setA(aVal); CHECK_EQUAL(aVal, lkb.getA()); lkb.setM(mVal); CHECK_EQUAL(mVal, lkb.getM()); lkb.setD50(d50Val); CHECK_EQUAL(d50Val, lkb.getD50()); if (lkb.init()) { CHECK_CLOSE(7.36294657754956700000e-007, lkb.getValue(), tolerance); } //test RTNTCPRSModel rs = rttb::models::NTCPRSModel(); gammaVal = 1.7; sVal = 1; CHECK_EQUAL(0, rs.getGamma()); CHECK_EQUAL(0, rs.getS()); CHECK_EQUAL(0, rs.getD50()); rs.setDVH(dvhPtr3); CHECK_EQUAL(dvhPtr3, rs.getDVH()); rs.setD50(d50Val); CHECK_EQUAL(d50Val, rs.getD50()); rs.setGamma(gammaVal); CHECK_EQUAL(gammaVal, rs.getGamma()); rs.setS(sVal); CHECK_EQUAL(sVal, rs.getS()); if (rs.init()) { CHECK_CLOSE(1.76778795490939440000e-007, rs.getValue(), tolerance); } //DVH HT 3 rttb::io::other::DVHTxtFileReader dvhReader4 = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_NT3); DVHPointer dvhPtr4 = dvhReader4.generateDVH(); CHECK_CLOSE(2.18212982041056310000e+001, models::getEUD(dvhPtr4, 10), toleranceEUD); //test RTNTCPLKBModel aVal = 10; mVal = 0.16; d50Val = 55; lkb.setDVH(dvhPtr4); CHECK_EQUAL(dvhPtr4, lkb.getDVH()); lkb.setA(aVal); CHECK_EQUAL(aVal, lkb.getA()); lkb.setM(mVal); CHECK_EQUAL(mVal, lkb.getM()); lkb.setD50(d50Val); CHECK_EQUAL(d50Val, lkb.getD50()); if (lkb.init()) { CHECK_CLOSE(8.15234192641929420000e-005, lkb.getValue(), tolerance); } //test RTNTCPRSModel rs = rttb::models::NTCPRSModel(); gammaVal = 1.7; sVal = 1; CHECK_EQUAL(0, rs.getGamma()); CHECK_EQUAL(0, rs.getS()); CHECK_EQUAL(0, rs.getD50()); rs.setDVH(dvhPtr4); CHECK_EQUAL(dvhPtr4, rs.getDVH()); rs.setD50(d50Val); CHECK_EQUAL(d50Val, rs.getD50()); rs.setGamma(gammaVal); CHECK_EQUAL(gammaVal, rs.getGamma()); rs.setS(sVal); CHECK_EQUAL(sVal, rs.getS()); if (rs.init()) { CHECK_CLOSE(2.02607985020919480000e-004, rs.getValue(), tolerance); } //test using Virtuos Pleuramesotheliom MPM_LR_ah //DVH PTV rttb::io::other::DVHTxtFileReader dR_Target = rttb::io::other::DVHTxtFileReader(DVH_Virtuos_Target); DVHPointer dvhPtrTarget = dR_Target.generateDVH(); rttb::io::other::DVHTxtFileReader dR_Lung = rttb::io::other::DVHTxtFileReader(DVH_Virtuos_Lung); DVHPointer dvhPtrLung = dR_Lung.generateDVH(); //test TCP LQ Model models::BioModelParamType alphaMean = 0.34; models::BioModelParamType alphaVarianz = 0.02; models::BioModelParamType alpha_beta = 28; models::BioModelParamType rho = 1200; int numFractionsVirtuos = 27; rttb::models::TCPLQModel tcplqVirtuos = rttb::models::TCPLQModel(dvhPtrTarget, rho, numFractionsVirtuos, alpha_beta, alphaMean, alphaVarianz); if (tcplqVirtuos.init()) { CHECK_CLOSE(0.8894, tcplqVirtuos.getValue(), 1e-4); } models::BioModelParamType d50Mean = 20; models::BioModelParamType d50Varianz = 2; models::BioModelParamType m = 0.36; models::BioModelParamType a = 1.06; rttb::models::NTCPLKBModel lkbVirtuos = rttb::models::NTCPLKBModel(dvhPtrLung, d50Mean, m, a); if (lkbVirtuos.init()) { CHECK_CLOSE(0.0397, lkbVirtuos.getValue(), 1e-4); } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb \ No newline at end of file diff --git a/testing/examples/RTBioModelScatterPlotExampleTest.cpp b/testing/examples/RTBioModelScatterPlotExampleTest.cpp index 4cefe2f..3f8b8fe 100644 --- a/testing/examples/RTBioModelScatterPlotExampleTest.cpp +++ b/testing/examples/RTBioModelScatterPlotExampleTest.cpp @@ -1,480 +1,479 @@ // ----------------------------------------------------------------------- // 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 "litCheckMacros.h" #include "rttbBioModel.h" #include "rttbDVHTxtFileReader.h" #include "rttbDVH.h" #include "rttbTCPLQModel.h" #include "rttbNTCPLKBModel.h" #include "rttbNTCPRSModel.h" #include "rttbBioModelScatterPlots.h" #include "rttbBioModelCurve.h" #include "rttbDvhBasedModels.h" #include "../models/rttbScatterTester.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace testing { /*! @brief RTBioModelScatterPlotExampleTest. calculating Curves and Scatterplots for TCP and NTCP models. The values on curve and scatterplot need to be similar for similar dose values. The range of difference is given by the variance used to generate the scatter. WARNING: The values for comparison need to be adjusted if the input files are changed! */ int RTBioModelScatterPlotExampleTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; typedef rttb::models::CurveDataType CurveDataType; typedef rttb::models::ScatterPlotType ScatterPlotType; typedef core::DVH::DVHPointer DVHPointer; //increased accuracy requires double values in the calculation (rttbBaseType.h) double toleranceEUD = 1e-5; - double tolerance = 1e-7; //ARGUMENTS: 1: ptv dvh file name // 2: normal tissue 1 dvh file name // 3: TV dvh file name std::string DVH_FILENAME_PTV; std::string DVH_FILENAME_NT1; std::string DVH_FILENAME_TV_TEST; if (argc > 1) { DVH_FILENAME_PTV = argv[1]; } if (argc > 2) { DVH_FILENAME_NT1 = argv[2]; } if (argc > 3) { DVH_FILENAME_TV_TEST = argv[3]; } //DVH PTV rttb::io::other::DVHTxtFileReader dvhReader = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_PTV); DVHPointer dvhPtr = dvhReader.generateDVH(); rttb::io::other::DVHTxtFileReader dvhReader_test_tv = rttb::io::other::DVHTxtFileReader( DVH_FILENAME_TV_TEST); DVHPointer dvh_test_tv = dvhReader_test_tv.generateDVH(); //test TCP LQ Model models::BioModelParamType alpha = 0.35; models::BioModelParamType beta = 0.023333333333333; models::BioModelParamType roh = 10000000; int numFractions = 2; DoseTypeGy normalizationDose = 68; rttb::models::TCPLQModel tcplq = rttb::models::TCPLQModel(dvhPtr, alpha, beta, roh, numFractions); CHECK_NO_THROW(tcplq.init()); CurveDataType curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq, 0, alpha, 0, normalizationDose, 100)); ScatterPlotType tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 0, alpha, 0, normalizationDose, 100); CHECK_EQUAL(100, tcpScatter.size()); ScatterTester scatterCompare(curve, tcpScatter); CHECK_TESTER(scatterCompare); //test also with other parameter tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 3, roh, 0, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); CHECK_TESTER(scatterCompare); std::vector paramIdVec; models::BioModel::ParamVectorType meanVec; models::BioModel::ParamVectorType meanVecTest; meanVecTest.push_back(alpha); models::BioModel::ParamVectorType varianceVec; //"alphaMean":0,"alphaVariance":1,"alpha_beta":2, "rho":3 paramIdVec.push_back(0); meanVec.push_back(tcplq.getAlphaMean()); varianceVec.push_back(0); //setting parameter 1 will change the resulting scatter plot dramatically - is it meant to? //this is unexpected since the value was taken from the original calculation //paramIdVec.push_back(1); meanVec.push_back(tcplq.getAlphaVariance()); varianceVec.push_back(0); paramIdVec.push_back(2); meanVec.push_back(tcplq.getAlphaBeta()); varianceVec.push_back(0); paramIdVec.push_back(3); meanVec.push_back(tcplq.getRho()); varianceVec.push_back(0); CHECK_THROW_EXPLICIT(models::getScatterPlotVaryParameters(tcplq, paramIdVec, meanVecTest, varianceVec, normalizationDose, 50), core::InvalidParameterException); ScatterPlotType scatterVary = models::getScatterPlotVaryParameters(tcplq, paramIdVec, meanVec, varianceVec, normalizationDose, 50); CHECK_EQUAL(50, scatterVary.size()); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); models::BioModelParamType variance = 0.00015; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq, 0, alpha, variance, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 0, alpha, variance, normalizationDose, 100); scatterCompare.setVariance(variance); scatterCompare.setCompareScatter(tcpScatter); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); //test also with other parameter tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 3, roh, variance, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(tcplq.getAlphaMean()); varianceVec.push_back(variance); //paramIdVec.push_back(1); meanVec.push_back(tcplq.getAlphaVariance()); varianceVec.push_back(variance); paramIdVec.push_back(2); meanVec.push_back(tcplq.getAlphaBeta()); varianceVec.push_back(variance); paramIdVec.push_back(3); meanVec.push_back(tcplq.getRho()); varianceVec.push_back(variance); scatterVary = models::getScatterPlotVaryParameters(tcplq, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); models::BioModelParamType alphaBeta = 10; tcplq.setParameters(alpha, alphaBeta, roh, 0.08); tcplq.init(); normalizationDose = 40; curve = models::getCurveDoseVSBioModel(tcplq, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq, 0, alpha, 0, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 0, alpha, 0, normalizationDose, 100); scatterCompare.setReferenceCurve(curve); scatterCompare.setVariance(0); //do not allow larger deviations scatterCompare.setAllowExceptions(false); scatterCompare.setCompareScatter(tcpScatter); CHECK_TESTER(scatterCompare); variance = 0.25; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq, 0, alpha, variance, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq, 0, alpha, variance, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); scatterCompare.setVariance(variance); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); /*TCP LQ Test*/ alpha = 0.3; beta = 0.03; roh = 10000000; numFractions = 20; rttb::models::TCPLQModel tcplq_test = rttb::models::TCPLQModel(dvh_test_tv, alpha, beta, roh, numFractions); CHECK_NO_THROW(tcplq_test.init()); normalizationDose = 60; curve = models::getCurveDoseVSBioModel(tcplq_test, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq_test, 0, alpha, 0, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq_test, 0, alpha, 0, normalizationDose, 100); scatterCompare.setReferenceCurve(curve); scatterCompare.setVariance(0); scatterCompare.setCompareScatter(tcpScatter); CHECK_TESTER(scatterCompare); //test also with other parameter tcpScatter = models::getScatterPlotVary1Parameter(tcplq_test, 3, roh, 0, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(tcplq_test.getAlphaMean()); varianceVec.push_back(0); //paramIdVec.push_back(1); meanVec.push_back(tcplq_test.getAlphaVariance()); varianceVec.push_back(0); paramIdVec.push_back(2); meanVec.push_back(tcplq_test.getAlphaBeta()); varianceVec.push_back(0); paramIdVec.push_back(3); meanVec.push_back(tcplq_test.getRho()); varianceVec.push_back(0); scatterVary = models::getScatterPlotVaryParameters(tcplq_test, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); variance = 0.00025; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(tcplq_test, 0, alpha, variance, normalizationDose, 100)); tcpScatter = models::getScatterPlotVary1Parameter(tcplq_test, 0, alpha, variance, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); scatterCompare.setVariance(variance); CHECK_TESTER(scatterCompare); //test also with other parameter tcpScatter = models::getScatterPlotVary1Parameter(tcplq_test, 3, roh, variance, normalizationDose, 100); scatterCompare.setCompareScatter(tcpScatter); scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); scatterCompare.setAllowExceptions(false); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(tcplq_test.getAlphaMean()); varianceVec.push_back(variance); //paramIdVec.push_back(1); meanVec.push_back(tcplq_test.getAlphaVariance()); varianceVec.push_back(variance); paramIdVec.push_back(2); meanVec.push_back(tcplq_test.getAlphaBeta()); varianceVec.push_back(variance); paramIdVec.push_back(3); meanVec.push_back(tcplq_test.getRho()); varianceVec.push_back(variance); scatterVary = models::getScatterPlotVaryParameters(tcplq_test, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); //allow 5% of the points to deviate more scatterCompare.setAllowExceptions(true); CHECK_TESTER(scatterCompare); //DVH HT 1 rttb::io::other::DVHTxtFileReader dvhReader2 = rttb::io::other::DVHTxtFileReader(DVH_FILENAME_NT1); DVHPointer dvhPtr2 = dvhReader2.generateDVH(); CHECK_CLOSE(1.07920836034015810000e+001, models::getEUD(dvhPtr2, 10), toleranceEUD); //test RTNTCPLKBModel rttb::models::NTCPLKBModel lkb = rttb::models::NTCPLKBModel(); models::BioModelParamType aVal = 10; models::BioModelParamType mVal = 0.16; models::BioModelParamType d50Val = 55; lkb.setDVH(dvhPtr2); lkb.setA(aVal); lkb.setM(mVal); lkb.setD50(d50Val); CHECK_NO_THROW(lkb.init()); normalizationDose = 60; curve = models::getCurveDoseVSBioModel(lkb, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(lkb, 2, aVal, 0, normalizationDose, 100)); ScatterPlotType scatter = models::getScatterPlotVary1Parameter(lkb, 2, aVal, 0, normalizationDose, 100); scatterCompare.setReferenceCurve(curve); scatterCompare.setVariance(0); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); //"d50":0,"m":1,"a":2 //test also with other parameter scatter = models::getScatterPlotVary1Parameter(lkb, 0, d50Val, 0, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(lkb.getD50()); varianceVec.push_back(0); paramIdVec.push_back(1); meanVec.push_back(lkb.getM()); varianceVec.push_back(0); paramIdVec.push_back(2); meanVec.push_back(lkb.getA()); varianceVec.push_back(0); scatterVary = models::getScatterPlotVaryParameters(lkb, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); variance = 0.00025; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(lkb, 2, aVal, variance, normalizationDose, 100)); scatter = models::getScatterPlotVary1Parameter(lkb, 2, aVal, variance, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); scatterCompare.setVariance(variance); CHECK_TESTER(scatterCompare); //test also with other parameter scatter = models::getScatterPlotVary1Parameter(lkb, 0, d50Val, variance, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(lkb.getD50()); varianceVec.push_back(variance); paramIdVec.push_back(1); meanVec.push_back(lkb.getM()); varianceVec.push_back(variance); paramIdVec.push_back(2); meanVec.push_back(lkb.getA()); varianceVec.push_back(variance); scatterVary = models::getScatterPlotVaryParameters(lkb, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); //test RTNTCPRSModel rttb::models::NTCPRSModel rs = rttb::models::NTCPRSModel(); models::BioModelParamType gammaVal = 1.7; models::BioModelParamType sVal = 1; rs.setDVH(dvhPtr2); rs.setD50(d50Val); rs.setGamma(gammaVal); rs.setS(sVal); CHECK_NO_THROW(rs.init()); normalizationDose = 60; curve = models::getCurveDoseVSBioModel(rs, normalizationDose); CHECK_NO_THROW(models::getScatterPlotVary1Parameter(rs, 0, d50Val, 0, normalizationDose, 100)); scatter = models::getScatterPlotVary1Parameter(rs, 0, d50Val, 0, normalizationDose, 100); scatterCompare.setReferenceCurve(curve); scatterCompare.setVariance(0); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); //"d50":0,"gamma":1,"s":2 //test also with other parameter scatter = models::getScatterPlotVary1Parameter(rs, 1, gammaVal, 0, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(rs.getD50()); varianceVec.push_back(0); paramIdVec.push_back(1); meanVec.push_back(rs.getGamma()); varianceVec.push_back(0); paramIdVec.push_back(2); meanVec.push_back(rs.getS()); varianceVec.push_back(0); scatterVary = models::getScatterPlotVaryParameters(rs, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); variance = 0.0075; CHECK_NO_THROW(models::getScatterPlotVary1Parameter(rs, 0, d50Val, variance, normalizationDose, 100)); scatter = models::getScatterPlotVary1Parameter(rs, 0, d50Val, variance, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); scatterCompare.setVariance(variance); CHECK_TESTER(scatterCompare); //test also with other parameter scatter = models::getScatterPlotVary1Parameter(rs, 2, sVal, variance, normalizationDose, 100); scatterCompare.setCompareScatter(scatter); CHECK_TESTER(scatterCompare); paramIdVec.clear(); meanVec.clear(); varianceVec.clear(); paramIdVec.push_back(0); meanVec.push_back(rs.getD50()); varianceVec.push_back(variance); paramIdVec.push_back(1); meanVec.push_back(rs.getGamma()); varianceVec.push_back(variance); paramIdVec.push_back(2); meanVec.push_back(rs.getS()); varianceVec.push_back(variance); scatterVary = models::getScatterPlotVaryParameters(rs, paramIdVec, meanVec, varianceVec, normalizationDose, 50); scatterCompare.setCompareScatter(scatterVary); CHECK_TESTER(scatterCompare); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb \ No newline at end of file diff --git a/testing/examples/RTDVHCalculatorTest.cpp b/testing/examples/RTDVHCalculatorTest.cpp deleted file mode 100644 index fe89538..0000000 --- a/testing/examples/RTDVHCalculatorTest.cpp +++ /dev/null @@ -1,280 +0,0 @@ -// ----------------------------------------------------------------------- -// RTToolbox - DKFZ radiotherapy quantitative evaluation library -// -// (c) Copyright 2007, DKFZ, Heidelberg, Germany -// ALL RIGHTS RESERVED -// -// THIS FILE CONTAINS CONFIDENTIAL AND PROPRIETARY INFORMATION OF DKFZ. -// ANY DUPLICATION, MODIFICATION, DISTRIBUTION, OR -// DISCLOSURE IN ANY FORM, IN WHOLE, OR IN PART, IS STRICTLY PROHIBITED -// WITHOUT THE PRIOR EXPRESS WRITTEN PERMISSION OF DKFZ. -// -//------------------------------------------------------------------------ -/*! -// @file -// @version $Revision$ (last changed revision) -// @date $Date$ (last change date) -// @author zhangl (last changed by) -// @author *none* (Reviewer) -// @author zhangl (Programmer) -// -// Subversion HeadURL: $HeadURL: http://sidt-hpc1/dkfz_repository/NotMeVisLab/SIDT/RTToolbox/branch/restructure/testing/core/DVHCalculatorTest.cpp $ -*/ - -// this file defines the rttbCoreTests for the test driver -// and all it expects is that you have a function called RegisterTests - - -#include "litCheckMacros.h" - - -#include "dcmtk/config/osconfig.h" /* make sure OS specific configuration is included first */ - -#include "dcmtk/dcmrt/drtdose.h" -#include "dcmtk/dcmrt/drtstruct.h" - -#include "../../code/core/rttbMaskedDoseIteratorInterface_NEW.h" -#include "../../code/core/rttbGenericMaskedDoseIterator_NEW.h" -#include "../../code/core/rttbGenericDoseIterator_NEW.h" -#include "../../code/core/rttbDVHCalculator.h" -#include "../../code/io/dicom/rttbDicomStructureSet.h" -#include "../../code/core/rttbStructure.h" -#include "../../code/core/rttbNullPointerException.h" -#include "../../code/core/rttbInvalidParameterException.h" -#include "../../code/core/rttbInvalidDoseException.h" -#include "../../code/core/rttbException.h" -#include "../../code/core/rttbBaseType_NEW.h" - -//#include "CreateTestStructure.h" - currently not used - -// sollen externe Methoden verwendet werden? -#include "../../code/io/dicom/rttbDicomDoseAccessor.h" -#include "../../code/masks/rttbOTBMaskAccessor.h" -#include "../../code/algorithms/rttbDoseStatistics.h" - -#include - - - - - -namespace rttb -{ - - namespace testing - { - - - - /*! @brief DVHCalculatorTest. - 1. DVH calculator test: use MaskedDicomDoseIterator(&rtdose, &rtstr) -> now use rttbGenericMaskedDoseIterator_NEW.h - 2. DVH Calculator test: use cache of MaskedDicomDoseIterator resetDose -> obsolete - 3. DVH Calculator test: use MaskedDicomDoseIteratorRandomSampling ->obsolete use dummy mask? - 4. DVH Calculator test: use DicomDoseIteratorUseMask -> obsolete/ ../io/dicom? - 5. DVH Calculator test: use manually generated structure and Mask(DoseIteratorInterface*, Structure* ) -> check if this should be done here or in ../masks - 6. DVHCalculation Test using constructor Mask( DoseIteratorInterface* aDoseIter , std::vector aDoseVoxelVector ); -> mask generation should be tested in mask tests use default mask here - */ - - - /*static std::string MASK_FILENAME="../../../RTToolbox/testing/data/DICOM/TestDose/TestMask.dcm"; - static std::string HITDOSE_FILENAME="d:/Dotmobi-SVN/rondo/RondoInterface/Resources/TestData/Dicom/HITDemoData/RTDOSE/HITDemoData_20100316T175500_RTDOSE_450.IMA"; - static std::string HITSTRUCT_FILENAME="d:/Dotmobi-SVN/rondo/RondoInterface/Resources/TestData/Dicom/HITDemoData/RTSTRUCT/HITDemoData_20100407T152606_RTSTRUCT_0.IMA"; - static std::string HITBEAMDOSE_FILENAME="d:/Dotmobi-SVN/rondo/RondoInterface/Resources/TestData/Dicom/HITDemoData/RTDOSE/HITDemoData_20100407T152314_RTDOSE_0.IMA"; - - static std::string DVH_FILENAME_WRITE="../../../RTToolbox/testing/data/TestDVH/dvh_write.txt"; - - static std::string DVH_FILENAME_PTV_HIT="../../../RTToolbox/testing/data/TestDVH/dvh_PTV_HIT.txt";*/ - - int RTDVHCalculatorTest(int argc, char* argv[]) - { - typedef core::GenericDoseIterator::DoseAcessorPointer DoseAcessorPointer; - typedef core::DoseIteratorInterface::DoseIteratorPointer DoseIteratorPointer; - typedef core::MaskedDoseIteratorInterface::MaskAcessorPointer MaskAcessorPointer; - typedef masks::OTBMaskAccessor::StructTypePointer StructTypePointer; - typedef algorithms::DoseStatistics::ResultListPointer ResultListPointer; - - - PREPARE_DEFAULT_TEST_REPORTING; - //ARGUMENTS: 1: structure file name - // 2: dose1 file name - // 3: dose2 file name - // 4: dose3 file name - - std::string STRUCT_FILENAME; - std::string DOSE_FILENAME; - std::string DOSE2_FILENAME; - std::string DOSE3_FILENAME; - - /*STRUCT_FILENAME="d:/RTToolboxProj/trunk/testing/data/DICOM/H000090/RTSTRUCT/H000090_20100610T144958_RTSTRUCT_479.IMA"; - DOSE_FILENAME="d:/RTToolboxProj/trunk/testing/data/DICOM/H000090/RTDOSE/H000090_20100610T144958_RTDOSE_396.IMA"; - DOSE2_FILENAME="d:/RTToolboxProj/trunk/testing/data/DICOM/H000090/RTDOSE/H000090_20100610T144958_RTDOSE_397.IMA"; - DOSE3_FILENAME="d:/RTToolboxProj/trunk/testing/data/DICOM/H000090/RTDOSE/H000090_20100610T144958_RTDOSE_398.IMA";*/ - - - if (argc > 1) - { - STRUCT_FILENAME = argv[1]; - } - - if (argc > 2) - { - DOSE_FILENAME = argv[2]; - } - - if (argc > 3) - { - DOSE2_FILENAME = argv[3]; - } - - if (argc > 4) - { - DOSE3_FILENAME = argv[4]; - } - - - // CHECK_THROW_EXPLICIT(core::DVHCalculator(NULL,"",""),core::NullPointerException); -> tested in individual test - - OFCondition status; - DcmFileFormat fileformat; - - /* read dicom-rt dose */ - ::DRTDoseIOD rtdose; - // rttb::core::Structure* rtstr; - double deltaD = 0; - double deltaV = 0.03; - // use new architecture accessor+iterator - // import dicom files - CHECK_NO_THROW(DoseAcessorPointer spDoseAccessor(new core::DicomDoseAccessor( - DOSE_FILENAME.c_str()))); - DoseAcessorPointer spDoseAccessor(new core::DicomDoseAccessor(DOSE_FILENAME.c_str())); - std::cout << "ProcessDose " << DOSE_FILENAME.c_str() << std::endl; - core::GenericDoseIterator ddit(spDoseAccessor); - - /* read dicom-rt structure set*/ - // use toolbox extensions from Andreas?!: DoseToolDCMStructLoader.h - ::DRTStructureSetIOD rtstruct; - status = fileformat.loadFile(STRUCT_FILENAME.c_str()); - CHECK(status.good()); - - if (!status.good()) - { - std::cerr << "Error: load rt structure set loadFile(STRUCT_FILENAME) failed!" << std::endl; - return -1; - } - - status = rtstruct.read(*fileformat.getDataset()); - - if (!status.good()) - { - std::cerr << "Error: read DRTStructureSetIOD failed!" << std::endl; - return -1; - } - - //get maximum of rtdose - - CHECK_NO_THROW(rttb::core::DicomStructureSet rtstrset = rttb::core::DicomStructureSet(&rtstruct)); - rttb::core::DicomStructureSet rtstrset = rttb::core::DicomStructureSet(&rtstruct); - - clock_t start(clock()); - //std::vector mditVector; - ResultListPointer resultListMax; - ResultListPointer resultListMin; - - if (rtstrset.getNumberOfStructures() > 0) - { - for (int j = 0; j < rtstrset.getNumberOfStructures(); j++) - { - clock_t startLoop(clock()); - StructTypePointer spRtstr = boost::make_shared(*(rtstrset.getStructure(j))); - rttb::core::GeometricInfo geoInf = ddit.getGeometricInfo(); - - boost::shared_ptr pOTBMaskAccessor = - boost::make_shared(spRtstr, geoInf, rtdose); - MaskAcessorPointer spMaskAccessor(pOTBMaskAccessor); - boost::shared_ptr spTestMaskedDoseIterator = - boost::make_shared(spMaskAccessor, - spDoseAccessor); //using Mask - DoseIteratorPointer mddit(spTestMaskedDoseIterator); - mddit->reset(); - std::cout << "Initialization time: " << (clock() - startLoop) / 1000 << " s" << std::endl; - - rttb::core::DVHCalculator calc(mddit, spRtstr->getUID(), mddit->getDoseUID()); - rttb::core::DVH dvh = *(calc.getDVH()); - std::cout << "DVH(mask) max: " << dvh.getMaximum() << "in " << spRtstr->getVolume() << " cm^3" << - std::endl; - } - } - - clock_t finish(clock()); - - std::cout << "DVH Calculation time: " << (finish - start) / 60000 << " min" << std::endl; - - DoseAcessorPointer spDoseAccessor2(new core::DicomDoseAccessor(DOSE2_FILENAME.c_str())); - std::cout << "ProcessDose " << DOSE2_FILENAME.c_str() << std::endl; - - start = clock(); - resultListMax.reset(); - resultListMin.reset(); - - if (rtstrset.getNumberOfStructures() > 0) - { - for (int j = 0; j < rtstrset.getNumberOfStructures(); j++) - { - clock_t startLoop(clock()); - StructTypePointer spRtstr = boost::make_shared(*(rtstrset.getStructure(j))); - rttb::core::GeometricInfo geoInf = ddit.getGeometricInfo(); - - boost::shared_ptr pOTBMaskAccessor = - boost::make_shared(spRtstr, geoInf, rtdose); - MaskAcessorPointer spMaskAccessor(pOTBMaskAccessor); - boost::shared_ptr spTestMaskedDoseIterator = - boost::make_shared(spMaskAccessor, - spDoseAccessor2); //using Mask - DoseIteratorPointer mddit(spTestMaskedDoseIterator); - boost::shared_ptr spTestDoseIterator = - boost::make_shared(spDoseAccessor2); - DoseIteratorPointer ddit2(spTestDoseIterator); - mddit->reset(); - std::cout << "Initialization time: " << (clock() - startLoop) / 1000 << " s" << std::endl; - - rttb::core::DVHCalculator calc2(mddit, spRtstr->getUID(), mddit->getDoseUID()); - rttb::core::DVH dvh2 = *(calc2.getDVH()); - std::cout << "DVH(mask) max: " << dvh2.getMaximum() << "in " << spRtstr->getVolume() << " cm^3" << - std::endl; - - rttb::core::DVHCalculator calc(ddit2, spRtstr->getUID(), ddit2->getDoseUID()); - rttb::core::DVH dvh = *(calc.getDVH()); - std::cout << "DVH(dose) max: " << dvh.getMaximum() << std::endl; - - - //test statistics here? -> ../algorithms - clock_t startStats(clock()); - std::cout << "Statistics: " << std::endl; - rttb::algorithms::DoseStatistics stat(mddit); - std::cout << "Statistics time: " << (clock() - startStats) << " s" << std::endl; - std::cout << "max: " << stat.getMaximum(resultListMax) << std::endl; - std::cout << "mean: " << stat.getMean() << std::endl; - std::cout << "min: " << stat.getMinimum(resultListMin) << std::endl; - std::cout << "std: " << stat.getStdDeviation() << std::endl; - std::cout << "var: " << stat.getVariance() << std::endl; - } - } - - finish = clock(); - - std::cout << "Calculation time: " << (finish - start) / 60000 << " min" << std::endl; - - - - - RETURN_AND_REPORT_TEST_SUCCESS; - - - } - - - - - }//testing -}//rttb - diff --git a/testing/examples/RTDVHTest.cpp b/testing/examples/RTDVHTest.cpp index 233d9b4..f0ea1a4 100644 --- a/testing/examples/RTDVHTest.cpp +++ b/testing/examples/RTDVHTest.cpp @@ -1,115 +1,114 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (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 "litCheckMacros.h" #include "rttbDVH.h" #include "rttbDVHTxtFileReader.h" #include "rttbBaseType.h" #include "rttbDvhBasedModels.h" namespace rttb { namespace testing { /*! @brief RTDVHTest. Max, min, mean, modal, median, Vx, Dx, EUD, BED, LQED2 are tested. Test if calculation in new architecture returns similar results to the original implementation. Comparison of actual DVH values is performed in DVHCalculatorComparisonTest.cpp. WARNING: The values for comparison need to be adjusted if the input files are changed! */ int RTDVHTest(int argc, char* argv[]) { PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: dvh file name std::string RTDVH_FILENAME_PTV; if (argc > 1) { RTDVH_FILENAME_PTV = argv[1]; } typedef core::DVH::DVHPointer DVHPointer; /*test RT dvh*/ rttb::io::other::DVHTxtFileReader dvhReader = rttb::io::other::DVHTxtFileReader(RTDVH_FILENAME_PTV); const DoseCalcType expectedValue = 0.01305; //dvhReader DVHPointer dvh = dvhReader.generateDVH(); CHECK_CLOSE(expectedValue, models::getEUD(dvh, 10), errorConstant); std::cout << models::getEUD(dvh, 10) << std::endl; - std::map bedMap = models::calcBEDDVH(dvh, 15, 15); - std::map LqedMap = models::calcLQED2DVH(dvh, 15, 10); + CHECK_NO_THROW(models::calcBEDDVH(dvh, 15, 15)); CHECK_NO_THROW(models::calcLQED2DVH(dvh, 15, 10)); CHECK_NO_THROW(dvh->getDataDifferential(true)); CHECK_EQUAL(1, dvh->calcCumulativeDVH(true).at(0)); CHECK_NO_THROW(models::calcBEDDVH(dvh, 15, 15, true)); CHECK_NO_THROW(models::calcLQED2DVH(dvh, 15, 10, true)); - //test statistics (relative comulative data) + //test statistics (relative cumulative data) CHECK_CLOSE(expectedValue, dvh->getMaximum(), errorConstant); CHECK_CLOSE(expectedValue, dvh->getMinimum(), errorConstant); CHECK_CLOSE(expectedValue, dvh->getMean(), errorConstant); CHECK_CLOSE(expectedValue, dvh->getMedian(), errorConstant); CHECK_CLOSE(expectedValue, dvh->getModal(), errorConstant); CHECK_EQUAL(0, dvh->getVx(0.014)); CHECK_EQUAL(0.125, dvh->getVx(0.01)); - rttb::DoseTypeGy dTest = dvh->getDx(100); + CHECK_CLOSE(0.0131, dvh->getDx(100), errorConstant + errorConstant * 10); CHECK_CLOSE(0.013, dvh->getDx(249), errorConstant); - dTest = dvh->getDx(249); + CHECK_EQUAL(0, dvh->getDx(251)); - dTest = dvh->getDx(251); - //test statistics (absolute comulative data) + + //test statistics (absolute cumulative data) CHECK_EQUAL(2000, dvh->calcCumulativeDVH(false).at(0)); CHECK_EQUAL(0, dvh->getVx(0.014)); CHECK_EQUAL(250, dvh->getVx(0.01)); - dTest = dvh->getDx(100); + CHECK_CLOSE(0.0131, dvh->getDx(100), errorConstant + errorConstant * 10); CHECK_CLOSE(0.013, dvh->getDx(249), errorConstant); - dTest = dvh->getDx(249); + CHECK_EQUAL(0, dvh->getDx(251)); - dTest = dvh->getDx(251); + RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/examples/VoxelizationValidationTest.cpp b/testing/examples/VoxelizationValidationTest.cpp index 8b33bb8..d8324a6 100644 --- a/testing/examples/VoxelizationValidationTest.cpp +++ b/testing/examples/VoxelizationValidationTest.cpp @@ -1,201 +1,198 @@ // ----------------------------------------------------------------------- // 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" 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; if (argc > 4) { RTSTRUCT_FILENAME = argv[1]; RTDOSE_FILENAME = argv[2]; BoostMask_DIRNAME = argv[3]; OTBMask_DIRNAME = argv[4]; } 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. 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());*/ //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) //7: Ref.Pkt, mean = -1.#IND { CHECK_CLOSE(dvh.getMean(), dvh2.getMean(), 0.1); } CHECK_CLOSE(dvh.getMedian(), dvh2.getMedian(), 0.1); CHECK_CLOSE(dvh.getModal(), dvh2.getModal(), 0.1); //0: Aussenkontur and 3: Niere li. failed. if (j != 0 && j != 3) { CHECK_CLOSE(dvh.getVx(0), dvh2.getVx(0), dvh.getVx(0) * 0.05); //check volume difference < 5% } } } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/interpolation/InterpolationMatchPointTransformation/simpleRegistrationWorkflow.cpp b/testing/interpolation/InterpolationMatchPointTransformation/simpleRegistrationWorkflow.cpp index e2c0502..5a41fbb 100644 --- a/testing/interpolation/InterpolationMatchPointTransformation/simpleRegistrationWorkflow.cpp +++ b/testing/interpolation/InterpolationMatchPointTransformation/simpleRegistrationWorkflow.cpp @@ -1,113 +1,110 @@ // ----------------------------------------------------------------------- // 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) */ #undef MAP_SEAL_ALGORITHMS #include "simpleRegistrationWorkflow.h" #include "registrationHelper.h" simpleRegistrationWorkflow::simpleRegistrationWorkflow(std::string targetFilename, std::string movingFilename, - bool isDirectory) + bool isDirectory) : _targetFilename(targetFilename), _movingFilename(movingFilename) { - _targetFilename = targetFilename; - _movingFilename = movingFilename; - setImageFileNames(_targetFilename, _movingFilename, isDirectory, globals); loadData(globals); _spAlgorithmEuler = NULL; } vnl_vector simpleRegistrationWorkflow::getRegistrationParameters( Registration3D3DPointer reg) { typedef map::core::ModelBasedRegistrationKernel<3, 3> ModelBasedRegistrationKernel3D3D; const ModelBasedRegistrationKernel3D3D* pModelBasedDirectKernel3D3D = dynamic_cast(&(reg->getDirectMapping())); if (pModelBasedDirectKernel3D3D) { ModelBasedRegistrationKernel3D3D::ParametersType params = pModelBasedDirectKernel3D3D->getTransformModel()->GetParameters(); return params; } else { return vnl_vector(); } } void simpleRegistrationWorkflow::initializeAndPerformRegistration() { _spAlgorithmEuler = AlgorithmTypeEuler::New(); _spAlgorithmEuler->setProperty("PreinitTransform", map::core::MetaProperty::New(true)); _spAlgorithmEuler->setMovingImage(globals.spMovingImage); _spAlgorithmEuler->setTargetImage(globals.spTargetImage); AlgorithmTypeEuler::RegistrationType::Pointer spRegistration; try { spRegistration = _spAlgorithmEuler->getRegistration(); } catch (const map::core::ExceptionObject& e) { std::cerr << "caught an MatchPoint exception:\n"; e.Print(std::cerr); std::cerr << "\n"; } catch (const itk::ExceptionObject& e) { std::cerr << "caught an ITK exception:\n"; std::cerr << e.GetFile() << ":" << e.GetLine() << ":\n" << e.GetDescription() << "\n"; } catch (const std::exception& e) { std::cerr << "caught an exception:\n"; std::cerr << e.what() << "\n"; } catch (...) { std::cerr << "caught an unknown exception!!!\n"; } } map::core::Registration<3, 3>::Pointer simpleRegistrationWorkflow::getRegistration() { if (_spAlgorithmEuler.IsNull()) { initializeAndPerformRegistration(); } return _spAlgorithmEuler->getRegistration(); }; const itk::Image* simpleRegistrationWorkflow::getTargetImage() { if (_spAlgorithmEuler.IsNull()) { initializeAndPerformRegistration(); } return _spAlgorithmEuler->getTargetImage(); }; diff --git a/testing/io/dicom/DicomStructureSetGeneratorTest.cpp b/testing/io/dicom/DicomStructureSetGeneratorTest.cpp index 65e4e41..5f43d25 100644 --- a/testing/io/dicom/DicomStructureSetGeneratorTest.cpp +++ b/testing/io/dicom/DicomStructureSetGeneratorTest.cpp @@ -1,106 +1,106 @@ // ----------------------------------------------------------------------- // 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 "rttbDicomFileStructureSetGenerator.h" #include "rttbDicomIODStructureSetGenerator.h" #include "rttbDcmrtException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace testing { /*!@brief DicomIOTest - test structure set generator for dicom data 1) test dicom file structure set generator 2) test dicom IOD structure set generator */ int DicomStructureSetGeneratorTest(int argc, char* argv[]) { typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; //typedef boost::shared_ptr DRTStrSetIODPtr; typedef io::dicom::DicomIODStructureSetGenerator::DRTStrSetIODPtr DRTStrSetIODPtr; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: structure file name std::string RTSTRUCT_FILENAME; if (argc > 1) { RTSTRUCT_FILENAME = argv[1]; } /* structure set */ //1) test dicom file structure set generator CHECK_NO_THROW(io::dicom::DicomFileStructureSetGenerator("")); CHECK_NO_THROW(io::dicom::DicomFileStructureSetGenerator("Test.test")); CHECK_THROW_EXPLICIT(io::dicom::DicomFileStructureSetGenerator("Test.test").generateStructureSet(), rttb::core::InvalidParameterException); CHECK_NO_THROW(io::dicom::DicomFileStructureSetGenerator(RTSTRUCT_FILENAME.c_str())); StructureSetPointer rtStructureSet = io::dicom::DicomFileStructureSetGenerator( RTSTRUCT_FILENAME.c_str()).generateStructureSet(); //2) test dicom IOD structure set generator OFCondition status; DcmFileFormat fileformat; boost::shared_ptr drtStrSetIOD = boost::make_shared(); CHECK_NO_THROW(io::dicom::DicomIODStructureSetGenerator generator1(drtStrSetIOD)); CHECK_THROW_EXPLICIT(io::dicom::DicomIODStructureSetGenerator(drtStrSetIOD).generateStructureSet(), core::InvalidParameterException); - status = fileformat.loadFile(RTSTRUCT_FILENAME.c_str()); - status = drtStrSetIOD->read(*fileformat.getDataset()); + fileformat.loadFile(RTSTRUCT_FILENAME.c_str()); + drtStrSetIOD->read(*fileformat.getDataset()); CHECK_NO_THROW(io::dicom::DicomIODStructureSetGenerator generator2(drtStrSetIOD)); CHECK_NO_THROW(io::dicom::DicomIODStructureSetGenerator(drtStrSetIOD).generateStructureSet()); StructureSetPointer rtStructureSet2 = io::dicom::DicomIODStructureSetGenerator( drtStrSetIOD).generateStructureSet(); CHECK_EQUAL(rtStructureSet2->getNumberOfStructures(), rtStructureSet->getNumberOfStructures()); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb