diff --git a/Modules/DICOMReader/Resources/configurations/3D/descriptionbasedgrouping.xml b/Modules/DICOMReader/Resources/configurations/3D/imageposition.xml similarity index 100% rename from Modules/DICOMReader/Resources/configurations/3D/descriptionbasedgrouping.xml rename to Modules/DICOMReader/Resources/configurations/3D/imageposition.xml diff --git a/Modules/DICOMReader/files.cmake b/Modules/DICOMReader/files.cmake index c106e0c49f..e61be227c3 100644 --- a/Modules/DICOMReader/files.cmake +++ b/Modules/DICOMReader/files.cmake @@ -1,51 +1,51 @@ set(H_FILES mitkDICOMFileReader.h mitkDICOMImageFrameInfo.h mitkDICOMImageBlockDescriptor.h mitkDICOMGDCMImageFrameInfo.h mitkDICOMITKSeriesGDCMReader.h mitkDICOMDatasetSorter.h mitkDICOMEnums.h mitkDICOMTagBasedSorter.h mitkDICOMSortCriterion.h mitkDICOMSortByTag.h mitkEquiDistantBlocksSorter.h mitkSortByImagePositionPatient.h mitkClassicDICOMSeriesReader.h mitkThreeDnTDICOMSeriesReader.h mitkDICOMTag.h mitkDICOMReaderConfigurator.h mitkDICOMFileReaderSelector.h ) set(CPP_FILES mitkDICOMFileReader.cpp mitkDICOMImageBlockDescriptor.cpp mitkDICOMITKSeriesGDCMReader.cpp mitkDICOMDatasetSorter.cpp mitkDICOMTagBasedSorter.cpp mitkDICOMGDCMImageFrameInfo.cpp mitkDICOMImageFrameInfo.cpp mitkDICOMSortCriterion.cpp mitkDICOMSortByTag.cpp mitkITKDICOMSeriesReaderHelper.cpp mitkEquiDistantBlocksSorter.cpp mitkSortByImagePositionPatient.cpp mitkGantryTiltInformation.cpp mitkClassicDICOMSeriesReader.cpp mitkThreeDnTDICOMSeriesReader.cpp mitkDICOMTag.cpp mitkDICOMEnums.cpp mitkDICOMReaderConfigurator.cpp mitkDICOMFileReaderSelector.cpp ) set(RESOURCE_FILES configurations/3D/classicreader.xml - configurations/3D/descriptionbasedgrouping.xml + configurations/3D/imageposition.xml configurations/3D/instancenumber.xml configurations/3D/slicelocation.xml configurations/3D/imagetime.xml configurations/3DnT/classicreader.xml ) diff --git a/Modules/DICOMReader/mitkDICOMFileReaderSelector.cpp b/Modules/DICOMReader/mitkDICOMFileReaderSelector.cpp index 5d9509d169..7ee1bf1581 100644 --- a/Modules/DICOMReader/mitkDICOMFileReaderSelector.cpp +++ b/Modules/DICOMReader/mitkDICOMFileReaderSelector.cpp @@ -1,184 +1,226 @@ /*=================================================================== 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 "mitkDICOMFileReaderSelector.h" #include "mitkDICOMReaderConfigurator.h" #include "mitkModuleContext.h" #include #include #include #include mitk::DICOMFileReaderSelector ::DICOMFileReaderSelector() { } mitk::DICOMFileReaderSelector ::~DICOMFileReaderSelector() { } void mitk::DICOMFileReaderSelector ::AddConfigsFromResources(const std::string& path) { std::vector configs = GetModuleContext()->GetModule()->FindResources(path, "*.xml", false); for (std::vector::iterator iter = configs.begin(); iter != configs.end(); ++iter) { ModuleResource& resource = *iter; if (resource.IsValid()) { ModuleResourceStream stream(resource); // read all into string s std::string s; stream.seekg(0, std::ios::end); s.reserve(stream.tellg()); stream.seekg(0, std::ios::beg); s.assign((std::istreambuf_iterator(stream)), std::istreambuf_iterator()); //MITK_INFO << "------------ One resource with content\n" << s << "\n------- end -------"; this->AddConfig(s); } } } +void +mitk::DICOMFileReaderSelector +::AddConfigFromResource(ModuleResource& resource) +{ + if (resource.IsValid()) + { + ModuleResourceStream stream(resource); + + // read all into string s + std::string s; + + stream.seekg(0, std::ios::end); + s.reserve(stream.tellg()); + stream.seekg(0, std::ios::beg); + + s.assign((std::istreambuf_iterator(stream)), + std::istreambuf_iterator()); + + //MITK_INFO << "------------ One resource with content\n" << s << "\n------- end -------"; + this->AddConfig(s); + } +} + +void +mitk::DICOMFileReaderSelector +::AddConfigFromResource(const std::string& resourcename) +{ + ModuleResource r = GetModuleContext()->GetModule()->GetResource(resourcename); + this->AddConfigFromResource(r); +} + + void mitk::DICOMFileReaderSelector ::LoadBuiltIn3DConfigs() { - this->AddConfigsFromResources("configurations/3D"); + //this->AddConfigsFromResources("configurations/3D"); + // in this order of preference... + this->AddConfigFromResource("configurations/3D/instancenumber.xml"); + this->AddConfigFromResource("configurations/3D/imageposition.xml"); + this->AddConfigFromResource("configurations/3D/slicelocation.xml"); + this->AddConfigFromResource("configurations/3D/imagetime.xml"); } void mitk::DICOMFileReaderSelector ::LoadBuiltIn3DnTConfigs() { this->AddConfigsFromResources("configurations/3DnT"); } void mitk::DICOMFileReaderSelector ::AddConfig(const std::string& xmlDescription) { DICOMReaderConfigurator::Pointer configurator = DICOMReaderConfigurator::New(); DICOMFileReader::Pointer reader = configurator->CreateFromUTF8ConfigString(xmlDescription); if (reader.IsNotNull()) { m_Readers.push_back( reader ); m_PossibleConfigurations.push_back(xmlDescription); } else { std::stringstream ss; ss << "Could not parse reader configuration. Ignoring it."; throw std::invalid_argument( ss.str() ); } } void mitk::DICOMFileReaderSelector ::AddConfigFile(const std::string& filename) { std::ifstream file(filename.c_str()); std::string s; file.seekg(0, std::ios::end); s.reserve(file.tellg()); file.seekg(0, std::ios::beg); s.assign((std::istreambuf_iterator(file)), std::istreambuf_iterator()); this->AddConfig(s); } void mitk::DICOMFileReaderSelector ::SetInputFiles(StringList filenames) { m_InputFilenames = filenames; } const mitk::StringList& mitk::DICOMFileReaderSelector ::GetInputFiles() const { return m_InputFilenames; } mitk::DICOMFileReader::Pointer mitk::DICOMFileReaderSelector ::GetFirstReaderWithMinimumNumberOfOutputImages() { ReaderList workingCandidates; // let all readers analyze the file set unsigned int readerIndex(0); for (ReaderList::iterator rIter = m_Readers.begin(); rIter != m_Readers.end(); ++readerIndex, ++rIter) { (*rIter)->SetInputFiles( m_InputFilenames ); try { (*rIter)->AnalyzeInputFiles(); workingCandidates.push_back( *rIter ); MITK_INFO << "Reader " << readerIndex << " (" << (*rIter)->GetConfigurationLabel() << ") suggests " << (*rIter)->GetNumberOfOutputs() << " 3D blocks"; + if ((*rIter)->GetNumberOfOutputs() == 1) + { + MITK_DEBUG << "Early out with reader #" << readerIndex << " (" << (*rIter)->GetConfigurationLabel() << "), less than 1 block is not possible"; + return *rIter; + } } catch (std::exception& e) { MITK_ERROR << "Reader " << readerIndex << " (" << (*rIter)->GetConfigurationLabel() << ") threw exception during file analysis, ignoring this reader. Exception: " << e.what(); } catch (...) { MITK_ERROR << "Reader " << readerIndex << " (" << (*rIter)->GetConfigurationLabel() << ") threw unknown exception during file analysis, ignoring this reader."; } } DICOMFileReader::Pointer bestReader; unsigned int minimumNumberOfOutputs = std::numeric_limits::max(); readerIndex = 0; unsigned int bestReaderIndex(0); // select the reader with the minimum number of mitk::Images as output for (ReaderList::iterator rIter = workingCandidates.begin(); rIter != workingCandidates.end(); ++readerIndex, ++rIter) { if ( (*rIter)->GetNumberOfOutputs() < minimumNumberOfOutputs ) { minimumNumberOfOutputs = (*rIter)->GetNumberOfOutputs(); bestReader = *rIter; bestReaderIndex = readerIndex; } } MITK_DEBUG << "Decided for reader #" << bestReaderIndex << " (" << bestReader->GetConfigurationLabel() << ")"; MITK_DEBUG << m_PossibleConfigurations[bestReaderIndex]; return bestReader; } diff --git a/Modules/DICOMReader/mitkDICOMFileReaderSelector.h b/Modules/DICOMReader/mitkDICOMFileReaderSelector.h index 712a952265..7133c5a7b3 100644 --- a/Modules/DICOMReader/mitkDICOMFileReaderSelector.h +++ b/Modules/DICOMReader/mitkDICOMFileReaderSelector.h @@ -1,91 +1,95 @@ /*=================================================================== 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 mitkDICOMFileReaderSelector_h #define mitkDICOMFileReaderSelector_h #include "mitkDICOMFileReader.h" +#include + namespace mitk { /** \ingroup DICOMReaderModule \brief Simple best-reader selection. This class implements a process of comparing different DICOMFileReader%s and selecting the reader with the minimal number of mitk::Image%s in its output. The code found in this class can - just be used to select a reader using this simple strategy - be taken as an example of how to use DICOMFileReader%s To create a selection of potential readers, the class makes use of mitk::DICOMReaderConfigurator, i.e. DICOMFileReaderSelector also expects the configuration files/strings to be in the format expected by mitk::DICOMReaderConfigurator. Two convenience methods load "default" configurations from compiled-in resources: LoadBuiltIn3DConfigs() and LoadBuiltIn3DnTConfigs(). */ class DICOMReader_EXPORT DICOMFileReaderSelector : public itk::LightObject { public: mitkClassMacro( DICOMFileReaderSelector, itk::LightObject ) itkNewMacro( DICOMFileReaderSelector ) /// \brief Add a configuration as expected by DICOMReaderConfigurator. /// Configs can only be reset by instantiating a new DICOMFileReaderSelector. void AddConfig(const std::string& xmlDescription); /// \brief Add a configuration as expected by DICOMReaderConfigurator. /// Configs can only be reset by instantiating a new DICOMFileReaderSelector. void AddConfigFile(const std::string& filename); /// \brief Load 3D image creating configurations from the MITK module system (see \ref mitk::Module::FindResources). /// For a default set of configurations, look into the directory Resources of the DICOMReader module. void LoadBuiltIn3DConfigs(); /// \brief Load 3D+t image creating configurations from the MITK module system (see \ref mitk::Module::FindResources). /// For a default set of configurations, look into the directory Resources of the DICOMReader module. void LoadBuiltIn3DnTConfigs(); /// Input files void SetInputFiles(StringList filenames); /// Input files const StringList& GetInputFiles() const; /// Execute the analysis and selection process. The first reader with a minimal number of outputs will be returned. DICOMFileReader::Pointer GetFirstReaderWithMinimumNumberOfOutputImages(); protected: DICOMFileReaderSelector(); virtual ~DICOMFileReaderSelector(); void AddConfigsFromResources(const std::string& path); + void AddConfigFromResource(const std::string& resourcename); + void AddConfigFromResource(ModuleResource& resource); private: StringList m_PossibleConfigurations; StringList m_InputFilenames; typedef std::list ReaderList; ReaderList m_Readers; }; } // namespace #endif // mitkDICOMFileReaderSelector_h diff --git a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx index a4b3bb9796..898060e373 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx @@ -1,297 +1,297 @@ /*=================================================================== 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, itk::GDCMImageIO::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 = itk::GDCMImageIO::New(); typename ReaderType::Pointer reader = ReaderType::New(); reader->SetImageIO(io); - reader->ReverseOrderOff(); + reader->ReverseOrderOff(); // TODO check this.. 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 = FixUpTiltedGeometry( 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; } #define MITK_DEBUG_OUTPUT_FILELIST(list)\ MITK_DEBUG << "-------------------------------------------"; \ for (StringContainer::const_iterator _iter = (list).begin(); _iter!=(list).end(); ++_iter) \ { \ MITK_DEBUG <<" file '" << *_iter<< "'"; \ } \ MITK_DEBUG << "-------------------------------------------"; template mitk::Image::Pointer mitk::ITKDICOMSeriesReaderHelper ::LoadDICOMByITK3DnT( const StringContainerList& filenamesForTimeSteps, bool correctTilt, const GantryTiltInformation& tiltInfo, itk::GDCMImageIO::Pointer& io, Image::Pointer preLoadedImageBlock ) { unsigned int numberOfTimeSteps = filenamesForTimeSteps.size(); mitk::Image::Pointer image = mitk::Image::New(); typedef itk::Image ImageType; typedef itk::ImageSeriesReader ReaderType; io = itk::GDCMImageIO::New(); typename ReaderType::Pointer reader = ReaderType::New(); reader->SetImageIO(io); reader->ReverseOrderOff(); if (preLoadedImageBlock.IsNull()) { unsigned int currentTimeStep = 0; MITK_DEBUG << "Start loading timestep " << currentTimeStep; MITK_DEBUG_OUTPUT_FILELIST( filenamesForTimeSteps.front() ) reader->SetFileNames(filenamesForTimeSteps.front()); 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 = FixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); } image->InitializeByItk(readVolume.GetPointer(), 1, numberOfTimeSteps); image->SetImportVolume(readVolume->GetBufferPointer(), currentTimeStep++); // timestep 0 // for other time-steps for (StringContainerList::const_iterator timestepsIter = ++(filenamesForTimeSteps.begin()); // start with SECOND entry timestepsIter != filenamesForTimeSteps.end(); ++currentTimeStep, ++timestepsIter) { MITK_DEBUG << "Start loading timestep " << currentTimeStep; MITK_DEBUG_OUTPUT_FILELIST( *timestepsIter ) reader->SetFileNames(*timestepsIter); reader->Update(); readVolume = reader->GetOutput(); if (correctTilt) { readVolume = FixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); } image->SetImportVolume(readVolume->GetBufferPointer(), currentTimeStep); } } else { // TODO check and fix image = preLoadedImageBlock; StringContainer fakeList; fakeList.push_back( filenamesForTimeSteps.front().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 ::FixUpTiltedGeometry( 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. 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(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; }