BSplineTransform::TransformPoint walks SpaceDimension independent
ImageScanlineConstIterators 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:
ImageScanlineConstIterator<ImageType> 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<ScalarType>(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.
BSplineTransform::TransformPointwalksSpaceDimensionindependentImageScanlineConstIterators in lockstep, so each iteration issuesSpaceDimensionscattered loads. The loop is memory-latency-bound rather thancompute-bound, and consequently regresses under wider SIMD instead of benefiting.
Current code — itkBSplineTransform.hxx:542-575
TransformPointstarts at line 510. The inner loop:ImageScanlineConstIterator<ImageType> 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<ScalarType>(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, whileAffine-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
SpaceDimensionimages into one small dense scratch buffer (the support regionis
(VSplineOrder+1)^VDimensionelements, small enough to stack-allocate), thenrun a flat fused-multiply-add over that buffer against
weights.Constraints:
indices[]population must be preserved exactly — callers depend on it.-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.