5 #include <itkImageRegionConstIterator.h> 6 #include <itkImageRegionIterator.h> 8 #include <itkComposeDisplacementFieldsImageFilter.h> 9 #include <itkVectorLinearInterpolateNearestNeighborExtrapolateImageFunction.h> 14 template <
typename TPixelType,
unsigned int Dimension>
16 SVFExponentialImageFilter <TPixelType, Dimension>
17 ::BeforeThreadedGenerateData()
19 this->Superclass::BeforeThreadedGenerateData();
21 if ((m_ExponentiationOrder != 0)&&(m_ExponentiationOrder != 1))
22 itkExceptionMacro(
"Exponentiation order not supported");
25 if (m_ExponentiationOrder > 0)
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);
37 m_FieldJacobian = jacFilter->GetOutput();
38 m_FieldJacobian->DisconnectPipeline();
42 typedef itk::ImageRegionConstIterator <InputImageType> IteratorType;
43 IteratorType inItr(this->GetInput(),this->GetInput()->GetLargestPossibleRegion());
47 while (!inItr.IsAtEnd())
50 inputValue = inItr.Get();
52 for (
unsigned int i = 0;i < Dimension;++i)
53 norm += inputValue[i] * inputValue[i];
62 double pixelSpacing = this->GetInput()->GetSpacing()[0];
63 for (
unsigned int i = 1;i < Dimension;++i)
65 if (this->GetInput()->GetSpacing()[i] < pixelSpacing)
66 pixelSpacing = this->GetInput()->GetSpacing()[i];
69 maxNorm = std::sqrt(maxNorm);
71 double numIter = std::log(maxNorm / (m_MaximalDisplacementAmplitude * pixelSpacing)) / std::log(2.0);
73 m_NumberOfSquarings = 0;
75 m_NumberOfSquarings =
static_cast<unsigned int>(numIter + 1.0);
78 template <
typename TPixelType,
unsigned int Dimension>
83 typedef itk::ImageRegionConstIterator <InputImageType> InputIteratorType;
84 typedef itk::ImageRegionConstIterator <JacobianImageType> JacobianIteratorType;
85 typedef itk::ImageRegionIterator <OutputImageType> OutIteratorType;
87 InputIteratorType inputItr(this->GetInput(),outputRegionForThread);
88 OutIteratorType outItr(this->GetOutput(),outputRegionForThread);
90 JacobianIteratorType jacItr;
91 if (m_ExponentiationOrder > 0)
92 jacItr = JacobianIteratorType(m_FieldJacobian,outputRegionForThread);
98 double scalingFactor = 1.0 / std::pow(2.0, m_NumberOfSquarings);
99 while (!outItr.IsAtEnd())
101 inputValue = inputItr.Get();
102 for (
unsigned int i = 0;i < Dimension;++i)
103 outputValue[i] = scalingFactor * inputValue[i];
105 if (m_ExponentiationOrder > 0)
107 jacValue = jacItr.Get();
108 for (
unsigned int i = 0;i < Dimension;++i)
110 for (
unsigned int j = 0;j < Dimension;++j)
111 outputValue[i] += 0.5 * scalingFactor * scalingFactor * jacValue[i * Dimension + j] * inputValue[j];
115 outItr.Set(outputValue);
119 if (m_ExponentiationOrder > 0)
124 template <
typename TPixelType,
unsigned int Dimension>
127 ::AfterThreadedGenerateData()
129 this->Superclass::AfterThreadedGenerateData();
132 typedef itk::ComposeDisplacementFieldsImageFilter <OutputImageType,OutputImageType> ComposeFilterType;
133 typedef itk::VectorLinearInterpolateNearestNeighborExtrapolateImageFunction <
OutputImageType,
134 typename OutputImageType::PixelType::ValueType> VectorInterpolateFunctionType;
136 typename OutputImageType::Pointer outputPtr = this->GetOutput();
138 for (
unsigned int i = 0;i < m_NumberOfSquarings;++i)
140 typename ComposeFilterType::Pointer composer = ComposeFilterType::New();
141 composer->SetWarpingField(outputPtr);
142 composer->SetDisplacementField(outputPtr);
143 composer->SetNumberOfWorkUnits(this->GetNumberOfWorkUnits());
145 typename VectorInterpolateFunctionType::Pointer interpolator = VectorInterpolateFunctionType::New();
147 composer->SetInterpolator(interpolator);
149 outputPtr = composer->GetOutput();
150 outputPtr->DisconnectPipeline();
153 if (m_NumberOfSquarings > 0)
154 this->GraftOutput(outputPtr);
itk::Image< itk::Vector< TPixelType, Dimension >, Dimension > OutputImageType
InputImageType::PixelType InputPixelType
Compute the Jacobian matrix in real coordinates of a displacement field.
OutputImageType::PixelType OutputPixelType
Computes the exponentiation of a stationary velocity field using sclaing and squaring and approximate...
JacobianImageType::PixelType JacobianPixelType
Superclass::OutputImageRegionType OutputImageRegionType