new grid2 sim

simtools/grid2.cpp

@@ -0,0 +1,87 @@
+#include "sim2.hpp"
+#include "../util/bmpwriter.hpp"
+#include <iostream>
+#include <string>
+#include <iomanip>
+#include <sstream>
+
+// Function to convert simulation grid to pixel data for BMP
+std::vector<std::vector<Vec3>> gridToPixels(const Sim2& sim, int width, int height) {
+    std::vector<std::vector<Vec3>> pixels(height, std::vector<Vec3>(width, Vec3(0, 0, 0))); // Black background
+    
+    // Get all pixel data from the simulation grid
+    const auto& positions = sim.getPositions();
+    const auto& colors = sim.getColors();
+    
+    for (size_t i = 0; i < positions.size(); ++i) {
+        const Vec2& pos = positions[i];
+        const Vec4& color = colors[i];
+        
+        int x = static_cast<int>(pos.x);
+        int y = static_cast<int>(pos.y);
+        
+        // Only draw if within bounds
+        if (x >= 0 && x < width && y >= 0 && y < height) {
+            // Convert Vec4 (RGBA) to Vec3 (RGB)
+            pixels[y][x] = Vec3(color.x, color.y, color.z);
+        }
+    }
+    
+    return pixels;
+}
+
+int main() {
+    const int FRAME_COUNT = 60;
+    const float TOTAL_TIME = 1.0f; // 1 second
+    const float TIME_STEP = TOTAL_TIME / FRAME_COUNT;
+    
+    // Create output directory
+    std::filesystem::create_directories("output");
+    
+    // Create a simulation with 50x30 grid and gravity
+    Sim2 sim(50, 30, Vec2(0, -9.8f));
+    
+    // Create some objects
+    auto ball1 = sim.createBall(Vec2(10, 25), 3.0f, Vec4(1.0f, 0.5f, 0.0f, 1.0f), 1.0f); // Orange ball
+    auto ball2 = sim.createBall(Vec2(30, 20), 2.0f, Vec4(0.5f, 0.8f, 1.0f, 1.0f), 0.5f); // Blue ball
+    
+    // Create ground
+    auto ground = sim.createGround(Vec2(0, 0), 50, Vec4(0.0f, 1.0f, 0.0f, 1.0f)); // Green ground
+    
+    // Create walls
+    auto leftWall = sim.createWall(Vec2(0, 1), 29, Vec4(0.3f, 0.3f, 0.7f, 1.0f)); // Blue walls
+    auto rightWall = sim.createWall(Vec2(49, 1), 29, Vec4(0.3f, 0.3f, 0.7f, 1.0f));
+    
+    // Set world bounds
+    sim.setWorldBounds(Vec2(0, 0), Vec2(49, 29));
+    
+    // Simulation parameters
+    sim.setElasticity(0.8f); // Bouncy
+    sim.setFriction(0.05f);  // Low friction
+    
+    std::cout << "Rendering " << FRAME_COUNT << " frames over " << TOTAL_TIME << " seconds..." << std::endl;
+    
+    // Run simulation and save frames
+    for (int frame = 0; frame < FRAME_COUNT; ++frame) {
+        // Update simulation
+        sim.update(TIME_STEP);
+        
+        // Convert simulation state to pixel data
+        auto pixels = gridToPixels(sim, 50, 30);
+        
+        // Create filename with zero-padded frame number
+        std::ostringstream filename;
+        filename << "output/bounce" << std::setw(3) << std::setfill('0') << frame << ".bmp";
+        
+        // Save as BMP
+        if (BMPWriter::saveBMP(filename.str(), pixels)) {
+            std::cout << "Saved frame " << frame << " to " << filename.str() << std::endl;
+        } else {
+            std::cerr << "Failed to save frame " << frame << std::endl;
+        }
+    }
+    
+    std::cout << "Animation complete! " << FRAME_COUNT << " frames saved to output/ directory." << std::endl;
+    
+    return 0;
+}

simtools/sim2.hpp

@@ -1,265 +1,336 @@
 #ifndef SIM2_HPP
 #define SIM2_HPP
 
-#include "../util/noise2.hpp"
+#include "../util/grid/grid2.hpp"
-#include "../util/grid2.hpp"
 #include "../util/vec2.hpp"
-#include "../util/vec4.hpp"
+#include <vector>
-#include "../util/timing_decorator.hpp"
 #include <memory>
-#include <string>
 #include <unordered_map>
 
