InsightSoftwareConsortium / InsightSoftwareConsortium/ITK

PERF: Restructure `BSplineTransform::TransformPoint` inner loop for a SIMD-friendly gather pattern

Open
#6,746 0 comments 1 reaction 0 assignees View on GitHub
area:Registration type:Performance
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

Open the contributing 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

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.