Raytracing: Theory & Implementation

Raytracing: Theory & Implementation:

Part 1, Introduction
Part 2, Phong, Mirrors and Shadows
Part 3, Refractions and Beer's Law
Part 4, Spatial Subdivisions
Part 5, Soft Shadows
Part 6, Textures, Cameras and Speed
Part 7, Kd-Trees and More Speed

Author: Jacco Bikker

Note: These articles have been previously published at flipCode. They have been re-published at DevMaster.net based on author's request.

Very nice, I'm glad to hear that Jacco is still alive and that he knows devmaster.net

just wanted to say thanks for making this available on devmaster.net

No the honour is mine, flipcode is hardly visible at the moment and this is a great new spot.

Happy New Year to everyone, by the way. Hope you had some time to code. I sure did.

I'm really happy to see articles posted on here, regardless of whether or not they're new.

I hope this can help rekindle the small community feeling the flipCode once had but with its own personality. I am surprised to see so many old flipCode members here that I wrote off as being gone forever. Good to see you guys around again.

Originally Posted by phantom

Happy New Year to everyone, by the way.

Happy new year to you, comrade.

Ah, tbp! Long time no see! Did you notice that the Intel paper was finally published?

For the unaware: I worked together with an intel dude to write a paper on the SSE/packet tracing details of the latest version of the ray tracer. This cooperation went, well, 'suboptimal' but at least the paper is online:

http://www.intel.com/cd/ids/develope...ing/245711.htm

Enjoy. It's basically 'part 8' of the tutorial. It doesn't include full source code though.

Originally Posted by phantom

Happy New Year to everyone, by the way. Hope you had some time to code. I sure did.

Happy new year to you too. I'm really glad to see your wonderful tutorials published again and I hope this will inspire many people as it did inspire me!
I also hope you will continue the serie one day or that at least you'll let us see your progress through IOTDs/Demos ?

Thanks again for an amazing ressource.

Regards,
jonjon.

Indeed, I would dearly love to see what happened with the integration of the photon mapping

Originally Posted by phantom

Ah, tbp! Long time no see! Did you notice that the Intel paper was finally published?

Found out recently about it.

For the unaware: I worked together with an intel dude to write a paper on the SSE/packet tracing details of the latest version of the ray tracer. This cooperation went, well, 'suboptimal' but at least the paper is online:

Heh, please tell me more about that suboptimality

Anyway i already knew about all those features.
I haven't touched my code for like a year, but as i'm sure you have some fancy stuff under your sleeve i'll have to retaliate.
I hate you.

PS: ISP switch = slow reply.

Nah I haven't touched mine either. I've been working on a PocketPC game:
http://www.bik5.com/nc2

But after that I wanna do this:
ftp://download.intel.com/technology/...load/mlrta.pdf

Originally Posted by phantom

Nah I haven't touched mine either. I've been working on a PocketPC game:
http://www.bik5.com/nc2

Seen it, and even if PocketPC aren't my thing it seems very nice.

But after that I wanna do this:
ftp://download.intel.com/technology/...load/mlrta.pdf

After a sub-10s scan it sounds like something we've discussed a bit back then, i will give it a proper read when my gray cells are back online.

When you say you "wanna do this", that means you have a prototype or something?

PS: you've conveniently forgotten to clue me about alluded suboptimalities

This is not the place to bash people.

No I don't have a prototype, and yes, it sounds like what we have discussed. Their approach is however quite a bit simpler to implement than the corner ray approach.

Originally Posted by phantom

But after that I wanna do this:
ftp://download.intel.com/technology/...load/mlrta.pdf

Wow, that is an increadible speed up... Implement it, then write about it

Hi,

Just to share a little bit of my tests.
After reading your article, I went back to my raytracer left in an hided corner of ppModeler. It was slow, and buggy, so I didn't work on it for ages... (other thing to do before).
I have implemented the KD-tree algorithm, and it is quite good.
But, for an animated scene, the KD-tree construction is too slow. The aim of my raytracer is to render animated scene, not static scene. So, the render time should take into account the KD-tree construction.
Here are the results of my tests:

• I stated with a simple KD-tree construction, what 'Vlastimil Havran' called a BSP tree in its thesis (each time, the bounding box (BB) is splitten along is dominant axis, and the position is the middle). The KD-tree construction is recursive and the stop criteria are the number of triangles and the tree depth. Here is the results for a test scene (a human body in T-pose with 15 point lights):
  ######## RENDERING STATS ########
  # triangles: 32240
  # lights: 15
  Kd Tree
  Max Depth: 14
  # nodes: 4693
  # leafs: 2347
  max tri: 1151
  
  dt: 22883 ms
  dtKd: 111 ms
  dtAa: 9541 ms
  # pixels: 480000
  # rays: 679908
  # shadow rays: 1974323
  # aa rays: 199908
  
  # inter. tested: 288809639
  # inter. found: 502498
  #################################
  
• Then, instead of the of a recursive contruction, at each step the program choose the node with the maximum number of triangles, stop critaria being the number of triangles in the node and a maximum number of nodes (not reached in the tests). Here are the results:
  ######## RENDERING STATS ########
  # triangles: 32240
  # lights: 15
  Kd Tree
  Max Depth: 0
  # nodes: 16395
  # leafs: 8198
  max tri: 24
  
  dt: 19718 ms
  dtKd: 2263 ms
  dtAa: 5451 ms
  # pixels: 480000
  # rays: 676830
  # shadow rays: 1952922
  # aa rays: 196830
  
  # inter. tested: 114398083
  # inter. found: 496678
  #################################
  
• Then, I implemented the original version, but with a few modifications. Insted of looking at all possible splitting plane position and axis, the program looks for 256 different values on the 3 axes (this speed up things a lot, but give a KD-tree that can be improved), end criteria being the maximum number of triangles per nodes and the total number of nodes (not reached in the tests). Here are the results:
  ######## RENDERING STATS ########
  # triangles: 32240
  # lights: 15
  Kd Tree
  Max Depth: 0
  # nodes: 3123
  # leafs: 1562
  max tri: 64
  
  dt: 25687 ms
  dtKd: 3235 ms
  dtAa: 9414 ms
  # pixels: 480000
  # rays: 720273
  # shadow rays: 2239234
  # aa rays: 240273
  
  # inter. tested: 298921855
  # inter. found: 570076
  #################################
  

The results shows that simple BSP is better than full KD-tree (taking into account the KD-tree construction).

hope this help,

Thanks a lot for these articles,


craouette

Hello craouette,
let me start with 2 general remarks.

You need to amend your kd-tree (and compiler) to allow for faster updates while animating, see Wald's thesis for some ideas.

Also, Havran has shown that SAH is a better heuristic, in the general case, by a large margin.
Granted, it's quite hard to write an efficient compiler.


From what i see in your report (if i correctly interpret those numbers), it seems that it's your traversal that is really slow, nullifying all partitionning improvments.
You're shooting about ~700k+~2M+~200k ~= 3M ray in ~20s or 150kray/s, and while intersection tests vary between 115M(!) and 300M(!!) rendering times don't change much.

As Wald reports, with a proper kd-tree, you should get a couple of intersection tests per ray: that's one of the reason why kd-tree are so great.

Of course, if you have to mitigate your kd-tree quality for building times it should be a bit worse, but certainly not in the 100 intersections/ray range.

At that point, you'd better try with a bvh

PS: Jacco, that Intel paper was really interesting... but implementing that efficiently with so little details... ouch. Still, really hmm fascinating.

I once got a nice tip to speed up BSP compilation: Evaluate only one in every 10 faces. It's similar to what you propose, but with an important difference: In nodes that have less than 256 tris, you will test every triangle, while 'my' approach will still test 1/10th.

About kd-tree compilation: There is a lot that you can do to improve the process. The compiler that came with the article was slow, as I didn't optimize it at all; I only optimized the ray tracer. After writing the last article, I improved the ray tracer greatly (beats Wald) and also the kd-tree compiler (more than a 10-fold speed increase). However, it can be done much faster. TBP and a collegue managed to compile trees with the same heuristics in half the time and less.

