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