Skip to content

Commit fcc0c0a

Browse files
authored
Merge pull request #6368 from hjmjohnson/ingest-FixedPointInverseDisplacementField
ENH: Ingest ITKFixedPointInverseDisplacementField into Modules/Filtering
2 parents c25479f + 3d1f13c commit fcc0c0a

14 files changed

Lines changed: 549 additions & 50 deletions
Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
project(FixedPointInverseDisplacementField)
2+
3+
itk_module_impl()
Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,16 @@
1+
# FixedPointInverseDisplacementField
2+
3+
In-tree ITK module that computes the inverse of a displacement field using the
4+
fixed-point iteration of Chen et al. (*A simple fixed-point approach to invert
5+
a deformation field*, Medical Physics 35(1):81, 2008). Given a field mapping
6+
space A → B, it produces the field mapping B → A.
7+
8+
Flagship class: `itk::FixedPointInverseDisplacementFieldImageFilter`.
9+
10+
## Origin
11+
12+
Ingested from the standalone remote module
13+
[**InsightSoftwareConsortium/ITKFixedPointInverseDisplacementField**](https://github.com/InsightSoftwareConsortium/ITKFixedPointInverseDisplacementField)
14+
on 2026-06-03 via the v4 ingestion pipeline. The upstream repository will be
15+
archived read-only after this PR merges; it remains reachable at the URL above
16+
for historical reference.
Lines changed: 175 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,175 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
#ifndef itkFixedPointInverseDisplacementFieldImageFilter_h
19+
#define itkFixedPointInverseDisplacementFieldImageFilter_h
20+
21+
22+
#include "itkImageToImageFilter.h"
23+
24+
#include "itkWarpVectorImageFilter.h"
25+
#include "itkVectorLinearInterpolateImageFunction.h"
26+
#include "itkImageRegionIterator.h"
27+
#include "itkTimeProbe.h"
28+
29+
namespace itk
30+
{
31+
32+
/** \class FixedPointInverseDisplacementFieldImageFilter
33+
*
34+
* \brief Computes the inverse of a Displacement field using a fixed point iteration scheme.
35+
*
36+
* FixedPointInverseDisplacementFieldImageFilter takes a Displacement field as input and
37+
* computes the Displacement field that is its inverse. If the input Displacement
38+
* field was mapping coordinates from a space A into a space B, the output of
39+
* this filter will map coordinates from the space B into the space A.
40+
*
41+
* To compute the inverse of the given Displacement field, the fixed point algorithm by
42+
* Mingli Chen, Weiguo Lu, Quan Chen, Knneth J. Ruchala and Gusavo H. Olivera
43+
* described in the paper
44+
*
45+
\verbatim
46+
"A simple fixed-point approach to invert a Displacement field",
47+
Medical Physics, vol. 35, issue 1, p. 81
48+
\endverbatim
49+
* is applied.
50+
*
51+
* \author Marcel Lüthi, Computer Science Department, University of Basel
52+
*
53+
* \ingroup ImageToImageFilter
54+
* \ingroup FixedPointInverseDisplacementField
55+
*/
56+
57+
template <typename TInputImage, typename TOutputImage>
58+
class ITK_TEMPLATE_EXPORT FixedPointInverseDisplacementFieldImageFilter
59+
: public ImageToImageFilter<TInputImage, TOutputImage>
60+
{
61+
public:
62+
ITK_DISALLOW_COPY_AND_MOVE(FixedPointInverseDisplacementFieldImageFilter);
63+
64+
/** Standard class type alias. */
65+
using Self = FixedPointInverseDisplacementFieldImageFilter;
66+
using Superclass = ImageToImageFilter<TInputImage, TOutputImage>;
67+
using Pointer = SmartPointer<Self>;
68+
using ConstPointer = SmartPointer<const Self>;
69+
70+
/** Method for creation through the object factory. */
71+
itkNewMacro(Self);
72+
73+
/** Run-time type information (and related methods). */
74+
itkTypeMacro(FixedPointInverseDisplacementFieldImageFilter, ImageToImageFilter);
75+
76+
77+
/** Some type alias. */
78+
using InputImageType = TInputImage;
79+
using InputImageConstPointer = typename InputImageType::ConstPointer;
80+
using InputImagePointer = typename InputImageType::Pointer;
81+
using InputImagePointType = typename InputImageType::PointType;
82+
using InputImageRegionType = typename InputImageType::RegionType;
83+
using InputImageSpacingType = typename InputImageType::SpacingType;
84+
using InputImageIndexType = typename InputImageType::IndexType;
85+
86+
87+
using OutputImageType = TOutputImage;
88+
using OutputImagePointer = typename OutputImageType::Pointer;
89+
using OutputImagePixelType = typename OutputImageType::PixelType;
90+
using OutputImagePointType = typename OutputImageType::PointType;
91+
using OutputImageIndexType = typename OutputImageType::IndexType;
92+
using OutputImageValueType = typename OutputImagePixelType::ValueType;
93+
using OutputImageSizeType = typename OutputImageType::SizeType;
94+
using OutputImageSpacingType = typename OutputImageType::SpacingType;
95+
using OutputImageOriginPointType = typename TOutputImage::PointType;
96+
using TimeType = TimeProbe;
97+
98+
99+
/** Number of dimensions. */
100+
static constexpr unsigned int ImageDimension = TOutputImage::ImageDimension;
101+
102+
using InputConstIterator = ImageRegionConstIterator<InputImageType>;
103+
using InputIterator = ImageRegionIterator<InputImageType>;
104+
using OutputConstIterator = ImageRegionConstIterator<OutputImageType>;
105+
using OutputIterator = ImageRegionIterator<OutputImageType>;
106+
107+
using VectorWarperType = WarpVectorImageFilter<TOutputImage, TInputImage, TOutputImage>;
108+
using InterpolatorVectorType = typename NumericTraits<typename TOutputImage::PixelType>::RealType;
109+
using OutputVectorType = typename TOutputImage::PixelType;
110+
using FieldInterpolatorType = VectorLinearInterpolateImageFunction<TInputImage, double>;
111+
using FieldInterpolatorPointer = typename FieldInterpolatorType::Pointer;
112+
using FieldInterpolatorOutputType = typename FieldInterpolatorType::OutputType;
113+
114+
115+
itkSetMacro(NumberOfIterations, unsigned int);
116+
itkGetConstMacro(NumberOfIterations, unsigned int);
117+
118+
119+
/** Set the size of the output image. */
120+
itkSetMacro(Size, OutputImageSizeType);
121+
/** Get the size of the output image. */
122+
itkGetConstReferenceMacro(Size, OutputImageSizeType);
123+
124+
/** Set the output image spacing. */
125+
itkSetMacro(OutputSpacing, OutputImageSpacingType);
126+
virtual void
127+
SetOutputSpacing(const double * values);
128+
129+
/** Get the output image spacing. */
130+
itkGetConstReferenceMacro(OutputSpacing, OutputImageSpacingType);
131+
132+
/** Set the output image origin. */
133+
itkSetMacro(OutputOrigin, OutputImageOriginPointType);
134+
virtual void
135+
SetOutputOrigin(const double * values);
136+
137+
#ifdef ITK_USE_CONCEPT_CHECKING
138+
/** Begin concept checking */
139+
itkConceptMacro(OutputHasNumericTraitsCheck, (Concept::HasNumericTraits<OutputImageValueType>));
140+
itkConceptMacro(SameDimensionCheck,
141+
(Concept::SameDimension<TInputImage::ImageDimension, TOutputImage::ImageDimension>));
142+
143+
/** End concept checking */
144+
#endif
145+
146+
protected:
147+
FixedPointInverseDisplacementFieldImageFilter();
148+
~FixedPointInverseDisplacementFieldImageFilter() override = default;
149+
150+
void
151+
PrintSelf(std::ostream & os, Indent indent) const override;
152+
153+
void
154+
GenerateData() override;
155+
void
156+
GenerateOutputInformation() override;
157+
158+
void
159+
GenerateInputRequestedRegion() override;
160+
161+
unsigned int m_NumberOfIterations{ 5 };
162+
163+
164+
private:
165+
OutputImageSizeType m_Size; // Size of the output image
166+
OutputImageSpacingType m_OutputSpacing; // output image spacing
167+
OutputImageOriginPointType m_OutputOrigin; // output image origin
168+
};
169+
170+
} // end namespace itk
171+
172+
#ifndef ITK_MANUAL_INSTANTIATION
173+
# include "itkFixedPointInverseDisplacementFieldImageFilter.hxx"
174+
#endif
175+
#endif
Lines changed: 192 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,192 @@
1+
/*=========================================================================
2+
*
3+
* Copyright NumFOCUS
4+
*
5+
* Licensed under the Apache License, Version 2.0 (the "License");
6+
* you may not use this file except in compliance with the License.
7+
* You may obtain a copy of the License at
8+
*
9+
* https://www.apache.org/licenses/LICENSE-2.0.txt
10+
*
11+
* Unless required by applicable law or agreed to in writing, software
12+
* distributed under the License is distributed on an "AS IS" BASIS,
13+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14+
* See the License for the specific language governing permissions and
15+
* limitations under the License.
16+
*
17+
*=========================================================================*/
18+
#ifndef itkFixedPointInverseDisplacementFieldImageFilter_hxx
19+
#define itkFixedPointInverseDisplacementFieldImageFilter_hxx
20+
21+
#include <iostream>
22+
23+
namespace itk
24+
{
25+
//----------------------------------------------------------------------------
26+
// Constructor
27+
template <typename TInputImage, typename TOutputImage>
28+
FixedPointInverseDisplacementFieldImageFilter<TInputImage,
29+
TOutputImage>::FixedPointInverseDisplacementFieldImageFilter()
30+
{
31+
32+
m_OutputSpacing.Fill(1.0);
33+
m_OutputOrigin.Fill(0.0);
34+
for (unsigned int i = 0; i < ImageDimension; i++)
35+
{
36+
m_Size[i] = 0;
37+
}
38+
}
39+
40+
41+
/**
42+
* Set the output image spacing.
43+
*/
44+
template <typename TInputImage, typename TOutputImage>
45+
void
46+
FixedPointInverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::SetOutputSpacing(const double * spacing)
47+
{
48+
OutputImageSpacingType s(spacing);
49+
this->SetOutputSpacing(s);
50+
}
51+
52+
/**
53+
* Set the output image origin.
54+
*/
55+
template <typename TInputImage, typename TOutputImage>
56+
void
57+
FixedPointInverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::SetOutputOrigin(const double * origin)
58+
{
59+
OutputImageOriginPointType p(origin);
60+
this->SetOutputOrigin(p);
61+
}
62+
63+
64+
//----------------------------------------------------------------------------
65+
template <typename TInputImage, typename TOutputImage>
66+
void
67+
FixedPointInverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::GenerateData()
68+
{
69+
70+
const InputImageType * inputPtr = this->GetInput(0);
71+
OutputImageType * outputPtr = this->GetOutput(0);
72+
73+
InputConstIterator InputIt = InputConstIterator(inputPtr, inputPtr->GetRequestedRegion());
74+
75+
// We allocate a Displacement field that holds the output
76+
// and initialize it to 0
77+
outputPtr->SetBufferedRegion(outputPtr->GetRequestedRegion());
78+
outputPtr->Allocate();
79+
80+
OutputIterator outputIt = OutputIterator(outputPtr, outputPtr->GetRequestedRegion());
81+
OutputImagePixelType zero_pt;
82+
zero_pt.Fill(0);
83+
outputPtr->FillBuffer(zero_pt);
84+
85+
86+
// In the fixed point iteration, we will need to access non-grid points.
87+
// Currently, the best interpolator that is supported by itk for vector
88+
// images is the linear interpolater.
89+
typename FieldInterpolatorType::Pointer vectorInterpolator = FieldInterpolatorType::New();
90+
vectorInterpolator->SetInputImage(inputPtr);
91+
92+
93+
// Finally, perform the fixed point iteration.
94+
InputImagePointType mappedPt;
95+
OutputImagePointType pt;
96+
97+
for (unsigned int i = 0; i <= m_NumberOfIterations; i++)
98+
{
99+
for (outputIt.GoToBegin(); !outputIt.IsAtEnd(); ++outputIt)
100+
{
101+
const OutputImageIndexType index = outputIt.GetIndex();
102+
outputPtr->TransformIndexToPhysicalPoint(index, pt);
103+
const OutputImagePixelType displacement = outputIt.Get();
104+
for (unsigned int j = 0; j < ImageDimension; j++)
105+
{
106+
mappedPt[j] = pt[j] + displacement[j];
107+
}
108+
109+
110+
if (vectorInterpolator->IsInsideBuffer(mappedPt))
111+
{
112+
InterpolatorVectorType val = vectorInterpolator->Evaluate(mappedPt);
113+
OutputVectorType outputVector;
114+
for (unsigned int j = 0; j < ImageDimension; j++)
115+
{
116+
outputVector[j] = -val[j];
117+
}
118+
outputIt.Set(outputVector);
119+
}
120+
}
121+
}
122+
}
123+
124+
125+
/**
126+
* Inform pipeline of required output region
127+
*/
128+
template <typename TInputImage, typename TOutputImage>
129+
void
130+
FixedPointInverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::GenerateOutputInformation()
131+
{
132+
// call the superclass' implementation of this method
133+
Superclass::GenerateOutputInformation();
134+
135+
// get pointers to the input and output
136+
OutputImageType * outputPtr = this->GetOutput();
137+
if (!outputPtr)
138+
{
139+
return;
140+
}
141+
142+
// Set the size of the output region
143+
typename TOutputImage::RegionType outputLargestPossibleRegion;
144+
outputLargestPossibleRegion.SetSize(m_Size);
145+
outputPtr->SetLargestPossibleRegion(outputLargestPossibleRegion);
146+
147+
// Set spacing and origin
148+
outputPtr->SetSpacing(m_OutputSpacing);
149+
outputPtr->SetOrigin(m_OutputOrigin);
150+
151+
return;
152+
}
153+
154+
155+
//----------------------------------------------------------------------------
156+
template <typename TInputImage, typename TOutputImage>
157+
void
158+
FixedPointInverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion()
159+
{
160+
// call the superclass's implementation of this method
161+
Superclass::GenerateInputRequestedRegion();
162+
163+
if (!this->GetInput())
164+
{
165+
return;
166+
}
167+
168+
// get pointers to the input and output
169+
auto * inputPtr = const_cast<InputImageType *>(this->GetInput());
170+
171+
// Request the entire input image
172+
InputImageRegionType inputRegion;
173+
inputRegion = inputPtr->GetLargestPossibleRegion();
174+
inputPtr->SetRequestedRegion(inputRegion);
175+
}
176+
177+
//----------------------------------------------------------------------------
178+
template <typename TInputImage, typename TOutputImage>
179+
void
180+
FixedPointInverseDisplacementFieldImageFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream & os,
181+
Indent indent) const
182+
{
183+
184+
Superclass::PrintSelf(os, indent);
185+
186+
os << indent << "Number of iterations: " << m_NumberOfIterations << std::endl;
187+
os << std::endl;
188+
}
189+
190+
} // end namespace itk
191+
192+
#endif

0 commit comments

Comments
 (0)