Once you have an optimal kd-tree compiler, you can start experimenting the compilation/rendering balance. There are many ways to make compilation cheaper: I never skip triangles, so that would be a nice start, you could go on all the way down to a straight octree compiler, which should compile 1M polies in a snap (far faster than 22secs). If you have 4 or 5 basic algorithms with some adjustable settings, you could even decide which one to use based on the scene at hand.

If you really have a lot of spare time, you could go for an algo that recalculates only part of the tree. Or, even better, you could insert bounding boxes for animating objects, optimize the tree root (the most expensive part of the compilation) for these, and only update the dynamic portions of the tree this way.

Lots of options.

So I was interested in your explanation of barycentric coordinates, and wrote a test harness to check out your code.

Vertex vertex1 = { {-1.0f,-1.0f, 2.0f }, NULL };
Vertex vertex2 = { { 0.0f, 1.0f, 3.0f }, NULL };
Vertex vertex3 = { { 1.0f,-1.0f, 2.0f }, NULL };
/*
          | (0,1)
       2 /|\ 4
        /1|3\
   ----/--8--\---
      /_5_|_7_\
(-1,-1)   6   (1,-1)
          |
*/
Ray ray1 = { {-0.2f, 0.2f, 1.0f }, { 0.0f, 0.0f, 1.0f } };
Ray ray2 = { {-1.0f, 1.0f, 1.0f }, { 0.0f, 0.0f, 1.0f } };
Ray ray3 = { { 0.2f, 0.2f, 1.0f }, { 0.0f, 0.0f, 1.0f } };
Ray ray4 = { { 1.0f, 1.0f, 1.0f }, { 0.0f, 0.0f, 1.0f } };
Ray ray5 = { {-0.5f,-0.5f, 1.0f }, { 0.0f, 0.0f, 1.0f } };
Ray ray6 = { { 0.0f,-1.5f, 1.0f }, { 0.0f, 0.0f, 1.0f } };
Ray ray7 = { { 0.5f,-0.5f, 1.0f }, { 0.0f, 0.0f, 1.0f } };
Ray ray8 = { { 0.0f, 0.0f, 2.4f }, { 0.0f, 0.0f, 1.0f } };
Triangle triangle;
real distance, coordinate[2];
int hit;
triangle_init(&triangle, NULL, &vertex1, &vertex2, &vertex3);
/* result should be yes, no, yes, no... */
distance = 3.0f; hit = triangle_intersects(&triangle, &ray1, &distance, coordinate);
distance = 3.0f; hit = triangle_intersects(&triangle, &ray2, &distance, coordinate);
distance = 3.0f; hit = triangle_intersects(&triangle, &ray3, &distance, coordinate);
distance = 3.0f; hit = triangle_intersects(&triangle, &ray4, &distance, coordinate);
distance = 3.0f; hit = triangle_intersects(&triangle, &ray5, &distance, coordinate);
distance = 3.0f; hit = triangle_intersects(&triangle, &ray6, &distance, coordinate);
distance = 3.0f; hit = triangle_intersects(&triangle, &ray7, &distance, coordinate);
distance = 3.0f; hit = triangle_intersects(&triangle, &ray8, &distance, coordinate);

Above is a snippet, I simply shoot eight rays at a triangle and tested the distance and coordinates returned against known values.

I am sad to report that the distance values your triangle intersection test returned were out by 10, scalar not linear. This probably still worked, because every test was out by the same amount.

So I went back to the sources you cited and implemented the barycentric intersection test as discribed. Below is a listing of what I came up with...

typedef struct Triangle {
	Vertex *vertex[3];
	real normal[3];
	Material *material;
	/* pre-calculated barycentric coordinate variables */
	int axis;
	real b[3], c[3];
} Triangle;

