BUG: Make itkBresenhamLine::BuildLine finish at the correct end point

[slavust]

Overload of the BuildLine(Index, Index) function could end at index which is different from specified one. The issue was that when reconstructing LastIndex from the direction vector we didn't use actual Euclidean line length, but used Chebyshev distance instead. This led to rounding to integers closer to the line start, leading to bigger deviation at the line end.

The alternative approach would be to move the algorithm implementation to the BuildLine(Index, Index) and work fully on integers there, but this would require more code changes.
UPDATE: AI agent actually suggested code with algorithm implementation in BuildLine(Index, Index) overload. Since AI is allowed to do such modification, I'll proceed with it.

Bresenham line algorithm works on integer indices, and we cannot avoid rounding errors while having floating-point direction. But it's still possible to always finish at the user-specified point in the second function overload.

PR Checklist

• No API changes were made (or the changes have been approved)
• No major design changes were made (or the changes have been approved)
• Added test (or behavior not changed)
• Updated API documentation (or API not changed)
• Added license to new files (if any)
• Added Python wrapping to new files (if any) as described in ITK Software Guide Section 9.5
• Added ITK examples for all new major features (if any)

Refer to the ITK Software Guide for
further development details if necessary.

[github-actions]

Thank you for contributing a pull request! 🙏

Welcome to the ITK community! 🤗👋☀️

We are glad you are here and appreciate your contribution. Please keep in mind our community participation guidelines. 📜
More support and guidance on the contribution process can be found in our contributing guide. 📖

This is an automatic message. Allow for time for the ITK community to be able to read the pull request and comment
on it.

[hjmjohnson]

Test-before-fix verification

Confirmed that the new FinishAtEndIndex test fails on the unfixed code (current main), demonstrating the bug:

Expected: { 250, 250, 1 }
  Actual: { 251, 251, 0 }

The BuildLine(start={0,0,0}, end={250,250,1}) overshoots in X/Y and lands on the wrong Z plane because LastIndex is computed using the projected length on the dominant axis rather than the actual Euclidean line length. The PR's fix (euclideanLineLen * Direction[i] with proper rounding) corrects this.

Existing tests (BuildLineFromVector, BuildLineFromIndices) continue to pass.

[hjmjohnson]

Algorithm review: comparison with reference implementations and correctness analysis

Classic Bresenham properties

The classic Bresenham line algorithm (1965) is integer-only and guarantees hitting both endpoints exactly. It tracks an integer error term, steps one pixel at a time along the dominant axis, and decides whether to step in secondary axes based on accumulated error. The N-dimensional generalization preserves this: identify the dominant axis (largest absolute displacement), step along it, accumulate integer error terms for each of the other N-1 axes.

How other libraries handle this

Library	Approach	Endpoint guarantee
scikit-image (line_nd)	Integer endpoints directly, np.linspace along dominant axis	Yes, by construction
OpenCV (LineIterator)	Classic integer Bresenham, two-endpoint formulation	Yes
VTK	Parametric line; rasterization done separately with integer endpoints	Yes

The common pattern: the integer-endpoint-to-endpoint formulation is primary. Direction+length overloads are convenience wrappers, not the other way around. ITK inverts this — BuildLine(p0, p1) delegates to BuildLine(direction, length) via floating-point conversion, which is the root cause of the endpoint accuracy issue.

Current ITK bug analysis

In BuildLine(p0, p1):

1. Computes maxDistance = max(Abs(p0[i] - p1[i]) + 1) — the +1 conflates pixel count (fence-post) with displacement (fence-rail)
2. Passes maxDistance + 1 as length to overload 1 — double-counting
3. Overload 1 computes LastIndex[i] = (IndexValueType)(length * Direction[i]) using C-style cast (truncation toward zero)
4. Relies on an early break when indices.back() == p1 — but floating-point rounding means the Bresenham path may never exactly hit p1, so the break never fires

Test results on current main confirm the bug in both directions:

• Forward {0,0,0}→{250,250,1}: ends at {251,251,0} (overshoots X/Y, wrong Z)
• Reverse {250,250,1}→{0,0,0}: ends at {-1,-1,1} (undershoots X/Y, wrong Z)

PR #6024 fix assessment

The fix is mathematically sound for the forward case:

