Note on "Ray Tracing The Next Week" - Part 5

Note on "Ray Tracing The Next Week" - Part 5

Ayano Kagurazaka Lv3
  2024-02-21 2024-02-21 Created   2025-04-28 21:40:37 2025-04-28 21:40:37 Updated      3k Words  18 Mins  

The code is compiled with cmake, please check out CMakeLists.txt

icon-padding

I made some changes to the code, please check out this link

# Quadrilaterals

A quadrilateral(quad for short) in R3\mathbb{R}^3 can be defined by 3 vectors in R3\mathbb{R}^3, Q,u,v\mathbf{Q, u, v}, as:

P(s,t)=Q+su+tvP(s, t) = \mathbf{Q} + s\mathbf{u} + t\mathbf{v}

Which Q\mathbf{Q} indicates one vertex of the quad, and the other two vectors spans a plane. If we constrain the range of ss and tt to be a finite set AA and BB, the set {P(s,t):s∈A,t∈B}\{P(s, t) : s\in A,t\in B\} contains all the points on the quad.

example quad

If we recall the way we generate AABB, we will notice a problem: when we generate AABB for a plane:

MathUtil.cpp 
1
2
3
4
5
AABB::AABB(const Point3 &a, const Point3 &b) {
	x = Interval(fmin(a[0], b[0]), fmax(a[0], b[0]));
	y = Interval(fmin(a[1], b[1]), fmax(a[1], b[1]));
	z = Interval(fmin(a[2], b[2]), fmax(a[2], b[2]));
}

If the plane lies on XY, YZ, or XZ plane, the AABB will have thickness of 0 and will result in numeric error in our intersection function:

MathUtil.cpp 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
bool AABB::hit(const Ray &r, Interval ray_int) const {
	for (int i = 0; i < 3; i++) {
		auto invD = 1.0 / r.dir()[i];
		auto orig = r.pos()[i];

		auto t0 = (axis(i).min - orig) * invD;
		auto t1 = (axis(i).max - orig) * invD;

		if (invD < 0.0) {
			std::swap(t0, t1);
		}

		if (t0 > ray_int.min) ray_int.min = t0;
		if (t1 < ray_int.max) ray_int.max = t1;

		if (ray_int.max < ray_int.min) {
			return false;
		}
	}
	return true;
}

To solve this, we can add some padding to our plane, to do this, we first need to be able to expand our Interval:

MathUtil.h 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
class Interval {
public:
	double min, max;
	Interval();

	Interval(double min, double max);

	Interval(const Interval &first, const Interval &second);

	bool within(double x) const;

	bool surround(double x) const;

	double clamp(double x) const;

	Interval expand(double delta);

	static const Interval empty, universe;
};
MathUtil.cpp 
1
2
3
4
Interval Interval::expand(double delta) {
	auto padding = delta / 2;
	return Interval(min - padding, max + padding);
}

Then we implement a padding function in AABB, which returns a padded AABB.

MathUtil.h 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
class AABB {
public:
	Interval x, y, z;

	AABB() = default;

	AABB(const Point3 &a, const Point3 &b);

	AABB(const Interval &x, const Interval &y, const Interval &z);

	AABB(const AABB &up, const AABB &down);

	Interval axis(int i) const;

	AABB pad();

	[[nodiscard]] bool hit(const Ray &r, Interval ray_int) const;
};
MathUtil.cpp 
1
2
3
4
5
6
7
8
9
10
AABB AABB::pad() {
	double delta = 0.0001;
	double size_x = x.max - x.min;
	double size_y = y.max - y.min;
	double size_z = z.max - z.min;
	Interval new_x = (size_x >= delta ? x : x.expand(delta));
	Interval new_y = (size_y >= delta ? y : y.expand(delta));
	Interval new_z = (size_z >= delta ? z : z.expand(delta));
	return AABB(new_x, new_y, new_z);
}

# Intersection

Now we have to find a way to check whether the ray hits the quad. This can be broken down into three steps:

  1. Find the plane containing the quad.
  2. Solve for the intersection.
  3. Check the intersection lies in the quad or not.

# Plane intersection

We know the scalar form of a plane is:

Ax+By+Cz=DAx+By+Cz=D

where

[ABC]\begin{bmatrix} A \\ B \\ C \end{bmatrix}

is a vector that is orthogonal to the plane.

