ANIMA  4.0
animaSVFExponentialImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
5 #include <itkImageRegionConstIterator.h>
6 #include <itkImageRegionIterator.h>
7 
8 #include <itkComposeDisplacementFieldsImageFilter.h>
9 #include <itkVectorLinearInterpolateNearestNeighborExtrapolateImageFunction.h>
10 
11 namespace anima
12 {
13 
14 template <typename TPixelType, unsigned int Dimension>
15 void
16 SVFExponentialImageFilter <TPixelType, Dimension>
17 ::BeforeThreadedGenerateData()
18 {
19  this->Superclass::BeforeThreadedGenerateData();
20 
21  if ((m_ExponentiationOrder != 0)&&(m_ExponentiationOrder != 1))
22  itkExceptionMacro("Exponentiation order not supported");
23 
24  // Precomputes jacobian if needed
25  if (m_ExponentiationOrder > 0)
26  {
28 
29  typename JacobianFilterType::Pointer jacFilter = JacobianFilterType::New();
30  jacFilter->SetInput(this->GetInput());
31  jacFilter->SetNoIdentity(true);
32  jacFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
33  jacFilter->SetNeighborhood(0);
34 
35  jacFilter->Update();
36 
37  m_FieldJacobian = jacFilter->GetOutput();
38  m_FieldJacobian->DisconnectPipeline();
39  }
40 
41  // Computes field maximal norm
42  typedef itk::ImageRegionConstIterator <InputImageType> IteratorType;
43  IteratorType inItr(this->GetInput(),this->GetInput()->GetLargestPossibleRegion());
44 
45  double maxNorm = 0;
46  InputPixelType inputValue;
47  while (!inItr.IsAtEnd())
48  {
49  double norm = 0;
50  inputValue = inItr.Get();
51 
52  for (unsigned int i = 0;i < Dimension;++i)
53  norm += inputValue[i] * inputValue[i];
54 
55  if (norm > maxNorm)
56  maxNorm = norm;
57 
58  ++inItr;
59  }
60 
61  // Taken from Vercauteren et al. smart initialization of number of squarings necessary
62  double pixelSpacing = this->GetInput()->GetSpacing()[0];
63  for (unsigned int i = 1;i < Dimension;++i)
64  {
65  if (this->GetInput()->GetSpacing()[i] < pixelSpacing)
66  pixelSpacing = this->GetInput()->GetSpacing()[i];
67  }
68 
69  maxNorm = std::sqrt(maxNorm);
70 
71  double numIter = std::log(maxNorm / (m_MaximalDisplacementAmplitude * pixelSpacing)) / std::log(2.0);
72 
73  m_NumberOfSquarings = 0;
74  if (numIter + 1 > 0)
75  m_NumberOfSquarings = static_cast<unsigned int>(numIter + 1.0);
76 }
77 
78 template <typename TPixelType, unsigned int Dimension>
79 void
81 ::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread)
82 {
83  typedef itk::ImageRegionConstIterator <InputImageType> InputIteratorType;
84  typedef itk::ImageRegionConstIterator <JacobianImageType> JacobianIteratorType;
85  typedef itk::ImageRegionIterator <OutputImageType> OutIteratorType;
86 
87  InputIteratorType inputItr(this->GetInput(),outputRegionForThread);
88  OutIteratorType outItr(this->GetOutput(),outputRegionForThread);
89 
90  JacobianIteratorType jacItr;
91  if (m_ExponentiationOrder > 0)
92  jacItr = JacobianIteratorType(m_FieldJacobian,outputRegionForThread);
93 
94  InputPixelType inputValue;
95  OutputPixelType outputValue;
96  JacobianPixelType jacValue;
97 
98  double scalingFactor = 1.0 / std::pow(2.0, m_NumberOfSquarings);
99  while (!outItr.IsAtEnd())
100  {
101  inputValue = inputItr.Get();
102  for (unsigned int i = 0;i < Dimension;++i)
103  outputValue[i] = scalingFactor * inputValue[i];
104 
105  if (m_ExponentiationOrder > 0)
106  {
107  jacValue = jacItr.Get();
108  for (unsigned int i = 0;i < Dimension;++i)
109  {
110  for (unsigned int j = 0;j < Dimension;++j)
111  outputValue[i] += 0.5 * scalingFactor * scalingFactor * jacValue[i * Dimension + j] * inputValue[j];
112  }
113  }
114 
115  outItr.Set(outputValue);
116 
117  ++inputItr;
118  ++outItr;
119  if (m_ExponentiationOrder > 0)
120  ++jacItr;
121  }
122 }
123 
124 template <typename TPixelType, unsigned int Dimension>
125 void
127 ::AfterThreadedGenerateData()
128 {
129  this->Superclass::AfterThreadedGenerateData();
130 
131  // Compute recursive squaring of the output
132  typedef itk::ComposeDisplacementFieldsImageFilter <OutputImageType,OutputImageType> ComposeFilterType;
133  typedef itk::VectorLinearInterpolateNearestNeighborExtrapolateImageFunction <OutputImageType,
134  typename OutputImageType::PixelType::ValueType> VectorInterpolateFunctionType;
135 
136  typename OutputImageType::Pointer outputPtr = this->GetOutput();
137 
138  for (unsigned int i = 0;i < m_NumberOfSquarings;++i)
139  {
140  typename ComposeFilterType::Pointer composer = ComposeFilterType::New();
141  composer->SetWarpingField(outputPtr);
142  composer->SetDisplacementField(outputPtr);
143  composer->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
144 
145  typename VectorInterpolateFunctionType::Pointer interpolator = VectorInterpolateFunctionType::New();
146 
147  composer->SetInterpolator(interpolator);
148  composer->Update();
149  outputPtr = composer->GetOutput();
150  outputPtr->DisconnectPipeline();
151  }
152 
153  if (m_NumberOfSquarings > 0)
154  this->GraftOutput(outputPtr);
155 }
156 
157 } // end namespace anima
itk::Image< itk::Vector< TPixelType, Dimension >, Dimension > OutputImageType
Compute the Jacobian matrix in real coordinates of a displacement field.
Computes the exponentiation of a stationary velocity field using sclaing and squaring and approximate...
JacobianImageType::PixelType JacobianPixelType
Superclass::OutputImageRegionType OutputImageRegionType