10 template <
typename TInputImageType>
12 SymmetricBMRegistrationMethod <TInputImageType>
15 itk::TimeProbe tmpTime;
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();
26 if (this->GetVerboseProgression())
27 std::cout <<
"Matching performed in " << tmpTime.GetTotal() << std::endl;
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());
36 itk::TimeProbe tmpTimeReverse;
37 tmpTimeReverse.Start();
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();
45 tmpTimeReverse.Stop();
47 if (this->GetVerboseProgression())
48 std::cout <<
"Matching performed in " << tmpTimeReverse.GetTotal() << std::endl;
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());
57 if (this->GetAgregator()->GetOutputTransformType() == AgregatorType::SVF)
62 typedef typename SVFTransformType::VectorFieldType VelocityFieldType;
64 typedef typename itk::ImageRegionConstIterator <VelocityFieldType> VelocityFieldConstIterator;
65 typedef typename itk::ImageRegionIterator <VelocityFieldType> VelocityFieldIterator;
70 typename VelocityFieldType::Pointer usualAddOnSVF = const_cast <VelocityFieldType *> (usualAddOnCast->GetParametersAsVectorField());
72 VelocityFieldIterator usualAddOnItr(usualAddOnSVF,usualAddOnSVF->GetLargestPossibleRegion());
74 VelocityFieldConstIterator reverseAddOnItr(reverseAddOnCast->GetParametersAsVectorField(),
75 reverseAddOnCast->GetParametersAsVectorField()->GetLargestPossibleRegion());
77 typedef typename VelocityFieldType::PixelType VectorType;
79 while (!usualAddOnItr.IsAtEnd())
81 tmpVec = 0.5 * (usualAddOnItr.Get() - reverseAddOnItr.Get());
82 usualAddOnItr.Set(tmpVec);
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();
99 for (
unsigned int i = 0;i < NDimensions;++i)
101 for (
unsigned int j = 0;j < NDimensions;++j)
103 usualAddOnMatrix(i,j) = usualAddOnCast->GetMatrix()(i,j);
104 reverseAddOnMatrix(i,j) = reverseAddOnCast->GetMatrix()(i,j);
107 usualAddOnMatrix(i,NDimensions) = usualAddOnCast->GetOffset()[i];
108 reverseAddOnMatrix(i,NDimensions) = reverseAddOnCast->GetOffset()[i];
114 usualAddOnMatrix -= reverseAddOnMatrix;
115 usualAddOnMatrix /= 2.0;
119 typename AffineTransformType::MatrixType trsfMatrix;
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);
125 usualAddOnCast->SetMatrix(trsfMatrix);
127 typename AffineTransformType::OffsetType trsfOffset;
128 for (
unsigned int i = 0;i < NDimensions;++i)
129 trsfOffset[i] = usualAddOnMatrix(i,NDimensions);
131 usualAddOnCast->SetOffset(trsfOffset);
TInputImageType InputImageType
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
TransformType::Pointer TransformPointer