diff --git a/CMakeExternals/MITKData.cmake b/CMakeExternals/MITKData.cmake index 78af65a027..d586133631 100644 --- a/CMakeExternals/MITKData.cmake +++ b/CMakeExternals/MITKData.cmake @@ -1,38 +1,38 @@ #----------------------------------------------------------------------------- # MITK Data #----------------------------------------------------------------------------- # Sanity checks if(DEFINED MITK_DATA_DIR AND NOT EXISTS ${MITK_DATA_DIR}) message(FATAL_ERROR "MITK_DATA_DIR variable is defined but corresponds to non-existing directory") endif() set(proj MITK-Data) set(proj_DEPENDENCIES) set(MITK-Data_DEPENDS ${proj}) if(BUILD_TESTING) - set(revision_tag d01e488b) + set(revision_tag ab6c5213c778de18f1600975e8818541fd988803) #if(${proj}_REVISION_TAG) # set(revision_tag ${${proj}_REVISION_TAG}) #endif() ExternalProject_Add(${proj} URL ${MITK_THIRDPARTY_DOWNLOAD_PREFIX_URL}/MITK-Data_${revision_tag}.tar.gz UPDATE_COMMAND "" CONFIGURE_COMMAND "" BUILD_COMMAND "" INSTALL_COMMAND "" DEPENDS ${proj_DEPENDENCIES} ) set(MITK_DATA_DIR ${ep_source_dir}/${proj}) else() mitkMacroEmptyExternalProject(${proj} "${proj_DEPENDENCIES}") endif(BUILD_TESTING) diff --git a/Core/Code/DataManagement/mitkDataNodeFactory.cpp b/Core/Code/DataManagement/mitkDataNodeFactory.cpp index 86df353629..8a6c733691 100644 --- a/Core/Code/DataManagement/mitkDataNodeFactory.cpp +++ b/Core/Code/DataManagement/mitkDataNodeFactory.cpp @@ -1,482 +1,482 @@ /*=================================================================== 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 #include #include #include #include // C-Standard library includes #include #include // STL-related includes #include #include #include #include #include // VTK-related includes #include #include #include #include #include #include #include #include #include #include #include // ITK-related includes #include #include #include #include #include #include #include #include #include #include #include #ifdef NOMINMAX # define DEF_NOMINMAX # undef NOMINMAX #endif #include #ifdef DEF_NOMINMAX # ifndef NOMINMAX # define NOMINMAX # endif # undef DEF_NOMINMAX #endif #include #include // MITK-related includes #include "mitkSurface.h" #include "mitkPointSet.h" #include "mitkStringProperty.h" #include "mitkProperties.h" //#include "mitkMaterialProperty.h" #include "mitkLevelWindowProperty.h" #include "mitkVtkRepresentationProperty.h" #include "mitkVtkInterpolationProperty.h" #include "mitkVtkScalarModeProperty.h" #include "mitkImage.h" #include "mitkLookupTableProperty.h" #include "mitkLookupTable.h" #include "mitkImageChannelSelector.h" #include "mitkImageSliceSelector.h" #include "mitkCoreObjectFactory.h" #include "mitkTransferFunctionProperty.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkProgressBar.h" #include bool mitk::DataNodeFactory::m_TextureInterpolationActive = false; // default value for texture interpolation if nothing is defined in global options (see QmitkMainTemplate.ui.h) mitk::DataNodeFactory::DataNodeFactory() { m_Serie = false; m_OldProgress = 0; this->Modified(); //ensure that a CoreObjectFactory has been instantiated mitk::CoreObjectFactory::GetInstance(); } mitk::DataNodeFactory::~DataNodeFactory() {} void mitk::DataNodeFactory::SetImageSerie(bool serie) { m_Serie = serie; } void mitk::DataNodeFactory::GenerateData() { // IF filename is something.pic, and something.pic does not exist, try to read something.pic.gz // if there are both, something.pic and something.pic.gz, only the requested file is read // not only for images, but for all formats std::ifstream exists(m_FileName.c_str()); if (!exists) { std::string testfilename = m_FileName + ".gz"; std::ifstream exists(testfilename.c_str()); if (exists.good()) { m_FileName += ".gz"; } else { testfilename = m_FileName + ".GZ"; std::ifstream exists(testfilename.c_str()); if (exists.good()) { m_FileName += ".GZ"; } else { std::string message("File does not exist, or cannot be read. Filename = "); message += m_FileName; MITK_ERROR << message; itkExceptionMacro( << message ); } } } // part for DICOM // const char *numbers = "0123456789."; // std::string::size_type first_non_number; // first_non_number = itksys::SystemTools::GetFilenameName(m_FileName).find_first_not_of ( numbers ); if (DicomSeriesReader::IsDicom(this->m_FileName) /*|| first_non_number == std::string::npos*/) { this->ReadFileSeriesTypeDCM(); } else { bool usedNewDTNF = false; // the mitkBaseDataIO class returns a pointer of a vector of BaseData objects std::vector baseDataVector = mitk::BaseDataIO::LoadBaseDataFromFile( m_FileName, m_FilePrefix, m_FilePattern, m_Serie ); if( !baseDataVector.empty() ) this->ResizeOutputs((unsigned int)baseDataVector.size()); for(int i=0; i<(int)baseDataVector.size(); i++) { mitk::BaseData::Pointer baseData = baseDataVector.at(i); if( baseData.IsNotNull() ) { usedNewDTNF = true; mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(baseData); this->SetDefaultCommonProperties( node ); this->SetOutput(i, node); } } if(!usedNewDTNF && ( m_FileName != "" ) && !(m_Serie == false)) ReadFileSeriesTypeITKImageSeriesReader(); } } void mitk::DataNodeFactory::ResizeOutputs( const unsigned int& num ) { unsigned int prevNum = this->GetNumberOfOutputs(); this->SetNumberOfOutputs( num ); for ( unsigned int i = prevNum; i < num; ++i ) { this->SetNthOutput( i, this->MakeOutput( i ).GetPointer() ); } } bool mitk::DataNodeFactory::FileNameEndsWith( const std::string& name ) { if (m_FileName.size() < name.size()) return false; return m_FileName.substr(m_FileName.size() - name.size()) == name; } bool mitk::DataNodeFactory::FilePatternEndsWith( const std::string& name ) { return m_FilePattern.find( name ) != std::string::npos; } std::string mitk::DataNodeFactory::GetBaseFileName() { return itksys::SystemTools::GetFilenameName( m_FileName ); } std::string mitk::DataNodeFactory::GetBaseFilePrefix() { return itksys::SystemTools::GetFilenameName( m_FilePrefix ); } std::string mitk::DataNodeFactory::GetDirectory() { if ( !m_FileName.empty() ) return itksys::SystemTools::GetFilenamePath( m_FileName ); if ( !m_FilePrefix.empty() ) return itksys::SystemTools::GetFilenamePath( m_FilePrefix ); return std::string(); } void mitk::DataNodeFactory::ReadFileSeriesTypeDCM() { const char* previousCLocale = setlocale(LC_NUMERIC, NULL); setlocale(LC_NUMERIC, "C"); std::locale previousCppLocale( std::cin.getloc() ); std::locale l( "C" ); std::cin.imbue(l); if ( DicomSeriesReader::IsPhilips3DDicom(this->GetFileName()) ) { MITK_INFO << "it is a Philips3D US Dicom file" << std::endl; this->ResizeOutputs(1); DataNode::Pointer node = this->GetOutput(0); mitk::DicomSeriesReader::StringContainer stringvec; stringvec.push_back(this->GetFileName()); if (DicomSeriesReader::LoadDicomSeries(stringvec, *node)) { node->SetName(this->GetBaseFileName()); } setlocale(LC_NUMERIC, previousCLocale); std::cin.imbue(previousCppLocale); return; } - DicomSeriesReader::UidFileNamesMap names_map = DicomSeriesReader::GetSeries(this->GetDirectory(), this->m_SeriesRestrictions); + DicomSeriesReader::UidFileNamesMap names_map = DicomSeriesReader::GetSeries(this->GetDirectory(), true, this->m_SeriesRestrictions); // true = group gantry tilt images const unsigned int size = names_map.size(); this->ResizeOutputs(size); ProgressBar::GetInstance()->AddStepsToDo(size); ProgressBar::GetInstance()->Progress(); unsigned int outputIndex = 0u; const DicomSeriesReader::UidFileNamesMap::const_iterator n_end = names_map.end(); for (DicomSeriesReader::UidFileNamesMap::const_iterator n_it = names_map.begin(); n_it != n_end; ++n_it) { const std::string &uid = n_it->first; DataNode::Pointer node = this->GetOutput(outputIndex); MITK_INFO << "Reading series " << outputIndex << ": " << uid << std::endl; - if (DicomSeriesReader::LoadDicomSeries(n_it->second, *node)) + if (DicomSeriesReader::LoadDicomSeries(n_it->second, *node, true, true, true)) { std::string nodeName(uid); std::string studyDescription; if ( node->GetStringProperty( "dicom.study.StudyDescription", studyDescription ) ) { nodeName = studyDescription; std::string seriesDescription; if ( node->GetStringProperty( "dicom.series.SeriesDescription", seriesDescription ) ) { nodeName += "/" + seriesDescription; } } node->SetName(nodeName); ++outputIndex; } else { MITK_ERROR << "Skipping series " << outputIndex << " due to exception" << std::endl; } ProgressBar::GetInstance()->Progress(); } setlocale(LC_NUMERIC, previousCLocale); std::cin.imbue(previousCppLocale); } void mitk::DataNodeFactory::ReadFileSeriesTypeITKImageSeriesReader() { typedef itk::Image ImageType; typedef itk::ImageSeriesReader< ImageType > ReaderType; typedef itk::NumericSeriesFileNames NameGenerator; if ( ! this->GenerateFileList() ) { itkWarningMacro( "Sorry, file list could not be generated!" ); return ; } if ( m_MatchedFileNames.size() == 0 ) { itkWarningMacro( "Sorry, no files matched the given filename ("<< m_FileName <<")!" ); return ; } // // Finally, initialize the ITK-reader and load the files! // ReaderType::Pointer reader = ReaderType::New(); reader->SetFileNames( m_MatchedFileNames ); try { reader->Update(); ResizeOutputs( reader->GetNumberOfOutputs() ); for ( unsigned int i = 0; i < reader->GetNumberOfOutputs(); ++i ) { //Initialize mitk image from itk mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk( reader->GetOutput( i ) ); image->SetVolume( reader->GetOutput( i )->GetBufferPointer() ); //add the mitk image to the node mitk::DataNode::Pointer node = this->GetOutput( i ); node->SetData( image ); mitk::StringProperty::Pointer nameProp = mitk::StringProperty::New( m_FileName ); node->SetProperty( "name", nameProp ); } } catch ( const std::exception & e ) { itkWarningMacro( << e.what() ); return ; } } mitk::ColorProperty::Pointer mitk::DataNodeFactory::DefaultColorForOrgan( const std::string& organ ) { static bool initialized = false; static std::map< std::string, std::string > s_ColorMap; if (!initialized) { // all lowercase here, please! s_ColorMap.insert( std::make_pair( "ankle", "0xe38686") ); s_ColorMap.insert( std::make_pair( "appendix", "0xe38686") ); s_ColorMap.insert( std::make_pair( "blood vessels", "0xff3131") ); s_ColorMap.insert( std::make_pair( "bronchial tree", "0x3168ff") ); s_ColorMap.insert( std::make_pair( "bone", "0xd5d5d5") ); s_ColorMap.insert( std::make_pair( "brain", "0xff9cca") ); s_ColorMap.insert( std::make_pair( "coccyx", "0xe38686") ); s_ColorMap.insert( std::make_pair( "colon", "0xe38686") ); s_ColorMap.insert( std::make_pair( "cyst", "0xe38686") ); s_ColorMap.insert( std::make_pair( "elbow", "0xe38686") ); s_ColorMap.insert( std::make_pair( "eye", "0xe38686") ); s_ColorMap.insert( std::make_pair( "fallopian tube", "0xe38686") ); s_ColorMap.insert( std::make_pair( "fat", "0xff2bee") ); s_ColorMap.insert( std::make_pair( "hand", "0xe38686") ); s_ColorMap.insert( std::make_pair( "gall bladder", "0x567f18") ); s_ColorMap.insert( std::make_pair( "heart", "0xeb1d32") ); s_ColorMap.insert( std::make_pair( "hip", "0xe38686") ); s_ColorMap.insert( std::make_pair( "kidney", "0xd33f00") ); s_ColorMap.insert( std::make_pair( "knee", "0xe38686") ); s_ColorMap.insert( std::make_pair( "larynx", "0xe38686") ); s_ColorMap.insert( std::make_pair( "liver", "0xffcc3d") ); s_ColorMap.insert( std::make_pair( "lung", "0x6bdcff") ); s_ColorMap.insert( std::make_pair( "lymph node", "0xff0000") ); s_ColorMap.insert( std::make_pair( "muscle", "0xff456a") ); s_ColorMap.insert( std::make_pair( "nerve", "0xffea4f") ); s_ColorMap.insert( std::make_pair( "nose", "0xe38686") ); s_ColorMap.insert( std::make_pair( "oesophagus", "0xe38686") ); s_ColorMap.insert( std::make_pair( "ovaries", "0xe38686") ); s_ColorMap.insert( std::make_pair( "pancreas", "0xf9ab3d") ); s_ColorMap.insert( std::make_pair( "pelvis", "0xe38686") ); s_ColorMap.insert( std::make_pair( "penis", "0xe38686") ); s_ColorMap.insert( std::make_pair( "pharynx", "0xe38686") ); s_ColorMap.insert( std::make_pair( "prostate", "0xe38686") ); s_ColorMap.insert( std::make_pair( "rectum", "0xe38686") ); s_ColorMap.insert( std::make_pair( "sacrum", "0xe38686") ); s_ColorMap.insert( std::make_pair( "seminal vesicle", "0xe38686") ); s_ColorMap.insert( std::make_pair( "shoulder", "0xe38686") ); s_ColorMap.insert( std::make_pair( "spinal cord", "0xf5f93d") ); s_ColorMap.insert( std::make_pair( "spleen", "0xf96c3d") ); s_ColorMap.insert( std::make_pair( "stomach", "0xf96c3d") ); s_ColorMap.insert( std::make_pair( "teeth", "0xfffcd8") ); s_ColorMap.insert( std::make_pair( "testicles", "0xe38686") ); s_ColorMap.insert( std::make_pair( "thyroid", "0xfff694") ); s_ColorMap.insert( std::make_pair( "tongue", "0xe38686") ); s_ColorMap.insert( std::make_pair( "tumor", "0x937011") ); s_ColorMap.insert( std::make_pair( "urethra", "0xf8ff32") ); s_ColorMap.insert( std::make_pair( "urinary bladder", "0xf8ff32") ); s_ColorMap.insert( std::make_pair( "uterus", "0xe38686") ); s_ColorMap.insert( std::make_pair( "vagina", "0xe38686") ); s_ColorMap.insert( std::make_pair( "vertebra", "0xe38686") ); s_ColorMap.insert( std::make_pair( "wrist", "0xe38686") ); initialized = true; } std::string lowercaseOrgan(organ); for(unsigned int i = 0; i < organ.length(); i++) { lowercaseOrgan[i] = tolower(lowercaseOrgan[i]); } std::map< std::string, std::string >::iterator iter = s_ColorMap.find( lowercaseOrgan ); if ( iter != s_ColorMap.end() ) { std::string hexColor = iter->second; std::string hexRed = std::string("0x") + hexColor.substr( 2, 2 ); std::string hexGreen = std::string("0x") + hexColor.substr( 4, 2 ); std::string hexBlue = std::string("0x") + hexColor.substr( 6, 2 ); long int red = strtol( hexRed.c_str(), NULL, 16 ); long int green = strtol( hexGreen.c_str(), NULL, 16 ); long int blue = strtol( hexBlue.c_str(), NULL, 16 ); return ColorProperty::New( (float)red/ 255.0, (float)green/ 255.0, (float)blue/ 255.0 ); } else { // a default color (green) return ColorProperty::New( 0.0, 1.0, 0.0 ); } } void mitk::DataNodeFactory::SetDefaultCommonProperties(mitk::DataNode::Pointer &node) { // path mitk::StringProperty::Pointer pathProp = mitk::StringProperty::New( itksys::SystemTools::GetFilenamePath( m_FileName ) ); node->SetProperty( StringProperty::PATH, pathProp ); // name already defined? mitk::StringProperty::Pointer nameProp = dynamic_cast(node->GetProperty("name")); if(nameProp.IsNull() || (strcmp(nameProp->GetValue(),"No Name!")==0)) { // name already defined in BaseData mitk::StringProperty::Pointer baseDataNameProp = dynamic_cast(node->GetData()->GetProperty("name").GetPointer() ); if(baseDataNameProp.IsNull() || (strcmp(baseDataNameProp->GetValue(),"No Name!")==0)) { // name neither defined in node, nor in BaseData -> name = filename if (FileNameEndsWith( ".gz" )) m_FileName = m_FileName.substr( 0, m_FileName.length()-3 ); nameProp = mitk::StringProperty::New( itksys::SystemTools::GetFilenameWithoutLastExtension( m_FileName ) ); node->SetProperty( "name", nameProp ); } else { // name defined in BaseData! nameProp = mitk::StringProperty::New( baseDataNameProp->GetValue() ); node->SetProperty( "name", nameProp ); } } // visibility if(!node->GetProperty("visible")) node->SetVisibility(true); } diff --git a/Core/Code/IO/mitkDicomSeriesReader.cpp b/Core/Code/IO/mitkDicomSeriesReader.cpp index 18a1b9b43e..af84db1c53 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.cpp +++ b/Core/Code/IO/mitkDicomSeriesReader.cpp @@ -1,1131 +1,1331 @@ /*=================================================================== 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. ===================================================================*/ // uncomment for learning more about the internal sorting mechanisms //#define MBILOG_ENABLE_DEBUG #include #include #include #include #include #include #include #include "mitkProperties.h" namespace mitk { typedef itk::GDCMSeriesFileNames DcmFileNamesGeneratorType; + +DicomSeriesReader::SliceGroupingAnalysisResult::SliceGroupingAnalysisResult() +:m_GantryTilt(false) +{ +} + +DicomSeriesReader::StringContainer DicomSeriesReader::SliceGroupingAnalysisResult::GetBlockFilenames() +{ + return m_GroupedFiles; +} + +DicomSeriesReader::StringContainer DicomSeriesReader::SliceGroupingAnalysisResult::GetUnsortedFilenames() +{ + return m_UnsortedFiles; +} + +bool DicomSeriesReader::SliceGroupingAnalysisResult::ContainsGantryTilt() +{ + return m_GantryTilt; +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::AddFileToSortedBlock(const std::string& filename) +{ + m_GroupedFiles.push_back( filename ); +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::AddFileToUnsortedBlock(const std::string& filename) +{ + m_UnsortedFiles.push_back( filename ); +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::FlagGantryTilt() +{ + m_GantryTilt = true; +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::UndoPrematureGrouping() +{ + assert( !m_GroupedFiles.empty() ); + m_UnsortedFiles.insert( m_UnsortedFiles.begin(), m_GroupedFiles.back() ); + m_GroupedFiles.pop_back(); +} + + + + + + + const DicomSeriesReader::TagToPropertyMapType& DicomSeriesReader::GetDICOMTagsToMITKPropertyMap() { static bool initialized = false; static TagToPropertyMapType dictionary; if (!initialized) { /* Selection criteria: - no sequences because we cannot represent that - nothing animal related (specied, breed registration number), MITK focusses on human medical image processing. - only general attributes so far When extending this, we should make use of a real dictionary (GDCM/DCMTK and lookup the tag names there) */ // Patient module dictionary["0010|0010"] = "dicom.patient.PatientsName"; dictionary["0010|0020"] = "dicom.patient.PatientID"; dictionary["0010|0030"] = "dicom.patient.PatientsBirthDate"; dictionary["0010|0040"] = "dicom.patient.PatientsSex"; dictionary["0010|0032"] = "dicom.patient.PatientsBirthTime"; dictionary["0010|1000"] = "dicom.patient.OtherPatientIDs"; dictionary["0010|1001"] = "dicom.patient.OtherPatientNames"; dictionary["0010|2160"] = "dicom.patient.EthnicGroup"; dictionary["0010|4000"] = "dicom.patient.PatientComments"; dictionary["0012|0062"] = "dicom.patient.PatientIdentityRemoved"; dictionary["0012|0063"] = "dicom.patient.DeIdentificationMethod"; // General Study module dictionary["0020|000d"] = "dicom.study.StudyInstanceUID"; dictionary["0008|0020"] = "dicom.study.StudyDate"; dictionary["0008|0030"] = "dicom.study.StudyTime"; dictionary["0008|0090"] = "dicom.study.ReferringPhysiciansName"; dictionary["0020|0010"] = "dicom.study.StudyID"; dictionary["0008|0050"] = "dicom.study.AccessionNumber"; dictionary["0008|1030"] = "dicom.study.StudyDescription"; dictionary["0008|1048"] = "dicom.study.PhysiciansOfRecord"; dictionary["0008|1060"] = "dicom.study.NameOfPhysicianReadingStudy"; // General Series module dictionary["0008|0060"] = "dicom.series.Modality"; dictionary["0020|000e"] = "dicom.series.SeriesInstanceUID"; dictionary["0020|0011"] = "dicom.series.SeriesNumber"; dictionary["0020|0060"] = "dicom.series.Laterality"; dictionary["0008|0021"] = "dicom.series.SeriesDate"; dictionary["0008|0031"] = "dicom.series.SeriesTime"; dictionary["0008|1050"] = "dicom.series.PerformingPhysiciansName"; dictionary["0018|1030"] = "dicom.series.ProtocolName"; dictionary["0008|103e"] = "dicom.series.SeriesDescription"; dictionary["0008|1070"] = "dicom.series.OperatorsName"; dictionary["0018|0015"] = "dicom.series.BodyPartExamined"; dictionary["0018|5100"] = "dicom.series.PatientPosition"; dictionary["0028|0108"] = "dicom.series.SmallestPixelValueInSeries"; dictionary["0028|0109"] = "dicom.series.LargestPixelValueInSeries"; // VOI LUT module dictionary["0028|1050"] = "dicom.voilut.WindowCenter"; dictionary["0028|1051"] = "dicom.voilut.WindowWidth"; dictionary["0028|1055"] = "dicom.voilut.WindowCenterAndWidthExplanation"; initialized = true; } return dictionary; } DataNode::Pointer -DicomSeriesReader::LoadDicomSeries(const StringContainer &filenames, bool sort, bool check_4d, UpdateCallBackMethod callback) +DicomSeriesReader::LoadDicomSeries(const StringContainer &filenames, bool sort, bool check_4d, bool correctTilt, UpdateCallBackMethod callback) { DataNode::Pointer node = DataNode::New(); - if (DicomSeriesReader::LoadDicomSeries(filenames, *node, sort, check_4d, callback)) + if (DicomSeriesReader::LoadDicomSeries(filenames, *node, sort, check_4d, correctTilt, callback)) { if( filenames.empty() ) { return NULL; } return node; } else { return NULL; } } bool -DicomSeriesReader::LoadDicomSeries(const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, UpdateCallBackMethod callback) +DicomSeriesReader::LoadDicomSeries(const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, bool correctTilt, UpdateCallBackMethod callback) { if( filenames.empty() ) { MITK_WARN << "Calling LoadDicomSeries with empty filename string container. Probably invalid application logic."; node.SetData(NULL); return true; // this is not actually an error but the result is very simple } DcmIoType::Pointer io = DcmIoType::New(); try { if (io->CanReadFile(filenames.front().c_str())) { io->SetFileName(filenames.front().c_str()); io->ReadImageInformation(); switch (io->GetComponentType()) { case DcmIoType::UCHAR: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::CHAR: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::USHORT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::SHORT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::UINT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::INT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::ULONG: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::LONG: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::FLOAT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::DOUBLE: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); + DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; default: MITK_ERROR << "Found unsupported DICOM pixel type: (enum value) " << io->GetComponentType(); } if (node.GetData()) { return true; } } } catch(itk::MemoryAllocationError& e) { MITK_ERROR << "Out of memory. Cannot load DICOM series: " << e.what(); } catch(std::exception& e) { MITK_ERROR << "Error encountered when loading DICOM series:" << e.what(); } catch(...) { MITK_ERROR << "Unspecified error encountered when loading DICOM series."; } return false; } bool DicomSeriesReader::IsDicom(const std::string &filename) { DcmIoType::Pointer io = DcmIoType::New(); return io->CanReadFile(filename.c_str()); } bool DicomSeriesReader::IsPhilips3DDicom(const std::string &filename) { DcmIoType::Pointer io = DcmIoType::New(); if (io->CanReadFile(filename.c_str())) { //Look at header Tag 3001,0010 if it is "Philips3D" gdcm::Reader reader; reader.SetFileName(filename.c_str()); reader.Read(); gdcm::DataSet &data_set = reader.GetFile().GetDataSet(); gdcm::StringFilter sf; sf.SetFile(reader.GetFile()); if (data_set.FindDataElement(gdcm::Tag(0x3001, 0x0010)) && (sf.ToString(gdcm::Tag(0x3001, 0x0010)) == "Philips3D ")) { return true; } } return false; } bool DicomSeriesReader::ReadPhilips3DDicom(const std::string &filename, mitk::Image::Pointer output_image) { // Now get PhilipsSpecific Tags gdcm::PixmapReader reader; reader.SetFileName(filename.c_str()); reader.Read(); gdcm::DataSet &data_set = reader.GetFile().GetDataSet(); gdcm::StringFilter sf; sf.SetFile(reader.GetFile()); gdcm::Attribute<0x0028,0x0011> dimTagX; // coloumns || sagittal gdcm::Attribute<0x3001,0x1001, gdcm::VR::UL, gdcm::VM::VM1> dimTagZ; //I have no idea what is VM1. // (Philips specific) // transversal gdcm::Attribute<0x0028,0x0010> dimTagY; // rows || coronal gdcm::Attribute<0x0028,0x0008> dimTagT; // how many frames gdcm::Attribute<0x0018,0x602c> spaceTagX; // Spacing in X , unit is "physicalTagx" (usually centimeter) gdcm::Attribute<0x0018,0x602e> spaceTagY; gdcm::Attribute<0x3001,0x1003, gdcm::VR::FD, gdcm::VM::VM1> spaceTagZ; // (Philips specific) gdcm::Attribute<0x0018,0x6024> physicalTagX; // if 3, then spacing params are centimeter gdcm::Attribute<0x0018,0x6026> physicalTagY; gdcm::Attribute<0x3001,0x1002, gdcm::VR::US, gdcm::VM::VM1> physicalTagZ; // (Philips specific) dimTagX.Set(data_set); dimTagY.Set(data_set); dimTagZ.Set(data_set); dimTagT.Set(data_set); spaceTagX.Set(data_set); spaceTagY.Set(data_set); spaceTagZ.Set(data_set); physicalTagX.Set(data_set); physicalTagY.Set(data_set); physicalTagZ.Set(data_set); unsigned int dimX = dimTagX.GetValue(), dimY = dimTagY.GetValue(), dimZ = dimTagZ.GetValue(), dimT = dimTagT.GetValue(), physicalX = physicalTagX.GetValue(), physicalY = physicalTagY.GetValue(), physicalZ = physicalTagZ.GetValue(); float spaceX = spaceTagX.GetValue(), spaceY = spaceTagY.GetValue(), spaceZ = spaceTagZ.GetValue(); if (physicalX == 3) // spacing parameter in cm, have to convert it to mm. spaceX = spaceX * 10; if (physicalY == 3) // spacing parameter in cm, have to convert it to mm. spaceY = spaceY * 10; if (physicalZ == 3) // spacing parameter in cm, have to convert it to mm. spaceZ = spaceZ * 10; // Ok, got all necessary Tags! // Now read Pixeldata (7fe0,0010) X x Y x Z x T Elements const gdcm::Pixmap &pixels = reader.GetPixmap(); gdcm::RAWCodec codec; codec.SetPhotometricInterpretation(gdcm::PhotometricInterpretation::MONOCHROME2); codec.SetPixelFormat(pixels.GetPixelFormat()); codec.SetPlanarConfiguration(0); gdcm::DataElement out; codec.Decode(data_set.GetDataElement(gdcm::Tag(0x7fe0, 0x0010)), out); const gdcm::ByteValue *bv = out.GetByteValue(); const char *new_pixels = bv->GetPointer(); // Create MITK Image + Geometry typedef itk::Image ImageType; //Pixeltype might be different sometimes? Maybe read it out from header ImageType::RegionType myRegion; ImageType::SizeType mySize; ImageType::IndexType myIndex; ImageType::SpacingType mySpacing; ImageType::Pointer imageItk = ImageType::New(); mySpacing[0] = spaceX; mySpacing[1] = spaceY; mySpacing[2] = spaceZ; mySpacing[3] = 1; myIndex[0] = 0; myIndex[1] = 0; myIndex[2] = 0; myIndex[3] = 0; mySize[0] = dimX; mySize[1] = dimY; mySize[2] = dimZ; mySize[3] = dimT; myRegion.SetSize( mySize); myRegion.SetIndex( myIndex ); imageItk->SetSpacing(mySpacing); imageItk->SetRegions( myRegion); imageItk->Allocate(); imageItk->FillBuffer(0); itk::ImageRegionIterator iterator(imageItk, imageItk->GetLargestPossibleRegion()); iterator.GoToBegin(); unsigned long pixCount = 0; unsigned long planeSize = dimX*dimY; unsigned long planeCount = 0; unsigned long timeCount = 0; unsigned long numberOfSlices = dimZ; while (!iterator.IsAtEnd()) { unsigned long adressedPixel = pixCount + (numberOfSlices-1-planeCount)*planeSize // add offset to adress the first pixel of current plane + timeCount*numberOfSlices*planeSize; // add time offset iterator.Set( new_pixels[ adressedPixel ] ); pixCount++; ++iterator; if (pixCount == planeSize) { pixCount = 0; planeCount++; } if (planeCount == numberOfSlices) { planeCount = 0; timeCount++; } if (timeCount == dimT) { break; } } mitk::CastToMitkImage(imageItk, output_image); return true; // actually never returns false yet.. but exception possible } + +DicomSeriesReader::GantryTiltInformation::GantryTiltInformation() +: m_ShiftUp(0.0) +, m_ShiftRight(0.0) +, m_ShiftNormal(0.0) +, m_NumberOfSlicesApart(1) +{ +} + +DicomSeriesReader::GantryTiltInformation::GantryTiltInformation( + const Point3D& origin1, const Point3D& origin2, + const Vector3D& right, const Vector3D& up, + unsigned int numberOfSlicesApart) +: m_ShiftUp(0.0) +, m_ShiftRight(0.0) +, m_ShiftNormal(0.0) +, m_NumberOfSlicesApart(numberOfSlicesApart) +{ + assert(numberOfSlicesApart); + // 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 -DicomSeriesReader::TwoStringContainers + squared distance = | (pointAlongNormal - origin2) x (origin2 - origin1) | ^ 2 + / + |pointAlongNormal - origin2| ^ 2 + + ( x meaning the cross product ) + */ + + Vector3D normal = itk::CrossProduct(right, up); + Point3D 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 right: " << right; + MITK_DEBUG << " v up: " << up; + MITK_DEBUG << " v normal: " << normal; + + Point3D projectionRight = projectPointOnLine( origin1, origin2, right ); + Point3D projectionUp = projectPointOnLine( origin1, origin2, up ); + Point3D projectionNormal = projectPointOnLine( origin1, origin2, normal ); + + m_ShiftRight = (projectionRight - origin2).GetNorm(); + m_ShiftUp = (projectionUp - origin2).GetNorm(); + m_ShiftNormal = (projectionNormal - origin2).GetNorm(); + + MITK_DEBUG << " shift normal: " << m_ShiftNormal; + MITK_DEBUG << " shift up: " << m_ShiftUp; + MITK_DEBUG << " shift right: " << m_ShiftRight; + + MITK_DEBUG << " tilt angle (rad): " << tanh( m_ShiftUp / m_ShiftNormal ); + MITK_DEBUG << " tilt angle (deg): " << tanh( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535; + } +} + +Point3D +DicomSeriesReader::GantryTiltInformation::projectPointOnLine( Point3D p, Point3D lineOrigin, Vector3D lineDirection ) +{ + /** + See illustration at http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage472/ + + vector(lineOrigin,p) = normal * ( innerproduct((p - lineOrigin),normal) / squared-length(normal) ) + */ + + Vector3D lineOriginToP = p - lineOrigin; + ScalarType innerProduct = lineOriginToP * lineDirection; + + ScalarType factor = innerProduct / lineDirection.GetSquaredNorm(); + Point3D projection = lineOrigin + factor * lineDirection; + + return projection; +} + +ScalarType +DicomSeriesReader::GantryTiltInformation::GetTiltCorrectedAdditionalSize() const +{ + // this seems to be a bit too much sometimes, but better too much than cutting off parts of the image + return int(m_ShiftUp + 1.0); // to next bigger int: plus 1, then cut off after point +} + +ScalarType +DicomSeriesReader::GantryTiltInformation::GetMatrixCoefficientForCorrectionInWorldCoordinates() const +{ + // so many mm shifted per slice! + return m_ShiftUp / static_cast(m_NumberOfSlicesApart); +} + +ScalarType +DicomSeriesReader::GantryTiltInformation::GetRealZSpacing() const +{ + return m_ShiftNormal / static_cast(m_NumberOfSlicesApart); +} + + +bool +DicomSeriesReader::GantryTiltInformation::IsSheared() const +{ + return ( m_ShiftRight > 0.001 + || m_ShiftUp > 0.001); +} + + +bool +DicomSeriesReader::GantryTiltInformation::IsRegularGantryTilt() const +{ + return ( m_ShiftRight < 0.001 + && m_ShiftUp > 0.001); +} + + +std::string +DicomSeriesReader::ConstCharStarToString(const char* s) +{ + return s ? std::string(s) : std::string(); +} + +Point3D +DicomSeriesReader::DICOMStringToPoint3D(const std::string& s, bool& successful) +{ + Point3D p; + successful = true; + + std::istringstream originReader(s); + std::string coordinate; + unsigned int dim(0); + while( std::getline( originReader, coordinate, '\\' ) ) + { + if (dim > 4) break; // otherwise we access invalid array index + + p[dim++] = atof(coordinate.c_str()); + } + + if (dim != 3) + { + successful = false; + MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0032). Found " << dim << " instead of 3 values."; + } + + return p; +} + +void +DicomSeriesReader::DICOMStringToOrientationVectors(const std::string& s, Vector3D& right, Vector3D& up, bool& successful) +{ + successful = true; + + std::istringstream orientationReader(s); + std::string coordinate; + unsigned int dim(0); + while( std::getline( orientationReader, coordinate, '\\' ) ) + { + if (dim > 6) break; // otherwise we access invalid array index + + if (dim<3) + { + right[dim++] = atof(coordinate.c_str()); + } + else + { + up[dim++ - 3] = atof(coordinate.c_str()); + } + } + + if (dim != 6) + { + successful = false; + MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0037). Found " << dim << " instead of 6 values."; + } +} + + +DicomSeriesReader::SliceGroupingAnalysisResult DicomSeriesReader::AnalyzeFileForITKImageSeriesReaderSpacingAssumption( const StringContainer& files, + bool groupImagesWithGantryTilt, const gdcm::Scanner::MappingType& tagValueMappings_) { // result.first = files that fit ITK's assumption // result.second = files that do not fit, should be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption() again - TwoStringContainers result; + SliceGroupingAnalysisResult result; // we const_cast here, because I could not use a map.at(), which would make the code much more readable gdcm::Scanner::MappingType& tagValueMappings = const_cast(tagValueMappings_); const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation Vector3D fromFirstToSecondOrigin; fromFirstToSecondOrigin.Fill(0.0); bool fromFirstToSecondOriginInitialized(false); Point3D thisOrigin; Point3D lastOrigin; lastOrigin.Fill(0.0f); Point3D lastDifferentOrigin; lastDifferentOrigin.Fill(0.0f); bool lastOriginInitialized(false); MITK_DEBUG << "--------------------------------------------------------------------------------"; - MITK_DEBUG << "Analyzing files for z-spacing assumption of ITK's ImageSeriesReader "; + MITK_DEBUG << "Analyzing files for z-spacing assumption of ITK's ImageSeriesReader (group tilted: " << groupImagesWithGantryTilt << ")"; unsigned int fileIndex(0); for (StringContainer::const_iterator fileIter = files.begin(); fileIter != files.end(); ++fileIter, ++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 - const char* value = tagValueMappings[fileIter->c_str()][tagImagePositionPatient]; - if (value) - { - thisOriginString = value; - } - - std::istringstream originReader(thisOriginString); - std::string coordinate; - unsigned int dim(0); - while( std::getline( originReader, coordinate, '\\' ) ) thisOrigin[dim++] = atof(coordinate.c_str()); + thisOriginString = ConstCharStarToString( tagValueMappings[fileIter->c_str()][tagImagePositionPatient] ); - if (dim != 3) - { - MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0032). Found " << dim << "instead of 3 values."; - } + bool ignoredConversionError(-42); // hard to get here, no graceful way to react + thisOrigin = DICOMStringToPoint3D( thisOriginString, ignoredConversionError ); MITK_DEBUG << " " << fileIndex << " " << *fileIter << " at " << thisOriginString << "(" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")"; if ( lastOriginInitialized && (thisOrigin == lastOrigin) ) { MITK_DEBUG << " ==> Sort away " << *fileIter << " for separate time step"; // we already have one occupying this position - result.second.push_back( *fileIter ); + result.AddFileToUnsortedBlock( *fileIter ); fileFitsIntoPattern = false; } else { if (!fromFirstToSecondOriginInitialized && lastOriginInitialized) // calculate vector as soon as possible when we get a new position { fromFirstToSecondOrigin = thisOrigin - lastDifferentOrigin; fromFirstToSecondOriginInitialized = true; - // Now make sure this direction is along the normal vector of the first slice + // 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 loaded into a single mitk::Image at the moment + // 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. - // Again ugly code to read tag Image Orientation into two vEctors 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 thisOrientationString; - const char* value = tagValueMappings[fileIter->c_str()][tagImageOrientation]; - if (value) - { - thisOrientationString = value; - } - std::istringstream orientationReader(thisOrientationString); - std::string coordinate; - unsigned int dim(0); - while( std::getline( orientationReader, coordinate, '\\' ) ) - if (dim<3) right[dim++] = atof(coordinate.c_str()); - else up[dim++ - 3] = atof(coordinate.c_str()); - - if (dim != 6) - { - MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0037). Found " << dim << "instead of 6 values."; - } - - - /* - Determine if line (thisOrigin + l * normal) contains lastDifferentOrigin. - - Done by calculating the distance of lastDifferentOrigin from line (thisOrigin + l *normal) + DICOMStringToOrientationVectors( tagValueMappings[fileIter->c_str()][tagImageOrientation], right, up, ignoredConversionError ); + + GantryTiltInformation tiltInfo( lastDifferentOrigin, thisOrigin, right, up, 1 ); - E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html - - squared distance = | (pointAlongNormal - thisOrign) x (thisOrigin - lastDifferentOrigin) | ^ 2 - / - |pointAlongNormal - thisOrigin| ^ 2 - - ( x meaning the cross product ) - - MITK_DEBUG << "Tilt check: right vector (" << right[0] << "," << right[1] << "," << right[2] << "), " - "up vector (" << up[0] << "," << up[1] << "," << up[2] << ")"; - */ - - Vector3D normal = itk::CrossProduct(right, up); - Point3D pointAlongNormal = thisOrigin + normal; - - double numerator = itk::CrossProduct( pointAlongNormal - thisOrigin , thisOrigin - lastDifferentOrigin ).GetSquaredNorm(); - double denominator = (pointAlongNormal - thisOrigin).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 + if ( tiltInfo.IsSheared() ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt { - MITK_DEBUG << " Series might contain a tilted geometry"; - MITK_DEBUG << " Distance of expected slice origin from actual slice origin: " << distance; - MITK_DEBUG << " ==> Sort away " << *fileIter << " for later analysis"; - - /* Pessimistic approach: split block right here - - result.first.assign( files.begin(), fileIter ); - result.second.insert( result.second.end(), fileIter, files.end() ); + /* 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 (*fileIter) - return result; // stop processing with first split - */ + // tell apart gantry tilt from overall skewedness + // sort out irregularly sheared slices, that IS NOT tilting - /* optimistic approach: save file for later, check all further files */ - - result.second.push_back(*fileIter); - fileFitsIntoPattern = false; + if ( groupImagesWithGantryTilt && tiltInfo.IsRegularGantryTilt() ) + { + result.FlagGantryTilt(); + result.AddFileToSortedBlock(*fileIter); // this file is good for current block + fileFitsIntoPattern = true; + } + else + { + result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis + fileFitsIntoPattern = false; + } } else { - result.first.push_back(*fileIter); // this file is good for current block + result.AddFileToSortedBlock(*fileIter); // 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(); double toleratedError(0.005); // max. 1/10mm error when measurement crosses 20 slices in z direction if (norm > toleratedError) { MITK_DEBUG << " File does not fit into the inter-slice distance pattern (diff = " << norm << ", allowed " << toleratedError << ")."; MITK_DEBUG << " Expected position (" << assumedOrigin[0] << "," << assumedOrigin[1] << "," << assumedOrigin[2] << "), got position (" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")"; MITK_DEBUG << " ==> Sort away " << *fileIter << " 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 - /* - Pessimistic approach: split right here: - - result.first.assign( files.begin(), fileIter ); - result.second.insert( result.second.end(), fileIter, files.end() ); - - return result; // stop processing with first split - */ - /* Optimistic approach: check if any of the remaining slices fits in */ - result.second.push_back( *fileIter ); // sort away for further analysis + result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis fileFitsIntoPattern = false; } else { - result.first.push_back(*fileIter); // this file is good for current block + result.AddFileToSortedBlock(*fileIter); // this file is good for current block fileFitsIntoPattern = true; } } else // this should be the very first slice { - result.first.push_back(*fileIter); // this file is good for current block + result.AddFileToSortedBlock(*fileIter); // this file is good for current block fileFitsIntoPattern = true; } } - // recored current origin for reference in later iterations + // 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.GetBlockFilenames().size() == 2 ) + { + result.UndoPrematureGrouping(); + } + } return result; } DicomSeriesReader::UidFileNamesMap -DicomSeriesReader::GetSeries(const StringContainer& files, const StringContainer &restrictions) +DicomSeriesReader::GetSeries(const StringContainer& files, bool groupImagesWithGantryTilt, const StringContainer &restrictions) { - return GetSeries(files, true, restrictions); + return GetSeries(files, true, groupImagesWithGantryTilt, restrictions); } DicomSeriesReader::UidFileNamesMap -DicomSeriesReader::GetSeries(const StringContainer& files, bool sortTo3DPlust, const StringContainer& /*restrictions*/) +DicomSeriesReader::GetSeries(const StringContainer& files, bool sortTo3DPlust, bool groupImagesWithGantryTilt, const StringContainer& /*restrictions*/) { /** assumption about this method: returns a map of uid-like-key --> list(filename) each entry should contain filenames that have images of same - series instance uid (automatically done by GDCMSeriesFileNames - 0020,0037 image orientation (patient) - 0028,0030 pixel spacing (x,y) - 0018,0050 slice thickness */ UidFileNamesMap groupsOfSimilarImages; // preliminary result, refined into the final result mapOf3DPlusTBlocks // use GDCM directly, itk::GDCMSeriesFileNames does not work with GDCM 2 // PART I: scan files for sorting relevant DICOM tags, // separate images that differ in any of those // attributes (they cannot possibly form a 3D block) // scan 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 ); // additional tags read in this scan to allow later analysis // THESE tag are not used for initial separating of files const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) scanner.AddTag( tagImagePositionPatient ); // TODO add further restrictions from arguments // let GDCM scan files if ( !scanner.Scan( files ) ) { MITK_ERROR << "gdcm::Scanner failed when scanning " << files.size() << " input files."; return groupsOfSimilarImages; } // 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 << "Scan file " << fileIter->first << std::endl; if ( std::string(fileIter->first).empty() ) continue; // TODO understand why Scanner has empty string entries // we const_cast here, because I could not use a map.at() function in CreateMoreUniqueSeriesIdentifier. // doing the same thing with find would make the code less readable. Since we forget the Scanner results // anyway after this function, we can simply tolerate empty map entries introduced by bad operator[] access std::string moreUniqueSeriesId = CreateMoreUniqueSeriesIdentifier( const_cast(fileIter->second) ); groupsOfSimilarImages [ moreUniqueSeriesId ].push_back( fileIter->first ); } - // PART III: sort slices spatially + // PART II: sort slices spatially for ( UidFileNamesMap::const_iterator groupIter = groupsOfSimilarImages.begin(); groupIter != groupsOfSimilarImages.end(); ++groupIter ) { try { groupsOfSimilarImages[ groupIter->first ] = SortSeriesSlices( groupIter->second ); // sort each slice group spatially } catch(...) { MITK_ERROR << "Catched something."; } } - // PART II: analyze pre-sorted images for valid blocks (i.e. blocks of equal z-spacing), + // PART III: analyze pre-sorted images for valid blocks (i.e. blocks of equal z-spacing), // separate into multiple blocks if necessary. // // Analysis performs the following steps: // * imitate itk::ImageSeriesReader: use the distance between the first two images as z-spacing // * check what images actually fulfill ITK's z-spacing assumption // * separate all images that fail the test into new blocks, re-iterate analysis for these blocks UidFileNamesMap mapOf3DPlusTBlocks; // final result of this function for ( UidFileNamesMap::const_iterator groupIter = groupsOfSimilarImages.begin(); groupIter != groupsOfSimilarImages.end(); ++groupIter ) { UidFileNamesMap mapOf3DBlocks; // intermediate result for only this group(!) + std::map mapOf3DBlockAnalysisResults; StringContainer filesStillToAnalyze = groupIter->second; std::string groupUID = groupIter->first; unsigned int subgroup(0); MITK_DEBUG << "Analyze group " << groupUID; while (!filesStillToAnalyze.empty()) // repeat until all files are grouped somehow { - TwoStringContainers analysisResult = AnalyzeFileForITKImageSeriesReaderSpacingAssumption( filesStillToAnalyze, scanner.GetMappings() ); + SliceGroupingAnalysisResult analysisResult = + AnalyzeFileForITKImageSeriesReaderSpacingAssumption( filesStillToAnalyze, + groupImagesWithGantryTilt, + scanner.GetMappings() ); // enhance the UID for additional groups std::stringstream newGroupUID; newGroupUID << groupUID << '.' << subgroup; - mapOf3DBlocks[ newGroupUID.str() ] = analysisResult.first; + mapOf3DBlocks[ newGroupUID.str() ] = analysisResult.GetBlockFilenames(); MITK_DEBUG << "Result: sorted 3D group " << newGroupUID.str() << " with " << mapOf3DBlocks[ newGroupUID.str() ].size() << " files"; ++subgroup; - filesStillToAnalyze = analysisResult.second; // remember what needs further analysis + filesStillToAnalyze = analysisResult.GetUnsortedFilenames(); // remember what needs further analysis } // end of grouping, now post-process groups // PART IV: attempt to group blocks to 3D+t blocks if requested // inspect entries of mapOf3DBlocks // - if number of files is identical to previous entry, collect for 3D+t block // - as soon as number of files changes from previous entry, record collected blocks as 3D+t block, start a new one, continue // decide whether or not to group 3D blocks into 3D+t blocks where possible if ( !sortTo3DPlust ) { // copy 3D blocks to output // TODO avoid collisions (or prove impossibility) mapOf3DPlusTBlocks.insert( mapOf3DBlocks.begin(), mapOf3DBlocks.end() ); } else { // sort 3D+t (as described in "PART IV") MITK_DEBUG << "================================================================================"; MITK_DEBUG << "3D+t analysis:"; unsigned int numberOfFilesInPreviousBlock(0); std::string previousBlockKey; for ( UidFileNamesMap::const_iterator block3DIter = mapOf3DBlocks.begin(); block3DIter != mapOf3DBlocks.end(); ++block3DIter ) { unsigned int numberOfFilesInThisBlock = block3DIter->second.size(); std::string thisBlockKey = block3DIter->first; if (numberOfFilesInPreviousBlock == 0) { numberOfFilesInPreviousBlock = numberOfFilesInThisBlock; mapOf3DPlusTBlocks[thisBlockKey].insert( mapOf3DPlusTBlocks[thisBlockKey].end(), block3DIter->second.begin(), block3DIter->second.end() ); MITK_DEBUG << " 3D+t group " << thisBlockKey << " started"; previousBlockKey = thisBlockKey; } else { bool identicalOrigins; try { // check whether this and the previous block share a comon origin // TODO should be safe, but a little try/catch or other error handling wouldn't hurt const char *origin_value = scanner.GetValue( mapOf3DBlocks[thisBlockKey].front().c_str(), tagImagePositionPatient ), *previous_origin_value = scanner.GetValue( mapOf3DBlocks[previousBlockKey].front().c_str(), tagImagePositionPatient ), *destination_value = scanner.GetValue( mapOf3DBlocks[thisBlockKey].back().c_str(), tagImagePositionPatient ), *previous_destination_value = scanner.GetValue( mapOf3DBlocks[previousBlockKey].back().c_str(), tagImagePositionPatient ); if (!origin_value || !previous_origin_value || !destination_value || !previous_destination_value) { identicalOrigins = false; } else { - std::string thisOriginString = origin_value; - std::string previousOriginString = previous_origin_value; + std::string thisOriginString = ConstCharStarToString( origin_value ); + std::string previousOriginString = ConstCharStarToString( previous_origin_value ); // also compare last origin, because this might differ if z-spacing is different - std::string thisDestinationString = destination_value; - std::string previousDestinationString = previous_destination_value; + std::string thisDestinationString = ConstCharStarToString( destination_value ); + std::string previousDestinationString = ConstCharStarToString( previous_destination_value ); identicalOrigins = ( (thisOriginString == previousOriginString) && (thisDestinationString == previousDestinationString) ); } } catch(...) { identicalOrigins = false; } if (identicalOrigins && (numberOfFilesInPreviousBlock == numberOfFilesInThisBlock)) { // group with previous block mapOf3DPlusTBlocks[previousBlockKey].insert( mapOf3DPlusTBlocks[previousBlockKey].end(), block3DIter->second.begin(), block3DIter->second.end() ); MITK_DEBUG << " --> group enhanced with another timestep"; } else { // start a new block mapOf3DPlusTBlocks[thisBlockKey].insert( mapOf3DPlusTBlocks[thisBlockKey].end(), block3DIter->second.begin(), block3DIter->second.end() ); MITK_DEBUG << " ==> group closed with " << mapOf3DPlusTBlocks[previousBlockKey].size() / numberOfFilesInPreviousBlock << " time steps"; previousBlockKey = thisBlockKey; MITK_DEBUG << " 3D+t group " << thisBlockKey << " started"; } } numberOfFilesInPreviousBlock = numberOfFilesInThisBlock; } } } MITK_DEBUG << "================================================================================"; MITK_DEBUG << "Summary: "; for ( UidFileNamesMap::const_iterator groupIter = mapOf3DPlusTBlocks.begin(); groupIter != mapOf3DPlusTBlocks.end(); ++groupIter ) { MITK_DEBUG << " Image volume " << groupIter->first << " with " << groupIter->second.size() << " files"; } MITK_DEBUG << "Done. "; MITK_DEBUG << "================================================================================"; return mapOf3DPlusTBlocks; } DicomSeriesReader::UidFileNamesMap -DicomSeriesReader::GetSeries(const std::string &dir, const StringContainer &restrictions) +DicomSeriesReader::GetSeries(const std::string &dir, bool groupImagesWithGantryTilt, const StringContainer &restrictions) { gdcm::Directory directoryLister; directoryLister.Load( dir.c_str(), false ); // non-recursive - return GetSeries(directoryLister.GetFilenames(), restrictions); + return GetSeries(directoryLister.GetFilenames(), groupImagesWithGantryTilt, restrictions); } std::string DicomSeriesReader::CreateSeriesIdentifierPart( gdcm::Scanner::TagToValue& tagValueMap, const gdcm::Tag& tag ) { std::string result; try { result = IDifyTagValue( tagValueMap[ tag ] ? tagValueMap[ tag ] : std::string("") ); } catch (std::exception& e) { MITK_WARN << "Could not access tag " << tag << ": " << e.what(); } return result; } std::string DicomSeriesReader::CreateMoreUniqueSeriesIdentifier( 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[ tagSeriesInstanceUID ]; } catch (std::exception& e) { MITK_ERROR << "CreateMoreUniqueSeriesIdentifier() could not access series instance UID. Something is seriously wrong with this image."; MITK_ERROR << "Error from exception: " << e.what(); } constructedID += CreateSeriesIdentifierPart( tagValueMap, tagNumberOfRows ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagNumberOfColumns ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagPixelSpacing ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagSliceThickness ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagImageOrientation ); constructedID.resize( constructedID.length() - 1 ); // cut of trailing '.' return constructedID; } std::string DicomSeriesReader::IDifyTagValue(const std::string& value) { std::string IDifiedValue( value ); if (value.empty()) throw std::logic_error("IDifyTagValue() illegaly called with empty tag 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; } DicomSeriesReader::StringContainer -DicomSeriesReader::GetSeries(const std::string &dir, const std::string &series_uid, const StringContainer &restrictions) +DicomSeriesReader::GetSeries(const std::string &dir, const std::string &series_uid, bool groupImagesWithGantryTilt, const StringContainer &restrictions) { - UidFileNamesMap allSeries = GetSeries(dir, restrictions); + UidFileNamesMap allSeries = GetSeries(dir, groupImagesWithGantryTilt, restrictions); StringContainer resultingFileList; for ( UidFileNamesMap::const_iterator idIter = allSeries.begin(); idIter != allSeries.end(); ++idIter ) { if ( idIter->first.find( series_uid ) == 0 ) // this ID starts with given series_uid { resultingFileList.insert( resultingFileList.end(), idIter->second.begin(), idIter->second.end() ); // append } } return resultingFileList; } DicomSeriesReader::StringContainer DicomSeriesReader::SortSeriesSlices(const StringContainer &unsortedFilenames) { gdcm::Sorter sorter; sorter.SetSortFunction(DicomSeriesReader::GdcmSortFunction); try { sorter.Sort(unsortedFilenames); return sorter.GetFilenames(); } catch(std::logic_error&) { MITK_WARN << "Sorting error. Leaving series unsorted."; return unsortedFilenames; } } bool DicomSeriesReader::GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::DataSet &ds2) { // make sure we habe Image Position and Orientation if ( ! ( ds1.FindDataElement(gdcm::Tag(0x0020,0x0032)) && ds1.FindDataElement(gdcm::Tag(0x0020,0x0037)) && ds2.FindDataElement(gdcm::Tag(0x0020,0x0032)) && ds2.FindDataElement(gdcm::Tag(0x0020,0x0037)) ) ) { MITK_WARN << "Dicom images are missing attributes for a meaningful sorting."; throw std::logic_error("Dicom images are missing attributes for a meaningful sorting."); } gdcm::Attribute<0x0020,0x0032> image_pos1; // Image Position (Patient) gdcm::Attribute<0x0020,0x0037> image_orientation1; // Image Orientation (Patient) image_pos1.Set(ds1); image_orientation1.Set(ds1); gdcm::Attribute<0x0020,0x0032> image_pos2; gdcm::Attribute<0x0020,0x0037> image_orientation2; image_pos2.Set(ds2); image_orientation2.Set(ds2); if (image_orientation1 != image_orientation2) { MITK_ERROR << "Dicom images have different orientations."; throw std::logic_error("Dicom images have different orientations. Call GetSeries() first to separate images."); } 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]; 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]; } if ( fabs(dist1 - dist2) < mitk::eps) { gdcm::Attribute<0x0008,0x0032> acq_time1; // Acquisition time (may be missing, so we check existence first) gdcm::Attribute<0x0008,0x0032> acq_time2; gdcm::Attribute<0x0020,0x0012> acq_number1; // Acquisition number (may also be missing, so we check existence first) gdcm::Attribute<0x0020,0x0012> acq_number2; if (ds1.FindDataElement(gdcm::Tag(0x0008,0x0032)) && ds2.FindDataElement(gdcm::Tag(0x0008,0x0032))) { acq_time1.Set(ds1); acq_time2.Set(ds2); return acq_time1 < acq_time2; } else if (ds1.FindDataElement(gdcm::Tag(0x0020,0x0012)) && ds2.FindDataElement(gdcm::Tag(0x0020,0x0012))) { acq_number1.Set(ds1); acq_number2.Set(ds2); return acq_number1 < acq_number2; } else { return true; } } else { // default: compare position return dist1 < dist2; } } std::string DicomSeriesReader::GetConfigurationString() { std::stringstream configuration; configuration << "MITK_USE_GDCMIO: "; configuration << "true"; configuration << "\n"; configuration << "GDCM_VERSION: "; #ifdef GDCM_MAJOR_VERSION configuration << GDCM_VERSION; #endif //configuration << "\n"; return configuration.str(); } void DicomSeriesReader::CopyMetaDataToImageProperties(StringContainer filenames, const gdcm::Scanner::MappingType &tagValueMappings_, DcmIoType *io, Image *image) { std::list imageBlock; imageBlock.push_back(filenames); CopyMetaDataToImageProperties(imageBlock, tagValueMappings_, io, image); } void DicomSeriesReader::CopyMetaDataToImageProperties( std::list imageBlock, const gdcm::Scanner::MappingType& tagValueMappings_, DcmIoType* io, Image* image) { if (!io || !image) return; StringLookupTable filesForSlices; StringLookupTable sliceLocationForSlices; StringLookupTable instanceNumberForSlices; StringLookupTable SOPInstanceNumberForSlices; gdcm::Scanner::MappingType& tagValueMappings = const_cast(tagValueMappings_); //DICOM tags which should be added to the image properties const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number unsigned int timeStep(0); std::string propertyKeySliceLocation = "dicom.image.0020.1041"; std::string propertyKeyInstanceNumber = "dicom.image.0020.0013"; std::string propertyKeySOPInstanceNumber = "dicom.image.0008.0018"; // tags for each image for ( std::list::iterator i = imageBlock.begin(); i != imageBlock.end(); i++, timeStep++ ) { const StringContainer& files = (*i); unsigned int slice(0); for ( StringContainer::const_iterator fIter = files.begin(); fIter != files.end(); ++fIter, ++slice ) { filesForSlices.SetTableValue( slice, *fIter ); gdcm::Scanner::TagToValue tagValueMapForFile = tagValueMappings[fIter->c_str()]; if(tagValueMapForFile.find(tagSliceLocation) != tagValueMapForFile.end()) sliceLocationForSlices.SetTableValue(slice, tagValueMapForFile[tagSliceLocation]); if(tagValueMapForFile.find(tagInstanceNumber) != tagValueMapForFile.end()) instanceNumberForSlices.SetTableValue(slice, tagValueMapForFile[tagInstanceNumber]); if(tagValueMapForFile.find(tagSOPInstanceNumber) != tagValueMapForFile.end()) SOPInstanceNumberForSlices.SetTableValue(slice, tagValueMapForFile[tagSOPInstanceNumber]); } image->SetProperty( "files", StringLookupTableProperty::New( filesForSlices ) ); //If more than one time step add postfix ".t" + timestep if(timeStep != 0) { propertyKeySliceLocation.append(".t" + timeStep); propertyKeyInstanceNumber.append(".t" + timeStep); propertyKeySOPInstanceNumber.append(".t" + timeStep); } image->SetProperty( propertyKeySliceLocation.c_str(), StringLookupTableProperty::New( sliceLocationForSlices ) ); image->SetProperty( propertyKeyInstanceNumber.c_str(), StringLookupTableProperty::New( instanceNumberForSlices ) ); image->SetProperty( propertyKeySOPInstanceNumber.c_str(), StringLookupTableProperty::New( SOPInstanceNumberForSlices ) ); } // Copy tags for series, study, patient level (leave interpretation to application). // These properties will be copied to the DataNode by DicomSeriesReader. // tags for the series (we just use the one that ITK copied to its dictionary (proably that of the last slice) const itk::MetaDataDictionary& dict = io->GetMetaDataDictionary(); const TagToPropertyMapType& propertyLookup = DicomSeriesReader::GetDICOMTagsToMITKPropertyMap(); itk::MetaDataDictionary::ConstIterator dictIter = dict.Begin(); while ( dictIter != dict.End() ) { - MITK_DEBUG << "Key " << dictIter->first; + //MITK_DEBUG << "Key " << dictIter->first; std::string value; if ( itk::ExposeMetaData( dict, dictIter->first, value ) ) { - MITK_DEBUG << "Value " << value; + //MITK_DEBUG << "Value " << value; TagToPropertyMapType::const_iterator valuePosition = propertyLookup.find( dictIter->first ); if ( valuePosition != propertyLookup.end() ) { std::string propertyKey = valuePosition->second; - MITK_DEBUG << "--> " << propertyKey; + //MITK_DEBUG << "--> " << propertyKey; image->SetProperty( propertyKey.c_str(), StringProperty::New(value) ); } } else { MITK_WARN << "Tag " << dictIter->first << " not read as string as expected. Ignoring..." ; } ++dictIter; } } } // end namespace mitk #include diff --git a/Core/Code/IO/mitkDicomSeriesReader.h b/Core/Code/IO/mitkDicomSeriesReader.h index 2e14dc9a14..d7484fb62d 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.h +++ b/Core/Code/IO/mitkDicomSeriesReader.h @@ -1,471 +1,687 @@ /*=================================================================== 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 mitkDicomSeriesReader_h #define mitkDicomSeriesReader_h #include "mitkDataNode.h" #include "mitkConfig.h" #include #include #include #ifdef NOMINMAX # define DEF_NOMINMAX # undef NOMINMAX #endif #include #ifdef DEF_NOMINMAX # ifndef NOMINMAX # define NOMINMAX # endif # undef DEF_NOMINMAX #endif #include #include #include #include #include #include namespace mitk { /** \brief Loading DICOM images as MITK images. - \ref DicomSeriesReader_purpose - \ref DicomSeriesReader_limitations - \ref DicomSeriesReader_usage - \ref DicomSeriesReader_sorting - \ref DicomSeriesReader_sorting1 - \ref DicomSeriesReader_sorting2 - \ref DicomSeriesReader_sorting3 - \ref DicomSeriesReader_sorting4 - \ref DicomSeriesReader_tests \section DicomSeriesReader_purpose Purpose DicomSeriesReader serves as a central class for loading DICOM images as mitk::Image. As the term "DICOM image" covers a huge variety of possible modalities and implementations, and since MITK assumes that 3D images are made up of continuous blocks of slices without any gaps or changes in orientation, the loading mechanism must implement a number of decisions and compromises. The main intention of this implementation is not efficiency but correcness of generated slice positions! \section DicomSeriesReader_limitations Assumptions and limitations The class is working only with GDCM 2.0.14 (or possibly newer). This version is the default of an MITK super-build. Support for other versions or ITK's DicomIO was dropped because of the associated complexity of DicomSeriesReader. \b Assumptions - expected to work for SOP Classes CT Image Storage and MR Image Storage (NOT for the "Enhanced" variants containing multi-frame images) - special treatment for a certain type of Philips 3D ultrasound (recogized by tag 3001,0010 set to "Philips3D") - loader will always attempt to read multiple single slices as a single 3D image volume (i.e. mitk::Image) - slices will be grouped by basic properties such as orientation, rows, columns, spacing and grouped into as large blocks as possible \b Options - images that cover the same piece of space (i.e. position, orientation, and dimensions are equal) can be interpreted as time-steps of the same image, i.e. a series will be loaded as 3D+t \b Limitations - the 3D+t assumption only works if all time-steps have an equal number of slices and if all have the Acquisition Time attribute set to meaningful values - Images from tilted CT gantries CAN ONLY be loaded as a series of single-slice images, since mitk::Image or the accompanying mapper are not (yet?) capable of representing such geometries - Secondary Capture images are expected to have the (0018,2010) tag describing the pixel spacing. If only the (0028,0030) tag is set, the spacing will be misinterpreted as (1,1) \section DicomSeriesReader_usage Usage The starting point for an application is a set of DICOM files that should be loaded. For convenience, DicomSeriesReader can also parse a whole directory for DICOM files, but an application should better know exactly what to load. Loading is then done in two steps: 1. Group the files into spatial blocks by calling GetSeries(). This method will sort all passed files into meaningful blocks that could fit into an mitk::Image. Sorting for 3D+t loading is optional but default. The \b return value of this function is a list of identifiers similar to DICOM UIDs, each associated to a sorted list of file names. 2. Load a sorted set of files by calling LoadDicomSeries(). This method expects go receive the sorting output of GetSeries(). The method will then invoke ITK methods to actually load the files into memory and put them into mitk::Images. Again, loading as 3D+t is optional. Example: \code // only a directory is known at this point: /home/who/dicom DicomSeriesReader::UidFileNamesMap allImageBlocks = DicomSeriesReader::GetSeries("/home/who/dicom/"); // file now divided into groups of identical image size, orientation, spacing, etc. // each of these lists should be loadable as an mitk::Image. DicomSeriesReader::StringContainer seriesToLoad = allImageBlocks[...]; // decide what to load // final step: load into DataNode (can result in 3D+t image) DataNode::Pointer node = DicomSeriesReader::LoadDicomSeries( oneBlockSorted ); Image::Pointer image = dynamic_cast( node->GetData() ); \endcode \section DicomSeriesReader_sorting Logic for sorting 2D slices from DICOM images into 3D+t blocks for mitk::Image The general sorting mechanism (implemented in GetSeries) groups and sorts a set of DICOM files, each assumed to contain a single CT/MR slice. In the following we refer to those file groups as "blocks", since this is what they are meant to become when loaded into an mitk::Image. \subsection DicomSeriesReader_sorting1 Step 1: Avoiding pure non-sense A first pass separates slices that cannot possibly be loaded together because of restrictions of mitk::Image. After this steps, each block contains only slices that match in all of the following DICOM tags: - (0020,0037) Image Orientation - (0028,0030) Pixel Spacing - (0018,0050) Slice Thickness - (0028,0010) Number Of Rows - (0028,0011) Number Of Columns - (0020,000e) Series Instance UID : could be argued about, might be dropped in the future (optionally) \subsection DicomSeriesReader_sorting2 Step 2: Sort slices spatially Before slices are further analyzed, they are sorted spatially. As implemented by GdcmSortFunction(), slices are sorted by 1. distance from origin (calculated using (0020,0032) Image Position Patient and (0020,0037) Image Orientation) 2. when distance is equal, (0008,0032) Acquisition Time is used as a backup criterion (necessary for meaningful 3D+t sorting) \subsection DicomSeriesReader_sorting3 Step 3: Ensure equal z spacing 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. 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. \subsection DicomSeriesReader_sorting4 Step 4 (optional): group 3D blocks as 3D+t when possible This last step depends on an option of GetSeries(). When requested, image blocks from the previous step are merged again whenever two blocks occupy the same portion of space (i.e. same origin, number of slices and z-spacing). + \section DicomSeriesReader_gantrytilt Handling of gantry tilt + + 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 in conjunction with mitk::Image, DicomSeriesReader performs a correction for such series + if the groupImagesWithGantryTilt or correctGantryTilt flag in GetSeries and LoadDicomSeries is set (default = on). + + The correction algorithms undoes 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 slices of a series, both described as a pair (origin/orientation): + - we calculate if the second origin is on a line along the normal of the first 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. + + \section DicomSeriesReader_whynotinitk Why is this not in ITK? + + Some of this code would probably be better located in ITK. It is just a matter of resources that this is not the + case yet. Any attempts into this direction are welcome and can be supported. At least the gantry tilt correction + should be a simple addition to itk::ImageSeriesReader. + \section DicomSeriesReader_tests Tests regarding DICOM loading A number of tests have been implemented to check our assumptions regarding DICOM loading. Please see \ref DICOMTesting + \todo refactor all the protected helper objects/methods into a separate header so we compile faster */ class MITK_CORE_EXPORT DicomSeriesReader { public: /** \brief Lists of filenames. */ typedef std::vector StringContainer; /** \brief For grouped lists of filenames, assigned an ID each. */ typedef std::map UidFileNamesMap; /** \brief Interface for the progress callback. */ typedef void (*UpdateCallBackMethod)(float); /** \brief Provide combination of preprocessor defines that was active during compilation. Since this class is a combination of several possible implementations, separated only by ifdef's, calling instances might want to know which flags were active at compile time. */ static std::string GetConfigurationString(); /** \brief Checks if a specific file contains DICOM data. */ static bool IsDicom(const std::string &filename); /** \brief see other GetSeries(). Find all series (and sub-series -- see details) in a particular directory. */ static UidFileNamesMap GetSeries(const std::string &dir, + bool groupImagesWithGantryTilt, const StringContainer &restrictions = StringContainer()); /** \brief see other GetSeries(). \warning Untested, could or could not work. This differs only by having an additional restriction to a single known DICOM series. Internally, it uses the other GetSeries() method. */ static StringContainer GetSeries(const std::string &dir, const std::string &series_uid, + bool groupImagesWithGantryTilt, const StringContainer &restrictions = StringContainer()); /** \brief PREFERRED version of this method - scan and sort DICOM files. Parse a list of files for images of DICOM series. For each series, an enumeration of the files contained in it is created. \return The resulting maps UID-like keys (based on Series Instance UID and slice properties) to sorted lists of file names. SeriesInstanceUID will be enhanced to be unique for each set of file names that is later loadable as a single mitk::Image. This implies that Image orientation, slice thickness, pixel spacing, rows, and columns must be the same for each file (i.e. the image slice contained in the file). If this separation logic requires that a SeriesInstanceUID must be made more specialized, it will follow the same logic as itk::GDCMSeriesFileNames to enhance the UID with more digits and dots. Optionally, more tags can be used to separate files into different logical series by setting the restrictions parameter. \warning Adding restrictions is not yet implemented! */ static UidFileNamesMap - GetSeries(const StringContainer& files, bool sortTo3DPlust, const StringContainer &restrictions = StringContainer()); + GetSeries(const StringContainer& files, + bool sortTo3DPlust, + bool groupImagesWithGantryTilt, + const StringContainer &restrictions = StringContainer()); /** \brief See other GetSeries(). Use GetSeries(const StringContainer& files, bool sortTo3DPlust, const StringContainer &restrictions) instead. */ static UidFileNamesMap - GetSeries(const StringContainer& files, const StringContainer &restrictions = StringContainer()); + GetSeries(const StringContainer& files, + bool groupImagesWithGantryTilt, + const StringContainer &restrictions = StringContainer()); /** Loads a DICOM series composed by the file names enumerated in the file names container. If a callback method is supplied, it will be called after every progress update with a progress value in [0,1]. \param filenames The filenames to load. \param sort Whether files should be sorted spatially (true) or not (false - maybe useful if presorted) \param load4D Whether to load the files as 3D+t (if possible) */ static DataNode::Pointer LoadDicomSeries(const StringContainer &filenames, bool sort = true, bool load4D = true, + bool correctGantryTilt = true, UpdateCallBackMethod callback = 0); /** \brief See LoadDicomSeries! Just a slightly different interface. */ static bool LoadDicomSeries(const StringContainer &filenames, DataNode &node, bool sort = true, bool load4D = true, + bool correctGantryTilt = true, UpdateCallBackMethod callback = 0); protected: + /** + \brief Return type of DicomSeriesReader::AnalyzeFileForITKImageSeriesReaderSpacingAssumption. + + Class contains the grouping result of method DicomSeriesReader::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. + */ + StringContainer GetBlockFilenames(); + + /** + \brief Remaining files, which could not be grouped. + */ + StringContainer GetUnsortedFilenames(); + + /** + \brief Wheter or not the grouped result contain a gantry tilt. + */ + bool ContainsGantryTilt(); + + /** + \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. + */ + void AddFileToSortedBlock(const std::string& filename); + + /** + \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. + */ + void AddFileToUnsortedBlock(const std::string& filename); + + /** + \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(); + + /** + \brief Only meaningful for use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption. + */ + void UndoPrematureGrouping(); + + protected: + + StringContainer m_GroupedFiles; + StringContainer m_UnsortedFiles; + + bool m_GantryTilt; + }; + + /** + \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. + + All calculations are done in the constructor, results can then + be read via the remaining methods. + */ + class GantryTiltInformation + { + public: + + /** + \brief Just so we can create empty instances for assigning results later. + */ + GantryTiltInformation(); + + /** + \brief THE constructor, which does all the calculations. + + See code comments for explanation. + */ + GantryTiltInformation( const Point3D& origin1, + const Point3D& origin2, + const Vector3D& right, + const Vector3D& up, + unsigned int numberOfSlicesApart); + + /** + \brief Whether the slices were sheared. + */ + bool IsSheared() const; + + /** + \brief Whether the shearing is a gantry tilt or more complicated. + */ + bool IsRegularGantryTilt() const; + + /** + \brief The offset distance in Y direction for each slice (describes the tilt result). + */ + ScalarType GetMatrixCoefficientForCorrectionInWorldCoordinates() const; + + + /** + \brief The z / inter-slice spacing. Needed to correct ImageSeriesReader's result. + */ + ScalarType GetRealZSpacing() const; + + /** + \brief The shift between first and last slice in mm. + + Needed to resize an orthogonal image volume. + */ + ScalarType GetTiltCorrectedAdditionalSize() const; + + protected: + + /** + \brief Projection of point p onto line through lineOrigin in direction of lineDirection. + */ + Point3D projectPointOnLine( Point3D p, Point3D lineOrigin, Vector3D lineDirection ); + + ScalarType m_ShiftUp; + ScalarType m_ShiftRight; + ScalarType m_ShiftNormal; + unsigned int m_NumberOfSlicesApart; + }; + /** \brief for internal sorting. */ typedef std::pair TwoStringContainers; /** \brief Maps DICOM tags to MITK properties. */ typedef std::map TagToPropertyMapType; /** \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) */ static - TwoStringContainers - AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const StringContainer& files, const gdcm::Scanner::MappingType& tagValueMappings_); + SliceGroupingAnalysisResult + AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const StringContainer& files, bool groupsOfSimilarImages, const gdcm::Scanner::MappingType& tagValueMappings_); + + static + std::string + ConstCharStarToString(const char* s); + + static + Point3D + DICOMStringToPoint3D(const std::string& s, bool& successful); + + static + void + DICOMStringToOrientationVectors(const std::string& s, Vector3D& right, Vector3D& up, bool& successful); + + template + static + typename ImageType::Pointer + InPlaceFixUpTiltedGeometry( ImageType* input, const GantryTiltInformation& tiltInfo ); + /** \brief Sort a set of file names in an order that is meaningful for loading them into an mitk::Image. \warning This method assumes that input files are similar in basic properties such as slice thicknes, image orientation, pixel spacing, rows, columns. It should always be ok to put the result of a call to GetSeries(..) into this method. Sorting order is determined by 1. image position along its normal (distance from world origin) 2. acquisition time If P denotes a position and T denotes a time step, this method will order slices from three timesteps like this: \verbatim P1T1 P1T2 P1T3 P2T1 P2T2 P2T3 P3T1 P3T2 P3T3 \endverbatim */ static StringContainer SortSeriesSlices(const StringContainer &unsortedFilenames); public: /** \brief Checks if a specific file is a Philips3D ultrasound DICOM file. */ static bool IsPhilips3DDicom(const std::string &filename); protected: /** \brief Read a Philips3D ultrasound DICOM file and put into an mitk::Image. */ static bool ReadPhilips3DDicom(const std::string &filename, mitk::Image::Pointer output_image); /** \brief Construct a UID that takes into account sorting criteria from GetSeries(). */ static std::string CreateMoreUniqueSeriesIdentifier( gdcm::Scanner::TagToValue& tagValueMap ); /** \brief Helper for CreateMoreUniqueSeriesIdentifier */ static std::string CreateSeriesIdentifierPart( gdcm::Scanner::TagToValue& tagValueMap, const gdcm::Tag& tag ); /** \brief Helper for CreateMoreUniqueSeriesIdentifier */ static std::string IDifyTagValue(const std::string& value); typedef itk::GDCMImageIO DcmIoType; /** \brief Progress callback for DicomSeriesReader. */ class CallbackCommand : public itk::Command { public: CallbackCommand(UpdateCallBackMethod callback) : m_Callback(callback) { } void Execute(const itk::Object *caller, const itk::EventObject&) { (*this->m_Callback)(static_cast(caller)->GetProgress()); } void Execute(itk::Object *caller, const itk::EventObject&) { (*this->m_Callback)(static_cast(caller)->GetProgress()); } protected: UpdateCallBackMethod m_Callback; }; /** \brief Scan for slice image information */ static void ScanForSliceInformation( const StringContainer &filenames, gdcm::Scanner& scanner ); /** \brief Performs actual loading of a series and creates an image having the specified pixel type. */ template static void - LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, UpdateCallBackMethod callback); + LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, bool correctTilt, UpdateCallBackMethod callback); /** \brief 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); + LoadDICOMByITK( const StringContainer&, bool correctTilt, const GantryTiltInformation& tiltInfo, CallbackCommand* command = NULL); /** \brief Sort files into time step blocks of a 3D+t image. 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. \todo We can probably remove this method if we somehow transfer 3D+t information from GetSeries to LoadDicomSeries. */ static std::list SortIntoBlocksFor3DplusT( const StringContainer& presortedFilenames, const gdcm::Scanner::MappingType& tagValueMappings_, bool sort, bool& canLoadAs4D); /** \brief Defines spatial sorting for sorting by GDCM 2. Sorts by image position along image normal (distance from world origin). In cases of conflict, acquisition time is used as a secondary sort criterium. */ static bool GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::DataSet &ds2); /** \brief Copy information about files and DICOM tags from ITK's MetaDataDictionary and from the list of input files to the PropertyList of mitk::Image. \todo Tag copy must follow; image level will cause some additional files parsing, probably. */ static void CopyMetaDataToImageProperties( StringContainer filenames, const gdcm::Scanner::MappingType& tagValueMappings_, DcmIoType* io, Image* image); static void CopyMetaDataToImageProperties( std::list imageBlock, const gdcm::Scanner::MappingType& tagValueMappings_, DcmIoType* io, Image* image); /** \brief Map between DICOM tags and MITK properties. Uses as a positive list for copying specified DICOM tags (from ITK's ImageIO) to MITK properties. ITK provides MetaDataDictionary entries of form "gggg|eeee" (g = group, e = element), e.g. "0028,0109" (Largest Pixel in Series), which we want to sort as dicom.series.largest_pixel_in_series". */ static const TagToPropertyMapType& GetDICOMTagsToMITKPropertyMap(); }; } #endif /* MITKDICOMSERIESREADER_H_ */ diff --git a/Core/Code/IO/mitkDicomSeriesReader.txx b/Core/Code/IO/mitkDicomSeriesReader.txx index a67ba660fc..80c95e54c9 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.txx +++ b/Core/Code/IO/mitkDicomSeriesReader.txx @@ -1,317 +1,489 @@ /*=================================================================== 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 MITKDICOMSERIESREADER_TXX_ #define MITKDICOMSERIESREADER_TXX_ #include #include #include +#include +#include +#include + +#include + namespace mitk { template -void DicomSeriesReader::LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool load4D, UpdateCallBackMethod callback) +void DicomSeriesReader::LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool load4D, bool correctTilt, UpdateCallBackMethod callback) { const char* previousCLocale = setlocale(LC_NUMERIC, NULL); setlocale(LC_NUMERIC, "C"); std::locale previousCppLocale( std::cin.getloc() ); std::locale l( "C" ); std::cin.imbue(l); + + const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) + const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation try { mitk::Image::Pointer image = mitk::Image::New(); CallbackCommand *command = callback ? new CallbackCommand(callback) : 0; bool initialize_node = false; /* special case for Philips 3D+t ultrasound images */ if ( DicomSeriesReader::IsPhilips3DDicom(filenames.front().c_str()) ) { ReadPhilips3DDicom(filenames.front().c_str(), image); initialize_node = true; } else { /* default case: assume "normal" image blocks, possibly 3D+t */ bool canLoadAs4D(true); gdcm::Scanner scanner; ScanForSliceInformation(filenames, scanner); + + // need non-const access for map + gdcm::Scanner::MappingType& tagValueMappings = const_cast(scanner.GetMappings()); - std::list imageBlocks = SortIntoBlocksFor3DplusT( filenames, scanner.GetMappings(), sort, canLoadAs4D ); + std::list imageBlocks = SortIntoBlocksFor3DplusT( filenames, tagValueMappings, sort, canLoadAs4D ); unsigned int volume_count = imageBlocks.size(); + GantryTiltInformation tiltInfo; + + // check possibility of a single slice with many timesteps. In this case, don't check for tilt, no second slice possible + if ( !imageBlocks.empty() && imageBlocks.front().size() > 1 && correctTilt) + { + // check tiltedness here, potentially fixup ITK's loading result by shifting slice contents + // check first and last position slice from tags, make some calculations to detect tilt + + std::string firstFilename(imageBlocks.front().front()); + // calculate from first and last slice to minimize rounding errors + std::string secondFilename(imageBlocks.front().back()); + + std::string imagePosition1( ConstCharStarToString( tagValueMappings[ firstFilename.c_str() ][ tagImagePositionPatient ] ) ); + std::string imageOrientation( ConstCharStarToString( tagValueMappings[ firstFilename.c_str() ][ tagImageOrientation ] ) ); + std::string imagePosition2( ConstCharStarToString( tagValueMappings[secondFilename.c_str() ][ tagImagePositionPatient ] ) ); + + bool ignoredConversionError(-42); // hard to get here, no graceful way to react + Point3D origin1( DICOMStringToPoint3D( imagePosition1, ignoredConversionError ) ); + Point3D origin2( DICOMStringToPoint3D( imagePosition2, ignoredConversionError ) ); + + 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 + DICOMStringToOrientationVectors( imageOrientation, right, up, ignoredConversionError ); + + tiltInfo = GantryTiltInformation ( origin1, origin2, right, up, filenames.size()-1 ); + correctTilt = tiltInfo.IsSheared() && tiltInfo.IsRegularGantryTilt(); + + MITK_DEBUG << "** Loading now: shear? " << tiltInfo.IsSheared(); + MITK_DEBUG << "** Loading now: normal tilt? " << tiltInfo.IsRegularGantryTilt(); + MITK_DEBUG << "** Loading now: perform tilt correction? " << correctTilt; + } + else + { + correctTilt = false; // we CANNOT do that + } + if (volume_count == 1 || !canLoadAs4D || !load4D) { - image = LoadDICOMByITK( imageBlocks.front() , command ); // load first 3D block + image = LoadDICOMByITK( imageBlocks.front(), correctTilt, tiltInfo, command ); // load first 3D block initialize_node = true; } else if (volume_count > 1) { // It is 3D+t! Read it and store into mitk image typedef itk::Image ImageType; typedef itk::ImageSeriesReader ReaderType; DcmIoType::Pointer io = DcmIoType::New(); typename ReaderType::Pointer reader = ReaderType::New(); reader->SetImageIO(io); reader->ReverseOrderOff(); if (command) { reader->AddObserver(itk::ProgressEvent(), command); } unsigned int act_volume = 1u; reader->SetFileNames(imageBlocks.front()); reader->Update(); - image->InitializeByItk(reader->GetOutput(), 1, volume_count); - image->SetImportVolume(reader->GetOutput()->GetBufferPointer(), 0u); + + typename ImageType::Pointer readVolume = reader->GetOutput(); + // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right position + if (correctTilt) + { + readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); + } + + image->InitializeByItk( readVolume.GetPointer(), 1, volume_count); + image->SetImportVolume( readVolume->GetBufferPointer(), 0u); gdcm::Scanner scanner; ScanForSliceInformation(filenames, scanner); DicomSeriesReader::CopyMetaDataToImageProperties( imageBlocks, scanner.GetMappings(), io, image); 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: // 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 ); #endif MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " << image->GetGeometry()->GetSpacing()[1] << ", " << image->GetGeometry()->GetSpacing()[2] << "]"; for (std::list::iterator df_it = ++imageBlocks.begin(); df_it != imageBlocks.end(); ++df_it) { reader->SetFileNames(*df_it); reader->Update(); - image->SetImportVolume(reader->GetOutput()->GetBufferPointer(), act_volume++); + readVolume = reader->GetOutput(); + + if (correctTilt) + { + readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); + } + + image->SetImportVolume(readVolume->GetBufferPointer(), act_volume++); } initialize_node = true; } + } - if (initialize_node) - { - // forward some image properties to node + if (initialize_node) + { + // forward some image properties to node node.GetPropertyList()->ConcatenatePropertyList( image->GetPropertyList(), true ); 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; } } template -Image::Pointer DicomSeriesReader::LoadDICOMByITK( const StringContainer& filenames, CallbackCommand* command ) +Image::Pointer DicomSeriesReader::LoadDICOMByITK( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, CallbackCommand* command ) { /******** Normal Case, 3D (also for GDCM < 2 usable) ***************/ mitk::Image::Pointer image = mitk::Image::New(); typedef itk::Image ImageType; typedef itk::ImageSeriesReader ReaderType; DcmIoType::Pointer io = DcmIoType::New(); typename ReaderType::Pointer reader = ReaderType::New(); 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()); + typename ImageType::Pointer readVolume = reader->GetOutput(); + + // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right position + if (correctTilt) + { + readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); + } + + image->InitializeByItk(readVolume.GetPointer()); + image->SetImportVolume(readVolume->GetBufferPointer()); gdcm::Scanner scanner; ScanForSliceInformation(filenames, scanner); DicomSeriesReader::CopyMetaDataToImageProperties( filenames, scanner.GetMappings(), io, image); 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 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] << "]"; return image; } void DicomSeriesReader::ScanForSliceInformation(const StringContainer &filenames, gdcm::Scanner& scanner) { const gdcm::Tag ippTag(0x0020,0x0032); //Image position (Patient) scanner.AddTag(ippTag); + + const gdcm::Tag tagImageOrientation(0x0020,0x0037); //Image orientation + scanner.AddTag(tagImageOrientation); // TODO what if tags don't exist? const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location scanner.AddTag( tagSliceLocation ); const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number scanner.AddTag( tagInstanceNumber ); const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number scanner.AddTag( tagSOPInstanceNumber ); scanner.Scan(filenames); // make available image position for each file } std::list DicomSeriesReader::SortIntoBlocksFor3DplusT( const StringContainer& presortedFilenames, const gdcm::Scanner::MappingType& tagValueMappings, bool /*sort*/, bool& canLoadAs4D ) { std::list imageBlocks; // ignore sort request, because most likely re-sorting is now needed due to changes in GetSeries(bug #8022) StringContainer sorted_filenames = DicomSeriesReader::SortSeriesSlices(presortedFilenames); std::string firstPosition; unsigned int numberOfBlocks(0); // number of 3D image blocks const gdcm::Tag ippTag(0x0020,0x0032); //Image position (Patient) // loop files to determine number of image blocks for (StringContainer::const_iterator fileIter = sorted_filenames.begin(); fileIter != sorted_filenames.end(); ++fileIter) { gdcm::Scanner::TagToValue tagToValueMap = tagValueMappings.find( fileIter->c_str() )->second; if(tagToValueMap.find(ippTag) == tagToValueMap.end()) { continue; } std::string position = tagToValueMap.find(ippTag)->second; MITK_DEBUG << " " << *fileIter << " at " << position; if (firstPosition.empty()) { firstPosition = position; } if ( position == firstPosition ) { ++numberOfBlocks; } else { break; // enough information to know the number of image blocks } } MITK_DEBUG << " ==> Assuming " << numberOfBlocks << " time steps"; 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) { StringContainer filesOfCurrentBlock; 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; } +template +typename ImageType::Pointer +DicomSeriesReader::InPlaceFixUpTiltedGeometry( ImageType* input, const GantryTiltInformation& tiltInfo ) +{ + typedef itk::ResampleImageFilter ResampleFilterType; + typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); + resampler->SetInput( input ); + + /* + Transform for a point is + - transform from actual position to index coordinates + - apply a shear that undoes the gantry tilt + - transform back into world coordinates + + Anybody who does this in a simpler way: don't forget to write up how and why your solution works + */ + typedef itk::AffineTransform< double, ImageType::ImageDimension > TransformType; + typename TransformType::Pointer transformShear = TransformType::New(); + + /** + What I think should be done here: + + - apply a shear and spacing correction to the image block that corrects the ITK reader's error + - ITK ignores the shear and loads slices into an orthogonal volume + - ITK calculates the spacing from the origin distance, which is more than the actual spacing with gantry tilt images + - to undo the effect + - we have calculated some information in tiltInfo: + - the shift in Y direction that is added with each additional slice is the most important information + - the Y-shift is calculated in mm world coordinates + - we apply a shearing transformation to the ITK-read image volume + - to do this locally, + - we transform the image volume back to origin and "normal" orientation by applying the inverse of its transform + - this should bring us into the image's "index coordinate" system + - we apply a shear with the Y-shift factor put into a unit transform at row 1, col 2 + - we probably would need the shift factor in index coordinates but we use mm, which appears strange, see below + - we transform the image volume back to its actual position + - presumably back from index to world coordinates + - we lastly apply modify the image spacing in z direction by replacing this number with the correctly calulcated inter-slice distance + + Here comes the unsolved PROBLEM: + - WHY is it correct to apply the shift factor in millimeters, when we assume we are in a index coordinate system? + - or why is it not millimeters but index coordinates? + - or what else is the wrong assumption here? + - This is a serious question to anybody very firm in transforms, ITK, etc. + - this is also a chocolate reward question. Provide a good explanation with corrected code as a git commit and get it. + */ + + // ScalarType factor = tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() / input->GetSpacing()[1]; + ScalarType factor = tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates(); + // row 1, column 2 corrects shear in parallel to Y axis, proportional to distance in Z direction + transformShear->Shear( 1, 2, factor ); + + typename TransformType::Pointer imageIndexToWorld = TransformType::New(); + imageIndexToWorld->SetOffset( input->GetOrigin().GetVectorFromOrigin() ); + imageIndexToWorld->SetMatrix( input->GetDirection() ); + + typename TransformType::Pointer imageWorldToIndex = TransformType::New(); + imageIndexToWorld->GetInverse( imageWorldToIndex ); + + typename TransformType::Pointer gantryTiltCorrection = TransformType::New(); + gantryTiltCorrection->Compose( imageWorldToIndex ); + gantryTiltCorrection->Compose( transformShear ); + gantryTiltCorrection->Compose( imageIndexToWorld ); + + resampler->SetTransform( gantryTiltCorrection ); + + typedef itk::LinearInterpolateImageFunction< ImageType, double > InterpolatorType; + typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); + resampler->SetInterpolator( interpolator ); + /* + This would be the right place to invent a meaningful value for positions outside of the image. + For CT, HU -1000 might be meaningful, but a general solution seems not possible. Even for CT, + -1000 would only look natural for many not all images. + */ + resampler->SetDefaultPixelValue( std::numeric_limits::min() ); + + // adjust size in Y direction! (maybe just transform the outer last pixel to see how much space we would need + + resampler->SetOutputParametersFromImage( input ); // we basically need the same image again, just sheared + + typename ImageType::SizeType largerSize = resampler->GetSize(); // now the resampler already holds the input image's size. + largerSize[1] += tiltInfo.GetTiltCorrectedAdditionalSize(); + resampler->SetSize( largerSize ); + + resampler->Update(); + typename ImageType::Pointer result = resampler->GetOutput(); + + // ImageSeriesReader calculates z spacing as the distance between the first two origins. + // This is not correct in case of gantry tilt, so we set our calculated spacing. + typename ImageType::SpacingType correctedSpacing = result->GetSpacing(); + correctedSpacing[2] = tiltInfo.GetRealZSpacing(); + result->SetSpacing( correctedSpacing ); + + return result; +} + + } #endif diff --git a/Core/Code/Testing/DICOMTesting/DICOMTesting.dox b/Core/Code/Testing/DICOMTesting/DICOMTesting.dox.h similarity index 94% rename from Core/Code/Testing/DICOMTesting/DICOMTesting.dox rename to Core/Code/Testing/DICOMTesting/DICOMTesting.dox.h index 15aea78acd..96e8bb67b3 100644 --- a/Core/Code/Testing/DICOMTesting/DICOMTesting.dox +++ b/Core/Code/Testing/DICOMTesting/DICOMTesting.dox.h @@ -1,144 +1,148 @@ /** \page DICOMTesting MITK DICOM testing \section DICOMTesting_introduction Introduction Reading DICOM data into mitk::Images is a complicated business since DICOM and MITK have very different ideas of images. DicomSeriesReader brings DICOM and MITK as close together as possible by offering methods to load DICOM CT and MR images into mitk::Images, optionally grouping slices to 3D+t images. Since there are so many possible sources for mistakes with any change to this loading process, testing the many assumptions implemented in DicomSeriesReader is worthwhile. This document describes what and how theses kind of tests are implemented. \section DICOMTesting_problem Problem description The task of loading DICOM files into mitk::Images is a challenge because of major differences in the way that DICOM and MITK represent images: - DICOM images - are mostly stored as one slice per file - do not describe how they can be arraged into image stacks with orthogonal axis - sometimes they cannot be arranged in image stacks as described above (e.g. tilted gantry) - mitk::Image (at least its mature areas) - represents image stacks with orthogonal axis (nothing like a tilted gantry) - have a concept of a fourth dimension (time steps) Because image processing applications based on MITK still need a way to load DICOM images into mitk::Image objects, MITK needs a way to build image stacks and this needs to be well tested. For more background information, see David Clunie's most valuable posts on comp.protocols.dicom, e.g.: - http://groups.google.com/group/comp.protocols.dicom/browse_thread/thread/6db91972e161f0d4/6e0304ac264a6eb5 - http://groups.google.com/group/comp.protocols.dicom/browse_thread/thread/e9bd1497bea3e66b/187a7dc8810613d2 - http://groups.google.com/group/comp.protocols.dicom/browse_thread/thread/5d80bb0b7fafcb81/cf96119e3b024ed8 - http://groups.google.com/group/comp.protocols.dicom/browse_thread/thread/4568635e083a3fba/e2a8ceec23032601 \section DICOMTesting_testidea Test principle The general idea for DICOM loaing tests is to run a set of known DICOM files through DicomSeriesReader's methods GetSeries() and LoadDicomSeries() to generate mitk::Images. These images are then compared to expected image properties, such as the number of individual mitk::Images, positions, orientations, spacings, etc. Stored expectations look like this (should be self-explanatory): \verbatim -- Image 1 Pixeltype: s BitsPerPixel: 16 Dimension: 4 Dimensions: 64 64 6 1 Geometry: Matrix: 5.25 0 0 0 5.2468 0.139598 0 -0.183222 3.99757 Offset: -159.672 -309.974 -69.0122 Center: 0 0 0 Translation: -159.672 -309.974 -69.0122 Scale: 1 1 1 Origin: -159.672 -309.974 -69.0122 Spacing: 5.25 5.25 4 TimeBounds: 0 1 -- Image 2 Pixeltype: s BitsPerPixel: 16 Dimension: 4 Dimensions: 64 64 41 1 Geometry: Matrix: 5.25 0 0 0 5.25 0 0 0 4 Offset: -160.672 -311.672 -285 Center: 0 0 0 Translation: -160.672 -311.672 -285 Scale: 1 1 1 Origin: -160.672 -311.672 -285 Spacing: 5.25 5.25 4 TimeBounds: 0 1 \endverbatim \section DICOMTesting_implementation Implementation \section DICOMTesting_implementation_utils Test helpers (applications and classes) Application DumpDICOMMitkImage Takes a list of DICOM images, loads them using TestDICOMLoading, then dumps information about the resulting mitk::Images to standard output. This application is helpful when defining reference data for tests. Application VerifyDICOMMitkImageDump Takes a list of DICOM images and loads them using TestDICOMLoading. Takes a dump file as generated by DumpDICOMMitkImage, parses it and compares it to actually generated mitk::Images. This application is used to implement the majority of test cases. They all load images, then verify the expected result structure. Class TestDICOMLoading \section PageDICOMLoadingTests_testcaseimplementation Test case implementation Individual test cases are stored in the MITK-Data repository and constructed by Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt The CMake file parses given directories for subdirectories containing specific test cases. Each directory contains two files: - File "input": lists DICOM files that should be loaded for this test case - File "expected.dump": contains the image properties in the above mentioned dump format Each test case is translated into a CTest test which evaluates the return value of a call to VerifyDICOMMitkImageDump. \section PageDICOMLoadingTests_testcases Implemented test cases -From test set TinyCTAbdomen (see description.txt in this directory for details on test images): +From test set \b TinyCTAbdomen (see description.txt in this directory for details on test images): - singleslice : just a single slice (should work and contain meaningful spacing) - two_slices : two slices, spacing should be calculated correctly - all : load a "normal" series as a single 3D block - 3D_and_T : load a small set of slices with multiple time-steps - - diff_orientation : load a set of files containing two differently oriented image blocks + - diff_orientation : load a set of files containing two differently oriented image blocks; at least two images (110,111) have minor errors in small decimal digits - diff_orientation_gaps : load a set of files containing two differently oriented image blocks, each missing slices, so blocks must be split - diff_spacing : load a set of files containint two set of slices with different spacings - gap : load slices that cannot form a single 3D block, because single slices are missing - gaps : slices missing in differnt locations, so multiple splits needed - unsorted_gaps : as before, just file names are unsorted - single_negative_spacing : from reported bug related to single MR images with misleading spacing information - tilted_gantry : slice origins do not align along first slice normal (happens with tilted gantries) - interleaved : two volumes of slices with origins along the same line. The volumes' slices interleave in their border region. This test is meant to correctly sort apart the two blocks instead of generating many two-slices groups in the interleaved region. +From test set \b TiltHead (see description.txt in this directory for details on test images): + + - head_ct_tilt : load a series with gantry tilt and check its orientation + \section DICOMTesting_othertests Specific other tests This list is meant to provide an up-to-date list of all implemented DICOM loading tests. If you ever find this outdated, please update it or make the persons who invalidated the list update it. mitkDICOMTestingSanityTest_* These tests implement basic testing of the implemented helper classes. The tests use DicomSeriesReader to load a number of DICOM image. They verify: - DicomSeriesReader recognizes all input files as DICOM images - DicomSeriesReader generates a number of mitk::Images from the DICOM images - the number of mitk::Images matches a number given on the command line or CTest's add_test() - helper methods in class TestDICOMLoading make minimal sense (comparison of an image dump to itself must be valid) */ diff --git a/Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt b/Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt index 50def10d0c..5a6abc3527 100644 --- a/Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt +++ b/Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt @@ -1,67 +1,77 @@ MITK_CREATE_MODULE_TESTS(LABELS MITK-Core) include(mitkFunctionAddTestLabel) # verify minimum expectations: # files are recognized as DICOM # loading files results in a given number of images mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_NoFiles mitkDICOMTestingSanityTest 0) mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_CTImage mitkDICOMTestingSanityTest 1 ${MITK_DATA_DIR}/spacing-ok-ct.dcm) mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_MRImage mitkDICOMTestingSanityTest 1 ${MITK_DATA_DIR}/spacing-ok-mr.dcm) mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_SCImage mitkDICOMTestingSanityTest 1 ${MITK_DATA_DIR}/spacing-ok-sc.dcm) mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_DefectImage mitkDICOMTestingSanityTest 0 ${MITK_DATA_DIR}/spacing-ok-sc-no2032.dcm) #see bug 8108 if(NOT APPLE) set(VERIFY_DUMP_CMD ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/VerifyDICOMMitkImageDump) set(CT_ABDOMEN_DIR ${MITK_DATA_DIR}/TinyCTAbdomen) +set(CT_TILT_DIR ${MITK_DATA_DIR}/TiltHead) # these variables could be passed as parameters to a generic test creation function -set(CURRENT_DATASET_DIR ${CT_ABDOMEN_DIR}) set(TESTS_DIR Tests) set(INPUT_LISTNAME input) set(EXPECTED_DUMP expected.dump) -# find all test input lists -file(GLOB_RECURSE allInputs ${CURRENT_DATASET_DIR}/${TESTS_DIR}/*/${INPUT_LISTNAME}) - -function(expectFileExists filename) - if(NOT EXISTS ${filename}) - message(SEND_ERROR "Test case expected file ${filename} which does not exist! Please fix your CMake code or file layout.") - endif(NOT EXISTS ${filename}) -endfunction(expectFileExists) - -foreach(input ${allInputs}) - # extract only the name of the directory of this very test case - string(REGEX REPLACE ".*${TESTS_DIR}/([^/]+)/.*" "\\1" input ${input}) - set(inputfilelist "${CURRENT_DATASET_DIR}/${TESTS_DIR}/${input}/${INPUT_LISTNAME}") - set(expecteddump "${CURRENT_DATASET_DIR}/${TESTS_DIR}/${input}/${EXPECTED_DUMP}") - set(testname "DICOM_Load_${input}") - - message(STATUS "DICOM loading test case '${input}'") - expectFileExists(${inputfilelist}) - expectFileExists(${expecteddump}) - - # TODO provide error messages if input not valid - - set(dicomFiles) # clear list - # read list of file names from file "input" - file(STRINGS ${inputfilelist} rawDicomFiles) - foreach(raw ${rawDicomFiles}) - # prepend each file with an abosolute path - set(fileWithFullPath ${CURRENT_DATASET_DIR}/${raw}) - list(APPEND dicomFiles ${fileWithFullPath}) - endforeach(raw ${rawDicomFiles}) - - #message(STATUS " Load ${dicomFiles}") - add_test(${testname} - ${VERIFY_DUMP_CMD} ${expecteddump} ${dicomFiles}) - mitkFunctionAddTestLabel(${testname}) - -endforeach(input allInputs) -endif() +function(AddDicomTestsFromDataRepository CURRENT_DATASET_DIR TESTS_DIR INPUT_LISTNAME EXPECTED_DUMP) + + # find all test input lists + file(GLOB_RECURSE allInputs ${CURRENT_DATASET_DIR}/${TESTS_DIR}/*/${INPUT_LISTNAME}) + + function(expectFileExists filename) + if(NOT EXISTS ${filename}) + message(SEND_ERROR "Test case expected file ${filename} which does not exist! Please fix your CMake code or file layout.") + endif(NOT EXISTS ${filename}) + endfunction(expectFileExists) + + foreach(input ${allInputs}) + # extract only the name of the directory of this very test case + string(REGEX REPLACE ".*${TESTS_DIR}/([^/]+)/.*" "\\1" input ${input}) + set(inputfilelist "${CURRENT_DATASET_DIR}/${TESTS_DIR}/${input}/${INPUT_LISTNAME}") + set(expecteddump "${CURRENT_DATASET_DIR}/${TESTS_DIR}/${input}/${EXPECTED_DUMP}") + set(testname "DICOM_Load_${input}") + + message(STATUS "DICOM loading test case '${input}'") + expectFileExists(${inputfilelist}) + expectFileExists(${expecteddump}) + + # TODO provide error messages if input not valid + + set(dicomFiles) # clear list + # read list of file names from file "input" + file(STRINGS ${inputfilelist} rawDicomFiles) + foreach(raw ${rawDicomFiles}) + # prepend each file with an abosolute path + set(fileWithFullPath ${CURRENT_DATASET_DIR}/${raw}) + list(APPEND dicomFiles ${fileWithFullPath}) + endforeach(raw ${rawDicomFiles}) + + #message(STATUS " Load ${dicomFiles}") + add_test(${testname} + ${VERIFY_DUMP_CMD} ${expecteddump} ${dicomFiles}) + mitkFunctionAddTestLabel(${testname}) + + endforeach(input allInputs) + +endfunction(AddDicomTestsFromDataRepository CURRENT_DATASET_DIR TESTS_DIR INPUT_LISTNAME EXPECTED_DUMP) + + +AddDicomTestsFromDataRepository(${CT_ABDOMEN_DIR} ${TESTS_DIR} ${INPUT_LISTNAME} ${EXPECTED_DUMP}) +AddDicomTestsFromDataRepository(${CT_TILT_DIR} ${TESTS_DIR} ${INPUT_LISTNAME} ${EXPECTED_DUMP}) + + +endif() # apple diff --git a/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.cpp b/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.cpp index 691d34c4bd..b74a1c613b 100644 --- a/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.cpp +++ b/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.cpp @@ -1,432 +1,432 @@ /*=================================================================== 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 "mitkTestDICOMLoading.h" #include mitk::TestDICOMLoading::TestDICOMLoading() :m_PreviousCLocale(NULL) { } void mitk::TestDICOMLoading::SetDefaultLocale() { // remember old locale only once if (m_PreviousCLocale == NULL) { m_PreviousCLocale = setlocale(LC_NUMERIC, NULL); // set to "C" setlocale(LC_NUMERIC, "C"); m_PreviousCppLocale = std::cin.getloc(); std::locale l( "C" ); std::cin.imbue(l); std::cout.imbue(l); } } void mitk::TestDICOMLoading::ResetUserLocale() { if (m_PreviousCLocale) { setlocale(LC_NUMERIC, m_PreviousCLocale); std::cin.imbue(m_PreviousCppLocale); std::cout.imbue(m_PreviousCppLocale); m_PreviousCLocale = NULL; } } mitk::TestDICOMLoading::ImageList mitk::TestDICOMLoading::LoadFiles( const StringContainer& files ) { for (StringContainer::const_iterator iter = files.begin(); iter != files.end(); ++iter) { MITK_DEBUG << "File " << *iter; } ImageList result; - DicomSeriesReader::UidFileNamesMap seriesInFiles = DicomSeriesReader::GetSeries( files ); + DicomSeriesReader::UidFileNamesMap seriesInFiles = DicomSeriesReader::GetSeries( files, true ); // TODO sort series UIDs, implementation of map iterator might differ on different platforms (or verify this is a standard topic??) for (DicomSeriesReader::UidFileNamesMap::const_iterator seriesIter = seriesInFiles.begin(); seriesIter != seriesInFiles.end(); ++seriesIter) { StringContainer files = seriesIter->second; DataNode::Pointer node = DicomSeriesReader::LoadDicomSeries( files ); if (node.IsNotNull()) { Image::Pointer image = dynamic_cast( node->GetData() ); result.push_back( image ); } else { } } return result; } std::string mitk::TestDICOMLoading::TypeIDToString(const std::type_info& ti) { if (ti == typeid(unsigned char)) return "UCHAR"; else if (ti == typeid(char)) return "CHAR"; else if (ti == typeid(unsigned short)) return "USHORT"; else if (ti == typeid(short)) return "SHORT"; else if (ti == typeid(unsigned int)) return "UINT"; else if (ti == typeid(int)) return "INT"; else if (ti == typeid(long unsigned int)) return "ULONG"; else if (ti == typeid(long int)) return "LONG"; else if (ti == typeid(float)) return "FLOAT"; else if (ti == typeid(double)) return "DOUBLE"; else return "UNKNOWN"; } // add a line to stringstream result (see DumpImageInformation #define DumpLine(field, data) DumpILine(0, field, data) // add an indented(!) line to stringstream result (see DumpImageInformation #define DumpILine(indent, field, data) \ { \ std::string DumpLine_INDENT; DumpLine_INDENT.resize(indent, ' ' ); \ result << DumpLine_INDENT << field << ": " << data << "\n"; \ } std::string mitk::TestDICOMLoading::DumpImageInformation( const Image* image ) { std::stringstream result; if (image == NULL) return result.str(); SetDefaultLocale(); // basic image data DumpLine( "Pixeltype", TypeIDToString( (image->GetPixelType().GetTypeId()) )); DumpLine( "BitsPerPixel", image->GetPixelType().GetBpe() ); DumpLine( "Dimension", image->GetDimension() ); result << "Dimensions: "; for (unsigned int dim = 0; dim < image->GetDimension(); ++dim) result << image->GetDimension(dim) << " "; result << "\n"; // geometry data result << "Geometry: \n"; Geometry3D* geometry = image->GetGeometry(); if (geometry) { AffineTransform3D* transform = geometry->GetIndexToWorldTransform(); if (transform) { result << " " << "Matrix: "; const AffineTransform3D::MatrixType& matrix = transform->GetMatrix(); for (unsigned int i = 0; i < 3; ++i) for (unsigned int j = 0; j < 3; ++j) result << matrix[i][j] << " "; result << "\n"; result << " " << "Offset: "; const AffineTransform3D::OutputVectorType& offset = transform->GetOffset(); for (unsigned int i = 0; i < 3; ++i) result << offset[i] << " "; result << "\n"; result << " " << "Center: "; const AffineTransform3D::InputPointType& center = transform->GetCenter(); for (unsigned int i = 0; i < 3; ++i) result << center[i] << " "; result << "\n"; result << " " << "Translation: "; const AffineTransform3D::OutputVectorType& translation = transform->GetTranslation(); for (unsigned int i = 0; i < 3; ++i) result << translation[i] << " "; result << "\n"; result << " " << "Scale: "; const double* scale = transform->GetScale(); for (unsigned int i = 0; i < 3; ++i) result << scale[i] << " "; result << "\n"; result << " " << "Origin: "; const Point3D& origin = geometry->GetOrigin(); for (unsigned int i = 0; i < 3; ++i) result << origin[i] << " "; result << "\n"; result << " " << "Spacing: "; const Vector3D& spacing = geometry->GetSpacing(); for (unsigned int i = 0; i < 3; ++i) result << spacing[i] << " "; result << "\n"; result << " " << "TimeBounds: "; const TimeBounds timeBounds = geometry->GetTimeBounds(); for (unsigned int i = 0; i < 2; ++i) result << timeBounds[i] << " "; result << "\n"; } } ResetUserLocale(); return result.str(); } std::string mitk::TestDICOMLoading::trim(const std::string& pString, const std::string& pWhitespace) { const size_t beginStr = pString.find_first_not_of(pWhitespace); if (beginStr == std::string::npos) { // no content return ""; } const size_t endStr = pString.find_last_not_of(pWhitespace); const size_t range = endStr - beginStr + 1; return pString.substr(beginStr, range); } std::string mitk::TestDICOMLoading::reduce(const std::string& pString, const std::string& pFill, const std::string& pWhitespace) { // trim first std::string result(trim(pString, pWhitespace)); // replace sub ranges size_t beginSpace = result.find_first_of(pWhitespace); while (beginSpace != std::string::npos) { const size_t endSpace = result.find_first_not_of(pWhitespace, beginSpace); const size_t range = endSpace - beginSpace; result.replace(beginSpace, range, pFill); const size_t newStart = beginSpace + pFill.length(); beginSpace = result.find_first_of(pWhitespace, newStart); } return result; } bool mitk::TestDICOMLoading::CompareSpacedValueFields( const std::string& reference, const std::string& test, double /*eps*/ ) { bool result(true); // tokenize string, compare each token, if possible by float comparison std::stringstream referenceStream(reduce(reference)); std::stringstream testStream(reduce(test)); std::string refToken; std::string testToken; while ( std::getline( referenceStream, refToken, ' ' ) && std::getline ( testStream, testToken, ' ' ) ) { float refNumber; float testNumber; if ( this->StringToNumber(refToken, refNumber) ) { MITK_DEBUG << "Reference Token '" << refToken << "'" << " value " << refNumber << ", test Token '" << refToken << "'" << " value " << refNumber; if ( this->StringToNumber(testToken, testNumber) ) { result &= ( fabs(refNumber - testNumber) < mitk::eps ); } else { MITK_ERROR << refNumber << " cannot be compared to '" << testToken << "'"; } } else { MITK_DEBUG << "Token '" << refToken << "'" << " handled as string"; result &= refToken == testToken; } } if ( std::getline( referenceStream, refToken, ' ' ) ) { MITK_ERROR << "Reference string still had values when test string was already parsed: ref '" << reference << "', test '" << test << "'"; result = false; } else if ( std::getline( testStream, testToken, ' ' ) ) { MITK_ERROR << "Test string still had values when reference string was already parsed: ref '" << reference << "', test '" << test << "'"; result = false; } return result; } bool mitk::TestDICOMLoading::CompareImageInformationDumps( const std::string& referenceDump, const std::string& testDump ) { KeyValueMap reference = ParseDump(referenceDump); KeyValueMap test = ParseDump(testDump); bool testResult(true); // verify all expected values for (KeyValueMap::const_iterator refIter = reference.begin(); refIter != reference.end(); ++refIter) { const std::string& refKey = refIter->first; const std::string& refValue = refIter->second; if ( test.find(refKey) != test.end() ) { const std::string& testValue = test[refKey]; bool thisTestResult = CompareSpacedValueFields( refValue, testValue ); testResult &= thisTestResult; MITK_DEBUG << refKey << ": '" << refValue << "' == '" << testValue << "' ? " << (thisTestResult?"YES":"NO"); } else { MITK_ERROR << "Reference dump contains a key'" << refKey << "' (value '" << refValue << "')." ; MITK_ERROR << "This key is expected to be generated for tests (but was not). Most probably you need to update your test data."; return false; } } // now check test dump does not contain any additional keys for (KeyValueMap::const_iterator testIter = test.begin(); testIter != test.end(); ++testIter) { const std::string& key = testIter->first; const std::string& value = testIter->second; if ( reference.find(key) == reference.end() ) { MITK_ERROR << "Test dump contains an unexpected key'" << key << "' (value '" << value << "')." ; MITK_ERROR << "This key is not expected. Most probably you need to update your test data."; return false; } } return testResult; } mitk::TestDICOMLoading::KeyValueMap mitk::TestDICOMLoading::ParseDump( const std::string& dump ) { KeyValueMap parsedResult; std::string shredder(dump); std::stack surroundingKeys; std::stack expectedIndents; expectedIndents.push(0); while (true) { std::string::size_type newLinePos = shredder.find( '\n' ); if (newLinePos == std::string::npos || newLinePos == 0) break; std::string line = shredder.substr( 0, newLinePos ); shredder = shredder.erase( 0, newLinePos+1 ); std::string::size_type keyPosition = line.find_first_not_of( ' ' ); std::string::size_type colonPosition = line.find( ':' ); std::string key = line.substr(keyPosition, colonPosition - keyPosition); std::string::size_type firstSpacePosition = key.find_first_of(" "); if (firstSpacePosition != std::string::npos) { key.erase(firstSpacePosition); } if ( keyPosition > expectedIndents.top() ) { // more indent than before expectedIndents.push(keyPosition); } else if (keyPosition == expectedIndents.top() ) { if (!surroundingKeys.empty()) { surroundingKeys.pop(); // last of same length } } else { // less indent than before do expectedIndents.pop(); while (expectedIndents.top() != keyPosition); // unwind until current indent is found } if (!surroundingKeys.empty()) { key = surroundingKeys.top() + "." + key; // construct current key name } surroundingKeys.push(key); // this is the new embracing key std::string value = line.substr(colonPosition+1); MITK_DEBUG << " Key: '" << key << "' value '" << value << "'" ; parsedResult[key] = value; // store parsing result } return parsedResult; } diff --git a/Core/Code/Testing/mitkDicomSeriesReaderTest.cpp b/Core/Code/Testing/mitkDicomSeriesReaderTest.cpp index 9462c9c97c..86dd838e3c 100644 --- a/Core/Code/Testing/mitkDicomSeriesReaderTest.cpp +++ b/Core/Code/Testing/mitkDicomSeriesReaderTest.cpp @@ -1,162 +1,162 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include #include "mitkDicomSeriesReader.h" #include "mitkProperties.h" static std::map > GetTagInformationFromFile(mitk::DicomSeriesReader::StringContainer files) { gdcm::Scanner scanner; std::map > tagInformations; const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location scanner.AddTag( tagSliceLocation ); const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number scanner.AddTag( tagInstanceNumber ); const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number scanner.AddTag( tagSOPInstanceNumber ); //unsigned int slice(0); scanner.Scan(files); // return const_cast(scanner.GetMappings()); gdcm::Scanner::MappingType& tagValueMappings = const_cast(scanner.GetMappings()); for(std::vector::const_iterator fIter = files.begin(); fIter != files.end(); ++fIter) { std::map tags; tags.insert(std::pair (tagSliceLocation, tagValueMappings[fIter->c_str()][tagSliceLocation])); tags.insert(std::pair (tagInstanceNumber, tagValueMappings[fIter->c_str()][tagInstanceNumber])); tags.insert(std::pair (tagSOPInstanceNumber, tagValueMappings[fIter->c_str()][tagSOPInstanceNumber])); tagInformations.insert(std::pair > (fIter->c_str(), tags)); } return tagInformations; } int mitkDicomSeriesReaderTest(int argc, char* argv[]) { // always start with this! MITK_TEST_BEGIN("DicomSeriesReader") if(argc < 1) { MITK_ERROR << "No directory given!"; return -1; } char* dir; dir = argv[1]; //check if DICOMTags have been set as property for mitk::Image - mitk::DicomSeriesReader::UidFileNamesMap seriesInFiles = mitk::DicomSeriesReader::GetSeries( dir ); + mitk::DicomSeriesReader::UidFileNamesMap seriesInFiles = mitk::DicomSeriesReader::GetSeries( dir, true ); std::list images; std::map fileMap; // TODO sort series UIDs, implementation of map iterator might differ on different platforms (or verify this is a standard topic??) for (mitk::DicomSeriesReader::UidFileNamesMap::const_iterator seriesIter = seriesInFiles.begin(); seriesIter != seriesInFiles.end(); ++seriesIter) { mitk::DicomSeriesReader::StringContainer files = seriesIter->second; mitk::DataNode::Pointer node = mitk::DicomSeriesReader::LoadDicomSeries( files ); MITK_TEST_CONDITION_REQUIRED(node.IsNotNull(),"Testing node") if (node.IsNotNull()) { mitk::Image::Pointer image = dynamic_cast( node->GetData() ); images.push_back( image ); fileMap.insert( std::pair(image,files)); } } //Test if DICOM tags have been added correctly to the mitk::image properties const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number for ( std::list::const_iterator imageIter = images.begin(); imageIter != images.end(); ++imageIter ) { const mitk::Image::Pointer image = *imageIter; //Get tag information for all dicom files of this image std::map > tagInformations = GetTagInformationFromFile((*fileMap.find(image)).second); mitk::StringLookupTableProperty* sliceLocation = dynamic_cast(image->GetProperty("dicom.image.0020.1041").GetPointer()); mitk::StringLookupTableProperty* instanceNumber = dynamic_cast(image->GetProperty("dicom.image.0020.0013").GetPointer()); mitk::StringLookupTableProperty* SOPInstnaceNumber = dynamic_cast(image->GetProperty("dicom.image.0008.0018").GetPointer()); mitk::StringLookupTableProperty* files = dynamic_cast(image->GetProperty("files").GetPointer()); MITK_TEST_CONDITION(sliceLocation != NULL, "Test if tag for slice location has been set to mitk image"); if(sliceLocation != NULL) { for(int i = 0; i < (int)sliceLocation->GetValue().GetLookupTable().size(); i++) { if(i < (int)files->GetValue().GetLookupTable().size()) { MITK_INFO << "Table value: " << sliceLocation->GetValue().GetTableValue(i) << " and File value: " << tagInformations[files->GetValue().GetTableValue(i).c_str()][tagSliceLocation] << std::endl; MITK_INFO << "Filename: " << files->GetValue().GetTableValue(i).c_str() << std::endl; MITK_TEST_CONDITION(sliceLocation->GetValue().GetTableValue(i) == tagInformations[files->GetValue().GetTableValue(i).c_str()][tagSliceLocation], "Test if value for slice location is correct"); } } } MITK_TEST_CONDITION(instanceNumber != NULL, "Test if tag for image instance number has been set to mitk image"); if(instanceNumber != NULL) { for(int i = 0; i < (int)instanceNumber->GetValue().GetLookupTable().size(); i++) { if(i < (int)files->GetValue().GetLookupTable().size()) { MITK_INFO << "Table value: " << instanceNumber->GetValue().GetTableValue(i) << " and File value: " << tagInformations[files->GetValue().GetTableValue(i).c_str()][tagInstanceNumber] << std::endl; MITK_INFO << "Filename: " << files->GetValue().GetTableValue(i).c_str() << std::endl; MITK_TEST_CONDITION(instanceNumber->GetValue().GetTableValue(i) == tagInformations[files->GetValue().GetTableValue(i).c_str()][tagInstanceNumber], "Test if value for instance number is correct"); } } } MITK_TEST_CONDITION(SOPInstnaceNumber != NULL, "Test if tag for SOP instance number has been set to mitk image"); if(SOPInstnaceNumber != NULL) { for(int i = 0; i < (int)SOPInstnaceNumber->GetValue().GetLookupTable().size(); i++) { if(i < (int)files->GetValue().GetLookupTable().size()) { MITK_INFO << "Table value: " << instanceNumber->GetValue().GetTableValue(i) << " and File value: " << tagInformations[files->GetValue().GetTableValue(i).c_str()][tagSOPInstanceNumber] << std::endl; MITK_INFO << "Filename: " << files->GetValue().GetTableValue(i).c_str() << std::endl; MITK_TEST_CONDITION(SOPInstnaceNumber->GetValue().GetTableValue(i) == tagInformations[files->GetValue().GetTableValue(i).c_str()][tagSOPInstanceNumber], "Test if value for SOP instance number is correct"); } } } } MITK_TEST_END() }