RayTracingTheNextWeek.html

<meta charset="utf-8">
<!-- Markdeep: https://casual-effects.com/markdeep/ -->



                                   **Ray Tracing: The Next Week**
                                           Peter Shirley
                                     Version 2.0.0, 2019-Oct-07
                    <br>Copyright 2018-2019. Peter Shirley. All rights reserved.



Overview
====================================================================================================

In Ray Tracing in One Weekend, you built a simple brute force path tracer. In this installment we’ll
add textures, volumes (like fog), rectangles, instances, lights, and support for lots of objects
using a BVH. When done, you’ll have a “real” ray tracer.

A heuristic in ray tracing that many people--including me--believe, is that most optimizations
complicate the code without delivering much speedup. What I will do in this mini-book is go with the
simplest approach in each design decision I make. Check https://in1weekend.blogspot.com/ for
readings and references to a more sophisticated approach. However, I strongly encourage you to do no
premature optimization; if it doesn’t show up high in the execution time profile, it doesn’t need
optimization until all the features are supported!

The two hardest parts of this book are the BVH and the Perlin textures. This is why the title
suggests you take a week rather than a weekend for this endeavor. But you can save those for last if
you want a weekend project. Order is not very important for the concepts presented in this book, and
without BVH and Perlin texture you will still get a Cornell Box!


Acknowledgments
---------------
Thanks to Becker for his many helpful comments on the draft and to Matthew Heimlich for spotting a
critical motion blur error. Thanks to Andrew Kensler, Thiago Ize, and Ingo Wald for advice on
ray-AABB tests. Thanks to David Hart and Grue Debry for help with a bunch of the details. Thanks to
Jean Buckley for editing. Thanks to Dan Drummond for code fixes. Thanks to Steve Hollasch and Trevor
David Black for getting the book translated to Markdeep and moved to the web.



Motion Blur
====================================================================================================

When you decided to ray trace, you decided visual quality was worth more run-time. In your fuzzy
reflection and defocus blur you needed multiple samples per pixel. Once you have taken a step down
that road, the good news is that almost all effects can be brute-forced. Motion blur is certainly
one of those. In a real camera, the shutter opens and stays open for a time interval, and the camera
and objects may move during that time. Its really an average of what the camera sees over that
interval that we want. We can get a random estimate by sending each ray at some random time when the
shutter is open. As long as the objects are where they should be at that time, we can get the right
average answer with a ray that is at exactly a single time. This is fundamentally why random ray
tracing tends to be simple.

The basic idea is to generate rays at random times while the shutter is open and intersect the model
at that one time. The way it is usually done is to have the camera move and the objects move, but
have each ray exist at exactly one time. This way the “engine” of the ray tracer can just make sure
the objects are where they need to be for the ray, and the intersection guts don’t change much.

<div class='together'>
For this we will first need to have a ray store the time it exists at:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class ray
    {
        public:
            ray() {}
            ray(const vec3& a, const vec3& b, float ti = 0.0) { A = a; B = b; _time = ti;}
            vec3 origin() const       { return A; }
            vec3 direction() const    { return B; }
            float time() const    { return _time; }
            vec3 point_at_parameter(float t) const { return A + t*B; }

            vec3 A;
            vec3 B;
            float _time;
    };

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

Now we need to modify the camera to generate rays at a random time between `time1` and `time2`.
Should the camera keep track of `time1` and `time2` or should that be up to the user of camera when
a ray is created? When in doubt, I like to make constructors complicated if it makes calls simple,
so I will make the camera keep track, but that’s a personal preference. Not many changes are needed
to camera because for now it is not allowed to move; it just sends out rays over a time period.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class camera {
        public:
            camera(vec3 lookfrom, vec3 lookat, vec3 vup,
                   float vfov /* vfov is top to bottom in degrees */,
                   float aspect, float aperture, float focus_dist,
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
                   float t0, float t1) {
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++

                time0 = t0;
                time1 = t1;
                lens_radius = aperture / 2;
                float theta = vfov*M_PI/180;
                float half_height = tan(theta/2);
                float half_width = aspect * half_height;
                origin = lookfrom;
                w = unit_vector(lookfrom - lookat);
                u = unit_vector(cross(vup, w));
                v = cross(w, u);
                lower_left_corner = origin
                                  - half_width*focus_dist*u
                                  - half_height*focus_dist*v
                                  - focus_dist*w;
                horizontal = 2*half_width*focus_dist*u;
                vertical = 2*half_height*focus_dist*v;
            }

            ray get_ray(float s, float t) {
                vec3 rd = lens_radius*random_in_unit_disk();
                vec3 offset = u * rd.x() + v * rd.y();
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
                float time = time0 + random_double()*(time1-time0);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
                return ray(
                    origin + offset,
                    lower_left_corner + s*horizontal + t*vertical - origin - offset,
                    time);
            }

            vec3 origin;
            vec3 lower_left_corner;
            vec3 horizontal;
            vec3 vertical;
            vec3 u, v, w;
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
            float time0, time1;
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
            float lens_radius;
    };

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We also need a moving object. I’ll create a sphere class that has its center move linearly from
`center0` at `time0` to `center1` at `time1`. Outside that time interval it continues on, so those
times need not match up with the camera aperture open close.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class moving_sphere: public hittable {
        public:
            moving_sphere() {}
            moving_sphere(vec3 cen0, vec3 cen1, float t0, float t1, float r, material *m)
              : center0(cen0), center1(cen1), time0(t0), time1(t1), radius(r), mat_ptr(m)
                {};
            virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const;
            vec3 center(float time) const;
            vec3 center0, center1;
            float time0, time1;
            float radius;
            material *mat_ptr;
    };

    vec3 moving_sphere::center(float time) const{
        return center0 + ((time - time0) / (time1 - time0))*(center1 - center0);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

<div class='together'>
An alternative to making a new moving sphere class is to just make them all move and have the
stationary ones have the same begin and end point. I’m on the fence about that trade-off between
fewer classes and more efficient stationary spheres, so let your design taste guide you. The
intersection code barely needs a change: `center` just needs to become a function `center(time)`:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool moving_sphere::hit(
        const ray& r, float t_min, float t_max, hit_record& rec) const {

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
        vec3 oc = r.origin() - center(r.time());
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
        float a = dot(r.direction(), r.direction());
        float b = dot(oc, r.direction());
        float c = dot(oc, oc) - radius*radius;
        float discriminant = b*b - a*c;
        if (discriminant > 0) {
            float temp = (-b - sqrt(discriminant))/a;
            if (temp < t_max && temp > t_min) {
                rec.t = temp;
                rec.p = r.point_at_parameter(rec.t);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
                rec.normal = (rec.p - center(r.time())) / radius;
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
                rec.mat_ptr = mat_ptr;
                return true;
            }
            temp = (-b + sqrt(discriminant))/a;
            if (temp < t_max && temp > t_min) {
                rec.t = temp;
                rec.p = r.point_at_parameter(rec.t);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
                rec.normal = (rec.p - center(r.time())) / radius;
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
                rec.mat_ptr = mat_ptr;
                return true;
            }
        }
        return false;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
Be sure that in the materials you have the scattered rays be at the time of the incident ray.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class lambertian : public material {
        public:
            lambertian(const vec3& a) : albedo(a) {}

            virtual bool scatter(const ray& r_in, const hit_record& rec,
                vec3& attenuation, ray& scattered) const {

                vec3 target = rec.p + rec.normal + random_in_unit_sphere();
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
                scattered = ray(rec.p, target-rec.p, r_in.time());
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
                attenuation = albedo;
                return true;
            }

        vec3 albedo;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
If we take the example diffuse spheres from scene at the end of the last book and make them move
from their centers at `time==0`, to `center + vec3(0, 0.5*random_double(), 0)` at `time==1`, with
the camera aperture open over that frame.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    hittable *random_scene() {
        int n = 50000;
        hittable **list = new hittable*[n+1];
        list[0] =  new sphere(vec3(0,-1000,0), 1000, new lambertian(checker));
        int i = 1;
        for (int a = -10; a < 10; a++) {
            for (int b = -10; b < 10; b++) {
                float choose_mat = random_double();
                vec3 center(a+0.9*random_double(),0.2,b+0.9*random_double());
                if ((center-vec3(4,0.2,0)).length() > 0.9) {
                    if (choose_mat < 0.8) {  // diffuse
                        list[i++] = new moving_sphere(
                            center,
                            center+vec3(0, 0.5*random_double(), 0),
                            0.0, 1.0, 0.2,
                            new lambertian(
                                vec3(random_double()*random_double(),
                                     random_double()*random_double(),
                                     random_double()*random_double())
                            )
                        );
                    }
                    else if (choose_mat < 0.95) { // metal
                        list[i++] = new sphere(
                            center, 0.2,
                            new metal(
                                vec3(0.5*(1 + random_double()),
                                     0.5*(1 + random_double()),
                                     0.5*(1 + random_double())),
                                0.5*random_double()
                            )
                        );
                    }
                    else {  // glass
                        list[i++] = new sphere(center, 0.2, new dielectric(1.5));
                    }
                }
            }
        }

        list[i++] = new sphere(vec3(0, 1, 0), 1.0, new dielectric(1.5));
        list[i++] = new sphere(vec3(-4, 1, 0), 1.0, new lambertian(vec3(0.4, 0.2, 0.1)));
        list[i++] = new sphere(vec3(4, 1, 0), 1.0, new metal(vec3(0.7, 0.6, 0.5), 0.0));

        return new hittable_list(list,i);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And with these viewing parameters gives:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    vec3 lookfrom(13,2,3);
    vec3 lookat(0,0,0);
    float dist_to_focus = 10.0;
    float aperture = 0.0;

    camera cam(
        lookfrom, lookat, vec3(0,1,0), 20, float(nx)/float(ny), aperture,
        dist_to_focus, 0.0, 1.0
    );
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  ![Image 2-1](../images/img-2-2-01.jpg)

</div>



Bounding Volume Hierarchies
====================================================================================================

This part is by far the most difficult and involved part of the ray tracer we are working on. I am
sticking it in this chapter so the code can run faster, and because it refactors `hittable` a
little, and when I add rectangles and boxes we won't have to go back and refactor them.

The ray-object intersection is the main time-bottleneck in a ray tracer, and the time is linear with
the number of objects. But it’s a repeated search on the same model, so we ought to be able to make
it a logarithmic search in the spirit of binary search. Because we are sending millions to billions
of rays on the same model, we can do an analog of sorting the model and then each ray intersection
can be a sublinear search. The two most common families of sorting are to 1) divide the space, and
2) divide the objects. The latter is usually much easier to code up and just as fast to run for most
models.

<div class='together'>
The key idea of a bounding volume over a set of primitives is to find a volume that fully encloses
(bounds) all the objects. For example, suppose you computed a bounding sphere of 10 objects. Any ray
that misses the bounding sphere definitely misses all ten objects. If the ray hits the bounding
sphere, then it might hit one of the ten objects. So the bounding code is always of the form:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if (ray hits bounding object)
        return whether ray hits bounded objects
    else
        return false
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

A key thing is we are dividing objects into subsets. We are not dividing the screen or the volume.
Any object is in just one bounding volume, but bounding volumes can overlap.

<div class='together'>
To make things sub-linear we need to make the bounding volumes hierarchical. For example, if we
divided a set of objects into two groups, red and blue, and used rectangular bounding volumes, we’d
have:

  ![Figure 3-1](../images/fig-2-3-01.jpg)

</div>

<div class='together'>
Note that the blue and red bounding volumes are contained in the purple one, but they might
overlap, and they are not ordered -- they are just both inside. So the tree shown on the right has
no concept of ordering in the left and right children; they are simply inside. The code would be:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if (hits purple)
        hit0 = hits blue enclosed objects
        hit1 = hits red enclosed objects
        if (hit0 or hit1)
            return true and info of closer hit
    return false
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

To get that all to work we need a way to make good divisions, rather than bad ones, and a way to
intersect a ray with a bounding volume. A ray bounding volume intersection needs to be fast, and
bounding volumes need to be pretty compact. In practice for most models, axis-aligned boxes work
better than the alternatives, but this design choice is always something to keep in mind if you
encounter unusual types of models.

From now on we will call axis-aligned bounding rectangular parallelepiped (really, that is what they
need to be called if precise) axis-aligned bounding boxes, or AABBs. Any method you want to use to
intersect a ray with an AABB is fine. And all we need to know is whether or not we hit it; we don’t
need hit points or normals or any of that stuff that we need for an object we want to display.

<div class='together'>
Most people use the “slab” method. This is based on the observation that an n-dimensional AABB is
just the intersection of n axis-aligned intervals, often called “slabs” An interval is just
the points between two endpoints, _e.g._, $x$ such that $3 < x < 5$, or more succinctly $x$ in
$(3,5)$. In 2D, two intervals overlapping makes a 2D AABB (a rectangle):

  ![Figure 3-2](../images/fig-2-3-02.jpg)

</div>

<div class='together'>
For a ray to hit one interval we first need to figure out whether the ray hits the boundaries. For
example, again in 2D, this is the ray parameters $t_0$ and $t_1$. (If the ray is parallel to the
plane those will be undefined.)

  ![Figure 3-3](../images/fig-2-3-03.jpg)

</div>

<div class='together'>
In 3D, those boundaries are planes. The equations for the planes are $x = x_0$, and $x = x_1$. Where
does the ray hit that plane? Recall that the ray can be thought of as just a function that given a
$t$ returns a location $p(t)$:

  $$ p(t) = A + t \cdot B $$
</div>

<div class='together'>
That equation applies to all three of the x/y/z coordinates. For example $x(t) = A_x + t \cdot B_x$.
This ray hits the plane $x = x_0$ at the $t$ that satisfies this equation:

  $$ x_0 = A_x + t_0 \cdot B_x $$
</div>

<div class='together'>
Thus $t$ at that hitpoint is:

  $$ t_0 = \frac{x_0 - A_x}{B_x} $$
</div>

<div class='together'>
We get the similar expression for $x_1$:

  $$ t_1 = \frac{x_1 - A_x}{B_x} $$
</div>

<div class='together'>
The key observation to turn that 1D math into a hit test is that for a hit, the $t$-intervals need
to overlap. For example, in 2D the green and blue overlapping only happens if there is a hit:

  ![Figure 3-4](../images/fig-2-3-04.jpg)

</div>

<div class='together'>
What “do the t intervals in the slabs overlap?” would like in code is something like:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    compute (tx0, tx1)
    compute (ty0, ty1)
    return overlap?( (tx0, tx1), (ty0, ty1))
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

That is awesomely simple, and the fact that the 3D version also works is why people love the
slab method:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    compute (tx0, tx1)
    compute (ty0, ty1)
    compute (tz0, tz1)
    return overlap?( (tx0, tx1), (ty0, ty1), (tz0, tz1))
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
There are some caveats that make this less pretty than it first appears. First, suppose the ray is
travelling in the negative $x$ direction. The interval $(t_{x0}, t_{x1})$ as computed above might be
reversed, _e.g._ something like $(7, 3)$. Second, the divide in there could give us infinities. And
if the ray origin is on one of the slab boundaries, we can get a `NaN`. There are many ways these
issues are dealt with in various ray tracers’ AABB. (There are also vectorization issues like SIMD
which we will not discuss here. Ingo Wald’s papers are a great place to start if you want to go the
extra mile in vectorization for speed.) For our purposes, this is unlikely to be a major bottleneck
as long as we make it reasonably fast, so let’s go for simplest, which is often fastest anyway!
First let’s look at computing the intervals:

  $$ t_{x0} = \frac{x_0 - A_x}{B_x} $$
  $$ t_{x1} = \frac{x_1 - A_x}{B_x} $$
</div>

<div class='together'>
One troublesome thing is that perfectly valid rays will have $B_x = 0$, causing division by zero.
Some of those rays are inside the slab, and some are not. Also, the zero will have a ± sign under
IEEE floating point. The good news for $B_x = 0$ is that $t_{x0}$ and $t_{x1}$ will both be +∞ or
both be -∞ if not between $x_0$ and $x_1$. So, using min and max should get us the right answers:

  $$ t_{x0} = min(\frac{x_0 - A_x}{B_x}, \frac{x_1 - A_x}{B_x}) $$
  $$ t_{x1} = max(\frac{x_0 - A_x}{B_x}, \frac{x_1 - A_x}{B_x}) $$
</div>

The remaining troublesome case if we do that is if $B_x = 0$ and either $x_0 - A_x = 0$ or $x_1-A_x
= 0$ so we get a `NaN`. In that case we can probably accept either hit or no hit answer, but we’ll
revisit that later.

<div class='together'>
Now, let’s look at that overlap function. Suppose we can assume the intervals are not reversed (so
the first value is less than the second value in the interval) and we want to return true in that
case. The boolean overlap that also computes the overlap interval $(f, F)$ of intervals $(d, D)$ and
$(e, E)$ would be:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    bool overlap(d, D, e, E, f, F)
        f = max(d, e)
        F = min(D, E)
        return (f < F)
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
If there are any `NaN`s running around there, the compare will return false so we need to be sure
our bounding boxes have a little padding if we care about grazing cases (and we probably should
because in a ray tracer all cases come up eventually). With all three dimensions in a loop and
passing in the interval $t_{min}$, $t_{max}$ we get:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    inline float ffmin(float a, float b) { return a < b ? a : b; }
    inline float ffmax(float a, float b) { return a > b ? a : b; }

    class aabb {
        public:
            aabb() {}
            aabb(const vec3& a, const vec3& b) { _min = a; _max = b;}

            vec3 min() const {return _min; }
            vec3 max() const {return _max; }

            bool hit(const ray& r, float tmin, float tmax) const {
                for (int a = 0; a < 3; a++) {
                    float t0 = ffmin((_min[a] - r.origin()[a]) / r.direction()[a],
                                     (_max[a] - r.origin()[a]) / r.direction()[a]);
                    float t1 = ffmax((_min[a] - r.origin()[a]) / r.direction()[a],
                                     (_max[a] - r.origin()[a]) / r.direction()[a]);
                    tmin = ffmax(t0, tmin);
                    tmax = ffmin(t1, tmax);
                    if (tmax <= tmin)
                        return false;
                }
                return true;
            }

            vec3 _min;
            vec3 _max;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
Note that the built-in `fmax()` is replaced by `ffmax()` which is quite a bit faster because it
doesn’t worry about `NaN`s and other exceptions. In reviewing this intersection method, Andrew
Kensler at Pixar tried some experiments and has proposed this version of the code which works
extremely well on many compilers, and I have adopted it as my go-to method:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    inline bool aabb::hit(const ray& r, float tmin, float tmax) const {
        for (int a = 0; a < 3; a++) {
            float invD = 1.0f / r.direction()[a];
            float t0 = (min()[a] - r.origin()[a]) * invD;
            float t1 = (max()[a] - r.origin()[a]) * invD;
            if (invD < 0.0f)
                std::swap(t0, t1);
            tmin = t0 > tmin ? t0 : tmin;
            tmax = t1 < tmax ? t1 : tmax;
            if (tmax <= tmin)
                return false;
        }
        return true;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
We now need to add a function to compute bounding boxes to all of the hittables. Then we will make a
hierarchy of boxes over all the primitives and the individual primitives, like spheres, will live at
the leaves. That function returns a bool because not all primitives have bounding boxes (_e.g._,
infinite planes). In addition, objects move so it takes `time1` and `time2` for the interval of the
frame and the bounding box will bound the object moving through that interval.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class hittable {
        public:
            virtual bool hit(
                const ray& r, float t_min, float t_max, hit_record& rec) const = 0;
            virtual bool bounding_box(float t0, float t1, aabb& box) const = 0;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
For a sphere, that `bounding_box` function is easy:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool sphere::bounding_box(float t0, float t1, aabb& box) const {
        box = aabb(center - vec3(radius, radius, radius),
                   center + vec3(radius, radius, radius));
        return true;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
For `moving sphere`, we can take the box of the sphere at $t_0$, and the box of the sphere at $t_1$,
and compute the box of those two boxes:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool moving_sphere::bounding_box(float t0, float t1, aabb& box) const {
        aabb box0(center(t0) - vec3(radius, radius, radius),
                  center(t0) + vec3(radius, radius, radius));
        aabb box1(center(t1) - vec3(radius, radius, radius),
                  center(t1) + vec3(radius, radius, radius));
        box = surrounding_box(box0, box1);
        return true;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
For lists you can store the bounding box at construction, or compute it on the fly. I like doing it
the fly because it is only usually called at BVH construction.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool hittable_list::bounding_box(float t0, float t1, aabb& box) const {
        if (list_size < 1) return false;
        aabb temp_box;
        bool first_true = list[0]->bounding_box(t0, t1, temp_box);
        if (!first_true)
            return false;
        else
            box = temp_box;
        for (int i = 1; i < list_size; i++) {
            if(list[i]->bounding_box(t0, t1, temp_box)) {
                box = surrounding_box(box, temp_box);
            }
            else
                return false;
        }
        return true;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
This requires the `surrounding_box` function for `aabb` which computes the bounding box of two
boxes:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    aabb surrounding_box(aabb box0, aabb box1) {
        vec3 small( ffmin(box0.min().x(), box1.min().x()),
                    ffmin(box0.min().y(), box1.min().y()),
                    ffmin(box0.min().z(), box1.min().z()));
        vec3 big  ( ffmax(box0.max().x(), box1.max().x()),
                    ffmax(box0.max().y(), box1.max().y()),
                    ffmax(box0.max().z(), box1.max().z()));
        return aabb(small,big);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
A BVH is also going to be a `hittable` -- just like lists of `hittable`s. It’s really a container,
but it can respond to the query “does this ray hit you?”. One design question is whether we have two
classes, one for the tree, and one for the nodes in the tree; or do we have just one class and have
the root just be a node we point to. I am a fan of the one class design when feasible. Here is such
a class:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class bvh_node : public hittable {
        public:
            bvh_node() {}
            bvh_node(hittable **l, int n, float time0, float time1);

            virtual bool hit(const ray& r, float tmin, float tmax, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const;

            hittable *left;
            hittable *right;
            aabb box;
    };

    bool bvh_node::bounding_box(float t0, float t1, aabb& b) const {
        b = box;
        return true;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

Note that the children pointers are to generic hittables. They can be other `bvh_nodes`, or
`spheres`, or any other `hittable`.

<div class='together'>
The `hit` function is pretty straightforward: check whether the box for the node is hit, and if so,
check the children and sort out any details:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool bvh_node::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
        if (box.hit(r, t_min, t_max)) {
            hit_record left_rec, right_rec;
            bool hit_left = left->hit(r, t_min, t_max, left_rec);
            bool hit_right = right->hit(r, t_min, t_max, right_rec);
            if (hit_left && hit_right) {
                if (left_rec.t < right_rec.t)
                    rec = left_rec;
                else
                    rec = right_rec;
                return true;
            }
            else if (hit_left) {
                rec = left_rec;
                return true;
            }
            else if (hit_right) {
                rec = right_rec;
                return true;
            }
            else
                return false;
        }
        else return false;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
The most complicated part of any efficiency structure, including the BVH, is building it. We do this
in the constructor. A cool thing about BVHs is that as long as the list of objects in a `bvh_node`
gets divided into two sub-lists, the hit function will work. It will work best if the division is
done well, so that the two children have smaller bounding boxes than their parent’s bounding box,
but that is for speed not correctness. I’ll choose the middle ground, and at each node split the
list along one axis. I’ll go for simplicity:

1. randomly choose an axis
2. sort the primitives using library qsort
3. put half in each subtree

</div>

I used the old-school C `qsort` rather than the C++ sort because I need a different compare operator
depending on axis, and `qsort` takes a compare function rather than using the less-than operator. I
pass in a pointer to pointer -- this is just C for “array of pointers” because a pointer in C can
also just be a pointer to the first element of an array.

<div class='together'>
When the list coming in is two elements, I put one in each subtree and end the recursion. The
traverse algorithm should be smooth and not have to check for null pointers, so if I just have one
element I duplicate it in each subtree. Checking explicitly for three elements and just following
one recursion would probably help a little, but I figure the whole method will get optimized later.
This yields:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bvh_node::bvh_node(hittable **l, int n, float time0, float time1) {
        int axis = int(3*random_double());

        if (axis == 0)
           qsort(l, n, sizeof(hittable *), box_x_compare);
        else if (axis == 1)
           qsort(l, n, sizeof(hittable *), box_y_compare);
        else
           qsort(l, n, sizeof(hittable *), box_z_compare);

        if (n == 1) {
            left = right = l[0];
        }
        else if (n == 2) {
            left = l[0];
            right = l[1];
        }
        else {
            left = new bvh_node(l, n/2, time0, time1);
            right = new bvh_node(l + n/2, n - n/2, time0, time1);
        }

        aabb box_left, box_right;

        if (!left->bounding_box(time0, time1, box_left) ||
            !right->bounding_box(time0, time1, box_right)) {

            std::cerr << "no bounding box in bvh_node constructor\n";
        }

        box = surrounding_box(box_left, box_right);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

The check for whether there is a bounding box at all is in case you sent in something like an
infinite plane that doesn’t have a bounding box. We don’t have any of those primitives, so it
shouldn’t happen until you add such a thing.

The compare function has to take void pointers which you cast. This is old-school C and reminded me
why C++ was invented. I had to really mess with this to get all the pointer junk right. If you like
this part, you have a future as a systems person!

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    int box_x_compare (const void * a, const void * b) {
        aabb box_left, box_right;
        hittable *ah = *(hittable**)a;
        hittable *bh = *(hittable**)b;

        if (!ah->bounding_box(0,0, box_left) || !bh->bounding_box(0,0, box_right))
            std::cerr << "no bounding box in bvh_node constructor\n";

        if (box_left.min().x() - box_right.min().x() < 0.0)
            return -1;
        else
            return 1;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~



Solid Textures
====================================================================================================

A texture in graphics usually means a function that makes the colors on a surface procedural. This
procedure can be synthesis code, or it could be an image lookup, or a combination of both. We will
first make all colors a texture. Most programs keep constant rgb colors and textures different
classes so feel free to do something different, but I am a big believer in this architecture because
being able to make any color a texture is great.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class texture (
        public:
            virtual vec3 value(float u, float v, const vec3& p) const = O;
    };

    class constant_texture : public texture {
        public:
            constant_texture() {}
            constant_texture(vec3 c) : color(c) {}
            virtual vec3 value(float u, float v, const vec3& p) const {
                return color;
            }
            vec3 color;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

<div class='together'>
Now we can make textured materials by replacing the vec3 color with a texture pointer:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class lambertian : public material {
        public:
            lambertian(texture *a) : albedo(a) {}
            virtual bool scatter(const ray& r_in, const hit_record& rec,
                                 vec3& attenuation, ray& scattered) const {
                vec3 target = rec.p + rec.normal + random_in_unit_sphere();
                scattered = ray(rec.p, target - rec.p);
                attenuation = albedo->value(0, 0, rec.p);
                return true;
            }
            texture *albedo;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
where you used to have

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    new lambertian(vec3(0.5, 0.5, 0.5)))
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

now you should replace the `vec3(...)` with `new constant_texture(vec3(...))`

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    new lambertian(new constant_texture(vec3(0.5, 0.5, 0.5))))
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
We can create a checker texture by noting that the sign of sine and cosine just alternates in a
regular way and if we multiply trig functions in all three dimensions, the sign of that product
forms a 3D checker pattern.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class checker_texture : public texture {
        public:
            checker_texture() {}
            checker_texture(texture *t0, texture *tl): even(t0), odd(t1) {}
            virtual vec3 value(float u, float v, const vec3& p) const {
                float sines = sin(10*p.x())*sin(10*p.y())*sin(10*p.z());
                if (sines < 0)
                    return odd—>value(u, v, p);
                else
                    return even->value(u, v, p);
            }
            texture *odd;
            texture *even;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

Those checker odd/even pointers can be to a constant texture or to some other procedural texture.
This is in the spirit of shader networks introduced by Pat Hanrahan back in the 1980s.

<div class='together'>
If we add this to our random_scene() function’s base sphere:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    texture *checker = new checker_texture(
        new constant_texture(vec3(0.2, 0.3, 0.1)),
        new constant_texture(vec3(0.9, 0.9, 0.9))
    );
    list[0] = new sphere(vec3(0,-1000,0), 1000, new lambertian(checker));
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We get:

  ![Image 4-1](../images/img-2-4-01.jpg)

</div>

<div class='together'>
If we add a new scene:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    hittable *two_spheres() {
        texture *checker = new checker_texture(
            new constant_texture(vec3(0.2, 0.3, 0.1)),
            new constant_texture(vec3(0.9, 0.9, 0.9))
        );
        int n = 50;
        hittable **list = new hitable*[n+1];
        list[0] = new sphere(vec3(0,-10, 0), 10, new lambertian(checker));
        1ist[1] = new sphere(vec3(0, 10, 0), 10, new lambertian(checker));
        return new hittable_list(list,2);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

With camera:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    vec3 lookfrom(13,2,3);
    vec3 lookat(0,0,0);
    float dist_to_focus = 10.0;
    float aperture = 0.0;

    camera cam(lookfrom, lookat, vec3(0,1,0), 20, float(nx)/float(ny),
               aperture, dist_to_focus, 0.0, 1.0);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We get:

  ![Image 4-2](../images/img-2-4-02.jpg)

</div>



Perlin Noise
====================================================================================================

<div class='together'>
To get cool looking solid textures most people use some form of Perlin noise. These are named after
their inventor Ken Perlin. Perlin texture doesn’t return white noise like this:

  ![Image 5-1](../images/img-2-5-01.jpg)

Instead it returns something similar to blurred white noise:

  ![Image 5-2](../images/img-2-5-02.jpg)

</div>

A key part of Perlin noise is that it is repeatable: it takes a 3D point as input and always returns
the same randomish number. Nearby points return similar numbers. Another important part of Perlin
noise is that it be simple and fast, so it’s usually done as a hack. I’ll build that hack up
incrementally based on Andrew Kensler’s description.

<div class='together'>
We could just tile all of space with a 3D array of random numbers and use them in blocks. You get
something blocky where the repeating is clear:

  ![Image 5-3](../images/img-2-5-03.jpg)

</div>

<div class='together'>
Let’s just use some sort of hashing to scramble this, instead of tiling. This has a bit of support
code to make it all happen:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class perlin {
        public:
            float noise(const vec3& p) const {
                float u = p.x() - floor(p.x());
                float v = p.y() - floor(p.y());
                float w = p.z() - floor(p.z());
                int i = floor(p.x());
                int j = floor(p.y());
                int k = floor(p.z());
                return ranfloat[perm_x[i] ^ perm_y[j] ^ perm_z[k]];
            }
            static float *ranfloat;
            static int *perm_x;
            static int *perm_y;
            static int *perm_z;
    };

    static vec3* perlin_generate() {
        float * p = new float[256];
        for (int i = 0; i < 256; ++i)
            p[i] = random_double();
        return p;
    }

    void permute(int *p, int n) {
        for (int i = n-1; i > 0; i--) {
            int target = int(random_double()*(i+1));
            int tmp = p[i];
            p[i] = p[target];
            p[target] = tmp;
        }
        return;
    }

    static int* perlin_generate_perm() {
        int * p = new int[256];
        for (int i = 0; i < 256; i++)
            p[i] = i;
        permute(p, 256);
        return p;
    }

    float *perlin::ranfloat = perlin_generate();
    int *perlin::perm_x = perlin_generate_perm();
    int *perlin::perm_y = perlin_generate_perm();
    int *perlin::perm_z = perlin_generate_perm();
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
Now if we create an actual texture that takes these floats between 0 and 1 and creates grey
colors:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class noise_texture : public texture {
        public:
            noise_texture() {}
            virtual vec3 value(float u, float v, const vec3& p) const {
                return vec3(1,1,1) * noise.noise(p);
            }
            perlin noise;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
We can use that texture on some spheres:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    hittable *two_perlin_spheres() {
        texture *pertext = new noise_texture();
        hittable **list = new hittable*[2];
        list[0] = new sphere(vec3(0,-1000, 0), 1000, new lambertian(pertext));
        list[1] = new sphere(vec3(0, 2, 0), 2, new lambertian(pertext));
        return new hittable_list(list, 2);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
With the same camera as before:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    vec3 lookfrom(13,2,3);
    vec3 lookat(0,0,0);
    float dist_to_focus = 10.0;
    float aperture = 0.0;
    camera cam(lookfrom, lookat, vec3(0,1,0), 20, float(nx)/float(ny),
               aperture, dist_to_focus, 0.0, 1.0);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
Add the hashing does scramble as hoped:

  ![Image 5-4](../images/img-2-5-04.jpg)

</div>

<div class='together'>
To make it smooth, we can linearly interpolate:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    inline float trilinear_interp(float c[2][2][2], float u, float v, float w) {
        float accum = 0;
        for (int i=0; i < 2; 1++)
            for (int j=0; j < 2; j++)
                for (int k=0; k < 2; k++)
                    accum += (i*u + (1-i)*(1-u))*
                             (j*v + (1-j)*(1-v))*
                             (k*w + (1-k)*(1-w))*c[i][j][k];

        return accum;
    }

    class perlin {
        public:
            float noise(const vec3& p) const {
                float u = p.x() - floor(p.x());
                float v = p.y() - floor(p.y());
                float w = p.z() - floor(p.z());
                int i = floor(p.x());
                int j = floor(p.y());
                int k = floor(p.z());
                float c[2][2][2];
                for (int di=0; di < 2; di++)
                    for (int dj=0; dj < 2; dj++)
                        for (int dk=0; dk < 2; dk++)
                            c[di][dj][dk] = ranfloat[
                                perm_x[(i+di) & 255] ^
                                perm_y[(j+dj) & 255] ^
                                perm_z[(k+dk) & 255]
                            ];
                return trilinear_interp(c, u, v, w);
            }
            static int *perm_x;
            static int *perm_y;
            static int *perm_z;
        }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And we get:

  ![Image 5-5](../images/img-2-5-05.jpg)

</div>

<div class='together'>
Better, but there are obvious grid features in there. Some of it is Mach bands, a known perceptual
artifact of linear interpolation of color. A standard trick is to use a hermite cubic to round off
the interpolation:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class perlin (
        public:
            float noise(const vec3& p) const {
                float u = p.x() - floor(p.x());
                float v = p.y() - floor(p.y());
                float w = p.z() - floor(p.z());

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
                u = u*u*(3-2*u);
                v = v*v*(3-2*v);
                w = w*w*(3-2*w);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++

                int i = floor(p.x());
                int j = floor(p.y());
                int k = floor(p.z());
                ...
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
This gives a smoother looking image:

  ![Image 5-6](../images/img-2-5-06.jpg)

</div>

<div class='together'>
It is also a bit low frequency. We can scale the input point to make it vary more quickly:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class noise_texture : public texture {
        public:
            noise_texture() {}
            noise_texture(float sc) : scale(sc) {}
            virtual vec3 value(float u, float v, const vec3& p) const {
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ none delete
                return vec3(1,1,1) * noise.noise(p);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
                return vec3(1,1,1) * noise.noise(scale * p);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
            }
            perlin noise;
            float scale;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

which gives:

  ![Image 5-7](../images/img-2-5-07.jpg)

</div>

<div class='together'>
This is still a bit grid blocky looking, probably because the min and max of the pattern always
lands exactly on the integer x/y/z. Ken Perlin’s very clever trick was to instead put random unit
vectors (instead of just floats) on the lattice points, and use a dot product to move the min and
max off the lattice. So, first we need to change the random floats to random vectors:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    vec3 *perlin::ranvec = perlin_generate();
    int *perlin::perm_x = perlin_generate_perm();
    int *perlin::perm_y = perlin_generate_perm();
    int *perlin::perm_z = perlin_generate_perm();
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
These vectors are any reasonable set of irregular directions, and I won't bother to make them
exactly uniform:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    static vec3* perlin_generate() {
        vec3 *p = new vec3[256];
        for (int i = 0; i < 256; ++i) {
            double x_random = 2*random_double() - 1;
            double y_random = 2*random_double() - 1;
            double z_random = 2*random_double() - 1;
            p[i] = unit_vector(vec3(x_random, y_random, z_random));
        }
        return p;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
The Perlin class is now:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class perlin {
        public:
            float noise(const vec3& p) const {
                float u = p.x() - floor(p.x());
                float v = p.y() - floor(p.y());
                float w = p.z() - floor(p.z());
                int i = floor(p.x());
                int j = floor(p.y());
                int k = floor(p.z());
                vec3 c[2][2][2];
                for (int di=0; di < 2; di++)
                    for (int dj=0; dj < 2; dj++)
                        for (int dk=0; dk < 2; dk++)
                            c[di][dj][dk] = ranvec[
                                perm_x[(i+di) & 255] ^
                                perm_y[(j+dj) & 255] ^
                                pexm_z[(k+dk) & 255]
                            ];
                return perlin_interp(c, u, v, w);
            }
            static vec3 *ranvec;
            static int *perm_x;
            static int *perm_y;
            static int *perm_z;
        }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And the interpolation becomes a bit more complicated:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    inline float perlin_interp(vec3 c[2][2][2], float u, float v, float w) {
        float uu = u*u*(3-2*u);
        float vv = v*v*(3-2*v);
        float ww = w*w*(3-2*w);
        float accum = 0;
        for (int i=0; i < 2; i++)
            for (int j=0; j < 2; j++)
                for (int k=0; k < 2; k++) {
                    vec3 weight_v(u-i, v-j, w-k);
                    accum += (i*uu + (1-i)*(1-uu))*
                        (j*vv + (1-j)*(1-vv))*
                        (k*ww + (1-k)*(1-ww))*dot(c[i][j][k], weight_v);
                }
        return accum;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
This finally gives something more reasonable looking:

  ![Image 5-8](../images/img-2-5-08.jpg)

</div>

<div class='together'>
Very often, a composite noise that has multiple summed frequencies is used. This is usually called
turbulence and is a sum of repeated calls to noise:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    float turb(const vec3& p, int depth=7) const {
        float accum = 0;
        vec3 temp_p = p;
        float weight = 1.0;
        for (int i = 0; i < depth; i++) {
            accum += weight*noise(temp_p);
            weight *= 0.5;
            temp_p *= 2;
        }
        return fabs(accum);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

Here `fabs()` is the `math.h` absolute value function.

<div class='together'>
Used directly, turbulence gives a sort of camouflage netting appearance:

  ![Image 5-9](../images/img-2-5-09.jpg)

</div>

<div class='together'>
However, usually turbulence is used indirectly. For example, the “hello world” of procedural solid
textures is a simple marble-like texture. The basic idea is to make color proportional to something
like a sine function, and use turbulence to adjust the phase (so it shifts $x$ in $sin( x )$) which
makes the stripes undulate. Commenting out straight noise and turbulence, and giving a marble-like
effect is:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class noise_texture : public texture {
        public:
            noise_texture() {}
            noise_texture(float sc) : scale(sc) {}
            virtual vec3 value(float u, float v, const vec3& p) const {
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ none delete
                return vec3(1,1,1) * noise.noise(scale * p);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++ highlight
                return vec3(1,1,1) * 0.5 * (1 + sin(scale*p.z() + 10*noise.turb(p)));
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
            }
            perlin noise;
            float scale;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Which yields:

  ![Image 5-10](../images/img-2-5-10.jpg)

</div>



Image Texture Mapping
====================================================================================================

We used the hitpoint p before to index a procedure solid texture like marble. We can also read in an
image and use a 2D $(u,v)$ texture coordinate to index into the image.

<div class='together'>
A direct way to use scaled $(u,v)$ in an image is to round the u and v to integers, and use that as
$(i,j)$ pixels. This is awkward, because we don’t want to have to change the code when we change
image resolution. So instead, one of the the most universal unofficial standards in graphics is to
use texture coordinates instead of image pixel coordinates. These are just some form of fractional
position in the image. For example, for pixel $(i,j)$ in an nx by ny image, the image texture
position is:

  $$ u = i/(nx-1) $$
  $$ v = j/(nx-1) $$
</div>

<div class='together'>
This is just a fractional position. For a hittable, we need to also return a u and v in the hit
record. For spheres, this is usually based on some form of longitude and latitude, _i.e._, spherical
coordinates. So if we have a $(\theta,\phi)$ in spherical coordinates we just need to scale $\theta$
and $\phi$ to fractions. If $\theta$ is the angle down from the pole, and $\phi$ is the angle around
the axis through the poles, the normalization to $[0,1]$ would be:

  $$ u = \phi / (2\pi) $$
  $$ v = \theta / \pi $$
</div>

<div class='together'>
To compute $\theta$ and $\phi$, for a given hitpoint, the formula for spherical coordinates of a
unit radius sphere on the origin is:

  $$ x = cos(\phi) cos(\theta) $$
  $$ y = sin(\phi) cos(\theta) $$
  $$ z = sin(\theta) $$
</div>

<div class='together'>
We need to invert that. Because of the lovely `math.h` function `atan2()` which takes any number
proportional to sine and cosine and returns the angle, we can pass in $x$ and $y$ (the $cos(\theta)$
cancel):

  $$ \phi = atan2(y, x) $$
</div>

<div class='together'>
The $atan2$ returns in the range $-\pi$ to $\pi$ so we need to take a little care there.
The $\theta$ is more straightforward:

  $$ \theta = asin(z) $$

which returns numbers in the range $-\pi/2$ to $\pi/2$.
</div>

<div class='together'>
So for a sphere, the $(u,v)$ coord computation is accomplished by a utility function that expects
things on the unit sphere centered at the origin. The call inside sphere::hit should be:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    get_sphere_uv((rec.p-center)/radius, rec.u, rec.v);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
The utility function is:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    void get_sphere_uv(const vec3& p, float& u, float& v) {
        float phi = atan2(p.z(), p.x());
        float theta = asin(p.y());
        u = 1-(phi + M_PI) / (2*M_PI);
        v = (theta + M_PI/2) / M_PI;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

Now we also need to create a texture class that holds an image. I am going to use my favorite image
utility `stb_image`. It reads in an image into a big array of unsigned char. These are just packed
RGBs that each range 0..255 for black to fully-on.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class image_texture : public texture {
        public:
            image_texture() {}
            image_texture(unsigned char *pixels, int A, int B)
                : data(pixels), nx(A), ny(B) {}
            virtual vec3 value(float u, float v, const vec3& p) const;
            unsigned char *data;
            int nx, ny;
    };

    vec3 image_texture::value(float u, float v, const vec3& p) const {
         int i = (  u) * nx;
         int j = (1-v) * ny - 0.001;
         if (i < 0) i = 0;
         if (j < 0) j = 0;
         if (i > nx-1) i = nx-1;
         if (j > ny-1) j = ny-1;
         float r = int(data[3*i + 3*nx*j]  ) / 255.0;
         float g = int(data[3*i + 3*nx*j+1]) / 255.0;
         float b = int(data[3*i + 3*nx*j+2]) / 255.0;
         return vec3(r, g, b);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

<div class='together'>
The representation of a packed array in that order is pretty standard. Thankfully, the `stb_image`
package makes that super simple -- just include the header in `main.h` with a `#define`:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    #define STB_IMAGE_IMPLEMENTATION
    #include "stb_image.h"
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>

  ![Image 6-1: earthmap.jpg](../images/earthmap.jpg)

To read an image from a file earthmap.jpg (I just grabbed a random earth map from the web -- any
standard projection will do for our purposes), and then assign it to a diffuse material, the code
is:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    int nx, ny, nn;
    unsigned char *tex_data = stbi_load("earthmap.jpg", &nx, &ny, &nn, 0);
    material *mat = new lambertian(new image_texture(tex_data, nx, ny));
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

We start to see some of the power of all colors being textures -- we can assign any kind of texture
to the lambertian material, and lambertian doesn’t need to be aware of it.

<div class='together'>
To test this, assign it to a sphere, and then temporarily cripple the color() function in main to
just return attenuation. You should get something like:

  ![Image 6-2](../images/img-2-6-02.jpg)

</div>



Rectangles and Lights
====================================================================================================

<div class='together'>
First, let’s make a light emitting material. We need to add an emitted function (we could also add
it to hit_record instead -- that’s a matter of design taste). Like the background, it just tells the
ray what color it is and performs no reflection. It’s very simple:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class diffuse_light : public material {
        public:
            diffuse_light(texture *a) : emit(a) {}
            virtual bool scatter(const ray& r_in, const hit_record& rec,
                vec3& attenuation, ray& scattered) const { return false; }
            virtual vec3 emitted(float u, float v, const vec3& p) const {
                return emit->value(u, v, p);
            }
            texture *emit;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
So that I don’t have to make all the non-emitting materials implement `emitted()`, I have the base
class return black:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class material {
        public:
            virtual bool scatter(const ray& r_in, const hit_record& rec,
                vec3& attenuation, ray& scattered) const = 0;
            virtual vec3 emitted(float u, float v, const vec3& p) const {
                return vec3(0,0,0);
            }
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
Next, let’s make the background black in our color function, and pay attention to emitted:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    vec3 color(const ray& r, hittable *world, int depth) {
        hit_record rec;
        if (world->hit(r, 0.001, MAXFLOAT, rec)) {
            ray scattered;
            vec3 attenuation;
            vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
            if (depth < 50 && rec.mat_ptr->scatter(r, rec, attenuation, scattered))
                 return emitted + attenuation*color(scattered, world, depth+1);
            else
                return emitted;
        }
        else
            return vec3(0,0,0);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

Now, let’s make some rectangles. Rectangles are often convenient for modelling man-made
environments. I’m a fan of doing axis-aligned rectangles because they are easy. (We’ll get to
instancing so we can rotate them later.)

<div class='together'>
First, here is a rectangle in an xy plane. Such a plane is defined by its z value. For example, $z =
k$. An axis-aligned rectangle is defined by lines $x=x_0$ , $x=x_1$ , $y=y_0$ , $y=y_1$.

  ![Figure 7-1](../images/fig-2-7-01.jpg)

</div>

<div class='together'>
To determine whether a ray hits such a rectangle, we first determine where the ray hits the plane.
Recall that a ray $p(t) = a + t \cdot b$ has its z component defined by $z(t) = a_z + t \cdot b_z$.
Rearranging those terms we can solve for what the t is where $z=k$.

  $$ t = (k-a_z) / b_z $$
</div>

<div class='together'>
Once we have t, we can plug that into the equations for x and y:

  $$ x = a_x + t*b_x $$
  $$ y = a_y + t*b_y $$
</div>

It is a hit if $x_0 < x < x_1$ and $y_0 < y < y_1$.

<div class='together'>
The actual `xy_rect` class is thus:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class xy_rect: public hittable {
        public:
            xy_rect() {}
            xy_rect(float _x0, float _x1, float _y0, float _y1, float _k, material *mat)
                : x0(_x0), x1(_x1), y0(_y0), y1(_y1), k(_k), mp(mat) {};
            virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const {
                box =  aabb(vec3(x0,y0, k-0.0001), vec3(x1, y1, k+0.0001));
                return true;
            }
            material *mp;
            float x0, x1, y0, y1, k;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And the hit function is:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool xy_rect::hit(const ray& r, float t0, float t1, hit_record& rec) const {
        float t = (k-r.origin().z()) / r.direction().z();
        if (t < t0 || t > t1)
            return false;
        float x = r.origin().x() + t*r.direction().x();
        float y = r.origin().y() + t*r.direction().y();
        if (x < x0 || x > x1 || y < y0 || y > y1)
            return false;
        rec.u = (x-x0)/(x1-x0);
        rec.v = (y-y0)/(y1-y0);
        rec.t = t;
        rec.mat_ptr = mp;
        rec.p = r.point_at_parameter(t);
        rec.normal = vec3(0, 0, 1);
        return true;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
If we set up a rectangle as a light:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    hittable *simple_light() {
        texture *pertext = new noise_texture(4);
        hittable **list = new hittable*[4];
        list[0] = new sphere(vec3(0,-1000, 0), 1000, new lambertian(pertext));
        list[1] = new sphere(vec3(0, 2, 0), 2, new lambertian(pertext));
        list[2] = new sphere(vec3(0, 7, 0), 2,
            new diffuse_light(new constant_texture(vec3(4,4,4))));
        list[3] = new xy_rect(3, 5, 1, 3, -2,
            new diffuse_light(new constant_texture(vec3(4,4,4))));
        return new hittable_list(list,4);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
We get:

  ![Image 7-1](../images/img-2-7-01.jpg)

</div>

Note that the light is brighter than $(1,1,1)$. This allows it to be bright enough to light things.

<div class='together'>
Fool around with making some spheres lights too.

  ![Image 7-2](../images/img-2-7-02.jpg)

</div>

Now let’s add the other two axes and the famous Cornell Box.

<div class='together'>
This is xz and yz:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class xz_rect: public hittable {
        public:
            xz_rect() {}
            xz_rect(float _x0, float _x1, float _z0, float _z1, float _k, material *mat)
                : x0(_x0), x1(_x1), z0(_z0), z1(_z1), k(_k), mp(mat) {};
            virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const {
                box =  aabb(vec3(x0,k-0.0001,z0), vec3(x1, k+0.0001, z1));
                return true;
            }
            material *mp;
            float x0, x1, z0, z1, k;
    };

    class yz_rect: public hittable {
        public:
            yz_rect() {}
            yz_rect(float _y0, float _y1, float _z0, float _z1, float _k, material *mat)
                : y0(_y0), y1(_y1), z0(_z0), z1(_z1), k(_k), mp(mat) {};
            virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const {
                box =  aabb(vec3(k-0.0001, y0, z0), vec3(k+0.0001, y1, z1));
                return true;
            }
            material  *mp;
            float y0, y1, z0, z1, k;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
With unsurprising hit functions:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool xz_rect::hit(const ray& r, float t0, float t1, hit_record& rec) const {
        float t = (k-r.origin().y()) / r.direction().y();
        if (t < t0 || t > t1)
            return false;
        float x = r.origin().x() + t*r.direction().x();
        float z = r.origin().z() + t*r.direction().z();
        if (x < x0 || x > x1 || z < z0 || z > z1)
            return false;
        rec.u = (x-x0)/(x1-x0);
        rec.v = (z-z0)/(z1-z0);
        rec.t = t;
        rec.mat_ptr = mp;
        rec.p = r.point_at_parameter(t);
        rec.normal = vec3(0, 1, 0);
        return true;
    }

    bool yz_rect::hit(const ray& r, float t0, float t1, hit_record& rec) const {
        float t = (k-r.origin().x()) / r.direction().x();
        if (t < t0 || t > t1)
            return false;
        float y = r.origin().y() + t*r.direction().y();
        float z = r.origin().z() + t*r.direction().z();
        if (y < y0 || y > y1 || z < z0 || z > z1)
            return false;
        rec.u = (y-y0)/(y1-y0);
        rec.v = (z-z0)/(z1-z0);
        rec.t = t;
        rec.mat_ptr = mp;
        rec.p = r.point_at_parameter(t);
        rec.normal = vec3(1, 0, 0);
        return true;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
Let’s make the 5 walls and the light of the box:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    hittable *cornell_box() {
        hittable **list = new hitable*[5];
        int i = 0;
        material *red = new lambertian(new constant_texture(vec3(0.65, 0.05, 0.05)));
        material *white = new lambertian(new constant_texture(vec3(0.73, 0.73, 0.73)));
        material *green = new lambertian(new constant_texture(vec3(0.12, 0.45, 0.15)));
        material *light = new diffuse_light(new constant_texture(vec3(15, 15, 15)));

        list[i++] = new yz_rect(0, 555, 0, 555, 555, green);
        list[i++] = new yz_rect(0, 555, 0, 555, 0, red);
        list[i++] = new xz_rect(213, 343, 227, 332, 554, light);
        list[i++] = new xz_rect(0, 555, 0, 555, 0, white);
        list[i++] = new xy_rect(0, 555, 0, 555, 555, white);

        return new hittable_list(list,i);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And the view info:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    vec3 lookfrom(278, 278, -800);
    vec3 lookat(278,278,0);
    float dist_to_focus = 10.0;
    float aperture = 0.0;
    float vfov = 40.0;

    camera cam(lookfrom, lookat, vec3(0,1,0), vfov, float(nx)/float(ny),
        aperture, dist_to_focus, 0.0, 1.0);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
We get:

  ![Image 7-3](../images/img-2-7-03.jpg)

</div>

<div class='together'>
This is very noisy because the light is small. But why are the other walls missing? They are facing
the wrong way. We need outward facing normals. Let’s make a hittable that does nothing but hold
another hittable, but reverses the normals:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class flip_normals : public hittable {
        public:
            flip_normals(hittable *p) : ptr(p) {}

            virtual bool hit(
                const ray& r, float t_min, float t_max, hit_record& rec) const {

                if (ptr->hit(r, t_min, t_max, rec)) {
                    rec.normal = -rec.normal;
                    return true;
                }
                else
                    return false;
            }

            virtual bool bounding_box(float t0, float t1, aabb& box) const {
                return ptr->bounding_box(t0, t1, box);
            }

            hittable *ptr;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
This makes Cornell:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    hittable *cornell_box() {
        hittable **list = new hittable*[6];
        int i = 0;
        material *red = new lambertian(new constant_texture(vec3(0.65, 0.05, 0.05)));
        material *white = new lambertian(new constant_texture(vec3(0.73, 0.73, 0.73)));
        material *green = new lambertian(new constant_texture(vec3(0.12, 0.45, 0.15)));
        material *light = new diffuse_light(new constant_texture(vec3(15, 15, 15)));

        list[i++] = new flip_normals(new yz_rect(0, 555, 0, 555, 555, green));
        list[i++] = new yz_rect(0, 555, 0, 555, 0, red);
        list[i++] = new xz_rect(213, 343, 227, 332, 554, light);
        list[i++] = new flip_normals(new xz_rect(0, 555, 0, 555, 555, white));
        list[i++] = new xz_rect(0, 555, 0, 555, 0, white);
        list[i++] = new flip_normals(new xy_rect(0, 555, 0, 555, 555, white));

        return new hittable_list(list,i);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And voila:

  ![Image 7-4](../images/img-2-7-04.jpg)

</div>



Instances
====================================================================================================

The Cornell Box usually has two blocks in it. These are rotated relative to the walls. First, let’s
make an axis-aligned block primitive that holds 6 rectangles:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class box: public hittable {
        public:
            box() {}
            box(const vec3& p0, const vec3& p1, material *ptr);
            virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const {
                box =  aabb(pmin, pmax);
                return true;
            }
            vec3 pmin, pmax;
            hittable *list_ptr;
    };

    box::box(const vec3& p0, const vec3& p1, material *ptr) {
        pmin = p0;
        pmax = p1;
        hittable **list = new hittable*[6];
        list[0] = new xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p1.z(), ptr);
        list[1] = new flip_normals(
            new xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p0.z(), ptr));
        list[2] = new xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p1.y(), ptr);
        list[3] = new flip_normals(
            new xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p0.y(), ptr));
        list[4] = new yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p1.x(), ptr);
        list[5] = new flip_normals(
            new yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p0.x(), ptr));
        list_ptr = new hittable_list(list,6);
    }

    bool box::hit(const ray& r, float t0, float t1, hit_record& rec) const {
        return list_ptr->hit(r, t0, t1, rec);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

<div class='together'>
Now we can add two blocks (but not rotated)

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    list[i++] = new box(vec3(130, 0, 65), vec3(295, 165, 230), white);
    list[i++] = new box(vec3(265, 0, 295), vec3(430, 330, 460), white);
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

</div>

<div class='together'>
This gives:

  ![Image 8-1](../images/img-2-8-01.jpg)

</div>

<div class='together'>
Now that we have boxes, we need to rotate them a bit to have them match the _real_ Cornell box. In
ray tracing, this is usually done with an _instance_. An instance is a geometric primitive that has
been moved or rotated somehow. This is especially easy in ray tracing because we don’t move
anything; instead we move the rays in the opposite direction. For example, consider a _translation_
(often called a _move_). We could take the pink box at the origin and add 2 to all its x components,
or (as we almost always do in ray tracing) leave the box where it is, but in its hit routine
subtract 2 off the x-component of the ray origin.

  ![Figure 8-1](../images/fig-2-8-01.jpg)

</div>

<div class='together'>
Whether you think of this as a move or a change of coordinates is up to you. The code for this, to
move any underlying hittable is a _translate_ instance.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class translate : public hittable {
        public:
            translate(hittable *p, const vec3& displacement)
                : ptr(p), offset(displacement) {}
            virtual bool hit(
                const ray& r, float t_min, float t_max, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const;
            hittable *ptr;
            vec3 offset;
    };

    bool translate::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
        ray moved_r(r.origin() - offset, r.direction(), r.time());
        if (ptr->hit(moved_r, t_min, t_max, rec)) {
            rec.p += offset;
            return true;
        }
        else
            return false;
    }

    bool translate::bounding_box(float t0, float t1, aabb& box) const {
        if (ptr->bounding_box(t0, t1, box)) {
            box = aabb(box.min() + offset, box.max() + offset);
            return true;
        }
        else
            return false;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
Rotation isn’t quite as easy to understand or generate the formulas for. A common graphics tactic is
to apply all rotations about the x, y, and z axes. These rotations are in some sense axis-aligned.
First, let’s rotate by theta about the z-axis. That will be changing only x and y, and in ways that
don’t depend on z.

  ![Figure 8-2](../images/fig-2-8-02.jpg)

</div>

<div class='together'>
This involves some basic trigonometry that uses formulas that I will not cover here. That gives you
the correct impression it’s a little involved, but it is straightforward, and you can find it in any
graphics text and in many lecture notes. The result for rotating counter-clockwise about z is:

  $$ x' = cos(\theta) \cdot x - sin(\theta) \cdot y $$
  $$ y' = sin(\theta) \cdot x + cos(\theta) \cdot y $$
</div>

The great thing is that it works for any $\theta$ and doesn’t need any cases for quadrants or
anything like that. The inverse transform is the opposite geometric operation: rotate by $-\theta$.
Here, recall that $cos(\theta) = cos(-\theta)$ and $sin(-\theta) = -sin(\theta)$, so the formulas
are very simple.

<div class='together'>
Similarly, for rotating about y (as we want to do for the blocks in the box) the formulas are:

  $$ x' =  cos(\theta) \cdot x + sin(\theta) \cdot z $$
  $$ z' = -sin(\theta) \cdot x + cos(\theta) \cdot z $$

And about the x-axis:

  $$ y' = cos(\theta) \cdot y - sin(\theta) \cdot z $$
  $$ z' = sin(\theta) \cdot y + cos(\theta) \cdot z $$
</div>

Unlike the situation with translations, the surface normal vector also changes, so we need to
transform directions too if we get a hit. Fortunately for rotations, the same formulas apply. If you
add scales, things get more complicated. See the web page https://in1weekend.blogspot.com/ for links
to that.

<div class='together'>
For a y-rotation class we have:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class rotate_y : public hittable {
        public:
            rotate_y(hittable *p, float angle);
            virtual bool hit(
                const ray& r, float t_min, float t_max, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const {
                box = bbox; return hasbox;
            }
            hittable *ptr;
            float sin_theta;
            float cos_theta;
            bool hasbox;
            aabb bbox;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
With constructor:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    rotate_y::rotate_y(hittable *p, float angle) : ptr(p) {
        float radians = (M_PI / 180.) * angle;
        sin_theta = sin(radians);
        cos_theta = cos(radians);
        hasbox = ptr->bounding_box(0, 1, bbox);
        vec3 min(FLT_MAX, FLT_MAX, FLT_MAX);
        vec3 max(-FLT_MAX, -FLT_MAX, -FLT_MAX);
        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < 2; j++) {
                for (int k = 0; k < 2; k++) {
                    float x = i*bbox.max().x() + (1-i)*bbox.min().x();
                    float y = j*bbox.max().y() + (1-j)*bbox.min().y();
                    float z = k*bbox.max().z() + (1-k)*bbox.min().z();
                    float newx = cos_theta*x + sin_theta*z;
                    float newz = -sin_theta*x + cos_theta*z;
                    vec3 tester(newx, y, newz);
                    for ( int c = 0; c < 3; c++ )
                    {
                        if ( tester[c] > max[c] )
                            max[c] = tester[c];
                        if ( tester[c] < min[c] )
                            min[c] = tester[c];
                    }
                }
            }
        }
        bbox = aabb(min, max);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And the hit function:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool rotate_y::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
        vec3 origin = r.origin();
        vec3 direction = r.direction();
        origin[0] = cos_theta*r.origin()[0] - sin_theta*r.origin()[2];
        origin[2] =  sin_theta*r.origin()[0] + cos_theta*r.origin()[2];
        direction[0] = cos_theta*r.direction()[0] - sin_theta*r.direction()[2];
        direction[2] = sin_theta*r.direction()[0] + cos_theta*r.direction()[2];
        ray rotated_r(origin, direction, r.time());
        if (ptr->hit(rotated_r, t_min, t_max, rec)) {
            vec3 p = rec.p;
            vec3 normal = rec.normal;
            p[0] = cos_theta*rec.p[0] + sin_theta*rec.p[2];
            p[2] = -sin_theta*rec.p[0] + cos_theta*rec.p[2];
            normal[0] = cos_theta*rec.normal[0] + sin_theta*rec.normal[2];
            normal[2] = -sin_theta*rec.normal[0] + cos_theta*rec.normal[2];
            rec.p = p;
            rec.normal = normal;
            return true;
        }
        else
            return false;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And the changes to Cornell are:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    list[i++] = new translate(
                    new rotate_y(new box(vec3(0,0,0), vec3(165,165,165), white), -18),
                    vec3(130,0,65)
                );
    list[i++] = new translate(
                    new rotate_y(new box(vec3(0,0,0), vec3(165,330,165), white), 15),
                    vec3(265,0,295)
                );
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
Which yields:

  ![Image 8-2](../images/img-2-8-02.jpg)

</div>



Volumes
====================================================================================================

One thing it’s nice to add to a ray tracer is smoke/fog/mist. These are sometimes called _volumes_
or _participating media_. Another feature that is nice to add is subsurface scattering, which is
sort of like dense fog inside an object. This usually adds software architectural mayhem because
volumes are a different animal than surfaces. But a cute technique is to make a volume a random
surface. A bunch of smoke can be replaced with a surface that probabilistically might or might not
be there at every point in the volume. This will make more sense when you see the code.

First, let’s start with a volume of constant density. A ray going through there can either scatter
inside the volume, or it can make it all the way through like the middle ray in the figure. More
thin transparent volumes, like a light fog, are more likely to have rays like the middle one. How
far the ray has to travel through the volume also determines how likely it is for the ray to make it
through.

  ![Figure 9-1](../images/fig-2-9-01.jpg)

<div class='together'>
As the ray passes through the volume, it may scatter at any point. The denser the volume, the more
likely that is. The probability that the ray scatters in any small distance $\delta L$ is:

  $$ probability = C \cdot \delta L $$
</div>

<div class='together'>
where $C$ is proportional to the optical density of the volume. If you go through all the
differential equations, for a random number you get a distance where the scattering occurs. If that
distance is outside the volume, then there is no “hit”. For a constant volume we just need the
density $C$ and the boundary. I’ll use another hittable for the boundary. The resulting class is:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class constant_medium : public hittable {
        public:
            constant_medium(hittable *b, float d, texture *a) : boundary(b), density(d) {
                phase_function = new isotropic(a);
            }
            virtual bool hit(
                const ray& r, float t_min, float t_max, hit_record& rec) const;
            virtual bool bounding_box(float t0, float t1, aabb& box) const {
                return boundary->bounding_box(t0, t1, box);
            }
            hittable *boundary;
            float density;
            material *phase_function;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
The scattering function of isotropic picks a uniform random direction:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    class isotropic : public material {
        public:
            isotropic(texture *a) : albedo(a) {}
            virtual bool scatter(
                const ray& r_in,
                const hit_record& rec,
                vec3& attenuation,
                ray& scattered) const {

                scattered = ray(rec.p, random_in_unit_sphere());
                attenuation = albedo->value(rec.u, rec.v, rec.p);
                return true;
            }
            texture *albedo;
    };
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
And the hit function is:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    bool constant_medium::hit(const ray& r, float t_min, float t_max, hit_record& rec)
    const {

        // Print occasional samples when debugging. To enable, set enableDebug true.
        const enableDebug = false;
        bool debugging = enableDebug && random_double() < 0.00001;

        hit_record rec1, rec2;

        if (boundary->hit(r, -FLT_MAX, FLT_MAX, rec1)) {
            if (boundary->hit(r, rec1.t+0.0001, FLT_MAX, rec2)) {

                if (debugging) std::cerr << "\nt0 t1 " << rec1.t << " " << rec2.t << '\n';

                if (rec1.t < t_min)
                    rec1.t = t_min;

                if (rec2.t > t_max)
                    rec2.t = t_max;

                if (rec1.t >= rec2.t)
                    return false;

                if (rec1.t < 0)
                    rec1.t = 0;

                float distance_inside_boundary = (rec2.t - rec1.t)*r.direction().length();
                float hit_distance = -(1/density) * log(random_double());

                if (hit_distance < distance_inside_boundary) {

                    rec.t = rec1.t + hit_distance / r.direction().length();
                    rec.p = r.point_at_parameter(rec.t);

                    if (debugging) {
                        std::cerr << "hit_distance = " <<  hit_distance << '\n'
                                  << "rec.t = " <<  rec.t << '\n'
                                  << "rec.p = " <<  rec.p << '\n';
                    }

                    rec.normal = vec3(1,0,0);  // arbitrary
                    rec.mat_ptr = phase_function;
                    return true;
                }
            }
        }
        return false;
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

The reason we have to be so careful about the logic around the boundary is we need to make sure this
works for ray origins inside the volume. In clouds, things bounce around a lot so that is a common
case.

<div class='together'>
If we replace the two blocks with smoke and fog (dark and light particles) and make the light bigger
(and dimmer so it doesn’t blow out the scene) for faster convergence:

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    hittable *cornell_smoke() {
        hittable **list = new hittable*[8];
        int i = 0;
        material *red = new lambertian(new constant_texture(vec3(0.65, 0.05, 0.05)));
        material *white = new lambertian(new constant_texture(vec3(0.73, 0.73, 0.73)));
        material *green = new lambertian(new constant_texture(vec3(0.12, 0.45, 0.15)));
        material *light = new diffuse_light(new constant_texture(vec3(7, 7, 7)));

        list[i++] = new flip_normals(new yz_rect(0, 555, 0, 555, 555, green));
        list[i++] = new yz_rect(0, 555, 0, 555, 0, red);
        list[i++] = new xz_rect(113, 443, 127, 432, 554, light);
        list[i++] = new flip_normals(new xz_rect(0, 555, 0, 555, 555, white));
        list[i++] = new xz_rect(0, 555, 0, 555, 0, white);
        list[i++] = new flip_normals(new xy_rect(0, 555, 0, 555, 555, white));

        hittable *b1 = new translate(
            new rotate_y(new box(vec3(0, 0, 0), vec3(165, 165, 165), white), -18),
            vec3(130,0,65));
        hittable *b2 = new translate(
            new rotate_y(new box(vec3(0, 0, 0), vec3(165, 330, 165), white),  15),
            vec3(265,0,295));

        list[i++] = new constant_medium(
            b1, 0.01, new constant_texture(vec3(1.0, 1.0, 1.0)));
        list[i++] = new constant_medium(
            b2, 0.01, new constant_texture(vec3(0.0, 0.0, 0.0)));

        return new hittable_list(list,i);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
</div>

<div class='together'>
We get:

  ![Image 9-1](../images/img-2-9-01.jpg)

</div>



A Scene Testing All New Features
====================================================================================================

Let’s put it all together, with a big thin mist covering everything, and a blue subsurface
reflection sphere (we didn’t implement that explicitly, but a volume inside a dielectric is what a
subsurface material is). The biggest limitation left in the renderer is no shadow rays, but that is
why we get caustics and subsurface for free. It’s a double-edged design decision.

    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ C++
    hittable *final() {
        int nb = 20;
        hittable **list = new hittable*[30];
        hittable **boxlist = new hittable*[10000];
        hittable **boxlist2 = new hittable*[10000];
        material *white = new lambertian( new constant_texture(vec3(0.73, 0.73, 0.73)));
        material *ground = new lambertian( new constant_texture(vec3(0.48, 0.83, 0.53)));
        int b = 0;
        for (int i = 0; i < nb; i++) {
            for (int j = 0; j < nb; j++) {
                float w = 100;
                float x0 = -1000 + i*w;
                float z0 = -1000 + j*w;
                float y0 = 0;
                float x1 = x0 + w;
                float y1 = 100*(random_double()+0.01);
                float z1 = z0 + w;
                boxlist[b++] = new box(vec3(x0,y0,z0), vec3(x1,y1,z1), ground);
            }
        }
        int l = 0;
        list[l++] = new bvh_node(boxlist, b, 0, 1);
        material *light = new diffuse_light( new constant_texture(vec3(7, 7, 7)));
        list[l++] = new xz_rect(123, 423, 147, 412, 554, light);
        vec3 center(400, 400, 200);
        list[l++] = new moving_sphere(center, center+vec3(30, 0, 0),
            0, 1, 50, new lambertian(new constant_texture(vec3(0.7, 0.3, 0.1))));
        list[l++] = new sphere(vec3(260, 150, 45), 50, new dielectric(1.5));
        list[l++] = new sphere(vec3(0, 150, 145), 50,
            new metal(vec3(0.8, 0.8, 0.9), 10.0));
        hittable *boundary = new sphere(vec3(360, 150, 145), 70, new dielectric(1.5));
        list[l++] = boundary;
        list[l++] = new constant_medium(boundary, 0.2,
            new constant_texture(vec3(0.2, 0.4, 0.9)));
        boundary = new sphere(vec3(0, 0, 0), 5000, new dielectric(1.5));
        list[l++] = new constant_medium(boundary, 0.0001,
            new constant_texture(vec3(1.0, 1.0, 1.0)));
        int nx, ny, nn;
        unsigned char *tex_data = stbi_load("earthmap.jpg", &nx, &ny, &nn, 0);
        material *emat =  new lambertian(new image_texture(tex_data, nx, ny));
        list[l++] = new sphere(vec3(400, 200, 400), 100, emat);
        texture *pertext = new noise_texture(0.1);
        list[l++] =  new sphere(vec3(220, 280, 300), 80, new lambertian( pertext ));
        int ns = 1000;
        for (int j = 0; j < ns; j++) {
            boxlist2[j] = new sphere(
                vec3(165*random_double(), 165*random_double(), 165*random_double()),
                10, white);
        }
        list[l++] = new translate(new rotate_y(
            new bvh_node(boxlist2, ns, 0.0, 1.0), 15), vec3(-100,270,395));
        return new hittable_list(list,l);
    }
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

<div class='together'>
Running it with 10,000 rays per pixel yields:

  ![Image 10-1](../images/img-2-10-01.jpg)

</div>

Now go off and make a really cool image of your own! See https://in1weekend.blogspot.com/ for
pointers to further reading and features, and feel free to email questions, comments, and cool
images to me at ptrshrl@gmail.com.



<!-- Markdeep: https://casual-effects.com/markdeep/ -->
<link rel='stylesheet' href='../style/book.css'>
<style class="fallback">body{visibility:hidden;white-space:pre;font-family:monospace}</style>
<script src="markdeep.min.js"></script>
<script src="https://casual-effects.com/markdeep/latest/markdeep.min.js"></script>
<script>window.alreadyProcessedMarkdeep||(document.body.style.visibility="visible")</script>