Note on “Ray tracing in one week” - Part 2

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

right now you should have the code from part 1 (and enhancement from part -1), and it should be able to render a gradient background. In this part, we will add a sphere to the scene.

Sphere

From math, we know that a sphere is defined by the equation:

$$x^2 + y^2 + z^2 = r^2$$

Now it’s only 4 variables. Ez right? Well, if we consider a sphere at a point , then the equation becomes:

$$(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2$$

Now it’s 7 variables and the equation starts to become ugly. But we can simplify it by using linear algebra, we can define two vectors $\mathbf{P}$ and $\mathbf{C}$, such that:

Therefore, the equation simply becomes the dot product:

$$(\mathbf{P} - \mathbf{C}) \cdot (\mathbf{P} - \mathbf{C}) = r^2$$

Intersection

We defined $\mathbf{P}$ as the points on the sphere. Therefore, if there exist point(s) on a ray, such that the point(s) satisfy the equation above, the ray intersects with the sphere.

We used to define ray $\mathbf{P}$ as:

$$\mathbf{P} = \mathbf{A} + t\mathbf{B}$$

Therefore, we can substitute $\mathbf{P}$ into the equation above:

If we do some moving around, we can get:

$$(t\mathbf{B} + (\mathbf{A} - \mathbf{C})) \cdot (t\mathbf{B} + (\mathbf{A} - \mathbf{C})) = r^2$$

Distribute, we get:

$$t^2\mathbf{B} \cdot \mathbf{B} + 2t\mathbf{B} \cdot (\mathbf{A} - \mathbf{C}) + (\mathbf{A} - \mathbf{C}) \cdot (\mathbf{A} - \mathbf{C}) = r^2$$

This equation seems scary with matrix multiplication, but actually it’s not. $\mathbf{B} \cdot \mathbf{B}$ and other dot products are just numbers, and we have an expression of $r$ with respect to $t$. Since there is a $t^2$ in this equation, we can use the quadratic formula to solve for $t$:

where constants $a$, $b$, and $c$ being:

But, actually, we don’t even need to solve it right now. We can just use the discriminant (or $b^2 - 4ac$) to determine if the ray intersects with the sphere or not. Recall the algebra class, negative discriminant means no real solution, zero discriminant means one real solution, and positive discriminant means two real solutions. Therefore, we can use the discriminant to determine if the ray intersects with the sphere or not.

Implementation

To detect whether a ray hits the sphere, we need to have a function that checks intersection:

bool sphereIntersectionCheck(double r, const MathUtil::Point3& pos, const MathUtil::Ray& ray) {
    double a = ray.dir().dot(ray.dir());
    double b = 2 * ray.dir().dot(ray.pos() - pos);
    double c = (ray.pos() - pos).dot(ray.pos() - pos) - r * r;
    return (b * b - 4 * a * c >= 0);
}

In this function there is a bug, can you find it?

and then we need to add it to our ray->color mapping function:

ImageUtil::Color rayColor(const MathUtil::Ray& ray) {
    auto unit = ray.dir().unit();
    ImageUtil::Color color;
    if (!sphereIntersectionCheck(2, MathUtil::Point3(0, 0, -5), ray)) {
        color = ImageUtil::Color(std::min(0.5 * unit.x + 1, 1.0), std::min(0.5 * unit.y + 1, 1.0), 1);
    }
    else {
        color = ImageUtil::Color(0.67, 0.84, 0.9);
    }
    return color;
}

run, and now we get:

Great! But… it lacked a looooot of details, like shading, reflection, blah blah blah. But it’s a good start, since we are actually drawing a sphere with ray tracing! And it’s time to make some shading.

Surface normal

A normal vector is a vector that is orthogonal to the plane(or sphere) at that point. When finding the normal vector of one point, we have to normalize it. However, when doing normal vector calculation, we don’t have to use square root, since the sphere normal vector can be obtained by simply dividing the vector by the radius.

For a sphere, let’s denote the vector from the camera to the hitpoint as $\mathbf{R}$, the position vector of the sphere as $\mathbf{C}$, and the vector from the center of the sphere to the hitpoint as $\mathbf{N}$. Then, we can get the normal vector by:

$$\mathbf{N} = \mathbf{R} - \mathbf{C}$$

to normalize it, we can simply divide it by the radius of the sphere:

$$\mathbf{N'} = \frac{\mathbf{R} - \mathbf{C}}{r}$$

Since we already normalized the normal vector, we can simply map the color value (r, g, b) to the component of the normal vector, which is a simple way to represent shading.

Implementation

We need to modify our sphereIntersectionCheck a little bit, since we need to return the solution to the equation to find the hit point location (calculated from our ray).

double sphereIntersectionCheck(double r, const MathUtil::Point3& pos, const MathUtil::Ray& ray) {
    auto oc = ray.pos() - pos;
    double a = ray.dir().dot(ray.dir());
    double b = 2 * oc.dot(ray.dir());
    double c = oc.dot(oc) - r * r;
    double discriminant = b * b - 4 * a * c;
    if (discriminant < 0) {
        return -1.0;
    }
    return (-b - std::sqrt(discriminant)) / (2.0 * a);
}

In this modified function, we returned the solution to the equation, which will be used to calculate the hit point position.

Now, we apply this function to the ray->color mapping function:

ImageUtil::Color rayColor(const MathUtil::Ray& ray) {
    auto t = sphereIntersectionCheck(2, MathUtil::Point3(0, 0, -5), ray);
    ImageUtil::Color color;
    if (t <= 0.0) {
        auto unit = ray.dir().unit();
        color = ImageUtil::Color(std::min(0.5 * unit.x + 1, 1.0), std::min(0.5 * unit.y + 1, 1.0), 1);
    }
    else {
        auto normal = (ray.at(t) - MathUtil::Point3(0, 0, -5)).unit();
        color = {normal.x + 1, normal.y + 1, normal.z + 1};
        color *= 0.5;
    }
    return color;
}

And now, we have…

a shaded sphere!

Simplification of code

We can simplify the code a bit by using some math, let’s take a look at the sphereIntersectionCheck function:

double sphereIntersectionCheck(double r, const MathUtil::Point3& pos, const MathUtil::Ray& ray) {
    auto oc = ray.pos() - pos;
    double a = ray.dir().dot(ray.dir());
    double b = 2 * oc.dot(ray.dir());
    double c = oc.dot(oc) - r * r;
    double discriminant = b * b - 4 * a * c;
    if (discriminant < 0) {
        return -1.0;
    }
    return (-b - std::sqrt(discriminant)) / (2.0 * a);
}

We can see there are several self dot products, and recall from linear algebra, that the dot product of a vector with itself is the square of the length of the vector.

Also, we can take a look at the quadratic formula. If we substitute $2h = b$, we can get:

which simplifies calculation.

Therefore, we can simplify the function to:

double sphereIntersectionCheck(double r, const MathUtil::Point3& pos, const MathUtil::Ray& ray) {
    auto oc = ray.pos() - pos;
    double a = ray.pos().lengthSq();
    double half_b = oc.dot(ray.dir());
    double c = oc.lengthSq() - r * r;
    double discriminant = half_b * half_b - a * c;
    if (discriminant < 0) {
        return -1.0;
    }
    return (-half_b - std::sqrt(discriminant)) / (2.0 * a);
}

Associated changed files

main.cpp

#include "ImageUtil.h"
#include "MathUtil.h"
#include "GraphicObjects.h"
#include <vector>
#include <cmath>



double sphereIntersectionCheck(double r, const MathUtil::Point3& pos, const MathUtil::Ray& ray) {
    auto oc = ray.pos() - pos;
    double a = ray.pos().lengthSq();
    double half_b = oc.dot(ray.dir());
    double c = oc.lengthSq() - r * r;
    double discriminant = half_b * half_b - a * c;
    if (discriminant < 0) {
        return -1.0;
    }
    return (-half_b - std::sqrt(discriminant)) / (2.0 * a);
}

ImageUtil::Color rayColor(const MathUtil::Ray& ray) {
    auto t = sphereIntersectionCheck(2, MathUtil::Point3(0, 0, -5), ray);
    ImageUtil::Color color;
    if (t <= 0.0) {
        auto unit = ray.dir().unit();
        color = ImageUtil::Color(std::min(0.5 * unit.x + 1, 1.0), std::min(0.5 * unit.y + 1, 1.0), 1);
    }
    else {
        auto normal = (ray.at(t) - MathUtil::Point3(0, 0, -5)).unit();
        color = {normal.x + 1, normal.y + 1, normal.z + 1};
        color *= 0.5;
    }
    return color;
}


int main() {
    auto camera = Camera(1920, 16.0 / 9.0, 2.0, 1, MathUtil::Point3(0, 0, 0));
    auto img = std::vector<std::vector<ImageUtil::Color>>();
    for (int i = 0; i < camera.getHeight(); ++i) {
        auto v = std::vector<ImageUtil::Color>();
        if (i % 10 == 0 || i > camera.getHeight() - 10) spdlog::info("line remaining: {}", (camera.getHeight() - i + 1));
        for (int j = 0; j < camera.getWidth(); ++j) {
            auto ray_dir =  camera.getPixRayDir(j, i);
            auto ray = MathUtil::Ray(camera.getPosition(), ray_dir);
            auto color = rayColor(ray);
            v.emplace_back(color);
        }
        img.emplace_back(v);
    }
    ImageUtil::makePPM(camera.getWidth(), camera.getHeight(), img, "out", "test.ppm");
    return 0;
}

Associated files

directory layout

.
|-- CMakeLists.txt
|-- README.md
|-- include
|   `-- OneWeekend
|       |-- GraphicObjects.h
|       |-- ImageUtil.h
|       `-- MathUtil.h
|-- src
|   |-- GraphicObjects.cpp
|   |-- ImageUtil.cpp
|   |-- MathUtil.cpp
|   `-- main.cpp
`-- thirdparties
    |-- include
    |   `-- spdlog  [spdlog include omitted]
    `-- lib
        `-- Darwin
            `-- libspdlog.a