Initial commit

.gitignore

@@ -0,0 +1,5 @@
+.DS_Store
+/build/*
+/dist/*
+!/dist/*.xyz
+!/dist/*.txt

CMakeLists.txt

@@ -0,0 +1,6 @@
+project(cloud_tracer)
+
+file(GLOB SOURCES "src/*.cpp")
+
+set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/dist)
+add_executable(cloud_tracer ${SOURCES})

src/defer.hpp

@@ -0,0 +1,33 @@
+#pragma once
+#include <functional>
+
+template<typename F>
+class defer_finalizer {
+    F f;
+    bool moved;
+public:
+    template<typename T>
+    defer_finalizer(T && f_) : f(std::forward<T>(f_)), moved(false) { }
+
+    defer_finalizer(const defer_finalizer &) = delete;
+
+    defer_finalizer(defer_finalizer && other) : f(std::move(other.f)), moved(other.moved) {
+        other.moved = true;
+    }
+
+    ~defer_finalizer() {
+        if (!moved) f();
+    }
+};
+
+static struct {
+    template<typename F>
+    defer_finalizer<F> operator<<(F && f) {
+        return defer_finalizer<F>(std::forward<F>(f));
+    }
+} deferrer;
+
+
+#define concat(x, y) x ## y
+#define label(x, y) concat(x, y)
+#define defer auto label(__deferred_lambda_call, __COUNTER__) = deferrer << [&]

src/dt_pathtrace.cpp

@@ -0,0 +1,307 @@
+#include "dt_pathtrace.hpp"
+#include <cstdio>
+#include <cmath>
+
+#define INFO(...) printf("INFO: "  __VA_ARGS__);
+
+#ifndef pi
+#  define pi 3.1415926535897932384626433832795
+#  define piOverTwo 1.5707963267948966192313216916398
+#  define inverseOfPi 0.31830988618379067153776752674503
+#  define inverseOfTwoPi 0.15915494309189533576888376337251
+#  define two_pi 6.283185307179586476925286766559
+#endif
+namespace Random {
+    struct rng_state {
+        uint32_t x;
+        uint32_t y;
+        uint32_t z;
+        uint32_t w;
+    };
+
+    static rng_state RNG_STATE;
+
+    uint TausStep(uint z, int S1, int S2, int S3, uint M) {
+        uint b = (((z << S1) ^ z) >> S2);
+        return ((z & M) << S3) ^ b;
+    }
+    uint LCGStep(uint z, uint A, uint C) {
+        return A * z + C;
+    }
+
+    float HybridTaus() {
+        RNG_STATE.x = TausStep(RNG_STATE.x, 13, 19, 12, 4294967294);
+        RNG_STATE.y = TausStep(RNG_STATE.y, 2, 25, 4, 4294967288);
+        RNG_STATE.z = TausStep(RNG_STATE.z, 3, 11, 17, 4294967280);
+        RNG_STATE.w = LCGStep(RNG_STATE.w, 1664525, 1013904223);
+
+        return 2.3283064365387e-10 * (RNG_STATE.x ^ RNG_STATE.y ^ RNG_STATE.z ^ RNG_STATE.w);
+    }
+
+    float random() {
+        return HybridTaus();
+    }
+
+    void InitRNG(uint32_t seed) {
+
+        RNG_STATE.x = seed;
+        RNG_STATE.y = seed;
+        RNG_STATE.z = seed;
+        RNG_STATE.w = seed;
+
+        for (int i = 0; i < 23 + seed % 13; i++)
+            random();
+    }
+
+}
+inline double max(double a, double b) {
+    return a > b ? a : b;
+}
+
+inline double min(double a, double b) {
+    return a < b ? a : b;
+}
+
+inline float max(float a, float b) {
+    return a > b ? a : b;
+}
+
+inline float min(float a, float b) {
+    return a < b ? a : b;
+}
+
+inline float max(float a, float b, float c) {
+    return max(a, max(b, c));
+}
+
+inline float min(float a, float b, float c) {
+    return min(a, min(b, c));
+}
+
+inline uint32_t max(uint32_t a, uint32_t b) {
+    return a > b ? a : b;
+}
+
+inline uint32_t min(uint32_t a, uint32_t b) {
+    return a < b ? a : b;
+}
+
+inline float Lerp(float a, float t, float b) {
+    return a + t * (b - a);
+}
+
+float Sample3DTexture(texture3D Grid, float3 P) {
+    // cubic interpolation
+
+#define IDX(x, y, z) ((x * Grid.h * Grid.d) + (y * Grid.d) + z)
+
+    float fw = Grid.w * P.x;
+    float fh = Grid.h * P.y;
+    float fd = Grid.d * P.z;
+
+    float tw = fw - (int)fw;
+    float th = fh - (int)fh;
+    float td = fd - (int)fd;
+
+    uint32_t idx_w_left = floor(fw);
+    uint32_t idx_h_low  = floor(fh);
+    uint32_t idx_d_near = floor(fd);
+
+    uint32_t idx_w_right = ceil(fw);
+    uint32_t idx_h_high  = ceil(fh);
+    uint32_t idx_d_far   = ceil(fd);
+
+    float left_low_near   = Grid.Data[IDX(idx_w_left,  idx_h_low,  idx_d_near)];
+    float left_low_far    = Grid.Data[IDX(idx_w_left,  idx_h_low,  idx_d_far)];
+    float left_high_near  = Grid.Data[IDX(idx_w_left,  idx_h_high, idx_d_near)];
+    float left_high_far   = Grid.Data[IDX(idx_w_left,  idx_h_high, idx_d_far)];
+    float right_low_near  = Grid.Data[IDX(idx_w_right, idx_h_low,  idx_d_near)];
+    float right_low_far   = Grid.Data[IDX(idx_w_right, idx_h_low,  idx_d_far)];
+    float right_high_near = Grid.Data[IDX(idx_w_right, idx_h_high, idx_d_near)];
+    float right_high_far  = Grid.Data[IDX(idx_w_right, idx_h_high, idx_d_far)];
+
+    float low_near = Lerp(left_low_near, tw, right_low_near);
+    float low_far  = Lerp(left_low_far,  tw, right_low_far);
+    float low      = Lerp(low_near,      td, low_far);
+
+    float high_near = Lerp(left_high_near, tw, right_high_near);
+    float high_far  = Lerp(left_high_far,  tw, right_high_far);
+    float high      = Lerp(high_near,      td, high_far);
+
+    return Lerp(low, th, high);
+
+#undef IDX
+}
+
+void CreateOrthonormalBasis(float3 D, float3& B, float3& T) {
+    float3 other = fabs(D.z) >= 0.999 ? float3{1, 0, 0} : float3{0, 0, 1};
+    B = Normalize(Cross(other, D));
+    T = Normalize(Cross(D, B));
+}
+
+void CreateOrthonormalBasis2(float3 D, float3& B, float3& T) {
+    float3 other = fabs(D.z) >= 0.999 ? float3{1, 0, 0} : float3{0, 0, 1};
+    B = Normalize(Cross(other, D));
+    T = Normalize(Cross(D, B));
+}
+
+float3 RandomDirection(float3 D) {
+    float r1 = Random::random();
+    float r2 = Random::random() * 2 - 1;
+    float sqrR2 = r2 * r2;
+    float two_pi_by_r1 = two_pi * r1;
+    float sqrt_of_one_minus_sqrR2 = sqrt(1.0 - sqrR2);
+    float x = cos(two_pi_by_r1) * sqrt_of_one_minus_sqrR2;
+    float y = sin(two_pi_by_r1) * sqrt_of_one_minus_sqrR2;
+    float z = r2;
+
+    float3 t0, t1;
+    CreateOrthonormalBasis2(D, t0, t1);
+
+    return t0 * x + t1 * y + D * z;
+}
+
+float invertcdf(float GFactor, float xi) {
+#define one_minus_g2 (1.0 - (GFactor) * (GFactor))
+#define one_plus_g2 (1.0 + (GFactor) * (GFactor))
+#define one_over_2g (0.5 / (GFactor))
+
+    float t = (one_minus_g2) / (1.0f - GFactor + 2.0f * GFactor * xi);
+    return one_over_2g * (one_plus_g2 - t * t);
+
+#undef one_minus_g2
+#undef one_plus_g2
+#undef one_over_2g
+}
+
+float3 ImportanceSamplePhase(float GFactor, float3 D) {
+    if (fabs(GFactor) < 0.001) {
+        return RandomDirection(-D);
+    }
+
+    float phi = Random::random() * 2 * pi;
+    float cosTheta = invertcdf(GFactor, Random::random());
+    float sinTheta = sqrt(max(0, 1.0f - cosTheta * cosTheta));
+
+    float3 t0, t1;
+    CreateOrthonormalBasis(D, t0, t1);
+
+    return sinTheta * sin(phi) * t0 + sinTheta * cos(phi) * t1 +
+        cosTheta * D;
+}
+
+void GetGridBox(texture3D grid, float3& minim, float3& maxim) {
+    float maxDim = max(grid.w, grid.h, grid.d);
+    maxim = float3{(float)grid.w, (float)grid.h, (float)grid.d} / maxDim * 0.5;
+    minim = -maxim;
+}
+
+bool BoxIntersect(float3 bMin, float3 bMax, float3 P, float3 D, float& tMin, float& tMax) {
+    float2x3 C  = float2x3{ bMin - P, bMax - P };
+    float2x3 D2 = float2x3{ D, D };
+
+    float2x3 T =
+        Abs(D2) <= 0.000001
+        ? float2x3{float3{-1000, -1000, -1000}, float3{1000, 1000, 1000}}
+        : C / D2;
+
+    tMin = max(min(T.v1.v[0], T.v2.v[0]),
+               min(T.v1.v[1], T.v2.v[1]),
+               min(T.v1.v[2], T.v2.v[2]));
+
+    tMin = max(0.0f, tMin);
+
+    tMax = min(max(T.v1.v[0], T.v2.v[0]),
+               max(T.v1.v[1], T.v2.v[1]),
+               max(T.v1.v[2], T.v2.v[2]));
+
+    if (tMax < tMin || tMax < 0) {
+        return false;
+    }
+    return true;
+}
+
+void DTPathtrace(path_info PathInfo, volume_info VolumeInfo, ray_pack* RayPack) {
+
+    int PassNumber = PathInfo.PassNumber;
+    float3 x = PathInfo.x;
+    float3 w = PathInfo.w;
+
+    float density = VolumeInfo.Extinction.v[PassNumber % 3]; // extinction coefficient multiplier.
+    INFO("  density: %f\n", density);
+
+    float3 bMin, bMax;
+    GetGridBox(VolumeInfo.Grid, bMin, bMax);
+
+    float3 W = { 0 }; // weight
+    W.v[PassNumber % 3] = 3; // Wavelenght dependent importance multiplier
+
+    float tMin, tMax;
+    if (!BoxIntersect(bMin, bMax, x, w, tMin, tMax))
+        return;
+
+    RayPack->start_new_ray();
+    RayPack->add_vertex(x);
+
+
+    float d = tMax - tMin;
+    x += w * tMin;
+
+    RayPack->add_vertex(x);
+
+    while (true) {
+
+        float t =
+            density <= 0.00001
+            ? 10000000
+            : -log(max(0.00000000001, 1.0 - Random::random())) / density; // majorant of all volume data
+
+        INFO("  t: %f\n", t);
+        INFO("  d: %f\n", d);
+
+        if (t >= d) {
+            INFO("->Ray left the volume\n");
+            RayPack->add_vertex(x + (w * t));
+            return;
+        }
+
+        x += w * t; // move to next event position or border
+        RayPack->add_vertex(x);
+
+        float3 tSamplePosition = (x - bMin) / (bMax - bMin);
+        float probExt = Sample3DTexture(VolumeInfo.Grid, tSamplePosition);
+        INFO("  sample pos: %f %f %f\n", tSamplePosition.x, tSamplePosition.y, tSamplePosition.z);
+        INFO("  density there: %f\n", probExt);
+
+        float m_t = probExt * density; // extinction coef
+        float m_s = m_t * VolumeInfo.ScatteringAlbedo.v[PassNumber % 3]; // scattering coef
+        float m_a = m_t - m_s; // absorption coef
+        float m_n = density - m_t; // null coef
+
+        float xi = Random::random();
+
+        float Pa = m_a / density;
+        float Ps = m_s / density;
+        float Pn = m_n / density;
+
+        if (xi < Pa) { // absorption
+            INFO("->absorbtion\n");
+            return;
+        }
+
+        if (xi < 1 - Pn) { // scattering
+            INFO("->scatter\n");
+            w = ImportanceSamplePhase(VolumeInfo.G.v[PassNumber % 3], w); // scattering event...
+
+            if (!BoxIntersect(bMin, bMax, x, w, tMin, tMax))
+                return;
+
+            d = tMax - tMin;
+            x += w * tMin;
+        } else {
+            INFO("->null collision\n");
+            // if no absorption and no scattering null collision occurred
+            d -= t;
+        }
+    }
+}

src/dt_pathtrace.hpp

@@ -0,0 +1,25 @@
+#pragma once
+#include "stdint.h"
+#include "vector.hpp"
+#include "texture3d.hpp"
+#include "ray_pack.hpp"
+
+struct volume_info {
+    texture3D Grid;
+    float3 Extinction;
+    float3 ScatteringAlbedo;
+    float3 G;
+};
+
+struct path_info {
+    float3 x;
+    float3 w;
+    int PassNumber;
+};
+
+void GetGridBox(texture3D grid, float3& minim, float3& maxim);
+void DTPathtrace(path_info PathInfo, volume_info VolumeInfo, ray_pack* RayPack);
+
+namespace Random {
+    void InitRNG(uint32_t seed);
+}

src/list.hpp

@@ -0,0 +1,46 @@
+#pragma once
+#include <cstdlib>
+#include "stdint.h"
+
+template <typename type>
+struct list {
+    uint32_t count;
+    uint32_t length;
+    type* data;
+
+    inline void alloc(uint32_t initial_length = 16) {
+        count = 0;
+        length = initial_length;
+        data = (type*)malloc(sizeof(type) * length);
+    }
+
+    inline void dealloc() {
+        free(data);
+        data = nullptr;
+    }
+
+    inline void clear() {
+        count = 0;
+    }
+
+    void append(type element) {
+        if (count == length) {
+            length *= 2;
+            data = (type*)realloc(data, length * sizeof(type));
+        }
+        data[count] = element;
+        count++;
+    }
+
+    inline type* begin() {
+        return data;
+    }
+
+    inline type* end() {
+        return data+(count);
+    }
+
+     type& operator[](uint32_t index) {
+        return data[index];
+    }
+};

src/main.cpp

@@ -0,0 +1,112 @@
+#include <cstdio>
+#include <cmath>
+
+#include "vector.hpp"
+#include "texture3d.hpp"
+#include "list.hpp"
+#include "defer.hpp"
+#include "dt_pathtrace.hpp"
+#include "ray_pack.hpp"
+
+
+void WriteDomainAndRayPackToOBJFile(texture3D Grid, ray_pack* RayPack,
+                                    const char* OBJFileName)
+{
+    FILE* out = fopen(OBJFileName, "w");
+    if (!out) {
+        fprintf(stderr, "ERROR: file %s could not be opened\n", OBJFileName);
+        return;
+    }
+    defer {
+        fclose(out);
+    };
+
+    for (float3 v : RayPack->vertices) {
+        fprintf(out, "v %f %f %f\n", v.x, v.y, v.z);
+    }
+
+    // Domain Box
+    float3 maxim;
+    float3 minim;
+    GetGridBox(Grid, minim, maxim);
+    fprintf(out, "v %f %f %f\n",  maxim.x,  maxim.y,  maxim.z);
+    fprintf(out, "v %f %f %f\n",  maxim.x,  maxim.y, -maxim.z);
+    fprintf(out, "v %f %f %f\n",  maxim.x, -maxim.y,  maxim.z);
+    fprintf(out, "v %f %f %f\n",  maxim.x, -maxim.y, -maxim.z);
+    fprintf(out, "v %f %f %f\n", -maxim.x,  maxim.y,  maxim.z);
+    fprintf(out, "v %f %f %f\n", -maxim.x,  maxim.y, -maxim.z);
+    fprintf(out, "v %f %f %f\n", -maxim.x, -maxim.y,  maxim.z);
+    fprintf(out, "v %f %f %f\n", -maxim.x, -maxim.y, -maxim.z);
+
+    fprintf(out, "g Rays\n");
+    uint32_t start = 0;
+    uint32_t end;
+    // all but the last sample:
+    uint32_t i = 0;
+    for (; i < RayPack->ray_starts.count-1; ++i) {
+        end = RayPack->ray_starts[i+1];
+        fprintf(out, "o Ray_%u\n", i);
+
+        for (uint32_t j = start; j < end-1; ++j) {
+            // NOTE(Felix): +1 because obj vertex ids start counting at 1
+            fprintf(out, "l %u %u\n", j+1, j+2);
+        }
+        
+        start = end;
+    }
+    // NOTE(Felix): the last sample;
+    end = RayPack->vertices.count;
+    fprintf(out, "o Ray_%u\n", i);
+    printf("start %u end %u\n", start, end);
+    for (uint32_t j = start; j < end-1; ++j) {
+        // NOTE(Felix): +1 because obj vertex ids start counting at 1
+        fprintf(out, "l %u %u\n", j+1, j+2);
+    }
+
+    end++;
+    fprintf(out, "g Domain\n");
+    fprintf(out, "o Domain\n");
+    fprintf(out, "l %u %u\n", end+0, end+1);
+    fprintf(out, "l %u %u\n", end+2, end+3);
+    fprintf(out, "l %u %u\n", end+4, end+5);
+    fprintf(out, "l %u %u\n", end+6, end+7);
+
+    fprintf(out, "l %u %u\n", end+0, end+2);
+    fprintf(out, "l %u %u\n", end+1, end+3);
+    fprintf(out, "l %u %u\n", end+4, end+6);
+    fprintf(out, "l %u %u\n", end+5, end+7);
+
+    fprintf(out, "l %u %u\n", end+0, end+4);
+    fprintf(out, "l %u %u\n", end+1, end+5);
+    fprintf(out, "l %u %u\n", end+2, end+6);
+    fprintf(out, "l %u %u\n", end+3, end+7);
+
+}
+
+int main() {
+    ray_pack RayPack;
+    RayPack.alloc();
+    defer {
+        RayPack.dealloc();
+    };
+
+    path_info PI = { };
+    PI.x = {-2, -2, -2};
+    PI.w = Normalize(float3{0.4,0.4,0} - PI.x);
+    PI.PassNumber = 0;
+
+    volume_info VI = { };
+    VI.Grid       = ReadXYZFile("cloud-049.xyz");
+    VI.Extinction = {20, 20, 20};
+    VI.G          = {0.3, 0.3, 0.3};
+    VI.ScatteringAlbedo = {0.9, 0.9, 0.9};
+
+    Random::InitRNG(1);
+    for (int i = 0; i < 30; ++i) {
+        DTPathtrace(PI, VI, &RayPack);
+    }
+
+    WriteDomainAndRayPackToOBJFile(VI.Grid, &RayPack, "RayPack.obj");
+
+    return 0;
+}

src/ray_pack.hpp

@@ -0,0 +1,27 @@
+#pragma once
+#include "list.hpp"
+#include "vector.hpp"
+
+struct ray_pack {
+    list<float3>   vertices;
+    list<uint32_t> ray_starts;
+
+    void alloc() {
+        vertices.alloc();
+        ray_starts.alloc();
+    }
+
+    void dealloc() {
+        vertices.dealloc();
+        ray_starts.dealloc();
+    }
+
+    void add_vertex(float3 vert) {
+        vertices.append(vert);
+    }
+
+    void start_new_ray() {
+        ray_starts.append(vertices.count);
+    }
+
+};

src/scene.hpp

@@ -0,0 +1,17 @@
+#pragma once
+
+#include "dt_pathtrace.hpp"
+
+struct scene {
+    const char* SceneFile;
+
+    volume_info VolInfo;
+
+    float3 CameraPos;
+    float3 CameraLookAt;
+    float  CameraFOV;
+    uint32_t ResX;
+    uint32_t ResY;
+};
+
+scene LoadSceneFromFile(const char* FilePath);

src/texture3d.cpp

@@ -0,0 +1,47 @@
+#include <cstdio>
+#include <cstdlib>
+#include "texture3d.hpp"
+#include "defer.hpp"
+
+
+texture3D ReadXYZFile(const char* Path) {
+    uint32_t Buffer[9];
+
+    FILE* XYZFile = fopen(Path, "rb");
+
+    if (!XYZFile) {
+        fprintf(stderr, "ERROR: File %s could not be opened\n", Path);
+        return {};
+    }
+    defer {
+        fclose(XYZFile);
+    };
+
+
+    size_t Read = fread(Buffer, sizeof(uint32_t)*9, 1, XYZFile);
+
+    if (!Read) {
+        fprintf(stderr, "ERROR: Header could not be read\n");
+        return {};
+    }
+
+    texture3D Result;
+    Result.w = Buffer[0];
+    Result.h = Buffer[1];
+    Result.d = Buffer[2];
+
+    printf("INFO: Reading density map %u x %u x %u\n",
+           Result.w, Result.h, Result.d);
+
+    uint32_t FloatsCount = Result.w * Result.h * Result.d;
+    Result.Data = (float*)malloc(FloatsCount * sizeof(float));
+
+    Read = fread(Result.Data, FloatsCount * sizeof(float), 1, XYZFile);
+
+    if (!Read) {
+        fprintf(stderr, "ERROR: Data could not be read\n");
+        return {};
+    }
+
+    return Result;
+}

src/texture3d.hpp

@@ -0,0 +1,11 @@
+#pragma once
+#include "stdint.h"
+
+struct texture3D {
+    uint32_t w;
+    uint32_t h;
+    uint32_t d;
+    float* Data;
+};
+
+texture3D ReadXYZFile(const char* Path);

src/vector.cpp

@@ -0,0 +1,144 @@
+#include "vector.hpp"
+#include "math.h"
+
+
+float3 operator/(float3 v, float s) {
+    float3 result;
+    result.x = v.x / s;
+    result.y = v.y / s;
+    result.z = v.z / s;
+    return result;
+}
+
+float3 operator*(float3 v, float s) {
+    float3 result;
+    result.x = v.x * s;
+    result.y = v.y * s;
+    result.z = v.z * s;
+    return result;
+}
+
+float3 operator*(float s, float3 v) {
+    float3 result;
+    result.x = v.x * s;
+    result.y = v.y * s;
+    result.z = v.z * s;
+    return result;
+}
+
+float3 operator+(float3 v1, float3 v2) {
+    float3 result;
+    result.x = v1.x + v2.x;
+    result.y = v1.y + v2.y;
+    result.z = v1.z + v2.z;
+    return result;
+}
+
+float3 operator-(float3 v1) {
+    float3 result;
+    result.x = -v1.x;
+    result.y = -v1.y;
+    result.z = -v1.z;
+    return result;
+}
+
+float3 operator-(float3 v1, float3 v2) {
+    float3 result;
+    result.x = v1.x - v2.x;
+    result.y = v1.y - v2.y;
+    result.z = v1.z - v2.z;
+    return result;
+}
+
+float3 operator/(float3 v1, float3 v2) {
+    float3 result;
+    result.x = v1.x / v2.x;
+    result.y = v1.y / v2.y;
+    result.z = v1.z / v2.z;
+    return result;
+}
+
+float3& operator*=(float3& v, float s) {
+    v = v * s;
+    return v;
+}
+
+float3& operator+=(float3& v1, float3 v2) {
+    v1 = v1 + v2;
+    return v1;
+}
+
+float3& operator-=(float3& v1, float3 v2) {
+    v1 = v1 - v2;
+    return v1;
+}
+
+float3& operator/=(float3& v1, float3 v2) {
+    v1 = v1 / v2;
+    return v1;
+}
+
+float3 Cross(float3 a, float3 b) {
+    float3 result;
+
+    result.x = a.y*b.z - a.z*b.y;
+    result.y = a.z*b.x - a.x*b.z;
+    result.z = a.x*b.y - a.y*b.x;
+
+    return result;
+}
+
+float Dot(float3 a, float3 b) {
+    return
+        a.x * b.x +
+        a.y * b.y +
+        a.z * b.z;
+}
+
+float SquaredLength(float3 v) {
+    return Dot(v, v);
+}
+
+float3 Normalize(float3 v) {
+    float3 result { 0 };
+
+    float squaredLength = SquaredLength(v);
+    if(squaredLength > 0.000001f) {
+        result = v * (1.0f / sqrt(squaredLength));
+    }
+
+    return result;
+}
+
+float3 Abs(float3 v) {
+    float3 result;
+    result.x = fabs(v.x);
+    result.y = fabs(v.y);
+    result.z = fabs(v.z);
+    return result;
+}
+
+float2x3 Abs(float2x3 M) {
+    float2x3 result;
+    result.v1 = Abs(M.v1);
+    result.v2 = Abs(M.v2);
+    return result;
+}
+
+bool operator<=(float3 v, float s) {
+    return
+        v.x <= s &&
+        v.y <= s &&
+        v.z <= s;
+}
+
+bool operator<=(float2x3 M, float s) {
+    return M.v1 <= s && M.v2 <= s;
+}
+
+float2x3 operator/(float2x3 M1, float2x3 M2) {
+    float2x3 result;
+    result.v1 = M1.v1 / M2.v1;
+    result.v2 = M1.v2 / M2.v2;
+    return result;
+}

src/vector.hpp

@@ -0,0 +1,40 @@
+#pragma once
+
+typedef union {
+    struct {
+        float x, y, z;
+    };
+    struct {
+        float r, g, b;
+    };
+    float v[3];
+} float3;
+
+typedef struct {
+    float3 v1;
+    float3 v2;
+} float2x3;
+
+
+float3 Cross(float3, float3);
+float3 Normalize(float3);
+
+
+float3   Abs(float3);
+float2x3 Abs(float2x3);
+
+float2x3 operator/(float2x3, float2x3);
+bool operator<=(float2x3, float);
+
+float3 operator-(float3);
+float3 operator*(float3, float);
+float3 operator*(float, float3);
+float3 operator/(float3, float);
+float3 operator+(float3, float3);
+float3 operator-(float3, float3);
+float3 operator/(float3, float3);
+
+float3& operator*=(float3&, float);
+float3& operator/=(float3&, float3);
+float3& operator+=(float3&, float3);
+float3& operator-=(float3&, float3);