Medical Imaging Interaction Toolkit  2018.4.99-389bf124
Medical Imaging Interaction Toolkit
mitkITKDICOMSeriesReaderHelper.txx
Go to the documentation of this file.
1 /*============================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 #include "mitkITKDICOMSeriesReaderHelper.h"
14 
15 #include <itkImageSeriesReader.h>
16 #include <itkResampleImageFilter.h>
17 //#include <itkAffineTransform.h>
18 //#include <itkLinearInterpolateImageFunction.h>
19 //#include <itkTimeProbesCollectorBase.h>
20 
21 #include "dcmtk/ofstd/ofdatime.h"
22 
23 template <typename PixelType>
24 mitk::Image::Pointer
25 mitk::ITKDICOMSeriesReaderHelper
26 ::LoadDICOMByITK(
27  const StringContainer& filenames,
28  bool correctTilt,
29  const GantryTiltInformation& tiltInfo,
30  itk::GDCMImageIO::Pointer& io)
31 {
32  /******** Normal Case, 3D (also for GDCM < 2 usable) ***************/
33  mitk::Image::Pointer image = mitk::Image::New();
34 
35  typedef itk::Image<PixelType, 3> ImageType;
36  typedef itk::ImageSeriesReader<ImageType> ReaderType;
37 
38  io = itk::GDCMImageIO::New();
39  typename ReaderType::Pointer reader = ReaderType::New();
40 
41  reader->SetImageIO(io);
42  reader->ReverseOrderOff(); // at this point we require an order of input images so that
43  // the direction between the origin of the first and the last slice
44  // is the same direction as the image normals! Otherwise we might
45  // see images upside down. Unclear whether this is a bug in MITK,
46  // see NormalDirectionConsistencySorter.
47 
48  reader->SetFileNames(filenames);
49  reader->Update();
50  typename ImageType::Pointer readVolume = reader->GetOutput();
51 
52  // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right position
53  if (correctTilt)
54  {
55  readVolume = FixUpTiltedGeometry( reader->GetOutput(), tiltInfo );
56  }
57 
58  image->InitializeByItk(readVolume.GetPointer());
59  image->SetImportVolume(readVolume->GetBufferPointer());
60 
61 #ifdef MBILOG_ENABLE_DEBUG
62 
63  MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", "
64  << image->GetDimension(1) << ", "
65  << image->GetDimension(2) << "]";
66 
67  MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", "
68  << image->GetGeometry()->GetSpacing()[1] << ", "
69  << image->GetGeometry()->GetSpacing()[2] << "]";
70 #endif //MBILOG_ENABLE_DEBUG
71 
72  return image;
73 }
74 
75 #define MITK_DEBUG_OUTPUT_FILELIST(list)\
76  MITK_DEBUG << "-------------------------------------------"; \
77  for (StringContainer::const_iterator _iter = (list).cbegin(); _iter!=(list).cend(); ++_iter) \
78  { \
79  MITK_DEBUG <<" file '" << *_iter<< "'"; \
80  } \
81  MITK_DEBUG << "-------------------------------------------";
82 
83 template <typename PixelType>
84 mitk::Image::Pointer
85 mitk::ITKDICOMSeriesReaderHelper
86 ::LoadDICOMByITK3DnT(
87  const StringContainerList& filenamesForTimeSteps,
88  bool correctTilt,
89  const GantryTiltInformation& tiltInfo,
90  itk::GDCMImageIO::Pointer& io)
91 {
92  unsigned int numberOfTimeSteps = filenamesForTimeSteps.size();
93 
94  MITK_DEBUG << "Start extracting time bounds of time steps";
95  const TimeBoundsList timeBoundsList = ExtractTimeBoundsOfTimeSteps(filenamesForTimeSteps);
96  if (numberOfTimeSteps!=timeBoundsList.size())
97  {
98  mitkThrow() << "Error while loading 3D+t. Inconsistent size of generated time bounds list. List size: "<< timeBoundsList.size() << "; number of steps: "<<numberOfTimeSteps;
99  }
100 
101  mitk::Image::Pointer image = mitk::Image::New();
102 
103  typedef itk::Image<PixelType, 4> ImageType;
104  typedef itk::ImageSeriesReader<ImageType> ReaderType;
105 
106  io = itk::GDCMImageIO::New();
107  typename ReaderType::Pointer reader = ReaderType::New();
108 
109  reader->SetImageIO(io);
110  reader->ReverseOrderOff(); // at this point we require an order of input images so that
111  // the direction between the origin of the first and the last slice
112  // is the same direction as the image normals! Otherwise we might
113  // see images upside down. Unclear whether this is a bug in MITK,
114  // see NormalDirectionConsistencySorter.
115 
116 
117  unsigned int currentTimeStep = 0;
118 
119 #ifdef MBILOG_ENABLE_DEBUG
120  MITK_DEBUG << "Start loading timestep " << currentTimeStep;
121  MITK_DEBUG_OUTPUT_FILELIST( filenamesForTimeSteps.front() )
122 #endif // MBILOG_ENABLE_DEBUG
123 
124  reader->SetFileNames(filenamesForTimeSteps.front());
125  reader->Update();
126  typename ImageType::Pointer readVolume = reader->GetOutput();
127 
128  // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right position
129  if (correctTilt)
130  {
131  readVolume = FixUpTiltedGeometry( reader->GetOutput(), tiltInfo );
132  }
133 
134  image->InitializeByItk(readVolume.GetPointer(), 1, numberOfTimeSteps);
135  image->SetImportVolume(readVolume->GetBufferPointer(), currentTimeStep++); // timestep 0
136 
137  // for other time-steps
138  for (auto timestepsIter = ++(filenamesForTimeSteps.cbegin()); // start with SECOND entry
139  timestepsIter != filenamesForTimeSteps.cend();
140  ++currentTimeStep, ++timestepsIter)
141  {
142 #ifdef MBILOG_ENABLE_DEBUG
143  MITK_DEBUG << "Start loading timestep " << currentTimeStep;
144  MITK_DEBUG_OUTPUT_FILELIST( *timestepsIter )
145 #endif // MBILOG_ENABLE_DEBUG
146 
147  reader->SetFileNames( *timestepsIter );
148  reader->Update();
149  readVolume = reader->GetOutput();
150 
151  if (correctTilt)
152  {
153  readVolume = FixUpTiltedGeometry( reader->GetOutput(), tiltInfo );
154  }
155 
156  image->SetImportVolume(readVolume->GetBufferPointer(), currentTimeStep);
157  }
158 
159 #ifdef MBILOG_ENABLE_DEBUG
160  MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", "
161  << image->GetDimension(1) << ", "
162  << image->GetDimension(2) << "]";
163 
164  MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", "
165  << image->GetGeometry()->GetSpacing()[1] << ", "
166  << image->GetGeometry()->GetSpacing()[2] << "]";
167 #endif // MBILOG_ENABLE_DEBUG
168 
169  //construct timegeometry
170  TimeGeometry::Pointer timeGeometry = GenerateTimeGeometry(image->GetGeometry(),timeBoundsList);
171  image->SetTimeGeometry(timeGeometry);
172  return image;
173 }
174 
175 
176 template <typename ImageType>
177 typename ImageType::Pointer
178 mitk::ITKDICOMSeriesReaderHelper
179 ::FixUpTiltedGeometry( ImageType* input, const GantryTiltInformation& tiltInfo )
180 {
181  typedef itk::ResampleImageFilter<ImageType,ImageType> ResampleFilterType;
182  typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
183  resampler->SetInput( input );
184 
185  /*
186  Transform for a point is
187  - transform from actual position to index coordinates
188  - apply a shear that undoes the gantry tilt
189  - transform back into world coordinates
190 
191  Anybody who does this in a simpler way: don't forget to write up how and why your solution works
192  */
193  typedef itk::ScalableAffineTransform< double, ImageType::ImageDimension > TransformType;
194  typename TransformType::Pointer transformShear = TransformType::New();
195 
196  /**
197  - apply a shear and spacing correction to the image block that corrects the ITK reader's error
198  - ITK ignores the shear and loads slices into an orthogonal volume
199  - ITK calculates the spacing from the origin distance, which is more than the actual spacing with gantry tilt images
200  - to undo the effect
201  - we have calculated some information in tiltInfo:
202  - the shift in Y direction that is added with each additional slice is the most important information
203  - the Y-shift is calculated in mm world coordinates
204  - we apply a shearing transformation to the ITK-read image volume
205  - to do this locally,
206  - we transform the image volume back to origin and "normal" orientation by applying the inverse of its transform
207  (this brings us into the image's "index coordinate" system)
208  - we apply a shear with the Y-shift factor put into a unit transform at row 1, col 2
209  - we transform the image volume back to its actual position (from index to world coordinates)
210  - we lastly apply modify the image spacing in z direction by replacing this number with the correctly calulcated inter-slice distance
211  */
212 
213  const ScalarType factor = tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() / input->GetSpacing()[1];
214  // row 1, column 2 corrects shear in parallel to Y axis, proportional to distance in Z direction
215  transformShear->Shear( 1, 2, factor );
216 
217  typename TransformType::Pointer imageIndexToWorld = TransformType::New();
218  imageIndexToWorld->SetOffset( input->GetOrigin().GetVectorFromOrigin() );
219 
220  typename TransformType::MatrixType indexToWorldMatrix;
221  indexToWorldMatrix = input->GetDirection();
222 
223  typename ImageType::DirectionType scale;
224  for ( unsigned int i = 0; i < ImageType::ImageDimension; i++ )
225  {
226  scale[i][i] = input->GetSpacing()[i];
227  }
228  indexToWorldMatrix *= scale;
229 
230  imageIndexToWorld->SetMatrix( indexToWorldMatrix );
231 
232  typename TransformType::Pointer imageWorldToIndex = TransformType::New();
233  imageIndexToWorld->GetInverse( imageWorldToIndex );
234 
235  typename TransformType::Pointer gantryTiltCorrection = TransformType::New();
236  gantryTiltCorrection->Compose( imageWorldToIndex );
237  gantryTiltCorrection->Compose( transformShear );
238  gantryTiltCorrection->Compose( imageIndexToWorld );
239 
240  resampler->SetTransform( gantryTiltCorrection );
241 
242  typedef itk::LinearInterpolateImageFunction< ImageType, double > InterpolatorType;
243  typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
244  resampler->SetInterpolator( interpolator );
245  /*
246  This would be the right place to invent a meaningful value for positions outside of the image.
247  For CT, HU -1000 might be meaningful, but a general solution seems not possible. Even for CT,
248  -1000 would only look natural for many not all images.
249  */
250 
251  // TODO use (0028,0120) Pixel Padding Value if present
252  resampler->SetDefaultPixelValue( itk::NumericTraits< typename ImageType::PixelType >::min() );
253 
254  // adjust size in Y direction! (maybe just transform the outer last pixel to see how much space we would need
255 
256  resampler->SetOutputParametersFromImage( input ); // we basically need the same image again, just sheared
257 
258 
259  // if tilt positive, then we need additional pixels BELOW origin, otherwise we need pixels behind the end of the block
260 
261  // in any case we need more size to accomodate shifted slices
262  typename ImageType::SizeType largerSize = resampler->GetSize(); // now the resampler already holds the input image's size.
263  const double imageSizeZ = largerSize[2];
264  //MITK_DEBUG <<"Calculate lager size = " << largerSize[1] << " + " << tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) << " / " << input->GetSpacing()[1] << "+ 2.0";
265  largerSize[1] += static_cast<typename ImageType::SizeType::SizeValueType>(tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) / input->GetSpacing()[1]+ 2.0);
266  resampler->SetSize( largerSize );
267  //MITK_DEBUG << "Fix Y size of image w/ spacing " << input->GetSpacing()[1] << " from " << input->GetLargestPossibleRegion().GetSize()[1] << " to " << largerSize[1];
268 
269  // in SOME cases this additional size is below/behind origin
270  if ( tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() > 0.0 )
271  {
272  typename ImageType::DirectionType imageDirection = input->GetDirection();
273  Vector3D yDirection;
274  yDirection[0] = imageDirection[0][1];
275  yDirection[1] = imageDirection[1][1];
276  yDirection[2] = imageDirection[2][1];
277  yDirection.Normalize();
278 
279  typename ImageType::PointType shiftedOrigin;
280  shiftedOrigin = input->GetOrigin();
281 
282  // add some pixels to make everything fit
283  shiftedOrigin[0] -= yDirection[0] * (tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) + 1.0 * input->GetSpacing()[1]);
284  shiftedOrigin[1] -= yDirection[1] * (tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) + 1.0 * input->GetSpacing()[1]);
285  shiftedOrigin[2] -= yDirection[2] * (tiltInfo.GetTiltCorrectedAdditionalSize(imageSizeZ) + 1.0 * input->GetSpacing()[1]);
286 
287  resampler->SetOutputOrigin( shiftedOrigin );
288  }
289 
290  resampler->Update();
291  typename ImageType::Pointer result = resampler->GetOutput();
292 
293  // ImageSeriesReader calculates z spacing as the distance between the first two origins.
294  // This is not correct in case of gantry tilt, so we set our calculated spacing.
295  typename ImageType::SpacingType correctedSpacing = result->GetSpacing();
296  correctedSpacing[2] = tiltInfo.GetRealZSpacing();
297  result->SetSpacing( correctedSpacing );
298 
299  return result;
300 }
301 
302