• The Idea:
  • Trace rays backward from the view plane to the light source(s)
  • Find the intersection between the ray and objects in the scene
  • Each trace follows a single stochastic path
  • Multiple rays are used per pixel to account for propagation differences
• Existing Elements:
  • A Ray class, an empty (for semantics only) Material interface, a Shape3D interface, an Intersection class, a CompositeShape (in the sense of the composite pattern) class, some classes that implement the Shape3D interface, and a Light class

• What's Needed:
  • A non-empty PathTracingMaterial interface
  • Some classes that implement the RayTracingMaterial interface
  • A PathTracer
• Some Simplifying "Assumptions":
  • A single Light
  • Use stocahstic sampling for antialiasing
  • Ignore depth of field (i.e., focus)

cppexamples/graphics3d/pathtracing/PathTracingMaterial.h

#ifndef edu_jmu_cs_PathTracingMaterial_h
#define edu_jmu_cs_PathTracingMaterial_h

#include "../Intersection.h"
#include "../Material.h"
#include "../Ray.h"
#include <glm/glm.hpp>

/**
 * Requirements of a PathTracingMaterial.
 * 
 * @author Prof. David Bernstein, James Madison University
 */
class PathTracingMaterial: public Material {
 public:
  /**
   * Propagate the given Ray.
   *
   * @param ray           The Ray to propagate
   * @param intersection  The intersection of the Ray and the Shape3D
   * @param attenuation   The (outbound) attenuation
   * @param propagated    The (outbound) propagated Ray
   * @retun               true if propagation occurs; false otherwise
   */
  virtual bool propagate(const Ray& ray, const Intersection& intersection,
                         glm::vec3* attenuation, Ray* propagated) const = 0;
};

#endif
        

cppexamples/graphics3d/pathtracing/Matte.h

#ifndef edu_jmu_cs_Matte_h
#define edu_jmu_cs_Matte_h

#include "PathTracingMaterial.h"

/**
 * An encapsulation of a path tracing Material that is a perfect diffuser
 * and, hence, appears to be matte. It will appear the same from any angle.
 *
 * @author Prof. David Bernstein, James Madison University
 */
class Matte : public PathTracingMaterial {
 private:
  glm::vec3 albedo;

 public:
  /**
   * Explicit Value Constructor.
   * 
   * @param a   The albedo
   */
  explicit Matte(const glm::vec3& a);

  /**
   * Propogate the incoming Ray in a random direction.
   */
  virtual bool propagate(const Ray& ray, const Intersection& intersection,
                         glm::vec3* attenuation, Ray* propagated) const;
};

#endif




        

cppexamples/graphics3d/pathtracing/Matte.cpp

#include "Matte.h"
#include "../graphics3dUtilities.h"

Matte::Matte(const glm::vec3& a) : albedo(a) {
}

bool Matte::propagate(const Ray& ray, const Intersection& intersection,
                      glm::vec3* attenuation, Ray* propagated) const  {
  glm::vec3 target = intersection.p + intersection.n +
      randomPointInUnitDisk();

  *propagated  = Ray(intersection.p, target - intersection.p);
  *attenuation = albedo;

  return true;
}
        

cppexamples/graphics3d/graphics3dUtilities.cpp (Fragment: randomPointInUnitDisk)

glm::vec3 randomPointInUnitDisk() {
  //Generate a point at random and see if it is in the unit disk.
  // If not, try again.

  // Another approach is to use the spherical Fibonacci sequence
  // as discussed in Sanger et al.

  glm::vec3 p;
  do {
    p = 2.0f * (glm::vec3(static_cast<float>(drand48()),
                          static_cast<float>(drand48()), 0.0f)
                - glm::vec3(1.0f, 1.0f, 0.0f));
  } while (dot(p, p) >= 1.0f);

  return p;
}
        

cppexamples/graphics3d/pathtracing/Glossy.h

#ifndef edu_jmu_cs_Glossy_h
#define edu_jmu_cs_Glossy_h