By observation, we see this is just a dot product of the normal vector of the plane and the position vector of any point on the plane. The point can be the intersection point of our ray and the plane. Therefore, denote our ray as R(t)=P+td\mathbf{R}(t) = \mathbf{P} + t\mathbf{d} we can rewrite the equation as:

n⋅v=Dn⋅(P+td)=D\begin{align*} \mathbf{n}\cdot\mathbf{v} &= D \\ \mathbf{n}\cdot(\mathbf{P} + t\mathbf{d}) &= D \end{align*}

Do some simplification:

n⋅(P+td)=Dn⋅P+t⋅n⋅d=Dt=D−n⋅Pn⋅d\begin{align*} \mathbf{n}\cdot(\mathbf{P} + t\mathbf{d}) &= D \\ \mathbf{n}\cdot\mathbf{P} + t \cdot \mathbf{n}\cdot\mathbf{d} &= D \\ t = \frac{D - \mathbf{n}\cdot\mathbf{P}}{\mathbf{n}\cdot\mathbf{d}} \end{align*}

From this equation, we can see that if n⋅d=0\mathbf{n}\cdot\mathbf{d} = 0, we don’t have an intersection.

On the steps above we assumed we know the equation of the plane, but it’s quite simple to obtain. We know Q\mathbf{Q} must be on the plane, and normalized u×v\mathbf{u} \times \mathbf{v} is the normal vector. Therefore, D=Q⋅(u×v)D = \mathbf{Q} \cdot (\mathbf{u} \times \mathbf{v}), and we obtained the scalar equation with the information known.

# Quad intersection

When doing the quad intersection, since we already know the intersection point is on the plane of quad lying on, we reduced the problem to happen on R2\mathbb{R^2}. There are two vectors, and we know 2 linear independent vectors can span the whole R2\mathbb{R}^2, we can treat u\mathbf{u} and v\mathbf{v} as basis of the space.

We know, for a plane defined by:

P=Q+αu+βvP=\mathbf{Q} + \alpha\mathbf{u} + \beta\mathbf{v}

We can find any point on that plane by finding its α\alpha and β\beta. And if we want to see if the point lays inside the quad created by u\mathbf{u} and v\mathbf{v}, we just need to check the vector originated from Q\mathbf{Q} to the hit point has both 0≤α≤10 \leq \alpha \leq 1 and 0≤β≤10 \leq \beta \leq 1.

Denote p\mathbf{p} as the vector on the uv plane, originated from the origin to the hit point, we have:

p=P−Q=αu+βv\begin{align*} \mathbf{p} &= \mathbf{P} - \mathbf{Q}\\ &= \alpha\mathbf{u} + \beta\mathbf{v} \end{align*}

quadIntersection

to solve for α\alpha and β\beta, we cross each side by u\mathbf{u} and v\mathbf{v}:

u×p=u×(αu+βv)=u×αu+u×βv=u×βvv×p=v×(αu+βv)=v×αu+v×βv=v×αu\begin{align*} \mathbf{u} \times\mathbf{p} &= \mathbf{u} \times (\alpha \mathbf{u} + \beta \mathbf{v}) \\ &=\mathbf{u}\times\alpha\mathbf{u} + \mathbf{u}\times\beta\mathbf{v}\\ &=\mathbf{u}\times\beta\mathbf{v} \\ \mathbf{v} \times\mathbf{p} &= \mathbf{v} \times (\alpha \mathbf{u} + \beta \mathbf{v}) \\ &=\mathbf{v}\times\alpha\mathbf{u} + \mathbf{v}\times\beta\mathbf{v}\\ &=\mathbf{v}\times\alpha\mathbf{u} \end{align*}

both side dot by n\mathbf{n}, the normal vector, and move the terms to one side of the equal sign, we get:

α=n⋅(v×p)n⋅(v×u)β=n⋅(u×p)n⋅(u×v)\begin{align*} \alpha = \frac{\mathbf{n} \cdot (\mathbf{v}\times\mathbf{p})}{\mathbf{n} \cdot (\mathbf{v}\times\mathbf{u})} \\ \beta = \frac{\mathbf{n} \cdot (\mathbf{u}\times\mathbf{p})}{\mathbf{n} \cdot (\mathbf{u}\times\mathbf{v})} \\ \end{align*}

To reduce the complexity of computation, we flip the cross order in the expression with alpha:

α=n⋅(p×v)n⋅(u×v)\alpha = \frac{\mathbf{n}\cdot(\mathbf{p}\times\mathbf{v})}{\mathbf{n}\cdot(\mathbf{u}\times\mathbf{v})}