• euclideanLineLen = mainDirectionLen / Direction[maxDistanceDimension] correctly recovers the Euclidean distance
• RoundHalfIntegerUp for LastIndex is more accurate than truncation
• Removing the +1 from distance calculation fixes the off-by-two
• Removing the early break is correct since the line should now always have the right length

⚠️ Potential issue: std::max_element without absolute values

The PR computes maxDistanceDimension using:

std::distance(Direction.Begin(), std::max_element(Direction.Begin(), Direction.End()));

This finds the largest signed component, not the largest absolute component. For a direction like (-0.707, 0.707, 0.01), max_element returns dimension 1 (0.707), but for (-0.9, 0.1, 0.01), it returns dimension 1 (0.1) instead of dimension 0 (|-0.9|). This would select the wrong dominant axis for negative-dominant directions and compute an incorrect euclideanLineLen.

The original code correctly uses Abs(LastIndex[i]) to find the dominant axis. The fix should use absolute values when finding maxDistanceDimension.

Recommendation

1. The fix addresses a real bug — the endpoint guarantee is a fundamental Bresenham property that the current code violates
2. The std::max_element issue should be fixed before merge — use absolute values for dominant axis selection
3. Long-term: consider implementing integer-only Bresenham directly in BuildLine(p0, p1) rather than delegating through floating-point — this is what every other major library does and eliminates the entire class of rounding issues
4. Add test for negative directions (reversed endpoints) to catch the max_element issue

[hjmjohnson]

Integer-only Bresenham for BuildLine(Index, Index)

Added commit 7efbaf5 which replaces the floating-point delegation through BuildLine(Direction, length) with a direct integer-only N-dimensional Bresenham implementation. This eliminates all floating-point rounding issues and guarantees exact endpoints by construction.

Reference implementations

This follows the standard approach used by other major imaging libraries:

• scikit-image line_nd — integer-endpoint formulation is primary
• OpenCV LineIterator — classic integer Bresenham, two-endpoint formulation
• VTK vtkLine — parametric line with integer rasterization

In all three, the integer-endpoint-to-endpoint path is the primary/correct formulation. The direction+length overload (BuildLine(Direction, length)) is left unchanged as a separate code path for callers who genuinely have a direction vector.

New tests (7 added, 0 existing tests modified)

Test	What it verifies
AxisAligned2D	Horizontal/vertical lines hit exact endpoints, stay on-axis
SinglePixelLine	p0 == p1 returns exactly 1 pixel
Connectivity2D	All consecutive pixels are 8-connected (Chebyshev ≤ 1)
Connectivity3D	All consecutive pixels are 26-connected (Chebyshev ≤ 1)
PixelCountMatchesChebyshev	Pixel count = Chebyshev distance + 1 for various cases
ReverseSymmetry	Forward/reverse lines have same length, correct endpoints
SteepDiagonal3D	The original bug case {0,0,0}→{250,250,1} in both directions

All 11 tests pass (4 pre-existing + 7 new).

[hjmjohnson]

Performance analysis: upstream/main vs integer-only Bresenham

Benchmarked with -O2 on the same machine, 1000 iterations per measurement.

Benchmark	main (ms/call)	PR (ms/call)	Speedup
BuildLine(Index,Index) 2D 1000×500	0.0071	0.0052	1.37×
BuildLine(Index,Index) 2D 10×10	0.000167	0.000100	1.67×
BuildLine(Index,Index) 3D 500×300×200	0.0046	0.0030	1.56×
BuildLine(Index,Index) 3D 250×250×1	0.0029	0.0021	1.38×
BuildLine(Index,Index) 3D reverse	0.0029	0.0021	1.38×
BuildLine(Dir,len) 2D len=1001	0.0050	0.0049	~same
BuildLine(Dir,len) 3D len=501	0.0034	0.0034	~same
1000 short 3D lines (batch)	0.1737	0.0845	2.05×

Summary

• BuildLine(Index, Index) is 1.4–2× faster due to eliminating float Point allocation, direction vector normalization (sqrt + division), and float→int rounding
• BuildLine(Direction, length) is unchanged (expected — same code path)
• Batch of 1000 short lines shows the largest gain (2×) since per-call overhead of float conversion dominates for short lines
• No regressions in any benchmark

[hjmjohnson]

Branch reorganization

Reorganized the 6 commits into 2 clean concept-based commits on BRAINSia/ITK:bresenham-line-accuracy:

1. 54d08d63e8 ENH: Add comprehensive tests — all 16 GTest cases (test-only, no code changes). The FinishAtEndIndex test demonstrates the bug on unfixed code.
2. 59e3c0690d BUG: Integer-only Bresenham — the algorithm fix (.hxx only). All tests pass after this commit.

Original 6 commits backed up at BRAINSia/ITK:bresenham-line-accuracy-backup.

Review comment resolutions

Comment	Status	Resolution
abs could be taken only once before the main loop	Resolved	absDeltas[i] is precomputed before the loop — no Absolute() call inside
call the second method here	Resolved	BuildLine(Direction, length) delegates to BuildLine(Index, Index)
check each voxel / actually there might be difference with odd directions	Resolved	Confirmed that pixel-by-pixel reversal symmetry is NOT guaranteed — the error accumulator tie-breaking produces different intermediate pixels for some directions (e.g., {10,20,3}→{200,50,80} differs at index 95). Added a comment documenting this property. The test checks size, endpoints, and connectivity — all of which ARE guaranteed.

Authorship

• Author: Viacheslav Ustymenko (original contributor)
• Committer: Hans J. Johnson
• Co-Authored-By: Hans J. Johnson

Test results

18/18 BresenhamLine GTest cases pass (0.14 sec total).

[hjmjohnson]

@slavust THANK YOU for your attention today. Now that the changes are all passing and meeting your requirements, I have rebased and reorganized to keep the git history as clean as possible ( and improve git bisect performance).

[greptile-apps]

Greptile Summary

This PR fixes BresenhamLine::BuildLine(Index, Index) to always terminate at the user-specified endpoint by replacing the original floating-point direction-based approach with a pure integer N-dimensional Bresenham algorithm. BuildLine(Direction, length) is refactored to compute a correct LastIndex via the Euclidean line length (rather than the Chebyshev distance used before) and then delegates to the rewritten integer overload. The integer algorithm is mathematically correct: the dominant axis advances exactly maxAbsDelta steps and each secondary axis accumulates exactly absDeltas[i] overflows, guaranteeing p1 is the final point.

Confidence Score: 4/5

PR is safe to merge with one minor precision concern worth addressing.

The integer Bresenham rewrite is mathematically correct and well-tested. The only finding is a P2 single-precision float issue in BuildLine(Direction, length) that could silently return length-1 elements instead of length for very large inputs — easy one-line fix with double.

itkBresenhamLine.hxx line 44 — float precision for euclideanLineLen

Vulnerabilities

No security concerns identified.

Important Files Changed

Filename	Overview
Modules/Core/Common/include/itkBresenhamLine.hxx	Rewrites BuildLine(p0,p1) as integer-only N-D Bresenham (correct) and refactors BuildLine(Direction,length) to delegate to it; single-precision float for euclideanLineLen may lose precision for very large lines
Modules/Core/Common/test/itkBresenhamLineGTest.cxx	Comprehensive new tests covering 8 octants, connectivity, uniqueness, axis-aligned lines, negative coordinates, steep/shallow slopes, and very long lines — good coverage of the bug fix and edge cases

Flowchart

%%{init: {'theme': 'neutral'}}%%
flowchart TD
    A["BuildLine(Direction, length)"] --> B["Normalize Direction\nFind maxDistanceDimension"]
    B --> C["euclideanLineLen = (length-1) / |Dir[max]|"]
    C --> D["LastIndex[i] = RoundHalfIntegerUp(euclideanLineLen * Dir[i])"]
    D --> E["BuildLine(StartIndex=0, LastIndex)"]
    E --> F

    G["BuildLine(p0, p1)"] --> F["Compute absDeltas, step\nFind mainAxis, maxAbsDelta"]
    F --> H["numPixels = maxAbsDelta + 1"]
    H --> I["Loop s = 0..numPixels-1"]
    I --> J["push currentIndex"]
    J --> K["Advance mainAxis by step[mainAxis]"]
    K --> L["For each secondary axis:\naccumError += 2*absDelta\nif accum >= maxAbsDelta: index++, accum -= 2*maxAbsDelta"]
    L --> M{s < numPixels-1?}
    M -- yes --> I
    M -- no --> N["Return indices\n(first = p0, last = p1)"]

Loading

_{Reviews (1): Last reviewed commit: "BUG: Implement integer-only Bresenham fo..." | Re-trigger Greptile}