Cylinder intersection

For the first time in a while I’ve added a new primitive type to the raytracer. The more I’ve messed about with procedural geometry the more I’ve become frustrated by the limit of building objects out of cubes and spheres. The remaining obvious primitive type is the cylinder. Once this is added to a raytracer you can make all kinds of fun rounded shapes. More on this later.

First of all though, lets look at the maths of tracing an arbitrary ray with an arbitrary cylinder. The first simplification is to get rid of the arbitrary bit about the cylinder. As discussed in this article about transformations a major simplification comes about by translating the ray so that you only need to calculate the intersection with a primitive based at the origin.

Intersecting a ray against a cylinder with a radius of 1, based at the origin, and one unit high, is a much simpler task. In fact, it simplifies right down to two dimensions. As described in Neil Dodgson’s raytracing cribsheet, the ray tracing algorithm just needs to check the ray in two dimensions against a circle centred at the origin. The formula for a circle is

x2 + y2 = 1

The formula for a ray is

P(t) = E + tD

Substituing the ray formula into the circle formula gives

(xE + t.xD)2 + (yE + t.yD)2 = 1

Which can be expanded out to a quadratic of the form

a.t2 + b.t + c = 0
where
a = xD2 + yD2
b = 2.xD.xE + 2.yD.yE
c = xE2 + yE2 - 1

And as you may know a quadratic has 0, 1 or 2 solutions specified by the formula

t = (-b +/- sqrt(b*b - 4*a*c)) / 2*a

Next up then, we need to convert these lovely mathematical formula into a bit of code that allows you to calculate the values of t. In C++ this code simplifies down quickly to.

float t0, t1;
const DVector3& start = rRayContext.m_Ray.GetStart();
const DVector3& direction = rRayContext.m_Ray.GetDirection();

// a=xD2+yD2, b=2xExD+2yEyD, and c=xE2+yE2-1. 
float a = direction.mComponent[0] * direction.mComponent[0]
	+ direction.mComponent[2] * direction.mComponent[2];

float b = 2 * start.mComponent[0] * direction.mComponent[0]
	+ 2 * start.mComponent[2] * direction.mComponent[2];

float c = start.mComponent[0] * start.mComponent[0]
	+ start.mComponent[2] * start.mComponent[2]
	- 1;

float b24ac = b*b - 4*a*c;
if (b24ac<0)
	return false;
	
float sqb24ac = sqrtf(b24ac);
t0 = (-b + sqb24ac) / (2 * a);
t1 = (-b - sqb24ac) / (2 * a);

Next up we order the values of t so that they're in ascending order and calculate the y value at the two intersection points with the circle. Note that the x and y discussed in the formulae are actually x and z when you convert everything into a 3d raytracer. Hence why y is the height of the intersection point.

if (t0>t1) {float tmp = t0;t0=t1;t1=tmp;}

float y0 = start.mComponent[1] + t0 * direction.mComponent[1];
float y1 = start.mComponent[1] + t1 * direction.mComponent[1];

So far we've worked out that the ray intersects a two dimensional circle if we eliminate the y coordinate. Now we can start to look at three dimensions to work out whether or not the ray intersects the cylinder.

We're raytracing against a capped cylinder with caps at y = 1, and y = -1, so there are now five different possibilities for intersections.

if y0 >1, and y1 >1, then the ray misses the cylinder entirely.

if y0 >1, and y1 <1, then the ray hits the cylinder cap placed at +1. if y0 <1 and y0 >-1, then the ray intersects the side of the cylinder.

if y0 <-1 and y1 >-1, then the ray hits the cylinder cap placed at -1

if y0 <-1 and y1 <-1 then the ray misses the cylinder entirely. The code for this gets relatively lengthy as there are quite a few conditions, there's also a little bit of fiddling to get the cylinder cap collision points.

if (y0<-1)
{
	if (y1<-1)
		return false;
	else
	{
		// hit the cap
		float th = t0 + (t1-t0) * (y0+1) / (y0-y1);
		if (th<=0) return false;
	
		out_Response.mHitPosition = start + (direction*th);
		out_Response.mNormal = DVector3(0, -1, 0);
		return true;
	}
}
else if (y0>=-1 && y0<=1)
{
	// hit the cylinder bit
	if (t0<=0) return false;
		
	out_Response.mHitPosition = start + (direction*t0);
	out_Response.mNormal = DVector3(out_Response.mHitPosition.mComponent[0], 0, out_Response.mHitPosition.mComponent[2]);
	out_Response.mNormal.Normalise();
return true;
}
else if (y0>1)
{
	if (y1>1)
		return false;
	else
	{
		// hit the cap
		float th = t0 + (t1-t0) * (y0-1) / (y0-y1);
		if (th<=0) return false;

		out_Response.mHitPosition = start + (direction*th);
		out_Response.mNormal = DVector3(0, 1, 0);
		return true;
	}
}

return false;

And after a little bit more messing around with the Wootracer I can now add cylinders into the raytraced scenes I'm building. First off lets start with the simplest cylinder scene we can raytrace.

rule main {
diff=vec(1,0.5,0)
cylinder
}
Raytraced cylinder primitive

Raytraced cylinder primitive

Now on to a capsule (sphere capped cylinder).

rule main {
diff=vec(1,0.5,0)
sphere
pos.y+=0.5
cylinder
pos.y+=0.5
sphere
}
Raytraced capsule

Raytraced capsule

And very swiftly we can take this on to raytracing square objects made out of cylinders and spheres, stacked one on top of the other.

rule main {
repeat(35) { 
 { roundsq }
 scale*=vec(0.93,0.93,0.93)
 pos.y += 1
 ry += 14
 }
}

rule roundsq {
diff = vec(0.5:0.7, 0.7:0.9, 0.9:1)

pos.x -= 4
pos.z += 4.5
rx += 90
rz -= 90

repeat (4) {
 sphere
 pos.y += 0.5
 rz += 90
 scale *= vec(1, 9, 1)
 cylinder
 scale *= vec(1, 1 / 9, 1)
 pos.y += 8.5
 }
}
Stack of rounded squares

Stack of rounded squares

If you'd like to have a go at building these objects yourself then why not have a go? Wooscripter is now available to download, just follow the link!

You may also like...

3 Responses

  1. BobDylan says:

    This is a really elegant implementation. How would you generalize this to all cones?

Leave a Reply

Your email address will not be published. Required fields are marked *

Spam Protection *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>