ANIMA  4.0
animaSymmetricBMRegistrationMethod.hxx
Go to the documentation of this file.
1 #pragma once
3 
6 
7 namespace anima
8 {
9 
10 template <typename TInputImageType>
11 void
12 SymmetricBMRegistrationMethod <TInputImageType>
13 ::PerformOneIteration(InputImageType *refImage, InputImageType *movingImage, TransformPointer &addOn)
14 {
15  itk::TimeProbe tmpTime;
16  tmpTime.Start();
17 
18  this->GetBlockMatcher()->SetForceComputeBlocks(false);
19  this->GetBlockMatcher()->SetReferenceImage(this->GetFixedImage());
20  this->GetBlockMatcher()->SetMovingImage(movingImage);
21  this->GetBlockMatcher()->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
22  this->GetBlockMatcher()->Update();
23 
24  tmpTime.Stop();
25 
26  if (this->GetVerboseProgression())
27  std::cout << "Matching performed in " << tmpTime.GetTotal() << std::endl;
28 
29  this->GetAgregator()->SetInputRegions(this->GetBlockMatcher()->GetBlockRegions());
30  this->GetAgregator()->SetInputOrigins(this->GetBlockMatcher()->GetBlockPositions());
31  this->GetAgregator()->SetInputWeights(this->GetBlockMatcher()->GetBlockWeights());
32  this->GetAgregator()->SetInputTransforms(this->GetBlockMatcher()->GetBlockTransformPointers());
33 
34  TransformPointer usualAddOn = this->GetAgregator()->GetOutput();
35 
36  itk::TimeProbe tmpTimeReverse;
37  tmpTimeReverse.Start();
38 
39  m_ReverseBlockMatcher->SetForceComputeBlocks(false);
40  m_ReverseBlockMatcher->SetReferenceImage(this->GetMovingImage());
41  m_ReverseBlockMatcher->SetMovingImage(refImage);
42  m_ReverseBlockMatcher->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
43  m_ReverseBlockMatcher->Update();
44 
45  tmpTimeReverse.Stop();
46 
47  if (this->GetVerboseProgression())
48  std::cout << "Matching performed in " << tmpTimeReverse.GetTotal() << std::endl;
49 
50  this->GetAgregator()->SetInputRegions(m_ReverseBlockMatcher->GetBlockRegions());
51  this->GetAgregator()->SetInputOrigins(m_ReverseBlockMatcher->GetBlockPositions());
52  this->GetAgregator()->SetInputWeights(m_ReverseBlockMatcher->GetBlockWeights());
53  this->GetAgregator()->SetInputTransforms(m_ReverseBlockMatcher->GetBlockTransformPointers());
54 
55  TransformPointer reverseAddOn = this->GetAgregator()->GetOutput();
56 
57  if (this->GetAgregator()->GetOutputTransformType() == AgregatorType::SVF)
58  {
59  // Add update to current velocity field (cf. Vercauteren et al, 2008)
60  // First compute the SVF from two asymmetric ones: S = 0.5 * (S_0 - S_1)
61  // It's 0.5 since we are computing the full transform between the two images
62  typedef typename SVFTransformType::VectorFieldType VelocityFieldType;
63 
64  typedef typename itk::ImageRegionConstIterator <VelocityFieldType> VelocityFieldConstIterator;
65  typedef typename itk::ImageRegionIterator <VelocityFieldType> VelocityFieldIterator;
66 
67  SVFTransformType *usualAddOnCast = dynamic_cast <SVFTransformType *> (usualAddOn.GetPointer());
68  SVFTransformType *reverseAddOnCast = dynamic_cast <SVFTransformType *> (reverseAddOn.GetPointer());
69 
70  typename VelocityFieldType::Pointer usualAddOnSVF = const_cast <VelocityFieldType *> (usualAddOnCast->GetParametersAsVectorField());
71 
72  VelocityFieldIterator usualAddOnItr(usualAddOnSVF,usualAddOnSVF->GetLargestPossibleRegion());
73 
74  VelocityFieldConstIterator reverseAddOnItr(reverseAddOnCast->GetParametersAsVectorField(),
75  reverseAddOnCast->GetParametersAsVectorField()->GetLargestPossibleRegion());
76 
77  typedef typename VelocityFieldType::PixelType VectorType;
78  VectorType tmpVec;
79  while (!usualAddOnItr.IsAtEnd())
80  {
81  tmpVec = 0.5 * (usualAddOnItr.Get() - reverseAddOnItr.Get());
82  usualAddOnItr.Set(tmpVec);
83 
84  ++usualAddOnItr;
85  ++reverseAddOnItr;
86  }
87  }
88  else
89  {
90  AffineTransformType *usualAddOnCast = dynamic_cast <AffineTransformType *> (usualAddOn.GetPointer());
91  AffineTransformType *reverseAddOnCast = dynamic_cast <AffineTransformType *> (reverseAddOn.GetPointer());
92 
93  unsigned int NDimensions = InputImageType::ImageDimension;
94  vnl_matrix <double> usualAddOnMatrix(NDimensions+1,NDimensions+1,0);
95  vnl_matrix <double> reverseAddOnMatrix(NDimensions+1,NDimensions+1,0);
96  usualAddOnMatrix.set_identity();
97  reverseAddOnMatrix.set_identity();
98 
99  for (unsigned int i = 0;i < NDimensions;++i)
100  {
101  for (unsigned int j = 0;j < NDimensions;++j)
102  {
103  usualAddOnMatrix(i,j) = usualAddOnCast->GetMatrix()(i,j);
104  reverseAddOnMatrix(i,j) = reverseAddOnCast->GetMatrix()(i,j);
105  }
106 
107  usualAddOnMatrix(i,NDimensions) = usualAddOnCast->GetOffset()[i];
108  reverseAddOnMatrix(i,NDimensions) = reverseAddOnCast->GetOffset()[i];
109  }
110 
111  usualAddOnMatrix = anima::GetLogarithm(usualAddOnMatrix);
112  reverseAddOnMatrix = anima::GetLogarithm(reverseAddOnMatrix);
113 
114  usualAddOnMatrix -= reverseAddOnMatrix;
115  usualAddOnMatrix /= 2.0;
116 
117  usualAddOnMatrix = anima::GetExponential(usualAddOnMatrix);
118 
119  typename AffineTransformType::MatrixType trsfMatrix;
120 
121  for (unsigned int i = 0;i < NDimensions;++i)
122  for (unsigned int j = 0;j < NDimensions;++j)
123  trsfMatrix(i,j) = usualAddOnMatrix(i,j);
124 
125  usualAddOnCast->SetMatrix(trsfMatrix);
126 
127  typename AffineTransformType::OffsetType trsfOffset;
128  for (unsigned int i = 0;i < NDimensions;++i)
129  trsfOffset[i] = usualAddOnMatrix(i,NDimensions);
130 
131  usualAddOnCast->SetOffset(trsfOffset);
132  }
133 
134  addOn = usualAddOn;
135 }
136 
137 }
itk::AffineTransform< typename AgregatorType::ScalarType, TInputImageType::ImageDimension > AffineTransformType
vnl_matrix< T > GetExponential(const vnl_matrix< T > &m, const int numApprox=1)
Computation of the matrix exponential. Algo: classical scaling and squaring, as in Matlab...
vnl_matrix< T > GetLogarithm(const vnl_matrix< T > &m, const double precision=0.00000000001, const int numApprox=1)
Computation of the matrix logarithm. Algo: inverse scaling and squaring, variant proposed by Cheng et...
itk::StationaryVelocityFieldTransform< AgregatorScalarType, TInputImageType::ImageDimension > SVFTransformType