Photon mapping can be executed now. Buggy as heck.

kd_tree.cpp

@@ -1,6 +1,16 @@
-#include "kd_tree.hpp"
+#ifndef NDEBUG
+#include <iostream>
+#include <fstream>
+#endif
 #include <queue>
 
+#include "kd_tree.hpp"
+
+#ifndef NDEBUG
+using std::ofstream;
+using std::ios;
+using std::endl;
+#endif
 
 treeNode::treeNode(Photon p, superKey plane)
 {
@@ -125,6 +135,7 @@ kdTree::~kdTree(){}
 
 void kdTree::addPhoton(Photon p)
 {
+  
   Photons.push_back(p);
 }
 
@@ -259,6 +270,16 @@ bool kdTree::buildKdTree()
 
   //printTree();
 
+#ifndef NDEBUG
+  ofstream ofs("photons.txt", ios::out);
+  float r, g, b;
+  for (std::vector<Photon>::iterator it = Photons.begin(); it != Photons.end(); it++) {
+    rgbe2float(r, g, b, (*it).radiance);
+    ofs << (*it).position.x << " " << (*it).position.y << " " << (*it).position.z << " " << r << " " << g << " " << b << endl;
+  }
+  ofs.close();
+#endif
+  
   return true;
 }
 
@@ -287,7 +308,7 @@ void kdTree::printNode(treeNode* node)
 {
 }
 
-std::vector<Photon> kdTree::findInRange (Vec3 min, Vec3 max)
+std::vector<Photon> kdTree::findInRange (Vec3 min, Vec3 max) const
 {
   std::vector<Photon> photons;
 
@@ -296,7 +317,7 @@ std::vector<Photon> kdTree::findInRange (Vec3 min, Vec3 max)
   return photons;
 }
 
-void kdTree::findInRange (Vec3 min, Vec3 max, std::vector<Photon> &photons, treeNode *node)
+void kdTree::findInRange (Vec3 min, Vec3 max, std::vector<Photon> &photons, treeNode *node) const
 {	
   if(node == NULL) return;
 

kd_tree.hpp

@@ -148,7 +148,7 @@ public:
   void addPhoton(Photon p);
   bool buildKdTree();
   void printTree();
-  std::vector<Photon> findInRange (Vec3 min, Vec3 max);
+  std::vector<Photon> findInRange (Vec3 min, Vec3 max) const;
 
 private:
   treeNode* root;
@@ -179,5 +179,5 @@ private:
 
   void printNode(treeNode* node);
 
-  void findInRange (Vec3 min, Vec3 max, std::vector<Photon> &photons, treeNode *node);
+  void findInRange (Vec3 min, Vec3 max, std::vector<Photon> &photons, treeNode *node) const;
 };

main.cpp

@@ -17,6 +17,7 @@
 #include "tracer.hpp"
 #include "path_tracer.hpp"
 #include "whitted_tracer.hpp"
+#include "photon_tracer.hpp"
 
 using namespace std;
 using namespace glm;
@@ -55,6 +56,8 @@ static tracer_t g_tracer = NONE;
 static unsigned int g_max_depth = 5;
 static float g_gamma = 2.2f;
 static float g_exposure = 0.0f;
+static size_t g_photons = 15000;
+static float g_p_sample_radius = 0.01f;
 
 ////////////////////////////////////////////
 // Main function.
