ANIMA  4.0
animaOrientedModelBaseResampleImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
4 #include <itkImageRegionConstIteratorWithIndex.h>
5 #include <itkImageRegionIteratorWithIndex.h>
7 
8 #include <itkTranslationTransform.h>
9 #include <itkMatrixOffsetTransformBase.h>
10 #include <itkCompositeTransform.h>
11 
12 #include <animaBaseTensorTools.h>
13 
14 namespace anima
15 {
16 
17 template <typename TImageType, typename TInterpolatorPrecisionType>
18 void
21 {
23 
24  typename InterpolatorType::Pointer tmpInterpolator = InterpolatorType::New();
25  this->SetInterpolator(tmpInterpolator.GetPointer());
26 }
27 
28 
29 template <typename TImageType, typename TInterpolatorPrecisionType>
30 void
33 {
34  Superclass::GenerateInputRequestedRegion();
35  if (!this->GetInput())
36  return;
37 
38  InputImagePointer inputPtr = const_cast <TImageType *> (this->GetInput());
39  inputPtr->SetRequestedRegionToLargestPossibleRegion();
40 }
41 
42 template <typename TImageType, typename TInterpolatorPrecisionType>
43 void
46 {
47  Superclass::GenerateOutputInformation();
48 
49  this->GetOutput()->SetSpacing(m_OutputSpacing);
50  this->GetOutput()->SetOrigin(m_OutputOrigin);
51  this->GetOutput()->SetDirection(m_OutputDirection);
52  this->GetOutput()->SetRegions(m_OutputLargestPossibleRegion);
53  this->GetOutput()->SetNumberOfComponentsPerPixel(this->GetOutputVectorLength());
54 }
55 
56 template <typename TImageType, typename TInterpolatorPrecisionType>
57 unsigned int
60 {
61  if (this->GetNumberOfIndexedInputs() > 0)
62  return this->GetInput(0)->GetNumberOfComponentsPerPixel();
63  else
64  return 0;
65 }
66 
67 template <typename TImageType, typename TInterpolatorPrecisionType>
68 void
71 {
72  Superclass::BeforeThreadedGenerateData();
73 
74  if (m_Transform.IsNull())
75  itkExceptionMacro("No valid transformation...");
76 
77  if (m_Interpolator.IsNull())
78  this->InitializeInterpolator();
79 
80  m_Interpolator->SetInputImage(this->GetInput(0));
81 
82  if (!m_Transform->IsLinear())
83  {
84  m_StartIndDef = m_OutputLargestPossibleRegion.GetIndex();
85  m_EndIndDef = m_StartIndDef + m_OutputLargestPossibleRegion.GetSize();
86  }
87 }
88 
89 template <typename TImageType, typename TInterpolatorPrecisionType>
90 void
93 {
94  unsigned int threadId = this->GetSafeThreadId();
95 
96  if (m_Transform->IsLinear())
97  this->LinearThreadedGenerateData(outputRegionForThread,threadId);
98  else
99  this->NonLinearThreadedGenerateData(outputRegionForThread,threadId);
100 
101  this->SafeReleaseThreadId(threadId);
102 }
103 
104 template <typename TImageType, typename TInterpolatorPrecisionType>
105 void
107 ::LinearThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, unsigned int threadId)
108 {
109  typedef itk::ImageRegionIteratorWithIndex <InputImageType> IteratorType;
110 
111  IteratorType outputItr(this->GetOutput(),outputRegionForThread);
112 
113  InputIndexType tmpInd;
114  PointType tmpPoint;
115  unsigned int vectorSize = this->GetOutputVectorLength();
116  InputPixelType tmpRes(vectorSize), resRotated(vectorSize);
117  ContinuousIndexType index;
118 
119  vnl_matrix <double> orientationMatrix = this->ComputeLinearJacobianMatrix();
120  vnl_matrix <double> parametersRotationMatrix;
121  this->ComputeRotationParametersFromReorientationMatrix(orientationMatrix,parametersRotationMatrix);
122 
123  bool lastDimensionUseless = (this->GetInput(0)->GetLargestPossibleRegion().GetSize()[ImageDimension - 1] <= 1);
124 
125  while (!outputItr.IsAtEnd())
126  {
127  tmpInd = outputItr.GetIndex();
128 
129  this->GetOutput()->TransformIndexToPhysicalPoint(tmpInd,tmpPoint);
130 
131  tmpPoint = m_Transform->TransformPoint(tmpPoint);
132 
133  this->GetInput(0)->TransformPhysicalPointToContinuousIndex(tmpPoint,index);
134 
135  if (lastDimensionUseless)
136  index[ImageDimension - 1] = 0;
137 
138  if (m_Interpolator->IsInsideBuffer(index))
139  tmpRes = m_Interpolator->EvaluateAtContinuousIndex(index);
140  else
141  this->InitializeZeroPixel(tmpRes);
142 
143  if (!isZero(tmpRes))
144  {
145  this->ReorientInterpolatedModel(tmpRes,parametersRotationMatrix,resRotated,threadId);
146  outputItr.Set(resRotated);
147  }
148  else
149  outputItr.Set(tmpRes);
150 
151  ++outputItr;
152  }
153 }
154 
155 
156 template <typename TImageType, typename TInterpolatorPrecisionType>
157 void
159 ::NonLinearThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, unsigned int threadId)
160 {
161  typedef itk::ImageRegionIteratorWithIndex <InputImageType> IteratorType;
162 
163  IteratorType outputItr(this->GetOutput(),outputRegionForThread);
164 
165  InputIndexType tmpInd;
166  PointType tmpPoint;
167  unsigned int vectorSize = this->GetOutputVectorLength();
168  InputPixelType tmpRes(vectorSize), resRotated(vectorSize);
169  ContinuousIndexType index;
170 
171  vnl_matrix <double> orientationMatrix(ImageDimension,ImageDimension);
172  vnl_matrix <double> parametersRotationMatrix;
173 
174  while (!outputItr.IsAtEnd())
175  {
176  tmpInd = outputItr.GetIndex();
177  this->GetOutput()->TransformIndexToPhysicalPoint(tmpInd,tmpPoint);
178 
179  tmpPoint = m_Transform->TransformPoint(tmpPoint);
180 
181  this->GetInput(0)->TransformPhysicalPointToContinuousIndex(tmpPoint,index);
182 
183  if (m_Interpolator->IsInsideBuffer(index))
184  tmpRes = m_Interpolator->EvaluateAtContinuousIndex(index);
185  else
186  this->InitializeZeroPixel(tmpRes);
187 
188  if (!isZero(tmpRes))
189  {
190  this->ComputeLocalJacobianMatrix(tmpInd,orientationMatrix);
191  this->ComputeRotationParametersFromReorientationMatrix(orientationMatrix,parametersRotationMatrix);
192  this->ReorientInterpolatedModel(tmpRes,parametersRotationMatrix,resRotated,threadId);
193  outputItr.Set(resRotated);
194  }
195  else
196  outputItr.Set(tmpRes);
197 
198  ++outputItr;
199  }
200 }
201 
202 template <typename TImageType, typename TInterpolatorPrecisionType>
203 vnl_matrix <double>
206 {
207  vnl_matrix <double> reorientationMatrix(ImageDimension,ImageDimension);
208  reorientationMatrix.set_identity();
209 
210  typedef itk::TranslationTransform <TInterpolatorPrecisionType,ImageDimension> TranslationType;
211  if (dynamic_cast <TranslationType *> (m_Transform.GetPointer()))
212  return reorientationMatrix;
213 
214  typedef itk::MatrixOffsetTransformBase <TInterpolatorPrecisionType,ImageDimension,ImageDimension> AffineTransformType;
215 
216  AffineTransformType *tmpTrsf = dynamic_cast <AffineTransformType *> (m_Transform.GetPointer());
217 
218  if (tmpTrsf)
219  {
220  reorientationMatrix = tmpTrsf->GetMatrix().GetVnlMatrix().as_matrix();
221  return reorientationMatrix;
222  }
223 
224  typedef itk::CompositeTransform <TInterpolatorPrecisionType,ImageDimension> CompositeTransformType;
225  CompositeTransformType *compositeTrsf = dynamic_cast <CompositeTransformType *> (m_Transform.GetPointer());
226 
227  for (unsigned int i = 0;i < compositeTrsf->GetNumberOfTransforms();++i)
228  {
229  tmpTrsf = dynamic_cast <AffineTransformType *> (compositeTrsf->GetNthTransform(i).GetPointer());
230  if (tmpTrsf)
231  reorientationMatrix *= tmpTrsf->GetMatrix().GetVnlMatrix().as_matrix();
232  }
233 
234  return reorientationMatrix;
235 }
236 
237 template <typename TImageType, typename TInterpolatorPrecisionType>
238 void
240 ::ComputeLocalJacobianMatrix(InputIndexType &index, vnl_matrix <double> &reorientationMatrix)
241 {
242  vnl_matrix <double> jacMatrix(ImageDimension,ImageDimension);
243 
244  vnl_matrix <double> deltaMatrix(ImageDimension,ImageDimension,0);
245  InputIndexType posBef, posAfter;
246  PointType tmpPosBef, tmpPosAfter;
247 
248  vnl_matrix <double> resDiff(ImageDimension,ImageDimension,0);
249  OutputImagePointer outputPtr = this->GetOutput();
250 
251  PointType tmpPos;
252 
253  for (unsigned int i = 0;i < ImageDimension;++i)
254  {
255  posBef = index;
256  posBef[i]--;
257 
258  if (posBef[i] < m_StartIndDef[i])
259  posBef[i] = m_StartIndDef[i];
260 
261  outputPtr->TransformIndexToPhysicalPoint(posBef,tmpPosBef);
262 
263  posAfter = index;
264  posAfter[i]++;
265 
266  if (posAfter[i] >= m_EndIndDef[i])
267  posAfter[i] = m_EndIndDef[i] - 1;
268 
269  outputPtr->TransformIndexToPhysicalPoint(posAfter,tmpPosAfter);
270 
271  if (posAfter[i] == posBef[i])
272  {
273  deltaMatrix(i,i) = 1;
274  continue;
275  }
276 
277  double normDelta = 0.0;
278  for (unsigned int j = 0;j < ImageDimension;++j)
279  {
280  deltaMatrix(j,i) = tmpPosAfter[j] - tmpPosBef[j];
281  normDelta += deltaMatrix(j,i) * deltaMatrix(j,i);
282  }
283 
284  normDelta = std::sqrt(normDelta);
285  for (unsigned int j = 0;j < ImageDimension;++j)
286  deltaMatrix(j,i) /= normDelta;
287 
288  tmpPos = m_Transform->TransformPoint(tmpPosAfter);
289  for (unsigned int j = 0;j < ImageDimension;++j)
290  resDiff(j,i) = tmpPos[j];
291 
292  tmpPos = m_Transform->TransformPoint(tmpPosBef);
293  for (unsigned int j = 0;j < ImageDimension;++j)
294  resDiff(j,i) = (resDiff(j,i) - tmpPos[j]) / normDelta;
295  }
296 
297  vnl_matrix_inverse <double> invDelta(deltaMatrix);
298  deltaMatrix = invDelta.inverse();
299 
300  for (unsigned int i = 0;i < ImageDimension;++i)
301  {
302  for (unsigned int j = 0;j < ImageDimension;++j)
303  {
304  jacMatrix(i,j) = 0;
305  for (unsigned int k = 0;k < ImageDimension;++k)
306  jacMatrix(i,j) += resDiff(i,k) * deltaMatrix(k,j);
307  }
308  }
309 
310  if (m_StartIndDef[2] == (m_EndIndDef[2]-1))
311  jacMatrix(2,2) = 1;
312 
313  reorientationMatrix = jacMatrix;
314 }
315 
316 template <typename TImageType, typename TInterpolatorPrecisionType>
317 void
319 ::ComputeRotationParametersFromReorientationMatrix(vnl_matrix <double> &reorientationMatrix,
320  vnl_matrix <double> &modelOrientationMatrix)
321 {
322  if (m_FiniteStrainReorientation)
323  {
324  vnl_matrix <double> tmpMat(ImageDimension,ImageDimension);
325  anima::ExtractRotationFromJacobianMatrix(reorientationMatrix,modelOrientationMatrix,tmpMat);
326  }
327  else
328  modelOrientationMatrix = reorientationMatrix;
329 }
330 
331 } // end of namespace anima
void DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread) ITK_OVERRIDE
void ComputeLocalJacobianMatrix(InputIndexType &index, vnl_matrix< double > &reorientationMatrix)
virtual void InitializeInterpolator()
Initializes the default interpolator, might change in derived classes.
void LinearThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, unsigned int threadId)
itk::InterpolateImageFunction< InputImageType, TInterpolatorPrecisionType > InterpolatorType
void ExtractRotationFromJacobianMatrix(vnl_matrix< RealType > &jacobianMatrix, vnl_matrix< RealType > &rotationMatrix, vnl_matrix< RealType > &tmpMat)
void NonLinearThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, unsigned int threadId)
virtual void ComputeRotationParametersFromReorientationMatrix(vnl_matrix< double > &reorientationMatrix, vnl_matrix< double > &modelOrientationMatrix)