diff --git a/Modules/DICOM/src/mitkDICOMTagBasedSorter.cpp b/Modules/DICOM/src/mitkDICOMTagBasedSorter.cpp index 1bef2f93ee..7c8977956b 100644 --- a/Modules/DICOM/src/mitkDICOMTagBasedSorter.cpp +++ b/Modules/DICOM/src/mitkDICOMTagBasedSorter.cpp @@ -1,635 +1,633 @@ /*============================================================================ 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 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, 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(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] = IOVolumeSplitReason::New(); } else { for (auto& [key, value] : listForGroupID) { auto reason = IOVolumeSplitReason::New(); reason->AddReason(IOVolumeSplitReason::ReasonType::ValueSplitDifference); splitReasons[key] = reason; } } return listForGroupID; } mitk::DICOMTagBasedSorter::GroupIDToListType& mitk::DICOMTagBasedSorter ::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[groupKeyStr].push_back(*dataset); if (consecutiveReasons.find(groupKeyStr) == consecutiveReasons.end()) { auto dsReason = splitReasons[gIter->first]->Clone(); if (splitted) dsReason->AddReason(IOVolumeSplitReason::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& [key, group] : consecutiveGroups) { 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 96e25a455d..abca9df4ee 100644 --- a/Modules/DICOM/src/mitkEquiDistantBlocksSorter.cpp +++ b/Modules/DICOM/src/mitkEquiDistantBlocksSorter.cpp @@ -1,623 +1,621 @@ /*============================================================================ 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() : m_SplitReason(IOVolumeSplitReason::New()) { } const mitk::DICOMDatasetList& mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::GetBlockDatasets() const { return m_GroupedFiles; } const mitk::DICOMDatasetList& mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::GetUnsortedDatasets() const { return m_UnsortedFiles; } const mitk::IOVolumeSplitReason* mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult::GetSplitReason() const { return m_SplitReason; } mitk::IOVolumeSplitReason* 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; - m_SliceGroupingResults.clear(); while (!remainingInput.empty()) // repeat until all files are grouped somehow { 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 remainingInput = regularBlock->GetUnsortedDatasets(); if (remainingInput.empty() && !m_SliceGroupingResults.empty() && m_SliceGroupingResults.back()->GetSplitReason()->HasReason(IOVolumeSplitReason::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(IOVolumeSplitReason::ReasonType::OverlappingSlices); } m_SliceGroupingResults.push_back( regularBlock ); } unsigned int numberOfOutputs = m_SliceGroupingResults.size(); this->SetNumberOfOutputs(numberOfOutputs); unsigned int outputIndex(0); for (auto oIter = m_SliceGroupingResults.cbegin(); oIter != m_SliceGroupingResults.cend(); ++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(); } std::shared_ptr mitk::EquiDistantBlocksSorter ::AnalyzeFileForITKImageSeriesReaderSpacingAssumption( const DICOMDatasetList& datasets, bool groupImagesWithGantryTilt) { 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 { // ==> this is a group of its own, stop processing, come back later result->AddFileToSortedBlock( *dsIter ); DICOMDatasetList remainingFiles; remainingFiles.insert( remainingFiles.end(), dsIter+1, datasets.end() ); result->AddFilesToUnsortedBlock( remainingFiles ); if (!remainingFiles.empty()) //if there are remaining files add a split reason result->GetSplitReason()->AddReason(IOVolumeSplitReason::ReasonType::ImagePostionMissing); fileFitsIntoPattern = false; break; // no files anymore } else { // ==> this does not match, consider later 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->GetSplitReason()->AddReason(IOVolumeSplitReason::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) // 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() ); 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->GetSplitReason()->AddReason(IOVolumeSplitReason::ReasonType::GantryTiltDifference); fileFitsIntoPattern = false; } } else // not sheared { 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 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(IOVolumeSplitReason::ReasonType::SliceDistanceInconsistency, std::to_string(fromLastToThisOriginDistance)); if (missingSlicesCount == 0) { result->GetSplitReason()->RemoveReason(IOVolumeSplitReason::ReasonType::MissingSlices); } else if (missingSlicesCount < 0) { result->GetSplitReason()->AddReason(IOVolumeSplitReason::ReasonType::OverlappingSlices); } else if (!result->GetSplitReason()->HasReason(IOVolumeSplitReason::ReasonType::OverlappingSlices)) { //If the missing slice count is positive, but no overlapping was flagged, add the missing slice reason. //We only do it if overlapping was not flagged, to avoid false positives, that could be triggered by slices //of the overlapping volume. result->GetSplitReason()->AddReason(IOVolumeSplitReason::ReasonType::MissingSlices, std::to_string(missingSlicesCount)); } fileFitsIntoPattern = false; } else { result->AddFileToSortedBlock( *dsIter ); // this file is good for current block result->GetSplitReason()->RemoveReason(IOVolumeSplitReason::ReasonType::SliceDistanceInconsistency); result->GetSplitReason()->RemoveReason(IOVolumeSplitReason::ReasonType::MissingSlices); fileFitsIntoPattern = true; missingSlicesCount = 0; } } else // this should be the very first slice { 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() ) { // 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 ) { result->UndoPrematureGrouping(); result->GetSplitReason()->AddReason(IOVolumeSplitReason::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 ) { try { 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 )); } catch (...) { // just do not flag anything, we are ok } } return result; }