Merge pull request #14 from vonhyou/global-illum

Global illumination by path tracing

src/Geometry.cc

@@ -40,6 +40,8 @@ Vector3f Sphere::normal(const Vector3f &p) const {
   return (p - center).normalized();
 }
 
+Vector3f Sphere::sample() const { return center + Vector3f(radius, 0, 0); }
+
 bool isInRectangle(const Vector3f &p, const Vector3f &a, const Vector3f &b,
                    const Vector3f &c, const Vector3f &d, const Vector3f &n) {
   float s1 = (b - a).cross(p - a).dot(n);
@@ -53,7 +55,7 @@ bool isInRectangle(const Vector3f &p, const Vector3f &a, const Vector3f &b,
 
 Optional<float> Rectangle::intersect(const Ray &r) const {
   float denom = normal_.dot(r.direction());
-  if (abs(denom) < 1e-6f)
+  if (std::fabs(denom) < 1e-6f)
     return Optional<float>::nullopt;
 
   float t = -normal_.dot(r.origin() - p1) / denom;
@@ -67,3 +69,5 @@ Optional<float> Rectangle::intersect(const Ray &r) const {
 }
 
 Vector3f Rectangle::normal(const Vector3f &p) const { return normal_; }
+
+Vector3f Rectangle::sample() const { return p1; }

src/Geometry.h

@@ -20,6 +20,7 @@ public:
   virtual ~Geometry() = default;
   virtual Optional<float> intersect(const Ray &) const = 0;
   virtual Vector3f normal(const Vector3f &) const = 0;
+  virtual Vector3f sample() const = 0;
 
 protected:
   Geometry(Type type, float ka, float kd, float ks, const Vector3f &ca,
@@ -54,6 +55,7 @@ public:
 
   Optional<float> intersect(const Ray &) const override;
   Vector3f normal(const Vector3f &) const override;
+  Vector3f sample() const override;
 
 private:
   float radius;
@@ -70,6 +72,7 @@ public:
 
   Optional<float> intersect(const Ray &) const override;
   Vector3f normal(const Vector3f &) const override;
+  Vector3f sample() const override;
 
 private:
   Vector3f p1, p2, p3, p4;

src/HitRecord.h

@@ -9,6 +9,8 @@ using Eigen::Vector3f;
 
 class HitRecord {
 public:
+  HitRecord()
+      : t(0), ray_(Ray()), normal_(Vector3f::Zero()), geometry_(nullptr) {}
   HitRecord(float t, const Ray &r, Geometry *g) : t(t), ray_(r), geometry_(g) {}
   bool operator<(const HitRecord &) const;
 

src/Light.cc

@@ -20,6 +20,12 @@ Vector3f Light::is() const { return is_; }
 
 bool Light::isUse() const { return use; }
 
+Vector3f PointLight::getCenter() const { return center; }
+
+bool lightOnSurface(Vector3f center, Geometry *g) {
+  return (g->sample() - center).dot(g->normal(center)) < 1e-5;
+}
+
 Vector3f PointLight::illumination(const HitRecord &hit,
                                   const vector<Geometry *> &geometries) const {
   Vector3f shadingPoint = hit.point();
@@ -28,9 +34,9 @@ Vector3f PointLight::illumination(const HitRecord &hit,
   Ray shadowRay(shadingPoint, rayDirection);
 
   for (auto g : geometries)
-    if (g != geometry && g->intersect(shadowRay).hasValue() &&
-        g->type() == Geometry::Type::SPHERE)
-      return Vector3f::Zero();
+    if (g != geometry && g->intersect(shadowRay).hasValue())
+      if (g->type() == Geometry::Type::SPHERE || !lightOnSurface(center, g))
+        return Vector3f::Zero();
 
   Vector3f ambient_ =
       geometry->ka() * geometry->ca().array() * Scene::current->ai().array();
@@ -46,16 +52,19 @@ Vector3f PointLight::illumination(const HitRecord &hit,
   return specular_ + ambient_ + diffuse_;
 }
 
+Vector3f AreaLight::getCenter() const {
+  return p1 + (p4 - p1) / 2 + (p2 - p1) / 2;
+}
+
 Vector3f AreaLight::illumination(const HitRecord &hit,
                                  const vector<Geometry *> &geometries) const {
   Vector3f u = p4 - p1;
   Vector3f v = p2 - p1;
 
-  Vector3f color = Vector3f::Zero();
-
   if (useCenter) {
-    color += PointLight(*this, p1 + (u + v) / 2).illumination(hit, geometries);
+    return PointLight(*this, getCenter()).illumination(hit, geometries);
   } else {
+    Vector3f color = Vector3f::Zero();
     for (int y = 0; y < gridSize; ++y)
       for (int x = 0; x < gridSize; ++x) {
         Vector3f contribution =
@@ -63,7 +72,6 @@ Vector3f AreaLight::illumination(const HitRecord &hit,
                 .illumination(hit, geometries);
         color += contribution;
       }
+    return color / gridSize / gridSize;
   }
-
-  return color / gridSize / gridSize;
 }

src/Light.h

@@ -18,6 +18,7 @@ public:
   virtual ~Light() = default;
   virtual Vector3f illumination(const HitRecord &,
                                 const vector<Geometry *> &) const = 0;
+  virtual Vector3f getCenter() const = 0;
 
 protected:
   Light(Type type, const Vector3f &id, const Vector3f &is)
@@ -36,6 +37,7 @@ public:
   void setGridSize(unsigned int);
   void setUseCenter(bool);
   void setIsUse(bool);
+  Type type() const;
   Vector3f id() const;
   Vector3f is() const;
   bool isUse() const;
@@ -52,6 +54,9 @@ public:
 
 private:
   Vector3f p1, p2, p3, p4;
+
+public:
+  Vector3f getCenter() const override;
 };
 
 class PointLight : public Light {
@@ -67,6 +72,9 @@ public:
 
 private:
   Vector3f center;
+
+public:
+  Vector3f getCenter() const override;
 };
 
 #endif // !LIGHT_H_

src/Optional.h

@@ -1,4 +1,5 @@
 #ifndef OPTIONAL_H_
+#define OPTIONAL_H_
 
 namespace utils {
 

src/Output.h

@@ -11,11 +11,9 @@ using std::vector;
 
 class Output {
 public:
-  Output(const Vector3f &bgc, string path, int w, int h)
-      : red(vector<float>(w * h + 1, bgc.x())),
-        green(vector<float>(w * h + 1, bgc.y())),
-        blue(vector<float>(w * h + 1, bgc.z())), path(path), width(w),
-        height(h) {}
+  Output(string path, int w, int h)
+      : red(vector<float>(w * h + 1)), green(vector<float>(w * h + 1)),
+        blue(vector<float>(w * h + 1)), path(path), width(w), height(h) {}
 
   void write();
 

src/Parser.cc

@@ -18,9 +18,8 @@ const Vector3f getVector3f(const nlohmann::json &j) {
 const VectorXi getRpp(const nlohmann::json &j) {
   VectorXi rpp(j.size());
 
-  for (int i = 0; i < j.size(); ++i) {
+  for (int i = 0; i < j.size(); ++i)
     rpp[i] = j[i].get<int>();
-  }
 
   return rpp;
 }
@@ -40,6 +39,8 @@ Scene *Parser::getScene(const nlohmann::json &j) {
   sc->setAntialiasing(j.value("antialiasing", false));
   sc->setTwoSideRender(j.value("twosiderender", false));
   sc->setGlobalIllum(j.value("globalillum", false));
+  sc->setMaxBounce(j.value("maxbounce", 3));
+  sc->setProbTerminate(j.value("probTerminate", 0.33f));
   if (j.contains("raysperpixel"))
     sc->setRaysPerPixel(getRpp(j["raysperpixel"]));
 

src/Progress.h

@@ -0,0 +1,29 @@
+#ifndef PROGRESS_H_
+#define PROGRESS_H_
+
+#include <iostream>
+
+namespace utils {
+
+#define BAR_WIDTH 70
+
+class Progress {
+public:
+  static void of(float p) {
+    std::cout << "[";
+    int pos = BAR_WIDTH * p;
+    for (int i = 0; i < BAR_WIDTH; ++i) {
+      if (i < pos)
+        std::cout << "=";
+      else if (i == pos)
+        std::cout << ">";
+      else
+        std::cout << " ";
+    }
+    std::cout << "] " << int(p * 100.0) << " %\r";
+    std::cout.flush();
+  }
+};
+} // namespace utils
+
+#endif // !PROGRESS_H_

src/Ray.h

@@ -7,6 +7,7 @@ using Eigen::Vector3f;
 
 class Ray {
 public:
+  Ray() : origin_(Vector3f::Zero()), direction_(Vector3f::Zero()) {}
   Ray(const Vector3f &o, const Vector3f &d) : origin_(o), direction_(d) {}
 
 private:

src/RayTracer.cc

@@ -1,16 +1,60 @@
 #include "RayTracer.h"
-#include "../external/simpleppm.h"
 #include "HitRecord.h"
+#include "Light.h"
+#include "Optional.h"
 #include "Output.h"
 #include "Parser.h"
+#include "Progress.h"
 #include "Ray.h"
 
 #include <Eigen/Core>
+#include <algorithm>
 #include <cmath>
+#include <iostream>
 #include <queue>
+#include <random>
 
+using Eigen::VectorXi;
 using std::priority_queue;
 
+// help function declarations
+Ray getRay(int, int);
+Ray getRay(int, int, int, int);
+void writeColor(int, const Vector3f &);
+utils::Optional<Vector3f> trace(Ray r);
+Vector3f clamp(const Vector3f &);
+
+namespace camera {
+int width, height, gridWidth, gridHeight, raysPerPixel;
+Vector3f pos, u, v, du, dv, vpUpperLeft, pxUpperLeft, gdu, gdv;
+
+void init();
+} // namespace camera
+
+/**
+ * Student solution starts here
+ *
+ * The main.cpp provided by instructor will invoke this function
+ */
+void RayTracer::run() {
+  parse();
+
+  for (auto scene : scenes) {
+    if (Scene::current != nullptr)
+      delete Scene::current;
+
+    Scene::current = scene;
+    render();
+
+    Output::current->write();
+  }
+}
+
+/**
+ * Parse the scene stored in the json file
+ *
+ * Example scene files are in `assets` folder
+ */
 void RayTracer::parse() {
   for (auto i = json["output"].begin(); i != json["output"].end(); ++i)
     scenes.push_back(Parser::getScene(*i));
@@ -22,71 +66,264 @@ void RayTracer::parse() {
     lights.push_back(Parser::getLight(*i));
 }
 
-Ray getRay(int x, int y, const Vector3f &camPos, const Vector3f &pxUpperLeft,
-           const Vector3f &du, const Vector3f &dv) {
-  return Ray(camPos, pxUpperLeft + x * du + y * dv - camPos);
+Vector3f gammaCorrection(Vector3f color, float gammaInv) {
+  return Vector3f(std::pow(color.x(), gammaInv), std::pow(color.y(), gammaInv),
+                  std::pow(color.z(), gammaInv));
 }
 
-void RayTracer::calculateColor(const HitRecord &hit, int i) {
+/**
+ * Render the current scene
+ *
+ * For direction illumination and anti-aliasing, render by phong model
+ * for global illumination, use path-tracing method.
+ *
+ * (*) Global illumination will not work with anti-aliasing by the requirement
+ */
+void RayTracer::render() {
+  camera::init();
+
+  using namespace camera;
+  Output::current = new Output(Scene::current->name(), width, height);
+
+  for (int y = 0; y < height; ++y) {
+    utils::Progress::of((y + 1.0f) / height);
+
+    for (int x = 0; x < width; ++x) {
+      Vector3f color = Scene::current->backgroundColor();
+
+      if (Scene::current->globalIllum()) {
+        int success = 0;
+        Vector3f accumulate = Vector3f::Zero();
+        for (int j = 0; j < gridHeight; ++j)
+          for (int i = 0; i < gridWidth; ++i) {
+
+            if (x != width / 2 || y != height / 2 || i || j)
+              ; // goto DEBUG_COLOR;
+            Ray ray = getRay(x, y, i, j);
+            for (int rayNum = 0; rayNum < raysPerPixel; ++rayNum) {
+              utils::Optional<Vector3f> result = trace(ray);
+              if (result.hasValue()) {
+                accumulate += result.value();
+                success++;
+              }
+            }
+
+            //  std::cout << accumulate.transpose() << " (" << success <<
+            //  std::endl;
+          }
+
+        if (success)
+          color = gammaCorrection(accumulate / success, 1.0f / 2.1f);
+      } else {
+        Ray ray = getRay(x, y);
+        Optional<HitRecord> hitRecord = getHitRecord(ray);
+
+        if (hitRecord.hasValue()) {
+          HitRecord hit = hitRecord.value();
+          color = calculateColor(hit, y * width + x);
+        }
+      }
+    DEBUG_COLOR:
+      writeColor(y * width + x, clamp(color));
+    }
+  }
+  std::cout << std::endl;
+}
+
+/**
+ * Calculate color using phong model
+ */
+Vector3f RayTracer::calculateColor(const HitRecord &hit, int i) const {
   Vector3f result(0, 0, 0);
   for (auto light : lights)
     result += light->isUse() ? light->illumination(hit, geometries)
                              : Vector3f::Zero();
 
-  result = result.cwiseMax(0.0f).cwiseMin(1.0f);
-  Output::current->r(i, result.x());
-  Output::current->g(i, result.y());
-  Output::current->b(i, result.z());
+  return result;
 }
 
-void RayTracer::render() {
-  int width = Scene::current->width();
-  int height = Scene::current->height();
-  Vector3f cameraPos = Scene::current->center();
+/**
+ * Find the nearest geometry to intersect
+ */
+Optional<HitRecord> RayTracer::getHitRecord(Ray r, const Geometry *self) const {
+  priority_queue<HitRecord> records;
+  for (auto g : geometries) {
+    Optional<float> t = g->intersect(r);
+    if (t.hasValue() && g != self)
+      records.push(HitRecord(t.value(), r, g));
+  }
+
+  if (!records.empty()) {
+    HitRecord result = records.top();
+    result.calcNormal();
+    return Optional<HitRecord>(result);
+  }
+
+  return Optional<HitRecord>::nullopt;
+}
+
+Optional<HitRecord> RayTracer::getHitRecord(Ray r) const {
+  return getHitRecord(r, nullptr);
+}
+
+Light *RayTracer::singleLightSource() const {
+  for (auto light : lights)
+    if (light->isUse())
+      return light;
+  return nullptr;
+}
+
+Ray getRay(int x, int y) {
+  using namespace camera;
+  return Ray(pos, pxUpperLeft + x * du + y * dv - pos);
+}
+
+Ray getRay(int x, int y, int i, int j) {
+  using namespace camera;
+  return Ray(pos, vpUpperLeft + x * du + i * gdu + y * dv + j * gdv - pos);
+}
+
+Vector3f clamp(const Vector3f &color) {
+  return color.cwiseMax(0.0f).cwiseMin(1.0f);
+}
+
+void writeColor(int i, const Vector3f &color) {
+  Output::current->r(i, color.x());
+  Output::current->g(i, color.y());
+  Output::current->b(i, color.z());
+}
+
+// This should generate a higher quality random number
+float getRandomNumber() {
+  static std::uniform_real_distribution<float> distribution(0.0, 1.0);
+  static std::mt19937 generator;
+  return distribution(generator);
+}
+
+// Generate a randon point on a unit hemisphere
+Vector3f getRandomDirection() {
+RETRY_RANDOM:
+  float x = getRandomNumber() * 2 - 1;
+  float y = getRandomNumber() * 2 - 1;
+  if (x * x + y * y > 1)
+    goto RETRY_RANDOM;
+
+  return Vector3f(x, y, std::sqrt(1 - x * x - y * y));
+}
+
+Vector3f getGlobalRandDirection(Vector3f normal) {
+  Vector3f tangent = normal.cross(Vector3f::UnitX());
+  if (tangent.norm() < 1e-6f)
+    tangent = normal.cross(Vector3f::UnitY());
+
+  tangent.normalize();
+  Vector3f binormal = normal.cross(tangent);
+  Eigen::Matrix3f local2World;
+  local2World.col(0) = tangent;
+  local2World.col(1) = binormal.normalized();
+  local2World.col(2) = normal.normalized();
+
+  return local2World * getRandomDirection();
+}
+
+// Check if a light source in on a surface (or really near)
+bool lightOnSurface(HitRecord hit, const Light *l) {
+  Vector3f center = l->getCenter();
+  Geometry *g = hit.geometry();
+  Geometry::Type type = g->type();
+  if (type == Geometry::Type::RECTANGLE)
+    return (g->sample() - center).dot(g->normal(center)) < 1e-5;
+
+  return false;
+}
+
+Vector3f RayTracer::trace(HitRecord hit, int bounce, float prob) const {
+RETRY_TRACING:
+  bool finish = !bounce || (getRandomNumber() < prob);
+  Vector3f point = hit.point();
+  Light *light = singleLightSource();
+  Vector3f direction;
+  Geometry *geometry = hit.geometry();
+
+  if (finish)
+    direction = light->getCenter() - point;
+  else
+    direction = getGlobalRandDirection(hit.normal());
+
+  direction.normalize();
+  Ray ray(point + hit.normal() * 1e-6, direction);
+
+  Optional<HitRecord> hitRecord = getHitRecord(ray, geometry);
+  Vector3f traceColor = Vector3f::Zero();
+  if (!finish && hitRecord.hasValue())
+    traceColor = trace(hitRecord.value(), bounce - 1, prob);
+  else if (!finish && !hitRecord.hasValue())
+    goto RETRY_TRACING;
+  else if (finish)
+    if (!hitRecord.hasValue() ||
+        (hitRecord.hasValue() && lightOnSurface(hitRecord.value(), light)))
+      traceColor = light->id();
+
+  return traceColor.array() * geometry->cd().array() *
+         std::max(0.0f, hit.normal().dot(direction));
+}
+
+utils::Optional<Vector3f> RayTracer::trace(Ray r) const {
+  Optional<HitRecord> hitRecord = getHitRecord(r);
+  if (hitRecord.hasValue()) {
+    Vector3f color = trace(hitRecord.value(), Scene::current->maxBounce(),
+                           Scene::current->probTerminate());
+
+    if (color != Vector3f::Zero())
+      return utils::Optional<Vector3f>(color);
+  }
+
+  return utils::Optional<Vector3f>::nullopt;
+}
+
+namespace camera {
+int getGridWidth(VectorXi data) {
+  return data.size() != 2 && data.size() != 3 ? 1 : data.x();
+}
+
+int getGridHeight(VectorXi data) {
+  return data.size() == 2 ? data.x() : (data.size() == 3 ? data.y() : 1);
+}
+
+int getRayNumber(VectorXi data) {
+  return data.size() == 2 ? data.y() : (data.size() == 3 ? data.z() : 1);
+}
+
+/**
+ * Initialize camera parameters
+ *
+ * Construct the surface for viewing the world
+ */
+void init() {
+  width = Scene::current->width();
+  height = Scene::current->height();
+  pos = Scene::current->center();
   Vector3f lookAt = Scene::current->lookAt();
   float vpHeight =
-      2 * tan(Scene ::current->fov() / 180 * M_PI / 2) * lookAt.norm();
+      2 * tan(Scene::current->fov() / 180 * M_PI / 2) * lookAt.norm();
   float vpWidth = vpHeight * width / height;
-  Vector3f vpU = Vector3f(vpWidth, 0, 0);
-  Vector3f vpV = Vector3f(0, -vpHeight, 0);
-  Vector3f du = vpU / width;
-  Vector3f dv = vpV / height;
+  u = Vector3f(vpWidth, 0, 0);
+  v = Vector3f(0, -vpHeight, 0);
+  du = u / width;
+  dv = v / height;
+  vpUpperLeft = pos + lookAt - u / 2.0 - v / 2.0;
+  pxUpperLeft = vpUpperLeft + (du + dv) / 2.0;
 
-  Vector3f vpUpperLeft = cameraPos + lookAt - vpU / 2.0 - vpV / 2.0;
-  Vector3f pxUpperLeft = vpUpperLeft + (du + dv) / 2.0;
+  VectorXi data = Scene::current->raysPerPixel();
+  gridWidth = getGridWidth(data);
+  gridHeight = getGridHeight(data);
+  raysPerPixel = getRayNumber(data);
 
-  Output::current = new Output(Scene::current->backgroundColor(),
-                               Scene::current->name(), width, height);
-
-  for (int y = 0; y < height; ++y)
-    for (int x = 0; x < width; ++x) {
-      Ray ray = getRay(x, y, cameraPos, pxUpperLeft, du, dv);
-      priority_queue<HitRecord> records;
-      for (auto g : geometries) {
-        Optional<float> t = g->intersect(ray);
-        if (t.hasValue())
-          records.push(HitRecord(t.value(), ray, g));
-      }
-
-      if (!records.empty()) {
-        HitRecord hit = records.top();
-        hit.calcNormal();
-        calculateColor(hit, y * width + x);
-      }
-    }
-}
-
-void RayTracer::output() {
-  for (auto output : outputs)
-    output->write();
-}
-
-void RayTracer::run() {
-  parse();
-
-  for (auto scene : scenes) {
-    Scene::current = scene;
-    render();
-    Output::current->write();
+  gdu = Vector3f::Zero();
+  gdv = Vector3f::Zero();
+  if (gridWidth > 1 || gridHeight > 1) {
+    gdu = du / gridWidth;
+    gdv = dv / gridHeight;
   }
 }
+} // namespace camera

src/RayTracer.h

@@ -3,12 +3,15 @@
 
 #include "../external/json.hpp"
 #include "Geometry.h"
+#include "HitRecord.h"
 #include "Light.h"
 #include "Output.h"
 #include "Scene.h"
 
 #include <vector>
 
+using std::vector;
+
 class RayTracer {
 public:
   RayTracer(const nlohmann::json &j) : json(j) {}
@@ -16,15 +19,19 @@ public:
 
 private:
   nlohmann::json json;
-  std::vector<Scene *> scenes;
-  std::vector<Light *> lights;
-  std::vector<Geometry *> geometries;
-  std::vector<Output *> outputs;
+  vector<Scene *> scenes;
+  vector<Light *> lights;
+  vector<Geometry *> geometries;
+  vector<Output *> outputs;
 
   void parse();
-  void calculateColor(const HitRecord &, int);
   void render();
-  void output();
+  Optional<HitRecord> getHitRecord(Ray, const Geometry *) const;
+  Optional<HitRecord> getHitRecord(Ray) const;
+  Vector3f calculateColor(const HitRecord &, int) const;
+  Light *singleLightSource() const;
+  Optional<Vector3f> trace(Ray) const;
+  Vector3f trace(HitRecord, int, float) const;
 };
 
 #endif // !RAY_TRACER_H_

src/Scene.cc

@@ -10,6 +10,14 @@ int Scene::height() { return height_; }
 
 float Scene::fov() { return fov_; }
 
+bool Scene::globalIllum() { return globalIllum_; }
+
+int Scene::maxBounce() { return maxBounce_; }
+
+float Scene::probTerminate() { return probTerminate_; }
+
+Eigen::VectorXi Scene::raysPerPixel() const { return raysPerPixel_; }
+
 Vector3f Scene::ai() const { return ai_; }
 
 Vector3f Scene::center() const { return center_; }
@@ -35,3 +43,9 @@ void Scene::setTwoSideRender(bool twoSideRender) {
 void Scene::setGlobalIllum(bool globalIllum) {
   this->globalIllum_ = globalIllum;
 }
+
+void Scene::setMaxBounce(int maxBounce) { this->maxBounce_ = maxBounce; }
+
+void Scene::setProbTerminate(float probTerminate) {
+  this->probTerminate_ = probTerminate;
+}

src/Scene.h

@@ -30,22 +30,31 @@ private:
   bool antialiasing_ = false;
   bool twoSideRender_ = false;
   bool globalIllum_ = false;
+  int maxBounce_ = 3;
+  float probTerminate_ = 0.33;
 
 public:
+  static Scene *current;
+
   string name() const;
   int width();
   int height();
   float fov();
+  bool globalIllum();
+  int maxBounce();
+  float probTerminate();
   Vector3f ai() const;
   Vector3f center() const;
   Vector3f up() const;
   Vector3f lookAt() const;
   Vector3f backgroundColor() const;
+  Eigen::VectorXi raysPerPixel() const;
   void setRaysPerPixel(const Eigen::VectorXi &);
   void setAntialiasing(bool);
   void setTwoSideRender(bool);
   void setGlobalIllum(bool);
-  static Scene *current;
+  void setMaxBounce(int);
+  void setProbTerminate(float);
 };
 
 #endif // !SCENE_H_