Use simd for vec3 computations

9 files changed, 154 insertions, 125 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index b8c6ee1..6aaa0f0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -8,7 +8,7 @@ set(CMAKE_TOOLCHAIN_FILE "${CMAKE_CURRENT_SOURCE_DIR}/vcpkg/scripts/buildsystems
 
 if(WIN32)
 add_compile_options(
-    -fcolor-diagnostics -fansi-escape-codes # colored output when using ninja
+    -msse4.1 -fcolor-diagnostics -fansi-escape-codes # colored output when using ninja
 )
 endif()
 
diff --git a/include/aabb.hpp b/include/aabb.hpp
index cdbb79a..5f819db 100644
--- a/include/aabb.hpp
+++ b/include/aabb.hpp
@@ -16,41 +16,47 @@ class aabb {
             auto &v2 = tri.v2.position;

             auto &v3 = tri.v3.position;

 

-            minimum.x = fmin(fmin(v1.x, v2.x), v3.x);

-            minimum.y = fmin(fmin(v1.y, v2.y), v3.y);

-            minimum.z = fmin(fmin(v1.z, v2.z), v3.z);

+            // minimum.x = fmin(fmin(v1.x, v2.x), v3.x);

+            // minimum.y = fmin(fmin(v1.y, v2.y), v3.y);

+            // minimum.z = fmin(fmin(v1.z, v2.z), v3.z);

 

-            maximum.x = fmax(fmax(v1.x, v2.x), v3.x);

-            maximum.y = fmax(fmax(v1.y, v2.y), v3.y);

-            maximum.z = fmax(fmax(v1.z, v2.z), v3.z);

+            minimum.v[0] = fmin(fmin(v1.v[0], v2.v[0]), v3.v[0]);

+            minimum.v[1] = fmin(fmin(v1.v[1], v2.v[1]), v3.v[1]);

+            minimum.v[2] = fmin(fmin(v1.v[2], v2.v[2]), v3.v[2]);

+

+            // maximum.x = fmax(fmax(v1.x, v2.x), v3.x);

+            // maximum.y = fmax(fmax(v1.y, v2.y), v3.y);

+            // maximum.z = fmax(fmax(v1.z, v2.z), v3.z);

+

+            maximum.v[0] = fmax(fmax(v1.v[0], v2.v[0]), v3.v[0]);

+            maximum.v[1] = fmax(fmax(v1.v[1], v2.v[1]), v3.v[1]);

+            maximum.v[2] = fmax(fmax(v1.v[2], v2.v[2]), v3.v[2]);

         }

 

         void union_with(const aabb &box) {

-            minimum.x = minimum.x < box.minimum.x ? minimum.x : box.minimum.x;

-            minimum.y = minimum.y < box.minimum.y ? minimum.y : box.minimum.y;

-            minimum.z = minimum.z < box.minimum.z ? minimum.z : box.minimum.z;

-

-            maximum.x = maximum.x > box.maximum.x ? maximum.x : box.maximum.x;

-            maximum.y = maximum.y > box.maximum.y ? maximum.y : box.maximum.y;

-            maximum.z = maximum.z > box.maximum.z ? maximum.z : box.maximum.z;

+            minimum.min(box.minimum);

+            maximum.max(box.maximum);

         }

 

         void union_with(const point3 &point) {

-            minimum.x = minimum.x < point.x ? minimum.x : point.x;

-            minimum.y = minimum.y < point.y ? minimum.y : point.y;

-            minimum.z = minimum.z < point.z ? minimum.z : point.z;

-

-            maximum.x = maximum.x > point.x ? maximum.x : point.x;

-            maximum.y = maximum.y > point.y ? maximum.y : point.y;

-            maximum.z = maximum.z > point.z ? maximum.z : point.z;

+            minimum.min(point);

+            maximum.max(point);

         }

 

         int32_t maximum_axis() {

             vec3 diag = diagonal();

 

-            if (diag.x > diag.y && diag.x > diag.z) {

+            // if (diag.x > diag.y && diag.x > diag.z) {

+            //     return 0;

+            // } else if (diag.y > diag.z) {

+            //     return 1;

+            // } else {

+            //     return 2;

+            // }

+

+            if (diag.v[0] > diag.v[1] && diag.v[0] > diag.v[2]) {

                 return 0;

-            } else if (diag.y > diag.z) {

+            } else if (diag.v[1] > diag.v[2]) {

                 return 1;

             } else {

                 return 2;

@@ -67,7 +73,8 @@ class aabb {
 

         float surface_area() {

             auto diag = diagonal();

-            return 2.f * (diag.x * diag.y + diag.x * diag.z + diag.y * diag.z);

+            // return 2.f * (diag.x * diag.y + diag.x * diag.z + diag.y * diag.z);

+            return 2.f * (diag.v[0] * diag.v[1] + diag.v[0] * diag.v[2] + diag.v[1] * diag.v[2]);

         }

 

         vec3 minimum = vec3(std::numeric_limits<float>::max());

diff --git a/include/bvh-node.hpp b/include/bvh-node.hpp
index 09c9420..698977a 100644
--- a/include/bvh-node.hpp
+++ b/include/bvh-node.hpp
@@ -91,13 +91,13 @@ public:
 

             node.primitive_id = old_node.primitive_id;

 

-            node.min[0] = old_node.bounding_box.minimum.x;

-            node.min[1] = old_node.bounding_box.minimum.y;

-            node.min[2] = old_node.bounding_box.minimum.z;

+            node.min[0] = old_node.bounding_box.minimum.v[0];

+            node.min[1] = old_node.bounding_box.minimum.v[1];

+            node.min[2] = old_node.bounding_box.minimum.v[2];

 

-            node.max[0] = old_node.bounding_box.maximum.x;

-            node.max[1] = old_node.bounding_box.maximum.y;

-            node.max[2] = old_node.bounding_box.maximum.z;

+            node.max[0] = old_node.bounding_box.maximum.v[0];

+            node.max[1] = old_node.bounding_box.maximum.v[1];

+            node.max[2] = old_node.bounding_box.maximum.v[2];

 

             if (old_node.next_id != -1) {

                 node.next_id = temp_nodes[old_node.next_id].df_id;

diff --git a/include/utils.hpp b/include/utils.hpp
index e7178d8..8548a5f 100644
--- a/include/utils.hpp
+++ b/include/utils.hpp
@@ -9,6 +9,7 @@ std::vector<uint8_t> read_file(const char* path);
 std::vector<uint8_t> read_file(std::string& path);
 
 #define PI 3.14159265359
+#define EPSILON 0.000001
 
 inline float deg_to_rad(float deg) {
     return deg * PI / 180.0;
diff --git a/include/vec3.hpp b/include/vec3.hpp
index 6291e30..429d747 100644
--- a/include/vec3.hpp
+++ b/include/vec3.hpp
@@ -2,6 +2,7 @@
 #define __VEC3_HPP_
 
 #include <iostream>
+#include <xmmintrin.h>
 
 class vec3 {
     public:
@@ -9,10 +10,15 @@ class vec3 {
 
     vec3(float x);
     vec3(float x, float y, float z);
+    vec3(__m128& v);
+    vec3(__m128&& v);
+    vec3(const vec3& vec);
 
     vec3 operator-() const;
 
-    vec3& operator+=(vec3 vec);
+    vec3& operator=(const vec3& vec);
+
+    vec3& operator+=(const vec3& vec);
 
     vec3& operator*=(float scale);
 
@@ -28,7 +34,11 @@ class vec3 {
 
     vec3 cross(const vec3& vec) const;
 
-    vec3 unit() const;
+    vec3& min(const vec3& vec);
+
+    vec3& max(const vec3& vec);
+
+    vec3& normalize();
 
     bool near_zero() const;
 
@@ -36,38 +46,29 @@ class vec3 {
 
     static vec3 random(float min, float max);
 
-    float x, y, z;
-    float padding;
+    __m128 v;
 };
 
-vec3 lerp(vec3 u, vec3 v, float t);
+vec3 lerp(const vec3& u, const vec3& v, float t);
 
-vec3 operator+(vec3 u, vec3 v);
+vec3 operator+(const vec3& u, const vec3& v);
 
-vec3 operator-(vec3 u, vec3 v);
+vec3 operator-(const vec3& u, const vec3& v);
 
-vec3 operator*(vec3 u, vec3 v);
+vec3 operator*(const vec3& u, const vec3& v);
 
-vec3 operator*(vec3 u, float scale);
+vec3 operator*(const vec3& u, float scale);
 
-vec3 operator*(float scale, vec3 u);
+vec3 operator*(float scale, const vec3& u);
 
-vec3 operator/(vec3 u, float scale);
+vec3 operator/(const vec3& u, float scale);
 
-vec3 reflect(vec3 v, vec3 n);
+vec3 reflect(const vec3& v, const vec3& n);
 
 vec3 refract(const vec3& uv, const vec3& n, float etai_over_etat);
 
 std::ostream& operator<<(std::ostream& out, vec3 vec);
 
-vec3 random_in_unit_sphere();
-
-vec3 random_unit_vector();
-
-vec3 random_in_hemisphere(const vec3& normal);
-
-vec3 random_in_unit_disk();
-
 using point3 = vec3;
 using color = vec3;
 
diff --git a/src/camera.cpp b/src/camera.cpp
index 917cb4e..908ba09 100644
--- a/src/camera.cpp
+++ b/src/camera.cpp
@@ -12,9 +12,9 @@ camera::camera(
         float focus_dist
     ) : position(position), lens_radius(aperture / 2.0), fov(v_fov), focus_distance(focus_dist), aspect_ratio(aspect_ratio) {
 
-    forward = (target - position).unit();
-    right = forward.cross(vec3{ 0.0, 1.0, 0.0 }).unit();
-    up = right.cross(forward).unit();
+    forward = (target - position).normalize();
+    right = forward.cross(vec3{ 0.0, 1.0, 0.0 }).normalize();
+    up = right.cross(forward).normalize();
 
     auto h = tan(deg_to_rad(fov) / 2.0);
     auto viewport_height = 2.0 * h;
@@ -45,11 +45,12 @@ void camera::move(vec3 v) {
 }
 
 void camera::rotate_y(float theta) {
-    forward.x = forward.x * cos(theta) - forward.z * sin(theta);
-    forward.z = forward.x * sin(theta) + forward.z * cos(theta);
+    forward.v[0] = forward.v[0] * cos(theta) - forward.v[2] * sin(theta);
+    forward.v[2] = forward.v[0] * sin(theta) + forward.v[2] * cos(theta);
+    forward.normalize();
 
-    right = forward.cross(vec3{ 0.0, 1.0, 0.0 }).unit();
-    up = right.cross(forward).unit();
+    right = forward.cross(vec3{ 0.0, 1.0, 0.0 }).normalize();
+    up = right.cross(forward).normalize();
 
     auto h = tan(deg_to_rad(fov) / 2.0);
     auto viewport_height = 2.0 * h;
diff --git a/src/color.cpp b/src/color.cpp
index 92195de..555ff8d 100644
--- a/src/color.cpp
+++ b/src/color.cpp
@@ -2,9 +2,9 @@
 #include "utils.hpp"
 
 void write_color(std::ostream& out, const color& pixel_color) {
-    auto r = pixel_color.x;
-    auto g = pixel_color.y;
-    auto b = pixel_color.z;
+    auto r = pixel_color.v[0];
+    auto g = pixel_color.v[1];
+    auto b = pixel_color.v[2];
 
     out << static_cast<uint32_t>(256.0 * clamp(r, 0.0, 0.999)) << ' '
         << static_cast<uint32_t>(256.0 * clamp(g, 0.0, 0.999)) << ' '
diff --git a/src/main.cpp b/src/main.cpp
index 8d79000..37ce060 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -48,6 +48,8 @@ int main() {
     uint32_t height = width / aspect_ratio;
     const uint32_t samples_per_pixel = 500;
 
+    std::cerr << sizeof(vec3) << std::endl;
+
     window wnd { width , height };
 
     float delta_time = 0.f;
@@ -124,19 +126,19 @@ int main() {
 
                     switch (event.key) {
                         case KEYS::W: {
-                            move_vec = ((scene*) renderer.scene_buffer_ptr())->meta.cam.forward.unit() * move_speed;
+                            move_vec = ((scene*) renderer.scene_buffer_ptr())->meta.cam.forward * move_speed;
                             break;
                         }
                         case KEYS::S: {
-                            move_vec = -((scene*) renderer.scene_buffer_ptr())->meta.cam.forward.unit() * move_speed;
+                            move_vec = -((scene*) renderer.scene_buffer_ptr())->meta.cam.forward * move_speed;
                             break;
                         }
                         case KEYS::D: {
-                            move_vec = ((scene*) renderer.scene_buffer_ptr())->meta.cam.right.unit() * move_speed;
+                            move_vec = ((scene*) renderer.scene_buffer_ptr())->meta.cam.right * move_speed;
                             break;
                         }
                         case KEYS::A: {
-                            move_vec = -((scene*) renderer.scene_buffer_ptr())->meta.cam.right.unit() * move_speed;
+                            move_vec = -((scene*) renderer.scene_buffer_ptr())->meta.cam.right * move_speed;
                             break;
                         }
                         case KEYS::E: {
@@ -150,11 +152,11 @@ int main() {
                             break;
                         }
                         case KEYS::SPACE: {
-                            move_vec = ((scene*) renderer.scene_buffer_ptr())->meta.cam.up.unit() * move_speed;
+                            move_vec = ((scene*) renderer.scene_buffer_ptr())->meta.cam.up * move_speed;
                             break;
                         }
                         case KEYS::LCTRL: {
-                            move_vec = -((scene*) renderer.scene_buffer_ptr())->meta.cam.up.unit() * move_speed;
+                            move_vec = -((scene*) renderer.scene_buffer_ptr())->meta.cam.up * move_speed;
                             break;
                         }
                         default:
diff --git a/src/vec3.cpp b/src/vec3.cpp
index 166f681..d075473 100644
--- a/src/vec3.cpp
+++ b/src/vec3.cpp
@@ -1,87 +1,101 @@
 #include <cmath>
 #include <limits.h>
+#include <smmintrin.h>
 
 #include "vec3.hpp"
 #include "utils.hpp"
 
-vec3::vec3(): x(0.0), y(0.0), z(0.0) {}
+vec3::vec3(): v(_mm_setzero_ps()) {}
+vec3::vec3(float x): v(_mm_set1_ps(x)) {}
+vec3::vec3(float x, float y, float z): v(_mm_set_ps(0.f, z, y, x)) {}
 
-vec3::vec3(float x): x(x), y(x), z(x) {}
-vec3::vec3(float x, float y, float z): x(x), y(y), z(z) {}
+vec3::vec3(__m128& v): v(v) {}
+vec3::vec3(__m128&& v): v(std::move(v)) {}
+vec3::vec3(const vec3& vec) {
+    v = _mm_load_ps((float*)&vec.v);
+}
 
 vec3 vec3::operator-() const {
-    return vec3 {
-        -x,
-        -y,
-        -z
-    };
+    return vec3(_mm_xor_ps(v, _mm_set1_ps(-0.0f)));
+}
+
+vec3& vec3::operator=(const vec3& vec) {
+    v = _mm_load_ps((float*)&vec.v);
+
+    return *this;
 }
 
-vec3& vec3::operator+=(vec3 vec) {
-    x += vec.x;
-    y += vec.y;
-    z += vec.z;
+vec3& vec3::operator+=(const vec3& vec) {
+    v = _mm_add_ps(v, vec.v);
 
     return *this;
 }
 
 vec3& vec3::operator*=(float scale) {
-    x *= scale;
-    y *= scale;
-    z *= scale;
+    __m128 s = _mm_set1_ps(scale);
+    v = _mm_mul_ps(v, s);
 
     return *this;
 }
 
 vec3& vec3::operator/=(float numerator) {
-    x /= numerator;
-    y /= numerator;
-    z /= numerator;
+    __m128 s = _mm_set1_ps(numerator);
+    v = _mm_div_ps(v, s);
 
     return *this;
 }
 
 float vec3::operator[](int axis) {
-    switch (axis) {
-        case 0:
-            return x;
-        case 1:
-            return y;
-        case 2:
-            return z;
-        default:
-           return std::numeric_limits<float>::infinity();
+    if (axis > 2) {
+        return std::numeric_limits<float>::infinity();
     }
+
+    return v[axis];
 }
 
 float vec3::length() const {
-    return sqrt(length_sq());
+    return std::sqrt((_mm_dp_ps(v, v, 0b10001111))[0]);
 }
 
 float vec3::length_sq() const {
-    return x * x + y * y + z * z;
+    return (_mm_dp_ps(v, v, 0b10001111))[0];
 }
 
 float vec3::dot(const vec3& vec) const {
-    return x * vec.x + y * vec.y + z * vec.z;
+    return (_mm_dp_ps(v, vec.v, 0b10001111))[0];
 }
 
 vec3 vec3::cross(const vec3& vec) const {
-    return vec3(
-        y * vec.z - z * vec.y,
-        z * vec.x - x * vec.z,
-        x * vec.y - y * vec.x
-    );
+    __m128 tmp0 = _mm_shuffle_ps(v,v,_MM_SHUFFLE(3,0,2,1));
+    __m128 tmp1 = _mm_shuffle_ps(vec.v,vec.v,_MM_SHUFFLE(3,1,0,2));
+    __m128 tmp2 = _mm_mul_ps(tmp0,vec.v);
+    __m128 tmp3 = _mm_mul_ps(tmp0,tmp1);
+    __m128 tmp4 = _mm_shuffle_ps(tmp2,tmp2,_MM_SHUFFLE(3,0,2,1));
+    return _mm_sub_ps(tmp3,tmp4);
 }
 
-vec3 vec3::unit() const {
-    return *this / length();
+vec3& vec3::min(const vec3& vec) {
+    v = _mm_min_ps(v, vec.v);
+
+    return *this;
 }
 
-bool vec3::near_zero() const {
-    const auto s = 1e-8;
+vec3& vec3::max(const vec3& vec) {
+    v = _mm_max_ps(v, vec.v);
 
-    return (fabs(x) < s) && (fabs(y) < s) && (fabs(z) < s);
+    return *this;
+}
+
+vec3& vec3::normalize() {
+    v = _mm_div_ps(v, _mm_sqrt_ps(_mm_dp_ps(v, v, 0xff)));
+
+    return *this;
+}
+
+bool vec3::near_zero() const {
+    // __m128 e = _mm_set1_ps(EPSILON);
+    // _mm_cmplt_ps(_mm_andnot_ps(_mm_set1_ps(-0.0f), v), e);
+    return (fabs(v[0]) < EPSILON) && (fabs(v[1]) < EPSILON) && (fabs(v[2]) < EPSILON);
 }
 
 vec3 vec3::random() {
@@ -92,42 +106,45 @@ vec3 vec3::random(float min, float max) {
     return vec3 { randd(min, max), randd(min, max), randd(min, max) };
 }
 
-vec3 lerp(vec3 u, vec3 v, float t) {
+
+// TODO: Use simd instructions
+vec3 lerp(const vec3& u, const vec3& v, float t) {
     return (1.0 - t) * u + t * v;
 }
 
-vec3 operator+(vec3 u, vec3 v) {
-    return vec3(u.x + v.x, u.y + v.y, u.z + v.z);
+vec3 operator+(const vec3& u, const vec3& v) {
+    return vec3(_mm_add_ps(u.v, v.v));
 }
 
-vec3 operator-(vec3 u, vec3 v) {
-    return vec3(u.x - v.x, u.y - v.y, u.z - v.z);
+vec3 operator-(const vec3& u, const vec3& v) {
+    return vec3(_mm_sub_ps(u.v, v.v));
 }
 
-vec3 operator*(vec3 u, vec3 v) {
-    return vec3(u.x * v.x, u.y * v.y, u.z * v.z);
+vec3 operator*(const vec3& u, const vec3& v) {
+    return vec3(_mm_mul_ps(u.v, v.v));
 }
 
-vec3 operator*(vec3 u, float scale) {
-    return vec3(
-        u.x * scale,
-        u.y * scale,
-        u.z * scale
-    );
+vec3 operator*(const vec3& u, float scale) {
+    __m128 s = _mm_set1_ps(scale);
+    return vec3(_mm_mul_ps(u.v, s));
 }
 
-vec3 operator*(float scale, vec3 u) {
-    return u * scale;
+vec3 operator*(float scale, const vec3& u) {
+    __m128 s = _mm_set1_ps(scale);
+    return vec3(_mm_mul_ps(u.v, s));
 }
 
-vec3 operator/(vec3 u, float scale) {
-    return (1.0 / scale) * u;
+vec3 operator/(const vec3& u, float scale) {
+    __m128 s = _mm_set1_ps(scale);
+    return vec3(_mm_div_ps(u.v, s));
 }
 
-vec3 reflect(vec3 v, vec3 n) {
+// TODO: Use simd instructions
+vec3 reflect(const vec3& v, const vec3& n) {
     return v - 2 * v.dot(n) * n;
 }
 
+// TODO: Use simd instructions
 vec3 refract(const vec3& uv, const vec3& n, float etai_over_etat) {
     auto cos_theta = fmin(n.dot(-uv), 1.0);
     vec3 r_out_perp =  etai_over_etat * (uv + cos_theta*n);
@@ -137,5 +154,5 @@ vec3 refract(const vec3& uv, const vec3& n, float etai_over_etat) {
 }
 
 std::ostream& operator<<(std::ostream& out, vec3 vec) {
-    return out << vec.x << ' ' << vec.y << ' ' << vec.z << std::endl;
+    return out << vec[0] << ' ' << vec[1] << ' ' << vec[2] << std::endl;
 }