Ray tracing is a rendering technique, which traces paths of lights for pixels in an image plane. Despite its expensive computational cost, it is widely used because it can produce realistic images. In this project, we render various three-dimensional fractal images via ray tracing. Since directly calculating an intersection between a ray and a fractal object is a nontrivial problem, we use distance estimation and adaptive ray marching to find intersections. Normal vectors at intersection points are estimated via the finite difference method rather than being calculated geometrically. From calculated intersections and normal vectors, we can apply multiple lighting techniques to make images more realistic. In this project, we adopt the Phong lighting, reflection, refraction, soft shadows, and ambient occlusion. In addition, we calculate multiple camera positions based on Bezier curves and create videos by combining multiple rendered images.
Theoretical backgrounds of our core goals are given by paper Ray Tracing Deterministic 3-D Fractals of John C. Hart, Daniel J. Sandin, and Louis H. Kauffman.
Hart et al. suggested the adaptive ray marching technique. It is useful for rendering objects, whose boundary is not directly known, like three-dimensional fractals. Instead of computing an intersection at once by solving equations, a ray advances multiple times until the ray becomes close enough to an object. For this purpose, we should define a distance estimation function. The function calculates the minimum distance to the object from a given point. By advancing no more than the return value of the function, the ray avoids passing through the object rather than stopping at an intersection.
The below intersect function judges whether a given ray intersects with an object in a scene. The de function is a distance estimator.
bool intersect(Ray ray) {
float distance = 0;
int steps;
for (steps = 0; steps < max_ray_steps; steps++) {
float d = de(ray(distance));
distance += d;
if (d < min_distance) break;
}
return steps != max_ray_steps;
}The code is simplified for conciseness and one can find real implementation in the Scene.cpp file.
As Hart et al. noted, multiple methods to calculate normal vectors at an intersection point exist. We use the finite difference method, which is easy to implement and precise enough. Since a normal vector at an intersection point has a direction, which increases the distance to the object most fast, we can estimate the normal vector as the gradient vector of a distance estimator function at the intersection. We cannot compute the derivative of each component so that we use the finite difference method to estimate the gradient vector numerically.
The below estimate_normal function estimate the normal vector at a given point.
vec3 estimate_normal(vec3 point) {
vec3 x(normal_distance, 0, 0), y(0, normal_distance, 0), z(0, 0, normal_distance);
return normalize(vec3(
de(point + x) - de(point - x),
de(point + y) - de(point - y),
de(point + z) - de(point - z)));
}The code is simplified for conciseness and one can find real implementation in the Scene.cpp file.
A Julia quaternion fractal contains every quaternion \(z\) such that \(f^n(z)\) does not diverges as \(n\) goes to the infinity where \(f(z)=z^2+d\), superscript \(n\) denotes \(n\)-fold application of the function, and \(d\) is an arbitrary fixed quaternion. Hart et al. noted that we can use \(d(z)=\frac{|f^n(z)|}{2|f'^n(z)|}log|f^n(z)|\) where \(n\) is large enough as a distance estimator for Julia fractals. We can find both \(f^n\) and \(f’^n\) iteratively.
Since Julia fractals exist at a four-dimensional hyperspace but we can render objects in a three-dimensional space, we should project the fractals into some three-dimensional space. In this project, we simply project them into the \(k=k’\) space for some real number \(k’\).
The below de function is the distance estimator for Julia fractals.
float de(vec3 point) {
quat q = quat(point, k);
quat dq = quat(1, 0, 0, 0);
double r;
for (int i = 0; i < iter; i++) {
dq = 2 * q * dq;
q = q * q + d;
r = size(q);
}
return 0.5 * log(r) * r / size(dq);
}The code is simplified for conciseness and one can find real implementation in the Objects.cpp file.
Our implementation basically follows the paper.
In the paper, Section 3 discussed a way to find intersections between rays and fractal objects. We use the ray marching as they suggested in Section 3.1 and 3.2. In addition, we introduce a maximum distance of an increment as mentioned in Section 3.4. On the other hand, they suggested a minimum distance of an increment for faster convergence in Section 3.3 but, if an estimated distance is shorter than the minimum distance, we simply determine that a ray intersects with an object instead of forwarding a ray by the minimum distance. Our distance estimator for the Julia fractals is identical to their suggestion in Section 3.1.
Section 4 was about estimating normal vectors and the clarity of images. They proposed two methods to estimated normal vectors: Z-buffer neighbors in Section 4.1.1 and Gradient Computation in Section 4.1.2. We implement only gradient computation since it fits better to the ray tracing. Section 4.2 showed how to control the clarity of images by adjusting minimum distances dynamically. However, because we do not use a minimum distance as the criterion of intersecting rather than the shortest distance to proceed a ray, we decide not to implement dynamic minimum distances.
We can easily add shadows into images by shooting a ray from an intersection point to a light source as done in the assignment.
The below lighting function calculates the color at a given point by considering shadows.
vec3 lighting(vec3 point, vec3 normal, vec3 view, Material material) {
vec3 color = ambience * material.ambient;
for (Light light: lights) {
vec3 l = normalize(light.position - point);
vec3 r = normalize(mirror(l, normal));
double dd = dot(normal, l), ds = dot(view, r);
Ray ray(point + margin * l, l);
float t;
if (!(intersect(ray, t) && t < norm(point - light.position)))
color += light.color *
(material.diffuse * max(dd, 0) +
material.specular * pow(max(ds, 0), material.shininess));
}
return color;
}The code is simplified for conciseness and one can find real implementation in the Scene.cpp file.
To create more realistic images, we can use soft shadows by treating a light source as a small circle instead of a single point. We simply assume that light source is a circle, which is parallel to the tangent plane of an object at an intersection point. The intensity of light reaches the intersection point is a proportion of the circle that is not hidden by the object. To calculate the proportion hidden, we use the Monte Carlo method. We shoot lots of rays from the intersection and decide the proportion stochastically. It is computationally expensive but easy to implement. Since performance is not the most important aim for off-line rendering, we choose it. The following is the mathematical validation of our implementation:
\(~~~~{\sf soft\_shadow}(P)\)
\(=\frac{1}{\pi R^2}\int\int_{C}{\sf Intersect}(P,L(r,\theta))dA\)
\(=\frac{1}{\pi R^2}\int_{0}^{2\pi}\int_{0}^{R}{\sf Intersect}(P,L(r,\theta))rdrd\theta\)
\(=\frac{2\pi R}{\pi R^2}\int_{0}^{2\pi}\int_{0}^{R}{\sf Intersect}(P,L(r,\theta))r\frac{1}{R}\frac{1}{2\pi}drd\theta\)
\(=\frac{2}{R}E[{\sf Intersect}(P,L(r,\theta))\times r]\) where \(r\sim U[0,R]\) and \(\theta\sim U[0,2\pi]\)
\(\approx\frac{2}{RN}\sum_{i=1}^{N}{\sf Intersect}(P,L(r_i,\theta_i))\times r_i\)
(\({\sf Intersect}\) returns one if a ray intersects and zero otherwise.)
We can implement it as the following function:
double soft_shadow(vec3 light, vec3 point) {
double sum = 0;
for (int i = 0; i < N; i++) {
double r = rand_uniform(0, R);
double theta = rand_uniform(0, 2 * M_PI);
vec3 rand_light = adjust(light, point, r, theta);
vec3 l = rand_light - point;
double distance = norm(l);
Ray ray(point + margin * l, l);
float t;
if (intersect(ray, t) && t < distance) sum += r;
}
return 2 * sum / N / R;
}The lighting function can be modified as the following:
vec3 lighting(vec3 point, vec3 normal, vec3 view, Material material) {
vec3 color = ambience * material.ambient;
for (Light light: lights) {
vec3 l = normalize(light.position - point);
vec3 r = normalize(mirror(l, normal));
double dd = dot(normal, l), ds = dot(view, r);
color += light.color * (1 - soft_shadow(light.position, point)) *
(material.diffuse * max(dd, 0) +
material.specular * pow(max(ds, 0), material.shininess));
}
return color;
}The code is simplified for conciseness and one can find real implementation in the Scene.cpp file.
Ambient occlusion calculates how much a point receives the ambient light. In the Phong lighting, the intensity of the ambient light is constant so that every point is exposed to an equal amount of the ambient light. However, this assumption is not realistic since if more surface surrounds a point, then less ambient light is able to reach the point. Therefore, ambient occlusion is an important technique to render realistic images and to express detail of images.
Like soft shadows in the previous section, we use the Monte Carlo method to implement ambient occlusion. We shoot rays from a point to random points above the tangent plane of an object at the point. As more rays intersect with objects, the point is more occluded and receives less ambient light. The following is the mathematical validation of our implementation:
\(~~~~{\sf occlusion}(P)\)
\(=\frac{1}{2\pi}\int\int_{H}{\sf Intersect}(\phi,\theta)dS\)
\(=\frac{1}{2\pi}\int_{0}^{\pi/2}\int_{0}^{2\pi}{\sf Intersect}(\phi,\theta){\sf sin}\theta d\phi d\theta\)
\(=\frac{2\pi\times\pi/2}{2\pi}\int_{0}^{\pi/2}\int_{0}^{2\pi}{\sf Intersect}(\phi,\theta){\sf sin}\theta\frac{1}{2\pi}\frac{1}{\pi/2}d\phi d\theta\)
\(=\frac{\pi}{2}E[{\sf Intersect}(\phi,\theta)\times {\sf sin}\theta]\) where \(\phi\sim U[0,2\pi]\) and \(\theta\sim U[0,\pi/2]\)
\(\approx\frac{\pi}{2N}\sum_{i=1}^{N}{\sf Intersect}(\phi_i,\theta_i)\times {\sf sin}\theta_i\)
(\({\sf Intersect}\) returns one if a ray intersects and zero otherwise.)
We can implement it as the following function:
double occlusion(vec3 point, vec3 normal) {
double sum = 0;
for (int i = 0; i < N; i++) {
double phi = rand_uniform(0, M_PI * 2.0);
double theta = rand_uniform(0, M_PI * 0.5);
vec3 dir = calculate_dir(normal, phi, theta);
Ray ray(point + margin * dir, dir);
if (intersect(ray)) sum += sin(theta);
}
return sum * M_PI / 2 / N;
}The lighting function can be modified as the following:
vec3 lighting(vec3 point, vec3 normal, vec3 view, Material material) {
vec3 color = ambience * (1 - occlusion(point, normal)) * material.ambient;
for (Light light: lights) {
vec3 l = normalize(light.position - point);
vec3 r = normalize(mirror(l, normal));
double dd = dot(normal, l), ds = dot(view, r);
color += light.color *
(material.diffuse * max(dd, 0) +
material.specular * pow(max(ds, 0), material.shininess));
}
return color;
}The code is simplified for conciseness and one can find real implementation in the Scene.cpp file.
To implement transparent materials, we refer to an article in www.scratchapixel.com. We consider only completely transparent materials. Absorption at a surface or inside material does not exist in our implementation. When light meets a transparent object, some portion of the light is refracted and the remaining is reflected. We need to calculate the direction of the refracted ray in addition to the reflected ray. Snell’s law determines the direction:
\(\frac{{\sf sin}\theta_1}{{\sf sin}\theta_2}=\frac{\eta_2}{\eta_1}\)
\(\theta_1\) and \(\theta_2\) are the angles of incidence and refraction, respectively. \(\eta_1\) and \(\eta_2\) are the indices of refraction of the air and the material of an object, respectively.
By using the above facts, we can write the direction of the refracted ray \(T\) as the following: \(T=\frac{\eta_1}{\eta_2}(I+{\sf cos}\theta_1 N)-N\sqrt{1-(\frac{\eta_1}{\eta_2})^2(1-{\sf cos}^2\theta_1)}\) where \(N\) is the normal vector. \(\eta_1\) is always one in the code. \(\eta_2\) is stored in variable ior in the code. One thing we have to notice is a critical angle. If an angle of incidence is larger than the critical angle, then light is fully reflected, not refracted. It happens when \(\sqrt{1-(\frac{\eta_1}{\eta_2})^2(1-{\sf cos}^2\theta_1)}\) is negative.
We can implement it as the following function:
vec3 refract(vec3 v, vec3 n, float ior) {
float cosi = max(min(dot(v, n), 1), -1);
float etai = 1, etat = ior;
vec3 normal_ref = n;
if (cosi < 0) cosi = -cosi;
else { etat = 1; etai = ior; normal_ref = -n; }
float eta = etai / etat;
float k = 1 - eta * eta * (1 - cosi * cosi);
return (k < 0) ? vec3(0) : normalize(eta * v + (eta * cosi - sqrt(k)) * n);
}The function returns the direction of the refracted ray. The code is simplified for conciseness and one can find real implementation in the vec3.h file.
The Fresnel equations calculate the amount of reflected and refracted light:
\(F_{R\parallel}=(\frac{\eta_2{\sf cos}\theta_1-\eta_1{\sf cos}\theta_2}{\eta_2{\sf cos}\theta_1+\eta_1{\sf cos}\theta_2})^2\)
\(F_{R\perp}=(\frac{\eta_1{\sf cos}\theta_1-\eta_2{\sf cos}\theta_2}{\eta_1{\sf cos}\theta_1+\eta_2{\sf cos}\theta_2})^2\)
By taking the average of the two values, we get the actual proportion of reflected light:
\(F_R=\frac{1}{2}(F_{R\parallel}+F_{R\perp})\)
If \({\sf sin}\theta_2\) (sint in the code) is larger than one, the function returns one because it means that the incidence angle is larger than the critical angle.
We can implement it as the following function:
float fresnel(vec3 v, vec3 n, float ior) {
float cosi = max(min(dot(v, n), 1), -1);
float etai = 1, etat = ior;
if (cosi > 0) { etat = 1; etai = ior; }
float sint = etai / etat * sqrt(max(1 - cosi * cosi, 0));
if (sint >= 1) return 1;
else {
float cost = sqrt(max(1 - sint * sint, 0));
cosi = abs(cosi);
float Rs = ((etat * cosi) - (etai * cost)) / ((etat * cosi) + (etai * cost));
float Rp = ((etai * cosi) - (etat * cost)) / ((etai * cosi) + (etat * cost));
return (Rs * Rs + Rp * Rp) / 2;
}
}The code is simplified for conciseness and one can find real implementation in the vec3.h file.
The remaining is the proportion refracted.
\(F_T=1-F_R\)
Finally, we calculate the colors of refracted and reflected rays using recurisive calls. The colors are mixed with weight \(F_R\) (kr in the code) from the fresnel function.
We add the following code into the trace function.
vec3 trace(Ray ray, int depth) {
...
auto recursive_trace = [&](vec3 direction) {
vec3 origin = point + margin * direction;
Ray reflected_ray(origin, direction);
return trace(reflected_ray, depth + 1);
};
if (material.transparent) {
vec3 color_reflection = recursive_trace(reflect(ray.direction, normal));
vec3 color_refraction = recursive_trace(refract(ray.direction, normal, material.ior));
float kr = fresnel(ray.direction, normal, material.ior);
return mix(color_refraction, color_reflection, kr);
} else {
...
}The code is simplified for conciseness and one can find real implementation in the Scene.cpp file.
We refer to a blog article written by Mikael Hvidtfeldt Christensen to render the Mandelbulb Fractals. The Mandelbulb fractals are similar to the Julia fractals. A fractal contains every \(z\) such that \(f^n(z)\) does not diverges as \(n\) goes to the infinity where \(f(x)=x^2+z\) and superscript \(n\) denotes \(n\)-fold application of the function. Christensen noted that we can use higher powers, \(f(x)=x^8+z\) for example, to create more intersting images. According to Christensen, we can interpret the product as the rotation and scaling of a three-dimensional vector. Precisely, squaring a vector corresponds to squaring the length of the vector and doubling the angle of the vector.
We can implement it as the following function:
float de(vec3 point) {
vec3 p = point;
float dr = 1, r;
for (int i = 0; i < iter; i++) {
r = norm(p);
float theta = acos(p[2] / r);
float phi = atan2(p[1], p[0]);
dr = power * pow(r, power - 1) * dr + 1;
theta *= power;
phi *= power;
float st = sin(theta), ct = cos(theta), sp = sin(phi), cp = cos(phi);
p = pow(r, power) * vec3(st * cp, st * sp, ct) + point;
}
return 0.5 * log(r) * r / dr;
}The code is simplified for conciseness and one can find real implementation in the Objects.cpp file.
We can render Sierpinski tetrahedrons by defining a proper distance estimator. We refer to a blog article written by Christensen. However, the article did not discuss how to find the closest distance from a point to a tetrahedron. Therefore, we implement a function calculating the distance by ourselves. The code is complex and long so that we omit it here. One can find real implementation in the Tetra.cpp file.
Bezier curves are useful to create natural movements of a camera. Since the concept of Bezier curves is discussed in detail in the lecture, we do not explain it here. Now, the position of a camera is not fixed but changes as time passes.
We add the following code into the compute_vecs function.
void compute_vecs() {
vector<vec3> points = control_polygon;
double t = time / duration;
int length = points.size();
while (length-- > 1)
for (int i = 0; i < length; i++)
points[i] = mix(points[i], points[i + 1], t);
eye = points[0];
...
}The code is simplified for conciseness and one can find real implementation in the Camera.cpp file.
Every image in this section was rendered using the following configuration:
| Property | Value |
|---|---|
| Eye | (-2, 0, 4) |
| Center | (0, 0, 0) |
| Up | (0, 1, 0) |
| Field of View (Y) | 45 deg |
| Width | 640 px |
| Height | 360 px |
| Property | Value |
|---|---|
| Background Color | (0, 0, 0) |
| Ambience Light Color | (0.2, 0.2, 0.2) |
| Property | Value |
|---|---|
| Position | (-3, 1, 10) |
| Color | (1, 1, 1) |
| Property | Value |
|---|---|
| Ambient Color | (0.6, 0.6, 0.6) |
| Diffuse Color | (0.6, 0.6, 0.6) |
| Specular Color | (1, 1, 1) |
| Shininess | 50 |
| Mirror | 0 |
Rendering times were measured in a Google Cloud Platform Compute Engine n1-highcpu-8 instance.
| Property | Value |
|---|---|
| d | (-1, 0.2, 0, 0) |
| k | 0 |
| No shadows | Simple shadows | Soft shadows | |
|---|---|---|---|
| Constant ambience | 168.054 | 193.242 | 3266.59 |
| Ambient occlusion | 3099.13 | 3147.36 | 6135.28 |
| Property | Value |
|---|---|
| power | 8 |
| No shadows | Simple shadows | Soft shadows | |
|---|---|---|---|
| Constant ambience | 691.276 | 938.761 | 23244.7 |
| Ambient occlusion | 23462.8 | 23666 | 46639.2 |
| Property | Value |
|---|---|
| size | 1 |
| iter | 4 |
| No shadows | Simple shadows | Soft shadows | |
|---|---|---|---|
| Constant ambience | 2800.13 | 3477.99 | 82592.9 |
| Ambient occlusion | 92200.2 | 92943.7 | 171737 |
| Property | Value |
|---|---|
| Control Polygon | (-1, 0, 4.5) |
| (-1.5, 0, 3) | |
| (8, -0.5, 4) | |
| (3.5, -1, -7) | |
| (-4, -1.5, -2) | |
| (-0.5, -2, 1) | |
| (-0.25, -4, 0.5) | |
| Width | 1280 px |
| Height | 720 px |
| FPS | 30 |
| Frames | 201 |
| Property | Value |
|---|---|
| Position 1 | (-5, 0, 10) |
| Color 1 | (0.8, 0.8, 0.8) |
| Position 2 | (5, 0, -10) |
| Color 2 | (0.8, 0.8, 0.8) |
| Property | Value |
|---|---|
| Ambient Color | (0.6, 0.6, 0.6) |
| Diffuse Color | (0.6, 0.6, 0.6) |
| Specular Color | (1, 1, 1) |
| Shininess | 50 |
| Property | Value |
|---|---|
| d | (-1, 0.2, 0, 0) |
| k | 0 |
| Property | Value |
|---|---|
| Control Polygon | (-1, 0, 4.5) |
| (-0.5, -1, 2) | |
| (0, 4, 0) | |
| (-2, 2, -0.5) | |
| (-1, 2, -1) | |
| (2, 0, -2) | |
| Width | 1280 px |
| Height | 720 px |
| FPS | 30 |
| Frames | 201 |
| Property | Value |
|---|---|
| Background Color | (0, 0, 0) |
| Ambience Light Color | (0.2, 0.2, 0.2) |
| Property | Value |
|---|---|
| Position 1 | (-5, 0, 10) |
| Color 1 | (0.8, 0.8, 0.8) |
| Position 2 | (5, 0, -10) |
| Color 2 | (0.8, 0.8, 0.8) |
| Property | Value |
|---|---|
| Ambient Color | (0.6, 0.6, 0.6) |
| Diffuse Color | (0.6, 0.6, 0.6) |
| Specular Color | (1, 1, 1) |
| Shininess | 50 |
| Property | Value |
|---|---|
| d | (-0.213,-0.0410,-0.563,-0.560) |
| k | 0 |
| Property | Value |
|---|---|
| Control Polygon | (-1, 0, 4.5) |
| (-1.5, 0, 3) | |
| (8, -0.5, 4) | |
| (3.5, -1, -7) | |
| (-4, -1.5, -2) | |
| (-0.5, -2, 1) | |
| (-0.25, -4, 0.5) | |
| Width | 1280 px |
| Height | 720 px |
| FPS | 30 |
| Frames | 201 |
| Property | Value |
|---|---|
| Position | (3, 3, 3) |
| Color | (1, 1, 1) |
| Property | Value |
|---|---|
| Transparent | true |
| Index of Refraction | 1.5 |
| Property | Value |
|---|---|
| d | (-1, 0.2, 0, 0) |
| k | 0 |
| Property | Value |
|---|---|
| Control Polygon | (-1, 0, 4.5) |
| (-1.5, 0, 3) | |
| (0, 4, 0) | |
| (6, 0, 3) | |
| (6, -1, -3) | |
| (-3, -2, -1.5) | |
| (-3, -3, -1.5) | |
| Width | 640 px |
| Height | 360 px |
| FPS | 30 |
| Frames | 201 |
| Property | Value |
|---|---|
| Position 1 | (-5, 0, 10) |
| Color 1 | (0.8, 0.8, 0.8) |
| Position 2 | (5, 0, -10) |
| Color 2 | (0.8, 0.8, 0.8) |
| Property | Value |
|---|---|
| Ambient Color | (0.6, 0.6, 0.6) |
| Diffuse Color | (0.6, 0.6, 0.6) |
| Specular Color | (1, 1, 1) |
| Shininess | 50 |
| Property | Value |
|---|---|
| power | 8 |
| Property | Value |
|---|---|
| Control Polygon | (-4, -4, 0) |
| (-2.5, -2.5, 0) | |
| (2.5, -2.5, 0) | |
| (2.2, 0, 2.2) | |
| (2.2, 0, -2.2) | |
| (2.2, 2.2, -2.2) | |
| (0.1, 0.1, -0.1) | |
| Width | 640 px |
| Height | 360 px |
| FPS | 30 |
| Frames | 201 |
| Property | Value |
|---|---|
| Position 1 | (-5, 0, 10) |
| Color 1 | (0.8, 0.8, 0.8) |
| Position 2 | (5, 0, -10) |
| Color 2 | (0.8, 0.8, 0.8) |
| Property | Value |
|---|---|
| Ambient Color | (0.8, 0.8, 0.8) |
| Diffuse Color | (0.8, 0.8, 0.8) |
| Specular Color | (1, 1, 1) |
| Shininess | 100 |
| Property | Value |
|---|---|
| size | 1 |
| iter | 4 |