InsightSoftwareConsortium / InsightSoftwareConsortium/ITK
PERF: Restructure `BSplineTransform::TransformPoint` inner loop for a SIMD-friendly gather pattern
- Dominant language
- C++
- Stars
- 1.7k
- Forks
- 748
- Avg merge
- 1d 1h
- Merged PRs (30d)
- 64
Description
`BSplineTransform::TransformPoint` walks `SpaceDimension` independent
`ImageScanlineConstIterator`s in lockstep, so each iteration issues
`SpaceDimension` scattered loads. The loop is memory-latency-bound rather than
compute-bound, and consequently regresses under wider SIMD instead of benefiting.
Current code — itkBSplineTransform.hxx:542-575
`TransformPoint` starts at line 510. The inner loop:
```cpp
ImageScanlineConstIterator coeffIterator[SpaceDimension];
unsigned long counter = 0;
const ParametersValueType * basePointer = this->m_CoefficientImages[0]->GetBufferPointer();
for (unsigned int j = 0; j < SpaceDimension; ++j)
{
coeffIterator[j] = ImageScanlineConstIterator(this->m_CoefficientImages[j], supportRegion);
}
while (!coeffIterator[0].IsAtEnd())
{
while (!coeffIterator[0].IsAtEndOfLine())
{
// Multiply weight with coefficient
for (unsigned int j = 0; j < SpaceDimension; ++j)
{
outputPoint[j] += static_cast(weights[counter] * coeffIterator[j].Get());
}
// Populate the indices array
indices[counter] = &(coeffIterator[0].Value()) - basePointer;
// Go to next coefficient in the support region
++counter;
for (unsigned int j = 0; j < SpaceDimension; ++j)
{
++(coeffIterator[j]);
}
} // end scanline
for (unsigned int j = 0; j < SpaceDimension; ++j)
{
coeffIterator[j].NextLine();
}
}
```
The coefficient images are *separate* allocations, so `coeffIterator[j].Get()`
for j = 0..SpaceDimension-1 touches SpaceDimension unrelated cache lines per
scalar output. The compiler cannot turn this into a contiguous load.
Measured symptom
In a Resample benchmark sweep on a Sapphire Rapids Xeon w7-3545 (GCC 13.3),
BSpline-dominated variants **regress 12–17%** under `-march=x86-64-v4`, while
Affine-only variants **gain 17–20%**. The AVX-512 license-based frequency
throttle lowers the clock; because this loop is latency-bound, it pays the clock
cost and collects none of the width benefit.
Proposal and constraints
Restructure so coefficient reads form a contiguous or predictably strided
access — for example, pre-pack the support-region coefficients for all
`SpaceDimension` images into one small dense scratch buffer (the support region
is `(VSplineOrder+1)^VDimension` elements, small enough to stack-allocate), then
run a flat fused-multiply-add over that buffer against `weights`.
Constraints:
- Pure performance refactor: no API change, no behavioral change.
- `indices[]` population must be preserved exactly — callers depend on it.
- Correctness pinned by the existing BSplineTransform tests before and after.
- Re-benchmark at `-march=x86-64`, `-march=x86-64-v3`, and `-march=x86-64-v4`;
the v4 regression disappearing is the success criterion.
Related but distinct: #1278 (Enoki adoption) proposes a vectorization *library*
rather than restructuring this loop.
Contributor guide
Research direction
Start in itkBSplineTransform.hxx at BSplineTransform::TransformPoint around lines 510-575, and run the existing BSplineTransform tests before changing the loop. Benchmark the relevant Resample variants at -march=x86-64, -march=x86-64-v3, and -march=x86-64-v4. Done means preserving indices[] and behavior while eliminating the reported v4 regression.
Written by the indexing model from the issue text.
Assessment
- Tech stack
- cpp
- Domain
- performance
- Issue type
- Refactor
- Difficulty
- 5/5
- Estimated time
- Over a week
- Activity status
- Quiet
- Clarity
- Mostly clear
- Newbie friendliness
- 42/100