Added photon map implementation from the book.

.clang_complete

@@ -4,3 +4,4 @@
 -DGLM_FORCE_RADIANS
 -DUSE_CPP11_RANDOM
 -fopenmp
+-fno-builtin

Makefile

@@ -4,9 +4,10 @@ PVDIR = PhotonViewer
 OBJECTS = main.o sampling.o camera.o environment.o disk.o plane.o sphere.o \
           phong_brdf.o hsa_brdf.o directional_light.o point_light.o \
           spot_light.o sphere_area_light.o disk_area_light.o scene.o tracer.o \
-          path_tracer.o whitted_tracer.o rgbe.o kd_tree.o photon_tracer.o
+          path_tracer.o whitted_tracer.o rgbe.o kd_tree.o photon_tracer.o \
+          photonmap.o
 DEPENDS = $(OBJECTS:.o=.d)
-CXXFLAGS = -std=c++11 -pedantic -Wall -DGLM_FORCE_RADIANS -fopenmp -DUSE_CPP11_RANDOM #-DENABLE_KD_TREE
+CXXFLAGS = -std=c++11 -pedantic -Wall -DGLM_FORCE_RADIANS -fopenmp -DUSE_CPP11_RANDOM -fno-builtin #-DENABLE_KD_TREE
 LDLIBS = -lfreeimage -ljson_spirit
 
 .PHONY: all

main.cpp

@@ -64,6 +64,8 @@ static float g_exposure = 0.0f;
 static size_t g_photons = 15000;
 static float g_p_sample_radius = 0.01f;
 static float g_cone_filter_k = 1.0f;
+static int g_max_photons = 7000000;
+static int g_max_search  = 5000;
 
 ////////////////////////////////////////////
 // Main function.
