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