And extract a common constant from each value, denote it as w\mathbf{w}

w=nn⋅(u×v)=nn⋅n\begin{align*} \mathbf{w} &= \frac{\mathbf{n} }{\mathbf{n} \cdot (\mathbf{u}\times\mathbf{v})}\\ &=\frac{\mathbf{n} }{\mathbf{n}\cdot\mathbf{n}} \end{align*}

and we get the simplified version of α\alpha and β\beta:

α=w⋅(p×v)β=w⋅(u×p)\begin{align*} \alpha = \mathbf{w} \cdot (\mathbf{p} \times \mathbf{v}) \\ \beta = \mathbf{w} \cdot (\mathbf{u} \times \mathbf{p}) \\ \end{align*}

# Implementation

First we have to implement a quad class to store those values:

GraphicObjects.h 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
class Quad : public IHittable {
public:
	Quad(const AppleMath::Vector3& Q, const AppleMath::Vector3& u, const AppleMath::Vector3& v, std::shared_ptr<IMaterial> mat);

	virtual ~Quad() = default;

	AABB boundingBox() const override;

	virtual void setBoundingBox();

	bool hit(const Ray& r, Interval interval, HitRecord& record) const override;

	bool inside(double a, double b, HitRecord& rec) const;

private:
	AppleMath::Vector3 Q, u, v;
	AppleMath::Vector3 normal;
	double D;
	AppleMath::Vector3 w;
	std::shared_ptr<IMaterial> mat;
	AABB bbox;

};

In the constructor, we initialize some of invariant fields of this plane:

GraphicObjects.cpp

1
2
3
4
5
6
7
8
9
10
Quad::Quad(const AppleMath::Vector3 &Q, const AppleMath::Vector3 &u, const AppleMath::Vector3 &v,
		   std::shared_ptr<IMaterial> mat) :
	Q(Q),
	u(u), v(v), mat(mat) {
	auto n = u.cross(v);
	normal = n.normalized();
	D = normal.dot(Q);
	w = n / n.dot(n);
	setBoundingBox();
}

When setting the bounding box, to avoid zero-thickness error, we do a padding around AABB:

GraphicObjects.cpp

1
void Quad::setBoundingBox() { bbox = AABB(Q, Q + u + v).pad(); }

and AABB getter:

GraphicObjects.cpp

1
AABB Quad::boundingBox() const { return bbox; }

Then, in our hit function, we first check whether the ray hits the plane, and then check to see whether the solved t is valid:

GraphicObjects.cpp

1
2
3
4
5
6
7
8
9
10
bool Quad::hit(const Ray &r, Interval interval, HitRecord &record) const {
	double denom = normal.dot(r.dir());
	if (fabs(denom) < 1e-8) {
		return false;
	}
	auto t = (D - normal.dot(r.pos())) / denom;
	if (!interval.surround(t)) {
		return false;
	}
}

Then, we check the intersection, and find out if the hit point is inside given quad:

GraphicObjects.cpp

1
2
3
4
5
6
7
8
9
10
11
12
bool Quad::hit(const Ray &r, Interval interval, HitRecord &record) const {

	// ...

	auto intersection = r.at(t);
	auto plane_hit_vect = intersection - Q;
	auto alpha = w.dot(plane_hit_vect.cross(v));
	auto beta = w.dot(u.cross(plane_hit_vect));
	if (!inside(alpha, beta, record)) {
		return false;
	}
}

and the inside function:

GraphicObjects.cpp

1
2
3
4
5
6
7
bool Quad::inside(double a, double b, HitRecord &rec) const {
	if (a < 0 || a > 1 || b < 0 || b > 1)
		return false;
	rec.u = a;
	rec.v = b;
	return true;
}

Finally, we initialize the record fields with the calculated values:

GraphicObjects.cpp

1
2
3
4
5
6
7
8
9
10
11
bool Quad::hit(const Ray &r, Interval interval, HitRecord &record) const {

	// ...

	record.t = t;
	record.p = intersection;
	record.material = mat;
	record.setFaceNormal(r, normal);

	return true;
}

Full code:

GraphicObjects.cpp 
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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
AABB Quad::boundingBox() const { return bbox; }

void Quad::setBoundingBox() { bbox = AABB(Q, Q + u + v).pad(); }

