Medical Imaging Interaction Toolkit  2018.4.99-bd7b41ba
Medical Imaging Interaction Toolkit
mitkImageMappingHelper.cpp
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 <itkInterpolateImageFunction.h>
14 #include <itkNearestNeighborInterpolateImageFunction.h>
15 #include <itkLinearInterpolateImageFunction.h>
16 #include <itkBSplineInterpolateImageFunction.h>
17 #include <itkWindowedSincInterpolateImageFunction.h>
18 
19 #include <mitkImageAccessByItk.h>
20 #include <mitkImageCast.h>
21 #include <mitkGeometry3D.h>
22 #include <mitkImageToItk.h>
23 #include <mitkImageTimeSelector.h>
24 
25 #include "mapRegistration.h"
26 
27 #include "mitkImageMappingHelper.h"
28 #include "mitkRegistrationHelper.h"
29 
30 template <typename TImage >
31 typename ::itk::InterpolateImageFunction< TImage >::Pointer generateInterpolator(mitk::ImageMappingInterpolator::Type interpolatorType)
32 {
33  typedef ::itk::InterpolateImageFunction< TImage > BaseInterpolatorType;
34  typename BaseInterpolatorType::Pointer result;
35 
36  switch (interpolatorType)
37  {
39  {
40  result = ::itk::NearestNeighborInterpolateImageFunction<TImage>::New();
41  break;
42  }
44  {
45  typename ::itk::BSplineInterpolateImageFunction<TImage>::Pointer spInterpolator = ::itk::BSplineInterpolateImageFunction<TImage>::New();
46  spInterpolator->SetSplineOrder(3);
47  result = spInterpolator;
48  break;
49  }
51  {
52  result = ::itk::WindowedSincInterpolateImageFunction<TImage,4>::New();
53  break;
54  }
56  {
57  result = ::itk::WindowedSincInterpolateImageFunction<TImage,4,::itk::Function::WelchWindowFunction<4> >::New();
58  break;
59  }
60  default:
61  {
62  result = ::itk::LinearInterpolateImageFunction<TImage>::New();
63  break;
64  }
65 
66  }
67 
68  return result;
69 };
70 
71 template <typename TPixelType, unsigned int VImageDimension >
72 void doMITKMap(const ::itk::Image<TPixelType,VImageDimension>* input, mitk::ImageMappingHelper::ResultImageType::Pointer& result, const mitk::ImageMappingHelper::RegistrationType*& registration,
73  bool throwOnOutOfInputAreaError, const double& paddingValue, const mitk::ImageMappingHelper::ResultImageGeometryType*& resultGeometry,
74  bool throwOnMappingError, const double& errorValue, mitk::ImageMappingInterpolator::Type interpolatorType)
75 {
76  typedef ::map::core::Registration<VImageDimension,VImageDimension> ConcreteRegistrationType;
77  typedef ::map::core::ImageMappingTask<ConcreteRegistrationType, ::itk::Image<TPixelType,VImageDimension>, ::itk::Image<TPixelType,VImageDimension> > MappingTaskType;
78  typename MappingTaskType::Pointer spTask = MappingTaskType::New();
79 
80  typedef typename MappingTaskType::ResultImageDescriptorType ResultImageDescriptorType;
81  typename ResultImageDescriptorType::Pointer resultDescriptor;
82 
83  //check if image and result geometry fits the passed registration
85  if (registration->getMovingDimensions()!=VImageDimension)
86  {
87  map::core::OStringStream str;
88  str << "Dimension of MITK image ("<<VImageDimension<<") does not equal the moving dimension of the registration object ("<<registration->getMovingDimensions()<<").";
89  throw mitk::AccessByItkException(str.str());
90  }
91 
92  if (registration->getTargetDimensions()!=VImageDimension)
93  {
94  map::core::OStringStream str;
95  str << "Dimension of MITK image ("<<VImageDimension<<") does not equal the target dimension of the registration object ("<<registration->getTargetDimensions()<<").";
96  throw mitk::AccessByItkException(str.str());
97  }
98 
99  const ConcreteRegistrationType* castedReg = dynamic_cast<const ConcreteRegistrationType*>(registration);
100 
101  if (registration->getTargetDimensions()==2 && resultGeometry)
102  {
104 
105  if (bounds[4]!=0 || bounds[5]!=0)
106  {
107  //array "bounds" is constructed as [min Dim1, max Dim1, min Dim2, max Dim2, min Dim3, max Dim3]
108  //therfore [4] and [5] must be 0
109 
110  map::core::OStringStream str;
111  str << "Dimension of defined result geometry does not equal the target dimension of the registration object ("<<registration->getTargetDimensions()<<").";
112  throw mitk::AccessByItkException(str.str());
113  }
114  }
115 
116  //check/create resultDescriptor
118  if (resultGeometry)
119  {
120  resultDescriptor = ResultImageDescriptorType::New();
121 
122  typename ResultImageDescriptorType::PointType origin;
123  typename ResultImageDescriptorType::SizeType size;
124  typename ResultImageDescriptorType::SpacingType fieldSpacing;
125  typename ResultImageDescriptorType::DirectionType matrix;
126 
128  mitk::Vector3D geoSpacing = resultGeometry->GetSpacing();
129  mitk::Point3D geoOrigin = resultGeometry->GetOrigin();
130  mitk::AffineTransform3D::MatrixType geoMatrix = resultGeometry->GetIndexToWorldTransform()->GetMatrix();
131 
132  for (unsigned int i = 0; i<VImageDimension; ++i)
133  {
134  origin[i] = static_cast<typename ResultImageDescriptorType::PointType::ValueType>(geoOrigin[i]);
135  fieldSpacing[i] = static_cast<typename ResultImageDescriptorType::SpacingType::ValueType>(geoSpacing[i]);
136  size[i] = static_cast<typename ResultImageDescriptorType::SizeType::SizeValueType>(geoBounds[(2*i)+1]-geoBounds[2*i])*fieldSpacing[i];
137  }
138 
139  //Matrix extraction
140  matrix.SetIdentity();
141  unsigned int i;
142  unsigned int j;
143 
149 
150  // the following loop devides by spacing now to normalize columns.
151  // counterpart of InitializeByItk in mitkImage.h line 372 of revision 15092.
152 
153  // Check if information is lost
154  if ( VImageDimension == 2)
155  {
156  if ( ( geoMatrix[0][2] != 0) ||
157  ( geoMatrix[1][2] != 0) ||
158  ( geoMatrix[2][0] != 0) ||
159  ( geoMatrix[2][1] != 0) ||
160  (( geoMatrix[2][2] != 1) && ( geoMatrix[2][2] != -1) ))
161  {
162  // The 2D MITK image contains 3D rotation information.
163  // This cannot be expressed in a 2D ITK image, so the ITK image will have no rotation
164  }
165  else
166  {
167  // The 2D MITK image can be converted to an 2D ITK image without information loss!
168  for ( i=0; i < 2; ++i)
169  {
170  for( j=0; j < 2; ++j )
171  {
172  matrix[i][j] = geoMatrix[i][j]/fieldSpacing[j];
173  }
174  }
175  }
176  }
177  else if (VImageDimension == 3)
178  {
179  // Normal 3D image. Conversion possible without problem!
180  for ( i=0; i < 3; ++i)
181  {
182  for( j=0; j < 3; ++j )
183  {
184  matrix[i][j] = geoMatrix[i][j]/fieldSpacing[j];
185  }
186  }
187  }
188  else
189  {
190  assert(0);
191  throw mitk::AccessByItkException("Usage of resultGeometry for 2D images is not yet implemented.");
195  }
196 
197  resultDescriptor->setOrigin(origin);
198  resultDescriptor->setSize(size);
199  resultDescriptor->setSpacing(fieldSpacing);
200  resultDescriptor->setDirection(matrix);
201  }
202 
203  //do the mapping
205  typedef ::itk::InterpolateImageFunction< ::itk::Image<TPixelType,VImageDimension> > BaseInterpolatorType;
206  typename BaseInterpolatorType::Pointer interpolator = generateInterpolator< ::itk::Image<TPixelType,VImageDimension> >(interpolatorType);
207  assert(interpolator.IsNotNull());
208  spTask->setImageInterpolator(interpolator);
209  spTask->setInputImage(input);
210  spTask->setRegistration(castedReg);
211  spTask->setResultImageDescriptor(resultDescriptor);
212  spTask->setThrowOnMappingError(throwOnMappingError);
213  spTask->setErrorValue(errorValue);
214  spTask->setThrowOnPaddingError(throwOnOutOfInputAreaError);
215  spTask->setPaddingValue(paddingValue);
216 
217  spTask->execute();
218  mitk::CastToMitkImage<>(spTask->getResultImage(),result);
219 }
220 
223  bool throwOnOutOfInputAreaError, const double& paddingValue, const ResultImageGeometryType* resultGeometry,
224  bool throwOnMappingError, const double& errorValue, mitk::ImageMappingInterpolator::Type interpolatorType)
225 {
226  if (!registration)
227  {
228  mitkThrow() << "Cannot map image. Passed registration wrapper pointer is nullptr.";
229  }
230  if (!input)
231  {
232  mitkThrow() << "Cannot map image. Passed image pointer is nullptr.";
233  }
234 
236 
237  if(input->GetTimeSteps()==1)
238  { //map the image and done
239  AccessByItk_n(input, doMITKMap, (result, registration, throwOnOutOfInputAreaError, paddingValue, resultGeometry, throwOnMappingError, errorValue, interpolatorType));
240  }
241  else
242  { //map every time step and compose
243 
244  mitk::TimeGeometry::ConstPointer timeGeometry = input->GetTimeGeometry();
245  mitk::TimeGeometry::Pointer mappedTimeGeometry = timeGeometry->Clone();
246 
247  for (unsigned int i = 0; i<input->GetTimeSteps(); ++i)
248  {
249  ResultImageGeometryType::Pointer mappedGeometry = resultGeometry->Clone();
250  mappedTimeGeometry->SetTimeStepGeometry(mappedGeometry,i);
251  }
252 
253  result = mitk::Image::New();
254  result->Initialize(input->GetPixelType(),*mappedTimeGeometry, 1, input->GetTimeSteps());
255 
256  for (unsigned int i = 0; i<input->GetTimeSteps(); ++i)
257  {
259  imageTimeSelector->SetInput(input);
260  imageTimeSelector->SetTimeNr(i);
261  imageTimeSelector->UpdateLargestPossibleRegion();
262 
263  InputImageType::Pointer timeStepInput = imageTimeSelector->GetOutput();
264  ResultImageType::Pointer timeStepResult;
265  AccessByItk_n(timeStepInput, doMITKMap, (timeStepResult, registration, throwOnOutOfInputAreaError, paddingValue, resultGeometry, throwOnMappingError, errorValue, interpolatorType));
266  mitk::ImageReadAccessor readAccess(timeStepResult);
267  result->SetVolume(readAccess.GetData(),i);
268  }
269  }
270 
271  return result;
272 }
273 
276  bool throwOnOutOfInputAreaError, const double& paddingValue, const ResultImageGeometryType* resultGeometry,
277  bool throwOnMappingError, const double& errorValue, mitk::ImageMappingInterpolator::Type)
278 {
279  if (!registration)
280  {
281  mitkThrow() << "Cannot map image. Passed registration wrapper pointer is nullptr.";
282  }
283  if (!registration->GetRegistration())
284  {
285  mitkThrow() << "Cannot map image. Passed registration wrapper containes no registration.";
286  }
287  if (!input)
288  {
289  mitkThrow() << "Cannot map image. Passed image pointer is nullptr.";
290  }
291 
292  ResultImageType::Pointer result = map(input, registration->GetRegistration(), throwOnOutOfInputAreaError, paddingValue, resultGeometry, throwOnMappingError, errorValue);
293  return result;
294 }
295 
296 
299  refineGeometry(const InputImageType* input, const RegistrationType* registration,
300  bool throwOnError)
301 {
303 
304  if (!registration)
305  {
306  mitkThrow() << "Cannot refine image geometry. Passed registration pointer is nullptr.";
307  }
308  if (!input)
309  {
310  mitkThrow() << "Cannot refine image geometry. Passed image pointer is nullptr.";
311  }
312 
313  mitk::MITKRegistrationHelper::Affine3DTransformType::Pointer spTransform = mitk::MITKRegistrationHelper::getAffineMatrix(registration,false);
314  if(spTransform.IsNull() && throwOnError)
315  {
316  mitkThrow() << "Cannot refine image geometry. Registration does not contain a suitable direct mapping kernel (3D affine transformation or compatible required).";
317  }
318 
319  if(spTransform.IsNotNull())
320  {
321  //copy input image
322  result = input->Clone();
323 
324  //refine geometries
325  for(unsigned int i = 0; i < result->GetTimeSteps(); ++i)
326  { //refine every time step
327  result->GetGeometry(i)->Compose(spTransform);
328  }
329  result->GetTimeGeometry()->Update();
330  }
331 
332  return result;
333 }
334 
337  refineGeometry(const InputImageType* input, const MITKRegistrationType* registration,
338  bool throwOnError)
339 {
340  if (!registration)
341  {
342  mitkThrow() << "Cannot refine image geometry. Passed registration wrapper pointer is nullptr.";
343  }
344  if (!registration->GetRegistration())
345  {
346  mitkThrow() << "Cannot refine image geometry. Passed registration wrapper containes no registration.";
347  }
348  if (!input)
349  {
350  mitkThrow() << "Cannot refine image geometry. Passed image pointer is nullptr.";
351  }
352 
353  ResultImageType::Pointer result = refineGeometry(input, registration->GetRegistration(), throwOnError);
354  return result;
355 }
356 
357 bool
360 {
361  bool result = true;
362 
363  if (!registration)
364  {
365  mitkThrow() << "Cannot check refine capability of registration. Passed registration pointer is nullptr.";
366  }
367 
368  //if the helper does not return null, we can refine the geometry.
369  result = mitk::MITKRegistrationHelper::getAffineMatrix(registration,false).IsNotNull();
370 
371  return result;
372 }
373 
374 bool
377 {
378  if (!registration)
379  {
380  mitkThrow() << "Cannot check refine capability of registration. Passed registration wrapper pointer is nullptr.";
381  }
382  if (!registration->GetRegistration())
383  {
384  mitkThrow() << "Cannot check refine capability of registration. Passed registration wrapper containes no registration.";
385  }
386 
387  return canRefineGeometry(registration->GetRegistration());
388 }
389 
Pointer Clone() const
static Affine3DTransformType::Pointer getAffineMatrix(const mitk::MAPRegistrationWrapper *wrapper, bool inverseKernel)
BoundingBoxType::BoundsArrayType BoundsArrayType
const mitk::PixelType GetPixelType(int n=0) const
Returns the PixelType of channel n.
Definition: mitkImage.cpp:101
Pointer Clone() const
MAPRegistrationWrapper Wrapper class to allow the handling of MatchPoint registration objects as mitk...
#define AccessByItk_n(mitkImage, itkImageTypeFunction, va_tuple)
Access a MITK image by an ITK image with one or more parameters.
const mitk::TimeGeometry * GetTimeGeometry() const
Return the TimeGeometry of the data as const pointer.
Definition: mitkBaseData.h:66
::map::core::RegistrationBase * GetRegistration()
MITKMATCHPOINTREGISTRATION_EXPORT bool canRefineGeometry(const RegistrationType *registration)
ValueType
Type of the value held by a Value object.
Definition: jsoncpp.h:345
#define mitkThrow()
Image class for storing images.
Definition: mitkImage.h:72
MITKMATCHPOINTREGISTRATION_EXPORT ResultImageType::Pointer refineGeometry(const InputImageType *input, const RegistrationType *registration, bool throwOnError=true)
void doMITKMap(const ::itk::Image< TPixelType, VImageDimension > *input, mitk::ImageMappingHelper::ResultImageType::Pointer &result, const mitk::ImageMappingHelper::RegistrationType *&registration, bool throwOnOutOfInputAreaError, const double &paddingValue, const mitk::ImageMappingHelper::ResultImageGeometryType *&resultGeometry, bool throwOnMappingError, const double &errorValue, mitk::ImageMappingInterpolator::Type interpolatorType)
const Point3D GetOrigin() const
Get the origin, e.g. the upper-left corner of the plane.
static Pointer New()
MITKMATCHPOINTREGISTRATION_EXPORT ResultImageType::Pointer map(const InputImageType *input, const RegistrationType *registration, bool throwOnOutOfInputAreaError=false, const double &paddingValue=0, const ResultImageGeometryType *resultGeometry=nullptr, bool throwOnMappingError=true, const double &errorValue=0, mitk::ImageMappingInterpolator::Type interpolatorType=mitk::ImageMappingInterpolator::Linear)
Exception class thrown in AccessByItk macros.
const mitk::Vector3D GetSpacing() const
Get the spacing (size of a pixel).
unsigned int GetTimeSteps() const
Get the number of time steps from the TimeGeometry As the base data has not a data vector given by it...
Definition: mitkBaseData.h:360
ImageReadAccessor class to get locked read access for a particular image part.
typename ::itk::InterpolateImageFunction< TImage >::Pointer generateInterpolator(mitk::ImageMappingInterpolator::Type interpolatorType)
::map::core::RegistrationBase RegistrationType
const BoundsArrayType GetBounds() const
BaseGeometry Describes the geometry of a data object.
const void * GetData() const
Gives const access to the data.
static Pointer New()
mitk::AffineTransform3D * GetIndexToWorldTransform()
Get the transformation used to convert from index to world coordinates.