diff --git a/Examples/mitkdump/mitkdump.cpp b/Examples/mitkdump/mitkdump.cpp index 370d04a33b..24a4cfb213 100644 --- a/Examples/mitkdump/mitkdump.cpp +++ b/Examples/mitkdump/mitkdump.cpp @@ -1,86 +1,91 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDICOMITKSeriesGDCMReader.h" #include "mitkDICOMTagBasedSorter.h" #include "mitkDICOMSortByTag.h" int main(int argc, char* argv[]) { mitk::StringList inputFiles; // TODO for (int a = 1; a < argc; ++a) { inputFiles.push_back( std::string(argv[a]) ); } // ----------------- Configure reader ------------------- mitk::DICOMITKSeriesGDCMReader::Pointer gdcmReader = mitk::DICOMITKSeriesGDCMReader::New(); mitk::DICOMTagBasedSorter::Pointer tagSorter = mitk::DICOMTagBasedSorter::New(); // all the things that split by tag in DicomSeriesReader tagSorter->AddDistinguishingTag( std::make_pair(0x0028, 0x0010) ); // Number of Rows tagSorter->AddDistinguishingTag( std::make_pair(0x0028, 0x0011) ); // Number of Columns tagSorter->AddDistinguishingTag( std::make_pair(0x0028, 0x0030) ); // Pixel Spacing tagSorter->AddDistinguishingTag( std::make_pair(0x0018, 0x1164) ); // Imager Pixel Spacing tagSorter->AddDistinguishingTag( std::make_pair(0x0020, 0x0037) ); // Image Orientation (Patient) // TODO add tolerance parameter (l. 1572 of original code) tagSorter->AddDistinguishingTag( std::make_pair(0x0020, 0x000e) ); // Series Instance UID tagSorter->AddDistinguishingTag( std::make_pair(0x0018, 0x0050) ); // Slice Thickness tagSorter->AddDistinguishingTag( std::make_pair(0x0028, 0x0008) ); // Number of Frames + tagSorter->AddDistinguishingTag( std::make_pair(0x0020, 0x0052) ); // Frame of Reference UID // a sorter... // TODO ugly syntax, improve.. mitk::DICOMSortCriterion::ConstPointer sorting = mitk::DICOMSortByTag::New( std::make_pair(0x0020, 0x0013), // instance number mitk::DICOMSortByTag::New( std::make_pair(0x0020, 0x0012), // aqcuisition number mitk::DICOMSortByTag::New( std::make_pair(0x0008, 0x0032), // aqcuisition time mitk::DICOMSortByTag::New( std::make_pair(0x0018, 0x1060), // trigger time mitk::DICOMSortByTag::New( std::make_pair(0x0008, 0x0018) // SOP instance UID (last resort, not really meaningful but decides clearly) ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer(); tagSorter->SetSortCriterion( sorting ); gdcmReader->AddSortingElement( tagSorter ); // ----------------- Load ------------------- gdcmReader->SetInputFiles( inputFiles ); MITK_INFO << "Analyzing " << inputFiles.size() << " file ..."; gdcmReader->AnalyzeInputFiles(); + gdcmReader->PrintOutputs(std::cout, true); MITK_INFO << "Loading " << inputFiles.size() << " file ..."; gdcmReader->LoadImages(); unsigned int numberOfOutputs = gdcmReader->GetNumberOfOutputs(); for (unsigned int o = 0; o < numberOfOutputs; ++o) { const mitk::DICOMImageBlockDescriptor block = gdcmReader->GetOutput(o); const mitk::DICOMImageFrameList& outputFiles = block.GetImageFrameList(); mitk::Image::Pointer mitkImage = block.GetMitkImage(); MITK_INFO << "-------------------------------------------"; MITK_INFO << "Output " << o << " at " << (void*) mitkImage.GetPointer(); MITK_INFO << " Number of files: " << outputFiles.size(); - MITK_INFO << " Dimensions: " << mitkImage->GetDimension(0) << " " << mitkImage->GetDimension(1) << " " << mitkImage->GetDimension(2); + if (mitkImage.IsNotNull()) + { + MITK_INFO << " Dimensions: " << mitkImage->GetDimension(0) << " " << mitkImage->GetDimension(1) << " " << mitkImage->GetDimension(2); + } } } diff --git a/Modules/DICOMReader/Testing/mitkDICOMITKSeriesGDCMReaderBasicsTest.cpp b/Modules/DICOMReader/Testing/mitkDICOMITKSeriesGDCMReaderBasicsTest.cpp index d4d0183c2f..69869d5941 100644 --- a/Modules/DICOMReader/Testing/mitkDICOMITKSeriesGDCMReaderBasicsTest.cpp +++ b/Modules/DICOMReader/Testing/mitkDICOMITKSeriesGDCMReaderBasicsTest.cpp @@ -1,85 +1,86 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDICOMITKSeriesGDCMReader.h" #include "mitkDICOMFileReaderTestHelper.h" #include "mitkDICOMFilenameSorter.h" #include "mitkDICOMTagBasedSorter.h" #include "mitkDICOMSortByTag.h" #include "mitkTestingMacros.h" int mitkDICOMITKSeriesGDCMReaderBasicsTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkDICOMITKSeriesGDCMReaderBasicsTest"); mitk::DICOMITKSeriesGDCMReader::Pointer gdcmReader = mitk::DICOMITKSeriesGDCMReader::New(); MITK_TEST_CONDITION_REQUIRED(gdcmReader.IsNotNull(), "DICOMITKSeriesGDCMReader can be instantiated."); mitk::DICOMFileReaderTestHelper::SetTestInputFilenames( argc,argv ); // check the Set/GetInput function mitk::DICOMFileReaderTestHelper::TestInputFilenames( gdcmReader ); // check that output is a good reproduction of input (no duplicates, no new elements) mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); // repeat test with filename based sorter in-between mitk::DICOMFilenameSorter::Pointer filenameSorter = mitk::DICOMFilenameSorter::New(); gdcmReader->AddSortingElement( filenameSorter ); mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); // repeat test with some more realistic sorting gdcmReader = mitk::DICOMITKSeriesGDCMReader::New(); // this also tests destruction mitk::DICOMTagBasedSorter::Pointer tagSorter = mitk::DICOMTagBasedSorter::New(); // all the things that split by tag in DicomSeriesReader tagSorter->AddDistinguishingTag( std::make_pair(0x0028, 0x0010) ); // Number of Rows tagSorter->AddDistinguishingTag( std::make_pair(0x0028, 0x0011) ); // Number of Columns tagSorter->AddDistinguishingTag( std::make_pair(0x0028, 0x0030) ); // Pixel Spacing tagSorter->AddDistinguishingTag( std::make_pair(0x0018, 0x1164) ); // Imager Pixel Spacing tagSorter->AddDistinguishingTag( std::make_pair(0x0020, 0x0037) ); // Image Orientation (Patient) // TODO add tolerance parameter (l. 1572 of original code) // TODO handle as real vectors! cluster with configurable errors! tagSorter->AddDistinguishingTag( std::make_pair(0x0020, 0x000e) ); // Series Instance UID tagSorter->AddDistinguishingTag( std::make_pair(0x0018, 0x0050) ); // Slice Thickness tagSorter->AddDistinguishingTag( std::make_pair(0x0028, 0x0008) ); // Number of Frames + tagSorter->AddDistinguishingTag( std::make_pair(0x0020, 0x0052) ); // Frame of Reference UID // a sorter... // TODO ugly syntax, improve.. mitk::DICOMSortCriterion::ConstPointer sorting = mitk::DICOMSortByTag::New( std::make_pair(0x0020, 0x0013), // instance number mitk::DICOMSortByTag::New( std::make_pair(0x0020, 0x0012), // aqcuisition number mitk::DICOMSortByTag::New( std::make_pair(0x0008, 0x0032), // aqcuisition time mitk::DICOMSortByTag::New( std::make_pair(0x0018, 0x1060), // trigger time mitk::DICOMSortByTag::New( std::make_pair(0x0008, 0x0018) // SOP instance UID (last resort, not really meaningful but decides clearly) ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer(); tagSorter->SetSortCriterion( sorting ); gdcmReader->AddSortingElement( tagSorter ); mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); gdcmReader->PrintOutputs(std::cout, true); // really load images mitk::DICOMFileReaderTestHelper::TestMitkImagesAreLoaded( gdcmReader ); MITK_TEST_END(); } diff --git a/Modules/DICOMReader/files.cmake b/Modules/DICOMReader/files.cmake index c79b2fb116..4c2aee0fef 100644 --- a/Modules/DICOMReader/files.cmake +++ b/Modules/DICOMReader/files.cmake @@ -1,27 +1,30 @@ set(H_FILES mitkDICOMFileReader.h mitkDICOMImageFrameInfo.h mitkDICOMImageBlockDescriptor.h mitkDICOMGDCMImageFrameInfo.h mitkDICOMITKSeriesGDCMReader.h mitkDICOMDatasetSorter.h mitkDICOMFilenameSorter.h mitkDICOMEnums.h mitkDICOMTagBasedSorter.h mitkDICOMSortCriterion.h mitkDICOMSortByTag.h + mitkEquiDistantBlocksSorter.h ) set(CPP_FILES mitkDICOMFileReader.cpp mitkDICOMImageBlockDescriptor.cpp mitkDICOMITKSeriesGDCMReader.cpp mitkDICOMDatasetSorter.cpp mitkDICOMFilenameSorter.cpp mitkDICOMTagBasedSorter.cpp mitkDICOMGDCMImageFrameInfo.cpp mitkDICOMImageFrameInfo.cpp mitkDICOMSortCriterion.cpp mitkDICOMSortByTag.cpp mitkITKDICOMSeriesReaderHelper.cpp + mitkEquiDistantBlocksSorter.cpp + mitkGantryTiltInformation.cpp ) diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp index 6241db97d1..73792f05b7 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp @@ -1,272 +1,283 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDICOMITKSeriesGDCMReader.h" #include "mitkITKDICOMSeriesReaderHelper.h" +#include "mitkEquiDistantBlocksSorter.h" #include #include mitk::DICOMITKSeriesGDCMReader ::DICOMITKSeriesGDCMReader() :DICOMFileReader() ,m_FixTiltByShearing(false) ,m_Group3DplusT(false) { } mitk::DICOMITKSeriesGDCMReader ::DICOMITKSeriesGDCMReader(const DICOMITKSeriesGDCMReader& other ) :DICOMFileReader(other) ,m_FixTiltByShearing(false) ,m_Group3DplusT(false) { } mitk::DICOMITKSeriesGDCMReader ::~DICOMITKSeriesGDCMReader() { } mitk::DICOMITKSeriesGDCMReader& mitk::DICOMITKSeriesGDCMReader ::operator=(const DICOMITKSeriesGDCMReader& other) { if (this != &other) { DICOMFileReader::operator=(other); } return *this; } mitk::DICOMGDCMImageFrameList mitk::DICOMITKSeriesGDCMReader ::FromDICOMDatasetList(DICOMDatasetList& input) { DICOMGDCMImageFrameList output; output.reserve(input.size()); for(DICOMDatasetList::iterator inputIter = input.begin(); inputIter != input.end(); ++inputIter) { DICOMGDCMImageFrameInfo* gfi = dynamic_cast(*inputIter); assert(gfi); output.push_back(gfi); } return output; } mitk::DICOMDatasetList mitk::DICOMITKSeriesGDCMReader ::ToDICOMDatasetList(DICOMGDCMImageFrameList& input) { DICOMDatasetList output; output.reserve(input.size()); for(DICOMGDCMImageFrameList::iterator inputIter = input.begin(); inputIter != input.end(); ++inputIter) { DICOMDatasetAccess* da = inputIter->GetPointer(); assert(da); output.push_back(da); } return output; } mitk::DICOMImageFrameList mitk::DICOMITKSeriesGDCMReader ::ToDICOMImageFrameList(DICOMGDCMImageFrameList& input) { DICOMImageFrameList output; output.reserve(input.size()); for(DICOMGDCMImageFrameList::iterator inputIter = input.begin(); inputIter != input.end(); ++inputIter) { DICOMImageFrameInfo::Pointer fi = (*inputIter)->GetFrameInfo(); assert(fi.IsNotNull()); output.push_back(fi); } return output; } void mitk::DICOMITKSeriesGDCMReader ::AnalyzeInputFiles() { // TODO at this point, make sure we have a sorting element at the end that splits geometrically separate blocks + this->EnsureMandatorySortersArePresent(); + itk::TimeProbesCollectorBase timer; timer.Start("Reset"); this->ClearOutputs(); timer.Stop("Reset"); // prepare initial sorting (== list of input files) StringList inputFilenames = this->GetInputFiles(); // scan files for sorting-relevant tags // TODO provide tagToValueMappings to items in initialFramelist / m_SortingResultInProgress timer.Start("Setup scanning"); gdcm::Scanner gdcmScanner; for(SorterList::iterator sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sorterIter) { assert(sorterIter->IsNotNull()); DICOMTagList tags = (*sorterIter)->GetTagsOfInterest(); for(DICOMTagList::const_iterator tagIter = tags.begin(); tagIter != tags.end(); ++tagIter) { MITK_DEBUG << "Sorting uses tag " << tagIter->first << "," << tagIter->second; gdcmScanner.AddTag( gdcm::Tag(tagIter->first, tagIter->second) ); } } timer.Stop("Setup scanning"); timer.Start("Tag scanning"); gdcmScanner.Scan( inputFilenames ); timer.Stop("Tag scanning"); timer.Start("Setup sorting"); DICOMGDCMImageFrameList initialFramelist; for (StringList::const_iterator inputIter = inputFilenames.begin(); inputIter != inputFilenames.end(); ++inputIter) { // TODO check DICOMness and non-multi-framedness initialFramelist.push_back( DICOMGDCMImageFrameInfo::New( DICOMImageFrameInfo::New(*inputIter, 0), gdcmScanner.GetMapping(inputIter->c_str()) ) ); } m_SortingResultInProgress.clear(); m_SortingResultInProgress.push_back( initialFramelist ); timer.Stop("Setup sorting"); // sort and split blocks as configured timer.Start("Sorting frames"); SortingBlockList nextStepSorting; // we should not modify our input list while processing it unsigned int sorterIndex = 0; for(SorterList::iterator sorterIter = m_Sorter.begin(); sorterIter != m_Sorter.end(); ++sorterIndex, ++sorterIter) { std::stringstream ss; ss << "Sorting step " << sorterIndex; timer.Start( ss.str().c_str() ); nextStepSorting.clear(); DICOMDatasetSorter::Pointer& sorter = *sorterIter; for(SortingBlockList::iterator blockIter = m_SortingResultInProgress.begin(); blockIter != m_SortingResultInProgress.end(); ++blockIter) { DICOMGDCMImageFrameList& gdcmInfoFrameList = *blockIter; DICOMDatasetList datasetList = ToDICOMDatasetList( gdcmInfoFrameList ); sorter->SetInput(datasetList); sorter->Sort(); unsigned int numberOfResultingBlocks = sorter->GetNumberOfOutputs(); for (unsigned int b = 0; b < numberOfResultingBlocks; ++b) { DICOMDatasetList blockResult = sorter->GetOutput(b); DICOMGDCMImageFrameList sortedGdcmInfoFrameList = FromDICOMDatasetList(blockResult); nextStepSorting.push_back( sortedGdcmInfoFrameList ); } } m_SortingResultInProgress = nextStepSorting; timer.Stop( ss.str().c_str() ); } timer.Stop("Sorting frames"); // provide final result as output timer.Start("Output"); this->SetNumberOfOutputs( m_SortingResultInProgress.size() ); unsigned int o = 0; for (SortingBlockList::iterator blockIter = m_SortingResultInProgress.begin(); blockIter != m_SortingResultInProgress.end(); ++o, ++blockIter) { DICOMGDCMImageFrameList& gdcmFrameInfoList = *blockIter; DICOMImageFrameList frameList = ToDICOMImageFrameList( gdcmFrameInfoList ); DICOMImageBlockDescriptor block; block.SetImageFrameList( frameList ); this->SetOutput( o, block ); } timer.Stop("Output"); std::cout << "---------------------------------------------------------------" << std::endl; timer.Report( std::cout ); std::cout << "---------------------------------------------------------------" << std::endl; } // void AllocateOutputImages(); bool mitk::DICOMITKSeriesGDCMReader ::LoadImages() { mitk::ITKDICOMSeriesReaderHelper helper; unsigned int numberOfOutputs = this->GetNumberOfOutputs(); for (unsigned int o = 0; o < numberOfOutputs; ++o) { DICOMImageBlockDescriptor& block = this->InternalGetOutput(o); const DICOMImageFrameList& frames = block.GetImageFrameList(); bool correctTilt = false; // TODO bool tiltInfo = false; // TODO ITKDICOMSeriesReaderHelper::StringContainer filenames; for (DICOMImageFrameList::const_iterator frameIter = frames.begin(); frameIter != frames.end(); ++frameIter) { filenames.push_back( (*frameIter)->Filename ); } mitk::Image::Pointer mitkImage = helper.Load( filenames, correctTilt, tiltInfo ); // TODO preloaded images, caching..? block.SetMitkImage( mitkImage ); } return true; } bool mitk::DICOMITKSeriesGDCMReader ::CanHandleFile(const std::string& itkNotUsed(filename)) { return true; // can handle all } void mitk::DICOMITKSeriesGDCMReader ::AddSortingElement(DICOMDatasetSorter* sorter) { assert(sorter); m_Sorter.push_back( sorter ); } + +void +mitk::DICOMITKSeriesGDCMReader +::EnsureMandatorySortersArePresent() +{ + mitk::EquiDistantBlocksSorter::Pointer distanceKeeper = mitk::EquiDistantBlocksSorter::New(); + this->AddSortingElement( distanceKeeper ); +} diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h index 96bf0f4195..ac06642e79 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h @@ -1,73 +1,75 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDICOMITKSeriesGDCMReader_h #define mitkDICOMITKSeriesGDCMReader_h #include "mitkDICOMFileReader.h" #include "mitkDICOMDatasetSorter.h" #include "mitkDICOMGDCMImageFrameInfo.h" #include "DICOMReaderExports.h" namespace mitk { class DICOMReader_EXPORT DICOMITKSeriesGDCMReader : public DICOMFileReader { public: mitkClassMacro( DICOMITKSeriesGDCMReader, DICOMFileReader ); mitkCloneMacro( DICOMITKSeriesGDCMReader ); itkNewMacro( DICOMITKSeriesGDCMReader ); virtual void AnalyzeInputFiles(); // void AllocateOutputImages(); virtual bool LoadImages(); virtual bool CanHandleFile(const std::string& filename); virtual void AddSortingElement(DICOMDatasetSorter* sorter); protected: DICOMITKSeriesGDCMReader(); virtual ~DICOMITKSeriesGDCMReader(); DICOMITKSeriesGDCMReader(const DICOMITKSeriesGDCMReader& other); DICOMITKSeriesGDCMReader& operator=(const DICOMITKSeriesGDCMReader& other); DICOMDatasetList ToDICOMDatasetList(DICOMGDCMImageFrameList& input); DICOMGDCMImageFrameList FromDICOMDatasetList(DICOMDatasetList& input); DICOMImageFrameList ToDICOMImageFrameList(DICOMGDCMImageFrameList& input); + void EnsureMandatorySortersArePresent(); + private: bool m_FixTiltByShearing; bool m_Group3DplusT; typedef std::list SortingBlockList; SortingBlockList m_SortingResultInProgress; typedef std::list SorterList; SorterList m_Sorter; }; } #endif diff --git a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp new file mode 100644 index 0000000000..d6c4fe9422 --- /dev/null +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp @@ -0,0 +1,453 @@ +/*=================================================================== 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 "mitkEquiDistantBlocksSorter.h" +#include "mitkGantryTiltInformation.h" + +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::SliceGroupingAnalysisResult() +:m_GantryTilt(false) +{ +} + +mitk::DICOMDatasetList +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::GetBlockFilenames() +{ + return m_GroupedFiles; +} + +mitk::DICOMDatasetList +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::GetUnsortedFilenames() +{ + return m_UnsortedFiles; +} + +bool +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::ContainsGantryTilt() +{ + return m_GantryTilt; +} + +void +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::AddFileToSortedBlock(DICOMDatasetAccess* dataset) +{ + m_GroupedFiles.push_back( dataset ); +} + +void +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::AddFileToUnsortedBlock(DICOMDatasetAccess* dataset) +{ + m_UnsortedFiles.push_back( dataset ); +} + +void +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::AddFilesToUnsortedBlock(const DICOMDatasetList& datasets) +{ + m_UnsortedFiles.insert( m_UnsortedFiles.end(), datasets.begin(), datasets.end() ); +} + +void +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::FlagGantryTilt() +{ + m_GantryTilt = true; +} + +void +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::UndoPrematureGrouping() +{ + assert( !m_GroupedFiles.empty() ); + m_UnsortedFiles.insert( m_UnsortedFiles.begin(), m_GroupedFiles.back() ); + m_GroupedFiles.pop_back(); +} + +// ------------------------ end helper class + +mitk::EquiDistantBlocksSorter +::EquiDistantBlocksSorter() +:DICOMDatasetSorter() +{ +} + +mitk::EquiDistantBlocksSorter +::~EquiDistantBlocksSorter() +{ +} + +mitk::EquiDistantBlocksSorter +::EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other ) +:DICOMDatasetSorter(other) +{ +} + +mitk::EquiDistantBlocksSorter& +mitk::EquiDistantBlocksSorter +::operator=(const EquiDistantBlocksSorter& other) +{ + if (this != &other) + { + DICOMDatasetSorter::operator=(other); + } + return *this; +} + +mitk::DICOMTagList +mitk::EquiDistantBlocksSorter +::GetTagsOfInterest() +{ + DICOMTagList tags; + tags.push_back( std::make_pair(0x0020, 0x0032) ); // ImagePositionPatient + tags.push_back( std::make_pair(0x0020, 0x0037) ); // ImageOrientationPatient + + return tags; +} + +void +mitk::EquiDistantBlocksSorter +::Sort() +{ + DICOMDatasetList remainingInput = GetInput(); // copy + + bool acceptGantryTilt = false; // TODO + + typedef std::list OutputListType; + OutputListType outputs; + + while (!remainingInput.empty()) // repeat until all files are grouped somehow + { + SliceGroupingAnalysisResult regularBlock = this->AnalyzeFileForITKImageSeriesReaderSpacingAssumption( remainingInput, acceptGantryTilt ); + + outputs.push_back( regularBlock.GetBlockFilenames() ); + remainingInput = regularBlock.GetUnsortedFilenames(); + } + + unsigned int numberOfOutputs = outputs.size(); + this->SetNumberOfOutputs(numberOfOutputs); + + unsigned int outputIndex(0); + for (OutputListType::iterator oIter = outputs.begin(); + oIter != outputs.end(); + ++outputIndex, ++oIter) + { + this->SetOutput(outputIndex, *oIter); + } +} + +std::string +mitk::EquiDistantBlocksSorter +::ConstCharStarToString(const char* s) +{ + return s ? std::string(s) : std::string(); +} + +mitk::Point3D +mitk::EquiDistantBlocksSorter +::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, '\\' ) && dim < 3) + { + p[dim++]= atof(coordinate.c_str()); + } + + if (dim && dim != 3) + { + successful = false; + MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0032). Found " << dim << " instead of 3 values."; + } + else if (dim == 0) + { + successful = false; + p.Fill(0.0); // assume default (0,0,0) + } + + return p; +} + + +void +mitk::EquiDistantBlocksSorter +::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, '\\' ) && dim < 6 ) + { + if (dim<3) + { + right[dim++] = atof(coordinate.c_str()); + } + else + { + up[dim++ - 3] = atof(coordinate.c_str()); + } + } + + if (dim && dim != 6) + { + successful = false; + MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0037). Found " << dim << " instead of 6 values."; + } + else if (dim == 0) + { + // fill with defaults + right.Fill(0.0); + right[0] = 1.0; + + up.Fill(0.0); + up[1] = 1.0; + + successful = false; + } +} + + +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +mitk::EquiDistantBlocksSorter +::AnalyzeFileForITKImageSeriesReaderSpacingAssumption( + const DICOMDatasetList& datasets, + bool groupImagesWithGantryTilt) +{ + // result.first = files that fit ITK's assumption + // result.second = files that do not fit, should be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption() again + SliceGroupingAnalysisResult result; + + // we const_cast here, because I could not use a map.at(), which would make the code much more readable + const DICOMTag tagImagePositionPatient = std::make_pair(0x0020,0x0032); // Image Position (Patient) + const DICOMTag tagImageOrientation = std::make_pair(0x0020, 0x0037); // Image Orientation + const DICOMTag tagGantryTilt = std::make_pair(0x0018, 0x1120); // gantry tilt + + Vector3D fromFirstToSecondOrigin; fromFirstToSecondOrigin.Fill(0.0); + bool fromFirstToSecondOriginInitialized(false); + Point3D thisOrigin; + thisOrigin.Fill(0.0f); + Point3D lastOrigin; + lastOrigin.Fill(0.0f); + Point3D lastDifferentOrigin; + lastDifferentOrigin.Fill(0.0f); + + bool lastOriginInitialized(false); + + MITK_DEBUG << "--------------------------------------------------------------------------------"; + MITK_DEBUG << "Analyzing files for z-spacing assumption of ITK's ImageSeriesReader (group tilted: " << groupImagesWithGantryTilt << ")"; + unsigned int fileIndex(0); + for (DICOMDatasetList::const_iterator fileIter = datasets.begin(); + fileIter != datasets.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 + thisOriginString = (*fileIter)->GetTagValueAsString( tagImagePositionPatient ); + + if (thisOriginString.empty()) + { + // don't let such files be in a common group. Everything without position information will be loaded as a single slice: + // with standard DICOM files this can happen to: CR, DX, SC + MITK_DEBUG << " ==> Sort away " << *fileIter << " for later analysis (no position information)"; // we already have one occupying this position + + if ( result.GetBlockFilenames().empty() ) // nothing WITH position information yet + { + // ==> this is a group of its own, stop processing, come back later + result.AddFileToSortedBlock( *fileIter ); // TODO + + // TODO change to DICOMDatasetList! + DICOMDatasetList remainingFiles; + remainingFiles.insert( remainingFiles.end(), fileIter+1, datasets.end() ); + result.AddFilesToUnsortedBlock( remainingFiles ); + + fileFitsIntoPattern = false; + break; // no files anymore + } + else + { + // ==> this does not match, consider later + result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis + fileFitsIntoPattern = false; + continue; // next file + } + } + + bool ignoredConversionError(-42); // hard to get here, no graceful way to react + thisOrigin = DICOMStringToPoint3D( thisOriginString, ignoredConversionError ); + + MITK_DEBUG << " " << fileIndex << " " << *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.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis + fileFitsIntoPattern = false; + } + else + { + if (!fromFirstToSecondOriginInitialized && lastOriginInitialized) // calculate vector as soon as possible when we get a new position + { + fromFirstToSecondOrigin = thisOrigin - lastDifferentOrigin; + fromFirstToSecondOriginInitialized = true; + + // Here we calculate if this slice and the previous one are well aligned, + // i.e. we test if the previous origin is on a line through the current + // origin, directed into the normal direction of the current slice. + + // If this is NOT the case, then we have a data set with a TILTED GANTRY geometry, + // which cannot be simply loaded into a single mitk::Image at the moment. + // For this case, we flag this finding in the result and DicomSeriesReader + // can correct for that later. + + Vector3D right; right.Fill(0.0); + Vector3D up; right.Fill(0.0); // might be down as well, but it is just a name at this point + std::string orientationValue = (*fileIter)->GetTagValueAsString( tagImageOrientation ); + DICOMStringToOrientationVectors( orientationValue, right, up, ignoredConversionError ); + + GantryTiltInformation tiltInfo( lastDifferentOrigin, thisOrigin, right, up, 1 ); + + if ( tiltInfo.IsSheared() ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt + { + /* 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) + + // tell apart gantry tilt from overall skewedness + // sort out irregularly sheared slices, that IS NOT tilting + + if ( groupImagesWithGantryTilt && tiltInfo.IsRegularGantryTilt() ) + { + bool todoIsDone(false); // TODO add GantryTilt reading + // check if this is at least roughly the same angle as recorded in DICOM tags + if ( todoIsDone /*&& tagValueMappings[fileIter->c_str()].find(tagGantryTilt) != tagValueMappings[fileIter->c_str()].end() */) + { + // read value, compare to calculated angle + std::string tiltStr = (*fileIter)->GetTagValueAsString( tagGantryTilt ); + double angle = atof(tiltStr.c_str()); + + MITK_DEBUG << "Comparing recorded tilt angle " << angle << " against calculated value " << tiltInfo.GetTiltAngleInDegrees(); + // TODO we probably want the signs correct, too (that depends: this is just a rough check, nothing serious) + if ( fabs(angle) - tiltInfo.GetTiltAngleInDegrees() > 0.25) + { + result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis + fileFitsIntoPattern = false; + } + else // tilt angle from header is less than 0.25 degrees different from what we calculated, assume this is fine + { + result.FlagGantryTilt(); + result.AddFileToSortedBlock( *fileIter ); // this file is good for current block + fileFitsIntoPattern = true; + } + } + else // we cannot check the calculated tilt angle against the one from the dicom header (so we assume we are right) + { + result.FlagGantryTilt(); + result.AddFileToSortedBlock( *fileIter ); // this file is good for current block + fileFitsIntoPattern = true; + } + } + else // caller does not want tilt compensation OR shearing is more complicated than tilt + { + result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis + fileFitsIntoPattern = false; + } + } + else // not sheared + { + 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 + + /* Optimistic approach: check if any of the remaining slices fits in */ + result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis + fileFitsIntoPattern = false; + } + else + { + result.AddFileToSortedBlock( *fileIter ); // this file is good for current block + fileFitsIntoPattern = true; + } + } + else // this should be the very first slice + { + result.AddFileToSortedBlock( *fileIter ); // this file is good for current block + fileFitsIntoPattern = true; + } + } + + // record current origin for reference in later iterations + if ( !lastOriginInitialized || ( fileFitsIntoPattern && (thisOrigin != lastOrigin) ) ) + { + lastDifferentOrigin = thisOrigin; + } + + lastOrigin = thisOrigin; + lastOriginInitialized = true; + } + + if ( result.ContainsGantryTilt() ) + { + // check here how many files were grouped. + // IF it was only two files AND we assume tiltedness (e.g. save "distance") + // THEN we would want to also split the two previous files (simple) because + // we don't have any reason to assume they belong together + + if ( result.GetBlockFilenames().size() == 2 ) + { + result.UndoPrematureGrouping(); + } + } + + return result; +} diff --git a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h new file mode 100644 index 0000000000..18865065d8 --- /dev/null +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h @@ -0,0 +1,159 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef mitkEquiDistantBlocksSorter_h +#define mitkEquiDistantBlocksSorter_h + +#include "mitkDICOMDatasetSorter.h" +#include "mitkDICOMSortCriterion.h" + +#include "mitkVector.h" + +#include + +namespace mitk +{ + +/** + \brief Split inputs into blocks of equidant slices. + + This kind of splitting is used as a check before loading a DICOM series with ITK ImageSeriesReader. +*/ +class DICOMReader_EXPORT EquiDistantBlocksSorter : public DICOMDatasetSorter +{ + public: + + mitkClassMacro( EquiDistantBlocksSorter, DICOMDatasetSorter ) + itkNewMacro( EquiDistantBlocksSorter ) + + virtual DICOMTagList GetTagsOfInterest(); + + virtual void Sort(); + + 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. + */ + DICOMDatasetList GetBlockFilenames(); + + /** + \brief Remaining files, which could not be grouped. + */ + DICOMDatasetList GetUnsortedFilenames(); + + /** + \brief Wheter or not the grouped result contain a gantry tilt. + */ + bool ContainsGantryTilt(); + + /** + \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. + */ + void AddFileToSortedBlock(DICOMDatasetAccess* dataset); // TODO Dataset! + + /** + \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. + */ + void AddFileToUnsortedBlock(DICOMDatasetAccess* dataset); + void AddFilesToUnsortedBlock(const DICOMDatasetList& datasets); + + /** + \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. + \todo Could make sense to enhance this with an instance of GantryTiltInformation to store the whole result! + */ + void FlagGantryTilt(); + + /** + \brief Only meaningful for use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption. + */ + void UndoPrematureGrouping(); + + protected: + + DICOMDatasetList m_GroupedFiles; + DICOMDatasetList m_UnsortedFiles; + + bool m_GantryTilt; + }; + + /** + \brief Ensure an equal z-spacing for a group of files. + + Takes as input a number of images, which are all equally oriented and spatially sorted along their normal direction. + + Internally used by GetSeries. Returns two lists: the first one contins slices of equal inter-slice spacing. + The second list contains remaining files, which need to be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption again. + + Relevant code that is matched here is in + itkImageSeriesReader.txx (ImageSeriesReader::GenerateOutputInformation(void)), lines 176 to 245 (as of ITK 3.20) + */ + SliceGroupingAnalysisResult + AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const DICOMDatasetList& files, bool groupsOfSimilarImages); + + /** + \brief Safely convert const char* to std::string. + */ + std::string + ConstCharStarToString(const char* s); + + /** + \brief Convert DICOM string describing a point to Point3D. + + DICOM tags like ImagePositionPatient contain a position as float numbers separated by backslashes: + \verbatim + 42.7131\13.77\0.7 + \endverbatim + */ + Point3D + DICOMStringToPoint3D(const std::string& s, bool& successful); + + /** + \brief Convert DICOM string describing a point two Vector3D. + + DICOM tags like ImageOrientationPatient contain two vectors as float numbers separated by backslashes: + \verbatim + 42.7131\13.77\0.7\137.76\0.3 + \endverbatim + */ + void + DICOMStringToOrientationVectors(const std::string& s, Vector3D& right, Vector3D& up, bool& successful); + + EquiDistantBlocksSorter(); + virtual ~EquiDistantBlocksSorter(); + + EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other); + EquiDistantBlocksSorter& operator=(const EquiDistantBlocksSorter& other); +}; + +} + +#endif diff --git a/Modules/DICOMReader/mitkGantryTiltInformation.cpp b/Modules/DICOMReader/mitkGantryTiltInformation.cpp new file mode 100644 index 0000000000..e0109bbb7b --- /dev/null +++ b/Modules/DICOMReader/mitkGantryTiltInformation.cpp @@ -0,0 +1,185 @@ +#include "mitkGantryTiltInformation.h" + +#include "mitkLogMacros.h" + +mitk::GantryTiltInformation::GantryTiltInformation() +: m_ShiftUp(0.0) +, m_ShiftRight(0.0) +, m_ShiftNormal(0.0) +, m_ITKAssumedSliceSpacing(0.0) +, m_NumberOfSlicesApart(1) +{ +} + + +#define doublepoint(x) \ + Point3Dd x; \ + x[0] = x ## f[0]; \ + x[1] = x ## f[1]; \ + x[2] = x ## f[2]; + + +#define doublevector(x) \ + Vector3Dd x; \ + x[0] = x ## f[0]; \ + x[1] = x ## f[1]; \ + x[2] = x ## f[2]; + +mitk::GantryTiltInformation::GantryTiltInformation( + const Point3D& origin1f, const Point3D& origin2f, + const Vector3D& rightf, const Vector3D& upf, + unsigned int numberOfSlicesApart) +: m_ShiftUp(0.0) +, m_ShiftRight(0.0) +, m_ShiftNormal(0.0) +, m_NumberOfSlicesApart(numberOfSlicesApart) +{ + assert(numberOfSlicesApart); + + doublepoint(origin1); + doublepoint(origin2); + doublevector(right); + doublevector(up); + + // determine if slice 1 (imagePosition1 and imageOrientation1) and slice 2 can be in one orthogonal slice stack: + // calculate a line from origin 1, directed along the normal of slice (calculated as the cross product of orientation 1) + // check if this line passes through origin 2 + + /* + Determine if line (imagePosition2 + l * normal) contains imagePosition1. + Done by calculating the distance of imagePosition1 from line (imagePosition2 + l *normal) + + E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html + + squared distance = | (pointAlongNormal - origin2) x (origin2 - origin1) | ^ 2 + / + |pointAlongNormal - origin2| ^ 2 + + ( x meaning the cross product ) + */ + + Vector3Dd normal = itk::CrossProduct(right, up); + Point3Dd pointAlongNormal = origin2 + normal; + + double numerator = itk::CrossProduct( pointAlongNormal - origin2 , origin2 - origin1 ).GetSquaredNorm(); + double denominator = (pointAlongNormal - origin2).GetSquaredNorm(); + + double distance = sqrt(numerator / denominator); + + if ( distance > 0.001 ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt + { + MITK_DEBUG << " Series seems to contain a tilted (or sheared) geometry"; + MITK_DEBUG << " Distance of expected slice origin from actual slice origin: " << distance; + MITK_DEBUG << " ==> storing this shift for later analysis:"; + MITK_DEBUG << " v right: " << right; + MITK_DEBUG << " v up: " << up; + MITK_DEBUG << " v normal: " << normal; + + Point3Dd projectionRight = projectPointOnLine( origin1, origin2, right ); + Point3Dd projectionNormal = projectPointOnLine( origin1, origin2, normal ); + + m_ShiftRight = (projectionRight - origin2).GetNorm(); + m_ShiftNormal = (projectionNormal - origin2).GetNorm(); + + /* + now also check to which side the image is shifted. + + Calculation e.g. from + http://mathworld.wolfram.com/Point-PlaneDistance.html + */ + + Point3Dd testPoint = origin1; + Vector3Dd planeNormal = up; + + double signedDistance = ( + planeNormal[0] * testPoint[0] + + planeNormal[1] * testPoint[1] + + planeNormal[2] * testPoint[2] + - ( + planeNormal[0] * origin2[0] + + planeNormal[1] * origin2[1] + + planeNormal[2] * origin2[2] + ) + ) + / + sqrt( planeNormal[0] * planeNormal[0] + + planeNormal[1] * planeNormal[1] + + planeNormal[2] * planeNormal[2] + ); + + m_ShiftUp = signedDistance; + + m_ITKAssumedSliceSpacing = (origin2 - origin1).GetNorm(); + // How do we now this is assumed? See header documentation for ITK code references + //double itkAssumedSliceSpacing = sqrt( m_ShiftUp * m_ShiftUp + m_ShiftNormal * m_ShiftNormal ); + + MITK_DEBUG << " shift normal: " << m_ShiftNormal; + MITK_DEBUG << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing; + MITK_DEBUG << " shift up: " << m_ShiftUp; + MITK_DEBUG << " shift right: " << m_ShiftRight; + + MITK_DEBUG << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535; + } +} + +mitk::Point3D +mitk::GantryTiltInformation::projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ) +{ + /** + See illustration at http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage472/ + + vector(lineOrigin,p) = normal * ( innerproduct((p - lineOrigin),normal) / squared-length(normal) ) + */ + + Vector3Dd lineOriginToP = p - lineOrigin; + double innerProduct = lineOriginToP * lineDirection; + + double factor = innerProduct / lineDirection.GetSquaredNorm(); + Point3Dd projection = lineOrigin + factor * lineDirection; + + return projection; +} + +double +mitk::GantryTiltInformation::GetTiltCorrectedAdditionalSize() const +{ + return fabs(m_ShiftUp); +} + +double +mitk::GantryTiltInformation::GetTiltAngleInDegrees() const +{ + return atan( fabs(m_ShiftUp) / m_ShiftNormal ) * 180.0 / 3.1415926535; +} + +double +mitk::GantryTiltInformation::GetMatrixCoefficientForCorrectionInWorldCoordinates() const +{ + // so many mm need to be shifted per slice! + return m_ShiftUp / static_cast(m_NumberOfSlicesApart); +} + +double +mitk::GantryTiltInformation::GetRealZSpacing() const +{ + return m_ShiftNormal / static_cast(m_NumberOfSlicesApart); +} + + +bool +mitk::GantryTiltInformation::IsSheared() const +{ + return ( fabs(m_ShiftRight) > 0.001 + || fabs(m_ShiftUp) > 0.001); +} + + +bool +mitk::GantryTiltInformation::IsRegularGantryTilt() const +{ + return ( fabs(m_ShiftRight) < 0.001 + && fabs(m_ShiftUp) > 0.001); +} + + + diff --git a/Modules/DICOMReader/mitkGantryTiltInformation.h b/Modules/DICOMReader/mitkGantryTiltInformation.h new file mode 100644 index 0000000000..ed9736b45e --- /dev/null +++ b/Modules/DICOMReader/mitkGantryTiltInformation.h @@ -0,0 +1,127 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef mitkGantryTiltInformation_h +#define mitkGantryTiltInformation_h + +#include "mitkVector.h" + +namespace mitk +{ + +/** + \brief Gantry tilt analysis result. + + Takes geometry information for two slices of a DICOM series and + calculates whether these fit into an orthogonal block or not. + If NOT, they can either be the result of an acquisition with + gantry tilt OR completly broken by some shearing transformation. + + Most calculations are done in the constructor, results can then + be read via the remaining methods. + */ +class GantryTiltInformation +{ + public: + + // two types to avoid any rounding errors + typedef itk::Point Point3Dd; + typedef itk::Vector Vector3Dd; + + /** + \brief Just so we can create empty instances for assigning results later. + */ + GantryTiltInformation(); + + /** + \brief THE constructor, which does all the calculations. + + Determining the amount of tilt is done by checking the distances + of origin1 from planes through origin2. Two planes are considered: + - normal vector along normal of slices (right x up): gives the slice distance + - normal vector along orientation vector "up": gives the shift parallel to the plane orientation + + The tilt angle can then be calculated from these distances + + \param origin1 origin of the first slice + \param origin2 origin of the second slice + \param right right/up describe the orientatation of borth slices + \param up right/up describe the orientatation of borth slices + \param numberOfSlicesApart how many slices are the given origins apart (1 for neighboring slices) + */ + GantryTiltInformation( const Point3D& origin1, + const Point3D& origin2, + const Vector3D& right, + const Vector3D& up, unsigned int numberOfSlicesApart); + + /** + \brief Whether the slices were sheared. + + True if any of the shifts along right or up vector are non-zero. + */ + bool IsSheared() const; + + /** + \brief Whether the shearing is a gantry tilt or more complicated. + + Gantry tilt will only produce shifts in ONE orientation, not in both. + + Since the correction code currently only coveres one tilt direction + AND we don't know of medical images with two tilt directions, the + loading code wants to check if our assumptions are true. + */ + bool IsRegularGantryTilt() const; + + /** + \brief The offset distance in Y direction for each slice in mm (describes the tilt result). + */ + double GetMatrixCoefficientForCorrectionInWorldCoordinates() const; + + + /** + \brief The z / inter-slice spacing. Needed to correct ImageSeriesReader's result. + */ + double GetRealZSpacing() const; + + /** + \brief The shift between first and last slice in mm. + + Needed to resize an orthogonal image volume. + */ + double GetTiltCorrectedAdditionalSize() const; + + /** + \brief Calculated tilt angle in degrees. + */ + double GetTiltAngleInDegrees() const; + + protected: + + /** + \brief Projection of point p onto line through lineOrigin in direction of lineDirection. + */ + Point3D projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ); + + double m_ShiftUp; + double m_ShiftRight; + double m_ShiftNormal; + double m_ITKAssumedSliceSpacing; + unsigned int m_NumberOfSlicesApart; +}; + +} // namespace + +#endif diff --git a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp index 4dbad05b3d..db8cdc6926 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp @@ -1,277 +1,97 @@ /*=================================================================== 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 "mitkITKDICOMSeriesReaderHelper.h" #include "mitkITKDICOMSeriesReaderHelper.txx" -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::GantryTiltInformation() -: m_ShiftUp(0.0) -, m_ShiftRight(0.0) -, m_ShiftNormal(0.0) -, m_ITKAssumedSliceSpacing(0.0) -, m_NumberOfSlicesApart(1) -{ -} - - -#define doublepoint(x) \ - Point3Dd x; \ - x[0] = x ## f[0]; \ - x[1] = x ## f[1]; \ - x[2] = x ## f[2]; - - -#define doublevector(x) \ - Vector3Dd x; \ - x[0] = x ## f[0]; \ - x[1] = x ## f[1]; \ - x[2] = x ## f[2]; - -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::GantryTiltInformation( - const Point3D& origin1f, const Point3D& origin2f, - const Vector3D& rightf, const Vector3D& upf, - unsigned int numberOfSlicesApart) -: m_ShiftUp(0.0) -, m_ShiftRight(0.0) -, m_ShiftNormal(0.0) -, m_NumberOfSlicesApart(numberOfSlicesApart) -{ - assert(numberOfSlicesApart); - - doublepoint(origin1); - doublepoint(origin2); - doublevector(right); - doublevector(up); - - // determine if slice 1 (imagePosition1 and imageOrientation1) and slice 2 can be in one orthogonal slice stack: - // calculate a line from origin 1, directed along the normal of slice (calculated as the cross product of orientation 1) - // check if this line passes through origin 2 - - /* - Determine if line (imagePosition2 + l * normal) contains imagePosition1. - Done by calculating the distance of imagePosition1 from line (imagePosition2 + l *normal) - - E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html - - squared distance = | (pointAlongNormal - origin2) x (origin2 - origin1) | ^ 2 - / - |pointAlongNormal - origin2| ^ 2 - - ( x meaning the cross product ) - */ - - Vector3Dd normal = itk::CrossProduct(right, up); - Point3Dd pointAlongNormal = origin2 + normal; - - double numerator = itk::CrossProduct( pointAlongNormal - origin2 , origin2 - origin1 ).GetSquaredNorm(); - double denominator = (pointAlongNormal - origin2).GetSquaredNorm(); - - double distance = sqrt(numerator / denominator); - - if ( distance > 0.001 ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt - { - MITK_DEBUG << " Series seems to contain a tilted (or sheared) geometry"; - MITK_DEBUG << " Distance of expected slice origin from actual slice origin: " << distance; - MITK_DEBUG << " ==> storing this shift for later analysis:"; - MITK_DEBUG << " v right: " << right; - MITK_DEBUG << " v up: " << up; - MITK_DEBUG << " v normal: " << normal; - - Point3Dd projectionRight = projectPointOnLine( origin1, origin2, right ); - Point3Dd projectionNormal = projectPointOnLine( origin1, origin2, normal ); - - m_ShiftRight = (projectionRight - origin2).GetNorm(); - m_ShiftNormal = (projectionNormal - origin2).GetNorm(); - - /* - now also check to which side the image is shifted. - - Calculation e.g. from - http://mathworld.wolfram.com/Point-PlaneDistance.html - */ - - Point3Dd testPoint = origin1; - Vector3Dd planeNormal = up; - - double signedDistance = ( - planeNormal[0] * testPoint[0] - + planeNormal[1] * testPoint[1] - + planeNormal[2] * testPoint[2] - - ( - planeNormal[0] * origin2[0] - + planeNormal[1] * origin2[1] - + planeNormal[2] * origin2[2] - ) - ) - / - sqrt( planeNormal[0] * planeNormal[0] - + planeNormal[1] * planeNormal[1] - + planeNormal[2] * planeNormal[2] - ); - - m_ShiftUp = signedDistance; - - m_ITKAssumedSliceSpacing = (origin2 - origin1).GetNorm(); - // How do we now this is assumed? See header documentation for ITK code references - //double itkAssumedSliceSpacing = sqrt( m_ShiftUp * m_ShiftUp + m_ShiftNormal * m_ShiftNormal ); - - MITK_DEBUG << " shift normal: " << m_ShiftNormal; - MITK_DEBUG << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing; - MITK_DEBUG << " shift up: " << m_ShiftUp; - MITK_DEBUG << " shift right: " << m_ShiftRight; - - MITK_DEBUG << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535; - } -} - -mitk::Point3D -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ) -{ - /** - See illustration at http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage472/ - - vector(lineOrigin,p) = normal * ( innerproduct((p - lineOrigin),normal) / squared-length(normal) ) - */ - - Vector3Dd lineOriginToP = p - lineOrigin; - double innerProduct = lineOriginToP * lineDirection; - - double factor = innerProduct / lineDirection.GetSquaredNorm(); - Point3Dd projection = lineOrigin + factor * lineDirection; - - return projection; -} - -double -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::GetTiltCorrectedAdditionalSize() const -{ - return fabs(m_ShiftUp); -} - -double -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::GetTiltAngleInDegrees() const -{ - return atan( fabs(m_ShiftUp) / m_ShiftNormal ) * 180.0 / 3.1415926535; -} - -double -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::GetMatrixCoefficientForCorrectionInWorldCoordinates() const -{ - // so many mm need to be shifted per slice! - return m_ShiftUp / static_cast(m_NumberOfSlicesApart); -} - -double -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::GetRealZSpacing() const -{ - return m_ShiftNormal / static_cast(m_NumberOfSlicesApart); -} - - -bool -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::IsSheared() const -{ - return ( fabs(m_ShiftRight) > 0.001 - || fabs(m_ShiftUp) > 0.001); -} - - -bool -mitk::ITKDICOMSeriesReaderHelper::GantryTiltInformation::IsRegularGantryTilt() const -{ - return ( fabs(m_ShiftRight) < 0.001 - && fabs(m_ShiftUp) > 0.001); -} - - mitk::Image::Pointer mitk::ITKDICOMSeriesReaderHelper ::Load( const StringContainer& filenames, bool correctTilt, bool itkNotUsed(tiltInfo) ) { if( filenames.empty() ) { MITK_DEBUG << "Calling LoadDicomSeries with empty filename string container. Probably invalid application logic."; // TODO return NULL; // this is not actually an error but the result is very simple } DcmIoType::Pointer io = DcmIoType::New(); GantryTiltInformation tiltInfo; // TODO Image::Pointer preLoadedImageBlock; // TODO try { if (io->CanReadFile(filenames.front().c_str())) { io->SetFileName(filenames.front().c_str()); io->ReadImageInformation(); if (io->GetPixelType() == itk::ImageIOBase::SCALAR) { switch (io->GetComponentType()) { case DcmIoType::UCHAR: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::CHAR: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::USHORT: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::SHORT: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::UINT: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::INT: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::ULONG: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::LONG: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::FLOAT: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::DOUBLE: return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); default: MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); } } else if (io->GetPixelType() == itk::ImageIOBase::RGB) { switch (io->GetComponentType()) { case DcmIoType::UCHAR: return LoadDICOMByITK< itk::RGBPixel >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::CHAR: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::USHORT: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::SHORT: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::UINT: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::INT: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::ULONG: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::LONG: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::FLOAT: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); case DcmIoType::DOUBLE: return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, preLoadedImageBlock); default: MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); } } MITK_ERROR << "Unsupported DICOM pixel type"; return NULL; } } 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 NULL; } diff --git a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h index 296f20b743..a2d5040980 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h @@ -1,158 +1,56 @@ /*=================================================================== 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 mitkDICOMSeriesReaderHelper_h #define mitkDICOMSeriesReaderHelper_h -//#include "mitkVector.h" #include "mitkImage.h" +#include "mitkGantryTiltInformation.h" #include namespace mitk { class ITKDICOMSeriesReaderHelper { public: typedef std::vector StringContainer; typedef itk::GDCMImageIO DcmIoType; // TODO remove, we are NOT flexible here - /** - \brief Gantry tilt analysis result. - - Takes geometry information for two slices of a DICOM series and - calculates whether these fit into an orthogonal block or not. - If NOT, they can either be the result of an acquisition with - gantry tilt OR completly broken by some shearing transformation. - - Most calculations are done in the constructor, results can then - be read via the remaining methods. - */ - class GantryTiltInformation - { - public: - - // two types to avoid any rounding errors - typedef itk::Point Point3Dd; - typedef itk::Vector Vector3Dd; - - /** - \brief Just so we can create empty instances for assigning results later. - */ - GantryTiltInformation(); - - /** - \brief THE constructor, which does all the calculations. - - Determining the amount of tilt is done by checking the distances - of origin1 from planes through origin2. Two planes are considered: - - normal vector along normal of slices (right x up): gives the slice distance - - normal vector along orientation vector "up": gives the shift parallel to the plane orientation - - The tilt angle can then be calculated from these distances - - \param origin1 origin of the first slice - \param origin2 origin of the second slice - \param right right/up describe the orientatation of borth slices - \param up right/up describe the orientatation of borth slices - \param numberOfSlicesApart how many slices are the given origins apart (1 for neighboring slices) - */ - GantryTiltInformation( const Point3D& origin1, - const Point3D& origin2, - const Vector3D& right, - const Vector3D& up, unsigned int numberOfSlicesApart); - - /** - \brief Whether the slices were sheared. - - True if any of the shifts along right or up vector are non-zero. - */ - bool IsSheared() const; - - /** - \brief Whether the shearing is a gantry tilt or more complicated. - - Gantry tilt will only produce shifts in ONE orientation, not in both. - - Since the correction code currently only coveres one tilt direction - AND we don't know of medical images with two tilt directions, the - loading code wants to check if our assumptions are true. - */ - bool IsRegularGantryTilt() const; - - /** - \brief The offset distance in Y direction for each slice in mm (describes the tilt result). - */ - double GetMatrixCoefficientForCorrectionInWorldCoordinates() const; - - - /** - \brief The z / inter-slice spacing. Needed to correct ImageSeriesReader's result. - */ - double GetRealZSpacing() const; - - /** - \brief The shift between first and last slice in mm. - - Needed to resize an orthogonal image volume. - */ - double GetTiltCorrectedAdditionalSize() const; - - /** - \brief Calculated tilt angle in degrees. - */ - double GetTiltAngleInDegrees() const; - - protected: - - /** - \brief Projection of point p onto line through lineOrigin in direction of lineDirection. - */ - Point3D projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ); - - double m_ShiftUp; - double m_ShiftRight; - double m_ShiftNormal; - double m_ITKAssumedSliceSpacing; - unsigned int m_NumberOfSlicesApart; - }; - Image::Pointer Load( const StringContainer& filenames, bool correctTilt, bool tiltInfo ); template typename ImageType::Pointer // TODO this is NOT inplace! InPlaceFixUpTiltedGeometry( ImageType* input, const GantryTiltInformation& tiltInfo ); - - template Image::Pointer LoadDICOMByITK( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, Image::Pointer preLoadedImageBlock ); }; } #endif