From c9818de11f1a5e7bfa1a383b30f0904c5af28fe9 Mon Sep 17 00:00:00 2001
From: Zack Buhman <zack@buhman.org>
Date: Wed, 25 Jan 2023 03:41:52 -0800
Subject: [PATCH] add specular

---
 math/fp.hpp    | 23 ++++++++++++++++++++++-
 math/vec.hpp   |  5 +++++
 raytracing.cpp | 27 +++++++++++++++++++++++----
 3 files changed, 50 insertions(+), 5 deletions(-)

diff --git a/math/fp.hpp b/math/fp.hpp
index 64f0023..df45583 100644
--- a/math/fp.hpp
+++ b/math/fp.hpp
@@ -28,6 +28,8 @@ struct fp
   inline constexpr fp<T, I, B>& operator+=(fp<T, I, B> const& v);
 
   inline constexpr fp<T, I, B>& operator-=(fp<T, I, B> const& v);
+
+  inline constexpr fp<T, I, B>& operator*=(fp<T, I, B> const& v);
 };
 
 template <typename T, typename I, int B>
@@ -51,6 +53,13 @@ inline constexpr fp<T, I, B>& fp<T, I, B>::operator-=(fp<T, I, B> const& v)
   return *this;
 }
 
+template <typename T, typename I, int B>
+inline constexpr fp<T, I, B>& fp<T, I, B>::operator*=(fp<T, I, B> const& v)
+{
+  *this = *this * v;
+  return *this;
+}
+
 template <typename T, typename I, int B>
 constexpr inline fp<T, I, B> operator+(const fp<T, I, B>& a, const fp<T, I, B>& b) noexcept
 {
@@ -151,11 +160,23 @@ struct fp_limits<fp<T, I, B>>
   }
 };
 
+// functions
+
+template <typename T, typename I, int B>
+constexpr inline fp<T, I, B> pow(fp<T, I, B> a, fp<T, I, B> b) noexcept
+{
+  while (b > fp<T, I, B>(1)) {
+    a *= a;
+    b -= fp<T, I, B>(1);
+  }
+  return a;
+}
+
 // specializations
 
 using fp16_16 = fp<int32_t, int64_t, 16>;
 
-constexpr fp16_16 sqrt(fp16_16 n) noexcept
+constexpr inline fp16_16 sqrt(const fp16_16& n) noexcept
 {
   int32_t x = n.value;
   int32_t c = 0;
diff --git a/math/vec.hpp b/math/vec.hpp
index 7587cfc..7ba21d5 100644
--- a/math/vec.hpp
+++ b/math/vec.hpp
@@ -25,6 +25,11 @@ struct vec<3, T>
 
   inline constexpr T const& operator[](int i) const;
 
+  constexpr inline vec<3, T> operator-() const noexcept
+  {
+    return vec<3, T>(-x, -y, -z);
+  }
+
   inline constexpr vec<3, T>& operator=(vec<3, T> const& v);
 
   inline constexpr vec<3, T>& operator+=(vec<3, T> const& v);
diff --git a/raytracing.cpp b/raytracing.cpp
index 7b324a1..61f8c00 100644
--- a/raytracing.cpp
+++ b/raytracing.cpp
@@ -22,6 +22,7 @@ struct sphere {
   vec3 center;
   fp16_16 radius;
   vec3 color;
+  fp16_16 specular;
 };
 
 enum class light_type {
@@ -48,23 +49,27 @@ constexpr scene scene {
   { // spheres
     {
       {0, -1, 3}, // center
-      fp16_16(1), // radius
+      1,          // radius
       {1, 0, 0},  // color
+      8,          // specular
     },
     {
       {2, 0, 4},
-      fp16_16(1),
+      1,
       {0, 0, 1},
+      10
     },
     {
       {-2, 0, 4},
       fp16_16(1),
       {0, 1, 0},
+      10,
     },
     {
       {0, -61, 0},
       fp16_16(60),
       {1, 1, 0},
+      0
     }
   },
   { // lights
@@ -95,7 +100,8 @@ struct t1_t2 {
   fp16_16 t2;
 };
 
-fp16_16 compute_lighting(const vec3& point, const vec3& normal)
+fp16_16 compute_lighting(const vec3& point, const vec3& normal,
+                         const vec3& viewer, fp16_16 specular)
 {
   fp16_16 intensity{0};
 
@@ -110,10 +116,22 @@ fp16_16 compute_lighting(const vec3& point, const vec3& normal)
       } else {
         light_vector = light.direction;
       }
+
+      // diffuse
       auto n_dot_l = dot(normal, light_vector);
       if (n_dot_l > fp16_16(0)) {
         intensity += light.intensity * n_dot_l * (fp16_16(1) / length(light_vector));
       }
+
+      // specular
+      if (specular > fp16_16(0)) {
+        auto reflected = normal * fp16_16(2) * dot(normal, light_vector) - light_vector;
+        auto r_dot_v = dot(reflected, viewer);
+        if (r_dot_v > fp16_16(0)) {
+          auto base = r_dot_v / (length(reflected) * length(viewer));
+          intensity += light.intensity * pow(base, specular);
+        }
+      }
     }
   }
   return intensity;
@@ -169,7 +187,8 @@ static vec3 trace_ray
     vec3 point = origin + direction * closest_t;
     vec3 normal = point - closest_sphere->center;
     normal = normal * (fp16_16(1) / length(normal));
-    return closest_sphere->color * compute_lighting(point, normal);
+    return closest_sphere->color * compute_lighting(point, normal,
+                                                    -direction, closest_sphere->specular);
   }
 }