diff --git a/include/dcmqi/ImageSEGConverter.h b/include/dcmqi/ImageSEGConverter.h index 14e3c52d..f4735916 100644 --- a/include/dcmqi/ImageSEGConverter.h +++ b/include/dcmqi/ImageSEGConverter.h @@ -18,9 +18,7 @@ #include // DCMQI includes -//#include "dcmqi/OverlapUtil.h" #include "OverlapUtil.h" -#include "SegmentAttributes.h" #include "dcmqi/ConverterBase.h" #include "dcmqi/OverlapUtil.h" #include "dcmqi/JSONSegmentationMetaInformationHandler.h" @@ -94,33 +92,73 @@ namespace dcmqi { OFCondition getNonOverlappingSegmentGroups(const bool mergeSegments, OverlapUtil::SegmentGroups& segmentGroups); + /** + * Extract basic image info like directions, origin, spacing and image region. + * @return EC_Normal if successful, error otherwise + */ OFCondition extractBasicSegmentationInfo(); + /** + * Allocate an ITK image template with the same size and spacing as the DICOM image. + * This is used as a template for each slice that is to be inserted into the ITK space. + * @return EC_Normal if successful, error otherwise + */ ShortImageType::Pointer allocateITKImageTemplate(); + /** + * Allocate an ITK image based on the given template. + * @param imageTemplate The template to use for the new image. + * @return EC_Normal if successful, error otherwise + */ ShortImageType::Pointer allocateITKImageDuplicate(ShortImageType::Pointer imageTemplate); + /** + * Get the ITK image origin for the given frame. + * @param frameNo The (DICOM) frame number to get the origin for. + * @param origin The resulting origin. + * @return EC_Normal if successful, error otherwise + */ OFCondition getITKImageOrigin(const Uint32 frameNo, ShortImageType::PointType& origin); + /** + * Collect segment metadata for the given frame that will later go into the + * accompanying JSON segmentation description. + * @param segmentGroup The segment group number + * (i.e. uniquely identifying the NRRD output file in the end). + * @param segmentNumber The DICOM segment number. + * @param origin The resulting origin. + * @return EC_Normal if successful, error otherwise + */ OFCondition addSegmentMetadata(const size_t segmentGroup, const Uint16 segmentNumber); + /// The segmentation object, guarded by unique pointer OFunique_ptr m_segDoc; + /// Image direction in ITK speak ShortImageType::DirectionType m_direction; - // Spacing and origin + /// Computed image slice spacing double m_computedSliceSpacing; + /// Computed volume extent double m_computedVolumeExtent; + /// Slice direction in ITK speak vnl_vector m_sliceDirection; + /// Image origin in ITK speak ShortImageType::PointType m_imageOrigin; + /// Image slice spacing in ITK speak ShortImageType::SpacingType m_imageSpacing; + /// Image size in ITK speak ShortImageType::SizeType m_imageSize; + /// Image region in ITK speak ShortImageType::RegionType m_imageRegion; + /// Segment meta information handler for the JSON segmentation description JSONSegmentationMetaInformationHandler m_metaInfo; + /// OverlapUtil instance used by this class, used in DICOM segmentation + /// to itk conversion OverlapUtil m_overlapUtil; }; diff --git a/include/dcmqi/OverlapUtil.h b/include/dcmqi/OverlapUtil.h index ffc029f9..f5bcb901 100644 --- a/include/dcmqi/OverlapUtil.h +++ b/include/dcmqi/OverlapUtil.h @@ -1,6 +1,6 @@ /* * - * Copyright (C) 1994-2023, Open Connections GmbH + * Copyright (C) 2023, Open Connections GmbH * All rights reserved. See COPYRIGHT file for details. * * This software and supporting documentation were developed by @@ -11,7 +11,7 @@ * D-26121 Oldenburg, Germany * * - * Module: dcmseg + * Module: dcmqi * * Author: Michael Onken * @@ -19,14 +19,13 @@ * */ -#ifndef SEGOVLAP_H -#define SEGOVLAP_H +#ifndef DCMQI_OVERLAPUTIL_H +#define DCMQI_OVERLAPUTIL_H -#include "dcmtk/config/osconfig.h" /* make sure OS specific configuration is included first */ +#include "dcmtk/dcmiod/iodtypes.h" +#include "dcmtk/ofstd/ofcond.h" #include "dcmtk/ofstd/oftypes.h" #include "dcmtk/ofstd/ofvector.h" -#include "dcmtk/ofstd/ofcond.h" -#include "dcmtk/dcmiod/iodtypes.h" #include class DcmSegmentation; @@ -34,267 +33,328 @@ class DcmSegmentation; /// error: frames are not parallel extern const OFConditionConst SG_EC_FramesNotParallel; -namespace dcmqi { +namespace dcmqi +{ +/** Class that analyzes the frame-based segment structures of a given segmentation object. + * It provides the following main functionality: + * - Grouping of "physical" frames by their position in space (called "logical frames"), i.e. + * frames that are at the same position in space are grouped together. + * - For every logical frame (i.e. position in space), it lists the segments found at this + * position together with their respective physical frame number + * - Return physical frame numbers for a given segment number + * - Building an overlap matrix that stores for each arbitrary segment pair whether + * they overlap or not. + * - Return groups of segments, that no two overlapping segments will be in the same group. + * This will not necessarily return the optimal solution, but a solution that should be good enough. + * For the method used see getNonOverlappingSegments(). + */ class OverlapUtil { public: + /// Image Position Patient tuple (x,y,z) + typedef OFVector ImagePosition; // might be defined in respective functional group - // Image Position Patient tuple (x,y,z) - typedef OFVector ImagePosition; // might be defined in respective functional group - - struct FramePositionAndNumber - { - /** Default constructor required for vector initialization + /// Physical Frame number with its respective position + struct FramePositionAndNumber + { + /** Default constructor required for vector initialization + */ + FramePositionAndNumber() + : m_position() + , m_frameNumber(0) + { + } + /** Constructor + * @param pos Image position + * @param num Pyhsical frame number + */ + FramePositionAndNumber(const ImagePosition& pos, const Uint32& num) + : m_position(pos) + , m_frameNumber(num) + { + } + /// Frame position in space + ImagePosition m_position; + /// Physical frame number (number of frame in DICOM object) + Uint32 m_frameNumber; + }; + + /// Physical Frame number with its respective position + typedef OFVector FramePositions; + + /// Logical Frame, represented and defined by various physical frames (numbers) at the same position + typedef OFVector LogicalFrame; + + /// All distinct positions and for each position the physical frame numbers that can be found at it + typedef OFVector DistinctFramePositions; + + /// Lists frames for each segment where segment with index i is represented by the vector at index i, + /// and index 0 is unused. I.e. index i is segment number, value is vector of physical frame numbers. + typedef OFVector> FramesForSegment; + + /// Implements comparision operator to be used for sorting of frame positions, + /// making the sorting order depend on the coordinate given in the constructor + struct ComparePositions + { + /** Constructor, used to configure coordinate position to be used for sorting + * @param c Coordinate position to be used for sorting + */ + ComparePositions(size_t c) + : m_coordinate(c) + { + } + + /** Compare two frames*/ + bool operator()(const FramePositionAndNumber& a, const FramePositionAndNumber& b) const + { + return a.m_position[m_coordinate] < b.m_position[m_coordinate]; + } + /// Coordinate position (0-2, i.e. x,x,z) to be used for sorting + size_t m_coordinate; + }; + + /// Matrix of N x N segment numbers, where N is the number of segments. + /// Value is 1 at x,y if x and y overlap, 0 if they don't overlap, and -1 if not initialized. + typedef OFVector> OverlapMatrix; + + /// Group of non-overlapping segments (each represented by its segment number) + typedef OFVector> SegmentGroups; + + /** Represents a segment number and a logical frame number it is found at */ - FramePositionAndNumber() : m_position(), m_frameNumber(0) {} - /** Constructor - * @param pos Image position - * @param num Pyhsical frame number - */ - FramePositionAndNumber(const ImagePosition& pos, const Uint32& num) : m_position(pos), m_frameNumber(num) {} - /// Frame position in space - ImagePosition m_position; - /// Physical frame number (number of frame in DICOM object) - Uint32 m_frameNumber; - }; - - /// Physical Frame number with its respective position - typedef OFVector FramePositions; - - /// Logical Frame, represented by various physical frames (numbers) at the same position - typedef OFVector LogicalFrame; - - /// All distinct positions and for each position the physical frame numbers that can be found at it - typedef OFVector DistinctFramePositions; - - /// Lookup position of a frame in DistinctFramePositions vector by providing the frame number, - /// i.e. index is frame number, value is position in DistinctFramePositions vector; - /// or: which physical frame number is at which logical frame number/position? - typedef OFVector LogicalFrameToPositionIndex; - - /// Lists frames for each segment where segment with index i is represented by the vector at index i, - /// and index 0 is unused. I.e. index i is segment number, value is vector of physical frame numbers. - typedef OFVector > FramesForSegment; - - // Implements comparision operator to be used for sorting of frame positions, - // making the sorting order depend on the coordinate given in the constructor - struct ComparePositions - { - /** Constructor, used to configure coordinate position to be used for sorting - * @param c Coordinate position to be used for sorting - */ - ComparePositions(size_t c) : m_coordinate(c) { } - - /** Compare two frames*/ - bool operator()(const FramePositionAndNumber& a, const FramePositionAndNumber& b) const - { - return a.m_position[m_coordinate] < b.m_position[m_coordinate]; - } - /// Coordinate position (0-2, i.e. x-z) to be used for sorting - size_t m_coordinate; - }; - - typedef OFVector > OverlapMatrix; - // TODO: Check whether a vector would also do - typedef OFVector< std::set > SegmentGroups; - - - /** Represents a segment number and a logical frame number it is found at - */ - struct SegNumAndFrameNum - { - /** Constructor - * @param s Segment number - * @param f Logical frame number + struct SegNumAndFrameNum + { + /** Constructor + * @param s Segment number + * @param f Logical frame number + */ + SegNumAndFrameNum(const Uint16& s, const Uint16 f) + : m_segmentNumber(s) + , m_frameNumber(f) + { + } + /// Segment number as used in DICOM segmentation object (1-n) + Uint16 m_segmentNumber; + /// Logical frame number (number of frame in DistinctFramePositions vector) + Uint16 m_frameNumber; + /** Comparison operator + * @param rhs Right-hand side of comparison + * @return OFTrue if left-hand side is smaller than right-hand side + */ + bool operator<(const SegNumAndFrameNum& rhs) const + { + return m_segmentNumber < rhs.m_segmentNumber; + } + }; + + /// Segments and their phyiscal frame number (inner set), grouped by their + /// respective logical frame number (outer vector) + typedef OFVector> SegmentsByPosition; + + // ------------------------------------------ Methods ------------------------------------------ + + /** Constructor. Use setSegmentationObject() to set the segmentation object to work with. */ - SegNumAndFrameNum(const Uint16& s, const Uint16 f) : m_segmentNumber(s), m_frameNumber(f) {} - /// Segment number as used in DICOM segmentation object (1-n) - Uint16 m_segmentNumber; - /// Logical frame number (number of frame in DistinctFramePositions vector) - Uint16 m_frameNumber; - /** Comparison operator - * @param rhs Right-hand side of comparison - * @return OFTrue if left-hand side is smaller than right-hand side + OverlapUtil(); + + /** Destructor */ + ~OverlapUtil(); + + /** Set the segmentation object to work with and clears all old data. + * TODO: In the future, maybe have DcmSegmentation->getOverlapUtil() to access + * the OverlapUtil object, so that the user does not have to care about + * feeding the segmentation object to the OverlapUtil object. */ - bool operator<(const SegNumAndFrameNum& rhs) const - { - return m_segmentNumber < rhs.m_segmentNumber; - } - }; - - /// Segments and their phyiscal frame number (inner set), grouped by their respective logical frame number (outer vector) - typedef OFVector > SegmentsByPosition; - - // ------------------------------------------ Methods ------------------------------------------ - - - /** Constructor. Use setSegmentationObject() to set the segmentation object to work with. - */ - OverlapUtil(); - - /** Destructor */ - ~OverlapUtil(); - - /** Set the segmentation object to work with - * TODO: In the future, maybe have DcmSegmentation->getOverlapUtil() to access - * the OverlapUtil object, so that the user does not have to care about - * feeding the segmentation object to the OverlapUtil object. - */ - void setSegmentationObject(DcmSegmentation* seg); - - /** Clears all internal data (except segmentation object reference) - */ - void clear(); - - /** Get all distinct frame positions and the physical frame numbers at this position - * @param result Resulting vector of distinct frame positions - * @return EC_Normal if successful, error otherwise - */ - OFCondition getFramesByPosition(DistinctFramePositions& result); - - /** Get all segments and their phyisical frame number, grouped by their respective logical frame number - * @param result Resulting vector of segments grouped by logical frame number - * @return EC_Normal if successful, error otherwise - */ - OFCondition getSegmentsByPosition(SegmentsByPosition& result) ; - - /** TODO Not yet implemented, nice extra functionality */ - OFCondition getSegmentsForFrame(const Uint32 frameNumber, OFVector& segments); - - /** Get phyiscal frames for a specific segment - * @param segmentNumber Segment number to get frames for (1..n) - * @param frames Resulting vector of physical frame numbers - * (1 is first used index, 0 is first frame number) - * @return EC_Normal if successful, error otherwise - */ - OFCondition getFramesForSegment(const Uint32 segmentNumber, OFVector& frames); - - /** Returns computed overlap matrix - * @param matrix Resulting overlap matrix - * @return EC_Normal if successful, error otherwise - */ - OFCondition getOverlapMatrix(OverlapMatrix& matrix); - - /** Returns segments grouped together in a way, that no two overlapping - * segments will be in the same group. This method does not necessarily - * returns the optimal solution, but a solution that should be good enough. - * It is guaranteed, that segments in the same group don't overlap. - * - * It is based on the idea of a greedy algorithm that creates a first group - * containing the first segment. Then it goes to the next segment, checks whether - * it fits into the first group with no overlaps (easily checked in overlap matrix) - * and inserts it into that group if no overlaps exists. Otherwise, it creates a - * new group and continues with the next segment (trying to insert it into - * the first group, then second group, otherwise creates third group, and so on). - * @param segmentGroups Resulting vector of segment groups, each listing the - * segment numbers that are in that group - * @return EC_Normal if successful, error otherwise - */ - OFCondition getNonOverlappingSegments(SegmentGroups& segmentGroups); - - /** Prints overlap matrix to given stream - * @param ss The stream to dump to - */ - void printOverlapMatrix(OFStringStream& ss); - - /** Prints groups of non-overlapping segments to given stream - * @param ss The stream to dump to - */ - void printNonOverlappingSegments(OFStringStream& ss); + void setSegmentationObject(DcmSegmentation* seg); -protected: + /** Clears all internal data (except segmentation object reference). + * This should be called whenever the input data (i.e. the underlying) + * DICOM segmentation object changes, before calling any other method. + */ + void clear(); - /** Builds the overlap matrix, if not already done. - * @return EC_Normal if successful or already existant, error otherwise - */ - OFCondition buildOverlapMatrix(); - - /** Checks if frames are parallel, i.e. if the image position patient - * of all frames are parallel to each other (i.e. found in the shared functional group) - * @return EC_Normal if parallel, SG_EC_FramesNotParallel if image orientation is not shared, - * error otherwise - */ - OFCondition ensureFramesAreParallel(); - - /** Groups all physical frames by their position. This also works if the physical frames - * have slightly differennt positions, i.e. if they are not exactly the same and are only - * "close enough" to be considered the same. Right now, the maximum distance treated equal - * is if distance is smaller than slice thickness * 0.01 (i.e. 1% of slice thickness). - * @return EC_Normal if successful, error otherwise - */ - OFCondition groupFramesByPosition(); - - /** Checks whether the given two frames overlap - * @param f1 Frame 1, provided by its physical frame number - * @param f2 Frame 2, provided by its physical frame number - * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) - * @return EC_Normal if successful, error otherwise - */ - OFCondition checkFramesOverlap(const Uint32& f1, const Uint32& f2, OFBool& overlap); - - /** Checks whether the given two frames overlap by using comparing their pixel data - * by bitwise "and". Thi should be pretty efficient, however, only works and is called (right now), - * if row*cols % 8 = 0, so we can easily extract frames as binary bitsets without unpacking them. - * @param f1 Frame 1, provided by its physical frame number - * @param f2 Frame 2, provided by its physical frame number - * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) - * @return EC_Normal if successful, error otherwise - */ - OFCondition checkFramesOverlapBinary(const Uint32& f1, const Uint32& f2, const DcmIODTypes::Frame* f1_data, const DcmIODTypes::Frame* f2_data, const Uint16& rows, const Uint16 cols, OFBool& overlap); - - /** Checks whether the given two frames overlap by using comparing their pixel data - * after unpacking, i.e. expanding every bit to a byte, and then comparing whether the two - * related bytes of each frame are both non-zero. This is less efficient than checkFramesOverlapBinary(), - * @param f1 Frame 1, provided by its physical frame number - * @param f2 Frame 2, provided by its physical frame number - * @param f1_data Pixel data of frame 1 - * @param f2_data Pixel data of frame 2 - * @param rows Number of rows of the frame(s) - * @param cols Number of columns of the frame(s) - * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) - * @return EC_Normal if successful, error otherwise - */ - OFCondition checkFramesOverlapUnpacked(const Uint32& f1, const Uint32& f2, const DcmIODTypes::Frame* f1_data, const DcmIODTypes::Frame* f2_data, const Uint16& rows, const Uint16 cols, OFBool& overlap); + /** Get all distinct frame positions and the physical frame numbers at this position + * @param result Resulting vector of distinct frame positions + * @return EC_Normal if successful, error otherwise + */ + OFCondition getFramesByPosition(DistinctFramePositions& result); -private: + /** Get all segments and their physical frame number, grouped by their respective logical frame number + * @param result Resulting vector of segments grouped by logical frame number + * @return EC_Normal if successful, error otherwise + */ + OFCondition getSegmentsByPosition(SegmentsByPosition& result); - /// Frames with their respective positions (IPP) - FramePositions m_FramePositions; + /** Get phyiscal frames for a specific segment + * @param segmentNumber Segment number to get frames for (1..n) + * @param frames Resulting vector of physical frame numbers + * (1 is first used index, 0 is first frame number) + * @return EC_Normal if successful, error otherwise + */ + OFCondition getFramesForSegment(const Uint32 segmentNumber, OFVector& frames); - /// Frame numbers (starting from 0) grouped by segment number (first segment is 1, - /// i.e. index 0 is unused) - FramesForSegment m_framesForSegment; + /** Returns computed overlap matrix + * @param matrix Resulting overlap matrix + * @return EC_Normal if successful, error otherwise + */ + OFCondition getOverlapMatrix(OverlapMatrix& matrix); + + /** Returns segments grouped together in a way, that no two overlapping + * segments will be in the same group. This method does not necessarily + * returns the optimal solution, but a solution that should be good enough. + * It is guaranteed, that segments in the same group don't overlap. + * + * It is based on the idea of a greedy algorithm that creates a first group + * containing the first segment. Then it goes to the next segment, checks whether + * it fits into the first group with no overlaps (easily checked in overlap matrix) + * and inserts it into that group if no overlaps exists. Otherwise, it creates a + * new group and continues with the next segment (trying to insert it into + * the first group, then second group, otherwise creates third group, and so on). + * @param segmentGroups Resulting vector of segment groups, each listing the + * segment numbers that are in that group + * @return EC_Normal if successful, error otherwise + */ + OFCondition getNonOverlappingSegments(SegmentGroups& segmentGroups); - /// Logical frames, ie. physical frames with the same position are - /// grouped together in a logical frame. For every logical frame, we - /// store the related physical frame numbers. The logical frame number - /// is implicitly given by the index in the vector. - DistinctFramePositions m_logicalFramePositions; + /** Prints segments by their position in space + * @param ss The stream to dump to + */ + void printSegmentsByPosition(OFStringStream& ss); - // TODO: Not needed? - // Allows reverse lookup of logical frames by position in m_logicalFramePositions - // LogicalFrameToPositionIndex m_frameNumToPositionIndex; + /** Prints segment overlap matrix to given stream + * @param ss The stream to dump to + */ + void printOverlapMatrix(OFStringStream& ss); - /// Stores for each logical frame a collection of (paired) segment and physical frame number, - /// that exists at that position. - SegmentsByPosition m_segmentsByPosition; + /** Prints groups of non-overlapping segments (identified by their numbers) + * to given stream + * @param ss The stream to dump to + */ + void printNonOverlappingSegments(OFStringStream& ss); - /// Matrix that stores for each segment pair whether they overlap or not. - /// I.e. Matrix has size n x n, where n is the number of segments. - /// The diagonal is always 0 (no overlap), i.e. a segment never overlaps with itself. - /// If there is an overlap, the value is 1. - /// If the field is not initialized, the value is -1. - OverlapMatrix m_segmentOverlapMatrix; +protected: + /** Collect all physical frame positions in m_FramePositions + * @return EC_Normal if successful, error otherwise + */ + OFCondition collectPhysicalFramePositions(); + + /** Group physical frame positions into logical positions. This is done by sorting + * frames after *that* position coordinate that in its mean position difference is + * larger than slice thickness * 0.9. Then those frames that are close enough to + * each other (i.e. distance is smaller than slice thickness * 0.01), end up at the + * same logical position (considered a "logical frame") + * TODO: This should probably not use mean values for the coordinates + * since in some cases, the mean difference in a slice coordinate might be close to 0 + * if many frames are at the same position. Instead, the maximum difference, variance or + * something else could be used? + * @return EC_Normal if successful, error otherwise + */ + OFCondition groupFramesByLogicalPosition(); - //// Groups of segments that do not overlap with each other - SegmentGroups m_nonOverlappingSegments; + /** Builds the overlap matrix, if not already done. + * @return EC_Normal if successful or already existant, error otherwise + */ + OFCondition buildOverlapMatrix(); + + /** Checks if frames are parallel, i.e. if DICOM Image Position Patient is present and + * all frames are parallel to each other (i.e. found in the shared functional group) + * @return EC_Normal if parallel, SG_EC_FramesNotParallel if image orientation is not shared, + * error otherwise + */ + OFCondition ensureFramesAreParallel(); + + /** Groups all physical frames by their position. This also works if the physical frames + * have slightly different positions, i.e. if they are not exactly the same and are only + * "close enough" to be considered the same. Right now, the maximum distance treated equal + * is if distance is smaller than slice thickness * 0.01 (i.e. 1% of slice thickness). + * Only performs the computation, if not done before. + * @return EC_Normal if successful, error otherwise + */ + OFCondition groupFramesByPosition(); - /// Reference to segmentation object to work with - DcmSegmentation* m_seg; + /** Checks whether the given two frames overlap + * @param f1 Frame 1, provided by its physical frame number + * @param f2 Frame 2, provided by its physical frame number + * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) + * @return EC_Normal if successful, error otherwise + */ + OFCondition checkFramesOverlap(const Uint32& f1, const Uint32& f2, OFBool& overlap); + + /** Checks whether the given two frames overlap by using comparing their pixel data + * by bitwise "and". This is very efficient, however, only works and is called (right now), + * if row*cols % 8 = 0, so we can easily extract frames as binary bitsets without unpacking them. + * TODO: Check whether this can be easily extended to other cases as well. + * @param f1 Frame 1, provided by its physical frame number + * @param f2 Frame 2, provided by its physical frame number + * @param f1_data Pixel data of frame 1 + * @param f2_data Pixel data of frame 2 + * @param rows Number of rows of the frame(s) + * @param cols Number of columns of the frame(s) + * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) + * @return EC_Normal if successful, error otherwise + */ + OFCondition checkFramesOverlapBinary(const Uint32& f1, + const Uint32& f2, + const DcmIODTypes::Frame* f1_data, + const DcmIODTypes::Frame* f2_data, + const Uint16& rows, + const Uint16 cols, + OFBool& overlap); + + /** Checks whether the given two frames overlap by using comparing their pixel data + * after unpacking, i.e. expanding every bit to a byte, and then comparing whether the two + * related bytes of each frame are both non-zero. This is less efficient than checkFramesOverlapBinary(), + * @param f1 Frame 1, provided by its physical frame number + * @param f2 Frame 2, provided by its physical frame number + * @param f1_data Pixel data of frame 1 + * @param f2_data Pixel data of frame 2 + * @param rows Number of rows of the frame(s) + * @param cols Number of columns of the frame(s) + * @param overlap Resulting overlap (overlaps if OFTrue, otherwise not) + * @return EC_Normal if successful, error otherwise + */ + OFCondition checkFramesOverlapUnpacked(const Uint32& f1, + const Uint32& f2, + const DcmIODTypes::Frame* f1_data, + const DcmIODTypes::Frame* f2_data, + const Uint16& rows, + const Uint16 cols, + OFBool& overlap); +private: + /// Phyiscal frames with their respective positions (IPP) + FramePositions m_framePositions; + + /// Frame numbers (starting from 0) grouped by segment number + /// (first segment number is 1, i.e. index 0 is unused) + FramesForSegment m_framesForSegment; + + /// Logical frames, ie. physical frames with the same position are + /// grouped together to a logical frame. For every logical frame, we + /// store the related physical frame numbers. The logical frame number + /// is implicitly given by the index in the vector. + DistinctFramePositions m_logicalFramePositions; + + /// Stores for each logical frame a collection of (paired) segment and + /// physical frame number, that exists at that position. + SegmentsByPosition m_segmentsByPosition; + + /// Matrix that stores for each segment pair whether they overlap or not. + /// I.e. Matrix has size N x N, where N is the number of segments. + /// The diagonal is always 0 (no overlap), i.e. a segment never overlaps with itself. + /// If there is an overlap, the value is 1. + /// If the field is not initialized, the value is -1. + OverlapMatrix m_segmentOverlapMatrix; + + //// Groups of segments that do not overlap with each other + SegmentGroups m_nonOverlappingSegments; + + /// Reference to segmentation object to work with + /// Must be freed outside this class. + DcmSegmentation* m_seg; }; } // namespace dcmqi -#endif // SEGOVLAP_H \ No newline at end of file +#endif // DCMQI_OVERLAPUTIL_H \ No newline at end of file diff --git a/libsrc/ImageSEGConverter.cpp b/libsrc/ImageSEGConverter.cpp index a3e062ef..170f1c27 100644 --- a/libsrc/ImageSEGConverter.cpp +++ b/libsrc/ImageSEGConverter.cpp @@ -27,9 +27,10 @@ namespace dcmqi { , m_metaInfo() , m_overlapUtil() { - }; + // ------------------------------------------------------------------------------------- + DcmDataset* ImageSEGConverter::itkimage2dcmSegmentation(vector dcmDatasets, vector segmentations, const string &metaData, @@ -415,7 +416,6 @@ namespace dcmqi { // clean up for the next frame fgder->clearData(); } - } } } @@ -527,7 +527,7 @@ namespace dcmqi { segment2image, m_metaInfo.getJSONOutputAsString()); } -// ------------------------------------------------------------------------------------- + // ------------------------------------------------------------------------------------- std::map ImageSEGConverter::dcmSegmentation2itkimage(const bool mergeSegments) { @@ -676,7 +676,6 @@ namespace dcmqi { // Otherwise, just copy the pixel value (fractional value) from the frame. if (m_segDoc->getSegmentationType() == DcmSegTypes::ST_BINARY) { - itkImage->SetPixel(index, *segNum); } else @@ -728,12 +727,11 @@ namespace dcmqi { m_metaInfo.setBodyPartExamined(bodyPartExamined.c_str()); } - // ------------------------------------------------------------------------------------- OFCondition ImageSEGConverter::extractBasicSegmentationInfo() { - // TODO: error handling + // TODO: Better error handling OFCondition result; // Directions @@ -808,8 +806,8 @@ namespace dcmqi { size_t numSegs = m_segDoc->getNumberOfSegments(); for (size_t i = 1; i <= numSegs; ++i) { - std::set segs; - segs.insert(i); + OFVector segs; + segs.push_back(i); segmentGroups.push_back(segs); } cout << "Will not merge segments: Splitting segments into " << segmentGroups.size() << " groups" << endl; @@ -843,6 +841,8 @@ namespace dcmqi { return newSegmentImage; } + // ------------------------------------------------------------------------------------- + OFCondition ImageSEGConverter::getITKImageOrigin(const Uint32 frameNo, ShortImageType::PointType& origin) { FGInterface& fgInterface = m_segDoc->getFunctionalGroups(); @@ -862,6 +862,7 @@ namespace dcmqi { return EC_Normal; } + // ------------------------------------------------------------------------------------- OFCondition ImageSEGConverter::addSegmentMetadata(const size_t segmentGroup, const Uint16 segmentNumber) diff --git a/libsrc/OverlapUtil.cpp b/libsrc/OverlapUtil.cpp index 49bd7663..bbae9339 100644 --- a/libsrc/OverlapUtil.cpp +++ b/libsrc/OverlapUtil.cpp @@ -1,16 +1,36 @@ -#include "dcmtk/config/osconfig.h" /* make sure OS specific configuration is included first */ +/* + * + * Copyright (C) 2023, Open Connections GmbH + * All rights reserved. See COPYRIGHT file for details. + * + * This software and supporting documentation were developed by + * + * OFFIS e.V. + * R&D Division Health + * Escherweg 2 + * D-26121 Oldenburg, Germany + * + * + * Module: dcmqi + * + * Author: Michael Onken + * + * Purpose: Interface of class OverlapUtil + * + */ + #include "dcmqi/OverlapUtil.h" #include "dcmtk/dcmdata/dcerror.h" +#include "dcmtk/dcmfg/fginterface.h" #include "dcmtk/dcmfg/fgpixmsr.h" +#include "dcmtk/dcmfg/fgplanor.h" +#include "dcmtk/dcmfg/fgplanpo.h" #include "dcmtk/dcmfg/fgseg.h" +#include "dcmtk/dcmfg/fgtypes.h" +#include "dcmtk/dcmiod/iodtypes.h" #include "dcmtk/dcmseg/segdoc.h" #include "dcmtk/dcmseg/segtypes.h" #include "dcmtk/dcmseg/segutils.h" -#include "dcmtk/dcmiod/iodtypes.h" -#include "dcmtk/dcmfg/fginterface.h" -#include "dcmtk/dcmfg/fgtypes.h" -#include "dcmtk/dcmfg/fgplanor.h" -#include "dcmtk/dcmfg/fgplanpo.h" #include "dcmtk/ofstd/ofcond.h" #include "dcmtk/ofstd/ofstream.h" #include "dcmtk/ofstd/oftimer.h" @@ -21,80 +41,81 @@ makeOFConditionConst(SG_EC_FramesNotParallel, OFM_dcmseg, 7, OF_error, "Frames are not parallel"); -namespace dcmqi { +namespace dcmqi +{ OverlapUtil::OverlapUtil() -: m_FramePositions() -, m_logicalFramePositions() -//, m_frameNumToPositionIndex() -, m_segmentsByPosition() -, m_segmentOverlapMatrix(0) -, m_nonOverlappingSegments() -, m_seg() + : m_framePositions() + , m_framesForSegment() + , m_logicalFramePositions() + , m_segmentsByPosition() + , m_segmentOverlapMatrix(0) + , m_nonOverlappingSegments() + , m_seg() { - } - OverlapUtil::~OverlapUtil() { // nothing to do } - void OverlapUtil::setSegmentationObject(DcmSegmentation* seg) { m_seg = seg; clear(); } - void OverlapUtil::clear() { - m_FramePositions.clear(); - m_logicalFramePositions.clear(); - m_segmentsByPosition.clear(); - m_segmentOverlapMatrix.clear(); - m_nonOverlappingSegments.clear(); - m_framesForSegment.clear(); - + m_framePositions.clear(); + m_framesForSegment.clear(); + m_logicalFramePositions.clear(); + m_segmentsByPosition.clear(); + m_segmentOverlapMatrix.clear(); + m_nonOverlappingSegments.clear(); } - OFCondition OverlapUtil::getFramesByPosition(DistinctFramePositions& result) { OFCondition cond = ensureFramesAreParallel(); if (cond.good()) { - cond = groupFramesByPosition(); + if (m_logicalFramePositions.empty()) + { + cond = groupFramesByPosition(); + } + } + if (cond.good()) + { + result = m_logicalFramePositions; } return cond; } - OFCondition OverlapUtil::getFramesForSegment(const Uint32 segmentNumber, OFVector& frames) { - if ( (segmentNumber == 0) || (segmentNumber > m_seg->getNumberOfSegments() + 1 ) ) + if ((segmentNumber == 0) || (segmentNumber > m_seg->getNumberOfSegments() + 1)) { DCMSEG_ERROR("getFramesForSegment(): Segment number " << segmentNumber << " is out of range"); return EC_IllegalParameter; } if (m_framesForSegment.empty()) { - FGInterface& fg = m_seg->getFunctionalGroups(); + FGInterface& fg = m_seg->getFunctionalGroups(); Uint32 numFrames = m_seg->getNumberOfFrames(); m_framesForSegment.resize(numFrames); // Get Segmentation FG for each frame and remember the segment number for each frame // in the vector m_segmentsForFrame for (size_t f = 0; f < numFrames; f++) { - FGBase* group = NULL; + FGBase* group = NULL; FGSegmentation* segFG = NULL; - group = fg.get(f, DcmFGTypes::EFG_SEGMENTATION); - segFG = OFstatic_cast(FGSegmentation*, group); + group = fg.get(f, DcmFGTypes::EFG_SEGMENTATION); + segFG = OFstatic_cast(FGSegmentation*, group); if (segFG) { - Uint16 segNum = 0; + Uint16 segNum = 0; OFCondition cond = segFG->getReferencedSegmentNumber(segNum); if (cond.good() && segNum > 0) { @@ -102,12 +123,15 @@ OFCondition OverlapUtil::getFramesForSegment(const Uint32 segmentNumber, OFVecto } else if (segNum == 0) { - DCMSEG_WARN("getSegmentsForFrame(): Referenced Segment Number is 0 (not permitted) for frame #" << f << ", ignoring"); + DCMSEG_WARN("getFramesForSegment(): Referenced Segment Number is 0 (not permitted) for frame #" + << f << ", ignoring"); return EC_InvalidValue; } else { - DCMSEG_ERROR("getSegmentsForFrame(): Referenced Segment Number not found (not permitted) for frame #" << f << ", cannot add segment"); + DCMSEG_ERROR( + "getFramesForSegment(): Referenced Segment Number not found (not permitted) for frame #" + << f << ", cannot add segment"); return EC_TagNotFound; } } @@ -117,12 +141,11 @@ OFCondition OverlapUtil::getFramesForSegment(const Uint32 segmentNumber, OFVecto return EC_Normal; } - OFCondition OverlapUtil::ensureFramesAreParallel() { FGInterface& fg = m_seg->getFunctionalGroups(); OFCondition cond; - OFBool perFrame = OFFalse; + OFBool perFrame = OFFalse; FGPlaneOrientationPatient* pop = NULL; // Ensure that Image Orientation Patient is shared, i.e. we have parallel frames OFVector iop(6); @@ -136,7 +159,8 @@ OFCondition OverlapUtil::ensureFramesAreParallel() } else { - DCMSEG_ERROR("getFramesByPosition(): Image Orientation Patient is per-frame, frames are probably not parallel"); + DCMSEG_ERROR( + "getFramesByPosition(): Image Orientation Patient is per-frame, frames are probably not parallel"); return SG_EC_FramesNotParallel; } } @@ -148,165 +172,25 @@ OFCondition OverlapUtil::ensureFramesAreParallel() return EC_Normal; } - OFCondition OverlapUtil::groupFramesByPosition() { - m_FramePositions.clear(); - m_logicalFramePositions.clear(); - //m_frameNumToPositionIndex.clear(); + if (!m_framePositions.empty()) + { + // Already computed + return EC_Normal; + } OFTimer tm; + // Group all frames by position into m_logicalFramePositions. // After that, all frames at the same position will be in the same vector // assigned to the same key (the frame's coordinates) in the map. - FGInterface& fg = m_seg->getFunctionalGroups(); - size_t numFrames = m_seg->getNumberOfFrames(); - OFBool perFrame = OFFalse; - OFCondition cond; - // Vector of frame numbers with their respective position - m_FramePositions.reserve(numFrames); - - // Put all frames into vector along with their Image Position Patient coordinates - for (size_t i = 0; i < numFrames; ++i) - { - FGPlanePosPatient* ppp = NULL; - FGBase* group = fg.get(i, DcmFGTypes::EFG_PLANEPOSPATIENT, perFrame); - if (group) ppp = OFstatic_cast(FGPlanePosPatient*, group); - if (ppp) - { - // Get image position patient for frame i - OFVector ipp(3); - // Only in later DCMTK version: - // cond = ppp->getImagePositionPatient(ipp); - cond = ppp->getImagePositionPatient(ipp[0], ipp[1], ipp[2]); - if (cond.good()) - { - // Insert frame into map - m_FramePositions.push_back( FramePositionAndNumber(ipp, i )); - //m_frameNumToPositionIndex.push_back(m_FramePositions.size() - 1); // current position in m_logicalFramePositions - } - else - { - DCMSEG_ERROR("groupFramesByPosition(): Image Position Patient not found for frame " << i << ", cannot sort frames by position"); - cond = EC_TagNotFound; - break; - } - } - else - { - DCMSEG_ERROR("groupFramesByPosition(): Image Position Patient not found for frame " << i << ", cannot sort frames by position"); - cond = EC_TagNotFound; - break; - } - } - // Find all distinct positions and for each position the actual frames that can be found at it + // Group all frames by position into m_logicalFramePositions. + OFCondition cond = collectPhysicalFramePositions(); if (cond.good()) { - // Get Slice Thickness - Float64 sliceThickness = 0.0; - FGPixelMeasures* pm = NULL; - Uint8 relevantCoordinate = 0; - FGBase* group = fg.get(0, DcmFGTypes::EFG_PIXELMEASURES, perFrame); - if (group && (pm = OFstatic_cast(FGPixelMeasures*, group))) - { - cond = pm->getSliceThickness(sliceThickness); - if (cond.good()) - { - DCMSEG_DEBUG("groupFramesByPosition(): Slice Thickness is " << sliceThickness); - // Compute mean distance of frame positions in x, y and z direction - Float64 means[3] = {0}; - Float64 diff[3] = {0}; - size_t count[3] = {0}; - for (size_t j = 0; j < m_FramePositions.size()-1; ++j) - { - for (size_t xyz = 0; xyz < 3; ++xyz) - { - diff[xyz] = fabs(m_FramePositions[j].m_position[xyz] - m_FramePositions[j+1].m_position[xyz]); - if (diff[xyz] > sliceThickness*0.5) - { - means[xyz] += diff[xyz]; - count[xyz]++; - } - } - } - // Compute mean value for each coordinate - for (size_t xyz = 0; xyz < 3; ++xyz) - { - if (count[xyz] > 0) means[xyz] = means[xyz] / count[xyz]; - } - - // Output variance for debug purposes - DCMSEG_DEBUG("groupFramesByPosition(): Mean distance of x coordinate is " << means[0]); - DCMSEG_DEBUG("groupFramesByPosition(): Mean distance of y coordinate is " << means[1]); - DCMSEG_DEBUG("groupFramesByPosition(): Mean distance of z coordinate is " << means[2]); - // Ensure that mean diestance of one coordinate is close to slice thickness (not smaller than 90& as a rule of thumb) - // or multiple of it - if (fabs(means[0]) > sliceThickness*0.9) relevantCoordinate = 0; - else if (fabs(means[1]) > sliceThickness*0.9) relevantCoordinate = 1; - else if (fabs(means[2]) > sliceThickness*0.9) relevantCoordinate = 2; - // else if variance of all three coordinates is clos to zero, all frames are at the same position - // and we can use any of the coordinates for sorting, so we arbitrarily choose x. - // The worst thing that could happen is that we map segments internally to the same frame - // position later, and therefore, create more distinct segments than necessary. - else if (fabs(means[0]) < 0.1*sliceThickness && fabs(means[1]) < 0.1*sliceThickness && fabs(means[2]) < 0.1*sliceThickness) - { - DCMSEG_DEBUG("Frames are at the same position, arbitrarily choosing first coordinate for sorting"); - relevantCoordinate = 0; - } - else relevantCoordinate = 100; // error - if (relevantCoordinate < 100) - { - DCMSEG_DEBUG("Using coordinate " << OFstatic_cast(Uint16, relevantCoordinate) << " for sorting frames by position"); - ComparePositions c(relevantCoordinate); - std::sort(m_FramePositions.begin(), m_FramePositions.end(), c); - // vec will contain all frame numbers that are at the same position - OFVector vec; - vec.push_back(m_FramePositions[0].m_frameNumber); - m_logicalFramePositions.push_back(vec); // Initialize for first logical frame - for (size_t j = 1; j < m_FramePositions.size(); ++j) - { - // If frame is close to previous frame, add it to the same vector. - // 2.5 is chosen since it means the frames are not further away if clearly less than half a slice thickness - Float64 diff = fabs(m_FramePositions[j].m_position[relevantCoordinate] - m_FramePositions[j-1].m_position[relevantCoordinate]); - DCMSEG_DEBUG("Coordinates of both frames:"); - DCMSEG_DEBUG("Frame " << j << ": " << m_FramePositions[j].m_position[0] << ", " << m_FramePositions[j].m_position[1] << ", " << m_FramePositions[j].m_position[2]); - DCMSEG_DEBUG("Frame " << j-1 << ": " << m_FramePositions[j-1].m_position[0] << ", " << m_FramePositions[j-1].m_position[1] << ", " << m_FramePositions[j-1].m_position[2]); - DCMSEG_DEBUG("groupFramesByPosition(): Frame " << j << " is " << diff << " mm away from previous frame"); - // 1% inaccuracy for slice thickness will be considered as same logical position - if (diff < sliceThickness*0.01) - { - // Add frame to last vector - DCMSEG_DEBUG("Assigning to same frame bucket as previous frame"); - m_logicalFramePositions.back().push_back(m_FramePositions[j].m_frameNumber); // physical frame number - } - else - { - DCMSEG_DEBUG("Assigning to same new frame bucket"); - // Create new vector - OFVector vec; - vec.push_back(m_FramePositions[j].m_frameNumber); - m_logicalFramePositions.push_back(vec); - } - } - } else - { - DCMSEG_ERROR("groupFramesByPosition(): Slice Thickness not represented in frame positions, cannot sort frames by position"); - cond = EC_TagNotFound; - } - - } - else - { - DCMSEG_ERROR("groupFramesByPosition(): Slice Thickness not found, cannot sort frames by position"); - cond = EC_TagNotFound; - } - } - else - { - DCMSEG_ERROR("groupFramesByPosition(): Pixel Measures FG not found, cannot sort frames by position"); - cond = EC_TagNotFound; - } - + cond = groupFramesByLogicalPosition(); } + // print frame groups if debug log level is enabled: if (cond.good() && DCM_dcmsegLogger.isEnabledFor(OFLogger::DEBUG_LOG_LEVEL)) { @@ -316,7 +200,8 @@ OFCondition OverlapUtil::groupFramesByPosition() OFStringStream ss; for (size_t j = 0; j < m_logicalFramePositions[i].size(); ++j) { - if (j > 0) ss << ", "; + if (j > 0) + ss << ", "; ss << m_logicalFramePositions[i][j]; } DCMSEG_DEBUG("groupFramesByPosition(): Logical frame #" << i << ": " << ss.str()); @@ -326,22 +211,33 @@ OFCondition OverlapUtil::groupFramesByPosition() if (cond.bad()) { - m_FramePositions.clear(); + m_framePositions.clear(); m_logicalFramePositions.clear(); } return cond; } - OFCondition OverlapUtil::getSegmentsByPosition(SegmentsByPosition& result) { + if (!m_segmentsByPosition.empty()) + { + // Already computed + result = m_segmentsByPosition; + return EC_Normal; + } + // Make sure prerequisites are met OFTimer tm; - OFCondition cond; + OFCondition cond = groupFramesByPosition(); + if (cond.bad()) + { + return cond; + } size_t numSegments = m_seg->getNumberOfSegments(); if (m_logicalFramePositions.empty()) { - cond = getFramesByPosition(m_logicalFramePositions); - if (cond.bad()) return cond; + cond = getFramesByPosition(m_logicalFramePositions); + if (cond.bad()) + return cond; } m_segmentsByPosition.clear(); m_segmentsByPosition.resize(m_logicalFramePositions.size()); @@ -352,34 +248,39 @@ OFCondition OverlapUtil::getSegmentsByPosition(SegmentsByPosition& result) { Uint32 frameNumber = m_logicalFramePositions[l][f]; OFVector segs; - FGBase* group = NULL; + FGBase* group = NULL; FGSegmentation* segFG = NULL; - group = m_seg->getFunctionalGroups().get(frameNumber, DcmFGTypes::EFG_SEGMENTATION); - segFG = OFstatic_cast(FGSegmentation*, group); + group = m_seg->getFunctionalGroups().get(frameNumber, DcmFGTypes::EFG_SEGMENTATION); + segFG = OFstatic_cast(FGSegmentation*, group); if (segFG) { Uint16 segNum = 0; - cond = segFG->getReferencedSegmentNumber(segNum); + cond = segFG->getReferencedSegmentNumber(segNum); if (cond.good() && segNum > 0 && (segNum <= numSegments)) { m_segmentsByPosition[l].insert(SegNumAndFrameNum(segNum, frameNumber)); } else if (segNum == 0) { - DCMSEG_ERROR("getSegmentsByPosition(): Referenced Segment Number is 0 (not permitted), cannot add segment"); + DCMSEG_ERROR( + "getSegmentsByPosition(): Referenced Segment Number is 0 (not permitted), cannot add segment"); cond = EC_InvalidValue; break; } else if (segNum > numSegments) { - DCMSEG_ERROR("getSegmentsByPosition(): Found Referenced Segment Number " << segNum << " but only " << numSegments << " segments are present, cannot add segment"); - DCMSEG_ERROR("getSegmentsByPosition(): Segments are not numbered consecutively, cannot add segment"); + DCMSEG_ERROR("getSegmentsByPosition(): Found Referenced Segment Number " + << segNum << " but only " << numSegments + << " segments are present, cannot add segment"); + DCMSEG_ERROR( + "getSegmentsByPosition(): Segments are not numbered consecutively, cannot add segment"); cond = EC_InvalidValue; break; } else { - DCMSEG_ERROR("getSegmentsByPosition(): Referenced Segment Number not found (not permitted) , cannot add segment"); + DCMSEG_ERROR("getSegmentsByPosition(): Referenced Segment Number not found (not permitted) , " + "cannot add segment"); cond = EC_TagNotFound; break; } @@ -393,32 +294,27 @@ OFCondition OverlapUtil::getSegmentsByPosition(SegmentsByPosition& result) // print segments per logical frame if debug log level is enabled if (cond.good() && DCM_dcmsegLogger.isEnabledFor(OFLogger::DEBUG_LOG_LEVEL)) { - DCMSEG_DEBUG("getSegmentsByPosition(): Segments grouped by logical frame positions, (seg#,frame#):"); - for (size_t i = 0; i < m_segmentsByPosition.size(); ++i) - { - OFStringStream ss; - for (std::set::iterator it = m_segmentsByPosition[i].begin(); it != m_segmentsByPosition[i].end(); ++it) - { - if (it != m_segmentsByPosition[i].begin()) ss << ","; - ss << "(" << (*it).m_segmentNumber << "," << (*it).m_frameNumber << ")"; - } - DCMSEG_DEBUG("getSegmentsByPosition(): Logical frame #" << i << ": " << ss.str()); - } + OFStringStream ss; + printSegmentsByPosition(ss); + DCMSEG_DEBUG(ss.str()); } DCMSEG_DEBUG("groupFramesByPosition(): Grouping segments by position took " << tm.getDiff() << " s"); return cond; } - OFCondition OverlapUtil::getOverlapMatrix(OverlapMatrix& matrix) { - OFTimer tm; - OFCondition result; - if (m_segmentsByPosition.empty()) + if (!m_segmentOverlapMatrix.empty()) { - result = getSegmentsByPosition(m_segmentsByPosition); + // Already computed + matrix = m_segmentOverlapMatrix; + return EC_Normal; } + // Make sure prerequisites are met + OFTimer tm; + SegmentsByPosition dontCare; + OFCondition result = getSegmentsByPosition(dontCare); if (result.good()) { result = buildOverlapMatrix(); @@ -431,7 +327,6 @@ OFCondition OverlapUtil::getOverlapMatrix(OverlapMatrix& matrix) return result; } - OFCondition OverlapUtil::getNonOverlappingSegments(SegmentGroups& segmentGroups) { OFTimer tm; @@ -439,12 +334,11 @@ OFCondition OverlapUtil::getNonOverlappingSegments(SegmentGroups& segmentGroups) if (!m_nonOverlappingSegments.empty()) { // Already computed + segmentGroups = m_nonOverlappingSegments; return EC_Normal; } - if (m_segmentOverlapMatrix.empty()) - { - result = getOverlapMatrix(m_segmentOverlapMatrix); - } + // Make sure prerequisites are met + result = getOverlapMatrix(m_segmentOverlapMatrix); if (result.good()) { // Group those segments from the overlap matrix together, that do not @@ -452,7 +346,7 @@ OFCondition OverlapUtil::getNonOverlappingSegments(SegmentGroups& segmentGroups) // Go through all segments and check if they overlap with any of the already // grouped segments. If not, add them to the same group. If yes, create a new group // and add them there. - m_nonOverlappingSegments.push_back(std::set()); + m_nonOverlappingSegments.push_back(OFVector()); for (size_t i = 0; i < m_segmentOverlapMatrix.size(); ++i) { // Loop over all groups and check whether the current segment overlaps with any of them @@ -460,7 +354,9 @@ OFCondition OverlapUtil::getNonOverlappingSegments(SegmentGroups& segmentGroups) for (size_t j = 0; j < m_nonOverlappingSegments.size(); ++j) { // Loop over all segments in the current group - for (std::set::iterator it = m_nonOverlappingSegments[j].begin(); it != m_nonOverlappingSegments[j].end(); ++it) + for (OFVector::iterator it = m_nonOverlappingSegments[j].begin(); + it != m_nonOverlappingSegments[j].end(); + ++it) { // Check if the current segment overlaps with the segment in the current group if (m_segmentOverlapMatrix[i][(*it) - 1] == 1) @@ -472,15 +368,15 @@ OFCondition OverlapUtil::getNonOverlappingSegments(SegmentGroups& segmentGroups) if (!overlaps) { // Add segment to current group - m_nonOverlappingSegments[j].insert(i + 1); + m_nonOverlappingSegments[j].push_back(i + 1); break; } } if (overlaps) { // Create new group and add segment to it - m_nonOverlappingSegments.push_back(std::set()); - m_nonOverlappingSegments.back().insert(i + 1); + m_nonOverlappingSegments.push_back(OFVector()); + m_nonOverlappingSegments.back().push_back(i + 1); } } } @@ -502,6 +398,23 @@ OFCondition OverlapUtil::getNonOverlappingSegments(SegmentGroups& segmentGroups) return result; } +void OverlapUtil::printSegmentsByPosition(OFStringStream& ss) +{ + ss << "printSegmentsByPosition(): Segments grouped by logical frame positions, (seg#,frame#):" << OFendl; + for (size_t i = 0; i < m_segmentsByPosition.size(); ++i) + { + OFStringStream tempSS; + for (std::set::iterator it = m_segmentsByPosition[i].begin(); + it != m_segmentsByPosition[i].end(); + ++it) + { + if (it != m_segmentsByPosition[i].begin()) + tempSS << ","; + tempSS << "(" << (*it).m_segmentNumber << "," << (*it).m_frameNumber << ")"; + } + ss << "getSegmentsByPosition(): Logical frame #" << i << ": " << tempSS.str(); + } +} void OverlapUtil::printOverlapMatrix(OFStringStream& ss) { @@ -520,23 +433,24 @@ void OverlapUtil::printOverlapMatrix(OFStringStream& ss) } } - void OverlapUtil::printNonOverlappingSegments(OFStringStream& ss) { ss << "printNonOverlappingSegments(): Non-overlapping segments:" << OFendl; for (size_t i = 0; i < m_nonOverlappingSegments.size(); ++i) { ss << "Group #" << i << ": "; - for (std::set::iterator it = m_nonOverlappingSegments[i].begin(); it != m_nonOverlappingSegments[i].end(); ++it) + for (OFVector::iterator it = m_nonOverlappingSegments[i].begin(); + it != m_nonOverlappingSegments[i].end(); + ++it) { - if (it != m_nonOverlappingSegments[i].begin()) ss << ", "; + if (it != m_nonOverlappingSegments[i].begin()) + ss << ", "; ss << (*it); } ss << OFendl; } } - OFCondition OverlapUtil::buildOverlapMatrix() { // Make 2 dimensional array matrix of Sint8 type for (segment numbers) X (segment numbers). @@ -552,7 +466,6 @@ OFCondition OverlapUtil::buildOverlapMatrix() { m_segmentOverlapMatrix[i][i] = 0; } - // Go through all logical frame positions, and compare all segments at each position size_t index1, index2; index1 = index2 = 0; @@ -560,22 +473,30 @@ OFCondition OverlapUtil::buildOverlapMatrix() { DCMSEG_DEBUG("getOverlappingSegments(): Comparing segments at logical frame position " << i); // Compare all segments at this position - for (std::set::iterator it = m_segmentsByPosition[i].begin(); it != m_segmentsByPosition[i].end(); ++it) + for (std::set::iterator it = m_segmentsByPosition[i].begin(); + it != m_segmentsByPosition[i].end(); + ++it) { index1++; - for (std::set::iterator it2 = m_segmentsByPosition[i].begin(); it2 != m_segmentsByPosition[i].end(); ++it2) + for (std::set::iterator it2 = m_segmentsByPosition[i].begin(); + it2 != m_segmentsByPosition[i].end(); + ++it2) { index2++; // Skip comparison of same segments in reverse order (index2 < index1) - if (index2 <= index1) continue; + if (index2 <= index1) + continue; // Skip self-comparison (diagonal is always 0); (index1==index2) if (it->m_segmentNumber != it2->m_segmentNumber) { // Check if we already have found an overlap on another logical frame, and if so, skip - Sint8 existing_result = m_segmentOverlapMatrix[(*it).m_segmentNumber - 1][(*it2).m_segmentNumber - 1]; + Sint8 existing_result + = m_segmentOverlapMatrix[(*it).m_segmentNumber - 1][(*it2).m_segmentNumber - 1]; if (existing_result == 1) { - DCMSEG_DEBUG("getOverlappingSegments(): Skipping frame comparison on pos #" << i << " for segments " << (*it).m_segmentNumber << " and " << (*it2).m_segmentNumber << " (already marked as overlapping)"); + DCMSEG_DEBUG("getOverlappingSegments(): Skipping frame comparison on pos #" + << i << " for segments " << (*it).m_segmentNumber << " and " + << (*it2).m_segmentNumber << " (already marked as overlapping)"); continue; } // Compare pixels of the frames referenced by each segments. @@ -602,7 +523,6 @@ OFCondition OverlapUtil::buildOverlapMatrix() } } } - // print overlap matrix if debug log level is enabled if (DCM_dcmsegLogger.isEnabledFor(OFLogger::DEBUG_LOG_LEVEL)) { @@ -613,8 +533,6 @@ OFCondition OverlapUtil::buildOverlapMatrix() return EC_Normal; } - - OFCondition OverlapUtil::checkFramesOverlap(const Uint32& f1, const Uint32& f2, OFBool& overlap) { if (f1 == f2) @@ -628,12 +546,11 @@ OFCondition OverlapUtil::checkFramesOverlap(const Uint32& f1, const Uint32& f2, const DcmIODTypes::Frame* f1_data = m_seg->getFrame(f1); const DcmIODTypes::Frame* f2_data = m_seg->getFrame(f2); Uint16 rows, cols; - rows=cols=0; - DcmIODImage > *ip = - static_cast > *>(m_seg); + rows = cols = 0; + DcmIODImage>* ip = static_cast>*>(m_seg); ip->getImagePixel().getRows(rows); ip->getImagePixel().getColumns(cols); - if ( rows*cols % 8 != 0 ) + if (rows * cols % 8 != 0) { // We must copmare pixel by pixel of the unpacked frames (for now) result = checkFramesOverlapUnpacked(f1, f2, f1_data, f2_data, rows, cols, overlap); @@ -650,8 +567,13 @@ OFCondition OverlapUtil::checkFramesOverlap(const Uint32& f1, const Uint32& f2, return result; } - -OFCondition OverlapUtil::checkFramesOverlapBinary(const Uint32& f1, const Uint32& f2, const DcmIODTypes::Frame* f1_data, const DcmIODTypes::Frame* f2_data, const Uint16& rows, const Uint16 cols, OFBool& overlap) +OFCondition OverlapUtil::checkFramesOverlapBinary(const Uint32& f1, + const Uint32& f2, + const DcmIODTypes::Frame* f1_data, + const DcmIODTypes::Frame* f2_data, + const Uint16& rows, + const Uint16 cols, + OFBool& overlap) { DCMSEG_DEBUG("checkFramesOverlap(): Comparing frames " << f1 << " and " << f2 << " for overlap (fast binary mode)"); if (!f1_data || !f2_data) @@ -661,23 +583,18 @@ OFCondition OverlapUtil::checkFramesOverlapBinary(const Uint32& f1, const Uint32 } if (f1_data->length != f2_data->length) { - DCMSEG_ERROR("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " have different length, cannot compare"); + DCMSEG_ERROR("checkFramesOverlap(): Frames " << f1 << " and " << f2 + << " have different length, cannot compare"); return EC_IllegalCall; } - + // Compare byte (8 pixels at once) using bitwise AND (if both have a 1 at the same position, they overlap) for (size_t n = 0; n < f1_data->length; ++n) { - // Compate byte (8 pixels at once) using bitwise AND (if both have a 1 at the same position, they overlap) - if (f1 == 9732) - { - if ( (OFstatic_cast(Uint16, f1_data->pixData[n]) > 0) || (OFstatic_cast(Uint16, f2_data->pixData[n]) > 0) ) - { - DCMSEG_DEBUG("checkFramesOverlap(): Comparing pixel value " << OFstatic_cast(Uint16, f1_data->pixData[n]) << " of frame " << f1 << " with pixel value " << OFstatic_cast(Uint16, f2_data->pixData[n]) << " of frame " << f2); - } - } if (f1_data->pixData[n] & f2_data->pixData[n]) { - DCMSEG_DEBUG("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " do overlap, pixel value " << OFstatic_cast(Uint16, f1_data->pixData[n]) << " at index " << n << " is the same"); + DCMSEG_DEBUG("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " do overlap, pixel value " + << OFstatic_cast(Uint16, f1_data->pixData[n]) << " at index " + << n << " is the same"); overlap = OFTrue; break; } @@ -685,10 +602,16 @@ OFCondition OverlapUtil::checkFramesOverlapBinary(const Uint32& f1, const Uint32 return EC_Normal; } - -OFCondition OverlapUtil::checkFramesOverlapUnpacked(const Uint32& f1, const Uint32& f2, const DcmIODTypes::Frame* f1_data, const DcmIODTypes::Frame* f2_data, const Uint16& rows, const Uint16 cols, OFBool& overlap) +OFCondition OverlapUtil::checkFramesOverlapUnpacked(const Uint32& f1, + const Uint32& f2, + const DcmIODTypes::Frame* f1_data, + const DcmIODTypes::Frame* f2_data, + const Uint16& rows, + const Uint16 cols, + OFBool& overlap) { - DCMSEG_DEBUG("checkFramesOverlap(): Comparing frames " << f1 << " and " << f2 << " for overlap (slow unpacked mode)"); + DCMSEG_DEBUG("checkFramesOverlap(): Comparing frames " << f1 << " and " << f2 + << " for overlap (slow unpacked mode)"); OFunique_ptr f1_unpacked(DcmSegUtils::unpackBinaryFrame(f1_data, rows, cols)); OFunique_ptr f2_unpacked(DcmSegUtils::unpackBinaryFrame(f2_data, rows, cols)); if (!f1_unpacked || !f2_unpacked) @@ -698,17 +621,19 @@ OFCondition OverlapUtil::checkFramesOverlapUnpacked(const Uint32& f1, const Uint } if (f1_unpacked->length != f2_unpacked->length) { - DCMSEG_ERROR("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " have different length, cannot compare"); + DCMSEG_ERROR("checkFramesOverlap(): Frames " << f1 << " and " << f2 + << " have different length, cannot compare"); return EC_IllegalCall; } - // Compare pixels of both frames and check whether at least one has the same value DCMSEG_DEBUG("checkFramesOverlap(): Comparing frames " << f1 << " and " << f2 << " for overlap"); for (size_t n = 0; n < f1_unpacked->length; ++n) { if (f1_unpacked->pixData[n] != 0 && (f1_unpacked->pixData[n] == f2_unpacked->pixData[n])) { - DCMSEG_DEBUG("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " do overlap, pixel value " << OFstatic_cast(Uint16, f1_unpacked->pixData[n]) << " at index " << n << " is the same"); + DCMSEG_DEBUG("checkFramesOverlap(): Frames " << f1 << " and " << f2 << " do overlap, pixel value " + << OFstatic_cast(Uint16, f1_unpacked->pixData[n]) + << " at index " << n << " is the same"); overlap = OFTrue; break; } @@ -716,4 +641,185 @@ OFCondition OverlapUtil::checkFramesOverlapUnpacked(const Uint32& f1, const Uint return EC_Normal; } +OFCondition OverlapUtil::collectPhysicalFramePositions() +{ + // Group all frames by position into m_logicalFramePositions. + FGInterface& fg = m_seg->getFunctionalGroups(); + size_t numFrames = m_seg->getNumberOfFrames(); + OFBool perFrame = OFFalse; + OFCondition cond; + // Vector of frame numbers with their respective position + m_framePositions.clear(); + m_framePositions.reserve(numFrames); + + // Put all frames into vector along with their Image Position Patient coordinates + for (size_t i = 0; i < numFrames; ++i) + { + FGPlanePosPatient* ppp = NULL; + FGBase* group = fg.get(i, DcmFGTypes::EFG_PLANEPOSPATIENT, perFrame); + if (group) + ppp = OFstatic_cast(FGPlanePosPatient*, group); + if (ppp) + { + // Get image position patient for frame i + OFVector ipp(3); + // Only in later DCMTK version: + // cond = ppp->getImagePositionPatient(ipp); + cond = ppp->getImagePositionPatient(ipp[0], ipp[1], ipp[2]); + if (cond.good()) + { + // Insert frame into map + m_framePositions.push_back(FramePositionAndNumber(ipp, i)); + // m_frameNumToPositionIndex.push_back(m_FramePositions.size() - 1); // current position in + // m_logicalFramePositions + } + else + { + DCMSEG_ERROR("groupFramesByPosition(): Image Position Patient not found for frame " + << i << ", cannot sort frames by position"); + cond = EC_TagNotFound; + break; + } + } + else + { + DCMSEG_ERROR("groupFramesByPosition(): Image Position Patient not found for frame " + << i << ", cannot sort frames by position"); + cond = EC_TagNotFound; + break; + } + } + return cond; +} + +OFCondition OverlapUtil::groupFramesByLogicalPosition() +{ + OFCondition cond; + FGInterface& fg = m_seg->getFunctionalGroups(); + OFBool perFrame = OFFalse; + + // Find all distinct positions and for each position the actual frames that can be found at it + Float64 sliceThickness = 0.0; + FGPixelMeasures* pm = NULL; + Uint8 relevantCoordinate = 0; + FGBase* group = fg.get(0, DcmFGTypes::EFG_PIXELMEASURES, perFrame); + // Get/compute Slice Thickness + if (group && (pm = OFstatic_cast(FGPixelMeasures*, group))) + { + cond = pm->getSliceThickness(sliceThickness); + if (cond.good()) + { + DCMSEG_DEBUG("groupFramesByPosition(): Slice Thickness is " << sliceThickness); + // Compute mean distance of frame positions in x, y and z direction + Float64 means[3] = { 0 }; + Float64 diff[3] = { 0 }; + size_t count[3] = { 0 }; + for (size_t j = 0; j < m_framePositions.size() - 1; ++j) + { + for (size_t xyz = 0; xyz < 3; ++xyz) + { + diff[xyz] = fabs(m_framePositions[j].m_position[xyz] - m_framePositions[j + 1].m_position[xyz]); + if (diff[xyz] > sliceThickness * 0.5) + { + means[xyz] += diff[xyz]; + count[xyz]++; + } + } + } + // Compute mean value for each coordinate + for (size_t xyz = 0; xyz < 3; ++xyz) + { + if (count[xyz] > 0) + means[xyz] = means[xyz] / count[xyz]; + } + + // Output mean value for debug purposes + DCMSEG_DEBUG("groupFramesByPosition(): Mean distance of x coordinate is " << means[0]); + DCMSEG_DEBUG("groupFramesByPosition(): Mean distance of y coordinate is " << means[1]); + DCMSEG_DEBUG("groupFramesByPosition(): Mean distance of z coordinate is " << means[2]); + // Ensure that mean distance of one coordinate is close to slice thickness (not smaller than 90& as a rule + // of thumb) or multiple of it + if (fabs(means[0]) > sliceThickness * 0.9) + relevantCoordinate = 0; + else if (fabs(means[1]) > sliceThickness * 0.9) + relevantCoordinate = 1; + else if (fabs(means[2]) > sliceThickness * 0.9) + relevantCoordinate = 2; + // else if deviation of all three coordinates is close to zero, all frames are at the same position + // and we can use any of the coordinates for sorting, so we arbitrarily choose x. + // The worst thing that could happen is that we map segments internally to the same frame + // position later, and therefore, create more distinct segments than necessary. + else if (fabs(means[0]) < 0.1 * sliceThickness && fabs(means[1]) < 0.1 * sliceThickness + && fabs(means[2]) < 0.1 * sliceThickness) + { + DCMSEG_DEBUG("Frames are at the same position, arbitrarily choosing first coordinate for sorting"); + relevantCoordinate = 0; + } + else + relevantCoordinate = 100; // error + if (relevantCoordinate < 100) + { + DCMSEG_DEBUG("Using coordinate " << OFstatic_cast(Uint16, relevantCoordinate) + << " for sorting frames by position"); + ComparePositions c(relevantCoordinate); + std::sort(m_framePositions.begin(), m_framePositions.end(), c); + // vec will contain all frame numbers that are at the same position + OFVector vec; + vec.push_back(m_framePositions[0].m_frameNumber); + m_logicalFramePositions.push_back(vec); // Initialize for first logical frame + for (size_t j = 1; j < m_framePositions.size(); ++j) + { + // If frame is close to previous frame, add it to the same vector. + // 2.5 is chosen since it means the frames are not further away if clearly less than half a slice + // thickness + Float64 diff = fabs(m_framePositions[j].m_position[relevantCoordinate] + - m_framePositions[j - 1].m_position[relevantCoordinate]); + DCMSEG_DEBUG("Coordinates of both frames:"); + DCMSEG_DEBUG("Frame " << j << ": " << m_framePositions[j].m_position[0] << ", " + << m_framePositions[j].m_position[1] << ", " + << m_framePositions[j].m_position[2]); + DCMSEG_DEBUG("Frame " << j - 1 << ": " << m_framePositions[j - 1].m_position[0] << ", " + << m_framePositions[j - 1].m_position[1] << ", " + << m_framePositions[j - 1].m_position[2]); + DCMSEG_DEBUG("groupFramesByPosition(): Frame " << j << " is " << diff + << " mm away from previous frame"); + // 1% inaccuracy for slice thickness will be considered as same logical position + if (diff < sliceThickness * 0.01) + { + // Add frame to last vector + DCMSEG_DEBUG("Assigning to same frame bucket as previous frame"); + m_logicalFramePositions.back().push_back( + m_framePositions[j].m_frameNumber); // physical frame number + } + else + { + DCMSEG_DEBUG("Assigning to same new frame bucket"); + // Create new vector + OFVector vec; + vec.push_back(m_framePositions[j].m_frameNumber); + m_logicalFramePositions.push_back(vec); + } + } + } + else + { + DCMSEG_ERROR("groupFramesByPosition(): Slice Thickness not represented in frame positions, cannot sort " + "frames by position"); + cond = EC_TagNotFound; + } + } + else + { + DCMSEG_ERROR("groupFramesByPosition(): Slice Thickness not found, cannot sort frames by position"); + cond = EC_TagNotFound; + } + } + else + { + DCMSEG_ERROR("groupFramesByPosition(): Pixel Measures FG not found, cannot sort frames by position"); + cond = EC_TagNotFound; + } + return cond; +} + }