ANIMA  4.0
animaSVFLieBracketImageFilter.hxx
Go to the documentation of this file.
1 #pragma once
3 
5 #include <itkImageRegionConstIterator.h>
6 #include <itkImageRegionIterator.h>
7 
8 namespace anima
9 {
10 
11 template <typename TPixelType, unsigned int Dimension>
12 void
13 SVFLieBracketImageFilter <TPixelType, Dimension>
14 ::BeforeThreadedGenerateData()
15 {
16  this->Superclass::BeforeThreadedGenerateData();
17 
19 
20  if (m_FirstFieldJacobian.IsNull())
21  {
22  typename JacobianFilterType::Pointer jacFilter = JacobianFilterType::New();
23  jacFilter->SetInput(this->GetInput(0));
24  jacFilter->SetNoIdentity(true);
25  jacFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
26  jacFilter->SetNeighborhood(0);
27 
28  jacFilter->Update();
29 
30  m_FirstFieldJacobian = jacFilter->GetOutput();
31  m_FirstFieldJacobian->DisconnectPipeline();
32  }
33 
34  if (m_SecondFieldJacobian.IsNull())
35  {
36  typename JacobianFilterType::Pointer jacFilter = JacobianFilterType::New();
37  jacFilter->SetInput(this->GetInput(1));
38  jacFilter->SetNoIdentity(true);
39  jacFilter->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
40  jacFilter->SetNeighborhood(0);
41 
42  jacFilter->Update();
43 
44  m_SecondFieldJacobian = jacFilter->GetOutput();
45  m_SecondFieldJacobian->DisconnectPipeline();
46  }
47 }
48 
49 template <typename TPixelType, unsigned int Dimension>
50 void
52 ::DynamicThreadedGenerateData(const OutputImageRegionType &outputRegionForThread)
53 {
54  typedef itk::ImageRegionConstIterator <InputImageType> InputIteratorType;
55  typedef itk::ImageRegionConstIterator <JacobianImageType> JacobianIteratorType;
56  typedef itk::ImageRegionIterator <OutputImageType> OutIteratorType;
57 
58  InputIteratorType firstInputItr(this->GetInput(0),outputRegionForThread);
59  InputIteratorType secondInputItr(this->GetInput(1),outputRegionForThread);
60  OutIteratorType outItr(this->GetOutput(),outputRegionForThread);
61 
62  JacobianIteratorType firstJacItr(m_FirstFieldJacobian,outputRegionForThread);
63  JacobianIteratorType secondJacItr(m_SecondFieldJacobian,outputRegionForThread);
64 
65  InputPixelType inputValue;
66  OutputPixelType outputValue;
67  JacobianPixelType jacValue;
68 
69  while (!outItr.IsAtEnd())
70  {
71  outputValue.Fill(0);
72 
73  inputValue = secondInputItr.Get();
74  jacValue = firstJacItr.Get();
75 
76  for (unsigned int i = 0;i < Dimension;++i)
77  {
78  for (unsigned int j = 0;j < Dimension;++j)
79  outputValue[i] += jacValue[i * Dimension + j] * inputValue[j];
80  }
81 
82  inputValue = firstInputItr.Get();
83  jacValue = secondJacItr.Get();
84 
85  for (unsigned int i = 0;i < Dimension;++i)
86  {
87  for (unsigned int j = 0;j < Dimension;++j)
88  outputValue[i] -= jacValue[i * Dimension + j] * inputValue[j];
89  }
90 
91  outItr.Set(outputValue);
92 
93  ++firstInputItr;
94  ++secondInputItr;
95  ++outItr;
96  ++firstJacItr;
97  ++secondJacItr;
98  }
99 }
100 
101 } // end namespace anima
Compute the Jacobian matrix in real coordinates of a displacement field.
Computes the Lie bracket between two fields u and v as expressed by Bossa et al.
Superclass::OutputImageRegionType OutputImageRegionType
JacobianImageType::PixelType JacobianPixelType