commit 5a0d81134ed88c2eebe1c64ad7a037b98892819c
Author: Yggdrasil75 <cblackburn7557@gmail.com>
Date:   Wed Nov 5 11:56:56 2025 -0500

    init

diff --git a/c.cpp b/c.cpp
new file mode 100644
index 0000000..26c9049
--- /dev/null
+++ b/c.cpp
@@ -0,0 +1,937 @@
+#include <cmath>
+#include <iostream>
+#include <vector>
+#include <random>
+#include <algorithm>
+#include <unordered_map>
+#include <functional>
+#include "timing_decorator.hpp"
+#include <cstdint>
+#include <fstream>
+#include <string>
+
+const float EPSILON = 0.00000000001;
+
+//classes and structs
+
+struct Bool3 {
+    bool x, y, z;
+};
+
+class Vec3 {
+public:
+    double x, y, z;
+    
+    Vec3(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
+    
+    //Vec3(const VoxelIndex& idx) : x(static_cast<double>(idx.x)), y(static_cast<double>(idx.y)), z(static_cast<double>(idx.z)) {}
+    
+    inline double norm() const {
+        return std::sqrt(x*x + y*y + z*z);
+    }
+    
+    inline Vec3 normalize() const {
+        double n = norm();
+        return Vec3(x/n, y/n, z/n);
+    }
+    
+    inline Vec3 cross(const Vec3& other) const {
+        return Vec3(
+            y * other.z - z * other.y,
+            z * other.x - x * other.z,
+            x * other.y - y * other.x
+        );
+    }
+    
+    inline double dot(const Vec3& other) const {
+        return x * other.x + y * other.y + z * other.z;
+    }
+    
+    inline Vec3 operator+(float scalar) const {
+        return Vec3(x + scalar, y + scalar, z + scalar);
+    }
+
+    inline Vec3 operator+(const Vec3& other) const {
+        return Vec3(x + other.x, y + other.y, z + other.z);
+    }
+
+    Vec3 operator+(const Bool3& other) const {
+        return {
+            x + (other.x ? 1.0f : 0.0f),
+            y + (other.y ? 1.0f : 0.0f),
+            z + (other.z ? 1.0f : 0.0f)
+        };
+    }
+    
+    inline Vec3 operator-(const Vec3& other) const {
+        return Vec3(x - other.x, y - other.y, z - other.z);
+    }
+    
+    inline Vec3 operator*(float scalar) const {
+        return Vec3(x * scalar, y * scalar, z * scalar);
+    }
+    
+    inline Vec3 operator*(const Vec3& scalar) const {
+        return Vec3(x * scalar.x, y * scalar.y, z * scalar.z);
+    }
+
+    inline friend Vec3 operator*(float scalar, const Vec3& vec) {
+        return Vec3(scalar * vec.x, scalar * vec.y, scalar * vec.z);
+    }
+
+    inline Vec3 operator/(float scalar) const {
+        return Vec3(x / scalar, y / scalar, z / scalar);
+    }
+    
+    inline Vec3 operator/(const Vec3& scalar) const {
+        return Vec3(x / scalar.x, y / scalar.y, z / scalar.z);
+    }
+
+    inline friend Vec3 operator/(float scalar, const Vec3& vec) {
+        return Vec3(vec.x / scalar, vec.y / scalar, vec.z / scalar);
+    }
+
+    bool operator==(const Vec3& other) const {
+        return x == other.x && y == other.y && z == other.z;
+    }
+    
+    bool operator<(const Vec3& other) const {
+        if (x != other.x) return x < other.x;
+        if (y != other.y) return y < other.y;
+        return z < other.z;
+    }
+
+    const float& operator[](int index) const {
+        if (index == 0) return x;
+        if (index == 1) return y;
+        if (index == 2) return z;
+        throw std::out_of_range("Index out of range");
+    }
+
+    double& operator[](int index) {
+        if (index == 0) return x;
+        if (index == 1) return y;
+        if (index == 2) return z;
+        throw std::out_of_range("Index out of range");
+    }
+    
+    struct Hash {
+        size_t operator()(const Vec3& v) const {
+            size_t h1 = std::hash<double>()(std::round(v.x * 1000.0));
+            size_t h2 = std::hash<double>()(std::round(v.y * 1000.0));
+            size_t h3 = std::hash<double>()(std::round(v.z * 1000.0));
+            return h1 ^ (h2 << 1) ^ (h3 << 2);
+        }
+    };
+
+    std::string toString() const {
+        return "Vec3(" + std::to_string(x) + ", " + 
+                             std::to_string(y) + ", " + 
+                             std::to_string(z) + ")";
+    }
+
+    Vec3& safe_inverse_dir(float epsilon = 1e-6f) {
+        x = (std::abs(x) > epsilon) ? x : std::copysign(epsilon, -x);
+        y = (std::abs(y) > epsilon) ? y : std::copysign(epsilon, -y);
+        z = (std::abs(z) > epsilon) ? z : std::copysign(epsilon, -z);
+        return *this;
+    }
+
+    Vec3 sign() const {
+        return Vec3(
+            (x > 0) ? 1 : ((x < 0) ? -1 : 0),
+            (y > 0) ? 1 : ((y < 0) ? -1 : 0),
+            (z > 0) ? 1 : ((z < 0) ? -1 : 0)
+        );
+    }
+
+    Vec3 abs() {
+        return Vec3(std::abs(x), std::abs(y), std::abs(z));
+    }
+};
+
+class Vec4 {
+public:
+    double x, y, z, w;
+    
+    Vec4(double x = 0, double y = 0, double z = 0, double w = 0) : x(x), y(y), z(z), w(w) {}
+
+    inline double norm() const {
+        return std::sqrt(x*x + y*y + z*z + w*w);
+    }
+    
+    inline Vec4 normalize() const {
+        double n = norm();
+        return Vec4(x/n, y/n, z/n, w/n);
+    }
+    
+    inline std::array<Vec4, 6> wedge(const Vec4& other) const {
+        return {
+            Vec4(0, x*other.y - y*other.x, 0, 0),           // xy-plane
+            Vec4(0, 0, x*other.z - z*other.x, 0),           // xz-plane  
+            Vec4(0, 0, 0, x*other.w - w*other.x),           // xw-plane
+            Vec4(0, 0, y*other.z - z*other.y, 0),           // yz-plane
+            Vec4(0, 0, 0, y*other.w - w*other.y),           // yw-plane
+            Vec4(0, 0, 0, z*other.w - w*other.z)            // zw-plane
+        };
+    }
+    inline double dot(const Vec4& other) const {
+        return x * other.x + y * other.y + z * other.z + w * other.w;
+    }
+    
+    inline Vec4 operator+(const Vec4& other) const {
+        return Vec4(x + other.x, y + other.y, z + other.z, w + other.w);
+    }
+    
+    inline Vec4 operator-(const Vec4& other) const {
+        return Vec4(x - other.x, y - other.y, z - other.z, w - other.w);
+    }
+    
+    inline Vec4 operator*(double scalar) const {
+        return Vec4(x * scalar, y * scalar, z * scalar, w * scalar);
+    }
+    
+    inline Vec4 operator/(double scalar) const {
+        return Vec4(x / scalar, y / scalar, z / scalar, w / scalar);
+    }
+    
+    bool operator==(const Vec4& other) const {
+        return x == other.x && y == other.y && z == other.z && w == other.w;
+    }
+    
+    bool operator<(const Vec4& other) const {
+        if (x != other.x) return x < other.x;
+        if (y != other.y) return y < other.y;
+        if (z != other.z) return z < other.z;
+        return w < other.w;
+    }
+    
+    // Additional useful methods for 4D vectors
+    inline Vec4 homogenize() const {
+        if (w == 0) return *this;
+        return Vec4(x/w, y/w, z/w, 1.0);
+    }
+    
+    inline Vec3 xyz() const {
+        return Vec3(x, y, z);
+    }
+    
+    struct Hash {
+        size_t operator()(const Vec4& v) const {
+            size_t h1 = std::hash<double>()(std::round(v.x * 1000.0));
+            size_t h2 = std::hash<double>()(std::round(v.y * 1000.0));
+            size_t h3 = std::hash<double>()(std::round(v.z * 1000.0));
+            size_t h4 = std::hash<double>()(std::round(v.w * 1000.0));
+            return h1 ^ (h2 << 1) ^ (h3 << 2) ^ (h4 << 3);
+        }
+    };
+};
+
+struct VoxelIndex {
+    int x, y, z;
+
+    // Constructor
+    VoxelIndex(int x = 0, int y = 0, int z = 0) : x(x), y(y), z(z) {}
+
+    // Array-like access
+    int& operator[](size_t index) {
+        switch(index) {
+            case 0: return x;
+            case 1: return y;
+            case 2: return z;
+            default: throw std::out_of_range("Index out of range");
+        }
+    }
+
+    const int& operator[](size_t index) const {
+        switch(index) {
+            case 0: return x;
+            case 1: return y;
+            case 2: return z;
+            default: throw std::out_of_range("Index out of range");
+        }
+    }
+
+    inline VoxelIndex operator+(const Vec3& other) const {
+        return VoxelIndex(std::floor(x + other.x), std::floor(y + other.y), std::floor(z + other.z));
+    }
+
+    inline VoxelIndex operator-(const Vec3& other) const {
+        return VoxelIndex(std::floor(x - other.x), std::floor(y - other.y), std::floor(z - other.z));
+    }
+
+    inline VoxelIndex operator+(float scalar) const {
+        return VoxelIndex(std::floor(x + scalar), std::floor(y + scalar), std::floor(z + scalar));
+    }
+    
+    VoxelIndex operator+(const Bool3& other) const {
+        return {
+            x + (other.x ? 1 : 0),
+            y + (other.y ? 1 : 0),
+            z + (other.z ? 1 : 0)
+        };
+    }
+
+    inline VoxelIndex operator*(const Vec3& other) const {
+        return VoxelIndex(std::floor(x * other.x), std::floor(y * other.y), std::floor(z * other.z));
+    }
+    inline VoxelIndex operator*(float scalar) const {
+        return VoxelIndex(std::floor(x * scalar), std::floor(y * scalar), std::floor(z * scalar));
+    }
+
+    operator Vec3() const {
+        return Vec3(static_cast<double>(x), static_cast<double>(y), static_cast<double>(z));
+    }
+
+    Vec3 toVec3() const {
+        return Vec3(static_cast<double>(x), static_cast<double>(y), static_cast<double>(z));
+    }
+
+    // Hash function
+    size_t hash() const {
+        return std::hash<int>{}(x) ^ 
+               (std::hash<int>{}(y) << 1) ^ 
+               (std::hash<int>{}(z) << 2);
+    }
+
+    // Convert to string
+    std::string toString() const {
+        return "VoxelIndex(" + std::to_string(x) + ", " + 
+                             std::to_string(y) + ", " + 
+                             std::to_string(z) + ")";
+    }
+
+    // Comparison operators for completeness
+    bool operator==(const VoxelIndex& other) const {
+        return x == other.x && y == other.y && z == other.z;
+    }
+
+    bool operator!=(const VoxelIndex& other) const {
+        return !(*this == other);
+    }
+};
+
+namespace std {
+    template<>
+    struct hash<Vec3> {
+        size_t operator()(const Vec3& p) const {
+            return hash<float>()(p.x) ^ hash<float>()(p.y) ^ hash<float>()(p.z);
+        }
+    };
+    template<>
+    struct hash<VoxelIndex> {
+        size_t operator()(const VoxelIndex& idx) const {
+            return idx.hash();
+        }
+    };
+}
+
+class VoxelGrid {
+public:
+    // New structure - using position-based indexing
+    std::unordered_map<Vec3, size_t> pos_index_map;  // Maps voxel center -> index in positions/colors
+    std::vector<Vec3> positions;                     // Voxel center positions
+    std::vector<Vec4> colors;                        // Average colors per voxel
+    float gridsize;                                  // Voxel size
+    
+    // Old structure - to be removed/replaced
+    float voxel_size;
+    Vec3 min_bounds;
+    Vec3 max_bounds;
+        
+    // Grid dimensions
+    int dim_x, dim_y, dim_z;
+    
+    // Grid arrays for fast access (like Python)
+    std::vector<bool> grid_array;
+    std::vector<Vec4> color_array;
+    std::vector<int> count_array;
+    
+    bool arrays_initialized = false;
+
+    void initializeArrays() {
+        if (arrays_initialized) return;
+        
+        size_t total_size = dim_x * dim_y * dim_z;
+        grid_array.resize(total_size, false);
+        color_array.resize(total_size, Vec4{0, 0, 0, 0});
+        count_array.resize(total_size, 0);
+        
+        // Populate arrays from the new data structure
+        for (const auto& [voxel_pos, index] : pos_index_map) {
+            VoxelIndex voxel_idx = getVoxelIndex(voxel_pos);
+            if (isValidIndex(voxel_idx)) {
+                size_t array_idx = getArrayIndex(voxel_idx);
+                grid_array[array_idx] = true;
+                color_array[array_idx] = colors[index];
+                count_array[array_idx] = 1; // Each entry represents one voxel
+            }
+        }
+        
+        arrays_initialized = true;
+    }
+
+    size_t getArrayIndex(const VoxelIndex& index) const {
+        return index[0] + dim_x * (index[1] + dim_y * index[2]);
+    }
+
+    bool isValidIndex(const VoxelIndex& idx) const {
+        return idx[0] >= 0 && idx[0] < dim_x &&
+               idx[1] >= 0 && idx[1] < dim_y &&
+               idx[2] >= 0 && idx[2] < dim_z;
+    }
+
+    VoxelGrid(float size = 0.1f) : gridsize(size), voxel_size(size), arrays_initialized(false) {}
+    
+    void addPoints(const std::vector<Vec3>& points, const std::vector<Vec4>& input_colors) {
+        if (points.size() != input_colors.size()) {
+            std::cerr << "Error: Points and colors vectors must have the same size" << std::endl;
+            return;
+        }
+        
+        if (points.empty()) return;
+        
+        // Calculate bounds
+        min_bounds = points[0];
+        max_bounds = points[0];
+        
+        for (const auto& point : points) {
+            min_bounds.x = std::min(min_bounds.x, point.x);
+            min_bounds.y = std::min(min_bounds.y, point.y);
+            min_bounds.z = std::min(min_bounds.z, point.z);
+            
+            max_bounds.x = std::max(max_bounds.x, point.x);
+            max_bounds.y = std::max(max_bounds.y, point.y);
+            max_bounds.z = std::max(max_bounds.z, point.z);
+        }
+        
+        // Calculate grid dimensions
+        Vec3 range = max_bounds - min_bounds;
+        dim_x = static_cast<int>(std::ceil(range.x / gridsize)) + 1;
+        dim_y = static_cast<int>(std::ceil(range.y / gridsize)) + 1;
+        dim_z = static_cast<int>(std::ceil(range.z / gridsize)) + 1;
+        
+        // Temporary storage for averaging colors per voxel
+        std::unordered_map<Vec3, std::vector<Vec4>, Vec3::Hash> voxel_color_map;
+        
+        // Group points by voxel and collect colors
+        for (size_t i = 0; i < points.size(); ++i) {
+            Vec3 voxel_center = getVoxelCenter(getVoxelIndex(points[i]));
+            voxel_color_map[voxel_center].push_back(input_colors[i]);
+        }
+        
+        // Convert to the new structure - store average colors per voxel
+        positions.clear();
+        colors.clear();
+        pos_index_map.clear();
+        
+        for (auto& [voxel_pos, color_list] : voxel_color_map) {
+            // Calculate average color
+            Vec4 avg_color{0, 0, 0, 0};
+            for (const auto& color : color_list) {
+                avg_color.w += color.w;
+                avg_color.x += color.x;
+                avg_color.y += color.y;
+                avg_color.z += color.z;
+            }
+            avg_color.w /= color_list.size();
+            avg_color.x /= color_list.size();
+            avg_color.y /= color_list.size();
+            avg_color.z /= color_list.size();
+            
+            // Store in new structure
+            size_t index = positions.size();
+            positions.push_back(voxel_pos);
+            colors.push_back(avg_color);
+            pos_index_map[voxel_pos] = index;
+        }
+        
+        // Initialize arrays for fast access
+        initializeArrays();
+    }
+    
+    // Get voxel data using voxel center position
+    bool hasVoxelAt(const Vec3& voxel_center) const {
+        return pos_index_map.find(voxel_center) != pos_index_map.end();
+    }
+    
+    Vec4 getVoxelColor(const Vec3& voxel_center) const {
+        auto it = pos_index_map.find(voxel_center);
+        if (it != pos_index_map.end()) {
+            return colors[it->second];
+        }
+        return Vec4{0, 0, 0, 0};
+    }
+    
+    // Fast array access using VoxelIndex (compatibility with existing code)
+    bool getVoxelOccupied(const VoxelIndex& voxel_idx) const {
+        if (!arrays_initialized || !isValidIndex(voxel_idx)) return false;
+        return grid_array[getArrayIndex(voxel_idx)];
+    }
+    
+    Vec4 getVoxelColor(const VoxelIndex& voxel_idx) const {
+        if (!arrays_initialized || !isValidIndex(voxel_idx)) 
+            return Vec4{0, 0, 0, 0};
+        return color_array[getArrayIndex(voxel_idx)];
+    }
+    
+    int getVoxelPointCount(const VoxelIndex& voxel_idx) const {
+        if (!arrays_initialized || !isValidIndex(voxel_idx)) return 0;
+        return count_array[getArrayIndex(voxel_idx)];
+    }
+    
+    std::vector<VoxelIndex> getVoxelIndices() const {
+        std::vector<VoxelIndex> indices;
+        for (const auto& voxel_pos : positions) {
+            indices.push_back(getVoxelIndex(voxel_pos));
+        }
+        return indices;
+    }
+    
+    std::vector<Vec3> getVoxelCenters() const {
+        return positions;
+    }
+    
+    size_t getNumVoxels() const { return positions.size(); }
+    size_t getNumPoints() const {
+        // Since we're storing averaged voxels, each position represents multiple original points
+        // For simplicity, return number of voxels
+        return positions.size();
+    }
+    
+    Vec3 getMinBounds() const { return min_bounds; }
+    Vec3 getMaxBounds() const { return max_bounds; }
+    std::array<int, 3> getDimensions() const { return {dim_x, dim_y, dim_z}; }
+    float getVoxelSize() const { return gridsize; }
+    
+    VoxelIndex getVoxelIndex(const Vec3& point) const {
+        Vec3 normalized = point - min_bounds;
+        return VoxelIndex{
+            static_cast<int>(std::floor(normalized.x / gridsize)),
+            static_cast<int>(std::floor(normalized.y / gridsize)),
+            static_cast<int>(std::floor(normalized.z / gridsize))
+        };
+    }
+
+    Vec3 getVoxelCenter(const VoxelIndex& voxel_idx) const {
+        return Vec3(
+            min_bounds.x + (voxel_idx[0] + 0.5f) * gridsize,
+            min_bounds.y + (voxel_idx[1] + 0.5f) * gridsize,
+            min_bounds.z + (voxel_idx[2] + 0.5f) * gridsize
+        );
+    }
+
+    void clear() {
+        pos_index_map.clear();
+        positions.clear();
+        colors.clear();
+        grid_array.clear();
+        color_array.clear();
+        count_array.clear();
+        arrays_initialized = false;
+    }
+    
+    void setVoxelSize(float size) {
+        gridsize = size;
+        voxel_size = size;
+        clear();
+    }
+    
+    void printStats() const {
+        std::cout << "Voxel Grid Statistics:" << std::endl;
+        std::cout << "  Voxel size: " << gridsize << std::endl;
+        std::cout << "  Grid dimensions: " << dim_x << " x " << dim_y << " x " << dim_z << std::endl;
+        std::cout << "  Number of voxels: " << getNumVoxels() << std::endl;
+        std::cout << "  Using new position-based storage" << std::endl;
+    }
+};
+
+struct Image {
+    int width;
+    int height;
+    std::vector<uint8_t> data; // RGBA format
+    
+    Image(int w, int h) : width(w), height(h), data(w * h * 4) {
+        // Initialize to white background
+        for (int i = 0; i < w * h * 4; i += 4) {
+            data[i] = 255;     // R
+            data[i + 1] = 255; // G  
+            data[i + 2] = 255; // B
+            data[i + 3] = 255; // A
+        }
+    }    
+
+    // Helper methods
+    uint8_t* pixel(int x, int y) {
+        return &data[(y * width + x) * 4];
+    }
+    
+    void setPixel(int x, int y, uint8_t r, uint8_t g, uint8_t b, uint8_t a = 255) {
+        uint8_t* p = pixel(x, y);
+        p[0] = r; p[1] = g; p[2] = b; p[3] = a;
+    }
+
+    bool saveAsBMP(const std::string& filename) {
+        #pragma pack(push, 1)
+
+        // BMP file header (14 bytes)
+        struct BMPFileHeader {
+            uint16_t file_type{0x4D42}; // "BM"
+            uint32_t file_size{0};
+            uint16_t reserved1{0};
+            uint16_t reserved2{0};
+            uint32_t offset_data{0};
+        } file_header;
+        
+        // BMP info header (40 bytes)
+        struct BMPInfoHeader {
+            uint32_t size{40};
+            int32_t width{0};
+            int32_t height{0};
+            uint16_t planes{1};
+            uint16_t bit_count{32};
+            uint32_t compression{0}; // BI_RGB
+            uint32_t size_image{0};
+            int32_t x_pixels_per_meter{0};
+            int32_t y_pixels_per_meter{0};
+            uint32_t colors_used{0};
+            uint32_t colors_important{0};
+        } info_header;
+        
+        // Calculate sizes
+        uint32_t row_stride = width * 4;
+        uint32_t padding_size = (4 - (row_stride % 4)) % 4;
+        uint32_t data_size = (row_stride + padding_size) * height;
+        
+        file_header.file_size = sizeof(BMPFileHeader) + sizeof(BMPInfoHeader) + data_size;
+        file_header.offset_data = sizeof(BMPFileHeader) + sizeof(BMPInfoHeader);
+        
+        info_header.width = width;
+        info_header.height = -height; // Negative for top-down bitmap (no flipping needed)
+        info_header.size_image = data_size;
+        
+        // Open file
+        std::ofstream file(filename, std::ios::binary);
+        if (!file.is_open()) {
+            return false;
+        }
+        
+        // Write headers
+        file.write(reinterpret_cast<const char*>(&file_header), sizeof(file_header));
+        file.write(reinterpret_cast<const char*>(&info_header), sizeof(info_header));
+        
+        // Write pixel data (BMP stores as BGR, we have RGB)
+        std::vector<uint8_t> row_buffer(row_stride + padding_size);
+        
+        for (int y = 0; y < height; ++y) {
+            uint8_t* row_start = &data[y * row_stride];
+            
+            // Convert RGB to BGR and copy to buffer
+            for (int x = 0; x < width; ++x) {
+                uint8_t* src_pixel = row_start + (x * 4);
+                uint8_t* dst_pixel = row_buffer.data() + (x * 4);
+                
+                // Swap R and B channels (RGBA -> BGRA)
+                dst_pixel[0] = src_pixel[2]; // B
+                dst_pixel[1] = src_pixel[1]; // G
+                dst_pixel[2] = src_pixel[0]; // R
+                //dst_pixel[3] = src_pixel[3]; // A
+            }
+            
+            // Write row with padding
+            file.write(reinterpret_cast<const char*>(row_buffer.data()), row_stride + padding_size);
+        }
+        
+        return file.good();
+    }
+};
+
+
+//noise functions
+float fade(const float& a) {
+    TIME_FUNCTION;
+    return a * a * a * (10 + a * (-15 + a * 6));
+}
+
+float clamp(float x, float lowerlimit = 0.0f, float upperlimit = 1.0f) {
+    TIME_FUNCTION;
+    if (x < lowerlimit) return lowerlimit;
+    if (x > upperlimit) return upperlimit;
+    return x;
+}
+
+float pascalTri(const float& a, const float& b) {
+    TIME_FUNCTION;
+    int result = 1;
+    for (int i = 0; i < b; ++i){
+        result *= (a - 1) / (i + 1);
+    }
+    return result;
+}
+
+float genSmooth(int N, float x) {
+    TIME_FUNCTION;
+    x = clamp(x, 0, 1);
+    float result = 0;
+    for (int n = 0; n <= N; ++n){
+        result += pascalTri(-N - 1, n) * pascalTri(2 * N + 1, N-1) * pow(x, N + n + 1);
+    }
+    return result;
+}
+
+float inverse_smoothstep(float x) {
+    TIME_FUNCTION;
+    return 0.5 - sin(asin(1.0 - 2.0 * x) / 3.0);
+}
+
+float lerp(const float& t, const float& a, const float& b) {
+    TIME_FUNCTION;
+    return a + t * (b - a);
+}
+
+float grad(const int& hash, const float& b, const float& c, const float& d) {
+    TIME_FUNCTION;
+    int h = hash & 15;
+    float u = (h < 8) ? c : b;
+    float v = (h < 4) ? b : ((h == 12 || h == 14) ? c : d);
+    return (((h & 1) == 0) ? u : -u) + (((h & 2) == 0) ? v : -v);
+}
+
+float pnoise3d(const int p[512], const float& xf, const float& yf, const float& zf) {
+    TIME_FUNCTION;
+    int floorx = std::floor(xf);
+    int floory = std::floor(yf);
+    int floorz = std::floor(zf);
+    int iX = floorx & 255;
+    int iY = floory & 255;
+    int iZ = floorz & 255;
+
+    int x = xf - floorx;
+    int y = yf - floory;
+    int z = zf - floorz;
+
+    float u = fade(x);
+    float v = fade(y);
+    float w = fade(z);
+
+    int A  = p[iX] + iY;
+    int AA = p[A] + iZ;
+    int AB = p[A+1] + iZ;
+
+    int B  = p[iX + 1] + iY;
+    int BA = p[B] + iZ;
+    int BB = p[B+1] + iZ;
+
+    float f = grad(p[BA], x-1, y, z);
+    float g = grad(p[AA], x, y, z);
+    float h = grad(p[BB], x-1, y-1, z);
+    float j = grad(p[BB], x-1, y-1, z);
+    float k = grad(p[AA+1], x, y, z-1);
+    float l = grad(p[BA+1], x-1, y, z-1);
+    float m = grad(p[AB+1], x, y-1, z-1);
+    float n = grad(p[BB+1], x-1, y-1, z-1);
+
+    float o = lerp(u, m, n);
+    float q = lerp(u, k, l);
+    float r = lerp(u, h, j);
+    float s = lerp(u, f, g);
+    float t = lerp(v, q, o);
+    float e = lerp(v, s, r);
+    float d = lerp(w, e, t);
+    return d;
+}
+
+std::tuple<std::vector<Vec3>, std::vector<Vec4>> noiseBatch(int num_points, float scale, int sp[]) {
+    TIME_FUNCTION;
+    std::vector<Vec3> points;
+    std::vector<Vec4> colors;
+    points.reserve(num_points);
+    colors.reserve(num_points);
+
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<float> dis(-scale, scale);
+    
+    for (int i = 0; i < num_points; ++i) {
+        float x = dis(gen);
+        float y = dis(gen);
+        float z = dis(gen);
+        
+        float noise1 = pnoise3d(sp, (x * 0.5f), (y * 0.5f), (z * 0.5f));
+        float noise2 = pnoise3d(sp, (x * 0.3f), (y * 0.3f), (z * 0.3f));
+        float noise3 = pnoise3d(sp, (x * 0.7f), (y * 0.7f), (z * 0.7f));
+        float noise4 = pnoise3d(sp, (x * 0.7f), (y * 0.7f), (z * 0.7f));
+        
+        if (noise1 > 0.1f) {
+            float rt = (noise1 + 1.0f) * 0.5f;
+            float gt = (noise2 + 1.0f) * 0.5f;
+            float bt = (noise3 + 1.0f) * 0.5f;
+            float at = (noise4 + 1.0f) * 0.5f;
+
+            float maxV = std::max({rt, gt, bt});
+            if (maxV > 0) {
+                float r = rt / maxV;
+                float g = gt / maxV;
+                float b = bt / maxV;
+                float a = at / maxV;
+                points.push_back({x, y, z});
+                colors.push_back({r, g, b, a});
+            }
+        }
+    }
+    return std::make_tuple(points, colors);
+}
+
+std::tuple<std::vector<Vec3>, std::vector<Vec4>> genPointCloud(int numP, float scale, int seed) {
+    TIME_FUNCTION;
+    int permutation[256];
+    for (int i = 0; i < 256; ++i) {
+        permutation[i] = i;
+    }
+    std::mt19937 rng(seed);
+    std::shuffle(permutation, permutation+256, rng);
+    int p[512];
+    for (int i = 0; i < 256; ++i) {
+        p[i] = permutation[i];
+        p[i + 256] = permutation[i];
+    }
+    return noiseBatch(numP, scale, p);
+}
+
+
+//voxel stuff
+Bool3 greaterThanZero(const Vec3& v) {
+    return {v.x > 0, v.y > 0, v.z > 0};
+}
+
+Image render(int height, int width, Vec3 forward, Vec3 right, Vec3 up, Vec3 rayOrigin, Vec3 vbound,
+            float vsize, VoxelGrid grid, int dims) {
+    TIME_FUNCTION;
+
+    Image img = Image(height, width);
+    
+    float max_t = 50.0f;
+    int max_steps = 123;
+    float max_dist = 25.0f;
+    float screen_height = height;
+    float screen_width = width;
+    
+    float inv_w = 1.0f / width;
+    float inv_h = 1.0f / height;
+    float scr_w_half = screen_width * 0.5f;
+    float scr_h_half = screen_height * 0.5f;
+    
+    std::cout << height << std::endl;
+    std::cout << width << std::endl;
+    for (int y = 0; y < height; y++) {
+        int sy = std::ceil(1.0 - ((2.0 * y) * inv_h) * scr_h_half);
+        for (int x = 0; x < width; x++) {
+            int sx = std::ceil((((2.0 * x) * inv_w) - 1.0) * scr_w_half);
+            std::cout << "working in: " << x << ", " << y << std::endl;
+            Vec3 ray_dir = (forward + (sx * right) + (sy * up)).normalize();
+            std::cout << "current ray direction: " << ray_dir.toString() << std::endl;
+            Vec3 cv1 = ((rayOrigin - vbound) / vsize);
+            VoxelIndex cv = VoxelIndex(static_cast<int>(cv1.x), static_cast<int>(cv1.y), static_cast<int>(cv1.z));
+            std::cout << "cell at: " << cv.toString() << std::endl;
+            Vec3 inv_dir =  Vec3(
+                ray_dir.x != 0 ? 1.0f / ray_dir.x : std::numeric_limits<float>::max(),
+                ray_dir.y != 0 ? 1.0f / ray_dir.y : std::numeric_limits<float>::max(),
+                ray_dir.z != 0 ? 1.0f / ray_dir.z : std::numeric_limits<float>::max()
+            );
+            std::cout << "inverse of the current ray: " << inv_dir.toString() << std::endl;
+            Vec3 step = ray_dir.sign();
+            std::cout << "current ray signs: " << step.toString() << std::endl;
+            Bool3 step_mask = greaterThanZero(step);
+            VoxelIndex next_voxel_bound = (cv + step_mask) * vsize + vbound;
+            std::cout << "next cell at: " << next_voxel_bound.toString() << std::endl;
+            Vec3 t_max = (Vec3(next_voxel_bound) - rayOrigin) * inv_dir;
+            std::cout << "t_max at: " << t_max.toString() << std::endl;
+            Vec3 t_delta = vsize / inv_dir.abs();
+            std::cout << "t_delta: " << t_delta.toString() << std::endl;
+            float t = 0.0f;
+
+            Vec4 accumulatedColor = Vec4(0,0,0,0);
+            //std::cout << x << "," << y << std::endl;
+            for (int iter = 0; iter < max_steps; iter++) {
+                if (max_t < t) {
+                    std::cout << "t failed " << std::endl;
+                    break;
+                }
+                if (accumulatedColor.z >= 1.0) {
+                    std::cout << "z failed " << std::endl;
+                    break;
+                }
+
+                if (cv.x >= 0 && cv.x < dims &&
+                    cv.y >= 0 && cv.y < dims && 
+                    cv.z >= 0 && cv.z < dims) {
+                    std::cout << "found cell at: " << cv.toString() << std::endl;
+                    if (grid.getVoxelOccupied(cv)) {
+                        std::cout << "found occupied cell at: " << cv.toString() << std::endl;
+                        Vec4 voxel_color = grid.getVoxelColor(cv);
+                        
+                        float weight = voxel_color.z * (1.0f - accumulatedColor.z);
+                        accumulatedColor.w += voxel_color.w * weight;
+                        accumulatedColor.x += voxel_color.x * weight;
+                        accumulatedColor.y += voxel_color.y * weight;
+                        accumulatedColor.z += voxel_color.z * weight;
+                    }
+
+                    int minAxis = 0;
+                    if (t_max.y < t_max.x) minAxis = 1;
+                    if (t_max.z < t_max[minAxis]) minAxis = 2;
+                    cv[minAxis] += step[minAxis];
+                    t = t_max[minAxis];
+                    t_max[minAxis] += t_delta[minAxis];
+                }
+            }
+            if (accumulatedColor.z > 0) {
+                std::cout << "setting a color at " << x << " and " << y << std::endl;
+                float r = accumulatedColor.w + (1.0f - accumulatedColor.z) * 1.0f;
+                float g = accumulatedColor.x + (1.0f - accumulatedColor.z) * 1.0f;
+                float b = accumulatedColor.y + (1.0f - accumulatedColor.z) * 1.0f;
+                
+                img.setPixel(x, y, 
+                    static_cast<uint8_t>(r * 255),
+                    static_cast<uint8_t>(g * 255), 
+                    static_cast<uint8_t>(b * 255),
+                    255);
+            }
+            
+        }
+    }
+    return img;
+}
+
+int main() {
+    std::cout << "Generating point cloud" << std::endl;
+    auto [points, colors] = genPointCloud(150000, 10.0, 43);
+    std::cout << "Generating done" << std::endl;
+    
+    
+    VoxelGrid voxel_grid(0.2f);
+    
+    std::cout << "Adding points to voxel grid..." << std::endl;
+    voxel_grid.addPoints(points, colors);
+    voxel_grid.printStats();
+    
+    // Use the voxel grid's actual bounds
+    Vec3 min_bounds = voxel_grid.getMinBounds();
+    Vec3 max_bounds = voxel_grid.getMaxBounds();
+    Vec3 grid_center = (min_bounds + max_bounds) * 0.5;
+
+    // Proper camera setup
+    Vec3 rayOrigin(0, 0, 15);
+    Vec3 forward = (grid_center - rayOrigin).normalize();
+    Vec3 up(0, 1, 0);
+    Vec3 right = forward.cross(up).normalize();
+    auto dims = voxel_grid.getDimensions();
+    int max_dim = std::max({dims[0], dims[1], dims[2]});
+
+    Image img = render(50, 50, forward, right, up, rayOrigin, min_bounds, 
+                    voxel_grid.getVoxelSize(), voxel_grid, max_dim);
+    img.saveAsBMP("cpp_voxel_render.bmp");
+    
+    FunctionTimer::printStats(FunctionTimer::Mode::ENHANCED);
+    return 0;
+}
diff --git a/d.cpp b/d.cpp
new file mode 100644
index 0000000..11382a3
--- /dev/null
+++ b/d.cpp
@@ -0,0 +1,659 @@
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <algorithm>
+#include <unordered_map>
+#include <cstdio>
+#include <cstring>
+#include <fstream>
+#include <iostream>
+#include <vector>
+#include <random>
+#include <functional>
+#include "timing_decorator.hpp"
+#include <cstdint>
+#include <string>
+
+
+// Forward declarations
+class Vec3;
+class Vec4;
+
+class Vec3 {
+public:
+    float x, y, z;
+    
+    // Constructors
+    Vec3() : x(0), y(0), z(0) {}
+    Vec3(float x, float y, float z) : x(x), y(y), z(z) {}
+    
+    // Arithmetic operations
+    Vec3 operator+(const Vec3& other) const {
+        return Vec3(x + other.x, y + other.y, z + other.z);
+    }
+    
+    Vec3 operator-(const Vec3& other) const {
+        return Vec3(x - other.x, y - other.y, z - other.z);
+    }
+    
+    Vec3 operator*(float scalar) const {
+        return Vec3(x * scalar, y * scalar, z * scalar);
+    }
+    
+    Vec3 operator/(float scalar) const {
+        return Vec3(x / scalar, y / scalar, z / scalar);
+    }
+    
+    // Dot product
+    float dot(const Vec3& other) const {
+        return x * other.x + y * other.y + z * other.z;
+    }
+    
+    // Cross product
+    Vec3 cross(const Vec3& other) const {
+        return Vec3(
+            y * other.z - z * other.y,
+            z * other.x - x * other.z,
+            x * other.y - y * other.x
+        );
+    }
+    
+    // Length and normalization
+    float length() const {
+        return std::sqrt(x * x + y * y + z * z);
+    }
+    
+    Vec3 normalized() const {
+        float len = length();
+        if (len > 0) {
+            return *this / len;
+        }
+        return *this;
+    }
+    
+    // Component-wise absolute value
+    Vec3 abs() const {
+        return Vec3(std::abs(x), std::abs(y), std::abs(z));
+    }
+    
+    // Component-wise floor
+    Vec3 floor() const {
+        return Vec3(std::floor(x), std::floor(y), std::floor(z));
+    }
+    
+    // Equality operator for use in hash maps
+    bool operator==(const Vec3& other) const {
+        return x == other.x && y == other.y && z == other.z;
+    }
+    
+    inline Vec3 operator*(const Vec3& scalar) const {
+        return Vec3(x * scalar.x, y * scalar.y, z * scalar.z);
+    }
+    
+    inline Vec3 operator/(const Vec3& scalar) const {
+        return Vec3(x / scalar.x, y / scalar.y, z / scalar.z);
+    }
+};
+
+// Hash function for Vec3 to use in unordered_map
+namespace std {
+    template<>
+    struct hash<Vec3> {
+        size_t operator()(const Vec3& v) const {
+            return hash<float>()(v.x) ^ (hash<float>()(v.y) << 1) ^ (hash<float>()(v.z) << 2);
+        }
+    };
+}
+
+class Vec4 {
+public:
+    float r, g, b, a;
+    
+    // Constructors
+    Vec4() : r(0), g(0), b(0), a(1.0f) {}
+    Vec4(float r, float g, float b, float a = 1.0f) : r(r), g(g), b(b), a(a) {}
+    
+    // Construct from Vec3 with alpha
+    Vec4(const Vec3& rgb, float a = 1.0f) : r(rgb.x), g(rgb.y), b(rgb.z), a(a) {}
+    
+    // Arithmetic operations
+    Vec4 operator+(const Vec4& other) const {
+        return Vec4(r + other.r, g + other.g, b + other.b, a + other.a);
+    }
+    
+    Vec4 operator-(const Vec4& other) const {
+        return Vec4(r - other.r, g - other.g, b - other.b, a - other.a);
+    }
+    
+    Vec4 operator*(float scalar) const {
+        return Vec4(r * scalar, g * scalar, b * scalar, a * scalar);
+    }
+    
+    Vec4 operator/(float scalar) const {
+        return Vec4(r / scalar, g / scalar, b / scalar, a / scalar);
+    }
+    
+    // Component-wise multiplication
+    Vec4 operator*(const Vec4& other) const {
+        return Vec4(r * other.r, g * other.g, b * other.b, a * other.a);
+    }
+    
+    // Clamp values between 0 and 1
+    Vec4 clamped() const {
+        return Vec4(
+            std::clamp(r, 0.0f, 1.0f),
+            std::clamp(g, 0.0f, 1.0f),
+            std::clamp(b, 0.0f, 1.0f),
+            std::clamp(a, 0.0f, 1.0f)
+        );
+    }
+    
+    // Convert to Vec3 (ignoring alpha)
+    Vec3 toVec3() const {
+        return Vec3(r, g, b);
+    }
+    
+    // Convert to 8-bit color values
+    void toUint8(uint8_t& red, uint8_t& green, uint8_t& blue, uint8_t& alpha) const {
+        red = static_cast<uint8_t>(std::clamp(r, 0.0f, 1.0f) * 255);
+        green = static_cast<uint8_t>(std::clamp(g, 0.0f, 1.0f) * 255);
+        blue = static_cast<uint8_t>(std::clamp(b, 0.0f, 1.0f) * 255);
+        alpha = static_cast<uint8_t>(std::clamp(a, 0.0f, 1.0f) * 255);
+    }
+    
+    void toUint8(uint8_t& red, uint8_t& green, uint8_t& blue) const {
+        red = static_cast<uint8_t>(std::clamp(r, 0.0f, 1.0f) * 255);
+        green = static_cast<uint8_t>(std::clamp(g, 0.0f, 1.0f) * 255);
+        blue = static_cast<uint8_t>(std::clamp(b, 0.0f, 1.0f) * 255);
+    }
+};
+
+class VoxelGrid {
+private:
+    std::unordered_map<Vec3, size_t> positionToIndex;
+    std::vector<Vec3> positions;
+    std::vector<Vec4> colors;
+    
+    Vec3 gridSize;
+
+public:
+    Vec3 voxelSize;
+    VoxelGrid(const Vec3& size, const Vec3& voxelSize = Vec3(1, 1, 1)) 
+        : gridSize(size), voxelSize(voxelSize) {}
+    
+    // Add a voxel at position with color
+    void addVoxel(const Vec3& position, const Vec4& color) {
+        Vec3 gridPos = worldToGrid(position);
+        
+        // Check if voxel already exists
+        auto it = positionToIndex.find(gridPos);
+        if (it == positionToIndex.end()) {
+            // New voxel
+            size_t index = positions.size();
+            positions.push_back(gridPos);
+            colors.push_back(color);
+            positionToIndex[gridPos] = index;
+        } else {
+            // Update existing voxel (you might want to blend colors instead)
+            colors[it->second] = color;
+        }
+    }
+    
+    // Get voxel color at position
+    Vec4 getVoxel(const Vec3& position) const {
+        Vec3 gridPos = worldToGrid(position);
+        auto it = positionToIndex.find(gridPos);
+        if (it != positionToIndex.end()) {
+            return colors[it->second];
+        }
+        return Vec4(0, 0, 0, 0); // Transparent black for empty voxels
+    }
+    
+    // Check if position is occupied
+    bool isOccupied(const Vec3& position) const {
+        Vec3 gridPos = worldToGrid(position);
+        return positionToIndex.find(gridPos) != positionToIndex.end();
+    }
+    
+    // Convert world coordinates to grid coordinates
+    Vec3 worldToGrid(const Vec3& worldPos) const {
+        return (worldPos / voxelSize).floor();
+    }
+    
+    // Convert grid coordinates to world coordinates
+    Vec3 gridToWorld(const Vec3& gridPos) const {
+        return gridPos * voxelSize;
+    }
+    
+    // Get all occupied positions
+    const std::vector<Vec3>& getOccupiedPositions() const {
+        return positions;
+    }
+    
+    // Get all colors
+    const std::vector<Vec4>& getColors() const {
+        return colors;
+    }
+    
+    // Get the mapping from position to index
+    const std::unordered_map<Vec3, size_t>& getPositionToIndexMap() const {
+        return positionToIndex;
+    }
+    
+    // Get grid size
+    const Vec3& getGridSize() const {
+        return gridSize;
+    }
+    
+    // Get voxel size
+    const Vec3& getVoxelSize() const {
+        return voxelSize;
+    }
+    
+    // Clear the grid
+    void clear() {
+        positions.clear();
+        colors.clear();
+        positionToIndex.clear();
+    }
+};
+
+class AmanatidesWooAlgorithm {
+public:
+    struct Ray {
+        Vec3 origin;
+        Vec3 direction;
+        float tMax;
+        
+        Ray(const Vec3& origin, const Vec3& direction, float tMax = 1000.0f)
+            : origin(origin), direction(direction.normalized()), tMax(tMax) {}
+    };
+    
+    struct TraversalState {
+        Vec3 currentVoxel;
+        Vec3 tMax;
+        Vec3 tDelta;
+        Vec3 step;
+        bool hit;
+        float t;
+        
+        TraversalState() : hit(false), t(0) {}
+    };
+    
+    // Initialize traversal state for a ray
+    static TraversalState initTraversal(const Ray& ray, const Vec3& voxelSize) {
+        TraversalState state;
+        
+        // Find the starting voxel
+        Vec3 startPos = ray.origin / voxelSize;
+        state.currentVoxel = startPos.floor();
+        
+        // Determine step directions and initialize tMax
+        Vec3 rayDir = ray.direction;
+        Vec3 invDir = Vec3(
+            rayDir.x != 0 ? 1.0f / rayDir.x : std::numeric_limits<float>::max(),
+            rayDir.y != 0 ? 1.0f / rayDir.y : std::numeric_limits<float>::max(),
+            rayDir.z != 0 ? 1.0f / rayDir.z : std::numeric_limits<float>::max()
+        );
+        
+        // Calculate step directions
+        state.step = Vec3(
+            rayDir.x > 0 ? 1 : (rayDir.x < 0 ? -1 : 0),
+            rayDir.y > 0 ? 1 : (rayDir.y < 0 ? -1 : 0),
+            rayDir.z > 0 ? 1 : (rayDir.z < 0 ? -1 : 0)
+        );
+        
+        // Calculate tMax for each axis
+        Vec3 nextVoxelBoundary = state.currentVoxel;
+        if (state.step.x > 0) nextVoxelBoundary.x += 1;
+        if (state.step.y > 0) nextVoxelBoundary.y += 1;
+        if (state.step.z > 0) nextVoxelBoundary.z += 1;
+        
+        state.tMax = Vec3(
+            (nextVoxelBoundary.x - startPos.x) * invDir.x,
+            (nextVoxelBoundary.y - startPos.y) * invDir.y,
+            (nextVoxelBoundary.z - startPos.z) * invDir.z
+        );
+        
+        // Calculate tDelta
+        state.tDelta = Vec3(
+            state.step.x * invDir.x,
+            state.step.y * invDir.y,
+            state.step.z * invDir.z
+        );
+        
+        state.hit = false;
+        state.t = 0;
+        
+        return state;
+    }
+    
+    // Traverse the grid along the ray
+    static bool traverse(const Ray& ray, const VoxelGrid& grid, 
+                        std::vector<Vec3>& hitVoxels, std::vector<float>& hitDistances,
+                        int maxSteps = 1000) {
+        TraversalState state = initTraversal(ray, grid.voxelSize);
+        Vec3 voxelSize = grid.voxelSize;
+        
+        hitVoxels.clear();
+        hitDistances.clear();
+        
+        for (int step = 0; step < maxSteps; step++) {
+            // Check if current voxel is occupied
+            Vec3 worldPos = grid.gridToWorld(state.currentVoxel);
+            if (grid.isOccupied(worldPos)) {
+                hitVoxels.push_back(state.currentVoxel);
+                hitDistances.push_back(state.t);
+            }
+            
+            // Find next voxel
+            if (state.tMax.x < state.tMax.y) {
+                if (state.tMax.x < state.tMax.z) {
+                    state.currentVoxel.x += state.step.x;
+                    state.t = state.tMax.x;
+                    state.tMax.x += state.tDelta.x;
+                } else {
+                    state.currentVoxel.z += state.step.z;
+                    state.t = state.tMax.z;
+                    state.tMax.z += state.tDelta.z;
+                }
+            } else {
+                if (state.tMax.y < state.tMax.z) {
+                    state.currentVoxel.y += state.step.y;
+                    state.t = state.tMax.y;
+                    state.tMax.y += state.tDelta.y;
+                } else {
+                    state.currentVoxel.z += state.step.z;
+                    state.t = state.tMax.z;
+                    state.tMax.z += state.tDelta.z;
+                }
+            }
+            
+            // Check if we've exceeded maximum distance
+            if (state.t > ray.tMax) {
+                break;
+            }
+        }
+        
+        return !hitVoxels.empty();
+    }
+};
+
+class BMPWriter {
+private:
+    #pragma pack(push, 1)
+    struct BMPHeader {
+        uint16_t signature = 0x4D42; // "BM"
+        uint32_t fileSize;
+        uint16_t reserved1 = 0;
+        uint16_t reserved2 = 0;
+        uint32_t dataOffset = 54;
+    };
+    
+    struct BMPInfoHeader {
+        uint32_t headerSize = 40;
+        int32_t width;
+        int32_t height;
+        uint16_t planes = 1;
+        uint16_t bitsPerPixel = 24;
+        uint32_t compression = 0;
+        uint32_t imageSize;
+        int32_t xPixelsPerMeter = 0;
+        int32_t yPixelsPerMeter = 0;
+        uint32_t colorsUsed = 0;
+        uint32_t importantColors = 0;
+    };
+    #pragma pack(pop)
+
+public:
+    static bool saveBMP(const std::string& filename, const std::vector<uint8_t>& pixels, int width, int height) {
+        BMPHeader header;
+        BMPInfoHeader infoHeader;
+        
+        int rowSize = (width * 3 + 3) & ~3; // 24-bit, padded to 4 bytes
+        int imageSize = rowSize * height;
+        
+        header.fileSize = sizeof(BMPHeader) + sizeof(BMPInfoHeader) + imageSize;
+        infoHeader.width = width;
+        infoHeader.height = height;
+        infoHeader.imageSize = imageSize;
+        
+        std::ofstream file(filename, std::ios::binary);
+        if (!file) {
+            return false;
+        }
+        
+        file.write(reinterpret_cast<const char*>(&header), sizeof(header));
+        file.write(reinterpret_cast<const char*>(&infoHeader), sizeof(infoHeader));
+        
+        // Write pixel data (BMP stores pixels bottom-to-top)
+        std::vector<uint8_t> row(rowSize);
+        for (int y = height - 1; y >= 0; --y) {
+            const uint8_t* src = &pixels[y * width * 3];
+            std::memcpy(row.data(), src, width * 3);
+            file.write(reinterpret_cast<const char*>(row.data()), rowSize);
+        }
+        
+        return true;
+    }
+    
+    static bool saveVoxelGridSlice(const std::string& filename, const VoxelGrid& grid, int sliceZ) {
+        Vec3 gridSize = grid.getGridSize();
+        int width = static_cast<int>(gridSize.x);
+        int height = static_cast<int>(gridSize.y);
+        
+        std::vector<uint8_t> pixels(width * height * 3, 0);
+        
+        // Render the slice
+        for (int y = 0; y < height; ++y) {
+            for (int x = 0; x < width; ++x) {
+                Vec3 worldPos = grid.gridToWorld(Vec3(x, y, sliceZ));
+                Vec4 color = grid.getVoxel(worldPos);
+                
+                int index = (y * width + x) * 3;
+                color.toUint8(pixels[index + 2], pixels[index + 1], pixels[index]); // BMP is BGR
+            }
+        }
+        
+        return saveBMP(filename, pixels, width, height);
+    }
+    
+    static bool saveRayTraceResults(const std::string& filename, const VoxelGrid& grid, 
+                                   const std::vector<Vec3>& hitVoxels, 
+                                   const AmanatidesWooAlgorithm::Ray& ray,
+                                   int width = 800, int height = 600) {
+        std::vector<uint8_t> pixels(width * height * 3, 0);
+        
+        // Background color (dark gray)
+        for (int i = 0; i < width * height * 3; i += 3) {
+            pixels[i] = 50;     // B
+            pixels[i + 1] = 50; // G
+            pixels[i + 2] = 50; // R
+        }
+        
+        // Draw the grid bounds
+        Vec3 gridSize = grid.getGridSize();
+        drawGrid(pixels, width, height, gridSize);
+        
+        // Draw hit voxels
+        for (const auto& voxel : hitVoxels) {
+            drawVoxel(pixels, width, height, voxel, gridSize, Vec4(1, 0, 0, 1)); // Red for hit voxels
+        }
+        
+        // Draw the ray
+        drawRay(pixels, width, height, ray, gridSize);
+        
+        return saveBMP(filename, pixels, width, height);
+    }
+    
+private:
+    static void drawGrid(std::vector<uint8_t>& pixels, int width, int height, const Vec3& gridSize) {
+        // Draw grid boundaries in white
+        for (int x = 0; x <= gridSize.x; ++x) {
+            int screenX = mapToScreenX(x, gridSize.x, width);
+            for (int y = 0; y < height; ++y) {
+                int index = (y * width + screenX) * 3;
+                if (y % 5 == 0) { // Dotted line
+                    pixels[index] = 255;     // B
+                    pixels[index + 1] = 255; // G
+                    pixels[index + 2] = 255; // R
+                }
+            }
+        }
+        
+        for (int y = 0; y <= gridSize.y; ++y) {
+            int screenY = mapToScreenY(y, gridSize.y, height);
+            for (int x = 0; x < width; ++x) {
+                int index = (screenY * width + x) * 3;
+                if (x % 5 == 0) { // Dotted line
+                    pixels[index] = 255;     // B
+                    pixels[index + 1] = 255; // G
+                    pixels[index + 2] = 255; // R
+                }
+            }
+        }
+    }
+    
+    static void drawVoxel(std::vector<uint8_t>& pixels, int width, int height, 
+                         const Vec3& voxel, const Vec3& gridSize, const Vec4& color) {
+        int screenX = mapToScreenX(voxel.x, gridSize.x, width);
+        int screenY = mapToScreenY(voxel.y, gridSize.y, height);
+        
+        uint8_t r, g, b;
+        color.toUint8(r, g, b);
+        
+        // Draw a 4x4 square for the voxel
+        for (int dy = -2; dy <= 2; ++dy) {
+            for (int dx = -2; dx <= 2; ++dx) {
+                int px = screenX + dx;
+                int py = screenY + dy;
+                if (px >= 0 && px < width && py >= 0 && py < height) {
+                    int index = (py * width + px) * 3;
+                    pixels[index] = b;
+                    pixels[index + 1] = g;
+                    pixels[index + 2] = r;
+                }
+            }
+        }
+    }
+    
+    static void drawRay(std::vector<uint8_t>& pixels, int width, int height,
+                       const AmanatidesWooAlgorithm::Ray& ray, const Vec3& gridSize) {
+        // Draw ray origin and direction
+        int originX = mapToScreenX(ray.origin.x, gridSize.x, width);
+        int originY = mapToScreenY(ray.origin.y, gridSize.y, height);
+        
+        // Draw ray origin (green)
+        for (int dy = -3; dy <= 3; ++dy) {
+            for (int dx = -3; dx <= 3; ++dx) {
+                int px = originX + dx;
+                int py = originY + dy;
+                if (px >= 0 && px < width && py >= 0 && py < height) {
+                    int index = (py * width + px) * 3;
+                    pixels[index] = 0;       // B
+                    pixels[index + 1] = 255; // G
+                    pixels[index + 2] = 0;   // R
+                }
+            }
+        }
+        
+        // Draw ray direction (yellow line)
+        Vec3 endPoint = ray.origin + ray.direction * 10.0f; // Extend ray
+        int endX = mapToScreenX(endPoint.x, gridSize.x, width);
+        int endY = mapToScreenY(endPoint.y, gridSize.y, height);
+        
+        drawLine(pixels, width, height, originX, originY, endX, endY, Vec4(1, 1, 0, 1));
+    }
+    
+    static void drawLine(std::vector<uint8_t>& pixels, int width, int height,
+                        int x0, int y0, int x1, int y1, const Vec4& color) {
+        uint8_t r, g, b;
+        color.toUint8(r, g, b);
+        
+        int dx = std::abs(x1 - x0);
+        int dy = std::abs(y1 - y0);
+        int sx = (x0 < x1) ? 1 : -1;
+        int sy = (y0 < y1) ? 1 : -1;
+        int err = dx - dy;
+        
+        while (true) {
+            if (x0 >= 0 && x0 < width && y0 >= 0 && y0 < height) {
+                int index = (y0 * width + x0) * 3;
+                pixels[index] = b;
+                pixels[index + 1] = g;
+                pixels[index + 2] = r;
+            }
+            
+            if (x0 == x1 && y0 == y1) break;
+            
+            int e2 = 2 * err;
+            if (e2 > -dy) {
+                err -= dy;
+                x0 += sx;
+            }
+            if (e2 < dx) {
+                err += dx;
+                y0 += sy;
+            }
+        }
+    }
+    
+    static int mapToScreenX(float x, float gridWidth, int screenWidth) {
+        return static_cast<int>((x / gridWidth) * screenWidth);
+    }
+    
+    static int mapToScreenY(float y, float gridHeight, int screenHeight) {
+        return static_cast<int>((y / gridHeight) * screenHeight);
+    }
+};
+
+// Example usage
+int main() {
+    // Create a voxel grid
+    VoxelGrid grid(Vec3(10, 10, 10), Vec3(1, 1, 1));
+    
+    // Add some voxels
+    grid.addVoxel(Vec3(1, 1, 1), Vec4(1, 0, 0, 1)); // Red
+    grid.addVoxel(Vec3(2, 2, 2), Vec4(0, 1, 0, 0.5f)); // Green with transparency
+    grid.addVoxel(Vec3(3, 3, 3), Vec4(0, 0, 1, 1)); // Blue
+    grid.addVoxel(Vec3(4, 4, 4), Vec4(1, 1, 0, 1)); // Yellow
+    grid.addVoxel(Vec3(5, 5, 5), Vec4(1, 0, 1, 1)); // Magenta
+    
+    // Create a ray
+    AmanatidesWooAlgorithm::Ray ray(Vec3(0, 0, 0), Vec3(1, 1, 1).normalized());
+    
+    // Traverse the grid
+    std::vector<Vec3> hitVoxels;
+    std::vector<float> hitDistances;
+    bool hit = AmanatidesWooAlgorithm::traverse(ray, grid, hitVoxels, hitDistances);
+    
+    if (hit) {
+        printf("Ray hit %zu voxels:\n", hitVoxels.size());
+        for (size_t i = 0; i < hitVoxels.size(); i++) {
+            Vec4 color = grid.getVoxel(grid.gridToWorld(hitVoxels[i]));
+            printf("  Voxel at (%.1f, %.1f, %.1f), distance: %.2f, color: (%.1f, %.1f, %.1f, %.1f)\n",
+                   hitVoxels[i].x, hitVoxels[i].y, hitVoxels[i].z,
+                   hitDistances[i],
+                   color.r, color.g, color.b, color.a);
+        }
+    }
+    
+    // Save results to BMP files
+    printf("\nSaving results to BMP files...\n");
+    
+    // Save a slice of the voxel grid
+    if (BMPWriter::saveVoxelGridSlice("voxel_grid_slice.bmp", grid, 1)) {
+        printf("Saved voxel grid slice to 'voxel_grid_slice.bmp'\n");
+    } else {
+        printf("Failed to save voxel grid slice\n");
+    }
+    
+    // Save ray tracing visualization
+    if (BMPWriter::saveRayTraceResults("ray_trace_results.bmp", grid, hitVoxels, ray)) {
+        printf("Saved ray trace results to 'ray_trace_results.bmp'\n");
+    } else {
+        printf("Failed to save ray trace results\n");
+    }
+    
+    return 0;
+}
\ No newline at end of file
diff --git a/e.cpp b/e.cpp
new file mode 100644
index 0000000..0397384
--- /dev/null
+++ b/e.cpp
@@ -0,0 +1,854 @@
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <algorithm>
+#include <unordered_map>
+#include <cstdio>
+#include <cstring>
+#include <fstream>
+#include <iostream>
+#include <vector>
+#include <random>
+#include <functional>
+#include <tuple>
+#include "timing_decorator.hpp"
+#include "util/vec3.hpp"
+#include "util/vec4.hpp"
+
+class VoxelGrid {
+private:
+    std::unordered_map<Vec3, size_t> positionToIndex;
+    std::vector<Vec3> positions;
+    std::vector<Vec4> colors;
+    
+    Vec3 gridSize;
+
+public:
+    Vec3 voxelSize;
+    VoxelGrid(const Vec3& size, const Vec3& voxelSize = Vec3(1, 1, 1)) 
+        : gridSize(size), voxelSize(voxelSize) {}
+    
+    // Add a voxel at position with color
+    void addVoxel(const Vec3& position, const Vec4& color) {
+        Vec3 gridPos = worldToGrid(position);
+        
+        // Check if voxel already exists
+        auto it = positionToIndex.find(gridPos);
+        if (it == positionToIndex.end()) {
+            // New voxel
+            size_t index = positions.size();
+            positions.push_back(gridPos);
+            colors.push_back(color);
+            positionToIndex[gridPos] = index;
+        } else {
+            // Update existing voxel (you might want to blend colors instead)
+            colors[it->second] = color;
+        }
+    }
+    
+    // Get voxel color at position
+    Vec4 getVoxel(const Vec3& position) const {
+        Vec3 gridPos = worldToGrid(position);
+        auto it = positionToIndex.find(gridPos);
+        if (it != positionToIndex.end()) {
+            return colors[it->second];
+        }
+        return Vec4(0, 0, 0, 0); // Transparent black for empty voxels
+    }
+    
+    // Check if position is occupied
+    bool isOccupied(const Vec3& position) const {
+        Vec3 gridPos = worldToGrid(position);
+        return positionToIndex.find(gridPos) != positionToIndex.end();
+    }
+    
+    // Convert world coordinates to grid coordinates
+    Vec3 worldToGrid(const Vec3& worldPos) const {
+        return (worldPos / voxelSize).floor();
+    }
+    
+    // Convert grid coordinates to world coordinates
+    Vec3 gridToWorld(const Vec3& gridPos) const {
+        return gridPos * voxelSize;
+    }
+    
+    // Get all occupied positions
+    const std::vector<Vec3>& getOccupiedPositions() const {
+        return positions;
+    }
+    
+    // Get all colors
+    const std::vector<Vec4>& getColors() const {
+        return colors;
+    }
+    
+    // Get the mapping from position to index
+    const std::unordered_map<Vec3, size_t>& getPositionToIndexMap() const {
+        return positionToIndex;
+    }
+    
+    // Get grid size
+    const Vec3& getGridSize() const {
+        return gridSize;
+    }
+    
+    // Get voxel size
+    const Vec3& getVoxelSize() const {
+        return voxelSize;
+    }
+    
+    // Clear the grid
+    void clear() {
+        positions.clear();
+        colors.clear();
+        positionToIndex.clear();
+    }
+};
+
+class AmanatidesWooAlgorithm {
+public:
+    struct Ray {
+        Vec3 origin;
+        Vec3 direction;
+        float tMax;
+        
+        Ray(const Vec3& origin, const Vec3& direction, float tMax = 1000.0f)
+            : origin(origin), direction(direction.normalized()), tMax(tMax) {}
+    };
+    
+    struct TraversalState {
+        Vec3 currentVoxel;
+        Vec3 tMax;
+        Vec3 tDelta;
+        Vec3 step;
+        bool hit;
+        float t;
+        
+        TraversalState() : hit(false), t(0) {}
+    };
+    
+    // Initialize traversal state for a ray
+    static TraversalState initTraversal(const Ray& ray, const Vec3& voxelSize) {
+        TraversalState state;
+        
+        // Find the starting voxel
+        Vec3 startPos = ray.origin / voxelSize;
+        state.currentVoxel = startPos.floor();
+        
+        // Determine step directions and initialize tMax
+        Vec3 rayDir = ray.direction;
+        Vec3 invDir = Vec3(
+            rayDir.x != 0 ? 1.0f / rayDir.x : std::numeric_limits<float>::max(),
+            rayDir.y != 0 ? 1.0f / rayDir.y : std::numeric_limits<float>::max(),
+            rayDir.z != 0 ? 1.0f / rayDir.z : std::numeric_limits<float>::max()
+        );
+        
+        // Calculate step directions
+        state.step = Vec3(
+            rayDir.x > 0 ? 1 : (rayDir.x < 0 ? -1 : 0),
+            rayDir.y > 0 ? 1 : (rayDir.y < 0 ? -1 : 0),
+            rayDir.z > 0 ? 1 : (rayDir.z < 0 ? -1 : 0)
+        );
+        
+        // Calculate tMax for each axis
+        Vec3 nextVoxelBoundary = state.currentVoxel;
+        if (state.step.x > 0) nextVoxelBoundary.x += 1;
+        if (state.step.y > 0) nextVoxelBoundary.y += 1;
+        if (state.step.z > 0) nextVoxelBoundary.z += 1;
+        
+        state.tMax = Vec3(
+            (nextVoxelBoundary.x - startPos.x) * invDir.x,
+            (nextVoxelBoundary.y - startPos.y) * invDir.y,
+            (nextVoxelBoundary.z - startPos.z) * invDir.z
+        );
+        
+        // Calculate tDelta
+        state.tDelta = Vec3(
+            state.step.x * invDir.x,
+            state.step.y * invDir.y,
+            state.step.z * invDir.z
+        );
+        
+        state.hit = false;
+        state.t = 0;
+        
+        return state;
+    }
+    
+    // Traverse the grid along the ray
+    static bool traverse(const Ray& ray, const VoxelGrid& grid, 
+                        std::vector<Vec3>& hitVoxels, std::vector<float>& hitDistances,
+                        int maxSteps = 1000) {
+        TraversalState state = initTraversal(ray, grid.voxelSize);
+        Vec3 voxelSize = grid.voxelSize;
+        
+        hitVoxels.clear();
+        hitDistances.clear();
+        
+        for (int step = 0; step < maxSteps; step++) {
+            // Check if current voxel is occupied
+            Vec3 worldPos = grid.gridToWorld(state.currentVoxel);
+            if (grid.isOccupied(worldPos)) {
+                hitVoxels.push_back(state.currentVoxel);
+                hitDistances.push_back(state.t);
+            }
+            
+            // Find next voxel
+            if (state.tMax.x < state.tMax.y) {
+                if (state.tMax.x < state.tMax.z) {
+                    state.currentVoxel.x += state.step.x;
+                    state.t = state.tMax.x;
+                    state.tMax.x += state.tDelta.x;
+                } else {
+                    state.currentVoxel.z += state.step.z;
+                    state.t = state.tMax.z;
+                    state.tMax.z += state.tDelta.z;
+                }
+            } else {
+                if (state.tMax.y < state.tMax.z) {
+                    state.currentVoxel.y += state.step.y;
+                    state.t = state.tMax.y;
+                    state.tMax.y += state.tDelta.y;
+                } else {
+                    state.currentVoxel.z += state.step.z;
+                    state.t = state.tMax.z;
+                    state.tMax.z += state.tDelta.z;
+                }
+            }
+            
+            // Check if we've exceeded maximum distance
+            if (state.t > ray.tMax) {
+                break;
+            }
+        }
+        
+        return !hitVoxels.empty();
+    }
+};
+
+class BMPWriter {
+private:
+    #pragma pack(push, 1)
+    struct BMPHeader {
+        uint16_t signature = 0x4D42; // "BM"
+        uint32_t fileSize;
+        uint16_t reserved1 = 0;
+        uint16_t reserved2 = 0;
+        uint32_t dataOffset = 54;
+    };
+    
+    struct BMPInfoHeader {
+        uint32_t headerSize = 40;
+        int32_t width;
+        int32_t height;
+        uint16_t planes = 1;
+        uint16_t bitsPerPixel = 24;
+        uint32_t compression = 0;
+        uint32_t imageSize;
+        int32_t xPixelsPerMeter = 0;
+        int32_t yPixelsPerMeter = 0;
+        uint32_t colorsUsed = 0;
+        uint32_t importantColors = 0;
+    };
+    #pragma pack(pop)
+
+public:
+    static bool saveBMP(const std::string& filename, const std::vector<uint8_t>& pixels, int width, int height) {
+        BMPHeader header;
+        BMPInfoHeader infoHeader;
+        
+        int rowSize = (width * 3 + 3) & ~3; // 24-bit, padded to 4 bytes
+        int imageSize = rowSize * height;
+        
+        header.fileSize = sizeof(BMPHeader) + sizeof(BMPInfoHeader) + imageSize;
+        infoHeader.width = width;
+        infoHeader.height = height;
+        infoHeader.imageSize = imageSize;
+        
+        std::ofstream file(filename, std::ios::binary);
+        if (!file) {
+            return false;
+        }
+        
+        file.write(reinterpret_cast<const char*>(&header), sizeof(header));
+        file.write(reinterpret_cast<const char*>(&infoHeader), sizeof(infoHeader));
+        
+        // Write pixel data (BMP stores pixels bottom-to-top)
+        std::vector<uint8_t> row(rowSize);
+        for (int y = height - 1; y >= 0; --y) {
+            const uint8_t* src = &pixels[y * width * 3];
+            std::memcpy(row.data(), src, width * 3);
+            file.write(reinterpret_cast<const char*>(row.data()), rowSize);
+        }
+        
+        return true;
+    }
+    
+    static bool saveVoxelGridSlice(const std::string& filename, const VoxelGrid& grid, int sliceZ) {
+        Vec3 gridSize = grid.getGridSize();
+        int width = static_cast<int>(gridSize.x);
+        int height = static_cast<int>(gridSize.y);
+        
+        std::vector<uint8_t> pixels(width * height * 3, 0);
+        
+        // Render the slice
+        for (int y = 0; y < height; ++y) {
+            for (int x = 0; x < width; ++x) {
+                Vec3 worldPos = grid.gridToWorld(Vec3(x, y, sliceZ));
+                Vec4 color = grid.getVoxel(worldPos);
+                
+                int index = (y * width + x) * 3;
+                color.toUint8(pixels[index + 2], pixels[index + 1], pixels[index]); // BMP is BGR
+            }
+        }
+        
+        return saveBMP(filename, pixels, width, height);
+    }
+    
+    static bool saveRayTraceResults(const std::string& filename, const VoxelGrid& grid, 
+                                   const std::vector<Vec3>& hitVoxels, 
+                                   const AmanatidesWooAlgorithm::Ray& ray,
+                                   int width = 800, int height = 600) {
+        std::vector<uint8_t> pixels(width * height * 3, 0);
+        
+        // Background color (dark gray)
+        for (int i = 0; i < width * height * 3; i += 3) {
+            pixels[i] = 50;     // B
+            pixels[i + 1] = 50; // G
+            pixels[i + 2] = 50; // R
+        }
+        
+        // Draw the grid bounds
+        Vec3 gridSize = grid.getGridSize();
+        drawGrid(pixels, width, height, gridSize);
+        
+        // Draw hit voxels
+        for (const auto& voxel : hitVoxels) {
+            drawVoxel(pixels, width, height, voxel, gridSize, Vec4(1, 0, 0, 1)); // Red for hit voxels
+        }
+        
+        // Draw the ray
+        drawRay(pixels, width, height, ray, gridSize);
+        
+        return saveBMP(filename, pixels, width, height);
+    }
+    
+private:
+    static void drawGrid(std::vector<uint8_t>& pixels, int width, int height, const Vec3& gridSize) {
+        // Draw grid boundaries in white
+        for (int x = 0; x <= gridSize.x; ++x) {
+            int screenX = mapToScreenX(x, gridSize.x, width);
+            for (int y = 0; y < height; ++y) {
+                int index = (y * width + screenX) * 3;
+                if (y % 5 == 0) { // Dotted line
+                    pixels[index] = 255;     // B
+                    pixels[index + 1] = 255; // G
+                    pixels[index + 2] = 255; // R
+                }
+            }
+        }
+        
+        for (int y = 0; y <= gridSize.y; ++y) {
+            int screenY = mapToScreenY(y, gridSize.y, height);
+            for (int x = 0; x < width; ++x) {
+                int index = (screenY * width + x) * 3;
+                if (x % 5 == 0) { // Dotted line
+                    pixels[index] = 255;     // B
+                    pixels[index + 1] = 255; // G
+                    pixels[index + 2] = 255; // R
+                }
+            }
+        }
+    }
+    
+    static void drawVoxel(std::vector<uint8_t>& pixels, int width, int height, 
+                         const Vec3& voxel, const Vec3& gridSize, const Vec4& color) {
+        int screenX = mapToScreenX(voxel.x, gridSize.x, width);
+        int screenY = mapToScreenY(voxel.y, gridSize.y, height);
+        
+        uint8_t r, g, b;
+        color.toUint8(r, g, b);
+        
+        // Draw a 4x4 square for the voxel
+        for (int dy = -2; dy <= 2; ++dy) {
+            for (int dx = -2; dx <= 2; ++dx) {
+                int px = screenX + dx;
+                int py = screenY + dy;
+                if (px >= 0 && px < width && py >= 0 && py < height) {
+                    int index = (py * width + px) * 3;
+                    pixels[index] = b;
+                    pixels[index + 1] = g;
+                    pixels[index + 2] = r;
+                }
+            }
+        }
+    }
+    
+    static void drawRay(std::vector<uint8_t>& pixels, int width, int height,
+                       const AmanatidesWooAlgorithm::Ray& ray, const Vec3& gridSize) {
+        // Draw ray origin and direction
+        int originX = mapToScreenX(ray.origin.x, gridSize.x, width);
+        int originY = mapToScreenY(ray.origin.y, gridSize.y, height);
+        
+        // Draw ray origin (green)
+        for (int dy = -3; dy <= 3; ++dy) {
+            for (int dx = -3; dx <= 3; ++dx) {
+                int px = originX + dx;
+                int py = originY + dy;
+                if (px >= 0 && px < width && py >= 0 && py < height) {
+                    int index = (py * width + px) * 3;
+                    pixels[index] = 0;       // B
+                    pixels[index + 1] = 255; // G
+                    pixels[index + 2] = 0;   // R
+                }
+            }
+        }
+        
+        // Draw ray direction (yellow line)
+        Vec3 endPoint = ray.origin + ray.direction * 10.0f; // Extend ray
+        int endX = mapToScreenX(endPoint.x, gridSize.x, width);
+        int endY = mapToScreenY(endPoint.y, gridSize.y, height);
+        
+        drawLine(pixels, width, height, originX, originY, endX, endY, Vec4(1, 1, 0, 1));
+    }
+    
+    static void drawLine(std::vector<uint8_t>& pixels, int width, int height,
+                        int x0, int y0, int x1, int y1, const Vec4& color) {
+        uint8_t r, g, b;
+        color.toUint8(r, g, b);
+        
+        int dx = std::abs(x1 - x0);
+        int dy = std::abs(y1 - y0);
+        int sx = (x0 < x1) ? 1 : -1;
+        int sy = (y0 < y1) ? 1 : -1;
+        int err = dx - dy;
+        
+        while (true) {
+            if (x0 >= 0 && x0 < width && y0 >= 0 && y0 < height) {
+                int index = (y0 * width + x0) * 3;
+                pixels[index] = b;
+                pixels[index + 1] = g;
+                pixels[index + 2] = r;
+            }
+            
+            if (x0 == x1 && y0 == y1) break;
+            
+            int e2 = 2 * err;
+            if (e2 > -dy) {
+                err -= dy;
+                x0 += sx;
+            }
+            if (e2 < dx) {
+                err += dx;
+                y0 += sy;
+            }
+        }
+    }
+    
+    static int mapToScreenX(float x, float gridWidth, int screenWidth) {
+        return static_cast<int>((x / gridWidth) * screenWidth);
+    }
+    
+    static int mapToScreenY(float y, float gridHeight, int screenHeight) {
+        return static_cast<int>((y / gridHeight) * screenHeight);
+    }
+};
+
+// Noise generation functions
+float fade(const float& a) {
+    TIME_FUNCTION;
+    return a * a * a * (10 + a * (-15 + a * 6));
+}
+
+float clamp(float x, float lowerlimit = 0.0f, float upperlimit = 1.0f) {
+    TIME_FUNCTION;
+    if (x < lowerlimit) return lowerlimit;
+    if (x > upperlimit) return upperlimit;
+    return x;
+}
+
+float pascalTri(const float& a, const float& b) {
+    TIME_FUNCTION;
+    int result = 1;
+    for (int i = 0; i < b; ++i){
+        result *= (a - 1) / (i + 1);
+    }
+    return result;
+}
+
+float genSmooth(int N, float x) {
+    TIME_FUNCTION;
+    x = clamp(x, 0, 1);
+    float result = 0;
+    for (int n = 0; n <= N; ++n){
+        result += pascalTri(-N - 1, n) * pascalTri(2 * N + 1, N-1) * pow(x, N + n + 1);
+    }
+    return result;
+}
+
+float inverse_smoothstep(float x) {
+    TIME_FUNCTION;
+    return 0.5 - sin(asin(1.0 - 2.0 * x) / 3.0);
+}
+
+float lerp(const float& t, const float& a, const float& b) {
+    TIME_FUNCTION;
+    return a + t * (b - a);
+}
+
+float grad(const int& hash, const float& b, const float& c, const float& d) {
+    TIME_FUNCTION;
+    int h = hash & 15;
+    float u = (h < 8) ? c : b;
+    float v = (h < 4) ? b : ((h == 12 || h == 14) ? c : d);
+    return (((h & 1) == 0) ? u : -u) + (((h & 2) == 0) ? v : -v);
+}
+
+float pnoise3d(const int p[512], const float& xf, const float& yf, const float& zf) {
+    TIME_FUNCTION;
+    int floorx = std::floor(xf);
+    int floory = std::floor(yf);
+    int floorz = std::floor(zf);
+    int iX = floorx & 255;
+    int iY = floory & 255;
+    int iZ = floorz & 255;
+
+    float x = xf - floorx;
+    float y = yf - floory;
+    float z = zf - floorz;
+
+    float u = fade(x);
+    float v = fade(y);
+    float w = fade(z);
+
+    int A  = p[iX] + iY;
+    int AA = p[A] + iZ;
+    int AB = p[A+1] + iZ;
+
+    int B  = p[iX + 1] + iY;
+    int BA = p[B] + iZ;
+    int BB = p[B+1] + iZ;
+
+    float f = grad(p[BA], x-1, y, z);
+    float g = grad(p[AA], x, y, z);
+    float h = grad(p[BB], x-1, y-1, z);
+    float j = grad(p[BB], x-1, y-1, z);
+    float k = grad(p[AA+1], x, y, z-1);
+    float l = grad(p[BA+1], x-1, y, z-1);
+    float m = grad(p[AB+1], x, y-1, z-1);
+    float n = grad(p[BB+1], x-1, y-1, z-1);
+
+    float o = lerp(u, m, n);
+    float q = lerp(u, k, l);
+    float r = lerp(u, h, j);
+    float s = lerp(u, f, g);
+    float t = lerp(v, q, o);
+    float e = lerp(v, s, r);
+    float d = lerp(w, e, t);
+    return d;
+}
+
+std::tuple<std::vector<Vec3>, std::vector<Vec4>> noiseBatch(int num_points, float scale, int sp[]) {
+    TIME_FUNCTION;
+    std::vector<Vec3> points;
+    std::vector<Vec4> colors;
+    points.reserve(num_points);
+    colors.reserve(num_points);
+
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<float> dis(-scale, scale);
+    
+    for (int i = 0; i < num_points; ++i) {
+        float x = dis(gen);
+        float y = dis(gen);
+        float z = dis(gen);
+        
+        float noise1 = pnoise3d(sp, (x * 0.5f), (y * 0.5f), (z * 0.5f));
+        float noise2 = pnoise3d(sp, (x * 0.3f), (y * 0.3f), (z * 0.3f));
+        float noise3 = pnoise3d(sp, (x * 0.7f), (y * 0.7f), (z * 0.7f));
+        float noise4 = pnoise3d(sp, (x * 0.7f), (y * 0.7f), (z * 0.7f));
+        
+        if (noise1 > 0.1f) {
+            float rt = (noise1 + 1.0f) * 0.5f;
+            float gt = (noise2 + 1.0f) * 0.5f;
+            float bt = (noise3 + 1.0f) * 0.5f;
+            float at = (noise4 + 1.0f) * 0.5f;
+
+            float maxV = std::max({rt, gt, bt});
+            if (maxV > 0) {
+                float r = rt / maxV;
+                float g = gt / maxV;
+                float b = bt / maxV;
+                float a = at / maxV;
+                points.push_back({x, y, z});
+                colors.push_back({r, g, b, a});
+            }
+        }
+    }
+    return std::make_tuple(points, colors);
+}
+
+// Generate points in a more spherical distribution
+std::tuple<std::vector<Vec3>, std::vector<Vec4>> genPointCloud(int num_points, float radius, int seed) {
+    TIME_FUNCTION;
+    int permutation[256];
+    for (int i = 0; i < 256; ++i) {
+        permutation[i] = i;
+    }
+    std::mt19937 rng(seed);
+    std::shuffle(permutation, permutation+256, rng);
+    int p[512];
+    for (int i = 0; i < 256; ++i) {
+        p[i] = permutation[i];
+        p[i + 256] = permutation[i];
+    }
+    
+    std::vector<Vec3> points;
+    std::vector<Vec4> colors;
+    points.reserve(num_points);
+    colors.reserve(num_points);
+
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
+    std::uniform_real_distribution<float> radius_dist(0.0f, radius);
+    
+    for (int i = 0; i < num_points; ++i) {
+        // Generate random direction
+        Vec3 direction;
+        do {
+            direction = Vec3(dist(gen), dist(gen), dist(gen));
+        } while (direction.lengthSquared() == 0);
+        direction = direction.normalized();
+        
+        // Generate random radius with some noise
+        float r = radius_dist(gen);
+        
+        Vec3 point = direction * r;
+        
+        // Use noise based on spherical coordinates for more natural distribution
+        float noise = pnoise3d(p, point.x * 0.5f, point.y * 0.5f, point.z * 0.5f);
+        
+        if (noise > 0.1f) {
+            // Color based on position and noise
+            float rt = (point.x / radius + 1.0f) * 0.5f;
+            float gt = (point.y / radius + 1.0f) * 0.5f;
+            float bt = (point.z / radius + 1.0f) * 0.5f;
+            float at = (noise + 1.0f) * 0.5f;
+
+            points.push_back(point);
+            colors.push_back({rt, gt, bt, at});
+        }
+    }
+    return std::make_tuple(points, colors);
+}
+
+// Function to populate voxel grid with point cloud data
+void populateVoxelGridWithPointCloud(VoxelGrid& grid, 
+                                   const std::vector<Vec3>& points, 
+                                   const std::vector<Vec4>& colors) {
+    TIME_FUNCTION;
+    printf("Populating voxel grid with %zu points...\n", points.size());
+    
+    for (size_t i = 0; i < points.size(); i++) {
+        grid.addVoxel(points[i], colors[i]);
+    }
+    
+    printf("Voxel grid populated with %zu voxels\n", grid.getOccupiedPositions().size());
+}
+
+// Enhanced visualization function
+void visualizePointCloud(const std::vector<Vec3>& points, const std::vector<Vec4>& colors, 
+                        const std::string& filename, int width = 800, int height = 600) {
+    TIME_FUNCTION;
+    std::vector<uint8_t> pixels(width * height * 3, 0);
+    
+    // Background color (dark blue)
+    for (int i = 0; i < width * height * 3; i += 3) {
+        pixels[i] = 30;     // B
+        pixels[i + 1] = 30; // G
+        pixels[i + 2] = 50; // R
+    }
+    
+    // Find bounds of point cloud
+    Vec3 minPoint( std::numeric_limits<float>::max(),  std::numeric_limits<float>::max(),  std::numeric_limits<float>::max());
+    Vec3 maxPoint(-std::numeric_limits<float>::max(), -std::numeric_limits<float>::max(), -std::numeric_limits<float>::max());
+    
+    for (const auto& point : points) {
+        minPoint.x = std::min(minPoint.x, point.x);
+        minPoint.y = std::min(minPoint.y, point.y);
+        minPoint.z = std::min(minPoint.z, point.z);
+        maxPoint.x = std::max(maxPoint.x, point.x);
+        maxPoint.y = std::max(maxPoint.y, point.y);
+        maxPoint.z = std::max(maxPoint.z, point.z);
+    }
+    
+    Vec3 cloudSize = maxPoint - minPoint;
+    float maxDim = std::max({cloudSize.x, cloudSize.y, cloudSize.z});
+    
+    // Draw points
+    for (size_t i = 0; i < points.size(); i++) {
+        const auto& point = points[i];
+        const auto& color = colors[i];
+        
+        // Map 3D point to 2D screen coordinates (orthographic projection)
+        int screenX = static_cast<int>(((point.x - minPoint.x) / maxDim) * (width - 20)) + 10;
+        int screenY = static_cast<int>(((point.y - minPoint.y) / maxDim) * (height - 20)) + 10;
+        
+        if (screenX >= 0 && screenX < width && screenY >= 0 && screenY < height) {
+            uint8_t r, g, b;
+            color.toUint8(r, g, b);
+            
+            // Draw a 2x2 pixel for each point
+            for (int dy = -1; dy <= 1; dy++) {
+                for (int dx = -1; dx <= 1; dx++) {
+                    int px = screenX + dx;
+                    int py = screenY + dy;
+                    if (px >= 0 && px < width && py >= 0 && py < height) {
+                        int index = (py * width + px) * 3;
+                        pixels[index] = b;
+                        pixels[index + 1] = g;
+                        pixels[index + 2] = r;
+                    }
+                }
+            }
+        }
+    }
+    
+    BMPWriter::saveBMP(filename, pixels, width, height);
+}
+
+// Replace your main function with this improved version:
+int main() {
+    printf("=== Point Cloud Generation and Visualization ===\n\n");
+    
+    // Generate point cloud using noise function
+    printf("Generating point cloud...\n");
+    float cloudScale = 5.0f;
+    auto [points, colors] = genPointCloud(500000, cloudScale, 42);
+    printf("Generated %zu points\n\n", points.size());
+    
+    // Calculate actual bounds of the point cloud
+    Vec3 minPoint(std::numeric_limits<float>::max());
+    Vec3 maxPoint(-std::numeric_limits<float>::max());
+    
+    for (const auto& point : points) {
+        minPoint.x = std::min(minPoint.x, point.x);
+        minPoint.y = std::min(minPoint.y, point.y);
+        minPoint.z = std::min(minPoint.z, point.z);
+        maxPoint.x = std::max(maxPoint.x, point.x);
+        maxPoint.y = std::max(maxPoint.y, point.y);
+        maxPoint.z = std::max(maxPoint.z, point.z);
+    }
+    
+    Vec3 cloudCenter = (minPoint + maxPoint) * 0.5f;
+    Vec3 cloudSize = maxPoint - minPoint;
+    
+    printf("Point cloud bounds:\n");
+    printf("  Min: (%.2f, %.2f, %.2f)\n", minPoint.x, minPoint.y, minPoint.z);
+    printf("  Max: (%.2f, %.2f, %.2f)\n", maxPoint.x, maxPoint.y, maxPoint.z);
+    printf("  Center: (%.2f, %.2f, %.2f)\n", cloudCenter.x, cloudCenter.y, cloudCenter.z);
+    printf("  Size: (%.2f, %.2f, %.2f)\n", cloudSize.x, cloudSize.y, cloudSize.z);
+    
+    // Create a voxel grid that properly contains the point cloud
+    // Add some padding around the cloud
+    float padding = 2.0f;
+    Vec3 gridWorldSize = cloudSize + Vec3(padding * 2);
+    Vec3 gridWorldMin = cloudCenter - gridWorldSize * 0.5f;
+    
+    // Use smaller voxels for better resolution
+    Vec3 voxelSize(0.1f, 0.1f, 0.1f);
+    Vec3 gridSize = (gridWorldSize / voxelSize).ceil();
+    
+    printf("\nVoxel grid configuration:\n");
+    printf("  World size: (%.2f, %.2f, %.2f)\n", gridWorldSize.x, gridWorldSize.y, gridWorldSize.z);
+    printf("  Grid dimensions: (%d, %d, %d)\n", (int)gridSize.x, (int)gridSize.y, (int)gridSize.z);
+    printf("  Voxel size: (%.2f, %.2f, %.2f)\n", voxelSize.x, voxelSize.y, voxelSize.z);
+    
+    VoxelGrid grid(gridSize, voxelSize);
+    
+    // Center the point cloud in the voxel grid
+    printf("\nPopulating voxel grid...\n");
+    size_t voxelsAdded = 0;
+    for (size_t i = 0; i < points.size(); i++) {
+        // Use the original point positions - the grid will handle world-to-grid conversion
+        grid.addVoxel(points[i], colors[i]);
+        voxelsAdded++;
+    }
+    
+    printf("Voxel grid populated with %zu voxels (out of %zu points)\n", 
+           grid.getOccupiedPositions().size(), points.size());
+    
+    // Test if the cloud is properly centered by checking voxel distribution
+    auto& occupied = grid.getOccupiedPositions();
+    Vec3 gridMin(std::numeric_limits<float>::max());
+    Vec3 gridMax(-std::numeric_limits<float>::max());
+    
+    for (const auto& pos : occupied) {
+        gridMin.x = std::min(gridMin.x, pos.x);
+        gridMin.y = std::min(gridMin.y, pos.y);
+        gridMin.z = std::min(gridMin.z, pos.z);
+        gridMax.x = std::max(gridMax.x, pos.x);
+        gridMax.y = std::max(gridMax.y, pos.y);
+        gridMax.z = std::max(gridMax.z, pos.z);
+    }
+    
+    printf("\nVoxel distribution in grid:\n");
+    printf("  Grid min: (%.2f, %.2f, %.2f)\n", gridMin.x, gridMin.y, gridMin.z);
+    printf("  Grid max: (%.2f, %.2f, %.2f)\n", gridMax.x, gridMax.y, gridMax.z);
+    printf("  Grid center: (%.2f, %.2f, %.2f)\n", 
+           (gridMin.x + gridMax.x) * 0.5f, 
+           (gridMin.y + gridMax.y) * 0.5f, 
+           (gridMin.z + gridMax.z) * 0.5f);
+    
+    // Visualizations
+    printf("\nCreating visualizations...\n");
+    visualizePointCloud(points, colors, "point_cloud_visualization.bmp");
+    printf("Saved point cloud visualization to 'point_cloud_visualization.bmp'\n");
+    
+    // Save slices at different heights through the cloud
+    int centerZ = static_cast<int>(gridSize.z * 0.5f);
+    for (int offset = -2; offset <= 2; offset++) {
+        int sliceZ = centerZ + offset;
+        std::string filename = "voxel_slice_z" + std::to_string(offset) + ".bmp";
+        if (BMPWriter::saveVoxelGridSlice(filename, grid, sliceZ)) {
+            printf("Saved voxel grid slice to '%s'\n", filename.c_str());
+        }
+    }
+    
+    // Test ray tracing through the center of the cloud
+    printf("\n=== Ray Tracing Test ===\n");
+    
+    // Create rays that go through the center of the cloud from different directions
+    std::vector<AmanatidesWooAlgorithm::Ray> testRays = {
+        AmanatidesWooAlgorithm::Ray(cloudCenter - Vec3(10, 0, 0), Vec3(1, 0, 0), 20.0f),  // X direction
+        AmanatidesWooAlgorithm::Ray(cloudCenter - Vec3(0, 10, 0), Vec3(0, 1, 0), 20.0f),  // Y direction  
+        AmanatidesWooAlgorithm::Ray(cloudCenter - Vec3(0, 0, 10), Vec3(0, 0, 1), 20.0f),  // Z direction
+        AmanatidesWooAlgorithm::Ray(cloudCenter - Vec3(8, 8, 0), Vec3(1, 1, 0).normalized(), 20.0f)  // Diagonal
+    };
+    
+    for (size_t i = 0; i < testRays.size(); i++) {
+        std::vector<Vec3> hitVoxels;
+        std::vector<float> hitDistances;
+        
+        bool hit = AmanatidesWooAlgorithm::traverse(testRays[i], grid, hitVoxels, hitDistances);
+        
+        printf("Ray %zu: %s (%zu hits)\n", i, hit ? "HIT" : "MISS", hitVoxels.size());
+        
+        // Save visualization for this ray
+        std::string rayFilename = "ray_trace_" + std::to_string(i) + ".bmp";
+        BMPWriter::saveRayTraceResults(rayFilename, grid, hitVoxels, testRays[i]);
+        printf("  Saved ray trace to '%s'\n", rayFilename.c_str());
+    }
+    
+    printf("\n=== Statistics ===\n");
+    printf("Total points generated: %zu\n", points.size());
+    printf("Voxels in grid: %zu\n", grid.getOccupiedPositions().size());
+    printf("Grid size: (%.1f, %.1f, %.1f)\n", gridSize.x, gridSize.y, gridSize.z);
+    printf("Voxel size: (%.1f, %.1f, %.1f)\n", voxelSize.x, voxelSize.y, voxelSize.z);
+
+    FunctionTimer::printStats(FunctionTimer::Mode::ENHANCED);
+
+    return 0;
+}
\ No newline at end of file
diff --git a/f.cpp b/f.cpp
new file mode 100644
index 0000000..3395af3
--- /dev/null
+++ b/f.cpp
@@ -0,0 +1,363 @@
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <algorithm>
+#include <unordered_map>
+#include <cstdio>
+#include <cstring>
+#include <fstream>
+#include <iostream>
+#include <vector>
+#include <random>
+#include <functional>
+#include <tuple>
+#include "util/timing_decorator.hpp"
+#include "util/vec3.hpp"
+#include "util/vec4.hpp"
+#include "util/bmpwriter.hpp"
+#include "util/voxelgrid.hpp"
+
+std::vector<Vec3> generateSphere(int numPoints, float radius = 1.0f, float wiggleAmount = 0.1f) {
+    
+    printf("Generating sphere with %d points using grid method...\n", numPoints);
+    printf("Wiggle amount: %.3f\n", wiggleAmount);
+    
+    std::vector<Vec3> points;
+    
+    // Random number generator for wiggling
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
+    
+    // Calculate grid resolution based on desired number of points
+    // For a sphere, we need to account that only about 52% of points in a cube will be inside the sphere
+    int gridRes = static_cast<int>(std::cbrt(numPoints / 0.52f)) + 1;
+    
+    printf("Using grid resolution: %d x %d x %d\n", gridRes, gridRes, gridRes);
+    
+    // Calculate voxel size for wiggling (based on average distance between points)
+    float voxelSize = 2.0f / gridRes;
+    float maxWiggle = wiggleAmount * voxelSize;
+    
+    printf("Voxel size: %.4f, Max wiggle: %.4f\n", voxelSize, maxWiggle);
+    
+    // Generate points in a cube from -1 to 1 in all dimensions
+    for (int x = 0; x < gridRes; ++x) {
+        for (int y = 0; y < gridRes; ++y) {
+            for (int z = 0; z < gridRes; ++z) {
+                // Convert grid coordinates to normalized cube coordinates [-1, 1]
+                Vec3 point(
+                    (2.0f * x / (gridRes - 1)) - 1.0f,
+                    (2.0f * y / (gridRes - 1)) - 1.0f,
+                    (2.0f * z / (gridRes - 1)) - 1.0f
+                );
+                
+                // Check if point is inside the unit sphere
+                if (point.lengthSquared() <= 1.0f) {
+                    // Apply randomized wiggling
+                    Vec3 wiggle(
+                        dist(gen) * maxWiggle,
+                        dist(gen) * maxWiggle,
+                        dist(gen) * maxWiggle
+                    );
+                    
+                    Vec3 wiggledPoint = point + wiggle;
+                    
+                    // Re-normalize to maintain spherical shape while preserving the wiggle
+                    // This ensures the point stays within the sphere while having natural variation
+                    float currentLength = wiggledPoint.length();
+                    if (currentLength > 1.0f) {
+                        // Scale back to unit sphere surface, but preserve the wiggle direction
+                        wiggledPoint = wiggledPoint * (1.0f / currentLength);
+                    }
+                    
+                    points.push_back(wiggledPoint * radius); // Scale by radius
+                }
+            }
+        }
+    }
+    
+    printf("Generated %zu points inside sphere\n", points.size());
+    
+    // If we have too many points, randomly sample down to the desired number
+    if (points.size() > static_cast<size_t>(numPoints)) {
+        printf("Sampling down from %zu to %d points...\n", points.size(), numPoints);
+        
+        // Shuffle and resize
+        std::shuffle(points.begin(), points.end(), gen);
+        points.resize(numPoints);
+    }
+    // If we have too few points, we'll use what we have
+    else if (points.size() < static_cast<size_t>(numPoints)) {
+        printf("Warning: Only generated %zu points (requested %d)\n", points.size(), numPoints);
+    }
+    
+    return points;
+}
+
+// Alternative sphere generation with perlin-like noise for more natural wiggling
+std::vector<Vec3> generateSphereWithNaturalWiggle(int numPoints, float radius = 1.0f, float noiseStrength = 0.15f) {
+    
+    printf("Generating sphere with natural wiggling using %d points...\n", numPoints);
+    printf("Noise strength: %.3f\n", noiseStrength);
+    
+    std::vector<Vec3> points;
+    
+    // Random number generators
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_real_distribution<float> dist(-1.0f, 1.0f);
+    
+    // Calculate grid resolution
+    int gridRes = static_cast<int>(std::cbrt(numPoints / 0.52f)) + 1;
+    printf("Using grid resolution: %d x %d x %d\n", gridRes, gridRes, gridRes);
+    
+    float voxelSize = 2.0f / gridRes;
+    float maxDisplacement = noiseStrength * voxelSize;
+    
+    // Pre-compute some random offsets for pseudo-perlin noise
+    std::vector<float> randomOffsets(gridRes * gridRes * gridRes);
+    for (size_t i = 0; i < randomOffsets.size(); ++i) {
+        randomOffsets[i] = dist(gen);
+    }
+    
+    auto getNoise = [&](int x, int y, int z) -> float {
+        int idx = (x * gridRes * gridRes) + (y * gridRes) + z;
+        return randomOffsets[idx % randomOffsets.size()];
+    };
+    
+    for (int x = 0; x < gridRes; ++x) {
+        for (int y = 0; y < gridRes; ++y) {
+            for (int z = 0; z < gridRes; ++z) {
+                Vec3 point(
+                    (2.0f * x / (gridRes - 1)) - 1.0f,
+                    (2.0f * y / (gridRes - 1)) - 1.0f,
+                    (2.0f * z / (gridRes - 1)) - 1.0f
+                );
+                
+                if (point.lengthSquared() <= 1.0f) {
+                    // Use smoother noise for more natural displacement
+                    float noiseX = getNoise(x, y, z);
+                    float noiseY = getNoise(y, z, x);
+                    float noiseZ = getNoise(z, x, y);
+                    
+                    // Apply displacement that preserves the overall spherical structure
+                    Vec3 displacement(
+                        noiseX * maxDisplacement,
+                        noiseY * maxDisplacement,
+                        noiseZ * maxDisplacement
+                    );
+                    
+                    Vec3 displacedPoint = point + displacement;
+                    
+                    // Ensure we don't push points too far from sphere surface
+                    float currentLength = displacedPoint.length();
+                    if (currentLength > 1.0f) {
+                        displacedPoint = displacedPoint * (0.95f + 0.05f * dist(gen)); // Small random scaling
+                    }
+                    
+                    points.push_back(displacedPoint * radius);
+                }
+            }
+        }
+    }
+    
+    printf("Generated %zu points with natural wiggling\n", points.size());
+    
+    // Sample down if needed
+    if (points.size() > static_cast<size_t>(numPoints)) {
+        printf("Sampling down from %zu to %d points...\n", points.size(), numPoints);
+        std::shuffle(points.begin(), points.end(), gen);
+        points.resize(numPoints);
+    } else if (points.size() < static_cast<size_t>(numPoints)) {
+        printf("Warning: Only generated %zu points (requested %d)\n", points.size(), numPoints);
+    }
+    
+    return points;
+}
+
+
+// Modified function to populate voxel grid with sphere points and assign layers
+void populateVoxelGridWithLayeredSphere(VoxelGrid& grid, const std::vector<Vec3>& points) {
+    printf("Populating voxel grid with %zu sphere points...\n", points.size());
+    
+    // First add all voxels with a default color
+    Vec4 defaultColor(1.0f, 1.0f, 1.0f, 1.0f); // Temporary color
+    for (const auto& point : points) {
+        grid.addVoxel(point, defaultColor);
+    }
+    
+    printf("Voxel grid populated with %zu voxels\n", grid.getOccupiedPositions().size());
+    
+    // Now assign the planetary layers
+    grid.assignPlanetaryLayers();
+}
+
+Vec3 rotateX(const Vec3& point, float angle) {
+    float cosA = std::cos(angle);
+    float sinA = std::sin(angle);
+    return Vec3(
+        point.x,
+        point.y * cosA - point.z * sinA,
+        point.y * sinA + point.z * cosA
+    );
+}
+
+Vec3 rotateY(const Vec3& point, float angle) {
+    float cosA = std::cos(angle);
+    float sinA = std::sin(angle);
+    return Vec3(
+        point.x * cosA + point.z * sinA,
+        point.y,
+        -point.x * sinA + point.z * cosA
+    );
+}
+
+Vec3 rotateZ(const Vec3& point, float angle) {
+    float cosA = std::cos(angle);
+    float sinA = std::sin(angle);
+    return Vec3(
+        point.x * cosA - point.y * sinA,
+        point.x * sinA + point.y * cosA,
+        point.z
+    );
+}
+
+void visualizePointCloud(const std::vector<Vec3>& points, const std::vector<Vec4>& colors, 
+                                 const std::string& filename, int width = 1000, int height = 1000,
+                                 float angleX = 0.0f, float angleY = 0.0f, float angleZ = 0.0f) {
+    TIME_FUNCTION;
+    std::vector<uint8_t> pixels(width * height * 3, 0);
+    
+    // Background color (dark blue)
+    for (int i = 0; i < width * height * 3; i += 3) {
+        pixels[i] = 30;     // B
+        pixels[i + 1] = 30; // G
+        pixels[i + 2] = 50; // R
+    }
+    
+    // Find bounds of point cloud
+    Vec3 minPoint( std::numeric_limits<float>::max(),  std::numeric_limits<float>::max(),  std::numeric_limits<float>::max());
+    Vec3 maxPoint(-std::numeric_limits<float>::max(), -std::numeric_limits<float>::max(), -std::numeric_limits<float>::max());
+    
+    // Apply rotation to all points and find new bounds
+    std::vector<Vec3> rotatedPoints;
+    rotatedPoints.reserve(points.size());
+    
+    for (const auto& point : points) {
+        Vec3 rotated = point;
+        rotated = rotateX(rotated, angleX);
+        rotated = rotateY(rotated, angleY);
+        rotated = rotateZ(rotated, angleZ);
+        rotatedPoints.push_back(rotated);
+        
+        minPoint.x = std::min(minPoint.x, rotated.x);
+        minPoint.y = std::min(minPoint.y, rotated.y);
+        minPoint.z = std::min(minPoint.z, rotated.z);
+        maxPoint.x = std::max(maxPoint.x, rotated.x);
+        maxPoint.y = std::max(maxPoint.y, rotated.y);
+        maxPoint.z = std::max(maxPoint.z, rotated.z);
+    }
+    
+    Vec3 cloudSize = maxPoint - minPoint;
+    float maxDim = std::max({cloudSize.x, cloudSize.y, cloudSize.z});
+    
+    // Draw points
+    for (size_t i = 0; i < rotatedPoints.size(); i++) {
+        const auto& point = rotatedPoints[i];
+        const auto& color = colors[i];
+        
+        // Map 3D point to 2D screen coordinates (orthographic projection)
+        int screenX = static_cast<int>(((point.x - minPoint.x) / maxDim) * (width - 20)) + 10;
+        int screenY = static_cast<int>(((point.y - minPoint.y) / maxDim) * (height - 20)) + 10;
+        
+        if (screenX >= 0 && screenX < width && screenY >= 0 && screenY < height) {
+            uint8_t r, g, b;
+            color.toUint8(r, g, b);
+            
+            // Draw a 2x2 pixel for each point
+            for (int dy = -1; dy <= 1; dy++) {
+                for (int dx = -1; dx <= 1; dx++) {
+                    int px = screenX + dx;
+                    int py = screenY + dy;
+                    if (px >= 0 && px < width && py >= 0 && py < height) {
+                        int index = (py * width + px) * 3;
+                        pixels[index] = b;
+                        pixels[index + 1] = g;
+                        pixels[index + 2] = r;
+                    }
+                }
+            }
+        }
+    }
+    
+    BMPWriter::saveBMP(filename, pixels, width, height);
+}
+
+int main() {
+    printf("=== Layered Sphere Generation and Visualization ===\n\n");
+    
+    const int numPoints = 100000000;
+    const float radius = 2.0f;
+    
+    printf("Generating layered spheres with %d points each, radius %.1f\n\n", numPoints, radius);
+    
+    // Create voxel grid
+    VoxelGrid grid(Vec3(10, 10, 10), Vec3(0.1f, 0.1f, 0.1f));
+    
+    // Generate sphere with different wiggle options:
+    
+    // Option 1: Simple random wiggling (small amount - 10% of voxel size)
+    printf("1. Generating sphere with simple wiggling...\n");
+    auto sphere1 = generateSphereWithNaturalWiggle(numPoints, radius, 0.1f);
+    
+    populateVoxelGridWithLayeredSphere(grid, sphere1);
+    
+    // Extract positions and colors for visualization
+    std::vector<Vec3> occupiedPositions = grid.getOccupiedPositions();
+    std::vector<Vec4> layerColors = grid.getColors();
+    
+    // Convert to world coordinates for visualization
+    std::vector<Vec3> worldPositions;
+    for (const auto& gridPos : occupiedPositions) {
+        worldPositions.push_back(grid.gridToWorld(gridPos));
+    }
+    
+    // Create multiple visualizations from different angles
+    printf("\nGenerating views from different angles...\n");
+    
+    // Front view (0, 0, 0)
+    visualizePointCloud(worldPositions, layerColors, "output/sphere_front.bmp", 1000, 1000, 0.0f, 0.0f, 0.0f);
+    printf("  - sphere_front.bmp (front view)\n");
+    
+    // 45 degree rotation around Y axis
+    visualizePointCloud(worldPositions, layerColors, "output/sphere_45y.bmp", 1000, 1000, 0.0f, M_PI/4, 0.0f);
+    printf("  - sphere_45y.bmp (45° Y rotation)\n");
+    
+    // 90 degree rotation around Y axis (side view)
+    visualizePointCloud(worldPositions, layerColors, "output/sphere_side.bmp", 1000, 1000, 0.0f, M_PI/2, 0.0f);
+    printf("  - sphere_side.bmp (side view)\n");
+    
+    // 45 degree rotation around X axis (top-down perspective)
+    visualizePointCloud(worldPositions, layerColors, "output/sphere_45x.bmp", 1000, 1000, M_PI/4, 0.0f, 0.0f);
+    printf("  - sphere_45x.bmp (45° X rotation)\n");
+    
+    // Combined rotation (30° X, 30° Y)
+    visualizePointCloud(worldPositions, layerColors, "output/sphere_30x_30y.bmp", 1000, 1000, M_PI/6, M_PI/6, 0.0f);
+    printf("  - sphere_30x_30y.bmp (30° X, 30° Y rotation)\n");
+    
+    // Top view (90° X rotation)
+    visualizePointCloud(worldPositions, layerColors, "output/sphere_top.bmp", 1000, 1000, M_PI/2, 0.0f, 0.0f);
+    printf("  - sphere_top.bmp (top view)\n");
+    
+    printf("\n=== Sphere generated successfully ===\n");
+    printf("Files created in output/ directory:\n");
+    printf("  - sphere_front.bmp     (Front view)\n");
+    printf("  - sphere_45y.bmp       (45° Y rotation)\n");
+    printf("  - sphere_side.bmp      (Side view)\n");
+    printf("  - sphere_45x.bmp       (45° X rotation)\n");
+    printf("  - sphere_30x_30y.bmp   (30° X, 30° Y rotation)\n");
+    printf("  - sphere_top.bmp       (Top view)\n");
+    
+    return 0;
+}
\ No newline at end of file
diff --git a/g.cpp b/g.cpp
new file mode 100644
index 0000000..fd6f221
--- /dev/null
+++ b/g.cpp
@@ -0,0 +1,17 @@
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <algorithm>
+#include <unordered_map>
+#include <cstdio>
+#include <cstring>
+#include <fstream>
+#include <iostream>
+#include <vector>
+#include <random>
+#include <functional>
+#include <tuple>
+#include "timing_decorator.hpp"
+#include "util/vec3.hpp"
+#include "util/vec4.hpp"
+#include "util/bmpwriter.hpp"
diff --git a/gradient.bmp b/gradient.bmp
new file mode 100644
index 0000000..dca969a
Binary files /dev/null and b/gradient.bmp differ
diff --git a/main.cpp b/main.cpp
new file mode 100644
index 0000000..12fb964
--- /dev/null
+++ b/main.cpp
@@ -0,0 +1,136 @@
+#include <iostream>
+#include <string>
+#include <vector>
+#include "util/grid2.hpp"
+#include "util/bmpwriter.hpp"
+
+// Function to convert hex color string to Vec4
+Vec4 hexToVec4(const std::string& hex) {
+    if (hex.length() != 6) {
+        return Vec4(0, 0, 0, 1); // Default to black if invalid
+    }
+    
+    int r, g, b;
+    sscanf(hex.c_str(), "%02x%02x%02x", &r, &g, &b);
+    
+    return Vec4(r / 255.0f, g / 255.0f, b / 255.0f, 1.0f);
+}
+
+int main(int argc, char* argv[]) {
+    // Check for gradient flag
+    bool createGradient = false;
+    for (int i = 1; i < argc; ++i) {
+        std::string arg = argv[i];
+        if (arg == "--gradient" || arg == "-g") {
+            createGradient = true;
+            break;
+        }
+    }
+    
+    if (!createGradient) {
+        std::cout << "Usage: " << argv[0] << " --gradient (-g)" << std::endl;
+        std::cout << "Creates a gradient image with red, green, and blue corners" << std::endl;
+        return 1;
+    }
+    
+    // Create a grid with points arranged in a gradient pattern
+    const int WIDTH = 512;
+    const int HEIGHT = 512;
+    const int POINTS_PER_DIM = 256; // Resolution of the gradient
+    
+    Grid2 grid;
+    
+    // Define our target colors at specific positions
+    Vec4 red = hexToVec4("ff0000");      // Top-left corner
+    Vec4 green = hexToVec4("00ff00");    // Center
+    Vec4 blue = hexToVec4("0000ff");     // Bottom-right corner
+    Vec4 white = hexToVec4("ffffff");    // Top-right corner
+    Vec4 black = hexToVec4("000000");    // Bottom-left corner
+    
+    // Create gradient points
+    for (int y = 0; y < POINTS_PER_DIM; ++y) {
+        for (int x = 0; x < POINTS_PER_DIM; ++x) {
+            // Normalize coordinates to [0, 1]
+            float nx = static_cast<float>(x) / (POINTS_PER_DIM - 1);
+            float ny = static_cast<float>(y) / (POINTS_PER_DIM - 1);
+            
+            // Create position in [-1, 1] range
+            Vec2 pos(nx * 2.0f - 1.0f, ny * 2.0f - 1.0f);
+            
+            // Calculate interpolated color based on position
+            Vec4 color;
+            
+            if (nx + ny <= 1.0f) {
+                // Lower triangle: interpolate between red, green, and black
+                if (nx <= 0.5f && ny <= 0.5f) {
+                    // Bottom-left quadrant: red to black to green
+                    float t1 = nx * 2.0f; // Horizontal interpolation
+                    float t2 = ny * 2.0f; // Vertical interpolation
+                    
+                    if (t1 + t2 <= 1.0f) {
+                        // Interpolate between red and black
+                        color = red * (1.0f - t1 - t2) + black * (t1 + t2);
+                    } else {
+                        // Interpolate between black and green
+                        color = black * (2.0f - t1 - t2) + green * (t1 + t2 - 1.0f);
+                    }
+                } else {
+                    // Use bilinear interpolation for other areas
+                    Vec4 topLeft = red;
+                    Vec4 topRight = white;
+                    Vec4 bottomLeft = black;
+                    Vec4 bottomRight = green;
+                    
+                    Vec4 top = topLeft * (1.0f - nx) + topRight * nx;
+                    Vec4 bottom = bottomLeft * (1.0f - nx) + bottomRight * nx;
+                    color = bottom * (1.0f - ny) + top * ny;
+                }
+            } else {
+                // Upper triangle: interpolate between green, blue, and white
+                if (nx >= 0.5f && ny >= 0.5f) {
+                    // Top-right quadrant: green to white to blue
+                    float t1 = (nx - 0.5f) * 2.0f; // Horizontal interpolation
+                    float t2 = (ny - 0.5f) * 2.0f; // Vertical interpolation
+                    
+                    if (t1 + t2 <= 1.0f) {
+                        // Interpolate between green and white
+                        color = green * (1.0f - t1 - t2) + white * (t1 + t2);
+                    } else {
+                        // Interpolate between white and blue
+                        color = white * (2.0f - t1 - t2) + blue * (t1 + t2 - 1.0f);
+                    }
+                } else {
+                    // Use bilinear interpolation for other areas
+                    Vec4 topLeft = red;
+                    Vec4 topRight = white;
+                    Vec4 bottomLeft = black;
+                    Vec4 bottomRight = blue;
+                    
+                    Vec4 top = topLeft * (1.0f - nx) + topRight * nx;
+                    Vec4 bottom = bottomLeft * (1.0f - nx) + bottomRight * nx;
+                    color = bottom * (1.0f - ny) + top * ny;
+                }
+            }
+            
+            grid.addPoint(pos, color);
+        }
+    }
+    
+    // Render to RGB image
+    std::vector<uint8_t> imageData = grid.renderToRGB(WIDTH, HEIGHT);
+    
+    // Save as BMP
+    if (BMPWriter::saveBMP("output/gradient.bmp", imageData, WIDTH, HEIGHT)) {
+        std::cout << "Gradient image saved as 'gradient.bmp'" << std::endl;
+        std::cout << "Colors: " << std::endl;
+        std::cout << "  Top-left: ff0000 (red)" << std::endl;
+        std::cout << "  Center: 00ff00 (green)" << std::endl;
+        std::cout << "  Bottom-right: 0000ff (blue)" << std::endl;
+        std::cout << "  Gradient between ffffff and 000000 throughout" << std::endl;
+    } else {
+        std::cerr << "Failed to save gradient image" << std::endl;
+        return 1;
+    }
+    
+    return 0;
+}
\ No newline at end of file
diff --git a/pointcloud_renderer b/pointcloud_renderer
new file mode 100644
index 0000000..86889f4
Binary files /dev/null and b/pointcloud_renderer differ
diff --git a/sim.cpp b/sim.cpp
new file mode 100644
index 0000000..a591532
--- /dev/null
+++ b/sim.cpp
@@ -0,0 +1,479 @@
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <thread>
+#include <chrono>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+#include <algorithm>
+#include <unordered_map>
+#include <functional>
+#include <map>
+
+#ifdef _WIN32
+#include <winsock2.h>
+#include <ws2tcpip.h>
+#pragma comment(lib, "ws2_32.lib")
+#else
+#include <sys/socket.h>
+#include <netinet/in.h>
+#include <unistd.h>
+#include <arpa/inet.h>
+#endif
+
+#define M_PI 3.14159265358979323846
+#define M_PHI 1.61803398874989484820458683436563811772030917980576286213544862270526046281890
+#define M_TAU 6.2831853071795864769252867665590057683943387987502116419498891846156328125724179972560696506842341359
+
+class Vec3 {
+public:
+    double x, y, z;
+    
+    Vec3(double x = 0, double y = 0, double z = 0) : x(x), y(y), z(z) {}
+    
+    inline double norm() const {
+        return std::sqrt(x*x + y*y + z*z);
+    }
+    
+    inline Vec3 normalize() const {
+        double n = norm();
+        return Vec3(x/n, y/n, z/n);
+    }
+    
+    inline Vec3 cross(const Vec3& other) const {
+        return Vec3(
+            y * other.z - z * other.y,
+            z * other.x - x * other.z,
+            x * other.y - y * other.x
+        );
+    }
+    
+    inline double dot(const Vec3& other) const {
+        return x * other.x + y * other.y + z * other.z;
+    }
+    
+    inline Vec3 operator+(const Vec3& other) const {
+        return Vec3(x + other.x, y + other.y, z + other.z);
+    }
+    
+    inline Vec3 operator-(const Vec3& other) const {
+        return Vec3(x - other.x, y - other.y, z - other.z);
+    }
+    
+    inline Vec3 operator*(double scalar) const {
+        return Vec3(x * scalar, y * scalar, z * scalar);
+    }
+    
+    inline Vec3 operator/(double scalar) const {
+        return Vec3(x / scalar, y / scalar, z / scalar);
+    }
+    
+    bool operator==(const Vec3& other) const {
+        return x == other.x && y == other.y && z == other.z;
+    }
+    
+    bool operator<(const Vec3& other) const {
+        if (x != other.x) return x < other.x;
+        if (y != other.y) return y < other.y;
+        return z < other.z;
+    }
+    
+    struct Hash {
+        size_t operator()(const Vec3& v) const {
+            size_t h1 = std::hash<double>()(std::round(v.x * 1000.0));
+            size_t h2 = std::hash<double>()(std::round(v.y * 1000.0));
+            size_t h3 = std::hash<double>()(std::round(v.z * 1000.0));
+            return h1 ^ (h2 << 1) ^ (h3 << 2);
+        }
+    };
+};
+
+class Triangle {
+public:
+    Vec3 v0, v1, v2;
+    
+    Triangle(const Vec3& v0, const Vec3& v1, const Vec3& v2) : v0(v0), v1(v1), v2(v2) {}
+    
+    Vec3 normal() const {
+        Vec3 edge1 = v1 - v0;
+        Vec3 edge2 = v2 - v0;
+        return edge1.cross(edge2).normalize();
+    }
+};
+
+std::vector<Vec3> fibsphere(int numPoints, float radius) {
+    std::vector<Vec3> points;
+    points.reserve(numPoints);
+        
+    for (int i = 0; i < numPoints; ++i) {
+        double y = 1.0 - (i / (double)(numPoints - 1)) * 2.0;
+        double radius_at_y = std::sqrt(1.0 - y * y);
+        
+        double theta = 2.0 * M_PI * i / M_PHI;
+        
+        double x = std::cos(theta) * radius_at_y;
+        double z = std::sin(theta) * radius_at_y;
+        
+        points.emplace_back(x * radius, y * radius, z * radius);
+    }
+    
+    return points;
+}
+
+// Create proper triangulation for Fibonacci sphere
+std::vector<Triangle> createFibonacciSphereMesh(const std::vector<Vec3>& points) {
+    std::vector<Triangle> triangles;
+    int n = points.size();
+    
+    // Create a map to quickly find points by their spherical coordinates
+    std::map<std::pair<double, double>, int> pointMap;
+    for (int i = 0; i < n; i++) {
+        double phi = std::acos(points[i].y / 3.0); // theta = acos(y/R)
+        double theta = std::atan2(points[i].z, points[i].x);
+        if (theta < 0) theta += 2.0 * M_PI;
+        pointMap[{phi, theta}] = i;
+    }
+    
+    // Create triangles by connecting neighboring points
+    for (int i = 0; i < n - 1; i++) {
+        // For each point, connect it to its neighbors
+        if (i > 0) {
+            triangles.emplace_back(points[i-1], points[i], points[(i+1) % n]);
+        }
+    }
+    
+    // Add cap triangles
+    for (int i = 1; i < n/2 - 1; i++) {
+        triangles.emplace_back(points[0], points[i], points[i+1]);
+        triangles.emplace_back(points[n-1], points[n-1-i], points[n-2-i]);
+    }
+    
+    return triangles;
+}
+
+// Alternative: Use Delaunay triangulation on the sphere (simplified)
+std::vector<Triangle> createSphereMeshDelaunay(const std::vector<Vec3>& points) {
+    std::vector<Triangle> triangles;
+    int n = points.size();
+    
+    // Simple approach: connect each point to its nearest neighbors
+    // This is a simplified version - for production use a proper Delaunay triangulation
+    
+    for (int i = 0; i < n; i++) {
+        // Find nearest neighbors (simplified)
+        std::vector<std::pair<double, int>> distances;
+        for (int j = 0; j < n; j++) {
+            if (i != j) {
+                double dist = (points[i] - points[j]).norm();
+                distances.emplace_back(dist, j);
+            }
+        }
+        
+        // Sort by distance and take 6 nearest neighbors
+        std::sort(distances.begin(), distances.end());
+        int numNeighbors = min(6, (int)distances.size());
+        
+        // Create triangles with nearest neighbors
+        for (int k = 0; k < numNeighbors - 1; k++) {
+            triangles.emplace_back(points[i], points[distances[k].second], points[distances[k+1].second]);
+        }
+    }
+    
+    return triangles;
+}
+
+Vec3 rotate(const Vec3& point, double angleX, double angleY, double angleZ) {
+    // Rotate around X axis
+    double y1 = point.y * cos(angleX) - point.z * sin(angleX);
+    double z1 = point.y * sin(angleX) + point.z * cos(angleX);
+    
+    // Rotate around Y axis
+    double x2 = point.x * cos(angleY) + z1 * sin(angleY);
+    double z2 = -point.x * sin(angleY) + z1 * cos(angleY);
+    
+    // Rotate around Z axis
+    double x3 = x2 * cos(angleZ) - y1 * sin(angleZ);
+    double y3 = x2 * sin(angleZ) + y1 * cos(angleZ);
+    
+    return Vec3(x3, y3, z2);
+}
+
+std::string generateSVG(const std::vector<Vec3>& points, const std::vector<Triangle>& mesh, 
+                       double angleX, double angleY, double angleZ) {
+    std::stringstream svg;
+    int width = 800;
+    int height = 600;
+    
+    svg << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n";
+    svg << "<svg width=\"" << width << "\" height=\"" << height << "\" xmlns=\"http://www.w3.org/2000/svg\">\n";
+    svg << "<rect width=\"100%\" height=\"100%\" fill=\"#000000\"/>\n";
+    
+    // Project 3D to 2D
+    auto project = [&](const Vec3& point) -> std::pair<double, double> {
+        Vec3 rotated = rotate(point, angleX, angleY, angleZ);
+        // Perspective projection
+        double scale = 300.0 / (5.0 + rotated.z);
+        double x = width / 2 + rotated.x * scale;
+        double y = height / 2 + rotated.y * scale;
+        return {x, y};
+    };
+    
+    // Draw triangles with shading
+    for (const auto& triangle : mesh) {
+        Vec3 normal = triangle.normal();
+        Vec3 lightDir = Vec3(0.5, 0.7, 1.0).normalize();
+        double intensity = max(0.0, normal.dot(lightDir));
+        
+        // Calculate color based on intensity
+        int r = static_cast<int>(50 + intensity * 200);
+        int g = static_cast<int>(100 + intensity * 150);
+        int b = static_cast<int>(200 + intensity * 55);
+        
+        auto [x0, y0] = project(triangle.v0);
+        auto [x1, y1] = project(triangle.v1);
+        auto [x2, y2] = project(triangle.v2);
+        
+        // Only draw triangles facing the camera
+        Vec3 viewDir(0, 0, 1);
+        if (normal.dot(viewDir) > 0.1) {
+            svg << "<polygon points=\"" 
+                << x0 << "," << y0 << " "
+                << x1 << "," << y1 << " "
+                << x2 << "," << y2
+                << "\" fill=\"rgb(" << r << "," << g << "," << b << ")\" "
+                << "stroke=\"rgba(0,0,0,0.3)\" stroke-width=\"1\"/>\n";
+        }
+    }
+    
+    // Draw points for debugging
+    for (const auto& point : points) {
+        auto [x, y] = project(point);
+        svg << "<circle cx=\"" << x << "\" cy=\"" << y << "\" r=\"2\" fill=\"white\"/>\n";
+    }
+    
+    svg << "</svg>";
+    return svg.str();
+}
+
+// HTTP server class (keep your existing server implementation)
+class SimpleHTTPServer {
+private:
+    int serverSocket;
+    int port;
+    
+public:
+    SimpleHTTPServer(int port) : port(port), serverSocket(-1) {}
+    
+    ~SimpleHTTPServer() {
+        stop();
+    }
+    
+    bool start() {
+#ifdef _WIN32
+        WSADATA wsaData;
+        if (WSAStartup(MAKEWORD(2, 2), &wsaData) != 0) {
+            std::cerr << "WSAStartup failed" << std::endl;
+            return false;
+        }
+#endif
+        
+        serverSocket = socket(AF_INET, SOCK_STREAM, 0);
+        if (serverSocket < 0) {
+            std::cerr << "Socket creation failed" << std::endl;
+            return false;
+        }
+        
+        int opt = 1;
+#ifdef _WIN32
+        if (setsockopt(serverSocket, SOL_SOCKET, SO_REUSEADDR, (char*)&opt, sizeof(opt)) < 0) {
+#else
+        if (setsockopt(serverSocket, SOL_SOCKET, SO_REUSEADDR, &opt, sizeof(opt)) < 0) {
+#endif
+            std::cerr << "Setsockopt failed" << std::endl;
+            return false;
+        }
+        
+        sockaddr_in serverAddr;
+        serverAddr.sin_family = AF_INET;
+        serverAddr.sin_addr.s_addr = INADDR_ANY;
+        serverAddr.sin_port = htons(port);
+        
+        if (bind(serverSocket, (sockaddr*)&serverAddr, sizeof(serverAddr)) < 0) {
+            std::cerr << "Bind failed" << std::endl;
+            return false;
+        }
+        
+        if (listen(serverSocket, 10) < 0) {
+            std::cerr << "Listen failed" << std::endl;
+            return false;
+        }
+        
+        std::cout << "Server started on port " << port << std::endl;
+        return true;
+    }
+    
+    void stop() {
+        if (serverSocket >= 0) {
+#ifdef _WIN32
+            closesocket(serverSocket);
+            WSACleanup();
+#else
+            close(serverSocket);
+#endif
+            serverSocket = -1;
+        }
+    }
+    
+    void handleRequests() {
+        // Generate proper Fibonacci sphere
+        std::vector<Vec3> spherePoints = fibsphere(200, 3.0); // Reduced for performance
+        std::vector<Triangle> sphereMesh = createSphereMeshDelaunay(spherePoints);
+        
+        while (true) {
+            sockaddr_in clientAddr;
+#ifdef _WIN32
+            int clientAddrLen = sizeof(clientAddr);
+#else
+            socklen_t clientAddrLen = sizeof(clientAddr);
+#endif
+            int clientSocket = accept(serverSocket, (sockaddr*)&clientAddr, &clientAddrLen);
+            
+            if (clientSocket < 0) {
+                std::cerr << "Accept failed" << std::endl;
+                continue;
+            }
+            
+            char buffer[4096] = {0};
+            recv(clientSocket, buffer, sizeof(buffer), 0);
+            
+            std::string request(buffer);
+            std::string response;
+            
+            if (request.find("GET / ") != std::string::npos || request.find("GET /index.html") != std::string::npos) {
+                response = "HTTP/1.1 200 OK\r\nContent-Type: text/html\r\n\r\n" + getHTML();
+            } else if (request.find("GET /mesh.svg") != std::string::npos) {
+                static double angle = 0.0;
+                angle += 0.02;
+                std::string svg = generateSVG(sphereMesh, angle, angle * 0.7, angle * 0.3);
+                response = "HTTP/1.1 200 OK\r\nContent-Type: image/svg+xml\r\n\r\n" + svg;
+            } else {
+                response = "HTTP/1.1 404 Not Found\r\nContent-Type: text/plain\r\n\r\n404 Not Found";
+            }
+            
+            send(clientSocket, response.c_str(), response.length(), 0);
+            
+#ifdef _WIN32
+            closesocket(clientSocket);
+#else
+            close(clientSocket);
+#endif
+            
+            std::this_thread::sleep_for(std::chrono::milliseconds(10));
+        }
+    }
+    
+    std::string getHTML() {
+        return R"(
+<!DOCTYPE html>
+<html>
+<head>
+    <title>3D Sphere Mesh Renderer</title>
+    <style>
+        body {
+            font-family: 'Arial', sans-serif;
+            margin: 0;
+            padding: 20px;
+            background: linear-gradient(135deg, #667eea 0%, #764ba2 100%);
+            color: white;
+            text-align: center;
+            min-height: 100vh;
+            display: flex;
+            flex-direction: column;
+            align-items: center;
+            justify-content: center;
+        }
+        .container {
+            max-width: 900px;
+            background: rgba(255, 255, 255, 0.1);
+            padding: 30px;
+            border-radius: 15px;
+            backdrop-filter: blur(10px);
+            box-shadow: 0 8px 32px rgba(0, 0, 0, 0.3);
+        }
+        h1 {
+            margin-bottom: 20px;
+            text-shadow: 2px 2px 4px rgba(0, 0, 0, 0.3);
+            font-size: 2.5em;
+        }
+        #sphereContainer {
+            margin: 20px 0;
+            display: flex;
+            justify-content: center;
+        }
+        .instructions {
+            margin-top: 20px;
+            padding: 20px;
+            background: rgba(255, 255, 255, 0.2);
+            border-radius: 10px;
+            text-align: left;
+        }
+        .features {
+            display: grid;
+            grid-template-columns: repeat(auto-fit, minmax(250px, 1fr));
+            gap: 15px;
+            margin-top: 20px;
+        }
+        .feature {
+            background: rgba(255, 255, 255, 0.15);
+            padding: 15px;
+            border-radius: 8px;
+        }
+        .footer {
+            margin-top: 30px;
+            font-size: 0.9em;
+            opacity: 0.8;
+        }
+    </style>
+</head>
+<body>
+    <div class="container">
+        <h1>3D Sphere Mesh Renderer</h1>
+        
+        <div id="sphereContainer">
+            <img id="sphereImage" src="mesh.svg" width="600" height="450" alt="3D Sphere">
+        </div>
+        <--- DO NOT EDIT --->
+        <div class="footer">
+            <p>Built with C++ HTTP Server and SVG Graphics</p>
+        </div>
+    </div>
+
+    <script>
+        setInterval(function() {
+            const img = document.getElementById('sphereImage');
+            const timestamp = new Date().getTime();
+            img.src = 'mesh.svg?' + timestamp;
+        }, 50);
+    </script>
+</body>
+</html>
+)";
+    }
+};
+
+int main() {
+    SimpleHTTPServer server(5101);
+    
+    if (!server.start()) {
+        std::cerr << "Failed to start server" << std::endl;
+        return 1;
+    }
+    
+    std::cout << "Open your browser and navigate to http://localhost:5101" << std::endl;
+    std::cout << "Press Ctrl+C to stop the server" << std::endl;
+    
+    server.handleRequests();
+    
+    return 0;
+}
\ No newline at end of file
diff --git a/util/Vec2.hpp b/util/Vec2.hpp
new file mode 100644
index 0000000..35c05a2
--- /dev/null
+++ b/util/Vec2.hpp
@@ -0,0 +1,286 @@
+#ifndef VEC2_HPP
+#define VEC2_HPP
+
+#include <cmath>
+#include <algorithm>
+#include <string>
+
+class Vec2 {
+    public:
+    float x, y;
+    Vec2() : x(0), y(0) {}
+    Vec2(float x, float y) : x(x), y(y) {}
+    
+    Vec2 operator+(const Vec2& other) const {
+        return Vec2(x + other.x, y + other.y);
+    }
+    
+    Vec2 operator-(const Vec2& other) const {
+        return Vec2(x - other.x, y - other.y);
+    }
+    
+    Vec2 operator*(const Vec2& other) const {
+        return Vec2(x * other.x, y * other.y);
+    }
+    
+    Vec2 operator/(const Vec2& other) const {
+        return Vec2(x / other.x, y / other.y);
+    }
+
+    Vec2 operator+(float scalar) const {
+        return Vec2(x + scalar, y + scalar);
+    }
+
+    Vec2 operator-(float scalar) const {
+        return Vec2(x - scalar, y - scalar);
+    }
+    
+    Vec2 operator-() const {
+        return Vec2(-x, -y);
+    }
+
+    Vec2 operator*(float scalar) const {
+        return Vec2(x * scalar, y * scalar);
+    }
+
+    Vec2 operator/(float scalar) const {
+        return Vec2(x / scalar, y / scalar);
+    }
+
+    Vec2& operator=(float scalar) {
+        x = y = scalar;
+        return *this;
+    }
+    
+    Vec2& operator+=(const Vec2& other) {
+        x += other.x;
+        y += other.y;
+        return *this;
+    }
+    
+    Vec2& operator-=(const Vec2& other) {
+        x -= other.x;
+        y -= other.y;
+        return *this;
+    }
+    
+    Vec2& operator*=(const Vec2& other) {
+        x *= other.x;
+        y *= other.y;
+        return *this;
+    }
+    
+    Vec2& operator/=(const Vec2& other) {
+        x /= other.x;
+        y /= other.y;
+        return *this;
+    }
+    
+    Vec2& operator+=(float scalar) {
+        x += scalar;
+        y += scalar;
+        return *this;
+    }
+    
+    Vec2& operator-=(float scalar) {
+        x -= scalar;
+        y -= scalar;
+        return *this;
+    }
+    
+    Vec2& operator*=(float scalar) {
+        x *= scalar;
+        y *= scalar;
+        return *this;
+    }
+    
+    Vec2& operator/=(float scalar) {
+        x /= scalar;
+        y /= scalar;
+        return *this;
+    }
+
+    float dot(const Vec2& other) const {
+        return x * other.x + y * other.y;
+    }
+
+    float length() const {
+        return std::sqrt(x * x + y * y);
+    }
+    
+    float lengthSquared() const {
+        return x * x + y * y;
+    }
+    
+    float distance(const Vec2& other) const {
+        return (*this - other).length();
+    }
+    
+    float distanceSquared(const Vec2& other) const {
+        Vec2 diff = *this - other;
+        return diff.x * diff.x + diff.y * diff.y;
+    }
+
+    Vec2 normalized() const {
+        float len = length();
+        if (len > 0) {
+            return *this / len;
+        }
+        return *this;
+    }
+    
+    bool operator==(const Vec2& other) const {
+        return x == other.x && y == other.y;
+    }
+    
+    bool operator!=(const Vec2& other) const {
+        return x != other.x || y != other.y;
+    }
+    
+    bool operator<(const Vec2& other) const {
+        return (x < other.x) || (x == other.x && y < other.y);
+    }
+    
+    bool operator<=(const Vec2& other) const {
+        return (x < other.x) || (x == other.x && y <= other.y);
+    }
+    
+    bool operator>(const Vec2& other) const {
+        return (x > other.x) || (x == other.x && y > other.y);
+    }
+    
+    bool operator>=(const Vec2& other) const {
+        return (x > other.x) || (x == other.x && y >= other.y);
+    }
+    
+    Vec2 abs() const {
+        return Vec2(std::abs(x), std::abs(y));
+    }
+    
+    Vec2 floor() const {
+        return Vec2(std::floor(x), std::floor(y));
+    }
+    
+    Vec2 ceil() const {
+        return Vec2(std::ceil(x), std::ceil(y));
+    }
+    
+    Vec2 round() const {
+        return Vec2(std::round(x), std::round(y));
+    }
+    
+    Vec2 min(const Vec2& other) const {
+        return Vec2(std::min(x, other.x), std::min(y, other.y));
+    }
+    
+    Vec2 max(const Vec2& other) const {
+        return Vec2(std::max(x, other.x), std::max(y, other.y));
+    }
+    
+    Vec2 clamp(const Vec2& minVal, const Vec2& maxVal) const {
+        return Vec2(
+            std::clamp(x, minVal.x, maxVal.x),
+            std::clamp(y, minVal.y, maxVal.y)
+        );
+    }
+    
+    Vec2 clamp(float minVal, float maxVal) const {
+        return Vec2(
+            std::clamp(x, minVal, maxVal),
+            std::clamp(y, minVal, maxVal)
+        );
+    }
+    
+    bool isZero(float epsilon = 1e-10f) const {
+        return std::abs(x) < epsilon && std::abs(y) < epsilon;
+    }
+    
+    bool equals(const Vec2& other, float epsilon = 1e-10f) const {
+        return std::abs(x - other.x) < epsilon && 
+               std::abs(y - other.y) < epsilon;
+    }
+    
+    friend Vec2 operator+(float scalar, const Vec2& vec) {
+        return Vec2(scalar + vec.x, scalar + vec.y);
+    }
+    
+    friend Vec2 operator-(float scalar, const Vec2& vec) {
+        return Vec2(scalar - vec.x, scalar - vec.y);
+    }
+    
+    friend Vec2 operator*(float scalar, const Vec2& vec) {
+        return Vec2(scalar * vec.x, scalar * vec.y);
+    }
+    
+    friend Vec2 operator/(float scalar, const Vec2& vec) {
+        return Vec2(scalar / vec.x, scalar / vec.y);
+    }
+
+    Vec2 perpendicular() const {
+        return Vec2(-y, x);
+    }
+    
+    Vec2 reflect(const Vec2& normal) const {
+        return *this - 2.0f * this->dot(normal) * normal;
+    }
+    
+    Vec2 lerp(const Vec2& other, float t) const {
+        t = std::clamp(t, 0.0f, 1.0f);
+        return *this + (other - *this) * t;
+    }
+    
+    Vec2 slerp(const Vec2& other, float t) const {
+        t = std::clamp(t, 0.0f, 1.0f);
+        float dot = this->dot(other);
+        dot = std::clamp(dot, -1.0f, 1.0f);
+        
+        float theta = std::acos(dot) * t;
+        Vec2 relative = other - *this * dot;
+        relative = relative.normalized();
+        
+        return (*this * std::cos(theta)) + (relative * std::sin(theta));
+    }
+    
+    Vec2 rotate(float angle) const {
+        float cosA = std::cos(angle);
+        float sinA = std::sin(angle);
+        return Vec2(x * cosA - y * sinA, x * sinA + y * cosA);
+    }
+    
+    float angle() const {
+        return std::atan2(y, x);
+    }
+    
+    float angleTo(const Vec2& other) const {
+        return std::acos(this->dot(other) / (this->length() * other.length()));
+    }
+    
+    float& operator[](int index) {
+        return (&x)[index];
+    }
+    
+    const float& operator[](int index) const {
+        return (&x)[index];
+    }
+    
+    std::string toString() const {
+        return "(" + std::to_string(x) + ", " + std::to_string(y) + ")";
+    }
+    
+};
+
+inline std::ostream& operator<<(std::ostream& os, const Vec2& vec) {
+    os << vec.toString();
+    return os;
+}
+
+namespace std {
+    template<>
+    struct hash<Vec2> {
+        size_t operator()(const Vec2& v) const {
+            return hash<float>()(v.x) ^ (hash<float>()(v.y) << 1);
+        }
+    };
+}
+
+#endif
\ No newline at end of file
diff --git a/util/bmpwriter.hpp b/util/bmpwriter.hpp
new file mode 100644
index 0000000..95a1830
--- /dev/null
+++ b/util/bmpwriter.hpp
@@ -0,0 +1,167 @@
+#ifndef BMP_WRITER_HPP
+#define BMP_WRITER_HPP
+
+#include <vector>
+#include <fstream>
+#include <cstring>
+#include <string>
+#include <algorithm>
+#include "vec3.hpp"
+
+class BMPWriter {
+private:
+    #pragma pack(push, 1)
+    struct BMPHeader {
+        uint16_t signature = 0x4D42; // "BM"
+        uint32_t fileSize;
+        uint16_t reserved1 = 0;
+        uint16_t reserved2 = 0;
+        uint32_t dataOffset = 54;
+    };
+    
+    struct BMPInfoHeader {
+        uint32_t headerSize = 40;
+        int32_t width;
+        int32_t height;
+        uint16_t planes = 1;
+        uint16_t bitsPerPixel = 24;
+        uint32_t compression = 0;
+        uint32_t imageSize;
+        int32_t xPixelsPerMeter = 0;
+        int32_t yPixelsPerMeter = 0;
+        uint32_t colorsUsed = 0;
+        uint32_t importantColors = 0;
+    };
+    #pragma pack(pop)
+
+public:
+    // Save a 2D vector of Vec3 (RGB) colors as BMP
+    // Vec3 components: x = red, y = green, z = blue (values in range [0,1])
+    static bool saveBMP(const std::string& filename, const std::vector<std::vector<Vec3>>& pixels) {
+        if (pixels.empty() || pixels[0].empty()) {
+            return false;
+        }
+        
+        int height = static_cast<int>(pixels.size());
+        int width = static_cast<int>(pixels[0].size());
+        
+        // Validate that all rows have the same width
+        for (const auto& row : pixels) {
+            if (row.size() != width) {
+                return false;
+            }
+        }
+        
+        return saveBMP(filename, pixels, width, height);
+    }
+    
+    // Alternative interface with width/height and flat vector (row-major order)
+    static bool saveBMP(const std::string& filename, const std::vector<Vec3>& pixels, int width, int height) {
+        if (pixels.size() != width * height) {
+            return false;
+        }
+        
+        // Convert to 2D vector format
+        std::vector<std::vector<Vec3>> pixels2D(height, std::vector<Vec3>(width));
+        for (int y = 0; y < height; ++y) {
+            for (int x = 0; x < width; ++x) {
+                pixels2D[y][x] = pixels[y * width + x];
+            }
+        }
+        
+        return saveBMP(filename, pixels2D, width, height);
+    }
+    
+    // Save from 1D vector of uint8_t pixels (BGR order: pixels[i]=b, pixels[i+1]=g, pixels[i+2]=r)
+    static bool saveBMP(const std::string& filename, const std::vector<uint8_t>& pixels, int width, int height) {
+        if (pixels.size() != width * height * 3) {
+            return false;
+        }
+        
+        BMPHeader header;
+        BMPInfoHeader infoHeader;
+        
+        int rowSize = (width * 3 + 3) & ~3; // 24-bit, padded to 4 bytes
+        int imageSize = rowSize * height;
+        
+        header.fileSize = sizeof(BMPHeader) + sizeof(BMPInfoHeader) + imageSize;
+        infoHeader.width = width;
+        infoHeader.height = height;
+        infoHeader.imageSize = imageSize;
+        
+        std::ofstream file(filename, std::ios::binary);
+        if (!file) {
+            return false;
+        }
+        
+        file.write(reinterpret_cast<const char*>(&header), sizeof(header));
+        file.write(reinterpret_cast<const char*>(&infoHeader), sizeof(infoHeader));
+        
+        // Write pixel data (BMP stores pixels bottom-to-top)
+        std::vector<uint8_t> row(rowSize, 0);
+        for (int y = height - 1; y >= 0; --y) {
+            const uint8_t* srcRow = pixels.data() + (y * width * 3);
+            
+            // Copy and rearrange if necessary (input is already in BGR order)
+            for (int x = 0; x < width; ++x) {
+                int srcOffset = x * 3;
+                int dstOffset = x * 3;
+                
+                // Input is already BGR: pixels[i]=b, pixels[i+1]=g, pixels[i+2]=r
+                // So we can copy directly
+                row[dstOffset] = srcRow[srcOffset];     // B
+                row[dstOffset + 1] = srcRow[srcOffset + 1]; // G
+                row[dstOffset + 2] = srcRow[srcOffset + 2]; // R
+            }
+            file.write(reinterpret_cast<const char*>(row.data()), rowSize);
+        }
+        
+        return true;
+    }
+    
+private:
+    static bool saveBMP(const std::string& filename, const std::vector<std::vector<Vec3>>& pixels, int width, int height) {
+        BMPHeader header;
+        BMPInfoHeader infoHeader;
+        
+        int rowSize = (width * 3 + 3) & ~3; // 24-bit, padded to 4 bytes
+        int imageSize = rowSize * height;
+        
+        header.fileSize = sizeof(BMPHeader) + sizeof(BMPInfoHeader) + imageSize;
+        infoHeader.width = width;
+        infoHeader.height = height;
+        infoHeader.imageSize = imageSize;
+        
+        std::ofstream file(filename, std::ios::binary);
+        if (!file) {
+            return false;
+        }
+        
+        file.write(reinterpret_cast<const char*>(&header), sizeof(header));
+        file.write(reinterpret_cast<const char*>(&infoHeader), sizeof(infoHeader));
+        
+        // Write pixel data (BMP stores pixels bottom-to-top)
+        std::vector<uint8_t> row(rowSize, 0);
+        for (int y = height - 1; y >= 0; --y) {
+            for (int x = 0; x < width; ++x) {
+                const Vec3& color = pixels[y][x];
+                
+                // Convert from [0,1] float to [0,255] uint8_t
+                uint8_t r = static_cast<uint8_t>(std::clamp(color.x * 255.0f, 0.0f, 255.0f));
+                uint8_t g = static_cast<uint8_t>(std::clamp(color.y * 255.0f, 0.0f, 255.0f));
+                uint8_t b = static_cast<uint8_t>(std::clamp(color.z * 255.0f, 0.0f, 255.0f));
+                
+                // BMP is BGR order
+                int pixelOffset = x * 3;
+                row[pixelOffset] = b;
+                row[pixelOffset + 1] = g;
+                row[pixelOffset + 2] = r;
+            }
+            file.write(reinterpret_cast<const char*>(row.data()), rowSize);
+        }
+        
+        return true;
+    }
+};
+
+#endif
\ No newline at end of file
diff --git a/util/grid2.hpp b/util/grid2.hpp
new file mode 100644
index 0000000..70684b4
--- /dev/null
+++ b/util/grid2.hpp
@@ -0,0 +1,231 @@
+#ifndef GRID2_HPP
+#define GRID2_HPP
+
+#include "Vec2.hpp"
+#include "Vec4.hpp"
+#include <vector>
+#include <cstdint>
+#include <algorithm>
+#include <stdexcept>
+
+class Grid2 {
+public:
+    std::vector<Vec2> positions;
+    std::vector<Vec4> colors;
+    
+    Grid2() = default;
+    
+    // Constructor with initial size
+    Grid2(size_t size) {
+        positions.resize(size);
+        colors.resize(size);
+    }
+    
+    // Add a point with position and color
+    void addPoint(const Vec2& position, const Vec4& color) {
+        positions.push_back(position);
+        colors.push_back(color);
+    }
+    
+    // Clear all points
+    void clear() {
+        positions.clear();
+        colors.clear();
+    }
+    
+    // Get number of points
+    size_t size() const {
+        return positions.size();
+    }
+    
+    // Check if grid is empty
+    bool empty() const {
+        return positions.empty();
+    }
+    
+    // Resize the grid
+    void resize(size_t newSize) {
+        positions.resize(newSize);
+        colors.resize(newSize);
+    }
+    
+    // Render to RGB image data
+    std::vector<uint8_t> renderToRGB(int width, int height, const Vec4& backgroundColor = Vec4(0, 0, 0, 1)) const {
+        if (width <= 0 || height <= 0) {
+            throw std::invalid_argument("Width and height must be positive");
+        }
+        
+        std::vector<uint8_t> imageData(width * height * 3);
+        
+        // Initialize with background color
+        uint8_t bgR, bgG, bgB;
+        backgroundColor.toUint8(bgR, bgG, bgB);
+        
+        for (int i = 0; i < width * height * 3; i += 3) {
+            imageData[i] = bgR;
+            imageData[i + 1] = bgG;
+            imageData[i + 2] = bgB;
+        }
+        
+        // Find the bounding box of all points to map to pixel coordinates
+        if (positions.empty()) {
+            return imageData;
+        }
+        
+        Vec2 minPos = positions[0];
+        Vec2 maxPos = positions[0];
+        
+        for (const auto& pos : positions) {
+            minPos = minPos.min(pos);
+            maxPos = maxPos.max(pos);
+        }
+        
+        // Add a small margin to avoid division by zero and edge issues
+        Vec2 size = maxPos - minPos;
+        if (size.x < 1e-10f) size.x = 1.0f;
+        if (size.y < 1e-10f) size.y = 1.0f;
+        
+        float margin = 0.05f; // 5% margin
+        minPos -= size * margin;
+        maxPos += size * margin;
+        size = maxPos - minPos;
+        
+        // Render each point
+        for (size_t i = 0; i < positions.size(); i++) {
+            const Vec2& pos = positions[i];
+            const Vec4& color = colors[i];
+            
+            // Convert world coordinates to pixel coordinates
+            float normalizedX = (pos.x - minPos.x) / size.x;
+            float normalizedY = 1.0f - (pos.y - minPos.y) / size.y; // Flip Y for image coordinates
+            
+            int pixelX = static_cast<int>(normalizedX * width);
+            int pixelY = static_cast<int>(normalizedY * height);
+            
+            // Clamp to image bounds
+            pixelX = std::clamp(pixelX, 0, width - 1);
+            pixelY = std::clamp(pixelY, 0, height - 1);
+            
+            // Convert color to RGB
+            uint8_t r, g, b;
+            color.toUint8(r, g, b);
+            
+            // Set pixel color
+            int index = (pixelY * width + pixelX) * 3;
+            imageData[index] = r;
+            imageData[index + 1] = g;
+            imageData[index + 2] = b;
+        }
+        
+        return imageData;
+    }
+    
+    // Render to RGBA image data (with alpha channel)
+    std::vector<uint8_t> renderToRGBA(int width, int height, const Vec4& backgroundColor = Vec4(0, 0, 0, 1)) const {
+        if (width <= 0 || height <= 0) {
+            throw std::invalid_argument("Width and height must be positive");
+        }
+        
+        std::vector<uint8_t> imageData(width * height * 4);
+        
+        // Initialize with background color
+        uint8_t bgR, bgG, bgB, bgA;
+        backgroundColor.toUint8(bgR, bgG, bgB, bgA);
+        
+        for (int i = 0; i < width * height * 4; i += 4) {
+            imageData[i] = bgR;
+            imageData[i + 1] = bgG;
+            imageData[i + 2] = bgB;
+            imageData[i + 3] = bgA;
+        }
+        
+        if (positions.empty()) {
+            return imageData;
+        }
+        
+        // Find the bounding box (same as RGB version)
+        Vec2 minPos = positions[0];
+        Vec2 maxPos = positions[0];
+        
+        for (const auto& pos : positions) {
+            minPos = minPos.min(pos);
+            maxPos = maxPos.max(pos);
+        }
+        
+        Vec2 size = maxPos - minPos;
+        if (size.x < 1e-10f) size.x = 1.0f;
+        if (size.y < 1e-10f) size.y = 1.0f;
+        
+        float margin = 0.05f;
+        minPos -= size * margin;
+        maxPos += size * margin;
+        size = maxPos - minPos;
+        
+        // Render each point
+        for (size_t i = 0; i < positions.size(); i++) {
+            const Vec2& pos = positions[i];
+            const Vec4& color = colors[i];
+            
+            float normalizedX = (pos.x - minPos.x) / size.x;
+            float normalizedY = 1.0f - (pos.y - minPos.y) / size.y;
+            
+            int pixelX = static_cast<int>(normalizedX * width);
+            int pixelY = static_cast<int>(normalizedY * height);
+            
+            pixelX = std::clamp(pixelX, 0, width - 1);
+            pixelY = std::clamp(pixelY, 0, height - 1);
+            
+            uint8_t r, g, b, a;
+            color.toUint8(r, g, b, a);
+            
+            int index = (pixelY * width + pixelX) * 4;
+            imageData[index] = r;
+            imageData[index + 1] = g;
+            imageData[index + 2] = b;
+            imageData[index + 3] = a;
+        }
+        
+        return imageData;
+    }
+    
+    // Get the bounding box of all positions
+    void getBoundingBox(Vec2& minPos, Vec2& maxPos) const {
+        if (positions.empty()) {
+            minPos = Vec2(0, 0);
+            maxPos = Vec2(0, 0);
+            return;
+        }
+        
+        minPos = positions[0];
+        maxPos = positions[0];
+        
+        for (const auto& pos : positions) {
+            minPos = minPos.min(pos);
+            maxPos = maxPos.max(pos);
+        }
+    }
+    
+    // Scale all positions to fit within a specified range
+    void normalizePositions(const Vec2& targetMin = Vec2(-1, -1), const Vec2& targetMax = Vec2(1, 1)) {
+        if (positions.empty()) return;
+        
+        Vec2 currentMin, currentMax;
+        getBoundingBox(currentMin, currentMax);
+        
+        Vec2 currentSize = currentMax - currentMin;
+        Vec2 targetSize = targetMax - targetMin;
+        
+        if (currentSize.x < 1e-10f) currentSize.x = 1.0f;
+        if (currentSize.y < 1e-10f) currentSize.y = 1.0f;
+        
+        for (auto& pos : positions) {
+            float normalizedX = (pos.x - currentMin.x) / currentSize.x;
+            float normalizedY = (pos.y - currentMin.y) / currentSize.y;
+            
+            pos.x = targetMin.x + normalizedX * targetSize.x;
+            pos.y = targetMin.y + normalizedY * targetSize.y;
+        }
+    }
+};
+
+#endif
\ No newline at end of file
diff --git a/util/ray.cpp b/util/ray.cpp
new file mode 100644
index 0000000..f50f0f8
--- /dev/null
+++ b/util/ray.cpp
@@ -0,0 +1,27 @@
+#ifndef RAY_HPP
+#define RAY_HPP
+
+#include "Vec2.hpp"
+#include "Vec3.hpp"
+#include "Vec4.hpp"
+#include "ray2.hpp"
+#include "ray3.hpp"
+#include "ray4.hpp"
+
+// Stream operators for rays
+inline std::ostream& operator<<(std::ostream& os, const Ray2& ray) {
+    os << ray.toString();
+    return os;
+}
+
+inline std::ostream& operator<<(std::ostream& os, const Ray3& ray) {
+    os << ray.toString();
+    return os;
+}
+
+inline std::ostream& operator<<(std::ostream& os, const Ray4& ray) {
+    os << ray.toString();
+    return os;
+}
+
+#endif
\ No newline at end of file
diff --git a/util/ray2.hpp b/util/ray2.hpp
new file mode 100644
index 0000000..645662d
--- /dev/null
+++ b/util/ray2.hpp
@@ -0,0 +1,59 @@
+#ifndef RAY2_HPP
+#define RAY2_HPP
+
+#include "Vec2.hpp"
+
+class Ray2 {
+public:
+    Vec2 origin;
+    Vec2 direction;
+    
+    Ray2() : origin(Vec2()), direction(Vec2(1, 0)) {}
+    Ray2(const Vec2& origin, const Vec2& direction) 
+        : origin(origin), direction(direction.normalized()) {}
+    
+    // Get point at parameter t along the ray
+    Vec2 at(float t) const {
+        return origin + direction * t;
+    }
+    
+    // Reflect ray off a surface with given normal
+    Ray2 reflect(const Vec2& point, const Vec2& normal) const {
+        Vec2 reflectedDir = direction.reflect(normal);
+        return Ray2(point, reflectedDir);
+    }
+    
+    // Check if ray intersects with a circle
+    bool intersectsCircle(const Vec2& center, float radius, float& t1, float& t2) const {
+        Vec2 oc = origin - center;
+        float a = direction.dot(direction);
+        float b = 2.0f * oc.dot(direction);
+        float c = oc.dot(oc) - radius * radius;
+        
+        float discriminant = b * b - 4 * a * c;
+        
+        if (discriminant < 0) {
+            return false;
+        }
+        
+        discriminant = std::sqrt(discriminant);
+        t1 = (-b - discriminant) / (2.0f * a);
+        t2 = (-b + discriminant) / (2.0f * a);
+        
+        return true;
+    }
+    
+    // Get the distance from a point to this ray
+    float distanceToPoint(const Vec2& point) const {
+        Vec2 pointToOrigin = point - origin;
+        float projection = pointToOrigin.dot(direction);
+        Vec2 closestPoint = origin + direction * projection;
+        return point.distance(closestPoint);
+    }
+    
+    std::string toString() const {
+        return "Ray2(origin: " + origin.toString() + ", direction: " + direction.toString() + ")";
+    }
+};
+
+#endif
\ No newline at end of file
diff --git a/util/ray3.hpp b/util/ray3.hpp
new file mode 100644
index 0000000..8b8cd67
--- /dev/null
+++ b/util/ray3.hpp
@@ -0,0 +1,73 @@
+#ifndef RAY3_HPP
+#define RAY3_HPP
+
+#include "Vec3.hpp"
+
+class Ray3 {
+public:
+    Vec3 origin;
+    Vec3 direction;
+    
+    Ray3() : origin(Vec3()), direction(Vec3(1, 0, 0)) {}
+    Ray3(const Vec3& origin, const Vec3& direction) 
+        : origin(origin), direction(direction.normalized()) {}
+    
+    // Get point at parameter t along the ray
+    Vec3 at(float t) const {
+        return origin + direction * t;
+    }
+    
+    // Reflect ray off a surface with given normal
+    Ray3 reflect(const Vec3& point, const Vec3& normal) const {
+        Vec3 reflectedDir = direction.reflect(normal);
+        return Ray3(point, reflectedDir);
+    }
+    
+    // Check if ray intersects with a sphere
+    bool intersectsSphere(const Vec3& center, float radius, float& t1, float& t2) const {
+        Vec3 oc = origin - center;
+        float a = direction.dot(direction);
+        float b = 2.0f * oc.dot(direction);
+        float c = oc.dot(oc) - radius * radius;
+        
+        float discriminant = b * b - 4 * a * c;
+        
+        if (discriminant < 0) {
+            return false;
+        }
+        
+        discriminant = std::sqrt(discriminant);
+        t1 = (-b - discriminant) / (2.0f * a);
+        t2 = (-b + discriminant) / (2.0f * a);
+        
+        return true;
+    }
+    
+    // Check if ray intersects with a plane (defined by point and normal)
+    bool intersectsPlane(const Vec3& planePoint, const Vec3& planeNormal, float& t) const {
+        float denom = planeNormal.dot(direction);
+        
+        if (std::abs(denom) < 1e-10f) {
+            return false; // Ray is parallel to plane
+        }
+        
+        t = planeNormal.dot(planePoint - origin) / denom;
+        return t >= 0;
+    }
+    
+    // Get the distance from a point to this ray
+    float distanceToPoint(const Vec3& point) const {
+        Vec3 pointToOrigin = point - origin;
+        Vec3 crossProduct = direction.cross(pointToOrigin);
+        return crossProduct.length() / direction.length();
+    }
+    
+    // Transform ray by a 4x4 matrix (for perspective/affine transformations)
+    Ray3 transform(const class Mat4& matrix) const;
+    
+    std::string toString() const {
+        return "Ray3(origin: " + origin.toString() + ", direction: " + direction.toString() + ")";
+    }
+};
+
+#endif
\ No newline at end of file
diff --git a/util/ray4.hpp b/util/ray4.hpp
new file mode 100644
index 0000000..646b4dc
--- /dev/null
+++ b/util/ray4.hpp
@@ -0,0 +1,52 @@
+#ifndef RAY4_HPP
+#define RAY4_HPP
+
+#include "vec4.hpp"
+
+class Ray4 {
+public:
+    Vec4 origin;
+    Vec4 direction;
+    
+    Ray4() : origin(Vec4()), direction(Vec4(1, 0, 0, 0)) {}
+    Ray4(const Vec4& origin, const Vec4& direction) 
+        : origin(origin), direction(direction.normalized()) {}
+    
+    // Get point at parameter t along the ray (in 4D space)
+    Vec4 at(float t) const {
+        return origin + direction * t;
+    }
+    
+    // Get 3D projection of the ray (homogeneous coordinates)
+    // Ray3 projectTo3D() const {
+    //     Vec3 projOrigin = origin.homogenized().xyz();
+    //     Vec3 projDirection = direction.homogenized().xyz().normalized();
+    //     return Ray3(projOrigin, projDirection);
+    // }
+    
+    // Get the distance from a point to this ray in 4D space
+    float distanceToPoint(const Vec4& point) const {
+        Vec4 pointToOrigin = point - origin;
+        float projection = pointToOrigin.dot(direction);
+        Vec4 closestPoint = origin + direction * projection;
+        return point.distance(closestPoint);
+    }
+    
+    // Check if this 4D ray intersects with a 3D hyperplane
+    bool intersectsHyperplane(const Vec4& planePoint, const Vec4& planeNormal, float& t) const {
+        float denom = planeNormal.dot(direction);
+        
+        if (std::abs(denom) < 1e-10f) {
+            return false; // Ray is parallel to hyperplane
+        }
+        
+        t = planeNormal.dot(planePoint - origin) / denom;
+        return true;
+    }
+    
+    std::string toString() const {
+        return "Ray4(origin: " + origin.toString() + ", direction: " + direction.toString() + ")";
+    }
+};
+
+#endif
\ No newline at end of file
diff --git a/util/timing_decorator.cpp b/util/timing_decorator.cpp
new file mode 100644
index 0000000..b8ef2af
--- /dev/null
+++ b/util/timing_decorator.cpp
@@ -0,0 +1,115 @@
+#include "timing_decorator.hpp"
+#include <cmath>
+
+std::unordered_map<std::string, FunctionTimer::TimingStats> FunctionTimer::stats_;
+
+void FunctionTimer::recordTiming(const std::string& func_name, double elapsed_seconds) {
+    auto& stat = stats_[func_name];
+    stat.call_count++;
+    stat.total_time += elapsed_seconds;
+    stat.timings.push_back(elapsed_seconds);
+}
+
+FunctionTimer::PercentileStats FunctionTimer::calculatePercentiles(const std::vector<double>& timings) {
+    PercentileStats result;
+    if (timings.empty()) return result;
+    
+    std::vector<double> sorted_timings = timings;
+    std::sort(sorted_timings.begin(), sorted_timings.end());
+    
+    auto percentile_index = [&](double p) -> size_t {
+        return std::min(sorted_timings.size() - 1, 
+                       static_cast<size_t>(p * sorted_timings.size() / 100.0));
+    };
+    
+    result.min = sorted_timings.front();
+    result.max = sorted_timings.back();
+    result.median = sorted_timings[percentile_index(50.0)];
+    result.p90 = sorted_timings[percentile_index(90.0)];
+    result.p95 = sorted_timings[percentile_index(95.0)];
+    result.p99 = sorted_timings[percentile_index(99.0)];
+    result.p99_9 = sorted_timings[percentile_index(99.9)];
+    
+    return result;
+}
+
+std::unordered_map<std::string, FunctionTimer::TimingStats> FunctionTimer::getStats(Mode mode) {
+    return stats_; // In C++ we return all data, filtering happens in print
+}
+
+void FunctionTimer::printStats(Mode mode) {
+    if (stats_.empty()) {
+        std::cout << "No timing statistics available." << std::endl;
+        return;
+    }
+
+    // Determine column widths
+    size_t func_col_width = 0;
+    for (const auto& [name, _] : stats_) {
+        func_col_width = std::max(func_col_width, name.length());
+    }
+    func_col_width = std::max(func_col_width, size_t(8)); // "Function"
+    
+    const int num_width = 12;
+
+    if (mode == Mode::BASIC) {
+        std::cout << "\nBasic Function Timing Statistics:" << std::endl;
+        std::cout << std::string(func_col_width + 3 * num_width + 8, '-') << std::endl;
+        
+        std::cout << std::left << std::setw(func_col_width) << "Function"
+                  << std::setw(num_width) << "Calls"
+                  << std::setw(num_width) << "Total (s)"
+                  << std::setw(num_width) << "Avg (s)" << std::endl;
+                  
+        std::cout << std::string(func_col_width + 3 * num_width + 8, '-') << std::endl;
+        
+        for (const auto& [func_name, data] : stats_) {
+            if (data.call_count > 0) {
+                std::cout << std::left << std::setw(func_col_width) << func_name
+                          << std::setw(num_width) << data.call_count
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << data.total_time
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << data.avg_time() << std::endl;
+            }
+        }
+        
+        std::cout << std::string(func_col_width + 3 * num_width + 8, '-') << std::endl;
+    } else { // ENHANCED mode
+        std::cout << "\nEnhanced Function Timing Statistics:" << std::endl;
+        size_t total_width = func_col_width + 8 * num_width + 8;
+        std::cout << std::string(total_width, '-') << std::endl;
+        
+        std::cout << std::left << std::setw(func_col_width) << "Function"
+                  << std::setw(num_width) << "Calls"
+                  << std::setw(num_width) << "Total (s)"
+                  << std::setw(num_width) << "Avg (s)"
+                  << std::setw(num_width) << "Min (s)"
+                  << std::setw(num_width) << "Median (s)"
+                  << std::setw(num_width) << "P99 (s)"
+                  << std::setw(num_width) << "P99.9 (s)"
+                  << std::setw(num_width) << "Max (s)" << std::endl;
+                  
+        std::cout << std::string(total_width, '-') << std::endl;
+        
+        for (const auto& [func_name, data] : stats_) {
+            if (data.call_count > 0) {
+                auto percentiles = calculatePercentiles(data.timings);
+                
+                std::cout << std::left << std::setw(func_col_width) << func_name
+                          << std::setw(num_width) << data.call_count
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << data.total_time
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << data.avg_time()
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << percentiles.min
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << percentiles.median
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << percentiles.p99
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << percentiles.p99_9
+                          << std::setw(num_width) << std::fixed << std::setprecision(6) << percentiles.max << std::endl;
+            }
+        }
+        
+        std::cout << std::string(total_width, '-') << std::endl;
+    }
+}
+
+void FunctionTimer::clearStats() {
+    stats_.clear();
+}
\ No newline at end of file
diff --git a/util/timing_decorator.hpp b/util/timing_decorator.hpp
new file mode 100644
index 0000000..218ae86
--- /dev/null
+++ b/util/timing_decorator.hpp
@@ -0,0 +1,100 @@
+#pragma once
+
+#include <chrono>
+#include <unordered_map>
+#include <string>
+#include <vector>
+#include <algorithm>
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+#include <memory>
+#include <functional>
+
+class FunctionTimer {
+public:
+    enum class Mode { BASIC, ENHANCED };
+    
+    struct TimingStats {
+        size_t call_count = 0;
+        double total_time = 0.0;
+        std::vector<double> timings;
+        
+        double avg_time() const {
+            return call_count > 0 ? total_time / call_count : 0.0;
+        }
+    };
+    
+    struct PercentileStats {
+        double p99_9 = 0.0;
+        double p99 = 0.0;
+        double p95 = 0.0;
+        double p90 = 0.0;
+        double max = 0.0;
+        double min = 0.0;
+        double median = 0.0;
+    };
+    
+    // Record timing for a function
+    static void recordTiming(const std::string& func_name, double elapsed_seconds);
+    
+    // Get statistics
+    static std::unordered_map<std::string, TimingStats> getStats(Mode mode = Mode::BASIC);
+    
+    // Print statistics
+    static void printStats(Mode mode = Mode::ENHANCED);
+    
+    // Clear all statistics
+    static void clearStats();
+
+private:
+    static std::unordered_map<std::string, TimingStats> stats_;
+    
+    static PercentileStats calculatePercentiles(const std::vector<double>& timings);
+};
+
+// Macro to easily time functions - similar to Python decorator
+#define TIME_FUNCTION auto function_timer_scoped_ = ScopedFunctionTimer(__func__)
+
+// Scoped timer for RAII-style timing
+class ScopedFunctionTimer {
+public:
+    ScopedFunctionTimer(const std::string& func_name) 
+        : func_name_(func_name), start_(std::chrono::steady_clock::now()) {}
+    
+    ~ScopedFunctionTimer() {
+        auto end = std::chrono::steady_clock::now();
+        auto elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end - start_);
+        FunctionTimer::recordTiming(func_name_, elapsed.count() / 1000000.0);
+    }
+
+private:
+    std::string func_name_;
+    std::chrono::steady_clock::time_point start_;
+};
+
+// Template decorator for functions (similar to Python)
+template<typename Func, typename... Args>
+auto time_function_decorator(const std::string& func_name, Func&& func, Args&&... args) {
+    auto start = std::chrono::steady_clock::now();
+    
+    if constexpr (std::is_void_v<std::invoke_result_t<Func, Args...>>) {
+        // Void return type
+        std::invoke(std::forward<Func>(func), std::forward<Args>(args)...);
+        auto end = std::chrono::steady_clock::now();
+        auto elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
+        FunctionTimer::recordTiming(func_name, elapsed.count() / 1000000.0);
+    } else {
+        // Non-void return type
+        auto result = std::invoke(std::forward<Func>(func), std::forward<Args>(args)...);
+        auto end = std::chrono::steady_clock::now();
+        auto elapsed = std::chrono::duration_cast<std::chrono::microseconds>(end - start);
+        FunctionTimer::recordTiming(func_name, elapsed.count() / 1000000.0);
+        return result;
+    }
+}
+
+// Macro to create decorated functions
+#define DECORATE_FUNCTION(func) [&](auto&&... args) { \
+    return time_function_decorator(#func, func, std::forward<decltype(args)>(args)...); \
+}
diff --git a/util/vec.cpp b/util/vec.cpp
new file mode 100644
index 0000000..4df7aaf
--- /dev/null
+++ b/util/vec.cpp
@@ -0,0 +1,19 @@
+#ifndef vec_hpp
+#define vec_hpp
+
+#include "Vec4.hpp"
+#include "Vec3.hpp"
+#include "Vec2.hpp"
+
+Vec4::Vec4(const Vec3& vec3, float w) : x(vec3.x), y(vec3.y), z(vec3.z), w(w) {}
+Vec3::Vec3(const Vec2& vec2, float z) : x(vec2.x), y(vec2.y), z(z) {}
+
+Vec3 Vec4::xyz() const {
+    return Vec3(x, y, z);
+}
+
+Vec3 Vec4::rgb() const {
+    return Vec3(r, g, b);
+}
+
+#endif
\ No newline at end of file
diff --git a/util/vec3.hpp b/util/vec3.hpp
new file mode 100644
index 0000000..d0d6020
--- /dev/null
+++ b/util/vec3.hpp
@@ -0,0 +1,322 @@
+#ifndef VEC3_HPP
+#define VEC3_HPP
+
+#include <cmath>
+#include <algorithm>
+#include <string>
+#include <ostream>
+
+class Vec3 {
+public:
+    float x, y, z;
+    
+    Vec3() : x(0), y(0), z(0) {}
+    Vec3(float x, float y, float z) : x(x), y(y), z(z) {}
+    Vec3(float scalar) : x(scalar), y(scalar), z(scalar) {}
+
+    Vec3(const class Vec2& vec2, float z = 0.0f);
+    
+    // Arithmetic operations
+    Vec3 operator+(const Vec3& other) const {
+        return Vec3(x + other.x, y + other.y, z + other.z);
+    }
+    
+    Vec3 operator-(const Vec3& other) const {
+        return Vec3(x - other.x, y - other.y, z - other.z);
+    }
+    
+    Vec3 operator*(const Vec3& other) const {
+        return Vec3(x * other.x, y * other.y, z * other.z);
+    }
+    
+    Vec3 operator/(const Vec3& other) const {
+        return Vec3(x / other.x, y / other.y, z / other.z);
+    }
+
+    Vec3 operator+(float scalar) const {
+        return Vec3(x + scalar, y + scalar, z + scalar);
+    }
+
+    Vec3 operator-(float scalar) const {
+        return Vec3(x - scalar, y - scalar, z - scalar);
+    }
+    
+    Vec3 operator-() const {
+        return Vec3(-x, -y, -z);
+    }
+
+    Vec3 operator*(float scalar) const {
+        return Vec3(x * scalar, y * scalar, z * scalar);
+    }
+
+    Vec3 operator/(float scalar) const {
+        return Vec3(x / scalar, y / scalar, z / scalar);
+    }
+
+    Vec3& operator=(float scalar) {
+        x = y = z = scalar;
+        return *this;
+    }
+    
+    Vec3& operator+=(const Vec3& other) {
+        x += other.x;
+        y += other.y;
+        z += other.z;
+        return *this;
+    }
+    
+    Vec3& operator-=(const Vec3& other) {
+        x -= other.x;
+        y -= other.y;
+        z -= other.z;
+        return *this;
+    }
+    
+    Vec3& operator*=(const Vec3& other) {
+        x *= other.x;
+        y *= other.y;
+        z *= other.z;
+        return *this;
+    }
+    
+    Vec3& operator/=(const Vec3& other) {
+        x /= other.x;
+        y /= other.y;
+        z /= other.z;
+        return *this;
+    }
+    
+    Vec3& operator+=(float scalar) {
+        x += scalar;
+        y += scalar;
+        z += scalar;
+        return *this;
+    }
+    
+    Vec3& operator-=(float scalar) {
+        x -= scalar;
+        y -= scalar;
+        z -= scalar;
+        return *this;
+    }
+    
+    Vec3& operator*=(float scalar) {
+        x *= scalar;
+        y *= scalar;
+        z *= scalar;
+        return *this;
+    }
+    
+    Vec3& operator/=(float scalar) {
+        x /= scalar;
+        y /= scalar;
+        z /= scalar;
+        return *this;
+    }
+
+    float dot(const Vec3& other) const {
+        return x * other.x + y * other.y + z * other.z;
+    }
+    
+    Vec3 cross(const Vec3& other) const {
+        return Vec3(
+            y * other.z - z * other.y,
+            z * other.x - x * other.z,
+            x * other.y - y * other.x
+        );
+    }
+
+    float length() const {
+        return std::sqrt(x * x + y * y + z * z);
+    }
+    
+    float lengthSquared() const {
+        return x * x + y * y + z * z;
+    }
+    
+    float distance(const Vec3& other) const {
+        return (*this - other).length();
+    }
+    
+    float distanceSquared(const Vec3& other) const {
+        Vec3 diff = *this - other;
+        return diff.x * diff.x + diff.y * diff.y + diff.z * diff.z;
+    }
+
+    Vec3 normalized() const {
+        float len = length();
+        if (len > 0) {
+            return *this / len;
+        }
+        return *this;
+    }
+    
+    bool operator==(const Vec3& other) const {
+        return x == other.x && y == other.y && z == other.z;
+    }
+    
+    bool operator!=(const Vec3& other) const {
+        return x != other.x || y != other.y || z != other.z;
+    }
+    
+    bool operator<(const Vec3& other) const {
+        return (x < other.x) || 
+               (x == other.x && y < other.y) || 
+               (x == other.x && y == other.y && z < other.z);
+    }
+    
+    bool operator<=(const Vec3& other) const {
+        return (x < other.x) || 
+               (x == other.x && y < other.y) || 
+               (x == other.x && y == other.y && z <= other.z);
+    }
+    
+    bool operator>(const Vec3& other) const {
+        return (x > other.x) || 
+               (x == other.x && y > other.y) || 
+               (x == other.x && y == other.y && z > other.z);
+    }
+    
+    bool operator>=(const Vec3& other) const {
+        return (x > other.x) || 
+               (x == other.x && y > other.y) || 
+               (x == other.x && y == other.y && z >= other.z);
+    }
+    
+    Vec3 abs() const {
+        return Vec3(std::abs(x), std::abs(y), std::abs(z));
+    }
+    
+    Vec3 floor() const {
+        return Vec3(std::floor(x), std::floor(y), std::floor(z));
+    }
+    
+    Vec3 ceil() const {
+        return Vec3(std::ceil(x), std::ceil(y), std::ceil(z));
+    }
+    
+    Vec3 round() const {
+        return Vec3(std::round(x), std::round(y), std::round(z));
+    }
+    
+    Vec3 min(const Vec3& other) const {
+        return Vec3(std::min(x, other.x), std::min(y, other.y), std::min(z, other.z));
+    }
+    
+    Vec3 max(const Vec3& other) const {
+        return Vec3(std::max(x, other.x), std::max(y, other.y), std::max(z, other.z));
+    }
+    
+    Vec3 clamp(const Vec3& minVal, const Vec3& maxVal) const {
+        return Vec3(
+            std::clamp(x, minVal.x, maxVal.x),
+            std::clamp(y, minVal.y, maxVal.y),
+            std::clamp(z, minVal.z, maxVal.z)
+        );
+    }
+    
+    Vec3 clamp(float minVal, float maxVal) const {
+        return Vec3(
+            std::clamp(x, minVal, maxVal),
+            std::clamp(y, minVal, maxVal),
+            std::clamp(z, minVal, maxVal)
+        );
+    }
+    
+    bool isZero(float epsilon = 1e-10f) const {
+        return std::abs(x) < epsilon && std::abs(y) < epsilon && std::abs(z) < epsilon;
+    }
+    
+    bool equals(const Vec3& other, float epsilon = 1e-10f) const {
+        return std::abs(x - other.x) < epsilon && 
+               std::abs(y - other.y) < epsilon &&
+               std::abs(z - other.z) < epsilon;
+    }
+    
+    friend Vec3 operator+(float scalar, const Vec3& vec) {
+        return Vec3(scalar + vec.x, scalar + vec.y, scalar + vec.z);
+    }
+    
+    friend Vec3 operator-(float scalar, const Vec3& vec) {
+        return Vec3(scalar - vec.x, scalar - vec.y, scalar - vec.z);
+    }
+    
+    friend Vec3 operator*(float scalar, const Vec3& vec) {
+        return Vec3(scalar * vec.x, scalar * vec.y, scalar * vec.z);
+    }
+    
+    friend Vec3 operator/(float scalar, const Vec3& vec) {
+        return Vec3(scalar / vec.x, scalar / vec.y, scalar / vec.z);
+    }
+
+    Vec3 reflect(const Vec3& normal) const {
+        return *this - 2.0f * this->dot(normal) * normal;
+    }
+    
+    Vec3 lerp(const Vec3& other, float t) const {
+        t = std::clamp(t, 0.0f, 1.0f);
+        return *this + (other - *this) * t;
+    }
+    
+    Vec3 slerp(const Vec3& other, float t) const {
+        t = std::clamp(t, 0.0f, 1.0f);
+        float dot = this->dot(other);
+        dot = std::clamp(dot, -1.0f, 1.0f);
+        
+        float theta = std::acos(dot) * t;
+        Vec3 relative = other - *this * dot;
+        relative = relative.normalized();
+        
+        return (*this * std::cos(theta)) + (relative * std::sin(theta));
+    }
+    
+    Vec3 rotateX(float angle) const {
+        float cosA = std::cos(angle);
+        float sinA = std::sin(angle);
+        return Vec3(x, y * cosA - z * sinA, y * sinA + z * cosA);
+    }
+    
+    Vec3 rotateY(float angle) const {
+        float cosA = std::cos(angle);
+        float sinA = std::sin(angle);
+        return Vec3(x * cosA + z * sinA, y, -x * sinA + z * cosA);
+    }
+    
+    Vec3 rotateZ(float angle) const {
+        float cosA = std::cos(angle);
+        float sinA = std::sin(angle);
+        return Vec3(x * cosA - y * sinA, x * sinA + y * cosA, z);
+    }
+    
+    float angleTo(const Vec3& other) const {
+        return std::acos(this->dot(other) / (this->length() * other.length()));
+    }
+    
+    float& operator[](int index) {
+        return (&x)[index];
+    }
+    
+    const float& operator[](int index) const {
+        return (&x)[index];
+    }
+    
+    std::string toString() const {
+        return "(" + std::to_string(x) + ", " + std::to_string(y) + ", " + std::to_string(z) + ")";
+    }
+};
+
+inline std::ostream& operator<<(std::ostream& os, const Vec3& vec) {
+    os << vec.toString();
+    return os;
+}
+
+namespace std {
+    template<>
+    struct hash<Vec3> {
+        size_t operator()(const Vec3& v) const {
+            return hash<float>()(v.x) ^ (hash<float>()(v.y) << 1) ^ (hash<float>()(v.z) << 2);
+        }
+    };
+}
+
+#endif
\ No newline at end of file
diff --git a/util/vec4.hpp b/util/vec4.hpp
new file mode 100644
index 0000000..0fb9c1c
--- /dev/null
+++ b/util/vec4.hpp
@@ -0,0 +1,384 @@
+#ifndef VEC4_HPP
+#define VEC4_HPP
+
+#include "vec3.hpp"
+#include <cmath>
+#include <algorithm>
+#include <string>
+#include <ostream>
+#include <cstdint>
+
+class Vec4 {
+public:
+    union {
+        struct { float x, y, z, w; };
+        struct { float r, g, b, a; };
+        struct { float s, t, p, q; }; // For texture coordinates
+    };
+    
+    // Constructors
+    Vec4() : x(0), y(0), z(0), w(0) {}
+    Vec4(float x, float y, float z, float w) : x(x), y(y), z(z), w(w) {}
+    Vec4(float scalar) : x(scalar), y(scalar), z(scalar), w(scalar) {}
+    
+    Vec4(const Vec3& rgb, float w = 1.0f) : x(rgb.x), y(rgb.y), z(rgb.z), w(w) {}
+    static Vec4 RGB(float r, float g, float b, float a = 1.0f) { return Vec4(r, g, b, a); }
+    static Vec4 RGBA(float r, float g, float b, float a) { return Vec4(r, g, b, a); }
+    
+    
+    Vec4 operator+(const Vec4& other) const {
+        return Vec4(x + other.x, y + other.y, z + other.z, w + other.w);
+    }
+    
+    Vec4 operator-(const Vec4& other) const {
+        return Vec4(x - other.x, y - other.y, z - other.z, w - other.w);
+    }
+    
+    Vec4 operator*(const Vec4& other) const {
+        return Vec4(x * other.x, y * other.y, z * other.z, w * other.w);
+    }
+    
+    Vec4 operator/(const Vec4& other) const {
+        return Vec4(x / other.x, y / other.y, z / other.z, w / other.w);
+    }
+
+    Vec4 operator+(float scalar) const {
+        return Vec4(x + scalar, y + scalar, z + scalar, w + scalar);
+    }
+
+    Vec4 operator-(float scalar) const {
+        return Vec4(x - scalar, y - scalar, z - scalar, w - scalar);
+    }
+    
+    Vec4 operator-() const {
+        return Vec4(-x, -y, -z, -w);
+    }
+
+    Vec4 operator*(float scalar) const {
+        return Vec4(x * scalar, y * scalar, z * scalar, w * scalar);
+    }
+
+    Vec4 operator/(float scalar) const {
+        return Vec4(x / scalar, y / scalar, z / scalar, w / scalar);
+    }
+
+    Vec4& operator=(float scalar) {
+        x = y = z = w = scalar;
+        return *this;
+    }
+    
+    Vec4& operator+=(const Vec4& other) {
+        x += other.x;
+        y += other.y;
+        z += other.z;
+        w += other.w;
+        return *this;
+    }
+    
+    Vec4& operator-=(const Vec4& other) {
+        x -= other.x;
+        y -= other.y;
+        z -= other.z;
+        w -= other.w;
+        return *this;
+    }
+    
+    Vec4& operator*=(const Vec4& other) {
+        x *= other.x;
+        y *= other.y;
+        z *= other.z;
+        w *= other.w;
+        return *this;
+    }
+    
+    Vec4& operator/=(const Vec4& other) {
+        x /= other.x;
+        y /= other.y;
+        z /= other.z;
+        w /= other.w;
+        return *this;
+    }
+    
+    Vec4& operator+=(float scalar) {
+        x += scalar;
+        y += scalar;
+        z += scalar;
+        w += scalar;
+        return *this;
+    }
+    
+    Vec4& operator-=(float scalar) {
+        x -= scalar;
+        y -= scalar;
+        z -= scalar;
+        w -= scalar;
+        return *this;
+    }
+    
+    Vec4& operator*=(float scalar) {
+        x *= scalar;
+        y *= scalar;
+        z *= scalar;
+        w *= scalar;
+        return *this;
+    }
+    
+    Vec4& operator/=(float scalar) {
+        x /= scalar;
+        y /= scalar;
+        z /= scalar;
+        w /= scalar;
+        return *this;
+    }
+
+    float dot(const Vec4& other) const {
+        return x * other.x + y * other.y + z * other.z + w * other.w;
+    }
+    
+    // 4D cross product (returns vector perpendicular to 3 given vectors in 4D space)
+    Vec4 cross(const Vec4& v1, const Vec4& v2, const Vec4& v3) const {
+        float a = v1.y * (v2.z * v3.w - v2.w * v3.z) -
+                  v1.z * (v2.y * v3.w - v2.w * v3.y) +
+                  v1.w * (v2.y * v3.z - v2.z * v3.y);
+                  
+        float b = -v1.x * (v2.z * v3.w - v2.w * v3.z) +
+                   v1.z * (v2.x * v3.w - v2.w * v3.x) -
+                   v1.w * (v2.x * v3.z - v2.z * v3.x);
+                   
+        float c = v1.x * (v2.y * v3.w - v2.w * v3.y) -
+                  v1.y * (v2.x * v3.w - v2.w * v3.x) +
+                  v1.w * (v2.x * v3.y - v2.y * v3.x);
+                  
+        float d = -v1.x * (v2.y * v3.z - v2.z * v3.y) +
+                   v1.y * (v2.x * v3.z - v2.z * v3.x) -
+                   v1.z * (v2.x * v3.y - v2.y * v3.x);
+        
+        return Vec4(a, b, c, d);
+    }
+
+    float length() const {
+        return std::sqrt(x * x + y * y + z * z + w * w);
+    }
+    
+    float lengthSquared() const {
+        return x * x + y * y + z * z + w * w;
+    }
+    
+    float distance(const Vec4& other) const {
+        return (*this - other).length();
+    }
+    
+    float distanceSquared(const Vec4& other) const {
+        Vec4 diff = *this - other;
+        return diff.x * diff.x + diff.y * diff.y + diff.z * diff.z + diff.w * diff.w;
+    }
+
+    Vec4 normalized() const {
+        float len = length();
+        if (len > 0) {
+            return *this / len;
+        }
+        return *this;
+    }
+    
+    // Homogeneous normalization (divide by w)
+    Vec4 homogenized() const {
+        if (w != 0.0f) {
+            return Vec4(x / w, y / w, z / w, 1.0f);
+        }
+        return *this;
+    }
+    
+    // Clamp values between 0 and 1
+    Vec4 clamped() const {
+        return Vec4(
+            std::clamp(r, 0.0f, 1.0f),
+            std::clamp(g, 0.0f, 1.0f),
+            std::clamp(b, 0.0f, 1.0f),
+            std::clamp(a, 0.0f, 1.0f)
+        );
+    }
+    
+    // Convert to Vec3 (ignoring alpha)
+    Vec3 toVec3() const {
+        return Vec3(r, g, b);
+    }
+    
+    // Convert to 8-bit color values
+    void toUint8(uint8_t& red, uint8_t& green, uint8_t& blue, uint8_t& alpha) const {
+        red = static_cast<uint8_t>(std::clamp(r, 0.0f, 1.0f) * 255);
+        green = static_cast<uint8_t>(std::clamp(g, 0.0f, 1.0f) * 255);
+        blue = static_cast<uint8_t>(std::clamp(b, 0.0f, 1.0f) * 255);
+        alpha = static_cast<uint8_t>(std::clamp(a, 0.0f, 1.0f) * 255);
+    }
+    
+    void toUint8(uint8_t& red, uint8_t& green, uint8_t& blue) const {
+        red = static_cast<uint8_t>(std::clamp(r, 0.0f, 1.0f) * 255);
+        green = static_cast<uint8_t>(std::clamp(g, 0.0f, 1.0f) * 255);
+        blue = static_cast<uint8_t>(std::clamp(b, 0.0f, 1.0f) * 255);
+    }
+    // Get XYZ components as Vec3
+    class Vec3 xyz() const;
+    
+    // Get RGB components as Vec3
+    class Vec3 rgb() const;
+    
+    bool operator==(const Vec4& other) const {
+        return x == other.x && y == other.y && z == other.z && w == other.w;
+    }
+    
+    bool operator!=(const Vec4& other) const {
+        return x != other.x || y != other.y || z != other.z || w != other.w;
+    }
+    
+    bool operator<(const Vec4& other) const {
+        return (x < other.x) || 
+               (x == other.x && y < other.y) || 
+               (x == other.x && y == other.y && z < other.z) ||
+               (x == other.x && y == other.y && z == other.z && w < other.w);
+    }
+    
+    bool operator<=(const Vec4& other) const {
+        return (x < other.x) || 
+               (x == other.x && y < other.y) || 
+               (x == other.x && y == other.y && z < other.z) ||
+               (x == other.x && y == other.y && z == other.z && w <= other.w);
+    }
+    
+    bool operator>(const Vec4& other) const {
+        return (x > other.x) || 
+               (x == other.x && y > other.y) || 
+               (x == other.x && y == other.y && z > other.z) ||
+               (x == other.x && y == other.y && z == other.z && w > other.w);
+    }
+    
+    bool operator>=(const Vec4& other) const {
+        return (x > other.x) || 
+               (x == other.x && y > other.y) || 
+               (x == other.x && y == other.y && z > other.z) ||
+               (x == other.x && y == other.y && z == other.z && w >= other.w);
+    }
+    
+    Vec4 abs() const {
+        return Vec4(std::abs(x), std::abs(y), std::abs(z), std::abs(w));
+    }
+    
+    Vec4 floor() const {
+        return Vec4(std::floor(x), std::floor(y), std::floor(z), std::floor(w));
+    }
+    
+    Vec4 ceil() const {
+        return Vec4(std::ceil(x), std::ceil(y), std::ceil(z), std::ceil(w));
+    }
+    
+    Vec4 round() const {
+        return Vec4(std::round(x), std::round(y), std::round(z), std::round(w));
+    }
+    
+    Vec4 min(const Vec4& other) const {
+        return Vec4(std::min(x, other.x), std::min(y, other.y), 
+                   std::min(z, other.z), std::min(w, other.w));
+    }
+    
+    Vec4 max(const Vec4& other) const {
+        return Vec4(std::max(x, other.x), std::max(y, other.y), 
+                   std::max(z, other.z), std::max(w, other.w));
+    }
+    
+    Vec4 clamp(float minVal, float maxVal) const {
+        return Vec4(
+            std::clamp(x, minVal, maxVal),
+            std::clamp(y, minVal, maxVal),
+            std::clamp(z, minVal, maxVal),
+            std::clamp(w, minVal, maxVal)
+        );
+    }
+    
+    // Color-specific clamping (clamps RGB between 0 and 1)
+    Vec4 clampColor() const {
+        return Vec4(
+            std::clamp(r, 0.0f, 1.0f),
+            std::clamp(g, 0.0f, 1.0f),
+            std::clamp(b, 0.0f, 1.0f),
+            std::clamp(a, 0.0f, 1.0f)
+        );
+    }
+    
+    bool isZero(float epsilon = 1e-10f) const {
+        return std::abs(x) < epsilon && std::abs(y) < epsilon && 
+               std::abs(z) < epsilon && std::abs(w) < epsilon;
+    }
+    
+    bool equals(const Vec4& other, float epsilon = 1e-10f) const {
+        return std::abs(x - other.x) < epsilon && 
+               std::abs(y - other.y) < epsilon &&
+               std::abs(z - other.z) < epsilon &&
+               std::abs(w - other.w) < epsilon;
+    }
+    
+    friend Vec4 operator+(float scalar, const Vec4& vec) {
+        return Vec4(scalar + vec.x, scalar + vec.y, scalar + vec.z, scalar + vec.w);
+    }
+    
+    friend Vec4 operator-(float scalar, const Vec4& vec) {
+        return Vec4(scalar - vec.x, scalar - vec.y, scalar - vec.z, scalar - vec.w);
+    }
+    
+    friend Vec4 operator*(float scalar, const Vec4& vec) {
+        return Vec4(scalar * vec.x, scalar * vec.y, scalar * vec.z, scalar * vec.w);
+    }
+    
+    friend Vec4 operator/(float scalar, const Vec4& vec) {
+        return Vec4(scalar / vec.x, scalar / vec.y, scalar / vec.z, scalar / vec.w);
+    }
+
+    Vec4 lerp(const Vec4& other, float t) const {
+        t = std::clamp(t, 0.0f, 1.0f);
+        return *this + (other - *this) * t;
+    }
+    
+    // Convert to grayscale using standard RGB weights
+    float grayscale() const {
+        return r * 0.299f + g * 0.587f + b * 0.114f;
+    }
+    
+    // Color inversion (1.0 - color)
+    Vec4 inverted() const {
+        return Vec4(1.0f - r, 1.0f - g, 1.0f - b, a);
+    }
+    
+    float& operator[](int index) {
+        return (&x)[index];
+    }
+    
+    const float& operator[](int index) const {
+        return (&x)[index];
+    }
+    
+    std::string toString() const {
+        return "(" + std::to_string(x) + ", " + std::to_string(y) + ", " + 
+               std::to_string(z) + ", " + std::to_string(w) + ")";
+    }
+    
+    std::string toColorString() const {
+        return "RGBA(" + std::to_string(r) + ", " + std::to_string(g) + ", " + 
+               std::to_string(b) + ", " + std::to_string(a) + ")";
+    }
+};
+
+inline std::ostream& operator<<(std::ostream& os, const Vec4& vec) {
+    os << vec.toString();
+    return os;
+}
+
+namespace std {
+    template<>
+    struct hash<Vec4> {
+        size_t operator()(const Vec4& v) const {
+            return hash<float>()(v.x) ^ (hash<float>()(v.y) << 1) ^ 
+                   (hash<float>()(v.z) << 2) ^ (hash<float>()(v.w) << 3);
+        }
+    };
+}
+
+#endif
\ No newline at end of file
diff --git a/util/voxelgrid.hpp b/util/voxelgrid.hpp
new file mode 100644
index 0000000..fa96836
--- /dev/null
+++ b/util/voxelgrid.hpp
@@ -0,0 +1,217 @@
+#ifndef VOXEL_HPP
+#define VOXEL_HPP
+
+#include <vector>
+#include <cmath>
+#include <limits>
+#include <algorithm>
+#include <unordered_map>
+#include <cstdio>
+#include <cstring>
+#include <fstream>
+#include <iostream>
+#include <vector>
+#include <random>
+#include <functional>
+#include <tuple>
+#include "timing_decorator.hpp"
+#include "vec3.hpp"
+#include "vec4.hpp"
+
+class VoxelGrid {
+private:
+    std::unordered_map<Vec3, size_t> positionToIndex;
+    std::vector<Vec3> positions;
+    std::vector<Vec4> colors;
+    std::vector<int> layers;
+    
+    Vec3 gridSize;
+
+public:
+    Vec3 voxelSize;
+    
+    enum LayerType {
+        ATMOSPHERE = 0,
+        CRUST = 1,
+        MANTLE = 2,
+        OUTER_CORE = 3,
+        INNER_CORE = 4,
+        EMPTY = -1
+    };
+    
+    VoxelGrid(const Vec3& size, const Vec3& voxelSize = Vec3(1, 1, 1)) : gridSize(size), voxelSize(voxelSize) {}
+    
+    void addVoxel(const Vec3& position, const Vec4& color) {
+        Vec3 gridPos = worldToGrid(position);
+        
+        auto it = positionToIndex.find(gridPos);
+        if (it == positionToIndex.end()) {
+            size_t index = positions.size();
+            positions.push_back(gridPos);
+            colors.push_back(color);
+            layers.push_back(EMPTY);
+            positionToIndex[gridPos] = index;
+        } else {
+            colors[it->second] = color;
+        }
+    }
+    
+    void addVoxelWithLayer(const Vec3& position, const Vec4& color, int layer) {
+        Vec3 gridPos = worldToGrid(position);
+        
+        auto it = positionToIndex.find(gridPos);
+        if (it == positionToIndex.end()) {
+            size_t index = positions.size();
+            positions.push_back(gridPos);
+            colors.push_back(color);
+            layers.push_back(layer);
+            positionToIndex[gridPos] = index;
+        } else {
+            colors[it->second] = color;
+            layers[it->second] = layer;
+        }
+    }
+    
+    Vec4 getVoxel(const Vec3& position) const {
+        Vec3 gridPos = worldToGrid(position);
+        auto it = positionToIndex.find(gridPos);
+        if (it != positionToIndex.end()) {
+            return colors[it->second];
+        }
+        return Vec4(0, 0, 0, 0);
+    }
+    
+    int getVoxelLayer(const Vec3& position) const {
+        Vec3 gridPos = worldToGrid(position);
+        auto it = positionToIndex.find(gridPos);
+        if (it != positionToIndex.end()) {
+            return layers[it->second];
+        }
+        return EMPTY;
+    }
+    
+    bool isOccupied(const Vec3& position) const {
+        Vec3 gridPos = worldToGrid(position);
+        return positionToIndex.find(gridPos) != positionToIndex.end();
+    }
+    
+    Vec3 worldToGrid(const Vec3& worldPos) const {
+        return (worldPos / voxelSize).floor();
+    }
+    
+    Vec3 gridToWorld(const Vec3& gridPos) const {
+        return gridPos * voxelSize;
+    }
+    
+    const std::vector<Vec3>& getOccupiedPositions() const {
+        return positions;
+    }
+    
+    const std::vector<Vec4>& getColors() const {
+        return colors;
+    }
+    
+    const std::vector<int>& getLayers() const {
+        return layers;
+    }
+    
+    const std::unordered_map<Vec3, size_t>& getPositionToIndexMap() const {
+        return positionToIndex;
+    }
+    
+    const Vec3& getGridSize() const {
+        return gridSize;
+    }
+    
+    const Vec3& getVoxelSize() const {
+        return voxelSize;
+    }
+    
+    void clear() {
+        positions.clear();
+        colors.clear();
+        layers.clear();
+        positionToIndex.clear();
+    }
+    
+    void assignPlanetaryLayers(const Vec3& center = Vec3(0, 0, 0)) {
+        TIME_FUNCTION;
+        printf("Assigning planetary layers...\n");
+        
+        const float atmospherePercent = 0.05f;
+        const float crustPercent = 0.01f;
+        const float mantlePercent = 0.10f;
+        const float outerCorePercent = 0.42f;
+        const float innerCorePercent = 0.42f;
+        
+        float maxDistance = 0.0f;
+        for (const auto& pos : positions) {
+            Vec3 worldPos = gridToWorld(pos);
+            float distance = (worldPos - center).length();
+            maxDistance = std::max(maxDistance, distance);
+        }
+        
+        printf("Maximum distance from center: %.2f\n", maxDistance);
+        
+        const float atmosphereStart = maxDistance * (1.0f - atmospherePercent);
+        const float crustStart = maxDistance * (1.0f - atmospherePercent - crustPercent);
+        const float mantleStart = maxDistance * (1.0f - atmospherePercent - crustPercent - mantlePercent);
+        const float outerCoreStart = maxDistance * (1.0f - atmospherePercent - crustPercent - mantlePercent - outerCorePercent);
+        
+        printf("Layer boundaries:\n");
+        printf("  Atmosphere: %.2f to %.2f\n", atmosphereStart, maxDistance);
+        printf("  Crust: %.2f to %.2f\n", crustStart, atmosphereStart);
+        printf("  Mantle: %.2f to %.2f\n", mantleStart, crustStart);
+        printf("  Outer Core: %.2f to %.2f\n", outerCoreStart, mantleStart);
+        printf("  Inner Core: 0.00 to %.2f\n", outerCoreStart);
+        
+        int atmosphereCount = 0, crustCount = 0, mantleCount = 0, outerCoreCount = 0, innerCoreCount = 0;
+        
+        for (size_t i = 0; i < positions.size(); ++i) {
+            Vec3 worldPos = gridToWorld(positions[i]);
+            float distance = (worldPos - center).length();
+            
+            Vec4 layerColor;
+            int layerType;
+            
+            if (distance >= atmosphereStart) {
+                // Atmosphere - transparent blue
+                layerColor = Vec4(0.2f, 0.4f, 1.0f, 0.3f); // Semi-transparent blue
+                layerType = ATMOSPHERE;
+                atmosphereCount++;
+            } else if (distance >= crustStart) {
+                // Crust - light brown
+                layerColor = Vec4(0.8f, 0.7f, 0.5f, 1.0f); // Light brown
+                layerType = CRUST;
+                crustCount++;
+            } else if (distance >= mantleStart) {
+                // Mantle - reddish brown
+                layerColor = Vec4(0.7f, 0.3f, 0.2f, 1.0f); // Reddish brown
+                layerType = MANTLE;
+                mantleCount++;
+            } else if (distance >= outerCoreStart) {
+                // Outer Core - orange/yellow
+                layerColor = Vec4(1.0f, 0.6f, 0.2f, 1.0f); // Orange
+                layerType = OUTER_CORE;
+                outerCoreCount++;
+            } else {
+                // Inner Core - bright yellow
+                layerColor = Vec4(1.0f, 0.9f, 0.1f, 1.0f); // Bright yellow
+                layerType = INNER_CORE;
+                innerCoreCount++;
+            }
+            
+            colors[i] = layerColor;
+            layers[i] = layerType;
+        }
+        
+        printf("Layer distribution:\n");
+        printf("  Atmosphere: %d voxels (%.1f%%)\n", atmosphereCount, (atmosphereCount * 100.0f) / positions.size());
+        printf("  Crust: %d voxels (%.1f%%)\n", crustCount, (crustCount * 100.0f) / positions.size());
+        printf("  Mantle: %d voxels (%.1f%%)\n", mantleCount, (mantleCount * 100.0f) / positions.size());
+        printf("  Outer Core: %d voxels (%.1f%%)\n", outerCoreCount, (outerCoreCount * 100.0f) / positions.size());
+        printf("  Inner Core: %d voxels (%.1f%%)\n", innerCoreCount, (innerCoreCount * 100.0f) / positions.size());
+    }
+};
+
+#endif
\ No newline at end of file