diff --git a/README.md b/README.md index 156074b..0956754 100644 --- a/README.md +++ b/README.md @@ -1,272 +1,271 @@ # RTToolbox RTToolbox is a software library to support quantitative analysis of treatment outcome for radiotherapy. The RTToolbox was designed following object-oriented design (OOD) principles and was implemented in the language C++. Features include: * import of radiotherapy data (e.g. dose distributions and structure sets) from DICOM-RT format and other standard image processing formats * DVH calculation * Dose statistic calculation * arithmetic operations on dose distributions * structure relationship analyses (e.g. fully-contained, partially-contained) * Calculation of dose comparison indices such as Conformity Index (CI), Homogeneity Index (HI) and Conformation Number (CN) * Calculation of biological models including TCP, NTCP, EUD and BED Also, the RTToolbox provides apps e.g. for DVH/Dose Statistic calculation or Dose accumulation that provides a convenient access of RT scenarios without computer-science knowledge. ## Getting Started These instructions will get you a copy of the project up and running on your local machine for development and testing purposes. ### Prerequisites #### Build system * [CMake](https://cmake.org), version 3.1 or higher #### Compiler * Visual Studio 2013 * Visual Studio 2015 * Visual Studio 2017 * GCC 5.4 * GCC 7.3 Other compilers may work as well, but are not tested. #### Linking Static/Dynamic library support Can be changed with advanced option `BUILD_SHARED_LIBS` :warning: building RTToolbox as dynamic library under Windows and as static library under Linux is an experimental feature. #### Third party libraries * [boost](http://www.boost.org ), version 1.64.0 or higher * [DCMTK](http://dicom.offis.de/dcmtk.php.en ), with RT support - 3.6.1_20121102 or newer * [ITK](https://itk.org ), version 4.4 or higher (*optional*) * for DoseInterpolation support with ITK transformation or ITK File IO support * [MatchPoint](http://mitk.org/download/thirdparty/MatchPoint_rev1610.tar.gz ), version 0.12 or higher (*optional*) * for DoseInterpolation support with MatchPoint registration objects :information_source: To make sure everything runs smoothly, please make sure that all libraries and the RTToolbox are either compiled with `/MD` or `/MT` flags. ##### Boost In case you work with Windows, we recommend using the [pre-build versions of boost](https://sourceforge.net/projects/boost/files/boost-binaries/). If you want to build the library yourself, consider the following: Build (using the same compiler options as RTToolbox, usually `STATIC LINKING` and `x64` architecture). The following components are needed: * `filesystem`, -* `system`, -* `thread` and +* `system` and * `program_options` * if you plan to build the apps (*optional*) :information_source: eventually, it might be needed to add the CMake variable `BOOST_LIBRARYDIR` and set it to the respective library path of boost. For Windows: To build Boost open a command prompt, change to your boost source directory and copy following command(s): Debug: `b2 -j12 --with-filesystem --with-system --with-thread --with-program_options --with-date_time --with-atomic --with-chrono toolset=msvc-14.1 address-model=64 variant=debug threading=multi link=shared define=_BIND_TO_CURRENT_VCLIBS_VERSION` Release: `b2 -j12 --with-filesystem --with-system --with-thread --with-program_options --with-date_time --with-atomic --with-chrono toolset=msvc-14.1 address-model=64 variant=release threading=multi link=shared` If you don´t require `program_options` delete `--with-program_options` from the command before executing it. ##### DCMTK For Windows: To compile DCMTK with `/MD` flags (standard for all other libs), you need to patch the CMAKE options of DCMTK (`PathToDCMTK\CMake\dcmtkPrepare.cmake`), either by replacing `"/MT"` with `"/MD"` or by explicitly replacing lines 135 to 171 with the following lines: ``` IF(DCMTK_OVERWRITE_WIN32_COMPILER_FLAGS AND NOT BUILD_SHARED_LIBS) # settings for Microsoft Visual Studio IF(CMAKE_GENERATOR MATCHES "Visual Studio .*") # get Visual Studio Version STRING(REGEX REPLACE "Visual Studio ([0-9]+).*" "\\1" VS_VERSION "${CMAKE_GENERATOR}") # these settings never change even for C or C++ SET(CMAKE_C_FLAGS_DEBUG "/MDd /Z7 /Od") SET(CMAKE_C_FLAGS_RELEASE "/DNDEBUG /MD /O2") SET(CMAKE_C_FLAGS_MINSIZEREL "/DNDEBUG /MD /O2") SET(CMAKE_C_FLAGS_RELWITHDEBINFO "/DNDEBUG /MDd /Z7 /Od") SET(CMAKE_CXX_FLAGS_DEBUG "/MDd /Z7 /Od") SET(CMAKE_CXX_FLAGS_RELEASE "/DNDEBUG /MD /O2") SET(CMAKE_CXX_FLAGS_MINSIZEREL "/DNDEBUG /MD /O2") SET(CMAKE_CXX_FLAGS_RELWITHDEBINFO "/DNDEBUG /MDd /Z7 /Od") # specific settings for the various Visual Studio versions IF(VS_VERSION EQUAL 6) SET(CMAKE_C_FLAGS "/nologo /W3 /GX /Gy /YX") SET(CMAKE_CXX_FLAGS "/nologo /W3 /GX /Gy /YX /Zm500") # /Zm500 increments heap size which is needed on some system to compile templates in dcmimgle ENDIF(VS_VERSION EQUAL 6) IF(VS_VERSION EQUAL 7) SET(CMAKE_C_FLAGS "/nologo /W3 /Gy") SET(CMAKE_CXX_FLAGS "/nologo /W3 /Gy") ENDIF(VS_VERSION EQUAL 7) IF(VS_VERSION GREATER 7) SET(CMAKE_C_FLAGS "/nologo /W3 /Gy /EHsc") SET(CMAKE_CXX_FLAGS "/nologo /W3 /Gy /EHsc") ENDIF(VS_VERSION GREATER 7) ENDIF(CMAKE_GENERATOR MATCHES "Visual Studio .*") ENDIF(DCMTK_OVERWRITE_WIN32_COMPILER_FLAGS AND NOT BUILD_SHARED_LIBS) ``` `BUILD_APPS` can be switched off. Then build DCMTK. For Linux: install required dependencies (Ubuntu 18.04 and newer): `sudo apt-get install libpng-dev libtiff5-dev libxml2-dev libjpeg8-dev zlib1g-dev libwrap0-dev libssl-dev` install required dependencies (Ubuntu 17.10 and older): `sudo apt-get install libpng12-dev libtiff5-dev libxml2-dev libjpeg8-dev zlib1g-dev libwrap0-dev libssl-dev` Enable `BUILD_SHARED_LIBS`. `BUILD_APPS` can be switched off. ##### ITK Build ITK with default options. :warning: ensure that compiler enables C++11 features by setting `CMAKE_CXX_STANDARD=11` (default for supported compilers) :warning: if you use ITK 5, turn on `RTTB_ITK5_SUPPORT` :warning: Only use one ITK version consistently throughout all libraries and RTToolbox! Otherwise, linker errors will occur. ##### MatchPoint Configure MatchPoint. Please disable `BUILD_TESTING` before building it. :warning: ensure that compiler enables C++11 features by setting `CMAKE_CXX_STANDARD=11` (default for supported compilers) :warning: Only use one ITK version consistently throughout all libraries and RTToolbox! Otherwise, linker errors will occur. :warning: MatchPoint currently does not support ITK 5 ### Building RT-Toolbox * Configure with CMake * Set `BOOST_INCLUDE_DIR` to the main boost directory. Eventually set `BOOST_LIBRARYDIR` to the respective path (e.g. `/lib64-msvc-14.1\` for Visual Studio 2017 and 64-bit) * Select all packages you like to build (Parameters `BUILD_*` ; e.g. `BUILD_IO_Dicom`). * `BUILD_IO_Dicom`: Reading and writing of DICOM-RT files * `BUILD_IO_HELAX`: Reading of Helax DICOM files * `BUILD_IO_ITK`: Generic reading/writing with ITK * `BUILD_Interpolation`: Dose Interpolation * `BUILD_InterpolationMatchPointTransformation`: Dose Interpolation with Match Point registration support. * `BUILD_Masks`: Voxelization support * `BUILD_Models`: Calculation of dosimetrical models like TCP, NTCP etc. * `BUILD_Apps`: To build the RTTB command line apps (five available) * `BioModelCalc`: calculate the radiobiological effect based on dose * `DoseAcc`: Do dose accumulation * `DoseMap`: Do dose mapping * `DoseTool`: Compute Dose statistics and DVH * `VoxelizerTool`: Voxelize an RTSTRUCT file Some modules of RT-Toolbox are mandatory (e.g. `RTTBCore`) and build automatically. :information_source: enabling `BUILD_All_Modules` builds all modules (except Apps and Testing modules). :information_source: if you build RTTB with VS dynamic, you must ensure that code that uses RTTB DLLs uses the same STL Set `DCMTK_DIR` to your dcmtk binary file directory and `DCMTK_SOURCE_DIR` to your dcmtk source directory. If you want to build RT-Toolbox with ITK and/or MatchPoint set your `ITK_DIR` to your itk binary file directory and/or `MatchPoint_DIR` to your binary matchpoint directory. All directory entries left empty do not require a manual input. Finally, Generate the compilation files for your environment and built it. ### Examples Some examples can be found in ´testing/examples´: * `RTBioModelExampleTest`: Computation of Biological model indices (TCP/NTCP) from a given DVH * `RTDoseIndexTest`: Computation of dose indices (Conformation Number, Conformal Index, Conformity index) from a given DVH * `RTDoseStatisticsDicomTest`: Computation of dose statistics (max dose/mean dose/min dose/Dx/Vx) based on dose data for a specified structure * `RTDVHTest`: Computation of statistics (max value/mean value/min value/Dx/Vx) based on a DVH Other examples include: * `DVHCalculatorTest` (`testing/core`): Computation of a DVH from dose and structure * `VoxelizationValidationTest` (`testing/validation`): Computation of a voxelization * `ITKDoseAccessorConverterTest`: (`testing/io/itk`): Saving image RTToolbox image data as an ITK file ## Running the tests [CTest](https://cmake.org/Wiki/CMake/Testing_With_CTest) is used as testing framework. See their documentation for general testing questions. :information_source: The used testing library Litmus is build automatically. :warning: currently, you have access to testing data only with ssh. That means that a [phabricator](https://phabricator.mitk.org/) account and access to `RTTB-data` repository is mandatory. Please contact rttb(at)dkfz.de for further information. Enabling testing is done as follows: * Enable `BUILD_TESTING` * Configure with CMake * Enable tests of interest * Generate CMake configuration * Build RT-Toolbox * Run tests (build `RUN_TESTS` project or call `ctest` in commandline) to ensure that everything is correct. :information_source: `BUILD_Tester_All` builds all test modules. ## Contributing Please add a github issue and send a pull request if you want to contribute. ## Versioning We use the Ubuntu Release versioning scheme. v2017.02 was released in February 2017. We aim at releasing stable versions once a year. For the versions available, see the [tags on this repository](https://github.com/MIC-DKFZ/RTTB/tags). ## Authors See the list of [contributors](https://github.com/MIC-DKFZ/RTTB/contributors) who participated in this project. ## License This project is licensed under the BSD License - see the [LICENSE](LICENSE) file for details ## Contact Software Development for Integrated Diagnostics and Therapy (SIDT), German Cancer Research Center (DKFZ), Heidelberg, Germany. Web: https://www.dkfz-heidelberg.de/en/mic/research/SIDT/sidt_projects.html E-mail: rttb(at)dkfz.de ## Acknowledgments * **Billie Thompson** - *Template of the readme* - [PurpleBooth](https://github.com/PurpleBooth) diff --git a/cmake/PackageDepends/RTTB_BoostBinaries_Config.cmake b/cmake/PackageDepends/RTTB_BoostBinaries_Config.cmake index d0f34ec..dc73c72 100644 --- a/cmake/PackageDepends/RTTB_BoostBinaries_Config.cmake +++ b/cmake/PackageDepends/RTTB_BoostBinaries_Config.cmake @@ -1,19 +1,19 @@ IF(NOT DEFINED RTTB_USE_MITK_BOOST) OPTION(RTTB_USE_MITK_BOOST "RTTB should use a boost which is available in the MITK superbuild external projects structure." OFF) MARK_AS_ADVANCED(RTTB_USE_MITK_BOOST) ENDIF(NOT DEFINED RTTB_USE_MITK_BOOST) IF(BUILD_SHARED_LIBS OR RTTB_USE_MITK_BOOST) SET(Boost_USE_STATIC_LIBS OFF) ADD_DEFINITIONS(-DBOOST_ALL_DYN_LINK) ELSE(BUILD_SHARED_LIBS OR RTTB_USE_MITK_BOOST) SET(Boost_USE_STATIC_LIBS ON) ENDIF(BUILD_SHARED_LIBS OR RTTB_USE_MITK_BOOST) SET(BOOST_MIN_VERSION "1.64.0") - FIND_PACKAGE(Boost ${BOOST_MIN_VERSION} REQUIRED COMPONENTS filesystem system thread ${RTTB_Boost_ADDITIONAL_COMPONENT} QUIET) + FIND_PACKAGE(Boost ${BOOST_MIN_VERSION} REQUIRED COMPONENTS filesystem system ${RTTB_Boost_ADDITIONAL_COMPONENT} QUIET) if(Boost_LIBRARIES) LIST(APPEND ALL_INCLUDE_DIRECTORIES ${Boost_INCLUDE_DIRS}) LIST(APPEND ALL_LIBRARIES ${Boost_LIBRARIES}) link_directories(${Boost_LIBRARY_DIRS}) endif() diff --git a/code/algorithms/rttbDoseStatisticsCalculator.cpp b/code/algorithms/rttbDoseStatisticsCalculator.cpp index 6969acb..5b42344 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.cpp +++ b/code/algorithms/rttbDoseStatisticsCalculator.cpp @@ -1,413 +1,410 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseStatisticsCalculator.h" #include #include #include -#include #include "rttbNullPointerException.h" #include "rttbInvalidDoseException.h" #include "rttbInvalidParameterException.h" -#include - namespace rttb { namespace algorithms { DoseStatisticsCalculator::DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator) { if (aDoseIterator == nullptr) { throw core::NullPointerException("DoseIterator must not be nullptr"); } else { _doseIterator = aDoseIterator; } _simpleDoseStatisticsCalculated = false; _complexDoseStatisticsCalculated = false; _multiThreading = false; - _mutex = ::boost::make_shared<::boost::shared_mutex>(); + _mutex = ::boost::make_shared(); } DoseStatisticsCalculator::~DoseStatisticsCalculator() = default; DoseStatisticsCalculator::DoseIteratorPointer DoseStatisticsCalculator::getDoseIterator() const { return _doseIterator; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( bool computeComplexMeasures, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (_doseIterator == nullptr) { throw core::NullPointerException("_doseIterator must not be nullptr!"); } //"simple" dose statistics are mandatory calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); if (computeComplexMeasures) { //more complex dose statistics are optional with default maximum dose and default relative x values calculateComplexDoseStatistics(_statistics->getMaximum(), std::vector(), std::vector()); } return _statistics; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (_doseIterator == nullptr) { throw core::NullPointerException("_doseIterator must not be nullptr!"); } if (referenceDose <= 0) { throw rttb::core::InvalidParameterException("Reference dose must be > 0 !"); } //simple dose statistics calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); //more complex dose statistics with given reference dose and default x values calculateComplexDoseStatistics(referenceDose, std::vector(), std::vector()); return _statistics; } DoseStatisticsCalculator::DoseStatisticsPointer DoseStatisticsCalculator::calculateDoseStatistics( const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues, DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { if (_doseIterator == nullptr) { throw core::NullPointerException("_doseIterator must not be nullptr!"); } //"simple" dose statistics calculateSimpleDoseStatistics(maxNumberMinimaPositions, maxNumberMaximaPositions); if (referenceDose <= 0) { //more complex dose statistics with default maximum dose and relative x values calculateComplexDoseStatistics(_statistics->getMaximum(), precomputeDoseValues, precomputeVolumeValues); } else { //more complex dose statistics with given reference dose and relative x values calculateComplexDoseStatistics(referenceDose, precomputeDoseValues, precomputeVolumeValues); } return _statistics; } void DoseStatisticsCalculator::calculateSimpleDoseStatistics(unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions) { _doseVector.clear(); _voxelProportionVector.clear(); std::multimap doseValueVSIndexMap; std::vector voxelProportionVectorTemp; DoseStatisticType maximumDose = 0; DoseStatisticType minimumDose = std::numeric_limits::max(); DoseStatisticType meanDose; DoseStatisticType stdDeviationDose; DoseTypeGy sum = 0; VolumeType numVoxels = 0.0; DoseTypeGy squareSum = 0; VolumeType volume = 0; _doseIterator->reset(); int i = 0; DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid()) { doseValue = _doseIterator->getCurrentDoseValue(); if (i == 0) { minimumDose = doseValue; volume = _doseIterator->getCurrentVoxelVolume(); } rttb::FractionType voxelProportion = _doseIterator->getCurrentRelevantVolumeFraction(); sum += doseValue * voxelProportion; numVoxels += voxelProportion; squareSum += doseValue * doseValue * voxelProportion; if (doseValue > maximumDose) { maximumDose = doseValue; } else if (doseValue < minimumDose) { minimumDose = doseValue; } voxelProportionVectorTemp.push_back(voxelProportion); doseValueVSIndexMap.insert(std::pair(doseValue, i)); i++; _doseIterator->next(); } if (numVoxels != 0) { meanDose = sum / numVoxels; //standard deviation is defined only for n>=2 if (numVoxels >= 2) { //uncorrected variance is calculated DoseStatisticType varianceDose = (squareSum / numVoxels - meanDose * meanDose); if (varianceDose < errorConstant) { stdDeviationDose = 0; } else { stdDeviationDose = pow(varianceDose, 0.5); } } else { stdDeviationDose = 0; } } //sort dose values and corresponding volume fractions in member variables for (auto & it : doseValueVSIndexMap) { _doseVector.push_back((float)it.first); _voxelProportionVector.push_back(voxelProportionVectorTemp.at(it.second)); } volume *= numVoxels; _statistics = boost::make_shared(minimumDose, maximumDose, meanDose, stdDeviationDose, numVoxels, volume); _simpleDoseStatisticsCalculated = true; ResultListPointer minimumVoxelPositions = computeMinimumPositions(maxNumberMinimaPositions); ResultListPointer maximumVoxelPositions = computeMaximumPositions(maxNumberMaximaPositions); _statistics->setMinimumVoxelPositions(minimumVoxelPositions); _statistics->setMaximumVoxelPositions(maximumVoxelPositions); } void DoseStatisticsCalculator::calculateComplexDoseStatistics(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues) { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call calculateComplexDoseStatistics()"); } std::vector precomputeDoseValuesNonConst = precomputeDoseValues; std::vector precomputeVolumeValuesNonConst = precomputeVolumeValues; //set default values if (precomputeDoseValues.empty()) { std::vector defaultPrecomputeDoseValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)( 0.95)(0.98); precomputeDoseValuesNonConst = defaultPrecomputeDoseValues; } if (precomputeVolumeValues.empty()) { std::vector defaultPrecomputeVolumeValues = boost::assign::list_of(0.02)(0.05)(0.1)(0.9)( 0.95)(0.98); precomputeVolumeValuesNonConst = defaultPrecomputeVolumeValues; } _Vx = ::boost::make_shared(precomputeDoseValuesNonConst, referenceDose, _doseIterator); _Vx->compute(); _Dx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), _statistics->getMinimum()); _Dx->compute(); _MOHx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume()); _MOHx->compute(); _MOCx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume()); _MOCx->compute(); _MaxOHx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume()); _MaxOHx->compute(); _MinOCx = ::boost::make_shared(precomputeVolumeValuesNonConst, _statistics->getVolume(), this->_doseVector, this->_voxelProportionVector, this->_doseIterator->getCurrentVoxelVolume(), _statistics->getMinimum(), _statistics->getMaximum()); _MinOCx->compute(); _statistics->setVx(_Vx->getMeasureCollection()); _statistics->setDx(_Dx->getMeasureCollection()); _statistics->setMOHx(_MOHx->getMeasureCollection()); _statistics->setMOCx(_MOCx->getMeasureCollection()); _statistics->setMaxOHx(_MaxOHx->getMeasureCollection()); _statistics->setMinOCx(_MinOCx->getMeasureCollection()); _statistics->setReferenceDose(referenceDose); _complexDoseStatisticsCalculated = true; } void DoseStatisticsCalculator::addPrecomputeValues(const std::vector& values) { if (!_complexDoseStatisticsCalculated) { throw core::InvalidDoseException("Complex DoseStatistics have to be computed in order to call addPrecomputeDoseValues()"); } _Vx->addPrecomputeDoseValues(values); _Dx->addPrecomputeVolumeValues(values); _MOHx->addPrecomputeVolumeValues(values); _MOCx->addPrecomputeVolumeValues(values); _MaxOHx->addPrecomputeVolumeValues(values); _MinOCx->addPrecomputeVolumeValues(values); } void DoseStatisticsCalculator::recalculateDoseStatistics() { if (!_complexDoseStatisticsCalculated) { throw core::InvalidDoseException("Complex DoseStatistics have to be computed in order to call recalculateDoseStatistics()"); } _Vx->compute(); _Dx->compute(); _MOHx->compute(); _MOCx->compute(); _MaxOHx->compute(); _MinOCx->compute(); } DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMaximumPositions( unsigned int maxNumberMaxima) const { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMaximumPositions()"); } ResultListPointer maxVoxelVector = boost::make_shared > >(); unsigned int count = 0; this->_doseIterator->reset(); DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid() && count < maxNumberMaxima) { doseValue = _doseIterator->getCurrentDoseValue(); if (doseValue == _statistics->getMaximum()) { VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID(); std::pair voxel(doseValue, currentID); maxVoxelVector->push_back(voxel); count++; } _doseIterator->next(); } return maxVoxelVector; } DoseStatisticsCalculator::ResultListPointer DoseStatisticsCalculator::computeMinimumPositions( unsigned int maxNumberMinima) const { if (!_simpleDoseStatisticsCalculated) { throw core::InvalidDoseException("simple DoseStatistics have to be computed in order to call computeMinimumPositions()"); } ResultListPointer minVoxelVector = boost::make_shared > >(); /*! @todo: Architecture Annotation: Finding the positions for the minimum only once reduces computation time, but will require sensible use by the programmers. To be save the output vector minVoxelVector will be always cleared here to garantee that no false values are presented. This change may be revoced to increase computation speed later on (only compute if(minVoxelVector->size()==0)). */ unsigned int count = 0; this->_doseIterator->reset(); DoseTypeGy doseValue = 0; while (_doseIterator->isPositionValid() && count < maxNumberMinima) { doseValue = _doseIterator->getCurrentDoseValue(); if (doseValue == _statistics->getMinimum()) { VoxelGridID currentID = _doseIterator->getCurrentVoxelGridID(); std::pair voxel(doseValue, currentID); minVoxelVector->push_back(voxel); count++; } _doseIterator->next(); } return minVoxelVector; } void DoseStatisticsCalculator::setMultiThreading(const bool choice) { _multiThreading = choice; } }//end namespace algorithms }//end namespace rttb diff --git a/code/algorithms/rttbDoseStatisticsCalculator.h b/code/algorithms/rttbDoseStatisticsCalculator.h index d88fb09..69812b0 100644 --- a/code/algorithms/rttbDoseStatisticsCalculator.h +++ b/code/algorithms/rttbDoseStatisticsCalculator.h @@ -1,201 +1,201 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #ifndef __DOSE_STATISTICS_CALCULATOR_H #define __DOSE_STATISTICS_CALCULATOR_H #include +#include #include -#include #include "rttbDoseIteratorInterface.h" #include "rttbDoseStatistics.h" #include "RTTBAlgorithmsExports.h" #include "rttbDxVolumeToDoseMeasureCollectionCalculator.h" #include "rttbVxDoseToVolumeMeasureCollectionCalculator.h" #include "rttbMOHxVolumeToDoseMeasureCollectionCalculator.h" #include "rttbMOCxVolumeToDoseMeasureCollectionCalculator.h" #include "rttbMaxOHxVolumeToDoseMeasureCollectionCalculator.h" #include "rttbMinOCxVolumeToDoseMeasureCollectionCalculator.h" #ifdef _MSC_VER #pragma warning(push) #pragma warning(disable: 4251) #endif namespace rttb { namespace algorithms { /*! @class DoseStatisticsCalculator @brief Class for calculating different statistical values from a RT dose distribution @details These values range from standard statistical values such as minimum, maximum and mean to more complex dose specific measures such as Vx (volume irradiated with a dose >=x), Dx (minimal dose delivered to x% of the VOI) or MOHx (mean in the hottest volume). For a complete list, see calculateDoseStatistics(). @note the complex dose statistics are precomputed and cannot be computed "on the fly" lateron! The doses/volumes that should be used for precomputation have to be set in calculateDoseStatistics() */ class RTTBAlgorithms_EXPORT DoseStatisticsCalculator { public: using DoseIteratorPointer = core::DoseIteratorInterface::Pointer; using ResultListPointer = DoseStatistics::ResultListPointer; using DoseStatisticsPointer = DoseStatistics::Pointer; private: DoseIteratorPointer _doseIterator; /*! @brief Contains relevant dose values sorted in descending order. */ std::vector _doseVector; /*! @brief Contains the corresponding voxel proportions to the values in doseVector. */ std::vector _voxelProportionVector; /*! @brief The doseStatistics are stored here. */ DoseStatisticsPointer _statistics; bool _simpleDoseStatisticsCalculated; bool _complexDoseStatisticsCalculated; bool _multiThreading; - ::boost::shared_ptr _mutex; + ::boost::shared_ptr _mutex; VxDoseToVolumeMeasureCollectionCalculator::Pointer _Vx; DxVolumeToDoseMeasureCollectionCalculator::Pointer _Dx; MOHxVolumeToDoseMeasureCollectionCalculator::Pointer _MOHx; MOCxVolumeToDoseMeasureCollectionCalculator::Pointer _MOCx; MaxOHxVolumeToDoseMeasureCollectionCalculator::Pointer _MaxOHx; MinOCxVolumeToDoseMeasureCollectionCalculator::Pointer _MinOCx; /*! @brief Calculates the positions where the dose has its maximum @param maxNumberMaximaPositions the maximal amount of computed positions @pre maximumDose must be defined in _statistics with the correct value */ ResultListPointer computeMaximumPositions(unsigned int maxNumberMaximaPositions) const; /*! @brief Calculates the positions where the dose has its minimum @param maxNumberMinimaPositions the maximal amount of computed positions (they are read sequentially using the iterator until maxNumberMinimaPositions have been read, other positions are not considered) @pre minimumDose must be defined in _statistics with the correct value */ ResultListPointer computeMinimumPositions(unsigned int maxNumberMinimaPositions) const; /*! @brief Calculates simple dose statistics (min, mean, max, stdDev, minDosePositions, maxDosePositions) @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed */ void calculateSimpleDoseStatistics(unsigned int maxNumberMinimaPositions, unsigned int maxNumberMaximaPositions); /*! @brief Calculates complex dose statistics (Dx, Vx, MOHx, MOCx, MaxOHx, MinOCx) @warning computations can take quite long (>1 min) for large structures as many statistics are precomputed */ void calculateComplexDoseStatistics(DoseTypeGy referenceDose, const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues); public: ~DoseStatisticsCalculator(); /*! @brief Constructor @param aDoseIterator the dose to be analyzed */ explicit DoseStatisticsCalculator(DoseIteratorPointer aDoseIterator); DoseIteratorPointer getDoseIterator() const; /*! @brief Compute simple or complex dose statistics with default relative x values and the maximum dose as default reference dose (for Vx computation) @details The following statistics are calculated always (i.e. also if computeComplexMeasures=false):
  • minimum dose
  • mean dose
  • maximum dose
  • standard deviation dose
  • voxel positions of minimum dose
  • voxel positions of maximum dose
