diff --git a/Modules/DICOMReader/Testing/mitkDICOMFileReaderTestHelper.h b/Modules/DICOMReader/Testing/mitkDICOMFileReaderTestHelper.h index 5518633900..7368d04f93 100644 --- a/Modules/DICOMReader/Testing/mitkDICOMFileReaderTestHelper.h +++ b/Modules/DICOMReader/Testing/mitkDICOMFileReaderTestHelper.h @@ -1,129 +1,153 @@ /*=================================================================== 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 mitkDICOMFileReaderTestHelper_h #define mitkDICOMFileReaderTestHelper_h #include "mitkDICOMFileReader.h" #include "mitkDICOMEnums.h" #include "mitkTestingMacros.h" namespace mitk { class DICOMFileReaderTestHelper { public: static StringList& GetInputFilenames() { static StringList inputs; return inputs; } static void SetTestInputFilenames(int argc, char* argv[]) { mitk::StringList inputFiles; for (int a = 1; a < argc; ++a) { inputFiles.push_back( argv[a] ); } GetInputFilenames() = inputFiles; } static void SetTestInputFilenames(const StringList& filenames) { GetInputFilenames() = filenames; } static void TestInputFilenames(DICOMFileReader* reader) { StringList inputFiles = GetInputFilenames(); reader->SetInputFiles( inputFiles ); const StringList& inputFilesReturned = reader->GetInputFiles(); MITK_TEST_CONDITION( inputFilesReturned.size() == inputFiles.size(), "Input file list is received") MITK_TEST_CONDITION( reader->GetNumberOfOutputs() == 0, "No outputs without analysis") // TODO: check that strings are actually contained } static void TestOutputsContainInputs(DICOMFileReader* reader) { StringList inputFiles = GetInputFilenames(); reader->SetInputFiles( inputFiles ); reader->AnalyzeInputFiles(); StringList allSortedInputsFiles; unsigned int numberOfOutputs = reader->GetNumberOfOutputs(); for (unsigned int o = 0; o < numberOfOutputs; ++o) { - DICOMImageBlockDescriptor block = reader->GetOutput(o); + const DICOMImageBlockDescriptor block = reader->GetOutput(o); const DICOMImageFrameList& outputFiles = block.GetImageFrameList(); for(DICOMImageFrameList::const_iterator iter = outputFiles.begin(); iter != outputFiles.end(); ++iter) { // check that output is part of input StringList::iterator inputPositionOfCurrentOutput = std::find( inputFiles.begin(), inputFiles.end(), (*iter)->Filename ); if (inputPositionOfCurrentOutput != inputFiles.end()) { // check that output is only part of ONE output StringList::iterator outputPositionOfCurrentOutput = std::find( allSortedInputsFiles.begin(), allSortedInputsFiles.end(), (*iter)->Filename ); if (outputPositionOfCurrentOutput == allSortedInputsFiles.end()) { // was not in list before allSortedInputsFiles.push_back( *inputPositionOfCurrentOutput ); } else { reader->PrintOutputs(std::cout); MITK_TEST_CONDITION_REQUIRED(false, "File '" << (*iter)->Filename << "' appears in TWO outputs. Readers are expected to use each frame only once." ) } } else { reader->PrintOutputs(std::cout); MITK_TEST_CONDITION_REQUIRED(false, "File '" << (*iter)->Filename << "' appears in output, but it was never part of the input list." ) } } } MITK_TEST_CONDITION( allSortedInputsFiles.size() == inputFiles.size(), "Output list size (" << allSortedInputsFiles.size() << ") equals input list size (" << inputFiles.size() << ")" ) try { - DICOMImageBlockDescriptor block = reader->GetOutput( inputFiles.size() ); + const DICOMImageBlockDescriptor block = reader->GetOutput( inputFiles.size() ); MITK_TEST_CONDITION(false, "Invalid indices for GetOutput() should throw exception") } catch( std::invalid_argument& ) { MITK_TEST_CONDITION(true, "Invalid indices for GetOutput() should throw exception") } } +static void TestMitkImagesAreLoaded(DICOMFileReader* reader) +{ + StringList inputFiles = GetInputFilenames(); + reader->SetInputFiles( inputFiles ); + + reader->AnalyzeInputFiles(); + reader->LoadImages(); + + unsigned int numberOfOutputs = reader->GetNumberOfOutputs(); + for (unsigned int o = 0; o < numberOfOutputs; ++o) + { + const DICOMImageBlockDescriptor block = reader->GetOutput(o); + + const 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); + } +} + + }; // end test class } // namespace #endif diff --git a/Modules/DICOMReader/Testing/mitkDICOMITKSeriesGDCMReaderBasicsTest.cpp b/Modules/DICOMReader/Testing/mitkDICOMITKSeriesGDCMReaderBasicsTest.cpp index 4644513886..5574e730d1 100644 --- a/Modules/DICOMReader/Testing/mitkDICOMITKSeriesGDCMReaderBasicsTest.cpp +++ b/Modules/DICOMReader/Testing/mitkDICOMITKSeriesGDCMReaderBasicsTest.cpp @@ -1,78 +1,84 @@ /*=================================================================== 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) 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 // a sorter... // TODO ugly syntax, improve.. - mitk::DICOMSortCriterion::Pointer sorting = + 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(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 7e77d13a28..c79b2fb116 100644 --- a/Modules/DICOMReader/files.cmake +++ b/Modules/DICOMReader/files.cmake @@ -1,26 +1,27 @@ 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 ) 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 ) diff --git a/Modules/DICOMReader/mitkDICOMFileReader.cpp b/Modules/DICOMReader/mitkDICOMFileReader.cpp index 741efa51d7..dd1e4780e2 100644 --- a/Modules/DICOMReader/mitkDICOMFileReader.cpp +++ b/Modules/DICOMReader/mitkDICOMFileReader.cpp @@ -1,148 +1,163 @@ /*=================================================================== 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 "mitkDICOMFileReader.h" mitk::DICOMFileReader ::DICOMFileReader() :itk::LightObject() { } mitk::DICOMFileReader ::~DICOMFileReader() { } mitk::DICOMFileReader ::DICOMFileReader(const DICOMFileReader& other ) :itk::LightObject() ,m_Outputs( other.m_Outputs ) // TODO copy instead of reference! { } mitk::DICOMFileReader& mitk::DICOMFileReader ::operator=(const DICOMFileReader& other) { if (this != &other) { m_InputFilenames = other.m_InputFilenames; m_Outputs = other.m_Outputs; // TODO copy instead of reference! } return *this; } void mitk::DICOMFileReader ::SetInputFiles(StringList filenames) { m_InputFilenames = filenames; } const mitk::StringList& mitk::DICOMFileReader ::GetInputFiles() const { return m_InputFilenames; } unsigned int mitk::DICOMFileReader ::GetNumberOfOutputs() const { return m_Outputs.size(); } void mitk::DICOMFileReader ::ClearOutputs() { m_Outputs.clear(); } void mitk::DICOMFileReader ::SetNumberOfOutputs(unsigned int numberOfOutputs) { m_Outputs.resize(numberOfOutputs); } void mitk::DICOMFileReader ::SetOutput(unsigned int index, const mitk::DICOMImageBlockDescriptor& output) { if (index < m_Outputs.size()) { m_Outputs[index] = output; } else { std::stringstream ss; ss << "Index " << index << " out of range (" << m_Outputs.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } void mitk::DICOMFileReader ::PrintOutputs(std::ostream& os, bool filenameDetails) { os << "---- Outputs of DICOMFilereader " << (void*)this << "----"<< std::endl; for (unsigned int o = 0; o < m_Outputs.size(); ++o) { os << "-- Output " << o << std::endl; const DICOMImageBlockDescriptor& block = m_Outputs[o]; const DICOMImageFrameList& frames = block.GetImageFrameList(); os << " Number of frames: " << frames.size() << std::endl; os << " Pixels interpolated: " << (block.GetPixelsInterpolated() ? "true" : "false") << std::endl; os << " Pixel spacing interpretation: " << (int)block.GetPixelSpacingInterpretation() << std::endl; os << " MITK image: " << (void*)block.GetMitkImage().GetPointer() << std::endl; if (filenameDetails) { for (DICOMImageFrameList::const_iterator frameIter = frames.begin(); frameIter != frames.end(); ++frameIter) { os << " " << (*frameIter)->Filename; if ((*frameIter)->FrameNo > 0) { os << ", " << (*frameIter)->FrameNo; } os << std::endl; } } } os << "---- End of output list ----" << std::endl; } const mitk::DICOMImageBlockDescriptor& mitk::DICOMFileReader ::GetOutput(unsigned int index) const { if (index < m_Outputs.size()) { return m_Outputs[index]; } else { std::stringstream ss; ss << "Index " << index << " out of range (" << m_Outputs.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } +mitk::DICOMImageBlockDescriptor& +mitk::DICOMFileReader +::InternalGetOutput(unsigned int index) +{ + if (index < m_Outputs.size()) + { + return m_Outputs[index]; + } + else + { + std::stringstream ss; + ss << "Index " << index << " out of range (" << m_Outputs.size() << " indices reserved)"; + throw std::invalid_argument( ss.str() ); + } +} diff --git a/Modules/DICOMReader/mitkDICOMFileReader.h b/Modules/DICOMReader/mitkDICOMFileReader.h index 657f2be15b..4699c0ae66 100644 --- a/Modules/DICOMReader/mitkDICOMFileReader.h +++ b/Modules/DICOMReader/mitkDICOMFileReader.h @@ -1,77 +1,79 @@ /*=================================================================== 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 mitkDICOMFileReader_h #define mitkDICOMFileReader_h #include "itkObjectFactory.h" #include "mitkCommon.h" #include "DICOMReaderExports.h" #include "mitkDICOMImageBlockDescriptor.h" namespace mitk { class DICOMReader_EXPORT DICOMFileReader : public itk::LightObject { public: enum LoadingConfidence { NoSupport = 0, FullSupport = 1, PartialSupport = 2, }; mitkClassMacro( DICOMFileReader, itk::LightObject ) void SetInputFiles(StringList filenames); const StringList& GetInputFiles() const; virtual void AnalyzeInputFiles() = 0; unsigned int GetNumberOfOutputs() const; const DICOMImageBlockDescriptor& GetOutput(unsigned int index) const; // void AllocateOutputImages(); virtual bool LoadImages() = 0; virtual bool CanHandleFile(const std::string& filename) = 0; void PrintOutputs(std::ostream& os, bool filenameDetails = false); protected: DICOMFileReader(); virtual ~DICOMFileReader(); DICOMFileReader(const DICOMFileReader& other); DICOMFileReader& operator=(const DICOMFileReader& other); void ClearOutputs(); void SetNumberOfOutputs(unsigned int numberOfOutputs); void SetOutput(unsigned int index, const DICOMImageBlockDescriptor& output); + DICOMImageBlockDescriptor& InternalGetOutput(unsigned int index); + private: StringList m_InputFilenames; std::vector< DICOMImageBlockDescriptor > m_Outputs; }; } #endif diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp index e11a7796bc..6241db97d1 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp @@ -1,247 +1,272 @@ /*=================================================================== 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 #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 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() { - // does nothing + 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 ); } diff --git a/Modules/DICOMReader/mitkDICOMSortByTag.cpp b/Modules/DICOMReader/mitkDICOMSortByTag.cpp index 6bb880ac24..0b2d3b8e24 100644 --- a/Modules/DICOMReader/mitkDICOMSortByTag.cpp +++ b/Modules/DICOMReader/mitkDICOMSortByTag.cpp @@ -1,133 +1,133 @@ /*=================================================================== 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 "mitkDICOMSortByTag.h" mitk::DICOMSortByTag ::DICOMSortByTag(const DICOMTag& tag, DICOMSortCriterion::Pointer secondaryCriterion) :DICOMSortCriterion(secondaryCriterion) ,m_Tag(tag) { } mitk::DICOMSortByTag ::~DICOMSortByTag() { } mitk::DICOMSortByTag ::DICOMSortByTag(const DICOMSortByTag& other ) :DICOMSortCriterion(other) ,m_Tag(other.m_Tag) { } mitk::DICOMSortByTag& mitk::DICOMSortByTag ::operator=(const DICOMSortByTag& other) { if (this != &other) { DICOMSortCriterion::operator=(other); m_Tag = other.m_Tag; } return *this; } mitk::DICOMTagList mitk::DICOMSortByTag ::GetTagsOfInterest() const { DICOMTagList list; list.push_back(m_Tag); return list; } bool mitk::DICOMSortByTag ::IsLeftBeforeRight(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right) const { return this->NumericCompare(left, right, m_Tag); } bool mitk::DICOMSortByTag ::StringCompare(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right, const DICOMTag& tag) const { assert(left); assert(right); std::string leftString = left->GetTagValueAsString(tag); std::string rightString = right->GetTagValueAsString(tag); if (leftString != rightString) { return leftString.compare(rightString) < 0; } else { if (m_SecondaryCriterion.IsNotNull()) { return m_SecondaryCriterion->IsLeftBeforeRight(left, right); } else { return leftString.compare(rightString) < 0; // TODO last resort needed } } } bool mitk::DICOMSortByTag ::NumericCompare(const mitk::DICOMDatasetAccess* left, const mitk::DICOMDatasetAccess* right, const DICOMTag& tag) const { assert(left); assert(right); std::string leftString = left->GetTagValueAsString(tag); std::string rightString = right->GetTagValueAsString(tag); const char* leftInput(leftString.c_str()); const char* rightInput(rightString.c_str()); char* leftEnd(NULL); char* rightEnd(NULL); double leftDouble = strtod(leftInput, &leftEnd); double rightDouble = strtod(rightInput, &rightEnd); - if (leftEnd == leftInput || rightEnd == rightInput) + if (leftEnd == leftInput || rightEnd == rightInput) // no numerical conversion.. { - return this->StringCompare(left,right, tag); + return this->StringCompare(left,right, tag); // fallback to string compare } else { - if (leftDouble != rightDouble) + if (leftDouble != rightDouble) // can we decide? { return leftDouble < rightDouble; } - else + else // ask secondary criterion { if (m_SecondaryCriterion.IsNotNull()) { return m_SecondaryCriterion->IsLeftBeforeRight(left, right); } else { return leftDouble < rightDouble; // TODO last resort } } } } diff --git a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp new file mode 100644 index 0000000000..4dbad05b3d --- /dev/null +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp @@ -0,0 +1,277 @@ +/*=================================================================== + +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 new file mode 100644 index 0000000000..296f20b743 --- /dev/null +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h @@ -0,0 +1,158 @@ +/*=================================================================== + +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 + +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 diff --git a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx new file mode 100644 index 0000000000..1952f4f311 --- /dev/null +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx @@ -0,0 +1,205 @@ +/*=================================================================== + +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 +#include +//#include +//#include +//#include + + +template +mitk::Image::Pointer +mitk::ITKDICOMSeriesReaderHelper +::LoadDICOMByITK( + const StringContainer& filenames, + bool correctTilt, + const GantryTiltInformation& tiltInfo, + DcmIoType::Pointer& io, + Image::Pointer preLoadedImageBlock ) +{ + /******** Normal Case, 3D (also for GDCM < 2 usable) ***************/ + mitk::Image::Pointer image = mitk::Image::New(); + + typedef itk::Image ImageType; + typedef itk::ImageSeriesReader ReaderType; + + io = DcmIoType::New(); + typename ReaderType::Pointer reader = ReaderType::New(); + + reader->SetImageIO(io); + reader->ReverseOrderOff(); + + if (preLoadedImageBlock.IsNull()) + { + reader->SetFileNames(filenames); + reader->Update(); + 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()); + } + else + { + image = preLoadedImageBlock; + StringContainer fakeList; + fakeList.push_back( filenames.front() ); + reader->SetFileNames( fakeList ); // we always need to load at least one file to get the MetaDataDictionary + reader->Update(); + } + + MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " + << image->GetDimension(1) << ", " + << image->GetDimension(2) << "]"; + + MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " + << image->GetGeometry()->GetSpacing()[1] << ", " + << image->GetGeometry()->GetSpacing()[2] << "]"; + + return image; +} + +template +typename ImageType::Pointer +mitk::ITKDICOMSeriesReaderHelper +::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::ScalableAffineTransform< double, ImageType::ImageDimension > TransformType; + typename TransformType::Pointer transformShear = TransformType::New(); + + /** + - 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 brings 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 transform the image volume back to its actual position (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 + */ + + ScalarType factor = tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() / input->GetSpacing()[1]; + // 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() ); + + typename TransformType::MatrixType indexToWorldMatrix; + indexToWorldMatrix = input->GetDirection(); + + typename ImageType::DirectionType scale; + for ( unsigned int i = 0; i < ImageType::ImageDimension; i++ ) + { + scale[i][i] = input->GetSpacing()[i]; + } + indexToWorldMatrix *= scale; + + imageIndexToWorld->SetMatrix( indexToWorldMatrix ); + + 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. + */ + + // TODO use (0028,0120) Pixel Padding Value if present + resampler->SetDefaultPixelValue( itk::NumericTraits< typename ImageType::PixelType >::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 + + + // if tilt positive, then we need additional pixels BELOW origin, otherwise we need pixels behind the end of the block + + // in any case we need more size to accomodate shifted slices + typename ImageType::SizeType largerSize = resampler->GetSize(); // now the resampler already holds the input image's size. + largerSize[1] += static_cast(tiltInfo.GetTiltCorrectedAdditionalSize() / input->GetSpacing()[1]+ 2.0); + resampler->SetSize( largerSize ); + + // in SOME cases this additional size is below/behind origin + if ( tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() > 0.0 ) + { + typename ImageType::DirectionType imageDirection = input->GetDirection(); + Vector3D yDirection; + yDirection[0] = imageDirection[0][1]; + yDirection[1] = imageDirection[1][1]; + yDirection[2] = imageDirection[2][1]; + yDirection.Normalize(); + + typename ImageType::PointType shiftedOrigin; + shiftedOrigin = input->GetOrigin(); + + // add some pixels to make everything fit + shiftedOrigin[0] -= yDirection[0] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]); + shiftedOrigin[1] -= yDirection[1] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]); + shiftedOrigin[2] -= yDirection[2] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]); + + resampler->SetOutputOrigin( shiftedOrigin ); + } + + 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; +} + +