Global diffuse lighting. (gamma not corrected)

Some operations on pixel<T>.
Make ray3 support copy semantic.
Fix vec3 operands does not filter out vec3-vec3 as parameters.
random_uv_gen generating random unit vectors.

6 files changed, 143 insertions, 45 deletions

diff --git a/aa.h b/aa.h
index 368a14a..ed3bc7f 100644
--- a/aa.h
+++ b/aa.h
@@ -44,13 +44,14 @@ public:
             remaining -= n;
             for (unsigned i = 0; i < n; ++i) {
                 workers.emplace_back(std::thread{
-                        [&](int tid, uint64_t seed, std::vector<basic_viewport<T>> *subs, vec3d viewpoint,
+                        [&](int tid, uint64_t seed, uint64_t diffuse_seed, std::vector<basic_viewport<T>> *subs, vec3d viewpoint,
                             uint16_t image_width, uint16_t image_height) {
                             bias_ctx bc{seed};
-                            auto image = (*subs)[tid].render(world, viewpoint, image_width, image_height, bc);
+                            auto image = (*subs)[tid].render(
+                                    world, viewpoint, image_width, image_height, bc, diffuse_seed);
                             images[base + tid] = image;
                         },
-                        i, seedgen(), subviews, viewpoint, image_width, image_height
+                        i, seedgen(), seedgen(), subviews, viewpoint, image_width, image_height
                 });
             }
             for (auto &th: workers) {
diff --git a/bitmap.h b/bitmap.h
index 1e35395..5b07914 100644
--- a/bitmap.h
+++ b/bitmap.h
@@ -44,8 +44,30 @@ struct pixel {
     static inline pixel<T> from_normalized(const vec3d &v3d) {
         return from_normalized(v3d.x / 2.0 + 0.5, v3d.y / 2.0 + 0.5, v3d.z / 2.0 + 0.5);
     }
+
+    static inline pixel<T> black() {
+        return pixel<T>{(T) 0, (T) 0, (T) 0};
+    }
 };
 
+template<
+        typename T,
+        typename S,
+        typename = typename std::enable_if<std::is_arithmetic<S>::value, S>::type
+>
+pixel<T> operator*(const pixel<T> &pixel, S scale) {
+    return ::pixel < T > {.r=(T) (pixel.r * scale), .g=(T) (pixel.g * scale), .b=(T) (pixel.b * scale)};
+}
+
+template<
+        typename T,
+        typename S,
+        typename = typename std::enable_if<std::is_arithmetic<S>::value, S>::type
+>
+pixel<T> operator*(S scale, const pixel <T> &pixel) {
+    return ::pixel < T > {.r=(T) (pixel.r * scale), .g=(T) (pixel.g * scale), .b=(T) (pixel.b * scale)};
+}
+
 // Mix two colors a and b. Returns a*u + b*v
 template<typename T>
 inline pixel<T> mix(const pixel<T> &a, const pixel<T> &b, double u, double v) {
@@ -93,7 +115,7 @@ public:
         bitmap<T> result{images[0].width, images[0].height};
         const auto m = images.size();
         const auto n = images[0].content.size();
-        uintmax_t acc_r, acc_g, acc_b;
+        uint_fast32_t acc_r, acc_g, acc_b;
         for (size_t i = 0; i < n; ++i) {
             acc_r = 0;
             acc_g = 0;
@@ -103,7 +125,7 @@ public:
                 acc_g += images[j].content[i].g;
                 acc_b += images[j].content[i].b;
             }
-            result.content[i] = pixel<T>{(T) (acc_r / m), (T) (acc_g / m), (T) (acc_b / m)};
+            result.content[i] = pixel<T>{(T) (1.0 * acc_r / m), (T) (1.0 * acc_g / m), (T) (1.0 * acc_b / m)};
         }
         return result;
     }
diff --git a/hitlist.h b/hitlist.h
index 0c65e1a..9b02819 100644
--- a/hitlist.h
+++ b/hitlist.h
@@ -18,6 +18,9 @@
 #include <iostream>
 #include <cstdint>
 
+//#define T_SIMPLE_COLOR
+//#define T_NORM_VISUAL
+#define T_DIFFUSE
 
 // A world, T is color depth
 template<typename T>
@@ -35,40 +38,60 @@ public:
     }
 
     // Given a ray, compute the color.
