From c97a2d4fe3be5cb447c23208791690f9fb1e9f29 Mon Sep 17 00:00:00 2001 From: Miguel Angel Astor Romero Date: Tue, 27 Jun 2017 10:38:12 -0400 Subject: [PATCH] Added photon map implementation from the book. --- .clang_complete | 1 + Makefile | 5 +- main.cpp | 30 ++- path_tracer.cpp | 3 +- photon_tracer.cpp | 89 ++++++--- photon_tracer.hpp | 90 ++++++++- photonmap.cpp | 464 ++++++++++++++++++++++++++++++++++++++++++++++ photonmap.hpp | 103 ++++++++++ 8 files changed, 747 insertions(+), 38 deletions(-) create mode 100644 photonmap.cpp create mode 100644 photonmap.hpp diff --git a/.clang_complete b/.clang_complete index 3936014..4e4f8df 100644 --- a/.clang_complete +++ b/.clang_complete @@ -4,3 +4,4 @@ -DGLM_FORCE_RADIANS -DUSE_CPP11_RANDOM -fopenmp +-fno-builtin diff --git a/Makefile b/Makefile index a3fd96a..b6dd867 100644 --- a/Makefile +++ b/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 diff --git a/main.cpp b/main.cpp index 3f97fb6..b7fade3 100644 --- a/main.cpp +++ b/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)); @@ -403,6 +409,26 @@ void parse_args(int argc, char ** const argv) { exit(EXIT_FAILURE); } 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(optopt) << "\" requires an argument." << endl; diff --git a/path_tracer.cpp b/path_tracer.cpp index 9a6f894..534ab3b 100644 --- a/path_tracer.cpp +++ b/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())) + - (_f->m_mat->m_specular * dir_spec_color)); + color += ((dir_diff_color + ind_color + amb_color) * (_f->m_mat->m_diffuse / pi())) + (_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) { diff --git a/photon_tracer.cpp b/photon_tracer.cpp index 2713b10..65cbd34 100644 --- a/photon_tracer.cpp +++ b/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 photons; - vector caustics; + vector photons; + vector caustics; t = numeric_limits::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() * (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() * (radius * radius); + c_contrib /= (1.0f - (2.0f / (3.0f * m_cone_filter_k))) * pi() * (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() * (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
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(n_photons_per_ligth)); + power = (al->m_figure->m_mat->m_emission /* / static_cast(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(n_photons_per_ligth)); + power = (pl->m_diffuse /* / static_cast(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((static_cast(current) / static_cast(total)) * 100.0) << "% done."; } cout << endl; - cout << "Generated " << ANSI_BOLD_YELLOW << m_photon_map.getNumPhotons() << ANSI_RESET_STYLE << " total photons." << endl; - m_photon_map.save_photon_list(specular ? "caustics.txt" : "photons.txt"); + 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"); + + 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); - m_photon_map.addPhoton(ph); + ph = PhotonAux(Vec3(x, y, z), Vec3(dx, dy, dz), r, g, b, rc); + //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) { - Photon photon; +void PhotonTracer::trace_photon(PhotonAux & ph, Scene * s, const unsigned int rec_level) { + 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); - m_photon_map.addPhoton(photon); + photon = PhotonAux(p_pos, p_dir, red, green, blue, ph.ref_index); + //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())); 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); } } diff --git a/photon_tracer.hpp b/photon_tracer.hpp index b9c6538..b68b664 100644 --- a/photon_tracer.hpp +++ b/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" -class PhotonTracer: public Tracer { +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::epsilon() <= y) && (x + std::numeric_limits::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(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) { }; + PhotonTracer(): + 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_caustics_map; - void trace_photon(Photon & ph, Scene * s, const unsigned int rec_level); + /*kdTree m_photon_map; + kdTree m_caustics_map;*/ + PhotonMap m_photon_map; + int m_max_s_photons; + void trace_photon(PhotonAux & ph, Scene * s, const unsigned int rec_level); }; #endif diff --git a/photonmap.cpp b/photonmap.cpp new file mode 100644 index 0000000..6c559a3 --- /dev/null +++ b/photonmap.cpp @@ -0,0 +1,464 @@ +//---------------------------------------------------------------------------- +// photonmap.cc +// An example implementation of the photon map data structure +// +// Henrik Wann Jensen - February 2001 +//---------------------------------------------------------------------------- + +#include +#include +#include +#include +#include + +#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 (indexpos[ 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 (jfound && np->dist2[j]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 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]; + } + } +} + diff --git a/photonmap.hpp b/photonmap.hpp new file mode 100644 index 0000000..1bd3eaa --- /dev/null +++ b/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