-class Sim2 {
+class PhysicsObject {
-private:
+protected:
-    std::unique_ptr<Noise2> noiseGenerator;
+    Vec2 position;
-    Grid2 terrainGrid;
+    Vec2 velocity;
-    int gridWidth;
+    Vec2 acceleration;
-    int gridHeight;
+    float mass;
-    
+    bool isStatic;
-    // Terrain generation parameters
+    std::shared_ptr<Grid2> shape;
-    float scale;
-    int octaves;
-    float persistence;
-    float lacunarity;
-    uint32_t seed;
-    Vec2 offset;
-    
-    // Terrain modification parameters
-    float elevationMultiplier;
-    float waterLevel;
-    Vec4 landColor;
-    Vec4 waterColor;
     
 public:
-    Sim2(int width = 512, int height = 512, uint32_t seed = 42, float scale = 4.0f, int octaves = 4,
+    PhysicsObject(const Vec2& pos, float mass = 1.0f, bool isStatic = false)
-        float persistence = 0.5f, float lacunarity = 2.0f, float waterlevel = 0.3f, float elevation = 1.0f) 
+        : position(pos), velocity(0, 0), acceleration(0, 0), mass(mass), isStatic(isStatic) {}
-        : gridWidth(width), gridHeight(height), scale(scale), octaves(octaves), 
+    
-          persistence(persistence), lacunarity(lacunarity), seed(seed), offset(0, 0),
+    virtual ~PhysicsObject() = default;
-          elevationMultiplier(elevation), waterLevel(waterlevel),
+    
-          landColor(0.2f, 0.8f, 0.2f, 1.0f),  // Green
+    virtual void update(float dt) {
-          waterColor(0.2f, 0.3f, 0.8f, 1.0f)  // Blue
+        if (isStatic) return;
-    {
+        
-        noiseGenerator = std::make_unique<Noise2>(seed,Noise2::PERLIN,Noise2::PRECOMPUTED);
+        // Update velocity and position using basic physics
-        generateTerrain();
+        velocity += acceleration * dt;
+        position += velocity * dt;
+        acceleration = Vec2(0, 0); // Reset acceleration for next frame
     }
     
-    // Generate initial terrain
+    virtual void applyForce(const Vec2& force) {
-    void generateTerrain() {
+        if (!isStatic) {
-        TIME_FUNCTION;
+            acceleration += force / mass;
-        terrainGrid = noiseGenerator->generateTerrainNoise(
-            gridWidth, gridHeight, scale, octaves, persistence, seed, offset);
-        
-        applyTerrainColors();
-    }
-    
-    // Regenerate terrain with current parameters
-    void regenerate() {
-        generateTerrain();
-    }
-    
-    // Basic parameter modifications
-    void setScale(float newScale) {
-        scale = std::max(0.1f, newScale);
-        generateTerrain();
-    }
-    
-    void setOctaves(int newOctaves) {
-        octaves = std::max(1, newOctaves);
-        generateTerrain();
-    }
-    
-    void setPersistence(float newPersistence) {
-        persistence = std::clamp(newPersistence, 0.0f, 1.0f);
-        generateTerrain();
-    }
-    
-    void setLacunarity(float newLacunarity) {
-        lacunarity = std::max(1.0f, newLacunarity);
-        generateTerrain();
-    }
-    
-    void setSeed(uint32_t newSeed) {
-        seed = newSeed;
-        noiseGenerator->setSeed(seed);
-        generateTerrain();
-    }
-    
-    void setOffset(const Vec2& newOffset) {
-        offset = newOffset;
-        generateTerrain();
-    }
-    
-    void setElevationMultiplier(float multiplier) {
-        elevationMultiplier = std::max(0.0f, multiplier);
-        applyElevationModification();
-    }
-    
-    void setWaterLevel(float level) {
-        waterLevel = std::clamp(level, 0.0f, 1.0f);
-        applyTerrainColors();
-    }
-    
-    void setLandColor(const Vec4& color) {
-        landColor = color;
-        applyTerrainColors();
-    }
-    
-    void setWaterColor(const Vec4& color) {
-        waterColor = color;
-        applyTerrainColors();
-    }
-    
-    // Get current parameters
-    float getScale() const { return scale; }
-    int getOctaves() const { return octaves; }
-    float getPersistence() const { return persistence; }
-    float getLacunarity() const { return lacunarity; }
-    uint32_t getSeed() const { return seed; }
-    Vec2 getOffset() const { return offset; }
-    float getElevationMultiplier() const { return elevationMultiplier; }
-    float getWaterLevel() const { return waterLevel; }
-    Vec4 getLandColor() const { return landColor; }
-    Vec4 getWaterColor() const { return waterColor; }
-    
-    // Get the terrain grid
-    const Grid2& getTerrainGrid() const {
-        return terrainGrid;
-    }
-    
-    // Get terrain dimensions
-    int getWidth() const { return gridWidth; }
-    int getHeight() const { return gridHeight; }
-    
-    // Get elevation at specific coordinates
-    float getElevation(int x, int y) const {
-        if (x < 0 || x >= gridWidth || y < 0 || y >= gridHeight) {
-            return 0.0f;
-        }
-        return terrainGrid.colors[y * gridWidth + x].x; // Elevation stored in red channel
-    }
-    
-    // Render to RGB image
-    std::vector<uint8_t> renderToRGB(int width, int height, 
-                                    const Vec4& backgroundColor = Vec4(0, 0, 0, 1)) const {
-        return terrainGrid.renderToRGB(width, height, backgroundColor);
-    }
-    
-    // Render to RGBA image
-    std::vector<uint8_t> renderToRGBA(int width, int height, 
-                                     const Vec4& backgroundColor = Vec4(0, 0, 0, 1)) const {
-        return terrainGrid.renderToRGBA(width, height, backgroundColor);
-    }
-    
-    // Export terrain data as heightmap (grayscale)
-    Grid2 exportHeightmap() const {
-        Grid2 heightmap(gridWidth * gridHeight);
-        
-        for (int y = 0; y < gridHeight; y++) {
-            for (int x = 0; x < gridWidth; x++) {
-                int index = y * gridWidth + x;
-                float elevation = terrainGrid.colors[index].x;
-                heightmap.positions[index] = Vec2(x, y);
-                heightmap.colors[index] = Vec4(elevation, elevation, elevation, 1.0f);
         }
     }
     
-        return heightmap;
+    // Getters
+    const Vec2& getPosition() const { return position; }
+    const Vec2& getVelocity() const { return velocity; }
+    float getMass() const { return mass; }
+    bool getIsStatic() const { return isStatic; }
+    std::shared_ptr<Grid2> getShape() const { return shape; }
+    
+    // Setters
+    void setPosition(const Vec2& pos) { position = pos; }
+    void setVelocity(const Vec2& vel) { velocity = vel; }
+    void setMass(float m) { mass = m; }
+    void setShape(std::shared_ptr<Grid2> newShape) { shape = newShape; }
+    
+    // Check if this object contains a world position
+    virtual bool contains(const Vec2& worldPos) const {
+        if (!shape) return false;
+        
+        Vec2 localPos = worldPos - position;
+        return shape->isOccupied(localPos.round());
     }
     
-    // Generate random seed and regenerate
+    // Get all occupied positions in world coordinates
-    void randomizeSeed() {
+    virtual std::vector<Vec2> getWorldPositions() const {
-        std::random_device rd;
+        std::vector<Vec2> worldPositions;
-        setSeed(rd());
+        if (!shape) return worldPositions;
-    }
         
-    // Reset all parameters to default
+        const auto& localPositions = shape->getPositions();
-    void reset() {
+        for (const auto& localPos : localPositions) {
-        scale = 4.0f;
+            worldPositions.push_back(position + localPos);
-        octaves = 4;
+        }
-        persistence = 0.5f;
+        return worldPositions;
-        lacunarity = 2.0f;
-        elevationMultiplier = 1.0f;
-        waterLevel = 0.3f;
-        landColor = Vec4(0.2f, 0.8f, 0.2f, 1.0f);
-        waterColor = Vec4(0.2f, 0.3f, 0.8f, 1.0f);
-        generateTerrain();
     }
-    
-    // Get terrain statistics
-    struct TerrainStats {
-        float minElevation;
-        float maxElevation;
-        float averageElevation;
-        float landPercentage;
-        int landArea;
-        int waterArea;
 };
 
-    TerrainStats getTerrainStats() const {
+class FallingObject : public PhysicsObject {
-        TerrainStats stats = {1.0f, 0.0f, 0.0f, 0.0f, 0, 0};
+private:
-        float totalElevation = 0.0f;
+    Vec2 gravityForce;
     
-        for (const auto& color : terrainGrid.colors) {
+public:
-            float elevation = color.x;
+    FallingObject(const Vec2& pos, float mass = 1.0f, const Vec2& gravity = Vec2(0, -9.8f))
-            stats.minElevation = std::min(stats.minElevation, elevation);
+        : PhysicsObject(pos, mass, false), gravityForce(gravity * mass) {}
-            stats.maxElevation = std::max(stats.maxElevation, elevation);
-            totalElevation += elevation;
     
-            if (elevation > waterLevel) {
+    void update(float dt) override {
-                stats.landArea++;
+        if (!isStatic) {
-            } else {
+            // Apply gravity
-                stats.waterArea++;
+            applyForce(gravityForce);
         }
+        PhysicsObject::update(dt);
     }
     
-        stats.averageElevation = totalElevation / terrainGrid.colors.size();
+    void setGravity(const Vec2& gravity) {
-        stats.landPercentage = static_cast<float>(stats.landArea) / 
+        gravityForce = gravity * mass;
-                              (stats.landArea + stats.waterArea) * 100.0f;
+    }
     
-        return stats;
+    const Vec2 getGravity() const {
+        return gravityForce / mass;
+    }
+};
+
+class SolidObject : public PhysicsObject {
+public:
+    SolidObject(const Vec2& pos, float mass = 1.0f)
+        : PhysicsObject(pos, mass, true) {} // Solids are static by default
+    
+    // Solids don't move, so override update to do nothing
+    void update(float dt) override {
+        // Solids don't move
+    }
+    
+    void applyForce(const Vec2& force) override {
+        // Solids don't respond to forces
+    }
+};
+
+class Sim2 : public Grid2 {
+private:
+    std::vector<std::shared_ptr<PhysicsObject>> objects;
+    Vec2 globalGravity;
+    float elasticity; // Bounciness factor (0.0 to 1.0)
+    float friction;   // Surface friction (0.0 to 1.0)
+    Vec2 worldBoundsMin, worldBoundsMax;
+    bool useWorldBounds;
+    
+public:
+    Sim2(int width, int height, const Vec2& gravity = Vec2(0, -9.8f))
+        : Grid2(width, height), globalGravity(gravity), elasticity(0.7f), 
+          friction(0.1f), useWorldBounds(false) {}
+    
+    // Add objects to simulation
+    void addObject(std::shared_ptr<PhysicsObject> object) {
+        objects.push_back(object);
+    }
+    
+    // Create a falling ball
+    std::shared_ptr<FallingObject> createBall(const Vec2& position, float radius, 
+                                            const Vec4& color = Vec4(1.0f, 0.5f, 0.0f, 1.0f), 
+                                            float mass = 1.0f) {
+        auto ball = std::make_shared<FallingObject>(position, mass, globalGravity);
+        auto shape = std::make_shared<Grid2>(static_cast<int>(radius * 2 + 2));
+        shape->fillCircle(Vec2(radius, radius), radius, color);
+        ball->setShape(shape);
+        objects.push_back(ball);
+        return ball;
+    }
+    
+    // Create a solid ground
+    std::shared_ptr<SolidObject> createGround(const Vec2& position, int width, 
+                                            const Vec4& color = Vec4(0.0f, 1.0f, 0.0f, 1.0f)) {
+        auto ground = std::make_shared<SolidObject>(position);
+        auto shape = std::make_shared<Grid2>(width, 1);
+        for (int x = 0; x < width; ++x) {
+            shape->addPixel(x, 0, color);
+        }
+        ground->setShape(shape);
+        objects.push_back(ground);
+        return ground;
+    }
+    
+    // Create a solid wall
+    std::shared_ptr<SolidObject> createWall(const Vec2& position, int height,
+                                          const Vec4& color = Vec4(0.3f, 0.3f, 0.7f, 1.0f)) {
+        auto wall = std::make_shared<SolidObject>(position);
+        auto shape = std::make_shared<Grid2>(1, height);
+        for (int y = 0; y < height; ++y) {
+            shape->addPixel(0, y, color);
+        }
+        wall->setShape(shape);
+        objects.push_back(wall);
+        return wall;
+    }
+    
+    // Update the simulation
+    void update(float dt) {
+        // Clear the main grid
+        clear();
+        
+        // Update all objects
+        for (auto& obj : objects) {
+            obj->update(dt);
+        }
+        
+        // Handle collisions
+        handleCollisions();
+        
+        // Handle world bounds
+        if (useWorldBounds) {
+            handleWorldBounds();
+        }
+        
+        // Render all objects to the main grid
+        renderObjects();
+    }
+    
+    // Set world bounds for containment
+    void setWorldBounds(const Vec2& min, const Vec2& max) {
+        worldBoundsMin = min;
+        worldBoundsMax = max;
+        useWorldBounds = true;
+    }
+    
+    // Simulation parameters
+    void setElasticity(float e) { elasticity = std::clamp(e, 0.0f, 1.0f); }
+    void setFriction(float f) { friction = std::clamp(f, 0.0f, 1.0f); }
+    void setGlobalGravity(const Vec2& gravity) { globalGravity = gravity; }
+    
+    float getElasticity() const { return elasticity; }
+    float getFriction() const { return friction; }
+    const Vec2& getGlobalGravity() const { return globalGravity; }
+    
+    // Get all objects
+    const std::vector<std::shared_ptr<PhysicsObject>>& getObjects() const {
+        return objects;
+    }
+    
+    // Clear all objects
+    void clearObjects() {
+        objects.clear();
+        clear();
     }
 
 private:
-    void applyTerrainColors() {
+    void handleCollisions() {
-        for (int y = 0; y < gridHeight; y++) {
+        // Simple collision detection between falling objects and solids
-            for (int x = 0; x < gridWidth; x++) {
+        for (size_t i = 0; i < objects.size(); ++i) {
-                int index = y * gridWidth + x;
+            auto obj1 = objects[i];
-                float elevation = terrainGrid.colors[index].x;
+            if (obj1->getIsStatic()) continue; // Skip static objects for collision detection
             
-                // Apply water level and color based on elevation
+            for (size_t j = 0; j < objects.size(); ++j) {
-                if (elevation <= waterLevel) {
+                if (i == j) continue;
-                    // Water - optionally add depth variation
+                
-                    float depth = (waterLevel - elevation) / waterLevel;
+                auto obj2 = objects[j];
-                    Vec4 water = waterColor * (0.7f + 0.3f * depth);
+                if (!obj2->getIsStatic()) continue; // Only check against static objects for now
-                    water.w = 1.0f; // Ensure full alpha
+                
-                    terrainGrid.colors[index] = water;
+                checkAndResolveCollision(obj1, obj2);
-                } else {
-                    // Land - optionally add elevation-based color variation
-                    float height = (elevation - waterLevel) / (1.0f - waterLevel);
-                    Vec4 land = landColor * (0.8f + 0.2f * height);
-                    land.w = 1.0f; // Ensure full alpha
-                    terrainGrid.colors[index] = land;
-                }
             }
         }
     }
     
-    void applyElevationModification() {
+    void checkAndResolveCollision(std::shared_ptr<PhysicsObject> dynamicObj, 
-        for (int y = 0; y < gridHeight; y++) {
+                                 std::shared_ptr<PhysicsObject> staticObj) {
-            for (int x = 0; x < gridWidth; x++) {
+        if (!dynamicObj->getShape() || !staticObj->getShape()) return;
-                int index = y * gridWidth + x;
-                float originalElevation = terrainGrid.colors[index].x;
         
-                // Apply elevation multiplier with clamping
+        // Get all world positions of the dynamic object
-                float newElevation = std::clamp(originalElevation * elevationMultiplier, 0.0f, 1.0f);
+        auto dynamicPositions = dynamicObj->getWorldPositions();
-                terrainGrid.colors[index].x = newElevation;
+        
-                terrainGrid.colors[index].y = newElevation; // Keep grayscale for heightmap
+        for (const auto& dynPos : dynamicPositions) {
-                terrainGrid.colors[index].z = newElevation;
+            // Check if this position collides with the static object
+            if (staticObj->contains(dynPos)) {
+                // Simple collision response - bounce
+                Vec2 velocity = dynamicObj->getVelocity();
+                
+                // Determine collision normal (simplified - always bounce up)
+                Vec2 normal(0, 1);
+                
+                // Reflect velocity with elasticity
+                Vec2 reflectedVel = velocity.reflect(normal) * elasticity;
+                
+                // Apply friction
+                reflectedVel.x *= (1.0f - friction);
+                
+                dynamicObj->setVelocity(reflectedVel);
+                
+                // Move object out of collision
+                Vec2 newPos = dynamicObj->getPosition();
+                newPos.y += 1.0f; // Move up by 1 pixel
+                dynamicObj->setPosition(newPos);
+                
+                break; // Only handle first collision
+            }
         }
     }
     
-        applyTerrainColors();
+    void handleWorldBounds() {
+        for (auto& obj : objects) {
+            if (obj->getIsStatic()) continue;
+            
+            Vec2 pos = obj->getPosition();
+            Vec2 vel = obj->getVelocity();
+            bool collided = false;
+            
+            // Check bounds and bounce
+            if (pos.x < worldBoundsMin.x) {
+                pos.x = worldBoundsMin.x;
+                vel.x = -vel.x * elasticity;
+                collided = true;
+            } else if (pos.x > worldBoundsMax.x) {
+                pos.x = worldBoundsMax.x;
+                vel.x = -vel.x * elasticity;
+                collided = true;
+            }
+            
+            if (pos.y < worldBoundsMin.y) {
+                pos.y = worldBoundsMin.y;
+                vel.y = -vel.y * elasticity;
+                collided = true;
+            } else if (pos.y > worldBoundsMax.y) {
+                pos.y = worldBoundsMax.y;
+                vel.y = -vel.y * elasticity;
+                collided = true;
+            }
+            
+            if (collided) {
+                obj->setPosition(pos);
+                obj->setVelocity(vel);
+            }
+        }
+    }
+    
+    void renderObjects() {
+        // Render all objects to the main grid
+        for (const auto& obj : objects) {
+            if (!obj->getShape()) continue;
+            
+            const auto& shape = obj->getShape();
+            const Vec2& objPos = obj->getPosition();
+            
+            // Copy all pixels from object's shape to main grid at object's position
+            const auto& positions = shape->getPositions();
+            const auto& colors = shape->getColors();
+            const auto& sizes = shape->getSizes();
+            
+            for (size_t i = 0; i < positions.size(); ++i) {
+                Vec2 worldPos = objPos + positions[i];
+                Vec2 gridPos = worldPos.round();
+                
+                // Only add if within grid bounds
+                if (gridPos.x >= 0 && gridPos.x < width && 
+                    gridPos.y >= 0 && gridPos.y < height) {
+                    if (!isOccupied(gridPos)) {
+                        addPixel(gridPos, colors[i], sizes[i]);
+                    }
+                }
+            }
+        }
     }
 };
 

simtools/sim2old.hpp

@@ -0,0 +1,266 @@
+#ifndef SIM2_HPP
+#define SIM2_HPP
+
+#include "../util/noise2.hpp"
+#include "../util/grid2.hpp"
+#include "../util/vec2.hpp"
+#include "../util/vec4.hpp"
+#include "../util/timing_decorator.hpp"
+#include <memory>
+#include <string>
+#include <unordered_map>
+
+class Sim2 {
+private:
+    std::unique_ptr<Noise2> noiseGenerator;
+    Grid2 terrainGrid;
+    int gridWidth;
+    int gridHeight;
+    
+    // Terrain generation parameters
+    float scale;
+    int octaves;
+    float persistence;
+    float lacunarity;
+    uint32_t seed;
+    Vec2 offset;
+    
+    // Terrain modification parameters
+    float elevationMultiplier;
+    float waterLevel;
+    Vec4 landColor;
+    Vec4 waterColor;
+    
+public:
+    Sim2(int width = 512, int height = 512, uint32_t seed = 42, float scale = 4.0f, int octaves = 4,
+        float persistence = 0.5f, float lacunarity = 2.0f, float waterlevel = 0.3f, float elevation = 1.0f) 
+        : gridWidth(width), gridHeight(height), scale(scale), octaves(octaves), 
+          persistence(persistence), lacunarity(lacunarity), seed(seed), offset(0, 0),
+          elevationMultiplier(elevation), waterLevel(waterlevel),
+          landColor(0.2f, 0.8f, 0.2f, 1.0f),  // Green
+          waterColor(0.2f, 0.3f, 0.8f, 1.0f)  // Blue
+    {
+        noiseGenerator = std::make_unique<Noise2>(seed,Noise2::PERLIN,Noise2::PRECOMPUTED);
+        generateTerrain();
+    }
+    
+    // Generate initial terrain
+    void generateTerrain() {
+        TIME_FUNCTION;
+        terrainGrid = noiseGenerator->generateTerrainNoise(
+            gridWidth, gridHeight, scale, octaves, persistence, seed, offset);
+        
+        applyTerrainColors();
+    }
+    
+    // Regenerate terrain with current parameters
+    void regenerate() {
+        generateTerrain();
+    }
+    
+    // Basic parameter modifications
+    void setScale(float newScale) {
+        scale = std::max(0.1f, newScale);
+        generateTerrain();
+    }
+    
+    void setOctaves(int newOctaves) {
+        octaves = std::max(1, newOctaves);
+        generateTerrain();
+    }
+    
+    void setPersistence(float newPersistence) {
+        persistence = std::clamp(newPersistence, 0.0f, 1.0f);
+        generateTerrain();
+    }
+    
+    void setLacunarity(float newLacunarity) {
+        lacunarity = std::max(1.0f, newLacunarity);
+        generateTerrain();
+    }
+    
+    void setSeed(uint32_t newSeed) {
+        seed = newSeed;
+        noiseGenerator->setSeed(seed);
+        generateTerrain();
+    }
+    
+    void setOffset(const Vec2& newOffset) {
+        offset = newOffset;
+        generateTerrain();
+    }
+    
+    void setElevationMultiplier(float multiplier) {
+        elevationMultiplier = std::max(0.0f, multiplier);
+        applyElevationModification();
+    }
+    
+    void setWaterLevel(float level) {
+        waterLevel = std::clamp(level, 0.0f, 1.0f);
+        applyTerrainColors();
+    }
+    
+    void setLandColor(const Vec4& color) {
+        landColor = color;
+        applyTerrainColors();
+    }
+    
+    void setWaterColor(const Vec4& color) {
+        waterColor = color;
+        applyTerrainColors();
+    }
+    
+    // Get current parameters
+    float getScale() const { return scale; }
+    int getOctaves() const { return octaves; }
+    float getPersistence() const { return persistence; }
+    float getLacunarity() const { return lacunarity; }
+    uint32_t getSeed() const { return seed; }
+    Vec2 getOffset() const { return offset; }
+    float getElevationMultiplier() const { return elevationMultiplier; }
+    float getWaterLevel() const { return waterLevel; }
+    Vec4 getLandColor() const { return landColor; }
+    Vec4 getWaterColor() const { return waterColor; }
+    
+    // Get the terrain grid
+    const Grid2& getTerrainGrid() const {
+        return terrainGrid;
+    }
+    
+    // Get terrain dimensions
+    int getWidth() const { return gridWidth; }
+    int getHeight() const { return gridHeight; }
+    
+    // Get elevation at specific coordinates
+    float getElevation(int x, int y) const {
+        if (x < 0 || x >= gridWidth || y < 0 || y >= gridHeight) {
+            return 0.0f;
+        }
+        return terrainGrid.colors[y * gridWidth + x].x; // Elevation stored in red channel
+    }
+    
+    // Render to RGB image
+    std::vector<uint8_t> renderToRGB(int width, int height, 
+                                    const Vec4& backgroundColor = Vec4(0, 0, 0, 1)) const {
+        return terrainGrid.renderToRGB(width, height, backgroundColor);
+    }
+    
+    // Render to RGBA image
+    std::vector<uint8_t> renderToRGBA(int width, int height, 
+                                     const Vec4& backgroundColor = Vec4(0, 0, 0, 1)) const {
+        return terrainGrid.renderToRGBA(width, height, backgroundColor);
+    }
+    
+    // Export terrain data as heightmap (grayscale)
+    Grid2 exportHeightmap() const {
+        Grid2 heightmap(gridWidth * gridHeight);
+        
+        for (int y = 0; y < gridHeight; y++) {
+            for (int x = 0; x < gridWidth; x++) {
+                int index = y * gridWidth + x;
+                float elevation = terrainGrid.colors[index].x;
+                heightmap.positions[index] = Vec2(x, y);
+                heightmap.colors[index] = Vec4(elevation, elevation, elevation, 1.0f);
+            }
+        }
+        
+        return heightmap;
+    }
+    
+    // Generate random seed and regenerate
+    void randomizeSeed() {
+        std::random_device rd;
+        setSeed(rd());
+    }
+    
+    // Reset all parameters to default
+    void reset() {
+        scale = 4.0f;
+        octaves = 4;
+        persistence = 0.5f;
+        lacunarity = 2.0f;
+        elevationMultiplier = 1.0f;
+        waterLevel = 0.3f;
+        landColor = Vec4(0.2f, 0.8f, 0.2f, 1.0f);
+        waterColor = Vec4(0.2f, 0.3f, 0.8f, 1.0f);
+        generateTerrain();
+    }
+    
+    // Get terrain statistics
+    struct TerrainStats {
+        float minElevation;
+        float maxElevation;
+        float averageElevation;
+        float landPercentage;
+        int landArea;
+        int waterArea;
+    };
+    
+    TerrainStats getTerrainStats() const {
+        TerrainStats stats = {1.0f, 0.0f, 0.0f, 0.0f, 0, 0};
+        float totalElevation = 0.0f;
+        
+        for (const auto& color : terrainGrid.colors) {
+            float elevation = color.x;
+            stats.minElevation = std::min(stats.minElevation, elevation);
+            stats.maxElevation = std::max(stats.maxElevation, elevation);
+            totalElevation += elevation;
+            
+            if (elevation > waterLevel) {
+                stats.landArea++;
+            } else {
+                stats.waterArea++;
+            }
+        }
+        
+        stats.averageElevation = totalElevation / terrainGrid.colors.size();
+        stats.landPercentage = static_cast<float>(stats.landArea) / 
+                              (stats.landArea + stats.waterArea) * 100.0f;
+        
+        return stats;
+    }
+
+private:
+    void applyTerrainColors() {
+        for (int y = 0; y < gridHeight; y++) {
+            for (int x = 0; x < gridWidth; x++) {
+                int index = y * gridWidth + x;
+                float elevation = terrainGrid.colors[index].x;
+                
+                // Apply water level and color based on elevation
+                if (elevation <= waterLevel) {
+                    // Water - optionally add depth variation
+                    float depth = (waterLevel - elevation) / waterLevel;
+                    Vec4 water = waterColor * (0.7f + 0.3f * depth);
+                    water.w = 1.0f; // Ensure full alpha
+                    terrainGrid.colors[index] = water;
+                } else {
+                    // Land - optionally add elevation-based color variation
+                    float height = (elevation - waterLevel) / (1.0f - waterLevel);
+                    Vec4 land = landColor * (0.8f + 0.2f * height);
+                    land.w = 1.0f; // Ensure full alpha
+                    terrainGrid.colors[index] = land;
+                }
+            }
+        }
+    }
+    
+    void applyElevationModification() {
+        for (int y = 0; y < gridHeight; y++) {
+            for (int x = 0; x < gridWidth; x++) {
+                int index = y * gridWidth + x;
+                float originalElevation = terrainGrid.colors[index].x;
+                
+                // Apply elevation multiplier with clamping
+                float newElevation = std::clamp(originalElevation * elevationMultiplier, 0.0f, 1.0f);
+                terrainGrid.colors[index].x = newElevation;
+                terrainGrid.colors[index].y = newElevation; // Keep grayscale for heightmap
+                terrainGrid.colors[index].z = newElevation;
+            }
+        }
+        
+        applyTerrainColors();
+    }
+};
+
+#endif

util/grid/grid2.hpp

@@ -0,0 +1,295 @@
+#ifndef GRID2_HPP
+#define GRID2_HPP
+
+#include "../vec2.hpp"
+#include "../vec4.hpp"
+#include <vector>
+#include <unordered_map>
+#include <string>
+#include <algorithm>
+
+class Grid2 {
+public:
+    int width, height;
+    
+    Grid2() : width(0), height(0) {}
+    Grid2(int size) : width(size), height(size) {
+        positions.reserve(size);
+        colors.reserve(size);
+        sizes.reserve(size);
+    }
+    Grid2(int width, int height) : width(width), height(height) {
+        positions.reserve(width * height);
+        colors.reserve(width * height);
+        sizes.reserve(width * height);
+    }
+    
+    // Add a pixel at specific position
+    int addPixel(const Vec2& position, const Vec4& color, float size = 1.0f) {
+        int index = positions.size();
+        positions.push_back(position);
+        colors.push_back(color);
+        sizes.push_back(size);
+        positionToIndex[position] = index;
+        return index;
+    }
+    
+    // Add a pixel with integer coordinates
+    int addPixel(int x, int y, const Vec4& color, float size = 1.0f) {
+        return addPixel(Vec2(static_cast<float>(x), static_cast<float>(y)), color, size);
+    }
+    
+    // Check if position is occupied
+    bool isOccupied(const Vec2& position) const {
+        return positionToIndex.find(position) != positionToIndex.end();
+    }
+    
+    bool isOccupied(int x, int y) const {
+        return isOccupied(Vec2(static_cast<float>(x), static_cast<float>(y)));
+    }
+    
+    // Get pixel index at position, returns -1 if not found
+    int getPixelIndex(const Vec2& position) const {
+        auto it = positionToIndex.find(position);
+        return (it != positionToIndex.end()) ? it->second : -1;
+    }
+    
+    int getPixelIndex(int x, int y) const {
+        return getPixelIndex(Vec2(static_cast<float>(x), static_cast<float>(y)));
+    }
+    
+    // Remove pixel at position
+    bool removePixel(const Vec2& position) {
+        int index = getPixelIndex(position);
+        if (index == -1) return false;
+        
+        // Swap with last element and update map
+        if (index != positions.size() - 1) {
+            positions[index] = positions.back();
+            colors[index] = colors.back();
+            sizes[index] = sizes.back();
+            
+            // Update mapping for the moved element
+            positionToIndex[positions[index]] = index;
+        }
+        
+        // Remove last element
+        positions.pop_back();
+        colors.pop_back();
+        sizes.pop_back();
+        positionToIndex.erase(position);
+        
+        return true;
+    }
+    
+    bool removePixel(int x, int y) {
+        return removePixel(Vec2(static_cast<float>(x), static_cast<float>(y)));
+    }
+    
+    // Clear all pixels
+    void clear() {
+        positions.clear();
+        colors.clear();
+        sizes.clear();
+        positionToIndex.clear();
+    }
+    
+    // Get pixel count
+    size_t getPixelCount() const {
+        return positions.size();
+    }
+    
+    // Check if grid is empty
+    bool isEmpty() const {
+        return positions.empty();
+    }
+    
+    // Get bounding box of occupied pixels
+    void getBoundingBox(Vec2& minCorner, Vec2& maxCorner) const {
+        if (positions.empty()) {
+            minCorner = Vec2(0, 0);
+            maxCorner = Vec2(0, 0);
+            return;
+        }
+        
+        minCorner = positions[0];
+        maxCorner = positions[0];
+        
+        for (const auto& pos : positions) {
+            minCorner = minCorner.min(pos);
+            maxCorner = maxCorner.max(pos);
+        }
+    }
+    
+    // Fill a rectangular region
+    void fillRectangle(const Vec2& start, const Vec2& end, const Vec4& color, float size = 1.0f) {
+        int startX = static_cast<int>(std::min(start.x, end.x));
+        int endX = static_cast<int>(std::max(start.x, end.x));
+        int startY = static_cast<int>(std::min(start.y, end.y));
+        int endY = static_cast<int>(std::max(start.y, end.y));
+        
+        for (int y = startY; y <= endY; ++y) {
+            for (int x = startX; x <= endX; ++x) {
+                if (!isOccupied(x, y)) {
+                    addPixel(x, y, color, size);
+                }
+            }
+        }
+    }
+    
+    // Create a circle pattern
+    void fillCircle(const Vec2& center, float radius, const Vec4& color, float size = 1.0f) {
+        int centerX = static_cast<int>(center.x);
+        int centerY = static_cast<int>(center.y);
+        int radiusInt = static_cast<int>(radius);
+        
+        for (int y = centerY - radiusInt; y <= centerY + radiusInt; ++y) {
+            for (int x = centerX - radiusInt; x <= centerX + radiusInt; ++x) {
+                float dx = x - center.x;
+                float dy = y - center.y;
+                if (dx * dx + dy * dy <= radius * radius) {
+                    if (!isOccupied(x, y)) {
+                        addPixel(x, y, color, size);
+                    }
+                }
+            }
+        }
+    }
+    
+    // Create a line between two points
+    void drawLine(const Vec2& start, const Vec2& end, const Vec4& color, float size = 1.0f) {
+        // Bresenham's line algorithm
+        int x0 = static_cast<int>(start.x);
+        int y0 = static_cast<int>(start.y);
+        int x1 = static_cast<int>(end.x);
+        int y1 = static_cast<int>(end.y);
+        
+        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 (!isOccupied(x0, y0)) {
+                addPixel(x0, y0, color, size);
+            }
+            
+            if (x0 == x1 && y0 == y1) break;
+            
+            int e2 = 2 * err;
+            if (e2 > -dy) {
+                err -= dy;
+                x0 += sx;
+            }
+            if (e2 < dx) {
+                err += dx;
+                y0 += sy;
+            }
+        }
+    }
+    
+    // Get neighbors of a pixel (4-connected)
+    std::vector<int> getNeighbors4(const Vec2& position) const {
+        std::vector<int> neighbors;
+        Vec2 offsets[] = {Vec2(1, 0), Vec2(-1, 0), Vec2(0, 1), Vec2(0, -1)};
+        
+        for (const auto& offset : offsets) {
+            Vec2 neighborPos = position + offset;
+            int index = getPixelIndex(neighborPos);
+            if (index != -1) {
+                neighbors.push_back(index);
+            }
+        }
+        
+        return neighbors;
+    }
+    
+    // Get neighbors of a pixel (8-connected)
+    std::vector<int> getNeighbors8(const Vec2& position) const {
+        std::vector<int> neighbors;
+        
+        for (int dy = -1; dy <= 1; ++dy) {
+            for (int dx = -1; dx <= 1; ++dx) {
+                if (dx == 0 && dy == 0) continue;
+                
+                Vec2 neighborPos = position + Vec2(dx, dy);
+                int index = getPixelIndex(neighborPos);
+                if (index != -1) {
+                    neighbors.push_back(index);
+                }
+            }
+        }
+        
+        return neighbors;
+    }
+    
+    // Find connected components
+    std::vector<std::vector<int>> findConnectedComponents() const {
+        std::vector<std::vector<int>> components;
+        std::unordered_map<Vec2, bool, std::hash<Vec2>> visited;
+        
+        for (size_t i = 0; i < positions.size(); ++i) {
+            const Vec2& pos = positions[i];
+            if (visited.find(pos) == visited.end()) {
+                std::vector<int> component;
+                floodFill(pos, visited, component);
+                components.push_back(component);
+            }
+        }
+        
+        return components;
+    }
+    
+    // Getters
+    const std::vector<Vec2>& getPositions() const { return positions; }
+    const std::vector<Vec4>& getColors() const { return colors; }
+    const std::vector<float>& getSizes() const { return sizes; }
+    
+    Vec2 getPosition(int index) const { return positions[index]; }
+    Vec4 getColor(int index) const { return colors[index]; }
+    float getSize(int index) const { return sizes[index]; }
+    
+    void setColor(int index, const Vec4& color) { colors[index] = color; }
+    void setSize(int index, float size) { sizes[index] = size; }
+    
+    int getWidth() const { return width; }
+    int getHeight() const { return height; }
+
+private:
+    std::vector<Vec2> positions;
+    std::vector<Vec4> colors;
+    std::vector<float> sizes;
+    std::unordered_map<Vec2, int, std::hash<Vec2>> positionToIndex;
+    
+    void floodFill(const Vec2& start, std::unordered_map<Vec2, bool, std::hash<Vec2>>& visited, 
+                   std::vector<int>& component) const {
+        std::vector<Vec2> stack;
+        stack.push_back(start);
+        
+        while (!stack.empty()) {
+            Vec2 current = stack.back();
+            stack.pop_back();
+            
+            if (visited.find(current) != visited.end()) continue;
+            
+            visited[current] = true;
+            int index = getPixelIndex(current);
+            if (index != -1) {
+                component.push_back(index);
+                
+                // Add 4-connected neighbors
+                Vec2 neighbors[] = {current + Vec2(1, 0), current + Vec2(-1, 0), 
+                                   current + Vec2(0, 1), current + Vec2(0, -1)};
+                
+                for (const auto& neighbor : neighbors) {
+                    if (isOccupied(neighbor) && visited.find(neighbor) == visited.end()) {
+                        stack.push_back(neighbor);
+                    }
+                }
+            }
+        }
+    }
+};
+
+#endif

util/grid/grid3.hpp

@@ -0,0 +1,386 @@
+#ifndef GRID3_HPP
+#define GRID3_HPP
+
+#include "../vec3.hpp"
+#include "../vec4.hpp"
+#include <vector>
+#include <unordered_map>
+#include <string>
+#include <algorithm>
+#include <functional>
+
+class Grid3 {
+public:
+    // Constructors
+    Grid3() : width(0), height(0), depth(0) {}
+    Grid3(int size) : width(size), height(size), depth(size) {
+        positions.reserve(size * size * size);
+        colors.reserve(size * size * size);
+        sizes.reserve(size * size * size);
+    }
+    Grid3(int width, int height, int depth) : width(width), height(height), depth(depth) {
+        positions.reserve(width * height * depth);
+        colors.reserve(width * height * depth);
+        sizes.reserve(width * height * depth);
+    }
+    
+    // Add a voxel at specific position
+    int addVoxel(const Vec3& position, const Vec4& color, float size = 1.0f) {
+        int index = positions.size();
+        positions.push_back(position);
+        colors.push_back(color);
+        sizes.push_back(size);
+        positionToIndex[position] = index;
+        return index;
+    }
+    
+    // Add a voxel with integer coordinates
+    int addVoxel(int x, int y, int z, const Vec4& color, float size = 1.0f) {
+        return addVoxel(Vec3(static_cast<float>(x), static_cast<float>(y), static_cast<float>(z)), color, size);
+    }
+    
+    // Check if position is occupied
+    bool isOccupied(const Vec3& position) const {
+        return positionToIndex.find(position) != positionToIndex.end();
+    }
+    
+    bool isOccupied(int x, int y, int z) const {
+        return isOccupied(Vec3(static_cast<float>(x), static_cast<float>(y), static_cast<float>(z)));
+    }
+    
+    // Get voxel index at position, returns -1 if not found
+    int getVoxelIndex(const Vec3& position) const {
+        auto it = positionToIndex.find(position);
+        return (it != positionToIndex.end()) ? it->second : -1;
+    }
+    
+    int getVoxelIndex(int x, int y, int z) const {
+        return getVoxelIndex(Vec3(static_cast<float>(x), static_cast<float>(y), static_cast<float>(z)));
+    }
+    
+    // Remove voxel at position
+    bool removeVoxel(const Vec3& position) {
+        int index = getVoxelIndex(position);
+        if (index == -1) return false;
+        
+        // Swap with last element and update map
+        if (index != positions.size() - 1) {
+            positions[index] = positions.back();
+            colors[index] = colors.back();
+            sizes[index] = sizes.back();
+            
+            // Update mapping for the moved element
+            positionToIndex[positions[index]] = index;
+        }
+        
+        // Remove last element
+        positions.pop_back();
+        colors.pop_back();
+        sizes.pop_back();
+        positionToIndex.erase(position);
+        
+        return true;
+    }
+    
+    bool removeVoxel(int x, int y, int z) {
+        return removeVoxel(Vec3(static_cast<float>(x), static_cast<float>(y), static_cast<float>(z)));
+    }
+    
+    // Clear all voxels
+    void clear() {
+        positions.clear();
+        colors.clear();
+        sizes.clear();
+        positionToIndex.clear();
+    }
+    
+    // Get voxel count
+    size_t getVoxelCount() const {
+        return positions.size();
+    }
+    
+    // Check if grid is empty
+    bool isEmpty() const {
+        return positions.empty();
+    }
+    
+    // Get bounding box of occupied voxels
+    void getBoundingBox(Vec3& minCorner, Vec3& maxCorner) const {
+        if (positions.empty()) {
+            minCorner = Vec3(0, 0, 0);
+            maxCorner = Vec3(0, 0, 0);
+            return;
+        }
+        
+        minCorner = positions[0];
+        maxCorner = positions[0];
+        
+        for (const auto& pos : positions) {
+            minCorner = minCorner.min(pos);
+            maxCorner = maxCorner.max(pos);
+        }
+    }
+    
+    // Fill a rectangular prism region
+    void fillCuboid(const Vec3& start, const Vec3& end, const Vec4& color, float size = 1.0f) {
+        int startX = static_cast<int>(std::min(start.x, end.x));
+        int endX = static_cast<int>(std::max(start.x, end.x));
+        int startY = static_cast<int>(std::min(start.y, end.y));
+        int endY = static_cast<int>(std::max(start.y, end.y));
+        int startZ = static_cast<int>(std::min(start.z, end.z));
+        int endZ = static_cast<int>(std::max(start.z, end.z));
+        
+        for (int z = startZ; z <= endZ; ++z) {
+            for (int y = startY; y <= endY; ++y) {
+                for (int x = startX; x <= endX; ++x) {
+                    if (!isOccupied(x, y, z)) {
+                        addVoxel(x, y, z, color, size);
+                    }
+                }
+            }
+        }
+    }
+    
+    // Create a sphere
+    void fillSphere(const Vec3& center, float radius, const Vec4& color, float size = 1.0f) {
+        int centerX = static_cast<int>(center.x);
+        int centerY = static_cast<int>(center.y);
+        int centerZ = static_cast<int>(center.z);
+        int radiusInt = static_cast<int>(radius);
+        
+        for (int z = centerZ - radiusInt; z <= centerZ + radiusInt; ++z) {
+            for (int y = centerY - radiusInt; y <= centerY + radiusInt; ++y) {
+                for (int x = centerX - radiusInt; x <= centerX + radiusInt; ++x) {
+                    float dx = x - center.x;
+                    float dy = y - center.y;
+                    float dz = z - center.z;
+                    if (dx * dx + dy * dy + dz * dz <= radius * radius) {
+                        if (!isOccupied(x, y, z)) {
+                            addVoxel(x, y, z, color, size);
+                        }
+                    }
+                }
+            }
+        }
+    }
+    
+    // Create a hollow sphere (just the surface)
+    void fillHollowSphere(const Vec3& center, float radius, const Vec4& color, float thickness = 1.0f, float size = 1.0f) {
+        int centerX = static_cast<int>(center.x);
+        int centerY = static_cast<int>(center.y);
+        int centerZ = static_cast<int>(center.z);
+        int radiusInt = static_cast<int>(radius);
+        
+        for (int z = centerZ - radiusInt; z <= centerZ + radiusInt; ++z) {
+            for (int y = centerY - radiusInt; y <= centerY + radiusInt; ++y) {
+                for (int x = centerX - radiusInt; x <= centerX + radiusInt; ++x) {
+                    float dx = x - center.x;
+                    float dy = y - center.y;
+                    float dz = z - center.z;
+                    float distance = std::sqrt(dx * dx + dy * dy + dz * dz);
+                    
+                    if (std::abs(distance - radius) <= thickness) {
+                        if (!isOccupied(x, y, z)) {
+                            addVoxel(x, y, z, color, size);
+                        }
+                    }
+                }
+            }
+        }
+    }
+    
+    // Create a cylinder
+    void fillCylinder(const Vec3& baseCenter, const Vec3& axis, float radius, float height, 
+                     const Vec4& color, float size = 1.0f) {
+        // Simplified cylinder aligned with Y-axis
+        Vec3 normalizedAxis = axis.normalized();
+        
+        for (int h = 0; h < static_cast<int>(height); ++h) {
+            Vec3 center = baseCenter + normalizedAxis * static_cast<float>(h);
+            
+            for (int y = -static_cast<int>(radius); y <= static_cast<int>(radius); ++y) {
+                for (int x = -static_cast<int>(radius); x <= static_cast<int>(radius); ++x) {
+                    if (x * x + y * y <= radius * radius) {
+                        Vec3 pos = center + Vec3(x, y, 0);
+                        if (!isOccupied(pos)) {
+                            addVoxel(pos, color, size);
+                        }
+                    }
+                }
+            }
+        }
+    }
+    
+    // Get neighbors of a voxel (6-connected)
+    std::vector<int> getNeighbors6(const Vec3& position) const {
+        std::vector<int> neighbors;
+        Vec3 offsets[] = {
+            Vec3(1, 0, 0), Vec3(-1, 0, 0),
+            Vec3(0, 1, 0), Vec3(0, -1, 0),
+            Vec3(0, 0, 1), Vec3(0, 0, -1)
+        };
+        
+        for (const auto& offset : offsets) {
+            Vec3 neighborPos = position + offset;
+            int index = getVoxelIndex(neighborPos);
+            if (index != -1) {
+                neighbors.push_back(index);
+            }
+        }
+        
+        return neighbors;
+    }
+    
+    // Get neighbors of a voxel (26-connected)
+    std::vector<int> getNeighbors26(const Vec3& position) const {
+        std::vector<int> neighbors;
+        
+        for (int dz = -1; dz <= 1; ++dz) {
+            for (int dy = -1; dy <= 1; ++dy) {
+                for (int dx = -1; dx <= 1; ++dx) {
+                    if (dx == 0 && dy == 0 && dz == 0) continue;
+                    
+                    Vec3 neighborPos = position + Vec3(dx, dy, dz);
+                    int index = getVoxelIndex(neighborPos);
+                    if (index != -1) {
+                        neighbors.push_back(index);
+                    }
+                }
+            }
+        }
+        
+        return neighbors;
+    }
+    
+    // Create a simple teapot model (simplified)
+    void createTeapot(const Vec3& position, float scale, const Vec4& color, float size = 1.0f) {
+        // This is a very simplified teapot representation
+        // In practice, you'd load a proper voxel model
+        
+        // Teapot body (ellipsoid)
+        fillEllipsoid(position + Vec3(0, scale * 0.3f, 0), 
+                     Vec3(scale * 0.4f, scale * 0.3f, scale * 0.4f), color, size);
+        
+        // Teapot lid (smaller ellipsoid on top)
+        fillEllipsoid(position + Vec3(0, scale * 0.6f, 0), 
+                     Vec3(scale * 0.3f, scale * 0.1f, scale * 0.3f), color, size);
+        
+        // Teapot spout (cylinder)
+        fillCylinder(position + Vec3(scale * 0.3f, scale * 0.2f, 0), 
+                    Vec3(1, 0.2f, 0), scale * 0.05f, scale * 0.3f, color, size);
+        
+        // Teapot handle (torus segment)
+        fillTorusSegment(position + Vec3(-scale * 0.3f, scale * 0.3f, 0), 
+                        Vec3(0, 1, 0), scale * 0.1f, scale * 0.2f, color, size);
+    }
+    
+    // Fill an ellipsoid
+    void fillEllipsoid(const Vec3& center, const Vec3& radii, const Vec4& color, float size = 1.0f) {
+        int radiusX = static_cast<int>(radii.x);
+        int radiusY = static_cast<int>(radii.y);
+        int radiusZ = static_cast<int>(radii.z);
+        
+        for (int z = -radiusZ; z <= radiusZ; ++z) {
+            for (int y = -radiusY; y <= radiusY; ++y) {
+                for (int x = -radiusX; x <= radiusX; ++x) {
+                    float normalizedX = static_cast<float>(x) / radii.x;
+                    float normalizedY = static_cast<float>(y) / radii.y;
+                    float normalizedZ = static_cast<float>(z) / radii.z;
+                    
+                    if (normalizedX * normalizedX + normalizedY * normalizedY + normalizedZ * normalizedZ <= 1.0f) {
+                        Vec3 pos = center + Vec3(x, y, z);
+                        if (!isOccupied(pos)) {
+                            addVoxel(pos, color, size);
+                        }
+                    }
+                }
+            }
+        }
+    }
+    
+    // Fill a torus segment
+    void fillTorusSegment(const Vec3& center, const Vec3& axis, float majorRadius, float minorRadius, 
+                         const Vec4& color, float size = 1.0f) {
+        Vec3 normalizedAxis = axis.normalized();
+        
+        // Simplified torus - in practice you'd use proper torus equation
+        for (float angle = 0; angle < 2 * M_PI; angle += 0.2f) {
+            Vec3 circleCenter = center + Vec3(std::cos(angle) * majorRadius, 0, std::sin(angle) * majorRadius);
+            fillSphere(circleCenter, minorRadius, color, size);
+        }
+    }
+    
+    // Find connected components in 3D
+    std::vector<std::vector<int>> findConnectedComponents() const {
+        std::vector<std::vector<int>> components;
+        std::unordered_map<Vec3, bool, std::hash<Vec3>> visited;
+        
+        for (size_t i = 0; i < positions.size(); ++i) {
+            const Vec3& pos = positions[i];
+            if (visited.find(pos) == visited.end()) {
+                std::vector<int> component;
+                floodFill3D(pos, visited, component);
+                components.push_back(component);
+            }
+        }
+        
+        return components;
+    }
+    
+    // Getters
+    const std::vector<Vec3>& getPositions() const { return positions; }
+    const std::vector<Vec4>& getColors() const { return colors; }
+    const std::vector<float>& getSizes() const { return sizes; }
+    
+    Vec3 getPosition(int index) const { return positions[index]; }
+    Vec4 getColor(int index) const { return colors[index]; }
+    float getSize(int index) const { return sizes[index]; }
+    
+    void setColor(int index, const Vec4& color) { colors[index] = color; }
+    void setSize(int index, float size) { sizes[index] = size; }
+    
+    int getWidth() const { return width; }
+    int getHeight() const { return height; }
+    int getDepth() const { return depth; }
+
+private:
+    std::vector<Vec3> positions;
+    std::vector<Vec4> colors;
+    std::vector<float> sizes;
+    std::unordered_map<Vec3, int, std::hash<Vec3>> positionToIndex;
+    int width, height, depth;
+    
+    void floodFill3D(const Vec3& start, std::unordered_map<Vec3, bool, std::hash<Vec3>>& visited, 
+                     std::vector<int>& component) const {
+        std::vector<Vec3> stack;
+        stack.push_back(start);
+        
+        while (!stack.empty()) {
+            Vec3 current = stack.back();
+            stack.pop_back();
+            
+            if (visited.find(current) != visited.end()) continue;
+            
+            visited[current] = true;
+            int index = getVoxelIndex(current);
+            if (index != -1) {
+                component.push_back(index);
+                
+                // Add 6-connected neighbors
+                Vec3 neighbors[] = {
+                    current + Vec3(1, 0, 0), current + Vec3(-1, 0, 0),
+                    current + Vec3(0, 1, 0), current + Vec3(0, -1, 0),
+                    current + Vec3(0, 0, 1), current + Vec3(0, 0, -1)
+                };
+                
+                for (const auto& neighbor : neighbors) {
+                    if (isOccupied(neighbor) && visited.find(neighbor) == visited.end()) {
+                        stack.push_back(neighbor);
+                    }
+                }
+            }
+        }
+    }
+};
+
+#endif