How do I generalise Bresenham's line algorithm to floating-point endpoints?

Question

I'm trying to combine two things. I'm writing a game and I need to determine the grid squares lying on a line with the floating-point endpoints.

Moreover I need it to include all the grid squares it touches (i.e. not just Bresenham's line but the blue one):

Can someone offer me any insight on how to do it? The obvious solution is to use naive line algorithm, but is there something more optimized (faster)?

Answer 1

You are looking for a grid traversal algorithm. This paper gives a good implementation;

Here's the basic implementation in 2D found on the paper:

loop {
    if (tMaxX < tMaxY) {
        tMaxX = tMaxX + tDeltaX;
        X = X + stepX;
    } else {
        tMaxY = tMaxY + tDeltaY;
        Y = Y + stepY;
    }
    NextVoxel(X, Y);
}

There's also a 3D ray-casting version on the paper.

In case the link rots, you can find many mirrors with its name: A faster voxel traversal algorithm for raytracing.

---

Relevant bits from the paper:

...it is a variant of the DDA line algorithm. However, instead of basing it on the simple DDA..., in which an unconditional step along one axis is required, ours has no preferred axis. This considerably simplifies the inner loop and allows for easy testing of an intersection point to see if it is in the current voxel.

...

The new traversal algorithm breaks down the ray into intervals of t, each of which spans one voxel. We start at the ray origin and visit each of these voxels in interval order. The traversal algorithm consists of two phases: initialization and incremental traversal.

...

The integer variables X and Y are initialized to the starting voxel coordinates. In addition, the variables stepX and stepY are initialized to either 1 or -1 indicating whether X and Y are incremented or decremented as the ray crosses voxel boundaries (this is determined by the sign of the x and y components...)

...

Next, we determine the value of t at which the ray crosses the first vertical voxel boundary and store it in variable tMaxX. We perform a similar computation in y and store the result in tMaxY. The minimum of these two values will indicate how much we can travel along the ray and still remain in the current voxel.

Finally, we compute tDeltaX and tDeltaY. tDeltaX indicates how far along the ray we must move (in units of t) for the horizontal component of such a movement to equal the width of a voxel. Similarly, we store in tDeltaY the amount of movement along the ray which has a vertical component equal to the height of a voxel.

Answer 2

Blue's idea is good, but the implementation is a bit clumsy. In fact, you can easily do it without sqrt. Let's assume for the moment that you exclude degenerate cases (BeginX==EndX || BeginY==EndY) and focus only on line directions in the first quadrant, so BeginX < EndX && BeginY < EndY. You'll have to implement a version for at least one other quadrant too, but that's very similar to the version for the first quadrant - you only check other edges. In C'ish pseudo code:

int cx = floor(BeginX); // Begin/current cell coords
int cy = floor(BeginY);
int ex = floor(EndX); // End cell coords
int ey = floor(EndY);

// Delta or direction
double dx = EndX-BeginX;
double dy = EndY-BeginY;

while (cx < ex && cy < ey)
{
  // find intersection "time" in x dir
  float t0 = (ceil(BeginX)-BeginX)/dx;
  float t1 = (ceil(BeginY)-BeginY)/dy;

  visit_cell(cx, cy);

  if (t0 < t1) // cross x boundary first=?
  {
    ++cx;
    BeginX += t0*dx;
    BeginY += t0*dy;
  }
  else
  {
    ++cy;
    BeginX += t1*dx;
    BeginY += t1*dy;
  }
}

Now for other quadrants, you just change the ++cx or ++cy and the loop condition. If you use this for collision, you probably have to implement all 4 versions, otherwise you can get away with two by appropriately swapping begin and end points.

Answer 3

float difX = end.x - start.x;
float difY = end.y - start.y;
float dist = abs(difX) + abs(difY);

float dx = difX / dist;
float dy = difY / dist;

for (int i = 0, int x, int y; i <= ceil(dist); i++) {
    x = floor(start.x + dx * i);
    y = floor(start.y + dy * i);
    draw(x,y);
}
return true;

JS Demo :

//C# stackoverflow

function verifyLoS(start, end) {
    var difX = end.x - start.x;
    var difY = end.y - start.y;
	var dist = Math.abs(difX) + Math.abs(difY);

    var dx = difX / dist;
    var dy = difY / dist;

    for (var i = 0, x, y; i <= Math.ceil(dist); i++) {
        x = Math.floor(start.x + dx * i);
        y = Math.floor(start.y + dy * i);
        draw(x,y);
    }
    return true;
}

var HEIGHT = 14,
    WIDTH = 14,
    PX = 32;

var canvas, ctx;

function main () {
  canvas = document.getElementById("canvas");
  ctx = canvas.getContext("2d");
  
  canvas.height = HEIGHT * PX;
  canvas.width = WIDTH * PX;

  loop();
}

function loop() {
  for (var x=0;x<HEIGHT; x++) {
    for (var y=0; y<WIDTH; y++) {
      ctx.fillStyle = ((x+y)%2===0)?"#AFA":"#AEA";
      ctx.fillRect(x*PX,y*PX,PX,PX);
    }
  }
  var a = {x:1.5, y:9.5},
      b = {x:12.5, y:1.5};
  verifyLoS(a, b);
  ctx.beginPath();
  ctx.moveTo(a.x*PX,a.y*PX);
  ctx.lineTo(b.x*PX,b.y*PX);
  ctx.stroke();
  //window.requestAnimationFrame(loop);
}

function draw(x, y) {
    ctx.fillStyle = "rgba(255,0,0,0.2)";
    ctx.fillRect(x*PX,y*PX,PX,PX);
}

main();

<!DOCTYPE html>
<html>
<head>
  <meta charset="utf-8">
  <title>JS Bin</title>
</head>
<body>
  <canvas id="canvas"></canvas>
</body>
</html>

Answer 4

Your assumption isn't necessarily to find the cells but the lines it cross on this grid.

For example taking your image we can highlight not the cells, but lines of the grid it crosses:

This then shows that if it crosses a grid line that the cells either side of this line are those that are filled.

You can use an intersection algorithm to find if your floating point line will cross these by scaling your points to pixels. If you have a 1.0:1 ratio of floating coordinates:pixels then you're sorted and you can just translate it directly. Using the Line segment intersection algorithm you can check if your lower left line (1,7)(2,7) intersects with your line (1.3,6.2)(6.51,2.9). http://alienryderflex.com/intersect/

Some translation from c to C# will be needed but you can get the idea from that paper. I'll put the code below in-case the link breaks.

//  public domain function by Darel Rex Finley, 2006

//  Determines the intersection point of the line defined by points A and B with the
//  line defined by points C and D.
//
//  Returns YES if the intersection point was found, and stores that point in X,Y.
//  Returns NO if there is no determinable intersection point, in which case X,Y will
//  be unmodified.

bool lineIntersection(
double Ax, double Ay,
double Bx, double By,
double Cx, double Cy,
double Dx, double Dy,
double *X, double *Y) {

  double  distAB, theCos, theSin, newX, ABpos ;

  //  Fail if either line is undefined.
  if (Ax==Bx && Ay==By || Cx==Dx && Cy==Dy) return NO;

  //  (1) Translate the system so that point A is on the origin.
  Bx-=Ax; By-=Ay;
  Cx-=Ax; Cy-=Ay;
  Dx-=Ax; Dy-=Ay;

  //  Discover the length of segment A-B.
  distAB=sqrt(Bx*Bx+By*By);

  //  (2) Rotate the system so that point B is on the positive X axis.
  theCos=Bx/distAB;
  theSin=By/distAB;
  newX=Cx*theCos+Cy*theSin;
  Cy  =Cy*theCos-Cx*theSin; Cx=newX;
  newX=Dx*theCos+Dy*theSin;
  Dy  =Dy*theCos-Dx*theSin; Dx=newX;

  //  Fail if the lines are parallel.
  if (Cy==Dy) return NO;

  //  (3) Discover the position of the intersection point along line A-B.
  ABpos=Dx+(Cx-Dx)*Dy/(Dy-Cy);

  //  (4) Apply the discovered position to line A-B in the original coordinate system.
  *X=Ax+ABpos*theCos;
  *Y=Ay+ABpos*theSin;

  //  Success.
  return YES; }

If you need to find out only when (and where) the line segments intersect, you can modify the function as follows:

//  public domain function by Darel Rex Finley, 2006  

//  Determines the intersection point of the line segment defined by points A and B
//  with the line segment defined by points C and D.
//
//  Returns YES if the intersection point was found, and stores that point in X,Y.
//  Returns NO if there is no determinable intersection point, in which case X,Y will
//  be unmodified.

bool lineSegmentIntersection(
double Ax, double Ay,
double Bx, double By,
double Cx, double Cy,
double Dx, double Dy,
double *X, double *Y) {

  double  distAB, theCos, theSin, newX, ABpos ;

  //  Fail if either line segment is zero-length.
  if (Ax==Bx && Ay==By || Cx==Dx && Cy==Dy) return NO;

  //  Fail if the segments share an end-point.
  if (Ax==Cx && Ay==Cy || Bx==Cx && By==Cy
  ||  Ax==Dx && Ay==Dy || Bx==Dx && By==Dy) {
    return NO; }

  //  (1) Translate the system so that point A is on the origin.
  Bx-=Ax; By-=Ay;
  Cx-=Ax; Cy-=Ay;
  Dx-=Ax; Dy-=Ay;

  //  Discover the length of segment A-B.
  distAB=sqrt(Bx*Bx+By*By);

  //  (2) Rotate the system so that point B is on the positive X axis.
  theCos=Bx/distAB;
  theSin=By/distAB;
  newX=Cx*theCos+Cy*theSin;
  Cy  =Cy*theCos-Cx*theSin; Cx=newX;
  newX=Dx*theCos+Dy*theSin;
  Dy  =Dy*theCos-Dx*theSin; Dx=newX;

  //  Fail if segment C-D doesn't cross line A-B.
  if (Cy<0. && Dy<0. || Cy>=0. && Dy>=0.) return NO;

  //  (3) Discover the position of the intersection point along line A-B.
  ABpos=Dx+(Cx-Dx)*Dy/(Dy-Cy);

  //  Fail if segment C-D crosses line A-B outside of segment A-B.
  if (ABpos<0. || ABpos>distAB) return NO;

  //  (4) Apply the discovered position to line A-B in the original coordinate system.
  *X=Ax+ABpos*theCos;
  *Y=Ay+ABpos*theSin;

  //  Success.
  return YES; }

Answer 5

I ran into the same problem today and made a pretty big mountain of spaghetti out of a mole hill but ended up with something that works: https://github.com/SnpM/Pan-Line-Algorithm.

From the ReadMe:

The core concept of this algorithm is similar to Bresenham's in that it increments by 1 unit on one axis and tests the increase on the other axis. Fractions make incrementing considerably more difficult, however, and a lot of pizzas had to be added in. For example, incrementing from X = .21 to X = 1.21 with a slope of 5 makes for a complex problem (coordinate patterns between those nasty numbers are hard to predict) but incrementing from 1 to 2 with a slope of 5 makes for an easy problem. Coordinate pattern between integers is very easy to solve (just a line perpendicular to the incrementing axis). To get the easy problem, the incrementing is offset to an integer number with all the calculations done seperately for the fractional piece. So instead of starting the incrementing on .21, the algorithm tries to start the incrementing on 1 and goes from there.

The ReadMe explains the solution much better than the code. I'm planning on revising it to be less headache-inducing.

I know I'm about a year late to this question but I hope this gets to others who are looking for a solution to this problem.

Answer 6

Got it working perfectly! I combined Gustavo Maciel's answer with A-312's (using algorithm from paper Gustavo cited). Thanks guys!

Solution (see intersect() function and comments): https://jsfiddle.net/6x7t4q1o/5

function intersect(start, end) {
  //Grid cells are 1.0 X 1.0.
  let x = Math.floor(start.x);
  let y = Math.floor(start.y);
  let diffX = end.x - start.x;
  let diffY = end.y - start.y;
  let stepX = Math.sign(diffX);
  let stepY = Math.sign(diffY);
  
  //Ray/Slope related maths.
  //Straight distance to the first vertical grid boundary.
  let xOffset = end.x > start.x ?
  	(Math.ceil(start.x) - start.x) :
    (start.x - Math.floor(start.x));
  //Straight distance to the first horizontal grid boundary.
  let yOffset = end.y > start.y ?
  	(Math.ceil(start.y) - start.y) :
    (start.y - Math.floor(start.y));
  //Angle of ray/slope.
  let angle = Math.atan2(-diffY, diffX);
  //NOTE: These can be divide by 0's, but JS just yields Infinity! :)
  //How far to move along the ray to cross the first vertical grid cell boundary.
  let tMaxX = xOffset / Math.cos(angle);
  //How far to move along the ray to cross the first horizontal grid cell boundary.
  let tMaxY = yOffset / Math.sin(angle);
  //How far to move along the ray to move horizontally 1 grid cell.
  let tDeltaX = 1.0 / Math.cos(angle);
  //How far to move along the ray to move vertically 1 grid cell.
  let tDeltaY = 1.0 / Math.sin(angle);
  
  //Travel one grid cell at a time.
  let manhattanDistance = Math.abs(Math.floor(end.x) - Math.floor(start.x)) +
  	Math.abs(Math.floor(end.y) - Math.floor(start.y));
  for (let t = 0; t <= manhattanDistance; ++t) {
    drawSquare(x, y);
    //Only move in either X or Y coordinates, not both.
  	if (Math.abs(tMaxX) < Math.abs(tMaxY)) {
    	tMaxX += tDeltaX;
      x += stepX;
    } else {
    	tMaxY += tDeltaY;
      y += stepY;
    }
  }
  
  infos.innerHTML = "";
  infos.innerHTML += "diffX:   " + diffX   + "<br/>diffY:   " + diffY + "<br/>";
  infos.innerHTML += "stepX:   " + stepX   + "<br/>stepY:   " + stepY + "<br/>";
  infos.innerHTML += "xOffset: " + xOffset + "<br/>yOffset: " + yOffset + "<br/>";
  infos.innerHTML += "tMaxX:   " + tMaxX   + "<br/>tMaxY:   " + tMaxY + "<br/>";
  infos.innerHTML += "tDeltaX: " + tDeltaX + "<br/>tDeltaY: " + tDeltaY + "<br/><br/>";
  infos.innerHTML += "angle: " + angle + "<br/>";
  infos.innerHTML += "manhattanDistance: " + manhattanDistance + "<br/>";
}

let WIDTH = 15,
    HEIGHT = 15,
    PX = 32;
let animationRate = 100;

let canvas, ctx;
let lastAnimationTime = -animationRate;
let animations = 0;

//Line coordinates.
let a = { x: WIDTH / 3,     y: HEIGHT * 2 / 3 };
let b = { x: WIDTH * 2 / 3, y: HEIGHT / 3 };
let resetA = { x: a.x, y: a.y };
let resetB = { x: b.x, y: b.y };

let infos;

function main () {
  canvas = document.getElementById("canvas");
  ctx = canvas.getContext("2d");
  canvas.height = HEIGHT * PX;
  canvas.width = WIDTH * PX;
  ctx.lineWidth = 1.5;
  
  infos = document.getElementById("infos");
  
  requestAnimationFrame(animate);
}

function animate(elapsedTime) {
	if (elapsedTime - lastAnimationTime > animationRate) {
  	lastAnimationTime += animationRate;
    ++animations;
    loop(elapsedTime);
  }
  if (animations > 500) {
  	a = { x: resetA.x, y: resetA. y };
    b = { x: resetB.x, y: resetB. y };
    animations = 0;
  }
	requestAnimationFrame(animate);
}

function loop(elapsedTime) {
  //Draw checkerboard.
  for (let x = 0; x < HEIGHT; x++) {
    for (let y = 0; y < WIDTH; y++) {
      ctx.fillStyle = ((x + y) % 2) ? "#90FF90": "#90F090";
      ctx.fillRect(x * PX, y * PX, PX, PX);
    }
  }
  
  let radius = Math.sin(elapsedTime * 0.0001) * 5.0;
  a.x = Math.cos(elapsedTime * 0.0002) * -radius + WIDTH * 0.5;
  a.y = Math.sin(elapsedTime * 0.0002) * radius + HEIGHT * 0.5;
  b.x = Math.cos(elapsedTime * 0.0002) * radius + WIDTH * 0.5;
  b.y = Math.sin(elapsedTime * 0.0002) * -radius + HEIGHT * 0.5;
  
  //Calculate and draw intersection squares.
  intersect(a, b);
  
  //Draw line.
  ctx.beginPath();
  ctx.moveTo(a.x * PX, a.y * PX);
  ctx.lineTo(b.x * PX, b.y * PX);
  ctx.stroke();
  
  //Draw S and E.
  ctx.fillStyle = "#00000050";
  ctx.fillText("S",
  	a.x * PX + 50.0 * Math.sign(a.x - b.x),
  	a.y * PX + 50.0 * Math.sign(a.y - b.y));
  ctx.fillText("E",
  	b.x * PX + 50.0 * Math.sign(b.x - a.x),
  	b.y * PX + 50.0 * Math.sign(b.y - a.y));
}

function drawSquare(x, y) {
  ctx.fillStyle = "rgba(255, 0, 0, 0.2)";
  ctx.fillRect(x * PX, y * PX, PX, PX);
}

main();

<!DOCTYPE html>
<html>
<head>
  <meta charset="utf-8">
  <title>JS Bin</title>
</head>
<body style="background-color: white;">
  <canvas id="canvas"></canvas>
  <pre id="infos"></pre>
</body>
</html>