From 2a0643adca2256f9c11cda632fe81df2ceeb22b2 Mon Sep 17 00:00:00 2001 From: Daniel Maleike Date: Fri, 28 Jan 2011 10:03:21 +0100 Subject: [PATCH] FIX (#6506): sort by position first, then by acquisition time; adapt 3D+t loading to this change; sort files into groups of similar image properties --- mitk/Core/Code/IO/mitkDicomSeriesReader.cpp | 161 ++++++++++++++++-- mitk/Core/Code/IO/mitkDicomSeriesReader.h | 26 +++ mitk/Core/Code/IO/mitkDicomSeriesReader.txx | 236 ++++++++++++++++++--------- 3 files changed, 326 insertions(+), 97 deletions(-) diff --git a/mitk/Core/Code/IO/mitkDicomSeriesReader.cpp b/mitk/Core/Code/IO/mitkDicomSeriesReader.cpp index 733630c..a290a66 100644 --- a/mitk/Core/Code/IO/mitkDicomSeriesReader.cpp +++ b/mitk/Core/Code/IO/mitkDicomSeriesReader.cpp @@ -23,6 +23,8 @@ PURPOSE. See the above copyright notices for more information. #include #include #include + #include + #include #endif namespace mitk @@ -283,7 +285,7 @@ bool DicomSeriesReader::ReadPhilips3DDicom(const std::string &filename, mitk::Im DicomSeriesReader::UidFileNamesMap DicomSeriesReader::GetSeries(const std::string &dir, const StringContainer &restrictions) { - DcmFileNamesGeneratorType::Pointer name_generator = DcmFileNamesGeneratorType::New(); + UidFileNamesMap map; // result variable /** @@ -296,6 +298,10 @@ DicomSeriesReader::UidFileNamesMap DicomSeriesReader::GetSeries(const std::strin - 0018,0050 slice thickness */ + +#if GDCM_MAJOR_VERSION < 2 + DcmFileNamesGeneratorType::Pointer name_generator = DcmFileNamesGeneratorType::New(); + name_generator->SetUseSeriesDetails(true); name_generator->AddSeriesRestriction("0020|0037"); // image orientation (patient) name_generator->AddSeriesRestriction("0028|0030"); // pixel spacing (x,y) @@ -312,7 +318,6 @@ DicomSeriesReader::UidFileNamesMap DicomSeriesReader::GetSeries(const std::strin name_generator->SetDirectory(dir.c_str()); - UidFileNamesMap map; const StringContainer &series_uids = name_generator->GetSeriesUIDs(); const StringContainer::const_iterator series_end = series_uids.end(); @@ -323,9 +328,123 @@ DicomSeriesReader::UidFileNamesMap DicomSeriesReader::GetSeries(const std::strin map[uid] = name_generator->GetFileNames(uid); } +#else + + // set directory + gdcm::Directory gDirectory; + gDirectory.Load( dir.c_str(), false ); // non-recursive + const gdcm::Directory::FilenamesType &gAllFiles = gDirectory.GetFilenames(); + + // scann for relevant tags in dicom files + gdcm::Scanner scanner; + const gdcm::Tag tagSeriesInstanceUID(0x0020,0x000e); // Series Instance UID + scanner.AddTag( tagSeriesInstanceUID ); + + const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // image orientation + scanner.AddTag( tagImageOrientation ); + + const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing + scanner.AddTag( tagPixelSpacing ); + + const gdcm::Tag tagSliceThickness(0x0018, 0x0050); // slice thickness + scanner.AddTag( tagSliceThickness ); + + const gdcm::Tag tagNumberOfRows(0x0028, 0x0010); // number rows + scanner.AddTag( tagNumberOfRows ); + + const gdcm::Tag tagNumberOfColumns(0x0028, 0x0011); // number cols + scanner.AddTag( tagNumberOfColumns ); + + // TODO add further restrictions from arguments + + // let GDCM scan files + if ( !scanner.Scan( gDirectory.GetFilenames() ) ) + { + MITK_ERROR << "gdcm::Scanner failed scanning " << dir; + return map; + } + + // assign files IDs that will separate them for loading into image blocks + for (gdcm::Scanner::ConstIterator fileIter = scanner.Begin(); + fileIter != scanner.End(); + ++fileIter) + { + MITK_DEBUG << "File " << fileIter->first << std::endl; + if ( std::string(fileIter->first).empty() ) continue; // TODO understand why Scanner has empty string entries + + std::string moreUniqueSeriesId = CreateMoreUniqueSeriesIdentifier( fileIter->second ); + map [ moreUniqueSeriesId ].push_back( fileIter->first ); + } + +#endif + + for ( UidFileNamesMap::const_iterator i = map.begin(); i != map.end(); ++i ) + { + MITK_INFO << "Entry " << i->first << " with " << i->second.size() << " files"; + } + return map; } +#if GDCM_MAJOR_VERSION >= 2 +std::string DicomSeriesReader::CreateMoreUniqueSeriesIdentifier( const gdcm::Scanner::TagToValue& tagValueMap ) +{ + const gdcm::Tag tagSeriesInstanceUID(0x0020,0x000e); // Series Instance UID + const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // image orientation + const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing + const gdcm::Tag tagSliceThickness(0x0018, 0x0050); // slice thickness + const gdcm::Tag tagNumberOfRows(0x0028, 0x0010); // number rows + const gdcm::Tag tagNumberOfColumns(0x0028, 0x0011); // number cols + + std::string constructedID; + + try + { + constructedID = tagValueMap.at( tagSeriesInstanceUID ); + + constructedID += IDifyTagValue( tagValueMap.at(tagNumberOfRows) ); + constructedID += IDifyTagValue( tagValueMap.at(tagNumberOfColumns) ); + constructedID += IDifyTagValue( tagValueMap.at(tagPixelSpacing) ); + constructedID += IDifyTagValue( tagValueMap.at(tagSliceThickness) ); + constructedID += IDifyTagValue( tagValueMap.at(tagImageOrientation) ); + + constructedID.resize( constructedID.length() - 1 ); // cut of trailing '.' + + MITK_DEBUG << "ID: " << constructedID; + return constructedID; + } + catch (std::exception& e) + { + MITK_ERROR << "CreateMoreUniqueSeriesIdentifier() could not access all required DICOM tags. You are calling it wrongly. Using partial ID."; + return constructedID; + } +} + +std::string DicomSeriesReader::IDifyTagValue(const std::string& value) +{ + std::string IDifiedValue( value ); + + // Eliminate non-alnum characters, including whitespace... + // that may have been introduced by concats. + for(std::size_t i=0; i= 'a' && IDifiedValue[i] <= 'z') + || (IDifiedValue[i] >= '0' && IDifiedValue[i] <= '9') + || (IDifiedValue[i] >= 'A' && IDifiedValue[i] <= 'Z'))) + { + IDifiedValue.erase(i, 1); + } + } + + + IDifiedValue += "."; + return IDifiedValue; +} + +#endif + DicomSeriesReader::StringContainer DicomSeriesReader::GetSeries(const std::string &dir, const std::string &series_uid, const StringContainer &restrictions) { DcmFileNamesGeneratorType::Pointer name_generator = DcmFileNamesGeneratorType::New(); @@ -381,31 +500,35 @@ bool DicomSeriesReader::GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::D if (image_orientation1 != image_orientation2) { MITK_ERROR << "Dicom images have different orientations."; - throw std::logic_error("Dicom images have different orientations."); + throw std::logic_error("Dicom images have different orientations. Call GetSeries() first to separate images."); } - if (acq_time1 == acq_time2) - { - double normal[3]; + double normal[3]; - normal[0] = image_orientation1[1] * image_orientation1[5] - image_orientation1[2] * image_orientation1[4]; - normal[1] = image_orientation1[2] * image_orientation1[3] - image_orientation1[0] * image_orientation1[5]; - normal[2] = image_orientation1[0] * image_orientation1[4] - image_orientation1[1] * image_orientation1[3]; + normal[0] = image_orientation1[1] * image_orientation1[5] - image_orientation1[2] * image_orientation1[4]; + normal[1] = image_orientation1[2] * image_orientation1[3] - image_orientation1[0] * image_orientation1[5]; + normal[2] = image_orientation1[0] * image_orientation1[4] - image_orientation1[1] * image_orientation1[3]; - double - dist1 = 0.0, - dist2 = 0.0; + double + dist1 = 0.0, + dist2 = 0.0; - for (unsigned char i = 0u; i < 3u; ++i) - { - dist1 += normal[i] * image_pos1[i]; - dist2 += normal[i] * image_pos2[i]; - } + for (unsigned char i = 0u; i < 3u; ++i) + { + dist1 += normal[i] * image_pos1[i]; + dist2 += normal[i] * image_pos2[i]; + } + if ( fabs(dist1 - dist2) < mitk::eps) + { + // exception: same position: compare by acquisition time + return acq_time1 < acq_time2; + } + else + { + // default: compare position return dist1 < dist2; } - - return acq_time1 < acq_time2; } #endif diff --git a/mitk/Core/Code/IO/mitkDicomSeriesReader.h b/mitk/Core/Code/IO/mitkDicomSeriesReader.h index d442395..45b7136 100644 --- a/mitk/Core/Code/IO/mitkDicomSeriesReader.h +++ b/mitk/Core/Code/IO/mitkDicomSeriesReader.h @@ -85,6 +85,12 @@ public: * Read a Philips3D Ultrasound DICOM file and put into an mitk image */ static bool ReadPhilips3DDicom(const std::string &filename, mitk::Image::Pointer output_image); + + /** + Construct a UID that takes into account sorting criteria from GetSeries(). + */ + static std::string CreateMoreUniqueSeriesIdentifier( const gdcm::Scanner::TagToValue& tagValueMap ); + static std::string IDifyTagValue(const std::string& value); #endif /** @@ -169,6 +175,26 @@ protected: template static void LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, UpdateCallBackMethod callback); + /** + Feed files into itk::ImageSeriesReader and retrieve a 3D MITK image. + + \param command can be used for progress reporting + */ + template + static Image::Pointer LoadDICOMByITK( const StringContainer&, CallbackCommand* command = NULL); + + + /** + Called by LoadDicom. Expects to be fed a single list of filenames that have been sorted by + GetSeries previously (one map entry). This method will check how many timestep can be filled + with given files. + + Assumption is that the number of time steps is determined by how often the first position in + space repeats. I.e. if the first three files in the input parameter all describe the same + location in space, we'll construct three lists of files. and sort the remaining files into them. + */ + static std::list SortIntoBlocksFor3DplusT( const StringContainer& presortedFilenames, bool sort, bool& canLoadAs4D ); + #if GDCM_MAJOR_VERSION >= 2 /* * Auxiliary sort function for Gdcm Dicom sorting. It is used for sorting diff --git a/mitk/Core/Code/IO/mitkDicomSeriesReader.txx b/mitk/Core/Code/IO/mitkDicomSeriesReader.txx index 989290e..94f9538 100644 --- a/mitk/Core/Code/IO/mitkDicomSeriesReader.txx +++ b/mitk/Core/Code/IO/mitkDicomSeriesReader.txx @@ -52,48 +52,17 @@ void DicomSeriesReader::LoadDicom(const StringContainer &filenames, DataNode &no } /******** For 4D data split in multiple files ***************/ - std::list decomposed_filenames; - unsigned int volume_count = 1u; + bool canLoadAs4D(true); + std::list imageBlocks = SortIntoBlocksFor3DplusT( filenames, sort, canLoadAs4D ); + unsigned int volume_count = imageBlocks.size(); - StringContainer sorted_filenames = sort - ? DicomSeriesReader::SortSeriesSlices(filenames) - : filenames; - if (check_4d) + if (!canLoadAs4D) { - gdcm::Tag tag(0x0020,0x0032); //Image position (Patient) - gdcm::Scanner scanner; - - scanner.AddTag(tag); - scanner.Scan(sorted_filenames); - - const StringContainer::const_iterator f_end = sorted_filenames.end(); - const char *act_value = scanner.GetValue(sorted_filenames.front().c_str(), tag); - - decomposed_filenames.push_back(StringContainer()); - decomposed_filenames.back().push_back(sorted_filenames.front()); - - for (StringContainer::const_iterator f_it = ++sorted_filenames.begin(); f_it != f_end; ++f_it) - { - const char *value = scanner.GetValue(f_it->c_str(), tag); - - if (!strcmp(act_value, value)) - { - decomposed_filenames.push_back(StringContainer()); - ++volume_count; - } - - decomposed_filenames.back().push_back(*f_it); - } + image = LoadDICOMByITK( imageBlocks.front() , command ); // load first 3D block } else { - decomposed_filenames.push_back(sorted_filenames); - } - - if (volume_count > 1u) - { - // It is 4D! Read it and store into mitk image - + // It is 3D+t! Read it and store into mitk image typedef itk::Image ImageType; typedef itk::ImageSeriesReader ReaderType; @@ -108,16 +77,18 @@ void DicomSeriesReader::LoadDicom(const StringContainer &filenames, DataNode &no reader->AddObserver(itk::ProgressEvent(), command); } - const std::list::const_iterator df_end = decomposed_filenames.end(); + const std::list::const_iterator df_end = imageBlocks.end(); unsigned int act_volume = 1u; - reader->SetFileNames(decomposed_filenames.front()); + reader->SetFileNames(imageBlocks.front()); reader->Update(); image->InitializeByItk(reader->GetOutput(), 1, volume_count); image->SetImportVolume(reader->GetOutput()->GetBufferPointer(), 0u); - MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " << image->GetDimension(1) << ", " << image->GetDimension(2) << ", " << image->GetDimension(3) << "]"; - MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " << image->GetGeometry()->GetSpacing()[1] << ", " << image->GetGeometry()->GetSpacing()[2] << "]"; + MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " + << image->GetDimension(1) << ", " + << image->GetDimension(2) << ", " + << image->GetDimension(3) << "]"; #if (GDCM_MAJOR_VERSION == 2) && (GDCM_MINOR_VERSION < 1) && (GDCM_BUILD_VERSION < 15) // workaround for a GDCM 2 bug until version 2.0.15: @@ -132,71 +103,180 @@ void DicomSeriesReader::LoadDicom(const StringContainer &filenames, DataNode &no image->GetGeometry()->SetSpacing( correctedImageSpacing ); #endif + MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " + << image->GetGeometry()->GetSpacing()[1] << ", " + << image->GetGeometry()->GetSpacing()[2] << "]"; - for (std::list::iterator df_it = ++decomposed_filenames.begin(); df_it != df_end; ++df_it) + for (std::list::iterator df_it = ++imageBlocks.begin(); df_it != df_end; ++df_it) { reader->SetFileNames(*df_it); reader->Update(); image->SetImportVolume(reader->GetOutput()->GetBufferPointer(), act_volume++); } - node.SetData(image); - setlocale(LC_NUMERIC, previousCLocale); - std::cin.imbue(previousCppLocale); - return; // finished } + +#else + image = LoadDICOMByITK( filenames, command ); #endif - /******** Normal Case, 3D (also for GDCM < 2 usable) ***************/ - typedef itk::Image ImageType; - typedef itk::ImageSeriesReader ReaderType; + node.SetData( image ); + setlocale(LC_NUMERIC, previousCLocale); + std::cin.imbue(previousCppLocale); + } + catch (std::exception& e) + { + // reset locale then throw up + setlocale(LC_NUMERIC, previousCLocale); + std::cin.imbue(previousCppLocale); + + throw e; + } +} - DcmIoType::Pointer io = DcmIoType::New(); - typename ReaderType::Pointer reader = ReaderType::New(); +template +Image::Pointer DicomSeriesReader::LoadDICOMByITK( const StringContainer& filenames, CallbackCommand* command ) +{ + /******** Normal Case, 3D (also for GDCM < 2 usable) ***************/ + mitk::Image::Pointer image = mitk::Image::New(); - reader->SetImageIO(io); - reader->ReverseOrderOff(); + typedef itk::Image ImageType; + typedef itk::ImageSeriesReader ReaderType; - if (command) - { - reader->AddObserver(itk::ProgressEvent(), command); - } + DcmIoType::Pointer io = DcmIoType::New(); + typename ReaderType::Pointer reader = ReaderType::New(); - reader->SetFileNames(filenames); - reader->Update(); - image->InitializeByItk(reader->GetOutput()); - image->SetImportVolume(reader->GetOutput()->GetBufferPointer()); + reader->SetImageIO(io); + reader->ReverseOrderOff(); + + /* + if (command) + { + reader->AddObserver(itk::ProgressEvent(), command); + } + */ + + reader->SetFileNames(filenames); + reader->Update(); + image->InitializeByItk(reader->GetOutput()); + image->SetImportVolume(reader->GetOutput()->GetBufferPointer()); + + + MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " + << image->GetDimension(1) << ", " + << image->GetDimension(2) << "]"; #if (GDCM_MAJOR_VERSION == 2) && (GDCM_MINOR_VERSION < 1) && (GDCM_BUILD_VERSION < 15) - // workaround for a GDCM 2 bug until version 2.0.15: - // GDCM read spacing vector wrongly. Instead of "row spacing, column spacing", it misinterprets the DICOM tag as "column spacing, row spacing". - // this is undone here, until we use a GDCM that has this issue fixed. - // From the commit comments, GDCM 2.0.15 fixed the spacing interpretation with bug 2901181 - // http://sourceforge.net/tracker/index.php?func=detail&aid=2901181&group_id=137895&atid=739587 + // workaround for a GDCM 2 bug until version 2.0.15: + // GDCM read spacing vector wrongly. Instead of "row spacing, column spacing", it misinterprets the DICOM tag as "column spacing, row spacing". + // this is undone here, until we use a GDCM that has this issue fixed. + // From the commit comments, GDCM 2.0.15 fixed the spacing interpretation with bug 2901181 + // http://sourceforge.net/tracker/index.php?func=detail&aid=2901181&group_id=137895&atid=739587 - Vector3D correctedImageSpacing = image->GetGeometry()->GetSpacing(); - std::swap( correctedImageSpacing[0], correctedImageSpacing[1] ); - image->GetGeometry()->SetSpacing( correctedImageSpacing ); + Vector3D correctedImageSpacing = image->GetGeometry()->GetSpacing(); + std::swap( correctedImageSpacing[0], correctedImageSpacing[1] ); + image->GetGeometry()->SetSpacing( correctedImageSpacing ); #endif + MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " + << image->GetGeometry()->GetSpacing()[1] << ", " + << image->GetGeometry()->GetSpacing()[2] << "]"; - MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " << image->GetDimension(1) << ", " << image->GetDimension(2) << "]"; - MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " << image->GetGeometry()->GetSpacing()[1] << ", " << image->GetGeometry()->GetSpacing()[2] << "]"; + return image; +} + +std::list +DicomSeriesReader::SortIntoBlocksFor3DplusT( const StringContainer& presortedFilenames, bool sort, bool& canLoadAs4D ) +{ + std::list imageBlocks; - node.SetData(image); - setlocale(LC_NUMERIC, previousCLocale); - std::cin.imbue(previousCppLocale); + // sort only if requested (default) + StringContainer sorted_filenames = sort + ? DicomSeriesReader::SortSeriesSlices(presortedFilenames) + : presortedFilenames; + + gdcm::Tag ippTag(0x0020,0x0032); //Image position (Patient) + gdcm::Scanner scanner; + + scanner.AddTag(ippTag); + scanner.Scan(sorted_filenames); // make available image position for each file + std::string firstPosition; + unsigned int numberOfBlocks(0); // number of 3D image blocks + + // loop files to determine number of image blocks + for (StringContainer::const_iterator fileIter = sorted_filenames.begin(); + fileIter != sorted_filenames.end(); + ++fileIter) + { + std::string position = scanner.GetValue( fileIter->c_str(), ippTag); + MITK_DEBUG << "File " << *fileIter << " at " << position; + if (firstPosition.empty()) + { + firstPosition = position; + } + + if ( position == firstPosition ) + { + ++numberOfBlocks; + } + else + { + //break; // enough information to know the number of image blocks + } } - catch (std::exception& e) + + MITK_DEBUG << "Assuming " << numberOfBlocks << " image blocks"; + + if (numberOfBlocks == 0) return imageBlocks; // only possible if called with no files + + + // loop files to sort them into image blocks + unsigned int numberOfExpectedSlices(0); + for (unsigned int block = 0; block < numberOfBlocks; ++block) { - setlocale(LC_NUMERIC, previousCLocale); - std::cin.imbue(previousCppLocale); + StringContainer filesOfCurrentBlock; - throw e; + for ( StringContainer::const_iterator fileIter = sorted_filenames.begin() + block; + fileIter != sorted_filenames.end(); + //fileIter += numberOfBlocks) // TODO shouldn't this work? give invalid iterators on first attempts + ) + { + filesOfCurrentBlock.push_back( *fileIter ); + for (unsigned int b = 0; b < numberOfBlocks; ++b) + { + if (fileIter != sorted_filenames.end()) + ++fileIter; + } + } + + imageBlocks.push_back(filesOfCurrentBlock); + + if (block == 0) + { + numberOfExpectedSlices = filesOfCurrentBlock.size(); + } + else + { + if (filesOfCurrentBlock.size() != numberOfExpectedSlices) + { + MITK_WARN << "DicomSeriesReader expected " << numberOfBlocks + << " image blocks of " + << numberOfExpectedSlices + << " images each. Block " + << block + << " got " + << filesOfCurrentBlock.size() + << " instead. Cannot load this as 3D+t"; // TODO implement recovery (load as many slices 3D+t as much as possible) + canLoadAs4D = false; + } + } } + + return imageBlocks; } + } #endif -- 1.7.0.4