25

How do I calculate the intersection between a ray and a plane?

Code

This produces the wrong results.

float denom = normal.dot(ray.direction);

if (denom > 0)
{
    float t = -((center - ray.origin).dot(normal)) / denom;

    if (t >= 0)
    {
        rec.tHit = t;
        rec.anyHit = true;
        computeSurfaceHitFields(ray, rec);
        return true;
    }
}

Parameters

ray represents the ray object.
ray.direction is the direction vector.
ray.origin is the origin vector.
rec represents the result object.
rec.tHit is the value of the hit.
rec.anyHit is a boolean.

My function has access to the plane:
center and normal defines the plane

Quonux
  • 2,975
  • 1
  • 24
  • 32
Pontus Magnusson
  • 632
  • 1
  • 10
  • 25
  • 3
    When you say this does not work, what specifically isn't working? Is it crashing, are the results wrong? does it fail to compile? – Matt Coubrough Jun 01 '14 at 00:18
  • 2
    Is your normal vector guaranteed to point away from the ray origin? Otherwise, *denom < 0* might very well still produce an intersection. – wonce Jun 01 '14 at 00:21

4 Answers4

36

As wonce commented, you want to also allow the denominator to be negative, otherwise you will miss intersections with the front face of your plane. However, you still want a test to avoid a division by zero, which would indicate the ray being parallel to the plane. You also have a superfluous negation in your computation of t. Overall, it should look like this:

float denom = normal.dot(ray.direction);
if (abs(denom) > 0.0001f) // your favorite epsilon
{
    float t = (center - ray.origin).dot(normal) / denom;
    if (t >= 0) return true; // you might want to allow an epsilon here too
}
return false;
Trillian
  • 6,207
  • 1
  • 26
  • 36
  • 1
    What is center and what is ray origin? – Lily Dec 06 '15 at 07:05
  • 1
    `center` here is a known point on the plane. This uses the "Point-normal" description for the plane. – Joan Charmant Mar 06 '16 at 10:34
  • 1
    I have to point out that you should use `fabs()` and not `abs()` because `abs()` converts the input parameter to `int`. – Tara Jun 24 '16 at 18:56
  • 5
    @Duckdoom5: I was obviously talking about C++, because the question is tagged as C++. I assume you are talking about `std::abs()` (which has different overloads) and not `abs()` (which only deals with integers). Since the answer contains `abs()` instead of `std::abs()` and doesn't use a `using` directive either, I assume this to be a mistake. See here for more information: http://stackoverflow.com/questions/21392627/abs-vs-stdabs-what-does-the-reference-say – Tara Oct 11 '16 at 01:56
  • 1
    @Duckdoom5 I hope you'll delete that comment. It's wrong and you'll confuse novices with wrong information. – Glenn Maynard May 21 '17 at 06:13
  • 1
    @GlennMaynard You are right, I was assuming `std::abs()` was used here, so I'll delete the comment. – Duckdoom5 Oct 05 '17 at 12:32
  • I used this code and it works fine. Math explanation: https://www.cs.princeton.edu/courses/archive/fall00/cs426/lectures/raycast/sld017.htm and t is the distance from ray to plane. – Ardiya Jan 18 '21 at 11:14
  • Are there no bounds to the plane? I don't understand where its size is determined. – Ryad Shaanbi Nov 30 '21 at 15:54
8

First consider the math of the ray-plane intersection:

In general one intersects the parametric form of the ray, with the implicit form of the geometry.

So given a ray of the form x = a * t + a0, y = b * t + b0, z = c * t + c0;

and a plane of the form: A x * B y * C z + D = 0;

now substitute the x, y and z ray equations into the plane equation and you will get a polynomial in t. you then solve that polynomial for the real values of t. With those values of t you can back substitute into the ray equation to get the real values of x, y and z. Here it is in Maxima:

enter image description here

Note that the answer looks like the quotient of two dot products! The normal to a plane is the first three coefficients of the plane equation A, B, and C. You still need D to uniquely determine the plane. Then you code that up in the language of your choice like so:

Point3D intersectRayPlane(Ray ray, Plane plane)
{
    Point3D point3D;

    //  Do the dot products and find t > epsilon that provides intersection.


    return (point3D);
}
vwvan
  • 417
  • 1
  • 4
  • 12
5

Math

Define:

  • Let the ray be given parametrically by q = p + t*v for initial point p and direction vector v for t >= 0.

  • Let the plane be the set of points r satisfying the equation dot(n, r) + d = 0 for normal vector n = (a, b, c) and constant d. Fully expanded, the plane equation may also be written in the familiar form ax + by + cz + d = 0.

The ray-plane intersection occurs when q satisfies the plane equation. Substituting, we have:

d = -dot(n, q)
  = -dot(n, p + t * v)
  = -dot(n, p) + t * dot(n, v)

Rearranging:

t = -(dot(n, p) + d) / dot(n, v)

This value of t can be used to determine the intersection by plugging it back into p + t*v.

Example implementation

std::optional<vec3> intersectRayWithPlane(
    vec3 p, vec3 v,  // ray
    vec3 n, float d  // plane
) {
    float denom = dot(n, v);

    // Prevent divide by zero:
    if (abs(denom) <= 1e-4f)
        return std::nullopt;

    // If you want to ensure the ray reflects off only
    // the "top" half of the plane, use this instead:
    //
    // if (-denom <= 1e-4f)
    //     return std::nullopt;

    float t = -(dot(n, p) + d) / dot(n, v);

    // Use pointy end of the ray.
    // It is technically correct to compare t < 0,
    // but that may be undesirable in a raytracer.
    if (t <= 1e-4)
        return std::nullopt;

    return p + t * v;
}
Mateen Ulhaq
  • 24,552
  • 19
  • 101
  • 135
  • I wrote this answer to provide a complete working example for the `ax + by + cz + d = 0` form. Additionally discusses what to do if we only wish to reflect rays off a normally oriented plane (i.e. the front side of the plane... not its back). – Mateen Ulhaq Nov 12 '19 at 13:41
  • Also, see [this answer](https://physics.stackexchange.com/questions/435512/snells-law-in-vector-form/436252#436252) for a derivation of Snell's Law in vector form, and its solution for the parameter `t`. – Mateen Ulhaq Nov 14 '19 at 13:11
  • t = - (dot(n, p) + d) / dot(n, v) should be t = - (dot(n, p0) + d) / dot(n, v) – KeithB Nov 18 '21 at 11:35
4

implementation of vwvan's answer

Vector3 Intersect(Vector3 planeP, Vector3 planeN, Vector3 rayP, Vector3 rayD)
{
    var d = Vector3.Dot(planeP, -planeN);
    var t = -(d + Vector3.Dot(rayP, planeN)) / Vector3.Dot(rayD, planeN);
    return rayP + t * rayD;
}
Bas Smit
  • 645
  • 1
  • 7
  • 19