Quad::Quad(const AppleMath::Vector3 &Q, const AppleMath::Vector3 &u, const AppleMath::Vector3 &v,
		   std::shared_ptr<IMaterial> mat) :
	Q(Q),
	u(u), v(v), mat(mat) {
	auto n = u.cross(v);
	normal = n.normalized();
	D = normal.dot(Q);
	w = n / n.dot(n);
	setBoundingBox();
}

bool Quad::inside(double a, double b, HitRecord &rec) const {
	if (a < 0 || a > 1 || b < 0 || b > 1)
		return false;
	rec.u = a;
	rec.v = b;
	return true;
}


bool Quad::hit(const Ray &r, Interval interval, HitRecord &record) const {
	double denom = normal.dot(r.dir());
	if (fabs(denom) < 1e-8) {
		return false;
	}
	auto t = (D - normal.dot(r.pos())) / denom;
	if (!interval.surround(t)) {
		return false;
	}
	auto intersection = r.at(t);
	auto plane_hit_vect = intersection - Q;
	auto alpha = w.dot(plane_hit_vect.cross(v));
	auto beta = w.dot(u.cross(plane_hit_vect));
	if (!inside(alpha, beta, record)) {
		return false;
	}
	record.t = t;
	record.p = intersection;
	record.material = mat;
	record.setFaceNormal(r, normal);

	return true;
}

To test things out, we write a sample scene:

scenes.h

1
2
3
4
5
//...

void quads();

//...

scenes.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
void quads() {
	HittableList world;

	auto blue = std::make_shared<Lambertian>(Color{0.36, 0.81, 0.98});
	auto pink = std::make_shared<Lambertian>(Color{0.96, 0.66, 0.72});
	auto white = std::make_shared<Lambertian>(Color{1, 1, 1});

	world.add(std::make_shared<Quad>(Point3{-3,-2, 5}, AppleMath::Vector3{0, 0,-4}, AppleMath::Vector3{0, 4, 0}, blue));
    world.add(std::make_shared<Quad>(Point3{-2,-2, 0}, AppleMath::Vector3{4, 0, 0}, AppleMath::Vector3{0, 4, 0}, white));
    world.add(std::make_shared<Quad>(Point3{ 3,-2, 1}, AppleMath::Vector3{0, 0, 4}, AppleMath::Vector3{0, 4, 0}, pink));
    world.add(std::make_shared<Quad>(Point3{-2, 3, 1}, AppleMath::Vector3{4, 0, 0}, AppleMath::Vector3{0, 0, 4}, blue));
    world.add(std::make_shared<Quad>(Point3{-2,-3, 5}, AppleMath::Vector3{4, 0, 0}, AppleMath::Vector3{0, 0,-4}, pink));
	Camera camera(1920, 9.0 / 9.0, 90, {0, 0, 9}, {0, 0, 0}, 0);
	camera.setSampleCount(100);
	camera.setShutterSpeed(1.0 / 24.0);
	camera.setRenderDepth(50);
	camera.setRenderThreadCount(12);
	camera.setChunkDimension(64);
	render(world, camera, "quads.ppm");
}

In our main function, we can test that out:

main.cpp

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
30
31
32
#include "scenes.h"

int main() {

	switch (6) {
		case 0:
			randomSpheres();
			break;
		case 1:
			twoSpheres();
			break;
		case 2:
			huajiSphere();
			break;
		case 3:
			perlinSpheres();
			break;
		case 4:
			terrain();
			break;
		case 5:
			rotationTest();
			break;
		case 6:
			quads();
			break;
		default:
			break;
	}
	return 0;
}

The result:

quads

# triangles

icon-padding

This part is not in The next week

For triangles, we can still use the plain-interior way to check intersection. Therefore, we only need to think of a way to check whether the point is inside the triangle or not.

If we connect the vertex of a triangle with the intersection point, we get a vector pointing towards the intersection point, and there are only two situations we need to consider:

  • The vertex-intersection vector is in the left of the vertex-(next vertex) vector
  • The vertex-intersection vector is in the right of the vertex-(next vertex) vector, or they are in the same direction

We can use cross product to differentiate the first two cases: since they have different directions, the cross product must have same magnitude but opposite direction. Therefore, we only need to check its relative direction to the normal vector, which can be done by checking the sign of dot product.

Similar to quad, our triangle class has similar definition:

GraphicObjects.h 
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
class Triangle : public IHittable {
public:

	Triangle(const AppleMath::Vector3& Q, const AppleMath::Vector3& u, const AppleMath::Vector3& v, std::shared_ptr<IMaterial> mat);

	virtual ~Triangle() = default;

	virtual void setBoundingBox();

