Fast Voxel Traversal 2D

Question

I'm trying traverse all the cells that a line goes through. I've found the Fast Voxel Traversal Algorithm that seems to fit my needs, but I'm currently finding to be inaccurate. Below is a graph with a red line and points as voxel coordinates that the algorithm gives. As you can see it is almost correct except for the (4, 7) point, as it should be (5,6). I'm not sure if i'm implementing the algorithm correctly either so I've included it in Python. So i guess my question is my implementation correct or is there a better algo to this?

Thanks

def getVoxelTraversalPts(strPt, endPt, geom):
  Max_Delta = 1000000.0
  #origin
  x0 = geom[0]
  y0 = geom[3]

  (sX, sY) = (strPt[0], strPt[1])
  (eX, eY) = (endPt[0], endPt[1])
  dx = geom[1]
  dy = geom[5]

  sXIndex = ((sX - x0) / dx)
  sYIndex = ((sY - y0) / dy)
  eXIndex = ((eX - sXIndex) / dx) + sXIndex
  eYIndex = ((eY - sYIndex) / dy) + sYIndex

  deltaX = float(eXIndex - sXIndex)
  deltaXSign = 1 if deltaX > 0 else -1 if deltaX < 0 else 0
  stepX = deltaXSign

  tDeltaX = min((deltaXSign / deltaX), Max_Delta) if deltaXSign != 0 else Max_Delta
  maxX = tDeltaX * (1 - sXIndex + int(sXIndex)) if deltaXSign > 0 else tDeltaX * (sXIndex - int(sXIndex))

  deltaY = float(eYIndex - sYIndex)
  deltaYSign = 1 if deltaY > 0 else -1 if deltaY < 0 else 0
  stepY = deltaYSign

  tDeltaY = min(deltaYSign / deltaY, Max_Delta) if deltaYSign != 0 else Max_Delta
  maxY = tDeltaY * (1 - sYIndex + int(sYIndex)) if deltaYSign > 0 else tDeltaY * (sYIndex - int(sYIndex))

  x = sXIndex
  y = sYIndex

  ptsIndexes = []
  pt = [round(x), round(y)]
  ptsIndexes.append(pt)
  prevPt = pt
  while True:
    if maxX < maxY:
        maxX += tDeltaX
        x += deltaXSign
    else:
        maxY += tDeltaY
        y += deltaYSign

    pt = [round(x), round(y)]
    if pt != prevPt:
        #print pt
        ptsIndexes.append(pt)
        prevPt = pt

    if maxX > 1 and maxY > 1:
        break

  return (ptsIndexes)

Answer 1

The voxels that you are walking start at 0.0, i.e. the first voxel spans space from 0.0 to 1.0, a not from -0.5 to 0.5 as you seem to be assuming. In other words, they are the ones marked with dashed line, and not the solid one.

If you want voxels to be your way, you will have to fix initial maxX and maxY calculations.

Answer 2

Ain't nobody got time to read the paper you posted and figure out if you've implemented it correctly.

Here's a question, though. Is the algorithm you've used (a) actually meant to determine all the cells that a line passes through or (b) form a decent voxel approximation of a straight line between two points?

I'm more familiar with Bresenham's line algorithm which performs (b). Here's a picture of it in action:

Note that the choice of cells is "aesthetic", but omits certain cells the line passes through. Including these would make the line "uglier".

I suspect a similar thing is going on with your voxel line algorithm. However, looking at your data and the Bresenham image suggests a simple solution. Walk along the line of discovered cells, but, whenever you have to make a diagonal step, consider the two intermediate cells. You can then use a line-rectangle intersection algorithm (see here) to determine which of the candidate cells should have, but wasn't, included.

Answer 3

I guess just to be complete, I decided to use a different algo. the one referenced here dtb's answer on another question.

here's the implementation

def getIntersectPts(strPt, endPt, geom=[0,1,0,0,0,1]):
'''
    Find intersections pts for every half cell size
    ** cell size has only been tested with 1

    Returns cell coordinates that the line passes through
'''

x0 = geom[0]
y0 = geom[3]

(sX, sY) = (strPt[0], strPt[1])
(eX, eY) = (endPt[0], endPt[1])
xSpace = geom[1]
ySpace = geom[5]

sXIndex = ((sX - x0) / xSpace)
sYIndex = ((sY - y0) / ySpace)
eXIndex = ((eX - sXIndex) / xSpace) + sXIndex
eYIndex = ((eY - sYIndex) / ySpace) + sYIndex


dx = (eXIndex - sXIndex)
dy = (eYIndex - sYIndex)
xHeading = 1.0 if dx > 0 else -1.0 if dx < 0 else 0.0
yHeading = 1.0 if dy > 0 else -1.0 if dy < 0 else 0.0

xOffset = (1 - (math.modf(sXIndex)[0]))
yOffset = (1 - (math.modf(sYIndex)[0]))

ptsIndexes = []
x = sXIndex
y = sYIndex
pt = (x, y) #1st pt

if dx != 0:
    m = (float(dy) / float(dx))
    b = float(sY - sX * m )

dx = abs(int(dx))
dy = abs(int(dy))

if dx == 0:
    for h in range(0, dy + 1):
        pt = (x, y + (yHeading *h))
        ptsIndexes.append(pt)

    return ptsIndexes

#print("m {}, dx {}, dy {}, b {}, xdir {}, ydir {}".format(m, dx, dy, b, xHeading, yHeading))
#print("x {}, y {}, {} {}".format(sXIndex, sYIndex, eXIndex, eYIndex))

#snap to half a cell size so we can find intersections on cell boundaries
sXIdxSp = round(2.0 * sXIndex) / 2.0
sYIdxSp = round(2.0 * sYIndex) / 2.0
eXIdxSp = round(2.0 * eXIndex) / 2.0
eYIdxSp = round(2.0 * eYIndex) / 2.0
# ptsIndexes.append(pt)
prevPt = False
#advance half grid size
for w in range(0, dx * 4):
    x = xHeading * (w / 2.0) + sXIdxSp
    y = (x * m + b)
    if xHeading < 0:
        if x < eXIdxSp:
            break
    else:
        if x > eXIdxSp:
            break

    pt = (round(x), round(y)) #snapToGrid
    # print(w, x, y)

    if prevPt != pt:
        ptsIndexes.append(pt)
        prevPt = pt
#advance half grid size
for h in range(0, dy * 4):
    y = yHeading * (h / 2.0) + sYIdxSp
    x = ((y - b) / m)
    if yHeading < 0:
        if y < eYIdxSp:
            break
    else:
        if y > eYIdxSp:
            break
    pt = (round(x), round(y)) # snapToGrid
    # print(h, x, y)

    if prevPt != pt:
        ptsIndexes.append(pt)
        prevPt = pt

return set(ptsIndexes) #elminate duplicates