Literate Ray Tracer in Clojure
by Abhirag
Table of Contents
Overview
A literate ray tracer implemented in Clojure, based on the book Ray Tracing in One Weekend by Peter Shirley. The orginal source code for the book can be found here.
Output an image
Whenever you start a renderer, you need a way to see an image. The most straightforward way is to write it to a file. The catch is, there are so many formats and many of those are complex. I always start with a plain text ppm file. Here’s a nice description from Wikipedia:
This is an example of a color RGB image stored in PPM format. There is a newline character at the end of each line.
P3 3 2 255 # The part above is the header # "P3" means this is a RGB color image in ASCII # "3 2" is the width and height of the image in pixels # "255" is the maximum value for each color # The part below is image data: RGB triplets 255 0 0 0 255 0 0 0 255 255 255 0 255 255 255 0 0 0
Let’s write some Clojure code to output such a thing:
(let [nx 200 ny 100 header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [r (/ i nx) g (/ j ny) b 0.2 ir (int (* 255.99 r)) ig (int (* 255.99 g)) ib (int (* 255.99 b))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_1.ppm" header :append true) (spit "./output_images/image_1.ppm" body :append true))
There are some things to note in the code above:
- The pixels are written out in rows with pixels left to right.
- The rows are written out from top to bottom.
- By convention, each of the red/green/blue components range from 0.0 to 1.0.
We will relax that later when we internally use high dynamic range, but before
output we will tone map to the zero to one range, so this code won’t change. - Red goes from black to fully on from left to right, and green goes from black at the
bottom to fully on at the top. Red and green together make yellow so we should expect
the upper right corner to be yellow.
Opening the output file shows: try it in your favorite image viewer and google “ppm viewer” if your viewer doesn’t support it
Hooray! This is the graphics “hello world”. If your image doesn’t look like that, open the output file in a text editor and see what it looks like. It should start something like this:
P3 200 100 255 0 253 51 1 253 51 2 253 51 3 253 51 5 253 51
Vectors
Almost all graphics programs have some way of representing geometric vectors and colors. In many systems these vectors are 4D (3D plus a homogeneous coordinate for geometry, and RGB plus an alpha transparency channel for colors). For our purposes, three coordinates suffices. We’ll use the same representation for colors, locations, directions, offsets, whatever. Some people don’t like this because it doesn’t prevent you from doing something silly, like adding a color to a location. They have a good point, but we’re going to always take the “less code” route when not obviously wrong.
Here’s my code for dealing with vectors:
Defining a vector:
(def vec3 [230 411 518])
Accessing a vector:
(println (get vec3 0)) (println (get vec3 1)) (println (get vec3 2))
230 411 518
Vector Operations:
(defn vec3-squared-length [[a1 a2 a3]] (+ (* a1 a1) (* a2 a2) (* a3 a3))) (defn vec3-length [vec3] (Math/sqrt (vec3-squared-length vec3))) (defn vec3-product [[a1 a2 a3] [b1 b2 b3]] [(* a1 b1) (* a2 b2) (* a3 b3)]) (defn vec3-plus [[a1 a2 a3] [b1 b2 b3]] [(+ a1 b1) (+ a2 b2) (+ a3 b3)]) (defn vec3-minus [[a1 a2 a3] [b1 b2 b3]] [(- a1 b1) (- a2 b2) (- a3 b3)]) (defn vec3-division [[a1 a2 a3] [b1 b2 b3]] [(/ a1 b1) (/ a2 b2) (/ a3 b3)]) (defn vec3-scalar-multiply [scalar vec3] (mapv #(* scalar %) vec3))
The vec3-division and vec3-product operations are for colors and it’s unlikely you want to use them for things like locations. Similarly, we’ll need some geometric operations occasionally:
(defn vec3-dot-product [[a1 a2 a3] [b1 b2 b3]] (+ (* a1 b1) (* a2 b2) (* a3 b3))) (defn vec3-cross-product [[a1 a2 a3] [b1 b2 b3]] [(- (* a2 b3) (* a3 b2)), (- (* a3 b1) (* a1 b3)) (- (* a1 b2) (* a2 b1))])
And to make a unit vector in the same direction as the input vector:
(defn vec3-unit-vector [vec3] (vec (for [elem vec3] (/ elem (vec3-length vec3)))))
(println (vec3-squared-length [3 4 5])) (println (vec3-length [3 4 5])) (println (vec3-product [1 2 3] [4 5 6])) (println (vec3-plus [1 2 3] [4 5 6])) (println (vec3-minus [1 2 3] [4 5 6])) (println (vec3-division [1 2 3] [4 5 6])) (println (vec3-dot-product [1 3 -5] [4 -2 -1])) (println (vec3-cross-product [1 2 3] [4 5 6])) (println (vec3-unit-vector [0 3 4])) (println (vec3-scalar-multiply 2 [1 2 3]))
50 7.0710678118654755 [4 10 18] [5 7 9] [-3 -3 -3] [1/4 2/5 1/2] 3 [-3 6 -3] [0.0 0.6 0.8] [2 4 6]
Now we can change the code we wrote to output a PPM image to use this:
(let [nx 200 ny 100 header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [color [(/ i nx) (/ j ny) 0.2] ir (int (* 255.99 (get color 0))) ig (int (* 255.99 (get color 1))) ib (int (* 255.99 (get color 2)))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_2.ppm" header :append true) (spit "./output_images/image_2.ppm" body :append true))
Rays, a simple camera, and background
The one thing that all ray tracers have is a ray representation, and a computation of what color is seen along a ray. Let’s think of a ray as a function \(p(t) = A + t \times B\). Here \(p\) is a 3D position along a line in 3D. \(A\) is the ray origin and \(B\) is the ray direction. The ray parameter \(t\) is a real number. Plug in a different \(t\) and \(p(t)\) moves the point along the ray. Add in negative \(t\) and you can go anywhere on the 3D line. For positive \(t\), you get only the parts in front of \(A\), and this is what is often called a half-line or ray. The example \(C = p(2)\) is shown here:
The function \(p(t)\) in more verbose code form I call \(point\_at\_parameter(t)\):
(defn create-ray [vec3-origin vec3-direction] {:origin vec3-origin :direction vec3-direction}) (defn ray-origin [ray] (:origin ray)) (defn ray-direction [ray] (:direction ray)) (defn point_at_parameter [{:keys [origin direction]} t] (vec3-plus origin (vec3-scalar-multiply t direction)))
Now we are ready to turn the corner and make a ray tracer. At the core of a ray tracer is to send rays through pixels and compute what color is seen in the direction of those rays. This is of the form calculate which ray goes from the eye to a pixel, compute what that ray intersects, and compute a color for that intersection point. When first developing a ray tracer, I always do a simple camera for getting the code up and running. I also make a simple \(\mathtt{color(ray)}\) function that returns the color of the background (a simple gradient).
I’ve often gotten into trouble using square images for debugging because I transpose x and
y too often, so I’ll stick with a 200 by 100 image. I’ll put the “eye” (or camera center if you
think of a camera) at (0,0,0). I will have the y-axis go up, and the x-axis to the right.
In
order to respect the convention of a right handed coordinate system, into the screen is the
negative z-axis. I will traverse the screen from the lower left hand corner and use two
offset vectors along the screen sides to move the ray endpoint across the screen. Note that
I do not make the ray direction a unit length vector because I think not doing that makes
for simpler and slightly faster code.
Below in code, the ray goes to approximately the pixel centers: I won’t worry about exactness for now because we’ll add antialiasing later
(defn color [ray] (let [direction (ray-direction ray) y (get direction 1) t (* 0.5 (+ 1 y)) white [1 1 1] blue [0.5 0.7 1]] (vec3-plus (vec3-scalar-multiply (- 1 t) white) (vec3-scalar-multiply t blue)))) (let [nx 200 ny 100 lower-left-corner [-2 -1 -1] origin [0 0 0] horizontal [4 0 0] vertical [0 2 0] header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [u (/ i nx) v (/ j ny) vector-on-plane (vec3-plus (vec3-scalar-multiply u horizontal) (vec3-scalar-multiply v vertical)) direction (vec3-plus lower-left-corner vector-on-plane) ray (create-ray origin direction) pixel-color (color ray) ir (int (* 255.99 (get pixel-color 0))) ig (int (* 255.99 (get pixel-color 1))) ib (int (* 255.99 (get pixel-color 2)))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_3.ppm" header :append true) (spit "./output_images/image_3.ppm" body :append true))
The \(\mathtt{color(ray)}\) function linearly blends white and blue depending on the up/downess of the \(y\) coordinate. \(-1.0 \lt y \lt 1.0\), so I did a standard graphics trick of scaling that to \(0.0 \lt t \lt 1.0\). When \(t = 1.0\), I want blue. When \(t = 0.0\), I want white. In between, I want a blend. This forms a “linear blend”, or “linear interpolation”, or “lerp” for short, between two things. A lerp is always of the form: \(blended\_value = (1- t) \times start\_value + t \times end\_value\), with \(t\) going from zero to one. In our case this produces:
Adding a sphere
Let’s add a single object to our ray tracer. People often use spheres in ray tracers because calculating whether a ray hits a sphere is pretty straightforward. Recall that the equation for a sphere centered at the origin of radius \(R\) is \(x^2 + y^2 + z^2 = R^2\). The way you can read that equation is “for any \((x, y, z)\), if \(x^2 + y^2 + z^2 = R^2\) then \((x,y,z)\) is on the sphere and otherwise it is not”. It gets uglier if the sphere center is at \((cx, cy, cz)\): \[(x-cx)^2 + (y-cy)^2 + (z-cz)^2 = R^2\] In graphics, you almost always want your formulas to be in terms of vectors. You might note that the vector from center \(C = (cx,cy,cz)\) to point \(p = (x,y,z)\) is \((p - C)\). And \(dot((p - C), (p - C)) = (x-cx)^2 + (y-cy)^2 + (z-cz)^2\). So the equation of the sphere in vector form is: \[dot((p - C), (p - C)) = R^2\] We can read this as “any point \(p\) that satisfies this equation is on the sphere”. We want to know if our ray \(p(t) = A + t \times B\) ever hits the sphere anywhere. If it does hit the sphere, there is some \(t\) for which \(p(t)\) satisfies the sphere equation. So we are looking for any \(t\) where this is true: \[dot((p(t) - C), (p(t) - C)) = R^2\] or expanding the full form of the ray \(p(t)\): \[dot((A + t \times B - C),(A + t \times B - C)) = R^2\] The rules of vector algebra are all that we would want here, and if we expand that equation and move all the terms to the left hand side we get: \[t^2 \times dot(B, B) + 2 \times t \times dot(A-C, B) + dot(A-C, A-C) - R^2 = 0\] The vectors and \(R\) in that equation are all constant and known. The unknown is \(t\), and the equation is a quadratic, like you probably saw in your high school math class. You can solve for \(t\) and there is a square root part that is either positive (meaning two real solutions), negative (meaning no real solutions), or zero (meaning one real solution). In graphics, the algebra almost always relates very directly to the geometry. What we have is:
If we take that math and hard-code it into our program, we can test it by coloring red any pixel that hits a small sphere we place at -1 on the z-axis:
(defn hit-sphere [center radius ray] (let [origin (ray-origin ray) direction (ray-direction ray) oc (vec3-minus origin center) a (vec3-dot-product direction direction) b (* 2 (vec3-dot-product oc direction)) c (- (vec3-dot-product oc oc) (* radius radius)) discriminant (- (* b b) (* 4 a c))] (> discriminant 0))) (defn color [ray] (if (hit-sphere [0 0 -1] 0.5 ray) (let [red [1 0 0]] red) (let [direction (ray-direction ray) y (get direction 1) t (* 0.5 (+ 1 y)) white [1 1 1] blue [0.5 0.7 1]] (vec3-plus (vec3-scalar-multiply (- 1 t) white) (vec3-scalar-multiply t blue)))))
What we get is this:
Now this lacks all sorts of things — like shading and reflection rays and more than one object — but we are closer to halfway done than we are to our start! One thing to be aware of is that we tested whether the ray hits the sphere at all, but \(t \lt 0\) solutions work fine. If you change your sphere center to \(z = 1\) you will get exactly the same picture because you see the things behind you. This is not a feature! We’ll fix those issues next.
Surface normals and multiple objects
First, let’s get ourselves a surface normal so we can shade. This is a vector that is perpendicular to the surface, and by convention, points out. One design decision is whether these normals (again by convention) are unit length. That is convenient for shading so I will say yes, but I won’t enforce that in the code. This could allow subtle bugs, so be aware this is personal preference as are most design decisions like that. For a sphere, the normal is in the direction of the hitpoint minus the center:
On the earth, this implies that the vector from the earth’s center to you points straight up. Let’s throw that into the code now, and shade it. We don’t have any lights or anything yet, so let’s just visualize the normals with a color map. A common trick used for visualizing normals (because it’s easy and somewhat intuitive to assume \(N\) is a unit length vector– so each component is between -1 and 1) is to map each component to the interval from 0 to 1, and then map x/y/z to r/g/b. For the normal we need the hit point, not just whether we hit or not. Let’s assume the closest hit point (smallest \(t\)). These changes in the code let us compute and visualize \(N\):
(defn hit-sphere [center radius ray] (let [origin (ray-origin ray) direction (ray-direction ray) oc (vec3-minus origin center) a (vec3-dot-product direction direction) b (* 2 (vec3-dot-product oc direction)) c (- (vec3-dot-product oc oc) (* radius radius)) discriminant (- (* b b) (* 4 a c))] (if (< discriminant 0) -1 (/ (- (- b) (Math/sqrt discriminant)) (* 2 a))))) (defn color [ray] (let [t (hit-sphere [0 0 -1] 0.5 ray)] (if (> t 0) (let [N (vec3-unit-vector (vec3-minus (point_at_parameter ray t) [0 0 -1])) color-normal (vec3-scalar-multiply 0.5 (mapv + N [1 1 1]))] color-normal) (let [direction (ray-direction ray) y (get direction 1) t (* 0.5 (+ 1 y)) white [1 1 1] blue [0.5 0.7 1]] (vec3-plus (vec3-scalar-multiply (- 1 t) white) (vec3-scalar-multiply t blue))))))
And that yields this picture:
Now, how about several spheres? While it is tempting to have an array of spheres, a very clean solution is to make a representation for anything a ray might hit and make both a sphere and a list of spheres just something you can hit. What that representation should be called is something of a quandary— calling it an “object” would be good if not for “object oriented” programming. “Surface” is often used, with the weakness being maybe we will want volumes. “Hitable” emphasizes the function that unites them. I don’t love any of these but I will go with “hitable”.
Most ray tracers have found it convenient to add a valid interval for hits \(t_{min}\) to \(t_{max}\), so the hit only “counts” if \(t_{min} \lt t \lt t_{max}\). For the initial rays this is positive \(t\), but as we will see, it can help some details in the code to have an interval \(t_{min}\) to \(t_{max}\). One design question is whether to do things like compute the normal if we hit something; we might end up hitting something closer as we do our search, and we will only need the normal of the closest thing. I will go with the simple solution and compute a bundle of stuff I will store in some hash-map. I know we’ll want motion blur at some point, so I’ll add a time input variable. Here’s the code template for such a thing:
(let [hit-record {:t 0.0 :p [0 0 0] :normal [0 0 0]}] (defn hit [ray t-min t-max hit-record]))
And here’s the sphere (note that I eliminated a bunch of redundant 2’s that cancel each other out):
(defn hit-sphere [ray sphere t-min t-max hit-record] (let [origin (ray-origin ray) direction (ray-direction ray) center (:center sphere) radius (:radius sphere) oc (vec3-minus origin center) a (vec3-dot-product direction direction) b (vec3-dot-product oc direction) c (- (vec3-dot-product oc oc) (* radius radius)) discriminant (- (* b b) (* a c)) root1 (/ (- (- b) (Math/sqrt discriminant)) a) root2 (/ (+ (- b) (Math/sqrt discriminant)) a)] (if (< discriminant 0) [false hit-record] (cond (and (< root1 t-max) (> root1 t-min)) [true (assoc hit-record :t root1 :p (point_at_parameter ray root1) :normal (vec3-unit-vector (vec3-minus (point_at_parameter ray root1) center)))] (and (< root2 t-max) (> root2 t-min)) [true (assoc hit-record :t root2 :p (point_at_parameter ray root2) :normal (vec3-unit-vector (vec3-minus (point_at_parameter ray root2) center)))] :else [false hit-record]))))
And a list of objects:
(defn hit-list [ray hitable-list t-min t-max hit-record] (loop [l hitable-list closest-so-far t-max hit-anything? false h-rec hit-record] (if (empty? l) [hit-anything? h-rec] (let [[head & tail] l [hit? {:keys [t p normal] :as rec}] (hit-sphere ray head t-min closest-so-far h-rec)] (if hit? (recur tail t true rec) (recur tail closest-so-far hit-anything? h-rec))))))
Let's create the output image:
(defn color [ray hitable-list] (let [[hit? {:keys [t p normal]}] (hit-list ray hitable-list 0 Float/POSITIVE_INFINITY {:t 0.0 :p [0 0 0] :normal [0 0 0]})] (if hit? (let [color-normal (vec3-scalar-multiply 0.5 (mapv + normal [1 1 1]))] color-normal) (let [direction (ray-direction ray) y (get direction 1) t (* 0.5 (+ 1 y)) white [1 1 1] blue [0.5 0.7 1]] (vec3-plus (vec3-scalar-multiply (- 1 t) white) (vec3-scalar-multiply t blue)))))) (let [nx 200 ny 100 lower-left-corner [-2 -1 -1] origin [0 0 0] horizontal [4 0 0] vertical [0 2 0] hitable-list [{:center [0 -100.5 -1] :radius 100} {:center [0 0 -1] :radius 0.5}] header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [u (/ i nx) v (/ j ny) vector-on-plane (vec3-plus (vec3-scalar-multiply u horizontal) (vec3-scalar-multiply v vertical)) direction (vec3-plus lower-left-corner vector-on-plane) ray (create-ray origin direction) pixel-color (color ray hitable-list) ir (int (* 255.99 (get pixel-color 0))) ig (int (* 255.99 (get pixel-color 1))) ib (int (* 255.99 (get pixel-color 2)))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_7.ppm" header :append true) (spit "./output_images/image_7.ppm" body :append true))
This yields a picture that is really just a visualization of where the spheres are along with their surface normal. This is often a great way to look at your model for flaws and characteristics.
Antialiasing
When a real camera takes a picture, there are usually no jaggies along edges because the edge pixels are a blend of some foreground and some background. We can get the same effect by averaging a bunch of samples inside each pixel. We will not bother with stratification, which is controversial but is usual for my programs. For some ray tracers it is critical, but the kind of general one we are writing doesn’t benefit very much from it and it makes the code uglier We abstract the camera representation a bit so we can make a cooler camera later.
One thing we need is a random number generator that returns real random numbers. Whatever your infrastructure, find a function that returns a canonical random number which by convention returns random real in the range \(0 \le ran \lt 1\). The “less than” before the 1 is important as we will sometimes take advantage of that.
For a given pixel we have several samples within that pixel and send rays through each of the samples. The colors of these rays are then averaged:
Putting that all together yields this:
(defn get-ray [camera u v] (let [origin (:origin camera) horizontal (:horizontal camera) vertical (:vertical camera) lower-left-corner (:lower-left-corner camera) vector-on-plane (vec3-plus (vec3-scalar-multiply u horizontal) (vec3-scalar-multiply v vertical)) direction (vec3-plus lower-left-corner vector-on-plane)] (create-ray origin direction))) (defn average-color [nx ny ns i j camera hitable-list] (loop [n 0 color-acc [0 0 0]] (if (= n ns) (mapv / color-acc [ns ns ns]) (let [u (/ (+ i (rand)) nx) v (/ (+ j (rand)) ny) ray (get-ray camera u v) col (color ray hitable-list)] (recur (inc n) (mapv + col color-acc))))))
Let's create the output image:
(let [nx 200 ny 100 ns 100 camera {:lower-left-corner [-2 -1 -1] :horizontal [4 0 0] :vertical [0 2 0] :origin [0 0 0]} hitable-list [{:center [0 -100.5 -1] :radius 100} {:center [0 0 -1] :radius 0.5}] header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [pixel-color (average-color nx ny ns i j camera hitable-list) ir (int (* 255.99 (get pixel-color 0))) ig (int (* 255.99 (get pixel-color 1))) ib (int (* 255.99 (get pixel-color 2)))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_8.ppm" header :append true) (spit "./output_images/image_8.ppm" body :append true))
Zooming into the image that is produced, the big change is in edge pixels that are part background and part foreground:
Diffuse Materials
Now that we have objects and multiple rays per pixel, we can make some realistic looking materials. We’ll start with diffuse (matte) materials. One question is whether we can mix and match shapes and materials (so we assign a sphere a material) or if it’s put together so the geometry and material are tightly bound (that could be useful for procedural objects where the geometry and material are linked). We’ll go with separate— which is usual in most renderers— but do be aware of the limitation.
Diffuse objects that don’t emit light merely take on the color of their surroundings, but they modulate that with their own intrinsic color. Light that reflects off a diffuse surface has its direction randomized. So, if we send three rays into a crack between two diffuse surfaces they will each have different random behavior:
They also might be absorbed rather than reflected. The darker the surface, the more likely absorption is. (That’s why it is dark!) Really any algorithm that randomizes direction will produce surfaces that look matte. One of the simplest ways to do this turns out to be exactly correct for ideal diffuse surfaces (mathematically ideal Lambertian).
Pick a random point \(s\) from the unit radius sphere that is tangent to the hitpoint, and send a ray from the hitpoint \(p\) to the random point \(s\). That sphere has center \((p+N)\):
We also need a way to pick a random point in a unit radius sphere centered at the origin. We’ll use what is usually the easiest algorithm: a rejection method. First, we pick a random point in the unit cube where \(x\), \(y\), and \(z\) all range from -1 to +1. We reject this point and try again if the point is outside the sphere:
(defn random-in-unit-sphere [] (loop [p (vec3-minus (vec3-scalar-multiply 2 [(rand) (rand) (rand)]) [1 1 1])] (if (< (vec3-dot-product p p) 1) p (recur (vec3-minus (vec3-scalar-multiply 2 [(rand) (rand) (rand)]) [1 1 1]))))) (defn color [ray hitable-list] (let [[hit? {:keys [t p normal]}] (hit-list ray hitable-list 0 Float/POSITIVE_INFINITY {:t 0.0 :p [0 0 0] :normal [0 0 0]})] (if hit? (let [center-transform (vec3-plus p normal) target (vec3-plus center-transform (random-in-unit-sphere)) reflection-direction (vec3-minus target p) reflected-ray (create-ray p reflection-direction)] (vec3-scalar-multiply 0.5 (color reflected-ray hitable-list))) (let [direction (ray-direction ray) y (get direction 1) t (* 0.5 (+ 1 y)) white [1 1 1] blue [0.5 0.7 1]] (vec3-plus (vec3-scalar-multiply (- 1 t) white) (vec3-scalar-multiply t blue))))))
This gives us:
Note the shadowing under the sphere. This picture is very dark, but our spheres only absorb half the energy on each bounce, so they are 50% reflectors. If you can’t see the shadow, don’t worry, we will fix that now. These spheres should look pretty light (in real life, a light grey). The reason for this is that almost all image viewers assume that the image is “gamma corrected”, meaning the 0 to 1 values have some transform before being stored as a byte. There are many good reasons for that, but for our purposes we just need to be aware of it. To a first approximation, we can use “gamma 2” which means raising the color to the power \(1 / gamma\), or in our simple case \(1 ⁄ 2\), which is just square-root:
(let [nx 200 ny 100 ns 100 camera {:lower-left-corner [-2 -1 -1] :horizontal [4 0 0] :vertical [0 2 0] :origin [0 0 0]} hitable-list [{:center [0 -100.5 -1] :radius 100} {:center [0 0 -1] :radius 0.5}] header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [pixel-color (average-color nx ny ns i j camera hitable-list) gamma-corrected (mapv #(Math/sqrt %) pixel-color) ir (int (* 255.99 (get gamma-corrected 0))) ig (int (* 255.99 (get gamma-corrected 1))) ib (int (* 255.99 (get gamma-corrected 2)))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_10.ppm" header :append true) (spit "./output_images/image_10.ppm" body :append true))
That yields light grey, as we desire:
Metal
If we want different objects to have different materials, we have a design decision. We could have:
- A universal material with lots of parameters
- Different material types, just zero out some of those parameters
I am a fan of the first approach. For our program the material needs to do two things:
- Produce a scattered ray (or say it absorbed the incident ray)
- If scattered, say how much the ray should be attenuated
Material representation:
(def material {:type :diffuse :attenuation [0.8 0.3 0.3]})
For the Lambertian (diffuse) case we already have, it can either scatter always and attenuate by its reflectance \(R\), or it can scatter with no attenuation but absorb the fraction \(1-R\) of the rays. Or it could be a mixture of those strategies. For Lambertian materials we get this function:
(defn scatter [ray-in hit-rec] (let [{:keys [t p normal material]} hit-rec {:keys [type attenuation]} material] (cond (= type :diffuse) (let [center-transform (vec3-plus p normal) target (vec3-plus center-transform (random-in-unit-sphere)) scatter-direction (vec3-minus target p) scattered-ray (create-ray p scatter-direction)] [true scattered-ray]) :else [false nil]))) (defn hit-sphere [ray sphere t-min t-max] (let [hit-record {:t 0.0 :p [0 0 0] :normal [0 0 0] :material (:material sphere)} origin (ray-origin ray) direction (ray-direction ray) center (:center sphere) radius (:radius sphere) oc (vec3-minus origin center) a (vec3-dot-product direction direction) b (vec3-dot-product oc direction) c (- (vec3-dot-product oc oc) (* radius radius)) discriminant (- (* b b) (* a c)) root1 (/ (- (- b) (Math/sqrt discriminant)) a) root2 (/ (+ (- b) (Math/sqrt discriminant)) a)] (if (< discriminant 0) [false hit-record] (cond (and (< root1 t-max) (> root1 t-min)) [true (assoc hit-record :t root1 :p (point_at_parameter ray root1) :normal (vec3-unit-vector (vec3-minus (point_at_parameter ray root1) center)))] (and (< root2 t-max) (> root2 t-min)) [true (assoc hit-record :t root2 :p (point_at_parameter ray root2) :normal (vec3-unit-vector (vec3-minus (point_at_parameter ray root2) center)))] :else [false hit-record])))) (defn hit-list [ray hitable-list t-min t-max] (loop [l hitable-list closest-so-far t-max hit-anything? false h-rec {:t 0.0 :p [0 0 0] :normal [0 0 0] :material nil}] (if (empty? l) [hit-anything? h-rec] (let [[head & tail] l [hit? {:keys [t p normal] :as rec}] (hit-sphere ray head t-min closest-so-far)] (if hit? (recur tail t true rec) (recur tail closest-so-far hit-anything? h-rec))))))
We need to modify the color function to use this:
(defn color ([ray hitable-list depth] (let [[hit? hit-rec] (hit-list ray hitable-list 0 Float/POSITIVE_INFINITY)] (if hit? (let [[scatters? scattered-ray] (scatter ray hit-rec) {:keys [material]} hit-rec {:keys [attenuation]} material] (if (and (< depth 50) scatters?) (vec3-product attenuation (color scattered-ray hitable-list (inc depth))) [0 0 0])) (let [direction (ray-direction ray) y (get direction 1) t (* 0.5 (+ 1 y)) white [1 1 1] blue [0.5 0.7 1]] (vec3-plus (vec3-scalar-multiply (- 1 t) white) (vec3-scalar-multiply t blue)))))) ([ray hitable-list] (color ray hitable-list 0)))
For smooth metals the ray won’t be randomly scattered. The key math is: how does a ray get reflected from a metal mirror? Vector math is our friend here:
The reflected ray direction in red is just \((v + 2B)\). In our design, \(N\) is a unit vector, but \(v\) may not be. The length of \(B\) should be \(dot(v,N)\). Because \(v\) points in, we will need a minus sign yielding:
(defn reflect [v n] (vec3-minus v (vec3-scalar-multiply (* 2 (vec3-dot-product v n)) n)))
The metal material just reflects rays using this formula:
(defn scatter [ray-in hit-rec] (let [{:keys [t p normal material]} hit-rec {:keys [type attenuation]} material] (cond (= type :diffuse) (let [center-transform (vec3-plus p normal) target (vec3-plus center-transform (random-in-unit-sphere)) scatter-direction (vec3-minus target p) scattered-ray (create-ray p scatter-direction)] [true scattered-ray]) (= type :metal) (let [incident-direction (vec3-unit-vector (:direction ray-in)) reflected-direction (reflect incident-direction normal) scattered-ray (create-ray p reflected-direction) scatter? (> (vec3-dot-product reflected-direction normal) 0)] [scatter? scattered-ray]) :else [false nil])))
And add some metal spheres:
(let [nx 200 ny 100 ns 100 camera {:lower-left-corner [-2 -1 -1] :horizontal [4 0 0] :vertical [0 2 0] :origin [0 0 0]} hitable-list [{:center [0 -100.5 -1] :radius 100 :material {:type :diffuse :attenuation [0.8 0.8 0.0]}} {:center [0 0 -1] :radius 0.5 :material {:type :diffuse :attenuation [0.8 0.3 0.3]}} {:center [1 0 -1] :radius 0.5 :material {:type :metal :attenuation [0.8 0.6 0.2]}} {:center [-1 0 -1] :radius 0.5 :material {:type :metal :attenuation [0.8 0.8 0.8]}}] header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [pixel-color (average-color nx ny ns i j camera hitable-list) gamma-corrected (mapv #(Math/sqrt %) pixel-color) ir (int (* 255.99 (get gamma-corrected 0))) ig (int (* 255.99 (get gamma-corrected 1))) ib (int (* 255.99 (get gamma-corrected 2)))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_12.ppm" header :append true) (spit "./output_images/image_12.ppm" body :append true))
Which gives:
Figure 16: Metal
We can also randomize the reflected direction by using a small sphere and choosing a new endpoint for the ray:
The bigger the sphere, the fuzzier the reflections will be. This suggests adding a fuzziness parameter that is just the radius of the sphere (so zero is no perturbation). The catch is that for big spheres or grazing rays, we may scatter below the surface. We can just have the surface absorb those. We’ll put a maximum of 1 on the radius of the sphere which yields:
(= type :metal) (let [fuzz (:fuzz material) fuzz-constrained (if (< fuzz 1) fuzz 1) incident-direction (vec3-unit-vector (:direction ray-in)) reflected-direction (reflect incident-direction normal) scattered-ray (create-ray p (vec3-plus (vec3-scalar-multiply fuzz-constrained (random-in-unit-sphere)) reflected-direction)) scatter? (> (vec3-dot-product reflected-direction normal) 0)] [scatter? scattered-ray])
We can try that out by adding fuzziness 0.3 and 1.0 to the metals:
(let [nx 200 ny 100 ns 100 camera {:lower-left-corner [-2 -1 -1] :horizontal [4 0 0] :vertical [0 2 0] :origin [0 0 0]} hitable-list [{:center [0 -100.5 -1] :radius 100 :material {:type :diffuse :attenuation [0.8 0.8 0.0]}} {:center [0 0 -1] :radius 0.5 :material {:type :diffuse :attenuation [0.8 0.3 0.3]}} {:center [1 0 -1] :radius 0.5 :material {:type :metal :attenuation [0.8 0.6 0.2] :fuzz 1}} {:center [-1 0 -1] :radius 0.5 :material {:type :metal :attenuation [0.8 0.8 0.8] :fuzz 0.3}}] header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [pixel-color (average-color nx ny ns i j camera hitable-list) gamma-corrected (mapv #(Math/sqrt %) pixel-color) ir (int (* 255.99 (get gamma-corrected 0))) ig (int (* 255.99 (get gamma-corrected 1))) ib (int (* 255.99 (get gamma-corrected 2)))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_13.ppm" header :append true) (spit "./output_images/image_13.ppm" body :append true))
Dielectrics
Clear materials such as water, glass, and diamonds are dielectrics. When a light ray hits them, it splits into a reflected ray and a refracted (transmitted) ray. We’ll handle that by randomly choosing between reflection or refraction and only generating one scattered ray per interaction.
The refraction is described by Snell’s law: \[n\,sin(\theta) = n’\,sin(\theta’)\] Where \(n\) and \(n’\) are the refractive indices (typically air = 1, glass = 1.3-1.7, diamond = 2.4)and the geometry is:
One troublesome practical issue is that when the ray is in the material with the higher refractive index, there is no real solution to Snell’s law and thus there is no refraction possible. Here all the light is reflected, and because in practice that is usually inside solid objects, it is called “total internal reflection”. This is why sometimes the water-air boundary acts as a perfect mirror when you are submerged. The code for refraction is thus a bit more complicated than for reflection:
(defn refract [v n ni_over_nt] (let [uv (vec3-unit-vector v) dt (vec3-dot-product uv n) discriminant (- 1 (* ni_over_nt ni_over_nt (- 1 (* dt dt))))] (if (> discriminant 0) (let [vec-term1 (vec3-scalar-multiply ni_over_nt (vec3-minus v (vec3-scalar-multiply dt n))) vec-term2 (vec3-scalar-multiply (Math/sqrt discriminant) n) refracted (vec3-minus vec-term1 vec-term2)] [true refracted]) [false nil])))
And the dielectric material that always refracts when possible is:
(= type :dielectric) (let [incident-direction (vec3-unit-vector (:direction ray-in)) reflected-direction (reflect incident-direction normal) dot-product (vec3-dot-product incident-direction normal) ref-idx (:ref-idx material) outward-normal (if (> dot-product 0) (vec3-scalar-multiply -1 normal) normal) ni_over_nt (if (> dot-product 0) ref-idx (/ 1 ref-idx)) [refract? refracted-direction] (refract incident-direction outward-normal ni_over_nt) refracted-ray (create-ray p refracted-direction) reflected-ray (create-ray p reflected-direction)] (if refract? [true refracted-ray] [true reflected-ray]))
Attenuation is always 1, the glass surface absorbs nothing.
If we try that out with these parameters:
hitable-list [{:center [0 -100.5 -1] :radius 100 :material {:type :diffuse :attenuation [0.8 0.8 0.0]}} {:center [0 0 -1] :radius 0.5 :material {:type :diffuse :attenuation [0.8 0.3 0.3]}} {:center [1 0 -1] :radius 0.5 :material {:type :metal :attenuation [0.8 0.6 0.2] :fuzz 1}} {:center [-1 0 -1] :radius 0.5 :material {:type :dielectric :attenuation [1.0 1.0 1.0] :ref-idx 1.5}}]
We get:
Now real glass has reflectivity that varies with angle, look at a window at a steep angle and it becomes a mirror. There is a big ugly equation for that, but almost everybody uses a simple and surprisingly simple polynomial approximation by Christophe Schlick:
(defn schlick [cosine ref_idx] (let [sqrt-r0 (/ (- 1 ref_idx) (+ 1 ref_idx)) r0 (* sqrt-r0 sqrt-r0)] (+ r0 (* (- 1 r0) (Math/pow (- 1 cosine) 5)))))
This yields our full glass material:
(= type :dielectric) (let [incident-direction (vec3-unit-vector (:direction ray-in)) reflected-direction (reflect incident-direction normal) dot-product (vec3-dot-product incident-direction normal) ref-idx (:ref-idx material) outward-normal (if (> dot-product 0) (vec3-scalar-multiply -1 normal) normal) ni_over_nt (if (> dot-product 0) ref-idx (/ 1 ref-idx)) cosine (if (> dot-product 0) (/ (* ref-idx dot-product) (vec3-length incident-direction)) (- (/ dot-product (vec3-length incident-direction)))) [refract? refracted-direction] (refract incident-direction outward-normal ni_over_nt) refracted-ray (create-ray p refracted-direction) reflected-ray (create-ray p reflected-direction) reflect-prob (schlick cosine ref-idx)] (if refract? (if (< (rand) reflect-prob) [true reflected-ray] [true refracted-ray]) [true reflected-ray]))
An interesting and easy trick with dielectric spheres is to note that if you use a negative radius, the geometry is unaffected but the surface normal points inward, so it can be used as a bubble to make a hollow glass sphere:
hitable-list [{:center [0 -100.5 -1] :radius 100 :material {:type :diffuse :attenuation [0.8 0.8 0.0]}} {:center [0 0 -1] :radius 0.5 :material {:type :diffuse :attenuation [0.1 0.2 0.5]}} {:center [1 0 -1] :radius 0.5 :material {:type :metal :attenuation [0.8 0.6 0.2] :fuzz 1}} {:center [-1 0 -1] :radius 0.5 :material {:type :dielectric :attenuation [1.0 1.0 1.0] :ref-idx 1.5}} {:center [-1 0 -1] :radius (- 0.45) :material {:type :dielectric :attenuation [1.0 1.0 1.0] :ref-idx 1.5}}]
This gives:
Positionable camera
Cameras, like dielectrics, are a pain to debug. So I always develop mine incrementally. First, let’s allow an adjustable field of view (fov). This is the angle you see through the portal. Since our image is not square, the fov is different horizontally and vertically. I always use vertical fov. I also usually specify it in degrees and change to radians inside the implementation— a matter of personal taste.
I first keep the rays coming from the origin and heading to the \(z = -1\) plane. We could make it the \(z = -2\) plane, or whatever, as long as we made \(h\) a ratio to that distance. Here is our setup:
This implies \(h = tan(\theta/2)\). Our camera now becomes:
(defn get-camera [vfov aspect] (let [theta (Math/toRadians vfov) half-height (Math/tan (/ theta 2)) half-width (* aspect half-height) lower-left-corner [(- half-width) (- half-height) -1] horizontal [(* 2 half-width) 0 0] vertical [0 (* 2 half-height) 0] origin [0 0 0]] {:lower-left-corner lower-left-corner :horizontal horizontal :vertical vertical :origin origin})) (defn get-ray [camera u v] (let [origin (:origin camera) horizontal (:horizontal camera) vertical (:vertical camera) lower-left-corner (:lower-left-corner camera) vector-on-plane (vec3-plus (vec3-scalar-multiply u horizontal) (vec3-scalar-multiply v vertical)) direction (vec3-plus lower-left-corner (vec3-minus vector-on-plane origin))] (create-ray origin direction))) (defn get-hitable-list [r] [{:center [(- r) 0 -1] :radius r :material {:type :diffuse :attenuation [0 0 1]}} {:center [r 0 -1] :radius r :material {:type :diffuse :attenuation [1 0 0]}}])
(let [nx 200 ny 100 ns 100 camera (get-camera 90 (/ nx ny)) hitable-list (get-hitable-list (Math/cos (/ Math/PI 4))) header (str "P3\n" nx " " ny "\n255\n") body (clojure.string/join (for [j (range (dec ny) -1 -1) i (range 0 nx) :let [pixel-color (average-color nx ny ns i j camera hitable-list) gamma-corrected (mapv #(Math/sqrt %) pixel-color) ir (int (* 255.99 (get gamma-corrected 0))) ig (int (* 255.99 (get gamma-corrected 1))) ib (int (* 255.99 (get gamma-corrected 2)))]] (str ir " " ig " " ib "\n")))] (spit "./output_images/image_17.ppm" header :append true) (spit "./output_images/image_17.ppm" body :append true))
gives:
To get an arbitrary viewpoint, let’s first name the points we care about. We’ll call the position where we place the camera lookfrom, and the point we look at lookat. (Later, if you want, you could define a direction to look in instead of a point to look at.)
We also need a way to specify the roll, or sideways tilt, of the camera; the rotation around the lookat-lookfrom axis. Another way to think about it is even if you keep lookfrom and lookat constant, you can still rotate your head around your nose. What we need is a way to specify an up vector for the camera. Notice we already we already have a plane that the up vector should be in, the plane orthogonal to the view direction.
We can actually use any up vector we want, and simply project it onto this plane to get an up vector for the camera. I use the common convention of naming a “view up” (vup) vector. A couple of cross products, and we now have a complete orthonormal basis \((u,v,w)\) to describe our camera’s orientation.
Remember that vup, \(v\), and \(w\) are all in the same plane. Note that, like before when our fixed camera faced \(-Z\), our arbitrary view camera faces \(-w\). And keep in mind that we can — but we don’t have to— use world up \((0,1,0)\) to specify vup. This is convenient and will naturally keep your camera horizontally level until you decide to experiment with crazy camera angles.
(defn get-camera [lookfrom lookat vup vfov aspect] (let [theta (Math/toRadians vfov) half-height (Math/tan (/ theta 2)) half-width (* aspect half-height) w (vec3-unit-vector (vec3-minus lookfrom lookat)) u (vec3-unit-vector (vec3-cross-product vup w)) v (vec3-cross-product w u) horizontal (vec3-scalar-multiply (* 2 half-width) u) vertical (vec3-scalar-multiply (* 2 half-height) v) origin lookfrom subtrahend (vec3-plus (vec3-plus (vec3-scalar-multiply half-width u) (vec3-scalar-multiply half-height v)) w) lower-left-corner (vec3-minus origin subtrahend)] {:lower-left-corner lower-left-corner :horizontal horizontal :vertical vertical :origin origin}))
This allows us to change the viewpoint:
(get-camera [-2 2 1] [0 0 -1] [0 1 0] 90 (/ nx ny))
to get:
And we can change field of view to get:
Defocus Blur
Now our final feature: defocus blur. Note, all photographers will call it “depth of field” so be aware of only using “defocus blur” among friends.
The reason we defocus blur in real cameras is because they need a big hole (rather than just a pinhole) to gather light. This would defocus everything, but if we stick a lens in the hole, there will be a certain distance where everything is in focus. The distance to that plane where things are in focus is controlled by the distance between the lens and the film/sensor. That is why you see the lens move relative to the camera when you change what is in focus (that may happen in your phone camera too, but the sensor moves). The “aperture” is a hole to control how big the lens is effectively. For a real camera, if you need more light you make the aperture bigger, and will get more defocus blur. For our virtual camera, we can have a perfect sensor and never need more light, so we only have an aperture when we want defocus blur.
A real camera has a compound lens that is complicated. For our code we could simulate the order: sensor, then lens, then aperture, and figure out where to send the rays and flip the image once computed (the image is projected upside down on the film). Graphics people usually use a thin lens approximation.
We also don’t need to simulate any of the inside of the camera. For the purposes of rendering an image outside the camera, that would be unnecessary complexity. Instead I usually start rays from the surface of the lens, and send them toward a virtual film plane, by finding the projection of the film on the plane that is in focus (at the distance \(focus\_dist\)).
For that we just need to have the ray origins be on a disk around lookfrom rather than from a point:
(defn random-in-unit-disk [] (loop [p (vec3-minus (vec3-scalar-multiply 2 [(rand) (rand) 0]) [1 1 0])] (if (< (vec3-dot-product p p) 1) p (recur (vec3-minus (vec3-scalar-multiply 2 [(rand) (rand) 0]) [1 1 0]))))) (defn get-camera [lookfrom lookat vup vfov aspect aperture focus_dist] (let [lens-radius (/ aperture 2) theta (Math/toRadians vfov) half-height (Math/tan (/ theta 2)) half-width (* aspect half-height) w (vec3-unit-vector (vec3-minus lookfrom lookat)) u (vec3-unit-vector (vec3-cross-product vup w)) v (vec3-cross-product w u) horizontal (vec3-scalar-multiply (* 2 half-width focus_dist) u) vertical (vec3-scalar-multiply (* 2 half-height focus_dist) v) origin lookfrom subtrahend (vec3-plus (vec3-plus (vec3-scalar-multiply (* focus_dist half-width) u) (vec3-scalar-multiply (* focus_dist half-height) v)) (vec3-scalar-multiply focus_dist w)) lower-left-corner (vec3-minus origin subtrahend)] {:lower-left-corner lower-left-corner :horizontal horizontal :vertical vertical :origin origin :lens-radius lens-radius :u u :v v})) (defn get-ray [camera s t] (let [origin (:origin camera) horizontal (:horizontal camera) vertical (:vertical camera) lower-left-corner (:lower-left-corner camera) lens-radius (:lens-radius camera) u (:u camera) v (:v camera) rd (vec3-scalar-multiply lens-radius (random-in-unit-disk)) offset (vec3-plus (vec3-scalar-multiply (get rd 0) u) (vec3-scalar-multiply (get rd 1) v)) vector-on-plane (vec3-plus (vec3-scalar-multiply s horizontal) (vec3-scalar-multiply t vertical)) direction (vec3-plus lower-left-corner (vec3-minus (vec3-minus vector-on-plane origin) offset))] (create-ray (vec3-plus origin offset) direction)))
Using a big aperture:
camera (let [lookfrom [3 3 2] lookat [0 0 -1] dist_to_focus (vec3-length (vec3-minus lookfrom lookat)) aperture 2] (get-camera lookfrom lookat [0 1 0] 20 (/ nx ny) aperture dist_to_focus))
We get:
