From 9b3d3615e3fbc1f4961432ce0d77276608ab3e96 Mon Sep 17 00:00:00 2001
From: kuhyx <147418882+kuhyx@users.noreply.github.com>
Date: Sat, 18 Jan 2025 15:16:46 +0100
Subject: [PATCH] Added cpp photonmapping (#3)

* feat: load exr map

* feat: scale and flip the hdr map

* feat: added camera

* feat: add simple slow but working photonmapping

* feat: added cpp version of photonmapping
---
 code/photonmapping/cpp/compile.sh |   1 +
 code/photonmapping/cpp/main.cpp   | 329 ++++++++++++++++++++++++++++++
 code/photonmapping/python/main.py | 221 ++++++++++++++++++++
 3 files changed, 551 insertions(+)
 create mode 100755 code/photonmapping/cpp/compile.sh
 create mode 100644 code/photonmapping/cpp/main.cpp
 create mode 100644 code/photonmapping/python/main.py

diff --git a/code/photonmapping/cpp/compile.sh b/code/photonmapping/cpp/compile.sh
new file mode 100755
index 00000000..5f567f3f
--- /dev/null
+++ b/code/photonmapping/cpp/compile.sh
@@ -0,0 +1 @@
+g++ -std=c++11 -O2 main.cpp -o photon_mapping
\ No newline at end of file
diff --git a/code/photonmapping/cpp/main.cpp b/code/photonmapping/cpp/main.cpp
new file mode 100644
index 00000000..6a6cd629
--- /dev/null
+++ b/code/photonmapping/cpp/main.cpp
@@ -0,0 +1,329 @@
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <random>
+#include <fstream>
+
+// Define basic vector operations
+class Vector3 {
+public:
+    double x, y, z;
+    
+    Vector3(double x_=0, double y_=0, double z_=0): x(x_), y(y_), z(z_) {}
+    
+    Vector3 operator + (const Vector3& v) const { return Vector3(x+v.x, y+v.y, z+v.z); }
+    Vector3 operator - (const Vector3& v) const { return Vector3(x-v.x, y-v.y, z-v.z); }
+    Vector3 operator * (double scalar) const { return Vector3(x*scalar, y*scalar, z*scalar); }
+    Vector3 operator / (double scalar) const { return Vector3(x/scalar, y/scalar, z/scalar); } // Added operator/
+    Vector3 operator - () const { return Vector3(-x, -y, -z); } // Added unary minus operator
+    Vector3 operator * (const Vector3& v) const { return Vector3(x*v.x, y*v.y, z*v.z); } // Component-wise multiplication
+    double dot(const Vector3& v) const { return x*v.x + y*v.y + z*v.z; }
+    double norm() const { return std::sqrt(x*x + y*y + z*z); }
+    Vector3 normalize() const { double n = norm(); return Vector3(x/n, y/n, z/n); }
+};
+
+typedef Vector3 Color; // Alias for RGB color
+
+// Define the photon
+struct Photon {
+    Vector3 position;
+    Vector3 direction;
+    Color power;
+    
+    Photon(const Vector3& pos, const Vector3& dir, const Color& pow)
+        : position(pos), direction(dir), power(pow) {}
+};
+
+// Define a simple sphere
+class Sphere {
+public:
+    Vector3 center;
+    double radius;
+    Color color;
+    
+    Sphere(const Vector3& c, double r, const Color& col)
+        : center(c), radius(r), color(col) {}
+    
+    bool intersect(const Vector3& ray_origin, const Vector3& ray_direction, double& t, Vector3& hit_point, Vector3& normal) const {
+        // Solve quadratic equation for intersection
+        Vector3 oc = ray_origin - center;
+        double a = ray_direction.dot(ray_direction);
+        double b = 2.0 * oc.dot(ray_direction);
+        double c = oc.dot(oc) - radius * radius;
+        double discriminant = b*b - 4*a*c;
+        if (discriminant < 0) {
+            return false; // No intersection
+        } else {
+            double sqrt_discriminant = std::sqrt(discriminant);
+            double t0 = (-b - sqrt_discriminant) / (2.0 * a);
+            double t1 = (-b + sqrt_discriminant) / (2.0 * a);
+            t = (t0 < t1 && t0 > 1e-4) ? t0 : t1;
+            if (t < 1e-4) {
+                return false;
+            }
+            hit_point = ray_origin + ray_direction * t;
+            normal = (hit_point - center).normalize();
+            return true;
+        }
+    }
+};
+
+// Define a simple plane
+class Plane {
+public:
+    Vector3 point;   // A point on the plane
+    Vector3 normal;  // Normal to the plane
+    Color color;
+    
+    Plane(const Vector3& p, const Vector3& n, const Color& col)
+        : point(p), normal(n.normalize()), color(col) {}
+    
+    bool intersect(const Vector3& ray_origin, const Vector3& ray_direction, double& t, Vector3& hit_point, Vector3& hit_normal) const {
+        double denom = normal.dot(ray_direction);
+        if (std::abs(denom) > 1e-6) {
+            t = (point - ray_origin).dot(normal) / denom;
+            if (t >= 1e-4) {
+                hit_point = ray_origin + ray_direction * t;
+                hit_normal = normal;
+                return true;
+            }
+        }
+        return false;
+    }
+};
+
+// Random number generators
+std::mt19937 rng;
+std::uniform_real_distribution<double> uni_dist(0.0, 1.0);
+
+// Random unit vector in a sphere
+Vector3 random_unit_vector() {
+    double theta = 2 * M_PI * uni_dist(rng);
+    double z = 2 * uni_dist(rng) - 1;
+    double r = std::sqrt(1 - z*z);
+    return Vector3(r * std::cos(theta), r * std::sin(theta), z);
+}
+
+// Random direction in hemisphere around normal
+Vector3 random_hemisphere_direction(const Vector3& normal) {
+    Vector3 dir = random_unit_vector();
+    if (dir.dot(normal) < 0) {
+        dir = dir * -1;
+    }
+    return dir;
+}
+
+// Global variables
+std::vector<Photon> photon_map;
+std::vector<void*> objects;  // Pointers to objects
+std::vector<int> object_types; // 0 for Sphere, 1 for Plane
+
+// Scene setup
+Sphere sphere(Vector3(0, 0, -5), 1.0, Color(1, 0, 0));         // Red sphere
+Plane plane(Vector3(0, -1, 0), Vector3(0, 1, 0), Color(0.5, 0.5, 0.5)); // Gray plane
+
+// Light source
+Vector3 light_position(-5, 5, -5);
+Color light_power(1000.0, 1000.0, 1000.0);  // Intense white light, will scale in code
+
+// Parameters
+int num_photons = 10000;   // Number of photons to emit
+int max_depth = 5;         // Maximum number of bounces
+double gather_radius = 0.5; // Radius for radiance estimation
+
+// Functions
+void emit_photons();
+void trace_photon(Photon photon, int depth);
+Color trace_ray(const Vector3& ray_origin, const Vector3& ray_direction);
+Color compute_direct_light(const Vector3& point, const Vector3& normal);
+Color estimate_radiance(const Vector3& point, const Vector3& normal);
+
+// Main execution
+int main() {
+    // Seed random number generator
+    rng.seed(std::random_device()());
+    
+    std::cout << "Emitting photons..." << std::endl;
+    // Add objects to the scene
+    objects.push_back(&sphere); object_types.push_back(0);
+    objects.push_back(&plane);  object_types.push_back(1);
+    
+    emit_photons();
+    std::cout << "Photons stored in photon map: " << photon_map.size() << std::endl;
+    
+    std::cout << "Rendering image..." << std::endl;
+    int width = 200;
+    int height = 100;
+    std::vector<Color> image(width * height);
+    
+    double aspect_ratio = double(width) / height;
+    double fov = M_PI / 3.0;  // 60 degrees field of view
+    
+    for (int y = 0; y < height; ++y) {
+        for (int x = 0; x < width; ++x) {
+            // Convert pixel coordinate to camera ray
+            double px = (2 * (x + 0.5) / double(width) - 1) * tan(fov / 2.0) * aspect_ratio;
+            double py = (1 - 2 * (y + 0.5) / double(height)) * tan(fov / 2.0);
+            Vector3 ray_origin(0, 0, 0);
+            Vector3 ray_direction = Vector3(px, py, -1).normalize();
+            Color color = trace_ray(ray_origin, ray_direction);
+            image[y * width + x] = Color(
+                std::min(color.x, 1.0),
+                std::min(color.y, 1.0),
+                std::min(color.z, 1.0)
+            );
+        }
+    }
+    
+    // Save image to PPM file
+    std::ofstream ofs("render.ppm");
+    ofs << "P3\n" << width << " " << height << "\n255\n";
+    double gamma = 1.0 / 2.2; // Corrected gamma division
+    for (const auto& pixel : image) {
+        int r = static_cast<int>(std::pow(pixel.x, gamma) * 255);
+        int g = static_cast<int>(std::pow(pixel.y, gamma) * 255);
+        int b = static_cast<int>(std::pow(pixel.z, gamma) * 255);
+        // Clamp values between 0 and 255
+        r = std::max(0, std::min(255, r));
+        g = std::max(0, std::min(255, g));
+        b = std::max(0, std::min(255, b));
+        ofs << r << " " << g << " " << b << "\n";
+    }
+    ofs.close();
+    std::cout << "Image saved to render.ppm" << std::endl;
+    return 0;
+}
+
+void emit_photons() {
+    Color per_photon_power = light_power / num_photons; // Scale light power per photon
+    for (int i = 0; i < num_photons; ++i) {
+        // Emit photons in random directions from the light source
+        Vector3 direction = random_unit_vector();
+        Photon photon(light_position, direction, per_photon_power);
+        trace_photon(photon, 0);
+    }
+}
+
+void trace_photon(Photon photon, int depth) {
+    if (depth > max_depth) {
+        return;
+    }
+    double closest_t = std::numeric_limits<double>::infinity();
+    void* hit_object = nullptr;
+    int hit_type = -1;
+    Vector3 hit_point, normal;
+    // Find the nearest intersection
+    for (size_t i = 0; i < objects.size(); ++i) {
+        double t;
+        Vector3 temp_hit_point, temp_normal;
+        bool hit = false;
+        if (object_types[i] == 0) {
+            Sphere* obj = static_cast<Sphere*>(objects[i]);
+            hit = obj->intersect(photon.position, photon.direction, t, temp_hit_point, temp_normal);
+        } else if (object_types[i] == 1) {
+            Plane* obj = static_cast<Plane*>(objects[i]);
+            hit = obj->intersect(photon.position, photon.direction, t, temp_hit_point, temp_normal);
+        }
+        if (hit && t < closest_t) {
+            closest_t = t;
+            hit_object = objects[i];
+            hit_type = object_types[i];
+            hit_point = temp_hit_point;
+            normal = temp_normal;
+        }
+    }
+    if (hit_object) {
+        photon_map.push_back(Photon(hit_point, photon.direction, photon.power));
+        // Diffuse reflection
+        Vector3 new_direction = random_hemisphere_direction(normal);
+        photon.position = hit_point;
+        photon.direction = new_direction;
+        // Absorb some power
+        photon.power = photon.power * 0.8; // Simple absorption
+        trace_photon(photon, depth + 1);
+    }
+}
+
+Color trace_ray(const Vector3& ray_origin, const Vector3& ray_direction) {
+    double closest_t = std::numeric_limits<double>::infinity();
+    void* hit_object = nullptr;
+    int hit_type = -1;
+    Vector3 hit_point, normal;
+    Color obj_color;
+    // Find the nearest intersection
+    for (size_t i = 0; i < objects.size(); ++i) {
+        double t;
+        Vector3 temp_hit_point, temp_normal;
+        bool hit = false;
+        if (object_types[i] == 0) {
+            Sphere* obj = static_cast<Sphere*>(objects[i]);
+            hit = obj->intersect(ray_origin, ray_direction, t, temp_hit_point, temp_normal);
+            if (hit) obj_color = obj->color;
+        } else if (object_types[i] == 1) {
+            Plane* obj = static_cast<Plane*>(objects[i]);
+            hit = obj->intersect(ray_origin, ray_direction, t, temp_hit_point, temp_normal);
+            if (hit) obj_color = obj->color;
+        }
+        if (hit && t < closest_t) {
+            closest_t = t;
+            hit_object = objects[i];
+            hit_type = object_types[i];
+            hit_point = temp_hit_point;
+            normal = temp_normal;
+        }
+    }
+    if (hit_object) {
+        Color direct_light = compute_direct_light(hit_point, normal);
+        Color indirect_light = estimate_radiance(hit_point, normal);
+        return obj_color * (direct_light + indirect_light); // Component-wise multiplication
+    } else {
+        return Color(0, 0, 0); // Background color
+    }
+}
+
+Color compute_direct_light(const Vector3& point, const Vector3& normal) {
+    // Simple Lambertian reflection from light source
+    Vector3 direction_to_light = (light_position - point).normalize();
+    // Shadow ray
+    Vector3 shadow_origin = point + normal * 1e-5;
+    Vector3 shadow_ray = direction_to_light;
+    bool in_shadow = false;
+    for (size_t i = 0; i < objects.size(); ++i) {
+        double t;
+        Vector3 temp_hit_point, temp_normal;
+        bool hit = false;
+        if (object_types[i] == 0) {
+            Sphere* obj = static_cast<Sphere*>(objects[i]);
+            hit = obj->intersect(shadow_origin, shadow_ray, t, temp_hit_point, temp_normal);
+        } else if (object_types[i] == 1) {
+            Plane* obj = static_cast<Plane*>(objects[i]);
+            hit = obj->intersect(shadow_origin, shadow_ray, t, temp_hit_point, temp_normal);
+        }
+        if (hit) {
+            in_shadow = true;
+            break;
+        }
+    }
+    if (in_shadow) {
+        return Color(0, 0, 0);
+    } else {
+        double intensity = std::max(0.0, normal.dot(direction_to_light));
+        double distance2 = (light_position - point).dot(light_position - point);
+        return (light_power * intensity) / (4 * M_PI * distance2);
+    }
+}
+
+Color estimate_radiance(const Vector3& point, const Vector3& normal) {
+    // Gather photons within the gather_radius
+    Color accumulated_power(0.0, 0.0, 0.0);
+    for (const auto& photon : photon_map) {
+        double distance = (photon.position - point).norm();
+        if (distance < gather_radius) {
+            double weight = std::max(0.0, normal.dot((-photon.direction).normalize()));
+            accumulated_power = accumulated_power + photon.power * weight;
+        }
+    }
+    double area = M_PI * gather_radius * gather_radius;
+    return accumulated_power / area; // Removed division by num_photons
+}
\ No newline at end of file
diff --git a/code/photonmapping/python/main.py b/code/photonmapping/python/main.py
new file mode 100644
index 00000000..9aec82d6
--- /dev/null
+++ b/code/photonmapping/python/main.py
@@ -0,0 +1,221 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Define basic vector operations
+class Vector3:
+    def __init__(self, x, y, z):
+        self.x = x
+        self.y = y
+        self.z = z
+        
+    def __add__(self, other):
+        return Vector3(self.x+other.x, self.y+other.y, self.z+other.z)
+    
+    def __sub__(self, other):
+        return Vector3(self.x-other.x, self.y-other.y, self.z-other.z)
+    
+    def __mul__(self, scalar):
+        return Vector3(self.x*scalar, self.y*scalar, self.z*scalar)
+    
+    def dot(self, other):
+        return self.x*other.x + self.y*other.y + self.z*other.z
+    
+    def norm(self):
+        return np.sqrt(self.dot(self))
+    
+    def normalize(self):
+        n = self.norm()
+        return Vector3(self.x/n, self.y/n, self.z/n)
+
+# Define the photon
+class Photon:
+    def __init__(self, position, direction, power):
+        self.position = position
+        self.direction = direction
+        self.power = power
+
+# Define a simple sphere
+class Sphere:
+    def __init__(self, center, radius, color):
+        self.center = center
+        self.radius = radius
+        self.color = color
+
+    def intersect(self, ray_origin, ray_direction):
+        # Solve quadratic equation for intersection
+        oc = ray_origin - self.center
+        a = ray_direction.dot(ray_direction)
+        b = 2.0 * oc.dot(ray_direction)
+        c = oc.dot(oc) - self.radius*self.radius
+        discriminant = b*b - 4*a*c
+        if discriminant < 0:
+            return None  # No intersection
+        else:
+            t = (-b - np.sqrt(discriminant)) / (2.0*a)
+            if t < 0:
+                t = (-b + np.sqrt(discriminant)) / (2.0*a)
+            if t < 0:
+                return None
+            hit_point = ray_origin + ray_direction * t
+            normal = (hit_point - self.center).normalize()
+            return (t, hit_point, normal)
+        
+# Define a simple plane
+class Plane:
+    def __init__(self, point, normal, color):
+        self.point = point
+        self.normal = normal.normalize()
+        self.color = color
+        
+    def intersect(self, ray_origin, ray_direction):
+        denom = self.normal.dot(ray_direction)
+        if abs(denom) > 1e-6:
+            t = (self.point - ray_origin).dot(self.normal) / denom
+            if t >= 0:
+                hit_point = ray_origin + ray_direction * t
+                return (t, hit_point, self.normal)
+        return None
+
+# Scene setup
+sphere = Sphere(Vector3(0, 0, -5), 1.0, np.array([1, 0, 0]))  # Red sphere
+plane = Plane(Vector3(0, -1, 0), Vector3(0, 1, 0), np.array([0.5, 0.5, 0.5]))  # Gray plane
+
+objects = [sphere, plane]
+
+# Light source
+light_position = Vector3(-5, 5, -5)
+light_power = np.array([1, 1, 1]) * 1000  # Intense white light
+
+# Photon map
+photon_map = []
+
+# Parameters
+num_photons = 10000  # Number of photons to emit
+max_depth = 5  # Maximum number of bounces
+gather_radius = 0.5  # Radius for radiance estimation
+
+def emit_photons():
+    for _ in range(num_photons):
+        # Emit photons in random directions from the light source
+        direction = random_unit_vector()
+        power = light_power / num_photons
+        photon = Photon(light_position, direction, power)
+        trace_photon(photon, 0)
+
+def trace_photon(photon, depth):
+    if depth > max_depth:
+        return
+    closest_t = np.inf
+    hit_object = None
+    hit_info = None
+    # Find the nearest intersection
+    for obj in objects:
+        result = obj.intersect(photon.position, photon.direction)
+        if result:
+            t, hit_point, normal = result
+            if t < closest_t:
+                closest_t = t
+                hit_object = obj
+                hit_info = (hit_point, normal)
+    if hit_object:
+        hit_point, normal = hit_info
+        photon_map.append((hit_point, photon.power))
+        # Diffuse reflection
+        new_direction = random_hemisphere_direction(normal)
+        photon.position = hit_point
+        photon.direction = new_direction
+        # Absorb some power
+        photon.power = photon.power * 0.8  # Simple absorption
+        trace_photon(photon, depth+1)
+
+def random_unit_vector():
+    theta = np.random.uniform(0, 2*np.pi)
+    z = np.random.uniform(-1, 1)
+    r = np.sqrt(1 - z*z)
+    return Vector3(r * np.cos(theta), r * np.sin(theta), z)
+
+def random_hemisphere_direction(normal):
+    dir = random_unit_vector()
+    if dir.dot(normal) < 0:
+        dir = Vector3(-dir.x, -dir.y, -dir.z)
+    return dir
+
+def render_image(width, height):
+    aspect_ratio = width / height
+    fov = np.pi / 3  # 60 degrees field of view
+    image = np.zeros((height, width, 3))
+    for y in range(height):
+        for x in range(width):
+            # Convert pixel coordinate to camera ray
+            px = (2 * (x + 0.5) / width - 1) * np.tan(fov / 2) * aspect_ratio
+            py = (1 - 2 * (y + 0.5) / height) * np.tan(fov / 2)
+            ray_origin = Vector3(0, 0, 0)
+            ray_direction = Vector3(px, py, -1).normalize()
+            color = trace_ray(ray_origin, ray_direction)
+            image[y, x, :] = np.clip(color, 0, 1)
+    return image
+
+def trace_ray(ray_origin, ray_direction):
+    closest_t = np.inf
+    hit_object = None
+    hit_info = None
+    # Find the nearest intersection
+    for obj in objects:
+        result = obj.intersect(ray_origin, ray_direction)
+        if result:
+            t, hit_point, normal = result
+            if t < closest_t:
+                closest_t = t
+                hit_object = obj
+                hit_info = (hit_point, normal, obj.color)
+    if hit_object:
+        hit_point, normal, color = hit_info
+        direct_light = compute_direct_light(hit_point, normal)
+        indirect_light = estimate_radiance(hit_point, normal)
+        return color * (direct_light + indirect_light)
+    else:
+        return np.array([0, 0, 0])  # Background color
+
+def compute_direct_light(point, normal):
+    # Simple Lambertian reflection from light source
+    direction_to_light = (light_position - point).normalize()
+    # Shadow ray
+    shadow_origin = point + normal * 1e-5
+    shadow_ray = direction_to_light
+    in_shadow = False
+    for obj in objects:
+        result = obj.intersect(shadow_origin, shadow_ray)
+        if result:
+            in_shadow = True
+            break
+    if in_shadow:
+        return np.array([0, 0, 0])
+    else:
+        intensity = max(0, normal.dot(direction_to_light))
+        return intensity * light_power / (4 * np.pi * (light_position - point).norm()**2)
+
+def estimate_radiance(point, normal):
+    # Gather photons within the gather_radius
+    accumulated_power = np.array([0.0, 0.0, 0.0])
+    for photon_pos, photon_power in photon_map:
+        distance = (photon_pos - point).norm()
+        if distance < gather_radius:
+            weight = max(0, normal.dot((photon_pos - point).normalize()))
+            accumulated_power += photon_power * weight
+    area = np.pi * gather_radius ** 2
+    return accumulated_power / (area * num_photons)
+
+# Main execution
+if __name__ == '__main__':
+    print("Emitting photons...")
+    emit_photons()
+    
+    print("Rendering image...")
+    width = 200
+    height = 100
+    image = render_image(width, height)
+    
+    # Display the image
+    plt.imshow(image)
+    plt.axis('off')
+    plt.show()