@@ -63,6 +66,7 @@ int main(int argc, char ** argv) {
   Ray r;
   vec2 sample;
   Tracer * tracer;
+  PhotonTracer * p_tracer;
   size_t total;
   size_t current = 0;
   FIBITMAP * input_bitmap;
@@ -108,9 +112,12 @@ int main(int argc, char ** argv) {
     tracer = static_cast<Tracer *>(new PathTracer(g_max_depth));
 
   } else if(g_tracer == JENSEN) {
-    cerr << "Photon mapping coming soon." << endl;
+    cout << "Using " << ANSI_BOLD_YELLOW << "Jensen's photon mapping" << ANSI_RESET_STYLE << " with ray tracing." << endl;
-    return EXIT_FAILURE;
+    p_tracer = new PhotonTracer(g_max_depth, g_p_sample_radius);
-
+    cout << "Building photon map with " << ANSI_BOLD_YELLOW << g_photons << ANSI_RESET_STYLE << " primary photons per light source." << endl;
+    p_tracer->build_photon_map(scn, g_photons, false);
+    tracer = static_cast<Tracer *>(p_tracer);
+    
   } else {
     cerr << "Must specify a ray tracer with \"-t\"." << endl;
     print_usage(argv);
@@ -208,9 +215,9 @@ void print_usage(char ** const argv) {
   cerr << "Renders the scene specified by the scene file FILE." << endl << endl;
   cerr << "Mandatory options: " << endl;
   cerr << "  -t\tRay tracing method to use. Valid values: " << endl;
-  cerr << "    \twhitted     Classic Whitted ray tracing." << endl;
+  cerr << "    \t" << ANSI_BOLD_YELLOW << "whitted" << ANSI_RESET_STYLE << "     Classic Whitted ray tracing." << endl;
-  cerr << "    \tmonte_carlo Monte Carlo path tracing." << endl;
+  cerr << "    \t" << ANSI_BOLD_YELLOW << "monte_carlo" << ANSI_RESET_STYLE << " Monte Carlo path tracing." << endl;
-  cerr << "    \tjensen      Photon mapping. " << endl << endl;
+  cerr << "    \t" << ANSI_BOLD_YELLOW << "jensen" << ANSI_RESET_STYLE << "      Photon mapping. " << endl << endl;
   cerr << "Extra options:" << endl;
   cerr << "  -o\tOutput image file name with extension." << endl;
   cerr << "    \tDefaults to \"output.png\"." << endl;
@@ -226,11 +233,16 @@ void print_usage(char ** const argv) {
   cerr << "    \tDefaults to 2.2" << endl;
   cerr << "  -e\tExposure scale factor (in [-8, 8])." << endl;
   cerr << "    \tDefaults to 0.0 (no correction)." << endl;
+  cerr << "  -p\tNumber of primary photons per light source." << endl;
+  cerr << "    \tDefaults to 15000." << endl;
+  cerr << "  -h\tHemisphere radius for photon map sampling (> 0)." << endl;
+  cerr << "    \tDefaults to 0.01f ." << endl;
 }
 
 void parse_args(int argc, char ** const argv) {
   int opt;
   int x_pos;
+  int photons;
 
   // Check command line arguments.
   if(argc == 1) {
@@ -238,7 +250,7 @@ void parse_args(int argc, char ** const argv) {
     exit(EXIT_FAILURE);
   }
 
-  while((opt = getopt(argc, argv, "-:t:s:w:f:o:r:g:e:")) != -1) {
+  while((opt = getopt(argc, argv, "-:t:s:w:f:o:r:g:e:p:h:")) != -1) {
     switch (opt) {
     case 1:
       g_input_file = (char *)malloc((strlen(optarg) + 1) * sizeof(char));
@@ -328,7 +340,28 @@ void parse_args(int argc, char ** const argv) {
       }
 
       break;
-      
+
+    case 'p':
+      photons = atoi(optarg);
+      if (photons <= 0) {
+	cerr << "The number of photons must be a positive integer." << endl;
+	print_usage(argv);
+	exit(EXIT_FAILURE);
+      }
+      g_photons = (size_t)photons;
+
+      break;
+
+    case 'h':
+      g_p_sample_radius = atof(optarg);
+      if (g_p_sample_radius <= 0.0f) {
+	cerr << "Photon map sampling radius must be greater than 0.0" << endl;
+	print_usage(argv);
+	exit(EXIT_FAILURE);
+      }
+
+      break;
+
     case ':':
       cerr << "Option \"-" << static_cast<char>(optopt) << "\" requires an argument." << endl;
       print_usage(argv);

photon_tracer.cpp

@@ -1,24 +1,39 @@
+#include <iostream>
+#include <iomanip>
 #include <limits>
+#include <vector>
 
 #include <glm/gtc/constants.hpp>
 
 #include "photon_tracer.hpp"
 #include "sampling.hpp"
 #include "area_light.hpp"
+#include "rgbe.hpp"
 
+using std::cout;
+using std::endl;
+using std::setw;
+using std::vector;
 using std::numeric_limits;
 using namespace glm;
 
+#define ANSI_BOLD_YELLOW "\x1b[1;33m"
+#define ANSI_RESET_STYLE "\x1b[m"
+
 PhotonTracer::~PhotonTracer() { }
 
 vec3 PhotonTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
   float t, _t;
   Figure * _f;
-  vec3 n, color, i_pos, ref, dir_diff_color, dir_spec_color;
+  vec3 n, color, i_pos, ref, dir_diff_color, dir_spec_color, p_contrib;
   Ray mv_r, sr, rr;
   bool vis, is_area_light;
   float kr;
   AreaLight * al;
+  Vec3 mn;
+  Vec3 mx;
+  vector<Photon> photons;
+  float red, green, blue;
 
   t = numeric_limits<float>::max();
   _f = NULL;
@@ -97,6 +112,17 @@ vec3 PhotonTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
       } else if (_f->m_mat->m_rho > 0.0f && rec_level >= m_max_depth)
 	  return vec3(0.0f);
 
+      // TODO: Change photon map search method for hemisphere search.
+      mn = Vec3(i_pos.x - m_h_radius, i_pos.y - m_h_radius, i_pos.z - m_h_radius);
+      mx = Vec3(i_pos.x + m_h_radius, i_pos.y + m_h_radius, i_pos.z + m_h_radius);
+      photons =  m_photon_map.findInRange(mn, mx);
+      for (vector<Photon>::iterator it = photons.begin(); it != photons.end(); it++) {
+	rgbe2float(red, green, blue, (*it).radiance);
+	p_contrib += vec3(red, green, blue);
+      }
+      p_contrib /= pi<float>() * (m_h_radius * m_h_radius);
+      color += p_contrib;
+      
     } else {
       // If the material has transmission enabled, calculate the Fresnel term.
       kr = fresnel(r.m_direction, n, r.m_ref_index, _f->m_mat->m_ref_index);
@@ -130,17 +156,26 @@ void PhotonTracer::build_photon_map(Scene * s, const size_t n_photons_per_ligth,
   vec3 l_sample, s_normal, h_sample;
   float r1, r2;
   Ray rr;
+  size_t total = 0, current = 0;
 
   for (vector<Light *>::iterator it = s->m_lights.begin(); it != s->m_lights.end(); it++) {
+    total += (*it)->light_type() == Light::AREA ? 1 : 0;
+  }
+  total *= n_photons_per_ligth;
+
+  cout << "Tracing a total of " << ANSI_BOLD_YELLOW << total << ANSI_RESET_STYLE << " primary photons:" << endl;
+  
+  for (vector<Light *>::iterator it = s->m_lights.begin(); it != s->m_lights.end(); it++) {
+    l = *it;
+
+    /* Only area lights supported right now. */
+    if (l->light_type() != Light::AREA)
+      continue;
+
+    al = static_cast<AreaLight *>(l);
+      
+#pragma omp parallel for schedule(dynamic, 1) private(l_sample, s_normal, h_sample, r1, r2, rr) shared(current)
     for (size_t p = 0; p < n_photons_per_ligth; p++) {
-      l = *it;
-
-      /* Only area lights supported right now. */
-      if (l->light_type() != Light::AREA)
-	continue;
-
-      al = static_cast<AreaLight *>(l);
-
       if (!specular) {
 	l_sample = al->sample_at_surface();
 	s_normal = al->normal_at_last_sample();
@@ -156,8 +191,17 @@ void PhotonTracer::build_photon_map(Scene * s, const size_t n_photons_per_ligth,
       }
       
       trace_photon(rr, s, 0, specular);
+
+#pragma omp atomic
+	current++;
     }
+
+    cout << "\r" << setw(3) << static_cast<size_t>((static_cast<double>(current) / static_cast<double>(total)) * 100.0) << "% done.";
   }
+  cout << endl;
+
+  cout << "Building photon map Kd-tree." << endl;
+  m_photon_map.buildKdTree();
 }
 
 vec3 PhotonTracer::trace_photon(Ray &r, Scene * s, const unsigned int rec_level, const bool specular) {
@@ -199,7 +243,11 @@ vec3 PhotonTracer::trace_photon(Ray &r, Scene * s, const unsigned int rec_level,
     if (is_area_light) {
       p_pos = Vec3(i_pos.x, i_pos.y, i_pos.z);
       photon = Photon(p_pos, _f->m_mat->m_emission.r, _f->m_mat->m_emission.g, _f->m_mat->m_emission.b);
-      m_photon_map.addPhoton(photon);
+
+#pragma omp critical
+      {
+	m_photon_map.addPhoton(photon);
+      }
       
       return _f->m_mat->m_emission;
 
@@ -312,7 +360,10 @@ vec3 PhotonTracer::trace_photon(Ray &r, Scene * s, const unsigned int rec_level,
     
     p_pos = Vec3(i_pos.x, i_pos.y, i_pos.z);
     photon = Photon(p_pos, color.r, color.g, color.b);
-    m_photon_map.addPhoton(photon);
+#pragma omp critical
+    {
+      m_photon_map.addPhoton(photon);
+    }
 
     // Return final color.
     return color;

photon_tracer.hpp

@@ -7,8 +7,8 @@
 
 class PhotonTracer: public Tracer {
 public:
-  PhotonTracer(): Tracer() { }
+  PhotonTracer(): Tracer(), m_h_radius(0.5f) { }
-  PhotonTracer(unsigned int max_depth): Tracer(max_depth) { };
+  PhotonTracer(unsigned int max_depth, float _r = 0.5f): Tracer(max_depth), m_h_radius(_r) { };
 
   virtual ~PhotonTracer();
   virtual vec3 trace_ray(Ray & r, Scene * s, unsigned int rec_level) const;
@@ -16,6 +16,7 @@ public:
   void build_photon_map(Scene * s, const size_t n_photons_per_ligth = 10000, const bool specular = false);
 
 private:
+  float m_h_radius;
   kdTree m_photon_map;
   vec3 trace_photon(Ray &r, Scene * s, const unsigned int rec_level, const bool specular = false);
 };