Object Intersection

Some notes for basic intersection tests. These are not necessarily the most optimized or clever algorithms, but they are general and simple. Much craftier solutions can be found for special cases, and there are links at the bottom for those.


Ray-Plane
The plane is defined by an arbitrary point on its surface,
 a, and a normal vector nThe ray is defined by its position, p, and its direction dThe point the ray intersects the plane is at xSince x is on the plane's surface, a vector between it and any other point on the plane should be orthogonal to the plane's normal; the dot product between these two vectors must be zero. We already know the point a, so:

n \cdot (x - a) = 0

It's also clear that the point x should be on the same line that the ray is on. We can get this point by adding the ray's direction times some scalar to the ray's position. If t is positive, the intersection point is in the direction of v; if it is negative, the ray doesn't intersect the plane and the ray is facing away from the plane's surface:

x = p+dt

With these two equations, we can solve for t to find the intersection:

n\cdot(p + dt - a) = 0
n\cdot(p-a) + n\cdot dt = 0


t = -n\cdot (p-a)/n\cdot d
x = p + d(-n\cdot(p-a)/n\cdot d))

Vec3 rayPlaneIntersect(Ray r, Plane p) {
    float nDotD = dot(p.n, r.dir);

    if (nDotD == 0)
        return null;

    return r.pos + r.dir * (-dot(p.n, r.pos - p.a) / nDotD);
}

Ray-Triangle
There are two parts to finding the intersection of a ray with a triangle:
1. Intersect the ray with the triangle's plane
2. Check if the intersection point is inside all three edges

The plane is defined by any vertex of the triangle and the normal vector
n. Using the algorithm above, it is easy to find out where the ray would intersect this plane. This gives the point on the triangle's plane x. While the ray may intersect the plane, it's not yet clear if it intersects the triangle, so we check if the intersection point is on the correct side of the edges of the triangle. 

Looking at the triangle from the front (the normal is facing toward you), the vertices should be ordered counter-clockwise. The three edges of the triangle are in red, below. For x to be inside the triangle's edges, it must be "left" of each edge: if you are standing on the starting vertex of the edge facing the end vertex, x should be on your left side. If this is the case, the cross product of the edge with the vector from the starting vertex to x should be in the same direction as the normal vector, n

If x is to the "right" side of the edge, the cross product would give a vector pointing in the opposite direction of n. If the cross products are all in the same direction as n, the intersection point x is inside the triangle and is valid. In other words, x is inside the triangle if these three conditions are true:

((b-a)\times(x-a))\cdot n> 0
((c-b)\times(x-b))\cdot n> 0
((a-c)\times(x-c))\cdot n> 0




Vec3 rayTriangleIntersect(Ray r, Triangle t) {
    Vec3 x = rayPlaneIntersect(r, new Plane(t.a, t.n));

    // ray doesn't intersect the plane, so it can't possibly intersect triangle
    if (x == null)
        return null;

    // check the 3 edges: if one fails, the whole test fails
    if (dot(cross(b-a, x-a), n) < 0)) return null;
    if (dot(cross(c-b, x-b), n) < 0)) return null;
    if (dot(cross(a-c, x-c), n) < 0)) return null;

    return x;
}

Ray-Polygon
Intersecting a ray with a convex polygon is just a generalized version of the ray-triangle intersection. Since all polygons can be split into convex polygons, this method will work in most situations.

Vec3 rayPolygonIntersect(Ray r, Polygon p) {
    Vec3 x = rayPlaneIntersect(r, new Plane(p.v[0], p.n));
    if (x == null)
        return;

    int n = p.v.length;
    for (int i=0; i<n; i++)
        if (dot(cross(v[(i+1)%n] - v[i], x - v[i]), p.n) < 0) return null;

    return x;
}



Ray-Sphere
Again, the ray is given by its position p and direction d, and the intersection with the sphere is at x. For x to be on the sphere, it's distance from the center of the sphere must be the same as the sphere's radius r:

|x-c| = r

The distance between x and c is the same as the length of the vector between the two points, and the length of a vector can be determined with the dot product:

\sqrt{(x-c)\cdot (x-c)}  = r
(x-c)\cdot (x-c) = r^2

Once again, x can be found using the ray's position and direction, so this can be substituted in:

p+dt=x
(p+dt - c) \cdot (p+dt-c) = r^2

This equation can be rearranged a bit and then expanded:

(dt+(p-c))\cdot (dt+(p-c)) = r^2
(d\cdot d)t^2 +2d\cdot (p-c)t + (p-c)\cdot (p-c) - r^2 = 0

The above is in the form of a quadratic equation:

at^2 + bt + c = 0

where 

a = d\cdot d
b = 2d\cdot (p-c)
c = (p-c)\cdot (p-c) - r^2

The quadratic formula can be used to find t:

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

t = \frac{-2d\cdot(p-c) \pm \sqrt{(2d\cdot(p-c))^2 -4(d\cdot d)((p-c)\cdot(p-c)-r^2)}}{2(d\cdot d)}


Vec3 intersect(Ray r, Sphere s) {
Vec3 pMinusC = r.pos - s.c;
float a = dot(r.dir, r.dir);
float b = 2 * dot(d, pMinusC);
float c = dot(pMinusC, pMinusC) - s.r * s.r;
float disc = b * b - 4 * a * c;
// no solution if discriminant is negative
if (disc < 0)
return null;
float t;
if (disc == 0) {
        // one solution, barely touches side
t = (float)(-b + Math.sqrt(disc)) / (2 * a);
    } else {
// two solutions: the smallest t is the closest intersection
        float temp = sqrt(disc) / (2 * a);
float t1 = -b + temp;
float t2 = -b - temp;
t = min(t1, t2);
}
return r.pos + r.dir * t;
}


The structures and methods referenced in the code above:

struct Vec3 {
    float x, y, z;
}

struct Ray {
    Vec3 pos;      // position   
    Vec3 dir;      // direction
}

struct Plane {
    Vec3 a;        // arbitrary point on the plane
    Vec3 n;        // normal
}

struct Triangle {
    Vec3 a, b, c;  // the three vertices ordered CCW
    Vec3 n;        // normal
}

struct Polygon {
    Vec3[] verts;  // vertices ordered CCW
    Vec3 n;        // normal
}

float dot(Vec3 a, Vec3 b) {
    return a.x * b.x + a.y * b.y + a.z * b.z
}

Vec3 cross(Vec3 a, Vec3 b) {
    Vec3 c;
    c.x = a.y * b.z - a.z * b.y;
    c.y = a.z * b.x - a.x * b.z;
    c.z = a.x * b.y - a.y * b.x;
    return c;
}


Useful Links:


Comments