diff --git a/Modules/DICOMReader/TODOs.txt b/Modules/DICOMReader/TODOs.txt index a6e2b3eee8..bd0f2ccdc9 100644 --- a/Modules/DICOMReader/TODOs.txt +++ b/Modules/DICOMReader/TODOs.txt @@ -1,37 +1,49 @@ Important -[ ] ONLY load a pre-sorted list +[...] ONLY load a pre-sorted list + - Should work now by: + - get the list of images, consturct FrameInfo list (order must be known from prior sorting) + - prepare a tagcache object that provides tag information (from some kind of database) + - create DICOMImageBlockDescriptor with frame list and tag cache + - Call SetFixTiltByShearing() and TiltInfo object as appropriate... (the second part should not be necessary..) + - call DICOMITKSeriesGDCMReader::LoadMitkImageForImageBlockDescriptor(block) [ ] Either NRRD cache option or slice-by-slice loading [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 +[ ] 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) [ ] Check input images fo DICOMness and non-multi-frameness [ ] Use CanReadFile() in selection! [x] Check TODO in mitkDICOMTag.cpp - operator< for DICOMTag has been re-implemented in a readable and safer way. [ ] 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) [ ] ... diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp index 2d523b04d7..096a58b6fc 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp @@ -1,571 +1,581 @@ /*=================================================================== 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() :DICOMFileReader() ,m_FixTiltByShearing(true) { this->EnsureMandatorySortersArePresent(); } mitk::DICOMITKSeriesGDCMReader ::DICOMITKSeriesGDCMReader(const DICOMITKSeriesGDCMReader& other ) :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() ) { this->EnsureMandatorySortersArePresent(); } 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(); } return *this; } void mitk::DICOMITKSeriesGDCMReader ::SetFixTiltByShearing(bool on) { m_FixTiltByShearing = on; } 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(); // 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); - m_SortingResultInProgress = this->InternalExecuteSortingStep(sorterIndex, m_NormalDirectionConsistencySorter.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; - DICOMImageFrameList frameList = ToDICOMImageFrameList( gdcmFrameInfoList ); 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 ); - - const GantryTiltInformation& tiltInfo = m_EquiDistantBlocksSorter->GetTiltInformation( (gdcmFrameInfoList.front())->GetFilenameIfAvailable() ); - block.SetFlag("gantryTilt", tiltInfo.IsRegularGantryTilt()); block.SetTiltInformation( tiltInfo ); static const DICOMTag tagPixelSpacing(0x0028,0x0030); static const DICOMTag tagImagerPixelSpacing(0x0018,0x1164); std::string pixelSpacingString = (gdcmFrameInfoList.front())->GetTagValueAsString( tagPixelSpacing ); std::string imagerPixelSpacingString = gdcmFrameInfoList.front()->GetTagValueAsString( tagImagerPixelSpacing ); block.SetPixelSpacingTagValues(pixelSpacingString, imagerPixelSpacingString); static const DICOMTag tagSOPClassUID(0x0008,0x0016); std::string sopClassUID = (gdcmFrameInfoList.front())->GetTagValueAsString( tagSOPClassUID ); block.SetSOPClassUID(sopClassUID); block.SetReaderImplementationLevel( this->GetReaderImplementationLevel(sopClassUID) ); 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; + 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 = block.GetFlag("gantryTilt", false); + 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 ); // TODO preloaded images, caching..? 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& itkNotUsed(filename)) { return true; // can handle all // TODO peek into file, check DCM // TODO nice-to-have: check multi-framedness } void mitk::DICOMITKSeriesGDCMReader ::AddSortingElement(DICOMDatasetSorter* sorter, bool atFront) { assert(sorter); if (atFront) { m_Sorter.push_front( sorter ); } else { m_Sorter.push_back( sorter ); } } void mitk::DICOMITKSeriesGDCMReader ::EnsureMandatorySortersArePresent() { // TODO do just once! Do it in c'tor and handle this step extra! 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(5) ); // Image Orientation (Patient) // TODO: configurable! 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(); } } 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 0ca1016332..aaf3c616d3 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h @@ -1,329 +1,327 @@ /*=================================================================== 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 // TODO should reject multi-frame // TODO tests seem to pass if reader creates NO output at all! // TODO preloaded volumes?? could be solved in a different way.. 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 - The output blocks described by DICOMImageBlockDescriptor will contains the following properties (compare DICOMImageBlockDescriptor): - - \b "gantryTilt": true if the image is has a gantry tilt that will be corrected during loading - \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(). 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 ); /** \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); /** \brief Controls whether to "fix" tilted acquisitions by shearing the output (see \ref DICOMITKSeriesGDCMReader_GantryTilt). */ void SetFixTiltByShearing(bool on); 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(); 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 Creates the required sorting steps described in \ref DICOMITKSeriesGDCMReader_ForcedConfiguration void EnsureMandatorySortersArePresent(); /// \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: 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 - 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; }; } #endif diff --git a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp index 52239c5da5..0cede4cb06 100644 --- a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp +++ b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp @@ -1,655 +1,655 @@ /*=================================================================== 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_PixelSpacing("") ,m_ImagerPixelSpacing("") ,m_SOPClassUID("") ,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_PixelSpacing( other.m_PixelSpacing ) ,m_ImagerPixelSpacing( other.m_ImagerPixelSpacing ) ,m_SOPClassUID( other.m_SOPClassUID ) ,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_PixelSpacing = other.m_PixelSpacing; m_ImagerPixelSpacing = other.m_ImagerPixelSpacing; m_SOPClassUID = other.m_SOPClassUID; 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; } 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) { 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; } void mitk::DICOMImageBlockDescriptor ::SetPixelSpacingTagValues(const std::string& pixelSpacing, const std::string& imagerPixelSpacing) { m_PixelSpacing = pixelSpacing; m_ImagerPixelSpacing = imagerPixelSpacing; } /* 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_PixelSpacing.empty()) { if (m_ImagerPixelSpacing.empty()) { return SpacingUnknown; } else { return SpacingAtDetector; } } else // Pixel Spacing defined { if (m_ImagerPixelSpacing.empty()) { return SpacingInPatient; } else if (m_PixelSpacing != m_ImagerPixelSpacing) { return SpacingInPatient; } else { return SpacingAtDetector; } } } void mitk::DICOMImageBlockDescriptor ::GetDesiredMITKImagePixelSpacing( ScalarType& spacingX, ScalarType& spacingY) const { // preference for "in patient" pixel spacing if ( !DICOMStringToSpacing( m_PixelSpacing, spacingX, spacingY ) ) { // fallback to "on detector" spacing if ( !DICOMStringToSpacing( m_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( m_SOPClassUID ) ); 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->GetFlag("gantryTilt",false) )); + 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; } void mitk::DICOMImageBlockDescriptor ::SetSOPClassUID(const std::string& uid) { m_SOPClassUID = uid; } std::string mitk::DICOMImageBlockDescriptor ::GetSOPClassUID() const { return m_SOPClassUID; } std::string mitk::DICOMImageBlockDescriptor ::GetSOPClassUIDAsName() const { gdcm::UIDs uidKnowledge; uidKnowledge.SetFromUID( m_SOPClassUID.c_str() ); const char* name = uidKnowledge.GetName(); if (name) { return std::string(name); } else { return std::string(""); } } void mitk::DICOMImageBlockDescriptor ::SetTagCache(DICOMTagCache* privateCache) { // TODO this must only be used during loading and never afterwards // TODO better: somehow make possible to have smartpointers here 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->GetFlag("gantryTilt",false)) + 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_TagCache && !m_ImageFrameList.empty()) { 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/mitkEquiDistantBlocksSorter.cpp b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp index 3288ebe75f..1bb3eb8035 100644 --- a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp @@ -1,501 +1,493 @@ /*=================================================================== 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) { } mitk::EquiDistantBlocksSorter ::EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other ) :DICOMDatasetSorter(other) ,m_AcceptTilt(false) { } mitk::EquiDistantBlocksSorter ::~EquiDistantBlocksSorter() { } void mitk::EquiDistantBlocksSorter ::PrintConfiguration(std::ostream& os, const std::string& indent) const { os << indent << "Sort into blocks of equidistant, well-aligned slices " << (m_AcceptTilt ? "(accepting a gantry tilt)" : "") << std::endl; } void mitk::EquiDistantBlocksSorter ::SetAcceptTilt(bool accept) { m_AcceptTilt = accept; } mitk::EquiDistantBlocksSorter& mitk::EquiDistantBlocksSorter ::operator=(const EquiDistantBlocksSorter& other) { if (this != &other) { DICOMDatasetSorter::operator=(other); m_AcceptTilt = other.m_AcceptTilt; } 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); } } -const mitk::GantryTiltInformation -mitk::EquiDistantBlocksSorter -::GetTiltInformation(const std::string& filename) -{ - for (ResultsList::iterator ri = m_SliceGroupingResults.begin(); - ri != m_SliceGroupingResults.end(); - ++ri) - { - SliceGroupingAnalysisResult& result = *ri; - - if (filename == result.GetFirstFilenameOfBlock()) - { - return result.GetTiltInfo(); - } - } - - return GantryTiltInformation(); // empty -} - 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; - 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) + // classic mode without tolerance! + bool adaptiveErrorTolerance = false; // TODO make an option + if (adaptiveErrorTolerance) + { + 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 = 0.005; + } 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() ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt + 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 ) { 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 ) { - DICOMDatasetList datasets = result.GetBlockDatasets(); - DICOMDatasetAccess* firstDataset = datasets.front(); - DICOMDatasetAccess* lastDataset = datasets.back(); - unsigned int numberOfSlicesApart = datasets.size() - 1; - - 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 = firstDataset->GetTagValueAsString( tagImageOrientation ); - bool orientationConversion(false); - DICOMStringToOrientationVectors( orientationValue, right, up, orientationConversion ); - - if (orientationConversion) + 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 ); - if (!firstOriginString.empty() && !lastOriginString.empty()) - { - bool firstOriginConversion(false); - bool lastOriginConversion(false); - - Point3D firstOrigin = DICOMStringToPoint3D( firstOriginString, firstOriginConversion ); - Point3D lastOrigin = DICOMStringToPoint3D( lastOriginString, lastOriginConversion ); - - if (firstOriginConversion && lastOriginConversion) - { - GantryTiltInformation updatedTiltInfo( firstOrigin, lastOrigin, right, up, numberOfSlicesApart ); - result.FlagGantryTilt(updatedTiltInfo); - } - } + 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 37e9083fbe..be110347bd 100644 --- a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h @@ -1,180 +1,182 @@ /*=================================================================== 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. 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); - const GantryTiltInformation GetTiltInformation(const std::string& filename); virtual void PrintConfiguration(std::ostream& os, const std::string& indent = "") 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; }; } #endif diff --git a/Modules/DICOMReader/mitkGantryTiltInformation.cpp b/Modules/DICOMReader/mitkGantryTiltInformation.cpp index d7e11d45d5..d4b1dea4f0 100644 --- a/Modules/DICOMReader/mitkGantryTiltInformation.cpp +++ b/Modules/DICOMReader/mitkGantryTiltInformation.cpp @@ -1,219 +1,258 @@ /*=================================================================== 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 "mitkGantryTiltInformation.h" +#include "mitkDICOMTag.h" + #include "mitkLogMacros.h" mitk::GantryTiltInformation::GantryTiltInformation() : m_ShiftUp(0.0) , m_ShiftRight(0.0) , m_ShiftNormal(0.0) , m_ITKAssumedSliceSpacing(0.0) , m_NumberOfSlicesApart(0) { } #define doublepoint(x) \ Point3Dd x; \ x[0] = x ## f[0]; \ x[1] = x ## f[1]; \ x[2] = x ## f[2]; #define doublevector(x) \ Vector3Dd x; \ x[0] = x ## f[0]; \ x[1] = x ## f[1]; \ x[2] = x ## f[2]; mitk::GantryTiltInformation::GantryTiltInformation( const Point3D& origin1f, const Point3D& origin2f, const Vector3D& rightf, const Vector3D& upf, unsigned int numberOfSlicesApart) : m_ShiftUp(0.0) , m_ShiftRight(0.0) , m_ShiftNormal(0.0) , m_NumberOfSlicesApart(numberOfSlicesApart) { assert(numberOfSlicesApart); doublepoint(origin1); doublepoint(origin2); doublevector(right); doublevector(up); // determine if slice 1 (imagePosition1 and imageOrientation1) and slice 2 can be in one orthogonal slice stack: // calculate a line from origin 1, directed along the normal of slice (calculated as the cross product of orientation 1) // check if this line passes through origin 2 /* Determine if line (imagePosition2 + l * normal) contains imagePosition1. Done by calculating the distance of imagePosition1 from line (imagePosition2 + l *normal) E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html squared distance = | (pointAlongNormal - origin2) x (origin2 - origin1) | ^ 2 / |pointAlongNormal - origin2| ^ 2 ( x meaning the cross product ) */ Vector3Dd normal = itk::CrossProduct(right, up); Point3Dd pointAlongNormal = origin2 + normal; double numerator = itk::CrossProduct( pointAlongNormal - origin2 , origin2 - origin1 ).GetSquaredNorm(); double denominator = (pointAlongNormal - origin2).GetSquaredNorm(); double distance = sqrt(numerator / denominator); if ( distance > 0.001 ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt { MITK_DEBUG << " Series seems to contain a tilted (or sheared) geometry"; MITK_DEBUG << " Distance of expected slice origin from actual slice origin: " << distance; MITK_DEBUG << " ==> storing this shift for later analysis:"; MITK_DEBUG << " v origin1: " << origin1; MITK_DEBUG << " v origin2: " << origin2; MITK_DEBUG << " v right: " << right; MITK_DEBUG << " v up: " << up; MITK_DEBUG << " v normal: " << normal; Point3Dd projectionRight = projectPointOnLine( origin1, origin2, right ); Point3Dd projectionNormal = projectPointOnLine( origin1, origin2, normal ); m_ShiftRight = (projectionRight - origin2).GetNorm(); m_ShiftNormal = (projectionNormal - origin2).GetNorm(); /* now also check to which side the image is shifted. Calculation e.g. from http://mathworld.wolfram.com/Point-PlaneDistance.html */ Point3Dd testPoint = origin1; Vector3Dd planeNormal = up; double signedDistance = ( planeNormal[0] * testPoint[0] + planeNormal[1] * testPoint[1] + planeNormal[2] * testPoint[2] - ( planeNormal[0] * origin2[0] + planeNormal[1] * origin2[1] + planeNormal[2] * origin2[2] ) ) / sqrt( planeNormal[0] * planeNormal[0] + planeNormal[1] * planeNormal[1] + planeNormal[2] * planeNormal[2] ); m_ShiftUp = signedDistance; m_ITKAssumedSliceSpacing = (origin2 - origin1).GetNorm(); // How do we now this is assumed? See header documentation for ITK code references //double itkAssumedSliceSpacing = sqrt( m_ShiftUp * m_ShiftUp + m_ShiftNormal * m_ShiftNormal ); MITK_DEBUG << " calculated from slices " << m_NumberOfSlicesApart << " slices apart"; MITK_DEBUG << " shift normal: " << m_ShiftNormal; MITK_DEBUG << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing; MITK_DEBUG << " shift up: " << m_ShiftUp; MITK_DEBUG << " shift right: " << m_ShiftRight; MITK_DEBUG << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535; } } +mitk::GantryTiltInformation +mitk::GantryTiltInformation +::MakeFromTagValues( + const std::string& origin1String, + const std::string& origin2String, + const std::string orientationString, + unsigned int numberOfSlicesApart) +{ + 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 + bool orientationConversion(false); + DICOMStringToOrientationVectors( orientationString, right, up, orientationConversion ); + + if (orientationConversion + && !origin1String.empty() && !origin2String.empty() + ) + { + bool firstOriginConversion(false); + bool lastOriginConversion(false); + + Point3D firstOrigin = DICOMStringToPoint3D( origin1String, firstOriginConversion ); + Point3D lastOrigin = DICOMStringToPoint3D( origin2String, lastOriginConversion ); + + if (firstOriginConversion && lastOriginConversion) + { + return GantryTiltInformation( firstOrigin, lastOrigin, right, up, numberOfSlicesApart ); + } + } + + std::stringstream ss; + ss << "Invalid tag values when constructing tilt information from origin1 '" << origin1String + << "', origin2 '" << origin2String + << "', and orientation '" << orientationString << "'"; + + throw std::invalid_argument(ss.str()); +} + void mitk::GantryTiltInformation ::Print(std::ostream& os) const { os << " calculated from slices " << m_NumberOfSlicesApart << " slices apart" << std::endl; os << " shift normal: " << m_ShiftNormal << std::endl; os << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing << std::endl; os << " shift up: " << m_ShiftUp << std::endl; os << " shift right: " << m_ShiftRight << std::endl; os << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535 << std::endl; } mitk::Point3D mitk::GantryTiltInformation::projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ) { /** See illustration at http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage472/ vector(lineOrigin,p) = normal * ( innerproduct((p - lineOrigin),normal) / squared-length(normal) ) */ Vector3Dd lineOriginToP = p - lineOrigin; double innerProduct = lineOriginToP * lineDirection; double factor = innerProduct / lineDirection.GetSquaredNorm(); Point3Dd projection = lineOrigin + factor * lineDirection; return projection; } double mitk::GantryTiltInformation::GetTiltCorrectedAdditionalSize(unsigned int imageSizeZ) const { return fabs(m_ShiftUp / static_cast(m_NumberOfSlicesApart) * static_cast(imageSizeZ-1)); } double mitk::GantryTiltInformation::GetTiltAngleInDegrees() const { return atan( fabs(m_ShiftUp) / m_ShiftNormal ) * 180.0 / 3.1415926535; } double mitk::GantryTiltInformation::GetMatrixCoefficientForCorrectionInWorldCoordinates() const { // so many mm need to be shifted per slice! return m_ShiftUp / static_cast(m_NumberOfSlicesApart); } double mitk::GantryTiltInformation::GetRealZSpacing() const { return m_ShiftNormal / static_cast(m_NumberOfSlicesApart); } bool mitk::GantryTiltInformation::IsSheared() const { return m_NumberOfSlicesApart && ( fabs(m_ShiftRight) > 0.001 || fabs(m_ShiftUp) > 0.001); } bool mitk::GantryTiltInformation::IsRegularGantryTilt() const { return m_NumberOfSlicesApart && ( fabs(m_ShiftRight) < 0.001 && fabs(m_ShiftUp) > 0.001); } - - diff --git a/Modules/DICOMReader/mitkGantryTiltInformation.h b/Modules/DICOMReader/mitkGantryTiltInformation.h index 425818cc4e..18b4312ac4 100644 --- a/Modules/DICOMReader/mitkGantryTiltInformation.h +++ b/Modules/DICOMReader/mitkGantryTiltInformation.h @@ -1,133 +1,145 @@ /*=================================================================== 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 mitkGantryTiltInformation_h #define mitkGantryTiltInformation_h #include "mitkVector.h" namespace mitk { /** \ingroup DICOMReaderModule \brief Gantry tilt analysis result. Takes geometry information for two slices of a DICOM series and calculates whether these fit into an orthogonal block or not. If NOT, they can either be the result of an acquisition with gantry tilt OR completly broken by some shearing transformation. Most calculations are done in the constructor, results can then be read via the remaining methods. This class is a helper to DICOMITKSeriesGDCMReader and can not be used outside of \ref DICOMReaderModule */ class GantryTiltInformation { public: // two types to avoid any rounding errors typedef itk::Point Point3Dd; typedef itk::Vector Vector3Dd; /** \brief Just so we can create empty instances for assigning results later. */ GantryTiltInformation(); void Print(std::ostream& os) const; /** \brief THE constructor, which does all the calculations. Determining the amount of tilt is done by checking the distances of origin1 from planes through origin2. Two planes are considered: - normal vector along normal of slices (right x up): gives the slice distance - normal vector along orientation vector "up": gives the shift parallel to the plane orientation The tilt angle can then be calculated from these distances \param origin1 origin of the first slice \param origin2 origin of the second slice \param right right/up describe the orientatation of borth slices \param up right/up describe the orientatation of borth slices \param numberOfSlicesApart how many slices are the given origins apart (1 for neighboring slices) */ GantryTiltInformation( const Point3D& origin1, const Point3D& origin2, const Vector3D& right, const Vector3D& up, unsigned int numberOfSlicesApart); + /** + \brief Factory method to create GantryTiltInformation from tag values (strings). + + Parameters as the regular c'tor. + */ + static + GantryTiltInformation + MakeFromTagValues( const std::string& origin1String, + const std::string& origin2String, + const std::string orientationString, + unsigned int numberOfSlicesApart); + /** \brief Whether the slices were sheared. True if any of the shifts along right or up vector are non-zero. */ bool IsSheared() const; /** \brief Whether the shearing is a gantry tilt or more complicated. Gantry tilt will only produce shifts in ONE orientation, not in both. Since the correction code currently only coveres one tilt direction AND we don't know of medical images with two tilt directions, the loading code wants to check if our assumptions are true. */ bool IsRegularGantryTilt() const; /** \brief The offset distance in Y direction for each slice in mm (describes the tilt result). */ double GetMatrixCoefficientForCorrectionInWorldCoordinates() const; /** \brief The z / inter-slice spacing. Needed to correct ImageSeriesReader's result. */ double GetRealZSpacing() const; /** \brief The shift between first and last slice in mm. Needed to resize an orthogonal image volume. */ double GetTiltCorrectedAdditionalSize(unsigned int imageSizeZ) const; /** \brief Calculated tilt angle in degrees. */ double GetTiltAngleInDegrees() const; private: /** \brief Projection of point p onto line through lineOrigin in direction of lineDirection. */ Point3D projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ); double m_ShiftUp; double m_ShiftRight; double m_ShiftNormal; double m_ITKAssumedSliceSpacing; unsigned int m_NumberOfSlicesApart; }; } // namespace #endif diff --git a/Modules/DICOMReader/mitkNormalDirectionConsistencySorter.cpp b/Modules/DICOMReader/mitkNormalDirectionConsistencySorter.cpp index 9507acbada..40666d8a96 100644 --- a/Modules/DICOMReader/mitkNormalDirectionConsistencySorter.cpp +++ b/Modules/DICOMReader/mitkNormalDirectionConsistencySorter.cpp @@ -1,133 +1,165 @@ /*=================================================================== 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 "mitkNormalDirectionConsistencySorter.h" #include mitk::NormalDirectionConsistencySorter ::NormalDirectionConsistencySorter() :DICOMDatasetSorter() { } mitk::NormalDirectionConsistencySorter ::NormalDirectionConsistencySorter(const NormalDirectionConsistencySorter& other ) :DICOMDatasetSorter(other) { } mitk::NormalDirectionConsistencySorter ::~NormalDirectionConsistencySorter() { } void mitk::NormalDirectionConsistencySorter ::PrintConfiguration(std::ostream& os, const std::string& indent) const { - os << indent << "NormalDirectionConsistencySorter TODO" << std::endl; + os << indent << "NormalDirectionConsistencySorter" << std::endl; } mitk::NormalDirectionConsistencySorter& mitk::NormalDirectionConsistencySorter ::operator=(const NormalDirectionConsistencySorter& other) { if (this != &other) { DICOMDatasetSorter::operator=(other); } return *this; } mitk::DICOMTagList mitk::NormalDirectionConsistencySorter ::GetTagsOfInterest() { DICOMTagList tags; tags.push_back( DICOMTag(0x0020, 0x0032) ); // ImagePositionPatient tags.push_back( DICOMTag(0x0020, 0x0037) ); // ImageOrientationPatient return tags; } void mitk::NormalDirectionConsistencySorter ::Sort() { DICOMDatasetList datasets = GetInput(); if (datasets.size() > 1) { // at some point in the code, there is the expectation that // the direction of the slice normals is the same as the direction between // first and last slice origin. We need to make this sure here, because // we want to feed the files into itk::ImageSeriesReader with the consistent // setting of ReverseOrderOff. static const DICOMTag tagImagePositionPatient = DICOMTag(0x0020,0x0032); // Image Position (Patient) static const DICOMTag tagImageOrientation = DICOMTag(0x0020, 0x0037); // Image Orientation DICOMDatasetAccess* firstDS = datasets.front(); DICOMDatasetAccess* lastDS = datasets.back(); // make sure here that the direction from slice to slice is the direction of // image normals... std::string imageOrientationString = firstDS->GetTagValueAsString( tagImageOrientation ); std::string imagePositionPatientFirst = firstDS->GetTagValueAsString( tagImagePositionPatient ); std::string imagePositionPatientLast = lastDS->GetTagValueAsString( tagImagePositionPatient ); static Vector3D right; right.Fill(0.0); static Vector3D up; up.Fill(0.0); static bool hasOrientation(false); DICOMStringToOrientationVectors( imageOrientationString, right, up, hasOrientation ); static Point3D firstOrigin; firstOrigin.Fill(0.0f); static bool firstHasOrigin(false); firstOrigin = DICOMStringToPoint3D( imagePositionPatientFirst, firstHasOrigin ); static Point3D lastOrigin; lastOrigin.Fill(0.0f); static bool lastHasOrigin(false); lastOrigin = DICOMStringToPoint3D( imagePositionPatientLast, lastHasOrigin ); static Vector3D normal; normal[0] = right[1] * up[2] - right[2] * up[1]; normal[1] = right[2] * up[0] - right[0] * up[2]; normal[2] = right[0] * up[1] - right[1] * up[0]; normal.Normalize(); static Vector3D directionOfSlices; directionOfSlices = lastOrigin - firstOrigin; directionOfSlices.Normalize(); + static double projection = 0.0; + projection = 0.0; + projection = normal * directionOfSlices; + MITK_DEBUG << "Making sense of \norientation '" << imageOrientationString - << "'\nfirst position '" << imagePositionPatientFirst - << "'\nlast position '" << imagePositionPatientLast << "'"; + << "'\nfirst position '" << imagePositionPatientFirst + << "'\nlast position '" << imagePositionPatientLast << "'"; MITK_DEBUG << "Normal: " << normal; MITK_DEBUG << "Direction of slices: " << directionOfSlices; + MITK_DEBUG << "Projection of direction onto slice normal: " << projection; - if ( (directionOfSlices - normal).GetNorm() > 0.5 ) + if ( projection < 0.0 ) { MITK_DEBUG << "Need to reverse filenames"; std::reverse( datasets.begin(), datasets.end() ); + + m_TiltInfo = GantryTiltInformation::MakeFromTagValues( + imagePositionPatientLast, + imagePositionPatientFirst, + imageOrientationString, + datasets.size() - 1 + ); } + else + { + m_TiltInfo = GantryTiltInformation::MakeFromTagValues( + imagePositionPatientFirst, + imagePositionPatientLast, + imageOrientationString, + datasets.size() - 1 + ); + } + } + else // just ONE dataset, do not forget to reset tilt information + { + m_TiltInfo = GantryTiltInformation(); // empty info } this->SetNumberOfOutputs(1); this->SetOutput(0, datasets); } + +mitk::GantryTiltInformation +mitk::NormalDirectionConsistencySorter +::GetTiltInformation() const +{ + return m_TiltInfo; +} diff --git a/Modules/DICOMReader/mitkNormalDirectionConsistencySorter.h b/Modules/DICOMReader/mitkNormalDirectionConsistencySorter.h index 17b6c21210..b0dbfb478f 100644 --- a/Modules/DICOMReader/mitkNormalDirectionConsistencySorter.h +++ b/Modules/DICOMReader/mitkNormalDirectionConsistencySorter.h @@ -1,63 +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 mitkNormalDirectionConsistencySorter_h #define mitkNormalDirectionConsistencySorter_h #include "mitkDICOMDatasetSorter.h" +#include "mitkGantryTiltInformation.h" namespace mitk { /** \ingroup DICOMReaderModule \brief Makes sure that the order of files is along the image plane normals. When loading with ImageSeriesReader and initializing an mitk::Image with the result we need to make sure that the order of inputs for the ImageSeriesReader is along the normal of the images. I.e. The direction of the normal needs to be the same direction as the vector from the first to the last image origin. + Since this class is used as a last sorting step before loading, it will also + calculate (and return) an updated GantryTiltInformation object. + \note This class might be a workaround for another bug in MITK, but until this issue is completely understood, the workaround fixes the problem of images that appear upside-down. */ class DICOMReader_EXPORT NormalDirectionConsistencySorter : public DICOMDatasetSorter { public: mitkClassMacro( NormalDirectionConsistencySorter, DICOMDatasetSorter ) itkNewMacro( NormalDirectionConsistencySorter ) virtual DICOMTagList GetTagsOfInterest(); /// See class description. virtual void Sort(); + /// See class description and DICOMITKSeriesGDCMReader. + GantryTiltInformation GetTiltInformation() const; + virtual void PrintConfiguration(std::ostream& os, const std::string& indent = "") const; protected: NormalDirectionConsistencySorter(); virtual ~NormalDirectionConsistencySorter(); NormalDirectionConsistencySorter(const NormalDirectionConsistencySorter& other); NormalDirectionConsistencySorter& operator=(const NormalDirectionConsistencySorter& other); + + GantryTiltInformation m_TiltInfo; }; } #endif diff --git a/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.cpp b/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.cpp index 99d6291f1d..cd6c6c5740 100644 --- a/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.cpp +++ b/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.cpp @@ -1,241 +1,248 @@ /*=================================================================== 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 "mitkThreeDnTDICOMSeriesReader.h" #include "mitkITKDICOMSeriesReaderHelper.h" mitk::ThreeDnTDICOMSeriesReader ::ThreeDnTDICOMSeriesReader() :DICOMITKSeriesGDCMReader() ,m_Group3DandT(true) { } mitk::ThreeDnTDICOMSeriesReader ::ThreeDnTDICOMSeriesReader(const ThreeDnTDICOMSeriesReader& other ) :DICOMITKSeriesGDCMReader(other) ,m_Group3DandT(true) { } mitk::ThreeDnTDICOMSeriesReader ::~ThreeDnTDICOMSeriesReader() { } mitk::ThreeDnTDICOMSeriesReader& mitk::ThreeDnTDICOMSeriesReader ::operator=(const ThreeDnTDICOMSeriesReader& other) { if (this != &other) { DICOMITKSeriesGDCMReader::operator=(other); this->m_Group3DandT = other.m_Group3DandT; } return *this; } void mitk::ThreeDnTDICOMSeriesReader ::SetGroup3DandT(bool on) { m_Group3DandT = on; } mitk::DICOMITKSeriesGDCMReader::SortingBlockList mitk::ThreeDnTDICOMSeriesReader ::Condense3DBlocks(SortingBlockList& resultOf3DGrouping) { if (!m_Group3DandT) { return resultOf3DGrouping; // don't work if nobody asks us to } SortingBlockList remainingBlocks = resultOf3DGrouping; SortingBlockList non3DnTBlocks; SortingBlockList true3DnTBlocks; std::vector true3DnTBlocksTimeStepCount; // TODO we should provide this tag as needed via a function // (however, we currently know that the superclass will use this tag) const DICOMTag tagImagePositionPatient(0x0020, 0x0032); while (!remainingBlocks.empty()) { // new block to fill up DICOMGDCMImageFrameList& firstBlock = remainingBlocks.front(); DICOMGDCMImageFrameList current3DnTBlock = firstBlock; int current3DnTBlockNumberOfTimeSteps = 1; // get block characteristics of first block unsigned int currentBlockNumberOfSlices = firstBlock.size(); std::string currentBlockFirstOrigin = firstBlock.front()->GetTagValueAsString( tagImagePositionPatient ); std::string currentBlockLastOrigin = firstBlock.back()->GetTagValueAsString( tagImagePositionPatient ); remainingBlocks.pop_front(); // compare all other blocks against the first one for (SortingBlockList::iterator otherBlockIter = remainingBlocks.begin(); otherBlockIter != remainingBlocks.end(); /*++otherBlockIter*/) // <-- inside loop { // get block characteristics from first block DICOMGDCMImageFrameList& otherBlock = *otherBlockIter; unsigned int otherBlockNumberOfSlices = otherBlock.size(); std::string otherBlockFirstOrigin = otherBlock.front()->GetTagValueAsString( tagImagePositionPatient ); std::string otherBlockLastOrigin = otherBlock.back()->GetTagValueAsString( tagImagePositionPatient ); // add matching blocks to current3DnTBlock // keep other blocks for later if ( otherBlockNumberOfSlices == currentBlockNumberOfSlices && otherBlockFirstOrigin == currentBlockFirstOrigin && otherBlockLastOrigin == currentBlockLastOrigin ) { // matching block ++current3DnTBlockNumberOfTimeSteps; current3DnTBlock.insert( current3DnTBlock.end(), otherBlock.begin(), otherBlock.end() ); // append // remove this block from remainingBlocks otherBlockIter = remainingBlocks.erase(otherBlockIter); // make sure iterator otherBlockIter is valid afterwards } else { ++otherBlockIter; } } // in any case, we now now all about the first block of our list ... // ... and we wither call it 3D o 3D+t if (current3DnTBlockNumberOfTimeSteps > 1) { true3DnTBlocks.push_back(current3DnTBlock); true3DnTBlocksTimeStepCount.push_back(current3DnTBlockNumberOfTimeSteps); } else { non3DnTBlocks.push_back(current3DnTBlock); } } // create output for real 3D+t blocks (other outputs will be created by superclass) // set 3D+t flag on output block this->SetNumberOfOutputs( true3DnTBlocks.size() ); unsigned int o = 0; for (SortingBlockList::iterator blockIter = true3DnTBlocks.begin(); blockIter != true3DnTBlocks.end(); ++o, ++blockIter) { + // bad copy&paste code from DICOMITKSeriesGDCMReader, should be handled in a better way DICOMGDCMImageFrameList& gdcmFrameInfoList = *blockIter; - DICOMImageFrameList frameList = ToDICOMImageFrameList( gdcmFrameInfoList ); 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 ); - - // bad copy&paste code, should be handled in a better way - - const GantryTiltInformation& tiltInfo = m_EquiDistantBlocksSorter->GetTiltInformation( (gdcmFrameInfoList.front())->GetFilenameIfAvailable() ); block.SetTiltInformation( tiltInfo ); // assume static const DICOMTag tagPixelSpacing(0x0028,0x0030); static const DICOMTag tagImagerPixelSpacing(0x0018,0x1164); std::string pixelSpacingString = (gdcmFrameInfoList.front())->GetTagValueAsString( tagPixelSpacing ); std::string imagerPixelSpacingString = gdcmFrameInfoList.front()->GetTagValueAsString( tagImagerPixelSpacing ); block.SetPixelSpacingTagValues(pixelSpacingString, imagerPixelSpacingString); block.SetFlag("3D+t", true); block.SetIntProperty("timesteps", true3DnTBlocksTimeStepCount[o]); MITK_DEBUG << "Found " << true3DnTBlocksTimeStepCount[o] << " timesteps"; this->SetOutput( o, block ); } return non3DnTBlocks; } bool mitk::ThreeDnTDICOMSeriesReader ::LoadImages() { bool success = true; unsigned int numberOfOutputs = this->GetNumberOfOutputs(); for (unsigned int o = 0; o < numberOfOutputs; ++o) { DICOMImageBlockDescriptor& block = this->InternalGetOutput(o); if (block.GetFlag("3D+t", false)) { success &= this->LoadMitkImageForOutput(o); } else { success &= DICOMITKSeriesGDCMReader::LoadMitkImageForOutput(o); // let superclass handle non-3D+t } } return success; } bool mitk::ThreeDnTDICOMSeriesReader ::LoadMitkImageForImageBlockDescriptor(DICOMImageBlockDescriptor& block) const { PushLocale(); const DICOMImageFrameList& frames = block.GetImageFrameList(); const GantryTiltInformation tiltInfo = block.GetTiltInformation(); - bool hasTilt = block.GetFlag("gantryTilt", false); + bool hasTilt = tiltInfo.IsRegularGantryTilt(); int numberOfTimesteps = block.GetIntProperty("timesteps", 1); if (numberOfTimesteps == 1) { return DICOMITKSeriesGDCMReader::LoadMitkImageForImageBlockDescriptor(block); } int numberOfFramesPerTimestep = frames.size() / numberOfTimesteps; assert( int(double((double)frames.size() / (double)numberOfTimesteps )) == numberOfFramesPerTimestep ); // this should hold ITKDICOMSeriesReaderHelper::StringContainerList filenamesPerTimestep; for (int timeStep = 0; timeStepFilename ); } filenamesPerTimestep.push_back( filenamesOfThisTimeStep ); } mitk::ITKDICOMSeriesReaderHelper helper; mitk::Image::Pointer mitkImage = helper.Load3DnT( filenamesPerTimestep, m_FixTiltByShearing && hasTilt, tiltInfo ); // TODO preloaded images, caching..? block.SetMitkImage( mitkImage ); PopLocale(); return true; }