From ac5652bb68de08697a0726e1b11d1ad56df6956c Mon Sep 17 00:00:00 2001 From: Viacheslav Ustymenko Date: Wed, 8 Apr 2026 16:36:24 -0500 Subject: [PATCH 1/2] BUG: Implement integer-only Bresenham for BuildLine(Index, Index) Replace the floating-point delegation through BuildLine(Direction, length) with a direct integer-only N-dimensional Bresenham that guarantees exact start and end points by construction. The previous implementation converted integer endpoints to a floating-point direction vector, computed Chebyshev distance for the line length, then reconstructed LastIndex via rounding. This introduced endpoint error for lines where the dominant-axis projection differed from the Euclidean length (e.g., {0,0,0} to {250,250,1} would overshoot to {251,251,0}). Reference: classic Bresenham extended to N-D as used by scikit-image (line_nd), OpenCV (LineIterator), and VTK (vtkLine). (cherry picked from commit 59e3c0690d) --- .../Core/Common/include/itkBresenhamLine.hxx | 154 ++++++++---------- 1 file changed, 69 insertions(+), 85 deletions(-) diff --git a/Modules/Core/Common/include/itkBresenhamLine.hxx b/Modules/Core/Common/include/itkBresenhamLine.hxx index 0aa6b1675ae..7c237c9405a 100644 --- a/Modules/Core/Common/include/itkBresenhamLine.hxx +++ b/Modules/Core/Common/include/itkBresenhamLine.hxx @@ -21,125 +21,109 @@ #include "itkPoint.h" #include "itkMath.h" +#include +#include + namespace itk { template auto BresenhamLine::BuildLine(LType Direction, IdentifierType length) -> OffsetArray { - // copied from the line iterator - /** Variables that drive the Bresenham-Algorithm */ - // The dimension with the largest difference between start and end - unsigned int m_MainDirection; - - // Accumulated error for the other dimensions - IndexType m_AccumulateError; - - // Increment for the error for each step. Two times the difference between - // start and end - IndexType m_IncrementError; - - // If enough is accumulated for a dimension, the index has to be - // incremented. Will be the number of pixels in the line - IndexType m_MaximalError; - - // Direction of increment. -1 or 1 - IndexType m_OverflowIncrement; - - // After an overflow, the accumulated error is reduced again. Will be - // two times the number of pixels in the line - IndexType m_ReduceErrorAfterIncrement; + Direction.Normalize(); - OffsetArray result(length); + // The dimension with the largest absolute component + const unsigned int maxDistanceDimension = std::distance( + Direction.Begin(), std::max_element(Direction.Begin(), Direction.End(), [](const auto a, const auto b) { + return itk::Math::Absolute(a) < itk::Math::Absolute(b); + })); - IndexType m_CurrentImageIndex, LastIndex; + // compute actual line length because the shorter distance + // the larger deviation due to rounding to integers + const IdentifierType mainDirectionLen = length - 1; + const float euclideanLineLen = mainDirectionLen / itk::Math::Absolute(Direction[maxDistanceDimension]); - Direction.Normalize(); // we are going to start at 0 - m_CurrentImageIndex.Fill(0); - constexpr IndexType StartIndex = { { 0 } }; + constexpr IndexType StartIndex{ 0 }; + IndexType LastIndex; for (unsigned int i = 0; i < VDimension; ++i) { - LastIndex[i] = (IndexValueType)(length * Direction[i]); + // round to closest voxel center to minimize the deviation from Direction + LastIndex[i] = Math::RoundHalfIntegerUp(euclideanLineLen * Direction[i]); } - // Find the dominant direction - IndexValueType maxDistance = 0; - unsigned int maxDistanceDimension = 0; - for (unsigned int i = 0; i < VDimension; ++i) + + const IndexArray indices = this->BuildLine(StartIndex, LastIndex); + OffsetArray offsets; + offsets.reserve(indices.size()); + for (const IndexType & index : indices) { - auto distance = static_cast(itk::Math::abs(LastIndex[i])); - if (distance > maxDistance) - { - maxDistance = distance; - maxDistanceDimension = i; - } - m_IncrementError[i] = 2 * distance; - m_OverflowIncrement[i] = (LastIndex[i] < 0 ? -1 : 1); + offsets.push_back(index - StartIndex); } - m_MainDirection = maxDistanceDimension; - m_MaximalError.Fill(maxDistance); - m_ReduceErrorAfterIncrement.Fill(2 * maxDistance); - m_AccumulateError.Fill(0); - unsigned int steps = 1; - result[0] = m_CurrentImageIndex - StartIndex; - while (steps < length) - { - // This part is from ++ in LineConstIterator - // We need to modify m_AccumulateError, m_CurrentImageIndex, m_IsAtEnd - for (unsigned int i = 0; i < VDimension; ++i) - { - if (i == m_MainDirection) - { - m_CurrentImageIndex[i] += m_OverflowIncrement[i]; - } - else - { - m_AccumulateError[i] += m_IncrementError[i]; - if (m_AccumulateError[i] >= m_MaximalError[i]) - { - m_CurrentImageIndex[i] += m_OverflowIncrement[i]; - m_AccumulateError[i] -= m_ReduceErrorAfterIncrement[i]; - } - } - } - - result[steps] = m_CurrentImageIndex - StartIndex; // produce an offset - ++steps; - } - return (result); + return offsets; } template auto BresenhamLine::BuildLine(IndexType p0, IndexType p1) -> IndexArray { - itk::Point point0; - itk::Point point1; - IdentifierType maxDistance = 0; + // Integer-only N-dimensional Bresenham, guaranteeing exact start and end + // points. This avoids floating-point direction conversion entirely. + // Reference: classic Bresenham extended to N-D as used by scikit-image, + // OpenCV, and VTK. + + // Compute displacements and find the dominant axis + IndexType absDeltas; + IndexType step; // +1 or -1 per dimension + IndexValueType maxAbsDelta = 0; + unsigned int mainAxis = 0; for (unsigned int i = 0; i < VDimension; ++i) { - point0[i] = p0[i]; - point1[i] = p1[i]; - IdentifierType distance = itk::Math::abs(p0[i] - p1[i]) + 1; - if (distance > maxDistance) + const IndexValueType delta = p1[i] - p0[i]; + step[i] = (delta >= 0) ? 1 : -1; + absDeltas[i] = static_cast(itk::Math::Absolute(delta)); + + if (absDeltas[i] > maxAbsDelta) { - maxDistance = distance; + maxAbsDelta = absDeltas[i]; + mainAxis = i; } } - OffsetArray offsets = this->BuildLine(point1 - point0, maxDistance + 1); + // Number of pixels = steps along dominant axis + 1 + const IdentifierType numPixels = static_cast(maxAbsDelta) + 1; IndexArray indices; - indices.reserve(offsets.size()); // we might not have to use the last one - for (unsigned int i = 0; i < offsets.size(); ++i) + indices.reserve(numPixels); + + // Error accumulators for each secondary axis (2 * absDelta[i] per step, + // overflow when accumulated error >= maxAbsDelta) + auto accumulateError = MakeFilled(0); + auto currentIndex = p0; + + for (IdentifierType s = 0; s < numPixels; ++s) { - indices.push_back(p0 + offsets[i]); - if (indices.back() == p1) + indices.push_back(currentIndex); + + // Advance along each dimension + for (unsigned int i = 0; i < VDimension; ++i) { - break; + if (i == mainAxis) + { + currentIndex[i] += step[i]; + } + else + { + accumulateError[i] += 2 * absDeltas[i]; + if (accumulateError[i] >= maxAbsDelta) + { + currentIndex[i] += step[i]; + accumulateError[i] -= 2 * maxAbsDelta; + } + } } } + return indices; } From 2003fb40050b9d4365609768440a16f6c39b3266 Mon Sep 17 00:00:00 2001 From: Viacheslav Ustymenko Date: Thu, 9 Apr 2026 01:05:55 +0300 Subject: [PATCH 2/2] BUG: BresenhamLine - use double precision floats for ending index computation MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Use explicit static_cast() on both operands of the division computing euclideanLineLen so the division evaluates in double precision, not float (LType = Vector). Assisted-by: Greptile — identified float-precision division in BuildLine (cherry picked from commit b7248b1ca4) --- Modules/Core/Common/include/itkBresenhamLine.hxx | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Modules/Core/Common/include/itkBresenhamLine.hxx b/Modules/Core/Common/include/itkBresenhamLine.hxx index 7c237c9405a..1cea4756e89 100644 --- a/Modules/Core/Common/include/itkBresenhamLine.hxx +++ b/Modules/Core/Common/include/itkBresenhamLine.hxx @@ -35,13 +35,14 @@ BresenhamLine::BuildLine(LType Direction, IdentifierType length) -> // The dimension with the largest absolute component const unsigned int maxDistanceDimension = std::distance( Direction.Begin(), std::max_element(Direction.Begin(), Direction.End(), [](const auto a, const auto b) { - return itk::Math::Absolute(a) < itk::Math::Absolute(b); + return itk::Math::abs(a) < itk::Math::abs(b); })); // compute actual line length because the shorter distance // the larger deviation due to rounding to integers const IdentifierType mainDirectionLen = length - 1; - const float euclideanLineLen = mainDirectionLen / itk::Math::Absolute(Direction[maxDistanceDimension]); + const double euclideanLineLen = + static_cast(mainDirectionLen) / static_cast(itk::Math::abs(Direction[maxDistanceDimension])); // we are going to start at 0 constexpr IndexType StartIndex{ 0 }; @@ -81,7 +82,7 @@ BresenhamLine::BuildLine(IndexType p0, IndexType p1) -> IndexArray { const IndexValueType delta = p1[i] - p0[i]; step[i] = (delta >= 0) ? 1 : -1; - absDeltas[i] = static_cast(itk::Math::Absolute(delta)); + absDeltas[i] = static_cast(itk::Math::abs(delta)); if (absDeltas[i] > maxAbsDelta) {