ANIMA  4.0
animaReorientation.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <itkImage.h>
4 #include <itkOrientImageFilter.h>
6 #include <itkExtractImageFilter.h>
7 
8 namespace anima
9 {
10 template <class ImageType>
11 typename itk::SmartPointer<ImageType>
12 reorient3DImage(typename itk::SmartPointer<ImageType> input,
13  itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
14 {
15  // Type redefinition and cast are only here for compilation purpose.
16  typedef itk::Image<typename ImageType::PixelType, 3> CheckedImageType;
17  typename CheckedImageType::Pointer checkedInput;
18  checkedInput = dynamic_cast<CheckedImageType *>(input.GetPointer());
19 
20  typedef itk::OrientImageFilter<CheckedImageType, CheckedImageType> OrientImageFilter;
21  typename OrientImageFilter::Pointer orienter = OrientImageFilter::New();
22  orienter->UseImageDirectionOn();
23  orienter->SetDesiredCoordinateOrientation(orientation);
24  orienter->SetInput(checkedInput);
25  orienter->Update();
26  return dynamic_cast<ImageType *>(orienter->GetOutput());
27 }
28 
29 template <class ImageType>
30 typename itk::SmartPointer<ImageType>
31 reorientImage(typename itk::SmartPointer<ImageType> input,
32  itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
33 {
34  if(ImageType::ImageDimension == 4)
35  {
36  typedef itk::Image<typename ImageType::PixelType, 3> ImageToReorientType;
37  typedef itk::ExtractImageFilter<ImageType, ImageToReorientType> ExtractImageFilterType;
38  typename ExtractImageFilterType::Pointer extractFilter = ExtractImageFilterType::New();
39  extractFilter->SetInput(input);
40  extractFilter->SetDirectionCollapseToGuess();
41 
42  typename ImageType::IndexType extractIndex;
43  typename ImageType::SizeType extractSize;
44  typename ImageType::RegionType extractRegion;
45 
46  for (unsigned int d = 0; d < 3; ++d)
47  {
48  extractIndex[d] = 0;
49  extractSize[d] = input->GetLargestPossibleRegion().GetSize()[d];
50  }
51  extractIndex[3] = 0;
52  extractRegion.SetIndex(extractIndex);
53  extractSize[3] = 0;
54  extractRegion.SetSize(extractSize);
55  extractFilter->SetExtractionRegion(extractRegion);
56 
57  extractFilter->Update();
58 
59  typename ImageToReorientType::Pointer reorientedBase;
60  reorientedBase = reorient3DImage<ImageToReorientType>(extractFilter->GetOutput(), orientation);
61  typename ImageType::PointType reorientedOrigin;
62  typename ImageType::SpacingType reorientedSpacing;
63  typename ImageType::DirectionType reorientedDirection;
64  typename ImageType::IndexType reorientedStart;
65  typename ImageType::SizeType reorientedEnd;
66 
67 
68  for(unsigned int d = 0; d < 3; ++d)
69  {
70  reorientedOrigin[d] = reorientedBase->GetOrigin()[d];
71  reorientedSpacing[d] = reorientedBase->GetSpacing()[d];
72  for(unsigned int e = 0; e < 3; ++e)
73  reorientedDirection[d][e] = reorientedBase->GetDirection()[d][e];
74  reorientedDirection[d][3] = input->GetDirection()[d][3];
75 
76  reorientedStart[d] = reorientedBase->GetLargestPossibleRegion().GetIndex()[d];
77  reorientedEnd[d] = reorientedBase->GetLargestPossibleRegion().GetSize()[d];
78  }
79  reorientedOrigin[3] = input->GetOrigin()[3];
80  reorientedSpacing[3] = input->GetSpacing()[3];
81  for(unsigned int f = 0; f < 4; ++f)
82  reorientedDirection[3][f] = input->GetDirection()[3][f];
83  reorientedStart[3] = input->GetLargestPossibleRegion().GetIndex()[3];
84  reorientedEnd[3] = input->GetLargestPossibleRegion().GetSize()[3];
85 
86  typename ImageType::Pointer reorientedImage = ImageType::New();
87  reorientedImage->SetRegions(typename ImageType::RegionType(reorientedStart, reorientedEnd));
88  reorientedImage->SetOrigin(reorientedOrigin);
89  reorientedImage->SetSpacing(reorientedSpacing);
90  reorientedImage->SetDirection(reorientedDirection);
91  reorientedImage->Allocate();
92 
93  typedef itk::ImageRegionIterator <ImageType> FillIteratorType;
94  FillIteratorType fillItr(reorientedImage, reorientedImage->GetLargestPossibleRegion());
95 
96  typedef itk::ImageRegionConstIterator <ImageToReorientType> ReorientedIteratorType;
97  ReorientedIteratorType reoItr(reorientedBase, reorientedBase->GetLargestPossibleRegion());
98  while (!reoItr.IsAtEnd())
99  {
100  fillItr.Set(reoItr.Get());
101  ++fillItr;
102  ++reoItr;
103  }
104 
105  for(unsigned int i = 1; i < input->GetLargestPossibleRegion().GetSize()[3]; ++i)
106  {
107  extractIndex[3] = i;
108  extractRegion.SetIndex(extractIndex);
109  extractFilter->SetExtractionRegion(extractRegion);
110  extractFilter->SetDirectionCollapseToGuess();
111  extractFilter->Update();
112 
113  reorientedBase = reorient3DImage<ImageToReorientType>(extractFilter->GetOutput(), orientation);
114  reoItr = ReorientedIteratorType(reorientedBase, reorientedBase->GetLargestPossibleRegion());
115  while (!reoItr.IsAtEnd())
116  {
117  fillItr.Set(reoItr.Get());
118  ++fillItr;
119  ++reoItr;
120  }
121  }
122  return reorientedImage;
123  }
124  else if(ImageType::ImageDimension == 3)
125  {
126  return reorient3DImage<ImageType>(input, orientation);
127  }
128  else
129  {
130  unsigned int dimension = ImageType::ImageDimension;
131  std::string msg = "Reorient an image of dimension "
132  + std::to_string(dimension)
133  + "is not supported yet.";
134  throw itk::ExceptionObject(__FILE__, __LINE__, msg , ITK_LOCATION);
135  }
136 }
137 
138 template <class ImageType, class GradientType>
139 void reorientGradients(typename itk::SmartPointer<ImageType> input, std::vector<GradientType> &gradients)
140 {
141  // What we want to do here is reorient a list of gradient
142  // vectors in order to have them in image coordinates from a set of gradients in real coordinates
143 
144  // To do so we apply 'reorientedGrad = D^-1 o grad'
145 
146  typename ImageType::DirectionType inversedD1 = input->GetInverseDirection();
147 
148  for(unsigned int g = 0; g < gradients.size(); g++)
149  {
150  GradientType reorientedGrad(0.0);
151  for(unsigned int i = 0; i < 3; ++i)
152  {
153  for(unsigned int j = 0; j < 3; ++j)
154  reorientedGrad[i] += inversedD1[i][j] * gradients[g][j];
155  }
156  gradients[g] = reorientedGrad;
157  }
158 }
159 
160 } // end namespace anima
161 
itk::SmartPointer< ImageType > reorient3DImage(typename itk::SmartPointer< ImageType > input, itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
itk::SmartPointer< ImageType > reorientImage(typename itk::SmartPointer< ImageType > input, itk::SpatialOrientation::ValidCoordinateOrientationFlags orientation)
void reorientGradients(typename itk::SmartPointer< ImageType > input, std::vector< GradientType > &gradients)