Help me with some ray-tracing math

Questions about the LÖVE API, installing LÖVE and other support related questions go here.
Forum rules
Before you make a thread asking for help, read this.
User avatar
kikito
Inner party member
Posts: 3153
Joined: Sat Oct 03, 2009 5:22 pm
Location: Madrid, Spain
Contact:

Help me with some ray-tracing math

Post by kikito »

Hi there,

I am trying to implement ray tracing in bump.lua

After looking for a while, I found this paper, which describes an algorithm which seems very appropiate for my needs. Unfortunately I don't understand this part:
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, store in tDeltaY the amount of movement along the ray which has a vertical component equal to the height of a voxel.
Can any of you math-magicians give it a look and tell me how am I supposed to initialize tMaxX, tMaxY, tDeltaX and tDeltaY? Math expression would be fine, Lua code would be perfect.

Also, I'm posting this in stackoverflow, if you want some rep points:

http://stackoverflow.com/questions/1236 ... hm-for-ray
When I write def I mean function.
User avatar
dreadkillz
Party member
Posts: 223
Joined: Sun Mar 04, 2012 2:04 pm
Location: USA

Re: Help me with some ray-tracing math

Post by dreadkillz »

Hello kikito. I recently implemented this in my sweep and prune when I learned how to raycast.

The paper is a bit hard to understand so I'll try to explain it another way. I'm not good at explaining things so bare with me.

Imagine a ray traversing a grid. the ray initially starts in a cell within the grid. What happens initially is that the ray touches either a vertical line or a horizontal line in the first cell. The point at which the ray hits this first line is the starting point of your NEXT cell (the cell your ray starts in need to be checked too). You'll first need to figure this first hit.

First you'll need to rasterize your ray's origin (x,y) into grid coordinates. In my case I just use:

cell_x = math.floor(x/cell_size)
cell_y = math.floor(y/cell_size)

The next part is to figure out which direction your ray is moving on either the x axis or y axis. If your ray is moving to the right on the x axis, you would add one to the x component of your cell coordinate b/c your first vertical line hit is to the right of you. If your ray is moving downward, you would add one to the y component of your cell coordinate cause first horizontal line hit is below you.

Next you'll have to figure out the smallest distance ratio of your ray to figure out which line it starts at. So let's say you describe your ray as x+dx,y+dy where dx and dy are the length of the ray's component on each axis. You'll need to figure out the distance from your first line hit to your starting point on each axis.

To do this, you'll just turn the coordinate of your first line hit into "real" coordinates and subtract your starting point. So if my ray is moving right and downward:

xd = ((cell_x+1) * cell_size) - x
yd = ((cell_y+1) * cell_size) - y

Now what happens next is you divide the distance to your first lines by the maximum lengths the ray travels on either axis to figure out the smallest distance ratio. The logic behind this is that the first line you hit (vertical or horizontal) is controlled by the shortest distance the ray has to travel.

So:
x_ratio = xd/dx
y_ratio = yd/dy

Say that the y distance ratio of the ray is greater than x distance ratio; if we multiply the ray's component by the y ratio (x+y_ratio*dx) (y+y_ratio*dy), we will miss our next horizontal line because we traveled too far to hit the next vertical line.

We will need to check the ratio's in a loop whenever the ray is moving onto the next line to ensure that it takes the shortest distance possible.

Next we need to calculate the ratio needed on each axis to reach the next vertical or horizontal line. Just divide the cell size on both axis by its total lengths
So:

x_delta_ratio = cell_size/dx
y_delta_ratio = cell_size/dy

Explaination:
So let's say your cell size is 5 and your ray's dx is 10 --> (5/10) = 1/2 . What this means is that you'll need half the x length of your ray to hit the next vertical line.

So you have your initial ratio's and you know which line that ray hits first. To figure out the next line and cell to check, Add your delta ratio to your smallest ratio (because smallest ratio dictates length of ray to the next cell) and increase the cell step.

So:

Say the first cell that you hit is (0,0). This doesn't include the cell that you start in.
The x ratio and y ratio of this first cell are 0.1,0.2 respectively.
The x delta ratio and y delta ratio are 0.5,0.4 respectively.

The smallest ratio is the x ratio so we add x delta to it:
0.1 + 0.5 = 0.6 --> new x ratio

We increase the cell coordinate by -1 or 1 depending on our direction of the ray. The next cell to check is either (-1,0) or (1,0).

The next loop our y ratio is now smaller (0.2). We increase it and the cell stepping:
new y ratio = (0.2+0.4) = 0.6. The next cell is (-1,1) or (1,1)

Now we have equal ratio's (0.6 and 0.6). This means our next cell is next diagonal cell. Just step on either axis as we will reach our desire cells So:

(0,0) cell and next is (1,1). Increase x first ---> (1,0), then next loop increase y (1,1).

Here's the relevant code:

Code: Select all

local raycast = function(self,x,y,dx,dy,isCoroutine)
	local set   = {}
	local x0,y0 = floor(x/self.width),floor(y/self.height)
	
	local dxRatio,dyRatio,xDelta,yDelta,xStep,yStep,smallest,xStart,yStart
	-- visualization of sides that we can check: [0 1]
	-- if we're moving to the right, +1 to cell-x to
	-- check the right side first as it's > our starting point
	-- positive steps when moving right, neg steps when moving left
	if dx > 0 then 
		xStep   = 1 
		xStart  = 1 
	else 
		xStep   = -1 
		xStart  = 0
	end
	if dy > 0 then 
		yStep   = 1 
		yStart  = y0 + 1
	else 
		yStep   = -1 
		yStart  = 0
	end
	
	-- dxRatio: (x0*width-x)/dx precalculations
	-- xDelta is the ratio of the ray's width to reach the next vertical line on the x-axis
	-- always take the shortest ratio to reach the nearest line on the grid
	-- zero hack
	if dx == 0 then
		dxRatio = math.huge
		xDelta  = 0
	else
		local a,b = self.width/dx,x/dx
		dxRatio   = a*(x0+xStart)-b
		xDelta    = a*xStep
	end
	if dy == 0 then
		dyRatio = math.huge
		yDelta  = 0
	else
		local a,b = self.height/dy,y/dy
		dyRatio   = a*(y0+yStart)-b
		yDelta    = a*yStep
	end
	
	-- Use a repeat loop so that the ray checks its starting cell first before moving.
	repeat
		local row = self.cells[x0]
		if row and row[y0] then
			-- if called as an iterator, iterate through all objects in the cell and the next cells
			-- otherwise, just look for the earliest hit and return
			if isCoroutine then
				for obj,hitx,hity in row[y0]:iterRay(x,y,dx,dy) do
						if not set[obj] then coroutine.yield(obj,hitx,hity); set[obj]=true end
				end
			else
				local obj,hitx,hity = row[y0]:rayQuery(x,y,dx,dy)
				if obj then return obj,hitx,hity end
			end
		end
		
		if dxRatio < dyRatio then
			smallest  = dxRatio
			dxRatio   = dxRatio + xDelta
			x0        = x0 + xStep
		else
			smallest  = dyRatio
			dyRatio   = dyRatio + yDelta
			y0        = y0 + yStep
		end
	until smallest > 1
end
User avatar
vrld
Party member
Posts: 917
Joined: Sun Apr 04, 2010 9:14 pm
Location: Germany
Contact:

Re: Help me with some ray-tracing math

Post by vrld »

First, some definitions they skipped in the paper (note: I write vectors as bold letters instead of \(\vec{l}etters\) with arrows on top):
A ray is a linear, vectorial function \(\mathbf{r}(t) = \mathbf{u} + \mathbf{v}\cdot t\) of some parameter \(t >= 0\) that returns a position \(\mathbf{x} = \mathbf{r}(t)\). You can interpret it as following: (In space) starting at point u, you throw some object in the direction v. The ray function returns the position the object is at time t.

tMaxX and tMaxY
... so, informal, you need to compute the times when the ray hits the first vertical and horizontal grid boundaries.

Formally, lets assume the grid has the origin \(\mathbf{g} = (g_x, g_y)\) (more often than not \(\mathbf{g} = (0,0)\)) and the grid size \(\mathbf{s} = (s_x, s_y)\), i.e. each cell has the width \(s_x\) and height \(s_y\).

With \(i\) being an integer, the vertical lines have an x-coordinate of \(x(i) = g_x + i \cdot s_x\). Likewise, for the horizontal lines \(y(k) = g_y + k\cdot s_y\). The \(i\) and \(k\) also uniquely denote each cell of the grid: The cell index \((i,k)\) has the same parameters as the left and upper bounding line. This is basically just a complicated way of describing the same thing you do with tile-based games when you calculate the tile coordinates for a given position.

You already know the cell where the ray starts. Lets denote this cell with \((i_\mathbf{u}, k_\mathbf{u})\).

The time \(tMaxX\) when the ray hits the first vertical line can be computed by solving \(u_x + v_x \cdot t = g_x + (i_\mathbf{u} + 1)\cdot s_x\) for \(t\),
\[tMaxX = \frac{g_x + (i_\mathbf{u} + 1)\cdot s_x - u_x}{v_x}.\]
The same has to be done to find \(tMaxY\).


tDeltaX and tDeltaY
These sound more complicated than they are. All they really ask is: How big is \(t\) for \(u_x + v_x \cdot t = u_x + s_x\) (and likewise for y), i.e.: \[tDeltaX = \frac{s_x}{v_x}.\]

Edit: Ninja'd by dreadkillz :shock:
Last edited by vrld on Tue Sep 11, 2012 12:32 pm, edited 1 time in total.
I have come here to chew bubblegum and kick ass... and I'm all out of bubblegum.

hump | HC | SUIT | moonshine
User avatar
dreadkillz
Party member
Posts: 223
Joined: Sun Mar 04, 2012 2:04 pm
Location: USA

Re: Help me with some ray-tracing math

Post by dreadkillz »

Yea I need to learn Latex so I can write it out all pretty with vector notations. That was painful for me (hopefully not for you to read) to type.
User avatar
kikito
Inner party member
Posts: 3153
Joined: Sat Oct 03, 2009 5:22 pm
Location: Madrid, Spain
Contact:

Re: Help me with some ray-tracing math

Post by kikito »

self.mind:blow() :D

EDIT:

I'm trying vrld's solution but there is something I don't understand.

It seems to me that tmaxX and tdeltaX (and their Y counterparts) can be either positive or negative, since \(v_x\) and \(v_y\) can be, and they are used in the calculations.

If that's the case, then the iterative part can not work:

Code: Select all

loop {
  if(tMaxX < tMaxY) {
    tMaxX= tMaxX + tDeltaX;
    X= X + stepX;
  } else {
    tMaxY= tMaxY + tDeltaY;
    Y= Y + stepY;
  }
  NextVoxel(X,Y);
}
Notice how it says tMaxX < tMaxY . Unless they are both positive numbers, this conditional surely will fail to do what it seems to be intending to do (in fact, right now I have some very nice infinite loops).
Must I understand that the paper is only dealing with the first quadrant on this part of the code? Or are tMaxX and tDeltaX supposed to be positive, and I'm missing some math.abs somewhere?
When I write def I mean function.
User avatar
vrld
Party member
Posts: 917
Joined: Sun Apr 04, 2010 9:14 pm
Location: Germany
Contact:

Re: Help me with some ray-tracing math

Post by vrld »

You're absolutely right. When calculating the initial \(tMaxX\), you also have to take the direction of the ray into account - when the ray travels to the left, (\(v_x < 0\)), the first vertical cell boundary is the one of cell \((i_u,\cdot)\), not of cell \((i_u+1,\cdot)\). In this case, the equation to solve becomes \[tMaxX = \frac{g_x + i_u \cdot s_x - u_x}{v_x}.\]
\(tMaxX\) will be positive: If \(v_x < 0\), i.e. the ray travels to the left, then \(g_x + i_u\cdot s_x \le u_x\) and thus \(g_x + i_u\cdot s_x - u_x \le 0\).

As you said, to make the iteration work \(tDeltaX\) has to be positive. The formulation in the paper ("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") is a bit unfortunate, but you have to in fact calculate \[tDeltaX = \left|\frac{s_x}{v_x}\right|.\]
The value is still "how far we must move for the horizontal component of such a movement to equal the width of a voxel" ;)
I have come here to chew bubblegum and kick ass... and I'm all out of bubblegum.

hump | HC | SUIT | moonshine
User avatar
dreadkillz
Party member
Posts: 223
Joined: Sun Mar 04, 2012 2:04 pm
Location: USA

Re: Help me with some ray-tracing math

Post by dreadkillz »

In addition to the sign change for tDeltaX and tMaxX, you'll also need to check for the case when 0/0 (indeterminate). In this case, just convert tMaxX or tMaxY to infinity (math.huge).
User avatar
kikito
Inner party member
Posts: 3153
Joined: Sat Oct 03, 2009 5:22 pm
Location: Madrid, Spain
Contact:

Re: Help me with some ray-tracing math

Post by kikito »

Hi, I'm happy to report that I got the broad phase working, thanks to vrld's indications. Here's the releavant code:

Code: Select all


local abs, floor = math.abs, math.floor

local __cellSize = 128

local function _toGrid(wx, wy)
  return floor(wx / __cellSize), floor(wy / __cellSize)
end

local function sign(n)
  return n < 0 and -1 or (n > 0  and 1 or 0)
end

function printSegment(x1,y1,x2,y2)

  local vx, vy = x2 - x1, y2 - y1
  local x, y = _toGrid(x1,y1)
  local ox, oy = _toGrid(x2,y2)

  print('==================')
  print(x1,y1, "     =>     ", x2, y2)
  print(x,y, "     =>     ", ox, oy)
  print(vx, vy)
  print('--------')


  local stepX, stepY   = sign(vx), sign(vy)
  local deltaX, deltaY = abs(__cellSize / vx), abs(__cellSize / vy)

  local maxX = ((x + (vx > 0 and 1 or 0)) * __cellSize - x1) / vx
  local maxY = ((y + (vy > 0 and 1 or 0)) * __cellSize - y1) / vy

  print(x,y,maxX,maxY)

  while x~=ox or y~=oy do

    if maxX < maxY then
      maxX = maxX + deltaX
      x = x + stepX
    else
      maxY = maxY + deltaY
      y = y + stepY
    end

    print(x, y, maxX, maxY)

  end
  print('==================')
end


printSegment(515,515,76,263)
Interestingly enough, I didn't need to do anything special to treat the vx = 0 or vy = 0 cases. Apparently dividing by 0 returns something very similar (I end up getting nan instead of inf in maxX or maxY, but they work the same way in the conditions I need).

Now that I have (with lots of unestimable help) finished the broad phase, it's time to program the narrow one: colliding boxes with a segment. I'll let you know if I need help with that one.
Last edited by kikito on Wed Nov 14, 2012 11:01 pm, edited 1 time in total.
When I write def I mean function.
User avatar
dreadkillz
Party member
Posts: 223
Joined: Sun Mar 04, 2012 2:04 pm
Location: USA

Re: Help me with some ray-tracing math

Post by dreadkillz »

Very nice. Looking forward to seeing it in action.
User avatar
kikito
Inner party member
Posts: 3153
Joined: Sat Oct 03, 2009 5:22 pm
Location: Madrid, Spain
Contact:

Re: Help me with some ray-tracing math

Post by kikito »

I'm afraid I have to resurrect this old thread.

I've been rewriting bump from scratch with tests and could not track this down until last night.

I did some interactive tests and realized that algorithm I developed with vrld's enters infinite loops at a seemingly random positions.

See attached .love:
gridSegment.love
This program has an infinite loop
(1007 Bytes) Downloaded 234 times
The program will draw a grid, a line going from 1 point to the mouse, and highlight the elements on the grid which "touch" the line. You may click around to change the origin of the line.

Keep doing it for a while and you will soon find that the program "blocks", and you will have to manually kill it :(

If you open it from the console you will see the lines it's calculating being printed out. The last one will be the one causing the hang.

The offending code must be somewhere in the gridSegment function:

Code: Select all

function gridSegment(cellSize, x1,y1, x2,y2)
  local x, y    = toGrid(cellSize, x1,y1)
  local ox, oy  = toGrid(cellSize, x2,y2)

  local vx, vy          = x2 - x1, y2 - y1
  local stepX, stepY    = sign(vx), sign(vy)
  local deltaX, deltaY  = abs(cellSize / vx), abs(cellSize / vy)

  local maxX = ((x + (vx > 0 and 0 or -1)) * cellSize - x1 + 1) / vx
  local maxY = ((y + (vy > 0 and 0 or -1)) * cellSize - y1 + 1) / vy

  local result, len = {x,y}, 2

  while x~=ox or y~=oy do
    if maxX < maxY then
      maxX = maxX + deltaX
      x = x + stepX
    else
      maxY = maxY + deltaY
      y = y + stepY
    end

    result[len + 1] = x
    result[len + 2] = y
    len = len + 2
  end

  return result, len
end
As a reference, this function enters an infinite loop when cellSize = 64, origin = 531,226 and destination = 464,226.

If anyone feels knowledgeable enough to trackle the problem, or knows a better knows a better I line-drawing algorithm, I'd be grateful. I'm open to suggestions.
When I write def I mean function.
Post Reply

Who is online

Users browsing this forum: No registered users and 203 guests