diff --git a/Modules/DICOMReader/Testing/mitkDICOMFileReaderTestHelper.h b/Modules/DICOMReader/Testing/mitkDICOMFileReaderTestHelper.h index 7368d04f93..42bd14be87 100644 --- a/Modules/DICOMReader/Testing/mitkDICOMFileReaderTestHelper.h +++ b/Modules/DICOMReader/Testing/mitkDICOMFileReaderTestHelper.h @@ -1,153 +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) { 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 { 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); + MITK_DEBUG << "-------------------------------------------"; + MITK_DEBUG << "Output " << o << " at " << (void*) mitkImage.GetPointer(); + MITK_DEBUG << " Number of files: " << outputFiles.size(); + MITK_DEBUG << " Dimensions: " << mitkImage->GetDimension(0) << " " << mitkImage->GetDimension(1) << " " << mitkImage->GetDimension(2); } } }; // end test class } // namespace #endif diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp index 1b2ed69505..33c023558a 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.cpp @@ -1,335 +1,361 @@ /*=================================================================== 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 "mitkGantryTiltInformation.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; } +void +mitk::DICOMITKSeriesGDCMReader +::SetFixTiltByShearing(bool on) +{ + m_FixTiltByShearing = on; +} + 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->GetGroup() << "," << tagIter->GetElement(); gdcmScanner.AddTag( gdcm::Tag(tagIter->GetGroup(), tagIter->GetElement()) ); } } // Add some of our own interest gdcmScanner.AddTag( gdcm::Tag(0x0018,0x1164) ); gdcmScanner.AddTag( gdcm::Tag(0x0028,0x0030) ); 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; MITK_DEBUG << "================================================================================"; MITK_DEBUG << "DICOMIKTSeriesGDCMReader: " << ss.str() << ": " << m_SortingResultInProgress.size() << " groups input"; unsigned int groupIndex = 0; for(SortingBlockList::iterator blockIter = m_SortingResultInProgress.begin(); blockIter != m_SortingResultInProgress.end(); ++groupIndex, ++blockIter) { DICOMGDCMImageFrameList& gdcmInfoFrameList = *blockIter; DICOMDatasetList datasetList = ToDICOMDatasetList( gdcmInfoFrameList ); MITK_DEBUG << "--------------------------------------------------------------------------------"; MITK_DEBUG << "DICOMIKTSeriesGDCMReader: " << ss.str() << ", dataset group " << groupIndex << " (" << datasetList.size() << " datasets): "; for (DICOMDatasetList::iterator oi = datasetList.begin(); oi != datasetList.end(); ++oi) { MITK_DEBUG << " INPUT : " << (*oi)->GetFilenameIfAvailable(); } sorter->SetInput(datasetList); sorter->Sort(); unsigned int numberOfResultingBlocks = sorter->GetNumberOfOutputs(); for (unsigned int b = 0; b < numberOfResultingBlocks; ++b) { DICOMDatasetList blockResult = sorter->GetOutput(b); for (DICOMDatasetList::iterator oi = blockResult.begin(); oi != blockResult.end(); ++oi) { MITK_DEBUG << " OUTPUT(" << b << ") :" << (*oi)->GetFilenameIfAvailable(); } DICOMGDCMImageFrameList sortedGdcmInfoFrameList = FromDICOMDatasetList(blockResult); nextStepSorting.push_back( sortedGdcmInfoFrameList ); } } 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 ); + assert(!gdcmFrameInfoList.empty()); + assert(!frameList.empty()); DICOMImageBlockDescriptor block; block.SetImageFrameList( frameList ); - if (!gdcmFrameInfoList.empty()) - { - // assume - static const DICOMTag tagPixelSpacing(0x0028,0x0030); - static const DICOMTag tagImagerPixelSpacing(0x0018,0x1164); - std::string pixelSpacingString = (gdcmFrameInfoList.front())->GetTagValueAsString( tagPixelSpacing ); - std::string imagerPixelSpacingString = gdcmFrameInfoList.front()->GetTagValueAsString( tagImagerPixelSpacing ); - block.SetPixelSpacingInformation(pixelSpacingString, imagerPixelSpacingString); - } + bool hasTilt = false; + const GantryTiltInformation& tiltInfo = m_EquiDistantBlocksSorter->GetTiltInformation( (gdcmFrameInfoList.front())->GetFilenameIfAvailable(), hasTilt ); + block.SetHasGantryTilt( hasTilt ); + block.SetTiltInformation( tiltInfo ); + + // assume + static const DICOMTag tagPixelSpacing(0x0028,0x0030); + static const DICOMTag tagImagerPixelSpacing(0x0018,0x1164); + std::string pixelSpacingString = (gdcmFrameInfoList.front())->GetTagValueAsString( tagPixelSpacing ); + std::string imagerPixelSpacingString = gdcmFrameInfoList.front()->GetTagValueAsString( tagImagerPixelSpacing ); + block.SetPixelSpacingInformation(pixelSpacingString, imagerPixelSpacingString); + 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 + const GantryTiltInformation tiltInfo = block.GetTiltInformation(); + bool hasTilt = block.HasGantryTilt(); + if (hasTilt) + { + MITK_DEBUG << "When loading image " << o << ": got tilt info:"; + tiltInfo.Print(std::cout); + } + else + { + MITK_DEBUG << "When loading image " << o << ": has NO info."; + } ITKDICOMSeriesReaderHelper::StringContainer filenames; for (DICOMImageFrameList::const_iterator frameIter = frames.begin(); frameIter != frames.end(); ++frameIter) { filenames.push_back( (*frameIter)->Filename ); } - mitk::Image::Pointer mitkImage = helper.Load( filenames, correctTilt, tiltInfo ); // TODO preloaded images, caching..? + mitk::Image::Pointer mitkImage = helper.Load( filenames, m_FixTiltByShearing && hasTilt, tiltInfo ); // TODO preloaded images, caching..? Vector3D imageSpacing = mitkImage->GetGeometry()->GetSpacing(); ScalarType desiredSpacingX = imageSpacing[0]; ScalarType desiredSpacingY = imageSpacing[1]; block.GetDesiredMITKImagePixelSpacing( desiredSpacingX, desiredSpacingY ); // prefer pixel spacing over imager pixel spacing MITK_DEBUG << "Loaded image with spacing " << imageSpacing[0] << ", " << imageSpacing[1]; MITK_DEBUG << "Found correct spacing info " << desiredSpacingX << ", " << desiredSpacingY; imageSpacing[0] = desiredSpacingX; imageSpacing[1] = desiredSpacingY; mitkImage->GetGeometry()->SetSpacing( imageSpacing ); 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() { // TODO: cols, rows, etc. are also not optional - mitk::EquiDistantBlocksSorter::Pointer distanceKeeper = mitk::EquiDistantBlocksSorter::New(); - this->AddSortingElement( distanceKeeper ); + if (m_EquiDistantBlocksSorter.IsNull()) + { + m_EquiDistantBlocksSorter = mitk::EquiDistantBlocksSorter::New(); + } + m_EquiDistantBlocksSorter->SetAcceptTilt( m_FixTiltByShearing ); + + // TODO CHECK + this->AddSortingElement( m_EquiDistantBlocksSorter ); } diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h index 371a23480e..0d3626d715 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h @@ -1,79 +1,83 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDICOMITKSeriesGDCMReader_h #define mitkDICOMITKSeriesGDCMReader_h #include "mitkDICOMFileReader.h" #include "mitkDICOMDatasetSorter.h" #include "mitkDICOMGDCMImageFrameInfo.h" +#include "mitkEquiDistantBlocksSorter.h" #include "DICOMReaderExports.h" // TODO ensure "C" locale!! -// TODO tilt // TODO 3D+t 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); + void SetFixTiltByShearing(bool on); + 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; + + mitk::EquiDistantBlocksSorter::Pointer m_EquiDistantBlocksSorter; }; } #endif diff --git a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp index 39f3ecf02a..ef72c0347c 100644 --- a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp +++ b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.cpp @@ -1,211 +1,244 @@ /*=================================================================== 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 "mitkDICOMImageBlockDescriptor.h" mitk::DICOMImageBlockDescriptor ::DICOMImageBlockDescriptor() :m_PixelsInterpolated(false) ,m_PixelSpacingInterpretation() ,m_PixelSpacing("") ,m_ImagerPixelSpacing("") +,m_HasGantryTilt( false ) { } mitk::DICOMImageBlockDescriptor ::~DICOMImageBlockDescriptor() { } mitk::DICOMImageBlockDescriptor ::DICOMImageBlockDescriptor(const DICOMImageBlockDescriptor& other) :m_ImageFrameList( other.m_ImageFrameList ) ,m_MitkImage( other.m_MitkImage ) ,m_SliceIsLoaded( other.m_SliceIsLoaded ) ,m_PixelsInterpolated( other.m_PixelsInterpolated ) ,m_PixelSpacingInterpretation( other.m_PixelSpacingInterpretation ) ,m_PixelSpacing( other.m_PixelSpacing ) ,m_ImagerPixelSpacing( other.m_ImagerPixelSpacing ) +,m_HasGantryTilt( other.m_HasGantryTilt ) +,m_TiltInformation( other.m_TiltInformation ) { if (m_MitkImage) { m_MitkImage = m_MitkImage->Clone(); } } mitk::DICOMImageBlockDescriptor& mitk::DICOMImageBlockDescriptor ::operator=(const DICOMImageBlockDescriptor& other) { if (this != &other) { m_ImageFrameList = other.m_ImageFrameList; m_MitkImage = other.m_MitkImage; m_SliceIsLoaded = other.m_SliceIsLoaded; m_PixelsInterpolated = other.m_PixelsInterpolated; m_PixelSpacingInterpretation = other.m_PixelSpacingInterpretation; m_PixelSpacing = other.m_PixelSpacing; m_ImagerPixelSpacing = other.m_ImagerPixelSpacing; + m_HasGantryTilt = other.m_HasGantryTilt; + m_TiltInformation = other.m_TiltInformation; if (m_MitkImage) { m_MitkImage = m_MitkImage->Clone(); } } return *this; } +bool +mitk::DICOMImageBlockDescriptor +::HasGantryTilt() const +{ + return m_HasGantryTilt; +} + +bool +mitk::DICOMImageBlockDescriptor +::SetHasGantryTilt(bool hasi) +{ + m_HasGantryTilt = hasi; +} + +void +mitk::DICOMImageBlockDescriptor +::SetTiltInformation(const GantryTiltInformation& info) +{ + m_TiltInformation = info; +} + +const mitk::GantryTiltInformation +mitk::DICOMImageBlockDescriptor +::GetTiltInformation() const +{ + return m_TiltInformation; +} + void mitk::DICOMImageBlockDescriptor ::SetImageFrameList(const DICOMImageFrameList& framelist) { m_ImageFrameList = framelist; m_SliceIsLoaded.resize(framelist.size()); m_SliceIsLoaded.assign(framelist.size(), false); } const mitk::DICOMImageFrameList& mitk::DICOMImageBlockDescriptor ::GetImageFrameList() const { return m_ImageFrameList; } void mitk::DICOMImageBlockDescriptor ::SetMitkImage(Image::Pointer image) { if (m_MitkImage != image) { m_MitkImage = image; } } mitk::Image::Pointer mitk::DICOMImageBlockDescriptor ::GetMitkImage() const { return m_MitkImage; } void mitk::DICOMImageBlockDescriptor ::SetSliceIsLoaded(unsigned int index, bool isLoaded) { if (index < m_SliceIsLoaded.size()) { m_SliceIsLoaded[index] = isLoaded; } else { std::stringstream ss; ss << "Index " << index << " out of range (" << m_SliceIsLoaded.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } bool mitk::DICOMImageBlockDescriptor ::IsSliceLoaded(unsigned int index) const { if (index < m_SliceIsLoaded.size()) { return m_SliceIsLoaded[index]; } else { std::stringstream ss; ss << "Index " << index << " out of range (" << m_SliceIsLoaded.size() << " indices reserved)"; throw std::invalid_argument( ss.str() ); } } bool mitk::DICOMImageBlockDescriptor ::AllSlicesAreLoaded() const { bool allLoaded = true; for (BoolList::const_iterator iter = m_SliceIsLoaded.begin(); iter != m_SliceIsLoaded.end(); ++iter) { allLoaded &= *iter; } return allLoaded; } void mitk::DICOMImageBlockDescriptor ::SetPixelsInterpolated(bool pixelsAreInterpolated) { if (m_PixelsInterpolated != pixelsAreInterpolated) { m_PixelsInterpolated = pixelsAreInterpolated; } } bool mitk::DICOMImageBlockDescriptor ::GetPixelsInterpolated() const { return m_PixelsInterpolated; } void mitk::DICOMImageBlockDescriptor ::SetPixelSpacingInterpretation( PixelSpacingInterpretation interpretation ) { if (m_PixelSpacingInterpretation != interpretation) { m_PixelSpacingInterpretation = interpretation; } } mitk::PixelSpacingInterpretation mitk::DICOMImageBlockDescriptor ::GetPixelSpacingInterpretation() const { return m_PixelSpacingInterpretation; } void mitk::DICOMImageBlockDescriptor ::SetPixelSpacingInformation(const std::string& pixelSpacing, const std::string& imagerPixelSpacing) { m_PixelSpacing = pixelSpacing; m_ImagerPixelSpacing = imagerPixelSpacing; } void mitk::DICOMImageBlockDescriptor ::GetDesiredMITKImagePixelSpacing( float& spacingX, float& spacingY) const { // preference for "in patient" pixel spacing if ( !DICOMStringToSpacing( m_PixelSpacing, spacingX, spacingY ) ) { // fallback to "on detector" spacing if ( !DICOMStringToSpacing( m_ImagerPixelSpacing, spacingX, spacingY ) ) { // last resort: invent something spacingX = spacingY = 1.0; } } } diff --git a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.h b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.h index 5952aa8ba5..7c8a827fca 100644 --- a/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.h +++ b/Modules/DICOMReader/mitkDICOMImageBlockDescriptor.h @@ -1,72 +1,83 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDICOMImageBlockDescriptor_h #define mitkDICOMImageBlockDescriptor_h #include "mitkDICOMEnums.h" #include "mitkDICOMImageFrameInfo.h" #include "mitkDICOMTag.h" #include "mitkImage.h" +#include "mitkGantryTiltInformation.h" + namespace mitk { class DICOMReader_EXPORT DICOMImageBlockDescriptor { public: DICOMImageBlockDescriptor(); ~DICOMImageBlockDescriptor(); DICOMImageBlockDescriptor(const DICOMImageBlockDescriptor& other); DICOMImageBlockDescriptor& operator=(const DICOMImageBlockDescriptor& other); void SetImageFrameList(const DICOMImageFrameList& framelist); const DICOMImageFrameList& GetImageFrameList() const; void SetMitkImage(Image::Pointer image); Image::Pointer GetMitkImage() const; void SetSliceIsLoaded(unsigned int index, bool isLoaded); bool IsSliceLoaded(unsigned int index) const; bool AllSlicesAreLoaded() const; void SetPixelsInterpolated(bool pixelsAreInterpolated); bool GetPixelsInterpolated() const; void SetPixelSpacingInterpretation( PixelSpacingInterpretation interpretation ); PixelSpacingInterpretation GetPixelSpacingInterpretation() const; void SetPixelSpacingInformation(const std::string& pixelSpacing, const std::string& imagerPixelSpacing); void GetDesiredMITKImagePixelSpacing( float& spacingX, float& spacingY) const; + bool HasGantryTilt() const; + bool SetHasGantryTilt(bool hasi); + + void SetTiltInformation(const GantryTiltInformation& info); + const GantryTiltInformation GetTiltInformation() const; + private: DICOMImageFrameList m_ImageFrameList; Image::Pointer m_MitkImage; BoolList m_SliceIsLoaded; bool m_PixelsInterpolated; PixelSpacingInterpretation m_PixelSpacingInterpretation; std::string m_PixelSpacing; std::string m_ImagerPixelSpacing; + + bool m_HasGantryTilt; + GantryTiltInformation m_TiltInformation; }; } #endif diff --git a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp index e9be5c1c85..8783460ebc 100644 --- a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.cpp @@ -1,390 +1,457 @@ /*=================================================================== 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() +::SetFirstFilenameOfBlock(const std::string& filename) +{ + m_FirstFilenameOfBlock = filename; +} + +std::string +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::GetFirstFilenameOfBlock() const +{ + return m_FirstFilenameOfBlock; +} + +void +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::FlagGantryTilt(const GantryTiltInformation& tiltInfo) { m_GantryTilt = true; + m_TiltInfo = tiltInfo; +} + +const mitk::GantryTiltInformation& +mitk::EquiDistantBlocksSorter::SliceGroupingAnalysisResult +::GetTiltInfo() const +{ + return m_TiltInfo; } 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() +,m_AcceptTilt(false) +{ +} + +mitk::EquiDistantBlocksSorter +::EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other ) +:DICOMDatasetSorter(other) +,m_AcceptTilt(false) { } mitk::EquiDistantBlocksSorter ::~EquiDistantBlocksSorter() { } +void mitk::EquiDistantBlocksSorter -::EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other ) -:DICOMDatasetSorter(other) +::SetAcceptTilt(bool accept) { + m_AcceptTilt = accept; } + 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( DICOMTag(0x0020, 0x0032) ); // ImagePositionPatient tags.push_back( DICOMTag(0x0020, 0x0037) ); // ImageOrientationPatient + tags.push_back( DICOMTag(0x0018, 0x1120) ); // GantryTilt return tags; } void mitk::EquiDistantBlocksSorter ::Sort() { DICOMDatasetList remainingInput = GetInput(); // copy - bool acceptGantryTilt = false; // TODO - typedef std::list OutputListType; OutputListType outputs; + m_SliceGroupingResults.clear(); + while (!remainingInput.empty()) // repeat until all files are grouped somehow { - SliceGroupingAnalysisResult regularBlock = this->AnalyzeFileForITKImageSeriesReaderSpacingAssumption( remainingInput, acceptGantryTilt ); + SliceGroupingAnalysisResult regularBlock = this->AnalyzeFileForITKImageSeriesReaderSpacingAssumption( remainingInput, m_AcceptTilt ); DICOMDatasetList inBlock = regularBlock.GetBlockFilenames(); DICOMDatasetList laterBlock = regularBlock.GetUnsortedFilenames(); MITK_DEBUG << "Result: sorted 3D group with " << inBlock.size() << " files"; for (DICOMDatasetList::const_iterator diter = inBlock.begin(); diter != inBlock.end(); ++diter) MITK_DEBUG << " IN " << (*diter)->GetFilenameIfAvailable(); for (DICOMDatasetList::const_iterator diter = laterBlock.begin(); diter != laterBlock.end(); ++diter) MITK_DEBUG << " OUT " << (*diter)->GetFilenameIfAvailable(); outputs.push_back( regularBlock.GetBlockFilenames() ); + m_SliceGroupingResults.push_back( regularBlock ); 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); } } +const mitk::GantryTiltInformation +mitk::EquiDistantBlocksSorter +::GetTiltInformation(const std::string& filename, bool& hasTiltInfo) +{ + for (ResultsList::iterator ri = m_SliceGroupingResults.begin(); + ri != m_SliceGroupingResults.end(); + ++ri) + { + SliceGroupingAnalysisResult& result = *ri; + + if (filename == result.GetFirstFilenameOfBlock()) + { + hasTiltInfo = result.ContainsGantryTilt(); // this is a returning statement, don't remove + if (hasTiltInfo) + { + return result.GetTiltInfo(); + } + } + } + + return GantryTiltInformation(); // empty +} + std::string mitk::EquiDistantBlocksSorter ::ConstCharStarToString(const char* s) { return s ? std::string(s) : std::string(); } 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 = DICOMTag(0x0020,0x0032); // Image Position (Patient) const DICOMTag tagImageOrientation = DICOMTag(0x0020, 0x0037); // Image Orientation const DICOMTag tagGantryTilt = DICOMTag(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 " << datasets.size() << " files for z-spacing assumption of ITK's ImageSeriesReader (group tilted: " << groupImagesWithGantryTilt << ")"; unsigned int fileIndex(0); for (DICOMDatasetList::const_iterator dsIter = datasets.begin(); dsIter != datasets.end(); ++dsIter, ++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 = (*dsIter)->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 " << *dsIter << " 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( *dsIter ); DICOMDatasetList remainingFiles; remainingFiles.insert( remainingFiles.end(), dsIter+1, datasets.end() ); result.AddFilesToUnsortedBlock( remainingFiles ); fileFitsIntoPattern = false; break; // no files anymore } else { // ==> this does not match, consider later result.AddFileToUnsortedBlock( *dsIter ); // 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 << " " << (*dsIter)->GetFilenameIfAvailable() << " at " /* << thisOriginString */ << "(" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")"; if ( lastOriginInitialized && (thisOrigin == lastOrigin) ) { MITK_DEBUG << " ==> Sort away " << *dsIter << " for separate time step"; // we already have one occupying this position result.AddFileToUnsortedBlock( *dsIter ); // 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 = (*dsIter)->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 (*dsIter) // 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[dsIter->c_str()].find(tagGantryTilt) != tagValueMappings[dsIter->c_str()].end() */) + double angle = 0.0; + std::string tiltStr = (*dsIter)->GetTagValueAsString( tagGantryTilt ); + const char* convertInput = tiltStr.c_str(); + char* convertEnd(NULL); + if (!tiltStr.empty()) { // read value, compare to calculated angle - std::string tiltStr = (*dsIter)->GetTagValueAsString( tagGantryTilt ); - double angle = atof(tiltStr.c_str()); + angle = strtod(convertInput, &convertEnd); // TODO check errors! + } + if (convertEnd != convertInput) + { 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( *dsIter ); // 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(); + assert(!datasets.empty()); + + result.FlagGantryTilt(tiltInfo); result.AddFileToSortedBlock( *dsIter ); // this file is good for current block + result.SetFirstFilenameOfBlock( datasets.front()->GetFilenameIfAvailable() ); 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(); + assert(!datasets.empty()); + + result.FlagGantryTilt(tiltInfo); result.AddFileToSortedBlock( *dsIter ); // this file is good for current block + result.SetFirstFilenameOfBlock( datasets.front()->GetFilenameIfAvailable() ); fileFitsIntoPattern = true; } } else // caller does not want tilt compensation OR shearing is more complicated than tilt { result.AddFileToUnsortedBlock( *dsIter ); // sort away for further analysis fileFitsIntoPattern = false; } } else // not sheared { result.AddFileToSortedBlock( *dsIter ); // 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 " << *dsIter << " 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( *dsIter ); // sort away for further analysis fileFitsIntoPattern = false; } else { result.AddFileToSortedBlock( *dsIter ); // this file is good for current block fileFitsIntoPattern = true; } } else // this should be the very first slice { result.AddFileToSortedBlock( *dsIter ); // 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 index cfdf1685c1..293943a73c 100644 --- a/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h +++ b/Modules/DICOMReader/mitkEquiDistantBlocksSorter.h @@ -1,137 +1,156 @@ /*=================================================================== 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 "mitkGantryTiltInformation.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(); + void SetAcceptTilt(bool accept); + const GantryTiltInformation GetTiltInformation(const std::string& filename, bool& hasTiltInfo); // TODO ugly, but bool flag will be removed + 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(); + void SetFirstFilenameOfBlock(const std::string& filename); + std::string GetFirstFilenameOfBlock() const; + /** \brief Remaining files, which could not be grouped. */ DICOMDatasetList GetUnsortedFilenames(); /** \brief Wheter or not the grouped result contain a gantry tilt. */ bool ContainsGantryTilt(); + /** + \brief Detailed description of gantry tilt. + */ + const GantryTiltInformation& GetTiltInfo() const; + /** \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(); + void FlagGantryTilt(const GantryTiltInformation& tiltInfo); /** \brief Only meaningful for use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption. */ void UndoPrematureGrouping(); protected: DICOMDatasetList m_GroupedFiles; DICOMDatasetList m_UnsortedFiles; - bool m_GantryTilt; + bool m_GantryTilt; // TODO make the flag part of GantryTiltInformation! + GantryTiltInformation m_TiltInfo; + std::string m_FirstFilenameOfBlock; }; /** \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); EquiDistantBlocksSorter(); virtual ~EquiDistantBlocksSorter(); EquiDistantBlocksSorter(const EquiDistantBlocksSorter& other); EquiDistantBlocksSorter& operator=(const EquiDistantBlocksSorter& other); + + bool m_AcceptTilt; + + typedef std::vector ResultsList; + ResultsList m_SliceGroupingResults; }; } #endif diff --git a/Modules/DICOMReader/mitkGantryTiltInformation.cpp b/Modules/DICOMReader/mitkGantryTiltInformation.cpp index e0109bbb7b..f4874cc972 100644 --- a/Modules/DICOMReader/mitkGantryTiltInformation.cpp +++ b/Modules/DICOMReader/mitkGantryTiltInformation.cpp @@ -1,185 +1,200 @@ #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 origin1: " << origin1; + MITK_DEBUG << " v origin2: " << origin2; 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; } } +void +mitk::GantryTiltInformation +::Print(std::ostream& os) const +{ + os << " calculated from slices " << m_NumberOfSlicesApart << " slices apart " << std::endl; + os << " shift normal: " << m_ShiftNormal << std::endl; + os << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing << std::endl; + os << " shift up: " << m_ShiftUp << std::endl; + os << " shift right: " << m_ShiftRight << std::endl; + + os << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535 << std::endl; +} + 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 +mitk::GantryTiltInformation::GetTiltCorrectedAdditionalSize(unsigned int imageSizeZ) const { - return fabs(m_ShiftUp); + return fabs(m_ShiftUp / static_cast(m_NumberOfSlicesApart) * static_cast(imageSizeZ-1)); // TODO care for SIZE of image } 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 index ed9736b45e..b43b5a5f48 100644 --- a/Modules/DICOMReader/mitkGantryTiltInformation.h +++ b/Modules/DICOMReader/mitkGantryTiltInformation.h @@ -1,127 +1,129 @@ /*=================================================================== 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(); + void Print(std::ostream& os) const; + /** \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; + double GetTiltCorrectedAdditionalSize(unsigned int imageSizeZ) const; // TODO: imageSizeZ is a bit of a loss of information: when loading, GantryTiltInformation *should* be initialized with the most distant slices /** \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 db8cdc6926..8579d2a6c6 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp @@ -1,97 +1,96 @@ /*=================================================================== 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::Image::Pointer mitk::ITKDICOMSeriesReaderHelper -::Load( const StringContainer& filenames, bool correctTilt, bool itkNotUsed(tiltInfo) ) +::Load( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& 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 a2d5040980..b5c0d53ab8 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h @@ -1,56 +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 "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 - Image::Pointer Load( const StringContainer& filenames, bool correctTilt, bool tiltInfo ); + Image::Pointer Load( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& 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 index e1e6b2c10f..172c85209b 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx @@ -1,205 +1,209 @@ /*=================================================================== 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->ReverseOrderOff(); 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 ) { + tiltInfo.Print(std::cout); 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); + double imageSizeZ = largerSize[2]; + MITK_DEBUG <<"Calculate lager size = " << largerSize[1] << " + " << tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) << " / " << input->GetSpacing()[1] << "+ 2.0"; + largerSize[1] += static_cast(tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) / input->GetSpacing()[1]+ 2.0); resampler->SetSize( largerSize ); + MITK_DEBUG << "Fix Y size of image w/ spacing " << input->GetSpacing()[1] << " from " << input->GetLargestPossibleRegion().GetSize()[1] << " to " << largerSize[1]; // 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]); + shiftedOrigin[0] -= yDirection[0] * (tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) + 1.0 * input->GetSpacing()[1]); + shiftedOrigin[1] -= yDirection[1] * (tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) + 1.0 * input->GetSpacing()[1]); + shiftedOrigin[2] -= yDirection[2] * (tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) + 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; } diff --git a/Modules/DICOMTesting/mitkTestDICOMLoading.cpp b/Modules/DICOMTesting/mitkTestDICOMLoading.cpp index 179481ed85..00491c74ec 100644 --- a/Modules/DICOMTesting/mitkTestDICOMLoading.cpp +++ b/Modules/DICOMTesting/mitkTestDICOMLoading.cpp @@ -1,436 +1,437 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "mitkTestDICOMLoading.h" #include mitk::TestDICOMLoading::TestDICOMLoading() :m_PreviousCLocale(NULL) { } void mitk::TestDICOMLoading::SetDefaultLocale() { // remember old locale only once if (m_PreviousCLocale == NULL) { m_PreviousCLocale = setlocale(LC_NUMERIC, NULL); // set to "C" setlocale(LC_NUMERIC, "C"); m_PreviousCppLocale = std::cin.getloc(); std::locale l( "C" ); std::cin.imbue(l); std::cout.imbue(l); } } void mitk::TestDICOMLoading::ResetUserLocale() { if (m_PreviousCLocale) { setlocale(LC_NUMERIC, m_PreviousCLocale); std::cin.imbue(m_PreviousCppLocale); std::cout.imbue(m_PreviousCppLocale); m_PreviousCLocale = NULL; } } mitk::TestDICOMLoading::ImageList mitk::TestDICOMLoading ::LoadFiles( const StringList& files ) { for (StringList::const_iterator iter = files.begin(); iter != files.end(); ++iter) { MITK_DEBUG << "File " << *iter; } ImageList result; ClassicDICOMSeriesReader::Pointer reader = ClassicDICOMSeriesReader::New(); + reader->SetFixTiltByShearing(true); reader->SetInputFiles( files ); reader->AnalyzeInputFiles(); - reader->PrintOutputs(std::cout); + reader->PrintOutputs(std::cout,true); reader->LoadImages(); unsigned int numberOfImages = reader->GetNumberOfOutputs(); for (unsigned imageIndex = 0; imageIndex < numberOfImages; ++imageIndex) { const DICOMImageBlockDescriptor& block = reader->GetOutput(imageIndex); result.push_back( block.GetMitkImage() ); } return result; } std::string mitk::TestDICOMLoading::ComponentTypeToString(int type) { if (type == itk::ImageIOBase::UCHAR) return "UCHAR"; else if (type == itk::ImageIOBase::CHAR) return "CHAR"; else if (type == itk::ImageIOBase::USHORT) return "USHORT"; else if (type == itk::ImageIOBase::SHORT) return "SHORT"; else if (type == itk::ImageIOBase::UINT) return "UINT"; else if (type == itk::ImageIOBase::INT) return "INT"; else if (type == itk::ImageIOBase::ULONG) return "ULONG"; else if (type == itk::ImageIOBase::LONG) return "LONG"; else if (type == itk::ImageIOBase::FLOAT) return "FLOAT"; else if (type == itk::ImageIOBase::DOUBLE) return "DOUBLE"; else return "UNKNOWN"; } // add a line to stringstream result (see DumpImageInformation #define DumpLine(field, data) DumpILine(0, field, data) // add an indented(!) line to stringstream result (see DumpImageInformation #define DumpILine(indent, field, data) \ { \ std::string DumpLine_INDENT; DumpLine_INDENT.resize(indent, ' ' ); \ result << DumpLine_INDENT << field << ": " << data << "\n"; \ } std::string mitk::TestDICOMLoading::DumpImageInformation( const Image* image ) { std::stringstream result; if (image == NULL) return result.str(); SetDefaultLocale(); // basic image data DumpLine( "Pixeltype", ComponentTypeToString(image->GetPixelType().GetComponentType()) ); DumpLine( "BitsPerPixel", image->GetPixelType().GetBpe() ); DumpLine( "Dimension", image->GetDimension() ); result << "Dimensions: "; for (unsigned int dim = 0; dim < image->GetDimension(); ++dim) result << image->GetDimension(dim) << " "; result << "\n"; // geometry data result << "Geometry: \n"; Geometry3D* geometry = image->GetGeometry(); if (geometry) { AffineTransform3D* transform = geometry->GetIndexToWorldTransform(); if (transform) { result << " " << "Matrix: "; const AffineTransform3D::MatrixType& matrix = transform->GetMatrix(); for (unsigned int i = 0; i < 3; ++i) for (unsigned int j = 0; j < 3; ++j) result << matrix[i][j] << " "; result << "\n"; result << " " << "Offset: "; const AffineTransform3D::OutputVectorType& offset = transform->GetOffset(); for (unsigned int i = 0; i < 3; ++i) result << offset[i] << " "; result << "\n"; result << " " << "Center: "; const AffineTransform3D::InputPointType& center = transform->GetCenter(); for (unsigned int i = 0; i < 3; ++i) result << center[i] << " "; result << "\n"; result << " " << "Translation: "; const AffineTransform3D::OutputVectorType& translation = transform->GetTranslation(); for (unsigned int i = 0; i < 3; ++i) result << translation[i] << " "; result << "\n"; result << " " << "Scale: "; const double* scale = transform->GetScale(); for (unsigned int i = 0; i < 3; ++i) result << scale[i] << " "; result << "\n"; result << " " << "Origin: "; const Point3D& origin = geometry->GetOrigin(); for (unsigned int i = 0; i < 3; ++i) result << origin[i] << " "; result << "\n"; result << " " << "Spacing: "; const Vector3D& spacing = geometry->GetSpacing(); for (unsigned int i = 0; i < 3; ++i) result << spacing[i] << " "; result << "\n"; result << " " << "TimeBounds: "; const TimeBounds timeBounds = geometry->GetTimeBounds(); for (unsigned int i = 0; i < 2; ++i) result << timeBounds[i] << " "; result << "\n"; } } ResetUserLocale(); return result.str(); } std::string mitk::TestDICOMLoading::trim(const std::string& pString, const std::string& pWhitespace) { const size_t beginStr = pString.find_first_not_of(pWhitespace); if (beginStr == std::string::npos) { // no content return ""; } const size_t endStr = pString.find_last_not_of(pWhitespace); const size_t range = endStr - beginStr + 1; return pString.substr(beginStr, range); } std::string mitk::TestDICOMLoading::reduce(const std::string& pString, const std::string& pFill, const std::string& pWhitespace) { // trim first std::string result(trim(pString, pWhitespace)); // replace sub ranges size_t beginSpace = result.find_first_of(pWhitespace); while (beginSpace != std::string::npos) { const size_t endSpace = result.find_first_not_of(pWhitespace, beginSpace); const size_t range = endSpace - beginSpace; result.replace(beginSpace, range, pFill); const size_t newStart = beginSpace + pFill.length(); beginSpace = result.find_first_of(pWhitespace, newStart); } return result; } bool mitk::TestDICOMLoading::CompareSpacedValueFields( const std::string& reference, const std::string& test, double /*eps*/ ) { bool result(true); // tokenize string, compare each token, if possible by float comparison std::stringstream referenceStream(reduce(reference)); std::stringstream testStream(reduce(test)); std::string refToken; std::string testToken; while ( std::getline( referenceStream, refToken, ' ' ) && std::getline ( testStream, testToken, ' ' ) ) { float refNumber; float testNumber; if ( this->StringToNumber(refToken, refNumber) ) { if ( this->StringToNumber(testToken, testNumber) ) { // print-out compared tokens if DEBUG output allowed MITK_DEBUG << "Reference Token '" << refToken << "'" << " value " << refNumber << ", test Token '" << testToken << "'" << " value " << testNumber; bool old_result = result; result &= ( fabs(refNumber - testNumber) < mitk::eps ); // log the token/number which causes the test to fail if( old_result != result) { MITK_ERROR << std::setprecision(16) << "Reference Token '" << refToken << "'" << " value " << refNumber << ", test Token '" << testToken << "'" << " value " << testNumber; MITK_ERROR << "[FALSE] - difference: " << std::setprecision(16) << fabs(refNumber - testNumber) << " EPS: " << mitk::eps; } } else { MITK_ERROR << refNumber << " cannot be compared to '" << testToken << "'"; } } else { MITK_DEBUG << "Token '" << refToken << "'" << " handled as string"; result &= refToken == testToken; } } if ( std::getline( referenceStream, refToken, ' ' ) ) { MITK_ERROR << "Reference string still had values when test string was already parsed: ref '" << reference << "', test '" << test << "'"; result = false; } else if ( std::getline( testStream, testToken, ' ' ) ) { MITK_ERROR << "Test string still had values when reference string was already parsed: ref '" << reference << "', test '" << test << "'"; result = false; } return result; } bool mitk::TestDICOMLoading::CompareImageInformationDumps( const std::string& referenceDump, const std::string& testDump ) { KeyValueMap reference = ParseDump(referenceDump); KeyValueMap test = ParseDump(testDump); bool testResult(true); // verify all expected values for (KeyValueMap::const_iterator refIter = reference.begin(); refIter != reference.end(); ++refIter) { const std::string& refKey = refIter->first; const std::string& refValue = refIter->second; if ( test.find(refKey) != test.end() ) { const std::string& testValue = test[refKey]; bool thisTestResult = CompareSpacedValueFields( refValue, testValue ); testResult &= thisTestResult; MITK_DEBUG << refKey << ": '" << refValue << "' == '" << testValue << "' ? " << (thisTestResult?"YES":"NO"); } else { MITK_ERROR << "Reference dump contains a key'" << refKey << "' (value '" << refValue << "')." ; MITK_ERROR << "This key is expected to be generated for tests (but was not). Most probably you need to update your test data."; return false; } } // now check test dump does not contain any additional keys for (KeyValueMap::const_iterator testIter = test.begin(); testIter != test.end(); ++testIter) { const std::string& key = testIter->first; const std::string& value = testIter->second; if ( reference.find(key) == reference.end() ) { MITK_ERROR << "Test dump contains an unexpected key'" << key << "' (value '" << value << "')." ; MITK_ERROR << "This key is not expected. Most probably you need to update your test data."; return false; } } return testResult; } mitk::TestDICOMLoading::KeyValueMap mitk::TestDICOMLoading::ParseDump( const std::string& dump ) { KeyValueMap parsedResult; std::string shredder(dump); std::stack surroundingKeys; std::stack expectedIndents; expectedIndents.push(0); while (true) { std::string::size_type newLinePos = shredder.find( '\n' ); if (newLinePos == std::string::npos || newLinePos == 0) break; std::string line = shredder.substr( 0, newLinePos ); shredder = shredder.erase( 0, newLinePos+1 ); std::string::size_type keyPosition = line.find_first_not_of( ' ' ); std::string::size_type colonPosition = line.find( ':' ); std::string key = line.substr(keyPosition, colonPosition - keyPosition); std::string::size_type firstSpacePosition = key.find_first_of(" "); if (firstSpacePosition != std::string::npos) { key.erase(firstSpacePosition); } if ( keyPosition > expectedIndents.top() ) { // more indent than before expectedIndents.push(keyPosition); } else if (keyPosition == expectedIndents.top() ) { if (!surroundingKeys.empty()) { surroundingKeys.pop(); // last of same length } } else { // less indent than before do expectedIndents.pop(); while (expectedIndents.top() != keyPosition); // unwind until current indent is found } if (!surroundingKeys.empty()) { key = surroundingKeys.top() + "." + key; // construct current key name } surroundingKeys.push(key); // this is the new embracing key std::string value = line.substr(colonPosition+1); MITK_DEBUG << " Key: '" << key << "' value '" << value << "'" ; parsedResult[key] = value; // store parsing result } return parsedResult; }