Additionally, these statistics are computed if computeComplexMeasures=true:
  • Dx (the minimal dose delivered to a volume >= x)
  • Vx (the volume irradiated with a dose >= x)
  • MOHx (mean dose of the hottest x volume)
  • MOCx (mean dose of the coldest x volume)
  • MaxOHx (Maximum outside of the hottest x volume)
  • MinOCx (Minimum outside of the coldest x volume)
Default x values for Vx are 0.02, 0.05, 0.1, 0.9, 0.95 and 0.98, with respect to maxDose. Default x values for Dx, MOHx, MOCx, MaxOHx and MinOCx are 0.02, 0.05, 0.1, 0.9, 0.95 and 0.98, with respect to volume. @param computeComplexMeasures should complex statistics be calculated? If it is true, the complex dose statistics will be calculated with default relative x values and the maximum dose as reference dose @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed @warning If computeComplexMeasures==true, computations can take quite long (>1 min) for large structures as many statistics are precomputed @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! Only the default x values can be requested in DoseStatistics! */ DoseStatisticsPointer calculateDoseStatistics(bool computeComplexMeasures = false, unsigned int maxNumberMinimaPositions = 10, unsigned int maxNumberMaximaPositions = 10); /*! @brief Compute complex dose statistics with given reference dose and default relative x values @param referenceDose the reference dose to compute Vx, normally it should be the prescribed dose @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed @exception InvalidParameterException thrown if referenceDose <= 0 @warning Computations can take quite long (>1 min) for large structures as many statistics are precomputed @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! Only the default x values can be requested in DoseStatistics! */ DoseStatisticsPointer calculateDoseStatistics(DoseTypeGy referenceDose, unsigned int maxNumberMinimaPositions = 10, unsigned int maxNumberMaximaPositions = 10); /*! @brief Compute complex dose statistics with given relative x values and reference dose @param precomputeDoseValues the relative dose values for Vx precomputation, e.g. 0.02, 0.05, 0.95... @param precomputeVolumeValues the relative volume values for Dx, MOHx, MOCx, MaxOHx and MinOCx precomputation, e.g. 0.02, 0.05, 0.95... @param referenceDose the reference dose to compute Vx, normally it should be the prescribed dose. Default value is the maximum dose. @param maxNumberMinimaPositions the maximal amount of computed positions where the dose has its minimum that is computed @param maxNumberMaximaPositions the maximal amount of computed positions where the dose has its maximum that is computed @warning Computations can take quite long (>1 min) for large structures as many statistics are precomputed @note The complex dose statistics are precomputed and cannot be computed "on the fly" lateron! The doses/volumes that should be used for precomputation have to be set by in precomputeDoseValues and precomputeVolumeValues. Only these values can be requested in DoseStatistics! */ DoseStatisticsPointer calculateDoseStatistics(const std::vector& precomputeDoseValues, const std::vector& precomputeVolumeValues, DoseTypeGy referenceDose = -1, unsigned int maxNumberMinimaPositions = 10, unsigned int maxNumberMaximaPositions = 10); /*! @brief Adds additonal precompute values for all complex Dose Statistics @exception InvalidDoseException if complexDoseStatistics are not already calculated */ void addPrecomputeValues(const std::vector& values); /*! @brief Recalculates the DoseMeasures for all complex Dose Statistics @exception InvalidDoseException if complexDoseStatistics are not already calculated */ void recalculateDoseStatistics(); void setMultiThreading(bool choice); }; } } #ifdef _MSC_VER #pragma warning(pop) #endif #endif diff --git a/code/algorithms/rttbDoseToVolumeMeasureCollectionCalculator.cpp b/code/algorithms/rttbDoseToVolumeMeasureCollectionCalculator.cpp index 419c861..5d9bde6 100644 --- a/code/algorithms/rttbDoseToVolumeMeasureCollectionCalculator.cpp +++ b/code/algorithms/rttbDoseToVolumeMeasureCollectionCalculator.cpp @@ -1,86 +1,85 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbDoseToVolumeMeasureCollectionCalculator.h" -#include + +#include + #include "rttbInvalidParameterException.h" #include "rttbUtils.h" #include -//#include - namespace rttb { namespace algorithms { DoseToVolumeMeasureCollectionCalculator::DoseToVolumeMeasureCollectionCalculator(const std::vector& precomputeDoseValues, const DoseTypeGy referenceDose, const core::DoseIteratorInterface::Pointer doseIterator, DoseToVolumeMeasureCollection::complexStatistics name, bool multiThreading) : _doseIterator(doseIterator), _referenceDose(referenceDose), _measureCollection(::boost::make_shared(name)), _multiThreading(multiThreading) { addPrecomputeDoseValues(precomputeDoseValues); } void DoseToVolumeMeasureCollectionCalculator::compute() { - std::vector threads; + std::vector threads; for (double _precomputeDoseValue : _precomputeDoseValues) { double xAbsolute = _precomputeDoseValue * _referenceDose; if (!rttb::core::isKey(_measureCollection->getAllValues(), xAbsolute)) { if (_multiThreading) { throw rttb::core::InvalidParameterException("MultiThreading is not implemented yet."); - //threads.push_back(boost::thread(&DoseToVolumeMeasureCollectionCalculator::insertIntoMeasureCollection, this , xAbsolute, computeSpecificValue(xAbsolute))); } else { insertIntoMeasureCollection(xAbsolute, this->computeSpecificValue(xAbsolute)); } } } for (auto & thread : threads) { thread.join(); } } void DoseToVolumeMeasureCollectionCalculator::addPrecomputeDoseValues(const std::vector& values) { for (double value : values) { if (value > 1 || value < 0) { throw rttb::core::InvalidParameterException("Values must be between 1 and 0!"); } if (!rttb::core::isKey(_precomputeDoseValues, value)) { _precomputeDoseValues.push_back(value); } } } DoseToVolumeMeasureCollection::Pointer DoseToVolumeMeasureCollectionCalculator::getMeasureCollection() { return _measureCollection; } void DoseToVolumeMeasureCollectionCalculator::insertIntoMeasureCollection(DoseTypeGy xAbsolute, VolumeType resultVolume) { _measureCollection->insertValue(xAbsolute, resultVolume); } } } diff --git a/code/algorithms/rttbVolumeToDoseMeasureCollectionCalculator.cpp b/code/algorithms/rttbVolumeToDoseMeasureCollectionCalculator.cpp index abc876b..dea50cf 100644 --- a/code/algorithms/rttbVolumeToDoseMeasureCollectionCalculator.cpp +++ b/code/algorithms/rttbVolumeToDoseMeasureCollectionCalculator.cpp @@ -1,84 +1,87 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbVolumeToDoseMeasureCollectionCalculator.h" -#include + +#include + #include "rttbInvalidParameterException.h" #include "rttbUtils.h" #include -//#include namespace rttb { namespace algorithms { VolumeToDoseMeasureCollectionCalculator::VolumeToDoseMeasureCollectionCalculator(const std::vector& precomputeVolumeValues, const VolumeType volume, const std::vector& doseVector, const std::vector& voxelProportionVector, const DoseVoxelVolumeType currentVoxelVolume, VolumeToDoseMeasureCollection::complexStatistics name, bool multiThreading) : _doseVector(doseVector), _currentVoxelVolume(currentVoxelVolume), _voxelProportionVector(voxelProportionVector), _volume(volume), _measureCollection(::boost::make_shared(name)), _multiThreading(multiThreading) { addPrecomputeVolumeValues(precomputeVolumeValues); } void VolumeToDoseMeasureCollectionCalculator::compute() { - std::vector threads; + std::vector threads; for (double _precomputeVolumeValue : _precomputeVolumeValues) { double xAbsolute = _precomputeVolumeValue * _volume; if (!rttb::core::isKey(_measureCollection->getAllValues(), xAbsolute)) { if (_multiThreading) { throw rttb::core::InvalidParameterException("MultiThreading is not implemented yet."); - //threads.push_back(boost::thread(&VolumeToDoseMeasureCollectionCalculator::insertIntoMeasureCollection, this, xAbsolute, computeSpecificValue(xAbsolute))); } else { insertIntoMeasureCollection(xAbsolute, this->computeSpecificValue(xAbsolute)); } } } for (auto & thread : threads) { - thread.join(); + if (thread.joinable()) + { + thread.join(); + } } } void VolumeToDoseMeasureCollectionCalculator::addPrecomputeVolumeValues(const std::vector& values) { for (double value : values) { if (value > 1 || value < 0) { throw rttb::core::InvalidParameterException("Values must be between 1 and 0!"); } if (!rttb::core::isKey(_precomputeVolumeValues, value)) { _precomputeVolumeValues.push_back(value); } } } VolumeToDoseMeasureCollection::Pointer VolumeToDoseMeasureCollectionCalculator::getMeasureCollection() { return _measureCollection; } void VolumeToDoseMeasureCollectionCalculator::insertIntoMeasureCollection(VolumeType xAbsolute, DoseTypeGy resultDose) { _measureCollection->insertValue(xAbsolute, resultDose); } } } diff --git a/code/masks/rttbBoostMask.cpp b/code/masks/rttbBoostMask.cpp index 696ddc5..f0115eb 100644 --- a/code/masks/rttbBoostMask.cpp +++ b/code/masks/rttbBoostMask.cpp @@ -1,558 +1,571 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include +#include #include #include #include -#include #include "rttbBoostMask.h" #include "rttbNullPointerException.h" #include "rttbInvalidParameterException.h" #include "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbBoostMaskVoxelizationThread.h" namespace rttb { namespace masks { namespace boost { BoostMask::BoostMask(core::GeometricInfo::Pointer aDoseGeoInfo, core::Structure::Pointer aStructure, bool strict, unsigned int numberOfThreads) : _geometricInfo(aDoseGeoInfo), _structure(aStructure), _strict(strict), _numberOfThreads(numberOfThreads), _voxelizationThickness(0.0), _voxelInStructure(::boost::make_shared()) { _isUpToDate = false; if (_geometricInfo == nullptr) { throw rttb::core::NullPointerException("Error: Geometric info is nullptr!"); } else if (_structure == nullptr) { throw rttb::core::NullPointerException("Error: Structure is nullptr!"); } if (_numberOfThreads == 0) { - _numberOfThreads = ::boost::thread::hardware_concurrency(); + _numberOfThreads = std::thread::hardware_concurrency(); if (_numberOfThreads == 0) { throw rttb::core::InvalidParameterException("Error: detection of the number of hardware threads is not possible. Please specify number of threads for voxelization explicitly as parameter in BoostMask."); } } } BoostMask::MaskVoxelListPointer BoostMask::getRelevantVoxelVector() { if (!_isUpToDate) { calcMask(); } return _voxelInStructure; } void BoostMask::calcMask() { preprocessing(); voxelization(); generateMaskVoxelList(); _isUpToDate = true; } void BoostMask::preprocessing() { rttb::PolygonSequenceType polygonSequence = _structure->getStructureVector(); //Convert world coordinate polygons to the polygons with geometry coordinate rttb::PolygonSequenceType geometryCoordinatePolygonVector; rttb::PolygonSequenceType::iterator it; rttb::ContinuousVoxelGridIndex3D globalMaxGridIndex(std::numeric_limits::min(), std::numeric_limits::min(), std::numeric_limits::min()); rttb::ContinuousVoxelGridIndex3D globalMinGridIndex(_geometricInfo->getNumColumns(), _geometricInfo->getNumRows(), 0); for (it = polygonSequence.begin(); it != polygonSequence.end(); ++it) { PolygonType rttbPolygon = *it; PolygonType geometryCoordinatePolygon; //1. convert polygon to geometry coordinate polygons //2. calculate global min/max //3. check if polygon is planar if (!preprocessingPolygon(rttbPolygon, geometryCoordinatePolygon, globalMinGridIndex, globalMaxGridIndex, errorConstant)) { throw rttb::core::Exception("TiltedMaskPlaneException"); } geometryCoordinatePolygonVector.push_back(geometryCoordinatePolygon); } rttb::VoxelGridIndex3D minIndex = VoxelGridIndex3D(GridIndexType(globalMinGridIndex(0) ), GridIndexType(globalMinGridIndex(1) ), GridIndexType(globalMinGridIndex(2) )); rttb::VoxelGridIndex3D maxIndex = VoxelGridIndex3D(GridIndexType(globalMaxGridIndex(0) ), GridIndexType(globalMaxGridIndex(1) ), GridIndexType(globalMaxGridIndex(2) )); _globalBoundingBox.push_back(minIndex); _globalBoundingBox.push_back(maxIndex); //convert rttb polygon sequence to a map of z index and a vector of boost ring 2d (without holes) _ringMap = convertRTTBPolygonSequenceToBoostRingMap(geometryCoordinatePolygonVector); } void BoostMask::voxelization() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } BoostRingMap::iterator itMap; size_t mapSizeInAThread = _ringMap.size() / _numberOfThreads; unsigned int count = 0; unsigned int countThread = 0; BoostPolygonMap polygonMap; std::vector polygonMapVector; //check donut and convert to a map of z index and a vector of boost polygon 2d (with or without holes) for (itMap = _ringMap.begin(); itMap != _ringMap.end(); ++itMap) { //the vector of all boost 2d polygons with the same z grid index(donut polygon is accepted). BoostPolygonVector polygonVector = checkDonutAndConvert((*itMap).second); if (count == mapSizeInAThread && countThread < (_numberOfThreads - 1)) { polygonMapVector.push_back(polygonMap); polygonMap.clear(); count = 0; countThread++; } polygonMap.insert(std::pair((*itMap).first, polygonVector)); count++; } _voxelizationMap = ::boost::make_shared >(); polygonMapVector.push_back(polygonMap); //insert the last one //generate voxelization map, multi-threading - ::boost::thread_group threads; - auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); + std::vector threads; + + auto aMutex = ::boost::make_shared(); for (const auto & i : polygonMapVector) { - BoostMaskVoxelizationThread t(i, _globalBoundingBox, - _voxelizationMap, aMutex, _strict); - threads.create_thread(t); + BoostMaskVoxelizationThread t(i, _globalBoundingBox, + _voxelizationMap, aMutex, _strict); + threads.emplace_back(t); } - threads.join_all(); + for (auto& thread : threads) + { + if (thread.joinable()) + { + thread.join(); + } + } } void BoostMask::generateMaskVoxelList() { if (_globalBoundingBox.size() < 2) { throw rttb::core::InvalidParameterException("Bounding box calculation failed! "); } //check homogeneous of the voxelization plane (the contours plane) if (!calcVoxelizationThickness(_voxelizationThickness)) { throw rttb::core::InvalidParameterException("Error: The contour plane should be homogeneous!"); } - ::boost::thread_group threads; - auto aMutex = ::boost::make_shared<::boost::shared_mutex>(); + std::vector threads; + auto aMutex = ::boost::make_shared(); unsigned int sliceNumberInAThread = _geometricInfo->getNumSlices() / _numberOfThreads; //generate mask voxel list, multi-threading for (unsigned int i = 0; i < _numberOfThreads; ++i) { unsigned int beginSlice = i * sliceNumberInAThread; unsigned int endSlice; if (i < _numberOfThreads - 1) { endSlice = (i + 1) * sliceNumberInAThread; } else { endSlice = _geometricInfo->getNumSlices(); } BoostMaskGenerateMaskVoxelListThread t(_globalBoundingBox, _geometricInfo, _voxelizationMap, _voxelizationThickness, beginSlice, endSlice, _voxelInStructure, _strict, aMutex); - threads.create_thread(t); + threads.emplace_back(t); } - threads.join_all(); + for (auto& thread : threads) + { + if (thread.joinable()) + { + thread.join(); + } + } } bool BoostMask::preprocessingPolygon(const rttb::PolygonType& aRTTBPolygon, rttb::PolygonType& geometryCoordinatePolygon, rttb::ContinuousVoxelGridIndex3D& minimum, rttb::ContinuousVoxelGridIndex3D& maximum, double aErrorConstant) const { double minZ = _geometricInfo->getNumSlices(); double maxZ = 0.0; for (auto worldCoordinatePoint : aRTTBPolygon) { //convert to geometry coordinate polygon rttb::ContinuousVoxelGridIndex3D geometryCoordinatePoint; _geometricInfo->worldCoordinateToContinuousIndex(worldCoordinatePoint, geometryCoordinatePoint); geometryCoordinatePolygon.push_back(geometryCoordinatePoint); //calculate the current global min/max //min and max for x if (geometryCoordinatePoint(0) < minimum(0)) { minimum(0) = geometryCoordinatePoint(0); } if (geometryCoordinatePoint(0) > maximum(0)) { maximum(0) = geometryCoordinatePoint(0); } //min and max for y if (geometryCoordinatePoint(1) < minimum(1)) { minimum(1) = geometryCoordinatePoint(1); } if (geometryCoordinatePoint(1) > maximum(1)) { maximum(1) = geometryCoordinatePoint(1); } //min and max for z if (geometryCoordinatePoint(2) < minimum(2)) { minimum(2) = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maximum(2)) { maximum(2) = geometryCoordinatePoint(2); } //check planar if (geometryCoordinatePoint(2) < minZ) { minZ = geometryCoordinatePoint(2); } if (geometryCoordinatePoint(2) > maxZ) { maxZ = geometryCoordinatePoint(2); } } return (std::abs(maxZ - minZ) <= aErrorConstant); } BoostMask::BoostRing2D BoostMask::convertRTTBPolygonToBoostRing(const rttb::PolygonType& aRTTBPolygon) const { BoostMask::BoostRing2D polygon2D; BoostPoint2D firstPoint; for (unsigned int i = 0; i < aRTTBPolygon.size(); i++) { rttb::WorldCoordinate3D rttbPoint = aRTTBPolygon.at(i); BoostPoint2D boostPoint(rttbPoint[0], rttbPoint[1]); if (i == 0) { firstPoint = boostPoint; } ::boost::geometry::append(polygon2D, boostPoint); } ::boost::geometry::append(polygon2D, firstPoint); return polygon2D; } BoostMask::BoostRingMap BoostMask::convertRTTBPolygonSequenceToBoostRingMap( const rttb::PolygonSequenceType& aRTTBPolygonVector) const { rttb::PolygonSequenceType::const_iterator it; BoostMask::BoostRingMap aRingMap; for (it = aRTTBPolygonVector.begin(); it != aRTTBPolygonVector.end(); ++it) { rttb::PolygonType rttbPolygon = *it; double zIndex = rttbPolygon.at(0)[2];//get the first z index of the polygon bool isFirstZ = true; if (!aRingMap.empty()) { auto findIt = findNearestKey(aRingMap, zIndex, errorConstant); //if the z index is found (same slice), add the polygon to vector if (findIt != aRingMap.end()) { //BoostRingVector ringVector = ; (*findIt).second.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); isFirstZ = false; } } //if it is the first z index in the map, insert vector with the polygon if (isFirstZ) { BoostRingVector ringVector; ringVector.push_back(convertRTTBPolygonToBoostRing(rttbPolygon)); aRingMap.insert(std::pair(zIndex, ringVector)); } } return aRingMap; } BoostMask::BoostRingMap::iterator BoostMask::findNearestKey(BoostMask::BoostRingMap& aBoostRingMap, double aIndex, double aErrorConstant) const { auto find = aBoostRingMap.find(aIndex); //if find a key equivalent to aIndex, found if (find != aBoostRingMap.end()) { return find; } else { auto lowerBound = aBoostRingMap.lower_bound(aIndex); //if all keys go before aIndex, check the last key if (lowerBound == aBoostRingMap.end()) { lowerBound = --aBoostRingMap.end(); } //if the lower bound very close to aIndex, found if (std::abs((*lowerBound).first - aIndex) <= aErrorConstant) { return lowerBound; } else { //if the lower bound is the beginning, not found if (lowerBound == aBoostRingMap.begin()) { return aBoostRingMap.end(); } else { auto lowerBound1 = --lowerBound;//the key before the lower bound //if the key before the lower bound very close to a Index, found if (std::abs((*lowerBound1).first - aIndex) <= aErrorConstant) { return lowerBound1; } //else, not found else { return aBoostRingMap.end(); } } } } } BoostMask::BoostPolygonVector BoostMask::checkDonutAndConvert(const BoostMask::BoostRingVector& aRingVector) const { //check donut BoostMask::BoostRingVector::const_iterator it1; BoostMask::BoostRingVector::const_iterator it2; BoostMask::BoostPolygonVector boostPolygonVector; std::vector donutIndexVector;//store the outer and inner ring index BoostMask::BoostPolygonVector donutVector;//store new generated donut polygon //Get donut index and donut polygon unsigned int index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++) { bool it1IsDonut = false; //check if the ring is already determined as a donut for (unsigned int i : donutIndexVector) { if (i == index1) { it1IsDonut = true; break; } } //if not jet, check now if (!it1IsDonut) { bool it2IsDonut = false; unsigned int index2 = 0; for (it2 = aRingVector.begin(); it2 != aRingVector.end(); ++it2, index2++) { if (it2 != it1) { BoostMask::BoostPolygon2D polygon2D; if (::boost::geometry::within(*it1, *it2)) { ::boost::geometry::append(polygon2D, *it2);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it1, 0);//append a ring to the interior ring it2IsDonut = true; } //if donut else if (::boost::geometry::within(*it2, *it1)) { ::boost::geometry::append(polygon2D, *it1);//append an outer ring to the polygon ::boost::geometry::interior_rings(polygon2D).resize(1);//create an interior ring ::boost::geometry::append(polygon2D, *it2, 0);//append a ring to the interior ring it2IsDonut = true; } if (it2IsDonut) { donutIndexVector.push_back(index1); donutIndexVector.push_back(index2); donutVector.push_back(polygon2D);//store donut polygon break;//Only store the first donut! } } } } } //Store no donut polygon to boostPolygonVector index1 = 0; for (it1 = aRingVector.begin(); it1 != aRingVector.end(); ++it1, index1++) { bool it1IsDonut = false; //check if the ring is the outer or inner of a donut for (unsigned int i : donutIndexVector) { if (i == index1) { it1IsDonut = true; break; } } if (!it1IsDonut) { BoostMask::BoostPolygon2D polygon2D; ::boost::geometry::append(polygon2D, *it1); boostPolygonVector.push_back(polygon2D);//insert the ring, which is not a part of donut } } //Append donut polygon to boostPolygonVector BoostMask::BoostPolygonVector::iterator itDonut; for (itDonut = donutVector.begin(); itDonut != donutVector.end(); ++itDonut) { boostPolygonVector.push_back(*itDonut);//append donuts } return boostPolygonVector; } bool BoostMask::calcVoxelizationThickness(double& aThickness) const { if (_voxelizationMap->size() <= 1) { aThickness = 1; return true; } double thickness = 0; auto it = _voxelizationMap->cbegin(); auto it2 = ++_voxelizationMap->cbegin(); for (; it != _voxelizationMap->cend() && it2 != _voxelizationMap->cend(); ++it, ++it2) { if (thickness == 0) { thickness = it2->first - it->first; } else { double curThickness = it2->first - it->first; //if not homogeneous (leave out double imprecisions), return false if (std::abs(thickness-curThickness)>errorConstant) { //return false; std::cout << "Two polygons are far from each other?" << std::endl; } } } if (thickness != 0) { aThickness = thickness; } else { aThickness = 1; } return true; } } } } diff --git a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp index 8bff3c7..818136c 100644 --- a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp +++ b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.cpp @@ -1,150 +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. // //------------------------------------------------------------------------ + #include "rttbBoostMaskGenerateMaskVoxelListThread.h" #include "rttbInvalidParameterException.h" -#include namespace rttb { namespace masks { namespace boost { BoostMaskGenerateMaskVoxelListThread::BoostMaskGenerateMaskVoxelListThread( const VoxelIndexVector& aGlobalBoundingBox, core::GeometricInfo::Pointer aGeometricInfo, BoostArrayMapPointer aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, MaskVoxelListPointer aMaskVoxelList, bool strictVoxelization, - ::boost::shared_ptr<::boost::shared_mutex> aMutex) : + ::boost::shared_ptr aMutex) : _globalBoundingBox(aGlobalBoundingBox), _geometricInfo(aGeometricInfo), _voxelizationMap(aVoxelizationMap), _voxelizationThickness(aVoxelizationThickness), _beginSlice(aBeginSlice), _endSlice(aEndSlice), _resultMaskVoxelList(aMaskVoxelList), _strictVoxelization(strictVoxelization), _mutex(aMutex) {} void BoostMaskGenerateMaskVoxelListThread::operator()() { rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); unsigned int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; unsigned int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; std::vector maskVoxelsInThread; for (unsigned int indexZ = _beginSlice; indexZ < _endSlice; ++indexZ) { //calculate weight vector std::map weightVectorForZ; calcWeightVector(indexZ, weightVectorForZ); //For each x,y, calc sum of all voxelization plane, use weight vector for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) { for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = x + minIndex[0]; currentIndex[1] = y + minIndex[1]; currentIndex[2] = indexZ; rttb::VoxelGridID gridID; _geometricInfo->convert(currentIndex, gridID); double volumeFraction = 0; auto it = weightVectorForZ.cbegin(); auto itMap = _voxelizationMap->cbegin(); for (; it != weightVectorForZ.cend() && itMap != _voxelizationMap->cend(); ++it, ++itMap) { double weight = it->second; if (weight > 0){ BoostArray2DPointer voxelizationArray = itMap->second; //calc sum of all voxelization plane, use weight volumeFraction += (*voxelizationArray)[x][y] * weight; } } if (volumeFraction > 1 && ((volumeFraction - 1) <= errorConstant || !_strictVoxelization)) { volumeFraction = 1; } else if (volumeFraction < 0 || (volumeFraction - 1) > errorConstant) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } //insert mask voxel if volumeFraction > 0 if (volumeFraction > 0) { core::MaskVoxel maskVoxelPtr = core::MaskVoxel(gridID, volumeFraction); maskVoxelsInThread.push_back(maskVoxelPtr); } } } } - ::boost::unique_lock<::boost::shared_mutex> lock(*_mutex); - _resultMaskVoxelList->insert(_resultMaskVoxelList->end(), maskVoxelsInThread.begin(), maskVoxelsInThread.end()); + std::unique_lock lock(*_mutex); + _resultMaskVoxelList->insert(_resultMaskVoxelList->end(), maskVoxelsInThread.begin(), maskVoxelsInThread.end()); } void BoostMaskGenerateMaskVoxelListThread::calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const { double indexZMin = aIndexZ - 0.5; double indexZMax = aIndexZ + 0.5; for (auto & it : *_voxelizationMap) { double voxelizationPlaneIndexMin = it.first - 0.5 * _voxelizationThickness; double voxelizationPlaneIndexMax = it.first + 0.5 * _voxelizationThickness; double weight = 0; if ((voxelizationPlaneIndexMin < indexZMin) && (voxelizationPlaneIndexMax > indexZMin)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = voxelizationPlaneIndexMax - indexZMin; } else { weight = 1; } } else if ((voxelizationPlaneIndexMin >= indexZMin) && (voxelizationPlaneIndexMin < indexZMax)) { if (voxelizationPlaneIndexMax < indexZMax) { weight = _voxelizationThickness; } else { weight = indexZMax - voxelizationPlaneIndexMin; } } weightVector.insert(std::pair(it.first, weight)); } } } } } diff --git a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h index 1f42945..ec1e7e4 100644 --- a/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h +++ b/code/masks/rttbBoostMaskGenerateMaskVoxelListThread.h @@ -1,84 +1,86 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #ifndef __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H #define __BOOST_MASK_GENERATE_MASK_VOXEL_LIST_H #include "rttbBaseType.h" #include "rttbGeometricInfo.h" #include "rttbMaskAccessorInterface.h" +#include +#include + #include #include -#include namespace rttb { namespace masks { namespace boost { /*! @class BoostMaskGenerateMaskVoxelListThread * */ class BoostMaskGenerateMaskVoxelListThread { public: typedef ::boost::multi_array BoostArray2D; using BoostArray2DPointer = ::boost::shared_ptr; typedef ::boost::shared_ptr > BoostArrayMapPointer; using VoxelIndexVector = std::vector; using MaskVoxelListPointer = core::MaskAccessorInterface::MaskVoxelListPointer; BoostMaskGenerateMaskVoxelListThread(const VoxelIndexVector& aGlobalBoundingBox, core::GeometricInfo::Pointer aGeometricInfo, BoostArrayMapPointer aVoxelizationMap, double aVoxelizationThickness, unsigned int aBeginSlice, unsigned int aEndSlice, - MaskVoxelListPointer aMaskVoxelList, bool strictVoxelization, ::boost::shared_ptr<::boost::shared_mutex> aMutex); + MaskVoxelListPointer aMaskVoxelList, bool strictVoxelization, ::boost::shared_ptr aMutex); void operator()(); private: VoxelIndexVector _globalBoundingBox; core::GeometricInfo::Pointer _geometricInfo; BoostArrayMapPointer _voxelizationMap; bool _strictVoxelization=true; //(for example, the first contour has the double grid index 0.1, the second 0.3, the third 0.5, then the thickness is 0.2) double _voxelizationThickness; unsigned int _beginSlice; /*! @brief _endSlice is the first index not to be processed (like a end iterator) */ unsigned int _endSlice; MaskVoxelListPointer _resultMaskVoxelList; - ::boost::shared_ptr<::boost::shared_mutex> _mutex; + ::boost::shared_ptr _mutex; /*! @brief For each dose grid index z, calculate the weight vector for each structure contour */ void calcWeightVector(const rttb::VoxelGridID& aIndexZ, std::map& weightVector) const; }; } } } #endif \ No newline at end of file diff --git a/code/masks/rttbBoostMaskVoxelizationThread.cpp b/code/masks/rttbBoostMaskVoxelizationThread.cpp index 693898b..399d406 100644 --- a/code/masks/rttbBoostMaskVoxelizationThread.cpp +++ b/code/masks/rttbBoostMaskVoxelizationThread.cpp @@ -1,168 +1,167 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #include "rttbBoostMaskVoxelizationThread.h" #include "rttbInvalidParameterException.h" #include -#include #include namespace rttb { namespace masks { namespace boost { BoostMaskVoxelizationThread::BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, - const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex, bool strict) : _geometryCoordinateBoostPolygonMap(APolygonMap), + const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr aMutex, bool strict) : _geometryCoordinateBoostPolygonMap(APolygonMap), _globalBoundingBox(aGlobalBoundingBox), _resultVoxelization(anArrayMap), _mutex(aMutex), _strict(strict) { } void BoostMaskVoxelizationThread::operator()() { rttb::VoxelGridIndex3D minIndex = _globalBoundingBox.at(0); rttb::VoxelGridIndex3D maxIndex = _globalBoundingBox.at(1); const unsigned int globalBoundingBoxSize0 = maxIndex[0] - minIndex[0] + 1; const unsigned int globalBoundingBoxSize1 = maxIndex[1] - minIndex[1] + 1; std::map > voxelizationMapInThread; for (auto & it : _geometryCoordinateBoostPolygonMap) { BoostArray2D maskArray(::boost::extents[globalBoundingBoxSize0][globalBoundingBoxSize1]); BoostPolygonVector boostPolygonVec = it.second; for (unsigned int x = 0; x < globalBoundingBoxSize0; ++x) { for (unsigned int y = 0; y < globalBoundingBoxSize1; ++y) { rttb::VoxelGridIndex3D currentIndex; currentIndex[0] = x + minIndex[0]; currentIndex[1] = y + minIndex[1]; currentIndex[2] = 0; //Get intersection polygons of the dose voxel and the structure BoostPolygonDeque polygons = getIntersections(currentIndex, boostPolygonVec); //Calc areas of all intersection polygons double volumeFraction = calcArea(polygons); volumeFraction = correctForErrorAndStrictness(volumeFraction, _strict); if (volumeFraction < 0 || volumeFraction > 1 ) { throw rttb::core::InvalidParameterException("Mask calculation failed! The volume fraction should >= 0 and <= 1!"); } maskArray[x][y] = volumeFraction; } } voxelizationMapInThread.insert(std::pair(it.first, ::boost::make_shared(maskArray))); } //insert gathered values into voxelization map - ::boost::unique_lock<::boost::shared_mutex> lock(*_mutex); + std::unique_lock lock(*_mutex); _resultVoxelization->insert(voxelizationMapInThread.begin(), voxelizationMapInThread.end()); } /*Get intersection polygons of the contour and a voxel polygon*/ BoostMaskVoxelizationThread::BoostPolygonDeque BoostMaskVoxelizationThread::getIntersections( const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons) { BoostMaskVoxelizationThread::BoostPolygonDeque polygonDeque; BoostRing2D voxelPolygon = get2DContour(aVoxelIndex3D); ::boost::geometry::correct(voxelPolygon); BoostPolygonVector::const_iterator it; for (it = intersectionSlicePolygons.begin(); it != intersectionSlicePolygons.end(); ++it) { BoostPolygon2D contour = *it; ::boost::geometry::correct(contour); BoostPolygonDeque intersection; ::boost::geometry::intersection(voxelPolygon, contour, intersection); polygonDeque.insert(polygonDeque.end(), intersection.begin(), intersection.end()); } return polygonDeque; } BoostMaskVoxelizationThread::BoostRing2D BoostMaskVoxelizationThread::get2DContour( const rttb::VoxelGridIndex3D& aVoxelGrid3D) { BoostRing2D polygon; BoostPoint2D point1(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] - 0.5); ::boost::geometry::append(polygon, point1); BoostPoint2D point2(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] - 0.5); ::boost::geometry::append(polygon, point2); BoostPoint2D point3(aVoxelGrid3D[0] + 0.5, aVoxelGrid3D[1] + 0.5); ::boost::geometry::append(polygon, point3); BoostPoint2D point4(aVoxelGrid3D[0] - 0.5, aVoxelGrid3D[1] + 0.5); ::boost::geometry::append(polygon, point4); ::boost::geometry::append(polygon, point1); return polygon; } /*Calculate the intersection area*/ double BoostMaskVoxelizationThread::calcArea(const BoostPolygonDeque& aPolygonDeque) { double area = 0; BoostPolygonDeque::const_iterator it; for (it = aPolygonDeque.begin(); it != aPolygonDeque.end(); ++it) { area += ::boost::geometry::area(*it); } return area; } double BoostMaskVoxelizationThread::correctForErrorAndStrictness(double volumeFraction, bool strict) const { if (strict){ if (volumeFraction > 1 && (volumeFraction - 1) <= errorConstant) { volumeFraction = 1; } } else { if (volumeFraction > 1){ volumeFraction = 1; } else if (volumeFraction < 0){ volumeFraction = 0; } } return volumeFraction; } } } } diff --git a/code/masks/rttbBoostMaskVoxelizationThread.h b/code/masks/rttbBoostMaskVoxelizationThread.h index b31079f..5b6def9 100644 --- a/code/masks/rttbBoostMaskVoxelizationThread.h +++ b/code/masks/rttbBoostMaskVoxelizationThread.h @@ -1,102 +1,103 @@ // ----------------------------------------------------------------------- // 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. // //------------------------------------------------------------------------ #ifndef __BOOST_MASK_VOXELIZATION_THREAD_H #define __BOOST_MASK_VOXELIZATION_THREAD_H #include +#include +#include #include "rttbBaseType.h" #include #include -#include #include #include namespace rttb { namespace masks { namespace boost { /*! @class BoostMaskGenerateMaskVoxelListThread * */ class BoostMaskVoxelizationThread { public: using BoostPolygon2D = ::boost::geometry::model::polygon< ::boost::geometry::model::d2::point_xy >; using BoostPolygonVector = std::vector;//polygon with or without holes typedef std::map BoostPolygonMap;//map of the z index with the vector of boost 2d polygon using VoxelIndexVector = std::vector; typedef ::boost::multi_array BoostArray2D; using BoostArray2DPointer = ::boost::shared_ptr; typedef ::boost::shared_ptr > BoostArrayMapPointer; /*! @brief Constructor * @param aMutex a mutex for thread-safe handling of the _resultVoxelization * @param strict true means that volumeFractions of <0 and >1 are NOT corrected. Otherwise, they are automatically corrected to 0 or 1, respectively. */ BoostMaskVoxelizationThread(const BoostPolygonMap& APolygonMap, - const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr<::boost::shared_mutex> aMutex, bool strict); + const VoxelIndexVector& aGlobalBoundingBox, BoostArrayMapPointer anArrayMap, ::boost::shared_ptr aMutex, bool strict); void operator()(); private: using BoostPolygonDeque = std::deque; using BoostRing2D = ::boost::geometry::model::ring< ::boost::geometry::model::d2::point_xy >; using BoostPoint2D = ::boost::geometry::model::d2::point_xy; BoostPolygonMap _geometryCoordinateBoostPolygonMap; VoxelIndexVector _globalBoundingBox; - BoostArrayMapPointer _resultVoxelization; - ::boost::shared_ptr<::boost::shared_mutex> _mutex; - bool _strict; + BoostArrayMapPointer _resultVoxelization; + ::boost::shared_ptr _mutex; + bool _strict; /*! @brief Get intersection polygons of the contour and a voxel polygon * @param aVoxelIndex3D The 3d grid index of the voxel * @param intersectionSlicePolygons The polygons of the slice intersecting the voxel * @return Return all intersection polygons of the structure and the voxel */ static BoostPolygonDeque getIntersections(const rttb::VoxelGridIndex3D& aVoxelIndex3D, const BoostPolygonVector& intersectionSlicePolygons); /*! @brief Get the voxel 2d contour polygon in geometry coordinate*/ static BoostRing2D get2DContour(const rttb::VoxelGridIndex3D& aVoxelGrid3D); /*! @brief Calculate the area of all polygons * @param aPolygonDeque The deque of polygons * @return Return the area of all polygons */ static double calcArea(const BoostPolygonDeque& aPolygonDeque); /*! @brief Corrects the volumeFraction * @details the volume fraction is corrected in case of strict=true. Otherwise, it's only corrected for double imprecision * @return The corrected volumeFraction */ double correctForErrorAndStrictness(double volumeFraction, bool strict) const; }; } } } #endif \ No newline at end of file diff --git a/documentation/RTToolbox_build.txt b/documentation/RTToolbox_build.txt index 189236b..3c1e069 100644 --- a/documentation/RTToolbox_build.txt +++ b/documentation/RTToolbox_build.txt @@ -1,209 +1,208 @@ /** \page Build Build instructions \tableofcontents \section tests Compiler tests RTToolbox is currently tested with the following compilers (only x64): - Visual Studio 2013 - Visual Studio 2015 - GCC 5.4 Other compiler versions or other compilers may or may not work. \note The apps can't be compiled with Visual Studio 2010. \section linking Static/Dynamic library support Can be changed with advanced option \c BUILD_SHARED_LIBS \warning building RTToolbox as dynamic library under windows and as static library under Linux is an experimental feature. \section thirdPartyLibraries Needed Third Party Libraries Please load and compile the following third party libraries/tools: \subsection required Required - CMake (version 3.1 or higher) - boost (version 1.64.0 or higher) \subsection optional Optional - DCMTK (with RT support - 3.6.1_20121102 or newer) - if you want DICOM support - ITK (version 4.4 or higher) - if you want DoseInterpolation support with \c itk::Transform or ITK File IO support - MatchPoint (version 0.12 or higher) - if you want DoseInterpolation support with MatchPoint Registration objects) - doxygen - if you want to generate a source code documentation \subsection installing Third party installation Instruction \note To make sure everything runs smoothly, please make sure that all libraries and the RTToolbox are either compiled with \c /MD or \c /MT flags. If third party library packages cannot be found automatically, CMake will ask for them. Please give the location of the binary folder, where the libraries where built for boost, DCMTK, ITK and MatchPoint. \subsubsection boost Boost Build (using the same compiler options as RTToolbox, usually \c STATIC LINKING and x64 architecture). The following components are needed: - \c filesystem, -- \c system, -- \c thread and +- \c system and - \c program_options if you plan to build the apps. \note eventually, it might be needed to add the CMake variable BOOST_LIBRARY dir and set it to the respective library. \subsubsection DCMTK DCMTK For Windows: To compile DCMTK with \c /MD flags (standard for all other libs), you need to patch the CMAKE options of DCMTK (PathToDCMTK\CMake\dcmtkPrepare.cmake), either by replacing \c "/MT" with \c "/MD" or by explicitly replacing lines 253 to 283 with the following lines: \code{.cmake} IF(DCMTK_OVERWRITE_WIN32_COMPILER_FLAGS AND NOT BUILD_SHARED_LIBS) # settings for Microsoft Visual Studio IF(CMAKE_GENERATOR MATCHES "Visual Studio .*") # get Visual Studio Version STRING(REGEX REPLACE "Visual Studio ([0-9]+).*" "\\1" VS_VERSION "${CMAKE_GENERATOR}") # these settings never change even for C or C++ SET(CMAKE_C_FLAGS_DEBUG "/MDd /Z7 /Od") SET(CMAKE_C_FLAGS_RELEASE "/DNDEBUG /MD /O2") SET(CMAKE_C_FLAGS_MINSIZEREL "/DNDEBUG /MD /O2") SET(CMAKE_C_FLAGS_RELWITHDEBINFO "/DNDEBUG /MDd /Z7 /Od") SET(CMAKE_CXX_FLAGS_DEBUG "/MDd /Z7 /Od") SET(CMAKE_CXX_FLAGS_RELEASE "/DNDEBUG /MD /O2") SET(CMAKE_CXX_FLAGS_MINSIZEREL "/DNDEBUG /MD /O2") SET(CMAKE_CXX_FLAGS_RELWITHDEBINFO "/DNDEBUG /MDd /Z7 /Od") # specific settings for the various Visual Studio versions IF(VS_VERSION EQUAL 6) SET(CMAKE_C_FLAGS "/nologo /W3 /GX /Gy /YX") SET(CMAKE_CXX_FLAGS "/nologo /W3 /GX /Gy /YX /Zm500") # /Zm500 increments heap size which is needed on some system to compile templates in dcmimgle ENDIF(VS_VERSION EQUAL 6) IF(VS_VERSION EQUAL 7) SET(CMAKE_C_FLAGS "/nologo /W3 /Gy") SET(CMAKE_CXX_FLAGS "/nologo /W3 /Gy") ENDIF(VS_VERSION EQUAL 7) IF(VS_VERSION GREATER 7) SET(CMAKE_C_FLAGS "/nologo /W3 /Gy /EHsc") SET(CMAKE_CXX_FLAGS "/nologo /W3 /Gy /EHsc") ENDIF(VS_VERSION GREATER 7) ENDIF(CMAKE_GENERATOR MATCHES "Visual Studio .*") ENDIF(DCMTK_OVERWRITE_WIN32_COMPILER_FLAGS AND NOT BUILD_SHARED_LIBS) \endcode Then build DCMTK. \c BUILD_APPS can be switched off. \subsubsection ITK ITK Build ITK. \warning Only use one ITK version consistently throughout all libraries and RTToolbox! Otherwise, linker errors will occur. \warning ensure that compiler enables C++11 features by setting CMAKE_CXX_STANDARD=11 (default for supported compilers) \subsubsection MatchPoint MatchPoint Build MatchPoint. \warning Only use one ITK version consistently throughout all libraries and RTToolbox! Otherwise, linker errors will occur. \warning ensure that compiler enables C++11 features by setting CMAKE_CXX_STANDARD=11 (default for supported compilers) \section installation RT-Toolbox installation Instruction \subsection build Configure/Build RTToolbox \subsubsection Modules Modules - Configure with CMake - Set \c BOOST_INCLUDE_DIR and \c BOOST_DIR to the main boost directory (where "boost_build.jam" is located). Eventually, you have to set BOOST_LIBRARYDIR - Select all packages you like to build (Parameters \c "BUILD_*" ; e.g. \c BUILD_IO_Dicom). - \c BUILD_IO_Dicom: Reading and writing of DICOM-RT files - \c BUILD_IO_HELAX: Reading of Helax DICOM files - \c BUILD_IO_ITK: Generic reading/writing with ITK - \c BUILD_Interpolation: Dose Interpolation - \c BUILD_InterpolationMatchPointTransformation: Dose Interpolation with Match Point registration support. - \c BUILD_Masks: Voxelization support - \c BUILD_Masks_Legacy: Legacy voxelization support. - \c BUILD_Models: Calculation of dosimetrical models like TCP, NTCP etc. Some modules of RT-Toolbox are mandatory (e.g. \c RTTBCore) and build automatically. \warning The \c BUILD_Masks_Legacy module is not supported anymore and may yield to wrong results. \note enabling \c BUILD_All_Modules builds all modules (except Apps, Testing modules and Masks_Legacy). \subsubsection Apps Apps \note The used library ArgumentParsingLib for parsing command line arguments is built automatically. If you want to use the RTTB command line apps, select \c BUILD_Apps. You can select out of five apps: - \c BioModelCalc: calculate the radiobiological effect based on dose - \c DoseAcc: Do dose accumulation - \c DoseMap: Do dose mapping - \c DoseTool: Compute Dose statistics and DVH - \c VoxelizerTool: Voxelize an RTSTRUCT file - Generate CMake configuration - Build RTToolbox and have fun! \subsubsection optionalDependencies Optional third party dependencies - If you select \c BUILD_IO_DICOM, the configuration will ask you for the DCMTK main directory (parameter \c DCMTK_DIR; where you have built DCMTK). - If you select \c BUILD_Interpolation or \c BUILD_IO_ITK, set \c ITK_DIR to the binary directory of ITK. - If you select \c BUILD_InterpolationMatchPointTransformation, set \c MatchPoint_DIR to the binary directory of MatchPoint. \note If you have built DCMTK as out source build you have to also set the CMake parameter \c DCMTK_SOURCE_DIR to the root directory of the DCMTK source. \subsection testRTTB Testing RTToolbox CTest is used as testing framework. See their documentation for general testing questions. \note The used testing library Litmus is build automatically. \warning currently, you have access to testing data only with ssh. That means that a phabricator account and access to RTTB-data repository is mandatory. - Enable \c BUILD_TESTING - Configure with CMake - Enable tests of interest - Generate CMake configuration - Run tests (build \c RUN_TESTS project or call \c ctest in commandline) to ensure that everything is correct. \note enabling \c BUILD_Tester_All builds all test modules. \subsection documentation Documentation Generate the documentation using doxygen and the configuration found in "RTTB_binary_dir/documentation/doxygen.config". \subsection examples Examples - Enable \c BUILD_TESTING - Enable \c BUILD_RTToolbox_Test_examples - Run tests (build \c RUN_TESTS project or call \c ctest in commandline) to also run the examples. These examples (located in testing/examples/) are a good starting point for own projects and show basic functionality of RTToolbox like DVH generation. **/