February 11, 2008

Raytracing: Intersection with a sphere

Follow-up to Raytracing: Intersection with a plane from Nick's blog

Last time I went through intersection of a ray with a plane. Plane’s are all well and good, but you can’t have a ray tracer without spheres everywhere :-)

This article’s another maths-heavy one I’m afraid – more vectors and dot products, this time with a quadratic equation thrown in at the end. Enjoy.

Equation for our ray:

$P = O + Dt$

(P is the position, O is the origin of the ray, D is the direction of the ray, t is the parametric variable)

Equation for a sphere:

$(x - C_x)^2 + (y - C_y)^2 + (z - C_z)^2 = r^2$

(x, y and z represent points around the surface of the sphere, C is the centre of the sphere, r is the radius of the sphere)

As with the plane intersection in the previous article, we solve the equations simultaneously, setting P = (x,y,z).

First we’ll rewrite the sphere equation in vector form. Summing the squares of the x, y and z components is the same as the dot product.

$(P - C) \cdot (P - C) = r^2$

Now substitute in the equation for our ray:

$(O + Dt - C) \cdot (O + Dt - C) = r^2$

We’re trying to solve for t, so let’s break up the above equation using the distributivity of the dot product:

$(O + Dt - C) \cdot O + (O + Dt - C) \cdot Dt + (O + Dt - C) \cdot (-C)= r^2$

Now extract the Dt terms:

$(O - C) \cdot O + (O - C) \cdot Dt + (O - C) \cdot (-C) + (Dt) \cdot O + (Dt) \cdot (Dt) + (Dt) \cdot(-C) = r^2$

Collect like terms:

$(O - C) \cdot O + 2(O - C) \cdot Dt + (O - C) \cdot (-C) + (Dt) \cdot (Dt) = r^2$

Now if we take the scalar t outside of the dot products, and move the r squared to the left hand side we get:

$(O - C) \cdot O + (O - C) \cdot (-C) + t(2(O - C) \cdot D) + t^2(D \cdot D) - r^2 = 0$

Which we can solve using the quadratic formula

Where
$a = D \cdot D = 1 \\ b = 2(O - C) \cdot D \\ c = (O-C) \cdot O + (O-C) \cdot (-C) - r^2 = (O-C) \cdot (O-C) - r^2$

This is where I notice that I’ve reused “C” as a variable name – don’t get confused between “C” for centre” and the “c” in the quadratic formula.

Here I’ve also simplified A and C. A is always 1 in this case, as the direction of our ray “D” is normalised. I’ve also combined C into a single dot product.

Here’s the quadratic formula, which we can simplify by removing the “a” multiplicand.

$t = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a} \\ = \frac{-b \pm \sqrt{b^2 - 4c}}{2}$

The rest is just a case of plugging in your values.

In order to solve this, you should check the value of the determinant (the bit inside the square root). If it’s < 0 then we have no real solutions, and therefore no intersections.

If the determinant is >= 0, then you need to check both solutions to the equation(-b + ... and -b – ...). We want the closest (lowest value) positive solution. A -ve value means the intersection is behind the ray. Also, if you have one +ve and one -ve solution, then your ray has been cast from inside the sphere.

One final thing worth noting is that there are special methods used to compute quadratic solutions which help to reduce the effect of rounding errors.

Ok, well hope that makes sense. Let me know if you have any questions or comments. Next post will probably cover lighting, so you can actually get planes and spheres in colour!

9 comments by 2 or more people

1. Mathew Mannion

Ok, well hope that makes sense.

No

11 Feb 2008, 14:31

2. Any bit in particular?

Luckily you don’t really need to understand all of it. You just need to know which values to plug in to b and c to get the two solutions out of the quadratic at the end. The next article on lighting should be a bit more easy going – with a pretty picture at the end :-)

11 Feb 2008, 14:46

3. Simon Claret

Makes sense to me. I would add one thing.

You might need to calculate the surface normal at the point of intersection (so that you can calculate the color at the point of intersection using the phong model). This can be done by:

n = (O-C)/r

23 Mar 2008, 21:05

4. Yeah, the normal calculation is needed for that. You also need the surface normal for things such as reflection & refraction calculations.

24 Mar 2008, 01:46

5. Nathan

Thanks for that Nick.
I’m just curious about one thing. Most tutorials that I see, their implementation omits the constants, do they all just cancel out?

04 Jun 2008, 03:46

6. Hi Nathan,

Which constants are you referring to?

Everything above the calculation of the quadratic parameters a, b and c is provided for explanation (so wouldn’t be included in an implementation). The calculation of the quadratic formula results in an implementation fairly similar to most others I’ve seen.

04 Jun 2008, 10:02

7. Nathan

Nick,
What I’m reffering to is the 2 when calculating b, the 4 that c is multiplied by, and the divisor in the quadratic formula. It looks like they cancel out, but I just wanted to make sure.

07 Jun 2008, 17:11

8. Rstockham

Just want to make sure with the math. All variables that are rays or positions multiplied with another ray or position use dot product, otherwise simple multiplication. In future examples to make it easier to read, could you differentiate dot product and simple multiplication.

14 Aug 2008, 05:52

9. Hello,

The dot product is indicated by the dot. Simple multiplication is used when multiplying with scalar values (indicated by lower case names). Perhaps I should’ve made that more clear, but I don’t think there’s any need to differentiate further. (using a cross to indicate multiplication would likely only be confused with the cross product)

17 Aug 2008, 12:21

You are not allowed to comment on this entry as it has restricted commenting permissions.

February 2008

Mo Tu We Th Fr Sa Su
Jan |  Today  | Mar
1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29