#include "PathTracingMaterial.h"

/**
 * A Material that has specular reflections.
 *
 * @author Prof. David Bernstein, James Madison University
 */
class Glossy: public PathTracingMaterial {
 private:
  glm::vec3 albedo;
  float fuzz;

 public:
  /**
   * Explicit Value Constructor.
   * 
   * @param a   The albedo
   * @param f   The fuzz factir
   */
  Glossy(const glm::vec3& a, float f);

  /**
   * Propogate the incoming Ray based on the angle of incidence.
   */
  virtual bool propagate(const Ray& ray, const Intersection& intersection,
                         glm::vec3* attenuation, Ray* propagated) const;
};

#endif




        

cppexamples/graphics3d/pathtracing/Glossy.cpp

#include "Glossy.h"
#include "../graphics3dUtilities.h"

Glossy::Glossy(const glm::vec3& a, float f) : albedo(a) {
  if (f < 1.0f) fuzz = f;
  else          fuzz = 1.0f;
}

bool Glossy::propagate(const Ray& ray, const Intersection& intersection,
                    glm::vec3* attenuation, Ray* propagated) const  {
    glm::vec3 reflected = reflect(normalize(ray.d),
                                  intersection.n);

    *propagated = Ray(intersection.p,
                      reflected + (fuzz*randomPointInUnitDisk()));

    *attenuation = albedo;
    return (dot(propagated->d, intersection.n) > 0);
  }
        

cppexamples/graphics3d/graphics3dUtilities.cpp (Fragment: reflect)

/**
 * Note: The orientation of l is reversed from some presentations.
 */
glm::vec3 reflect(const glm::vec3& l, const glm::vec3& n) {
  return l - 2*dot(l, n)*n;
}
        

cppexamples/graphics3d/pathtracing/Transparent.h

#ifndef edu_jmu_cs_Transparent_h
#define edu_jmu_cs_Transparent_h

#include "PathTracingMaterial.h"

/**
 * A Material that is Transparent and, hence, both reflects and
 * transmits the ray (after refraction).
 *
 * @author Prof. David Bernstein, James Madison University
 */
class Transparent : public PathTracingMaterial {
 public:
  float refractionIndex;

  /**
   * Explicit Value Constructor.
   *
   * @param ri   The index of refraction
   */
  explicit Transparent(float ri);

  /**
   * Randomly either reflect the incoming Ray (based on the angle of
   * incidence) or trasnmit it (based on the index of refraction).
   */
  virtual bool propagate(const Ray& ray, const Intersection& intersection,
                         glm::vec3* attenuation, Ray* propagated) const;
};

#endif




        

cppexamples/graphics3d/pathtracing/Transparent.cpp

#include "Transparent.h"
#include "../graphics3dUtilities.h"

#include <iostream>

Transparent::Transparent(float ri) : refractionIndex(ri) {
}

bool Transparent::propagate(const Ray& ray, const Intersection& intersection,
                         glm::vec3* attenuation, Ray* propagated) const  {
  *attenuation = glm::vec3(1.0f, 1.0f, 1.0f);
  glm::vec3 reflected = reflect(ray.d, intersection.n);

  glm::vec3 refracted;
  float     kappaR;    // The proportion of the energy that is reflected
  if (refract(ray.d, intersection.n, refractionIndex, &refracted))
    fresnel(ray.d, intersection.n, refractionIndex, &kappaR);
  else
    kappaR = 1.0;

  kappaR = clamp(0.0f, 1.0f, kappaR);

  // In this case, treat kappaR as the probability that the Ray is reflected
  if (static_cast<float>(drand48()) < kappaR)
    *propagated = Ray(intersection.p, reflected);
  else
    *propagated = Ray(intersection.p, refracted);
  return true;
}
        

cppexamples/graphics3d/graphics3dUtilities.cpp (Fragment: refract)