-    pixel<T> color(const ray3d &r) const {
-        // Detect hits
-        bool hit = false;
-        double hit_t = std::numeric_limits<double>::infinity();
-        std::shared_ptr<object> hit_obj;
-        // Check the nearest object we hit
-        for (const auto &obj: objects) {
-            double t_;
-            if (obj->hit(r, t_, 0.0) && t_ < hit_t) {
-                hit = true;
-                hit_t = t_;
-                hit_obj = obj;
+    pixel<T> color(ray3d r, random_uv_gen_3d &ruvg, uint_fast32_t max_recursion_depth = 64) const {
+        double decay = 1;
+        while (max_recursion_depth-- > 0) {
+            // Detect hits
+            bool hit = false;
+            double hit_t = std::numeric_limits<double>::infinity();
+            std::shared_ptr<object> hit_obj;
+            // Check the nearest object we hit
+            for (const auto &obj: objects) {
+                double t_;
+                if (obj->hit(r, t_, 0.0) && t_ < hit_t) {
+                    hit = true;
+                    hit_t = t_;
+                    hit_obj = obj;
+                }
+            }
+            if (hit) {
+#ifdef T_SIMPLE_COLOR
+                // simply returns color of the object
+                return hit_obj->color();
+#endif
+#ifdef T_NORM_VISUAL
+                // normal vector on hit point
+                const auto nv = hit_obj->normal_vector(r.at(hit_t));
+                // visualize normal vector at hit point
+                return pixel<T>::from_normalized(nv);
+#endif
+#ifdef T_DIFFUSE
+                const auto hit_point = r.at(hit_t); // hit point, on the surface
+                auto nv = hit_obj->normal_vector(hit_point);
+                if (dot(nv, r.direction()) > 0) return pixel<T>::black(); // discard rays from inner (or invert nv)
+                vec3d diffuse_target = hit_point + nv + ruvg();
+                decay *= 0.5; // lose 50% light when diffused
+                r = ray3d{hit_point, diffuse_target - hit_point}; // the new diffused ray we trace on
+                continue;
+#endif
             }
-        }
-        if (hit) {
-            // normal vector on hit point
-            const auto nv = hit_obj->normal_vector(r.at(hit_t));
-//                return obj->color();
-            // visualize normal vector at hit point
-            return pixel<T>::from_normalized(nv);
-        }
-
 
-        // Does not hit anything. Get background color (infinity)
-        const auto u = (r.direction().y + 1.0) * 0.5;
-        return mix(
-                pixel<T>::from_normalized(1.0, 1.0, 1.0),
-                pixel<T>::from_normalized(0.5, 0.7, 1.0),
-                1.0 - u,
-                u
-        );
+            // Does not hit anything. Get background color (infinity)
+            const auto u = (r.direction().y + 1.0) * 0.5;
+            const auto c = mix(
+                    pixel<T>::from_normalized(1.0, 1.0, 1.0),
+                    pixel<T>::from_normalized(0.5, 0.7, 1.0),
+                    1.0 - u,
+                    u
+            );
+#ifdef T_DIFFUSE
+            return decay * c;
+#else
+            return c;
+#endif
+        }
+        return pixel<T>::black(); // reached recursion time limit, very little light
     }
-
-
 };
 
 using hitlist8b = hitlist<uint8_t>;
diff --git a/ray.h b/ray.h
index 5a68557..db28d0f 100644
--- a/ray.h
+++ b/ray.h
@@ -10,11 +10,22 @@
 // A ray is a half-line, starts from a 3d point, along the direction of a unit vector, to the infinity
 template<typename T>
 class ray3 {
-    const vec3<T> source_;
-    const vec3<T> direction_; // unit vector
+    vec3<T> source_;
+    vec3<T> direction_; // unit vector
 
 public:
+    ~ray3() = default;
+
     ray3() = delete;
+
+    ray3(const ray3<T> &other) = default;
+
+    ray3<T> &operator=(const ray3<T> &other) {
+        source_ = other.source_;
+        direction_ = other.direction_;
+        return *this;
+    }
+
     ray3(const vec3<T> &source, const vec3<T> &direction) : source_(source), direction_(direction.unit_vec()) {}
 
     // Get the source point from where the ray emits.
@@ -32,6 +43,11 @@ public:
     vec3<T> at(U t) const {
         return source_ + direction_ * t;
     }
+
+    // Get a ray starts from zero and directs to undefined direction.
+    static ray3<T> null() {
+        return ray3<T>{vec3<T>::zero(), vec3<T>::zero()};
+    }
 };
 
 using ray3d = ray3<double>;
diff --git a/vec.h b/vec.h
index bf40e7d..2c8c2f4 100644
--- a/vec.h
+++ b/vec.h
@@ -6,6 +6,7 @@
 #define RT_VEC_H
 
 #include <cmath>
+#include <random>
 
 static inline bool eq(int a, int b) {
     return a == b;
@@ -88,19 +89,31 @@ inline std::ostream &operator<<(std::ostream &out, const vec3<T> &vec) {
 }
 
 // product vec3 by a scalar
-template<typename T, typename S>
+template<
+        typename T,
+        typename S,
+        typename = typename std::enable_if<std::is_arithmetic<S>::value, S>::type
+>
 inline vec3<T> operator*(const vec3<T> &vec, const S &b) {
     return vec3<T>{.x=(T) (vec.x * b), .y=(T) (vec.y * b), .z=(T) (vec.z * b)};
 }
 
 // product vec3 by a scalar
-template<typename T, typename S>
+template<
+        typename T,
+        typename S,
+        typename = typename std::enable_if<std::is_arithmetic<S>::value, S>::type
+>
 inline vec3<T> operator*(const S &b, const vec3<T> &vec) {
     return vec3<T>{.x=(T) (vec.x * b), .y=(T) (vec.y * b), .z=(T) (vec.z * b)};
 }
 
 // product vec3 by the inversion of a scalar (div by a scalar)
-template<typename T, typename S>
+template<
+        typename T,
+        typename S,
+        typename = typename std::enable_if<std::is_arithmetic<S>::value, S>::type
+>
 inline vec3<T> operator/(const vec3<T> &vec, const S &b) {
     return vec3<T>{.x=(T) (vec.x / b), .y=(T) (vec.y / b), .z=(T) (vec.z / b)};
 }
@@ -129,4 +142,24 @@ using vec3f = vec3<float>;
 // 3-dim vector (double)
 using vec3d = vec3<double>;
 
+// random unit vector generator
+template<typename T>
+class random_uv_gen {
+    std::mt19937_64 mt;
+    std::uniform_real_distribution<T> uni{-1.0, 1.0};
+
+public:
+    explicit random_uv_gen(uint64_t seed) : mt{seed} {}
+
+    vec3<T> operator()() {
+        while (true) {
+            const auto x = uni(mt), y = uni(mt), z = uni(mt);
+            const auto vec = vec3<T>{.x=x, .y=y, .z=z};
+            if (vec.mod2() <= 1.0) return vec;
+        }
+    }
+};
+
+using random_uv_gen_3d = random_uv_gen<double>;
+
 #endif //RT_VEC_H
\ No newline at end of file
diff --git a/viewport.h b/viewport.h
index 74bdfcd..30a22dd 100644
--- a/viewport.h
+++ b/viewport.h
@@ -44,6 +44,7 @@ template<typename T>
 class viewport {
 public:
     virtual bitmap<T> render(const hitlist<T> &world, vec3d viewpoint, uint16_t image_width, uint16_t image_height) = 0;
+
     virtual ~viewport() = default;
 };
 
@@ -64,7 +65,8 @@ public:
     virtual bitmap<T>
     render(const hitlist<T> &world, vec3d viewpoint, uint16_t image_width, uint16_t image_height) override {
         bias_ctx bc{};
-        return render(world, viewpoint, image_width, image_height, bc);
+        static constexpr uint64_t default_diffuse_seed = 123456789012345678ULL;
+        return render(world, viewpoint, image_width, image_height, bc, default_diffuse_seed);
     }
 
     /**
@@ -75,8 +77,9 @@ public:
      */
     virtual bitmap<T> render(const hitlist<T> &world, vec3d viewpoint,
                              uint16_t image_width, uint16_t image_height,
-                             bias_ctx &bias) const {
+                             bias_ctx &bias, uint64_t diffuse_seed) const {
         bitmap<T> image{image_width, image_height};
+        random_uv_gen_3d ruvg{diffuse_seed};
         double bx, by;
         const auto r = center - viewpoint;
         const int img_hw = image_width / 2, img_hh = image_height / 2;
@@ -94,8 +97,8 @@ public:
                         .z=0.0
                 }; // offset on screen plane
                 const auto dir = r + off; // direction vector from camera to current pixel on screen
-                const ray3d ray{viewpoint, dir}; // from camera to pixel (on the viewport)
-                const auto pixel = world.color(ray);
+                ray3d ray{viewpoint, dir}; // from camera to pixel (on the viewport)
+                const auto pixel = world.color(ray, ruvg);
                 image.set(i + img_hw, -j + img_hh, pixel);
             }
         }