diff --git a/Modules/DICOMReader/TODOs.txt b/Modules/DICOMReader/TODOs.txt index ccb8d3696d..b62e34cfeb 100644 --- a/Modules/DICOMReader/TODOs.txt +++ b/Modules/DICOMReader/TODOs.txt @@ -1,60 +1,64 @@ Important [x] Make sure that ..Selector does tag scanning only once -[ ] Check Gantry Tilt TODO/Bug in old reader -[ ] - and in new reader, how do we WANT to handle this? +[x] Check Gantry Tilt TODO/Bug in old reader +[x] - and in new reader, how do we WANT to handle this? + - old and new code both did not correctly check the stored value against the calculated value + - it SEEMS like tilts below 45 degrees were never checked against the tag, and bigger values are probably unrealistic(?) + - so, nothing wrong. It just seems that we calculated a DIFFERENT angle than the one recorded + - in the new reader we remove the check. We are right as far as our shearing is affected. [x] Keep OLD DICOMTesting module, in order to test code until it is removed - restored this module as DCMTesting module [x] ONLY load a pre-sorted list - Should work now by: - get the list of images, construct FrameInfo list (order must be known from prior sorting) - prepare a tagcache object that provides tag information (from some kind of database, or run AnalyzeInputFiles() on file reader) - create DICOMImageBlockDescriptor with frame list and tag cache - Call SetFixTiltByShearing() and TiltInfo object as appropriate... - call DICOMITKSeriesGDCMReader::LoadMitkImageForImageBlockDescriptor(block) [x] Option to just provide properties to a pre-loaded mitk::Image - see TestDICOMLoading::DecorateVerifyCachedImage() [x] Performance of mitkdump / DICOMReaderSelector (TagCache?) - the note in GDCM..Reader::GetTagValue() is correct. This piece of code is taking lots of time. Approx. half of the time is tag scanning, the other half is construction of block describing properties, which accesses GetTagValue. [x] - maybe we can evaluate these properties in a lazy way (only when asked for). Solution: Yes, this helps, implemented. [x] Gantry tilt broken during implementation - Was a wrong implementation of NormalDirectionConsistencySorter [x] Accepted image origin error during EquiDistantBlocksSorter must be configurable. To support legacy behavior of the reader, the tolerance must be fixable to a constant value. For newer readers, it should be adapted to the actual slice distance. [x] Sorting by "Image Time" seems undeterministic (clarkson test data) [x] - Numeric conversion of "1.2.3.4" yields 1.2 on Windows, seems to fail on Linux(?) - need to verify that the whole string (except whitespace at the begin/end) was converted Solution: GetTagValue as String does a right-trim plus we check after conversion! Works. [x] Implement Configurator::GetBuiltIn3DReaders() (don't ignore "classic reader") [x] Complete MITK properties of images: level/window and color type (MONOCHROME1/2) [x] Check input images fo DICOMness and non-multi-frameness [x] Use CanHandleFile() in selection! - implicitly now, the reader checks itself and produces no output if it cannot handle the files [x] Check ToDo in mitkDICOMTag.cpp - operator< for DICOMTag has been re-implemented in a readable and safer way. [x] Configurable precision for tag value comparisons [x] Images are upside-down in some cases - error was hidden assumption somewhere: filenames for ImageSeriesReader need to be in an order that goes along the image normals, not in the opposite direction Questionable [ ] Fallback to filename instead of instance UID? Nice-to-have [ ] Multi-Frame images (DCMTK) [ ] ... -[ ] Add the tags used by DICOMImageBlockDescriptor dynamically instead of hard-coded +[x] Add the tags used by DICOMImageBlockDescriptor dynamically instead of hard-coded (look for ToDo around m_GDCMScanner.AddTag() in mitkDICOMITKSeriesGDCMReader.cpp) [ ] Let DICOMImageBlockDescriptor describe attributes common and different for all blocks of a series [ ] Let DICOMImageBlockDescriptor check image properties when mitk::Image is set (pre-loaded case) diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp index 84738d9d5e..b3032c0aaf 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp @@ -1,733 +1,721 @@ /*=================================================================== 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 //#define ENABLE_TIMING #include "mitkDICOMITKSeriesGDCMReader.h" #include "mitkITKDICOMSeriesReaderHelper.h" #include "mitkGantryTiltInformation.h" #include "mitkDICOMTagBasedSorter.h" #include "mitkDICOMGDCMTagScanner.h" #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) ,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 } #if defined(MBILOG_ENABLE_DEBUG) || defined (ENABLE_TIMING) #define timeStart(part) timer.Start(part); #define timeStop(part) timer.Stop(part); #else #define timeStart(part) #define timeStop(part) #endif void mitk::DICOMITKSeriesGDCMReader ::AnalyzeInputFiles() { itk::TimeProbesCollectorBase timer; timeStart("Reset"); this->ClearOutputs(); timeStop("Reset"); // prepare initial sorting (== list of input files) StringList inputFilenames = this->GetInputFiles(); timeStart("Check input for DCM"); if ( inputFilenames.empty() || !this->CanHandleFile( inputFilenames.front() ) // first || !this->CanHandleFile( inputFilenames.back() ) // last || !this->CanHandleFile( inputFilenames[ inputFilenames.size() / 2] ) // roughly central file ) { // TODO a read-as-many-as-possible fallback could be implemented here MITK_DEBUG << "Reader unable to process files.."; return; } timeStop("Check input for DCM"); // scan files for sorting-relevant tags if (m_TagCache.IsNull()) { timeStart("Tag scanning"); DICOMGDCMTagScanner::Pointer filescanner = DICOMGDCMTagScanner::New(); m_TagCache = filescanner.GetPointer(); // keep alive and make accessible to sub-classes filescanner->SetInputFiles(inputFilenames); filescanner->AddTags( this->GetTagsOfInterest() ); PushLocale(); filescanner->Scan(); PopLocale(); timeStop("Tag scanning"); } else { // ensure that the tag cache contains our required tags AND files and has scanned! } m_SortingResultInProgress.clear(); // TODO We should remove the following cast // DICOMImageFrameInfo would need to inherit DICOMDatasetAccess! // - then the DICOMGDCMTagScanner class could create a DICOMGDCMImageFrameList internally // - and return it as a DICOMImageFrameList // - like this, DICOMITKSeriesGDCMReader could use the DICOMImageFrameInfoList to feed its sorters // - problem: // - DICOMImageFrameInfo is also part of DICOMImageBlockDescriptor, which is meant // to describe the scanner output, even after the reader (and its scanner) is deleted. // - if DICOMImageFrameInfo now inherits DICOMDatasetAccess, it would also need to implement // GetTagValueAsString(). // - so this could all work if we implement a default response in DICOMImageFrameInfo::GetTagValueAsString() (like in GetFilenameIfAvailable) // and overwrite it in DICOMGDCMImageFrameInfo, which also knows about a specific GDCM scanner result // (which again COULD (no need to?) be hidden as a point to a DICOMGDCMTagScanner class) // if ( DICOMGDCMTagScanner* tagCache = dynamic_cast(m_TagCache.GetPointer()) ) { m_SortingResultInProgress.push_back( tagCache->GetFrameInfoList() ); } else { throw std::logic_error("Bad implementation error: DICOMITKSeriesGDCMReader now unable to find dataset/tag information for its input."); } // sort and split blocks as configured timeStart("Sorting frames"); unsigned int sorterIndex = 0; for(SorterList::iterator sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sorterIndex, ++sorterIter) { std::stringstream ss; ss << "Sorting step " << sorterIndex; timeStart( ss.str().c_str() ); m_SortingResultInProgress = this->InternalExecuteSortingStep(sorterIndex, *sorterIter, m_SortingResultInProgress); timeStop( ss.str().c_str() ); } // a last extra-sorting step: ensure equidistant slices timeStart( "EquiDistantBlocksSorter" ); m_SortingResultInProgress = this->InternalExecuteSortingStep(sorterIndex++, m_EquiDistantBlocksSorter.GetPointer(), m_SortingResultInProgress); timeStop( "EquiDistantBlocksSorter" ); timeStop("Sorting frames"); timeStart("Condensing 3D blocks"); m_SortingResultInProgress = this->Condense3DBlocks( m_SortingResultInProgress ); timeStop("Condensing 3D blocks"); // provide final result as output timeStart("Output"); unsigned int o = this->GetNumberOfOutputs(); this->SetNumberOfOutputs( o + m_SortingResultInProgress.size() ); // Condense3DBlocks may already have added outputs! for (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->GetTagCache() ); // 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 ); } timeStop("Output"); #if defined(MBILOG_ENABLE_DEBUG) || defined (ENABLE_TIMING) std::cout << "---------------------------------------------------------------" << std::endl; timer.Report( std::cout ); std::cout << "---------------------------------------------------------------" << std::endl; #endif } mitk::DICOMITKSeriesGDCMReader::SortingBlockList mitk::DICOMITKSeriesGDCMReader ::InternalExecuteSortingStep( unsigned int sortingStepIndex, DICOMDatasetSorter::Pointer sorter, const SortingBlockList& input) { SortingBlockList nextStepSorting; // we should not modify our input list while processing it std::stringstream ss; ss << "Sorting step " << sortingStepIndex << " '"; #if defined(MBILOG_ENABLE_DEBUG) sorter->PrintConfiguration(ss); #endif ss << "'"; nextStepSorting.clear(); MITK_DEBUG << "================================================================================"; MITK_DEBUG << "DICOMITKSeriesGDCMReader: " << ss.str() << ": " << input.size() << " groups input"; 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 ); } } 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; } mitk::DICOMTagCache::Pointer mitk::DICOMITKSeriesGDCMReader ::GetTagCache() const { return m_TagCache; } void mitk::DICOMITKSeriesGDCMReader ::SetTagCache(DICOMTagCache::Pointer tagCache) { m_TagCache = tagCache; } mitk::DICOMTagList mitk::DICOMITKSeriesGDCMReader ::GetTagsOfInterest() const { DICOMTagList completeList; + // check all configured sorters for(SorterList::const_iterator sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sorterIter) { assert(sorterIter->IsNotNull()); DICOMTagList tags = (*sorterIter)->GetTagsOfInterest(); completeList.insert( completeList.end(), tags.begin(), tags.end() ); } - // 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) - completeList.push_back( DICOMTag(0x0018,0x1164) ); // pixel spacing - completeList.push_back( DICOMTag(0x0028,0x0030) ); // imager pixel spacing - - completeList.push_back( DICOMTag(0x0028,0x1050) ); // window center - completeList.push_back( DICOMTag(0x0028,0x1051) ); // window width - completeList.push_back( DICOMTag(0x0008,0x0008) ); // image type - completeList.push_back( DICOMTag(0x0028,0x0004) ); // photometric interpretation - - completeList.push_back( DICOMTag(0x0020,0x1041) ); // slice location - completeList.push_back( DICOMTag(0x0020,0x0013) ); // instance number - completeList.push_back( DICOMTag(0x0008,0x0016) ); // sop class UID - completeList.push_back( DICOMTag(0x0008,0x0018) ); // sop instance UID - - completeList.push_back( DICOMTag(0x0020,0x0011) ); // series number - completeList.push_back( DICOMTag(0x0008,0x1030) ); // study description - completeList.push_back( DICOMTag(0x0008,0x103e) ); // series description - completeList.push_back( DICOMTag(0x0008,0x0060) ); // modality - completeList.push_back( DICOMTag(0x0020,0x0012) ); // acquisition number - completeList.push_back( DICOMTag(0x0018,0x0024) ); // sequence name - completeList.push_back( DICOMTag(0x0020,0x0037) ); // image orientation - completeList.push_back( DICOMTag(0x0020,0x0032) ); // ipp + // check our own forced sorters + DICOMTagList tags = m_EquiDistantBlocksSorter->GetTagsOfInterest(); + completeList.insert( completeList.end(), tags.begin(), tags.end() ); + + tags = m_NormalDirectionConsistencySorter->GetTagsOfInterest(); + completeList.insert( completeList.end(), tags.begin(), tags.end() ); + + // add the tags for DICOMImageBlockDescriptor + tags = DICOMImageBlockDescriptor::GetTagsOfInterest(); + completeList.insert( completeList.end(), tags.begin(), tags.end() ); return completeList; } diff --git a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp index 7fddbfc90c..a9e145e0b6 100644 --- a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp +++ b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp @@ -1,711 +1,742 @@ /*=================================================================== 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 "mitkDICOMImageBlockDescriptor.h" #include "mitkStringProperty.h" #include "mitkLevelWindowProperty.h" #include mitk::DICOMImageBlockDescriptor ::DICOMImageBlockDescriptor() :m_ReaderImplementationLevel(SOPClassUnknown) ,m_PropertyList(PropertyList::New()) ,m_TagCache(NULL) ,m_PropertiesOutOfDate(true) { } mitk::DICOMImageBlockDescriptor ::~DICOMImageBlockDescriptor() { } mitk::DICOMImageBlockDescriptor ::DICOMImageBlockDescriptor(const DICOMImageBlockDescriptor& other) :m_ImageFrameList( other.m_ImageFrameList ) ,m_MitkImage( other.m_MitkImage ) ,m_SliceIsLoaded( other.m_SliceIsLoaded ) ,m_ReaderImplementationLevel( other.m_ReaderImplementationLevel ) ,m_TiltInformation( other.m_TiltInformation ) ,m_PropertyList( other.m_PropertyList->Clone() ) ,m_TagCache( other.m_TagCache ) ,m_PropertiesOutOfDate( other.m_PropertiesOutOfDate ) { if (m_MitkImage) { m_MitkImage = m_MitkImage->Clone(); } } mitk::DICOMImageBlockDescriptor& mitk::DICOMImageBlockDescriptor ::operator=(const DICOMImageBlockDescriptor& other) { if (this != &other) { m_ImageFrameList = other.m_ImageFrameList; m_MitkImage = other.m_MitkImage; m_SliceIsLoaded = other.m_SliceIsLoaded; m_ReaderImplementationLevel = other.m_ReaderImplementationLevel; m_TiltInformation = other.m_TiltInformation; if (other.m_PropertyList) { m_PropertyList = other.m_PropertyList->Clone(); } if (other.m_MitkImage) { m_MitkImage = other.m_MitkImage->Clone(); } m_TagCache = other.m_TagCache; m_PropertiesOutOfDate = other.m_PropertiesOutOfDate; } return *this; } +mitk::DICOMTagList +mitk::DICOMImageBlockDescriptor::GetTagsOfInterest() +{ + DICOMTagList completeList; + + completeList.push_back( DICOMTag(0x0018,0x1164) ); // pixel spacing + completeList.push_back( DICOMTag(0x0028,0x0030) ); // imager pixel spacing + + completeList.push_back( DICOMTag(0x0008,0x0018) ); // sop instance UID + completeList.push_back( DICOMTag(0x0008,0x0016) ); // sop class UID + + completeList.push_back( DICOMTag(0x0020,0x0011) ); // series number + completeList.push_back( DICOMTag(0x0008,0x1030) ); // study description + completeList.push_back( DICOMTag(0x0008,0x103e) ); // series description + completeList.push_back( DICOMTag(0x0008,0x0060) ); // modality + completeList.push_back( DICOMTag(0x0018,0x0024) ); // sequence name + completeList.push_back( DICOMTag(0x0020,0x0037) ); // image orientation + + completeList.push_back( DICOMTag(0x0020,0x1041) ); // slice location + completeList.push_back( DICOMTag(0x0020,0x0012) ); // acquisition number + completeList.push_back( DICOMTag(0x0020,0x0013) ); // instance number + completeList.push_back( DICOMTag(0x0020,0x0032) ); // image position patient + + completeList.push_back( DICOMTag(0x0028,0x1050) ); // window center + completeList.push_back( DICOMTag(0x0028,0x1051) ); // window width + completeList.push_back( DICOMTag(0x0008,0x0008) ); // image type + completeList.push_back( DICOMTag(0x0028,0x0004) ); // photometric interpretation + + return completeList; +} + void mitk::DICOMImageBlockDescriptor ::SetTiltInformation(const GantryTiltInformation& info) { m_TiltInformation = info; } const mitk::GantryTiltInformation mitk::DICOMImageBlockDescriptor ::GetTiltInformation() const { return m_TiltInformation; } void mitk::DICOMImageBlockDescriptor ::SetImageFrameList(const DICOMImageFrameList& framelist) { m_ImageFrameList = framelist; m_SliceIsLoaded.resize(framelist.size()); m_SliceIsLoaded.assign(framelist.size(), false); m_PropertiesOutOfDate = true; } const mitk::DICOMImageFrameList& mitk::DICOMImageBlockDescriptor ::GetImageFrameList() const { return m_ImageFrameList; } void mitk::DICOMImageBlockDescriptor ::SetMitkImage(Image::Pointer image) { if (m_MitkImage != image) { if (m_TagCache.IsNull()) { MITK_ERROR << "Unable to describe MITK image with properties without a tag-cache object!"; m_MitkImage = NULL; return; } if (m_ImageFrameList.empty()) { MITK_ERROR << "Unable to describe MITK image with properties without a frame list!"; m_MitkImage = NULL; return; } // Should verify that the image matches m_ImageFrameList and m_TagCache // however, this is hard to do without re-analyzing all // TODO we should at least make sure that the number of frames is identical (plus rows/columns, orientation) // without gantry tilt correction, we can also check image origin m_MitkImage = this->DescribeImageWithProperties( this->FixupSpacing(image) ); } } mitk::Image::Pointer mitk::DICOMImageBlockDescriptor ::GetMitkImage() const { return m_MitkImage; } mitk::Image::Pointer mitk::DICOMImageBlockDescriptor ::FixupSpacing(Image* mitkImage) { if (mitkImage) { Vector3D imageSpacing = mitkImage->GetGeometry()->GetSpacing(); ScalarType desiredSpacingX = imageSpacing[0]; ScalarType desiredSpacingY = imageSpacing[1]; this->GetDesiredMITKImagePixelSpacing( desiredSpacingX, desiredSpacingY ); // prefer pixel spacing over imager pixel spacing MITK_DEBUG << "Loaded image with spacing " << imageSpacing[0] << ", " << imageSpacing[1]; MITK_DEBUG << "Found correct spacing info " << desiredSpacingX << ", " << desiredSpacingY; imageSpacing[0] = desiredSpacingX; imageSpacing[1] = desiredSpacingY; mitkImage->GetGeometry()->SetSpacing( imageSpacing ); } return mitkImage; } void mitk::DICOMImageBlockDescriptor ::SetSliceIsLoaded(unsigned int index, bool isLoaded) { if (index < m_SliceIsLoaded.size()) { m_SliceIsLoaded[index] = isLoaded; } else { std::stringstream ss; ss << "Index " << index << " out of range (" << m_SliceIsLoaded.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } bool mitk::DICOMImageBlockDescriptor ::IsSliceLoaded(unsigned int index) const { if (index < m_SliceIsLoaded.size()) { return m_SliceIsLoaded[index]; } else { std::stringstream ss; ss << "Index " << index << " out of range (" << m_SliceIsLoaded.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } bool mitk::DICOMImageBlockDescriptor ::AllSlicesAreLoaded() const { bool allLoaded = true; for (BoolList::const_iterator iter = m_SliceIsLoaded.begin(); iter != m_SliceIsLoaded.end(); ++iter) { allLoaded &= *iter; } return allLoaded; } /* PS defined IPS defined PS==IPS 0 0 --> UNKNOWN spacing, loader will invent 0 1 --> spacing as at detector surface 1 0 --> spacing as in patient 1 1 0 --> detector surface spacing CORRECTED for geometrical magnifications: spacing as in patient 1 1 1 --> detector surface spacing NOT corrected for geometrical magnifications: spacing as at detector */ mitk::PixelSpacingInterpretation mitk::DICOMImageBlockDescriptor ::GetPixelSpacingInterpretation() const { if ( m_ImageFrameList.empty() || m_TagCache.IsNull() ) { MITK_ERROR << "Invalid call to GetPixelSpacingInterpretation. Need to have initialized tag-cache!"; return SpacingUnknown; } std::string pixelSpacing = this->GetPixelSpacing(); std::string imagerPixelSpacing = this->GetImagerPixelSpacing(); if (pixelSpacing.empty()) { if (imagerPixelSpacing.empty()) { return SpacingUnknown; } else { return SpacingAtDetector; } } else // Pixel Spacing defined { if (imagerPixelSpacing.empty()) { return SpacingInPatient; } else if (pixelSpacing != imagerPixelSpacing) { return SpacingInPatient; } else { return SpacingAtDetector; } } } std::string mitk::DICOMImageBlockDescriptor ::GetPixelSpacing() const { if ( m_ImageFrameList.empty() || m_TagCache.IsNull() ) { MITK_ERROR << "Invalid call to GetPixelSpacing. Need to have initialized tag-cache!"; return std::string(""); } static const DICOMTag tagPixelSpacing(0x0028,0x0030); return m_TagCache->GetTagValue( m_ImageFrameList.front(), tagPixelSpacing ); } std::string mitk::DICOMImageBlockDescriptor ::GetImagerPixelSpacing() const { if ( m_ImageFrameList.empty() || m_TagCache.IsNull() ) { MITK_ERROR << "Invalid call to GetImagerPixelSpacing. Need to have initialized tag-cache!"; return std::string(""); } static const DICOMTag tagImagerPixelSpacing(0x0018,0x1164); return m_TagCache->GetTagValue( m_ImageFrameList.front(), tagImagerPixelSpacing ); } void mitk::DICOMImageBlockDescriptor ::GetDesiredMITKImagePixelSpacing( ScalarType& spacingX, ScalarType& spacingY) const { std::string pixelSpacing = this->GetPixelSpacing(); // preference for "in patient" pixel spacing if ( !DICOMStringToSpacing( pixelSpacing, spacingX, spacingY ) ) { std::string imagerPixelSpacing = this->GetImagerPixelSpacing(); // fallback to "on detector" spacing if ( !DICOMStringToSpacing( imagerPixelSpacing, spacingX, spacingY ) ) { // last resort: invent something spacingX = spacingY = 1.0; } } } void mitk::DICOMImageBlockDescriptor ::SetProperty(const std::string& key, BaseProperty* value) { m_PropertyList->SetProperty(key, value); } mitk::BaseProperty* mitk::DICOMImageBlockDescriptor ::GetProperty(const std::string& key) const { this->UpdateImageDescribingProperties(); return m_PropertyList->GetProperty(key); } std::string mitk::DICOMImageBlockDescriptor ::GetPropertyAsString(const std::string& key) const { this->UpdateImageDescribingProperties(); mitk::BaseProperty::Pointer property = m_PropertyList->GetProperty(key); if (property.IsNotNull()) { return property->GetValueAsString(); } else { return std::string(""); } } void mitk::DICOMImageBlockDescriptor ::SetFlag(const std::string& key, bool value) { m_PropertyList->ReplaceProperty(key, BoolProperty::New(value)); } bool mitk::DICOMImageBlockDescriptor ::GetFlag(const std::string& key, bool defaultValue) const { this->UpdateImageDescribingProperties(); BoolProperty::ConstPointer boolProp = dynamic_cast( this->GetProperty(key) ); if (boolProp.IsNotNull()) { return boolProp->GetValue(); } else { return defaultValue; } } void mitk::DICOMImageBlockDescriptor ::SetIntProperty(const std::string& key, int value) { m_PropertyList->ReplaceProperty(key, IntProperty::New(value)); } int mitk::DICOMImageBlockDescriptor ::GetIntProperty(const std::string& key, int defaultValue) const { this->UpdateImageDescribingProperties(); IntProperty::ConstPointer intProp = dynamic_cast( this->GetProperty(key) ); if (intProp.IsNotNull()) { return intProp->GetValue(); } else { return defaultValue; } } double mitk::DICOMImageBlockDescriptor ::stringtodouble(const std::string& str) const { double d; std::string trimmedstring(str); trimmedstring = trimmedstring.erase(trimmedstring.find_last_not_of(" \n\r\t")+1); std::istringstream converter(trimmedstring); if ( !trimmedstring.empty() && (converter >> d) && converter.eof() ) { return d; } else { throw std::invalid_argument("Argument is not a convertable number"); } } mitk::Image::Pointer mitk::DICOMImageBlockDescriptor ::DescribeImageWithProperties(Image* mitkImage) { // TODO: this is a collection of properties that have been provided by the // legacy DicomSeriesReader. // We should at some point clean up this collection and name them in a more // consistent way! if (!mitkImage) return mitkImage; // first part: add some tags that describe individual slices // these propeties are defined at analysis time (see UpdateImageDescribingProperties()) std::string propertyKeySliceLocation = "dicom.image.0020.1041"; std::string propertyKeyInstanceNumber = "dicom.image.0020.0013"; std::string propertyKeySOPInstanceUID = "dicom.image.0008.0018"; mitkImage->SetProperty( propertyKeySliceLocation.c_str(), this->GetProperty("sliceLocationForSlices") ); mitkImage->SetProperty( propertyKeyInstanceNumber.c_str(), this->GetProperty("instanceNumberForSlices") ); mitkImage->SetProperty( propertyKeySOPInstanceUID.c_str(), this->GetProperty("SOPInstanceUIDForSlices") ); // second part: add properties that describe the whole image block mitkImage->SetProperty("dicomseriesreader.SOPClassUID", StringProperty::New( this->GetSOPClassUID() ) ); mitkImage->SetProperty("dicomseriesreader.SOPClass", StringProperty::New( this->GetSOPClassUIDAsName() )); mitkImage->SetProperty("dicomseriesreader.PixelSpacingInterpretationString", StringProperty::New( PixelSpacingInterpretationToString( this->GetPixelSpacingInterpretation() )) ); mitkImage->SetProperty("dicomseriesreader.PixelSpacingInterpretation", GenericProperty::New( this->GetPixelSpacingInterpretation() )); mitkImage->SetProperty("dicomseriesreader.ReaderImplementationLevelString", StringProperty::New( ReaderImplementationLevelToString( m_ReaderImplementationLevel ) )); mitkImage->SetProperty("dicomseriesreader.ReaderImplementationLevel", GenericProperty::New( m_ReaderImplementationLevel )); mitkImage->SetProperty("dicomseriesreader.GantyTiltCorrected", BoolProperty::New( this->GetTiltInformation().IsRegularGantryTilt() )); mitkImage->SetProperty("dicomseriesreader.3D+t", BoolProperty::New( this->GetFlag("3D+t",false) )); // level window std::string windowCenter = this->GetPropertyAsString("windowCenter"); std::string windowWidth = this->GetPropertyAsString("windowWidth"); try { double level = stringtodouble( windowCenter ); double window = stringtodouble( windowWidth ); mitkImage->SetProperty("levelwindow", LevelWindowProperty::New(LevelWindow(level,window)) ); } catch (...) { // nothing, no levelwindow to be predicted... } mitkImage->SetProperty("dicom.pixel.PhotometricInterpretation", this->GetProperty("photometricInterpretation") ); mitkImage->SetProperty("dicom.image.imagetype", this->GetProperty("imagetype") ); mitkImage->SetProperty("dicom.study.StudyDescription", this->GetProperty("studyDescription") ); mitkImage->SetProperty("dicom.series.SeriesDescription", this->GetProperty("seriesDescription") ); // third part: get something from ImageIO. BUT this needs to be created elsewhere. or not at all! return mitkImage; } void mitk::DICOMImageBlockDescriptor ::SetReaderImplementationLevel(const ReaderImplementationLevel& level) { m_ReaderImplementationLevel = level; } mitk::ReaderImplementationLevel mitk::DICOMImageBlockDescriptor ::GetReaderImplementationLevel() const { return m_ReaderImplementationLevel; } std::string mitk::DICOMImageBlockDescriptor ::GetSOPClassUID() const { if ( !m_ImageFrameList.empty() && m_TagCache.IsNotNull() ) { static const DICOMTag tagSOPClassUID(0x0008,0x0016); return m_TagCache->GetTagValue( m_ImageFrameList.front(), tagSOPClassUID ); } else { MITK_ERROR << "Invalid call to DICOMImageBlockDescriptor::GetSOPClassUID(). Need to have initialized tag-cache!"; return std::string(""); } } std::string mitk::DICOMImageBlockDescriptor ::GetSOPClassUIDAsName() const { if ( !m_ImageFrameList.empty() && m_TagCache.IsNotNull() ) { gdcm::UIDs uidKnowledge; uidKnowledge.SetFromUID( this->GetSOPClassUID().c_str() ); const char* name = uidKnowledge.GetName(); if (name) { return std::string(name); } else { return std::string(""); } } else { MITK_ERROR << "Invalid call to DICOMImageBlockDescriptor::GetSOPClassUIDAsName(). Need to have initialized tag-cache!"; return std::string(""); } } void mitk::DICOMImageBlockDescriptor ::SetTagCache(DICOMTagCache* privateCache) { // this must only be used during loading and never afterwards m_TagCache = privateCache; } #define printPropertyRange(label, property_name) \ { \ std::string first = this->GetPropertyAsString( #property_name "First"); \ std::string last = this->GetPropertyAsString( #property_name "Last"); \ if (!first.empty() || !last.empty()) \ { \ if (first == last) \ { \ os << " " label ": '" << first << "'" << std::endl; \ } \ else \ { \ os << " " label ": '" << first << "' - '" << last << "'" << std::endl; \ } \ } \ } #define printProperty(label, property_name) \ { \ std::string first = this->GetPropertyAsString( #property_name ); \ if (!first.empty()) \ { \ os << " " label ": '" << first << "'" << std::endl; \ } \ } #define printBool(label, commands) \ { \ os << " " label ": '" << (commands?"yes":"no") << "'" << std::endl; \ } void mitk::DICOMImageBlockDescriptor ::Print(std::ostream& os, bool filenameDetails) const { os << " Number of Frames: '" << m_ImageFrameList.size() << "'" << std::endl; os << " SOP class: '" << this->GetSOPClassUIDAsName() << "'" << std::endl; printProperty("Series Number", seriesNumber); printProperty("Study Description", studyDescription); printProperty("Series Description", seriesDescription); printProperty("Modality", modality); printProperty("Sequence Name", sequenceName); printPropertyRange("Slice Location", sliceLocation); printPropertyRange("Acquisition Number", acquisitionNumber); printPropertyRange("Instance Number", instanceNumber); printPropertyRange("Image Position", imagePositionPatient); printProperty("Image Orientation", orientation); os << " Pixel spacing interpretation: '" << PixelSpacingInterpretationToString(this->GetPixelSpacingInterpretation()) << "'" << std::endl; printBool("Gantry Tilt", this->GetTiltInformation().IsRegularGantryTilt()) //printBool("3D+t", this->GetFlag("3D+t",false)) //os << " MITK image loaded: '" << (this->GetMitkImage().IsNotNull() ? "yes" : "no") << "'" << std::endl; if (filenameDetails) { os << " Files in this image block:" << std::endl; for (DICOMImageFrameList::const_iterator frameIter = m_ImageFrameList.begin(); frameIter != m_ImageFrameList.end(); ++frameIter) { os << " " << (*frameIter)->Filename; if ((*frameIter)->FrameNo > 0) { os << ", " << (*frameIter)->FrameNo; } os << std::endl; } } } #define storeTagValueToProperty(tag_name, tag_g, tag_e) \ { \ const DICOMTag t(tag_g, tag_e); \ std::string tagValue = m_TagCache->GetTagValue( firstFrame, t ); \ const_cast(this)->SetProperty(#tag_name, StringProperty::New( tagValue ) ); \ } #define storeTagValueRangeToProperty(tag_name, tag_g, tag_e) \ { \ const DICOMTag t(tag_g, tag_e); \ std::string tagValueFirst = m_TagCache->GetTagValue( firstFrame, t ); \ std::string tagValueLast = m_TagCache->GetTagValue( lastFrame, t ); \ const_cast(this)->SetProperty(#tag_name "First", StringProperty::New( tagValueFirst ) ); \ const_cast(this)->SetProperty(#tag_name "Last", StringProperty::New( tagValueLast ) ); \ } void mitk::DICOMImageBlockDescriptor ::UpdateImageDescribingProperties() const { if (!m_PropertiesOutOfDate) return; if (!m_ImageFrameList.empty()) { if (m_TagCache.IsNull()) { MITK_ERROR << "Invalid call to DICOMImageBlockDescriptor::UpdateImageDescribingProperties(). Need to have initialized tag-cache!"; return; } DICOMImageFrameInfo::Pointer firstFrame = m_ImageFrameList.front(); DICOMImageFrameInfo::Pointer lastFrame = m_ImageFrameList.back(); // see macros above storeTagValueToProperty(seriesNumber,0x0020,0x0011) storeTagValueToProperty(studyDescription,0x0008,0x1030) storeTagValueToProperty(seriesDescription,0x0008,0x103e) storeTagValueToProperty(modality,0x0008,0x0060) storeTagValueToProperty(sequenceName,0x0018,0x0024) storeTagValueToProperty(orientation,0x0020,0x0037) storeTagValueRangeToProperty(sliceLocation,0x0020,0x1041) storeTagValueRangeToProperty(acquisitionNumber,0x0020,0x0012) storeTagValueRangeToProperty(instanceNumber,0x0020,0x0013) storeTagValueRangeToProperty(imagePositionPatient,0x0020,0x0032) storeTagValueToProperty(windowCenter,0x0028,0x1050) storeTagValueToProperty(windowWidth,0x0028,0x1051) storeTagValueToProperty(imageType,0x0008,0x0008) storeTagValueToProperty(photometricInterpretation,0x0028,0x0004) // some per-image attributes // frames are just numbered starting from 0. timestep 1 (the second time-step) has frames starting at (number-of-frames-per-timestep) std::string propertyKeySliceLocation = "dicom.image.0020.1041"; std::string propertyKeyInstanceNumber = "dicom.image.0020.0013"; std::string propertyKeySOPInstanceNumber = "dicom.image.0008.0018"; StringLookupTable sliceLocationForSlices; StringLookupTable instanceNumberForSlices; StringLookupTable SOPInstanceUIDForSlices; const DICOMTag tagSliceLocation(0x0020,0x1041); const DICOMTag tagInstanceNumber(0x0020,0x0013); const DICOMTag tagSOPInstanceNumber(0x0008,0x0018); unsigned int slice(0); for (DICOMImageFrameList::const_iterator frameIter = m_ImageFrameList.begin(); frameIter != m_ImageFrameList.end(); ++slice, ++frameIter) { std::string sliceLocation = m_TagCache->GetTagValue( *frameIter, tagSliceLocation ); sliceLocationForSlices.SetTableValue(slice, sliceLocation); std::string instanceNumber = m_TagCache->GetTagValue( *frameIter, tagInstanceNumber ); instanceNumberForSlices.SetTableValue(slice, instanceNumber); std::string sopInstanceUID = m_TagCache->GetTagValue( *frameIter, tagSOPInstanceNumber ); SOPInstanceUIDForSlices.SetTableValue(slice, sopInstanceUID); MITK_DEBUG << "Tag info for slice " << slice << ": SL '" << sliceLocation << "' IN '" << instanceNumber << "' SOP instance UID '" << sopInstanceUID << "'"; // add property or properties with proper names const_cast(this)->SetProperty( "sliceLocationForSlices", StringLookupTableProperty::New( sliceLocationForSlices ) ); const_cast(this)->SetProperty( "instanceNumberForSlices", StringLookupTableProperty::New( instanceNumberForSlices ) ); const_cast(this)->SetProperty( "SOPInstanceUIDForSlices", StringLookupTableProperty::New( SOPInstanceUIDForSlices ) ); } m_PropertiesOutOfDate = false; } } diff --git a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.h b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.h index f5a30f4e4c..e3e060b6ba 100644 --- a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.h +++ b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.h @@ -1,161 +1,163 @@ /*=================================================================== 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 mitkDICOMImageBlockDescriptor_h #define mitkDICOMImageBlockDescriptor_h #include "mitkDICOMEnums.h" #include "mitkDICOMImageFrameInfo.h" #include "mitkDICOMTag.h" #include "mitkDICOMTagCache.h" #include "mitkImage.h" #include "mitkProperties.h" #include "mitkWeakPointer.h" #include "mitkGantryTiltInformation.h" namespace mitk { /** \ingroup DICOMReaderModule \brief Output descriptor for DICOMFileReader. As a result of analysis by a mitk::DICOMFileReader, this class describes the properties of a single mitk::Images that could be loaded by the file reader. The descriptor contains the following information: - the mitk::Image itself. This will be NULL after analysis and only be present after actual loading. - a list of frames (mostly: filenames) that went into composition of the mitk::Image. - an assessment of the reader's ability to load this set of files (ReaderImplementationLevel) - this can be used for reader selection when one reader is able to load an image with correct colors and the other is able to produce only gray values, for example - description of aspects of the image. Mostly a key-value list implemented by means of mitk::PropertyList. - for specific keys and possible values, see documentation of specific readers. \note an mitk::Image may both consist of multiple files (the "old" DICOM way) or a mitk::Image may be described by a single DICOM file or even only parts of a DICOM file (the newer multi-frame DICOM classes). To reflect this DICOMImageFrameList describes a list of frames from different or a single file. Described aspects of an image are: - whether pixel spacing is meant to be in-patient or on-detector (mitk::PixelSpacingInterpretation) - details about a possible gantry tilt (intended for use by file readers, may be hidden later) */ class DICOMReader_EXPORT DICOMImageBlockDescriptor { public: DICOMImageBlockDescriptor(); ~DICOMImageBlockDescriptor(); DICOMImageBlockDescriptor(const DICOMImageBlockDescriptor& other); DICOMImageBlockDescriptor& operator=(const DICOMImageBlockDescriptor& other); + static DICOMTagList GetTagsOfInterest(); + /// List of frames that constitute the mitk::Image (DICOMImageFrame%s) void SetImageFrameList(const DICOMImageFrameList& framelist); /// List of frames that constitute the mitk::Image (DICOMImageFrame%s) const DICOMImageFrameList& GetImageFrameList() const; /// The 3D mitk::Image that is loaded from the DICOM files of a DICOMImageFrameList void SetMitkImage(Image::Pointer image); /// the 3D mitk::Image that is loaded from the DICOM files of a DICOMImageFrameList Image::Pointer GetMitkImage() const; /// Reader's capability to appropriately load this set of frames ReaderImplementationLevel GetReaderImplementationLevel() const; /// Reader's capability to appropriately load this set of frames void SetReaderImplementationLevel(const ReaderImplementationLevel& level); /// Key-value store describing aspects of the image to be loaded void SetProperty(const std::string& key, BaseProperty* value); /// Key-value store describing aspects of the image to be loaded BaseProperty* GetProperty(const std::string& key) const; /// Convenience function around GetProperty() std::string GetPropertyAsString(const std::string&) const; /// Convenience function around SetProperty() void SetFlag(const std::string& key, bool value); /// Convenience function around GetProperty() bool GetFlag(const std::string& key, bool defaultValue) const; /// Convenience function around SetProperty() void SetIntProperty(const std::string& key, int value); /// Convenience function around GetProperty() int GetIntProperty(const std::string& key, int defaultValue) const; private: // For future implementation: load slice-by-slice, mark this using these methods void SetSliceIsLoaded(unsigned int index, bool isLoaded); // For future implementation: load slice-by-slice, mark this using these methods bool IsSliceLoaded(unsigned int index) const; // For future implementation: load slice-by-slice, mark this using these methods bool AllSlicesAreLoaded() const; public: /// Describe how the mitk::Image's pixel spacing should be interpreted PixelSpacingInterpretation GetPixelSpacingInterpretation() const; /// Describe the correct x/y pixel spacing of the mitk::Image (which some readers might need to adjust after loading) void GetDesiredMITKImagePixelSpacing( ScalarType& spacingXinMM, ScalarType& spacingYinMM) const; /// Describe the gantry tilt of the acquisition void SetTiltInformation(const GantryTiltInformation& info); /// Describe the gantry tilt of the acquisition const GantryTiltInformation GetTiltInformation() const; /// SOP Class UID of this set of frames void SetSOPClassUID(const std::string& uid); /// SOP Class UID of this set of frames std::string GetSOPClassUID() const; /// SOP Class as human readable name (e.g. "CT Image Storage") std::string GetSOPClassUIDAsName() const; void SetTagCache(DICOMTagCache* privateCache); /// Print information about this image block to given stream void Print(std::ostream& os, bool filenameDetails) const; private: // read values from tag cache std::string GetPixelSpacing() const; std::string GetImagerPixelSpacing() const; Image::Pointer FixupSpacing(Image* mitkImage); Image::Pointer DescribeImageWithProperties(Image* mitkImage); void UpdateImageDescribingProperties() const; double stringtodouble(const std::string& str) const; DICOMImageFrameList m_ImageFrameList; Image::Pointer m_MitkImage; BoolList m_SliceIsLoaded; ReaderImplementationLevel m_ReaderImplementationLevel; GantryTiltInformation m_TiltInformation; PropertyList::Pointer m_PropertyList; mitk::WeakPointer m_TagCache; mutable bool m_PropertiesOutOfDate; }; } #endif diff --git a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp index b65efa5ff0..c840004922 100644 --- a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp @@ -1,599 +1,571 @@ /*=================================================================== 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 + tags.push_back( DICOMTag(0x0018, 0x1120) ); // GantryDetectorTilt return tags; } void mitk::EquiDistantBlocksSorter ::Sort() { DICOMDatasetList remainingInput = GetInput(); // copy typedef std::list OutputListType; OutputListType outputs; m_SliceGroupingResults.clear(); while (!remainingInput.empty()) // repeat until all files are grouped somehow { SliceGroupingAnalysisResult regularBlock = this->AnalyzeFileForITKImageSeriesReaderSpacingAssumption( remainingInput, m_AcceptTilt ); 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; - } + 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 // 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; }