diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp index 9939060689..2d523b04d7 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp @@ -1,563 +1,571 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #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() +::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() +::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()); 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; 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 -::LoadMitkImageForOutput(unsigned int o) +::LoadMitkImageForImageBlockDescriptor(DICOMImageBlockDescriptor& block) const { PushLocale(); - DICOMImageBlockDescriptor& block = this->InternalGetOutput(o); const DICOMImageFrameList& frames = block.GetImageFrameList(); const GantryTiltInformation tiltInfo = block.GetTiltInformation(); bool hasTilt = block.GetFlag("gantryTilt", false); - if (hasTilt) - { - MITK_DEBUG << "When loading image " << o << ": got tilt info:"; - //tiltInfo.Print(std::cout); - } - else - { - MITK_DEBUG << "When loading image " << o << ": has NO info."; - } ITKDICOMSeriesReaderHelper::StringContainer filenames; for (DICOMImageFrameList::const_iterator frameIter = frames.begin(); frameIter != frames.end(); ++frameIter) { filenames.push_back( (*frameIter)->Filename ); } mitk::ITKDICOMSeriesReaderHelper helper; - mitk::Image::Pointer mitkImage = helper.Load( filenames, m_FixTiltByShearing && hasTilt, tiltInfo ); // TODO preloaded images, caching..? + 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(); + } - block.SetMitkImage( mitkImage ); PopLocale(); - return true; + 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 d61bc7ae4f..0ca1016332 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h @@ -1,327 +1,329 @@ /*=================================================================== 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(); + 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(); + 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; protected: // NOT nice, made available to ThreeDnTDICOMSeriesReader due to lack of time mitk::EquiDistantBlocksSorter::Pointer m_EquiDistantBlocksSorter; mitk::NormalDirectionConsistencySorter::Pointer m_NormalDirectionConsistencySorter; private: - std::stack m_ReplacedCLocales; - std::stack m_ReplacedCinLocales; + mutable std::stack m_ReplacedCLocales; + mutable std::stack m_ReplacedCinLocales; DICOMGDCMImageFrameList m_InputFrameList; gdcm::Scanner m_GDCMScanner; }; } #endif diff --git a/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.cpp b/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.cpp index 0117b8621c..4145ce7acd 100644 --- a/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.cpp +++ b/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.cpp @@ -1,246 +1,235 @@ /*=================================================================== 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) { DICOMGDCMImageFrameList& gdcmFrameInfoList = *blockIter; DICOMImageFrameList frameList = ToDICOMImageFrameList( gdcmFrameInfoList ); assert(!gdcmFrameInfoList.empty()); 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; } -// TODO why not handle 3D as a special case of 3D+t?? Late insight?.. .. because of ITK VImageDimension? .. could be a parameter bool mitk::ThreeDnTDICOMSeriesReader -::LoadMitkImageForOutput(unsigned int o) +::LoadMitkImageForImageBlockDescriptor(DICOMImageBlockDescriptor& block) const { PushLocale(); - DICOMImageBlockDescriptor& block = this->InternalGetOutput(o); const DICOMImageFrameList& frames = block.GetImageFrameList(); const GantryTiltInformation tiltInfo = block.GetTiltInformation(); bool hasTilt = block.GetFlag("gantryTilt", false); - if (hasTilt) - { - MITK_DEBUG << "When loading image " << o << ": got tilt info:"; - //tiltInfo.Print(std::cout); - } - else - { - MITK_DEBUG << "When loading image " << o << ": has NO info."; - } int numberOfTimesteps = block.GetIntProperty("timesteps", 1); 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; } diff --git a/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.h b/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.h index c10b908a1b..6e8ada9dc0 100644 --- a/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.h +++ b/Modules/DICOMReader/mitkThreeDnTDICOMSeriesReader.h @@ -1,83 +1,83 @@ /*=================================================================== 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 mitkThreeDnTDICOMSeriesReader_h #define mitkThreeDnTDICOMSeriesReader_h #include "mitkDICOMITKSeriesGDCMReader.h" #include "DICOMReaderExports.h" namespace mitk { /** \ingroup DICOMReader \brief Extends DICOMITKSeriesGDCMReader by sorting/grouping into 3D+t image blocks. This class reuses the DICOMITKSeriesGDCMReader class and adds the option of grouping 3D blocks at the same spatial position into a single block, which is loaded as a 3D+t mitk::Image (compare to \ref DICOMITKSeriesGDCMReader_Condensing). To group two output blocks into a single 3D+t block, this class tests a number of requirements that the two blocks must fulfill: - the origin of the first slice must be identical - the origin of the last slice must be identical - the number of slices must be identical The output blocks described by DICOMImageBlockDescriptor will contains the following properties: - \b "3D+t": true if the image is 3D+t - \b "timesteps": number of timesteps of an image (only defined if "3D+t" is true) */ class DICOMReader_EXPORT ThreeDnTDICOMSeriesReader : public DICOMITKSeriesGDCMReader { public: mitkClassMacro( ThreeDnTDICOMSeriesReader, DICOMITKSeriesGDCMReader ); mitkCloneMacro( ThreeDnTDICOMSeriesReader ); itkNewMacro( ThreeDnTDICOMSeriesReader ); /// \brief Control whether 3D+t grouping shall actually be attempted. void SetGroup3DandT(bool on); // void AllocateOutputImages(); /// \brief Load via multiple calls to itk::ImageSeriesReader. virtual bool LoadImages(); protected: ThreeDnTDICOMSeriesReader(); virtual ~ThreeDnTDICOMSeriesReader(); ThreeDnTDICOMSeriesReader(const ThreeDnTDICOMSeriesReader& other); ThreeDnTDICOMSeriesReader& operator=(const ThreeDnTDICOMSeriesReader& other); /** \brief Analyze the groups produced by DICOMITKSeriesGDCMReader for 3D+t properties. This method tests whether some blocks are at the same spatial position and groups them into single blocks. */ virtual SortingBlockList Condense3DBlocks(SortingBlockList&); - bool LoadMitkImageForOutput(unsigned int o); + bool LoadMitkImageForImageBlockDescriptor(DICOMImageBlockDescriptor& block) const; bool m_Group3DandT; }; } #endif