diff --git a/CMakeExternals/DCMTK.cmake b/CMakeExternals/DCMTK.cmake index bb31a8de..9a510f85 100644 --- a/CMakeExternals/DCMTK.cmake +++ b/CMakeExternals/DCMTK.cmake @@ -101,6 +101,7 @@ if(NOT DEFINED DCMTK_DIR AND NOT ${CMAKE_PROJECT_NAME}_USE_SYSTEM_${proj}) -DDCMTK_ENABLE_BUILTIN_DICTIONARY:BOOL=ON -DDCMTK_ENABLE_PRIVATE_TAGS:BOOL=ON -DDCMTK_COMPILE_WIN32_MULTITHREADED_DLL:BOOL=ON + -DDCMTK_ENABLE_STL:BOOL=ON DEPENDS ${${proj}_DEPENDENCIES} ) diff --git a/apps/seg/Testing/CMakeLists.txt b/apps/seg/Testing/CMakeLists.txt index 1ead4de9..fccdab32 100644 --- a/apps/seg/Testing/CMakeLists.txt +++ b/apps/seg/Testing/CMakeLists.txt @@ -30,6 +30,10 @@ dcmqi_add_test( --outputDICOM ${MODULE_TEMP_DIR}/liver.dcm ) +# Creates a DICOM segmentation file that has 3 segments: +# - segment for liver (DICOM Segment Number 1) +# - segment for spine (DICOM Segment Number 2) +# - segment for heart (DICOM Segment Number 3) dcmqi_add_test( NAME ${itk2dcm}_makeSEG_multiple_segment_files MODULE_NAME ${MODULE_NAME} @@ -50,6 +54,7 @@ dcmqi_add_test( --outputDICOM ${MODULE_TEMP_DIR}/liver_heart_seg_reordered.dcm ) + find_program(DCIODVFY_EXECUTABLE dciodvfy) if(EXISTS ${DCIODVFY_EXECUTABLE}) @@ -139,6 +144,38 @@ dcmqi_add_test( ${itk2dcm}_makeSEG_multiple_segment_files_reordered ) + # Reads a DICOM segmentation file that has 3 segments (liver, spine, heart - in this order). + # Heart and liver segments overlap. + # The goal is to export these segments to NRRD+JSON. Since the liver and heart segments overlap, + # one output file will contain 2 segments (liver and spine) and the other will contain only + # the heart. This will be reflected in the JSON metadata files with two entries in the outer + # array (one for each NRRD file), one with 2 segments (liver,spine) and one with 1 segment (heart). + dcmqi_add_test( + NAME ${dcm2itk}_makeNRRD_merged_segment_file + MODULE_NAME ${MODULE_NAME} + COMMAND ${SEM_LAUNCH_COMMAND} $ + --compare ${BASELINE}/liver_spine_seg.nrrd ${MODULE_TEMP_DIR}/makeNRRD_merged_segment_file-1.nrrd + --compare ${BASELINE}/heart_seg.nrrd ${MODULE_TEMP_DIR}/makeNRRD_merged_segment_file-2.nrrd + ${dcm2itk}Test + --inputDICOM ${MODULE_TEMP_DIR}/liver_heart_seg.dcm + --outputDirectory ${MODULE_TEMP_DIR} + --prefix makeNRRD_merged_segment_file + --mergeSegments + TEST_DEPENDS + ${itk2dcm}_makeSEG_multiple_segment_files + ) + + # Compare expected JSON output coming from makeNRRD_merged_segment_file test. + dcmqi_add_test( + NAME ${dcm2itk}_makeNRRD_merged_segment_file_JSON + MODULE_NAME ${MODULE_NAME} + COMMAND python ${CMAKE_SOURCE_DIR}/util/comparejson.py + ${CMAKE_SOURCE_DIR}/doc/examples/seg-example_multiple_segments_merged.json + ${MODULE_TEMP_DIR}/makeNRRD_merged_segment_file-meta.json + TEST_DEPENDS + ${dcm2itk}_makeNRRD_merged_segment_file + ) + dcmqi_add_test( NAME seg_meta_roundtrip MODULE_NAME ${MODULE_NAME} @@ -160,8 +197,6 @@ dcmqi_add_test( ${dcm2itk}_makeNRRD_multiple_segment_files ) - - set(TEST_SEG_SIZES 24x38x3 23x38x3) foreach(seg_size ${TEST_SEG_SIZES}) diff --git a/apps/seg/itkimage2segimage.cxx b/apps/seg/itkimage2segimage.cxx index e3e0b329..8c3c7fbf 100644 --- a/apps/seg/itkimage2segimage.cxx +++ b/apps/seg/itkimage2segimage.cxx @@ -1,9 +1,9 @@ // CLP includes +#include "dcmqi/Itk2DicomConverter.h" #include "itkimage2segimageCLP.h" // DCMQI includes #undef HAVE_SSTREAM // Avoid redefinition warning -#include "dcmqi/ImageSEGConverter.h" #include "dcmqi/internal/VersionConfigure.h" // DCMTK includes @@ -122,7 +122,7 @@ int main(int argc, char *argv[]) } try { - DcmDataset* result = dcmqi::ImageSEGConverter::itkimage2dcmSegmentation(dcmDatasets, segmentations, metadata, skipEmptySlices); + DcmDataset* result = dcmqi::Itk2DicomConverter::itkimage2dcmSegmentation(dcmDatasets, segmentations, metadata, skipEmptySlices); if (result == NULL){ std::cerr << "ERROR: Conversion failed." << std::endl; diff --git a/apps/seg/segimage2itkimage.cxx b/apps/seg/segimage2itkimage.cxx index ff076f51..0cae329a 100644 --- a/apps/seg/segimage2itkimage.cxx +++ b/apps/seg/segimage2itkimage.cxx @@ -1,15 +1,17 @@ // CLP includes #include "segimage2itkimageCLP.h" +#include // DCMQI includes #undef HAVE_SSTREAM // Avoid redefinition warning -#include "dcmqi/ImageSEGConverter.h" +#include "dcmqi/Dicom2ItkConverter.h" #include "dcmqi/internal/VersionConfigure.h" // DCMTK includes #include typedef dcmqi::Helper helper; +typedef itk::ImageFileWriter WriterType; int main(int argc, char *argv[]) @@ -23,6 +25,11 @@ int main(int argc, char *argv[]) dcmtk::log4cplus::BasicConfigurator::doConfigure(); } + if (mergeSegments && outputType != "nrrd") { + std::cerr << "ERROR: mergeSegments option is only supported when output format is NRRD!" << std::endl; + return EXIT_FAILURE; + } + if(helper::isUndefinedOrPathDoesNotExist(inputSEGFileName, "Input DICOM file") || helper::isUndefinedOrPathDoesNotExist(outputDirName, "Output directory")) return EXIT_FAILURE; @@ -34,23 +41,39 @@ int main(int argc, char *argv[]) DcmDataset* dataset = sliceFF.getDataset(); try { - pair , string> result = dcmqi::ImageSEGConverter::dcmSegmentation2itkimage(dataset); + dcmqi::Dicom2ItkConverter converter; + std::string metaInfo; + OFCondition result = converter.dcmSegmentation2itkimage(dataset, metaInfo, mergeSegments); + if (result.bad()) + { + std::cerr << "ERROR: Failed to convert DICOM SEG to ITK image: " << result.text() << std::endl; + return EXIT_FAILURE; + } string outputPrefix = prefix.empty() ? "" : prefix + "-"; string fileExtension = dcmqi::Helper::getFileExtensionFromType(outputType); - - for(map::const_iterator sI=result.first.begin();sI!=result.first.end();++sI){ - typedef itk::ImageFileWriter WriterType; + itk::SmartPointer itkImage = converter.begin(); + size_t fileIndex = 1; + while (itkImage) + { stringstream imageFileNameSStream; - - imageFileNameSStream << outputDirName << "/" << outputPrefix << sI->first << fileExtension; - - WriterType::Pointer writer = WriterType::New(); - writer->SetFileName(imageFileNameSStream.str().c_str()); - writer->SetInput(sI->second); - writer->SetUseCompression(1); - writer->Update(); + cout << "Writing itk image to " << outputDirName << "/" << outputPrefix << fileIndex << fileExtension; + imageFileNameSStream << outputDirName << "/" << outputPrefix << fileIndex << fileExtension; + + try { + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName(imageFileNameSStream.str().c_str()); + writer->SetInput(itkImage); + writer->SetUseCompression(1); + writer->Update(); + cout << " ... done" << endl; + } catch (itk::ExceptionObject & error) { + std::cerr << "fatal ITK error: " << error << std::endl; + return EXIT_FAILURE; + } + itkImage = converter.next(); + fileIndex++; } stringstream jsonOutput; @@ -58,7 +81,7 @@ int main(int argc, char *argv[]) ofstream outputFile; outputFile.open(jsonOutput.str().c_str()); - outputFile << result.second; + outputFile << metaInfo; outputFile.close(); return EXIT_SUCCESS; diff --git a/apps/seg/segimage2itkimage.xml b/apps/seg/segimage2itkimage.xml index 80ca3f5f..af3a72ea 100644 --- a/apps/seg/segimage2itkimage.xml +++ b/apps/seg/segimage2itkimage.xml @@ -66,6 +66,15 @@ Display more verbose output, useful for troubleshooting. + + mergeSegments + + input + mergeSegments + false + Save all segments into a single file. When segments are non-overlapping, output is a single 3D file. If overlapping, single 4D following conventions of 3D Slicer segmentations format. Only supported when the output is NRRD for now. + + diff --git a/data/segmentations/liver_spine_seg.nrrd b/data/segmentations/liver_spine_seg.nrrd new file mode 100644 index 00000000..d536c5cf Binary files /dev/null and b/data/segmentations/liver_spine_seg.nrrd differ diff --git a/doc/examples/seg-example_multiple_segments_merged.json b/doc/examples/seg-example_multiple_segments_merged.json new file mode 100644 index 00000000..36588411 --- /dev/null +++ b/doc/examples/seg-example_multiple_segments_merged.json @@ -0,0 +1,68 @@ +{ + "BodyPartExamined" : "", + "ClinicalTrialCoordinatingCenterName" : "BWH", + "ClinicalTrialSeriesID" : "Session1", + "ClinicalTrialTimePointID" : "1", + "ContentCreatorName" : "Doe^John", + "InstanceNumber" : "1", + "SeriesDescription" : "Segmentation", + "SeriesNumber" : "300", + "segmentAttributes" : [ + [ + { + "SegmentAlgorithmName" : "SlicerEditor", + "SegmentAlgorithmType" : "SEMIAUTOMATIC", + "SegmentDescription" : "Liver Segmentation", + "SegmentLabel" : "Liver", + "SegmentedPropertyCategoryCodeSequence" : { + "CodeMeaning" : "Tissue", + "CodeValue" : "85756007", + "CodingSchemeDesignator" : "SCT" + }, + "SegmentedPropertyTypeCodeSequence" : { + "CodeMeaning" : "Liver", + "CodeValue" : "10200004", + "CodingSchemeDesignator" : "SCT" + }, + "labelID" : 1, + "recommendedDisplayRGBValue" : [ 220, 129, 101 ] + }, + { + "SegmentAlgorithmType" : "MANUAL", + "SegmentDescription" : "Anatomical Structure", + "SegmentLabel" : "Thoracic spine", + "SegmentedPropertyCategoryCodeSequence" : { + "CodeMeaning" : "Anatomical Structure", + "CodeValue" : "123037004", + "CodingSchemeDesignator" : "SCT" + }, + "SegmentedPropertyTypeCodeSequence" : { + "CodeMeaning" : "Thoracic spine", + "CodeValue" : "122495006", + "CodingSchemeDesignator" : "SCT" + }, + "labelID" : 2, + "recommendedDisplayRGBValue" : [ 226, 202, 134 ] + } + ], + [ + { + "SegmentAlgorithmType" : "MANUAL", + "SegmentDescription" : "Anatomical Structure", + "SegmentLabel" : "Heart", + "SegmentedPropertyCategoryCodeSequence" : { + "CodeMeaning" : "Anatomical Structure", + "CodeValue" : "123037004", + "CodingSchemeDesignator" : "SCT" + }, + "SegmentedPropertyTypeCodeSequence" : { + "CodeMeaning" : "Heart", + "CodeValue" : "80891009", + "CodingSchemeDesignator" : "SCT" + }, + "labelID" : 3, + "recommendedDisplayRGBValue" : [ 206, 110, 84 ] + } + ] + ] +} diff --git a/docker/Makefile b/docker/Makefile index daf7773c..59ad31cd 100644 --- a/docker/Makefile +++ b/docker/Makefile @@ -44,17 +44,19 @@ prereq.build_dicom3tools: prereq.pull_dockcross && ( $(TMP)/dockcross make -j$(grep -c processor /proc/cpuinfo) World > /dev/null 2>&1 )\ || echo "Reusing " -# Download and install "ajv" +# Download and install "ajv" prereq.install_npm_packages: prereq.pull_dockcross mkdir -p $(ROOT_DIR)/$(BUILD_DIR) + echo "Root dir: $(ROOT_DIR)" + echo "Build dir: $(BUILD_DIR)" # If needed, download node cd $(ROOT_DIR)/$(BUILD_DIR) && \ - [ ! -e node-v6.9.5-linux-x64 ] && \ - wget --no-check-certificate https://nodejs.org/dist/v6.9.5/node-v6.9.5-linux-x64.tar.xz && \ - tar xf node-v6.9.5-linux-x64.tar.xz || true + [ ! -e node-v12.19.1-linux-x64 ] && \ + wget --no-check-certificate https://nodejs.org/dist/v12.19.1/node-v12.19.1-linux-x64.tar.xz && \ + tar xf node-v12.19.1-linux-x64.tar.xz || true # Install tools required to run DCMQI "doc" tests cd $(ROOT_DIR) && \ - $(TMP)/dockcross bash -c "export PATH=/work/build/node-v6.9.5-linux-x64/bin:$$PATH && sudo npm install ajv-cli@3.3.0 -g" + $(TMP)/dockcross bash -c "export PATH=/work/build/node-v12.19.1-linux-x64/bin:$$PATH && sudo npm install ajv-cli@3.3.0 -g" # Configure, build and package dcmqi.generate_package: prereq.pull_dockcross diff --git a/include/dcmqi/ConverterBase.h b/include/dcmqi/ConverterBase.h index 37c7abfe..ed9df404 100644 --- a/include/dcmqi/ConverterBase.h +++ b/include/dcmqi/ConverterBase.h @@ -37,7 +37,9 @@ using namespace std; typedef short ShortPixelType; +typedef unsigned char CharPixelType; typedef itk::Image ShortImageType; +typedef itk::Image CharImageType; typedef itk::ImageFileReader ShortReaderType; namespace dcmqi { @@ -219,6 +221,21 @@ namespace dcmqi { sort(originDistances.begin(), originDistances.end()); sliceSpacing = fabs(originDistances[0]-originDistances[1]); + if (sliceSpacing == 0) + { + cout << "Slice spacing is zero, trying to get/use it from DICOM file instead" << endl; + // Get Slice Spacing as defined in Pixel Measures FG + OFBool isPerFrame; + FGPixelMeasures *pixelMeasures = OFstatic_cast(FGPixelMeasures*, + fgInterface.get(0, DcmFGTypes::EFG_PIXELMEASURES, isPerFrame)); + if(!pixelMeasures){ + cerr << "Canont get Slice Spacing, Pixel Measures FG is missing!" << endl; + return EXIT_FAILURE; + } + if(pixelMeasures->getSpacingBetweenSlices(sliceSpacing,0).good()){ + cout << "Using Slice Spacing (from DICOM file): " << sliceSpacing << endl; + } + } // for(size_t i=1;isecond>1) overlappingFramesCnt++; } - cout << "Total frames with unique IPP: " << originDistances.size() << endl; cout << "Total overlapping frames: " << overlappingFramesCnt << endl; } @@ -244,6 +260,9 @@ namespace dcmqi { sliceSpacing = 1.0; } cout << "Origin: " << imageOrigin << endl; + cout << "Slice extent: " << sliceExtent << endl; + cout << "Slice spacing: " << sliceSpacing << endl; + return 0; } diff --git a/include/dcmqi/Dicom2ItkConverter.h b/include/dcmqi/Dicom2ItkConverter.h new file mode 100644 index 00000000..ba26814f --- /dev/null +++ b/include/dcmqi/Dicom2ItkConverter.h @@ -0,0 +1,188 @@ +#ifndef DCMQI_SEGMENTATION_CONVERTER_H +#define DCMQI_SEGMENTATION_CONVERTER_H + +// DCMTK includes +#include +#include +#include +#include +#include +#include + +// ITK includes +#include +#include +#include +#include +#include +#include + +// DCMQI includes +#include "OverlapUtil.h" +#include "dcmqi/ConverterBase.h" +#include "dcmqi/JSONSegmentationMetaInformationHandler.h" +#include "dcmqi/OverlapUtil.h" + +using namespace std; + +typedef itk::LabelImageToLabelMapFilter LabelToLabelMapFilterType; + +namespace dcmqi +{ + +/** + * @brief The ImageSEGConverter class provides methods to convert from DICOM Segmentation objects to itk images. + * It should be used as follows: + * - Call dcmSegmentation2itkimage() to start the conversion of DICOM + * Segmentation object into a map of itk images. This will return to you + * the JSON meta information of the conversion. + * - Call begin() to get the first ITK image result of the conversion. + * Note that if conversion fails, this can return a null pointer. + * - Call next() to get the next ITK image result of the conversion, until + * it returns a null pointer. + */ +class Dicom2ItkConverter : public ConverterBase +{ + +public: + Dicom2ItkConverter(); + + /** + * @brief Converts a DICOM Segmentation object into a map of itk images and metadata. + * + * @param segDataset A pointer to the DICOM Segmentation object to be converted. + * @param mergeSegments A boolean indicating whether to merge segments during the conversion. + * Defaults to false. + * @return A pair containing the resulting map of itk images and the metadata. + */ + OFCondition dcmSegmentation2itkimage(DcmDataset* segDataset, std::string& metaInfo, const bool mergeSegments = false); + + /** Get first ITK image result of conversion, or null pointer if conversion failed + * @return Shared pointer to first ITK image resulting from the conversion + */ + itk::SmartPointer begin(); + + /** Get next ITK image result of conversion, or nuĺl pointer if no more results + * @return Shared pointer to next ITK image resulting from the conversion, or null + */ + itk::SmartPointer next(); + + /** Get JSON the metadata of the conversion + * @return Metadata of the conversion + */ + JSONSegmentationMetaInformationHandler getMetaInformation(); + +protected: + /** Internal result loop, produced one result at a time (or null) + * @return Shared pointer to first/next ITK image resulting from the conversion + */ + itk::SmartPointer nextResult(); + + /** + * @brief Converts a DICOM Segmentation object into a map of itk images. + * + * @param mergeSegments A boolean indicating whether to merge segments during the conversion. + * Defaults to false. + * @return A map of itk images resulting from the conversion. + */ + + OFCondition dcmSegmentation2itkimage(const bool mergeSegments = false); + + /** + * @brief Populates the metadata of a DICOM Segmentation object from a DICOM dataset. + * + * @param segDataset A pointer to the DICOM dataset containing the metadata. + */ + void populateMetaInformationFromDICOM(DcmDataset* segDataset); + + /** + * Helper method that uses the OverlapUtil class to retrieve non-overlapping segment + * groups (by their segment numbers) from the DICOM Segmentation object. + * If mergeSegments is false, all segments are returned within a single group. + * + * @param mergeSegments Whether or not to merge overlapping segments into the same group. + * @param segmentGroups The resulting segment groups. + * @return EC_Normal if successful, error otherwise + */ + OFCondition getNonOverlappingSegmentGroups(const bool mergeSegments, OverlapUtil::SegmentGroups& segmentGroups); + + /** + * Extract basic image info like directions, origin, spacing and image region. + * @return EC_Normal if successful, error otherwise + */ + OFCondition extractBasicSegmentationInfo(); + + /** + * Allocate an ITK image with the same size and spacing as the DICOM image. + * This is used for each slice that is to be inserted into the ITK space. + * @return EC_Normal if successful, error otherwise + */ + ShortImageType::Pointer allocateITKImage(); + + /** + * Allocate an ITK image based on the given template. + * @param imageTemplate The template to use for the new image. + * @return EC_Normal if successful, error otherwise + */ + ShortImageType::Pointer allocateITKImageDuplicate(ShortImageType::Pointer imageTemplate); + + /** + * Get the ITK image origin for the given frame. + * @param frameNo The (DICOM) frame number to get the origin for. + * @param origin The resulting origin. + * @return EC_Normal if successful, error otherwise + */ + OFCondition getITKImageOrigin(const Uint32 frameNo, ShortImageType::PointType& origin); + + /** + * Collect segment metadata for the given frame that will later go into the + * accompanying JSON segmentation description. + * @param segmentGroup The segment group number + * (i.e. uniquely identifying the NRRD output file in the end). + * @param segmentNumber The DICOM segment number. + * @param origin The resulting origin. + * @return EC_Normal if successful, error otherwise + */ + OFCondition addSegmentMetadata(const size_t segmentGroup, const Uint16 segmentNumber); + + OFCondition createMetaInfo(); + + /// The segmentation object, guarded by unique pointer + OFunique_ptr m_segDoc; + + /// Image direction in ITK speak + ShortImageType::DirectionType m_direction; + /// Computed image slice spacing + double m_computedSliceSpacing; + /// Computed volume extent + double m_computedVolumeExtent; + /// Slice direction in ITK speak + vnl_vector m_sliceDirection; + + /// Image origin in ITK speak + ShortImageType::PointType m_imageOrigin; + /// Image slice spacing in ITK speak + ShortImageType::SpacingType m_imageSpacing; + /// Image size in ITK speak + ShortImageType::SizeType m_imageSize; + /// Image region in ITK speak + ShortImageType::RegionType m_imageRegion; + + /// Segment meta information handler for the JSON segmentation description + JSONSegmentationMetaInformationHandler m_metaInfo; + + /// Segment groups, used to iterate over results + OverlapUtil::SegmentGroups m_segmentGroups; + + /// Internal iterator to current segment group, used while iterating + /// over results using begin() and next() + OverlapUtil::SegmentGroups::iterator m_groupIterator; + + /// OverlapUtil instance used by this class, used in DICOM segmentation + /// to itk conversion + OverlapUtil m_overlapUtil; +}; + +} + +#endif // DCMQI_SEGMENTATION_CONVERTER_H diff --git a/include/dcmqi/ImageSEGConverter.h b/include/dcmqi/ImageSEGConverter.h deleted file mode 100644 index 01bde85b..00000000 --- a/include/dcmqi/ImageSEGConverter.h +++ /dev/null @@ -1,50 +0,0 @@ -#ifndef DCMQI_SEGMENTATION_CONVERTER_H -#define DCMQI_SEGMENTATION_CONVERTER_H - -// DCMTK includes -#include -#include -#include -#include -#include -#include - -// ITK includes -#include -#include -#include -#include -#include - -// DCMQI includes -#include "dcmqi/ConverterBase.h" -#include "dcmqi/JSONSegmentationMetaInformationHandler.h" - -using namespace std; - - -typedef itk::LabelImageToLabelMapFilter LabelToLabelMapFilterType; - -namespace dcmqi { - - class ImageSEGConverter : public ConverterBase { - - public: - static DcmDataset* itkimage2dcmSegmentation(vector dcmDatasets, - vector segmentations, - const string &metaData, - bool skipEmptySlices=true); - - - static pair, string> dcmSegmentation2itkimage(DcmDataset *segDataset); - static map dcmSegmentation2itkimage( - DcmSegmentation *segdoc, - JSONSegmentationMetaInformationHandler *metaInfo=NULL); - - static void populateMetaInformationFromDICOM(DcmDataset *segDataset, DcmSegmentation *segdoc, - JSONSegmentationMetaInformationHandler &metaInfo); - }; - -} - -#endif //DCMQI_SEGMENTATION_CONVERTER_H diff --git a/include/dcmqi/Itk2DicomConverter.h b/include/dcmqi/Itk2DicomConverter.h new file mode 100644 index 00000000..90f86387 --- /dev/null +++ b/include/dcmqi/Itk2DicomConverter.h @@ -0,0 +1,56 @@ +#ifndef DCMQI_ITK2DICOM_CONVERTER_H +#define DCMQI_ITK2DICOM_CONVERTER_H + +// DCMTK includes +#include +#include +#include +#include +#include +#include + +// ITK includes +#include +#include +#include +#include +#include +#include + +// DCMQI includes +#include "dcmqi/ConverterBase.h" + + +using namespace std; + +typedef itk::LabelImageToLabelMapFilter LabelToLabelMapFilterType; + +namespace dcmqi { + + /** + * @brief The Itk2DicomConverter class provides methods to convert from itk images to DICOM Segmentation objects. + */ + class Itk2DicomConverter : public ConverterBase { + + public: + + Itk2DicomConverter(); + + /** + * @brief Converts itk images data into a DICOM Segmentation object. + * + * @param dcmDatasets A vector of DICOM datasets with the images that the segmentation is based on. + * @param segmentations A vector of itk images to be converted. + * @param metaData A string containing the metadata to be used for the DICOM Segmentation object. + * @param skipEmptySlices A boolean indicating whether to skip empty slices during the conversion. + * @return A pointer to the resulting DICOM Segmentation object. + */ + static DcmDataset* itkimage2dcmSegmentation(vector dcmDatasets, + vector segmentations, + const string &metaData, + bool skipEmptySlices=true); + }; + +} + +#endif //DCMQI_ITK2DICOM_CONVERTER_H diff --git a/include/dcmqi/JSONSegmentationMetaInformationHandler.h b/include/dcmqi/JSONSegmentationMetaInformationHandler.h index 219b0087..c7d312a6 100644 --- a/include/dcmqi/JSONSegmentationMetaInformationHandler.h +++ b/include/dcmqi/JSONSegmentationMetaInformationHandler.h @@ -35,7 +35,10 @@ namespace dcmqi { void read(); bool write(string filename); - SegmentAttributes* createAndGetNewSegment(unsigned labelID); + // segGroupNumber starts with 0 and refers to the item in segmentsAttributesMappingList. + // if segGroupNumber is invalid, create segment within a newly created segmentation group + // otherwise add segment to the existing segmentation group identified by segGroupNumber + SegmentAttributes* createOrGetSegment(const unsigned int segGroupNumber, unsigned labelID); protected: @@ -46,7 +49,7 @@ namespace dcmqi { void readSegmentAttributes(); - Json::Value createAndGetSegmentAttributes(); + Json::Value createAndGetSegmentAttributesJSON(); }; } diff --git a/include/dcmqi/OverlapUtil.h b/include/dcmqi/OverlapUtil.h new file mode 100644 index 00000000..3116d222 --- /dev/null +++ b/include/dcmqi/OverlapUtil.h @@ -0,0 +1,373 @@ +/* + * + * Copyright (C) 2023, Open Connections GmbH + * All rights reserved. See COPYRIGHT file for details. + * + * This software and supporting documentation were developed by + * + * OFFIS e.V. + * R&D Division Health + * Escherweg 2 + * D-26121 Oldenburg, Germany + * + * + * Module: dcmqi + * + * Author: Michael Onken + * + * Purpose: Interface of class OverlapUtil + * + */ + +#ifndef DCMQI_OVERLAPUTIL_H +#define DCMQI_OVERLAPUTIL_H + +#include "dcmtk/dcmiod/iodtypes.h" +#include "dcmtk/ofstd/ofcond.h" +#include "dcmtk/ofstd/oftypes.h" +#include "dcmtk/ofstd/ofvector.h" +#include + +class DcmSegmentation; + +/// error: frames are not parallel +extern const OFConditionConst SG_EC_FramesNotParallel; + +namespace dcmqi +{ + +/** Class that analyzes the frame-based segment structures of a given segmentation object. + * It provides the following main functionality: + * - Grouping of "physical" frames by their position in space (called "logical frames"), i.e. + * frames that are at the same position in space are grouped together. + * - For every logical frame (i.e. position in space), it lists the segments found at this + * position together with their respective physical frame number + * - Return physical frame numbers for a given segment number + * - Building an overlap matrix that stores for each arbitrary segment pair whether + * they overlap or not. + * - Return groups of segments, that no two overlapping segments will be in the same group. + * This will not necessarily return the optimal solution, but a solution that should be good enough. + * For the method used see getNonOverlappingSegments(). + */ +class OverlapUtil +{ +public: + /// Image Position Patient tuple (x,y,z) + typedef OFVector ImagePosition; // might be defined in respective functional group + + /// Physical Frame number with its respective position + struct FramePositionAndNumber + { + /** Default constructor required for vector initialization + */ + FramePositionAndNumber() + : m_position() + , m_frameNumber(0) + { + } + /** Constructor + * @param pos Image position + * @param num Physical frame number + */ + FramePositionAndNumber(const ImagePosition& pos, const Uint32& num) + : m_position(pos) + , m_frameNumber(num) + { + } + /// Frame position in space + ImagePosition m_position; + /// Physical frame number (number of frame in DICOM object) + Uint32 m_frameNumber; + }; + + /// Physical Frame number with its respective position + typedef OFVector FramePositions; + + /// Logical Frame, represented and defined by various physical frames (numbers) at the same position + typedef OFVector LogicalFrame; + + /// All distinct positions and for each position the physical frame numbers that can be found at it + typedef OFVector DistinctFramePositions; + + /// Lists frames for each segment where segment with index i is represented by the vector at index i, + /// and index 0 is unused. I.e. index i is segment number, value is vector of physical frame numbers. + typedef OFVector> FramesForSegment; + + /// Implements comparision operator to be used for sorting of frame positions, + /// making the sorting order depend on the coordinate given in the constructor + struct ComparePositions + { + /** Constructor, used to configure coordinate position to be used for sorting + * @param c Coordinate position to be used for sorting + */ + ComparePositions(size_t c) + : m_coordinate(c) + { + } + + /** Compare two frames*/ + bool operator()(const FramePositionAndNumber& a, const FramePositionAndNumber& b) const + { + return a.m_position[m_coordinate] < b.m_position[m_coordinate]; + } + /// Coordinate position (0-2, i.e. x,x,z) to be used for sorting + size_t m_coordinate; + }; + + /// Matrix of N x N segment numbers, where N is the number of segments. + /// Value is 1 at x,y if x and y overlap, 0 if they don't overlap, and -1 if not initialized. + typedef OFVector> OverlapMatrix; + + /// Group of non-overlapping segments (each represented by its segment number) + typedef OFVector> SegmentGroups; + + /** Represents a segment number and a logical frame number it is found at + */ + struct SegNumAndFrameNum + { + /** Constructor + * @param s Segment number + * @param f Logical frame number + */ + SegNumAndFrameNum(const Uint16& s, const Uint16 f) + : m_segmentNumber(s) + , m_frameNumber(f) + { + } + /// Segment number as used in DICOM segmentation object (1-n) + Uint16 m_segmentNumber; + /// Logical frame number (number of frame in DistinctFramePositions vector) + Uint16 m_frameNumber; + /** Comparison operator + * @param rhs Right-hand side of comparison + * @return OFTrue if left-hand side is smaller than right-hand side + */ + bool operator<(const SegNumAndFrameNum& rhs) const + { + return m_segmentNumber < rhs.m_segmentNumber; + } + }; + + /// Segments and their phyiscal frame number (inner set), grouped by their + /// respective logical frame number (outer vector) + typedef OFVector> SegmentsByPosition; + + // ------------------------------------------ Methods ------------------------------------------ + + /** Constructor. Use setSegmentationObject() to set the segmentation object to work with. + */ + OverlapUtil(); + + /** Destructor */ + ~OverlapUtil(); + + /** Set the segmentation object to work with and clears all old data. + * TODO: In the future, maybe have DcmSegmentation->getOverlapUtil() to access + * the OverlapUtil object, so that the user does not have to care about + * feeding the segmentation object to the OverlapUtil object. + */ + void setSegmentationObject(DcmSegmentation* seg); + + /** Clears all internal data (except segmentation object reference). + * This should be called whenever the input data (i.e. the underlying) + * DICOM segmentation object changes, before calling any other method. + */ + void clear(); + + /** Get all distinct frame positions and the physical frame numbers at this position + * @param result Resulting vector of distinct frame positions + * @return EC_Normal if successful, error otherwise + */ + OFCondition getFramesByPosition(DistinctFramePositions& result); + + /** Get all segments and their physical frame number, grouped by their respective logical frame number + * @param result Resulting vector of segments grouped by logical frame number + * @return EC_Normal if successful, error otherwise + */ + OFCondition getSegmentsByPosition(SegmentsByPosition& result); + + /** Get phyiscal frames for a specific segment by its segment number + * @param segmentNumber Segment number to get frames for (1..n) + * @param frames Resulting vector of physical frame numbers (first frame is frame 0) + * @return EC_Normal if successful, error otherwise + */ + OFCondition getFramesForSegment(const Uint32 segmentNumber, OFVector& frames); + + /** Returns computed overlap matrix + * @param matrix Resulting overlap matrix + * @return EC_Normal if successful, error otherwise + */ + OFCondition getOverlapMatrix(OverlapMatrix& matrix); + + /** Returns segments grouped together in a way, that no two overlapping + * segments will be in the same group. This method does not necessarily + * returns the optimal solution, but a solution that should be good enough. + * It is guaranteed, that segments in the same group don't overlap. + * + * It is based on the idea of a greedy algorithm that creates a first group + * containing the first segment. Then it goes to the next segment, checks whether + * it fits into the first group with no overlaps (easily checked in overlap matrix) + * and inserts it into that group if no overlaps exists. Otherwise, it creates a + * new group and continues with the next segment (trying to insert it into + * the first group, then second group, otherwise creates third group, and so on). + * @param segmentGroups Resulting vector of segment groups, each listing the + * segment numbers that are in that group + * @return EC_Normal if successful, error otherwise + */ + OFCondition getNonOverlappingSegments(SegmentGroups& segmentGroups); + + /** Prints segments by their position in space + * @param ss The stream to dump to + */ + void printSegmentsByPosition(OFStringStream& ss); + + /** Prints segment overlap matrix to given stream + * @param ss The stream to dump to + */ + void printOverlapMatrix(OFStringStream& ss); + + /** Prints groups of non-overlapping segments (identified by their numbers) + * to given stream + * @param ss The stream to dump to + */ + void printNonOverlappingSegments(OFStringStream& ss); + +protected: + /** Collect all physical frame positions in m_FramePositions + * @return EC_Normal if successful, error otherwise + */ + OFCondition collectPhysicalFramePositions(); + + /** Group physical frame positions into logical positions. This is done by sorting + * frames after *that* position coordinate that in its mean position difference is + * larger than slice thickness * 0.9. Then those frames that are close enough to + * each other (i.e. distance is smaller than slice thickness * 0.01), end up at the + * same logical position (considered a "logical frame") + * TODO: This should probably not use mean values for the coordinates + * since in some cases, the mean difference in a slice coordinate might be close to 0 + * if many frames are at the same position. Instead, the maximum difference, variance or + * something else could be used? + * @return EC_Normal if successful, error otherwise + */ + OFCondition groupFramesByLogicalPosition(); + + /** Builds the overlap matrix, if not already done. + * @return EC_Normal if successful or already existant, error otherwise + */ + OFCondition buildOverlapMatrix(); + + /** Checks if frames are parallel, i.e. if DICOM Image Position Patient is present and + * all frames are parallel to each other (i.e. found in the shared functional group) + * @return EC_Normal if parallel, SG_EC_FramesNotParallel if image orientation is not shared, + * error otherwise + */ + OFCondition ensureFramesAreParallel(); + + /** Groups all physical frames by their position. This also works if the physical frames + * have slightly different positions, i.e. if they are not exactly the same and are only + * "close enough" to be considered the same. Right now, the maximum distance treated equal + * is if distance is smaller than slice thickness * 0.01 (i.e. 1% of slice thickness). + * Only performs the computation, if not done before. + * @return EC_Normal if successful, error otherwise + */ + OFCondition groupFramesByPosition(); + + /** Checks whether the given two frames overlap + * @param f1 Frame 1, provided by its physical frame number + * @param f2 Frame 2, provided by its physical frame number + * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) + * @return EC_Normal if successful, error otherwise + */ + OFCondition checkFramesOverlap(const Uint32& f1, const Uint32& f2, OFBool& overlap); + + /** Checks whether the given two frames overlap by using comparing their pixel data + * by bitwise "and". This is very efficient, however, only works and is called (right now), + * if row*cols % 8 = 0, so we can easily extract frames as binary bitsets without unpacking them. + * TODO: Check whether this can be easily extended to other cases as well. + * @param f1 Frame 1, provided by its physical frame number + * @param f2 Frame 2, provided by its physical frame number + * @param f1_data Pixel data of frame 1 + * @param f2_data Pixel data of frame 2 + * @param rows Number of rows of the frame(s) + * @param cols Number of columns of the frame(s) + * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) + * @return EC_Normal if successful, error otherwise + */ + OFCondition checkFramesOverlapBinary(const Uint32& f1, + const Uint32& f2, + const DcmIODTypes::Frame* f1_data, + const DcmIODTypes::Frame* f2_data, + const Uint16& rows, + const Uint16 cols, + OFBool& overlap); + + /** Checks whether the given two frames overlap by using comparing their pixel data + * after unpacking, i.e. expanding every bit to a byte, and then comparing whether the two + * related bytes of each frame are both non-zero. This is less efficient than checkFramesOverlapBinary(), + * @param f1 Frame 1, provided by its physical frame number + * @param f2 Frame 2, provided by its physical frame number + * @param f1_data Pixel data of frame 1 + * @param f2_data Pixel data of frame 2 + * @param rows Number of rows of the frame(s) + * @param cols Number of columns of the frame(s) + * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) + * @return EC_Normal if successful, error otherwise + */ + OFCondition checkFramesOverlapUnpacked(const Uint32& f1, + const Uint32& f2, + const DcmIODTypes::Frame* f1_data, + const DcmIODTypes::Frame* f2_data, + const Uint16& rows, + const Uint16 cols, + OFBool& overlap); + + /** Return the most relevant (changing) coordinate, computed by multiplying + * x and y vectors of the image orientation and selecting the coordinate + * with the largest absolute value. + * @param imageOrientation Image orientation patient (3 coordinates for x vector, + * 3 coordinates for y vector ) + * @return 0 if x, 1 if y, 2 if z, 3 if not determinable + */ + static Uint8 identifyChangingCoordinate(const OFVector& imageOrientation); + +private: + /// Image Orientation Patient + OFVector m_imageOrientation; + + /// Phyiscal frames with their respective positions (IPP) + FramePositions m_framePositions; + + /// Outer vector with one entry per segment. Index is the DICOM segment + /// number where segment 1 goes to index 0, segment 2 to index 1, and so on. + /// Inner vector contains the physical frame numbers that represent the + /// segment. + FramesForSegment m_framesForSegment; + + /// Logical frames, ie. physical frames with the same position are + /// grouped together to a logical frame. For every logical frame, we + /// store the related physical frame numbers. The logical frame number + /// is implicitly given by the index in the vector. + DistinctFramePositions m_logicalFramePositions; + + /// Stores for each logical frame a collection of (paired) segment and + /// physical frame number, that exists at that position. + SegmentsByPosition m_segmentsByPosition; + + /// Matrix that stores for each segment pair whether they overlap or not. + /// I.e. Matrix has size N x N, where N is the number of segments. + /// The diagonal is always 0 (no overlap), i.e. a segment never overlaps with itself. + /// If there is an overlap, the value is 1. + /// If the field is not initialized, the value is -1. + OverlapMatrix m_segmentOverlapMatrix; + + //// Groups of segments that do not overlap with each other + SegmentGroups m_nonOverlappingSegments; + + /// Reference to segmentation object to work with + /// Must be freed outside this class. + DcmSegmentation* m_seg; +}; + +} // namespace dcmqi + +#endif // DCMQI_OVERLAPUTIL_H \ No newline at end of file diff --git a/include/dcmqi/itkTestMain.h b/include/dcmqi/itkTestMain.h index 9940967e..fc6f7cf7 100644 --- a/include/dcmqi/itkTestMain.h +++ b/include/dcmqi/itkTestMain.h @@ -41,8 +41,6 @@ #include "itkImage.h" #include "itkImageFileReader.h" #include "itkImageFileWriter.h" -#include "itkImageRegionConstIterator.h" -#include "itkSubtractImageFilter.h" #include "itkRescaleIntensityImageFilter.h" #include "itkExtractImageFilter.h" #include "itkFloatingPointExceptions.h" @@ -257,6 +255,9 @@ int RegressionTestImage(const char *testImageFilename, ::itk::SizeValueType numberOfPixelsTolerance, unsigned int radiusTolerance) { + // Print filenames + std::cout << "Test image: " << testImageFilename << std::endl; + std::cout << "Baseline image: " << baselineImageFilename << std::endl; // Use the factory mechanism to read the test and baseline files and convert // them to double typedef itk::Image ImageType; @@ -314,12 +315,25 @@ int RegressionTestImage(const char *testImageFilename, diff->SetDifferenceThreshold(intensityTolerance); diff->SetToleranceRadius(radiusTolerance); diff->UpdateLargestPossibleRegion(); + std::cout << "Number of differences: " << diff->GetNumberOfPixelsWithDifferences() << std::endl; + // Write difference image to file "testoutput" if there are discrepencies + if( diff->GetNumberOfPixelsWithDifferences() > 0 ) + { + std::cout << "Difference image was saved as " << "testoutput" << std::endl; + typedef itk::ImageFileWriter WriterType; + WriterType::Pointer writer = WriterType::New(); + writer->SetFileName("testoutput.nrrd"); + writer->SetInput(diff->GetOutput()); + writer->UpdateLargestPossibleRegion(); + } const itk::SizeValueType status = diff->GetNumberOfPixelsWithDifferences(); // if there are discrepencies, create an diff image if( ( status > numberOfPixelsTolerance ) && reportErrors ) { + std::cout << "RegressionTestimage: Number of pixels with differences = " + << status << ", threshold = " << numberOfPixelsTolerance << std::endl; typedef itk::RescaleIntensityImageFilter RescaleType; typedef itk::ExtractImageFilter ExtractType; typedef itk::ImageFileWriter WriterType; @@ -465,9 +479,11 @@ int RegressionTestImage(const char *testImageFilename, std::cout << testName.str(); std::cout << "" << std::endl; } + std::cout << "RegressionTestimage: Number of pixels with differences = " + << status << ", threshold = " << numberOfPixelsTolerance << std::endl; return ( status > numberOfPixelsTolerance ) ? 1 : 0; -} +} // // Generate all of the possible baselines // The possible baselines are generated fromn the baselineFilename using the diff --git a/libsrc/CMakeLists.txt b/libsrc/CMakeLists.txt index 2faa2b74..9307b538 100644 --- a/libsrc/CMakeLists.txt +++ b/libsrc/CMakeLists.txt @@ -10,28 +10,32 @@ set(HDRS ${INCLUDE_DIR}/QIICRConstants.h ${INCLUDE_DIR}/QIICRUIDs.h ${INCLUDE_DIR}/ConverterBase.h + ${INCLUDE_DIR}/Dicom2ItkConverter.h ${INCLUDE_DIR}/Exceptions.h ${INCLUDE_DIR}/framesorter.h - ${INCLUDE_DIR}/ImageSEGConverter.h + ${INCLUDE_DIR}/Itk2DicomConverter.h ${INCLUDE_DIR}/ParaMapConverter.h ${INCLUDE_DIR}/Helper.h ${INCLUDE_DIR}/ColorUtilities.h ${INCLUDE_DIR}/JSONMetaInformationHandlerBase.h ${INCLUDE_DIR}/JSONParametricMapMetaInformationHandler.h ${INCLUDE_DIR}/JSONSegmentationMetaInformationHandler.h + ${INCLUDE_DIR}/OverlapUtil.h ${INCLUDE_DIR}/SegmentAttributes.h ${INCLUDE_DIR}/TID1500Reader.h ) set(SRCS ConverterBase.cpp - ImageSEGConverter.cpp + Dicom2ItkConverter.cpp ParaMapConverter.cpp Helper.cpp ColorUtilities.cpp + Itk2DicomConverter.cpp JSONMetaInformationHandlerBase.cpp JSONParametricMapMetaInformationHandler.cpp JSONSegmentationMetaInformationHandler.cpp + OverlapUtil.cpp SegmentAttributes.cpp TID1500Reader.cpp ) @@ -77,8 +81,8 @@ set(${lib_name}_INCLUDE_DIRS ${_dcmtk_includes} ${ITK_INCLUDE_DIRS} ${JsonCpp_INCLUDE_DIR} - $ - $ + $ + $ ) if(ITK_INSTALL_PREFIX) if(EXISTS ${ITK_INSTALL_PREFIX}/include/vxl/core AND EXISTS ${ITK_INSTALL_PREFIX}/include/vxl/vcl) diff --git a/libsrc/Dicom2ItkConverter.cpp b/libsrc/Dicom2ItkConverter.cpp new file mode 100644 index 00000000..16fa8da0 --- /dev/null +++ b/libsrc/Dicom2ItkConverter.cpp @@ -0,0 +1,542 @@ + +// DCMQI includes +#include "dcmqi/Dicom2ItkConverter.h" +#include "dcmqi/ColorUtilities.h" +#include "dcmqi/OverlapUtil.h" + +// DCMTK includes +#include +#include +#include +#include +#include +#include + +namespace dcmqi +{ + +Dicom2ItkConverter::Dicom2ItkConverter() + : m_segDoc() + , m_direction() + , m_computedSliceSpacing() + , m_computedVolumeExtent() + , m_sliceDirection(3) + , m_imageOrigin() + , m_imageSpacing() + , m_imageSize() + , m_imageRegion() + , m_metaInfo() + , m_groupIterator() + , m_overlapUtil() {}; + +// ------------------------------------------------------------------------------------- + +OFCondition +Dicom2ItkConverter::dcmSegmentation2itkimage(DcmDataset* segDataset, std::string& metaInfo, const bool mergeSegments) +{ + DcmSegmentation* segdoc = NULL; + + // Make sure RLE-compressed images can be decompressed + DcmRLEDecoderRegistration::registerCodecs(); + + // Load the DICOM segmentation dataset into DcmSegmentation member + OFCondition cond = DcmSegmentation::loadDataset(*segDataset, segdoc); + if (!segdoc) + { + cerr << "ERROR: Failed to load segmentation dataset! " << cond.text() << endl; + throw -1; + } + m_segDoc.reset(segdoc); + + // Populate DICOM series information into accompanying JSON metainfo member + populateMetaInformationFromDICOM(segDataset); + + OFCondition result = dcmSegmentation2itkimage(mergeSegments); + if (result.good()) + { + metaInfo = m_metaInfo.getJSONOutputAsString(); + } + return result; +} + +// ------------------------------------------------------------------------------------- + +OFCondition Dicom2ItkConverter::dcmSegmentation2itkimage(const bool mergeSegments) +{ + // Setup logging (not used so far?) + OFLogger dcemfinfLogger = OFLog::getLogger("qiicr.apps"); + dcemfinfLogger.setLogLevel(dcmtk::log4cplus::OFF_LOG_LEVEL); + + // Extract directions, origin, spacing and image region + OFCondition result = extractBasicSegmentationInfo(); + + // Find groups of segments that can go into the same ITK image (i.e. that are non-overlapping) + m_overlapUtil.setSegmentationObject(m_segDoc.get()); + result = getNonOverlappingSegmentGroups(mergeSegments, m_segmentGroups); + if (result.bad()) + { + return result; + } + + // Create JSON meta info for all segments first since this is returned + // immediately from this call, while the ITK result images are + // made accessible through result iterators only. + result = createMetaInfo(); + if (result.bad()) + { + return result; + } + + return result; +} + +// ------------------------------------------------------------------------------------- + +itk::SmartPointer Dicom2ItkConverter::begin() +{ + // Set result iterator to first group, i.e. make sure that the first call to nextResult() + // will return the ITK image for the first group. + m_groupIterator = m_segmentGroups.begin(); + return nextResult(); +} + +// ------------------------------------------------------------------------------------- + +itk::SmartPointer Dicom2ItkConverter::next() +{ + return nextResult(); +} + +// ------------------------------------------------------------------------------------- + +itk::SmartPointer Dicom2ItkConverter::nextResult() +{ + OFCondition result; + ShortImageType::Pointer itkImage = nullptr; + if (m_groupIterator != m_segmentGroups.end()) + { + // Create target ITK image for this group + itkImage = allocateITKImage(); + + // Loop over segments belonging to this segment group + auto segNum = m_groupIterator->begin(); + while (segNum != m_groupIterator->end()) + { + // Iterate over frames for this segment, and copy the data into the ITK image. + // Afterwards, the ITK image will have the complete data belonging to that segment + OverlapUtil::FramesForSegment::value_type framesForSegment; + m_overlapUtil.getFramesForSegment(*segNum, framesForSegment); + for (size_t frameIndex = 0; frameIndex < framesForSegment.size(); frameIndex++) + { + // Copy the data from the frame into the ITK image + ShortImageType::PointType frameOriginPoint; + ShortImageType::IndexType frameOriginIndex; + result = getITKImageOrigin(framesForSegment[frameIndex], frameOriginPoint); + if (result.bad()) + { + cerr << "ERROR: Failed to get origin for frame " << framesForSegment[frameIndex] << " of segment " + << *segNum << endl; + m_groupIterator = m_segmentGroups.end(); + return nullptr; + } + if (!itkImage->TransformPhysicalPointToIndex(frameOriginPoint, frameOriginIndex)) + { + cerr << "ERROR: Frame " << framesForSegment[frameIndex] << " origin " << frameOriginPoint + << " is outside image geometry!" << frameOriginIndex << endl; + cerr << "Image size: " << itkImage->GetBufferedRegion().GetSize() << endl; + m_groupIterator = m_segmentGroups.end(); + return nullptr; + } + // Handling differs depending on whether the segmentation is binary or fractional + // (we have to unpack binary frames before copying them into the ITK image) + const DcmIODTypes::Frame* rawFrame = m_segDoc->getFrame(framesForSegment[frameIndex]); + const DcmIODTypes::Frame* unpackedFrame = NULL; + if (m_segDoc->getSegmentationType() == DcmSegTypes::ST_BINARY) + { + unpackedFrame = DcmSegUtils::unpackBinaryFrame(rawFrame, + m_imageSize[1], // Rows + m_imageSize[0]); // Cols + } + else + { + // fractional segmentation frames can be used as is + unpackedFrame = rawFrame; + } + unsigned slice = frameOriginIndex[2]; + + { + /* WIP */ + // use ConstIterator, as in this example: + // https://github.com/InsightSoftwareConsortium/ITK/blob/633f84548311600845d54ab2463d3412194690a8/Examples/Iterators/ImageRegionIterator.cxx#L201-L212 + // CharImageType::RegionType frameRegion; + // CharImageType::RegionType::IndexType regionStart; + // CharImageType::RegionType::SizeType regionSize; + + // regionStart[0] = 0; + // regionStart[1] = 0; + // regionStart[2] = slice; + + // regionSize[0] = imageSize[0]; + // regionSize[1] = imageSize[1]; + // regionSize[2] = 1; + + // frameRegion.SetSize(regionSize); + // frameRegion.SetIndex(regionStart); + + // const bool importImageFilterWillOwnTheBuffer = false; + + // // itkImportImageFilter example: + // // + // https://itk.org/Doxygen/html/SphinxExamples_2src_2Core_2Common_2ConvertArrayToImage_2Code_8cxx-example.html#_a1 + // const unsigned int numberOfPixels = imageSize[0]*imageSize[1]; + // using ImportFilterType = itk::ImportImageFilter; + // ImportFilterType::Pointer importFilter = ImportFilterType::New(); + // importFilter->SetImportPointer((unsigned char*)unpackedFrame->pixData, numberOfPixels, + // importImageFilterWillOwnTheBuffer); importFilter->SetRegion(frameRegion); + // importFilter->Update(); + + // using ConstIteratorType = itk::ImageRegionConstIterator; + // using IteratorType = itk::ImageRegionIterator; + // ConstIteratorType inputIt(importFilter->GetOutput(), frameRegion); + // IteratorType outputIt(segment2image[targetSegmentsImageIndex], frameRegion); + + // inputIt.GoToBegin(); + // outputIt.GoToBegin(); + + // while(!inputIt.IsAtEnd()){ + // if(inputIt.Get() != 0) + // outputIt.Set(segmentIdLabel); + // ++outputIt; + // ++inputIt; + // } + /* WIP */ + } + + for (unsigned row = 0; row < m_imageSize[1]; row++) + { + for (unsigned col = 0; col < m_imageSize[0]; col++) + { + ShortImageType::PixelType pixel; + unsigned bitCnt = row * m_imageSize[0] + col; + pixel = unpackedFrame->pixData[bitCnt]; + ShortImageType::IndexType index; + if (pixel != 0) + { + index[0] = col; + index[1] = row; + index[2] = slice; + // Use the segment number as the pixel value in case of binary segmentations. + // Otherwise, just copy the pixel value (fractional value) from the frame. + if (m_segDoc->getSegmentationType() == DcmSegTypes::ST_BINARY) + { + itkImage->SetPixel(index, *segNum); + } + else + { + itkImage->SetPixel(index, pixel); + } + } + } + } + if (m_segDoc->getSegmentationType() == DcmSegTypes::ST_BINARY) + { + delete unpackedFrame; + unpackedFrame = NULL; + } + } + segNum++; + } + m_groupIterator++; + } + return itk::SmartPointer(itkImage); +} + +// ------------------------------------------------------------------------------------- + +void Dicom2ItkConverter::populateMetaInformationFromDICOM(DcmDataset* segDataset) +{ + OFString creatorName, sessionID, timePointID, seriesDescription, seriesNumber, instanceNumber, bodyPartExamined, + coordinatingCenter; + + segDataset->findAndGetOFString(DCM_InstanceNumber, instanceNumber); + m_segDoc->getContentIdentification().getContentCreatorName(creatorName); + + segDataset->findAndGetOFString(DCM_ClinicalTrialTimePointID, timePointID); + segDataset->findAndGetOFString(DCM_ClinicalTrialSeriesID, sessionID); + segDataset->findAndGetOFString(DCM_ClinicalTrialCoordinatingCenterName, coordinatingCenter); + + m_segDoc->getSeries().getBodyPartExamined(bodyPartExamined); + m_segDoc->getSeries().getSeriesNumber(seriesNumber); + m_segDoc->getSeries().getSeriesDescription(seriesDescription); + + m_metaInfo.setClinicalTrialCoordinatingCenterName(coordinatingCenter.c_str()); + m_metaInfo.setContentCreatorName(creatorName.c_str()); + m_metaInfo.setClinicalTrialSeriesID(sessionID.c_str()); + m_metaInfo.setSeriesNumber(seriesNumber.c_str()); + m_metaInfo.setClinicalTrialTimePointID(timePointID.c_str()); + m_metaInfo.setSeriesDescription(seriesDescription.c_str()); + m_metaInfo.setInstanceNumber(instanceNumber.c_str()); + m_metaInfo.setBodyPartExamined(bodyPartExamined.c_str()); +} + +// ------------------------------------------------------------------------------------- + +OFCondition Dicom2ItkConverter::extractBasicSegmentationInfo() +{ + // TODO: Better error handling + OFCondition result; + + // Directions + FGInterface& fgInterface = m_segDoc->getFunctionalGroups(); + if (getImageDirections(fgInterface, m_direction)) + { + cerr << "ERROR: Failed to get image directions!" << endl; + throw -1; + } + + // Origin + m_sliceDirection[0] = m_direction[0][2]; + m_sliceDirection[1] = m_direction[1][2]; + m_sliceDirection[2] = m_direction[2][2]; + if (computeVolumeExtent( + fgInterface, m_sliceDirection, m_imageOrigin, m_computedSliceSpacing, m_computedVolumeExtent)) + { + cerr << "ERROR: Failed to compute origin and/or slice spacing!" << endl; + throw -1; + } + + // Spacing + m_imageSpacing.Fill(0); + if (getDeclaredImageSpacing(fgInterface, m_imageSpacing)) + { + cerr << "ERROR: Failed to get image spacing from DICOM!" << endl; + throw -1; + } + + if (!m_imageSpacing[2]) + { + cerr << "ERROR: No sufficient information to derive slice spacing! Unable to interpret the data." << endl; + throw -1; + } + + // Region, defined by DICOM rows and columns + { + DcmIODImage>* ip + = static_cast>*>(m_segDoc.get()); + Uint16 value; + if (ip->getImagePixel().getRows(value).good()) + { + m_imageSize[1] = value; + } + if (ip->getImagePixel().getColumns(value).good()) + { + m_imageSize[0] = value; + } + } + // Number of slices should be computed, since segmentation may have empty frames + m_imageSize[2] = round(m_computedVolumeExtent / m_imageSpacing[2]) + 1; + + // Initialize the image template. This will only serve as a template for the individual + // frames, which will be created via ImageDuplicator later on. + m_imageRegion.SetSize(m_imageSize); + return result; +} + +// ------------------------------------------------------------------------------------- + +OFCondition Dicom2ItkConverter::getNonOverlappingSegmentGroups(const bool mergeSegments, + OverlapUtil::SegmentGroups& segmentGroups) +{ + OFCondition result; + if (mergeSegments) + { + result = m_overlapUtil.getNonOverlappingSegments(segmentGroups); + if (result.bad()) + { + cout << "WARNING: Failed to compute non-overlapping segments (Error: " << result.text() << "), " + << "falling back to one group per segment instead." << endl; + } + cout << "Identified " << segmentGroups.size() << " groups of non-overlapping segments" << endl; + } + // Otherwise, use single group containing all segments (which might overlap) + if (!mergeSegments || result.bad()) + { + size_t numSegs = m_segDoc->getNumberOfSegments(); + for (size_t i = 1; i <= numSegs; ++i) + { + OFVector segs; + segs.push_back(i); + segmentGroups.push_back(segs); + } + cout << "Will not merge segments: Splitting segments into " << segmentGroups.size() << " groups" << endl; + } + return result; + ; +} + +// ------------------------------------------------------------------------------------- + +ShortImageType::Pointer Dicom2ItkConverter::allocateITKImage() +{ + // Initialize the image + ShortImageType::Pointer segImage = ShortImageType::New(); + segImage->SetRegions(m_imageRegion); + segImage->SetOrigin(m_imageOrigin); + segImage->SetSpacing(m_imageSpacing); + segImage->SetDirection(m_direction); + segImage->Allocate(); + segImage->FillBuffer(0); + return segImage; +} + +ShortImageType::Pointer Dicom2ItkConverter::allocateITKImageDuplicate(ShortImageType::Pointer imageTemplate) +{ + typedef itk::ImageDuplicator DuplicatorType; + DuplicatorType::Pointer dup = DuplicatorType::New(); + dup->SetInputImage(imageTemplate); + dup->Update(); + ShortImageType::Pointer newSegmentImage = dup->GetOutput(); + newSegmentImage->FillBuffer(0); + return newSegmentImage; +} + +// ------------------------------------------------------------------------------------- + +OFCondition Dicom2ItkConverter::getITKImageOrigin(const Uint32 frameNo, ShortImageType::PointType& origin) +{ + FGInterface& fgInterface = m_segDoc->getFunctionalGroups(); + + FGPlanePosPatient* planposfg + = OFstatic_cast(FGPlanePosPatient*, fgInterface.get(frameNo, DcmFGTypes::EFG_PLANEPOSPATIENT)); + assert(planposfg); + + for (int j = 0; j < 3; j++) + { + OFString planposStr; + if (planposfg->getImagePositionPatient(planposStr, j).good()) + { + origin[j] = atof(planposStr.c_str()); + } + } + return EC_Normal; +} + +// ------------------------------------------------------------------------------------- + +OFCondition Dicom2ItkConverter::createMetaInfo() +{ + OFCondition result; + OverlapUtil::SegmentGroups::const_iterator group = m_segmentGroups.begin(); + size_t groupNumber = 1; + while (group != m_segmentGroups.end()) + { + auto segNum = group->begin(); + while (result.good() && (segNum != group->end())) + { + result = addSegmentMetadata(groupNumber, *segNum); + segNum++; + } + group++; + groupNumber++; + } + if (result.bad()) + { + cerr << "ERROR: Failed to create segment metadata: " << result.text() << endl; + } + return result; +} + +// ------------------------------------------------------------------------------------- + +OFCondition Dicom2ItkConverter::addSegmentMetadata(const size_t segmentGroup, const Uint16 segmentNumber) +{ + SegmentAttributes* segmentAttributes = m_metaInfo.createOrGetSegment(segmentGroup, segmentNumber); + // NOTE: Segment numbers in DICOM start with 1 + DcmSegment* segment = m_segDoc->getSegment(segmentNumber); + if (segment == NULL) + { + cerr << "Failed to get segment for segment ID " << segmentNumber << endl; + return EC_IllegalParameter; + } + + // get CIELab color for the segment + Uint16 ciedcm[3]; + int rgb[3]; + if (segment->getRecommendedDisplayCIELabValue(ciedcm[0], ciedcm[1], ciedcm[2]).bad()) + { + // NOTE: if the call above fails, it overwrites the values anyway, + // not sure if this is a dcmtk bug or not + ciedcm[0] = 43803; + ciedcm[1] = 26565; + ciedcm[2] = 37722; + cerr << "Failed to get CIELab values - initializing to default " << ciedcm[0] << "," << ciedcm[1] << "," + << ciedcm[2] << endl; + } + + ColorUtilities::getSRGBFromIntegerScaledCIELabPCS(rgb[0], rgb[1], rgb[2], ciedcm[0], ciedcm[1], ciedcm[2]); + + if (segmentAttributes) + { + segmentAttributes->setLabelID(segmentNumber); + DcmSegTypes::E_SegmentAlgoType algorithmType = segment->getSegmentAlgorithmType(); + string readableAlgorithmType = DcmSegTypes::algoType2OFString(algorithmType).c_str(); + segmentAttributes->setSegmentAlgorithmType(readableAlgorithmType); + + if (algorithmType == DcmSegTypes::SAT_UNKNOWN) + { + cerr << "ERROR: AlgorithmType is not valid with value " << readableAlgorithmType << endl; + throw -1; + } + if (algorithmType != DcmSegTypes::SAT_MANUAL) + { + OFString segmentAlgorithmName; + segment->getSegmentAlgorithmName(segmentAlgorithmName); + if (segmentAlgorithmName.length() > 0) + segmentAttributes->setSegmentAlgorithmName(segmentAlgorithmName.c_str()); + } + + OFString segmentDescription, segmentLabel, trackingIdentifier, trackingUniqueIdentifier; + + segment->getSegmentDescription(segmentDescription); + segmentAttributes->setSegmentDescription(segmentDescription.c_str()); + + segment->getSegmentLabel(segmentLabel); + segmentAttributes->setSegmentLabel(segmentLabel.c_str()); + + segment->getTrackingID(trackingIdentifier); + segment->getTrackingUID(trackingUniqueIdentifier); + + if (trackingIdentifier.length() > 0) + { + segmentAttributes->setTrackingIdentifier(trackingIdentifier.c_str()); + } + if (trackingUniqueIdentifier.length() > 0) + { + segmentAttributes->setTrackingUniqueIdentifier(trackingUniqueIdentifier.c_str()); + } + + segmentAttributes->setRecommendedDisplayRGBValue(rgb[0], rgb[1], rgb[2]); + segmentAttributes->setSegmentedPropertyCategoryCodeSequence(segment->getSegmentedPropertyCategoryCode()); + segmentAttributes->setSegmentedPropertyTypeCodeSequence(segment->getSegmentedPropertyTypeCode()); + + if (segment->getSegmentedPropertyTypeModifierCode().size() > 0) + { + segmentAttributes->setSegmentedPropertyTypeModifierCodeSequence( + segment->getSegmentedPropertyTypeModifierCode()[0]); + } + + GeneralAnatomyMacro& anatomyMacro = segment->getGeneralAnatomyCode(); + CodeSequenceMacro& anatomicRegionSequence = anatomyMacro.getAnatomicRegion(); + if (anatomicRegionSequence.check(true).good()) + { + segmentAttributes->setAnatomicRegionSequence(anatomyMacro.getAnatomicRegion()); + } + if (anatomyMacro.getAnatomicRegionModifier().size() > 0) + { + segmentAttributes->setAnatomicRegionModifierSequence(*(anatomyMacro.getAnatomicRegionModifier()[0])); + } + } + return EC_Normal; +} + +} diff --git a/libsrc/ImageSEGConverter.cpp b/libsrc/Itk2DicomConverter.cpp similarity index 61% rename from libsrc/ImageSEGConverter.cpp rename to libsrc/Itk2DicomConverter.cpp index dc6734e1..6b011050 100644 --- a/libsrc/ImageSEGConverter.cpp +++ b/libsrc/Itk2DicomConverter.cpp @@ -1,18 +1,24 @@ // DCMQI includes -#include "dcmqi/ImageSEGConverter.h" +#include "dcmqi/Itk2DicomConverter.h" #include "dcmqi/ColorUtilities.h" +#include "dcmqi/JSONSegmentationMetaInformationHandler.h" - -//DCMTK includes -#include +// DCMTK includes #include -#include + namespace dcmqi { - DcmDataset* ImageSEGConverter::itkimage2dcmSegmentation(vector dcmDatasets, + + Itk2DicomConverter::Itk2DicomConverter() + { + }; + + // ------------------------------------------------------------------------------------- + + DcmDataset* Itk2DicomConverter::itkimage2dcmSegmentation(vector dcmDatasets, vector segmentations, const string &metaData, bool skipEmptySlices) { @@ -397,7 +403,6 @@ namespace dcmqi { // clean up for the next frame fgder->clearData(); } - } } } @@ -421,10 +426,10 @@ namespace dcmqi { CHECK_COND(segdoc->getFrameOfReference().setFrameOfReferenceUID(frameOfRefUIDchar)); } + std::cout << "Writing DICOM segmentation dataset " << std::endl; // Don't check functional groups since its very time consuming and we trust // ourselves to put together valid datasets segdoc->setCheckFGOnWrite(OFFalse); - OFCondition writeResult = segdoc->writeDataset(segdocDataset); if(writeResult.bad()){ cerr << "FATAL ERROR: Writing of the SEG dataset failed!"; @@ -436,6 +441,7 @@ namespace dcmqi { } // Set reader/session/timepoint information + std::cout << "Patching in extra meta information into DICOM dataset" << std::endl; CHECK_COND(segdocDataset.putAndInsertString(DCM_SeriesDescription, metaInfo.getSeriesDescription().c_str())); CHECK_COND(segdocDataset.putAndInsertString(DCM_ContentCreatorName, metaInfo.getContentCreatorName().c_str())); CHECK_COND(segdocDataset.putAndInsertString(DCM_ClinicalTrialSeriesID, metaInfo.getClinicalTrialSeriesID().c_str())); @@ -483,328 +489,4 @@ namespace dcmqi { return new DcmDataset(segdocDataset); } - - pair , string> ImageSEGConverter::dcmSegmentation2itkimage(DcmDataset *segDataset) { - DcmSegmentation *segdoc = NULL; - - DcmRLEDecoderRegistration::registerCodecs(); - - OFCondition cond = DcmSegmentation::loadDataset(*segDataset, segdoc); - if(!segdoc){ - cerr << "ERROR: Failed to load segmentation dataset! " << cond.text() << endl; - throw -1; - } - OFunique_ptr segdocguard(segdoc); - - JSONSegmentationMetaInformationHandler metaInfo; - populateMetaInformationFromDICOM(segDataset, segdoc, metaInfo); - - const std::map segment2image = - dcmSegmentation2itkimage(segdoc, &metaInfo); - return pair, string>( - segment2image, metaInfo.getJSONOutputAsString()); - } - - std::map - ImageSEGConverter::dcmSegmentation2itkimage( - DcmSegmentation *segdoc, - JSONSegmentationMetaInformationHandler *metaInfo) { - - OFLogger dcemfinfLogger = OFLog::getLogger("qiicr.apps"); - dcemfinfLogger.setLogLevel(dcmtk::log4cplus::OFF_LOG_LEVEL); - - // Directions - FGInterface &fgInterface = segdoc->getFunctionalGroups(); - ShortImageType::DirectionType direction; - if(getImageDirections(fgInterface, direction)){ - cerr << "ERROR: Failed to get image directions!" << endl; - throw -1; - } - - // Spacing and origin - double computedSliceSpacing, computedVolumeExtent; - vnl_vector sliceDirection(3); - sliceDirection[0] = direction[0][2]; - sliceDirection[1] = direction[1][2]; - sliceDirection[2] = direction[2][2]; - - ShortImageType::PointType imageOrigin; - if(computeVolumeExtent(fgInterface, sliceDirection, imageOrigin, computedSliceSpacing, computedVolumeExtent)){ - cerr << "ERROR: Failed to compute origin and/or slice spacing!" << endl; - throw -1; - } - - ShortImageType::SpacingType imageSpacing; - imageSpacing.Fill(0); - if(getDeclaredImageSpacing(fgInterface, imageSpacing)){ - cerr << "ERROR: Failed to get image spacing from DICOM!" << endl; - throw -1; - } - - if(!imageSpacing[2]){ - cerr << "ERROR: No sufficient information to derive slice spacing! Unable to interpret the data." << endl; - throw -1; - } - - // Region size - ShortImageType::SizeType imageSize; - { - DcmIODImage > *ip = - static_cast > *>(segdoc); - Uint16 value; - if (ip->getImagePixel().getRows(value).good()) { - imageSize[1] = value; - } - if (ip->getImagePixel().getColumns(value).good()) { - imageSize[0] = value; - } - } - // number of slices should be computed, since segmentation may have empty frames - imageSize[2] = round(computedVolumeExtent/imageSpacing[2])+1; - - // Initialize the image - ShortImageType::RegionType imageRegion; - imageRegion.SetSize(imageSize); - ShortImageType::Pointer segImage = ShortImageType::New(); - segImage->SetRegions(imageRegion); - segImage->SetOrigin(imageOrigin); - segImage->SetSpacing(imageSpacing); - segImage->SetDirection(direction); - segImage->Allocate(); - segImage->FillBuffer(0); - - // ITK images corresponding to the individual segments - map segment2image; - - // Iterate over frames, find the matching slice for each of the frames based on - // ImagePositionPatient, set non-zero pixels to the segment number. Notify - // about pixels that are initialized more than once. - - const DcmIODTypes::Frame *unpackedFrame = NULL; - - for(size_t frameId=0;frameIdgetFrame(frameId); - bool isPerFrame; - - FGPlanePosPatient *planposfg = - OFstatic_cast(FGPlanePosPatient*,fgInterface.get(frameId, DcmFGTypes::EFG_PLANEPOSPATIENT, isPerFrame)); - assert(planposfg); - -#ifndef NDEBUG - FGFrameContent *fracon = - OFstatic_cast(FGFrameContent*,fgInterface.get(frameId, DcmFGTypes::EFG_FRAMECONTENT, isPerFrame)); - assert(fracon); -#endif - - FGSegmentation *fgseg = - OFstatic_cast(FGSegmentation*,fgInterface.get(frameId, DcmFGTypes::EFG_SEGMENTATION, isPerFrame)); - assert(fgseg); - - Uint16 segmentId = -1; - if(fgseg->getReferencedSegmentNumber(segmentId).bad()){ - cerr << "ERROR: Failed to get ReferencedSegmentNumber!"; - throw -1; - } - - // WARNING: this is needed only for David's example, which numbers - // (incorrectly!) segments starting from 0, should start from 1 - if(segmentId == 0){ - cerr << "ERROR: ReferencedSegmentNumber value of 0 was encountered. Segment numbers and references to segment numbers should both start from 1!" << endl; - throw -1; - } - - if(segment2image.find(segmentId) == segment2image.end()){ - typedef itk::ImageDuplicator DuplicatorType; - DuplicatorType::Pointer dup = DuplicatorType::New(); - dup->SetInputImage(segImage); - dup->Update(); - ShortImageType::Pointer newSegmentImage = dup->GetOutput(); - newSegmentImage->FillBuffer(0); - segment2image[segmentId] = newSegmentImage; - } - - // populate meta information needed for Slicer ScalarVolumeNode initialization - // (for example) - { - // NOTE: according to the standard, segment numbering should start from 1, - // not clear if this is intentional behavior or a bug in DCMTK expecting - // it to start from 0 - - DcmSegment* segment = segdoc->getSegment(segmentId); - if(segment == NULL){ - cerr << "Failed to get segment for segment ID " << segmentId << endl; - continue; - } - - // get CIELab color for the segment - Uint16 ciedcm[3]; - int rgb[3]; - if(segment->getRecommendedDisplayCIELabValue( - ciedcm[0], ciedcm[1], ciedcm[2] - ).bad()) { - // NOTE: if the call above fails, it overwrites the values anyway, - // not sure if this is a dcmtk bug or not - ciedcm[0] = 43803; - ciedcm[1] = 26565; - ciedcm[2] = 37722; - cerr << "Failed to get CIELab values - initializing to default " << - ciedcm[0] << "," << ciedcm[1] << "," << ciedcm[2] << endl; - } - - /* Debug CIELab 2 RGB conversion - double lab[3]; - - cout << "DICOM Lab " << ciedcm[0] << ", " << ciedcm[1] << ", " << ciedcm[2] << endl; - - IODCIELabUtil::dicomlab2Lab(lab[0], lab[1], lab[2], ciedcm[0], ciedcm[1], ciedcm[2]); - - cout << "Lab " << ciedcm[0] << ", " << ciedcm[1] << ", " << ciedcm[2] << endl; - - IODCIELabUtil::lab2Rgb(rgb[0], rgb[1], rgb[2], lab[0], lab[1], lab[2]); - - cout << "RGB " << unsigned(rgb[0]*256) << ", " << unsigned(rgb[1]*256) << ", " << unsigned(rgb[2]*256) << endl; - */ - - ColorUtilities::getSRGBFromIntegerScaledCIELabPCS(rgb[0], rgb[1], rgb[2], ciedcm[0], ciedcm[1], ciedcm[2]); - //IODCIELabUtil::dicomLab2RGB(rgb[0], rgb[1], rgb[2], ciedcm[0], ciedcm[1], ciedcm[2]); - - //rgb[0] = unsigned(rgb[0]*256); - //rgb[1] = unsigned(rgb[1]*256); - //rgb[2] = unsigned(rgb[2]*256); - - SegmentAttributes* segmentAttributes = metaInfo?metaInfo->createAndGetNewSegment(segmentId):NULL; - - if (segmentAttributes) { - segmentAttributes->setLabelID(segmentId); - DcmSegTypes::E_SegmentAlgoType algorithmType = segment->getSegmentAlgorithmType(); - string readableAlgorithmType = DcmSegTypes::algoType2OFString(algorithmType).c_str(); - segmentAttributes->setSegmentAlgorithmType(readableAlgorithmType); - - if (algorithmType == DcmSegTypes::SAT_UNKNOWN) { - cerr << "ERROR: AlgorithmType is not valid with value " << readableAlgorithmType << endl; - throw -1; - } - if (algorithmType != DcmSegTypes::SAT_MANUAL) { - OFString segmentAlgorithmName; - segment->getSegmentAlgorithmName(segmentAlgorithmName); - if(segmentAlgorithmName.length() > 0) - segmentAttributes->setSegmentAlgorithmName(segmentAlgorithmName.c_str()); - } - - OFString segmentDescription, segmentLabel, trackingIdentifier, trackingUniqueIdentifier; - - segment->getSegmentDescription(segmentDescription); - segmentAttributes->setSegmentDescription(segmentDescription.c_str()); - - segment->getSegmentLabel(segmentLabel); - segmentAttributes->setSegmentLabel(segmentLabel.c_str()); - - segment->getTrackingID(trackingIdentifier); - segment->getTrackingUID(trackingUniqueIdentifier); - - if (trackingIdentifier.length() > 0) { - segmentAttributes->setTrackingIdentifier(trackingIdentifier.c_str()); - } - if (trackingUniqueIdentifier.length() > 0) { - segmentAttributes->setTrackingUniqueIdentifier(trackingUniqueIdentifier.c_str()); - } - - segmentAttributes->setRecommendedDisplayRGBValue(rgb[0], rgb[1], rgb[2]); - segmentAttributes->setSegmentedPropertyCategoryCodeSequence(segment->getSegmentedPropertyCategoryCode()); - segmentAttributes->setSegmentedPropertyTypeCodeSequence(segment->getSegmentedPropertyTypeCode()); - - if (segment->getSegmentedPropertyTypeModifierCode().size() > 0) { - segmentAttributes->setSegmentedPropertyTypeModifierCodeSequence( - segment->getSegmentedPropertyTypeModifierCode()[0]); - } - - GeneralAnatomyMacro &anatomyMacro = segment->getGeneralAnatomyCode(); - CodeSequenceMacro& anatomicRegionSequence = anatomyMacro.getAnatomicRegion(); - if (anatomicRegionSequence.check(true).good()) { - segmentAttributes->setAnatomicRegionSequence(anatomyMacro.getAnatomicRegion()); - } - if (anatomyMacro.getAnatomicRegionModifier().size() > 0) { - segmentAttributes->setAnatomicRegionModifierSequence(*(anatomyMacro.getAnatomicRegionModifier()[0])); - } - } - } - - // get string representation of the frame origin - ShortImageType::PointType frameOriginPoint; - ShortImageType::IndexType frameOriginIndex; - for(int j=0;j<3;j++){ - OFString planposStr; - if(planposfg->getImagePositionPatient(planposStr, j).good()){ - frameOriginPoint[j] = atof(planposStr.c_str()); - } - } - - if(!segment2image[segmentId]->TransformPhysicalPointToIndex(frameOriginPoint, frameOriginIndex)){ - cerr << "ERROR: Frame " << frameId << " origin " << frameOriginPoint << - " is outside image geometry!" << frameOriginIndex << endl; - cerr << "Image size: " << segment2image[segmentId]->GetBufferedRegion().GetSize() << endl; - throw -1; - } - - unsigned slice = frameOriginIndex[2]; - - bool deleteFrame(true); - if(segdoc->getSegmentationType() == DcmSegTypes::ST_BINARY) - unpackedFrame = DcmSegUtils::unpackBinaryFrame(frame, - imageSize[1], // Rows - imageSize[0]); // Cols - else { - unpackedFrame = frame; - deleteFrame = false; - } - - // initialize slice with the frame content - for(unsigned row=0;rowpixData[bitCnt]; - - if(pixel!=0){ - ShortImageType::IndexType index; - index[0] = col; - index[1] = row; - index[2] = slice; - segment2image[segmentId]->SetPixel(index, segmentId); - } - } - } - - if(deleteFrame) - delete unpackedFrame; - } - - return segment2image; - } - - void ImageSEGConverter::populateMetaInformationFromDICOM(DcmDataset *segDataset, DcmSegmentation *segdoc, - JSONSegmentationMetaInformationHandler &metaInfo) { - OFString creatorName, sessionID, timePointID, seriesDescription, seriesNumber, instanceNumber, bodyPartExamined, coordinatingCenter; - - segDataset->findAndGetOFString(DCM_InstanceNumber, instanceNumber); - segdoc->getContentIdentification().getContentCreatorName(creatorName); - - segDataset->findAndGetOFString(DCM_ClinicalTrialTimePointID, timePointID); - segDataset->findAndGetOFString(DCM_ClinicalTrialSeriesID, sessionID); - segDataset->findAndGetOFString(DCM_ClinicalTrialCoordinatingCenterName, coordinatingCenter); - - segdoc->getSeries().getBodyPartExamined(bodyPartExamined); - segdoc->getSeries().getSeriesNumber(seriesNumber); - segdoc->getSeries().getSeriesDescription(seriesDescription); - - metaInfo.setClinicalTrialCoordinatingCenterName(coordinatingCenter.c_str()); - metaInfo.setContentCreatorName(creatorName.c_str()); - metaInfo.setClinicalTrialSeriesID(sessionID.c_str()); - metaInfo.setSeriesNumber(seriesNumber.c_str()); - metaInfo.setClinicalTrialTimePointID(timePointID.c_str()); - metaInfo.setSeriesDescription(seriesDescription.c_str()); - metaInfo.setInstanceNumber(instanceNumber.c_str()); - metaInfo.setBodyPartExamined(bodyPartExamined.c_str()); - } - } diff --git a/libsrc/JSONParametricMapMetaInformationHandler.cpp b/libsrc/JSONParametricMapMetaInformationHandler.cpp index 59bee849..f380bea4 100644 --- a/libsrc/JSONParametricMapMetaInformationHandler.cpp +++ b/libsrc/JSONParametricMapMetaInformationHandler.cpp @@ -204,9 +204,11 @@ namespace dcmqi { data["SourceImageDiffusionBValues"].append(*it); } - Json::StyledWriter styledWriter; - - ss << styledWriter.write(data); + Json::StreamWriterBuilder styledBuilder; + styledBuilder["indentation"] = " "; + std::unique_ptr writer(styledBuilder.newStreamWriter()); + writer->write(data, &ss); + ss << std::endl; return ss.str(); } diff --git a/libsrc/JSONSegmentationMetaInformationHandler.cpp b/libsrc/JSONSegmentationMetaInformationHandler.cpp index a9ad3a0a..8780e7ad 100644 --- a/libsrc/JSONSegmentationMetaInformationHandler.cpp +++ b/libsrc/JSONSegmentationMetaInformationHandler.cpp @@ -86,7 +86,7 @@ namespace dcmqi { data["InstanceNumber"] = this->instanceNumber; data["BodyPartExamined"] = this->bodyPartExamined; - data["segmentAttributes"] = createAndGetSegmentAttributes(); + data["segmentAttributes"] = createAndGetSegmentAttributesJSON(); Json::StyledWriter styledWriter; ss << styledWriter.write(data); @@ -94,11 +94,15 @@ namespace dcmqi { return ss.str(); } - Json::Value JSONSegmentationMetaInformationHandler::createAndGetSegmentAttributes() { + Json::Value JSONSegmentationMetaInformationHandler::createAndGetSegmentAttributesJSON() { // return a list of lists, where each inner list contains just one item (segment) Json::Value values(Json::arrayValue); + for (vector >::const_iterator vIt = this->segmentsAttributesMappingList.begin(); vIt != this->segmentsAttributesMappingList.end(); ++vIt) { + + Json::Value innerList(Json::arrayValue); + for(map::const_iterator mIt=vIt->begin();mIt!=vIt->end();++mIt){ Json::Value segment; SegmentAttributes* segmentAttributes = mIt->second; @@ -139,28 +143,30 @@ namespace dcmqi { rgb.append(segmentAttributes->getRecommendedDisplayRGBValue()[1]); rgb.append(segmentAttributes->getRecommendedDisplayRGBValue()[2]); segment["recommendedDisplayRGBValue"] = rgb; - Json::Value innerList(Json::arrayValue); innerList.append(segment); - values.append(innerList); } + + values.append(innerList); } return values; } - SegmentAttributes *JSONSegmentationMetaInformationHandler::createAndGetNewSegment(unsigned labelID) { - for (vector >::const_iterator vIt = this->segmentsAttributesMappingList.begin(); - vIt != this->segmentsAttributesMappingList.end(); ++vIt) { - for(map::const_iterator mIt = vIt->begin();mIt!=vIt->end();++mIt){ - SegmentAttributes *segmentAttributes = mIt->second; - if (segmentAttributes->getLabelID() == labelID) - return NULL; - } + SegmentAttributes *JSONSegmentationMetaInformationHandler::createOrGetSegment(const unsigned int segGroupNumber, const unsigned labelID) { + if (segGroupNumber < 1) + { + std::cerr << "ERROR: Segment group number must be >= 1" << std::endl; + return NULL; } SegmentAttributes *segment = new SegmentAttributes(labelID); - map tempMap; - tempMap[labelID] = segment; - this->segmentsAttributesMappingList.push_back(tempMap); + if (this->segmentsAttributesMappingList.size() < segGroupNumber) { + map tempMap; + tempMap[labelID] = segment; + this->segmentsAttributesMappingList.push_back(tempMap); + } else { + this->segmentsAttributesMappingList[segGroupNumber-1][labelID] = segment; + } + return segment; } @@ -237,7 +243,6 @@ namespace dcmqi { Json::Value elem = segment["TrackingUniqueIdentifier"]; segmentAttribute->setTrackingUniqueIdentifier(elem.asString()); } - } segmentsAttributesMappingList.push_back(labelID2SegmentAttributes); } diff --git a/libsrc/OverlapUtil.cpp b/libsrc/OverlapUtil.cpp new file mode 100644 index 00000000..a5b77f6c --- /dev/null +++ b/libsrc/OverlapUtil.cpp @@ -0,0 +1,830 @@ +/* + * + * Copyright (C) 2023, Open Connections GmbH + * All rights reserved. See COPYRIGHT file for details. + * + * This software and supporting documentation were developed by + * + * OFFIS e.V. + * R&D Division Health + * Escherweg 2 + * D-26121 Oldenburg, Germany + * + * + * Module: dcmqi + * + * Author: Michael Onken + * + * Purpose: Interface of class OverlapUtil + * + */ + +#include "dcmqi/OverlapUtil.h" +#include "dcmtk/dcmdata/dcerror.h" +#include "dcmtk/dcmfg/fginterface.h" +#include "dcmtk/dcmfg/fgpixmsr.h" +#include "dcmtk/dcmfg/fgplanor.h" +#include "dcmtk/dcmfg/fgplanpo.h" +#include "dcmtk/dcmfg/fgseg.h" +#include "dcmtk/dcmfg/fgtypes.h" +#include "dcmtk/dcmiod/iodtypes.h" +#include "dcmtk/dcmseg/segdoc.h" +#include "dcmtk/dcmseg/segtypes.h" +#include "dcmtk/dcmseg/segutils.h" +#include "dcmtk/ofstd/ofcond.h" +#include "dcmtk/ofstd/ofstream.h" +#include "dcmtk/ofstd/oftimer.h" +#include "dcmtk/ofstd/oftypes.h" + +#include +#include + +makeOFConditionConst(SG_EC_FramesNotParallel, OFM_dcmseg, 7, OF_error, "Frames are not parallel"); + +namespace dcmqi +{ + +OverlapUtil::OverlapUtil() + : m_imageOrientation() + , m_framePositions() + , m_framesForSegment() + , m_logicalFramePositions() + , m_segmentsByPosition() + , m_segmentOverlapMatrix(0) + , m_nonOverlappingSegments() + , m_seg() +{ +} + +OverlapUtil::~OverlapUtil() +{ + // nothing to do +} + +void OverlapUtil::setSegmentationObject(DcmSegmentation* seg) +{ + m_seg = seg; + clear(); +} + +void OverlapUtil::clear() +{ + m_imageOrientation.clear(); + m_framePositions.clear(); + m_framesForSegment.clear(); + m_logicalFramePositions.clear(); + m_segmentsByPosition.clear(); + m_segmentOverlapMatrix.clear(); + m_nonOverlappingSegments.clear(); +} + +OFCondition OverlapUtil::getFramesByPosition(DistinctFramePositions& result) +{ + OFCondition cond; + if (m_logicalFramePositions.empty()) + { + cond = groupFramesByPosition(); + } + if (cond.good()) + { + result = m_logicalFramePositions; + } + return cond; +} + +OFCondition OverlapUtil::getFramesForSegment(const Uint32 segmentNumber, OFVector& frames) +{ + if ((segmentNumber == 0) || (segmentNumber > m_seg->getNumberOfSegments() + 1)) + { + DCMSEG_ERROR("getFramesForSegment(): Segment number " << segmentNumber << " is out of range"); + return EC_IllegalParameter; + } + if (m_framesForSegment.empty()) + { + FGInterface& fg = m_seg->getFunctionalGroups(); + size_t tempNum = m_seg->getNumberOfFrames(); + if (tempNum > 4294967295) + { + DCMSEG_ERROR("getFramesForSegment(): Number of frames " << tempNum << " exceeds maximum number of possible frames (2^32-1)"); + return EC_IllegalParameter; + } + Uint32 numFrames = static_cast(m_seg->getNumberOfFrames()); + m_framesForSegment.clear(); + m_framesForSegment.resize(m_seg->getNumberOfSegments()); + // Get Segmentation FG for each frame and remember the segment number for each frame + // in the vector m_segmentsForFrame + for (Uint32 f = 0; f < numFrames; f++) + { + FGBase* group = NULL; + FGSegmentation* segFG = NULL; + group = fg.get(f, DcmFGTypes::EFG_SEGMENTATION); + segFG = OFstatic_cast(FGSegmentation*, group); + if (segFG) + { + Uint16 segNum = 0; + OFCondition cond = segFG->getReferencedSegmentNumber(segNum); + if (cond.good() && segNum > 0) + { + m_framesForSegment[segNum - 1].push_back(f); // physical frame number for segment + } + else if (segNum == 0) + { + DCMSEG_WARN("getFramesForSegment(): Referenced Segment Number is 0 (not permitted) for frame #" + << f << ", ignoring"); + return EC_InvalidValue; + } + else + { + DCMSEG_ERROR( + "getFramesForSegment(): Referenced Segment Number not found (not permitted) for frame #" + << f << ", cannot add segment"); + return EC_TagNotFound; + } + } + } + } + frames = m_framesForSegment[segmentNumber - 1]; + return EC_Normal; +} + +OFCondition OverlapUtil::ensureFramesAreParallel() +{ + FGInterface& fg = m_seg->getFunctionalGroups(); + OFCondition cond; + OFBool perFrame = OFFalse; + FGPlaneOrientationPatient* pop = NULL; + // Ensure that Image Orientation Patient is shared, i.e. we have parallel frames + m_imageOrientation.clear(); + m_imageOrientation.resize(6); + FGBase* group = fg.get(0, DcmFGTypes::EFG_PLANEORIENTPATIENT, perFrame); + if (group && (pop = OFstatic_cast(FGPlaneOrientationPatient*, group))) + { + if (perFrame == OFFalse) + { + DCMSEG_DEBUG("ensureFramesAreParallel(): Image Orientation Patient is shared, frames are parallel"); + m_imageOrientation.resize(6); + cond = pop->getImageOrientationPatient(m_imageOrientation[0], + m_imageOrientation[1], + m_imageOrientation[2], + m_imageOrientation[3], + m_imageOrientation[4], + m_imageOrientation[5]); + std::cout << "Image Orientation Patient set to : " << m_imageOrientation[0] << ", " << m_imageOrientation[1] + << ", " << m_imageOrientation[2] << ", " << m_imageOrientation[3] << ", " << m_imageOrientation[4] + << ", " << m_imageOrientation[5] << std::endl; + return cond; + } + else + { + DCMSEG_ERROR( + "ensureFramesAreParallel(): Image Orientation Patient is per-frame, frames are probably not parallel"); + return SG_EC_FramesNotParallel; + } + } + else + { + DCMSEG_ERROR( + "ensureFramesAreParallel(): Plane Orientation (Patient) FG not found, cannot check for parallel frames"); + return EC_TagNotFound; + } + return EC_Normal; +} + +OFCondition OverlapUtil::groupFramesByPosition() +{ + if (!m_framePositions.empty()) + { + // Already computed + return EC_Normal; + } + + OFCondition cond = ensureFramesAreParallel(); + if (cond.bad()) + { + return cond; + } + + OFTimer tm; + + // Group all frames by position into m_logicalFramePositions. + // After that, all frames at the same position will be in the same vector + // assigned to the same key (the frame's coordinates) in the map. + // Group all frames by position into m_logicalFramePositions. + cond = collectPhysicalFramePositions(); + if (cond.good()) + { + cond = groupFramesByLogicalPosition(); + } + + // print frame groups if debug log level is enabled: + if (cond.good() && DCM_dcmsegLogger.isEnabledFor(OFLogger::DEBUG_LOG_LEVEL)) + { + DCMSEG_DEBUG("groupFramesByPosition(): Frames grouped by position:"); + for (size_t i = 0; i < m_logicalFramePositions.size(); ++i) + { + OFStringStream ss; + for (size_t j = 0; j < m_logicalFramePositions[i].size(); ++j) + { + if (j > 0) + ss << ", "; + ss << m_logicalFramePositions[i][j]; + } + DCMSEG_DEBUG("groupFramesByPosition(): Logical frame #" << i << ": " << ss.str()); + } + } + DCMSEG_DEBUG("groupFramesByPosition(): Grouping frames by position took " << tm.getDiff() << " s"); + + if (cond.bad()) + { + m_framePositions.clear(); + m_logicalFramePositions.clear(); + } + return cond; +} + +OFCondition OverlapUtil::getSegmentsByPosition(SegmentsByPosition& result) +{ + if (!m_segmentsByPosition.empty()) + { + // Already computed + result = m_segmentsByPosition; + return EC_Normal; + } + // Make sure prerequisites are met + OFTimer tm; + OFCondition cond = groupFramesByPosition(); + if (cond.bad()) + { + return cond; + } + size_t numSegments = m_seg->getNumberOfSegments(); + if (m_logicalFramePositions.empty()) + { + cond = getFramesByPosition(m_logicalFramePositions); + if (cond.bad()) + return cond; + } + m_segmentsByPosition.clear(); + m_segmentsByPosition.resize(m_logicalFramePositions.size()); + for (size_t l = 0; l < m_logicalFramePositions.size(); ++l) + { + OFVector segments; + for (size_t f = 0; f < m_logicalFramePositions[l].size(); ++f) + { + Uint32 frameNumber = m_logicalFramePositions[l][f]; + OFVector segs; + FGBase* group = NULL; + FGSegmentation* segFG = NULL; + group = m_seg->getFunctionalGroups().get(frameNumber, DcmFGTypes::EFG_SEGMENTATION); + segFG = OFstatic_cast(FGSegmentation*, group); + if (segFG) + { + Uint16 segNum = 0; + cond = segFG->getReferencedSegmentNumber(segNum); + if (cond.good() && segNum > 0 && (segNum <= numSegments)) + { + m_segmentsByPosition[l].insert(SegNumAndFrameNum(segNum, frameNumber)); + } + else if (segNum == 0) + { + DCMSEG_ERROR( + "getSegmentsByPosition(): Referenced Segment Number is 0 (not permitted), cannot add segment"); + cond = EC_InvalidValue; + break; + } + else if (segNum > numSegments) + { + DCMSEG_ERROR("getSegmentsByPosition(): Found Referenced Segment Number " + << segNum << " but only " << numSegments + << " segments are present, cannot add segment"); + DCMSEG_ERROR( + "getSegmentsByPosition(): Segments are not numbered consecutively, cannot add segment"); + cond = EC_InvalidValue; + break; + } + else + { + DCMSEG_ERROR("getSegmentsByPosition(): Referenced Segment Number not found (not permitted) , " + "cannot add segment"); + cond = EC_TagNotFound; + break; + } + } + } + if (cond.bad()) + { + break; + } + } + // print segments per logical frame if debug log level is enabled + if (cond.good() && DCM_dcmsegLogger.isEnabledFor(OFLogger::DEBUG_LOG_LEVEL)) + { + OFStringStream ss; + printSegmentsByPosition(ss); + DCMSEG_DEBUG(ss.str()); + } + DCMSEG_DEBUG("groupFramesByPosition(): Grouping segments by position took " << tm.getDiff() << " s"); + return cond; +} + +OFCondition OverlapUtil::getOverlapMatrix(OverlapMatrix& matrix) +{ + if (!m_segmentOverlapMatrix.empty()) + { + // Already computed + matrix = m_segmentOverlapMatrix; + return EC_Normal; + } + // Make sure prerequisites are met + OFTimer tm; + SegmentsByPosition dontCare; + OFCondition result = getSegmentsByPosition(dontCare); + if (result.good()) + { + result = buildOverlapMatrix(); + } + if (result.good()) + { + matrix = m_segmentOverlapMatrix; + } + DCMSEG_DEBUG("getOverlappingSegments(): Building overlap matrix took " << tm.getDiff() << " s"); + return result; +} + +OFCondition OverlapUtil::getNonOverlappingSegments(SegmentGroups& segmentGroups) +{ + OFTimer tm; + OFCondition result; + if (!m_nonOverlappingSegments.empty()) + { + // Already computed + segmentGroups = m_nonOverlappingSegments; + return EC_Normal; + } + // Make sure prerequisites are met + result = getOverlapMatrix(m_segmentOverlapMatrix); + if (result.good()) + { + // Group those segments from the overlap matrix together, that do not + // overlap with each other. + // Go through all segments and check if they overlap with any of the already + // grouped segments. If not, add them to the same group. If yes, create a new group + // and add them there. + m_nonOverlappingSegments.push_back(OFVector()); + for (size_t i = 0; i < m_segmentOverlapMatrix.size(); ++i) + { + // Loop over all groups and check whether the current segment overlaps with any of them + OFBool overlaps = OFFalse; + for (size_t j = 0; j < m_nonOverlappingSegments.size(); ++j) + { + // Loop over all segments in the current group + for (OFVector::iterator it = m_nonOverlappingSegments[j].begin(); + it != m_nonOverlappingSegments[j].end(); + ++it) + { + // Check if the current segment overlaps with the segment in the current group + if (m_segmentOverlapMatrix[i][(*it) - 1] == 1) + { + overlaps = OFTrue; + break; + } + } + if (!overlaps) + { + // Add segment to current group + m_nonOverlappingSegments[j].push_back(i + 1); + break; + } + } + if (overlaps) + { + // Create new group and add segment to it + m_nonOverlappingSegments.push_back(OFVector()); + m_nonOverlappingSegments.back().push_back(i + 1); + } + } + } + DCMSEG_DEBUG("getNonOverlappingSegments(): Grouping non-overlapping segments took " << tm.getDiff() << " s"); + if (result.good()) + { + // print non-overlapping segments if debug log level is enabled + if (DCM_dcmsegLogger.isEnabledFor(OFLogger::DEBUG_LOG_LEVEL)) + { + OFStringStream ss; + printNonOverlappingSegments(ss); + DCMSEG_DEBUG(ss.str()); + } + } + if (result.good()) + { + segmentGroups = m_nonOverlappingSegments; + } + return result; +} + +void OverlapUtil::printSegmentsByPosition(OFStringStream& ss) +{ + ss << "printSegmentsByPosition(): Segments grouped by logical frame positions, (seg#,frame#):" << OFendl; + for (size_t i = 0; i < m_segmentsByPosition.size(); ++i) + { + OFStringStream tempSS; + for (std::set::iterator it = m_segmentsByPosition[i].begin(); + it != m_segmentsByPosition[i].end(); + ++it) + { + if (it != m_segmentsByPosition[i].begin()) + tempSS << ","; + tempSS << "(" << (*it).m_segmentNumber << "," << (*it).m_frameNumber << ")"; + } + ss << "printSegmentsByPosition(): Logical frame #" << i << ": " << tempSS.str(); + } +} + +void OverlapUtil::printOverlapMatrix(OFStringStream& ss) +{ + ss << "printOverlapMatrix(): Overlap matrix:" << OFendl; + for (size_t i = 0; i < m_segmentOverlapMatrix.size(); ++i) + { + for (size_t j = 0; j < m_segmentOverlapMatrix[i].size(); ++j) + { + if (m_segmentOverlapMatrix[i][j] >= 0) + ss << OFstatic_cast(Uint32, m_segmentOverlapMatrix[i][j]); + else + ss << "1"; + ss << " "; + } + ss << OFendl; + } +} + +void OverlapUtil::printNonOverlappingSegments(OFStringStream& ss) +{ + ss << "printNonOverlappingSegments(): Non-overlapping segments:" << OFendl; + for (size_t i = 0; i < m_nonOverlappingSegments.size(); ++i) + { + ss << "Group #" << i << ": "; + for (OFVector::iterator it = m_nonOverlappingSegments[i].begin(); + it != m_nonOverlappingSegments[i].end(); + ++it) + { + if (it != m_nonOverlappingSegments[i].begin()) + ss << ", "; + ss << (*it); + } + ss << OFendl; + } +} + +OFCondition OverlapUtil::buildOverlapMatrix() +{ + // Make 2 dimensional array matrix of Sint8 type for (segment numbers) X (segment numbers). + // Initialize with -1 (not checked yet) + m_segmentOverlapMatrix.clear(); + m_segmentOverlapMatrix.resize(m_seg->getNumberOfSegments()); + for (size_t i = 0; i < m_segmentOverlapMatrix.size(); ++i) + { + m_segmentOverlapMatrix[i].resize(m_seg->getNumberOfSegments(), -1); + } + // Diagonal is always 0 (segment does not interfere/overlap with itself) + for (size_t i = 0; i < m_segmentOverlapMatrix.size(); ++i) + { + m_segmentOverlapMatrix[i][i] = 0; + } + // Go through all logical frame positions, and compare all segments at each position + size_t index1, index2; + index1 = index2 = 0; + for (size_t i = 0; i < m_segmentsByPosition.size(); ++i) + { + DCMSEG_DEBUG("getOverlappingSegments(): Comparing segments at logical frame position " << i); + // Compare all segments at this position + for (std::set::iterator it = m_segmentsByPosition[i].begin(); + it != m_segmentsByPosition[i].end(); + ++it) + { + index1++; + for (std::set::iterator it2 = m_segmentsByPosition[i].begin(); + it2 != m_segmentsByPosition[i].end(); + ++it2) + { + index2++; + // Skip comparison of same segments in reverse order (index2 < index1) + if (index2 <= index1) + continue; + // Skip self-comparison (diagonal is always 0); (index1==index2) + if (it->m_segmentNumber != it2->m_segmentNumber) + { + // Check if we already have found an overlap on another logical frame, and if so, skip + Sint8 existing_result + = m_segmentOverlapMatrix[(*it).m_segmentNumber - 1][(*it2).m_segmentNumber - 1]; + if (existing_result == 1) + { + DCMSEG_DEBUG("getOverlappingSegments(): Skipping frame comparison on pos #" + << i << " for segments " << (*it).m_segmentNumber << " and " + << (*it2).m_segmentNumber << " (already marked as overlapping)"); + continue; + } + // Compare pixels of the frames referenced by each segments. + // If they overlap, mark as overlapping + OFBool overlap = OFFalse; + checkFramesOverlap(it->m_frameNumber, it2->m_frameNumber, overlap); + + // Enter result into overlap matrix + m_segmentOverlapMatrix[(*it).m_segmentNumber - 1][(*it2).m_segmentNumber - 1] = overlap ? 1 : 0; + m_segmentOverlapMatrix[(*it2).m_segmentNumber - 1][(*it).m_segmentNumber - 1] = overlap ? 1 : 0; + } + } + } + } + // Since we don't compare all segments (since not all are showing up together on a single logical frame), + // we set all remaining entries that are still not initialized (-1) to 0 (no overlap) + for (size_t i = 0; i < m_segmentOverlapMatrix.size(); ++i) + { + for (size_t j = 0; j < m_segmentOverlapMatrix[i].size(); ++j) + { + if (m_segmentOverlapMatrix[i][j] == -1) + { + m_segmentOverlapMatrix[i][j] = 0; + } + } + } + // print overlap matrix if debug log level is enabled + if (DCM_dcmsegLogger.isEnabledFor(OFLogger::DEBUG_LOG_LEVEL)) + { + OFStringStream ss; + printOverlapMatrix(ss); + DCMSEG_DEBUG(ss.str()); + } + return EC_Normal; +} + +OFCondition OverlapUtil::checkFramesOverlap(const Uint32& f1, const Uint32& f2, OFBool& overlap) +{ + if (f1 == f2) + { + // The same frame should not be considered overlapping at all + overlap = OFFalse; + return EC_Normal; + } + overlap = OFFalse; + OFCondition result; + const DcmIODTypes::Frame* f1_data = m_seg->getFrame(f1); + const DcmIODTypes::Frame* f2_data = m_seg->getFrame(f2); + Uint16 rows, cols; + rows = cols = 0; + DcmIODImage>* ip = static_cast>*>(m_seg); + ip->getImagePixel().getRows(rows); + ip->getImagePixel().getColumns(cols); + if (rows * cols % 8 != 0) + { + // We must compare pixel by pixel of the unpacked frames (for now) + result = checkFramesOverlapUnpacked(f1, f2, f1_data, f2_data, rows, cols, overlap); + } + else + { + // We can compare byte by byte using bitwise AND (if both have a 1 at the same position, they overlap) + result = checkFramesOverlapBinary(f1, f2, f1_data, f2_data, rows, cols, overlap); + } + if (result.good() && !overlap) + { + DCMSEG_DEBUG("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " don't overlap"); + } + return result; +} + +OFCondition OverlapUtil::checkFramesOverlapBinary(const Uint32& f1, + const Uint32& f2, + const DcmIODTypes::Frame* f1_data, + const DcmIODTypes::Frame* f2_data, + const Uint16& rows, + const Uint16 cols, + OFBool& overlap) +{ + DCMSEG_DEBUG("checkFramesOverlap(): Comparing frames " << f1 << " and " << f2 << " for overlap (fast binary mode)"); + if (!f1_data || !f2_data) + { + DCMSEG_ERROR("checkFramesOverlap(): Cannot access binary frames " << f1 << " and " << f2 << " for comparison"); + return EC_IllegalCall; + } + if (f1_data->length != f2_data->length) + { + DCMSEG_ERROR("checkFramesOverlap(): Frames " << f1 << " and " << f2 + << " have different length, cannot compare"); + return EC_IllegalCall; + } + // Compare byte (8 pixels at once) using bitwise AND (if both have a 1 at the same position, they overlap) + for (size_t n = 0; n < f1_data->length; ++n) + { + if (f1_data->pixData[n] & f2_data->pixData[n]) + { + DCMSEG_DEBUG("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " do overlap, pixel value " + << OFstatic_cast(Uint16, f1_data->pixData[n]) << " at index " + << n << " is the same"); + overlap = OFTrue; + break; + } + } + return EC_Normal; +} + +OFCondition OverlapUtil::checkFramesOverlapUnpacked(const Uint32& f1, + const Uint32& f2, + const DcmIODTypes::Frame* f1_data, + const DcmIODTypes::Frame* f2_data, + const Uint16& rows, + const Uint16 cols, + OFBool& overlap) +{ + DCMSEG_DEBUG("checkFramesOverlap(): Comparing frames " << f1 << " and " << f2 + << " for overlap (slow unpacked mode)"); + OFunique_ptr f1_unpacked(DcmSegUtils::unpackBinaryFrame(f1_data, rows, cols)); + OFunique_ptr f2_unpacked(DcmSegUtils::unpackBinaryFrame(f2_data, rows, cols)); + if (!f1_unpacked || !f2_unpacked) + { + DCMSEG_ERROR("checkFramesOverlap(): Cannot unpack frames " << f1 << " and " << f2 << " for comparison"); + return EC_IllegalCall; + } + if (f1_unpacked->length != f2_unpacked->length) + { + DCMSEG_ERROR("checkFramesOverlap(): Frames " << f1 << " and " << f2 + << " have different length, cannot compare"); + return EC_IllegalCall; + } + // Compare pixels of both frames and check whether at least one has the same value + DCMSEG_DEBUG("checkFramesOverlap(): Comparing frames " << f1 << " and " << f2 << " for overlap"); + for (size_t n = 0; n < f1_unpacked->length; ++n) + { + if (f1_unpacked->pixData[n] != 0 && (f1_unpacked->pixData[n] == f2_unpacked->pixData[n])) + { + DCMSEG_DEBUG("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " do overlap, pixel value " + << OFstatic_cast(Uint16, f1_unpacked->pixData[n]) + << " at index " << n << " is the same"); + overlap = OFTrue; + break; + } + } + return EC_Normal; +} + +OFCondition OverlapUtil::collectPhysicalFramePositions() +{ + // Group all frames by position into m_logicalFramePositions. + FGInterface& fg = m_seg->getFunctionalGroups(); + size_t numFrames = m_seg->getNumberOfFrames(); + OFBool perFrame = OFFalse; + OFCondition cond; + // Vector of frame numbers with their respective position + m_framePositions.clear(); + m_framePositions.reserve(numFrames); + + // Put all frames into vector along with their Image Position Patient coordinates + for (size_t i = 0; i < numFrames; ++i) + { + FGPlanePosPatient* ppp = NULL; + FGBase* group = fg.get(i, DcmFGTypes::EFG_PLANEPOSPATIENT, perFrame); + if (group) + ppp = OFstatic_cast(FGPlanePosPatient*, group); + if (ppp) + { + // Get image position patient for frame i + OFVector ipp(3); + // Only in later DCMTK version: + // cond = ppp->getImagePositionPatient(ipp); + cond = ppp->getImagePositionPatient(ipp[0], ipp[1], ipp[2]); + if (cond.good()) + { + // Insert frame into map + m_framePositions.push_back(FramePositionAndNumber(ipp, i)); + } + else + { + DCMSEG_ERROR("groupFramesByPosition(): Image Position Patient not found for frame " + << i << ", cannot sort frames by position"); + cond = EC_TagNotFound; + break; + } + } + else + { + DCMSEG_ERROR("groupFramesByPosition(): Image Position Patient not found for frame " + << i << ", cannot sort frames by position"); + cond = EC_TagNotFound; + break; + } + } + return cond; +} + +OFCondition OverlapUtil::groupFramesByLogicalPosition() +{ + OFCondition cond; + FGInterface& fg = m_seg->getFunctionalGroups(); + OFBool perFrame = OFFalse; + + // Find all distinct positions and for each position the actual frames that can be found at it + Float64 sliceThickness = 0.0; + FGPixelMeasures* pm = NULL; + Uint8 relevantCoordinate = 3; // not determined yet + FGBase* group = fg.get(0, DcmFGTypes::EFG_PIXELMEASURES, perFrame); + if (group && (pm = OFstatic_cast(FGPixelMeasures*, group))) + { + // Get/compute Slice Thickness + cond = pm->getSliceThickness(sliceThickness); + if (cond.good()) + { + DCMSEG_DEBUG("groupFramesByPosition(): Slice Thickness is " << sliceThickness); + // Identify coordinate to be used for frame sorting + relevantCoordinate = identifyChangingCoordinate(m_imageOrientation); + if (relevantCoordinate < 3) + { + DCMSEG_DEBUG("Using coordinate " << OFstatic_cast(Uint16, relevantCoordinate) + << " for sorting frames by position"); + ComparePositions c(relevantCoordinate); + std::sort(m_framePositions.begin(), m_framePositions.end(), c); + // vec will contain all frame numbers that are at the same position + OFVector vec; + vec.push_back(m_framePositions[0].m_frameNumber); + m_logicalFramePositions.push_back(vec); // Initialize for first logical frame + for (size_t j = 1; j < m_framePositions.size(); ++j) + { + // If frame is close to previous frame, add it to the same vector. + // 2.5 is chosen since it means the frames are not further away if clearly less than half a slice + // thickness + Float64 diff = fabs(m_framePositions[j].m_position[relevantCoordinate] + - m_framePositions[j - 1].m_position[relevantCoordinate]); + DCMSEG_DEBUG("Coordinates of both frames:"); + DCMSEG_DEBUG("Frame " << j << ": " << m_framePositions[j].m_position[0] << ", " + << m_framePositions[j].m_position[1] << ", " + << m_framePositions[j].m_position[2]); + DCMSEG_DEBUG("Frame " << j - 1 << ": " << m_framePositions[j - 1].m_position[0] << ", " + << m_framePositions[j - 1].m_position[1] << ", " + << m_framePositions[j - 1].m_position[2]); + DCMSEG_DEBUG("groupFramesByPosition(): Frame " << j << " is " << diff + << " mm away from previous frame"); + // 1% inaccuracy for slice thickness will be considered as same logical position + if (diff < sliceThickness * 0.01) + { + // Add frame to last vector + DCMSEG_DEBUG("Assigning to same frame bucket as previous frame"); + m_logicalFramePositions.back().push_back( + m_framePositions[j].m_frameNumber); // physical frame number + } + else + { + DCMSEG_DEBUG("Assigning to same new frame bucket"); + // Create new vector + OFVector vec; + vec.push_back(m_framePositions[j].m_frameNumber); + m_logicalFramePositions.push_back(vec); + } + } + } + else + { + std::cerr + << "groupFramesByPosition(): Cannot identify coordinate relevant for sorting frames by position" + << std::endl; + cond = EC_InvalidValue; + } + } + else + { + std::cerr << "groupFramesByPosition(): Slice Thickness not found, cannot sort frames by position" + << std::endl; + cond = EC_TagNotFound; + } + } + else + { + std::cerr << "groupFramesByPosition(): Pixel Measures FG not found, cannot sort frames by position" + << std::endl; + cond = EC_TagNotFound; + } + return cond; +} + +Uint8 OverlapUtil::identifyChangingCoordinate(const OFVector& imageOrientation) +{ + Float64 cross_product[3]; + // Compute cross product of image orientation vectors. + // We are only interested into the absolute values for later comparison + cross_product[0] = fabs(imageOrientation[1] * imageOrientation[5] - imageOrientation[2] * imageOrientation[4]); + cross_product[1] = fabs(imageOrientation[2] * imageOrientation[3] - imageOrientation[0] * imageOrientation[5]); + cross_product[2] = fabs(imageOrientation[0] * imageOrientation[4] - imageOrientation[1] * imageOrientation[3]); + // Find out which coordinate is changing the most (biggest absolute coordinate value of cross product) + if ((cross_product[0] > cross_product[1]) && (cross_product[0] > cross_product[2])) + { + return 0; + } + if ((cross_product[1] > cross_product[0]) && (cross_product[1] > cross_product[2])) + { + return 1; + } + if ((cross_product[2] > cross_product[0]) && (cross_product[2] > cross_product[1])) + { + return 2; + } + // No clear winner + return 3; +} + +} diff --git a/libsrc/ParaMapConverter.cpp b/libsrc/ParaMapConverter.cpp index ef8da16f..4185b1ea 100644 --- a/libsrc/ParaMapConverter.cpp +++ b/libsrc/ParaMapConverter.cpp @@ -5,7 +5,6 @@ // DCMQI includes #include "dcmqi/ParaMapConverter.h" -#include "dcmqi/ImageSEGConverter.h" // DCMTK includes #include