diff --git a/Modules/DICOM/include/mitkDICOMDatasetSorter.h b/Modules/DICOM/include/mitkDICOMDatasetSorter.h index c2b8ad4066..98c294b556 100644 --- a/Modules/DICOM/include/mitkDICOMDatasetSorter.h +++ b/Modules/DICOM/include/mitkDICOMDatasetSorter.h @@ -1,91 +1,95 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkDICOMDatasetSorter_h #define mitkDICOMDatasetSorter_h #include "itkObjectFactory.h" #include "mitkCommon.h" #include "mitkDICOMDatasetAccess.h" +#include "mitkDICOMSplitReason.h" namespace mitk { /** \ingroup DICOMModule \brief The sorting/splitting building-block of DICOMITKSeriesGDCMReader. This class describes the interface of the sorting/splitting process described as part of DICOMITKSeriesGDCMReader::AnalyzeInputFiles() (see \ref DICOMITKSeriesGDCMReader_LoadingStrategy). The procedure is simple: - take a list of input datasets (DICOMDatasetAccess) - sort them (to be defined by sub-classes, based on specific tags) - return the sorting result as outputs (the single input might be distributed into multiple outputs) The simplest and most generic form of sorting is implemented in sub-class DICOMTagBasedSorter. */ class MITKDICOM_EXPORT DICOMDatasetSorter : public itk::LightObject { public: mitkClassMacroItkParent( DICOMDatasetSorter, itk::LightObject ); /** \brief Return the tags of interest (to facilitate scanning) */ virtual DICOMTagList GetTagsOfInterest() = 0; /// \brief Input for sorting void SetInput(DICOMDatasetList filenames); /// \brief Input for sorting const DICOMDatasetList& GetInput() const; /// \brief Sort input datasets into one or multiple outputs. virtual void Sort() = 0; /// \brief Output of the sorting process. unsigned int GetNumberOfOutputs() const; /// \brief Output of the sorting process. const DICOMDatasetList& GetOutput(unsigned int index) const; /// \brief Output of the sorting process. DICOMDatasetList& GetOutput(unsigned int index); + const DICOMSplitReason* GetSplitReason(unsigned int index) const; + /// \brief Print configuration details into stream. virtual void PrintConfiguration(std::ostream& os, const std::string& indent = "") const = 0; virtual bool operator==(const DICOMDatasetSorter& other) const = 0; protected: DICOMDatasetSorter(); ~DICOMDatasetSorter() override; DICOMDatasetSorter(const DICOMDatasetSorter& other); DICOMDatasetSorter& operator=(const DICOMDatasetSorter& other); void ClearOutputs(); void SetNumberOfOutputs(unsigned int numberOfOutputs); - void SetOutput(unsigned int index, const DICOMDatasetList& output); + void SetOutput(unsigned int index, const DICOMDatasetList& output, const DICOMSplitReason* splitReason = nullptr); private: DICOMDatasetList m_Input; std::vector< DICOMDatasetList > m_Outputs; + std::vector< DICOMSplitReason::Pointer > m_SplitReasons; }; } #endif diff --git a/Modules/DICOM/include/mitkDICOMITKSeriesGDCMReader.h b/Modules/DICOM/include/mitkDICOMITKSeriesGDCMReader.h index 62b242f8d3..1abb94f66f 100644 --- a/Modules/DICOM/include/mitkDICOMITKSeriesGDCMReader.h +++ b/Modules/DICOM/include/mitkDICOMITKSeriesGDCMReader.h @@ -1,377 +1,378 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkDICOMITKSeriesGDCMReader_h #define mitkDICOMITKSeriesGDCMReader_h #include #include #include "mitkDICOMFileReader.h" #include "mitkDICOMDatasetSorter.h" #include "mitkDICOMGDCMImageFrameInfo.h" #include "mitkEquiDistantBlocksSorter.h" #include "mitkNormalDirectionConsistencySorter.h" #include "MitkDICOMExports.h" namespace itk { class TimeProbesCollectorBase; } namespace mitk { /** \ingroup DICOMModule \brief Flexible reader based on itk::ImageSeriesReader and GDCM, for single-slice modalities like CT, MR, PET, CR, etc. Implements the loading processed as structured by DICOMFileReader offers configuration of its loading strategy. Documentation sections: - \ref DICOMITKSeriesGDCMReader_LoadingStrategy - \ref DICOMITKSeriesGDCMReader_ForcedConfiguration - \ref DICOMITKSeriesGDCMReader_UserConfiguration - \ref DICOMITKSeriesGDCMReader_GantryTilt - \ref DICOMITKSeriesGDCMReader_Testing - \ref DICOMITKSeriesGDCMReader_Internals - \ref DICOMITKSeriesGDCMReader_RelatedClasses - \ref DICOMITKSeriesGDCMReader_TiltInternals - \ref DICOMITKSeriesGDCMReader_Condensing \section DICOMITKSeriesGDCMReader_LoadingStrategy Loading strategy The set of input files is processed by a number of DICOMDatasetSorter objects which may do two sort of things: 1. split a list of input frames into multiple lists, based on DICOM tags such as "Rows", "Columns", which cannot be mixed within a single mitk::Image 2. sort the frames within the input lists, based on the values of DICOM tags such as "Image Position Patient" When the DICOMITKSeriesGDCMReader is configured with DICOMDatasetSorter%s, the list of input files is processed as follows: 1. build an initial set of output groups, simply by grouping all input files. 2. for each configured DICOMDatasetSorter, process: - for each output group: 1. set this group's files as input to the sorter 2. let the sorter sort (and split) 3. integrate the sorter's output groups with our own output groups \section DICOMITKSeriesGDCMReader_ForcedConfiguration Forced Configuration In all cases, the reader will add two DICOMDatasetSorter objects that are required to load mitk::Images properly via itk::ImageSeriesReader: 1. As a \b first step, the input files will be split into groups that are not compatible because they differ in essential aspects: - (0028,0010) Number of Rows - (0028,0011) Number of Columns - (0028,0030) Pixel Spacing - (0018,1164) Imager Pixel Spacing - (0020,0037) %Image Orientation (Patient) - (0018,0050) Slice Thickness - (0028,0008) Number of Frames 2. As are two forced \b last steps: 1. There will always be an instance of EquiDistantBlocksSorter, which ensures that there is an equal distance between all the frames of an Image. This is required to achieve correct geometrical positions in the mitk::Image, i.e. it is essential to be able to make measurements in images. - whether or not the distance is required to be orthogonal to the image planes is configured by SetFixTiltByShearing(). - during this check, we need to tolerate some minor errors in documented vs. calculated image origins. The amount of tolerance can be adjusted by SetToleratedOriginOffset() and SetToleratedOriginOffsetToAdaptive(). Please see EquiDistantBlocksSorter for more details. The default should be good for most cases. 2. There is always an instance of NormalDirectionConsistencySorter, which makes the order of images go along the image normals (see NormalDirectionConsistencySorter) \section DICOMITKSeriesGDCMReader_UserConfiguration User Configuration The user of this class can add more sorting steps (similar to the one described in above section) by calling AddSortingElement(). Usually, an application will add sorting by "Image Position Patient", by "Instance Number", and by other relevant tags here. \section DICOMITKSeriesGDCMReader_GantryTilt Gantry tilt handling When CT gantry tilt is used, the gantry plane (= X-Ray source and detector ring) and the vertical plane do not align anymore. This scanner feature is used for example to reduce metal artifacts (e.g. Lee C , Evaluation of Using CT Gantry Tilt Scan on Head and Neck Cancer Patients with Dental Structure: Scans Show Less Metal Artifacts. Presented at: Radiological Society of North America 2011 Scientific Assembly and Annual Meeting; November 27- December 2, 2011 Chicago IL.). The acquired planes of such CT series do not match the expectations of a orthogonal geometry in mitk::Image: if you stack the slices, they show a small shift along the Y axis: \verbatim without tilt with tilt |||||| ////// |||||| ////// -- |||||| --------- ////// -------- table orientation |||||| ////// |||||| ////// Stacked slices: without tilt with tilt -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- \endverbatim As such gemetries do not "work" in conjunction with mitk::Image, DICOMITKSeriesGDCMReader is able to perform a correction for such series. Whether or not such correction should be attempted is controlled by SetFixTiltByShearing(), the default being correction. For details, see "Internals" below. \section DICOMITKSeriesGDCMReader_Testing Testing A number of tests is implemented in module DICOMTesting, which is documented at \ref DICOMTesting. \section DICOMITKSeriesGDCMReader_Internals Class internals Internally, the class is based on GDCM and it depends heavily on the gdcm::Scanner class. Since the sorting elements (see DICOMDatasetSorter and DICOMSortCriterion) can access tags only via the DICOMDatasetAccess interface, BUT DICOMITKSeriesGDCMReader holds a list of more specific classes DICOMGDCMImageFrameInfo, we must convert between the two types sometimes. This explains the methods ToDICOMDatasetList(), FromDICOMDatasetList(). The intermediate result of all the sorting efforts is held in m_SortingResultInProgress, which is modified through InternalExecuteSortingStep(). \subsection DICOMITKSeriesGDCMReader_RelatedClasses Overview of related classes The following diagram gives an overview of the related classes: \image html implementeditkseriesgdcmreader.jpg \subsection DICOMITKSeriesGDCMReader_TiltInternals Details about the tilt correction The gantry tilt "correction" algorithm fixes two errors introduced by ITK's ImageSeriesReader: - the plane shift that is ignored by ITK's reader is recreated by applying a shearing transformation using itk::ResampleFilter. - the spacing is corrected (it is calculated by ITK's reader from the distance between two origins, which is NOT the slice distance in this special case) Both errors are introduced in itkImageSeriesReader.txx (ImageSeriesReader::GenerateOutputInformation(void)), lines 176 to 245 (as of ITK 3.20) For the correction, we examine two consecutive slices of a series, both described as a pair (origin/orientation): - we calculate if the first origin is on a line along the normal of the second slice - if this is not the case, the geometry will not fit a normal mitk::Image/mitk::%Geometry3D - we then project the second origin into the first slice's coordinate system to quantify the shift - both is done in class GantryTiltInformation with quite some comments. The geometry of image stacks with tilted geometries is illustrated below: - green: the DICOM images as described by their tags: origin as a point with the line indicating the orientation - red: the output of ITK ImageSeriesReader: wrong, larger spacing, no tilt - blue: how much a shear must correct \image html Modules/DICOM/doc/Doxygen/tilt-correction.jpg \subsection DICOMITKSeriesGDCMReader_Condensing Sub-classes can condense multiple blocks into a single larger block The sorting/splitting process described above is helpful for at least two more DICOM readers, which either try to load 3D+t images or which load diffusion data. In both cases, a single pixel of the mitk::Image is made up of multiple values, in one case values over time, in the other case multiple measurements of a single point. The specialized readers for these cases (e.g. ThreeDnTDICOMSeriesReader) can reuse most of the methods in DICOMITKSeriesGDCMReader, except that they need an extra step after the usual sorting, in which they can merge already grouped 3D blocks. What blocks are merged depends on the specialized reader's understanding of these images. To allow for such merging, a method Condense3DBlocks() is called as an absolute last step of AnalyzeInputFiles(). Given this, a sub-class could implement only LoadImages() and Condense3DBlocks() instead repeating most of AnalyzeInputFiles(). */ class MITKDICOM_EXPORT DICOMITKSeriesGDCMReader : public DICOMFileReader { public: mitkClassMacro( DICOMITKSeriesGDCMReader, DICOMFileReader ); mitkCloneMacro( DICOMITKSeriesGDCMReader ); itkFactorylessNewMacro( DICOMITKSeriesGDCMReader ); mitkNewMacro1Param( DICOMITKSeriesGDCMReader, unsigned int ); mitkNewMacro2Param( DICOMITKSeriesGDCMReader, unsigned int, bool ); /** \brief Runs the sorting / splitting process described in \ref DICOMITKSeriesGDCMReader_LoadingStrategy. Method required by DICOMFileReader. */ void AnalyzeInputFiles() override; // void AllocateOutputImages(); /** \brief Loads images using itk::ImageSeriesReader, potentially applies shearing to correct gantry tilt. */ bool LoadImages() override; // re-implemented from super-class bool CanHandleFile(const std::string& filename) override; /** \brief Add an element to the sorting procedure described in \ref DICOMITKSeriesGDCMReader_LoadingStrategy. */ virtual void AddSortingElement(DICOMDatasetSorter* sorter, bool atFront = false); typedef const std::list ConstSorterList; ConstSorterList GetFreelyConfiguredSortingElements() const; /** \brief Controls whether to "fix" tilted acquisitions by shearing the output (see \ref DICOMITKSeriesGDCMReader_GantryTilt). */ void SetFixTiltByShearing(bool on); bool GetFixTiltByShearing() const; /** \brief Controls whether groups of only two images are accepted when ensuring consecutive slices via EquiDistantBlocksSorter. */ void SetAcceptTwoSlicesGroups(bool accept) const; bool GetAcceptTwoSlicesGroups() const; /** \brief See \ref DICOMITKSeriesGDCMReader_ForcedConfiguration. */ void SetToleratedOriginOffsetToAdaptive(double fractionOfInterSliceDistanct = 0.3) const; /** \brief See \ref DICOMITKSeriesGDCMReader_ForcedConfiguration. */ void SetToleratedOriginOffset(double millimeters = 0.005) const; /** \brief Ignore all dicom tags that are non-essential for simple 3D volume import. */ void SetSimpleVolumeReading(bool read) { m_SimpleVolumeReading = read; }; /** \brief Ignore all dicom tags that are non-essential for simple 3D volume import. */ bool GetSimpleVolumeReading() { return m_SimpleVolumeReading; }; double GetToleratedOriginError() const; bool IsToleratedOriginOffsetAbsolute() const; double GetDecimalPlacesForOrientation() const; bool operator==(const DICOMFileReader& other) const override; DICOMTagPathList GetTagsOfInterest() const override; static int GetDefaultDecimalPlacesForOrientation() { return m_DefaultDecimalPlacesForOrientation; } static bool GetDefaultSimpleVolumeImport() { return m_DefaultSimpleVolumeImport; } static bool GetDefaultFixTiltByShearing() { return m_DefaultFixTiltByShearing; } protected: void InternalPrintConfiguration(std::ostream& os) const override; /// \brief Return active C locale - static std::string GetActiveLocale(); + static std::string GetActiveLocale(); /** \brief Remember current locale on stack, activate "C" locale. "C" locale is required for correct parsing of numbers by itk::ImageSeriesReader */ void PushLocale() const; /** \brief Activate last remembered locale from locale stack "C" locale is required for correct parsing of numbers by itk::ImageSeriesReader */ void PopLocale() const; const static int m_DefaultDecimalPlacesForOrientation = 5; const static bool m_DefaultSimpleVolumeImport = false; const static bool m_DefaultFixTiltByShearing = true; DICOMITKSeriesGDCMReader(unsigned int decimalPlacesForOrientation = m_DefaultDecimalPlacesForOrientation, bool simpleVolumeImport = m_DefaultSimpleVolumeImport); ~DICOMITKSeriesGDCMReader() override; DICOMITKSeriesGDCMReader(const DICOMITKSeriesGDCMReader& other); DICOMITKSeriesGDCMReader& operator=(const DICOMITKSeriesGDCMReader& other); - typedef std::vector SortingBlockList; + using SortingBlockListItemType = std::pair; + using SortingBlockList = std::vector ; /** \brief "Hook" for sub-classes, see \ref DICOMITKSeriesGDCMReader_Condensing \return REMAINING blocks */ virtual SortingBlockList Condense3DBlocks(SortingBlockList& resultOf3DGrouping); virtual DICOMTagCache::Pointer GetTagCache() const; void SetTagCache( const DICOMTagCache::Pointer& ) override; /// \brief Sorting step as described in \ref DICOMITKSeriesGDCMReader_LoadingStrategy - static SortingBlockList InternalExecuteSortingStep( + static SortingBlockList InternalExecuteSortingStep( unsigned int sortingStepIndex, const DICOMDatasetSorter::Pointer& sorter, const SortingBlockList& input); /// \brief Loads the mitk::Image by means of an itk::ImageSeriesReader virtual bool LoadMitkImageForOutput(unsigned int o); virtual bool LoadMitkImageForImageBlockDescriptor(DICOMImageBlockDescriptor& block) const; /// \brief Describe this reader's confidence for given SOP class UID - static ReaderImplementationLevel GetReaderImplementationLevel(const std::string sopClassUID); + static ReaderImplementationLevel GetReaderImplementationLevel(const std::string sopClassUID); private: /// \brief Creates the required sorting steps described in \ref DICOMITKSeriesGDCMReader_ForcedConfiguration void EnsureMandatorySortersArePresent(unsigned int decimalPlacesForOrientation, bool simpleVolumeImport = false); protected: // NOT nice, made available to ThreeDnTDICOMSeriesReader due to lack of time bool m_FixTiltByShearing; // could be removed by ITKDICOMSeriesReader NOT flagging tilt unless requested to fix it! bool m_SimpleVolumeReading; private: SortingBlockList m_SortingResultInProgress; typedef std::list SorterList; SorterList m_Sorter; protected: // NOT nice, made available to ThreeDnTDICOMSeriesReader and ClassicDICOMSeriesReader due to lack of time mitk::EquiDistantBlocksSorter::Pointer m_EquiDistantBlocksSorter; mitk::NormalDirectionConsistencySorter::Pointer m_NormalDirectionConsistencySorter; private: static std::mutex s_LocaleMutex; mutable std::stack m_ReplacedCLocales; mutable std::stack m_ReplacedCinLocales; double m_DecimalPlacesForOrientation; DICOMTagCache::Pointer m_TagCache; bool m_ExternalCache; }; } #endif diff --git a/Modules/DICOM/include/mitkDICOMImageBlockDescriptor.h b/Modules/DICOM/include/mitkDICOMImageBlockDescriptor.h index db2a634911..be839c3861 100644 --- a/Modules/DICOM/include/mitkDICOMImageBlockDescriptor.h +++ b/Modules/DICOM/include/mitkDICOMImageBlockDescriptor.h @@ -1,237 +1,247 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkDICOMImageBlockDescriptor_h #define mitkDICOMImageBlockDescriptor_h #include "mitkDICOMEnums.h" #include "mitkDICOMImageFrameInfo.h" #include "mitkDICOMTag.h" #include "mitkDICOMTagCache.h" +#include "mitkDICOMSplitReason.h" #include "mitkImage.h" #include "mitkProperties.h" #include "mitkWeakPointer.h" #include "mitkIPropertyProvider.h" #include "mitkGantryTiltInformation.h" #include namespace mitk { struct DICOMCachedValueInfo { unsigned int TimePoint; unsigned int SliceInTimePoint; std::string Value; }; class DICOMCachedValueLookupTable : public GenericLookupTable< DICOMCachedValueInfo > { public: typedef DICOMCachedValueLookupTable Self; typedef GenericLookupTable< DICOMCachedValueInfo > Superclass; const char *GetNameOfClass() const override { return "DICOMCachedValueLookupTable"; } DICOMCachedValueLookupTable() {} Superclass& operator=(const Superclass& other) override { return Superclass::operator=(other); } ~DICOMCachedValueLookupTable() override {} }; /** \ingroup DICOMModule \brief Output descriptor for DICOMFileReader. As a result of analysis by a mitk::DICOMFileReader, this class describes the properties of a single mitk::Images that could be loaded by the file reader. The descriptor contains the following information: - the mitk::Image itself. This will be nullptr after analysis and only be present after actual loading. - a list of frames (mostly: filenames) that went into composition of the mitk::Image. - an assessment of the reader's ability to load this set of files (ReaderImplementationLevel) - this can be used for reader selection when one reader is able to load an image with correct colors and the other is able to produce only gray values, for example - description of aspects of the image. Mostly a key-value list implemented by means of mitk::PropertyList. - for specific keys and possible values, see documentation of specific readers. \note an mitk::Image may both consist of multiple files (the "old" DICOM way) or a mitk::Image may be described by a single DICOM file or even only parts of a DICOM file (the newer multi-frame DICOM classes). To reflect this DICOMImageFrameList describes a list of frames from different or a single file. Described aspects of an image are: - whether pixel spacing is meant to be in-patient or on-detector (mitk::PixelSpacingInterpretation) - details about a possible gantry tilt (intended for use by file readers, may be hidden later) */ class MITKDICOM_EXPORT DICOMImageBlockDescriptor: public IPropertyProvider { public: DICOMImageBlockDescriptor(); ~DICOMImageBlockDescriptor() override; DICOMImageBlockDescriptor(const DICOMImageBlockDescriptor& other); DICOMImageBlockDescriptor& operator=(const DICOMImageBlockDescriptor& other); static DICOMTagList GetTagsOfInterest(); /// List of frames that constitute the mitk::Image (DICOMImageFrame%s) void SetImageFrameList(const DICOMImageFrameList& framelist); /// List of frames that constitute the mitk::Image (DICOMImageFrame%s) const DICOMImageFrameList& GetImageFrameList() const; /// The 3D mitk::Image that is loaded from the DICOM files of a DICOMImageFrameList void SetMitkImage(Image::Pointer image); /// the 3D mitk::Image that is loaded from the DICOM files of a DICOMImageFrameList Image::Pointer GetMitkImage() const; /// Reader's capability to appropriately load this set of frames ReaderImplementationLevel GetReaderImplementationLevel() const; /// Reader's capability to appropriately load this set of frames void SetReaderImplementationLevel(const ReaderImplementationLevel& level); /// Key-value store describing aspects of the image to be loaded void SetProperty(const std::string& key, BaseProperty* value); /// Key-value store describing aspects of the image to be loaded BaseProperty* GetProperty(const std::string& key) const; /// Convenience function around GetProperty() std::string GetPropertyAsString(const std::string&) const; + /** Returns the pointer to the split reason of this block descriptor.*/ + const DICOMSplitReason* GetSplitReason() const; + /** Returns the pointer to the split reason of this block descriptor.*/ + DICOMSplitReason* GetSplitReason(); + /** Sets the split reason for the block descriptor.*/ + void SetSplitReason(DICOMSplitReason* reason); + /// Convenience function around SetProperty() void SetFlag(const std::string& key, bool value); /// Convenience function around GetProperty() bool GetFlag(const std::string& key, bool defaultValue) const; /// Convenience function around SetProperty() void SetIntProperty(const std::string& key, int value); /// Convenience function around GetProperty() int GetIntProperty(const std::string& key, int defaultValue) const; BaseProperty::ConstPointer GetConstProperty(const std::string &propertyKey, const std::string &contextName = "", bool fallBackOnDefaultContext = true) const override; std::vector GetPropertyKeys(const std::string &contextName = "", bool includeDefaultContext = false) const override; std::vector GetPropertyContextNames() const override; private: // For future implementation: load slice-by-slice, mark this using these methods void SetSliceIsLoaded(unsigned int index, bool isLoaded); // For future implementation: load slice-by-slice, mark this using these methods bool IsSliceLoaded(unsigned int index) const; // For future implementation: load slice-by-slice, mark this using these methods bool AllSlicesAreLoaded() const; public: /// Describe how the mitk::Image's pixel spacing should be interpreted PixelSpacingInterpretation GetPixelSpacingInterpretation() const; /// Describe the correct x/y pixel spacing of the mitk::Image (which some readers might need to adjust after loading) void GetDesiredMITKImagePixelSpacing(ScalarType& spacingXinMM, ScalarType& spacingYinMM) const; /// Describe the gantry tilt of the acquisition void SetTiltInformation(const GantryTiltInformation& info); /// Describe the gantry tilt of the acquisition const GantryTiltInformation GetTiltInformation() const; /// SOP Class UID of this set of frames void SetSOPClassUID(const std::string& uid); /// SOP Class UID of this set of frames std::string GetSOPClassUID() const; /// SOP Class as human readable name (e.g. "CT Image Storage") std::string GetSOPClassUIDAsName() const; /**Convenience method that returns the property timesteps*/ int GetNumberOfTimeSteps() const; /**return the number of frames that constitute one timestep.*/ int GetNumberOfFramesPerTimeStep() const; void SetTagCache(DICOMTagCache* privateCache); /** Type specifies additional tags of interest. Key is the tag path of interest. * The value is an optional user defined name for the property that should be used to store the tag value(s). * Empty value is default and will imply to use the found DICOMTagPath as property name.*/ typedef std::map AdditionalTagsMapType; /** * \brief Set a list of DICOMTagPaths that specify all DICOM-Tags that will be copied into the property of the mitk::Image. * * This method can be used to specify a list of DICOM-tags that shall be available after the loading. * The value in the tagMap is an optional user defined name for the property key that should be used * when storing the property). Empty value is default and will imply to use the found DICOMTagPath * as property key. * By default the content of the DICOM tags will be stored in a StringLookupTable on the mitk::Image. * This behaviour can be changed by setting a different TagLookupTableToPropertyFunctor via * SetTagLookupTableToPropertyFunctor(). */ void SetAdditionalTagsOfInterest(const AdditionalTagsMapType& tagMap); typedef std::function TagLookupTableToPropertyFunctor; /** * \brief Set a functor that defines how the slice-specific tag-values are stored in a Property. * * This method sets a functor that is given a StringLookupTable that contains the values of one DICOM tag * mapped to the slice index. * The functor is supposed to store these values in an mitk Property. * * By default, the StringLookupTable is stored in a StringLookupTableProperty except if all values are * identical. In this case, the unique value is stored only once in a StringProperty. */ void SetTagLookupTableToPropertyFunctor(TagLookupTableToPropertyFunctor); /// Print information about this image block to given stream void Print(std::ostream& os, bool filenameDetails) const; private: // read values from tag cache std::string GetPixelSpacing() const; std::string GetImagerPixelSpacing() const; Image::Pointer FixupSpacing(Image* mitkImage); Image::Pointer DescribeImageWithProperties(Image* mitkImage); void UpdateImageDescribingProperties() const; static mitk::BaseProperty::Pointer GetPropertyForDICOMValues(const DICOMCachedValueLookupTable& cacheLookupTable); double stringtodouble(const std::string& str) const; DICOMImageFrameList m_ImageFrameList; Image::Pointer m_MitkImage; BoolList m_SliceIsLoaded; ReaderImplementationLevel m_ReaderImplementationLevel; GantryTiltInformation m_TiltInformation; PropertyList::Pointer m_PropertyList; - mitk::WeakPointer m_TagCache; + DICOMSplitReason::Pointer m_SplitReason; + + WeakPointer m_TagCache; mutable bool m_PropertiesOutOfDate; AdditionalTagsMapType m_AdditionalTagMap; std::set m_FoundAdditionalTags; TagLookupTableToPropertyFunctor m_PropertyFunctor; }; } #endif diff --git a/Modules/DICOM/include/mitkDICOMTagBasedSorter.h b/Modules/DICOM/include/mitkDICOMTagBasedSorter.h index b708bf53c1..0a3d0509de 100644 --- a/Modules/DICOM/include/mitkDICOMTagBasedSorter.h +++ b/Modules/DICOM/include/mitkDICOMTagBasedSorter.h @@ -1,202 +1,206 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkDICOMTagBasedSorter_h #define mitkDICOMTagBasedSorter_h #include "mitkDICOMDatasetSorter.h" #include "mitkDICOMSortCriterion.h" namespace mitk { /** \ingroup DICOMModule \brief Sort DICOM datasets based on configurable tags. This class implements sorting of input DICOM datasets into multiple outputs as described in \ref DICOMITKSeriesGDCMReader_LoadingStrategy. The logic of sorting and splitting is most simple and most generic: 1. Datasets will be put into different groups, if they differ in their value of specific tags (defined by AddDistinguishingTag()) - there might be multiple distinguishing tags defined - tag values might be processed before comparison by means of TagValueProcessor (e.g. round to a number of decimal places) 2. Each of the groups will be sorted by comparing their tag values using multiple DICOMSortCriterion - DICOMSortCriterion might evaluate a single tag (e.g. Instance Number) or multiple values (as in SortByImagePositionPatient) - only a single DICOMSortCriterion is defined for DICOMTagBasedSorter, because each DICOMSortCriterion holds a "secondary sort criterion", i.e. an application can define multiple tags for sorting by chaining \link DICOMSortCriterion DICOMSortCriteria \endlink - applications should make sure that sorting is always defined (to avoid problems with standard containers), e.g. by adding a comparison of filenames or instance UIDs as a last sorting fallback. */ class MITKDICOM_EXPORT DICOMTagBasedSorter : public DICOMDatasetSorter { public: /** \brief Processes tag values before they are compared. These classes could do some kind of normalization such as rounding, lower case formatting, etc. */ class MITKDICOM_EXPORT TagValueProcessor { public: /// \brief Implements the "processing". virtual std::string operator()(const std::string&) const = 0; virtual TagValueProcessor* Clone() const = 0; virtual ~TagValueProcessor() {} }; /** \brief Cuts a number after configured number of decimal places. An instance of this class can be used to avoid errors when comparing minimally different image orientations. */ class MITKDICOM_EXPORT CutDecimalPlaces : public TagValueProcessor { public: CutDecimalPlaces(unsigned int precision); CutDecimalPlaces(const CutDecimalPlaces& other); unsigned int GetPrecision() const; std::string operator()(const std::string&) const override; TagValueProcessor* Clone() const override; private: unsigned int m_Precision; }; mitkClassMacro( DICOMTagBasedSorter, DICOMDatasetSorter ); itkNewMacro( DICOMTagBasedSorter ); /** \brief Datasets that differ in given tag's value will be sorted into separate outputs. */ void AddDistinguishingTag( const DICOMTag&, TagValueProcessor* tagValueProcessor = nullptr ); DICOMTagList GetDistinguishingTags() const; const TagValueProcessor* GetTagValueProcessorForDistinguishingTag(const DICOMTag&) const; /** \brief Define the sorting criterion (which holds seconardy criteria) */ void SetSortCriterion( DICOMSortCriterion::ConstPointer criterion ); DICOMSortCriterion::ConstPointer GetSortCriterion() const; /** \brief A list of all the tags needed for processing (facilitates scanning). */ DICOMTagList GetTagsOfInterest() override; /** \brief Whether or not groups should be checked for consecutive tag values. When this flag is set (default in constructor=off), the sorter will not only sort in a way that the values of a configured tag are ascending BUT in addition the sorter will enforce a constant numerical distance between values. Having this flag is useful for handling of series with missing slices, e.g. Instance Numbers 1 2 3 5 6 7 8. With the flag set to true, the sorter would split this group into two, because the initial distance of 1 is not kept between Instance Numbers 3 and 5. A special case of this behavior can be configured by SetExpectDistanceOne(). When this additional flag is set to true, the sorter will expect distance 1 exactly. This can help if the second slice is missing already. Without this additional flag, we would "learn" about a wrong distance of 2 (or similar) and then sort completely wrong. */ void SetStrictSorting(bool strict); bool GetStrictSorting() const; /** \brief Flag for a special case in "strict sorting". Please see documentation of SetStrictSorting(). \sa SetStrictSorting */ void SetExpectDistanceOne(bool strict); bool GetExpectDistanceOne() const; /** \brief Actually sort as described in the Detailed Description. */ void Sort() override; /** \brief Print configuration details into given stream. */ void PrintConfiguration(std::ostream& os, const std::string& indent = "") const override; bool operator==(const DICOMDatasetSorter& other) const override; static bool GetDefaultStrictSorting() { return m_DefaultStrictSorting; } static bool GetDefaultExpectDistanceOne() { return m_DefaultExpectDistanceOne; } protected: /** \brief Helper struct to feed into std::sort, configured via DICOMSortCriterion. */ struct ParameterizedDatasetSort { ParameterizedDatasetSort(DICOMSortCriterion::ConstPointer); bool operator() (const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right); DICOMSortCriterion::ConstPointer m_SortCriterion; }; DICOMTagBasedSorter(); ~DICOMTagBasedSorter() override; DICOMTagBasedSorter(const DICOMTagBasedSorter& other); DICOMTagBasedSorter& operator=(const DICOMTagBasedSorter& other); /** \brief Helper for SplitInputGroups(). */ std::string BuildGroupID( DICOMDatasetAccess* dataset ); - typedef std::map GroupIDToListType; + using GroupIDToListType = std::map; + + using SplitReasonListType = std::map; /** \brief Implements the "distiguishing tags". To sort datasets into different groups, a long string will be built for each dataset. The string concatenates all tags and their respective values. Datasets that match in all values will end up with the same string. + @param splitReasons Reference to the split reasons vector. It will be also updated by the method to reflect the reasons for the returned groups. */ - GroupIDToListType SplitInputGroups(); + GroupIDToListType SplitInputGroups(SplitReasonListType& splitReasons); /** \brief Implements the sorting step. Relatively simple implementation thanks to std::sort and a parameterization via DICOMSortCriterion. + @param splitReasons Reference to the split reasons vector. It will be also updated by the method to reflect the reasons for the returned groups. */ - GroupIDToListType& SortGroups(GroupIDToListType& groups); + GroupIDToListType& SortGroups(GroupIDToListType& groups, SplitReasonListType& splitReasons); DICOMTagList m_DistinguishingTags; typedef std::map TagValueProcessorMap; TagValueProcessorMap m_TagValueProcessor; DICOMSortCriterion::ConstPointer m_SortCriterion; bool m_StrictSorting; bool m_ExpectDistanceOne; const static bool m_DefaultStrictSorting = false; const static bool m_DefaultExpectDistanceOne = false; }; } #endif diff --git a/Modules/DICOM/include/mitkEquiDistantBlocksSorter.h b/Modules/DICOM/include/mitkEquiDistantBlocksSorter.h index 26e1b0da35..f9636861e7 100644 --- a/Modules/DICOM/include/mitkEquiDistantBlocksSorter.h +++ b/Modules/DICOM/include/mitkEquiDistantBlocksSorter.h @@ -1,211 +1,222 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #ifndef mitkEquiDistantBlocksSorter_h #define mitkEquiDistantBlocksSorter_h #include "mitkDICOMDatasetSorter.h" #include "mitkDICOMSortCriterion.h" #include "mitkGantryTiltInformation.h" #include "mitkVector.h" namespace mitk { /** \ingroup DICOMModule \brief Split inputs into blocks of equidistant slices (for use in DICOMITKSeriesGDCMReader). Since inter-slice distance is not recorded in DICOM tags, we must ensure that blocks are made up of slices that have equal distances between neighboring slices. This is especially necessary because itk::ImageSeriesReader is later used for the actual loading, and this class expects (and does nocht verify) equal inter-slice distance (see \ref DICOMITKSeriesGDCMReader_ForcedConfiguration). To achieve such grouping, the inter-slice distance is calculated from the first two different slice positions of a block. Following slices are added to a block as long as they can be added by adding the calculated inter-slice distance to the last slice of the block. Slices that do not fit into the expected distance pattern, are set aside for further analysis. This grouping is done until each file has been assigned to a group. Slices that share a position in space are also sorted into separate blocks during this step. So the result of this step is a set of blocks that contain only slices with equal z spacing and unique slices at each position. During sorting, the origins (documented in tag image position patient) are compared against expected origins (from former origin plus moving direction). As there will be minor differences in numbers (from both calculations and imprecise tag values), we must be a bit tolerant here. The default behavior is to expect that an origin is not further away from the expected position than 30% of the inter-slice distance. To support a legacy behavior of a former loader (DicomSeriesReader), this default can be restricted to a constant number of millimeters by calling SetToleratedOriginOffset(mm). + REMARK: The EquiDistantBlocksSorter assumes that the order of the provided input + is sorted by image position like it is preferred by the reader and does not sort + it again. This assumption can lead to splittings even for complete volumes if + the input is not sorted by image position can (Reason: gaps will be detected + because the next slice will not have the assumed distance and will be sorted out.) + Detailed implementation in AnalyzeFileForITKImageSeriesReaderSpacingAssumption(). */ class MITKDICOM_EXPORT EquiDistantBlocksSorter : public DICOMDatasetSorter { public: mitkClassMacro( EquiDistantBlocksSorter, DICOMDatasetSorter ); itkNewMacro( EquiDistantBlocksSorter ); DICOMTagList GetTagsOfInterest() override; /** \brief Delegates work to AnalyzeFileForITKImageSeriesReaderSpacingAssumption(). AnalyzeFileForITKImageSeriesReaderSpacingAssumption() is called until it does not create multiple blocks anymore. */ void Sort() override; /** \brief Whether or not to accept images from a tilted acquisition in a single output group. */ void SetAcceptTilt(bool accept); bool GetAcceptTilt() const; /** \brief See class description and SetToleratedOriginOffset(). */ void SetToleratedOriginOffsetToAdaptive(double fractionOfInterSliceDistanct = 0.3); /** \brief See class description and SetToleratedOriginOffsetToAdaptive(). Default value of 0.005 is calculated so that we get a maximum of 1/10mm error when having a measurement crosses 20 slices in z direction (too strict? we don't know better..). */ void SetToleratedOriginOffset(double millimeters = 0.005); double GetToleratedOriginOffset() const; bool IsToleratedOriginOffsetAbsolute() const; void SetAcceptTwoSlicesGroups(bool accept); bool GetAcceptTwoSlicesGroups() const; void PrintConfiguration(std::ostream& os, const std::string& indent = "") const override; bool operator==(const DICOMDatasetSorter& other) const override; protected: /** \brief Return type of AnalyzeFileForITKImageSeriesReaderSpacingAssumption(). Class contains the grouping result of method AnalyzeFileForITKImageSeriesReaderSpacingAssumption(), which takes as input a number of images, which are all equally oriented and spatially sorted along their normal direction. The result contains of two blocks: a first one is the grouping result, all of those images can be loaded into one image block because they have an equal origin-to-origin distance without any gaps in-between. */ class SliceGroupingAnalysisResult { public: SliceGroupingAnalysisResult(); /** \brief Grouping result, all same origin-to-origin distance w/o gaps. */ - DICOMDatasetList GetBlockDatasets(); + const DICOMDatasetList& GetBlockDatasets() const; void SetFirstFilenameOfBlock(const std::string& filename); std::string GetFirstFilenameOfBlock() const; void SetLastFilenameOfBlock(const std::string& filename); std::string GetLastFilenameOfBlock() const; /** \brief Remaining files, which could not be grouped. */ - DICOMDatasetList GetUnsortedDatasets(); + const DICOMDatasetList& GetUnsortedDatasets() const; + + const DICOMSplitReason* GetSplitReason() const; + DICOMSplitReason* GetSplitReason(); /** \brief Whether or not the grouped result contain a gantry tilt. */ bool ContainsGantryTilt(); /** \brief Detailed description of gantry tilt. */ const GantryTiltInformation& GetTiltInfo() const; /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. */ void AddFileToSortedBlock(DICOMDatasetAccess* dataset); /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. */ void AddFileToUnsortedBlock(DICOMDatasetAccess* dataset); void AddFilesToUnsortedBlock(const DICOMDatasetList& datasets); /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. \todo Could make sense to enhance this with an instance of GantryTiltInformation to store the whole result! */ void FlagGantryTilt(const GantryTiltInformation& tiltInfo); /** \brief Only meaningful for use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption. */ void UndoPrematureGrouping(); protected: DICOMDatasetList m_GroupedFiles; DICOMDatasetList m_UnsortedFiles; + DICOMSplitReason::Pointer m_SplitReason; + GantryTiltInformation m_TiltInfo; std::string m_FirstFilenameOfBlock; std::string m_LastFilenameOfBlock; }; /** \brief Ensure an equal z-spacing for a group of files. Takes as input a number of images, which are all equally oriented and spatially sorted along their normal direction. Internally used by GetSeries. Returns two lists: the first one contains slices of equal inter-slice spacing. The second list contains remaining files, which need to be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption again. Relevant code that is matched here is in itkImageSeriesReader.txx (ImageSeriesReader::GenerateOutputInformation(void)), lines 176 to 245 (as of ITK 3.20) */ - SliceGroupingAnalysisResult + std::shared_ptr AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const DICOMDatasetList& files, bool groupsOfSimilarImages); /** \brief Safely convert const char* to std::string. */ std::string ConstCharStarToString(const char* s); EquiDistantBlocksSorter(); ~EquiDistantBlocksSorter() override; EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other); EquiDistantBlocksSorter& operator=(const EquiDistantBlocksSorter& other); bool m_AcceptTilt; - typedef std::vector ResultsList; + typedef std::vector > ResultsList; ResultsList m_SliceGroupingResults; double m_ToleratedOriginOffset; bool m_ToleratedOriginOffsetIsAbsolute; bool m_AcceptTwoSlicesGroups; }; } #endif diff --git a/Modules/DICOM/src/mitkDICOMDatasetSorter.cpp b/Modules/DICOM/src/mitkDICOMDatasetSorter.cpp index d6af270fe5..e9c3a8c754 100644 --- a/Modules/DICOM/src/mitkDICOMDatasetSorter.cpp +++ b/Modules/DICOM/src/mitkDICOMDatasetSorter.cpp @@ -1,117 +1,137 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkDICOMDatasetSorter.h" mitk::DICOMDatasetSorter ::DICOMDatasetSorter() :itk::LightObject() { } mitk::DICOMDatasetSorter ::~DICOMDatasetSorter() { } mitk::DICOMDatasetSorter ::DICOMDatasetSorter(const DICOMDatasetSorter& other ) :itk::LightObject() ,m_Outputs( other.m_Outputs ) { } mitk::DICOMDatasetSorter& mitk::DICOMDatasetSorter ::operator=(const DICOMDatasetSorter& other) { if (this != &other) { m_Input = other.m_Input; m_Outputs = other.m_Outputs; } return *this; } void mitk::DICOMDatasetSorter ::SetInput(DICOMDatasetList datasets) { m_Input = datasets; } const mitk::DICOMDatasetList& mitk::DICOMDatasetSorter ::GetInput() const { return m_Input; } unsigned int mitk::DICOMDatasetSorter ::GetNumberOfOutputs() const { return m_Outputs.size(); } void mitk::DICOMDatasetSorter ::ClearOutputs() { m_Outputs.clear(); + m_SplitReasons.clear(); } void mitk::DICOMDatasetSorter ::SetNumberOfOutputs(unsigned int numberOfOutputs) { m_Outputs.resize(numberOfOutputs); + m_SplitReasons.resize(numberOfOutputs); } void mitk::DICOMDatasetSorter -::SetOutput(unsigned int index, const DICOMDatasetList& output) +::SetOutput(unsigned int index, const DICOMDatasetList& output, const DICOMSplitReason* splitReason) { if (index < m_Outputs.size()) { m_Outputs[index] = output; + if (nullptr == splitReason) + m_SplitReasons[index] = DICOMSplitReason::New(); + else + m_SplitReasons[index] = splitReason->Clone(); } else { std::stringstream ss; - ss << "Index " << index << " out of range (" << m_Outputs.size() << " indices reserved)"; + ss << "Cannot get output. Index " << index << " out of range (" << m_Outputs.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } const mitk::DICOMDatasetList& mitk::DICOMDatasetSorter ::GetOutput(unsigned int index) const { return const_cast(this)->GetOutput(index); } mitk::DICOMDatasetList& mitk::DICOMDatasetSorter ::GetOutput(unsigned int index) { if (index < m_Outputs.size()) { return m_Outputs[index]; } else { std::stringstream ss; - ss << "Index " << index << " out of range (" << m_Outputs.size() << " indices reserved)"; + ss << "Cannot get output. Index " << index << " out of range (" << m_Outputs.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } + +const mitk::DICOMSplitReason* +mitk::DICOMDatasetSorter +::GetSplitReason(unsigned int index) const +{ + if (index >= m_Outputs.size()) + { + std::stringstream ss; + ss << "Cannot get split reason. Index " << index << " out of range (" << m_Outputs.size() << " indices reserved)"; + throw std::invalid_argument(ss.str()); + } + + return m_SplitReasons[index]; +} diff --git a/Modules/DICOM/src/mitkDICOMITKSeriesGDCMReader.cpp b/Modules/DICOM/src/mitkDICOMITKSeriesGDCMReader.cpp index c7c4c24179..c2ab71e52b 100644 --- a/Modules/DICOM/src/mitkDICOMITKSeriesGDCMReader.cpp +++ b/Modules/DICOM/src/mitkDICOMITKSeriesGDCMReader.cpp @@ -1,630 +1,635 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ //#define MBILOG_ENABLE_DEBUG #define ENABLE_TIMING #include #include #include "mitkDICOMITKSeriesGDCMReader.h" #include "mitkITKDICOMSeriesReaderHelper.h" #include "mitkGantryTiltInformation.h" #include "mitkDICOMTagBasedSorter.h" #include "mitkDICOMGDCMTagScanner.h" std::mutex mitk::DICOMITKSeriesGDCMReader::s_LocaleMutex; mitk::DICOMITKSeriesGDCMReader::DICOMITKSeriesGDCMReader( unsigned int decimalPlacesForOrientation, bool simpleVolumeImport ) : DICOMFileReader() , m_FixTiltByShearing(m_DefaultFixTiltByShearing) , m_SimpleVolumeReading( simpleVolumeImport ) , m_DecimalPlacesForOrientation( decimalPlacesForOrientation ) , m_ExternalCache(false) { this->EnsureMandatorySortersArePresent( decimalPlacesForOrientation, simpleVolumeImport ); } mitk::DICOMITKSeriesGDCMReader::DICOMITKSeriesGDCMReader( const DICOMITKSeriesGDCMReader& other ) : DICOMFileReader( other ) , m_FixTiltByShearing( other.m_FixTiltByShearing) , m_SimpleVolumeReading( other.m_SimpleVolumeReading) , m_SortingResultInProgress( other.m_SortingResultInProgress ) , m_Sorter( other.m_Sorter ) , m_EquiDistantBlocksSorter( other.m_EquiDistantBlocksSorter->Clone() ) , m_NormalDirectionConsistencySorter( other.m_NormalDirectionConsistencySorter->Clone() ) , m_ReplacedCLocales( other.m_ReplacedCLocales ) , m_ReplacedCinLocales( other.m_ReplacedCinLocales ) , m_DecimalPlacesForOrientation( other.m_DecimalPlacesForOrientation ) , m_TagCache( other.m_TagCache ) , m_ExternalCache(other.m_ExternalCache) { } mitk::DICOMITKSeriesGDCMReader::~DICOMITKSeriesGDCMReader() { } mitk::DICOMITKSeriesGDCMReader& mitk::DICOMITKSeriesGDCMReader:: operator=( const DICOMITKSeriesGDCMReader& other ) { if ( this != &other ) { DICOMFileReader::operator =( other ); this->m_FixTiltByShearing = other.m_FixTiltByShearing; this->m_SimpleVolumeReading = other.m_SimpleVolumeReading; this->m_SortingResultInProgress = other.m_SortingResultInProgress; this->m_Sorter = other.m_Sorter; // TODO should clone the list items this->m_EquiDistantBlocksSorter = other.m_EquiDistantBlocksSorter->Clone(); this->m_NormalDirectionConsistencySorter = other.m_NormalDirectionConsistencySorter->Clone(); this->m_ReplacedCLocales = other.m_ReplacedCLocales; this->m_ReplacedCinLocales = other.m_ReplacedCinLocales; this->m_DecimalPlacesForOrientation = other.m_DecimalPlacesForOrientation; this->m_TagCache = other.m_TagCache; } return *this; } bool mitk::DICOMITKSeriesGDCMReader::operator==( const DICOMFileReader& other ) const { if ( const auto* otherSelf = dynamic_cast( &other ) ) { if ( this->m_FixTiltByShearing == otherSelf->m_FixTiltByShearing && *( this->m_EquiDistantBlocksSorter ) == *( otherSelf->m_EquiDistantBlocksSorter ) && ( fabs( this->m_DecimalPlacesForOrientation - otherSelf->m_DecimalPlacesForOrientation ) < eps ) ) { // test sorters for equality if ( this->m_Sorter.size() != otherSelf->m_Sorter.size() ) return false; auto mySorterIter = this->m_Sorter.cbegin(); auto oSorterIter = otherSelf->m_Sorter.cbegin(); for ( ; mySorterIter != this->m_Sorter.cend() && oSorterIter != otherSelf->m_Sorter.cend(); ++mySorterIter, ++oSorterIter ) { if ( !( **mySorterIter == **oSorterIter ) ) return false; // this sorter differs } // nothing differs ==> all is equal return true; } else { return false; } } else { return false; } } void mitk::DICOMITKSeriesGDCMReader::SetFixTiltByShearing( bool on ) { this->Modified(); m_FixTiltByShearing = on; } bool mitk::DICOMITKSeriesGDCMReader::GetFixTiltByShearing() const { return m_FixTiltByShearing; } void mitk::DICOMITKSeriesGDCMReader::SetAcceptTwoSlicesGroups( bool accept ) const { this->Modified(); m_EquiDistantBlocksSorter->SetAcceptTwoSlicesGroups( accept ); } bool mitk::DICOMITKSeriesGDCMReader::GetAcceptTwoSlicesGroups() const { return m_EquiDistantBlocksSorter->GetAcceptTwoSlicesGroups(); } void mitk::DICOMITKSeriesGDCMReader::InternalPrintConfiguration( std::ostream& os ) const { unsigned int sortIndex( 1 ); for ( auto sorterIter = m_Sorter.cbegin(); sorterIter != m_Sorter.cend(); ++sortIndex, ++sorterIter ) { os << "Sorting step " << sortIndex << ":" << std::endl; ( *sorterIter )->PrintConfiguration( os, " " ); } os << "Sorting step " << sortIndex << ":" << std::endl; m_EquiDistantBlocksSorter->PrintConfiguration( os, " " ); } std::string mitk::DICOMITKSeriesGDCMReader::GetActiveLocale() { return setlocale( LC_NUMERIC, nullptr ); } void mitk::DICOMITKSeriesGDCMReader::PushLocale() const { s_LocaleMutex.lock(); std::string currentCLocale = setlocale( LC_NUMERIC, nullptr ); m_ReplacedCLocales.push( currentCLocale ); setlocale( LC_NUMERIC, "C" ); std::locale currentCinLocale( std::cin.getloc() ); m_ReplacedCinLocales.push( currentCinLocale ); std::locale l( "C" ); std::cin.imbue( l ); s_LocaleMutex.unlock(); } void mitk::DICOMITKSeriesGDCMReader::PopLocale() const { s_LocaleMutex.lock(); if ( !m_ReplacedCLocales.empty() ) { setlocale( LC_NUMERIC, m_ReplacedCLocales.top().c_str() ); m_ReplacedCLocales.pop(); } else { MITK_WARN << "Mismatched PopLocale on DICOMITKSeriesGDCMReader."; } if ( !m_ReplacedCinLocales.empty() ) { std::cin.imbue( m_ReplacedCinLocales.top() ); m_ReplacedCinLocales.pop(); } else { MITK_WARN << "Mismatched PopLocale on DICOMITKSeriesGDCMReader."; } s_LocaleMutex.unlock(); } mitk::DICOMITKSeriesGDCMReader::SortingBlockList mitk::DICOMITKSeriesGDCMReader::Condense3DBlocks( SortingBlockList& input ) { return input; // to be implemented differently by sub-classes } #if defined( MBILOG_ENABLE_DEBUG ) || defined( ENABLE_TIMING ) #define timeStart( part ) timer.Start( part ); #define timeStop( part ) timer.Stop( part ); #else #define timeStart( part ) #define timeStop( part ) #endif void mitk::DICOMITKSeriesGDCMReader::AnalyzeInputFiles() { itk::TimeProbesCollectorBase timer; timeStart( "Reset" ); this->ClearOutputs(); timeStop( "Reset" ); // prepare initial sorting (== list of input files) const StringList inputFilenames = this->GetInputFiles(); timeStart( "Check input for DCM" ); if ( inputFilenames.empty() || !this->CanHandleFile( inputFilenames.front() ) // first || !this->CanHandleFile( inputFilenames.back() ) // last || !this->CanHandleFile( inputFilenames[inputFilenames.size() / 2] ) // roughly central file ) { // TODO a read-as-many-as-possible fallback could be implemented here MITK_DEBUG << "Reader unable to process files.."; return; } timeStop( "Check input for DCM" ); // scan files for sorting-relevant tags if ( m_TagCache.IsNull() || ( m_TagCache->GetMTime()GetMTime() && !m_ExternalCache )) { timeStart( "Tag scanning" ); DICOMGDCMTagScanner::Pointer filescanner = DICOMGDCMTagScanner::New(); filescanner->SetInputFiles( inputFilenames ); filescanner->AddTagPaths( this->GetTagsOfInterest() ); PushLocale(); filescanner->Scan(); PopLocale(); m_TagCache = filescanner->GetScanCache(); // keep alive and make accessible to sub-classes timeStop("Tag scanning"); } else { // ensure that the tag cache contains our required tags AND files and has scanned! } m_SortingResultInProgress.clear(); - m_SortingResultInProgress.push_back(m_TagCache->GetFrameInfoList()); + m_SortingResultInProgress.push_back(std::make_pair(m_TagCache->GetFrameInfoList(), DICOMSplitReason::New())); // sort and split blocks as configured timeStart( "Sorting frames" ); unsigned int sorterIndex = 0; for ( auto sorterIter = m_Sorter.cbegin(); sorterIter != m_Sorter.cend(); ++sorterIndex, ++sorterIter ) { std::stringstream ss; ss << "Sorting step " << sorterIndex; timeStart( ss.str().c_str() ); m_SortingResultInProgress = this->InternalExecuteSortingStep( sorterIndex, *sorterIter, m_SortingResultInProgress ); timeStop( ss.str().c_str() ); } if ( !m_SimpleVolumeReading ) { // a last extra-sorting step: ensure equidistant slices timeStart( "EquiDistantBlocksSorter" ); m_SortingResultInProgress = this->InternalExecuteSortingStep( sorterIndex++, m_EquiDistantBlocksSorter.GetPointer(), m_SortingResultInProgress ); timeStop( "EquiDistantBlocksSorter" ); } timeStop( "Sorting frames" ); timeStart( "Condensing 3D blocks" ); m_SortingResultInProgress = this->Condense3DBlocks( m_SortingResultInProgress ); timeStop( "Condensing 3D blocks" ); // provide final result as output timeStart( "Output" ); unsigned int o = this->GetNumberOfOutputs(); this->SetNumberOfOutputs( o + m_SortingResultInProgress.size() ); // Condense3DBlocks may already have added outputs! for ( auto blockIter = m_SortingResultInProgress.cbegin(); blockIter != m_SortingResultInProgress.cend(); ++o, ++blockIter ) { - const DICOMDatasetAccessingImageFrameList& gdcmFrameInfoList = *blockIter; + const auto& gdcmFrameInfoList = blockIter->first; + auto& splitReason = blockIter->second; + assert( !gdcmFrameInfoList.empty() ); // reverse frames if necessary // update tilt information from absolute last sorting const DICOMDatasetList datasetList = ConvertToDICOMDatasetList( gdcmFrameInfoList ); m_NormalDirectionConsistencySorter->SetInput( datasetList ); m_NormalDirectionConsistencySorter->Sort(); const DICOMDatasetAccessingImageFrameList sortedGdcmInfoFrameList = ConvertToDICOMDatasetAccessingImageFrameList( m_NormalDirectionConsistencySorter->GetOutput( 0 ) ); const GantryTiltInformation& tiltInfo = m_NormalDirectionConsistencySorter->GetTiltInformation(); // set frame list for current block const DICOMImageFrameList frameList = ConvertToDICOMImageFrameList( sortedGdcmInfoFrameList ); assert( !frameList.empty() ); DICOMImageBlockDescriptor block; block.SetTagCache( this->GetTagCache() ); // important: this must be before SetImageFrameList(), because // SetImageFrameList will trigger reading of lots of interesting // tags! block.SetAdditionalTagsOfInterest( GetAdditionalTagsOfInterest() ); block.SetTagLookupTableToPropertyFunctor( GetTagLookupTableToPropertyFunctor() ); block.SetImageFrameList( frameList ); block.SetTiltInformation( tiltInfo ); + block.SetSplitReason(splitReason); block.SetReaderImplementationLevel( this->GetReaderImplementationLevel( block.GetSOPClassUID() ) ); this->SetOutput( o, block ); } timeStop( "Output" ); #if defined( MBILOG_ENABLE_DEBUG ) || defined( ENABLE_TIMING ) std::cout << "---------------------------------------------------------------" << std::endl; timer.Report( std::cout ); std::cout << "---------------------------------------------------------------" << std::endl; #endif } mitk::DICOMITKSeriesGDCMReader::SortingBlockList mitk::DICOMITKSeriesGDCMReader::InternalExecuteSortingStep( unsigned int sortingStepIndex, const DICOMDatasetSorter::Pointer& sorter, const SortingBlockList& input ) { SortingBlockList nextStepSorting; // we should not modify our input list while processing it std::stringstream ss; ss << "Sorting step " << sortingStepIndex << " '"; #if defined( MBILOG_ENABLE_DEBUG ) sorter->PrintConfiguration( ss ); #endif ss << "'"; nextStepSorting.clear(); MITK_DEBUG << "================================================================================"; MITK_DEBUG << "DICOMITKSeriesGDCMReader: " << ss.str() << ": " << input.size() << " groups input"; #if defined( MBILOG_ENABLE_DEBUG ) unsigned int groupIndex = 0; #endif for ( auto blockIter = input.cbegin(); blockIter != input.cend(); #if defined( MBILOG_ENABLE_DEBUG ) ++groupIndex, #endif ++blockIter ) { - const DICOMDatasetAccessingImageFrameList& gdcmInfoFrameList = *blockIter; + const auto& gdcmInfoFrameList = blockIter->first; + const auto& inputSplitReason = blockIter->second; const DICOMDatasetList datasetList = ConvertToDICOMDatasetList( gdcmInfoFrameList ); #if defined( MBILOG_ENABLE_DEBUG ) MITK_DEBUG << "--------------------------------------------------------------------------------"; MITK_DEBUG << "DICOMITKSeriesGDCMReader: " << ss.str() << ", dataset group " << groupIndex << " (" << datasetList.size() << " datasets): "; for ( auto oi = datasetList.cbegin(); oi != datasetList.cend(); ++oi ) { MITK_DEBUG << " INPUT : " << ( *oi )->GetFilenameIfAvailable(); } #endif sorter->SetInput( datasetList ); sorter->Sort(); unsigned int numberOfResultingBlocks = sorter->GetNumberOfOutputs(); for ( unsigned int b = 0; b < numberOfResultingBlocks; ++b ) { const DICOMDatasetList blockResult = sorter->GetOutput( b ); for ( auto oi = blockResult.cbegin(); oi != blockResult.cend(); ++oi ) { MITK_DEBUG << " OUTPUT(" << b << ") :" << ( *oi )->GetFilenameIfAvailable(); } DICOMDatasetAccessingImageFrameList sortedGdcmInfoFrameList = ConvertToDICOMDatasetAccessingImageFrameList( blockResult ); - nextStepSorting.push_back( sortedGdcmInfoFrameList ); + + nextStepSorting.push_back( std::make_pair(sortedGdcmInfoFrameList, inputSplitReason->ExtendReason(sorter->GetSplitReason(b))) ); } } return nextStepSorting; } mitk::ReaderImplementationLevel mitk::DICOMITKSeriesGDCMReader::GetReaderImplementationLevel( const std::string sopClassUID ) { if ( sopClassUID.empty() ) { return SOPClassUnknown; } gdcm::UIDs uidKnowledge; uidKnowledge.SetFromUID( sopClassUID.c_str() ); gdcm::UIDs::TSName gdcmType = static_cast((gdcm::UIDs::TSType)uidKnowledge); switch ( gdcmType ) { case gdcm::UIDs::CTImageStorage: case gdcm::UIDs::MRImageStorage: case gdcm::UIDs::PositronEmissionTomographyImageStorage: case gdcm::UIDs::ComputedRadiographyImageStorage: case gdcm::UIDs::DigitalXRayImageStorageForPresentation: case gdcm::UIDs::DigitalXRayImageStorageForProcessing: return SOPClassSupported; case gdcm::UIDs::NuclearMedicineImageStorage: return SOPClassPartlySupported; case gdcm::UIDs::SecondaryCaptureImageStorage: return SOPClassImplemented; default: return SOPClassUnsupported; } } // void AllocateOutputImages(); bool mitk::DICOMITKSeriesGDCMReader::LoadImages() { bool success = true; unsigned int numberOfOutputs = this->GetNumberOfOutputs(); for ( unsigned int o = 0; o < numberOfOutputs; ++o ) { success &= this->LoadMitkImageForOutput( o ); } return success; } bool mitk::DICOMITKSeriesGDCMReader::LoadMitkImageForImageBlockDescriptor( DICOMImageBlockDescriptor& block ) const { PushLocale(); const DICOMImageFrameList& frames = block.GetImageFrameList(); const GantryTiltInformation tiltInfo = block.GetTiltInformation(); bool hasTilt = tiltInfo.IsRegularGantryTilt(); ITKDICOMSeriesReaderHelper::StringContainer filenames; filenames.reserve( frames.size() ); for ( auto frameIter = frames.cbegin(); frameIter != frames.cend(); ++frameIter ) { filenames.push_back( ( *frameIter )->Filename ); } mitk::ITKDICOMSeriesReaderHelper helper; bool success( true ); try { mitk::Image::Pointer mitkImage = helper.Load( filenames, m_FixTiltByShearing && hasTilt, tiltInfo ); block.SetMitkImage( mitkImage ); } catch ( const std::exception& e ) { success = false; MITK_ERROR << "Exception during image loading: " << e.what(); } PopLocale(); return success; } bool mitk::DICOMITKSeriesGDCMReader::LoadMitkImageForOutput( unsigned int o ) { DICOMImageBlockDescriptor& block = this->InternalGetOutput( o ); return this->LoadMitkImageForImageBlockDescriptor( block ); } bool mitk::DICOMITKSeriesGDCMReader::CanHandleFile( const std::string& filename ) { return ITKDICOMSeriesReaderHelper::CanHandleFile( filename ); } void mitk::DICOMITKSeriesGDCMReader::AddSortingElement( DICOMDatasetSorter* sorter, bool atFront ) { assert( sorter ); if ( atFront ) { m_Sorter.push_front( sorter ); } else { m_Sorter.push_back( sorter ); } this->Modified(); } mitk::DICOMITKSeriesGDCMReader::ConstSorterList mitk::DICOMITKSeriesGDCMReader::GetFreelyConfiguredSortingElements() const { std::list result; unsigned int sortIndex( 0 ); for ( auto sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sortIndex, ++sorterIter ) { if ( sortIndex > 0 ) // ignore first element (see EnsureMandatorySortersArePresent) { result.push_back( ( *sorterIter ).GetPointer() ); } } return result; } void mitk::DICOMITKSeriesGDCMReader::EnsureMandatorySortersArePresent( unsigned int decimalPlacesForOrientation, bool simpleVolumeImport ) { DICOMTagBasedSorter::Pointer splitter = DICOMTagBasedSorter::New(); splitter->AddDistinguishingTag( DICOMTag(0x0028, 0x0010) ); // Number of Rows splitter->AddDistinguishingTag( DICOMTag(0x0028, 0x0011) ); // Number of Columns splitter->AddDistinguishingTag( DICOMTag(0x0028, 0x0030) ); // Pixel Spacing splitter->AddDistinguishingTag( DICOMTag(0x0018, 0x1164) ); // Imager Pixel Spacing splitter->AddDistinguishingTag( DICOMTag(0x0020, 0x0037), new mitk::DICOMTagBasedSorter::CutDecimalPlaces(decimalPlacesForOrientation) ); // Image Orientation (Patient) splitter->AddDistinguishingTag( DICOMTag(0x0018, 0x0050) ); // Slice Thickness if ( simpleVolumeImport ) { MITK_DEBUG << "Simple volume reading: ignoring number of frames"; } else { splitter->AddDistinguishingTag( DICOMTag(0x0028, 0x0008) ); // Number of Frames } this->AddSortingElement( splitter, true ); // true = at front if ( m_EquiDistantBlocksSorter.IsNull() ) { m_EquiDistantBlocksSorter = mitk::EquiDistantBlocksSorter::New(); } m_EquiDistantBlocksSorter->SetAcceptTilt( m_FixTiltByShearing ); if ( m_NormalDirectionConsistencySorter.IsNull() ) { m_NormalDirectionConsistencySorter = mitk::NormalDirectionConsistencySorter::New(); } } void mitk::DICOMITKSeriesGDCMReader::SetToleratedOriginOffsetToAdaptive( double fractionOfInterSliceDistance ) const { assert( m_EquiDistantBlocksSorter.IsNotNull() ); m_EquiDistantBlocksSorter->SetToleratedOriginOffsetToAdaptive( fractionOfInterSliceDistance ); this->Modified(); } void mitk::DICOMITKSeriesGDCMReader::SetToleratedOriginOffset( double millimeters ) const { assert( m_EquiDistantBlocksSorter.IsNotNull() ); m_EquiDistantBlocksSorter->SetToleratedOriginOffset( millimeters ); this->Modified(); } double mitk::DICOMITKSeriesGDCMReader::GetToleratedOriginError() const { assert( m_EquiDistantBlocksSorter.IsNotNull() ); return m_EquiDistantBlocksSorter->GetToleratedOriginOffset(); } bool mitk::DICOMITKSeriesGDCMReader::IsToleratedOriginOffsetAbsolute() const { assert( m_EquiDistantBlocksSorter.IsNotNull() ); return m_EquiDistantBlocksSorter->IsToleratedOriginOffsetAbsolute(); } double mitk::DICOMITKSeriesGDCMReader::GetDecimalPlacesForOrientation() const { return m_DecimalPlacesForOrientation; } mitk::DICOMTagCache::Pointer mitk::DICOMITKSeriesGDCMReader::GetTagCache() const { return m_TagCache; } void mitk::DICOMITKSeriesGDCMReader::SetTagCache( const DICOMTagCache::Pointer& tagCache ) { m_TagCache = tagCache; m_ExternalCache = tagCache.IsNotNull(); } mitk::DICOMTagPathList mitk::DICOMITKSeriesGDCMReader::GetTagsOfInterest() const { DICOMTagPathList completeList; // check all configured sorters for ( auto sorterIter = m_Sorter.cbegin(); sorterIter != m_Sorter.cend(); ++sorterIter ) { assert( sorterIter->IsNotNull() ); const DICOMTagList tags = ( *sorterIter )->GetTagsOfInterest(); completeList.insert( completeList.end(), tags.cbegin(), tags.cend() ); } // check our own forced sorters DICOMTagList tags = m_EquiDistantBlocksSorter->GetTagsOfInterest(); completeList.insert( completeList.end(), tags.cbegin(), tags.cend() ); tags = m_NormalDirectionConsistencySorter->GetTagsOfInterest(); completeList.insert( completeList.end(), tags.cbegin(), tags.cend() ); // add the tags for DICOMImageBlockDescriptor tags = DICOMImageBlockDescriptor::GetTagsOfInterest(); completeList.insert( completeList.end(), tags.cbegin(), tags.cend() ); const AdditionalTagsMapType tagList = GetAdditionalTagsOfInterest(); for ( auto iter = tagList.cbegin(); iter != tagList.cend(); ++iter ) { completeList.push_back( iter->first ) ; } return completeList; } diff --git a/Modules/DICOM/src/mitkDICOMImageBlockDescriptor.cpp b/Modules/DICOM/src/mitkDICOMImageBlockDescriptor.cpp index 334d8da6fb..be10aa8af6 100644 --- a/Modules/DICOM/src/mitkDICOMImageBlockDescriptor.cpp +++ b/Modules/DICOM/src/mitkDICOMImageBlockDescriptor.cpp @@ -1,912 +1,934 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkDICOMImageBlockDescriptor.h" #include "mitkStringProperty.h" #include "mitkLevelWindowProperty.h" #include "mitkPropertyKeyPath.h" #include "mitkDICOMIOMetaInformationPropertyConstants.h" #include #include #include #include mitk::DICOMImageBlockDescriptor::DICOMImageBlockDescriptor() : m_ReaderImplementationLevel( SOPClassUnknown ) , m_PropertyList( PropertyList::New() ) +, m_SplitReason(DICOMSplitReason::New()) , m_TagCache( nullptr ) , m_PropertiesOutOfDate( true ) { m_PropertyFunctor = &mitk::DICOMImageBlockDescriptor::GetPropertyForDICOMValues; } mitk::DICOMImageBlockDescriptor::~DICOMImageBlockDescriptor() { } mitk::DICOMImageBlockDescriptor::DICOMImageBlockDescriptor( const DICOMImageBlockDescriptor& other ) : m_ImageFrameList( other.m_ImageFrameList ) , m_MitkImage( other.m_MitkImage ) , m_SliceIsLoaded( other.m_SliceIsLoaded ) , m_ReaderImplementationLevel( other.m_ReaderImplementationLevel ) , m_TiltInformation( other.m_TiltInformation ) , m_PropertyList( other.m_PropertyList->Clone() ) +, m_SplitReason( other.m_SplitReason->Clone() ) , m_TagCache( other.m_TagCache ) , m_PropertiesOutOfDate( other.m_PropertiesOutOfDate ) , m_AdditionalTagMap(other.m_AdditionalTagMap) , m_FoundAdditionalTags(other.m_FoundAdditionalTags) , m_PropertyFunctor(other.m_PropertyFunctor) { if ( m_MitkImage ) { m_MitkImage = m_MitkImage->Clone(); } } mitk::DICOMImageBlockDescriptor& mitk::DICOMImageBlockDescriptor:: operator=( const DICOMImageBlockDescriptor& other ) { if ( this != &other ) { m_ImageFrameList = other.m_ImageFrameList; m_MitkImage = other.m_MitkImage; m_SliceIsLoaded = other.m_SliceIsLoaded; m_ReaderImplementationLevel = other.m_ReaderImplementationLevel; m_TiltInformation = other.m_TiltInformation; m_AdditionalTagMap = other.m_AdditionalTagMap; m_FoundAdditionalTags = other.m_FoundAdditionalTags; m_PropertyFunctor = other.m_PropertyFunctor; if ( other.m_PropertyList ) { m_PropertyList = other.m_PropertyList->Clone(); } + if (other.m_SplitReason) + { + m_SplitReason = other.m_SplitReason->Clone(); + } + if ( other.m_MitkImage ) { m_MitkImage = other.m_MitkImage->Clone(); } m_TagCache = other.m_TagCache; m_PropertiesOutOfDate = other.m_PropertiesOutOfDate; } return *this; } mitk::DICOMTagList mitk::DICOMImageBlockDescriptor::GetTagsOfInterest() { DICOMTagList completeList; completeList.push_back( DICOMTag( 0x0018, 0x1164 ) ); // pixel spacing completeList.push_back( DICOMTag( 0x0028, 0x0030 ) ); // imager pixel spacing completeList.push_back( DICOMTag( 0x0008, 0x0018 ) ); // sop instance UID completeList.push_back( DICOMTag( 0x0008, 0x0016 ) ); // sop class UID completeList.push_back( DICOMTag( 0x0020, 0x0011 ) ); // series number completeList.push_back( DICOMTag( 0x0008, 0x1030 ) ); // study description completeList.push_back( DICOMTag( 0x0008, 0x103e ) ); // series description completeList.push_back( DICOMTag( 0x0008, 0x0060 ) ); // modality completeList.push_back( DICOMTag( 0x0018, 0x0024 ) ); // sequence name completeList.push_back( DICOMTag( 0x0020, 0x0037 ) ); // image orientation completeList.push_back( DICOMTag( 0x0020, 0x1041 ) ); // slice location completeList.push_back( DICOMTag( 0x0020, 0x0012 ) ); // acquisition number completeList.push_back( DICOMTag( 0x0020, 0x0013 ) ); // instance number completeList.push_back( DICOMTag( 0x0020, 0x0032 ) ); // image position patient completeList.push_back( DICOMTag( 0x0028, 0x1050 ) ); // window center completeList.push_back( DICOMTag( 0x0028, 0x1051 ) ); // window width completeList.push_back( DICOMTag( 0x0008, 0x0008 ) ); // image type completeList.push_back( DICOMTag( 0x0028, 0x0004 ) ); // photometric interpretation return completeList; } void mitk::DICOMImageBlockDescriptor::SetAdditionalTagsOfInterest( const AdditionalTagsMapType& tagMap) { m_AdditionalTagMap = tagMap; } void mitk::DICOMImageBlockDescriptor::SetTiltInformation( const GantryTiltInformation& info ) { m_TiltInformation = info; } const mitk::GantryTiltInformation mitk::DICOMImageBlockDescriptor::GetTiltInformation() const { return m_TiltInformation; } void mitk::DICOMImageBlockDescriptor::SetImageFrameList( const DICOMImageFrameList& framelist ) { m_ImageFrameList = framelist; m_SliceIsLoaded.resize( framelist.size() ); m_SliceIsLoaded.assign( framelist.size(), false ); m_PropertiesOutOfDate = true; } const mitk::DICOMImageFrameList& mitk::DICOMImageBlockDescriptor::GetImageFrameList() const { return m_ImageFrameList; } void mitk::DICOMImageBlockDescriptor::SetMitkImage( Image::Pointer image ) { if ( m_MitkImage != image ) { if ( m_TagCache.IsExpired() ) { MITK_ERROR << "Unable to describe MITK image with properties without a tag-cache object!"; m_MitkImage = nullptr; return; } if ( m_ImageFrameList.empty() ) { MITK_ERROR << "Unable to describe MITK image with properties without a frame list!"; m_MitkImage = nullptr; return; } // Should verify that the image matches m_ImageFrameList and m_TagCache // however, this is hard to do without re-analyzing all // TODO we should at least make sure that the number of frames is identical (plus rows/columns, // orientation) // without gantry tilt correction, we can also check image origin m_MitkImage = this->DescribeImageWithProperties( this->FixupSpacing( image ) ); } } mitk::Image::Pointer mitk::DICOMImageBlockDescriptor::GetMitkImage() const { return m_MitkImage; } mitk::Image::Pointer mitk::DICOMImageBlockDescriptor::FixupSpacing( Image* mitkImage ) { if ( mitkImage ) { Vector3D imageSpacing = mitkImage->GetGeometry()->GetSpacing(); ScalarType desiredSpacingX = imageSpacing[0]; ScalarType desiredSpacingY = imageSpacing[1]; this->GetDesiredMITKImagePixelSpacing( desiredSpacingX, desiredSpacingY ); // prefer pixel spacing over imager pixel spacing if ( desiredSpacingX <= 0 || desiredSpacingY <= 0 ) { return mitkImage; } MITK_DEBUG << "Loaded image with spacing " << imageSpacing[0] << ", " << imageSpacing[1]; MITK_DEBUG << "Found correct spacing info " << desiredSpacingX << ", " << desiredSpacingY; imageSpacing[0] = desiredSpacingX; imageSpacing[1] = desiredSpacingY; mitkImage->GetGeometry()->SetSpacing( imageSpacing ); } return mitkImage; } void mitk::DICOMImageBlockDescriptor::SetSliceIsLoaded( unsigned int index, bool isLoaded ) { if ( index < m_SliceIsLoaded.size() ) { m_SliceIsLoaded[index] = isLoaded; } else { std::stringstream ss; ss << "Index " << index << " out of range (" << m_SliceIsLoaded.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } bool mitk::DICOMImageBlockDescriptor::IsSliceLoaded( unsigned int index ) const { if ( index < m_SliceIsLoaded.size() ) { return m_SliceIsLoaded[index]; } else { std::stringstream ss; ss << "Index " << index << " out of range (" << m_SliceIsLoaded.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } bool mitk::DICOMImageBlockDescriptor::AllSlicesAreLoaded() const { bool allLoaded = true; for ( auto iter = m_SliceIsLoaded.cbegin(); iter != m_SliceIsLoaded.cend(); ++iter ) { allLoaded &= *iter; } return allLoaded; } /* PS defined IPS defined PS==IPS 0 0 --> UNKNOWN spacing, loader will invent 0 1 --> spacing as at detector surface 1 0 --> spacing as in patient 1 1 0 --> detector surface spacing CORRECTED for geometrical magnifications: spacing as in patient 1 1 1 --> detector surface spacing NOT corrected for geometrical magnifications: spacing as at detector */ mitk::PixelSpacingInterpretation mitk::DICOMImageBlockDescriptor::GetPixelSpacingInterpretation() const { if ( m_ImageFrameList.empty() || m_TagCache.IsExpired() ) { MITK_ERROR << "Invalid call to GetPixelSpacingInterpretation. Need to have initialized tag-cache!"; return SpacingUnknown; } const std::string pixelSpacing = this->GetPixelSpacing(); const std::string imagerPixelSpacing = this->GetImagerPixelSpacing(); if ( pixelSpacing.empty() ) { if ( imagerPixelSpacing.empty() ) { return SpacingUnknown; } else { return SpacingAtDetector; } } else // Pixel Spacing defined { if ( imagerPixelSpacing.empty() ) { return SpacingInPatient; } else if ( pixelSpacing != imagerPixelSpacing ) { return SpacingInPatient; } else { return SpacingAtDetector; } } } std::string mitk::DICOMImageBlockDescriptor::GetPixelSpacing() const { auto tagCache = m_TagCache.Lock(); if ( m_ImageFrameList.empty() || tagCache.IsNull() ) { MITK_ERROR << "Invalid call to GetPixelSpacing. Need to have initialized tag-cache!"; return std::string( "" ); } static const DICOMTag tagPixelSpacing( 0x0028, 0x0030 ); return tagCache->GetTagValue( m_ImageFrameList.front(), tagPixelSpacing ).value; } std::string mitk::DICOMImageBlockDescriptor::GetImagerPixelSpacing() const { auto tagCache = m_TagCache.Lock(); if ( m_ImageFrameList.empty() || tagCache.IsNull() ) { MITK_ERROR << "Invalid call to GetImagerPixelSpacing. Need to have initialized tag-cache!"; return std::string( "" ); } static const DICOMTag tagImagerPixelSpacing( 0x0018, 0x1164 ); return tagCache->GetTagValue( m_ImageFrameList.front(), tagImagerPixelSpacing ).value; } void mitk::DICOMImageBlockDescriptor::GetDesiredMITKImagePixelSpacing( ScalarType& spacingX, ScalarType& spacingY ) const { const std::string pixelSpacing = this->GetPixelSpacing(); // preference for "in patient" pixel spacing if ( !DICOMStringToSpacing( pixelSpacing, spacingX, spacingY ) ) { const std::string imagerPixelSpacing = this->GetImagerPixelSpacing(); // fallback to "on detector" spacing if ( !DICOMStringToSpacing( imagerPixelSpacing, spacingX, spacingY ) ) { // at this point we have no hints whether the spacing is correct // do a quick sanity check and either trust in the input or set both to 1 // We assume neither spacing to be negative, zero or unexpectedly large for // medical images if (spacingX < mitk::eps || spacingX > 1000 || spacingY < mitk::eps || spacingY > 1000) { spacingX = spacingY = 1.0; } } } } void mitk::DICOMImageBlockDescriptor::SetProperty( const std::string& key, BaseProperty* value ) { m_PropertyList->SetProperty( key, value ); } mitk::BaseProperty* mitk::DICOMImageBlockDescriptor::GetProperty( const std::string& key ) const { this->UpdateImageDescribingProperties(); return m_PropertyList->GetProperty( key ); } std::string mitk::DICOMImageBlockDescriptor::GetPropertyAsString( const std::string& key ) const { this->UpdateImageDescribingProperties(); const mitk::BaseProperty::Pointer property = m_PropertyList->GetProperty( key ); if ( property.IsNotNull() ) { return property->GetValueAsString(); } else { return std::string( "" ); } } +const mitk::DICOMSplitReason* mitk::DICOMImageBlockDescriptor::GetSplitReason() const +{ + return m_SplitReason; +} + +mitk::DICOMSplitReason* mitk::DICOMImageBlockDescriptor::GetSplitReason() +{ + return m_SplitReason; +} + +void mitk::DICOMImageBlockDescriptor::SetSplitReason(DICOMSplitReason* reason) +{ + m_SplitReason = reason; +} + void mitk::DICOMImageBlockDescriptor::SetFlag( const std::string& key, bool value ) { m_PropertyList->ReplaceProperty( key, BoolProperty::New( value ) ); } bool mitk::DICOMImageBlockDescriptor::GetFlag( const std::string& key, bool defaultValue ) const { this->UpdateImageDescribingProperties(); BoolProperty::ConstPointer boolProp = dynamic_cast( this->GetProperty( key ) ); if ( boolProp.IsNotNull() ) { return boolProp->GetValue(); } else { return defaultValue; } } void mitk::DICOMImageBlockDescriptor::SetIntProperty( const std::string& key, int value ) { m_PropertyList->ReplaceProperty( key, IntProperty::New( value ) ); } int mitk::DICOMImageBlockDescriptor::GetIntProperty( const std::string& key, int defaultValue ) const { this->UpdateImageDescribingProperties(); IntProperty::ConstPointer intProp = dynamic_cast( this->GetProperty( key ) ); if ( intProp.IsNotNull() ) { return intProp->GetValue(); } else { return defaultValue; } } double mitk::DICOMImageBlockDescriptor::stringtodouble( const std::string& str ) const { double d; std::string trimmedstring( str ); try { trimmedstring = trimmedstring.erase( trimmedstring.find_last_not_of( " \n\r\t" ) + 1 ); } catch ( ... ) { // no last not of } std::string firstcomponent( trimmedstring ); try { firstcomponent = trimmedstring.erase( trimmedstring.find_first_of( "\\" ) ); } catch ( ... ) { // no last not of } std::istringstream converter( firstcomponent ); if ( !firstcomponent.empty() && ( converter >> d ) && converter.eof() ) { return d; } else { throw std::invalid_argument( "Argument is not a convertible number" ); } } mitk::Image::Pointer mitk::DICOMImageBlockDescriptor::DescribeImageWithProperties( Image* mitkImage ) { // TODO: this is a collection of properties that have been provided by the // legacy DicomSeriesReader. // We should at some point clean up this collection and name them in a more // consistent way! if ( !mitkImage ) return mitkImage; mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_FILES()), this->GetProperty("filenamesForSlices")); mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_PIXEL_SPACING_INTERPRETATION_STRING()), StringProperty::New(PixelSpacingInterpretationToString(this->GetPixelSpacingInterpretation()))); mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_PIXEL_SPACING_INTERPRETATION()), GenericProperty::New(this->GetPixelSpacingInterpretation())); mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_IMPLEMENTATION_LEVEL_STRING()), StringProperty::New(ReaderImplementationLevelToString(m_ReaderImplementationLevel))); mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_IMPLEMENTATION_LEVEL()), GenericProperty::New(m_ReaderImplementationLevel)); mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_GANTRY_TILT_CORRECTED()), BoolProperty::New(this->GetTiltInformation().IsRegularGantryTilt())); mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_3D_plus_t()), BoolProperty::New(this->GetFlag("3D+t", false))); mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_GDCM()), StringProperty::New(gdcm::Version::GetVersion())); mitkImage->SetProperty(PropertyKeyPathToPropertyName(DICOMIOMetaInformationPropertyConstants::READER_DCMTK()), StringProperty::New(PACKAGE_VERSION)); // get all found additional tags of interest for (const auto &tag : m_FoundAdditionalTags) { BaseProperty* prop = this->GetProperty(tag); if (prop) { mitkImage->SetProperty(tag.c_str(), prop); } } ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// //// Deprecated properties should be removed sooner then later (see above) ///////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////// // first part: add some tags that describe individual slices // these properties are defined at analysis time (see UpdateImageDescribingProperties()) const char* propertyKeySliceLocation = "dicom.image.0020.1041"; const char* propertyKeyInstanceNumber = "dicom.image.0020.0013"; const char* propertyKeySOPInstanceUID = "dicom.image.0008.0018"; mitkImage->SetProperty( propertyKeySliceLocation, this->GetProperty( "sliceLocationForSlices" ) ); mitkImage->SetProperty( propertyKeyInstanceNumber, this->GetProperty( "instanceNumberForSlices" ) ); mitkImage->SetProperty( propertyKeySOPInstanceUID, this->GetProperty( "SOPInstanceUIDForSlices" ) ); mitkImage->SetProperty( "files", this->GetProperty( "filenamesForSlices_deprecated" ) ); // second part: add properties that describe the whole image block mitkImage->SetProperty( "dicomseriesreader.SOPClassUID", StringProperty::New( this->GetSOPClassUID() ) ); mitkImage->SetProperty( "dicomseriesreader.SOPClass", StringProperty::New( this->GetSOPClassUIDAsName() ) ); mitkImage->SetProperty( "dicomseriesreader.PixelSpacingInterpretationString", StringProperty::New( PixelSpacingInterpretationToString( this->GetPixelSpacingInterpretation() ) ) ); mitkImage->SetProperty( "dicomseriesreader.PixelSpacingInterpretation", GenericProperty::New( this->GetPixelSpacingInterpretation() ) ); mitkImage->SetProperty( "dicomseriesreader.ReaderImplementationLevelString", StringProperty::New( ReaderImplementationLevelToString( m_ReaderImplementationLevel ) ) ); mitkImage->SetProperty( "dicomseriesreader.ReaderImplementationLevel", GenericProperty::New( m_ReaderImplementationLevel ) ); mitkImage->SetProperty( "dicomseriesreader.GantyTiltCorrected", BoolProperty::New( this->GetTiltInformation().IsRegularGantryTilt() ) ); mitkImage->SetProperty( "dicomseriesreader.3D+t", BoolProperty::New( this->GetFlag( "3D+t", false ) ) ); // level window const std::string windowCenter = this->GetPropertyAsString( "windowCenter" ); const std::string windowWidth = this->GetPropertyAsString( "windowWidth" ); try { const double level = stringtodouble( windowCenter ); const double window = stringtodouble( windowWidth ); mitkImage->SetProperty( "levelwindow", LevelWindowProperty::New( LevelWindow( level, window ) ) ); } catch ( ... ) { // nothing, no levelwindow to be predicted... } const std::string modality = this->GetPropertyAsString( "modality" ); mitkImage->SetProperty( "modality", StringProperty::New( modality ) ); mitkImage->SetProperty( "dicom.pixel.PhotometricInterpretation", this->GetProperty( "photometricInterpretation" ) ); mitkImage->SetProperty( "dicom.image.imagetype", this->GetProperty( "imagetype" ) ); mitkImage->SetProperty( "dicom.study.StudyDescription", this->GetProperty( "studyDescription" ) ); mitkImage->SetProperty( "dicom.series.SeriesDescription", this->GetProperty( "seriesDescription" ) ); mitkImage->SetProperty( "dicom.pixel.Rows", this->GetProperty( "rows" ) ); mitkImage->SetProperty( "dicom.pixel.Columns", this->GetProperty( "columns" ) ); // fourth part: get something from ImageIO. BUT this needs to be created elsewhere. or not at all! return mitkImage; } void mitk::DICOMImageBlockDescriptor::SetReaderImplementationLevel( const ReaderImplementationLevel& level ) { m_ReaderImplementationLevel = level; } mitk::ReaderImplementationLevel mitk::DICOMImageBlockDescriptor::GetReaderImplementationLevel() const { return m_ReaderImplementationLevel; } std::string mitk::DICOMImageBlockDescriptor::GetSOPClassUID() const { auto tagCache = m_TagCache.Lock(); if ( !m_ImageFrameList.empty() && tagCache.IsNotNull() ) { static const DICOMTag tagSOPClassUID( 0x0008, 0x0016 ); return tagCache->GetTagValue( m_ImageFrameList.front(), tagSOPClassUID ).value; } else { MITK_ERROR << "Invalid call to DICOMImageBlockDescriptor::GetSOPClassUID(). Need to have initialized tag-cache!"; return std::string( "" ); } } std::string mitk::DICOMImageBlockDescriptor::GetSOPClassUIDAsName() const { if ( !m_ImageFrameList.empty() && !m_TagCache.IsExpired() ) { gdcm::UIDs uidKnowledge; uidKnowledge.SetFromUID( this->GetSOPClassUID().c_str() ); const char* name = uidKnowledge.GetName(); if ( name ) { return std::string( name ); } else { return std::string( "" ); } } else { MITK_ERROR << "Invalid call to DICOMImageBlockDescriptor::GetSOPClassUIDAsName(). Need to have " "initialized tag-cache!"; return std::string( "" ); } } int mitk::DICOMImageBlockDescriptor::GetNumberOfTimeSteps() const { int result = 1; this->m_PropertyList->GetIntProperty("timesteps", result); return result; }; int mitk::DICOMImageBlockDescriptor::GetNumberOfFramesPerTimeStep() const { const int numberOfTimesteps = this->GetNumberOfTimeSteps(); int numberOfFramesPerTimestep = this->m_ImageFrameList.size() / numberOfTimesteps; assert(int(double((double)this->m_ImageFrameList.size() / (double)numberOfTimesteps)) == numberOfFramesPerTimestep); // this should hold return numberOfFramesPerTimestep; }; void mitk::DICOMImageBlockDescriptor::SetTagCache( DICOMTagCache* privateCache ) { // this must only be used during loading and never afterwards m_TagCache = privateCache; } #define printPropertyRange( label, property_name ) \ \ { \ const std::string first = this->GetPropertyAsString( #property_name "First" ); \ const std::string last = this->GetPropertyAsString( #property_name "Last" ); \ if ( !first.empty() || !last.empty() ) \ { \ if ( first == last ) \ { \ os << " " label ": '" << first << "'" << std::endl; \ } \ else \ { \ os << " " label ": '" << first << "' - '" << last << "'" << std::endl; \ } \ } \ \ } #define printProperty( label, property_name ) \ \ { \ const std::string first = this->GetPropertyAsString( #property_name ); \ if ( !first.empty() ) \ { \ os << " " label ": '" << first << "'" << std::endl; \ } \ \ } #define printBool( label, commands ) \ \ { \ os << " " label ": '" << ( commands ? "yes" : "no" ) << "'" << std::endl; \ \ } void mitk::DICOMImageBlockDescriptor::Print(std::ostream& os, bool filenameDetails) const { os << " Number of Frames: '" << m_ImageFrameList.size() << "'" << std::endl; os << " SOP class: '" << this->GetSOPClassUIDAsName() << "'" << std::endl; printProperty( "Series Number", seriesNumber ); printProperty( "Study Description", studyDescription ); printProperty( "Series Description", seriesDescription ); printProperty( "Modality", modality ); printProperty( "Sequence Name", sequenceName ); printPropertyRange( "Slice Location", sliceLocation ); printPropertyRange( "Acquisition Number", acquisitionNumber ); printPropertyRange( "Instance Number", instanceNumber ); printPropertyRange( "Image Position", imagePositionPatient ); printProperty( "Image Orientation", orientation ); os << " Pixel spacing interpretation: '" << PixelSpacingInterpretationToString( this->GetPixelSpacingInterpretation() ) << "'" << std::endl; printBool( "Gantry Tilt", this->GetTiltInformation().IsRegularGantryTilt() ) // printBool("3D+t", this->GetFlag("3D+t",false)) // os << " MITK image loaded: '" << (this->GetMitkImage().IsNotNull() ? "yes" : "no") << "'" << // std::endl; if ( filenameDetails ) { os << " Files in this image block:" << std::endl; for ( auto frameIter = m_ImageFrameList.begin(); frameIter != m_ImageFrameList.end(); ++frameIter ) { os << " " << ( *frameIter )->Filename; if ( ( *frameIter )->FrameNo > 0 ) { os << ", " << ( *frameIter )->FrameNo; } os << std::endl; } } } #define storeTagValueToProperty( tag_name, tag_g, tag_e ) \ \ { \ const DICOMTag t( tag_g, tag_e ); \ const std::string tagValue = tagCache->GetTagValue( firstFrame, t ).value; \ const_cast( this ) \ ->SetProperty( #tag_name, StringProperty::New( tagValue ) ); \ \ } #define storeTagValueRangeToProperty( tag_name, tag_g, tag_e ) \ \ { \ const DICOMTag t( tag_g, tag_e ); \ const std::string tagValueFirst = tagCache->GetTagValue( firstFrame, t ).value; \ const std::string tagValueLast = tagCache->GetTagValue( lastFrame, t ).value; \ const_cast( this ) \ ->SetProperty( #tag_name "First", StringProperty::New( tagValueFirst ) ); \ const_cast( this ) \ ->SetProperty( #tag_name "Last", StringProperty::New( tagValueLast ) ); \ \ } void mitk::DICOMImageBlockDescriptor::UpdateImageDescribingProperties() const { if ( !m_PropertiesOutOfDate ) return; if ( !m_ImageFrameList.empty() ) { auto tagCache = m_TagCache.Lock(); if (tagCache.IsNull()) { MITK_ERROR << "Invalid call to DICOMImageBlockDescriptor::UpdateImageDescribingProperties(). Need to " "have initialized tag-cache!"; return; } const DICOMImageFrameInfo::Pointer firstFrame = m_ImageFrameList.front(); const DICOMImageFrameInfo::Pointer lastFrame = m_ImageFrameList.back(); // see macros above storeTagValueToProperty( seriesNumber, 0x0020, 0x0011 ); storeTagValueToProperty( studyDescription, 0x0008, 0x1030 ); storeTagValueToProperty( seriesDescription, 0x0008, 0x103e ); storeTagValueToProperty( modality, 0x0008, 0x0060 ); storeTagValueToProperty( sequenceName, 0x0018, 0x0024 ); storeTagValueToProperty( orientation, 0x0020, 0x0037 ); storeTagValueToProperty( rows, 0x0028, 0x0010 ); storeTagValueToProperty( columns, 0x0028, 0x0011 ); storeTagValueRangeToProperty( sliceLocation, 0x0020, 0x1041 ); storeTagValueRangeToProperty( acquisitionNumber, 0x0020, 0x0012 ); storeTagValueRangeToProperty( instanceNumber, 0x0020, 0x0013 ); storeTagValueRangeToProperty( imagePositionPatient, 0x0020, 0x0032 ); storeTagValueToProperty( windowCenter, 0x0028, 0x1050 ); storeTagValueToProperty( windowWidth, 0x0028, 0x1051 ); storeTagValueToProperty( imageType, 0x0008, 0x0008 ); storeTagValueToProperty( photometricInterpretation, 0x0028, 0x0004 ); // some per-image attributes // frames are just numbered starting from 0. timestep 1 (the second time-step) has frames starting at // (number-of-frames-per-timestep) // std::string propertyKeySliceLocation = "dicom.image.0020.1041"; // std::string propertyKeyInstanceNumber = "dicom.image.0020.0013"; // std::string propertyKeySOPInstanceNumber = "dicom.image.0008.0018"; StringLookupTable sliceLocationForSlices; StringLookupTable instanceNumberForSlices; StringLookupTable SOPInstanceUIDForSlices; StringLookupTable filenamesForSlices_deprecated; DICOMCachedValueLookupTable filenamesForSlices; const DICOMTag tagSliceLocation( 0x0020, 0x1041 ); const DICOMTag tagInstanceNumber( 0x0020, 0x0013 ); const DICOMTag tagSOPInstanceNumber( 0x0008, 0x0018 ); std::unordered_map additionalTagResultList; unsigned int slice(0); int timePoint(-1); const int framesPerTimeStep = this->GetNumberOfFramesPerTimeStep(); for ( auto frameIter = m_ImageFrameList.begin(); frameIter != m_ImageFrameList.end(); ++slice, ++frameIter ) { unsigned int zSlice = slice%framesPerTimeStep; if ( zSlice == 0) { timePoint++; } const std::string sliceLocation = tagCache->GetTagValue( *frameIter, tagSliceLocation ).value; sliceLocationForSlices.SetTableValue( slice, sliceLocation ); const std::string instanceNumber = tagCache->GetTagValue( *frameIter, tagInstanceNumber ).value; instanceNumberForSlices.SetTableValue( slice, instanceNumber ); const std::string sopInstanceUID = tagCache->GetTagValue( *frameIter, tagSOPInstanceNumber ).value; SOPInstanceUIDForSlices.SetTableValue( slice, sopInstanceUID ); const std::string filename = ( *frameIter )->Filename; filenamesForSlices_deprecated.SetTableValue( slice, filename ); filenamesForSlices.SetTableValue(slice, { static_cast(timePoint), zSlice, filename }); MITK_DEBUG << "Tag info for slice " << slice << ": SL '" << sliceLocation << "' IN '" << instanceNumber << "' SOP instance UID '" << sopInstanceUID << "'"; for (const auto& tag : m_AdditionalTagMap) { const DICOMTagCache::FindingsListType findings = tagCache->GetTagValue( *frameIter, tag.first ); for (const auto& finding : findings) { if (finding.isValid) { std::string propKey = (tag.second.empty()) ? DICOMTagPathToPropertyName(finding.path) : tag.second; DICOMCachedValueInfo info{ static_cast(timePoint), zSlice, finding.value }; additionalTagResultList[propKey].SetTableValue(slice, info); } } } } // add property or properties with proper names auto* thisInstance = const_cast( this ); thisInstance->SetProperty( "sliceLocationForSlices", StringLookupTableProperty::New( sliceLocationForSlices ) ); thisInstance->SetProperty( "instanceNumberForSlices", StringLookupTableProperty::New( instanceNumberForSlices ) ); thisInstance->SetProperty( "SOPInstanceUIDForSlices", StringLookupTableProperty::New( SOPInstanceUIDForSlices ) ); thisInstance->SetProperty( "filenamesForSlices_deprecated", StringLookupTableProperty::New( filenamesForSlices_deprecated ) ); thisInstance->SetProperty("filenamesForSlices", m_PropertyFunctor(filenamesForSlices)); //add properties for additional tags of interest for ( auto iter = additionalTagResultList.cbegin(); iter != additionalTagResultList.cend(); ++iter ) { thisInstance->SetProperty( iter->first, m_PropertyFunctor( iter->second ) ); thisInstance->m_FoundAdditionalTags.insert(m_FoundAdditionalTags.cend(),iter->first); } m_PropertiesOutOfDate = false; } } mitk::BaseProperty::Pointer mitk::DICOMImageBlockDescriptor::GetPropertyForDICOMValues(const DICOMCachedValueLookupTable& cacheLookupTable) { const auto& lookupTable = cacheLookupTable.GetLookupTable(); typedef std::pair PairType; if ( std::adjacent_find( lookupTable.cbegin(), lookupTable.cend(), []( const PairType& lhs, const PairType& rhs ) { return lhs.second.Value != rhs.second.Value; } ) == lookupTable.cend() ) { return static_cast( mitk::StringProperty::New(cacheLookupTable.GetTableValue(0).Value).GetPointer()); } StringLookupTable stringTable; for (const auto &element : lookupTable) { stringTable.SetTableValue(element.first, element.second.Value); } return static_cast( mitk::StringLookupTableProperty::New(stringTable).GetPointer()); } void mitk::DICOMImageBlockDescriptor::SetTagLookupTableToPropertyFunctor( TagLookupTableToPropertyFunctor functor ) { if ( functor != nullptr ) { m_PropertyFunctor = functor; } } mitk::BaseProperty::ConstPointer mitk::DICOMImageBlockDescriptor::GetConstProperty(const std::string &propertyKey, const std::string &/*contextName*/, bool /*fallBackOnDefaultContext*/) const { this->UpdateImageDescribingProperties(); return m_PropertyList->GetConstProperty(propertyKey); }; std::vector mitk::DICOMImageBlockDescriptor::GetPropertyKeys(const std::string &/*contextName*/, bool /*includeDefaultContext*/) const { this->UpdateImageDescribingProperties(); return m_PropertyList->GetPropertyKeys(); }; std::vector mitk::DICOMImageBlockDescriptor::GetPropertyContextNames() const { return std::vector(); }; diff --git a/Modules/DICOM/src/mitkDICOMTagBasedSorter.cpp b/Modules/DICOM/src/mitkDICOMTagBasedSorter.cpp index fc354a7bd6..7b43b41df4 100644 --- a/Modules/DICOM/src/mitkDICOMTagBasedSorter.cpp +++ b/Modules/DICOM/src/mitkDICOMTagBasedSorter.cpp @@ -1,601 +1,635 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkDICOMTagBasedSorter.h" #include #include mitk::DICOMTagBasedSorter::CutDecimalPlaces ::CutDecimalPlaces(unsigned int precision) :m_Precision(precision) { } mitk::DICOMTagBasedSorter::CutDecimalPlaces ::CutDecimalPlaces(const CutDecimalPlaces& other) :m_Precision(other.m_Precision) { } std::string mitk::DICOMTagBasedSorter::CutDecimalPlaces ::operator()(const std::string& input) const { // be a bit tolerant for tags such as image orientation orientation, let only the first few digits matter (https://phabricator.mitk.org/T12263) // iterate all fields, convert each to a number, cut this number as configured, then return a concatenated string with all cut-off numbers std::ostringstream resultString; resultString.str(std::string()); resultString.clear(); resultString.setf(std::ios::fixed, std::ios::floatfield); resultString.precision(m_Precision); std::stringstream ss(input); ss.str(input); ss.clear(); std::string item; double number(0); std::istringstream converter(item); while (std::getline(ss, item, '\\')) { converter.str(item); converter.clear(); if (converter >> number && converter.eof()) { // converted to double resultString << number; } else { // did not convert to double resultString << item; // just paste the unmodified string } if (!ss.eof()) { resultString << "\\"; } } return resultString.str(); } mitk::DICOMTagBasedSorter::TagValueProcessor* mitk::DICOMTagBasedSorter::CutDecimalPlaces ::Clone() const { return new CutDecimalPlaces(*this); } unsigned int mitk::DICOMTagBasedSorter::CutDecimalPlaces ::GetPrecision() const { return m_Precision; } mitk::DICOMTagBasedSorter ::DICOMTagBasedSorter() :DICOMDatasetSorter() ,m_StrictSorting(m_DefaultStrictSorting) ,m_ExpectDistanceOne(m_DefaultExpectDistanceOne) { } mitk::DICOMTagBasedSorter ::~DICOMTagBasedSorter() { for(auto ti = m_TagValueProcessor.cbegin(); ti != m_TagValueProcessor.cend(); ++ti) { delete ti->second; } } mitk::DICOMTagBasedSorter ::DICOMTagBasedSorter(const DICOMTagBasedSorter& other ) :DICOMDatasetSorter(other) ,m_DistinguishingTags( other.m_DistinguishingTags ) ,m_SortCriterion( other.m_SortCriterion ) ,m_StrictSorting( other.m_StrictSorting ) ,m_ExpectDistanceOne( other.m_ExpectDistanceOne ) { for(auto ti = other.m_TagValueProcessor.cbegin(); ti != other.m_TagValueProcessor.cend(); ++ti) { m_TagValueProcessor[ti->first] = ti->second->Clone(); } } mitk::DICOMTagBasedSorter& mitk::DICOMTagBasedSorter ::operator=(const DICOMTagBasedSorter& other) { if (this != &other) { DICOMDatasetSorter::operator=(other); m_DistinguishingTags = other.m_DistinguishingTags; m_SortCriterion = other.m_SortCriterion; m_StrictSorting = other.m_StrictSorting; m_ExpectDistanceOne = other.m_ExpectDistanceOne; for(auto ti = other.m_TagValueProcessor.cbegin(); ti != other.m_TagValueProcessor.cend(); ++ti) { m_TagValueProcessor[ti->first] = ti->second->Clone(); } } return *this; } bool mitk::DICOMTagBasedSorter ::operator==(const DICOMDatasetSorter& other) const { if (const auto* otherSelf = dynamic_cast(&other)) { if (this->m_StrictSorting != otherSelf->m_StrictSorting) return false; if (this->m_ExpectDistanceOne != otherSelf->m_ExpectDistanceOne) return false; bool allTagsPresentAndEqual(true); if (this->m_DistinguishingTags.size() != otherSelf->m_DistinguishingTags.size()) return false; for (auto myTag = this->m_DistinguishingTags.cbegin(); myTag != this->m_DistinguishingTags.cend(); ++myTag) { allTagsPresentAndEqual &= (std::find( otherSelf->m_DistinguishingTags.cbegin(), otherSelf->m_DistinguishingTags.cend(), *myTag ) != otherSelf->m_DistinguishingTags.cend()); // other contains this tags // since size is equal, we don't need to check the inverse } if (!allTagsPresentAndEqual) return false; if (this->m_SortCriterion.IsNotNull() && otherSelf->m_SortCriterion.IsNotNull()) { return *(this->m_SortCriterion) == *(otherSelf->m_SortCriterion); } else { return this->m_SortCriterion.IsNull() && otherSelf->m_SortCriterion.IsNull(); } } else { return false; } } void mitk::DICOMTagBasedSorter ::PrintConfiguration(std::ostream& os, const std::string& indent) const { os << indent << "Tag based sorting " << "(strict=" << (m_StrictSorting?"true":"false") << ", expectDistanceOne=" << (m_ExpectDistanceOne?"true":"false") << "):" << std::endl; for (auto tagIter = m_DistinguishingTags.begin(); tagIter != m_DistinguishingTags.end(); ++tagIter) { os << indent << " Split on "; tagIter->Print(os); os << std::endl; } DICOMSortCriterion::ConstPointer crit = m_SortCriterion.GetPointer(); while (crit.IsNotNull()) { os << indent << " Sort by "; crit->Print(os); os << std::endl; crit = crit->GetSecondaryCriterion(); } } void mitk::DICOMTagBasedSorter ::SetStrictSorting(bool strict) { m_StrictSorting = strict; } bool mitk::DICOMTagBasedSorter ::GetStrictSorting() const { return m_StrictSorting; } void mitk::DICOMTagBasedSorter ::SetExpectDistanceOne(bool strict) { m_ExpectDistanceOne = strict; } bool mitk::DICOMTagBasedSorter ::GetExpectDistanceOne() const { return m_ExpectDistanceOne; } mitk::DICOMTagList mitk::DICOMTagBasedSorter ::GetTagsOfInterest() { DICOMTagList allTags = m_DistinguishingTags; if (m_SortCriterion.IsNotNull()) { const DICOMTagList sortingRelevantTags = m_SortCriterion->GetAllTagsOfInterest(); allTags.insert( allTags.end(), sortingRelevantTags.cbegin(), sortingRelevantTags.cend() ); // append } return allTags; } mitk::DICOMTagList mitk::DICOMTagBasedSorter ::GetDistinguishingTags() const { return m_DistinguishingTags; } const mitk::DICOMTagBasedSorter::TagValueProcessor* mitk::DICOMTagBasedSorter ::GetTagValueProcessorForDistinguishingTag(const DICOMTag& tag) const { auto loc = m_TagValueProcessor.find(tag); if (loc != m_TagValueProcessor.cend()) { return loc->second; } else { return nullptr; } } void mitk::DICOMTagBasedSorter ::AddDistinguishingTag( const DICOMTag& tag, TagValueProcessor* tagValueProcessor ) { m_DistinguishingTags.push_back(tag); m_TagValueProcessor[tag] = tagValueProcessor; } void mitk::DICOMTagBasedSorter ::SetSortCriterion( DICOMSortCriterion::ConstPointer criterion ) { m_SortCriterion = criterion; } mitk::DICOMSortCriterion::ConstPointer mitk::DICOMTagBasedSorter ::GetSortCriterion() const { return m_SortCriterion; } void mitk::DICOMTagBasedSorter ::Sort() { + SplitReasonListType splitReasons; + // 1. split - // 2. sort each group - GroupIDToListType groups = this->SplitInputGroups(); - GroupIDToListType& sortedGroups = this->SortGroups( groups ); + GroupIDToListType groups = this->SplitInputGroups(splitReasons); + + // 2. sort each group (can also lead to a split due to distance) + GroupIDToListType& sortedGroups = this->SortGroups( groups, splitReasons); // 3. define output this->SetNumberOfOutputs(sortedGroups.size()); unsigned int outputIndex(0); for (auto groupIter = sortedGroups.cbegin(); groupIter != sortedGroups.cend(); ++outputIndex, ++groupIter) { - this->SetOutput(outputIndex, groupIter->second); + this->SetOutput(outputIndex, groupIter->second, splitReasons[groupIter->first]); } } std::string mitk::DICOMTagBasedSorter ::BuildGroupID( DICOMDatasetAccess* dataset ) { // just concatenate all tag values assert(dataset); std::stringstream groupID; groupID << "g"; for (auto tagIter = m_DistinguishingTags.cbegin(); tagIter != m_DistinguishingTags.cend(); ++tagIter) { groupID << tagIter->GetGroup() << tagIter->GetElement(); // make group/element part of the id to cover empty tags DICOMDatasetFinding rawTagValue = dataset->GetTagValueAsString(*tagIter); std::string processedTagValue; if ( m_TagValueProcessor[*tagIter] != nullptr && rawTagValue.isValid) { processedTagValue = (*m_TagValueProcessor[*tagIter])(rawTagValue.value); } else { processedTagValue = rawTagValue.value; } groupID << "#" << processedTagValue; } // shorten ID? return groupID.str(); } mitk::DICOMTagBasedSorter::GroupIDToListType mitk::DICOMTagBasedSorter -::SplitInputGroups() +::SplitInputGroups(SplitReasonListType& splitReasons) { DICOMDatasetList input = GetInput(); // copy GroupIDToListType listForGroupID; for (auto dsIter = input.cbegin(); dsIter != input.cend(); ++dsIter) { DICOMDatasetAccess* dataset = *dsIter; assert(dataset); const std::string groupID = this->BuildGroupID( dataset ); MITK_DEBUG << "Group ID for for " << dataset->GetFilenameIfAvailable() << ": " << groupID; listForGroupID[groupID].push_back(dataset); } MITK_DEBUG << "After tag based splitting: " << listForGroupID.size() << " groups"; + splitReasons.clear(); + if (listForGroupID.size() == 1) + { + //no split -> no reason + splitReasons[listForGroupID.begin()->first] = DICOMSplitReason::New(); + } + else + { + for (auto& [key, value] : listForGroupID) + { + auto reason = DICOMSplitReason::New(); + reason->AddReason(DICOMSplitReason::ReasonType::ValueSplitDifference); + splitReasons[key] = reason; + } + } + return listForGroupID; } mitk::DICOMTagBasedSorter::GroupIDToListType& mitk::DICOMTagBasedSorter -::SortGroups(GroupIDToListType& groups) +::SortGroups(GroupIDToListType& groups, SplitReasonListType& splitReasons) { if (m_SortCriterion.IsNotNull()) { /* Three steps here: 1. sort within each group - this may result in orders such as 1 2 3 4 6 7 8 10 12 13 14 2. create new groups by enforcing consecutive order within each group - resorts above example like 1 2 3 4 ; 6 7 8 ; 10 ; 12 13 14 3. sort all of the groups (not WITHIN each group) by their first frame - if earlier "distinguish" steps created groups like 6 7 8 ; 1 2 3 4 ; 10, then this step would sort them like 1 2 3 4 ; 6 7 8 ; 10 */ // Step 1: sort within the groups // for each output // sort by all configured tags, use secondary tags when equal or empty // make configurable: // - sorting order (ascending, descending) // - sort numerically // - ... ? #ifdef MBILOG_ENABLE_DEBUG unsigned int groupIndex(0); #endif for (auto gIter = groups.begin(); gIter != groups.end(); #ifdef MBILOG_ENABLE_DEBUG ++groupIndex, #endif ++gIter) { DICOMDatasetList& dsList = gIter->second; #ifdef MBILOG_ENABLE_DEBUG MITK_DEBUG << " --------------------------------------------------------------------------------"; MITK_DEBUG << " DICOMTagBasedSorter before sorting group : " << groupIndex; for (auto oi = dsList.begin(); oi != dsList.cend(); ++oi) { MITK_DEBUG << " INPUT : " << (*oi)->GetFilenameIfAvailable(); } #endif // #ifdef MBILOG_ENABLE_DEBUG std::sort( dsList.begin(), dsList.end(), ParameterizedDatasetSort( m_SortCriterion ) ); #ifdef MBILOG_ENABLE_DEBUG MITK_DEBUG << " --------------------------------------------------------------------------------"; MITK_DEBUG << " DICOMTagBasedSorter after sorting group : " << groupIndex; for (auto oi = dsList.cbegin(); oi != dsList.cend(); ++oi) { MITK_DEBUG << " OUTPUT : " << (*oi)->GetFilenameIfAvailable(); } MITK_DEBUG << " --------------------------------------------------------------------------------"; #endif // MBILOG_ENABLE_DEBUG } GroupIDToListType consecutiveGroups; + SplitReasonListType consecutiveReasons; if (m_StrictSorting) { // Step 2: create new groups by enforcing consecutive order within each group unsigned int groupIndex(0); for (auto gIter = groups.begin(); gIter != groups.end(); ++gIter) { std::stringstream groupKey; groupKey << std::setfill('0') << std::setw(6) << groupIndex++; + std::string groupKeyStr = groupKey.str(); DICOMDatasetList& dsList = gIter->second; + DICOMDatasetAccess* previousDS(nullptr); unsigned int dsIndex(0); double constantDistance(0.0); bool constantDistanceInitialized(false); for (auto dataset = dsList.cbegin(); dataset != dsList.cend(); ++dsIndex, ++dataset) { + bool splitted = false; if (dsIndex >0) // ignore the first dataset, we cannot check any distances yet.. { // for the second and every following dataset: // let the sorting criterion calculate a "distance" // if the distance is not 1, split off a new group! const double currentDistance = m_SortCriterion->NumericDistance(previousDS, *dataset); if (constantDistanceInitialized) { if (fabs(currentDistance - constantDistance) < fabs(constantDistance * 0.01)) // ok, deviation of up to 1% of distance is tolerated { // nothing to do, just ok MITK_DEBUG << "Checking currentDistance==" << currentDistance << ": small enough"; } //else if (currentDistance < mitk::eps) // close enough to 0 else { MITK_DEBUG << "Split consecutive group at index " << dsIndex << " (current distance " << currentDistance << ", constant distance " << constantDistance << ")"; // split! this is done by simply creating a new group (key) groupKey.str(std::string()); groupKey.clear(); groupKey << std::setfill('0') << std::setw(6) << groupIndex++; + groupKeyStr = groupKey.str(); + splitted = true; } } else { // second slice: learn about the expected distance! // heuristic: if distance is an integer, we check for a special case: // if the distance is integer and not 1/-1, then we assume // a missing slice right after the first slice // ==> split off slices // in all other cases: second dataset at this position, no need to split already, we are still learning about the images // addition to the above: when sorting by imagepositions, a distance other than 1 between the first two slices is // not unusual, actually expected... then we should not split if (m_ExpectDistanceOne) { if ((currentDistance - (int)currentDistance == 0.0) && fabs(currentDistance) != 1.0) // exact comparison. An integer should not be expressed as 1.000000000000000000000000001! { MITK_DEBUG << "Split consecutive group at index " << dsIndex << " (special case: expected distance 1 exactly)"; groupKey.str(std::string()); groupKey.clear(); groupKey << std::setfill('0') << std::setw(6) << groupIndex++; + groupKeyStr = groupKey.str(); + splitted = true; } } MITK_DEBUG << "Initialize strict distance to currentDistance=" << currentDistance; constantDistance = currentDistance; constantDistanceInitialized = true; } } - consecutiveGroups[groupKey.str()].push_back(*dataset); + consecutiveGroups[groupKeyStr].push_back(*dataset); + + if (consecutiveReasons.find(groupKeyStr) == consecutiveReasons.end()) + { + auto dsReason = splitReasons[gIter->first]->Clone(); + if (splitted) + dsReason->AddReason(DICOMSplitReason::ReasonType::ValueSortDistance); + consecutiveReasons[groupKeyStr] = dsReason; + } + previousDS = *dataset; } } } else { consecutiveGroups = groups; } // Step 3: sort all of the groups (not WITHIN each group) by their first frame /* build a list-1 of datasets with the first dataset one of each group sort this list-1 build a new result list-2: - iterate list-1, for each dataset - find the group that contains this dataset - add this group as the next element to list-2 return list-2 as the sorted output */ DICOMDatasetList firstSlices; for (auto gIter = consecutiveGroups.cbegin(); gIter != consecutiveGroups.cend(); ++gIter) { assert(!gIter->second.empty()); firstSlices.push_back(gIter->second.front()); } std::sort( firstSlices.begin(), firstSlices.end(), ParameterizedDatasetSort( m_SortCriterion ) ); GroupIDToListType sortedResultBlocks; + SplitReasonListType sortedResultsReasons; + unsigned int groupKeyValue(0); - for (auto firstSlice = firstSlices.cbegin(); - firstSlice != firstSlices.cend(); - ++firstSlice) + + for (auto& [key, group] : consecutiveGroups) { - for (auto gIter = consecutiveGroups.cbegin(); - gIter != consecutiveGroups.cend(); - ++groupKeyValue, ++gIter) - { - if (gIter->second.front() == *firstSlice) - { - std::stringstream groupKey; - groupKey << std::setfill('0') << std::setw(6) << groupKeyValue; // try more than 999,999 groups and you are doomed (your application already is) - sortedResultBlocks[groupKey.str()] = gIter->second; - } - } + auto findSliceIterator = std::find(firstSlices.begin(), firstSlices.end(), group.front()); + + std::stringstream groupKey; + groupKey << std::setfill('0') << std::setw(6) << std::distance(firstSlices.begin(),findSliceIterator); // try more than 999,999 groups and you are doomed (your application already is) + const auto groupKeyStr = groupKey.str(); + sortedResultBlocks[groupKeyStr] = group; + sortedResultsReasons[groupKeyStr] = consecutiveReasons[key]; } groups = sortedResultBlocks; + splitReasons = sortedResultsReasons; } #ifdef MBILOG_ENABLE_DEBUG unsigned int groupIndex( 0 ); for ( auto gIter = groups.begin(); gIter != groups.end(); ++groupIndex, ++gIter ) { DICOMDatasetList& dsList = gIter->second; MITK_DEBUG << " --------------------------------------------------------------------------------"; MITK_DEBUG << " DICOMTagBasedSorter after sorting group : " << groupIndex; for ( auto oi = dsList.begin(); oi != dsList.end(); ++oi ) { MITK_DEBUG << " OUTPUT : " << ( *oi )->GetFilenameIfAvailable(); } MITK_DEBUG << " --------------------------------------------------------------------------------"; } #endif // MBILOG_ENABLE_DEBUG return groups; } mitk::DICOMTagBasedSorter::ParameterizedDatasetSort ::ParameterizedDatasetSort(DICOMSortCriterion::ConstPointer criterion) :m_SortCriterion(criterion) { } bool mitk::DICOMTagBasedSorter::ParameterizedDatasetSort ::operator() (const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right) { assert(left); assert(right); assert(m_SortCriterion.IsNotNull()); return m_SortCriterion->IsLeftBeforeRight(left, right); } diff --git a/Modules/DICOM/src/mitkEquiDistantBlocksSorter.cpp b/Modules/DICOM/src/mitkEquiDistantBlocksSorter.cpp index e40fc2ad99..842fa67ad7 100644 --- a/Modules/DICOM/src/mitkEquiDistantBlocksSorter.cpp +++ b/Modules/DICOM/src/mitkEquiDistantBlocksSorter.cpp @@ -1,568 +1,603 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "mitkEquiDistantBlocksSorter.h" mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult -::SliceGroupingAnalysisResult() +::SliceGroupingAnalysisResult() : m_SplitReason(DICOMSplitReason::New()) { } -mitk::DICOMDatasetList +const mitk::DICOMDatasetList& mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult -::GetBlockDatasets() +::GetBlockDatasets() const { return m_GroupedFiles; } -mitk::DICOMDatasetList +const mitk::DICOMDatasetList& mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult -::GetUnsortedDatasets() +::GetUnsortedDatasets() const { return m_UnsortedFiles; } +const mitk::DICOMSplitReason* +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult::GetSplitReason() const +{ + return m_SplitReason; +} + +mitk::DICOMSplitReason* +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult::GetSplitReason() +{ + return m_SplitReason; +} + bool mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::ContainsGantryTilt() { return m_TiltInfo.IsRegularGantryTilt(); } void mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::AddFileToSortedBlock(DICOMDatasetAccess* dataset) { m_GroupedFiles.push_back( dataset ); } void mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::AddFileToUnsortedBlock(DICOMDatasetAccess* dataset) { m_UnsortedFiles.push_back( dataset ); } void mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::AddFilesToUnsortedBlock(const DICOMDatasetList& datasets) { m_UnsortedFiles.insert( m_UnsortedFiles.end(), datasets.begin(), datasets.end() ); } void mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::SetFirstFilenameOfBlock(const std::string& filename) { m_FirstFilenameOfBlock = filename; } std::string mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::GetFirstFilenameOfBlock() const { return m_FirstFilenameOfBlock; } void mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::SetLastFilenameOfBlock(const std::string& filename) { m_LastFilenameOfBlock = filename; } std::string mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::GetLastFilenameOfBlock() const { return m_LastFilenameOfBlock; } void mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::FlagGantryTilt(const GantryTiltInformation& tiltInfo) { m_TiltInfo = tiltInfo; } const mitk::GantryTiltInformation& mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::GetTiltInfo() const { return m_TiltInfo; } void mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::UndoPrematureGrouping() { assert( !m_GroupedFiles.empty() ); m_UnsortedFiles.insert( m_UnsortedFiles.begin(), m_GroupedFiles.back() ); m_GroupedFiles.pop_back(); m_TiltInfo = GantryTiltInformation(); } // ------------------------ end helper class mitk::EquiDistantBlocksSorter ::EquiDistantBlocksSorter() :DICOMDatasetSorter() ,m_AcceptTilt(false) ,m_ToleratedOriginOffset(0.3) ,m_ToleratedOriginOffsetIsAbsolute(false) ,m_AcceptTwoSlicesGroups(true) { } mitk::EquiDistantBlocksSorter ::EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other ) :DICOMDatasetSorter(other) ,m_AcceptTilt(other.m_AcceptTilt) ,m_ToleratedOriginOffset(other.m_ToleratedOriginOffset) ,m_ToleratedOriginOffsetIsAbsolute(other.m_ToleratedOriginOffsetIsAbsolute) ,m_AcceptTwoSlicesGroups(other.m_AcceptTwoSlicesGroups) { } mitk::EquiDistantBlocksSorter ::~EquiDistantBlocksSorter() { } bool mitk::EquiDistantBlocksSorter ::operator==(const DICOMDatasetSorter& other) const { if (const auto* otherSelf = dynamic_cast(&other)) { return this->m_AcceptTilt == otherSelf->m_AcceptTilt && this->m_ToleratedOriginOffsetIsAbsolute == otherSelf->m_ToleratedOriginOffsetIsAbsolute && this->m_AcceptTwoSlicesGroups == otherSelf->m_AcceptTwoSlicesGroups && (fabs(this->m_ToleratedOriginOffset - otherSelf->m_ToleratedOriginOffset) < eps); } else { return false; } } void mitk::EquiDistantBlocksSorter ::PrintConfiguration(std::ostream& os, const std::string& indent) const { std::stringstream ts; if (!m_ToleratedOriginOffsetIsAbsolute) { ts << "adaptive"; } else { ts << m_ToleratedOriginOffset << "mm"; } os << indent << "Sort into blocks of equidistant, well-aligned (tolerance " << ts.str() << ") slices " << (m_AcceptTilt ? "(accepting a gantry tilt)" : "") << std::endl; } void mitk::EquiDistantBlocksSorter ::SetAcceptTilt(bool accept) { m_AcceptTilt = accept; } bool mitk::EquiDistantBlocksSorter ::GetAcceptTilt() const { return m_AcceptTilt; } void mitk::EquiDistantBlocksSorter ::SetAcceptTwoSlicesGroups(bool accept) { m_AcceptTwoSlicesGroups = accept; } bool mitk::EquiDistantBlocksSorter ::GetAcceptTwoSlicesGroups() const { return m_AcceptTwoSlicesGroups; } mitk::EquiDistantBlocksSorter& mitk::EquiDistantBlocksSorter ::operator=(const EquiDistantBlocksSorter& other) { if (this != &other) { DICOMDatasetSorter::operator=(other); m_AcceptTilt = other.m_AcceptTilt; m_ToleratedOriginOffset = other.m_ToleratedOriginOffset; m_ToleratedOriginOffsetIsAbsolute = other.m_ToleratedOriginOffsetIsAbsolute; m_AcceptTwoSlicesGroups = other.m_AcceptTwoSlicesGroups; } return *this; } mitk::DICOMTagList mitk::EquiDistantBlocksSorter ::GetTagsOfInterest() { DICOMTagList tags; tags.push_back( DICOMTag(0x0020, 0x0032) ); // ImagePositionPatient tags.push_back( DICOMTag(0x0020, 0x0037) ); // ImageOrientationPatient tags.push_back( DICOMTag(0x0018, 0x1120) ); // GantryDetectorTilt return tags; } void mitk::EquiDistantBlocksSorter ::Sort() { DICOMDatasetList remainingInput = GetInput(); // copy typedef std::list OutputListType; - OutputListType outputs; m_SliceGroupingResults.clear(); while (!remainingInput.empty()) // repeat until all files are grouped somehow { - SliceGroupingAnalysisResult regularBlock = this->AnalyzeFileForITKImageSeriesReaderSpacingAssumption( remainingInput, m_AcceptTilt ); + auto regularBlock = this->AnalyzeFileForITKImageSeriesReaderSpacingAssumption( remainingInput, m_AcceptTilt ); #ifdef MBILOG_ENABLE_DEBUG DICOMDatasetList inBlock = regularBlock.GetBlockDatasets(); DICOMDatasetList laterBlock = regularBlock.GetUnsortedDatasets(); MITK_DEBUG << "Result: sorted 3D group with " << inBlock.size() << " files"; for (DICOMDatasetList::const_iterator diter = inBlock.cbegin(); diter != inBlock.cend(); ++diter) MITK_DEBUG << " IN " << (*diter)->GetFilenameIfAvailable(); for (DICOMDatasetList::const_iterator diter = laterBlock.cbegin(); diter != laterBlock.cend(); ++diter) MITK_DEBUG << " OUT " << (*diter)->GetFilenameIfAvailable(); #endif // MBILOG_ENABLE_DEBUG - outputs.push_back( regularBlock.GetBlockDatasets() ); + remainingInput = regularBlock->GetUnsortedDatasets(); + + if (remainingInput.empty() && !m_SliceGroupingResults.empty() && m_SliceGroupingResults.back()->GetSplitReason()->ReasonExists(DICOMSplitReason::ReasonType::OverlappingSlices)) + { + //if all inputs are processed and there is already a preceding grouping result that has overlapping as split reason, add also overlapping as split reason for the current block + regularBlock->GetSplitReason()->AddReason(DICOMSplitReason::ReasonType::OverlappingSlices); + } + m_SliceGroupingResults.push_back( regularBlock ); - remainingInput = regularBlock.GetUnsortedDatasets(); } - unsigned int numberOfOutputs = outputs.size(); + unsigned int numberOfOutputs = m_SliceGroupingResults.size(); this->SetNumberOfOutputs(numberOfOutputs); unsigned int outputIndex(0); - for (auto oIter = outputs.cbegin(); - oIter != outputs.cend(); + for (auto oIter = m_SliceGroupingResults.cbegin(); + oIter != m_SliceGroupingResults.cend(); ++outputIndex, ++oIter) { - this->SetOutput(outputIndex, *oIter); + this->SetOutput(outputIndex, (*oIter)->GetBlockDatasets(), (*oIter)->GetSplitReason()); } } void mitk::EquiDistantBlocksSorter ::SetToleratedOriginOffsetToAdaptive(double fractionOfInterSliceDistance) { m_ToleratedOriginOffset = fractionOfInterSliceDistance; m_ToleratedOriginOffsetIsAbsolute = false; if (m_ToleratedOriginOffset < 0.0) { MITK_WARN << "Call SetToleratedOriginOffsetToAdaptive() only with positive numbers between 0.0 and 1.0, read documentation!"; } if (m_ToleratedOriginOffset > 0.5) { MITK_WARN << "EquiDistantBlocksSorter is now accepting large errors, take care of measurements, they could appear at imprecise locations!"; } } void mitk::EquiDistantBlocksSorter ::SetToleratedOriginOffset(double millimeters) { m_ToleratedOriginOffset = millimeters; m_ToleratedOriginOffsetIsAbsolute = true; if (m_ToleratedOriginOffset < 0.0) { MITK_WARN << "Negative tolerance set to SetToleratedOriginOffset()!"; } } double mitk::EquiDistantBlocksSorter ::GetToleratedOriginOffset() const { return m_ToleratedOriginOffset; } bool mitk::EquiDistantBlocksSorter ::IsToleratedOriginOffsetAbsolute() const { return m_ToleratedOriginOffsetIsAbsolute; } std::string mitk::EquiDistantBlocksSorter ::ConstCharStarToString(const char* s) { return s ? std::string(s) : std::string(); } -mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +std::shared_ptr mitk::EquiDistantBlocksSorter ::AnalyzeFileForITKImageSeriesReaderSpacingAssumption( const DICOMDatasetList& datasets, bool groupImagesWithGantryTilt) { - SliceGroupingAnalysisResult result; + auto result = std::make_shared(); const DICOMTag tagImagePositionPatient = DICOMTag(0x0020,0x0032); // Image Position (Patient) const DICOMTag tagImageOrientation = DICOMTag(0x0020, 0x0037); // Image Orientation Vector3D fromFirstToSecondOrigin; fromFirstToSecondOrigin.Fill(0.0); bool fromFirstToSecondOriginInitialized(false); Point3D thisOrigin; thisOrigin.Fill(0.0f); Point3D lastOrigin; lastOrigin.Fill(0.0f); Point3D lastDifferentOrigin; lastDifferentOrigin.Fill(0.0f); bool lastOriginInitialized(false); MITK_DEBUG << "--------------------------------------------------------------------------------"; MITK_DEBUG << "Analyzing " << datasets.size() << " files for z-spacing assumption of ITK's ImageSeriesReader (group tilted: " << groupImagesWithGantryTilt << ")"; unsigned int fileIndex(0); double toleratedOriginError(0.005); // default: max. 1/10mm error when measurement crosses 20 slices in z direction (too strict? we don't know better) for (auto dsIter = datasets.cbegin(); dsIter != datasets.cend(); ++dsIter, ++fileIndex) { bool fileFitsIntoPattern(false); std::string thisOriginString; // Read tag value into point3D. PLEASE replace this by appropriate GDCM code if you figure out how to do that thisOriginString = (*dsIter)->GetTagValueAsString(tagImagePositionPatient).value; if (thisOriginString.empty()) { // don't let such files be in a common group. Everything without position information will be loaded as a single slice: // with standard DICOM files this can happen to: CR, DX, SC MITK_DEBUG << " ==> Sort away " << *dsIter << " for later analysis (no position information)"; // we already have one occupying this position - if ( result.GetBlockDatasets().empty() ) // nothing WITH position information yet + if ( result->GetBlockDatasets().empty() ) // nothing WITH position information yet { // ==> this is a group of its own, stop processing, come back later - result.AddFileToSortedBlock( *dsIter ); + result->AddFileToSortedBlock( *dsIter ); DICOMDatasetList remainingFiles; remainingFiles.insert( remainingFiles.end(), dsIter+1, datasets.end() ); - result.AddFilesToUnsortedBlock( remainingFiles ); - + result->AddFilesToUnsortedBlock( remainingFiles ); + result->GetSplitReason()->AddReason(DICOMSplitReason::ReasonType::ImagePostionMissing); fileFitsIntoPattern = false; break; // no files anymore } else { // ==> this does not match, consider later - result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis + result->AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis fileFitsIntoPattern = false; continue; // next file } } bool ignoredConversionError(-42); // hard to get here, no graceful way to react thisOrigin = DICOMStringToPoint3D( thisOriginString, ignoredConversionError ); + int missingSlicesCount = 0; + MITK_DEBUG << " " << fileIndex << " " << (*dsIter)->GetFilenameIfAvailable() << " at " /* << thisOriginString */ << "(" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")"; if ( lastOriginInitialized && (thisOrigin == lastOrigin) ) { MITK_DEBUG << " ==> Sort away " << *dsIter << " for separate time step"; // we already have one occupying this position - result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis + result->AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis + result->GetSplitReason()->AddReason(DICOMSplitReason::ReasonType::OverlappingSlices); fileFitsIntoPattern = false; } else { if (!fromFirstToSecondOriginInitialized && lastOriginInitialized) // calculate vector as soon as possible when we get a new position { fromFirstToSecondOrigin = thisOrigin - lastDifferentOrigin; fromFirstToSecondOriginInitialized = true; // classic mode without tolerance! if (!m_ToleratedOriginOffsetIsAbsolute) { MITK_DEBUG << "Distance of two slices: " << fromFirstToSecondOrigin.GetNorm() << "mm"; toleratedOriginError = fromFirstToSecondOrigin.GetNorm() * 0.3; // a third of the slice distance // (less than half, which would mean that a slice is displayed where another slice should actually be) } else { toleratedOriginError = m_ToleratedOriginOffset; } MITK_DEBUG << "Accepting errors in actual versus expected origin up to " << toleratedOriginError << "mm"; // Here we calculate if this slice and the previous one are well aligned, // i.e. we test if the previous origin is on a line through the current // origin, directed into the normal direction of the current slice. // If this is NOT the case, then we have a data set with a TILTED GANTRY geometry, // which cannot be simply loaded into a single mitk::Image at the moment. // For this case, we flag this finding in the result and DicomSeriesReader // can correct for that later. Vector3D right; right.Fill(0.0); Vector3D up; right.Fill(0.0); // might be down as well, but it is just a name at this point std::string orientationValue = (*dsIter)->GetTagValueAsString( tagImageOrientation ).value; DICOMStringToOrientationVectors( orientationValue, right, up, ignoredConversionError ); GantryTiltInformation tiltInfo( lastDifferentOrigin, thisOrigin, right, up, 1 ); if ( tiltInfo.IsSheared() ) { /* optimistic approach, accepting gantry tilt: save file for later, check all further files */ // at this point we have TWO slices analyzed! if they are the only two files, we still split, because there is no third to verify our tilting assumption. // later with a third being available, we must check if the initial tilting vector is still valid. if yes, continue. - // if NO, we need to split the already sorted part (result.first) and the currently analyzed file (*dsIter) + // if NO, we need to split the already sorted part (result->first) and the currently analyzed file (*dsIter) // tell apart gantry tilt from overall skewedness // sort out irregularly sheared slices, that IS NOT tilting if ( groupImagesWithGantryTilt && tiltInfo.IsRegularGantryTilt() ) { assert(!datasets.empty()); - result.FlagGantryTilt(tiltInfo); - result.AddFileToSortedBlock( *dsIter ); // this file is good for current block - result.SetFirstFilenameOfBlock( datasets.front()->GetFilenameIfAvailable() ); - result.SetLastFilenameOfBlock( datasets.back()->GetFilenameIfAvailable() ); + result->FlagGantryTilt(tiltInfo); + result->AddFileToSortedBlock( *dsIter ); // this file is good for current block + result->SetFirstFilenameOfBlock( datasets.front()->GetFilenameIfAvailable() ); + result->SetLastFilenameOfBlock( datasets.back()->GetFilenameIfAvailable() ); fileFitsIntoPattern = true; } else // caller does not want tilt compensation OR shearing is more complicated than tilt { - result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis + result->AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis + result->GetSplitReason()->AddReason(DICOMSplitReason::ReasonType::GantryTiltDifference); fileFitsIntoPattern = false; } } else // not sheared { - result.AddFileToSortedBlock( *dsIter ); // this file is good for current block + result->AddFileToSortedBlock( *dsIter ); // this file is good for current block fileFitsIntoPattern = true; } } else if (fromFirstToSecondOriginInitialized) // we already know the offset between slices { Point3D assumedOrigin = lastDifferentOrigin + fromFirstToSecondOrigin; Vector3D originError = assumedOrigin - thisOrigin; double norm = originError.GetNorm(); if (norm > toleratedOriginError) { MITK_DEBUG << " File does not fit into the inter-slice distance pattern (diff = " << norm << ", allowed " << toleratedOriginError << ")."; MITK_DEBUG << " Expected position (" << assumedOrigin[0] << "," << assumedOrigin[1] << "," << assumedOrigin[2] << "), got position (" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")"; MITK_DEBUG << " ==> Sort away " << *dsIter << " for later analysis"; // At this point we know we deviated from the expectation of ITK's ImageSeriesReader // We split the input file list at this point, i.e. all files up to this one (excluding it) // are returned as group 1, the remaining files (including the faulty one) are group 2 /* Optimistic approach: check if any of the remaining slices fits in */ - result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis + result->AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis + + const auto fromLastToThisOriginDistance = (lastDifferentOrigin - thisOrigin).GetNorm(); + const auto fromFirstToSecondOriginDistance = fromFirstToSecondOrigin.GetNorm(); + + auto currentMissCount = static_cast(std::round((fromLastToThisOriginDistance / fromFirstToSecondOriginDistance)-1)); + if (missingSlicesCount == 0 || missingSlicesCount > currentMissCount) + { + missingSlicesCount = currentMissCount; + } + + result->GetSplitReason()->AddReason(DICOMSplitReason::ReasonType::SliceDistanceInconsistency, "missing slices: " + std::to_string(missingSlicesCount)); fileFitsIntoPattern = false; } else { - result.AddFileToSortedBlock( *dsIter ); // this file is good for current block + result->AddFileToSortedBlock( *dsIter ); // this file is good for current block + result->GetSplitReason()->RemoveReason(DICOMSplitReason::ReasonType::SliceDistanceInconsistency); fileFitsIntoPattern = true; + missingSlicesCount = 0; } } else // this should be the very first slice { - result.AddFileToSortedBlock( *dsIter ); // this file is good for current block + result->AddFileToSortedBlock( *dsIter ); // this file is good for current block fileFitsIntoPattern = true; } } // record current origin for reference in later iterations if ( !lastOriginInitialized || ( fileFitsIntoPattern && (thisOrigin != lastOrigin) ) ) { lastDifferentOrigin = thisOrigin; } lastOrigin = thisOrigin; lastOriginInitialized = true; } - if ( result.ContainsGantryTilt() ) + if ( result->ContainsGantryTilt() ) { // check here how many files were grouped. // IF it was only two files AND we assume tiltedness (e.g. save "distance") // THEN we would want to also split the two previous files (simple) because // we don't have any reason to assume they belong together // Above behavior can be configured via m_AcceptTwoSlicesGroups, the default being "do accept" - if ( result.GetBlockDatasets().size() == 2 && !m_AcceptTwoSlicesGroups ) + if ( result->GetBlockDatasets().size() == 2 && !m_AcceptTwoSlicesGroups ) { - result.UndoPrematureGrouping(); + result->UndoPrematureGrouping(); + result->GetSplitReason()->AddReason(DICOMSplitReason::ReasonType::GantryTiltDifference); } } // update tilt info to get maximum precision // earlier, tilt was only calculated from first and second slice. // now that we know the whole range, we can re-calculate using the very first and last slice - if ( result.ContainsGantryTilt() && result.GetBlockDatasets().size() > 1 ) + if ( result->ContainsGantryTilt() && result->GetBlockDatasets().size() > 1 ) { try { - DICOMDatasetList datasets = result.GetBlockDatasets(); + DICOMDatasetList datasets = result->GetBlockDatasets(); DICOMDatasetAccess* firstDataset = datasets.front(); DICOMDatasetAccess* lastDataset = datasets.back(); unsigned int numberOfSlicesApart = datasets.size() - 1; std::string orientationString = firstDataset->GetTagValueAsString( tagImageOrientation ).value; std::string firstOriginString = firstDataset->GetTagValueAsString( tagImagePositionPatient ).value; std::string lastOriginString = lastDataset->GetTagValueAsString( tagImagePositionPatient ).value; - result.FlagGantryTilt( GantryTiltInformation::MakeFromTagValues( firstOriginString, lastOriginString, orientationString, numberOfSlicesApart )); + result->FlagGantryTilt( GantryTiltInformation::MakeFromTagValues( firstOriginString, lastOriginString, orientationString, numberOfSlicesApart )); } catch (...) { // just do not flag anything, we are ok } } return result; } diff --git a/Modules/DICOM/src/mitkThreeDnTDICOMSeriesReader.cpp b/Modules/DICOM/src/mitkThreeDnTDICOMSeriesReader.cpp index 599719fd36..0119f4a016 100644 --- a/Modules/DICOM/src/mitkThreeDnTDICOMSeriesReader.cpp +++ b/Modules/DICOM/src/mitkThreeDnTDICOMSeriesReader.cpp @@ -1,295 +1,305 @@ /*============================================================================ The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center (DKFZ) All rights reserved. Use of this source code is governed by a 3-clause BSD license that can be found in the LICENSE file. ============================================================================*/ #include "mitkThreeDnTDICOMSeriesReader.h" #include "mitkITKDICOMSeriesReaderHelper.h" mitk::ThreeDnTDICOMSeriesReader ::ThreeDnTDICOMSeriesReader(unsigned int decimalPlacesForOrientation) :DICOMITKSeriesGDCMReader(decimalPlacesForOrientation) ,m_Group3DandT(m_DefaultGroup3DandT), m_OnlyCondenseSameSeries(m_DefaultOnlyCondenseSameSeries) { } mitk::ThreeDnTDICOMSeriesReader ::ThreeDnTDICOMSeriesReader(const ThreeDnTDICOMSeriesReader& other ) :DICOMITKSeriesGDCMReader(other) ,m_Group3DandT(m_DefaultGroup3DandT), m_OnlyCondenseSameSeries(m_DefaultOnlyCondenseSameSeries) { } mitk::ThreeDnTDICOMSeriesReader ::~ThreeDnTDICOMSeriesReader() { } mitk::ThreeDnTDICOMSeriesReader& mitk::ThreeDnTDICOMSeriesReader ::operator=(const ThreeDnTDICOMSeriesReader& other) { if (this != &other) { DICOMITKSeriesGDCMReader::operator=(other); this->m_Group3DandT = other.m_Group3DandT; } return *this; } bool mitk::ThreeDnTDICOMSeriesReader ::operator==(const DICOMFileReader& other) const { if (const auto* otherSelf = dynamic_cast(&other)) { return DICOMITKSeriesGDCMReader::operator==(other) && this->m_Group3DandT == otherSelf->m_Group3DandT; } else { return false; } } void mitk::ThreeDnTDICOMSeriesReader ::SetGroup3DandT(bool on) { m_Group3DandT = on; } bool mitk::ThreeDnTDICOMSeriesReader ::GetGroup3DandT() const { return m_Group3DandT; } /** Helper function to make the code in mitk::ThreeDnTDICOMSeriesReader ::Condense3DBlocks(SortingBlockList& resultOf3DGrouping) more readable.*/ bool BlockShouldBeCondensed(bool onlyCondenseSameSeries, unsigned int currentBlockNumberOfSlices, unsigned int otherBlockNumberOfSlices, const mitk::DICOMDatasetFinding& currentBlockFirstOrigin, const mitk::DICOMDatasetFinding& currentBlockLastOrigin, const mitk::DICOMDatasetFinding& otherBlockFirstOrigin, const mitk::DICOMDatasetFinding& otherBlockLastOrigin, const mitk::DICOMDatasetFinding& currentBlockSeriesInstanceUID, const mitk::DICOMDatasetFinding& otherBlockSeriesInstanceUID) { if (otherBlockNumberOfSlices != currentBlockNumberOfSlices) return false; //don't condense blocks that have unequal slice count if (!otherBlockFirstOrigin.isValid || !otherBlockLastOrigin.isValid) return false; //don't condense blocks that have invalid origins if (!currentBlockFirstOrigin.isValid || !currentBlockLastOrigin.isValid) return false; //don't condense blocks that have invalid origins const bool sameSeries = otherBlockSeriesInstanceUID.isValid && currentBlockSeriesInstanceUID.isValid && otherBlockSeriesInstanceUID.value == currentBlockSeriesInstanceUID.value; if (onlyCondenseSameSeries && !sameSeries) return false; //don't condense blocks if it is only allowed to condense same series and series are not defined or not equal. if (otherBlockFirstOrigin.value != currentBlockFirstOrigin.value) return false; //don't condense blocks that have unequal first origins if (otherBlockLastOrigin.value != currentBlockLastOrigin.value) return false; //don't condense blocks that have unequal last origins return true; } mitk::DICOMITKSeriesGDCMReader::SortingBlockList mitk::ThreeDnTDICOMSeriesReader ::Condense3DBlocks(SortingBlockList& resultOf3DGrouping) { if (!m_Group3DandT) { return resultOf3DGrouping; // don't work if nobody asks us to } SortingBlockList remainingBlocks = resultOf3DGrouping; SortingBlockList non3DnTBlocks; SortingBlockList true3DnTBlocks; std::vector true3DnTBlocksTimeStepCount; // we should describe our need for this tag as needed via a function // (however, we currently know that the superclass will always need this tag) const DICOMTag tagImagePositionPatient(0x0020, 0x0032); const DICOMTag tagSeriesInstaceUID(0x0020, 0x000e); while (!remainingBlocks.empty()) { // new block to fill up - const DICOMDatasetAccessingImageFrameList& firstBlock = remainingBlocks.front(); + const DICOMDatasetAccessingImageFrameList& firstBlock = remainingBlocks.front().first; DICOMDatasetAccessingImageFrameList current3DnTBlock = firstBlock; + auto currentSplitReason = remainingBlocks.front().second; + int current3DnTBlockNumberOfTimeSteps = 1; // get block characteristics of first block const unsigned int currentBlockNumberOfSlices = firstBlock.size(); const auto currentBlockFirstOrigin = firstBlock.front()->GetTagValueAsString( tagImagePositionPatient ); const auto currentBlockLastOrigin = firstBlock.back()->GetTagValueAsString( tagImagePositionPatient ); const auto currentBlockSeriesInstanceUID = firstBlock.back()->GetTagValueAsString(tagSeriesInstaceUID); remainingBlocks.erase( remainingBlocks.begin() ); // compare all other blocks against the first one for (auto otherBlockIter = remainingBlocks.begin(); otherBlockIter != remainingBlocks.cend(); /*++otherBlockIter*/) // <-- inside loop { // get block characteristics from first block - const DICOMDatasetAccessingImageFrameList otherBlock = *otherBlockIter; + const DICOMDatasetAccessingImageFrameList otherBlock = otherBlockIter->first; const unsigned int otherBlockNumberOfSlices = otherBlock.size(); const auto otherBlockFirstOrigin = otherBlock.front()->GetTagValueAsString( tagImagePositionPatient ); const auto otherBlockLastOrigin = otherBlock.back()->GetTagValueAsString( tagImagePositionPatient ); const auto otherBlockSeriesInstanceUID = otherBlock.back()->GetTagValueAsString(tagSeriesInstaceUID); // add matching blocks to current3DnTBlock // keep other blocks for later if ( BlockShouldBeCondensed(m_OnlyCondenseSameSeries, currentBlockNumberOfSlices, otherBlockNumberOfSlices, currentBlockFirstOrigin, currentBlockLastOrigin, otherBlockFirstOrigin, otherBlockLastOrigin, currentBlockSeriesInstanceUID, otherBlockSeriesInstanceUID)) { // matching block ++current3DnTBlockNumberOfTimeSteps; current3DnTBlock.insert( current3DnTBlock.end(), otherBlock.begin(), otherBlock.end() ); // append // remove this block from remainingBlocks otherBlockIter = remainingBlocks.erase(otherBlockIter); // make sure iterator otherBlockIter is valid afterwards } else { ++otherBlockIter; } } // in any case, we now know all about the first block of our list ... // ... and we either call it 3D o 3D+t if (current3DnTBlockNumberOfTimeSteps > 1) { - true3DnTBlocks.push_back(current3DnTBlock); + true3DnTBlocks.push_back(std::make_pair(current3DnTBlock,currentSplitReason)); true3DnTBlocksTimeStepCount.push_back(current3DnTBlockNumberOfTimeSteps); } else { - non3DnTBlocks.push_back(current3DnTBlock); + non3DnTBlocks.push_back(std::make_pair(current3DnTBlock, currentSplitReason)); } } // create output for real 3D+t blocks (other outputs will be created by superclass) // set 3D+t flag on output block this->SetNumberOfOutputs( true3DnTBlocks.size() ); unsigned int o = 0; for (auto blockIter = true3DnTBlocks.cbegin(); blockIter != true3DnTBlocks.cend(); ++o, ++blockIter) { // bad copy&paste code from DICOMITKSeriesGDCMReader, should be handled in a better way - DICOMDatasetAccessingImageFrameList gdcmFrameInfoList = *blockIter; + const auto& gdcmFrameInfoList = blockIter->first; assert(!gdcmFrameInfoList.empty()); // reverse frames if necessary // update tilt information from absolute last sorting const DICOMDatasetList datasetList = ConvertToDICOMDatasetList( gdcmFrameInfoList ); m_NormalDirectionConsistencySorter->SetInput( datasetList ); m_NormalDirectionConsistencySorter->Sort(); const DICOMDatasetAccessingImageFrameList sortedGdcmInfoFrameList = ConvertToDICOMDatasetAccessingImageFrameList( m_NormalDirectionConsistencySorter->GetOutput(0) ); const GantryTiltInformation& tiltInfo = m_NormalDirectionConsistencySorter->GetTiltInformation(); // set frame list for current block const DICOMImageFrameList frameList = ConvertToDICOMImageFrameList( sortedGdcmInfoFrameList ); assert(!frameList.empty()); DICOMImageBlockDescriptor block; block.SetTagCache( this->GetTagCache() ); // important: this must be before SetImageFrameList(), because SetImageFrameList will trigger reading of lots of interesting tags! block.SetAdditionalTagsOfInterest(GetAdditionalTagsOfInterest()); block.SetTagLookupTableToPropertyFunctor(GetTagLookupTableToPropertyFunctor()); block.SetImageFrameList( frameList ); block.SetTiltInformation( tiltInfo ); + block.SetSplitReason(blockIter->second->Clone()); + + if (true3DnTBlocks.size() == 1 && non3DnTBlocks.empty()) + { + //if we have condensed everything into just on 3DnT block, we can remove the overlap reason, + //because no real overlap is existent any more. + block.GetSplitReason()->RemoveReason(DICOMSplitReason::ReasonType::OverlappingSlices); + } block.SetFlag("3D+t", true); block.SetIntProperty("timesteps", true3DnTBlocksTimeStepCount[o]); MITK_DEBUG << "Found " << true3DnTBlocksTimeStepCount[o] << " timesteps"; this->SetOutput( o, block ); } return non3DnTBlocks; } bool mitk::ThreeDnTDICOMSeriesReader ::LoadImages() { bool success = true; unsigned int numberOfOutputs = this->GetNumberOfOutputs(); for (unsigned int o = 0; o < numberOfOutputs; ++o) { const DICOMImageBlockDescriptor& block = this->InternalGetOutput(o); if (block.GetFlag("3D+t", false)) { success &= this->LoadMitkImageForOutput(o); } else { success &= DICOMITKSeriesGDCMReader::LoadMitkImageForOutput(o); // let superclass handle non-3D+t } } return success; } bool mitk::ThreeDnTDICOMSeriesReader ::LoadMitkImageForImageBlockDescriptor(DICOMImageBlockDescriptor& block) const { PushLocale(); const DICOMImageFrameList& frames = block.GetImageFrameList(); const GantryTiltInformation tiltInfo = block.GetTiltInformation(); const bool hasTilt = tiltInfo.IsRegularGantryTilt(); const int numberOfTimesteps = block.GetNumberOfTimeSteps(); if (numberOfTimesteps == 1) { return DICOMITKSeriesGDCMReader::LoadMitkImageForImageBlockDescriptor(block); } const int numberOfFramesPerTimestep = block.GetNumberOfFramesPerTimeStep(); ITKDICOMSeriesReaderHelper::StringContainerList filenamesPerTimestep; for (int timeStep = 0; timeStepFilename ); } filenamesPerTimestep.push_back( filenamesOfThisTimeStep ); } mitk::ITKDICOMSeriesReaderHelper helper; mitk::Image::Pointer mitkImage = helper.Load3DnT( filenamesPerTimestep, m_FixTiltByShearing && hasTilt, tiltInfo ); block.SetMitkImage( mitkImage ); PopLocale(); return true; }