Ray-bounded plane intersection

Question

I'm trying to write a ray tracer in my freetime. Currently trying to do ray - bounded plane intersections.

My program is already working with infinite planes. I'm trying to work out the math for non-infinite planes. Tried to google, but all of the resources talk only about infinite planes.

My plane has a corner point (called position), from which two vectors (u and v) extend (their length correspond to the length of the sides). The ray has an origin and a direction.

First I calculate the intersection point with an infinite plane with the formula

t = normal * (position - origin) / (normal * direction)

The normal is calculated as a cross product of u and v. Then with the formula

origin + direction * t

I get the intersection point itself.

The next step is checking if this point is in the bounds of the rectangle, and this is where I'm having trouble.

My idea was to get the relative vector intersection - position that is extending from the corner of the plane to the intersection point, then transform it to a new basis of u, normal and v then check if the lengths of the transformed vectors are shorter than the u and v vectors.

bool BoundedPlane::intersect(const Vec3f &origin, const Vec3f &direction, float &t) const {
    t = normal * (position - origin) / (normal * direction);
    Vec3f relative = (origin + direction * t) - position;

    Mat3f transform{
        Vec3f(u.x, normal.x, v.x),
        Vec3f(u.y, normal.y, v.y),
        Vec3f(u.z, normal.z, v.z)
    };

    Vec3f local = transform.mul(relative);
    return t > 0 && local.x >= 0 && local.x <= u.x && local.z <= 0 && local.z <= v.z;

}

At the end I check if t is larger than 0, meaning the intersection is in front of the camera, and if the lengths of the vectors are inside the bounds. This gives me a weird line:

.

The plane should appear below the spheres like this:

(this used manual checking to see if it appears correctly if the numbers are right).

I'm not sure what I'm doing wrong, and if there's an easier way to check the bounds. Thanks in advance.

Edit1:

I moved the transformation matrix calculations into the constructor, so now the intersection test is:

bool BoundedPlane::intersect(const Vec3f &origin, const Vec3f &direction, float &t) const {
    if (!InfinitePlane::intersect(origin, direction, t)) {
        return false;
    }

    Vec3f local = transform.mul((origin + direction * t) - position);

    return local.x >= 0 && local.x <= 1 && local.z >= 0 && local.z <= 1;
}

The transform member is the inverse of the transformation matrix.

Answer 1

Could I suggest another approach? Consider the frame with origin position and basis vectors

u = { u.x, u.y, u.z }  
v = { v.x, v.y, v.z }  
direction = { direction.x, direction.y, direction.z} 

Step 1: Form the matrix

M = {
      {u.x,  v.x,  direction.x},
      {u.y,  v.y,  direction.y},
      {u.z,  v.z,  direction.z} 
    }

Step 2: Calculate the vector w, which is a solution to the 3 x 3 system of liner equations

M * w = origin - position, i.e.

w = inverse(M) * (origin - position);

Make sure that direction is not coplanar with u, v, otherwise there is no intersection and inverse(M) does not exist.

Step 3: if 0.0 <= w.x && w.x <= 1.0 && 0.0 <= w.y && w.y <= 1.0 then the line intersects the parallelogram spanned by the vectors u, v and the point of intersection is

w0 = { w.x, w.y , 0 };

intersection = position + M * w0; 

else, the line does not intersect the parallelogram spanned by the vectors u, v

The idea of this algorithm is to consider the (non-orthonormal) frame position, u, v, direction. Then the matrix M changes everything in the coordinates of this new frame. In this frame, the line is vertical, parallel to the "z-"axis, the point origin has coordinates w, and the vertical line through w intersects the plane at w0.

Edit 1: Here is a templet formula for the inverse of a 3x3 matrix:

If original matrix M is

a  b  c

d  e  f

g  h  i

inverse is

(1 / det(M)) * { 
                 {e*i - f*h,     c*h - b*i,     b*f - c*e},

                 {f*g - d*i,     a*i - c*g,     c*d - a*f},

                 {d*h - e*g,     b*g - a*h,     a*e - b*d},
                } 

where

det(M) = a*(e*i - f*h) + b*(f*g - d*i) + c*(d*h - e*h)

is the determinant of M.

So the inversion algorithm can be as follows:

Given

M = {
      {a,  b,  c},
      {d,  e,  f},
      {g,  h,  i},
    }

1. Calculate

inv_M = { 
           {e*i - f*h,     c*h - b*i,     b*f - c*e},

           {f*g - d*i,     a*i - c*g,     c*d - a*f},

           {d*h - e*g,     b*g - a*h,     a*e - b*d},
         };  

2. Calculate

det_M = a*inv_M[1][1] + b*inv_M[2][1] + c*inv_M[3][1]; 

3. Return inverse matrix of M

inv_M = (1/det_M) * inv_M;

Edit 2: Let's try another approach in order to speed things up.

Step 1: For each plane, determined by the point position and the two vectors u and v, precompute the following quatntities:

normal = cross(u, v); 

u_dot_u = dot(u, u);

u_dot_v = dot(u, v);

v_dot_v = dot(v, v); // all these need to be computed only once for the u and v vectors

det = u_dot_u * v_dot_v - u_dot_v * u_dot_v; // again only once per u and v 

Step 2: Now, for a given line with point origin and direction direction, as before, calculate the intersection point int_point with the plane spanned by u and v:

t = dot(normal,  position - origin) / dot(normal, direction);

int_point = origin + t * direction;

rhs = int_point - position;

Step 3: Calcualte

u_dot_rhs = dot(u, rhs);

v_dot_rhs = dot(v, rhs);

w1 = (v_dot_v * u_dot_rhs - u_dot_v * v_dot_rhs) / det;

w2 = (- u_dot_v * u_dot_rhs + u_dot_u * v_dot_rhs) / det;

Step 4:

if (0 < = w1 && w1 <= 1 && 0 < = w2 && w2 <= 1 ){
   int_point is in the parallelogram;
}
else{
   int_point is not in the parallelogram;
}

So what I am doing here is basically finding the intersection point of the line origin, direction with the plane given by position, u, v and restricting myself to the plane, which allows me to work in 2D rather than 3D. I am representing

int_point = position + w1 * u + w2 * v;

rhs = int_point - position = w1 * u + w2 * v

and finding w1 and w2 by dot-multiplying of this vector expression with the basis vectors u and v, which results in a 2x2 linear system, which I am solving directly.