diff --git a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h index 0d3626d715..22e6be0984 100644 --- a/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h +++ b/Modules/DICOMReader/mitkDICOMITKSeriesGDCMReader.h @@ -1,83 +1,85 @@ /*=================================================================== 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 3D+t +// TODO providing tags as properties! +// TODO preloaded volumes?? could be solved in a different way.. 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/mitkGantryTiltInformation.h b/Modules/DICOMReader/mitkGantryTiltInformation.h index b43b5a5f48..00a23bdd40 100644 --- a/Modules/DICOMReader/mitkGantryTiltInformation.h +++ b/Modules/DICOMReader/mitkGantryTiltInformation.h @@ -1,129 +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(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 + double GetTiltCorrectedAdditionalSize(unsigned int imageSizeZ) const; /** \brief Calculated tilt angle in degrees. */ double GetTiltAngleInDegrees() const; protected: /** \brief Projection of point p onto line through lineOrigin in direction of lineDirection. */ Point3D projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ); double m_ShiftUp; double m_ShiftRight; double m_ShiftNormal; double m_ITKAssumedSliceSpacing; unsigned int m_NumberOfSlicesApart; }; } // namespace #endif diff --git a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp index 8579d2a6c6..1ffba9b26d 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.cpp @@ -1,96 +1,97 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkITKDICOMSeriesReaderHelper.h" #include "mitkITKDICOMSeriesReaderHelper.txx" mitk::Image::Pointer mitk::ITKDICOMSeriesReaderHelper ::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 + MITK_DEBUG << "Calling LoadDicomSeries with empty filename string container. Probably invalid application logic."; return NULL; // this is not actually an error but the result is very simple } + typedef itk::GDCMImageIO DcmIoType; DcmIoType::Pointer io = DcmIoType::New(); 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 b5c0d53ab8..c1042cf29f 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.h @@ -1,56 +1,53 @@ /*=================================================================== 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, 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, + itk::GDCMImageIO::Pointer& io, Image::Pointer preLoadedImageBlock ); }; } #endif diff --git a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx index 172c85209b..df477f5359 100644 --- a/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx +++ b/Modules/DICOMReader/mitkITKDICOMSeriesReaderHelper.txx @@ -1,209 +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, + 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 = DcmIoType::New(); + io = itk::GDCMImageIO::New(); typename ReaderType::Pointer reader = ReaderType::New(); reader->SetImageIO(io); reader->ReverseOrderOff(); if (preLoadedImageBlock.IsNull()) { reader->SetFileNames(filenames); reader->Update(); typename ImageType::Pointer readVolume = reader->GetOutput(); // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right position if (correctTilt) { readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); } image->InitializeByItk(readVolume.GetPointer()); image->SetImportVolume(readVolume->GetBufferPointer()); } else { image = preLoadedImageBlock; StringContainer fakeList; fakeList.push_back( filenames.front() ); reader->SetFileNames( fakeList ); // we always need to load at least one file to get the MetaDataDictionary reader->Update(); } MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " << image->GetDimension(1) << ", " << image->GetDimension(2) << "]"; MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " << image->GetGeometry()->GetSpacing()[1] << ", " << image->GetGeometry()->GetSpacing()[2] << "]"; return image; } template typename ImageType::Pointer mitk::ITKDICOMSeriesReaderHelper ::InPlaceFixUpTiltedGeometry( ImageType* input, const GantryTiltInformation& tiltInfo ) { 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. 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; }