Medical Imaging Interaction Toolkit  2018.4.99-389bf124
Medical Imaging Interaction Toolkit
mitkDicomSeriesReader.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 #ifndef MITKDICOMSERIESREADER_TXX_
14 #define MITKDICOMSERIESREADER_TXX_
15 
16 #include <legacy/mitkDicomSeriesReader.h>
17 
18 #include <itkImageSeriesReader.h>
19 #include <mitkProperties.h>
20 
21 #include <itkAffineTransform.h>
22 #include <itkLinearInterpolateImageFunction.h>
23 #include <itkResampleImageFilter.h>
24 #include <itkTimeProbesCollectorBase.h>
25 
26 #include <limits>
27 
28 #include <mitkImage.h>
29 
30 namespace mitk
31 {
32  template <typename PixelType>
33  Image::Pointer DicomSeriesReader::LoadDICOMByITK4D(std::list<StringContainer> &imageBlocks,
34  ImageBlockDescriptor imageBlockDescriptor,
35  bool correctTilt,
36  const GantryTiltInformation &tiltInfo,
37  DcmIoType::Pointer &io,
38  CallbackCommand *command,
39  Image::Pointer preLoadedImageBlock)
40  {
41  mitk::Image::Pointer image = mitk::Image::New();
42 
43  // It is 3D+t! Read it and store into mitk image
44  typedef itk::Image<PixelType, 4> ImageType;
45  typedef itk::ImageSeriesReader<ImageType> ReaderType;
46 
47  typename ReaderType::Pointer reader = ReaderType::New();
48 
49  reader->SetImageIO(io);
50  reader->ReverseOrderOff();
51 
52  if (command)
53  {
54  reader->AddObserver(itk::ProgressEvent(), command);
55  }
56 
57  if (preLoadedImageBlock.IsNull())
58  {
59  reader->SetFileNames(imageBlocks.front());
60 
61  reader->Update();
62 
63  typename ImageType::Pointer readVolume = reader->GetOutput();
64  // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right
65  // position
66  if (correctTilt)
67  {
68  readVolume = InPlaceFixUpTiltedGeometry(reader->GetOutput(), tiltInfo);
69  }
70 
71  unsigned int volume_count = imageBlocks.size();
72  image->InitializeByItk(readVolume.GetPointer(), 1, volume_count);
73  image->SetImportVolume(readVolume->GetBufferPointer(), 0u);
74 
75  FixSpacingInformation(image, imageBlockDescriptor);
76  }
77  else
78  {
79  StringContainer fakeList;
80  fakeList.push_back(imageBlocks.front().front());
81  reader->SetFileNames(fakeList); // only ONE first filename to get MetaDataDictionary
82 
83  image = preLoadedImageBlock;
84  }
85 
86  gdcm::Scanner scanner;
87  ScanForSliceInformation(imageBlockDescriptor.GetFilenames(), scanner);
88  CopyMetaDataToImageProperties(imageBlocks, scanner.GetMappings(), io, imageBlockDescriptor, image);
89 
90  MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " << image->GetDimension(1) << ", "
91  << image->GetDimension(2) << ", " << image->GetDimension(3) << "]";
92 
93  MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", "
94  << image->GetGeometry()->GetSpacing()[1] << ", " << image->GetGeometry()->GetSpacing()[2] << "]";
95 
96  if (preLoadedImageBlock.IsNull())
97  {
98  unsigned int act_volume = 1u;
99  for (auto df_it = ++imageBlocks.begin(); df_it != imageBlocks.end(); ++df_it)
100  {
101  reader->SetFileNames(*df_it);
102  reader->Update();
103  typename ImageType::Pointer readVolume = reader->GetOutput();
104 
105  if (correctTilt)
106  {
107  readVolume = InPlaceFixUpTiltedGeometry(reader->GetOutput(), tiltInfo);
108  }
109 
110  image->SetImportVolume(readVolume->GetBufferPointer(), act_volume++);
111  }
112  }
113 
114  return image;
115  }
116 
117  template <typename PixelType>
118  Image::Pointer DicomSeriesReader::LoadDICOMByITK(const StringContainer &filenames,
119  bool correctTilt,
120  const GantryTiltInformation &tiltInfo,
121  DcmIoType::Pointer &io,
122  CallbackCommand *command,
123  Image::Pointer preLoadedImageBlock)
124  {
125  /******** Normal Case, 3D (also for GDCM < 2 usable) ***************/
126  mitk::Image::Pointer image = mitk::Image::New();
127 
128  typedef itk::Image<PixelType, 3> ImageType;
129  typedef itk::ImageSeriesReader<ImageType> ReaderType;
130 
131  io = DcmIoType::New();
132  typename ReaderType::Pointer reader = ReaderType::New();
133 
134  reader->SetImageIO(io);
135  reader->ReverseOrderOff();
136 
137  if (command)
138  {
139  reader->AddObserver(itk::ProgressEvent(), command);
140  }
141 
142  if (preLoadedImageBlock.IsNull())
143  {
144  reader->SetFileNames(filenames);
145  reader->Update();
146  typename ImageType::Pointer readVolume = reader->GetOutput();
147 
148  // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right
149  // position
150  if (correctTilt)
151  {
152  readVolume = InPlaceFixUpTiltedGeometry(reader->GetOutput(), tiltInfo);
153  }
154 
155  image->InitializeByItk(readVolume.GetPointer());
156  image->SetImportVolume(readVolume->GetBufferPointer());
157  }
158  else
159  {
160  image = preLoadedImageBlock;
161  StringContainer fakeList;
162  fakeList.push_back(filenames.front());
163  reader->SetFileNames(fakeList); // we always need to load at least one file to get the MetaDataDictionary
164  reader->Update();
165  }
166 
167  MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " << image->GetDimension(1) << ", "
168  << image->GetDimension(2) << "]";
169 
170  MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", "
171  << image->GetGeometry()->GetSpacing()[1] << ", " << image->GetGeometry()->GetSpacing()[2] << "]";
172 
173  return image;
174  }
175 
176  template <typename ImageType>
177  typename ImageType::Pointer DicomSeriesReader::InPlaceFixUpTiltedGeometry(ImageType *input,
178  const GantryTiltInformation &tiltInfo)
179  {
180  typedef itk::ResampleImageFilter<ImageType, ImageType> ResampleFilterType;
181  typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
182  resampler->SetInput(input);
183 
184  /*
185  Transform for a point is
186  - transform from actual position to index coordinates
187  - apply a shear that undoes the gantry tilt
188  - transform back into world coordinates
189 
190  Anybody who does this in a simpler way: don't forget to write up how and why your solution works
191  */
192  typedef itk::ScalableAffineTransform<double, ImageType::ImageDimension> TransformType;
193  typename TransformType::Pointer transformShear = TransformType::New();
194 
195  /**
196  - apply a shear and spacing correction to the image block that corrects the ITK reader's error
197  - ITK ignores the shear and loads slices into an orthogonal volume
198  - ITK calculates the spacing from the origin distance, which is more than the actual spacing with gantry tilt
199  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
207  transform
208  (this brings us into the image's "index coordinate" system)
209  - we apply a shear with the Y-shift factor put into a unit transform at row 1, col 2
210  - we transform the image volume back to its actual position (from index to world coordinates)
211  - we lastly apply modify the image spacing in z direction by replacing this number with the correctly calulcated
212  inter-slice distance
213  */
214 
215  ScalarType factor = tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() / input->GetSpacing()[1];
216  // row 1, column 2 corrects shear in parallel to Y axis, proportional to distance in Z direction
217  transformShear->Shear(1, 2, factor);
218 
219  typename TransformType::Pointer imageIndexToWorld = TransformType::New();
220  imageIndexToWorld->SetOffset(input->GetOrigin().GetVectorFromOrigin());
221 
222  typename TransformType::MatrixType indexToWorldMatrix;
223  indexToWorldMatrix = input->GetDirection();
224 
225  typename ImageType::DirectionType scale;
226  for (unsigned int i = 0; i < ImageType::ImageDimension; i++)
227  {
228  scale[i][i] = input->GetSpacing()[i];
229  }
230  indexToWorldMatrix *= scale;
231 
232  imageIndexToWorld->SetMatrix(indexToWorldMatrix);
233 
234  typename TransformType::Pointer imageWorldToIndex = TransformType::New();
235  imageIndexToWorld->GetInverse(imageWorldToIndex);
236 
237  typename TransformType::Pointer gantryTiltCorrection = TransformType::New();
238  gantryTiltCorrection->Compose(imageWorldToIndex);
239  gantryTiltCorrection->Compose(transformShear);
240  gantryTiltCorrection->Compose(imageIndexToWorld);
241 
242  resampler->SetTransform(gantryTiltCorrection);
243 
244  typedef itk::LinearInterpolateImageFunction<ImageType, double> InterpolatorType;
245  typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
246  resampler->SetInterpolator(interpolator);
247  /*
248  This would be the right place to invent a meaningful value for positions outside of the image.
249  For CT, HU -1000 might be meaningful, but a general solution seems not possible. Even for CT,
250  -1000 would only look natural for many not all images.
251  */
252 
253  // TODO use (0028,0120) Pixel Padding Value if present
254  resampler->SetDefaultPixelValue(itk::NumericTraits<typename ImageType::PixelType>::min());
255 
256  // adjust size in Y direction! (maybe just transform the outer last pixel to see how much space we would need
257 
258  resampler->SetOutputParametersFromImage(input); // we basically need the same image again, just sheared
259 
260  // if tilt positive, then we need additional pixels BELOW origin, otherwise we need pixels behind the end of the
261  // block
262 
263  // in any case we need more size to accomodate shifted slices
264  typename ImageType::SizeType largerSize =
265  resampler->GetSize(); // now the resampler already holds the input image's size.
266  largerSize[1] += static_cast<typename ImageType::SizeType::SizeValueType>(
267  tiltInfo.GetTiltCorrectedAdditionalSize() / input->GetSpacing()[1] + 2.0);
268  resampler->SetSize(largerSize);
269 
270  // in SOME cases this additional size is below/behind origin
271  if (tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() > 0.0)
272  {
273  typename ImageType::DirectionType imageDirection = input->GetDirection();
274  Vector3D yDirection;
275  yDirection[0] = imageDirection[0][1];
276  yDirection[1] = imageDirection[1][1];
277  yDirection[2] = imageDirection[2][1];
278  yDirection.Normalize();
279 
280  typename ImageType::PointType shiftedOrigin;
281  shiftedOrigin = input->GetOrigin();
282 
283  // add some pixels to make everything fit
284  shiftedOrigin[0] -= yDirection[0] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]);
285  shiftedOrigin[1] -= yDirection[1] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]);
286  shiftedOrigin[2] -= yDirection[2] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]);
287 
288  resampler->SetOutputOrigin(shiftedOrigin);
289  }
290 
291  resampler->Update();
292  typename ImageType::Pointer result = resampler->GetOutput();
293 
294  // ImageSeriesReader calculates z spacing as the distance between the first two origins.
295  // This is not correct in case of gantry tilt, so we set our calculated spacing.
296  typename ImageType::SpacingType correctedSpacing = result->GetSpacing();
297  correctedSpacing[2] = tiltInfo.GetRealZSpacing();
298  result->SetSpacing(correctedSpacing);
299 
300  return result;
301  }
302 }
303 
304 #endif