From 5a9cc2f680706dfbe3acd9be90ccfd765ea78767 Mon Sep 17 00:00:00 2001 From: Michael Onken Date: Tue, 5 Dec 2023 10:12:31 +0100 Subject: [PATCH 1/3] Ensure to also handle more than 65535 frames. In order to make the class more generic, also handle more than 65535 frames (Number of Frames can be 2^32-1 max). --- include/dcmqi/framesorter.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/dcmqi/framesorter.h b/include/dcmqi/framesorter.h index ae2eaf87..919454cf 100644 --- a/include/dcmqi/framesorter.h +++ b/include/dcmqi/framesorter.h @@ -69,7 +69,7 @@ class FrameSorter /// The error code should be set in any case (default: EC_Normal) OFCondition errorCode; /// The frame numbers, in sorted order (default: empty) - OFVector frameNumbers; + OFVector frameNumbers; /// Tag key that contains the information that was crucial for sorting. /// This is especially useful for creating dimension indices. Should be /// set to (0xffff,0xfff) if none was used (default). From 05b36fbf716a7097ad425fc2705c007c0f12b2f8 Mon Sep 17 00:00:00 2001 From: Michael Onken Date: Tue, 5 Dec 2023 11:50:18 +0100 Subject: [PATCH 2/3] Remember frame positions in Result, small changes. Also remember frame positions in sorting result since it is usually accessed anyway and interesting for library users that otherwise call the underlying method on FGInterface for all frames again. The DummySorter does not provide frame positions to keep it as simple as possible (i.e. keep it "zero" cost). Allow 2^32-1 frames as permitted in in DICOM (i.e. use Uint32 instead of Uint16) for frame number index. Changed precision from Float32 to Float64 where applicable. Added some documentation. --- include/dcmqi/framesorter.h | 58 +++++++++++++++++++++++++++++-------- 1 file changed, 46 insertions(+), 12 deletions(-) diff --git a/include/dcmqi/framesorter.h b/include/dcmqi/framesorter.h index 919454cf..8fb8aab0 100644 --- a/include/dcmqi/framesorter.h +++ b/include/dcmqi/framesorter.h @@ -35,6 +35,9 @@ #include #include +typedef OFVector ImagePosition; // Image Position Patient +typedef OFPair FrameWithPos; // DICOM frame number and Image Position Patient + /** Abstract class for sorting a set of frames in a functional group. The * sorting criteria are up to the actual implementation classes. */ @@ -70,6 +73,10 @@ class FrameSorter OFCondition errorCode; /// The frame numbers, in sorted order (default: empty) OFVector frameNumbers; + /// The frame positions, in sorted order, if provided + /// by the sorter (default: empty). If not empty, contains the same number + /// of items as frameNumbers. + OFVector framePositions; /// Tag key that contains the information that was crucial for sorting. /// This is especially useful for creating dimension indices. Should be /// set to (0xffff,0xfff) if none was used (default). @@ -87,6 +94,11 @@ class FrameSorter OFString fgPrivateCreator; }; + struct FramePositions + { + OFVector positions; + }; + /** Default constructor, does nothing */ FrameSorter(){}; @@ -126,6 +138,12 @@ class FrameSorter FGInterface* m_fg; }; +/** Dummy sorter implementing the FrameSorter interface, + * but not doing any sorting at all. As a result it provides + * a list of frames in their natural order, as found in the underlying + * DICOM dataset. The results will not contain frame position + * to make this implementation as lightweight and "stupid" as possible. + */ class FrameSorterIdentity : public FrameSorter { @@ -143,7 +161,11 @@ class FrameSorterIdentity : public FrameSorter return "Returns frames in the order defined in the functional group, i.e. as defined in the image file"; } - + /** Performs actual sorting. Does only set Results.frameNumbers and errorCode, + * leaving the rest untouched. + * @param results The results produced by dummy sorter (list of frame numbers as + * found in the underlying DICOM dataset, and EC_Normal as error code) + */ virtual void sort(Results& results) { if (m_fg == NULL) @@ -168,6 +190,11 @@ class FrameSorterIdentity : public FrameSorter }; +/** Sorter implementing the FrameSorter interface that sorts frames + * by their Image Position (Patient) attribute. The sorting is done + * by projecting the Image Position (Patient) on the slice direction + * (as defined by the Image Orientation (Patient) attribute). + */ class FrameSorterIPP : public FrameSorter { public: @@ -179,10 +206,15 @@ class FrameSorterIPP : public FrameSorter frameId() {} - Float32 key; - unsigned frameId; + /// The key relevant for sorting + Float64 key; + /// The DICOM frame number + Uint32 frameId; + /// The frame position as found in Image Position (Patient) + ImagePosition framePos; }; + FrameSorterIPP(){}; ~FrameSorterIPP(){}; @@ -201,7 +233,7 @@ class FrameSorterIPP : public FrameSorter return; } - OFVector dirX, dirY; + OFVector dirX, dirY; OFString orientStr; for(int i=0;i<3;i++){ if(planorfg->getImageOrientationPatient(orientStr, i).good()){ @@ -253,7 +285,7 @@ class FrameSorterIPP : public FrameSorter return; } - OFVector sOrigin; + ImagePosition sOrigin; for(int j=0;j<3;j++){ OFString planposStr; if(planposfg->getImagePositionPatient(planposStr, j).good()){ @@ -264,10 +296,11 @@ class FrameSorterIPP : public FrameSorter } } - Float32 dist; + Float64 dist; dist = dot(sliceDirection, sOrigin); orderedFrameItems[frameId].key = dist; orderedFrameItems[frameId].frameId = frameId; + orderedFrameItems[frameId].framePos = sOrigin; } qsort(&orderedFrameItems[0], numFrames, sizeof(OrderedFrameItem), &compareIPPKeys); @@ -275,6 +308,7 @@ class FrameSorterIPP : public FrameSorter for(Uint32 count=0;count cross_3d(OFVector v1, OFVector v2){ - OFVector result; + OFVector cross_3d(OFVector v1, OFVector v2){ + OFVector result; result.push_back(v1[1]*v2[2]-v1[2]*v2[1]); result.push_back(v1[2]*v2[0]-v1[0]*v2[2]); result.push_back(v1[0]*v2[1]-v1[1]*v2[0]); @@ -293,11 +327,11 @@ class FrameSorterIPP : public FrameSorter return result; } - Float32 dot(OFVector v1, OFVector v2){ - return Float32(v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]); + Float64 dot(OFVector v1, OFVector v2){ + return Float64(v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]); } - void normalize(OFVector &v){ + void normalize(OFVector &v){ double norm = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); v[0] = v[0]/norm; v[1] = v[1]/norm; @@ -306,7 +340,7 @@ class FrameSorterIPP : public FrameSorter static int compareIPPKeys(const void *a, const void *b); - OFVector sliceDirection; + OFVector sliceDirection; }; From 3614ed352e17ac2b0aa0eba08f9fc4f9fe1bd1f3 Mon Sep 17 00:00:00 2001 From: Michael Onken Date: Tue, 5 Dec 2023 11:51:42 +0100 Subject: [PATCH 3/3] Use framesorter class instead homegrown sorting. --- include/dcmqi/OverlapUtil.h | 5 +- libsrc/OverlapUtil.cpp | 179 ++++++++++++------------------------ 2 files changed, 61 insertions(+), 123 deletions(-) diff --git a/include/dcmqi/OverlapUtil.h b/include/dcmqi/OverlapUtil.h index 3116d222..dcc37a59 100644 --- a/include/dcmqi/OverlapUtil.h +++ b/include/dcmqi/OverlapUtil.h @@ -233,10 +233,7 @@ class OverlapUtil void printNonOverlappingSegments(OFStringStream& ss); protected: - /** Collect all physical frame positions in m_FramePositions - * @return EC_Normal if successful, error otherwise - */ - OFCondition collectPhysicalFramePositions(); + /** Group physical frame positions into logical positions. This is done by sorting * frames after *that* position coordinate that in its mean position difference is diff --git a/libsrc/OverlapUtil.cpp b/libsrc/OverlapUtil.cpp index a5b77f6c..69fd37ee 100644 --- a/libsrc/OverlapUtil.cpp +++ b/libsrc/OverlapUtil.cpp @@ -20,6 +20,7 @@ */ #include "dcmqi/OverlapUtil.h" +#include "dcmqi/framesorter.h" #include "dcmtk/dcmdata/dcerror.h" #include "dcmtk/dcmfg/fginterface.h" #include "dcmtk/dcmfg/fgpixmsr.h" @@ -210,11 +211,23 @@ OFCondition OverlapUtil::groupFramesByPosition() // After that, all frames at the same position will be in the same vector // assigned to the same key (the frame's coordinates) in the map. // Group all frames by position into m_logicalFramePositions. - cond = collectPhysicalFramePositions(); - if (cond.good()) + FrameSorterIPP sorter; + sorter.setSorterInput(&(m_seg->getFunctionalGroups())); + FrameSorterIPP::Results results; + sorter.sort(results); + if (results.errorCode != EC_Normal) + { + DCMSEG_ERROR("groupFramesByPosition(): Cannot sort frames by position: " << results.errorCode.text()); + return results.errorCode; + } + // Copy results from frame sorter to overlap util framePositions member + m_framePositions.clear(); + m_framePositions.reserve(results.framePositions.size()); + for (size_t i = 0; i < results.framePositions.size(); ++i) { - cond = groupFramesByLogicalPosition(); + m_framePositions.push_back(FramePositionAndNumber(results.framePositions[i], results.frameNumbers[i])); } + cond = groupFramesByLogicalPosition(); // print frame groups if debug log level is enabled: if (cond.good() && DCM_dcmsegLogger.isEnabledFor(OFLogger::DEBUG_LOG_LEVEL)) @@ -665,140 +678,67 @@ OFCondition OverlapUtil::checkFramesOverlapUnpacked(const Uint32& f1, return EC_Normal; } -OFCondition OverlapUtil::collectPhysicalFramePositions() -{ - // Group all frames by position into m_logicalFramePositions. - FGInterface& fg = m_seg->getFunctionalGroups(); - size_t numFrames = m_seg->getNumberOfFrames(); - OFBool perFrame = OFFalse; - OFCondition cond; - // Vector of frame numbers with their respective position - m_framePositions.clear(); - m_framePositions.reserve(numFrames); - - // Put all frames into vector along with their Image Position Patient coordinates - for (size_t i = 0; i < numFrames; ++i) - { - FGPlanePosPatient* ppp = NULL; - FGBase* group = fg.get(i, DcmFGTypes::EFG_PLANEPOSPATIENT, perFrame); - if (group) - ppp = OFstatic_cast(FGPlanePosPatient*, group); - if (ppp) - { - // Get image position patient for frame i - OFVector ipp(3); - // Only in later DCMTK version: - // cond = ppp->getImagePositionPatient(ipp); - cond = ppp->getImagePositionPatient(ipp[0], ipp[1], ipp[2]); - if (cond.good()) - { - // Insert frame into map - m_framePositions.push_back(FramePositionAndNumber(ipp, i)); - } - else - { - DCMSEG_ERROR("groupFramesByPosition(): Image Position Patient not found for frame " - << i << ", cannot sort frames by position"); - cond = EC_TagNotFound; - break; - } - } - else - { - DCMSEG_ERROR("groupFramesByPosition(): Image Position Patient not found for frame " - << i << ", cannot sort frames by position"); - cond = EC_TagNotFound; - break; - } - } - return cond; -} OFCondition OverlapUtil::groupFramesByLogicalPosition() { OFCondition cond; FGInterface& fg = m_seg->getFunctionalGroups(); OFBool perFrame = OFFalse; - - // Find all distinct positions and for each position the actual frames that can be found at it Float64 sliceThickness = 0.0; FGPixelMeasures* pm = NULL; - Uint8 relevantCoordinate = 3; // not determined yet FGBase* group = fg.get(0, DcmFGTypes::EFG_PIXELMEASURES, perFrame); if (group && (pm = OFstatic_cast(FGPixelMeasures*, group))) { // Get/compute Slice Thickness cond = pm->getSliceThickness(sliceThickness); - if (cond.good()) + if (cond.bad()) { - DCMSEG_DEBUG("groupFramesByPosition(): Slice Thickness is " << sliceThickness); - // Identify coordinate to be used for frame sorting - relevantCoordinate = identifyChangingCoordinate(m_imageOrientation); - if (relevantCoordinate < 3) - { - DCMSEG_DEBUG("Using coordinate " << OFstatic_cast(Uint16, relevantCoordinate) - << " for sorting frames by position"); - ComparePositions c(relevantCoordinate); - std::sort(m_framePositions.begin(), m_framePositions.end(), c); - // vec will contain all frame numbers that are at the same position - OFVector vec; - vec.push_back(m_framePositions[0].m_frameNumber); - m_logicalFramePositions.push_back(vec); // Initialize for first logical frame - for (size_t j = 1; j < m_framePositions.size(); ++j) - { - // If frame is close to previous frame, add it to the same vector. - // 2.5 is chosen since it means the frames are not further away if clearly less than half a slice - // thickness - Float64 diff = fabs(m_framePositions[j].m_position[relevantCoordinate] - - m_framePositions[j - 1].m_position[relevantCoordinate]); - DCMSEG_DEBUG("Coordinates of both frames:"); - DCMSEG_DEBUG("Frame " << j << ": " << m_framePositions[j].m_position[0] << ", " - << m_framePositions[j].m_position[1] << ", " - << m_framePositions[j].m_position[2]); - DCMSEG_DEBUG("Frame " << j - 1 << ": " << m_framePositions[j - 1].m_position[0] << ", " - << m_framePositions[j - 1].m_position[1] << ", " - << m_framePositions[j - 1].m_position[2]); - DCMSEG_DEBUG("groupFramesByPosition(): Frame " << j << " is " << diff - << " mm away from previous frame"); - // 1% inaccuracy for slice thickness will be considered as same logical position - if (diff < sliceThickness * 0.01) - { - // Add frame to last vector - DCMSEG_DEBUG("Assigning to same frame bucket as previous frame"); - m_logicalFramePositions.back().push_back( - m_framePositions[j].m_frameNumber); // physical frame number - } - else - { - DCMSEG_DEBUG("Assigning to same new frame bucket"); - // Create new vector - OFVector vec; - vec.push_back(m_framePositions[j].m_frameNumber); - m_logicalFramePositions.push_back(vec); - } - } - } - else - { - std::cerr - << "groupFramesByPosition(): Cannot identify coordinate relevant for sorting frames by position" - << std::endl; - cond = EC_InvalidValue; - } + DCMSEG_ERROR("groupFramesByPosition(): Cannot get Slice Thickness from Pixel Measures FG: " + << cond.text()); + return cond; + } + } + + Uint8 relevantCoordinate = identifyChangingCoordinate(m_imageOrientation); + + // 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 { - std::cerr << "groupFramesByPosition(): Slice Thickness not found, cannot sort frames by position" - << std::endl; - cond = EC_TagNotFound; + DCMSEG_DEBUG("Assigning to same new frame bucket"); + // Create new vector + OFVector vec; + vec.push_back(m_framePositions[j].m_frameNumber); + m_logicalFramePositions.push_back(vec); } } - else - { - std::cerr << "groupFramesByPosition(): Pixel Measures FG not found, cannot sort frames by position" - << std::endl; - cond = EC_TagNotFound; - } + return cond; } @@ -827,4 +767,5 @@ Uint8 OverlapUtil::identifyChangingCoordinate(const OFVector& imageOrie return 3; } -} + +} // namespace dcmqi