bool refract(const glm::vec3 &I, const glm::vec3 &N, const float &ior,
             glm::vec3* refracted) {
  glm::vec3 i = normalize(I);
  glm::vec3 n = N;
  float etai = 1;  // Assuming the ray comes from air
  float etat = ior;

  float cosi = clamp(-1, 1, dot(i, n));
  if (cosi < 0) {  // The ray is from the air to the object
    cosi = -cosi;
  } else {         // The Ray is from the object to the air
    std::swap(etai, etat);
    n = -N;
  }

  float ratio = etai / etat;
  float discriminant = 1 - ratio*ratio * (1 - cosi * cosi);

  if (discriminant > 0.0f) {
    *refracted = ratio * (i - cosi*n) - n * sqrtf(discriminant);
    return true;
  } else {
    return false;
  }
}
        

cppexamples/graphics3d/graphics3dUtilities.cpp (Fragment: fresnel)

void fresnel(const glm::vec3 &I, const glm::vec3 &N,
             const float &ior, float* kR) {
  float cosi = clamp(-1, 1, dot(I, N));
  float etai = 1, etat = ior;
  if (cosi > 0) {
    std::swap(etai, etat);
  }

  // Compute sini using Snell's law
  float sint = etai / etat * sqrtf(std::max(0.f, 1 - cosi * cosi));


  if (sint >= 1) {  // Total internal reflection
    *kR = 1;
  } else {
    float cost = sqrtf(std::max(0.f, 1 - sint * sint));
    cosi = fabsf(cosi);
    float Rs = ((etat * cosi) - (etai * cost))
        / ((etat * cosi) + (etai * cost));
    float Rp = ((etai * cosi) - (etat * cost))
        / ((etai * cosi) + (etat * cost));

    *kR = (Rs * Rs + Rp * Rp) / 2;
    *kR = clamp(0.0f, 1.0f, *kR);
  }
  //  kT = 1 - kR because of the conservation of energy
}

glm::vec3 lightIntensityAt(const glm::vec3 &p, const glm::vec3 n,
                           const Tracer* tracer) {
  glm::vec3 lightAmt(0.0f);
  glm::vec3 oShadow = p;

  for (int i = 0; i < tracer->lightsSize; ++i) {
    glm::vec3 l = tracer->light[i]->location - p;

    float lightDistanceSquared = dot(l, l);

    l = normalize(l);
    float cosl = std::max(0.0f, dot(l, n));

    // Determine if the point is in the shadow of a relevant object
    Intersection shadowIntersection;
    bool inShadow = tracer->scene->intersectsWith(Ray(oShadow, l),
        0.001, MAXFLOAT, &shadowIntersection)
        && shadowIntersection.t * shadowIntersection.t < lightDistanceSquared;

    if (!inShadow)
      lightAmt += glm::vec3(1.0f) * tracer->light[i]->color * cosl;
  }
  return lightAmt;
}


glm::vec3 randomPointInUnitDisk() {
  //Generate a point at random and see if it is in the unit disk.
  // If not, try again.

  // Another approach is to use the spherical Fibonacci sequence
  // as discussed in Sanger et al.

  glm::vec3 p;
  do {
    p = 2.0f * (glm::vec3(static_cast<float>(drand48()),
                          static_cast<float>(drand48()), 0.0f)
                - glm::vec3(1.0f, 1.0f, 0.0f));
  } while (dot(p, p) >= 1.0f);

  return p;
}


/**
 * Note: The orientation of l is reversed from some presentations.
 */
glm::vec3 reflect(const glm::vec3& l, const glm::vec3& n) {
  return l - 2*dot(l, n)*n;
}


bool refract(const glm::vec3 &I, const glm::vec3 &N, const float &ior,
             glm::vec3* refracted) {
  glm::vec3 i = normalize(I);
  glm::vec3 n = N;
  float etai = 1;  // Assuming the ray comes from air
  float etat = ior;

  float cosi = clamp(-1, 1, dot(i, n));
  if (cosi < 0) {  // The ray is from the air to the object
    cosi = -cosi;
  } else {         // The Ray is from the object to the air
    std::swap(etai, etat);
    n = -N;
  }

  float ratio = etai / etat;
  float discriminant = 1 - ratio*ratio * (1 - cosi * cosi);

  if (discriminant > 0.0f) {
    *refracted = ratio * (i - cosi*n) - n * sqrtf(discriminant);
    return true;
  } else {
    return false;
  }
}


