diff --git a/code/core/rttbStructureSet.cpp b/code/core/rttbStructureSet.cpp index 0cc4a3d..fd3306e 100644 --- a/code/core/rttbStructureSet.cpp +++ b/code/core/rttbStructureSet.cpp @@ -1,89 +1,89 @@ // ----------------------------------------------------------------------- // 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 "rttbStructureSet.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace core { StructureSet::StructureSet() {} StructureSet::StructureSet(std::vector aStructureVector, IDType aPatientUID, IDType aUID) { _structureSetVector = aStructureVector; _patientUID = aPatientUID; _UID = aUID; if (_UID == "") { boost::uuids::uuid id; boost::uuids::random_generator generator; id = generator(); std::stringstream ss; ss << id; _UID = ss.str(); } } StructureSet::StructTypePointer StructureSet::getStructure(int aStructureNo) const { - int size = this->getNumberOfStructures() - 1; + int size = static_cast(this->getNumberOfStructures() - 1); if (aStructureNo < 0 || aStructureNo > size) { std::stringstream sstr; sstr << "aStructureNo must between 0 and " << size; throw InvalidParameterException(sstr.str()); } return _structureSetVector.at(aStructureNo); } StructureSet::NumberOfStructuresType StructureSet::getNumberOfStructures() const { return _structureSetVector.size(); } IDType StructureSet::getUID() { return _UID; } IDType StructureSet::getPatientUID() { return _patientUID; } }//end namespace core }//end namespace rttb diff --git a/code/io/dicom/rttbDicomFileDoseAccessorWriter.cpp b/code/io/dicom/rttbDicomFileDoseAccessorWriter.cpp index a7ec2a5..6629e73 100644 --- a/code/io/dicom/rttbDicomFileDoseAccessorWriter.cpp +++ b/code/io/dicom/rttbDicomFileDoseAccessorWriter.cpp @@ -1,259 +1,255 @@ // ----------------------------------------------------------------------- // 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 "rttbDicomFileDoseAccessorWriter.h" #include "rttbInvalidDoseException.h" #include "rttbGeometricInfo.h" #include "rttbGeometricInfo.h" #include "rttbGenericDoseIterator.h" #include "rttbDoseStatisticsCalculator.h" namespace rttb { namespace io { namespace dicom { DicomFileDoseAccessorWriter::DicomFileDoseAccessorWriter() { _doseIOD = boost::make_shared(); _dataset = _fileformat.getDataset(); } void DicomFileDoseAccessorWriter::setFileName(DICOMRTFileNameString aFileName) { _fileName = aFileName; } bool DicomFileDoseAccessorWriter::process() { OFCondition status; /* Prepare dcmtk */ DcmItem* dcm_item = 0; //get geometric info rttb::core::GeometricInfo geometricInfo = _doseAccessor->getGeometricInfo(); /* ----------------------------------------------------------------- */ /* Part 1 -- General header */ /* ----------------------------------------------------------------- */ OFString CreationUID(_doseAccessor->getUID().c_str()); _dataset->putAndInsertString(DCM_ImageType, "DERIVED\\SECONDARY\\REFORMATTED"); _dataset->putAndInsertOFStringArray(DCM_InstanceCreationDate, "");//Creation Date _dataset->putAndInsertOFStringArray(DCM_InstanceCreationTime, "");//Creation Time _dataset->putAndInsertOFStringArray(DCM_InstanceCreatorUID, CreationUID); _dataset->putAndInsertString(DCM_SOPClassUID, UID_RTDoseStorage); _dataset->putAndInsertString(DCM_SOPInstanceUID, _doseAccessor->getUID().c_str()); _dataset->putAndInsertOFStringArray(DCM_StudyDate, ""); _dataset->putAndInsertOFStringArray(DCM_StudyTime, ""); _dataset->putAndInsertOFStringArray(DCM_AccessionNumber, ""); _dataset->putAndInsertOFStringArray(DCM_Modality, "RTDOSE"); _dataset->putAndInsertString(DCM_Manufacturer, "RTToolbox"); _dataset->putAndInsertString(DCM_InstitutionName, ""); _dataset->putAndInsertString(DCM_ReferringPhysicianName, ""); _dataset->putAndInsertString(DCM_StationName, ""); _dataset->putAndInsertString(DCM_ManufacturerModelName, "RTToolbox"); /* (0008,1140) DCM_ReferencedImageSequence -- MIM likes this */ dcm_item = 0; _dataset->findOrCreateSequenceItem( DCM_ReferencedImageSequence, dcm_item, -2); dcm_item->putAndInsertString(DCM_ReferencedSOPClassUID, UID_CTImageStorage); dcm_item->putAndInsertString(DCM_ReferencedSOPInstanceUID, ""); _dataset->putAndInsertString(DCM_PatientName, ""); _dataset->putAndInsertString(DCM_PatientID, ""); _dataset->putAndInsertString(DCM_PatientBirthDate, ""); _dataset->putAndInsertString(DCM_PatientSex, "O"); _dataset->putAndInsertString(DCM_SliceThickness, boost::lexical_cast(geometricInfo.getSliceThickness()).c_str()); _dataset->putAndInsertString(DCM_SoftwareVersions, ""); _dataset->putAndInsertString(DCM_StudyInstanceUID, ""); _dataset->putAndInsertString(DCM_SeriesInstanceUID, ""); _dataset->putAndInsertString(DCM_StudyID, "10001"); _dataset->putAndInsertString(DCM_SeriesNumber, ""); _dataset->putAndInsertString(DCM_InstanceNumber, "1"); /* GCS FIX: PatientOrientation */ std::ostringstream sstr; sstr << geometricInfo.getImagePositionPatient().x() << "\\" << geometricInfo.getImagePositionPatient().y() << "\\" << geometricInfo.getImagePositionPatient().z(); _dataset->putAndInsertString(DCM_PatientOrientation, "L/P"); _dataset->putAndInsertString(DCM_ImagePositionPatient, sstr.str().c_str()); sstr.str(""); sstr << geometricInfo.getImageOrientationRow().x() << "\\" << geometricInfo.getImageOrientationRow().y() << "\\" << geometricInfo.getImageOrientationRow().z() << "\\" << geometricInfo.getImageOrientationColumn().x() << "\\" << geometricInfo.getImageOrientationColumn().y() << "\\" << geometricInfo.getImageOrientationColumn().z(); _dataset->putAndInsertString(DCM_ImageOrientationPatient, sstr.str().c_str()); _dataset->putAndInsertString(DCM_FrameOfReferenceUID, ""); _dataset->putAndInsertString(DCM_SamplesPerPixel, "1"); _dataset->putAndInsertString(DCM_PhotometricInterpretation, "MONOCHROME2"); sstr.str(""); sstr << geometricInfo.getNumSlices(); _dataset->putAndInsertString(DCM_NumberOfFrames, sstr.str().c_str()); /* GCS FIX: Add FrameIncrementPointer */ _dataset->putAndInsertString(DCM_FrameIncrementPointer, "(3004,000c)"); _dataset->putAndInsertUint16(DCM_Rows, geometricInfo.getNumRows()); _dataset->putAndInsertUint16(DCM_Columns, geometricInfo.getNumColumns()); sstr.str(""); sstr << geometricInfo.getSpacing()(1) << "\\" << geometricInfo.getSpacing()(0); _dataset->putAndInsertString(DCM_PixelSpacing, sstr.str().c_str()); _dataset->putAndInsertString(DCM_BitsAllocated, "32"); _dataset->putAndInsertString(DCM_BitsStored, "32"); _dataset->putAndInsertString(DCM_HighBit, "31"); _dataset->putAndInsertString(DCM_DoseUnits, "GY"); _dataset->putAndInsertString(DCM_DoseSummationType, "PLAN"); sstr.str("0"); for (int i = 1; i < geometricInfo.getNumSlices(); i++) { sstr << "\\" << i* geometricInfo.getSpacing()(2); } _dataset->putAndInsertString(DCM_GridFrameOffsetVector, sstr.str().c_str()); /* We need to convert image to uint16_t, but first we need to scale it so that the maximum dose fits in a 16-bit unsigned integer. Compute an appropriate scaling factor based on the maximum dose. */ /* Find the maximum value in the image */ boost::shared_ptr spTestDoseIterator = boost::make_shared(_doseAccessor); rttb::core::GenericDoseIterator::DoseIteratorPointer spDoseIterator(spTestDoseIterator); rttb::algorithms::DoseStatisticsCalculator myDoseStatsCalculator(spDoseIterator); rttb::algorithms::DoseStatistics::DoseStatisticsPointer doseStatistics = myDoseStatsCalculator.calculateDoseStatistics(); - //rttb::algorithms::DoseStatistics doseStat(spDoseIterator); - //boost::shared_ptr< std::vector > > myResultPairs = - // boost::make_shared< std::vector > >(); - //rttb::algorithms::DoseStatistics::ResultListPointer spMyResultPairs(myResultPairs); double maxDose = doseStatistics->getMaximum(); /* Find scale factor */ - float dose_scale; + double dose_scale; dose_scale = maxDose / PixelDataMaxValue; /* Scale the image and add scale factor to _dataset */ sstr.str(""); sstr << dose_scale; _dataset->putAndInsertString(DCM_DoseGridScaling, sstr.str().c_str()); /* (300c,0002) ReferencedRTPlanSequence -- for future expansion */ dcm_item = 0; _dataset->findOrCreateSequenceItem( DCM_ReferencedRTPlanSequence, dcm_item, -2); dcm_item->putAndInsertString(DCM_ReferencedSOPClassUID, UID_RTPlanStorage); dcm_item->putAndInsertString(DCM_ReferencedSOPInstanceUID, ""); /* (300c,0060) DCM_ReferencedStructureSetSequence -- MIM likes this */ dcm_item = 0; _dataset->findOrCreateSequenceItem( DCM_ReferencedStructureSetSequence, dcm_item, -2); dcm_item->putAndInsertString(DCM_ReferencedSOPClassUID, UID_RTStructureSetStorage); dcm_item->putAndInsertString(DCM_ReferencedSOPInstanceUID, ""); /* Convert image bytes to integer, then add to _dataset */ Uint16* pixelData; - int pixelCount = geometricInfo.getNumRows() * geometricInfo.getNumColumns() * - geometricInfo.getNumSlices(); + unsigned int pixelCount = static_cast(geometricInfo.getNumRows() * geometricInfo.getNumColumns() * + geometricInfo.getNumSlices()); pixelData = new Uint16[pixelCount]; for (unsigned int i = 0; i < pixelCount; ++i) { double doseValue = _doseAccessor->getValueAt(i); double pixelValue = doseValue / dose_scale; if (pixelValue > PixelDataMaxValue) { pixelValue = PixelDataMaxValue; } pixelData[i] = boost::numeric_cast(pixelValue); } status = _dataset->putAndInsertUint16Array(DCM_PixelData, pixelData, pixelCount); if (!status.good()) { throw core::InvalidDoseException("Error: put and insert pixel data failed!"); } //Write dose to file status = _fileformat.saveFile(_fileName.c_str(), EXS_LittleEndianExplicit); if (status.bad()) { std::cerr << "Error: cannot write DICOM RTDOSE!" << std::endl; } return true; } }//end namespace itk }//end namespace io }//end namespace rttb diff --git a/code/io/helax/rttbDicomHelaxDoseAccessor.cpp b/code/io/helax/rttbDicomHelaxDoseAccessor.cpp index 850c9c0..13c1015 100644 --- a/code/io/helax/rttbDicomHelaxDoseAccessor.cpp +++ b/code/io/helax/rttbDicomHelaxDoseAccessor.cpp @@ -1,328 +1,329 @@ // ----------------------------------------------------------------------- // 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 "boost/filesystem/operations.hpp" #include "boost/filesystem/path.hpp" #include "boost/progress.hpp" #include #include "rttbDicomHelaxDoseAccessor.h" #include "rttbInvalidDoseException.h" #include "rttbDcmrtException.h" #include "rttbIndexOutOfBoundsException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace io { namespace helax { DicomHelaxDoseAccessor::~DicomHelaxDoseAccessor() { } DicomHelaxDoseAccessor::DicomHelaxDoseAccessor(std::vector aDICOMRTDoseVector) { for (int i = 0; i < aDICOMRTDoseVector.size(); i++) { _doseVector.push_back(aDICOMRTDoseVector.at(i)); } this->begin(); } bool DicomHelaxDoseAccessor::begin() { if (_doseVector.size() == 0) { throw core::InvalidParameterException(" The size of aDICOMRTDoseVector is 0!"); } assembleGeometricInfo(); _doseData.clear(); OFString doseGridScalingStr; _doseVector.at(0)->getDoseGridScaling( doseGridScalingStr);//get the first dose grid scaling as _doseGridScaling try { _doseGridScaling = boost::lexical_cast(doseGridScalingStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Dose grid scaling not readable or = 0!") ; } for (int i = 0; i < _doseVector.size(); i++) { DRTDoseIODPtr dose = _doseVector.at(i); OFString currentDoseGridScalingStr; dose->getDoseGridScaling(currentDoseGridScalingStr); double currentDoseGridScaling; try { currentDoseGridScaling = boost::lexical_cast(currentDoseGridScalingStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Dose grid scaling not readable or = 0!") ; } OFCondition status; DcmFileFormat fileformat; DcmItem doseitem; status = dose->write(doseitem); if (status.good()) { unsigned long count; const Uint16* pixelData; status = doseitem.findAndGetUint16Array(DcmTagKey(0x7fe0, 0x0010), pixelData, &count); if (status.good()) { for (unsigned int i = 0; i < static_cast(_geoInfo.getNumColumns()*_geoInfo.getNumRows()); i++) { - this->_doseData.push_back(pixelData[i]*currentDoseGridScaling / - _doseGridScaling); //recalculate dose data + Uint16 data = static_cast(pixelData[i] * currentDoseGridScaling / + _doseGridScaling); + this->_doseData.push_back(data); //recalculate dose data } } else { throw dicom::DcmrtException("Read Pixel Data (7FE0,0010) failed!"); } } else { throw dicom::DcmrtException("Read DICOM-RT Dose file failed!"); } } return true; } bool DicomHelaxDoseAccessor::assembleGeometricInfo() { DRTDoseIODPtr dose = _doseVector.at(0); Uint16 temp = 0; dose->getColumns(temp); _geoInfo.setNumColumns(temp); temp = 0; dose->getRows(temp); _geoInfo.setNumRows(temp); OFString numberOfFramesStr; dose->getNumberOfFrames(numberOfFramesStr); if (!numberOfFramesStr.empty()) { _geoInfo.setNumSlices(boost::lexical_cast(numberOfFramesStr.c_str())); } else { _geoInfo.setNumSlices((VoxelGridDimensionType)_doseVector.size()); } if (_geoInfo.getNumColumns() == 0 || _geoInfo.getNumRows() == 0 || _geoInfo.getNumSlices() == 0) { throw core::InvalidDoseException("Empty dicom dose!") ; } OFString imageOrientationRowX; dose->getImageOrientationPatient(imageOrientationRowX, 0); OFString imageOrientationRowY; dose->getImageOrientationPatient(imageOrientationRowY, 1); OFString imageOrientationRowZ; dose->getImageOrientationPatient(imageOrientationRowZ, 2); OFString imageOrientationColumnX; dose->getImageOrientationPatient(imageOrientationColumnX, 3); OFString imageOrientationColumnY; dose->getImageOrientationPatient(imageOrientationColumnY, 4); OFString imageOrientationColumnZ; dose->getImageOrientationPatient(imageOrientationColumnZ, 5); WorldCoordinate3D imageOrientationRow; WorldCoordinate3D imageOrientationColumn; try { imageOrientationRow(0) = boost::lexical_cast(imageOrientationRowX.c_str()); imageOrientationRow(1) = boost::lexical_cast(imageOrientationRowY.c_str()); imageOrientationRow(2) = boost::lexical_cast(imageOrientationRowZ.c_str()); imageOrientationColumn(0) = boost::lexical_cast(imageOrientationColumnX.c_str()); imageOrientationColumn(1) = boost::lexical_cast(imageOrientationColumnY.c_str()); imageOrientationColumn(2) = boost::lexical_cast(imageOrientationColumnZ.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("boost::lexical_cast ImageOrientation failed! Can not read image orientation X/Y/Z!") ; } OrientationMatrix orientation; orientation(0, 0) = imageOrientationRow.x(); orientation(1, 0) = imageOrientationRow.y(); orientation(2, 0) = imageOrientationRow.z(); orientation(0, 1) = imageOrientationColumn.x(); orientation(1, 1) = imageOrientationColumn.y(); orientation(2, 1) = imageOrientationColumn.z(); WorldCoordinate3D perpendicular = imageOrientationRow.cross(imageOrientationColumn); orientation(0, 2) = perpendicular.x(); orientation(1, 2) = perpendicular.y(); orientation(2, 2) = perpendicular.z(); _geoInfo.setOrientationMatrix(orientation); OFString imagePositionX; dose->getImagePositionPatient(imagePositionX, 0); OFString imagePositionY; dose->getImagePositionPatient(imagePositionY, 1); OFString imagePositionZ; dose->getImagePositionPatient(imagePositionZ, 2); WorldCoordinate3D imagePositionPatient; try { imagePositionPatient(0) = boost::lexical_cast(imagePositionX.c_str()); imagePositionPatient(1) = boost::lexical_cast(imagePositionY.c_str()); imagePositionPatient(2) = boost::lexical_cast(imagePositionZ.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("boost::lexical_cast ImagePosition failed! Can not read image position X/Y/Z!") ; } _geoInfo.setImagePositionPatient(imagePositionPatient); SpacingVectorType3D spacingVector; OFString pixelSpacingRowStr; dose->getPixelSpacing(pixelSpacingRowStr, 0); OFString pixelSpacingColumnStr; dose->getPixelSpacing(pixelSpacingColumnStr, 1); try { spacingVector(1) = boost::lexical_cast(pixelSpacingRowStr.c_str()); spacingVector(0) = boost::lexical_cast(pixelSpacingColumnStr.c_str()); } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Can not read Pixel Spacing Row/Column!") ; } _geoInfo.setSpacing(spacingVector); if (_geoInfo.getPixelSpacingRow() == 0 || _geoInfo.getPixelSpacingColumn() == 0) { throw core::InvalidDoseException("Pixel spacing not readable or = 0!"); } OFString sliceThicknessStr; dose->getSliceThickness(sliceThicknessStr); try { spacingVector(2) = boost::lexical_cast(sliceThicknessStr.c_str()); } catch (boost::bad_lexical_cast&) { spacingVector(2) = 0 ;//if no information about slice thickness, set to 0 and calculate it using z coordinate difference between 1. and 2. dose } if (spacingVector(2) == 0) { if (_doseVector.size() > 1) { DRTDoseIODPtr dose2 = _doseVector.at(1);//get the 2. dose OFString imagePositionZ2; dose2->getImagePositionPatient(imagePositionZ2, 2); try { spacingVector(2) = boost::lexical_cast(imagePositionZ2.c_str()) - imagePositionPatient( 2); //caculate slicethickness } catch (boost::bad_lexical_cast&) { throw core::InvalidDoseException("Can not read image position Z of the 2. dose!"); } } else { std::cerr << "sliceThickness == 0! It will be replaced with pixelSpacingRow=" << _geoInfo.getPixelSpacingRow() << "!" << std::endl; spacingVector(2) = spacingVector(0); } } _geoInfo.setSpacing(spacingVector); return true; } GenericValueType DicomHelaxDoseAccessor::getValueAt(const VoxelGridID aID) const { return _doseData.at(aID) * _doseGridScaling; } GenericValueType DicomHelaxDoseAccessor::getValueAt(const VoxelGridIndex3D& aIndex) const { VoxelGridID aVoxelGridID; if (_geoInfo.convert(aIndex, aVoxelGridID)) { return getValueAt(aVoxelGridID); } else { return -1; } } } } } diff --git a/code/io/itk/rttbITKImageMaskAccessor.cpp b/code/io/itk/rttbITKImageMaskAccessor.cpp index be03864..0a7f4c2 100644 --- a/code/io/itk/rttbITKImageMaskAccessor.cpp +++ b/code/io/itk/rttbITKImageMaskAccessor.cpp @@ -1,175 +1,176 @@ // ----------------------------------------------------------------------- // 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 "rttbITKImageMaskAccessor.h" #include "rttbInvalidDoseException.h" namespace rttb { namespace io { namespace itk { ITKImageMaskAccessor::ITKImageMaskAccessor(ITKMaskImageType::ConstPointer aMaskImage) : _mask(aMaskImage) { if (_mask.IsNull()) { throw core::InvalidDoseException("Mask image = 0!") ; } assembleGeometricInfo(); } ITKImageMaskAccessor::~ITKImageMaskAccessor() { }; bool ITKImageMaskAccessor::assembleGeometricInfo() { _geoInfo = boost::make_shared(); _geoInfo->setSpacing(SpacingVectorType3D(_mask->GetSpacing()[0], _mask->GetSpacing()[1], _mask->GetSpacing()[2])); if (_geoInfo->getSpacing()[0] == 0 || _geoInfo->getSpacing()[1] == 0 || _geoInfo->getSpacing()[2] == 0) { throw core::InvalidDoseException("Pixel spacing = 0!"); } _geoInfo->setImagePositionPatient(WorldCoordinate3D(_mask->GetOrigin()[0], _mask->GetOrigin()[1], _mask->GetOrigin()[2])); OrientationMatrix OM(0); for (int col = 0; col < 3; ++col) { for (int row = 0; row < 3; ++row) { OM(col, row) = _mask->GetDirection()(col, row); } } _geoInfo->setOrientationMatrix(OM); _geoInfo->setNumColumns(_mask->GetLargestPossibleRegion().GetSize()[0]); _geoInfo->setNumRows(_mask->GetLargestPossibleRegion().GetSize()[1]); _geoInfo->setNumSlices(_mask->GetLargestPossibleRegion().GetSize()[2]); if (_geoInfo->getNumColumns() == 0 || _geoInfo->getNumRows() == 0 || _geoInfo->getNumSlices() == 0) { throw core::InvalidDoseException("Empty mask!") ; } return true; } void ITKImageMaskAccessor::updateMask() { return; } ITKImageMaskAccessor::MaskVoxelListPointer ITKImageMaskAccessor::getRelevantVoxelVector() { updateMask(); _relevantVoxelVector = getRelevantVoxelVector(0); return _relevantVoxelVector; } ITKImageMaskAccessor::MaskVoxelListPointer ITKImageMaskAccessor::getRelevantVoxelVector( float lowerThreshold) { MaskVoxelListPointer filteredVoxelVectorPointer(new MaskVoxelList); updateMask(); int size = _geoInfo->getNumColumns() * _geoInfo->getNumRows() * _geoInfo->getNumSlices(); filteredVoxelVectorPointer->reserve(size); for (int gridIndex = 0 ; gridIndex < size; gridIndex++) { core::MaskVoxel currentVoxel = core::MaskVoxel(gridIndex); if (getMaskAt(gridIndex, currentVoxel)) { if (currentVoxel.getRelevantVolumeFraction() > lowerThreshold) { filteredVoxelVectorPointer->push_back(currentVoxel); } } } return filteredVoxelVectorPointer; } bool ITKImageMaskAccessor::getMaskAt(VoxelGridID aID, core::MaskVoxel& voxel) const { VoxelGridIndex3D aVoxelGridIndex; if (_geoInfo->convert(aID, aVoxelGridIndex)) { return getMaskAt(aVoxelGridIndex, voxel); } else { return false; } } bool ITKImageMaskAccessor::getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const { voxel.setRelevantVolumeFraction(0); if (_geoInfo->validIndex(aIndex)) { const ITKMaskImageType::IndexType pixelIndex = {{aIndex[0], aIndex[1], aIndex[2]}}; double value = _mask->GetPixel(pixelIndex); VoxelGridID gridId; _geoInfo->convert(aIndex, gridId); if (value >= 0 && value <= 1) { voxel.setRelevantVolumeFraction(value); + return true; } else { std::cerr << "The pixel value of the mask should be >=0 and <=1!" << std::endl; return false; } } else { return false; } } const core::GeometricInfo& ITKImageMaskAccessor::getGeometricInfo() const { return *_geoInfo; }; } } } diff --git a/code/io/itk/rttbITKImageMaskAccessor.h b/code/io/itk/rttbITKImageMaskAccessor.h index 2eeebf0..f3c7cb2 100644 --- a/code/io/itk/rttbITKImageMaskAccessor.h +++ b/code/io/itk/rttbITKImageMaskAccessor.h @@ -1,114 +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) */ #ifndef __ITK_IMAGE_MASK_ACCESSOR_H #define __ITK_IMAGE_MASK_ACCESSOR_H #include "rttbMaskAccessorInterface.h" #include "rttbGeometricInfo.h" #include "rttbBaseType.h" #include "itkImage.h" namespace rttb { namespace io { namespace itk { /*! @class ITKImageMaskAccessor @brief This class gives access to mask information stored in an itk image */ class ITKImageMaskAccessor: public core::MaskAccessorInterface { public: typedef ::itk::Image ITKMaskImageType; typedef ::itk::ImageBase<3> ITKImageBaseType; typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef boost::shared_ptr GeometricInfoPointer; private: /** @brief The mask as itkImage */ ITKMaskImageType::ConstPointer _mask; IDType _maskUID; GeometricInfoPointer _geoInfo; /*! vector containing list of mask voxels*/ MaskVoxelListPointer _relevantVoxelVector; /*! @brief get all required data from the itk image contained in _Mask @exception InvalidDoseException if PixelSpacing is 0 or size in any dimension is 0. */ bool assembleGeometricInfo(); public: ~ITKImageMaskAccessor(); ITKImageMaskAccessor(ITKMaskImageType::ConstPointer aMaskImage); /*! @brief voxelization of the given structures according to the original RTToolbox algorithm*/ void updateMask(); /*! @brief get vector conatining al relevant voxels that are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(); /*! @brief get vector conatining al relevant voxels that have a relevant volume above the given threshold and are inside the given structure*/ MaskVoxelListPointer getRelevantVoxelVector(float lowerThreshold); /*!@brief determine how a given voxel on the dose grid is masked * @param aID ID of the voxel in grid. * @param voxel Reference to the voxel. - * @post after a valid call voxel containes the information of the specified grid voxel. If aID is not valid, voxel values are undefined. + * @post after a valid call voxel contains the information of the specified grid voxel. If aID is not valid, voxel values are undefined. * The relevant volume fraction will be set to zero. - * @return Indicates of the voxel exists and therefore if parameter voxel containes valid values.*/ + * @return Indicates of the voxel exists and therefore if parameter voxel contains valid values.*/ bool getMaskAt(const VoxelGridID aID, core::MaskVoxel& voxel) const; bool getMaskAt(const VoxelGridIndex3D& aIndex, core::MaskVoxel& voxel) const; /*! @brief give access to GeometricInfo*/ const core::GeometricInfo& getGeometricInfo() const; /* @ brief is true if Mask is on a homogeneous grid */ // Inhomogeneous grids are not supported at the moment, but if they will // be supported in the future the interface does not need to change. bool isGridHomogeneous() const { return true; }; IDType getMaskUID() const { return _maskUID; }; }; } } } #endif diff --git a/code/io/other/rttbDVHTxtFileReader.cpp b/code/io/other/rttbDVHTxtFileReader.cpp index ba3590d..366a631 100644 --- a/code/io/other/rttbDVHTxtFileReader.cpp +++ b/code/io/other/rttbDVHTxtFileReader.cpp @@ -1,226 +1,226 @@ // ----------------------------------------------------------------------- // 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; - int numberOfBins; + 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 = std::max(dataDifferential.size(), dataCumulative.size()); + 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 namepsace other + }//end namespace other }//end namespace io }//end namespace rttb diff --git a/code/io/other/rttbDVHTxtFileWriter.cpp b/code/io/other/rttbDVHTxtFileWriter.cpp index ac6f8cd..87651bd 100644 --- a/code/io/other/rttbDVHTxtFileWriter.cpp +++ b/code/io/other/rttbDVHTxtFileWriter.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) */ #include #include #include "rttbDVHTxtFileWriter.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbException.h" namespace rttb { namespace io { namespace other { DVHTxtFileWriter::DVHTxtFileWriter(FileNameString aFileName, DVHType aDVHType) { this->setFileName(aFileName); this->setDVHType(aDVHType); } void DVHTxtFileWriter::setDVHType(DVHType aDVHType) { _dvhType = aDVHType; } FileNameString DVHTxtFileWriter::getFileName() const { return _fileName; } void DVHTxtFileWriter::setFileName(FileNameString aFileName) { _fileName = aFileName; } DVHType DVHTxtFileWriter::getDVHType() const { return _dvhType; } void DVHTxtFileWriter::writeDVH(DVHPointer aDvh) { if (!aDvh) { throw core::NullPointerException("aDvh must not be NULL! "); } std::ofstream out_dvh_ofstream(this->_fileName.c_str(), std::ios::out); if (!out_dvh_ofstream.is_open() || !out_dvh_ofstream.good()) { throw core::InvalidParameterException("Invalid dvh file name: could not open write file"); } else { //setting string stream precission explicitly is mandatory to guarantee that tests //run sucessfully on different systems! out_dvh_ofstream.precision(10); if (_dvhType.Type != DVHType::Differential && _dvhType.Type != DVHType::Cumulative) { throw core::InvalidParameterException("DVH Type not acceptable: Only: DIFFERENTIAL/CUMULATIVE!"); } if (_dvhType.Type == DVHType::Differential) { out_dvh_ofstream << "DVH Type: DIFFERENTIAL\n"; } else if (_dvhType.Type == DVHType::Cumulative) { out_dvh_ofstream << "DVH Type: CUMULATIVE\n"; } out_dvh_ofstream << "DeltaD: " << aDvh->getDeltaD() << "\n"; out_dvh_ofstream << "DeltaV: " << aDvh->getDeltaV() << "\n"; out_dvh_ofstream << "StructureID: " << aDvh->getStructureID() << "\n"; out_dvh_ofstream << "DoseID: " << aDvh->getDoseID() << "\n"; out_dvh_ofstream << "DVH Data: " << "\n"; if (_dvhType.Type == DVHType::Differential) { DataDifferentialType dataDifferential = aDvh->getDataDifferential(); - int numberOfBins = dataDifferential.size(); + unsigned int numberOfBins = static_cast(dataDifferential.size()); - for (int i = 0; i < numberOfBins; i++) + for (unsigned int i = 0; i < numberOfBins; i++) { out_dvh_ofstream << i << "," << dataDifferential[i] << "\n"; } } else if (_dvhType.Type == DVHType::Cumulative) { DataDifferentialType dataCumulative = aDvh->calcCumulativeDVH(); - int numberOfBins = dataCumulative.size(); + unsigned int numberOfBins = static_cast(dataCumulative.size()); - for (int i = 0; i < numberOfBins; i++) + for (unsigned int i = 0; i < numberOfBins; i++) { out_dvh_ofstream << i << "," << dataCumulative[i] << "\n"; } } } } } } } \ No newline at end of file diff --git a/code/io/other/rttbDVHXMLFileReader.cpp b/code/io/other/rttbDVHXMLFileReader.cpp index e9e8126..020b799 100644 --- a/code/io/other/rttbDVHXMLFileReader.cpp +++ b/code/io/other/rttbDVHXMLFileReader.cpp @@ -1,148 +1,147 @@ // ----------------------------------------------------------------------- // 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 #include "rttbDVHXMLFileReader.h" #include "rttbInvalidParameterException.h" #include "rttbBaseType.h" namespace rttb { namespace io { namespace other { DVHXMLFileReader::DVHXMLFileReader(FileNameString aFileName) { _fileName = aFileName; _resetFile = true; } DVHXMLFileReader::~DVHXMLFileReader() {} void DVHXMLFileReader::resetFileName(FileNameString aFileName) { _fileName = aFileName; _resetFile = true; } void DVHXMLFileReader::createDVH() { using boost::property_tree::ptree; ptree pt; // Load the XML file into the property tree. If reading fails // (cannot open file, parse error), an exception is thrown. try { read_xml(_fileName, pt); } catch (boost::property_tree::xml_parser_error& /*e*/) { throw rttb::core::InvalidParameterException("DVH file name invalid: could not open the xml file!"); } std::string dvhType; - int numberOfBins; std::deque dataDifferential, dataCumulative; DoseTypeGy deltaD = 0; DoseVoxelVolumeType deltaV = 0; IDType strID; IDType doseID; dvhType = pt.get("dvh.type"); deltaD = pt.get("dvh.deltaD"); deltaV = pt.get("dvh.deltaV"); strID = pt.get("dvh.structureID"); doseID = pt.get("dvh.doseID"); if (dvhType != "DIFFERENTIAL" && dvhType != "CUMULATIVE") { throw core::InvalidParameterException("DVH Type invalid! Only: DIFFERENTIAL/CUMULATIVE!"); } int count = 0; BOOST_FOREACH(boost::property_tree::ptree::value_type & v, pt.get_child("dvh.data")) { if (count % 2 == 1) { if (dvhType == "DIFFERENTIAL") { dataDifferential.push_back(boost::lexical_cast(v.second.data())); } else if (dvhType == "CUMULATIVE") { dataCumulative.push_back(boost::lexical_cast(v.second.data())); } } count++; } - numberOfBins = std::max(dataDifferential.size(), dataCumulative.size()); + unsigned int numberOfBins = static_cast(std::max(dataDifferential.size(), dataCumulative.size())); if (dvhType == "CUMULATIVE") //dataDifferential should be calculated { 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); } } _dvh = boost::make_shared(dataDifferential, deltaD, deltaV, strID, doseID); _resetFile = false; } DVHXMLFileReader::DVHPointer DVHXMLFileReader::generateDVH() { if (_resetFile) { this->createDVH(); } return _dvh; } - }//end namepsace other + }//end namespace other }//end namespace io }//end namespace rttb diff --git a/code/io/virtuos/rttbVirtuosFileStructureSetGenerator.cpp b/code/io/virtuos/rttbVirtuosFileStructureSetGenerator.cpp index 6dcecd9..9aab674 100644 --- a/code/io/virtuos/rttbVirtuosFileStructureSetGenerator.cpp +++ b/code/io/virtuos/rttbVirtuosFileStructureSetGenerator.cpp @@ -1,241 +1,241 @@ // ----------------------------------------------------------------------- // 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 "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbStructure.h" #include "rttbVirtuosFileStructureSetGenerator.h" namespace rttb { namespace io { namespace virtuos { VirtuosFileStructureSetGenerator::VirtuosFileStructureSetGenerator(FileNameType aVirtuosVDXFileName, FileNameType aVirtuosCTXFileName) : _pPointerOnVirtuosCube(new Cubeinfo*), _patient(NULL) { _VDXFileName = aVirtuosVDXFileName; _CTXFileName = aVirtuosCTXFileName; //check if file names are valid if (!boost::filesystem::exists(_VDXFileName)) { throw core::InvalidParameterException("invalid vdx file name"); } if (!boost::filesystem::exists(_CTXFileName)) { throw core::InvalidParameterException("invalid ctx file name"); } *_pPointerOnVirtuosCube = create_cubeinfo(0); this->initializeVirtuosCube(_CTXFileName); } VirtuosFileStructureSetGenerator::~VirtuosFileStructureSetGenerator() { if (this->_patient != NULL) { delete this->_patient; } closecube((*(this->_pPointerOnVirtuosCube))); nc_free_cubeinfo((*(this->_pPointerOnVirtuosCube))); delete this->_pPointerOnVirtuosCube; } void VirtuosFileStructureSetGenerator::importStructureSet(FileNameType aVirtuosVDXFileName, FileNameType aVirtuosCTXFileName) { //check file name if (aVirtuosCTXFileName.empty() || aVirtuosVDXFileName.empty()) { throw core::InvalidParameterException("Virtuos VDX/CTX file name must not be empty!"); } - int vdxPosition = aVirtuosVDXFileName.find(".vdx"); + size_t vdxPosition = aVirtuosVDXFileName.find(".vdx"); if (vdxPosition == std::string::npos) { throw core::InvalidParameterException("Virtuos VDX file name must be *.vdx!"); } //get patientFileName, patientDataPath for Virtuos function std::string patientFileName, patientName, patientDataPath; patientFileName.assign(aVirtuosVDXFileName, aVirtuosVDXFileName.find_last_of("/") + 1, aVirtuosVDXFileName.length()); patientName.assign(patientFileName, 0, patientFileName.find_first_of(".")); patientDataPath.assign(aVirtuosVDXFileName, 0, aVirtuosVDXFileName.find_last_of("/") + 1); //Virtuos: voi create voi model int errorcode = voi_create_voi_model_dirolab(patientName.c_str(), patientDataPath.c_str(), 0, this->_patient); if (errorcode != 0) { //throw std::string ("Virtuos Routines unable to create VOI Model! "); std::cerr << "voi_create_voi_model_dirolab error: error code " << errorcode << std::endl; } //Virtuos: voi read vdx errorcode = voi_read_vdx_version_2_for_DIROlab(patientFileName.c_str(), patientDataPath.c_str(), this->_patient); if (errorcode != 0) { //throw std::string ("voi_read_vdx_version_2_for_DIROlab failed! "); std::cerr << "voi_read_vdx_version_2_for_DIROlab error: error code " << errorcode << std::endl; } int numberOfVois = this->_patient->getNumberOfVois(); float firstSliceInFrame = (*_pPointerOnVirtuosCube)->pos_list[0].position; - double sliceThickness = (*_pPointerOnVirtuosCube)->slicedist; + float sliceThickness = (*_pPointerOnVirtuosCube)->slicedist; - float lastSliceInFrame = ((*_pPointerOnVirtuosCube)->dimz - 1) * sliceThickness + firstSliceInFrame; + float lastSliceInFrame = static_cast(((*_pPointerOnVirtuosCube)->dimz - 1) * sliceThickness + firstSliceInFrame); for (int currentVoiNumber = 0; currentVoiNumber < numberOfVois; currentVoiNumber++) { std::string voiName = ""; char tmpVoiName[1024]; voi_get_voi_name_by_index_dirolab(currentVoiNumber, 1024, tmpVoiName, this->_patient); voiName.assign(tmpVoiName); /* prepare contour extraction */ D3PointList* contours = NULL; contours = d3_list_create(1000000); D3Point origin = {0, 0, 0}, y_axis_point = {0, 0, 0}, x_axis_point = {0, 0, 0}; int maxNumberOfContours = 100000; int* pNoOFContours = &maxNumberOfContours; PolygonSequenceType polygonSequence; for (float z = firstSliceInFrame; z <= lastSliceInFrame; z += sliceThickness) { origin.x = 0.0; origin.y = 0.0; origin.z = z ; x_axis_point.x = 1.0; x_axis_point.y = 0.0; x_axis_point.z = z ; y_axis_point.x = 0.0; y_axis_point.y = 1.0; y_axis_point.z = z; *pNoOFContours = 100000; //<-- reason is the next function call voi_get_CT_contours_dirolab(voiName.c_str(), origin, x_axis_point, y_axis_point, pNoOFContours, &contours, 1, this->_patient); for (int i = 0; i < *pNoOFContours; i++) { PolygonType polygon; for (int j = 0; j < contours[i].used - 1; j++) //Virtuos polygon the last point is the same as the first point { WorldCoordinate3D point; point(0) = contours[i].points[j].x; point(1) = contours[i].points[j].y; point(2) = contours[i].points[j].z; polygon.push_back(point); }//end for j polygonSequence.push_back(polygon); }//end for i }//end for z boost::shared_ptr spStruct = boost::make_shared(polygonSequence); spStruct->setLabel(voiName); _strVector.push_back(spStruct); }//end for currentVoiNumber } void VirtuosFileStructureSetGenerator::initializeVirtuosCube(FileNameType aVirtuosCTXFileName) { if (aVirtuosCTXFileName.empty()) { throw core::InvalidParameterException("Empty File Name"); } - int gzPosition = aVirtuosCTXFileName.find(".gz"); + size_t gzPosition = aVirtuosCTXFileName.find(".gz"); if (gzPosition != std::string::npos) { aVirtuosCTXFileName.erase(gzPosition, aVirtuosCTXFileName.length()); } nc_init_cubeinfo(*_pPointerOnVirtuosCube); int opencubeErrorCode = opencube(aVirtuosCTXFileName.c_str() , *_pPointerOnVirtuosCube); if (opencubeErrorCode != 0) { std::stringstream opencubeErrorCodeAsStringStream; opencubeErrorCodeAsStringStream << opencubeErrorCode; throw core::InvalidParameterException(std::string( std::string("VirtuosIO::openCube returned error Code: ") + opencubeErrorCodeAsStringStream.str())); } if ((*_pPointerOnVirtuosCube)->dimx == 0 || (*_pPointerOnVirtuosCube)->dimy == 0 || (*_pPointerOnVirtuosCube)->dimz == 0) { throw core::InvalidParameterException("Invalid Cube dimension: dimX/dimY/dimZ must not be zero! "); } } VirtuosFileStructureSetGenerator::StructureSetPointer VirtuosFileStructureSetGenerator::generateStructureSet() { this->importStructureSet(_VDXFileName, _CTXFileName); return boost::make_shared(_strVector); } }//end namespace virtuos }//end namespace io }//end namespace rttb diff --git a/code/masks/boost/rttbBoostMask.cpp b/code/masks/boost/rttbBoostMask.cpp index a9b44a5..186b08f 100644 --- a/code/masks/boost/rttbBoostMask.cpp +++ b/code/masks/boost/rttbBoostMask.cpp @@ -1,490 +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 < _geometricInfo->getNumSlices(); indexZ++) + 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 < _geometricInfo->getNumSlices() - 1; indexZ++) + 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 < _geometricInfo->getNumSlices(); indexZ++) + 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; - for (unsigned int indexZ = 0; indexZ < _geometricInfo->getNumSlices(); indexZ++) + 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) { - if (sliceNumber < 0 || sliceNumber >= _geometricInfo->getNumSlices()) + 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++) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } //if not jet, check now if (!it1IsDonut) { bool it2IsDonut = false; unsigned int index2 = 0; for (it2 = aRingVector.begin(); it2 != aRingVector.end(); it2++, index2++) { if (it2 != it1) { BoostMask::BoostPolygon2D polygon2D; if (::boost::geometry::within(*it1, *it2)) { ::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring it2IsDonut = true; } //if donut else if (::boost::geometry::within(*it2, *it1)) { ::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring it2IsDonut = true; } if (it2IsDonut) { donutIndexVector.push_back(index1); donutIndexVector.push_back(index2); donutVector.push_back(polygon2D);//store donut polygon break;//Only store the first donut! } } } } } //Store no donut polygon to boostPolygonVector index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); it1++, index1++) { bool it1IsDonut = false; //check if the ring is the outer or inner of a donut for (unsigned int i = 0; i < donutIndexVector.size(); i++) { if (donutIndexVector.at(i) == index1) { it1IsDonut = true; break; } } if (!it1IsDonut) { BoostMask::BoostPolygon2D polygon2D; ::boost::geometry::append(polygon2D, *it1); boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut } } //Append donut polygon to boostPolygonVector BoostMask::BoostPolygonVector::iterator itDonut; for (itDonut = donutVector.begin(); itDonut != donutVector.end(); itDonut++) { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } } } } \ No newline at end of file diff --git a/code/models/rttbBioModelScatterPlots.cpp b/code/models/rttbBioModelScatterPlots.cpp index 0918a14..e3430e2 100644 --- a/code/models/rttbBioModelScatterPlots.cpp +++ b/code/models/rttbBioModelScatterPlots.cpp @@ -1,219 +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 "rttbBioModelScatterPlots.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace models { /* Initiate Random Number generator with current time */ boost::random::mt19937 rng(static_cast(time(0))); /* Generate random number between 0 and 1 */ boost::random::uniform_01 uniDist(rng); ScatterPlotType getScatterPlotVary1Parameter(BioModel& aModel, int aParamId, BioModelParamType aMean, BioModelParamType aVariance, DoseTypeGy aNormalisationDose, int numberOfPoints, DoseTypeGy aMinDose, DoseTypeGy aMaxDose) { ScatterPlotType scatterPlotData; if (aVariance == 0) { //set to small positive value to avoid negative infinity! aVariance = 1e-30; } if (aMaxDose <= aMinDose) { throw core::InvalidParameterException("Parameter invalid: aMaxDose must be > aMinDose!"); } if (aNormalisationDose <= 0) { throw core::InvalidParameterException("Parameter invalid: aNormalisationDose must be > 0!"); } /* Choose Normal Distribution */ boost::random::normal_distribution gaussian_dist(0, aVariance); /* Create a Gaussian Random Number generator * by binding with previously defined * normal distribution object */ boost::random::variate_generator > generator( rng, gaussian_dist); int i = 0; while (i < numberOfPoints) { double paramValue, probability; double randomValue = generator(); paramValue = randomValue + aMean; probability = normal_pdf(randomValue, aVariance); if (probability > 1) { probability = 1; } //randomly select a dose between aMinDose and aMaxDose double dose = uniDist() * (aMaxDose - aMinDose) + aMinDose; if (probability > 0) { try { aModel.setParameterByID(aParamId, paramValue); aModel.init(dose / aNormalisationDose); double value = aModel.getValue(); std::pair modelProbPair = std::make_pair(value, probability); scatterPlotData.insert(std::pair >(dose, modelProbPair)); i++; } catch (core::InvalidParameterException& /*e*/) { //repeat evaluation to guarantee the correct number of scatter values continue; } } } return scatterPlotData; } double normal_pdf(double aValue, double aVariance) { - static const float inv_sqrt_2pi = 0.3989422804014327; + static const double inv_sqrt_2pi = 0.3989422804014327; double a = (aValue) / aVariance; return inv_sqrt_2pi / aVariance * std::exp(-0.5f * a * a); } ScatterPlotType getScatterPlotVaryParameters(BioModel& aModel, std::vector aParamIdVec, BioModel::ParamVectorType aMeanVec, BioModel::ParamVectorType aVarianceVec, DoseTypeGy aNormalisationDose, int numberOfPoints, DoseTypeGy aMinDose, DoseTypeGy aMaxDose) { ScatterPlotType scatterPlotData; if (aMaxDose <= aMinDose) { throw core::InvalidParameterException("Parameter invalid: aMaxDose must be > aMinDose!"); } if (aNormalisationDose <= 0) { throw core::InvalidParameterException("Parameter invalid: aNormalisationDose must be > 0!"); } //all input vectors need to have the same size if (((aVarianceVec.size() != aMeanVec.size()) || (aVarianceVec.size() != aParamIdVec.size()))) { throw core::InvalidParameterException("Parameter vectors have different sizes!"); } for (GridIndexType v = 0; v < aVarianceVec.size(); v++) { //set to small positive value to avoid negative infinity! if (aVarianceVec.at(v) == 0) { aVarianceVec.at(v) = 1e-30; } } double paramValue; // vary all parameters for each scattered point int i = 0; while (i < numberOfPoints) { double probability = 1; for (GridIndexType j = 0; j < aParamIdVec.size(); j++) { /* Choose Normal Distribution */ boost::random::normal_distribution gaussian_dist(0, aVarianceVec.at(j)); /* Create a Gaussian Random Number generator * by binding with previously defined * normal distribution object */ boost::random::variate_generator > generator( rng, gaussian_dist); double randomValue = generator(); paramValue = randomValue + aMeanVec.at(j); if (aVarianceVec.at(j) != 0) { /* calculate combined probability */ probability = probability * normal_pdf(randomValue, aVarianceVec.at(j)); } else { throw core::InvalidParameterException("Parameter invalid: Variance should not be 0!"); } aModel.setParameterByID(aParamIdVec.at(j), paramValue); } //randomly select a dose between aMinDose and aMaxDose double dose = uniDist() * (aMaxDose - aMinDose) + aMinDose; if (probability > 0) { try { aModel.init(dose / aNormalisationDose); double value = aModel.getValue(); std::pair modelProbPair = std::make_pair(value, probability); scatterPlotData.insert(std::pair >(dose, modelProbPair)); i++; } catch (core::InvalidParameterException& /*e*/) { //repeat evaluation to guarantee the correct number of scatter values continue; } } } return scatterPlotData; } }//end namespace models }//end namespace rttb diff --git a/code/models/rttbIntegration.h b/code/models/rttbIntegration.h index d9dc45b..54e2b49 100644 --- a/code/models/rttbIntegration.h +++ b/code/models/rttbIntegration.h @@ -1,117 +1,117 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __INTEGRATION_H #define __INTEGRATION_H #include "rttbBaseTypeModels.h" namespace rttb { namespace models { - typedef float integrationType; + typedef double integrationType; /*! @class LkbModelFunctor @brief A FunctorType: calculate the transformed LKB-Model using the transformation x = b - (1-t)/t. \int_{-\infty}^{b} dx f(x) = \int_0^1 dt f(b - (1-t)/t)/t^2 */ class LkbModelFunctor { public: /*!b: upper bound of the lkb integration*/ double b; /*! calculate the transformed LKB-Model using the transformation x = b - (1-t)/t.. */ double calculate(double x) const; }; /*! @class tcpModelFunctor @brief .A FunctorType: calculate the transformed TCP-LQ-Model using the transformation x = a + (1-t)/t. \int_{a}^{+\infty} dx f(x) =\int_0^1 dt f(a + (1-t)/t)/t^2 */ class tcpModelFunctor { public: //TCP parameters including alpha mean, alpha variance, rho, volume vector and bed vector TcpParams params; /*! a: lower bound of the tcp integration*/ double a; /*! calculate the transformed TCP-LQ-Model using the transformation x = a + (1-t)/t. */ double calculate(double x) const; }; /*! Function to be integrated for TCP LQ model. @param x: variable of the TCP LQ model @param params: TCP parameters including alpha mean, alpha variance, rho, volume vector and bed vector @return Return the function value */ double tcpFunction(double x, const TcpParams& tcp_params); /*! Compute integration step for f(x) = exp(-pow(t,2)/2). @param x: variable of the lkb function @return Return the function value */ double lkbFunction(double x); /*! Calculate LKB Integration over(-infinity,b). The integral is mapped onto the semi-open interval (0,1] using the transformation x = b - (1-t)/t @param b: upper bound of the lkb integration */ double integrateLKB(double b); /*! Calculate TCP integration over (a, infinity). The integral is mapped onto the semi-open interval (0,1] using the transformation x = a + (1-t)/t @param a: lower bound of the tcp integration @param params: TCP parameters including alpha mean, alpha variance, rho, volume vector and bed vector */ double integrateTCP(double a, const TcpParams& params); /*This function returns the nth stage of refinement of the extended trapezoidal rule. @param BMfunction: function to be integrated, for example a LkbModelFunctor or a tcpModelFunctor @param a: lower bound of the integral @param b: upper bound of the integral @param stepNum: the nth stage @param result: the current result */ template integrationType trapzd(const FunctorType& BMfunction, integrationType a, integrationType b, int stepNum); /*! Iterative integration routine @param BMfunction: function to be integrated, for example a LkbModelFunctor or a tcpModelFunctor @param a: lower bound of the integral @param b: upper bound of the integral @exception throw InvalidParameterException if integral calculation failed. */ template integrationType iterativeIntegration(const FunctorType& BMfunction, integrationType a, integrationType b); } } #endif \ No newline at end of file diff --git a/testing/io/dicom/DicomDoseAccessorConverterTest.cpp b/testing/io/dicom/DicomDoseAccessorConverterTest.cpp index 8907aec..833f982 100644 --- a/testing/io/dicom/DicomDoseAccessorConverterTest.cpp +++ b/testing/io/dicom/DicomDoseAccessorConverterTest.cpp @@ -1,158 +1,158 @@ // ----------------------------------------------------------------------- // 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 "rttbDicomDoseAccessor.h" #include "rttbDicomFileDoseAccessorGenerator.h" #include "rttbDicomFileDoseAccessorWriter.h" #include "rttbInvalidDoseException.h" #include "rttbDcmrtException.h" #include "rttbInvalidParameterException.h" namespace rttb { namespace testing { /*!@brief DicomDoseAccessorGeneratorTest - test the generators for dicom data 1) test dicom file generator generateDoseAccessor() 2) test dicom IOD generator generateDoseAccessor() */ int DicomDoseAccessorConverterTest(int argc, char* argv[]) { typedef boost::shared_ptr DRTDoseIODPointer; typedef rttb::io::dicom::DicomDoseAccessor::DoseAccessorPointer DoseAccessorPointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: // 1: the file name of the dose to read // 2: the file name of the dose to write std::string RTDOSE_FILENAME_R; std::string RTDOSE_FILENAME_W; if (argc > 2) { RTDOSE_FILENAME_R = argv[1]; RTDOSE_FILENAME_W = argv[2]; } double errorConstant = 1e-3; DoseAccessorPointer doseAccessor_r = io::dicom::DicomFileDoseAccessorGenerator( RTDOSE_FILENAME_R.c_str()).generateDoseAccessor(); CHECK_NO_THROW(io::dicom::DicomFileDoseAccessorWriter()); io::dicom::DicomFileDoseAccessorWriter doseConverter; CHECK_NO_THROW(doseConverter.setDoseAccessor(doseAccessor_r)); CHECK_NO_THROW(doseConverter.setFileName(RTDOSE_FILENAME_W)); CHECK_EQUAL(doseConverter.process(), true); DoseAccessorPointer doseAccessor_w = io::dicom::DicomFileDoseAccessorGenerator( RTDOSE_FILENAME_W).generateDoseAccessor(); //Check geometricinfo CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImagePositionPatient().x(), doseAccessor_w->getGeometricInfo().getImagePositionPatient().x(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImagePositionPatient().y(), doseAccessor_w->getGeometricInfo().getImagePositionPatient().y(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImagePositionPatient().z(), doseAccessor_w->getGeometricInfo().getImagePositionPatient().z(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImageOrientationColumn().x(), doseAccessor_w->getGeometricInfo().getImageOrientationColumn().x(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImageOrientationColumn().y(), doseAccessor_w->getGeometricInfo().getImageOrientationColumn().y(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImageOrientationColumn().z(), doseAccessor_w->getGeometricInfo().getImageOrientationColumn().z(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImageOrientationRow().x(), doseAccessor_w->getGeometricInfo().getImageOrientationRow().x(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImageOrientationRow().y(), doseAccessor_w->getGeometricInfo().getImageOrientationRow().y(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getImageOrientationRow().z(), doseAccessor_w->getGeometricInfo().getImageOrientationRow().z(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getSpacing().x(), doseAccessor_w->getGeometricInfo().getSpacing().x(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getSpacing().y(), doseAccessor_w->getGeometricInfo().getSpacing().y(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getSpacing().z(), doseAccessor_w->getGeometricInfo().getSpacing().z(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getNumColumns(), doseAccessor_w->getGeometricInfo().getNumColumns(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getNumRows(), doseAccessor_w->getGeometricInfo().getNumRows(), errorConstant); CHECK_CLOSE(doseAccessor_r->getGeometricInfo().getNumSlices(), doseAccessor_w->getGeometricInfo().getNumSlices(), errorConstant); //Check pixel data - int size = doseAccessor_r->getGeometricInfo().getNumColumns() * - doseAccessor_r->getGeometricInfo().getNumRows() * - doseAccessor_r->getGeometricInfo().getNumSlices() ; + unsigned int size = doseAccessor_r->getGeometricInfo().getNumColumns() * + doseAccessor_r->getGeometricInfo().getNumRows() * + doseAccessor_r->getGeometricInfo().getNumSlices() ; for (unsigned int index = 0; index < 30; index++) { CHECK_CLOSE(doseAccessor_r->getValueAt(index), doseAccessor_w->getValueAt(index), errorConstant); if (size / 2 - index >= 0 && size / 2 - index < size) { CHECK_CLOSE(doseAccessor_r->getValueAt(size / 2 - index), doseAccessor_w->getValueAt(size / 2 - index), errorConstant); } CHECK_CLOSE(doseAccessor_r->getValueAt(size - index - 1), doseAccessor_w->getValueAt(size - index - 1), errorConstant); } RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/io/rttbDoseAccessorTester.cpp b/testing/io/rttbDoseAccessorTester.cpp index 5085f7d..3b05604 100644 --- a/testing/io/rttbDoseAccessorTester.cpp +++ b/testing/io/rttbDoseAccessorTester.cpp @@ -1,152 +1,152 @@ // ----------------------------------------------------------------------- // 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 "rttbGeometricInfo.h" #include "rttbDoseAccessorTester.h" namespace rttb { namespace testing { DoseAccessorTester::DoseAccessorTester(DoseAccessorPointer aReferenceDose, DoseAccessorPointer aCompareDose) : _referenceDose(aReferenceDose), _compareDose(aCompareDose), _geometryIsSimilar(false), _sameGridSize(false), _conversionFailed(false), _numDifference(0) {} void DoseAccessorTester::setReferenceDose(const DoseAccessorPointer aReferenceDose) { _referenceDose = aReferenceDose; } void DoseAccessorTester::setCompareDose(const DoseAccessorPointer aCompareDose) { _compareDose = aCompareDose; } lit::StringType DoseAccessorTester::getTestDescription(void) const { return "Compare two DoseAccessors and determine if the contained doses are equal."; }; bool DoseAccessorTester::doCheck(void) const { _pResults->onTestStart(getCurrentTestLabel()); _geometryIsSimilar = (_referenceDose->getGeometricInfo() == _compareDose->getGeometricInfo()); if (!_geometryIsSimilar) { return false; } _sameGridSize = (_referenceDose->getGridSize() == _compareDose->getGridSize()); if (!_sameGridSize) { return false; } _numDifference = 0; _maxDifference = 0; VoxelGridIndex3D id3D; for (VoxelGridID id = 0; id < _referenceDose->getGridSize(); id++) { _compareDose->getGeometricInfo().convert(id, id3D); if (!_compareDose->getGeometricInfo().validIndex(id3D)) { _conversionFailed = true; _failedID = id; return false; } if ((_referenceDose-> getValueAt(id) != _referenceDose-> getValueAt(id3D)) || (_compareDose->getValueAt(id) != _compareDose-> getValueAt(id3D))) { _conversionFailed = true; _failedID = id; return false; } if ((_referenceDose->getValueAt(id) != _compareDose->getValueAt(id)) || (_referenceDose->getValueAt(id3D) != _compareDose->getValueAt(id3D))) { - float tmpDifference = abs(_referenceDose->getValueAt(id) - _compareDose->getValueAt(id)); + double tmpDifference = abs(_referenceDose->getValueAt(id) - _compareDose->getValueAt(id)); if (tmpDifference > _maxDifference) { _maxDifference = tmpDifference; } _numDifference++; } }//end for(VoxelGridID id = 0; id < _referenceDose->getGridSize(); id++) return (_numDifference == 0); } void DoseAccessorTester::handleSuccess(void) const { std::ostringstream stream; stream << "Both doses are equal" << std::endl; _pResults->onTestSuccess(getCurrentTestLabel(), stream.str()); } void DoseAccessorTester::handleFailure(void) const { std::ostringstream stream; stream << "Doses are different" << std::endl; if (_geometryIsSimilar) { if (_sameGridSize) { stream << std::endl << "Error voxel count: " << _numDifference << std::endl; stream << std::endl << "Maximum difference: " << _maxDifference << std::endl; if (_conversionFailed) { stream << std::endl << "Index conversion failed in: " << _failedID << std::endl; } } else { stream << "Doses have different grid sizes" << std::endl; stream << "Reference dose contains " << _referenceDose->getGridSize() << " voxels and comparison dose " << _compareDose->getGridSize() << std::endl; } } else { stream << "Doses have different geometry" << std::endl; } _pResults->onTestFailure(getCurrentTestLabel(), stream.str()); } } } \ No newline at end of file diff --git a/testing/io/rttbDoseAccessorTester.h b/testing/io/rttbDoseAccessorTester.h index c77fae0..80e8e40 100644 --- a/testing/io/rttbDoseAccessorTester.h +++ b/testing/io/rttbDoseAccessorTester.h @@ -1,97 +1,97 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __DOSE_ACCESSOR_TESTER_H #define __DOSE_ACCESSOR_TESTER_H #include #include "litTesterBase.h" #include "litString.h" #include "rttbBaseType.h" #include "../../code/core/rttbDoseAccessorInterface.h" namespace rttb { namespace testing { /*! class DoseAccessorTester @brief Tester class for io classes. Compares two given DoseAccessors for similarity. These DoseAccessors are not similar if their geometry or grid size are not similar, - if the conversion of a given ID is invalid for one of the acessors, or if the dose + if the conversion of a given ID is invalid for one of the accessors, or if the dose at a given ID is not the same for both accessors. */ class DoseAccessorTester: public lit::TesterBase { public: typedef core::DoseAccessorInterface::DoseAccessorPointer DoseAccessorPointer; private: DoseAccessorPointer _referenceDose; DoseAccessorPointer _compareDose; - mutable float _maxDifference; - mutable float _numDifference; + mutable double _maxDifference; + mutable double _numDifference; mutable bool _geometryIsSimilar; mutable bool _sameGridSize; mutable bool _conversionFailed; mutable VoxelGridID _failedID; public: DoseAccessorTester(DoseAccessorPointer aReferenceDose, DoseAccessorPointer aCompareDose); /*! Set the dose accessor pointer for the dose comparison. */ void setReferenceDose(const DoseAccessorPointer aReferenceDose); void setCompareDose(const DoseAccessorPointer aCompareDose); /*! Returns a string that specifies the test the tester currently performs. */ lit::StringType getTestDescription(void) const; lit::StringType getTestName(void) const { return "DoseAccessorTester"; }; protected: /*! performs the test and checks the results. @result Indicates if the test was successful (true) or if it failed (false) */ bool doCheck(void) const; /*! Function will be called by check() if test was succesfull. */ void handleSuccess(void) const; /*! Function will be called by check() if test was a failure. */ void handleFailure(void) const; }; } } #endif \ No newline at end of file diff --git a/testing/io/virtuos/VirtuosStructureSetGeneratorTest.cpp b/testing/io/virtuos/VirtuosStructureSetGeneratorTest.cpp index 011476f..b539af4 100644 --- a/testing/io/virtuos/VirtuosStructureSetGeneratorTest.cpp +++ b/testing/io/virtuos/VirtuosStructureSetGeneratorTest.cpp @@ -1,119 +1,119 @@ // ----------------------------------------------------------------------- // 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 "rttbGeometricInfo.h" #include "rttbDoseIteratorInterface.h" #include "rttbVirtuosDoseAccessor.h" #include "rttbVirtuosFileStructureSetGenerator.h" #include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" #include "rttbIndexOutOfBoundsException.h" namespace rttb { namespace testing { /*!@brief VirtuosStructureSetGeneratorTest - test structure set generator for virtuos data 1) test constructor 2) test generateStructureSet */ int VirtuosStructureSetGeneratorTest(int argc, char* argv[]) { typedef core::StructureSetGeneratorInterface::StructureSetPointer StructureSetPointer; PREPARE_DEFAULT_TEST_REPORTING; //ARGUMENTS: 1: structure file name ".../testing/data/Virtuos/NHH030000.vdx" // 2: ctx file name (virtuos) ".../testing/data/Virtuos/NHH030000.ctx.gz" std::string RTSTRUCT_FILENAME; std::string CT_FILENAME; if (argc > 1) { RTSTRUCT_FILENAME = argv[1]; } if (argc > 2) { CT_FILENAME = argv[2]; } /* structure set */ //1) check constructor CHECK_THROW_EXPLICIT(io::virtuos::VirtuosFileStructureSetGenerator("aStructureFileName.vdx", CT_FILENAME.c_str()), core::InvalidParameterException); CHECK_THROW_EXPLICIT(io::virtuos::VirtuosFileStructureSetGenerator(RTSTRUCT_FILENAME.c_str(), "aCTFileName.ctx"), core::InvalidParameterException); CHECK_NO_THROW(io::virtuos::VirtuosFileStructureSetGenerator(RTSTRUCT_FILENAME.c_str(), CT_FILENAME.c_str())); //2) check generateStructureSet CHECK_THROW_EXPLICIT(io::virtuos::VirtuosFileStructureSetGenerator("aStructureFileName.vdx", CT_FILENAME.c_str()).generateStructureSet(), core::InvalidParameterException); CHECK_THROW_EXPLICIT(io::virtuos::VirtuosFileStructureSetGenerator(RTSTRUCT_FILENAME.c_str(), "aCTFileName.ctx").generateStructureSet(), core::InvalidParameterException); std::string testFileName = RTSTRUCT_FILENAME; - int gzPosition = testFileName.find(".vdx"); + size_t gzPosition = testFileName.find(".vdx"); if (gzPosition != std::string::npos) { testFileName.erase(gzPosition, testFileName.length()); } CHECK_THROW_EXPLICIT(io::virtuos::VirtuosFileStructureSetGenerator(testFileName, CT_FILENAME.c_str()).generateStructureSet(), core::InvalidParameterException); testFileName = CT_FILENAME; gzPosition = testFileName.find(".ctx"); if (gzPosition != std::string::npos) { testFileName.erase(gzPosition, testFileName.length()); } CHECK_THROW_EXPLICIT(io::virtuos::VirtuosFileStructureSetGenerator(RTSTRUCT_FILENAME.c_str(), testFileName).generateStructureSet(), core::InvalidParameterException); CHECK_NO_THROW(io::virtuos::VirtuosFileStructureSetGenerator(RTSTRUCT_FILENAME.c_str(), CT_FILENAME.c_str()).generateStructureSet()); RETURN_AND_REPORT_TEST_SUCCESS; } }//testing }//rttb diff --git a/testing/masks/rttbMaskRectStructTester.h b/testing/masks/rttbMaskRectStructTester.h index 60b2ceb..4b11ac3 100644 --- a/testing/masks/rttbMaskRectStructTester.h +++ b/testing/masks/rttbMaskRectStructTester.h @@ -1,102 +1,102 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __MASK_RECT_TESTER_H #define __MASK_RECT_TESTER_H #include #include "litTesterBase.h" #include "litString.h" #include "rttbBaseType.h" #include "../../code/core/rttbMaskAccessorInterface.h" namespace rttb { namespace testing { /*! class MaskRectStructTester @brief Tests if masked voxel are all inside the given boundaries The boundaries are defined in DummyStructure::CreateRectangularStructureCentered @see DummyStructure */ class MaskRectStructTester: public lit::TesterBase { private: typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; typedef core::MaskAccessorInterface::MaskAccessorPointer MaskAccessorPointer; MaskVoxelListPointer _referenceList; MaskAccessorPointer _maskAccessor; GridIndexType _zIndex; mutable bool _conversionFailed; mutable bool _differentSlice; mutable bool _outsideStructure; mutable bool _wrongRVF; mutable VoxelGridID _failedListIndex; - mutable float _wrongVal; + mutable double _wrongVal; public: MaskRectStructTester(MaskAccessorPointer aMaskAccessor, GridIndexType z); /*! Set the mask accessor pointer for the dose comparison. */ void setAccessor(const MaskAccessorPointer aMaskAccessor); /*! Set index of slice containing the dummy structure. */ void setZ(const GridIndexType z); /*! Returns a string that specifies the test the tester currently performs. */ lit::StringType getTestDescription(void) const; lit::StringType getTestName(void) const { return "MaskRectStructTester"; }; protected: /*! performes the test and checks the results. @result Indicates if the test was successfull (true) or if it failed (false) */ bool doCheck(void) const; /*! Function will be called be check() if test was succesfull. */ void handleSuccess(void) const; /*! Function will be called be check() if test was a failure. */ void handleFailure(void) const; }; } } #endif \ No newline at end of file diff --git a/testing/masks/rttbMaskVoxelListTester.cpp b/testing/masks/rttbMaskVoxelListTester.cpp index 6534588..be7533b 100644 --- a/testing/masks/rttbMaskVoxelListTester.cpp +++ b/testing/masks/rttbMaskVoxelListTester.cpp @@ -1,153 +1,150 @@ // ----------------------------------------------------------------------- // 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 "rttbMaskVoxelListTester.h" namespace rttb { namespace testing { MaskVoxelListTester::MaskVoxelListTester(MaskVoxelListPointer aReferenceList, MaskVoxelListPointer aCompareList) { _referenceList = aReferenceList; _compareList = aCompareList; _sameSize = false; _masVoxelsDiffer = false; _numDifference = 0; } void MaskVoxelListTester::setReferenceList(const MaskVoxelListPointer aReferenceList) { _referenceList = aReferenceList; } void MaskVoxelListTester::setCompareList(const MaskVoxelListPointer aCompareList) { _compareList = aCompareList; } lit::StringType MaskVoxelListTester::getTestDescription(void) const { return "Compare two MaskVoxelLists and determine if the contained content is equal."; }; bool MaskVoxelListTester::doCheck(void) const { _pResults->onTestStart(getCurrentTestLabel()); _masVoxelsDiffer = false; if (_referenceList->size() == _compareList->size()) { _sameSize = true; } else { _sameSize = false; return false; } _numDifference = 0; _maxDifference = 0; MaskVoxelList::iterator iterR, iterC; iterC = _compareList->begin(); VoxelGridID index = 0; for (iterR = _referenceList->begin(); iterR != _referenceList->end(); ++iterR, ++iterC, ++index) { if (iterR->getVoxelGridID() == iterC->getVoxelGridID()) { if (iterR->getRelevantVolumeFraction() == iterC->getRelevantVolumeFraction()) { continue; } else { - float diff = iterR->getRelevantVolumeFraction() - iterC->getRelevantVolumeFraction(); + double diff = iterR->getRelevantVolumeFraction() - iterC->getRelevantVolumeFraction(); if (diff > _maxDifference) { _maxDifference = diff; } - /*if(diff > 0.001){ - std::cout <getVoxelGridID()<< ": ("<< iterR->getRelevantVolumeFraction() << ","<getRelevantVolumeFraction()<<"); "; - }*/ _numDifference++; } } else { _failedListIndex = index; _masVoxelsDiffer = true; return false; } - }//end for(VoxelGridID id = 0; id < _referenceList->getGridSi... + } if (_numDifference > 0) { return false; } return true; } void MaskVoxelListTester::handleSuccess(void) const { std::ostringstream stream; stream << "Both Lists are equal" << std::endl; _pResults->onTestSuccess(getCurrentTestLabel(), stream.str()); } void MaskVoxelListTester::handleFailure(void) const { std::ostringstream stream; stream << "Lists were different" << std::endl; if (_sameSize) { stream << std::endl << "Error of volume fraction voxel count: " << _numDifference << std::endl; stream << std::endl << "Maximum difference in volume fraction: " << _maxDifference << std::endl; if (_masVoxelsDiffer) { stream << std::endl << "Mask points to different grid position in: " << _failedListIndex << std::endl; } } else { stream << "Lists have different size" << std::endl; stream << "Reference List is " << _referenceList->size() << " voxels long and comparison List " << _compareList->size() << std::endl; } _pResults->onTestFailure(getCurrentTestLabel(), stream.str()); } } } \ No newline at end of file diff --git a/testing/masks/rttbMaskVoxelListTester.h b/testing/masks/rttbMaskVoxelListTester.h index efda9e5..b07d312 100644 --- a/testing/masks/rttbMaskVoxelListTester.h +++ b/testing/masks/rttbMaskVoxelListTester.h @@ -1,95 +1,95 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __MASK_VL_TESTER_H #define __MASK_VL_TESTER_H #include #include "litTesterBase.h" #include "litString.h" #include "rttbBaseType.h" #include "../../code/core/rttbMaskAccessorInterface.h" namespace rttb { namespace testing { /*! class MaskVoxelListTester @brief Tester class for lists. Compares two given lists for similarity. Both lists need to have the same length and have equal values in each element to be considered similar. */ class MaskVoxelListTester: public lit::TesterBase { private: typedef core::MaskAccessorInterface::MaskVoxelListPointer MaskVoxelListPointer; typedef core::MaskAccessorInterface::MaskVoxelList MaskVoxelList; MaskVoxelListPointer _referenceList; MaskVoxelListPointer _compareList; - mutable float _maxDifference; - mutable float _numDifference; + mutable double _maxDifference; + mutable double _numDifference; mutable bool _sameSize; mutable bool _masVoxelsDiffer; mutable VoxelGridID _failedListIndex; public: MaskVoxelListTester(MaskVoxelListPointer aReferenceList, MaskVoxelListPointer aCompareList); /*! Set the dose accessor pointer for the dose comparison. */ void setReferenceList(const MaskVoxelListPointer aReferenceList); void setCompareList(const MaskVoxelListPointer aCompareList); /*! Returns a string that specifies the test the tester currently performs. */ lit::StringType getTestDescription(void) const; lit::StringType getTestName(void) const { return "MaskVoxelListTester"; }; protected: /*! performes the test and checks the results. @result Indicates if the test was successfull (true) or if it failed (false) */ bool doCheck(void) const; /*! Function will be called be check() if test was succesfull. */ void handleSuccess(void) const; /*! Function will be called be check() if test was a failure. */ void handleFailure(void) const; }; } } #endif \ No newline at end of file diff --git a/testing/models/rttbScatterTester.h b/testing/models/rttbScatterTester.h index 8dcef69..ba746df 100644 --- a/testing/models/rttbScatterTester.h +++ b/testing/models/rttbScatterTester.h @@ -1,111 +1,111 @@ // ----------------------------------------------------------------------- // RTToolbox - DKFZ radiotherapy quantitative evaluation library // // Copyright (c) German Cancer Research Center (DKFZ), // Software development for Integrated Diagnostics and Therapy (SIDT). // ALL RIGHTS RESERVED. // See rttbCopyright.txt or // http://www.dkfz.de/en/sidt/projects/rttb/copyright.html [^] // // This software is distributed WITHOUT ANY WARRANTY; without even // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notices for more information. // //------------------------------------------------------------------------ /*! // @file // @version $Revision$ (last changed revision) // @date $Date$ (last change date) // @author $Author$ (last changed by) */ #ifndef __SCATTER_TESTER_H #define __SCATTER_TESTER_H #include "litTesterBase.h" #include "litString.h" #include "rttbBaseType.h" #include "rttbBioModelScatterPlots.h" #include "rttbBioModelCurve.h" namespace rttb { namespace testing { /*! class ScatterTester @brief Tester class for for model scatter plots. Compares a given scatter plot with a given model curve. The values should be similar for similar doses. Variations should only depend on the variance given for the scatter plot calculation. An additional deviation of 1e-4+givenVariance is ignored. */ class ScatterTester: public lit::TesterBase { private: models::CurveDataType _referenceCurve; models::ScatterPlotType _compareScatter; - /*! Additional variance that is allowed in the comparison. This value ususally corresponds to the value used in + /*! Additional variance that is allowed in the comparison. This value usually corresponds to the value used in the generation of the scatter plot. */ models::BioModelParamType _variance; mutable float _maxDifference; mutable int _numDifference; mutable float _meanDifference; - /*! If true allows up tp 5% of the scatter points to deviate without failing. + /*! If true allows up to 5% of the scatter points to deviate without failing. If false, all points have to correspond. */ bool _allowExceptions; public: ScatterTester(models::CurveDataType aReferenceCurve, models::ScatterPlotType aCompareScatter, models::BioModelParamType aVariance = 0); /*! Set the reference curve used in comparison. */ void setReferenceCurve(const models::CurveDataType aReferenceCurve); /*! Set the scatter data used in comparison. */ void setCompareScatter(const models::ScatterPlotType aCompareScatter); /*! Set the variance that is allowed for the scatter plot. Usually this matches the parameter used in the computation of the scattered values. */ void setVariance(const models::BioModelParamType aVariance); /*! If true allows up tp 5% of the scatter points to deviate without failing. If false, all points have to correspond. */ void setAllowExceptions(const bool allow); /*! Returns a string that specifies the test the tester currently performs. */ lit::StringType getTestDescription(void) const; lit::StringType getTestName(void) const { return "ScatterTester"; }; protected: /*! performes the test and checks the results. @result Indicates if the test was successfull (true) or if it failed (false) */ bool doCheck(void) const; /*! Function will be called be check() if test was succesfull. */ void handleSuccess(void) const; /*! Function will be called be check() if test was a failure. */ void handleFailure(void) const; }; } } #endif \ No newline at end of file