@@ -118,7 +120,7 @@ int main(int argc, char ** argv) {
 
   case JENSEN:
     cout << "Using " << ANSI_BOLD_YELLOW << "Jensen's photon mapping" << ANSI_RESET_STYLE << " with ray tracing." << endl;
-    p_tracer = new PhotonTracer(g_max_depth, g_p_sample_radius, g_cone_filter_k);
+    p_tracer = new PhotonTracer(g_max_depth, g_p_sample_radius, g_cone_filter_k, g_max_photons, g_max_search);
     if (g_photons_file == NULL && g_caustics_file == NULL) {
       cout << "Building global photon map with " << ANSI_BOLD_YELLOW << g_photons / 2 << ANSI_RESET_STYLE << " primary photons per light source." << endl;
       p_tracer->photon_tracing(scn, g_photons / 2);
@@ -260,6 +262,10 @@ void print_usage(char ** const argv) {
   cerr << "    \tthe photons defined in the specified file." << endl;
   cerr << "  -l\tCone filter constant." << endl;
   cerr << "    \tDefaults to 1.0f." << endl;
+  cerr << "  -m\tMax number of photons in the photon map." << endl;
+  cerr << "    \tDefaults to 7000000." << endl;
+  cerr << "  -z\tMax number of photons for radiance estimate." << endl;
+  cerr << "    \tDefaults to 5000." << endl;
 }
 
 void parse_args(int argc, char ** const argv) {
@@ -273,7 +279,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:p:h:k:c:l:")) != -1) {
+  while((opt = getopt(argc, argv, "-:t:s:w:f:o:r:g:e:p:h:k:c:l:m:z:")) != -1) {
     switch (opt) {
     case 1:
       g_input_file = (char *)malloc((strlen(optarg) + 1) * sizeof(char));
@@ -404,6 +410,26 @@ void parse_args(int argc, char ** const argv) {
       }
       break;
 
+      case 'm':
+      g_max_photons = atoi(optarg);
+      if (g_max_photons <= 0) {
+	cerr << "Need to trace at least 1 photon." << endl;
+	print_usage(argv);
+	exit(EXIT_FAILURE);
+      }
+
+      break;
+
+      case 'z':
+      g_max_search = atoi(optarg);
+      if (g_max_search <= 0) {
+	cerr << "Need to search at least 1 photon." << endl;
+	print_usage(argv);
+	exit(EXIT_FAILURE);
+      }
+
+      break;
+      
     case ':':
       cerr << "Option \"-" << static_cast<char>(optopt) << "\" requires an argument." << endl;
       print_usage(argv);

path_tracer.cpp

@@ -122,8 +122,7 @@ vec3 PathTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
       amb_color = vis ? s->m_env->get_color(rr) * max(dot(n, rr.m_direction), 0.0f) / PDF : vec3(0.0f);
 
       // Add lighting.
-      color += (1.0f - _f->m_mat->m_rho) * (((dir_diff_color + ind_color + amb_color) * (_f->m_mat->m_diffuse / pi<float>())) +
+      color += ((dir_diff_color + ind_color + amb_color) * (_f->m_mat->m_diffuse / pi<float>())) + (_f->m_mat->m_specular * dir_spec_color);
-					    (_f->m_mat->m_specular * dir_spec_color));
 
       // Determine the specular reflection color.
       if (_f->m_mat->m_rho > 0.0f && rec_level < m_max_depth) {

photon_tracer.cpp

@@ -19,6 +19,7 @@ using std::cout;
 using std::cerr;
 using std::endl;
 using std::ifstream;
+using std::ofstream;
 using std::ios;
 using std::setw;
 using std::vector;
@@ -32,15 +33,16 @@ using namespace glm;
 PhotonTracer::~PhotonTracer() { }
 
 vec3 PhotonTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
-  float t, _t, red, green, blue, kr, radius, r1, r2;
+  const float radius = m_h_radius * m_h_radius;
+  float t, _t, /*red, green, blue,*/ kr, r1, r2;
   Figure * _f;
   vec3 n, color, i_pos, ref, dir_spec_color, p_contrib, c_contrib, sample, amb_color;
   Ray mv_r, sr, rr;
   bool vis, is_area_light;
   AreaLight * al;
   Vec3 mn, mx;
-  vector<Photon> photons;
+  vector<PhotonAux> photons;
-  vector<Photon> caustics;
+  vector<PhotonAux> caustics;
 
   t = numeric_limits<float>::max();
   _f = NULL;
@@ -110,7 +112,7 @@ vec3 PhotonTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
       }
 
       // Calculate photon map contribution
-      radius = m_h_radius;
+      /*      radius = m_h_radius;
 #ifdef ENABLE_KD_TREE
       Vec3 vmin(i_pos.x - m_h_radius, i_pos.y - m_h_radius, i_pos.z - m_h_radius);
       Vec3 vmax(i_pos.x + m_h_radius, i_pos.y + m_h_radius, i_pos.z + m_h_radius);
@@ -146,7 +148,7 @@ vec3 PhotonTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
 #else
 	m_caustics_map.find_by_distance(caustics, i_pos, n, m_h_radius, 1000);
 #endif
-      }
+}
 
       for (Photon p : photons) {
 	p.getColor(red, green, blue);
@@ -154,11 +156,19 @@ vec3 PhotonTracer::trace_ray(Ray & r, Scene * s, unsigned int rec_level) const {
       }
       p_contrib /= (1.0f - (2.0f / (3.0f * m_cone_filter_k))) * pi<float>() * (radius * radius);
 
-      for (Photon p : caustics) {
+      for (PhotonAux p : caustics) {
 	p.getColor(red, green, blue);
 	c_contrib += vec3(red, green, blue);
       }
-      c_contrib /= (1.0f - (2.0f / (3.0f * m_cone_filter_k))) * pi<float>() * (radius * radius);
+      c_contrib /= (1.0f - (2.0f / (3.0f * m_cone_filter_k))) * pi<float>() * (radius * radius); */
+
+      float irrad[3];
+      float pos[3] {i_pos.x, i_pos.y, i_pos.z};
+      float normal[3] {n.x, n.y, n.z};
+      
+      m_photon_map.irradiance_estimate(irrad, pos, normal, m_h_radius, m_max_s_photons);
+      c_contrib = vec3(irrad[0], irrad[1], irrad[2]);
+      c_contrib /= (1.0f - (2.0f / (3.0f * m_cone_filter_k))) * pi<float>() * (radius);
       
       // Calculate environment light contribution
       vis = true;
@@ -221,7 +231,7 @@ void PhotonTracer::photon_tracing(Scene * s, const size_t n_photons_per_ligth, c
   vec3 l_sample, s_normal, h_sample, power;
   Vec3 ls, dir;
   float r1, r2;
-  Photon ph;
+  PhotonAux ph;
   uint64_t total = 0, current = 0;
   vector<Figure *> spec_figures;
 
@@ -282,7 +292,7 @@ void PhotonTracer::photon_tracing(Scene * s, const size_t n_photons_per_ligth, c
 	}
 
 	// Create the primary photon.
-	power = (al->m_figure->m_mat->m_emission / static_cast<float>(n_photons_per_ligth));
+	power = (al->m_figure->m_mat->m_emission /* / static_cast<float>(n_photons_per_ligth) */);
 	
       } else if (pl != NULL) {
 	l_sample = glm::vec3(pl->m_position.x, pl->m_position.y, pl->m_position.z);
@@ -294,12 +304,12 @@ void PhotonTracer::photon_tracing(Scene * s, const size_t n_photons_per_ligth, c
 	  h_sample = normalize(spec_figures[p % spec_figures.size()]->sample_at_surface() - l_sample);
 	}
 
-	power = (pl->m_diffuse / static_cast<float>(n_photons_per_ligth));
+	power = (pl->m_diffuse /* / static_cast<float>(n_photons_per_ligth)*/ );
       }
 
       ls = Vec3(l_sample.x, l_sample.y, l_sample.z);
       dir = Vec3(h_sample.x, h_sample.y, h_sample.z);
-      ph = Photon(ls, dir, power.r, power.g, power.b, 1.0f);
+      ph = PhotonAux(ls, dir, power.r, power.g, power.b, 1.0f);
       
       trace_photon(ph, s, 0);
 
@@ -307,16 +317,33 @@ void PhotonTracer::photon_tracing(Scene * s, const size_t n_photons_per_ligth, c
 	current++;
     }
 
+    m_photon_map.scale_photon_power(1.0f / n_photons_per_ligth);
+
     cout << "\r" << setw(3) << static_cast<size_t>((static_cast<double>(current) / static_cast<double>(total)) * 100.0) << "% done.";
   }
   cout << endl;
 
-  cout << "Generated " << ANSI_BOLD_YELLOW << m_photon_map.getNumPhotons() << ANSI_RESET_STYLE << " total photons." << endl;
+  cout << "Generated " << ANSI_BOLD_YELLOW << m_photon_map.stored_photons << ANSI_RESET_STYLE << " total photons." << endl;
-  m_photon_map.save_photon_list(specular ? "caustics.txt" : "photons.txt");
+  //m_photon_map.save_photon_list(specular ? "caustics.txt" : "photons.txt");
+
+  string file_name = specular ? "caustics.txt" : "photons.txt";
+  
+  cout << "Writing photons to \x1b[1;33m" << file_name << "\x1b[m" << endl;
+  ofstream ofs(file_name, ios::out);
+  for (int i = 0; i < m_photon_map.stored_photons; i++) {
+    float r, g, b;
+    float dir[3];
+    rgbe2float(r, g, b, m_photon_map.photons[i].power);
+    m_photon_map.photon_dir(dir, &m_photon_map.photons[i]);
+    ofs << m_photon_map.photons[i].pos[0] << " " << m_photon_map.photons[i].pos[1] << " " << m_photon_map.photons[i].pos[2] << " " <<
+      dir[0] << " " << dir[1] << " " << dir[2] << " " <<
+      r << " " << g << " " << b << " " << m_photon_map.photons[i].ref_index << endl;
+  }
+  ofs.close();
 }
 
 void PhotonTracer::build_photon_map(const char * photons_file, const bool caustics) {
-  Photon ph;
+  PhotonAux ph;
   float x, y, z, dx, dy, dz, r, g, b, rc;
   ifstream ifs;
 
@@ -333,10 +360,16 @@ void PhotonTracer::build_photon_map(const char * photons_file, const bool causti
   cout << "Reading photon definitions from " << ANSI_BOLD_YELLOW << photons_file << ANSI_RESET_STYLE << "." << endl;
   while (!ifs.eof()) {
     ifs >> x >> y >> z >> dx >> dy >> dz >> r >> g >> b >> rc;
-    ph = Photon(Vec3(x, y, z), Vec3(dx, dy, dz), r, g, b, rc);
+    ph = PhotonAux(Vec3(x, y, z), Vec3(dx, dy, dz), r, g, b, rc);
-    m_photon_map.addPhoton(ph);
+    //m_photon_map.addPhoton(ph);
+
+    float power[3] {r, g, b};
+    float pos[3] {x, y, z};
+    float dir[3] {dx, dy, dz};
+    
+    m_photon_map.store(power, pos, dir, rc);
   }
-  cout << "Read " << ANSI_BOLD_YELLOW << m_photon_map.getNumPhotons() << ANSI_RESET_STYLE << " photons from the file." << endl;
+  cout << "Read " << ANSI_BOLD_YELLOW << m_photon_map.stored_photons << ANSI_RESET_STYLE << " photons from the file." << endl;
 
   ifs.close();
 
@@ -351,10 +384,12 @@ void PhotonTracer::build_photon_map(const bool caustics) {
   else
   m_caustics_map.buildKdTree();
 #endif
+
+  m_photon_map.balance();
 }
 
-void PhotonTracer::trace_photon(Photon & ph, Scene * s, const unsigned int rec_level) {
+void PhotonTracer::trace_photon(PhotonAux & ph, Scene * s, const unsigned int rec_level) {
-  Photon photon;
+  PhotonAux photon;
   float t, _t, red, green, blue;
   Figure * _f;
   vec3 n, color, i_pos, sample, ph_dir, ph_pos;
@@ -387,8 +422,12 @@ void PhotonTracer::trace_photon(Photon & ph, Scene * s, const unsigned int rec_l
       {
 	p_pos = Vec3(i_pos.x, i_pos.y, i_pos.z);
 	p_dir = Vec3(-ph.direction.x, -ph.direction.y, -ph.direction.z);
-	photon = Photon(p_pos, p_dir, red, green, blue, ph.ref_index);
+	photon = PhotonAux(p_pos, p_dir, red, green, blue, ph.ref_index);
-	m_photon_map.addPhoton(photon);
+	//m_photon_map.addPhoton(photon);
+	float power[3] {red, green, blue};
+	float pos[3] {p_pos.x, p_pos.y, p_pos.z};
+	float dir[3] {p_dir.x, p_dir.y, p_dir.z};
+	m_photon_map.store(power, pos, dir, ph.ref_index);
       }
 
       // Generate a photon for diffuse reflection.
@@ -400,7 +439,7 @@ void PhotonTracer::trace_photon(Photon & ph, Scene * s, const unsigned int rec_l
       color = (1.0f - _f->m_mat->m_rho) * (vec3(red, green, blue) * (_f->m_mat->m_diffuse / pi<float>()));
       p_pos = Vec3(i_pos.x, i_pos.y, i_pos.z);
       p_dir = Vec3(sample.x, sample.y, sample.z);
-      photon = Photon(p_pos, p_dir, color.r, color.g, color.b, ph.ref_index);
+      photon = PhotonAux(p_pos, p_dir, color.r, color.g, color.b, ph.ref_index);
 
       // Trace diffuse-reflected photon.
       if (rec_level < m_max_depth)
@@ -413,7 +452,7 @@ void PhotonTracer::trace_photon(Photon & ph, Scene * s, const unsigned int rec_l
       	p_pos = Vec3(i_pos.x, i_pos.y, i_pos.z);
       	ph_dir = normalize(reflect(vec3(ph.direction.x, ph.direction.y, ph.direction.z), n));
       	p_dir = Vec3(ph_dir.x, ph_dir.y, ph_dir.z);
-      	photon = Photon(p_pos, p_dir, color.r, color.g, color.b, ph.ref_index);
+      	photon = PhotonAux(p_pos, p_dir, color.r, color.g, color.b, ph.ref_index);
       	trace_photon(photon, s, rec_level + 1);
       }
 
@@ -429,7 +468,7 @@ void PhotonTracer::trace_photon(Photon & ph, Scene * s, const unsigned int rec_l
       	p_pos = Vec3(i_pos.x, i_pos.y, i_pos.z);
       	ph_dir = normalize(reflect(vec3(ph.direction.x, ph.direction.y, ph.direction.z), n));
       	p_dir = Vec3(ph_dir.x, ph_dir.y, ph_dir.z);
-      	photon = Photon(p_pos, p_dir, color.r, color.g, color.b, ph.ref_index);
+      	photon = PhotonAux(p_pos, p_dir, color.r, color.g, color.b, ph.ref_index);
       	trace_photon(photon, s, rec_level + 1);
       }
 
@@ -440,7 +479,7 @@ void PhotonTracer::trace_photon(Photon & ph, Scene * s, const unsigned int rec_l
       	p_pos = Vec3(i_pos.x, i_pos.y, i_pos.z);
       	ph_dir = normalize(refract(vec3(ph.direction.x, ph.direction.y, ph.direction.z), n, ph.ref_index / _f->m_mat->m_ref_index));
       	p_dir = Vec3(ph_dir.x, ph_dir.y, ph_dir.z);
-      	photon = Photon(p_pos, p_dir, color.r, color.g, color.b, _f->m_mat->m_ref_index);
+      	photon = PhotonAux(p_pos, p_dir, color.r, color.g, color.b, _f->m_mat->m_ref_index);
       	trace_photon(photon, s, rec_level + 1);
       }
     }

photon_tracer.hpp

@@ -3,12 +3,85 @@
 #define PHOTON_TRACER_HPP
 
 #include "tracer.hpp"
-#include "kd_tree.hpp"
+//#include "kd_tree.hpp"
+#include "photonmap.hpp"
+#include "rgbe.hpp"
+
+struct Vec3
+{
+  float x;
+  float y;
+  float z;
+
+  Vec3(float _x = 0.0f, float _y = 0.0f, float _z = 0.0f): x(_x), y(_y), z(_z) { }
+
+  Vec3(const Vec3 & other) = default;
+
+  glm::vec3 toVec3() {
+    return glm::vec3(x, y, z);
+  }
+  
+  inline bool equalFloat(const float x, const float y)
+  {
+    return (x - std::numeric_limits<float>::epsilon() <= y) && (x + std::numeric_limits<float>::epsilon() >= y);
+  }
+
+  inline bool operator<=(const Vec3 b)
+  {
+    return (x < b.x || (equalFloat(x, b.x) && (y < b.y || (equalFloat(y, b.y) && z <= b.z))));
+  }
+
+  inline bool operator>=(const Vec3 b)
+  {
+    return (x > b.x || (equalFloat(x, b.x) && (y > b.y || (equalFloat(y, b.y) && z >= b.z))));
+  }
+
+  inline friend std::ostream& operator<<(std::ostream& out, const Vec3& v)
+  {
+    return out << "X:" << v.x << " Y:" << v.y << " Z:" << v.z;
+  }
+};
+
+struct PhotonAux
+{
+  Vec3 position;
+  Vec3 direction;
+  float ref_index;
+  unsigned char radiance[4];
+  float r, g, b;
+
+  PhotonAux(Vec3 _p = Vec3(), Vec3 _d = Vec3(), float red = 0.0f, float green = 0.0f, float blue = 0.0f, float _r = 1.0f):
+    position(_p),
+    direction(_d),
+    ref_index(_r),
+    r(red),
+    g(green),
+    b(blue)
+  {
+    float2rgbe(radiance, red, green, blue);
+  }
+
+  inline void getColor(float & red, float & green, float & blue) {
+    rgbe2float(red, green, blue, radiance);
+  }
+};
 
 class PhotonTracer: public Tracer {  
 public:
-  PhotonTracer(): Tracer(), m_h_radius(0.5f), m_cone_filter_k(1.0f) { }
+  PhotonTracer():
-  PhotonTracer(unsigned int max_depth, float _r = 0.5f, float _k = 1.0f): Tracer(max_depth), m_h_radius(_r), m_cone_filter_k(_k < 1.0f ? 1.0f : _k) { };
+    Tracer(), m_h_radius(0.5f),
+    m_cone_filter_k(1.0f),
+    m_photon_map(7000000),
+    m_max_s_photons(5000)
+  { }
+  
+  PhotonTracer(unsigned int max_depth, float _r = 0.5f, float _k = 1.0f, const int max_photons = 7000000, const int max_search = 5000):
+    Tracer(max_depth),
+    m_h_radius(_r),
+    m_cone_filter_k(_k < 1.0f ? 1.0f : _k),
+    m_photon_map(max_photons),
+    m_max_s_photons(max_search)
+  { };
 
   virtual ~PhotonTracer();
   virtual vec3 trace_ray(Ray & r, Scene * s, unsigned int rec_level) const;
@@ -16,12 +89,15 @@ public:
   void photon_tracing(Scene * s, const size_t n_photons_per_ligth = 10000, const bool specular = false);
   void build_photon_map(const char * photons_file, const bool caustics = false);
   void build_photon_map(const bool caustics = false);
+
 private:
   float m_h_radius;
   float m_cone_filter_k;
-  kdTree m_photon_map;
+  /*kdTree m_photon_map;
-  kdTree m_caustics_map;
+    kdTree m_caustics_map;*/
-  void trace_photon(Photon & ph, Scene * s, const unsigned int rec_level);
+  PhotonMap m_photon_map;
+  int m_max_s_photons;
+  void trace_photon(PhotonAux & ph, Scene * s, const unsigned int rec_level);
 };
 
 #endif

photonmap.cpp

@@ -0,0 +1,464 @@
+//----------------------------------------------------------------------------
+// photonmap.cc
+// An example implementation of the photon map data structure
+//
+// Henrik Wann Jensen - February 2001
+//----------------------------------------------------------------------------
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <alloca.h>
+#include <string.h>
+#include <math.h>
+
+#include "photonmap.hpp"
+#include "rgbe.hpp"
+
+/* This is the constructor for the photon map.
+ * To create the photon map it is necessary to specify the
+ * maximum number of photons that will be stored
+*/
+//************************************************
+PhotonMap :: PhotonMap( const int max_phot )
+//************************************************
+{
+  stored_photons = 0;
+  prev_scale = 1;
+  max_photons = max_phot;
+
+  photons = (Photon*)malloc( sizeof( Photon ) * ( max_photons+1 ) );
+
+  if (photons == NULL) {
+    fprintf(stderr,"Out of memory initializing photon map\n");
+    exit(-1);
+  }
+
+  bbox_min[0] = bbox_min[1] = bbox_min[2] = 1e8f;
+  bbox_max[0] = bbox_max[1] = bbox_max[2] = -1e8f;
+  
+  //----------------------------------------
+  // initialize direction conversion tables
+  //----------------------------------------
+
+  for (int i=0; i<256; i++) {
+    double angle = double(i)*(1.0/256.0)*M_PI;
+    costheta[i] = cos( angle );
+    sintheta[i] = sin( angle );
+    cosphi[i]   = cos( 2.0*angle );
+    sinphi[i]   = sin( 2.0*angle );
+  }
+}
+
+
+//*************************
+PhotonMap :: ~PhotonMap()
+//*************************
+{
+  free( photons );
+}
+
+
+/* photon_dir returns the direction of a photon
+ */
+//*****************************************************************
+void PhotonMap :: photon_dir( float *dir, const Photon *p ) const
+//*****************************************************************
+{
+  dir[0] = sintheta[p->theta]*cosphi[p->phi];
+  dir[1] = sintheta[p->theta]*sinphi[p->phi];
+  dir[2] = costheta[p->theta];
+}
+
+
+/* irradiance_estimate computes an irradiance estimate
+ * at a given surface position
+*/
+//**********************************************
+void PhotonMap :: irradiance_estimate(
+  float irrad[3],                // returned irradiance
+  const float pos[3],            // surface position
+  const float normal[3],         // surface normal at pos
+  const float max_dist,          // max distance to look for photons
+  const int nphotons ) const     // number of photons to use
+//**********************************************
+{
+  irrad[0] = irrad[1] = irrad[2] = 0.0;
+
+  NearestPhotons np;
+  np.dist2 = (float*)alloca( sizeof(float)*(nphotons+1) );
+  np.index = (const Photon**)alloca( sizeof(Photon*)*(nphotons+1) );
+
+  np.pos[0] = pos[0]; np.pos[1] = pos[1]; np.pos[2] = pos[2];
+  np.max = nphotons;
+  np.found = 0;
+  np.got_heap = 0;
+  np.dist2[0] = max_dist*max_dist;
+
+  // locate the nearest photons
+  locate_photons( &np, 1 );
+
+  // if less than 8 photons return
+  if (np.found<8)
+    return;
+
+  float pdir[3];
+
+  // sum irradiance from all photons
+  for (int i=1; i<=np.found; i++) {
+    const Photon *p = np.index[i];
+    // the photon_dir call and following if can be omitted (for speed)
+    // if the scene does not have any thin surfaces
+    photon_dir( pdir, p );
+    if ( (pdir[0]*normal[0]+pdir[1]*normal[1]+pdir[2]*normal[2]) < 0.0f ) {
+      float red, green, blue;
+
+      rgbe2float(red, green, blue, p->power);
+      
+      irrad[0] += red;
+      irrad[1] += green;
+      irrad[2] += blue;
+    }
+  }
+
+  const float tmp=(1.0f/M_PI)/(np.dist2[0]);	// estimate of density
+
+  irrad[0] *= tmp;
+  irrad[1] *= tmp;
+  irrad[2] *= tmp;
+}
+
+
+/* locate_photons finds the nearest photons in the
+ * photon map given the parameters in np
+*/
+//******************************************
+void PhotonMap :: locate_photons(
+  NearestPhotons *const np,
+  const int index ) const
+//******************************************
+{
+  const Photon *p = &photons[index];
+  float dist1;
+
+  if (index<half_stored_photons) {
+    dist1 = np->pos[ p->plane ] - p->pos[ p->plane ];
+
+    if (dist1>0.0) { // if dist1 is positive search right plane
+      locate_photons( np, 2*index+1 );
+      if ( dist1*dist1 < np->dist2[0] )
+        locate_photons( np, 2*index );
+    } else {         // dist1 is negative search left first
+      locate_photons( np, 2*index );
+      if ( dist1*dist1 < np->dist2[0] )
+        locate_photons( np, 2*index+1 );
+    }
+  }
+
+  // compute squared distance between current photon and np->pos
+
+  dist1 = p->pos[0] - np->pos[0];
+  float dist2 = dist1*dist1;
+  dist1 = p->pos[1] - np->pos[1];
+  dist2 += dist1*dist1;
+  dist1 = p->pos[2] - np->pos[2];
+  dist2 += dist1*dist1;
+  
+  if ( dist2 < np->dist2[0] ) {
+    // we found a photon :) Insert it in the candidate list
+
+    if ( np->found < np->max ) {
+      // heap is not full; use array
+      np->found++;
+      np->dist2[np->found] = dist2;
+      np->index[np->found] = p;
+    } else {
+      int j,parent;
+
+      if (np->got_heap==0) { // Do we need to build the heap?
+        // Build heap
+        float dst2;
+        const Photon *phot;
+        int half_found = np->found>>1;
+        for ( int k=half_found; k>=1; k--) {
+          parent=k;
+          phot = np->index[k];
+          dst2 = np->dist2[k];
+          while ( parent <= half_found ) {
+            j = parent+parent;
+            if (j<np->found && np->dist2[j]<np->dist2[j+1])
+              j++;
+            if (dst2>=np->dist2[j])
+              break;
+            np->dist2[parent] = np->dist2[j];
+            np->index[parent] = np->index[j];
+            parent=j;
+          }
+          np->dist2[parent] = dst2;
+          np->index[parent] = phot;
+        }
+        np->got_heap = 1;
+      }
+
+      // insert new photon into max heap
+      // delete largest element, insert new and reorder the heap
+
+      parent=1;
+      j = 2;
+      while ( j <= np->found ) {
+        if ( j < np->found && np->dist2[j] < np->dist2[j+1] )
+          j++;
+        if ( dist2 > np->dist2[j] )
+          break;
+        np->dist2[parent] = np->dist2[j];
+        np->index[parent] = np->index[j];
+        parent = j;
+        j += j;
+      }
+      np->index[parent] = p;
+      np->dist2[parent] = dist2;
+
+      np->dist2[0] = np->dist2[1];
+    }
+  }
+}
+
+
+/* store puts a photon into the flat array that will form
+ * the final kd-tree.
+ *
+ * Call this function to store a photon.
+*/
+//***************************
+void PhotonMap :: store(
+  const float power[3],
+  const float pos[3],
+  const float dir[3],
+  const float ref_index)
+//***************************
+{
+  if (stored_photons>=max_photons)
+    return;
+
+  stored_photons++;
+  Photon *const node = &photons[stored_photons];
+
+  node->ref_index = ref_index;
+  
+  for (int i=0; i<3; i++) {
+    node->pos[i] = pos[i];
+
+    if (node->pos[i] < bbox_min[i])
+      bbox_min[i] = node->pos[i];
+    if (node->pos[i] > bbox_max[i])
+      bbox_max[i] = node->pos[i];
+
+    //node->power[i] = power[i];
+  }
+
+  float2rgbe(node->power, power[0], power[1], power[2]);
+
+  int theta = int( acos(dir[2])*(256.0/M_PI) );
+  if (theta>255)
+    node->theta = 255;
+  else
+   node->theta = (unsigned char)theta;
+
+  int phi = int( atan2(dir[1],dir[0])*(256.0/(2.0*M_PI)) );
+  if (phi>255)
+    node->phi = 255;
+  else if (phi<0)
+    node->phi = (unsigned char)(phi+256);
+  else
+    node->phi = (unsigned char)phi;
+}
+
+
+/* scale_photon_power is used to scale the power of all
+ * photons once they have been emitted from the light
+ * source. scale = 1/(#emitted photons).
+ * Call this function after each light source is processed.
+*/
+//********************************************************
+void PhotonMap :: scale_photon_power( const float scale )
+//********************************************************
+{
+  for (int i=prev_scale; i<=stored_photons; i++) {
+    float red, green, blue;
+    rgbe2float(red, green, blue, photons[i].power);
+    
+    red *= scale;
+    green *= scale;
+    blue *= scale;
+
+    float2rgbe(photons[i].power, red, green, blue);
+  }
+  prev_scale = stored_photons;
+}
+
+
+/* balance creates a left balanced kd-tree from the flat photon array.
+ * This function should be called before the photon map
+ * is used for rendering.
+ */
+//******************************
+void PhotonMap :: balance(void)
+//******************************
+{
+  if (stored_photons>1) {
+    // allocate two temporary arrays for the balancing procedure
+    Photon **pa1 = (Photon**)malloc(sizeof(Photon*)*(stored_photons+1));
+    Photon **pa2 = (Photon**)malloc(sizeof(Photon*)*(stored_photons+1));
+
+    for (int i=0; i<=stored_photons; i++)
+      pa2[i] = &photons[i];
+
+    balance_segment( pa1, pa2, 1, 1, stored_photons );
+    free(pa2);
+
+    // reorganize balanced kd-tree (make a heap)
+    int d, j=1, foo=1;
+    Photon foo_photon = photons[j];
+
+    for (int i=1; i<=stored_photons; i++) {
+      d=pa1[j]-photons;
+      pa1[j] = NULL;
+      if (d != foo)
+        photons[j] = photons[d];
+      else {
+        photons[j] = foo_photon;
+
+        if (i<stored_photons) {
+          for (;foo<=stored_photons; foo++)
+            if (pa1[foo] != NULL)
+              break;
+          foo_photon = photons[foo];
+          j = foo;
+        }
+        continue;
+      }
+      j = d;
+    }
+    free(pa1);
+  }
+
+  half_stored_photons = stored_photons/2-1;
+}
+
+
+#define swap(ph,a,b) { Photon *ph2=ph[a]; ph[a]=ph[b]; ph[b]=ph2; }
+
+// median_split splits the photon array into two separate
+// pieces around the median with all photons below the
+// the median in the lower half and all photons above
+// than the median in the upper half. The comparison
+// criteria is the axis (indicated by the axis parameter)
+// (inspired by routine in "Algorithms in C++" by Sedgewick)
+//*****************************************************************
+void PhotonMap :: median_split(
+  Photon **p,
+  const int start,               // start of photon block in array
+  const int end,                 // end of photon block in array
+  const int median,              // desired median number
+  const int axis )               // axis to split along
+//*****************************************************************
+{
+  int left = start;
+  int right = end;
+
+  while ( right > left ) {
+    const float v = p[right]->pos[axis];
+    int i=left-1;
+    int j=right;
+    for (;;) {
+      while ( p[++i]->pos[axis] < v )
+        ;
+      while ( p[--j]->pos[axis] > v && j>left )
+        ;
+      if ( i >= j )
+        break;
+      swap(p,i,j);
+    }
+
+    swap(p,i,right);
+    if ( i >= median )
+      right=i-1;
+    if ( i <= median )
+      left=i+1;
+  }
+}
+
+  
+// See "Realistic image synthesis using Photon Mapping" chapter 6
+// for an explanation of this function
+//****************************
+void PhotonMap :: balance_segment(
+  Photon **pbal,
+  Photon **porg,
+  const int index,
+  const int start,
+  const int end )
+//****************************
+{
+  //--------------------
+  // compute new median
+  //--------------------
+
+  int median=1;
+  while ((4*median) <= (end-start+1))
+    median += median;
+
+  if ((3*median) <= (end-start+1)) {
+    median += median;
+    median += start-1;
+  } else	
+    median = end-median+1;
+
+  //--------------------------
+  // find axis to split along
+  //--------------------------
+
+  int axis=2;
+  if ((bbox_max[0]-bbox_min[0])>(bbox_max[1]-bbox_min[1]) &&
+      (bbox_max[0]-bbox_min[0])>(bbox_max[2]-bbox_min[2]))
+    axis=0;
+  else if ((bbox_max[1]-bbox_min[1])>(bbox_max[2]-bbox_min[2]))
+    axis=1;
+
+  //------------------------------------------
+  // partition photon block around the median
+  //------------------------------------------
+
+  median_split( porg, start, end, median, axis );
+
+  pbal[ index ] = porg[ median ];
+  pbal[ index ]->plane = axis;
+
+  //----------------------------------------------
+  // recursively balance the left and right block
+  //----------------------------------------------
+
+  if ( median > start ) {
+    // balance left segment
+    if ( start < median-1 ) {
+      const float tmp=bbox_max[axis];
+      bbox_max[axis] = pbal[index]->pos[axis];
+      balance_segment( pbal, porg, 2*index, start, median-1 );
+      bbox_max[axis] = tmp;
+    } else {
+      pbal[ 2*index ] = porg[start];
+    }
+  }
+
+  if ( median < end ) {
+    // balance right segment
+    if ( median+1 < end ) {
+      const float tmp = bbox_min[axis];		
+      bbox_min[axis] = pbal[index]->pos[axis];
+      balance_segment( pbal, porg, 2*index+1, median+1, end );
+      bbox_min[axis] = tmp;
+    } else {
+      pbal[ 2*index+1 ] = porg[end];
+    }
+  }	
+}
+

photonmap.hpp

@@ -0,0 +1,103 @@
+#ifndef PHOTONMAP_H
+#define PHOTONMAP_H
+
+/* This is the photon
+ * The power is not compressed so the
+ * size is 28 bytes
+*/
+//**********************
+typedef struct Photon {
+//**********************
+  float pos[3];                  // photon position
+  short plane;                   // splitting plane for kd-tree
+  unsigned char theta, phi;      // incoming direction
+  //float power[3];                // photon power (uncompressed)
+  unsigned char power[4];
+  float ref_index;
+} Photon;
+
+
+/* This structure is used only to locate the
+ * nearest photons
+*/
+//**********************
+typedef struct NearestPhotons {
+//**********************
+    int max; 
+    int found; 
+    int got_heap; 
+    float pos[3]; 
+    float *dist2; 
+    const Photon **index; 
+} NearestPhotons;
+
+
+/* This is the Photon_map class
+ */
+//****************
+class PhotonMap {
+//****************
+ public:
+    PhotonMap(int max_phot);
+    ~PhotonMap();
+  
+    void store(
+      const float power[3],        // photon power
+      const float pos[3],          // photon position
+      const float dir[3],
+      const float ref_index);        // photon direction
+
+    void scale_photon_power(
+      const float scale);          // 1/(number of emitted photons)
+
+    void balance(void);            // balance the kd-tree (before use!)
+
+    void irradiance_estimate(
+      float irrad[3],              // returned irradiance
+      const float pos[3],          // surface position
+      const float normal[3],       // surface normal at pos
+      const float max_dist,        // max distance to look for photons
+      const int nphotons ) const;  // number of photons to use
+
+    void locate_photons(
+      NearestPhotons *const np,    // np is used to locate the photons
+      const int index) const;      // call with index = 1
+
+    void photon_dir(
+      float *dir,                  // direction of photon (returned)
+      const Photon *p) const;      // the photon
+
+ private:
+    friend class PhotonTracer;
+
+    void balance_segment(
+      Photon **pbal, 
+      Photon **porg, 
+      const int index,
+      const int start, 
+      const int end );
+
+    void median_split(
+      Photon **p, 
+      const int start, 
+      const int end,
+      const int median, 
+      const int axis );
+  
+    Photon *photons; 
+
+    int stored_photons; 
+    int half_stored_photons; 
+    int max_photons; 
+    int prev_scale; 
+
+    float costheta[256]; 
+    float sintheta[256]; 
+    float cosphi[256]; 
+    float sinphi[256]; 
+  
+    float bbox_min[3];     // use bbox_min
+    float bbox_max[3];     // use bbox_max
+};
+
+#endif