Intersection of Ray and Plane in C++
In this tutorial, first, we will get complete guidance on the ray-plane intersection in C++.
First, we will discuss vector operations with their implementation. Next, we will discuss the concept of ray-plane intersection and its implementation in C++.
Use Vector Operations in C++
Here, we will cover related vector operations and their implementation in C++. We need the following functions:
- Minus operator: To calculate the distance between two points.
- Dot product: A function to calculate the dot product between two vectors, which results in a real number.
- Multiplication operator: To calculate the cross product of two vectors.
- Multiplication operator with float parameter: To calculate scalar multiplication between vector and scalar value.
- Magnitude function: To calculate the vector’s magnitude.
- Normalise function: To calculate the normal of a vector.
- Finally, the stream operator is overloaded to display the vector.
Here is the complete vector 3D class. The main function checks/demonstrates member functions of the vector 3D class.
#include <math.h>
#include <iostream>
using namespace std;
class Vec3 {
float x;
float y;
float z;
public:
Vec3() { x, y, z = 0; }
Vec3(float x, float y, float z) {
this->x = x;
this->y = y;
this->z = z;
}
Vec3 &operator+=(const Vec3 &b) {
x = x + b.x;
y = y + b.y;
z = z + b.z;
return *this;
}
Vec3 &operator-=(const Vec3 &b) {
x = x - b.x;
y = y - b.y;
z = z - b.z;
return *this;
}
Vec3 operator+(const Vec3 &b) {
Vec3 newV = *this;
newV += b;
return newV;
}
Vec3 operator-(const Vec3 &b) {
Vec3 newV = *this;
newV -= b;
return newV;
}
Vec3 operator*(const Vec3 &b) { // cross operator
Vec3 newV;
newV.x = y * b.z - z * b.y;
newV.y = x * b.z - z * b.x;
newV.z = x * b.y - y * b.x;
return newV;
}
Vec3 &operator*=(const float s) { // Dot equal operator
x = x * s;
y = y * s;
z = z * s;
return *this;
}
Vec3 operator*(const float s) { // Dot equal operator
Vec3 newV = *this;
return newV *= s;
}
float mag() { return sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); }
Vec3 &normalize() {
double mag = this->mag();
x /= mag;
y /= mag;
z /= mag;
return *this;
}
// Dot operator
float dot(const Vec3 &b) { return x * b.x + y * b.y + z * b.z; }
friend ostream &operator<<(ostream &out, const Vec3 &v) {
out << "(" << v.x << ", " << v.y << ", " << v.z << ")\n";
return out;
}
};
int main() {
Vec3 v1(-5, 7, 2);
Vec3 v2(4, 12, 1);
Vec3 v3(1, 1, 1);
cout << "Vec1 = " << v1 << endl;
v1 += v2;
cout << "Result of adding Vec2 in Vec1 = " << v1 << endl;
v1 -= v3;
cout << "Result of subtracting Vec3 from vec1 = " << v1 << endl;
v1 *= 5; // dot operation
cout << "Resultant after the scaling of Vec1 with a value of 5 = " << v1
<< endl;
cout << "Magnitude for Vec1 = " << v1.mag() << endl;
v1.normalize();
cout << "Vec1 after normalization = " << v1 << endl;
cout << "Dot product of Vec1 and Vec2 = " << v1 * v2 << endl;
return 0;
}
Here is the output of the main function, showing accurate implementation of the vector 3D class.
Vec1 = (-5, 7, 2)
Result of adding Vec2 in Vec1 = (-1, 19, 3)
Result of subtracting Vec3 from vec1 = (-2, 18, 2)
Resultant after the scaling of Vec1 with a value of 5 = (-10, 90, 10)
Magnitude for Vec1 = 91.1043
Vec1 after normalization = (-0.109764, 0.987878, 0.109764)
Dot product of Vec1 and Vec2 = (-0.329293, -0.548821, -5.26868)
Intersection of Ray and Plane in C++
From here onwards, we assume that the reader is quite comfortable with the concepts of vector operations. We further assume that the reader has a basic idea of planes.
If you are uncomfortable with these concepts, please follow this link to learn more about vector operations and planes.
Let’s look at some formal mathematical facts and details to find an intersection of rays in the plane.
We know the dot product of two orthogonal or perpendicular vectors is always 0. Suppose a and b are two perpendicular or orthogonal vectors, then a.b=0.
Now, consider point p0 on the plane representing the distance of the plane from the origin and vector n, normal to the plane. We can compute vector p by subtracting any point on the plane from point p0.
The resultant vector lies in the plane and perpendicular to the plane normal.
This fact gives us the equation:
(p-p0) . N = 0 (i)
Consider l0 and l as the ray’s starting point and the ray’s direction, respectively. We can reach the plane (means intersect at point p) using the parametric form:
l0 + l * t = p (ii)
If the ray isn’t parallel to the plane, then t times in the direction of the ray; the ray will interest the plane. Next, we can put the value of p from equation (ii) into the equation (i), and we get:
(l0 + l * t – p0) . n = 0
We want to compute t, which can help us to compute the position of the intersection point using the parametric equation. Let’s solve the equation:
l * t . n + (l0 – p0) . n = 0
l * t . n = - (l0 – p0) . n
t = - (l0 – p0) . n / l . n
We need to determine the dot product of the line and plane normal coming in the denominator because if the value is 0 or close to 0, the result will be infinite, which means the ray and plane are parallel. Therefore, we will check the denominator and return false, which means there is no solution (that is, the intersection point).
The following function can check whether the ray is intersecting the plane or not:
bool intersection plane(Vec3 &n, Vec3 &p0, Vec3 &lo, Vec3 &l, float &t) {
// assuming vectors are all normalized
float denom = n.dot(l);
if (denom < 0.00005 && denom > -0.00005) / denom is near to 0 return false;
Vec3 p010 = p0 – l0;
t = p010.dot(n);
if t>=0) return true;
return false;
}
Using the if condition, we check if the denominator is close to 0, which means there will be infinite results. The dot product of plane normal with ray direction gives 0 if the ray is parallel to the plane; therefore, we return false.
Otherwise, we will calculate t and return true, meaning the ray intersects the plane. We can find a point of intersection using t.
Next, we have below the complete code to find the intersection point of a ray and a plane.
#include <math.h>
#include <iostream>
using namespace std;
class Vec3 {
float x;
float y;
float z;
public:
Vec3() { x, y, z = 0; }
Vec3(float x, float y, float z) {
this->x = x;
this->y = y;
this->z = z;
}
Vec3 &operator+=(const Vec3 &b) {
x = x + b.x;
y = y + b.y;
z = z + b.z;
return *this;
}
Vec3 &operator+=(const float s) {
x = x + s;
y = y + s;
z = z + s;
return *this;
}
Vec3 operator+(const float s) {
Vec3 newV = *this;
return newV += s;
}
Vec3 &operator-=(const Vec3 &b) {
x = x - b.x;
y = y - b.y;
z = z - b.z;
return *this;
}
Vec3 operator+(const Vec3 &b) {
Vec3 newV = *this;
newV += b;
return newV;
}
Vec3 operator-(const Vec3 &b) {
Vec3 newV = *this;
newV -= b;
return newV;
}
Vec3 operator*(const Vec3 &b) { // cross operator
Vec3 newV;
newV.x = y * b.z - z * b.y;
newV.y = x * b.z - z * b.x;
newV.z = x * b.y - y * b.x;
return newV;
}
Vec3 &operator*(const float s) { // Dot equal operator
x = x * s;
y = y * s;
z = z * s;
return *this;
}
float mag() { return sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2)); }
Vec3 &normalize() {
double mag = this->mag();
x /= mag;
y /= mag;
z /= mag;
}
// Dot operator
float dot(const Vec3 &b) { return x * b.x + y * b.y + z * b.z; }
friend ostream &operator<<(ostream &out, const Vec3 &v) {
out << "(" << v.x << ", " << v.y << ", " << v.z << ")\n";
return out;
}
};
bool intersectPlane(Vec3 &n, Vec3 &p0, Vec3 &l0, Vec3 &l, float &t) {
// considering vectors are normalized
float denom = n.dot(l); // dot product n.l
if (denom > 1e-6) {
Vec3 p0l0 = p0 - l0;
t = n.dot(p0l0) / denom;
return (t >= 0);
}
return false;
}
int main() {
Vec3 n1(3, -9, 1); // Normal
Vec3 p01(-4, 2, 2);
Vec3 l01(1, 1, 1);
Vec3 l1(1, 2, 1);
float t;
n1.normalize();
p01.normalize();
l01.normalize();
l1.normalize();
if (intersectPlane(n1, p01, l01, l1, t))
cout << "T:" << t << '\n';
else
cout << "Ray is not intersecting the plane\n";
cout << "------------------------\n";
Vec3 n2(2, 2, -2); // Normal
Vec3 p02(2, 4, 1);
Vec3 l02(1, 1, 1);
Vec3 l2(1, 2, 1);
n2.normalize();
p02.normalize();
l02.normalize();
l2.normalize();
if (intersectPlane(n2, p02, l02, l2, t))
cout << "T:" << t << '\n';
else
cout << "Ray is not intersecting the plane\n";
Vec3 intersectionPoint = l02 + l2 * t;
cout << intersectionPoint;
return 0;
}
The output of this code is:
Ray is not intersecting the plane
------------------------
T:0.629199
(0.83422, 1.09109, 0.83422)
As you can see, in our first data point, the ray is parallel to the plane. Therefore, the denominator is 0.
However, in the second set of data points, the denominator is greater than 0, so we can find the intersection point.