float schlick(float cosine, float etat) {
  // Assume the other material is air (with an IOR of 1)
  float r0 = (1.0f - etat) / (1.0f + etat);
  r0 = r0 * r0;
  return r0 + (1.0f - r0)*pow((1.0f - cosine), 5.0f);
}

        

cppexamples/graphics3d/pathtracing/PathTracer.h

#ifndef edu_jmu_cs_PathTracer_h
#define edu_jmu_cs_PathTracer_h

#include "../Tracer.h"

/**
 * An encapsulation of a simple PathTracer.
 *
 * @author Prof. David Bernstein, James Madison University
 */
class PathTracer: public Tracer {
 private:
  int ns;

 public:
  /**
   * Explicit Value Constructor.
   *
   * @param s    The scene
   * @param l    The lights
   * @param ls   The number of lights (only one is used)
   * @param ms   The maximum number of steps/propagations
   * @param samples  The number of samples
   */
  PathTracer(Shape3D* s, Light** l, int ls, int ms, int samples);

  virtual glm::vec3 trace(const Ray& r, int steps) const;

  virtual void render(Camera cam, FrameBuffer *fb, int width, int height) const;
};

#endif
        

cppexamples/graphics3d/pathtracing/PathTracer.cpp

#include "../graphics3dUtilities.h"
#include "PathTracer.h"
#include "PathTracingMaterial.h"

#include <iostream>


PathTracer::PathTracer(Shape3D* s, Light** l, int ls, int ms, int samples) {
  scene = s;
  light = l;
  lightsSize = ls;
  maxSteps = ms;
  ns = samples;
}

glm::vec3 PathTracer::trace(const Ray& r, int steps) const {
  Intersection intersection;

  if (scene->intersectsWith(r, 0.001, MAXFLOAT, &intersection)) {
    Ray propagated;
    glm::vec3 attenuation;

    PathTracingMaterial* material =
        ((PathTracingMaterial*)(intersection.material));
        
    if (steps < maxSteps && material->propagate(r, intersection,
                                                &attenuation, &propagated)) {
      // Recurse
      return attenuation * trace(propagated, steps + 1);
    } else {
      // Done with the recursion so return black
      return glm::vec3(0.0f, 0.0f, 0.0f);
    }
  } else {
    // Return the color of the Ray from the Light
    glm::vec3 d = normalize(r.d);
    float t = 0.5*(d.y + 1.0);
    return (1.0f - t) * glm::vec3(1.0f, 1.0f, 1.0f) + t * light[0]->color;
  }
}


void PathTracer::render(Camera cam, FrameBuffer *fb,
                        int width, int height) const {
  for (int j = 0; j < height; j++) {
    for (int i = 0; i < width; i++) {
      glm::vec3 col(0.0f, 0.0f, 0.0f);
      for (int sample=0; sample < ns; sample++) {
        float s = static_cast<float>(i + static_cast<float>(drand48()))
                                     / static_cast<float>(width);
        float t = static_cast<float>(j + static_cast<float>(drand48()))
                                     / static_cast<float>(height);

        Ray r = cam.getRay(s, t);

        col += trace(r, 0);
      }

      col /= static_cast<float>(ns);
      col = glm::vec3(sqrt(col[0]), sqrt(col[1]), sqrt(col[2]));
      int ir = clamp(0, 255, static_cast<int>(255.99*col[0]));
      int ig = clamp(0, 255, static_cast<int>(255.99*col[1]));
      int ib = clamp(0, 255, static_cast<int>(255.99*col[2]));

      fb->setPixel(i, j, ir, ig, ib);
    }
    std::cout << "Row " << j << " of " << height << "\n";
  }
}