	AABB boundingBox() const override;

	bool hit(const Ray& r, Interval interval, HitRecord& record) const override;

	bool inside(const AppleMath::Vector3& intersection) const;

private:

	AppleMath::Vector3 Q, u, v;
	AppleMath::Vector3 normal;
	AppleMath::Vector3 w;
	double D;

	std::shared_ptr<IMaterial> mat;
	AABB box;
};

and similar constructor and bounding box definition:

(code below are inside GraphicObjects.cpp )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
Triangle::Triangle(const AppleMath::Vector3 &Q, const AppleMath::Vector3 &u, const AppleMath::Vector3 &v,
				   std::shared_ptr<IMaterial> mat) :
	Q(Q),
	u(u), v(v), mat(mat) {
	auto n = v.cross(u);
	normal = n.normalized();
	w = n / n.dot(n);
	setBoundingBox();
}

void Triangle::setBoundingBox() { box = AABB(Q, Q + u + v).pad(); }

AABB Triangle::boundingBox() const {
	return box;
}

We also have to create a helper function to detect containment:

1
2
3
4
5
6
7
8
9
10
bool Triangle::inside(const AppleMath::Vector3 &intersection) const {
	auto side1 = v;
	auto side2 = u - v;
	auto side3 = -u;
	auto dis1 = intersection - Q;
	auto dis2 = intersection - (Q + side1);
	auto dis3 = intersection - (Q + side1 + side2);
	return (side1.cross(dis1).dot(normal) > 0 && side2.cross(dis2).dot(normal) > 0 &&
			side3.cross(dis3).dot(normal) > 0);
}

and a hit function which is similar to the hit function of quad:

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
bool Triangle::hit(const Ray& r, Interval interval, HitRecord& record) const {
	double denom = normal.dot(r.dir());
	if (fabs(denom) < 1e-8) {
		return false;
	}

	auto t = (D - normal.dot(r.pos())) / denom;

	if (!interval.surround(t)) {
		return false;
	}

	auto intersection = r.at(t);
	if (!inside(intersection)) {
		return false;
	}
	auto plane_hit_vec = intersection - Q;
	auto alpha = w.dot(plane_hit_vec.cross(u));
	auto beta = w.dot(u.cross(plane_hit_vec));
	record.t = t;
	record.p = intersection;
	record.material = mat;
	record.setFaceNormal(r, normal);
	record.u = alpha;
	record.v = beta;
	return true;
}

Then, we make a demo scene to test the triangle out:

scenes.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#ifndef RAYTRACING_SCENES_H
#define RAYTRACING_SCENES_H



void randomSpheres();

void twoSpheres();

void huajiSphere();

void perlinSpheres();

void terrain();

void rotationTest();

void quads();

void triangles();

#endif // RAYTRACING_SCENES_H

scene.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void triangles() {
	HittableList world;

	auto blue = std::make_shared<Lambertian>(Color{0.36, 0.81, 0.98});
	auto pink = std::make_shared<Lambertian>(Color{0.96, 0.66, 0.72});
	auto white = std::make_shared<Lambertian>(Color{1, 1, 1});

	Point3 p1{1.5, 1.5, 0};
	Point3 p2{-1.5, -1.5, 0};
	AppleMath::Vector3 vertical{3, 0, 0};
	AppleMath::Vector3 horizontal{0, 3, 0};


	world.add(std::make_shared<Triangle>(p1, -vertical, -horizontal, blue));
	world.add(std::make_shared<Triangle>(p2, vertical, horizontal, pink));
	Camera camera(1920, 9.0 / 9.0, 90, {0, 0, 5}, {0, 0, 0}, 0);
	camera.setSampleCount(100);
	camera.setShutterSpeed(1.0 / 24.0);
	camera.setRenderDepth(50);
	camera.setRenderThreadCount(12);
	camera.setChunkDimension(64);
	render(world, camera, "triangle.ppm");
}

The result:

triangle

# Associated code:

The complete code till this point can be found in this link

Recommended reading~
Note on Note on "Ray tracing in one week" - Part 2 Note on Note on "Ray tracing in one week" - Part 7 Note on Note on "Ray Tracing The Next Week" - Part 3
Recommended reading~
Note on Note on "Ray tracing in one week" - Part 2 Note on Note on "Ray tracing in one week" - Part 7
Comments
On this page
Note on "Ray Tracing The Next Week" - Part 5
  1. Quadrilaterals
    1. Intersection
  2. triangles
  3. Associated code: