Sunday, January 31, 2016

Ray Tracing: the Next Week

There's been a lot of interest and positive feedback on my mini-book on ray tracing.

This page is for the sequel Ray Tracing: the Next Week, available on Kindle.

This page also gives links and pointers for each chapter.   The features covered are those in this picture:

Chapter 1: Motion Blur

This method was introduced in 1984 by Rob Cook.

Chapter 2: A Bounding Volume Hierarchy (BVH)

The construction method in the book can be improved by using the surface area heuristic (SAH).     When evaluating potential partitions, the one that minimized the surface area of the sum of volumes of the sub-trees is almost always good.  Here is a SAH-based build that cuts on the longest axis.

Chapter 3: Solid Texture Mapping

Chapter 4: Perlin Noise 

A fantastic tool explaining how it works by Andrew Kensler

Chapter 5: Image Texture Mapping

Chapter 6: Rectangles and Lights

The program in the book implicitly samples lights so there are no shadow rays.   If you want to get more efficient direct lighting you can either send shadow rays, or importance sample by sending more rays toward the lights.

Chapter 7: Instances

A general instance usually stores transformation matrices.   Composite transforms can all be in one node.   When scales are allowed handling the surface normals must be done with care.

Chapter 8: Volumes

Here's a derivation of the ray-constant-medium intersection.

It is straightforward to add nonuniform densities by adding a more sophisticated intersection method.   This is covered in this blog post.   It's pretty common knowledge in the ray tracing community, but not really in the intersection.


  1. In chapter 4 when you say you scale the point in the perlin constructor and you show the result images ... by how much do you scale it?

  2. I think this is it, but I will have to run it again to be sure. I'll do that.

    list[0] = new sphere(vec3(0,-1000, 0), 1000, new lambertian( pertext ));
    list[1] = new sphere(vec3(0, 2, 0), 2, new lambertian( pertext ));

  3. For the final scene (Chapter 9), you did not provide the parameters for the camera. Would you mind sharing so that I can reproduce a similar image?


    1. I **think** it is this (I will do a run and check when I get back to my computer):
      vec3 lookfrom(478, 278, -600);
      vec3 lookat(278,278,0);
      float dist_to_focus = 10.0;
      float aperture = 0.0;
      float vfov = 40.0;

  4. At some point in the book you say: "This is very noisy because the light is small". I have noticed right away after switching to lights (chapter 6) and rendering the very first image, that everything is incredibly noisy (and way more than your screenshots). I am trying to understand the actual reason: is it because when rays get scattered (randomly), because the light is small, it requires a lot more of them to hit the light? But then why is it that the result is noisy? I am not sure I understand what is actually happening.


    1. To get some color around 0.5 (just for example), if the light has a color 100.0, then about half a percent need to hit the light to get the 0.5. If you send 1000 rays per pixel, that is about 5 rays. You can easily due to luck get 3, 4, 5, 6, 7, 8 rays. Sometimes, because there are a million pixels, you might get 0 or 20 rays. As you add more rays per pixel you get more stability in the fraction of outliers. I'll do a blog post with more details at sometime today

  5. Thank you for your previous answers. In chapter 8 (Volumes) the code use a "isotropic" material (phase_function = new isotropic(a)). It does not seem that you defined this previously. Can you explain a bit more?


    1. It's a terrible term actually, but standard. It means constant. Wouldn't constant be a better term :P

      Geometrically it means scattering in every direction is equally likely.

    2. But what does the isotropic class look like?

    3. Here it is:

    4. This is not what I am asking. The constant_medium class does this in the constructor phase_function = new isotropic(a);

      What is this isotropic class? It is not defined in the book nor in the link you just sent.

    5. Nor is what you asked :) Sorry poor reading comprehension. Here it is:

      class isotropic : public material {
      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;

  6. In chapter 3, I do not understand how to write the bounding_box method for the hitable_list class and I can't find anything about it in the book. Any help would be much appreciated.

    1. Here is my version-- I do it at construction time.

      bool hitable_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;
      box = temp_box;
      for (int i = 1; i < list_size; i++) {
      if(list[0]->bounding_box(t0, t1, temp_box)) {
      box = surrounding_box(box, temp_box);
      return false;
      return true;

    2. And the utility function

      aabb surrounding_box(aabb box0, aabb box1) {
      vec3 small( fmin(box0.min().x(), box1.min().x()),
      fmin(box0.min().y(), box1.min().y()),
      fmin(box0.min().z(), box1.min().z()));
      vec3 big ( fmax(box0.max().x(), box1.max().x()),
      fmax(box0.max().y(), box1.max().y()),
      fmax(box0.max().z(), box1.max().z()));
      return aabb(small,big);

    3. In the for loop above you have list[0]. Should that be list[i]? Also, thanks for the awesome series. A section/week on PBR would be nice :)

  7. In chapter 4 when you introduce Ken Perlin's trick to put random unit vectors on the lattice points, I was unable to reproduce the image shown in the book before turbulence and the "turb()" method is discussed.

    I spent a lot of time and tracked down what seems to be the problem in the "value(float, float, vec3)" method of the noise_texture. Using the code for return value used when introducing "scale" was the problem:

    return vec3(1, 1, 1) * noise.noise(scale * p);

    This returned an undecipherable image of noise.

    I looked ahead and noticed that the value calculation code changed when using turbulence. So I tried the simplest return calculation and was able to get the results that I was seeking:

    return vec3(1, 1, 1)*0.5*(1 + noise.noise(scale * p));

    This code produced the correct image as seen in the book before continuing onto reading about turbulence.

    Thank you so much for these book Mr.Shirley. :) I'm having a great time stretching my brain while learning about ray tracing.

    1. Thanks for pointing this out! I imagine it's related to the perlin_generate function's...

      p[i] = unit_vector(Vec3(-1 + 2 * drand48(), -1 + 2 * drand48(), -1 + 2 * drand48()));

      ...which looks like it throws some lovely negative numbers into the mix sometimes, and negative numbers which become colors is not likely to end well.

  8. After i implement BVH my render time increased more than twice. According to profiler bvh.hit() /aabb.hit()/ vec3[] operator functions creates significant overhead. This is intended or somthing is wrong?

  9. not intended but could be right-- that is most definitely not a great BVH implementation. There are two main places that need improvement in the BVH presented-- the build (it would work better with a mildly more complicated SAH), and the ray-box intersection. The ray-box intersection that is reasonably easy to drop in is by Amy Williams and the code is online. But view this BVH as a placeholder with log(N) behavior and a bad time constant.

  10. This comment has been removed by the author.

  11. Thanks Peter. Hit function change in aabb greatly helped.
    434 seconds with optimized multi box hit function from Amy Williams pdf.
    2015 seconds with out BVH.
    4950 seconds with the basic function from the book.
    (last scene from first book with 50 samples. compiled with vs2015 update 3 community edition).

    1. Awesome! I am guessing that improved BVH build would help a lot too, but maybe not as much.

    2. As a quick and dirty 'hack' I dropped in OpenMP support to see how the path tracer would work over multiple threads (6 in my hardcoded case). Right before the 'for' loop that iterates through the per pixel samples, I dropped in:

      #pragma omp parallel for

      Coupled with BVH and Amy Williams optimisation one of my scenes went from 32 hours down to 5.5 hours! Quite an improvement. However, the image quality was compromised in places due to (I believe) the fact that the random number generator is not happy running across multiple threads. I guess each thread would needs it's own generator.

    3. I implemented bucket based multithreading using std lib only, did not see any image degradation so far (I use drand48 implementation for windows). I also use a little optimized sah based bvh from link above.
      P.S numbers I provided above was based only on multithread implementation.
      P.P.S By the way in book stated that there is some links about scaling here. But i don't see any. May be I missed something.

    4. Probably I have dropped the balls and don't have them yet! But if you want massive parallelism running a ton of full images and then averaging the images themselves in a tree fashion (so for example for 1024 images first average pairs to get 512 images and repeat until you have one).

  12. Hi Peter. A quick question about the motion blur in Chapter 1... I notice that the motion-blurred spheres in the scene are not blurred when looking at their reflections in the larger metal sphere. Is this a difficult thing to address? Many thanks, really enjoying the introduction to ray/path tracing! :)

    1. Good eye! You are right. The shadow and reflection rays need to keep the time of the camera ray.

    2. Ah yes, that makes sense! Thanks for your time and reply.

  13. Hi Peter. I am a little stuck when trying to implement non uniform volume acording your blog post Here is my code for class and hit function
    class nonuniform_medium : public hitable {
    nonuniform_medium(hitable *b, texture *d, double ov, texture *a) : boundary(b), density(d), overall_density(ov) {
    phase_function = new isotropic(a);
    ~nonuniform_medium() {
    delete density;
    delete boundary;
    delete phase_function;
    virtual bool hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const;
    virtual bool bounding_box(double t0, double t1, aabb& box) const {
    return boundary->bounding_box(t0, t1, box);

    double overall_density;
    hitable *boundary;
    texture *density;
    material *phase_function;
    bool nonuniform_medium::hit(const ray& r, const double t_min, const double t_max, hit_record& rec) const {
    hit_record rec1, rec2;
    if (boundary->hit(r, -DBL_MAX, DBL_MAX, rec1)) {
    if (boundary->hit(r, rec1.t + 0.0001, DBL_MAX, rec2)) {
    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;
    double distance_inside_boundary = (rec2.t - rec1.t)*r.direction().length();
    double hit_distance = 0;
    while (1) {
    double x = drand48();
    hit_distance += -log(1 - x) / overall_density;
    double y = drand48();
    vec3 new_pos = r.origin() + hit_distance*r.direction();
    if ((density->value(0, 0, new_pos)) / overall_density > y)
    if (hit_distance < distance_inside_boundary) {
    rec.t = rec1.t + hit_distance / r.direction().length();
    rec.p = r.point_at_parameter(rec.t);
    rec.normal = vec3(1, 0, 0);
    rec.mat_ptr = phase_function;
    return true;
    return false;

    I now i miss something but i dont understand what exactly.

    1. Hey Dimitry. This LOOKS right to me. So that suggests a subtle bug. I would try using this code for a constant density and seeing what happens-- that might provide some debugging case where you know the right answer.

    2. This comment has been removed by the author.

    3. Thanks for help Peter. This code, when I use it with perlin noise as density function produce result visually similar to constant density. (I setup low frequency perlin noise function on solid cube before use it on volume to get good variations). I also used same noise as color variation function for same volume and all works as i expected. Only problem is in density variation.

    4. Looks like i make it work.
      Here is some images:
      I post my (not sure right) solution tommorow.

    5. Excellent! That looks plausible. One issue with media is it's hard to know what's right. A debugging case I like is to make the medium **almost** constant. It should then look constant.

    6. Here is some of my thoughts and explanations about non-uniform medium. As I promised in previous post.
      As I understand probability density function are same for non-uniform and uniform medium.
      The only problem was this statement,
      if((density->value(0, 0, new_pos)) / overall_density > drand48())
      we need to check our density function at new position against random number given by drand48()
      As far as I understand we get random number from drand48() and check if density function returns density higher than this random number. Higher the density at this point higher the probability that this check will return true, and we have a hit. According to this observation output range of density function and output range of drand48() need to be same. Output of random texture function that I use as density function lies between zero and one. Output of drand48() lies between zero and one too. If I scale density function by max density of my non-uniform medium (for example smax = 0.1) according to blog post
      if (s(newposition) / smax > Y) break
      I will lose 90% of probability during check against drand48(). So I come to this solution:
      density->value(0, 0, new_pos)) > drand48)
      Of course this will work only if density function and drand48() outputs are normalized against each other. And if I need to scale overall density I need to do it in basic probability formula
      -log(1 - x) / overall_density;
      P.S Looks like this work but I am not sure that there are no mistakes in my conclusions.

  14. This comment has been removed by the author.

  15. Hi Peter, thanks for your series of books. I have really enjoyed following the code examples and putting together the path tracer.

    I had an issue whilst rendering the final image. The noise_texture was not rendering properly.

    In the book the code, the value method of the noise_class had the following return statement: -

    return vec3(1,1,1)*0.5*(1 + sin(scale*p.z() + 10*noise.turb(p)));

    on GitHub the return code statement is: -

    return vec3(1, 1, 1) * 0.5 * (1 + sin(scale * p.x() + 5 * noise.turb(scale * p)));

    I have tried the GitHub version and it produces a rendering as in the Image of the final scene, the book version does not work

    I hope this helps



  16. Thanks for great booklets.
    I think there is an error in the bounding_box method of hitable_list in the second chapter of the second booklet "ray tracing: the next week."
    I am reading a kindle version from Amazon.
    In the method, there is a for-loop iterating all the boxes of the elements of the list.
    However the iterating index i is not being used within the loop.
    A correction for this error seems to be changing list[0] to list[i].
    Thank you again for the good booklets series.