void triangle_init(Triangle *triangle, Material *material,
	Vertex* vertex1, Vertex* vertex2, Vertex* vertex3)
{
	real *A, *B, *C, *N, b[3], c[3], reciprocal;
	int axis, x, y;
	triangle->material = material;
	triangle->vertex[0] = vertex1;
	triangle->vertex[1] = vertex2;
	triangle->vertex[2] = vertex3;
	A = vertex1->point;
	B = vertex2->point;
	C = vertex3->point;
	N = triangle->normal;
	/* First we calculate the surface normal of the plane
		b = B - A
		c = C - A
		N = b X c
	*/
	vector_sub(b, B, A);
	vector_sub(c, C, A);
	vector_cross(N, b, c);
	/* Next we find the axis to project the barycentric coordinates onto. We
	want the axis with the highest numerical resolution, which corresponds to
	the axis of the normal with the largest absolute dimension. */
	if(fabs_(N[X]) > fabs_(N[Y])) {
		if(fabs_(N[X]) > fabs_(N[Z])) axis = X;
		else axis = Z;
	}
	else {
		if(fabs_(N[Y]) > fabs_(N[Z])) axis = Y;
		else axis = Z;
	}
	triangle->axis = axis;
	/* Here we look up x and y from the axis */
	x = xLut[axis];
	y = yLut[axis];
	/* Next we calculate the first line equation */
	reciprocal = 1.0f / ((b[y] * c[x]) - (b[x] * c[y]));
	triangle->b[X] =  b[x] * -reciprocal;
	triangle->b[Y] =  b[y] *  reciprocal; /* minus is cancelled out */
	triangle->b[Z] = (b[y] * A[x] - b[x] * A[y]) * -reciprocal;
	/* Then we calculate the second line equation */
	triangle->c[X] =  c[x] * reciprocal;
	triangle->c[Y] = -c[y] * reciprocal;
	triangle->c[Z] = (c[y] * A[x] - c[x] * A[y]) * reciprocal;
	/* Finally we make our normal a unit vector */
	vector_normalize(N);
	vertex1->normal = triangle->normal;
	vertex2->normal = triangle->normal;
	vertex3->normal = triangle->normal;
}

int triangle_intersects(Triangle *triangle,
	Ray *ray, real *distance, real coordinate[2])
{
	real *b, *c, *A, *N, *O, *D, u, v, Px, Py;
	real t[3], newDistance;
	int axis, x, y;
	O = ray->origin;
	D = ray->direction;
	A = triangle->vertex[0]->point;
	N = triangle->normal;
	b = triangle->b;
	c = triangle->c;
	axis = triangle->axis;
	/* First we look up x and y from the axis */
	x = xLut[axis];
	y = yLut[axis];
	/* Now we calculate the distance from the ray origin to the triangle
	plane */
	vector_sub(t, O, A);
	newDistance = -vector_dot(t, N) / vector_dot(D, N);
	if(newDistance <= 0 || newDistance > *distance) {
		return FALSE;
	}
	/* Calculate barycentric point */
	Px = O[x] + newDistance * D[x];
	Py = O[y] + newDistance * D[y];
	/* Compute intersection coordinates */
	u = Py * c[X] + Px * c[Y] + c[Z];
	if(!almost_greater_equal(u, 0, MAX_ULPS)) {
		return FALSE;
	}
	v = Py * b[X] + Px * b[Y] + b[Z];
	if(!almost_greater_equal(v, 0, MAX_ULPS)) {
		return FALSE;
	}
	if(!almost_less_equal(u + v, 1, MAX_ULPS)) {
		return FALSE;
	}
	coordinate[X] = u;
	coordinate[Y] = v;
	*distance = newDistance;
	return vector_dot(D, N) > 0 ? -1 : TRUE;
}

This code seems to work as expected. So I was wondering if you could explain the rationale behind the differences between this and your implementation.

Primitive::Primitive( int a_Type, Vertex* a_V1, Vertex* a_V2, Vertex* a_V3 )
{
	m_Type = a_Type;
	m_Material = 0;
	m_Vertex[0] = a_V1;
	m_Vertex[1] = a_V2;
	m_Vertex[2] = a_V3;
	// init precomp
	vector3 A = m_Vertex[0]->GetPos();
	vector3 B = m_Vertex[1]->GetPos();
	vector3 C = m_Vertex[2]->GetPos();
	vector3 c = B - A;
	vector3 b = C - A;
	m_N = b.Cross( c );
	int u, v;
	if (_fabs( m_N.x ) > _fabs( m_N.y))
	{
		if (_fabs( m_N.x ) > _fabs( m_N.z )) k = 0; else k = 2;
	}
	else
	{
		if (_fabs( m_N.y ) > _fabs( m_N.z )) k = 1; else k = 2;
	}
	u = (k + 1) % 3;
	v = (k + 2) % 3;
	// precomp
	real krec = 1.0f / m_N.cell[k];
	nu = m_N.cell[u] * krec;
	nv = m_N.cell[v] * krec;
	nd = m_N.Dot( A ) * krec;
	// first line equation
	real reci = 1.0f / (b.cell[u] * c.cell[v] - b.cell[v] * c.cell[u]);
	bnu = b.cell[u] * reci;
	bnv = -b.cell[v] * reci;
	// second line equation
	cnu = c.cell[v] * reci;
	cnv = -c.cell[u] * reci;
	// finalize normal
	m_N.Normalize();
	m_Vertex[0]->SetNormal( m_N );
	m_Vertex[1]->SetNormal( m_N );
	m_Vertex[2]->SetNormal( m_N );
}
int Primitive::Intersect( Ray& a_Ray, real& a_Dist )
{
	#define ku modulo[k + 1]
	#define kv modulo[k + 2]
	vector3 O = a_Ray.GetOrigin(), D = a_Ray.GetDirection(), A = m_Vertex[0]->GetPos();
	const real lnd = 1.0f / (D.cell[k] + nu * D.cell[ku] + nv * D.cell[kv]);
	const real t = (nd - O.cell[k] - nu * O.cell[ku] - nv * O.cell[kv]) * lnd;
	if (!(a_Dist > t && t > 0)) return MISS;
	real hu = O.cell[ku] + t * D.cell[ku] - A.cell[ku];
	real hv = O.cell[kv] + t * D.cell[kv] - A.cell[kv];
	real beta = m_U = hv * bnu + hu * bnv;
	if (beta < 0) return MISS;
	real gamma = m_V = hu * cnu + hv * cnv;
	if (gamma < 0) return MISS;
	if ((m_U + m_V) > 1) return MISS;
	a_Dist = t;
	return (DOT( D, m_N ) > 0)?INPRIM:HIT;
}

At first I thought you were trying to eliminate the divide, but it is still there?
Cheers, -h

Did you check in the source code if it matters at all if a correct distance is returned? Considering the research you put into it, you're probably right that it does return wrong distances, but I believe the test is only used to accept or reject rays. In that case only the sign of the distance is useful, probably. I'll check it in the code later today.

I once got a nice tip to speed up BSP compilation: Evaluate only one in every 10 faces. It's similar to what you propose, but with an important difference: In nodes that have less than 256 tris, you will test every triangle, while 'my' approach will still test 1/10th.

Yes, that's what I did: min(max(256, numTris/10), 5).
One more question: in your sources, you compute the node surface before splitting. I think it is unnecessary, because you minimise a function, so all contant (added or multiplied) just shift the absolute values, but do not change the one that minimized... Do I miss something ?

From what i see in your report (if i correctly interpret those numbers), it seems that it's your traversal that is really slow, nullifying all partitionning improvments.
You're shooting about ~700k+~2M+~200k ~= 3M ray in ~20s or 150kray/s, and while intersection tests vary between 115M(!) and 300M(!!) rendering times don't change much..

rendering time is dtKd - dt, and it changes (not that much so...).

As Wald reports, with a proper kd-tree, you should get a couple of intersection tests per ray: that's one of the reason why kd-tree are so great.

In this tests, the minimum number of triangles per leaf nodes is 24 => min tests is 24.

I change the kdTree construction, but not the rendering.
The kdtree use the following split strategy:
if bb of triangles included in bb of node -> split empty space
if num of triangles > 100 -> split middle largest axis
else use SAH

Looking at the kdtree generated before, I found the the SAH can be locked into one node: split 18 triangles node => 14 triangles node and 18 triangles node => ...
But I think that is linked to the approximation (test of splitting plane on 5 values for each axis).


Here are the results:

# triangles: 32240
# lights: 15
Kd Tree
Max Depth: 53
# nodes: 11263
# leafs: 5632
max tri: 32
# empty leafs: 2163

dt: 16263 ms
dtKd: 1169 ms
dtAa: 2868 ms
# pixels: 480000
# rays: 589776
# shadow rays: 1262269
# aa rays: 109776

# inter. search: 1852045
# inter. tested: 69198245
# inter. found: 324353

=> 37.3 intersections tested per search => 1.16 node per search.

craouette

Did you check in the source code if it matters at all if a correct distance is returned?

Actually that is your code in the last box (from raytrace7.zip scene.cpp). I though it would be a good idea to include it for comparison, but forgot to mention it.

... but I believe the test is only used to accept or reject rays. In that case only the sign of the distance is useful, probably.

I respectfully disagree. Depth matters when your scene has different primitives like cones, disks and spheres, each with different intersection tests. If things are close together, the ray will register a hit on the wrong thing. Also if you want to have depth queues in your shaders, obviously depth matters. Still I may have buggered up your code somehow while I was testing it, so who knows. Your demo worked fine.

Anyway the bits of your code that still have me scratching my head are as follows...

// precomp
real krec = 1.0f / m_N.cell[k];
nu = m_N.cell[u] * krec;
nv = m_N.cell[v] * krec;
nd = m_N.Dot( A ) * krec;

So here you have nu, nv and nd. If our axis (k) is z then this could be written as

nu = normal[x] / normal[z]
nv = normal[y] / normal[z]
nd = DOT(normal, A) / normal[z]

So these coordinates being projected onto the z axis. But then the next bit in the intersection test I just don't get.

const real lnd = 1.0f / (D.cell[k] + nu * D.cell[ku] + nv * D.cell[kv]);
const real t = (nd - O.cell[k] - nu * O.cell[ku] - nv * O.cell[kv]) * lnd;

Anyway I would be most greatful for your insight. BTW: I really enjoyed your articles on ray tracing and all the other stuff over the years on flipcode. Thanks.

-h

Thinking about this a bit more, I think the real answer I am after here is whether anyone has found a way to get rid of the divide in the barycentric triangle intersection test.

Is there a way to calculate the distance from the ray origin to the triangle plane using a pre calculated reciprical that would replace the following...

Vector *O, *D, *N, *A;
O = ray->origin;
D = ray->direction;
N = triangle->normal
A = triangle->vertex[0]->point;

distance = -DOT(O - A, N) / DOT(D, N);

With a bunch of multiplies?

-h

Kd-trees partition space and their traversal is cheap.
If you don't exploit those 2 properties, you're better using another hierarchy.
Like i said it seems you're missing those 2 points.

First you bundle lots of primitives in leaves.

Second, even in your last report, you're shooting ~2M rays in ~15s or 150k ray/s.
That's off by an order of magnitude (for mono ray in a straight fp implementation), see http://homepages.paradise.net.nz/nickamy/benchmark.html for various aproaches/numbers.
Take note that a (box) bvh for such isolated model is more or less on par with a kd-tree and is much more animation friendly.

On top of that, kd-tree compilation still only takes 7% of your total rendering time. It seems to me that you're barking at the wrong tree

There's many solutions to your problem.
I'd suggest using a SAH bvh as it's way easier to write in a efficient fashion and amend for animation while still giving interesting speedups, but i have to admit i'm biased
The animated kd-tree road is way bumpier.

But first things first, please fix your traversal ASAP.

yeah, you're right... I'm working on it, but without much succes I must say.
any direction ? i suppose reducing the size of the main structures mainly...


craouette

Craouette: Check out Wald's thesis. It has some suggestions for fast kd-tree traversal and triangle intersections, even without SIMD.