diff --git a/Modules/DICOMReader/Resources/configurations/3D/imagetime.xml b/Modules/DICOMReader/Resources/configurations/3D/imagetime.xml index 4de03a0d7a..d705ba630b 100644 --- a/Modules/DICOMReader/Resources/configurations/3D/imagetime.xml +++ b/Modules/DICOMReader/Resources/configurations/3D/imagetime.xml @@ -1,24 +1,24 @@ diff --git a/Modules/DICOMReader/mitkClassicDICOMSeriesReader.cpp b/Modules/DICOMReader/mitkClassicDICOMSeriesReader.cpp index 786e4d9ef0..4a0959f260 100644 --- a/Modules/DICOMReader/mitkClassicDICOMSeriesReader.cpp +++ b/Modules/DICOMReader/mitkClassicDICOMSeriesReader.cpp @@ -1,96 +1,99 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkClassicDICOMSeriesReader.h" #include "mitkDICOMTagBasedSorter.h" #include "mitkDICOMSortByTag.h" #include "mitkSortByImagePositionPatient.h" mitk::ClassicDICOMSeriesReader ::ClassicDICOMSeriesReader() :ThreeDnTDICOMSeriesReader() { mitk::DICOMTagBasedSorter::Pointer tagSorter = mitk::DICOMTagBasedSorter::New(); // all the things that split by tag in mitk::DicomSeriesReader tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x000e) ); // Series Instance UID //tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0052) ); // Frame of Reference UID // a sorter... mitk::DICOMSortCriterion::ConstPointer sorting = mitk::SortByImagePositionPatient::New( // image position patient and image orientation mitk::DICOMSortByTag::New( DICOMTag(0x0020, 0x0012), // aqcuisition number mitk::DICOMSortByTag::New( DICOMTag(0x0008, 0x0032), // aqcuisition time mitk::DICOMSortByTag::New( DICOMTag(0x0018, 0x1060), // trigger time mitk::DICOMSortByTag::New( DICOMTag(0x0008, 0x0018) // SOP instance UID (last resort, not really meaningful but decides clearly) ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer(); tagSorter->SetSortCriterion( sorting ); + tagSorter->SetStrictSorting(false); // nothing did enforce constant distances before, there was just the EquiDistantBlocksSorter logic // define above sorting for this class this->AddSortingElement( tagSorter ); - this->SetFixTiltByShearing(true); - this->SetToleratedOriginOffset(0.005); - this->SetGroup3DandT(true); + this->SetFixTiltByShearing(true); // TODO that was configurable, default was true + this->SetToleratedOriginOffset(0.005); // was hard-coded + this->SetGroup3DandT(true); // TODO that was configurable, default was true + + m_EquiDistantBlocksSorter->SetAcceptTwoSlicesGroups(false); // old reader did not accept that this->SetConfigurationLabel("2013 sorting logic"); this->SetConfigurationDescription("Sort by Image Position, then Acquisition Number, Time, Trigger time, group by 3D+t, group tilted images"); } mitk::ClassicDICOMSeriesReader ::ClassicDICOMSeriesReader(const ClassicDICOMSeriesReader& other ) :itk::Object() ,ThreeDnTDICOMSeriesReader(other) { } mitk::ClassicDICOMSeriesReader ::~ClassicDICOMSeriesReader() { } mitk::ClassicDICOMSeriesReader& mitk::ClassicDICOMSeriesReader ::operator=(const ClassicDICOMSeriesReader& other) { if (this != &other) { ThreeDnTDICOMSeriesReader::operator=(other); } return *this; } bool mitk::ClassicDICOMSeriesReader ::operator==(const DICOMFileReader& other) const { if (dynamic_cast(&other)) { return true; } else { return false; } } diff --git a/Modules/DICOMReader/mitkDICOMFileReaderSelector.cpp b/Modules/DICOMReader/mitkDICOMFileReaderSelector.cpp index e45b9daec5..d084d40bb5 100644 --- a/Modules/DICOMReader/mitkDICOMFileReaderSelector.cpp +++ b/Modules/DICOMReader/mitkDICOMFileReaderSelector.cpp @@ -1,244 +1,244 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDICOMFileReaderSelector.h" #include "mitkDICOMReaderConfigurator.h" #include "mitkModuleContext.h" #include #include #include #include mitk::DICOMFileReaderSelector ::DICOMFileReaderSelector() { } mitk::DICOMFileReaderSelector ::~DICOMFileReaderSelector() { } std::list mitk::DICOMFileReaderSelector ::GetAllConfiguredReaders() const { return m_Readers; } void mitk::DICOMFileReaderSelector ::AddConfigsFromResources(const std::string& path) { std::vector configs = GetModuleContext()->GetModule()->FindResources(path, "*.xml", false); for (std::vector::iterator iter = configs.begin(); iter != configs.end(); ++iter) { ModuleResource& resource = *iter; if (resource.IsValid()) { ModuleResourceStream stream(resource); // read all into string s std::string s; stream.seekg(0, std::ios::end); s.reserve(stream.tellg()); stream.seekg(0, std::ios::beg); s.assign((std::istreambuf_iterator(stream)), std::istreambuf_iterator()); this->AddConfig(s); } } } void mitk::DICOMFileReaderSelector ::AddConfigFromResource(ModuleResource& resource) { if (resource.IsValid()) { ModuleResourceStream stream(resource); // read all into string s std::string s; stream.seekg(0, std::ios::end); s.reserve(stream.tellg()); stream.seekg(0, std::ios::beg); s.assign((std::istreambuf_iterator(stream)), std::istreambuf_iterator()); this->AddConfig(s); } } void mitk::DICOMFileReaderSelector ::AddConfigFromResource(const std::string& resourcename) { ModuleResource r = GetModuleContext()->GetModule()->GetResource(resourcename); this->AddConfigFromResource(r); } void mitk::DICOMFileReaderSelector ::AddFileReaderCanditate(DICOMFileReader::Pointer reader) { if (reader.IsNotNull()) { m_Readers.push_back( reader ); } } void mitk::DICOMFileReaderSelector ::LoadBuiltIn3DConfigs() { //this->AddConfigsFromResources("configurations/3D"); // in this order of preference... this->AddConfigFromResource("configurations/3D/instancenumber.xml"); + this->AddConfigFromResource("configurations/3D/slicelocation.xml"); this->AddConfigFromResource("configurations/3D/imageposition.xml"); this->AddConfigFromResource("configurations/3D/imageposition_byacquisition.xml"); - this->AddConfigFromResource("configurations/3D/slicelocation.xml"); - this->AddConfigFromResource("configurations/3D/imagetime.xml"); - this->AddConfigFromResource("configurations/3D/classicreader.xml"); + //this->AddConfigFromResource("configurations/3D/imagetime.xml"); // no sense in this one? want to see a real-world example first + //this->AddConfigFromResource("configurations/3D/classicreader.xml"); // currently is 3D+t actually } void mitk::DICOMFileReaderSelector ::LoadBuiltIn3DnTConfigs() { this->AddConfigsFromResources("configurations/3DnT"); } void mitk::DICOMFileReaderSelector ::AddConfig(const std::string& xmlDescription) { DICOMReaderConfigurator::Pointer configurator = DICOMReaderConfigurator::New(); DICOMFileReader::Pointer reader = configurator->CreateFromUTF8ConfigString(xmlDescription); if (reader.IsNotNull()) { m_Readers.push_back( reader ); m_PossibleConfigurations.push_back(xmlDescription); } else { std::stringstream ss; ss << "Could not parse reader configuration. Ignoring it."; throw std::invalid_argument( ss.str() ); } } void mitk::DICOMFileReaderSelector ::AddConfigFile(const std::string& filename) { std::ifstream file(filename.c_str()); std::string s; file.seekg(0, std::ios::end); s.reserve(file.tellg()); file.seekg(0, std::ios::beg); s.assign((std::istreambuf_iterator(file)), std::istreambuf_iterator()); this->AddConfig(s); } void mitk::DICOMFileReaderSelector ::SetInputFiles(StringList filenames) { m_InputFilenames = filenames; } const mitk::StringList& mitk::DICOMFileReaderSelector ::GetInputFiles() const { return m_InputFilenames; } mitk::DICOMFileReader::Pointer mitk::DICOMFileReaderSelector ::GetFirstReaderWithMinimumNumberOfOutputImages() { ReaderList workingCandidates; // let all readers analyze the file set unsigned int readerIndex(0); for (ReaderList::iterator rIter = m_Readers.begin(); rIter != m_Readers.end(); ++readerIndex, ++rIter) { (*rIter)->SetInputFiles( m_InputFilenames ); try { (*rIter)->AnalyzeInputFiles(); workingCandidates.push_back( *rIter ); MITK_INFO << "Reader " << readerIndex << " (" << (*rIter)->GetConfigurationLabel() << ") suggests " << (*rIter)->GetNumberOfOutputs() << " 3D blocks"; if ((*rIter)->GetNumberOfOutputs() == 1) { MITK_DEBUG << "Early out with reader #" << readerIndex << " (" << (*rIter)->GetConfigurationLabel() << "), less than 1 block is not possible"; return *rIter; } } catch (std::exception& e) { MITK_ERROR << "Reader " << readerIndex << " (" << (*rIter)->GetConfigurationLabel() << ") threw exception during file analysis, ignoring this reader. Exception: " << e.what(); } catch (...) { MITK_ERROR << "Reader " << readerIndex << " (" << (*rIter)->GetConfigurationLabel() << ") threw unknown exception during file analysis, ignoring this reader."; } } DICOMFileReader::Pointer bestReader; unsigned int minimumNumberOfOutputs = std::numeric_limits::max(); readerIndex = 0; unsigned int bestReaderIndex(0); // select the reader with the minimum number of mitk::Images as output for (ReaderList::iterator rIter = workingCandidates.begin(); rIter != workingCandidates.end(); ++readerIndex, ++rIter) { unsigned int thisReadersNumberOfOutputs = (*rIter)->GetNumberOfOutputs(); if ( thisReadersNumberOfOutputs > 0 // we don't count readers that don't actually produce output && thisReadersNumberOfOutputs < minimumNumberOfOutputs ) { minimumNumberOfOutputs = (*rIter)->GetNumberOfOutputs(); bestReader = *rIter; bestReaderIndex = readerIndex; } } MITK_DEBUG << "Decided for reader #" << bestReaderIndex << " (" << bestReader->GetConfigurationLabel() << ")"; MITK_DEBUG << m_PossibleConfigurations[bestReaderIndex]; return bestReader; } diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp index 00553c0449..ea4598ccac 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp @@ -1,690 +1,705 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "mitkDICOMITKSeriesGDCMReader.h" #include "mitkITKDICOMSeriesReaderHelper.h" #include "mitkGantryTiltInformation.h" #include "mitkDICOMTagBasedSorter.h" #include #include #include mitk::DICOMITKSeriesGDCMReader ::DICOMITKSeriesGDCMReader(unsigned int decimalPlacesForOrientation) :DICOMFileReader() ,m_FixTiltByShearing(true) ,m_DecimalPlacesForOrientation(decimalPlacesForOrientation) { this->EnsureMandatorySortersArePresent(decimalPlacesForOrientation); } mitk::DICOMITKSeriesGDCMReader ::DICOMITKSeriesGDCMReader(const DICOMITKSeriesGDCMReader& other ) :itk::Object() ,DICOMFileReader(other) ,DICOMTagCache(other) ,m_FixTiltByShearing(false) ,m_Sorter( other.m_Sorter ) // TODO should clone the list items ,m_EquiDistantBlocksSorter( other.m_EquiDistantBlocksSorter->Clone() ) ,m_NormalDirectionConsistencySorter( other.m_NormalDirectionConsistencySorter->Clone() ) ,m_DecimalPlacesForOrientation(other.m_DecimalPlacesForOrientation) { } mitk::DICOMITKSeriesGDCMReader ::~DICOMITKSeriesGDCMReader() { } mitk::DICOMITKSeriesGDCMReader& mitk::DICOMITKSeriesGDCMReader ::operator=(const DICOMITKSeriesGDCMReader& other) { if (this != &other) { DICOMFileReader::operator=(other); this->m_FixTiltByShearing = other.m_FixTiltByShearing; this->m_Sorter = other.m_Sorter; // TODO should clone the list items this->m_EquiDistantBlocksSorter = other.m_EquiDistantBlocksSorter->Clone(); this->m_NormalDirectionConsistencySorter = other.m_NormalDirectionConsistencySorter->Clone(); this->m_DecimalPlacesForOrientation = other.m_DecimalPlacesForOrientation; } return *this; } bool mitk::DICOMITKSeriesGDCMReader ::operator==(const DICOMFileReader& other) const { if (const Self* otherSelf = dynamic_cast(&other)) { if ( this->m_FixTiltByShearing == otherSelf->m_FixTiltByShearing && *(this->m_EquiDistantBlocksSorter) == *(otherSelf->m_EquiDistantBlocksSorter) && (fabs(this->m_DecimalPlacesForOrientation - otherSelf->m_DecimalPlacesForOrientation) < eps) ) { // test sorters for equality if (this->m_Sorter.size() != otherSelf->m_Sorter.size()) return false; SorterList::const_iterator mySorterIter = this->m_Sorter.begin(); SorterList::const_iterator oSorterIter = otherSelf->m_Sorter.begin(); for(; mySorterIter != this->m_Sorter.end() && oSorterIter != otherSelf->m_Sorter.end(); ++mySorterIter, ++oSorterIter) { if ( ! (**mySorterIter == **oSorterIter ) ) return false; // this sorter differs } - + // nothing differs ==> all is equal return true; } else { return false; } } else { return false; } } void mitk::DICOMITKSeriesGDCMReader ::SetFixTiltByShearing(bool on) { m_FixTiltByShearing = on; } bool mitk::DICOMITKSeriesGDCMReader ::GetFixTiltByShearing() const { return m_FixTiltByShearing; } +void +mitk::DICOMITKSeriesGDCMReader +::SetAcceptTwoSlicesGroups(bool accept) +{ + m_EquiDistantBlocksSorter->SetAcceptTwoSlicesGroups(accept); +} + +bool +mitk::DICOMITKSeriesGDCMReader +::GetAcceptTwoSlicesGroups() const +{ + return m_EquiDistantBlocksSorter->GetAcceptTwoSlicesGroups(); +} + + mitk::DICOMGDCMImageFrameList mitk::DICOMITKSeriesGDCMReader ::FromDICOMDatasetList(const DICOMDatasetList& input) { DICOMGDCMImageFrameList output; output.reserve(input.size()); for(DICOMDatasetList::const_iterator inputIter = input.begin(); inputIter != input.end(); ++inputIter) { DICOMGDCMImageFrameInfo* gfi = dynamic_cast(*inputIter); assert(gfi); output.push_back(gfi); } return output; } mitk::DICOMDatasetList mitk::DICOMITKSeriesGDCMReader ::ToDICOMDatasetList(const DICOMGDCMImageFrameList& input) { DICOMDatasetList output; output.reserve(input.size()); for(DICOMGDCMImageFrameList::const_iterator inputIter = input.begin(); inputIter != input.end(); ++inputIter) { DICOMDatasetAccess* da = inputIter->GetPointer(); assert(da); output.push_back(da); } return output; } mitk::DICOMImageFrameList mitk::DICOMITKSeriesGDCMReader ::ToDICOMImageFrameList(const DICOMGDCMImageFrameList& input) { DICOMImageFrameList output; output.reserve(input.size()); for(DICOMGDCMImageFrameList::const_iterator inputIter = input.begin(); inputIter != input.end(); ++inputIter) { DICOMImageFrameInfo::Pointer fi = (*inputIter)->GetFrameInfo(); assert(fi.IsNotNull()); output.push_back(fi); } return output; } void mitk::DICOMITKSeriesGDCMReader ::InternalPrintConfiguration(std::ostream& os) const { unsigned int sortIndex(1); for(SorterList::const_iterator sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sortIndex, ++sorterIter) { os << "Sorting step " << sortIndex << ":" << std::endl; (*sorterIter)->PrintConfiguration(os, " "); } os << "Sorting step " << sortIndex << ":" << std::endl; m_EquiDistantBlocksSorter->PrintConfiguration(os, " "); } std::string mitk::DICOMITKSeriesGDCMReader ::GetActiveLocale() const { return setlocale(LC_NUMERIC, NULL); } void mitk::DICOMITKSeriesGDCMReader ::PushLocale() const { std::string currentCLocale = setlocale(LC_NUMERIC, NULL); m_ReplacedCLocales.push( currentCLocale ); setlocale(LC_NUMERIC, "C"); std::locale currentCinLocale( std::cin.getloc() ); m_ReplacedCinLocales.push( currentCinLocale ); std::locale l( "C" ); std::cin.imbue(l); } void mitk::DICOMITKSeriesGDCMReader ::PopLocale() const { if (!m_ReplacedCLocales.empty()) { setlocale(LC_NUMERIC, m_ReplacedCLocales.top().c_str()); m_ReplacedCLocales.pop(); } else { MITK_WARN << "Mismatched PopLocale on DICOMITKSeriesGDCMReader."; } if (!m_ReplacedCinLocales.empty()) { std::cin.imbue( m_ReplacedCinLocales.top() ); m_ReplacedCinLocales.pop(); } else { MITK_WARN << "Mismatched PopLocale on DICOMITKSeriesGDCMReader."; } } mitk::DICOMITKSeriesGDCMReader::SortingBlockList mitk::DICOMITKSeriesGDCMReader ::Condense3DBlocks(SortingBlockList& input) { return input; // to be implemented differently by sub-classes } void mitk::DICOMITKSeriesGDCMReader ::AnalyzeInputFiles() { itk::TimeProbesCollectorBase timer; timer.Start("Reset"); this->ClearOutputs(); m_InputFrameList.clear(); m_GDCMScanner.ClearTags(); timer.Stop("Reset"); // prepare initial sorting (== list of input files) StringList inputFilenames = this->GetInputFiles(); timer.Start("Check appropriateness of input files"); if ( inputFilenames.empty() || !this->CanHandleFile( inputFilenames.front() ) // first || !this->CanHandleFile( inputFilenames.back() ) // last || !this->CanHandleFile( inputFilenames[ inputFilenames.size() / 2] ) // roughly central file ) { // TODO a read-as-many-as-possible fallback could be implemented here MITK_DEBUG << "Reader unable to process files.."; return; } timer.Stop("Check appropriateness of input files"); // scan files for sorting-relevant tags timer.Start("Setup scanning"); for(SorterList::iterator sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sorterIter) { assert(sorterIter->IsNotNull()); DICOMTagList tags = (*sorterIter)->GetTagsOfInterest(); for(DICOMTagList::const_iterator tagIter = tags.begin(); tagIter != tags.end(); ++tagIter) { MITK_DEBUG << "Sorting uses tag " << tagIter->GetGroup() << "," << tagIter->GetElement(); m_GDCMScanner.AddTag( gdcm::Tag(tagIter->GetGroup(), tagIter->GetElement()) ); } } // Add some of our own interest // TODO all tags that are needed in DICOMImageBlockDescriptor should be added by DICOMFileReader (this code location here should query all superclasses for tags) m_GDCMScanner.AddTag( gdcm::Tag(0x0018,0x1164) ); // pixel spacing m_GDCMScanner.AddTag( gdcm::Tag(0x0028,0x0030) ); // imager pixel spacing m_GDCMScanner.AddTag( gdcm::Tag(0x0028,0x1050) ); // window center m_GDCMScanner.AddTag( gdcm::Tag(0x0028,0x1051) ); // window width m_GDCMScanner.AddTag( gdcm::Tag(0x0008,0x0008) ); // image type m_GDCMScanner.AddTag( gdcm::Tag(0x0028,0x0004) ); // photometric interpretation m_GDCMScanner.AddTag( gdcm::Tag(0x0020,0x1041) ); // slice location m_GDCMScanner.AddTag( gdcm::Tag(0x0020,0x0013) ); // instance number m_GDCMScanner.AddTag( gdcm::Tag(0x0008,0x0016) ); // sop class UID m_GDCMScanner.AddTag( gdcm::Tag(0x0008,0x0018) ); // sop instance UID m_GDCMScanner.AddTag( gdcm::Tag(0x0020,0x0011) ); // series number m_GDCMScanner.AddTag( gdcm::Tag(0x0008,0x1030) ); // study description m_GDCMScanner.AddTag( gdcm::Tag(0x0008,0x103e) ); // series description m_GDCMScanner.AddTag( gdcm::Tag(0x0008,0x0060) ); // modality m_GDCMScanner.AddTag( gdcm::Tag(0x0020,0x0012) ); // acquisition number m_GDCMScanner.AddTag( gdcm::Tag(0x0018,0x0024) ); // sequence name m_GDCMScanner.AddTag( gdcm::Tag(0x0020,0x0037) ); // image orientation m_GDCMScanner.AddTag( gdcm::Tag(0x0020,0x0032) ); // ipp timer.Stop("Setup scanning"); timer.Start("Tag scanning"); PushLocale(); m_GDCMScanner.Scan( inputFilenames ); PopLocale(); timer.Stop("Tag scanning"); timer.Start("Setup sorting"); for (StringList::const_iterator inputIter = inputFilenames.begin(); inputIter != inputFilenames.end(); ++inputIter) { m_InputFrameList.push_back( DICOMGDCMImageFrameInfo::New( DICOMImageFrameInfo::New(*inputIter, 0), m_GDCMScanner.GetMapping(inputIter->c_str()) ) ); } m_SortingResultInProgress.clear(); m_SortingResultInProgress.push_back( m_InputFrameList ); timer.Stop("Setup sorting"); // sort and split blocks as configured timer.Start("Sorting frames"); unsigned int sorterIndex = 0; for(SorterList::iterator sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sorterIndex, ++sorterIter) { m_SortingResultInProgress = this->InternalExecuteSortingStep(sorterIndex, *sorterIter, m_SortingResultInProgress, &timer); } // a last extra-sorting step: ensure equidistant slices m_SortingResultInProgress = this->InternalExecuteSortingStep(sorterIndex++, m_EquiDistantBlocksSorter.GetPointer(), m_SortingResultInProgress, &timer); timer.Stop("Sorting frames"); timer.Start("Condensing 3D blocks (3D+t or vector values)"); m_SortingResultInProgress = this->Condense3DBlocks( m_SortingResultInProgress ); timer.Stop("Condensing 3D blocks (3D+t or vector values)"); // provide final result as output timer.Start("Output"); unsigned int o = this->GetNumberOfOutputs(); this->SetNumberOfOutputs( o + m_SortingResultInProgress.size() ); // Condense3DBlocks may already have added outputs! for (SortingBlockList::iterator blockIter = m_SortingResultInProgress.begin(); blockIter != m_SortingResultInProgress.end(); ++o, ++blockIter) { DICOMGDCMImageFrameList& gdcmFrameInfoList = *blockIter; assert(!gdcmFrameInfoList.empty()); // reverse frames if necessary // update tilt information from absolute last sorting DICOMDatasetList datasetList = ToDICOMDatasetList( gdcmFrameInfoList ); m_NormalDirectionConsistencySorter->SetInput( datasetList ); m_NormalDirectionConsistencySorter->Sort(); DICOMGDCMImageFrameList sortedGdcmInfoFrameList = FromDICOMDatasetList( m_NormalDirectionConsistencySorter->GetOutput(0) ); const GantryTiltInformation& tiltInfo = m_NormalDirectionConsistencySorter->GetTiltInformation(); // set frame list for current block DICOMImageFrameList frameList = ToDICOMImageFrameList( sortedGdcmInfoFrameList ); assert(!frameList.empty()); DICOMImageBlockDescriptor block; block.SetTagCache( this ); // important: this must be before SetImageFrameList(), because SetImageFrameList will trigger reading of lots of interesting tags! block.SetImageFrameList( frameList ); block.SetTiltInformation( tiltInfo ); block.SetReaderImplementationLevel( this->GetReaderImplementationLevel( block.GetSOPClassUID() ) ); this->SetOutput( o, block ); } timer.Stop("Output"); #ifdef MBILOG_ENABLE_DEBUG std::cout << "---------------------------------------------------------------" << std::endl; timer.Report( std::cout ); std::cout << "---------------------------------------------------------------" << std::endl; #endif } mitk::DICOMITKSeriesGDCMReader::SortingBlockList mitk::DICOMITKSeriesGDCMReader ::InternalExecuteSortingStep( unsigned int sortingStepIndex, DICOMDatasetSorter::Pointer sorter, const SortingBlockList& input, itk::TimeProbesCollectorBase* timer) { SortingBlockList nextStepSorting; // we should not modify our input list while processing it std::stringstream ss; ss << "Sorting step " << sortingStepIndex << " '"; sorter->PrintConfiguration(ss); ss << "'"; timer->Start( ss.str().c_str() ); nextStepSorting.clear(); MITK_DEBUG << "================================================================================"; MITK_DEBUG << "DICOMITKSeriesGDCMReader: " << ss.str() << ": " << input.size() << " groups input"; unsigned int groupIndex = 0; for(SortingBlockList::const_iterator blockIter = input.begin(); blockIter != input.end(); ++groupIndex, ++blockIter) { const DICOMGDCMImageFrameList& gdcmInfoFrameList = *blockIter; DICOMDatasetList datasetList = ToDICOMDatasetList( gdcmInfoFrameList ); MITK_DEBUG << "--------------------------------------------------------------------------------"; MITK_DEBUG << "DICOMITKSeriesGDCMReader: " << ss.str() << ", dataset group " << groupIndex << " (" << datasetList.size() << " datasets): "; for (DICOMDatasetList::iterator oi = datasetList.begin(); oi != datasetList.end(); ++oi) { MITK_DEBUG << " INPUT : " << (*oi)->GetFilenameIfAvailable(); } sorter->SetInput(datasetList); sorter->Sort(); unsigned int numberOfResultingBlocks = sorter->GetNumberOfOutputs(); for (unsigned int b = 0; b < numberOfResultingBlocks; ++b) { DICOMDatasetList blockResult = sorter->GetOutput(b); for (DICOMDatasetList::iterator oi = blockResult.begin(); oi != blockResult.end(); ++oi) { MITK_DEBUG << " OUTPUT(" << b << ") :" << (*oi)->GetFilenameIfAvailable(); } DICOMGDCMImageFrameList sortedGdcmInfoFrameList = FromDICOMDatasetList(blockResult); nextStepSorting.push_back( sortedGdcmInfoFrameList ); } } timer->Stop( ss.str().c_str() ); return nextStepSorting; } mitk::ReaderImplementationLevel mitk::DICOMITKSeriesGDCMReader ::GetReaderImplementationLevel(const std::string sopClassUID) const { if (sopClassUID.empty()) { return SOPClassUnknown; } gdcm::UIDs uidKnowledge; uidKnowledge.SetFromUID( sopClassUID.c_str() ); gdcm::UIDs::TSType gdcmType = uidKnowledge; switch (gdcmType) { case gdcm::UIDs::CTImageStorage: case gdcm::UIDs::MRImageStorage: case gdcm::UIDs::PositronEmissionTomographyImageStorage: case gdcm::UIDs::ComputedRadiographyImageStorage: case gdcm::UIDs::DigitalXRayImageStorageForPresentation: case gdcm::UIDs::DigitalXRayImageStorageForProcessing: return SOPClassSupported; case gdcm::UIDs::NuclearMedicineImageStorage: return SOPClassPartlySupported; case gdcm::UIDs::SecondaryCaptureImageStorage: return SOPClassImplemented; default: return SOPClassUnsupported; } } // void AllocateOutputImages(); bool mitk::DICOMITKSeriesGDCMReader ::LoadImages() { bool success = true; unsigned int numberOfOutputs = this->GetNumberOfOutputs(); for (unsigned int o = 0; o < numberOfOutputs; ++o) { success &= this->LoadMitkImageForOutput(o); } return success; } bool mitk::DICOMITKSeriesGDCMReader ::LoadMitkImageForImageBlockDescriptor(DICOMImageBlockDescriptor& block) const { PushLocale(); const DICOMImageFrameList& frames = block.GetImageFrameList(); const GantryTiltInformation tiltInfo = block.GetTiltInformation(); bool hasTilt = tiltInfo.IsRegularGantryTilt(); ITKDICOMSeriesReaderHelper::StringContainer filenames; for (DICOMImageFrameList::const_iterator frameIter = frames.begin(); frameIter != frames.end(); ++frameIter) { filenames.push_back( (*frameIter)->Filename ); } mitk::ITKDICOMSeriesReaderHelper helper; bool success(true); try { mitk::Image::Pointer mitkImage = helper.Load( filenames, m_FixTiltByShearing && hasTilt, tiltInfo ); block.SetMitkImage( mitkImage ); } catch (std::exception& e) { success = false; MITK_ERROR << "Exception during image loading: " << e.what(); } PopLocale(); return success; } bool mitk::DICOMITKSeriesGDCMReader ::LoadMitkImageForOutput(unsigned int o) { DICOMImageBlockDescriptor& block = this->InternalGetOutput(o); return this->LoadMitkImageForImageBlockDescriptor(block); } bool mitk::DICOMITKSeriesGDCMReader ::CanHandleFile(const std::string& filename) { return ITKDICOMSeriesReaderHelper::CanHandleFile(filename); } void mitk::DICOMITKSeriesGDCMReader ::AddSortingElement(DICOMDatasetSorter* sorter, bool atFront) { assert(sorter); if (atFront) { m_Sorter.push_front( sorter ); } else { m_Sorter.push_back( sorter ); } } mitk::DICOMITKSeriesGDCMReader::ConstSorterList mitk::DICOMITKSeriesGDCMReader ::GetFreelyConfiguredSortingElements() const { std::list result; unsigned int sortIndex(0); for(SorterList::const_iterator sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sortIndex, ++sorterIter) { if (sortIndex > 0) // ignore first element (see EnsureMandatorySortersArePresent) { result.push_back( (*sorterIter).GetPointer() ); } } return result; } void mitk::DICOMITKSeriesGDCMReader ::EnsureMandatorySortersArePresent(unsigned int decimalPlacesForOrientation) { DICOMTagBasedSorter::Pointer splitter = DICOMTagBasedSorter::New(); splitter->AddDistinguishingTag( DICOMTag(0x0028, 0x0010) ); // Number of Rows splitter->AddDistinguishingTag( DICOMTag(0x0028, 0x0011) ); // Number of Columns splitter->AddDistinguishingTag( DICOMTag(0x0028, 0x0030) ); // Pixel Spacing splitter->AddDistinguishingTag( DICOMTag(0x0018, 0x1164) ); // Imager Pixel Spacing splitter->AddDistinguishingTag( DICOMTag(0x0020, 0x0037), new mitk::DICOMTagBasedSorter::CutDecimalPlaces(decimalPlacesForOrientation) ); // Image Orientation (Patient) splitter->AddDistinguishingTag( DICOMTag(0x0018, 0x0050) ); // Slice Thickness splitter->AddDistinguishingTag( DICOMTag(0x0028, 0x0008) ); // Number of Frames this->AddSortingElement( splitter, true ); // true = at front if (m_EquiDistantBlocksSorter.IsNull()) { m_EquiDistantBlocksSorter = mitk::EquiDistantBlocksSorter::New(); } m_EquiDistantBlocksSorter->SetAcceptTilt( m_FixTiltByShearing ); if (m_NormalDirectionConsistencySorter.IsNull()) { m_NormalDirectionConsistencySorter = mitk::NormalDirectionConsistencySorter::New(); } } void mitk::DICOMITKSeriesGDCMReader ::SetToleratedOriginOffsetToAdaptive(double fractionOfInterSliceDistance) { assert( m_EquiDistantBlocksSorter.IsNotNull() ); m_EquiDistantBlocksSorter->SetToleratedOriginOffsetToAdaptive(fractionOfInterSliceDistance); } void mitk::DICOMITKSeriesGDCMReader ::SetToleratedOriginOffset(double millimeters) { assert( m_EquiDistantBlocksSorter.IsNotNull() ); m_EquiDistantBlocksSorter->SetToleratedOriginOffset(millimeters); } double mitk::DICOMITKSeriesGDCMReader ::GetToleratedOriginError() const { assert( m_EquiDistantBlocksSorter.IsNotNull() ); return m_EquiDistantBlocksSorter->GetToleratedOriginOffset(); } bool mitk::DICOMITKSeriesGDCMReader ::IsToleratedOriginOffsetAbsolute() const { assert( m_EquiDistantBlocksSorter.IsNotNull() ); return m_EquiDistantBlocksSorter->IsToleratedOriginOffsetAbsolute(); } double mitk::DICOMITKSeriesGDCMReader ::GetDecimalPlacesForOrientation() const { return m_DecimalPlacesForOrientation; } std::string mitk::DICOMITKSeriesGDCMReader ::GetTagValue(DICOMImageFrameInfo* frame, const DICOMTag& tag) const { // TODO inefficient. if (m_InputFrameList.contains(frame)) return frame->GetTagValueAsString(tag); for(DICOMGDCMImageFrameList::const_iterator frameIter = m_InputFrameList.begin(); frameIter != m_InputFrameList.end(); ++frameIter) { if ( (*frameIter)->GetFrameInfo().IsNotNull() && (*((*frameIter)->GetFrameInfo()) == *frame ) ) { return (*frameIter)->GetTagValueAsString(tag); } } return ""; } diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h index ee3975c3af..42a8a46adf 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h @@ -1,351 +1,357 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDICOMITKSeriesGDCMReader_h #define mitkDICOMITKSeriesGDCMReader_h #include "mitkDICOMFileReader.h" #include "mitkDICOMDatasetSorter.h" #include "mitkDICOMGDCMImageFrameInfo.h" #include "mitkEquiDistantBlocksSorter.h" #include "mitkNormalDirectionConsistencySorter.h" #include "mitkITKDICOMSeriesReaderHelper.h" #include "DICOMReaderExports.h" #include namespace itk { class TimeProbesCollectorBase; } namespace mitk { /** \ingroup DICOMReaderModule \brief Flexible reader based on itk::ImageSeriesReader and GDCM, for single-slice modalities like CT, MR, PET, CR, etc. Implements the loading processed as structured by DICOMFileReader offers configuration of its loading strategy. Documentation sections: - \ref DICOMITKSeriesGDCMReader_LoadingStrategy - \ref DICOMITKSeriesGDCMReader_ForcedConfiguration - \ref DICOMITKSeriesGDCMReader_UserConfiguration - \ref DICOMITKSeriesGDCMReader_GantryTilt - \ref DICOMITKSeriesGDCMReader_Testing - \ref DICOMITKSeriesGDCMReader_Internals - \ref DICOMITKSeriesGDCMReader_RelatedClasses - \ref DICOMITKSeriesGDCMReader_TiltInternals - \ref DICOMITKSeriesGDCMReader_Condensing \section DICOMITKSeriesGDCMReader_LoadingStrategy Loading strategy The set of input files is processed by a number of DICOMDatasetSorter objects which may do two sort of things: 1. split a list of input frames into multiple lists, based on DICOM tags such as "Rows", "Columns", which cannot be mixed within a single mitk::Image 2. sort the frames within the input lists, based on the values of DICOM tags such as "Image Position Patient" When the DICOMITKSeriesGDCMReader is configured with DICOMDatasetSorter%s, the list of input files is processed as follows: 1. build an initial set of output groups, simply by grouping all input files. 2. for each configured DICOMDatasetSorter, process: - for each output group: 1. set this group's files as input to the sorter 2. let the sorter sort (and split) 3. integrate the sorter's output groups with our own output groups \section DICOMITKSeriesGDCMReader_ForcedConfiguration Forced Configuration In all cases, the reader will add two DICOMDatasetSorter objects that are required to load mitk::Images properly via itk::ImageSeriesReader: 1. As a \b first step, the input files will be split into groups that are not compatible because they differ in essential aspects: - (0028,0010) Number of Rows - (0028,0011) Number of Columns - (0028,0030) Pixel Spacing - (0018,1164) Imager Pixel Spacing - (0020,0037) %Image Orientation (Patient) - (0018,0050) Slice Thickness - (0028,0008) Number of Frames 2. As are two forced \b last steps: 1. There will always be an instance of EquiDistantBlocksSorter, which ensures that there is an equal distance between all the frames of an Image. This is required to achieve correct geometrical positions in the mitk::Image, i.e. it is essential to be able to make measurements in images. - whether or not the distance is required to be orthogonal to the image planes is configured by SetFixTiltByShearing(). - during this check, we need to tolerate some minor errors in documented vs. calculated image origins. The amount of tolerance can be adjusted by SetToleratedOriginOffset() and SetToleratedOriginOffsetToAdaptive(). Please see EquiDistantBlocksSorter for more details. The default should be good for most cases. 2. There is always an instance of NormalDirectionConsistencySorter, which makes the order of images go along the image normals (see NormalDirectionConsistencySorter) \section DICOMITKSeriesGDCMReader_UserConfiguration User Configuration The user of this class can add more sorting steps (similar to the one described in above section) by calling AddSortingElement(). Usually, an application will add sorting by "Image Position Patient", by "Instance Number", and by other relevant tags here. \section DICOMITKSeriesGDCMReader_GantryTilt Gantry tilt handling When CT gantry tilt is used, the gantry plane (= X-Ray source and detector ring) and the vertical plane do not align anymore. This scanner feature is used for example to reduce metal artifacs (e.g. Lee C , Evaluation of Using CT Gantry Tilt Scan on Head and Neck Cancer Patients with Dental Structure: Scans Show Less Metal Artifacts. Presented at: Radiological Society of North America 2011 Scientific Assembly and Annual Meeting; November 27- December 2, 2011 Chicago IL.). The acquired planes of such CT series do not match the expectations of a orthogonal geometry in mitk::Image: if you stack the slices, they show a small shift along the Y axis: \verbatim without tilt with tilt |||||| ////// |||||| ////// -- |||||| --------- ////// -------- table orientation |||||| ////// |||||| ////// Stacked slices: without tilt with tilt -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- \endverbatim As such gemetries do not "work" in conjunction with mitk::Image, DICOMITKSeriesGDCMReader is able to perform a correction for such series. Whether or not such correction should be attempted is controlled by SetFixTiltByShearing(), the default being correction. For details, see "Internals" below. \section DICOMITKSeriesGDCMReader_Testing Testing A number of tests is implemented in module DICOMTesting, which is documented at \ref DICOMTesting. \section DICOMITKSeriesGDCMReader_Internals Class internals Internally, the class is based on GDCM and it depends heavily on the gdcm::Scanner class. Since the sorting elements (see DICOMDatasetSorter and DICOMSortCriterion) can access tags only via the DICOMDatasetAccess interface, BUT DICOMITKSeriesGDCMReader holds a list of more specific classes DICOMGDCMImageFrameInfo, we must convert between the two types sometimes. This explains the methods ToDICOMDatasetList(), FromDICOMDatasetList(). The intermediate result of all the sorting efforts is held in m_SortingResultInProgress, which is modified through InternalExecuteSortingStep(). \subsection DICOMITKSeriesGDCMReader_RelatedClasses Overview of related classes The following diagram gives an overview of the related classes: \image html implementeditkseriesgdcmreader.jpg \subsection DICOMITKSeriesGDCMReader_TiltInternals Details about the tilt correction The gantry tilt "correction" algorithm fixes two errors introduced by ITK's ImageSeriesReader: - the plane shift that is ignored by ITK's reader is recreated by applying a shearing transformation using itk::ResampleFilter. - the spacing is corrected (it is calculated by ITK's reader from the distance between two origins, which is NOT the slice distance in this special case) Both errors are introduced in itkImageSeriesReader.txx (ImageSeriesReader::GenerateOutputInformation(void)), lines 176 to 245 (as of ITK 3.20) For the correction, we examine two consecutive slices of a series, both described as a pair (origin/orientation): - we calculate if the first origin is on a line along the normal of the second slice - if this is not the case, the geometry will not fit a normal mitk::Image/mitk::Geometry3D - we then project the second origin into the first slice's coordinate system to quantify the shift - both is done in class GantryTiltInformation with quite some comments. The geometry of image stacks with tilted geometries is illustrated below: - green: the DICOM images as described by their tags: origin as a point with the line indicating the orientation - red: the output of ITK ImageSeriesReader: wrong, larger spacing, no tilt - blue: how much a shear must correct \image html tilt-correction.jpg \subsection DICOMITKSeriesGDCMReader_Condensing Sub-classes can condense multiple blocks into a single larger block The sorting/splitting process described above is helpful for at least two more DICOM readers, which either try to load 3D+t images or which load diffusion data. In both cases, a single pixel of the mitk::Image is made up of multiple values, in one case values over time, in the other case multiple measurements of a single point. The specialized readers for these cases (e.g. ThreeDnTDICOMSeriesReader) can reuse most of the methods in DICOMITKSeriesGDCMReader, except that they need an extra step after the usual sorting, in which they can merge already grouped 3D blocks. What blocks are merged depends on the specialized reader's understanding of these images. To allow for such merging, a method Condense3DBlocks() is called as an absolute last step of AnalyzeInputFiles(). Given this, a sub-class could implement only LoadImages() and Condense3DBlocks() instead repeating most of AnalyzeInputFiles(). */ class DICOMReader_EXPORT DICOMITKSeriesGDCMReader : public DICOMFileReader, protected DICOMTagCache { public: mitkClassMacro( DICOMITKSeriesGDCMReader, DICOMFileReader ); mitkCloneMacro( DICOMITKSeriesGDCMReader ); itkNewMacro( DICOMITKSeriesGDCMReader ); mitkNewMacro1Param( DICOMITKSeriesGDCMReader, unsigned int ); /** \brief Runs the sorting / splitting process described in \ref DICOMITKSeriesGDCMReader_LoadingStrategy. Method required by DICOMFileReader. */ virtual void AnalyzeInputFiles(); // void AllocateOutputImages(); /** \brief Loads images using itk::ImageSeriesReader, potentially applies shearing to correct gantry tilt. */ virtual bool LoadImages(); // re-implemented from super-class virtual bool CanHandleFile(const std::string& filename); /** \brief Add an element to the sorting procedure described in \ref DICOMITKSeriesGDCMReader_LoadingStrategy. */ virtual void AddSortingElement(DICOMDatasetSorter* sorter, bool atFront = false); - + typedef const std::list ConstSorterList; ConstSorterList GetFreelyConfiguredSortingElements() const; /** \brief Controls whether to "fix" tilted acquisitions by shearing the output (see \ref DICOMITKSeriesGDCMReader_GantryTilt). */ void SetFixTiltByShearing(bool on); bool GetFixTiltByShearing() const; + /** + \brief Controls whether groups of only two images are accepted when ensuring consecutive slices via EquiDistantBlocksSorter. + */ + void SetAcceptTwoSlicesGroups(bool accept); + bool GetAcceptTwoSlicesGroups() const; + /** \brief See \ref DICOMITKSeriesGDCMReader_ForcedConfiguration. */ void SetToleratedOriginOffsetToAdaptive(double fractionOfInterSliceDistanct = 0.3); /** \brief See \ref DICOMITKSeriesGDCMReader_ForcedConfiguration. */ void SetToleratedOriginOffset(double millimeters = 0.005); double GetToleratedOriginError() const; bool IsToleratedOriginOffsetAbsolute() const; double GetDecimalPlacesForOrientation() const; virtual bool operator==(const DICOMFileReader& other) const; protected: virtual void InternalPrintConfiguration(std::ostream& os) const; // From DICOMTagCache virtual std::string GetTagValue(DICOMImageFrameInfo* frame, const DICOMTag& tag) const; /// \brief Return active C locale std::string GetActiveLocale() const; /** \brief Remember current locale on stack, activate "C" locale. "C" locale is required for correct parsing of numbers by itk::ImageSeriesReader */ void PushLocale() const; /** \brief Activate last remembered locale from locale stack "C" locale is required for correct parsing of numbers by itk::ImageSeriesReader */ void PopLocale() const; DICOMITKSeriesGDCMReader(unsigned int decimalPlacesForOrientation = 5); virtual ~DICOMITKSeriesGDCMReader(); DICOMITKSeriesGDCMReader(const DICOMITKSeriesGDCMReader& other); DICOMITKSeriesGDCMReader& operator=(const DICOMITKSeriesGDCMReader& other); /// \brief See \ref DICOMITKSeriesGDCMReader_Internals DICOMDatasetList ToDICOMDatasetList(const DICOMGDCMImageFrameList& input); /// \brief See \ref DICOMITKSeriesGDCMReader_Internals DICOMGDCMImageFrameList FromDICOMDatasetList(const DICOMDatasetList& input); /// \brief See \ref DICOMITKSeriesGDCMReader_Internals DICOMImageFrameList ToDICOMImageFrameList(const DICOMGDCMImageFrameList& input); typedef std::list SortingBlockList; /** \brief "Hook" for sub-classes, see \ref DICOMITKSeriesGDCMReader_Condensing \return REMAINING blocks */ virtual SortingBlockList Condense3DBlocks(SortingBlockList& resultOf3DGrouping); /// \brief Sorting step as described in \ref DICOMITKSeriesGDCMReader_LoadingStrategy SortingBlockList InternalExecuteSortingStep( unsigned int sortingStepIndex, DICOMDatasetSorter::Pointer sorter, const SortingBlockList& input, itk::TimeProbesCollectorBase* timer); /// \brief Loads the mitk::Image by means of an itk::ImageSeriesReader virtual bool LoadMitkImageForOutput(unsigned int o); virtual bool LoadMitkImageForImageBlockDescriptor(DICOMImageBlockDescriptor& block) const; /** \brief Shear the loaded mitk::Image to "correct" a spatial error introduced by itk::ImageSeriesReader See \ref DICOMITKSeriesGDCMReader_GantryTilt for details. */ Image::Pointer FixupSpacing(Image* mitkImage, const DICOMImageBlockDescriptor& block) const; /// \brief Describe this reader's confidence for given SOP class UID ReaderImplementationLevel GetReaderImplementationLevel(const std::string sopClassUID) const; private: /// \brief Creates the required sorting steps described in \ref DICOMITKSeriesGDCMReader_ForcedConfiguration void EnsureMandatorySortersArePresent(unsigned int decimalPlacesForOrientation); protected: // NOT nice, made available to ThreeDnTDICOMSeriesReader due to lack of time bool m_FixTiltByShearing; // could be removed by ITKDICOMSeriesReader NOT flagging tilt unless requested to fix it! private: SortingBlockList m_SortingResultInProgress; typedef std::list SorterList; SorterList m_Sorter; - mitk::EquiDistantBlocksSorter::Pointer m_EquiDistantBlocksSorter; - protected: - // NOT nice, made available to ThreeDnTDICOMSeriesReader due to lack of time + // NOT nice, made available to ThreeDnTDICOMSeriesReader and ClassicDICOMSeriesReader due to lack of time + mitk::EquiDistantBlocksSorter::Pointer m_EquiDistantBlocksSorter; + mitk::NormalDirectionConsistencySorter::Pointer m_NormalDirectionConsistencySorter; private: mutable std::stack m_ReplacedCLocales; mutable std::stack m_ReplacedCinLocales; DICOMGDCMImageFrameList m_InputFrameList; gdcm::Scanner m_GDCMScanner; - + double m_DecimalPlacesForOrientation; }; } #endif diff --git a/Modules/DICOMReader/mitkDICOMReaderConfigurator.cpp b/Modules/DICOMReader/mitkDICOMReaderConfigurator.cpp index f0b6da45ed..915404e70f 100644 --- a/Modules/DICOMReader/mitkDICOMReaderConfigurator.cpp +++ b/Modules/DICOMReader/mitkDICOMReaderConfigurator.cpp @@ -1,637 +1,661 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDICOMReaderConfigurator.h" #include "mitkDICOMSortByTag.h" #include "mitkSortByImagePositionPatient.h" mitk::DICOMReaderConfigurator ::DICOMReaderConfigurator() { } mitk::DICOMReaderConfigurator ::~DICOMReaderConfigurator() { } mitk::DICOMFileReader::Pointer mitk::DICOMReaderConfigurator ::CreateFromConfigFile(const std::string& filename) const { TiXmlDocument doc (filename); if (doc.LoadFile()) { return this->CreateFromTiXmlDocument( doc ); } else { MITK_ERROR << "Unable to load file at '" << filename <<"'"; return DICOMFileReader::Pointer(); } } mitk::DICOMFileReader::Pointer mitk::DICOMReaderConfigurator ::CreateFromUTF8ConfigString(const std::string& xmlContents) const { TiXmlDocument doc; doc.Parse(xmlContents.c_str(), 0, TIXML_ENCODING_UTF8); return this->CreateFromTiXmlDocument( doc ); } mitk::DICOMFileReader::Pointer mitk::DICOMReaderConfigurator ::CreateFromTiXmlDocument(TiXmlDocument& doc) const { TiXmlHandle root(doc.RootElement()); if (TiXmlElement* rootElement = root.ToElement()) { if (strcmp(rootElement->Value(), "DICOMFileReader")) // :-( no std::string methods { MITK_ERROR << "File should contain a tag at top-level! Found '" << (rootElement->Value() ? std::string(rootElement->Value()) : std::string("!nothing!")) << "' instead"; return NULL; } const char* classnameC = rootElement->Attribute("class"); if (!classnameC) { MITK_ERROR << "File should name a reader class in the class attribute: . Found nothing instead"; return NULL; } int version(1); if ( rootElement->QueryIntAttribute("version", &version) == TIXML_SUCCESS) { if (version != 1) { MITK_WARN << "This reader is only capable of creating DICOMFileReaders of version 1. " << "Will not continue, because given configuration is meant for version " << version << "."; return NULL; } } else { MITK_ERROR << "File should name the version of the reader class in the version attribute: ." << " Found nothing instead, assuming version 1!"; version = 1; } std::string classname(classnameC); double decimalPlacesForOrientation(5); bool useDecimalPlacesForOrientation(false); useDecimalPlacesForOrientation = rootElement->QueryDoubleAttribute("decimalPlacesForOrientation", &decimalPlacesForOrientation) == TIXML_SUCCESS; // attribute present and a double value if (classname == "ClassicDICOMSeriesReader") { mitk::ClassicDICOMSeriesReader::Pointer reader = mitk::ClassicDICOMSeriesReader::New(); return reader.GetPointer(); } if (classname == "ThreeDnTDICOMSeriesReader") { mitk::ThreeDnTDICOMSeriesReader::Pointer reader; if (useDecimalPlacesForOrientation) reader = mitk::ThreeDnTDICOMSeriesReader::New(decimalPlacesForOrientation); else reader = mitk::ThreeDnTDICOMSeriesReader::New(); return ConfigureThreeDnTDICOMSeriesReader(reader, rootElement).GetPointer(); } else if (classname == "DICOMITKSeriesGDCMReader") { mitk::DICOMITKSeriesGDCMReader::Pointer reader; if (useDecimalPlacesForOrientation) reader = mitk::DICOMITKSeriesGDCMReader::New(decimalPlacesForOrientation); else reader = mitk::DICOMITKSeriesGDCMReader::New(); return ConfigureDICOMITKSeriesGDCMReader(reader, rootElement).GetPointer(); } else { MITK_ERROR << "DICOMFileReader tag names unknown class '" << classname << "'"; return NULL; } } else { MITK_ERROR << "Great confusion: no root element in XML document. Expecting a DICOMFileReader tag at top-level."; return NULL; } } #define boolStringTrue(s) \ ( s == "true" || s == "on" || s == "1" \ || s == "TRUE" || s == "ON") mitk::ThreeDnTDICOMSeriesReader::Pointer mitk::DICOMReaderConfigurator ::ConfigureThreeDnTDICOMSeriesReader(ThreeDnTDICOMSeriesReader::Pointer reader, TiXmlElement* element) const { assert(element); // use all the base class configuration if (this->ConfigureDICOMITKSeriesGDCMReader( reader.GetPointer(), element ).IsNull()) { return NULL; } // add the "group3DnT" flag bool group3DnT(true); const char* group3DnTC = element->Attribute("group3DnT"); if (group3DnTC) { std::string group3DnTS(group3DnTC); group3DnT = boolStringTrue(group3DnTS); } reader->SetGroup3DandT( group3DnT ); return reader; } mitk::DICOMITKSeriesGDCMReader::Pointer mitk::DICOMReaderConfigurator ::ConfigureDICOMITKSeriesGDCMReader(DICOMITKSeriesGDCMReader::Pointer reader, TiXmlElement* element) const { assert(element); const char* configLabelC = element->Attribute("label"); if (configLabelC) { std::string configLabel(configLabelC); reader->SetConfigurationLabel(configLabel); } const char* configDescriptionC = element->Attribute("description"); if (configDescriptionC) { std::string configDescription(configDescriptionC); reader->SetConfigurationDescription(configDescriptionC); } // "fixTiltByShearing" flag bool fixTiltByShearing(false); const char* fixTiltByShearingC = element->Attribute("fixTiltByShearing"); if (fixTiltByShearingC) { std::string fixTiltByShearingS(fixTiltByShearingC); fixTiltByShearing = boolStringTrue(fixTiltByShearingS); } reader->SetFixTiltByShearing( fixTiltByShearing ); + // "acceptTwoSlicesGroups" flag + bool acceptTwoSlicesGroups(true); + const char* acceptTwoSlicesGroupsC = element->Attribute("acceptTwoSlicesGroups"); + if (acceptTwoSlicesGroupsC) + { + std::string acceptTwoSlicesGroupsS(acceptTwoSlicesGroupsC); + acceptTwoSlicesGroups = boolStringTrue(acceptTwoSlicesGroupsS); + } + + reader->SetAcceptTwoSlicesGroups( acceptTwoSlicesGroups ); + // "toleratedOriginError" attribute (double) bool toleratedOriginErrorIsAbsolute(false); const char* toleratedOriginErrorIsAbsoluteC = element->Attribute("toleratedOriginErrorIsAbsolute"); if (toleratedOriginErrorIsAbsoluteC) { std::string toleratedOriginErrorIsAbsoluteS(toleratedOriginErrorIsAbsoluteC); toleratedOriginErrorIsAbsolute = boolStringTrue(toleratedOriginErrorIsAbsoluteS); } double toleratedOriginError(0.3); if (element->QueryDoubleAttribute("toleratedOriginError", &toleratedOriginError) == TIXML_SUCCESS) // attribute present and a double value { if (toleratedOriginErrorIsAbsolute) { reader->SetToleratedOriginOffset( toleratedOriginError ); } else { reader->SetToleratedOriginOffsetToAdaptive( toleratedOriginError ); } } // DICOMTagBasedSorters are the only thing we create at this point // TODO for-loop over all child elements of type DICOMTagBasedSorter, BUT actually a single sorter of this type is enough. TiXmlElement* dElement = element->FirstChildElement("DICOMDatasetSorter"); if (dElement) { const char* classnameC = dElement->Attribute("class"); if (!classnameC) { MITK_ERROR << "File should name a DICOMDatasetSorter class in the class attribute of . Found nothing instead"; return NULL; } std::string classname(classnameC); if (classname == "DICOMTagBasedSorter") { DICOMTagBasedSorter::Pointer tagSorter = CreateDICOMTagBasedSorter(dElement); if (tagSorter.IsNotNull()) { reader->AddSortingElement( tagSorter ); } } else { MITK_ERROR << "DICOMDatasetSorter tag names unknown class '" << classname << "'"; return NULL; } } return reader; } mitk::DICOMTagBasedSorter::Pointer mitk::DICOMReaderConfigurator ::CreateDICOMTagBasedSorter(TiXmlElement* element) const { mitk::DICOMTagBasedSorter::Pointer tagSorter = mitk::DICOMTagBasedSorter::New(); + + // "strictSorting" parameter! + bool strictSorting(true); + const char* strictSortingC = element->Attribute("strictSorting"); + if (strictSortingC) + { + std::string strictSortingS(strictSortingC); + strictSorting = boolStringTrue(strictSortingS); + } + tagSorter->SetStrictSorting(strictSorting); + TiXmlElement* dElement = element->FirstChildElement("Distinguishing"); if (dElement) { for ( TiXmlElement* tChild = dElement->FirstChildElement(); tChild != NULL; tChild = tChild->NextSiblingElement() ) { try { mitk::DICOMTag tag = tagFromXMLElement(tChild); int i(5); if (tChild->QueryIntAttribute("cutDecimalPlaces", &i) == TIXML_SUCCESS) { tagSorter->AddDistinguishingTag( tag, new mitk::DICOMTagBasedSorter::CutDecimalPlaces(i) ); } else { tagSorter->AddDistinguishingTag( tag ); } } catch(...) { return NULL; } } } // "sorting tags" TiXmlElement* sElement = element->FirstChildElement("Sorting"); if (sElement) { DICOMSortCriterion::Pointer previousCriterion; DICOMSortCriterion::Pointer currentCriterion; for ( TiXmlNode* tChildNode = sElement->LastChild(); tChildNode != NULL; tChildNode = tChildNode->PreviousSibling() ) { TiXmlElement* tChild = tChildNode->ToElement(); if (!tChild) continue; if (!strcmp(tChild->Value(), "Tag")) { try { currentCriterion = this->CreateDICOMSortByTag(tChild, previousCriterion); } catch(...) { std::stringstream ss; ss << "Could not parse element at (input line " << tChild->Row() << ", col. " << tChild->Column() << ")!"; MITK_ERROR << ss.str(); return NULL; } } else if (!strcmp(tChild->Value(), "ImagePositionPatient")) { try { currentCriterion = this->CreateSortByImagePositionPatient(tChild, previousCriterion); } catch(...) { std::stringstream ss; ss << "Could not parse element at (input line " << tChild->Row() << ", col. " << tChild->Column() << ")!"; MITK_ERROR << ss.str(); return NULL; } } else { MITK_ERROR << "File contain unknown tag <" << tChild->Value() << "> tag as child to ! Cannot interpret..."; } previousCriterion = currentCriterion; } tagSorter->SetSortCriterion( currentCriterion.GetPointer() ); } return tagSorter; } std::string mitk::DICOMReaderConfigurator ::requiredStringAttribute(TiXmlElement* xmlElement, const std::string& key) const { assert(xmlElement); const char* gC = xmlElement->Attribute(key.c_str()); if (gC) { std::string gS(gC); return gS; } else { std::stringstream ss; ss << "Expected an attribute '" << key << "' at this position " "(input line " << xmlElement->Row() << ", col. " << xmlElement->Column() << ")!"; MITK_ERROR << ss.str(); throw std::invalid_argument( ss.str() ); } } unsigned int mitk::DICOMReaderConfigurator ::hexStringToUInt(const std::string& s) const { std::stringstream converter(s); unsigned int ui; converter >> std::hex >> ui; MITK_DEBUG << "Converted string '" << s << "' to unsigned int " << ui; return ui; } mitk::DICOMTag mitk::DICOMReaderConfigurator ::tagFromXMLElement(TiXmlElement* xmlElement) const { assert(xmlElement); if (strcmp(xmlElement->Value(), "Tag")) // :-( no std::string methods { std::stringstream ss; ss << "Expected a tag at this position " "(input line " << xmlElement->Row() << ", col. " << xmlElement->Column() << ")!"; MITK_ERROR << ss.str(); throw std::invalid_argument( ss.str() ); } std::string name = requiredStringAttribute(xmlElement, "name"); std::string groupS = requiredStringAttribute(xmlElement, "group"); std::string elementS = requiredStringAttribute(xmlElement, "element"); try { // convert string to int (assuming string is in hex format with leading "0x" like "0x0020") unsigned int group = hexStringToUInt(groupS); unsigned int element = hexStringToUInt(elementS); return DICOMTag(group, element); } catch(...) { std::stringstream ss; ss << "Expected group and element values in to be hexadecimal with leading 0x, e.g. '0x0020'" "(input line " << xmlElement->Row() << ", col. " << xmlElement->Column() << ")!"; MITK_ERROR << ss.str(); throw std::invalid_argument( ss.str() ); } } mitk::DICOMSortCriterion::Pointer mitk::DICOMReaderConfigurator ::CreateDICOMSortByTag(TiXmlElement* xmlElement, DICOMSortCriterion::Pointer secondaryCriterion) const { mitk::DICOMTag tag = tagFromXMLElement(xmlElement); return DICOMSortByTag::New(tag, secondaryCriterion).GetPointer(); } mitk::DICOMSortCriterion::Pointer mitk::DICOMReaderConfigurator ::CreateSortByImagePositionPatient(TiXmlElement*, DICOMSortCriterion::Pointer secondaryCriterion) const { return SortByImagePositionPatient::New(secondaryCriterion).GetPointer(); } std::string mitk::DICOMReaderConfigurator ::CreateConfigStringFromReader(DICOMFileReader::ConstPointer reader) const { // check possible sub-classes from the most-specific one up to the most generic one const DICOMFileReader* cPointer = reader; TiXmlElement* root; if (const ClassicDICOMSeriesReader* specificReader = dynamic_cast(cPointer)) { root = this->CreateConfigStringFromReader(specificReader); } else if (const ThreeDnTDICOMSeriesReader* specificReader = dynamic_cast(cPointer)) { root = this->CreateConfigStringFromReader(specificReader); } else if (const DICOMITKSeriesGDCMReader* specificReader = dynamic_cast(cPointer)) { root = this->CreateConfigStringFromReader(specificReader); } else { MITK_WARN << "Unknown reader class passed to DICOMReaderConfigurator::CreateConfigStringFromReader(). Cannot serialize."; return ""; // no serialization, what a pity } if (root) { TiXmlDocument document; document.LinkEndChild( root ); TiXmlPrinter printer; printer.SetIndent( " " ); document.Accept( &printer ); std::string xmltext = printer.CStr(); return xmltext; } else { MITK_WARN << "DICOMReaderConfigurator::CreateConfigStringFromReader() created empty serialization. Problem?"; return ""; } } TiXmlElement* mitk::DICOMReaderConfigurator ::CreateConfigStringFromReader(const DICOMITKSeriesGDCMReader* reader) const { TiXmlElement* root = this->CreateDICOMFileReaderTag(reader); assert(root); root->SetAttribute("fixTiltByShearing", toString(reader->GetFixTiltByShearing())); + root->SetAttribute("acceptTwoSlicesGroups", toString(reader->GetAcceptTwoSlicesGroups())); root->SetDoubleAttribute("toleratedOriginError", reader->GetToleratedOriginError()); root->SetAttribute("toleratedOriginErrorIsAbsolute", toString(reader->IsToleratedOriginOffsetAbsolute())); root->SetDoubleAttribute("decimalPlacesForOrientation", reader->GetDecimalPlacesForOrientation()); // iterate DICOMDatasetSorter objects DICOMITKSeriesGDCMReader::ConstSorterList sorterList = reader->GetFreelyConfiguredSortingElements(); for(DICOMITKSeriesGDCMReader::ConstSorterList::const_iterator sorterIter = sorterList.begin(); sorterIter != sorterList.end(); ++sorterIter) { const DICOMDatasetSorter* sorter = *sorterIter; if (const DICOMTagBasedSorter* specificSorter = dynamic_cast(sorter)) { TiXmlElement* sorterTag = this->CreateConfigStringFromDICOMDatasetSorter(specificSorter); root->LinkEndChild(sorterTag); } else { MITK_WARN << "Unknown DICOMDatasetSorter class passed to DICOMReaderConfigurator::CreateConfigStringFromReader(). Cannot serialize."; return NULL; } } return root; } TiXmlElement* mitk::DICOMReaderConfigurator ::CreateConfigStringFromDICOMDatasetSorter(const DICOMTagBasedSorter* sorter) const { assert(sorter); TiXmlElement* sorterTag = new TiXmlElement("DICOMDatasetSorter"); sorterTag->SetAttribute("class", sorter->GetNameOfClass()); + sorterTag->SetAttribute("strictSorting", toString(sorter->GetStrictSorting())); TiXmlElement* distinguishingTagsElement = new TiXmlElement("Distinguishing"); sorterTag->LinkEndChild(distinguishingTagsElement); mitk::DICOMTagList distinguishingTags = sorter->GetDistinguishingTags(); for (DICOMTagList::iterator tagIter = distinguishingTags.begin(); tagIter != distinguishingTags.end(); ++tagIter) { TiXmlElement* tag = this->CreateConfigStringFromDICOMTag(*tagIter); distinguishingTagsElement->LinkEndChild(tag); const DICOMTagBasedSorter::TagValueProcessor* processor = sorter->GetTagValueProcessorForDistinguishingTag(*tagIter); if (const DICOMTagBasedSorter::CutDecimalPlaces* specificProcessor = dynamic_cast(processor)) { tag->SetDoubleAttribute("cutDecimalPlaces", specificProcessor->GetPrecision()); } } TiXmlElement* sortingElement = new TiXmlElement("Sorting"); sorterTag->LinkEndChild(sortingElement); mitk::DICOMSortCriterion::ConstPointer sortCriterion = sorter->GetSortCriterion(); while (sortCriterion.IsNotNull()) { std::string classname = sortCriterion->GetNameOfClass(); if (classname == "SortByImagePositionPatient") { sortingElement->LinkEndChild( new TiXmlElement("ImagePositionPatient") ); // no parameters } else if (classname == "DICOMSortByTag") { DICOMTagList pseudoTagList = sortCriterion->GetTagsOfInterest(); if (pseudoTagList.size() == 1) { DICOMTag firstTag = pseudoTagList.front(); TiXmlElement* tagElement = this->CreateConfigStringFromDICOMTag(firstTag); sortingElement->LinkEndChild( tagElement ); } else { MITK_ERROR << "Encountered SortByTag class with MULTIPLE tag in CreateConfigStringFromDICOMDatasetSorter. Cannot serialize."; return NULL; } } else { MITK_ERROR << "Encountered unknown class '" << classname << "' in CreateConfigStringFromDICOMDatasetSorter. Cannot serialize."; return NULL; } sortCriterion = sortCriterion->GetSecondaryCriterion(); } return sorterTag; } TiXmlElement* mitk::DICOMReaderConfigurator ::CreateConfigStringFromDICOMTag(const DICOMTag& tag) const { TiXmlElement* tagElement = new TiXmlElement("Tag"); // name group element tagElement->SetAttribute("name", tag.GetName().c_str()); tagElement->SetAttribute("group", toHexString(tag.GetGroup())); tagElement->SetAttribute("element", toHexString(tag.GetElement())); return tagElement; } std::string mitk::DICOMReaderConfigurator ::toHexString(unsigned int i) const { std::stringstream ss; ss << "0x" << std::setfill ('0') << std::setw(4) << std::hex << i; return ss.str(); } TiXmlElement* mitk::DICOMReaderConfigurator ::CreateConfigStringFromReader(const ThreeDnTDICOMSeriesReader* reader) const { TiXmlElement* root = this->CreateConfigStringFromReader(static_cast(reader)); assert(root); root->SetAttribute("group3DnT", toString(reader->GetGroup3DandT())); return root; } const char* mitk::DICOMReaderConfigurator ::toString(bool b) const { return b ? "true" : "false"; } TiXmlElement* mitk::DICOMReaderConfigurator ::CreateConfigStringFromReader(const ClassicDICOMSeriesReader* reader) const { return this->CreateDICOMFileReaderTag(reader); } TiXmlElement* mitk::DICOMReaderConfigurator ::CreateDICOMFileReaderTag(const DICOMFileReader* reader) const { TiXmlElement* readerTag = new TiXmlElement("DICOMFileReader"); readerTag->SetAttribute("class", reader->GetNameOfClass()); readerTag->SetAttribute("label", reader->GetConfigurationLabel().c_str()); readerTag->SetAttribute("description", reader->GetConfigurationDescription().c_str()); readerTag->SetAttribute("version", "1"); return readerTag; } diff --git a/Modules/DICOMReader/mitkDICOMSortByTag.cpp b/Modules/DICOMReader/mitkDICOMSortByTag.cpp index 7c2e8e8e9a..1f56d826d1 100644 --- a/Modules/DICOMReader/mitkDICOMSortByTag.cpp +++ b/Modules/DICOMReader/mitkDICOMSortByTag.cpp @@ -1,146 +1,175 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDICOMSortByTag.h" mitk::DICOMSortByTag ::DICOMSortByTag(const DICOMTag& tag, DICOMSortCriterion::Pointer secondaryCriterion) :DICOMSortCriterion(secondaryCriterion) ,m_Tag(tag) { } mitk::DICOMSortByTag ::~DICOMSortByTag() { } mitk::DICOMSortByTag ::DICOMSortByTag(const DICOMSortByTag& other ) :DICOMSortCriterion(other) ,m_Tag(other.m_Tag) { } mitk::DICOMSortByTag& mitk::DICOMSortByTag ::operator=(const DICOMSortByTag& other) { if (this != &other) { DICOMSortCriterion::operator=(other); m_Tag = other.m_Tag; } return *this; } - + bool mitk::DICOMSortByTag ::operator==(const DICOMSortCriterion& other) const { if (const DICOMSortByTag* otherSelf = dynamic_cast(&other)) { if (!(this->m_Tag == otherSelf->m_Tag)) return false; if (this->m_SecondaryCriterion.IsNull() && otherSelf->m_SecondaryCriterion.IsNull()) return true; if (this->m_SecondaryCriterion.IsNull() || otherSelf->m_SecondaryCriterion.IsNull()) return false; return *(this->m_SecondaryCriterion) == *(otherSelf->m_SecondaryCriterion); } else { return false; } } void mitk::DICOMSortByTag ::Print(std::ostream& os) const { m_Tag.Print(os); } mitk::DICOMTagList mitk::DICOMSortByTag ::GetTagsOfInterest() const { DICOMTagList list; list.push_back(m_Tag); return list; } bool mitk::DICOMSortByTag ::IsLeftBeforeRight(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right) const { return this->NumericCompare(left, right, m_Tag); } bool mitk::DICOMSortByTag ::StringCompare(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right, const DICOMTag& tag) const { assert(left); assert(right); std::string leftString = left->GetTagValueAsString(tag); std::string rightString = right->GetTagValueAsString(tag); if (leftString != rightString) { return leftString.compare(rightString) < 0; } else { return this->NextLevelIsLeftBeforeRight(left, right); } } bool mitk::DICOMSortByTag ::NumericCompare(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right, const DICOMTag& tag) const { assert(left); assert(right); std::string leftString = left->GetTagValueAsString(tag); std::string rightString = right->GetTagValueAsString(tag); std::istringstream lefti(leftString); std::istringstream righti(rightString); double leftDouble(0); double rightDouble(0); if ( (lefti >> leftDouble) && (righti >> rightDouble) && lefti.eof() && righti.eof() ) { if (leftDouble != rightDouble) // can we decide? { return leftDouble < rightDouble; } else // ask secondary criterion { return this->NextLevelIsLeftBeforeRight(left, right); } } else // no numerical conversion.. { return this->StringCompare(left,right, tag); // fallback to string compare } } + +double +mitk::DICOMSortByTag +::NumericDistance(const mitk::DICOMDatasetAccess* from, const mitk::DICOMDatasetAccess* to) const +{ + assert(from); + assert(to); + + std::string fromString = from->GetTagValueAsString(m_Tag); + std::string toString = to->GetTagValueAsString(m_Tag); + + std::istringstream fromi(fromString); + std::istringstream toi(toString); + + double fromDouble(0); + double toDouble(0); + + if ( (fromi >> fromDouble) && (toi >> toDouble) + && fromi.eof() && toi.eof() ) + { + return toDouble - fromDouble; + } + else + { + return 0; + } + + // TODO second-level compare? +} diff --git a/Modules/DICOMReader/mitkDICOMSortByTag.h b/Modules/DICOMReader/mitkDICOMSortByTag.h index 984a115a4a..27784420ba 100644 --- a/Modules/DICOMReader/mitkDICOMSortByTag.h +++ b/Modules/DICOMReader/mitkDICOMSortByTag.h @@ -1,70 +1,72 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDICOMSortByTag_h #define mitkDICOMSortByTag_h #include "mitkDICOMSortCriterion.h" namespace mitk { /** \ingroup DICOMReaderModule \brief Compare two datasets by the value of a single tag (for use in DICOMTagBasedSorter). The class will compare the tag values by 1. numerical value if possible (i.e. both datasets have a value that is numerical) 2. alphabetical order otherwise If the comparison results in equalness, it is refered to the secondary criterion, see DICOMSortByTag::NextLevelIsLeftBeforeRight(). */ class DICOMReader_EXPORT DICOMSortByTag : public DICOMSortCriterion { public: mitkClassMacro( DICOMSortByTag, DICOMSortCriterion ); mitkNewMacro1Param( DICOMSortByTag, const DICOMTag& ); mitkNewMacro2Param( DICOMSortByTag, const DICOMTag&, DICOMSortCriterion::Pointer ); virtual DICOMTagList GetTagsOfInterest() const; virtual bool IsLeftBeforeRight(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right) const; + virtual double NumericDistance(const mitk::DICOMDatasetAccess* from, const mitk::DICOMDatasetAccess* to) const; + virtual void Print(std::ostream& os) const; - + virtual bool operator==(const DICOMSortCriterion& other) const; protected: DICOMSortByTag( const DICOMTag& tag, DICOMSortCriterion::Pointer secondaryCriterion = NULL ); virtual ~DICOMSortByTag(); DICOMSortByTag(const DICOMSortByTag& other); DICOMSortByTag& operator=(const DICOMSortByTag& other); bool StringCompare(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right, const DICOMTag& tag) const; bool NumericCompare(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right, const DICOMTag& tag) const; private: DICOMTag m_Tag; }; } #endif diff --git a/Modules/DICOMReader/mitkDICOMSortCriterion.h b/Modules/DICOMReader/mitkDICOMSortCriterion.h index 079f7be8f0..0364a44b47 100644 --- a/Modules/DICOMReader/mitkDICOMSortCriterion.h +++ b/Modules/DICOMReader/mitkDICOMSortCriterion.h @@ -1,78 +1,82 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDICOMSortCriterion_h #define mitkDICOMSortCriterion_h #include "itkObjectFactory.h" #include "mitkCommon.h" #include "mitkDICOMDatasetAccess.h" namespace mitk { /** \ingroup DICOMReaderModule \brief A tag based sorting criterion for use in DICOMTagBasedSorter. This class is used within std::sort (see DICOMTagBasedSorter::Sort()) and has to answer a simple question by implementing IsLeftBeforeRight(). Each time IsLeftBeforeRight() is called, the method should return whether the left dataset should be sorted before the right dataset. Because there are identical tags values quite oftenly, a DICOMSortCriterion will always hold a secondary DICOMSortCriterion. In cases of equal tag values, the decision is refered to the secondary criterion. */ class DICOMReader_EXPORT DICOMSortCriterion : public itk::LightObject { public: mitkClassMacro( DICOMSortCriterion, itk::LightObject ); /// \brief Tags used for comparison (includes seconary criteria). DICOMTagList GetAllTagsOfInterest() const; /// \brief Tags used for comparison. virtual DICOMTagList GetTagsOfInterest() const = 0; /// \brief Answer the sorting question. virtual bool IsLeftBeforeRight(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right) const = 0; + /// \brief Calculate a distance between two datasets. + /// This ansers the question of consecutive datasets. + virtual double NumericDistance(const mitk::DICOMDatasetAccess* from, const mitk::DICOMDatasetAccess* to) const = 0; + /// \brief The fallback criterion. DICOMSortCriterion::ConstPointer GetSecondaryCriterion() const; /// brief describe this class in given stream. virtual void Print(std::ostream& os) const = 0; virtual bool operator==(const DICOMSortCriterion& other) const = 0; protected: DICOMSortCriterion( DICOMSortCriterion::Pointer secondaryCriterion ); virtual ~DICOMSortCriterion(); bool NextLevelIsLeftBeforeRight(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right) const; DICOMSortCriterion(const DICOMSortCriterion& other); DICOMSortCriterion& operator=(const DICOMSortCriterion& other); DICOMSortCriterion::Pointer m_SecondaryCriterion; }; } #endif diff --git a/Modules/DICOMReader/mitkDICOMTagBasedSorter.cpp b/Modules/DICOMReader/mitkDICOMTagBasedSorter.cpp index b3a217e051..df1007768d 100644 --- a/Modules/DICOMReader/mitkDICOMTagBasedSorter.cpp +++ b/Modules/DICOMReader/mitkDICOMTagBasedSorter.cpp @@ -1,359 +1,509 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDICOMTagBasedSorter.h" #include +#include mitk::DICOMTagBasedSorter::CutDecimalPlaces ::CutDecimalPlaces(unsigned int precision) :m_Precision(precision) { } std::string mitk::DICOMTagBasedSorter::CutDecimalPlaces ::operator()(const std::string& input) const { // be a bit tolerant for tags such as image orientation orienatation, let only the first few digits matter (http://bugs.mitk.org/show_bug.cgi?id=12263) // iterate all fields, convert each to a number, cut this number as configured, then return a concatenated string with all cut-off numbers static std::ostringstream resultString; resultString.str(std::string()); resultString.clear(); resultString.setf(std::ios::fixed, std::ios::floatfield); resultString.precision(m_Precision); static std::stringstream ss(input); ss.str(input); ss.clear(); static std::string item; while (std::getline(ss, item, '\\')) { static std::istringstream converter(item); converter.str(item); converter.clear(); static double number(0); 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(); } unsigned int mitk::DICOMTagBasedSorter::CutDecimalPlaces ::GetPrecision() const { return m_Precision; } mitk::DICOMTagBasedSorter ::DICOMTagBasedSorter() :DICOMDatasetSorter() +,m_StrictSorting(true) { } mitk::DICOMTagBasedSorter ::~DICOMTagBasedSorter() { for(TagValueProcessorMap::iterator ti = m_TagValueProcessor.begin(); ti != m_TagValueProcessor.end(); ++ti) { delete ti->second; } } mitk::DICOMTagBasedSorter ::DICOMTagBasedSorter(const DICOMTagBasedSorter& other ) :DICOMDatasetSorter(other) +,m_StrictSorting(other.m_StrictSorting) { } mitk::DICOMTagBasedSorter& mitk::DICOMTagBasedSorter ::operator=(const DICOMTagBasedSorter& other) { if (this != &other) { DICOMDatasetSorter::operator=(other); + m_StrictSorting = other.m_StrictSorting; } return *this; } bool mitk::DICOMTagBasedSorter ::operator==(const DICOMDatasetSorter& other) const { if (const DICOMTagBasedSorter* otherSelf = dynamic_cast(&other)) { + if (this->m_StrictSorting != otherSelf->m_StrictSorting) return false; + bool allTagsPresentAndEqual(true); if (this->m_DistinguishingTags.size() != otherSelf->m_DistinguishingTags.size()) return false; for (DICOMTagList::const_iterator myTag = this->m_DistinguishingTags.begin(); myTag != this->m_DistinguishingTags.end(); ++myTag) { allTagsPresentAndEqual &= (std::find( otherSelf->m_DistinguishingTags.begin(), otherSelf->m_DistinguishingTags.end(), *myTag ) != otherSelf->m_DistinguishingTags.end()); // 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:" << std::endl; for (DICOMTagList::const_iterator 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; +} mitk::DICOMTagList mitk::DICOMTagBasedSorter ::GetTagsOfInterest() { DICOMTagList allTags = m_DistinguishingTags; DICOMTagList sortingRelevantTags = m_SortCriterion->GetAllTagsOfInterest(); allTags.insert( allTags.end(), sortingRelevantTags.begin(), sortingRelevantTags.end() ); // append return allTags; } mitk::DICOMTagList mitk::DICOMTagBasedSorter ::GetDistinguishingTags() const { return m_DistinguishingTags; } const mitk::DICOMTagBasedSorter::TagValueProcessor* mitk::DICOMTagBasedSorter ::GetTagValueProcessorForDistinguishingTag(const DICOMTag& tag) const { TagValueProcessorMap::const_iterator loc = m_TagValueProcessor.find(tag); if (loc != m_TagValueProcessor.end()) { return loc->second; } else { return NULL; } } 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() { // 1. split // 2. sort each group GroupIDToListType groups = this->SplitInputGroups(); GroupIDToListType& sortedGroups = this->SortGroups( groups ); // 3. define output this->SetNumberOfOutputs(sortedGroups.size()); unsigned int outputIndex(0); for (GroupIDToListType::iterator groupIter = sortedGroups.begin(); groupIter != sortedGroups.end(); ++outputIndex, ++groupIter) { this->SetOutput(outputIndex, groupIter->second); } } std::string mitk::DICOMTagBasedSorter ::BuildGroupID( DICOMDatasetAccess* dataset ) { // just concatenate all tag values assert(dataset); std::stringstream groupID; groupID << "g"; for (DICOMTagList::iterator tagIter = m_DistinguishingTags.begin(); tagIter != m_DistinguishingTags.end(); ++tagIter) { groupID << tagIter->GetGroup() << tagIter->GetElement(); // make group/element part of the id to cover empty tags std::string rawTagValue = dataset->GetTagValueAsString(*tagIter); std::string processedTagValue; if ( m_TagValueProcessor[*tagIter] != NULL ) { processedTagValue = (*m_TagValueProcessor[*tagIter])(rawTagValue); } else { processedTagValue = rawTagValue; } groupID << processedTagValue; } // shorten ID? return groupID.str(); } mitk::DICOMTagBasedSorter::GroupIDToListType mitk::DICOMTagBasedSorter ::SplitInputGroups() { DICOMDatasetList input = GetInput(); // copy GroupIDToListType listForGroupID; for (DICOMDatasetList::iterator dsIter = input.begin(); dsIter != input.end(); ++dsIter) { DICOMDatasetAccess* dataset = *dsIter; assert(dataset); 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"; return listForGroupID; } mitk::DICOMTagBasedSorter::GroupIDToListType& mitk::DICOMTagBasedSorter ::SortGroups(GroupIDToListType& groups) { 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 // - ... ? unsigned int groupIndex(0); for (GroupIDToListType::iterator gIter = groups.begin(); gIter != groups.end(); ++groupIndex, ++gIter) { DICOMDatasetList& dsList = gIter->second; MITK_DEBUG << " --------------------------------------------------------------------------------"; MITK_DEBUG << " DICOMTagBasedSorter before sorting group : " << groupIndex; for (DICOMDatasetList::iterator oi = dsList.begin(); oi != dsList.end(); ++oi) { MITK_DEBUG << " INPUT : " << (*oi)->GetFilenameIfAvailable(); } std::sort( dsList.begin(), dsList.end(), ParameterizedDatasetSort( m_SortCriterion ) ); MITK_DEBUG << " --------------------------------------------------------------------------------"; MITK_DEBUG << " DICOMTagBasedSorter after sorting group : " << groupIndex; for (DICOMDatasetList::iterator oi = dsList.begin(); oi != dsList.end(); ++oi) { MITK_DEBUG << " OUTPUT : " << (*oi)->GetFilenameIfAvailable(); } MITK_DEBUG << " --------------------------------------------------------------------------------"; } + + GroupIDToListType consecutiveGroups; + if (m_StrictSorting) + { + // Step 2: create new groups by enforcing consecutive order within each group + unsigned int groupIndex(0); + for (GroupIDToListType::iterator gIter = groups.begin(); + gIter != groups.end(); + ++gIter) + { + std::stringstream groupKey; + groupKey << std::setfill('0') << std::setw(6) << groupIndex++; + + DICOMDatasetList& dsList = gIter->second; + DICOMDatasetAccess* previousDS(NULL); + unsigned int dsIndex(0); + double constantDistance(0.0); + bool constantDistanceInitialized(false); + for (DICOMDatasetList::iterator dataset = dsList.begin(); + dataset != dsList.end(); + ++dsIndex, ++dataset) + { + 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! + double currentDistance = m_SortCriterion->NumericDistance(previousDS, *dataset); + if (constantDistanceInitialized) + { + if (fabs(currentDistance - constantDistance) < constantDistance * 0.01) // ok, deviation of up to 1% of distance is tolerated + { + // nothing to do, just ok + } + else if (currentDistance < mitk::eps) // close enough to 0 + { + // no numeric comparison possible? + // rare case(?), just accept + } + else + { + // 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++; + } + } + 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 + if ((currentDistance - (int)currentDistance == 0.0) && fabs(currentDistance) != 1.0) + // exact comparison. An integer should not be expressed as 1.000000000000000000000000001! + { + groupKey.str(std::string()); + groupKey.clear(); + groupKey << std::setfill('0') << std::setw(6) << groupIndex++; + } + + constantDistance = currentDistance; + constantDistanceInitialized = true; + } + } + consecutiveGroups[groupKey.str()].push_back(*dataset); + 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 (GroupIDToListType::iterator gIter = consecutiveGroups.begin(); + gIter != consecutiveGroups.end(); + ++gIter) + { + assert(!gIter->second.empty()); + firstSlices.push_back(gIter->second.front()); + } + + std::sort( firstSlices.begin(), firstSlices.end(), ParameterizedDatasetSort( m_SortCriterion ) ); + + GroupIDToListType sortedResultBlocks; + unsigned int groupKeyValue(0); + for (DICOMDatasetList::iterator firstSlice = firstSlices.begin(); + firstSlice != firstSlices.end(); + ++firstSlice) + { + for (GroupIDToListType::iterator gIter = consecutiveGroups.begin(); + gIter != consecutiveGroups.end(); + ++groupKeyValue, ++gIter) + { + if (gIter->second.front() == *firstSlice) + { + std::stringstream groupKey; + groupKey << std::setfill('0') << std::setw(6) << groupKeyValue; // try more than 999,999 groups and you are doomed (your application already is) + sortedResultBlocks[groupKey.str()] = gIter->second; + } + } + } + + groups = sortedResultBlocks; } 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/DICOMReader/mitkDICOMTagBasedSorter.h b/Modules/DICOMReader/mitkDICOMTagBasedSorter.h index 176e0d49c0..32722b48e1 100644 --- a/Modules/DICOMReader/mitkDICOMTagBasedSorter.h +++ b/Modules/DICOMReader/mitkDICOMTagBasedSorter.h @@ -1,154 +1,162 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDICOMTagBasedSorter_h #define mitkDICOMTagBasedSorter_h #include "mitkDICOMDatasetSorter.h" #include "mitkDICOMSortCriterion.h" namespace mitk { /** \ingroup DICOMReaderModule \brief Sort DICOM datasets based on configurable tags. This class implements sorting of input DICOM datasets into multiple outputs as described in \ref DICOMITKSeriesGDCMReader_LoadingStrategy. The logic of sorting and splitting is most simple and most generic: 1. Datasets will be put into different groups, if they differ in their value of specific tags (defined by AddDistinguishingTag()) - there might be multiple distinguishing tags defined - tag values might be processed before comparison by means of TagValueProcessor (e.g. round to a number of decimal places) 2. Each of the groups will be sorted by comparing their tag values using multiple DICOMSortCriterion - DICOMSortCriterion might evaluate a single tag (e.g. Instance Number) or multiple values (as in SortByImagePositionPatient) - only a single DICOMSortCriterion is defined for DICOMTagBasedSorter, because each DICOMSortCriterion holds a "secondary sort criterion", i.e. an application can define multiple tags for sorting by chaining \link DICOMSortCriterion DICOMSortCriteria \endlink - applications should make sure that sorting is always defined (to avoid problems with standard containers), e.g. by adding a comparison of filenames or instance UIDs as a last sorting fallback. */ class DICOMReader_EXPORT DICOMTagBasedSorter : public DICOMDatasetSorter { public: /** \brief Processes tag values before they are compared. These classes could do some kind of normalization such as rounding, lower case formatting, etc. */ class DICOMReader_EXPORT TagValueProcessor { public: /// \brief Implements the "processing". virtual std::string operator()(const std::string&) const = 0; }; /** \brief Cuts a number after configured number of decimal places. An instance of this class can be used to avoid errors when comparing minimally different image orientations. */ class DICOMReader_EXPORT CutDecimalPlaces : public TagValueProcessor { public: CutDecimalPlaces(unsigned int precision); unsigned int GetPrecision() const; virtual std::string operator()(const std::string&) const; private: unsigned int m_Precision; }; mitkClassMacro( DICOMTagBasedSorter, DICOMDatasetSorter ) itkNewMacro( DICOMTagBasedSorter ) /** \brief Datasets that differ in given tag's value will be sorted into separate outputs. */ void AddDistinguishingTag( const DICOMTag&, TagValueProcessor* tagValueProcessor = NULL ); DICOMTagList GetDistinguishingTags() const; const TagValueProcessor* GetTagValueProcessorForDistinguishingTag(const DICOMTag&) const; /** \brief Define the sorting criterion (which holds seconardy criteria) */ void SetSortCriterion( DICOMSortCriterion::ConstPointer criterion ); DICOMSortCriterion::ConstPointer GetSortCriterion() const; /** \brief A list of all the tags needed for processing (facilitates scanning). */ virtual DICOMTagList GetTagsOfInterest(); + /** + \brief Whether or not groups should be checked for consecutive tag values. + */ + void SetStrictSorting(bool strict); + bool GetStrictSorting() const; + /** \brief Actually sort as described in the Detailed Description. */ virtual void Sort(); /** \brief Print configuration details into given stream. */ virtual void PrintConfiguration(std::ostream& os, const std::string& indent = "") const; - + virtual bool operator==(const DICOMDatasetSorter& other) const; protected: /** \brief Helper struct to feed into std::sort, configured via DICOMSortCriterion. */ struct ParameterizedDatasetSort { ParameterizedDatasetSort(DICOMSortCriterion::ConstPointer); bool operator() (const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right); DICOMSortCriterion::ConstPointer m_SortCriterion; }; DICOMTagBasedSorter(); virtual ~DICOMTagBasedSorter(); DICOMTagBasedSorter(const DICOMTagBasedSorter& other); DICOMTagBasedSorter& operator=(const DICOMTagBasedSorter& other); /** \brief Helper for SplitInputGroups(). */ std::string BuildGroupID( DICOMDatasetAccess* dataset ); typedef std::map GroupIDToListType; /** \brief Implements the "distiguishing tags". To sort datasets into different groups, a long string will be built for each dataset. The string concatenates all tags and their respective values. Datasets that match in all values will end up with the same string. */ GroupIDToListType SplitInputGroups(); /** \brief Implements the sorting step. Relatively simple implementation thanks to std::sort and a parameterization via DICOMSortCriterion. */ GroupIDToListType& SortGroups(GroupIDToListType& groups); DICOMTagList m_DistinguishingTags; typedef std::map TagValueProcessorMap; TagValueProcessorMap m_TagValueProcessor; DICOMSortCriterion::ConstPointer m_SortCriterion; + + bool m_StrictSorting; }; } #endif diff --git a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp index 5c7846997b..b65efa5ff0 100644 --- a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp @@ -1,572 +1,599 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "mitkEquiDistantBlocksSorter.h" mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::SliceGroupingAnalysisResult() { } mitk::DICOMDatasetList mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::GetBlockDatasets() { return m_GroupedFiles; } mitk::DICOMDatasetList mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult ::GetUnsortedDatasets() { return m_UnsortedFiles; } 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 EquiDistantBlocksSorter* 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) ); // GantryTilt return tags; } void mitk::EquiDistantBlocksSorter ::Sort() { DICOMDatasetList remainingInput = GetInput(); // copy typedef std::list OutputListType; OutputListType outputs; m_SliceGroupingResults.clear(); while (!remainingInput.empty()) // repeat until all files are grouped somehow { SliceGroupingAnalysisResult regularBlock = this->AnalyzeFileForITKImageSeriesReaderSpacingAssumption( remainingInput, m_AcceptTilt ); DICOMDatasetList inBlock = regularBlock.GetBlockDatasets(); DICOMDatasetList laterBlock = regularBlock.GetUnsortedDatasets(); MITK_DEBUG << "Result: sorted 3D group with " << inBlock.size() << " files"; for (DICOMDatasetList::const_iterator diter = inBlock.begin(); diter != inBlock.end(); ++diter) MITK_DEBUG << " IN " << (*diter)->GetFilenameIfAvailable(); for (DICOMDatasetList::const_iterator diter = laterBlock.begin(); diter != laterBlock.end(); ++diter) MITK_DEBUG << " OUT " << (*diter)->GetFilenameIfAvailable(); outputs.push_back( regularBlock.GetBlockDatasets() ); m_SliceGroupingResults.push_back( regularBlock ); remainingInput = regularBlock.GetUnsortedDatasets(); } unsigned int numberOfOutputs = outputs.size(); this->SetNumberOfOutputs(numberOfOutputs); unsigned int outputIndex(0); for (OutputListType::iterator oIter = outputs.begin(); oIter != outputs.end(); ++outputIndex, ++oIter) { this->SetOutput(outputIndex, *oIter); } } 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 unprecise locations!"; } } void mitk::EquiDistantBlocksSorter ::SetToleratedOriginOffset(double millimeters) { m_ToleratedOriginOffset = millimeters; m_ToleratedOriginOffsetIsAbsolute = true; if (m_ToleratedOriginOffset < 0.0) { MITK_WARN << "Negative tolerance set to SetToleratedOriginOffset()!"; } } double mitk::EquiDistantBlocksSorter ::GetToleratedOriginOffset() const { return m_ToleratedOriginOffset; } bool mitk::EquiDistantBlocksSorter ::IsToleratedOriginOffsetAbsolute() const { return m_ToleratedOriginOffsetIsAbsolute; } std::string mitk::EquiDistantBlocksSorter ::ConstCharStarToString(const char* s) { return s ? std::string(s) : std::string(); } mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult mitk::EquiDistantBlocksSorter ::AnalyzeFileForITKImageSeriesReaderSpacingAssumption( const DICOMDatasetList& datasets, bool groupImagesWithGantryTilt) { // result.first = files that fit ITK's assumption // result.second = files that do not fit, should be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption() again SliceGroupingAnalysisResult result; // we const_cast here, because I could not use a map.at(), which would make the code much more readable const DICOMTag tagImagePositionPatient = DICOMTag(0x0020,0x0032); // Image Position (Patient) const DICOMTag tagImageOrientation = DICOMTag(0x0020, 0x0037); // Image Orientation const DICOMTag tagGantryTilt = DICOMTag(0x0018, 0x1120); // gantry tilt 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 (DICOMDatasetList::const_iterator dsIter = datasets.begin(); dsIter != datasets.end(); ++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 ); 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 ); 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 ); 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 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 ); 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() ) { // check if this is at least roughly the same angle as recorded in DICOM tags double angle = 0.0; std::string tiltStr = (*dsIter)->GetTagValueAsString( tagGantryTilt ); std::istringstream i(tiltStr); if (i >> angle) { MITK_DEBUG << "Comparing recorded tilt angle " << angle << " against calculated value " << tiltInfo.GetTiltAngleInDegrees(); // TODO we probably want the signs correct, too (that depends: this is just a rough check, nothing serious) if ( fabs(angle) - tiltInfo.GetTiltAngleInDegrees() > 0.25) { result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis fileFitsIntoPattern = false; } else // tilt angle from header is less than 0.25 degrees different from what we calculated, assume this is fine { 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 // we cannot check the calculated tilt angle against the one from the dicom header (so we assume we are right) { 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 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 fileFitsIntoPattern = false; } else { result.AddFileToSortedBlock( *dsIter ); // this file is good for current block fileFitsIntoPattern = true; } } 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 - if ( result.GetBlockDatasets().size() == 2 ) + // Above behavior can be configured via m_AcceptTwoSlicesGroups, the default being "do accept" + if ( result.GetBlockDatasets().size() == 2 && !m_AcceptTwoSlicesGroups ) { result.UndoPrematureGrouping(); } } // 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 ); std::string firstOriginString = firstDataset->GetTagValueAsString( tagImagePositionPatient ); std::string lastOriginString = lastDataset->GetTagValueAsString( tagImagePositionPatient ); result.FlagGantryTilt( GantryTiltInformation::MakeFromTagValues( firstOriginString, lastOriginString, orientationString, numberOfSlicesApart )); } catch (...) { // just do not flag anything, we are ok } } return result; } diff --git a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h index 92a4790ae0..9de12893ee 100644 --- a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h @@ -1,210 +1,216 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkEquiDistantBlocksSorter_h #define mitkEquiDistantBlocksSorter_h #include "mitkDICOMDatasetSorter.h" #include "mitkDICOMSortCriterion.h" #include "mitkGantryTiltInformation.h" #include "mitkVector.h" #include namespace mitk { /** \ingroup DICOMReaderModule \brief Split inputs into blocks of equidistant slices (for use in DICOMITKSeriesGDCMReader). Since inter-slice distance is not recorded in DICOM tags, we must ensure that blocks are made up of slices that have equal distances between neighboring slices. This is especially necessary because itk::ImageSeriesReader is later used for the actual loading, and this class expects (and does nocht verify) equal inter-slice distance (see \ref DICOMITKSeriesGDCMReader_ForcedConfiguration). To achieve such grouping, the inter-slice distance is calculated from the first two different slice positions of a block. Following slices are added to a block as long as they can be added by adding the calculated inter-slice distance to the last slice of the block. Slices that do not fit into the expected distance pattern, are set aside for further analysis. This grouping is done until each file has been assigned to a group. Slices that share a position in space are also sorted into separate blocks during this step. So the result of this step is a set of blocks that contain only slices with equal z spacing and uniqe slices at each position. During sorting, the origins (documented in tag image position patient) are compared against expected origins (from former origin plus moving direction). As there will be minor differences in numbers (from both calculations and unprecise tag values), we must be a bit tolerant here. The default behavior is to expect that an origin is not further away from the expected position than 30% of the inter-slice distance. To support a legacy behavior of a former loader (DicomSeriesReader), this default can be restricted to a constant number of millimeters by calling SetToleratedOriginOffset(mm). Detailed implementation in AnalyzeFileForITKImageSeriesReaderSpacingAssumption(). */ class DICOMReader_EXPORT EquiDistantBlocksSorter : public DICOMDatasetSorter { public: mitkClassMacro( EquiDistantBlocksSorter, DICOMDatasetSorter ) itkNewMacro( EquiDistantBlocksSorter ) virtual DICOMTagList GetTagsOfInterest(); /** \brief Delegates work to AnalyzeFileForITKImageSeriesReaderSpacingAssumption(). AnalyzeFileForITKImageSeriesReaderSpacingAssumption() is called until it does not create multiple blocks anymore. */ virtual void Sort(); /** \brief Whether or not to accept images from a tilted acquisition in a single output group. */ void SetAcceptTilt(bool accept); + bool GetAcceptTilt() const; /** \brief See class description and SetToleratedOriginOffset(). */ void SetToleratedOriginOffsetToAdaptive(double fractionOfInterSliceDistanct = 0.3); /** \brief See class description and SetToleratedOriginOffsetToAdaptive(). Default value of 0.005 is calculated so that we get a maximum of 1/10mm error when having a measurement crosses 20 slices in z direction (too strict? we don't know better..). */ void SetToleratedOriginOffset(double millimeters = 0.005); double GetToleratedOriginOffset() const; bool IsToleratedOriginOffsetAbsolute() const; + void SetAcceptTwoSlicesGroups(bool accept); + bool GetAcceptTwoSlicesGroups() const; + virtual void PrintConfiguration(std::ostream& os, const std::string& indent = "") const; virtual bool operator==(const DICOMDatasetSorter& other) const; protected: /** \brief Return type of AnalyzeFileForITKImageSeriesReaderSpacingAssumption(). Class contains the grouping result of method AnalyzeFileForITKImageSeriesReaderSpacingAssumption(), which takes as input a number of images, which are all equally oriented and spatially sorted along their normal direction. The result contains of two blocks: a first one is the grouping result, all of those images can be loaded into one image block because they have an equal origin-to-origin distance without any gaps in-between. */ class SliceGroupingAnalysisResult { public: SliceGroupingAnalysisResult(); /** \brief Grouping result, all same origin-to-origin distance w/o gaps. */ DICOMDatasetList GetBlockDatasets(); void SetFirstFilenameOfBlock(const std::string& filename); std::string GetFirstFilenameOfBlock() const; void SetLastFilenameOfBlock(const std::string& filename); std::string GetLastFilenameOfBlock() const; /** \brief Remaining files, which could not be grouped. */ DICOMDatasetList GetUnsortedDatasets(); /** \brief Wheter or not the grouped result contain a gantry tilt. */ bool ContainsGantryTilt(); /** \brief Detailed description of gantry tilt. */ const GantryTiltInformation& GetTiltInfo() const; /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. */ void AddFileToSortedBlock(DICOMDatasetAccess* dataset); /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. */ void AddFileToUnsortedBlock(DICOMDatasetAccess* dataset); void AddFilesToUnsortedBlock(const DICOMDatasetList& datasets); /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. \todo Could make sense to enhance this with an instance of GantryTiltInformation to store the whole result! */ void FlagGantryTilt(const GantryTiltInformation& tiltInfo); /** \brief Only meaningful for use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption. */ void UndoPrematureGrouping(); protected: DICOMDatasetList m_GroupedFiles; DICOMDatasetList m_UnsortedFiles; GantryTiltInformation m_TiltInfo; std::string m_FirstFilenameOfBlock; std::string m_LastFilenameOfBlock; }; /** \brief Ensure an equal z-spacing for a group of files. Takes as input a number of images, which are all equally oriented and spatially sorted along their normal direction. Internally used by GetSeries. Returns two lists: the first one contins slices of equal inter-slice spacing. The second list contains remaining files, which need to be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption again. Relevant code that is matched here is in itkImageSeriesReader.txx (ImageSeriesReader::GenerateOutputInformation(void)), lines 176 to 245 (as of ITK 3.20) */ SliceGroupingAnalysisResult AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const DICOMDatasetList& files, bool groupsOfSimilarImages); /** \brief Safely convert const char* to std::string. */ std::string ConstCharStarToString(const char* s); EquiDistantBlocksSorter(); virtual ~EquiDistantBlocksSorter(); EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other); EquiDistantBlocksSorter& operator=(const EquiDistantBlocksSorter& other); bool m_AcceptTilt; typedef std::vector ResultsList; ResultsList m_SliceGroupingResults; double m_ToleratedOriginOffset; bool m_ToleratedOriginOffsetIsAbsolute; + + bool m_AcceptTwoSlicesGroups; }; } #endif diff --git a/Modules/DICOMReader/mitkSortByImagePositionPatient.cpp b/Modules/DICOMReader/mitkSortByImagePositionPatient.cpp index da30f6be08..a3b141811e 100644 --- a/Modules/DICOMReader/mitkSortByImagePositionPatient.cpp +++ b/Modules/DICOMReader/mitkSortByImagePositionPatient.cpp @@ -1,142 +1,170 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkSortByImagePositionPatient.h" #include "mitkDICOMTag.h" mitk::SortByImagePositionPatient ::SortByImagePositionPatient(DICOMSortCriterion::Pointer secondaryCriterion) :DICOMSortCriterion(secondaryCriterion) { } mitk::SortByImagePositionPatient ::~SortByImagePositionPatient() { } mitk::SortByImagePositionPatient ::SortByImagePositionPatient(const SortByImagePositionPatient& other ) :DICOMSortCriterion(other) { } mitk::SortByImagePositionPatient& mitk::SortByImagePositionPatient ::operator=(const SortByImagePositionPatient& other) { if (this != &other) { DICOMSortCriterion::operator=(other); } return *this; } - + bool mitk::SortByImagePositionPatient ::operator==(const DICOMSortCriterion& other) const { return dynamic_cast(&other) != NULL; // same class } void mitk::SortByImagePositionPatient ::Print(std::ostream& os) const { os << "(0020,0032) Image Position (Patient) along normal of (0020,0037) Image Orientation (Patient)"; } mitk::DICOMTagList mitk::SortByImagePositionPatient ::GetTagsOfInterest() const { DICOMTagList tags; tags.push_back( DICOMTag(0x0020, 0x0032) ); // ImagePositionPatient tags.push_back( DICOMTag(0x0020, 0x0037) ); // ImageOrientationPatient return tags; } bool mitk::SortByImagePositionPatient ::IsLeftBeforeRight(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right) const +{ + bool possible(false); + double distance = InternalNumericDistance(left, right, possible); // returns 0.0 if not possible + if (possible) + { + return distance > 0.0; + } + else + { + return this->NextLevelIsLeftBeforeRight(left, right); + } +} + +double +mitk::SortByImagePositionPatient +::InternalNumericDistance(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right, bool& possible) const { // sort by distance to world origin, assuming (almost) equal orientation static const DICOMTag tagImagePositionPatient = DICOMTag(0x0020,0x0032); // Image Position (Patient) static const DICOMTag tagImageOrientation = DICOMTag(0x0020, 0x0037); // Image Orientation static Vector3D leftRight; leftRight.Fill(0.0); static Vector3D leftUp; leftUp.Fill(0.0); static bool leftHasOrientation(false); DICOMStringToOrientationVectors( left->GetTagValueAsString( tagImageOrientation ), leftRight, leftUp, leftHasOrientation ); static Vector3D rightRight; rightRight.Fill(0.0); static Vector3D rightUp; rightUp.Fill(0.0); static bool rightHasOrientation(false); DICOMStringToOrientationVectors( right->GetTagValueAsString( tagImageOrientation ), rightRight, rightUp, rightHasOrientation ); static Point3D leftOrigin; leftOrigin.Fill(0.0f); static bool leftHasOrigin(false); leftOrigin = DICOMStringToPoint3D( left->GetTagValueAsString( tagImagePositionPatient ), leftHasOrigin ); static Point3D rightOrigin; rightOrigin.Fill(0.0f); static bool rightHasOrigin(false); rightOrigin = DICOMStringToPoint3D( right->GetTagValueAsString( tagImagePositionPatient ), rightHasOrigin ); // we tolerate very small differences in image orientation, since we got to know about // acquisitions where these values change across a single series (7th decimal digit) // (http://bugs.mitk.org/show_bug.cgi?id=12263) // still, we want to check if our assumption of 'almost equal' orientations is valid for (unsigned int dim = 0; dim < 3; ++dim) { if ( fabs(leftRight[dim] - rightRight[dim]) > 0.0001 || fabs(leftUp[dim] - rightUp[dim]) > 0.0001) { MITK_ERROR << "Dicom images have different orientations."; throw std::logic_error("Dicom images have different orientations. Call GetSeries() first to separate images."); } } static Vector3D normal; normal[0] = leftRight[1] * leftUp[2] - leftRight[2] * leftUp[1]; normal[1] = leftRight[2] * leftUp[0] - leftRight[0] * leftUp[2]; normal[2] = leftRight[0] * leftUp[1] - leftRight[1] * leftUp[0]; static double leftDistance = 0.0; static double rightDistance = 0.0; leftDistance = 0.0; rightDistance = 0.0; // this computes the distance from world origin (0,0,0) ALONG THE NORMAL of the image planes for (unsigned int dim = 0; dim < 3; ++dim) { leftDistance += normal[dim] * leftOrigin[dim]; rightDistance += normal[dim] * rightOrigin[dim]; } // if we can sort by just comparing the distance, we do exactly that if ( fabs(leftDistance - rightDistance) >= mitk::eps) { + possible = true; // default: compare position - return leftDistance < rightDistance; + return rightDistance - leftDistance; // if (left < right> ==> diff > 0 } else { - return this->NextLevelIsLeftBeforeRight(left, right); + possible = false; + return 0.0; } } + + +double +mitk::SortByImagePositionPatient +::NumericDistance(const mitk::DICOMDatasetAccess* from, const mitk::DICOMDatasetAccess* to) const +{ + bool possible(false); + double retVal = InternalNumericDistance(from, to, possible); // returns 0.0 if not possible + return possible ? retVal : 0.0; +} diff --git a/Modules/DICOMReader/mitkSortByImagePositionPatient.h b/Modules/DICOMReader/mitkSortByImagePositionPatient.h index bc9fdad546..1542edaa05 100644 --- a/Modules/DICOMReader/mitkSortByImagePositionPatient.h +++ b/Modules/DICOMReader/mitkSortByImagePositionPatient.h @@ -1,64 +1,68 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkSortByImagePositionPatient_h #define mitkSortByImagePositionPatient_h #include "mitkDICOMSortCriterion.h" #include "mitkVector.h" namespace mitk { /** \ingroup DICOMReaderModule \brief Sort by distance of image origin along image normal (for use in DICOMTagBasedSorter). To compare two datasets, their distance to the world origin is calculated. This distance is calculated along the image normals because we do not know the image orientation in advance, to any of the three coordinates could be identical for all datasets. \note This class assumes that the datasets have identical orientations! */ class DICOMReader_EXPORT SortByImagePositionPatient : public DICOMSortCriterion { public: mitkClassMacro( SortByImagePositionPatient, DICOMSortCriterion ); mitkNewMacro1Param( SortByImagePositionPatient, DICOMSortCriterion::Pointer ); virtual DICOMTagList GetTagsOfInterest() const; virtual bool IsLeftBeforeRight(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right) const; + virtual double NumericDistance(const mitk::DICOMDatasetAccess* from, const mitk::DICOMDatasetAccess* to) const; + virtual void Print(std::ostream& os) const; - + virtual bool operator==(const DICOMSortCriterion& other) const; protected: SortByImagePositionPatient( DICOMSortCriterion::Pointer secondaryCriterion = NULL ); virtual ~SortByImagePositionPatient(); SortByImagePositionPatient(const SortByImagePositionPatient& other); SortByImagePositionPatient& operator=(const SortByImagePositionPatient& other); + double InternalNumericDistance(const mitk::DICOMDatasetAccess* from, const mitk::DICOMDatasetAccess* to, bool& possible) const; + private: }; } #endif