How do I initialize the t-variables in "A Fast Voxel Traversal Algorithm for Ray Tracing"?

Question

I am trying to implement the algorithm explained on this paper, used to traverse grid cells in order following a straight line, which is useful for ray tracing:

http://www.cse.yorku.ca/~amana/research/grid.pdf

The paper describes the algorithm as two parts: initialisation and iterative traversal. I can undersand the iterative traversal part, but I'm having trouble understanding how some of the variables in the initialisation part are calculated.

I need help initialising tMaxX, tMaxY, tDeltaX & tDeltaY. Their initialisation procedure is explained as follows:

Next, we determine the value of t at which the ray crosses the ﬁrst 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, store in tDeltaY the amount of movement along the ray which has a vertical component equal to the height of a voxel.

I'm not able to deduce the code I need form the English description given above. Can someone translate it to a math/pseudocode expression for me?

MBo · Accepted Answer · 2016-04-14T15:42:31.280

11

Initialization for X-coordinate variables (the same for Y)

  DX = X2 - X1
  tDeltaX = GridCellWidth / DX
  tMaxX = tDeltaX * (1.0 - Frac(X1 / GridCellWidth)) 
  //Frac if fractional part of float, for example, Frac(1.3) = 0.3, Frac(-1.7)=0.3

Example:

  GridCellWidth, Height = 20
  X1 = 5, X2 = 105 
  Y1 = 5, Y2 = 55 
  DX = 100, DY  = 50
  tDeltaX = 0.2, tDeltaY = 0.4 
  tMaxX = 0.2 * (1.0 - 0.25) = 0.15  //ray will meet first vertical line at this param
  tMaxY = 0.4 * (1.0 - 0.25) = 0.3   //ray will meet first horizontal line at this param

We can see that first cell border will be met at parameter t = 0.15

edited Apr 14 '16 at 15:42

answered Sep 11 '12 at 13:09

MBo

77,366
5
53
86

It must be noted, that this Frac() function must return a *positive* fraction for negative numbers in contrast to what is implemented in some standard libraries (which return a negative fraction for negative numbers). – josch Apr 14 '16 at 15:02
It doesn't work for me with negative values, example: x1 = 0, x2 = 0; x1 = 1, x2 = -1; width = 1 height = 1 dx = 1; dy = -1 tDeltaX = 1; tDeltaY = -1; tMaxX = 1 * (1 - 0) = 1 tMaxY = -1 * (1 - 0) = -1; In this case tMaxY will alway be smaller than tMaxX, and adding tDeltaY won't change much – SnowyCoder Mar 19 '17 at 10:32
@SnowyCoder Transform direction (using mirroring) to the first quadrant of coordinate plane. – MBo Mar 19 '17 at 12:12

score 10 · Answer 2 · edited Dec 27 '21 at 15:21

10

The one that worked for me in 3D for both positive and negative directions (in CUDA C):

#define SIGN(x) (x > 0 ? 1 : (x < 0 ? -1 : 0))
#define FRAC0(x) (x - floorf(x))
#define FRAC1(x) (1 - x + floorf(x))

float tMaxX, tMaxY, tMaxZ, tDeltaX, tDeltaY, tDeltaZ;
int3 voxel;

float x1, y1, z1; // start point   
float x2, y2, z2; // end point   

int dx = SIGN(x2 - x1);
if (dx != 0) tDeltaX = fmin(dx / (x2 - x1), 10000000.0f); else tDeltaX = 10000000.0f;
if (dx > 0) tMaxX = tDeltaX * FRAC1(x1); else tMaxX = tDeltaX * FRAC0(x1);
voxel.x = (int) x1;

int dy = SIGN(y2 - y1);
if (dy != 0) tDeltaY = fmin(dy / (y2 - y1), 10000000.0f); else tDeltaY = 10000000.0f;
if (dy > 0) tMaxY = tDeltaY * FRAC1(y1); else tMaxY = tDeltaY * FRAC0(y1);
voxel.y = (int) y1;

int dz = SIGN(z2 - z1);
if (dz != 0) tDeltaZ = fmin(dz / (z2 - z1), 10000000.0f); else tDeltaZ = 10000000.0f;
if (dz > 0) tMaxZ = tDeltaZ * FRAC1(z1); else tMaxZ = tDeltaZ * FRAC0(z1);
voxel.z = (int) z1;

while (true) {
    if (tMaxX < tMaxY) {
        if (tMaxX < tMaxZ) {
            voxel.x += dx;
            tMaxX += tDeltaX;
        } else {
            voxel.z += dz;
            tMaxZ += tDeltaZ;
        }
    } else {
        if (tMaxY < tMaxZ) {
            voxel.y += dy;
            tMaxY += tDeltaY;
        } else {
            voxel.z += dz;
            tMaxZ += tDeltaZ;
        }
    }
    if (tMaxX > 1 && tMaxY > 1 && tMaxZ > 1) break;
    // process voxel here
}

Note, grid cell's width/height/depth are all equal to 1 in my case, but you can easily modify it for your needs.

edited Dec 27 '21 at 15:21

ProjectPhysX

4,535
2
14
34

answered Jul 24 '16 at 13:42

a5kin

1,335
16
20

Excellent algorithm/code. Though some suggestions to make it flawless: 1. Add process voxel right after loop statement. 2. Add process voxel before breaking out to process the last voxel. These two changes make sure that the first and last voxel are processed as well. – oOo Mar 12 '22 at 06:36
It may have some slight floating point instabilities at end of ray, exit of grid, could use some more work there, maybe subtract tiny value of tMaxX and tMaxY in comparison – oOo Mar 17 '22 at 22:36
It has ending condition issues, can go slightly outside. – oOo Mar 17 '22 at 22:56
It's paradox... must be solved specially, can't have it both ways. – oOo Mar 17 '22 at 23:05
Intersection either belongs to begin or end, can't have it both ways, most be solved specially. Possibly be ending loop and process on extra voxel for end cell. – oOo Mar 17 '22 at 23:06
Works in 2D too , just drop the Z stuff. I found I needed to add `process voxel` for the actual start and end points. – c z Apr 24 '23 at 13:24

How do I initialize the t-variables in "A Fast Voxel Traversal Algorithm for Ray Tracing"?

2 Answers2

Linked