pushing some junk home. adding back meshing

tests/g3etest.cpp

@@ -10,6 +10,7 @@
 #include <algorithm>
 
 #include "../util/grid/grid3eigen.hpp"
+#include "../util/grid/gridmesh.hpp"
 #include "../util/output/bmpwriter.hpp"
 #include "../util/output/frame.hpp"
 #include "../util/timing_decorator.cpp"
@@ -36,6 +37,9 @@ struct defaults {
     int lodDist = 500;
     float lodDropoff = 0.1;
     PNoise2 noise = PNoise2(42);
+    
+    float meshResolution = 0.5f;
+    float meshIsoLevel = 0.4f;
 };
 
 struct spheredefaults {
@@ -85,6 +89,9 @@ const std::chrono::seconds STATS_UPDATE_INTERVAL(60);
 std::string cachedStats;
 bool statsNeedUpdate = true;
 
+CachedGridMesh<int> planetMesh(1); 
+bool meshNeedsUpdate = false;
+
 void createSphere(const defaults& config, const spheredefaults& sconfig, Octree<int>& grid) {
     if (!grid.empty()) grid.clear();
 
@@ -167,6 +174,8 @@ void createSphere(const defaults& config, const spheredefaults& sconfig, Octree<
             }
         }
     }
+    std::cout << "voxel grid ready" << std::endl;
+    meshNeedsUpdate = true;
 }
 
 void addStar(const defaults& config, const stardefaults& starconf, Octree<int>& grid) {
@@ -176,6 +185,8 @@ void addStar(const defaults& config, const stardefaults& starconf, Octree<int>&
     PointType pos(starconf.x, starconf.y, starconf.z);
 
     grid.set(2, pos, true, colorVec, starconf.size, true, 2, true, starconf.emittance, 0.0f, 0.0f);
+    
+    meshNeedsUpdate = true;
 }
 
 void updateStatsCache(Octree<int>& grid) {
@@ -190,16 +201,26 @@ void livePreview(Octree<int>& grid, defaults& config, const Camera& cam) {
     std::lock_guard<std::mutex> lock(PreviewMutex);
     updatePreview = true;
     
+    if (meshNeedsUpdate) {
+        planetMesh.setResolution(config.meshResolution);
+        planetMesh.setIsoLevel(config.meshIsoLevel);
+        planetMesh.update(grid);
+        meshNeedsUpdate = false;
+        statsNeedUpdate = true;
+    }
+
     auto renderStart = std::chrono::high_resolution_clock::now();
     frame currentPreviewFrame;
-    grid.setLODMinDistance(config.lodDist);
+
-    grid.setLODFalloff(config.lodDropoff);
+    currentPreviewFrame = planetMesh.render(cam, config.outWidth, config.outHeight, 0.1f, 10000.0f, frame::colormap::RGB);
-    if (config.slowRender) {
-        currentPreviewFrame = grid.renderFrame(cam, config.outWidth, config.outHeight, frame::colormap::RGB, config.rayCount, config.reflectCount, config.globalIllumination, config.useLod);
-    } else {
-        currentPreviewFrame = grid.fastRenderFrame(cam, config.outWidth, config.outHeight, frame::colormap::RGB);
-    }
     
+    // grid.setLODMinDistance(config.lodDist);
+    // grid.setLODFalloff(config.lodDropoff);
+    // if (config.slowRender) {
+    //     currentPreviewFrame = grid.renderFrame(cam, config.outWidth, config.outHeight, frame::colormap::RGB, config.rayCount, config.reflectCount, config.globalIllumination, config.useLod);
+    // } else {
+    //     currentPreviewFrame = grid.fastRenderFrame(cam, config.outWidth, config.outHeight, frame::colormap::RGB);
+    // }
     
     auto renderEnd = std::chrono::high_resolution_clock::now();
     renderFrameTime = std::chrono::duration<double>(renderEnd - renderStart).count();
@@ -294,7 +315,7 @@ int main() {
     #endif
 
     bool application_not_closed = true;
-    GLFWwindow* window = glfwCreateWindow((int)(1280), (int)(800), "voxelgrid live renderer", nullptr, nullptr);
+    GLFWwindow* window = glfwCreateWindow((int)(1280), (int)(800), "StupidSim", nullptr, nullptr);
     if (window == nullptr) {
         glfwTerminate();
         return 1;
@@ -355,7 +376,7 @@ int main() {
     float rotationSpeedZ = 0.05f;
     float autoRotationTime = 0.0f;
     PointType initialViewDir(0, 0, 1);
-    float rotationRadius = 255.0f;
+    float rotationRadius = 2000.0f;
     float yawSpeed = 0.5f;
     float pitchSpeed = 0.3f;
     float rollSpeed = 0.2f;
@@ -384,6 +405,7 @@ int main() {
         gridInitialized = true;
         updateStatsCache(grid);
         resetView(cam, config.gridSizecube);
+        meshNeedsUpdate = true;
     }
     
     while (!glfwWindowShouldClose(window)) {
@@ -490,20 +512,29 @@ int main() {
                 sphereConf.centerZ = pos[2];
             }
             ImGui::DragFloat("Radius", &sphereConf.radius, 0.5f, 1.0f, 250.0f);
-            ImGui::DragFloat("Density (Overlap)", &sphereConf.voxelSize, 0.05f, 0.1f, 5.0f);
+            ImGui::DragFloat("Density (Overlap)", &sphereConf.voxelSize, 0.1f, 1.0f, 100.0f);
             if (ImGui::IsItemHovered()) {
                 ImGui::SetTooltip("Multiplies calculated point size. >1.0 ensures solid surface.");
             }
             ImGui::ColorEdit3("Color", sphereConf.color);
             
             ImGui::Separator();
-            ImGui::Checkbox("Is Light", &sphereConf.light);
+            ImGui::Text("Marching Cubes Config");
-            if(sphereConf.light) {
+            if (ImGui::SliderFloat("Mesh Resolution", &config.meshResolution, 0.1f, 2.0f)) {
-                ImGui::DragFloat("Emittance", &sphereConf.emittance, 0.1f, 0.0f, 100.0f);
+                meshNeedsUpdate = true;
             }
-            ImGui::SliderFloat("Reflection", &sphereConf.reflection, 0.0f, 1.0f);
+            if (ImGui::SliderFloat("Iso Level", &config.meshIsoLevel, 0.01f, 1.0f)) {
-            ImGui::SliderFloat("Refraction", &sphereConf.refraction, 0.0f, 1.0f);
+                meshNeedsUpdate = true;
-            ImGui::Checkbox("Fill Inside", &sphereConf.fillInside);
+            }
+
+            ImGui::Separator();
+            ImGui::Checkbox("Is Light", &sphereConf.light);
+            // if(sphereConf.light) {
+            //     ImGui::DragFloat("Emittance", &sphereConf.emittance, 0.1f, 0.0f, 100.0f);
+            // }
+            // ImGui::SliderFloat("Reflection", &sphereConf.reflection, 0.0f, 1.0f);
+            // ImGui::SliderFloat("Refraction", &sphereConf.refraction, 0.0f, 1.0f);
+            // ImGui::Checkbox("Fill Inside", &sphereConf.fillInside);
 
             if (ImGui::CollapsingHeader("Star/Sun Parameters", ImGuiTreeNodeFlags_DefaultOpen)) {
                 ImGui::Checkbox("Enable Star", &starConf.enabled);
@@ -516,7 +547,7 @@ int main() {
                     starConf.z = starPos[2];
                 }
 
-                ImGui::DragFloat("Size (Radius)", &starConf.size, 1.0f, 1.0f, 500.0f);
+                ImGui::DragFloat("Size (Radius)", &starConf.size, 1.0f, 1.0f, 1000.0f);
                 ImGui::DragFloat("Brightness", &starConf.emittance, 1.0f, 0.0f, 1000.0f);
                 ImGui::ColorEdit3("Light Color", starConf.color);
             }
@@ -604,16 +635,15 @@ int main() {
             float maxSliderValueX = config.gridSizecube;
             float maxSliderValueY = config.gridSizecube;
             float maxSliderValueZ = config.gridSizecube;
-            float maxSliderValueRotation = 360.0f;
             
-            ImGui::Text("Position (0 to grid size²):");
+            ImGui::Text("Position (0 to grid size):");
-            if (ImGui::SliderFloat("Camera X", &camX, 0.0f, maxSliderValueX)) {
+            if (ImGui::SliderFloat("Camera X", &camX, -maxSliderValueX, maxSliderValueX)) {
                 cam.origin[0] = camX;
             }
-            if (ImGui::SliderFloat("Camera Y", &camY, 0.0f, maxSliderValueY)) {
+            if (ImGui::SliderFloat("Camera Y", &camY, -maxSliderValueY, maxSliderValueY)) {
                 cam.origin[1] = camY;
             }
-            if (ImGui::SliderFloat("Camera Z", &camZ, 0.0f, maxSliderValueZ)) {
+            if (ImGui::SliderFloat("Camera Z", &camZ, -maxSliderValueZ, maxSliderValueZ)) {
                 cam.origin[2] = camZ;
             }
             
@@ -629,13 +659,12 @@ int main() {
                 cam.direction[2] = camvZ;
             }
 
-            if (ImGui::SliderFloat("Camera Speed", &camspeed, 1, 100)) {
+            if (ImGui::SliderFloat("Camera Speed", &camspeed, 1, 500)) {
                 cam.movementSpeed = camspeed;
             }
             
             ImGui::Separator();
             
-            // Focus Button
             if (ImGui::Button("Focus on Planet")) {
                 PointType target(sphereConf.centerX, sphereConf.centerY, sphereConf.centerZ);
                 PointType newDir = (target - cam.origin).normalized();
@@ -672,10 +701,7 @@ int main() {
 
             ImGui::Separator();
             ImGui::Text("Current Camera Position:");
-            ImGui::Text("X: %.2f, Y: %.2f, Z: %.2f", 
+            ImGui::Text("X: %.2f, Y: %.2f, Z: %.2f", cam.origin[0], cam.origin[1], cam.origin[2]);
-                        cam.origin[0], 
-                        cam.origin[1], 
-                        cam.origin[2]);
 
             ImGui::Text("Auto-Rotation:");
             
@@ -727,7 +753,8 @@ int main() {
                 // Update camera
                 cam.origin = PointType(camX, camY, camZ);
                 cam.direction = PointType(camvX, camvY, camvZ);
-                
+                cam.lookAt(PointType(sphereConf.centerX, sphereConf.centerY, sphereConf.centerZ));
+
                 // Sliders to control rotation speeds
                 ImGui::SliderFloat("X Speed", &rotationSpeedX, 0.01f, 1.0f);
                 ImGui::SliderFloat("Y Speed", &rotationSpeedY, 0.01f, 1.0f);
@@ -735,7 +762,7 @@ int main() {
             }
             
             // Slider for orbit radius
-            ImGui::SliderFloat("Orbit Radius", &rotationRadius, 10.0f, 2000.0f);
+            ImGui::SliderFloat("Orbit Radius", &rotationRadius, 10.0f, 5000.0f);
 
             if (autoRotateView) {
                 ImGui::SameLine();
@@ -791,16 +818,19 @@ int main() {
 
             ImGui::Separator();
 
-            ImGui::Checkbox("update Preview", &worldPreview);
+            ImGui::Checkbox("Continuous Preview", &worldPreview);
-            ImGui::Checkbox("Use Slower renderer", &config.slowRender);
+
-            if (config.slowRender) {
+            // Removed Raytracing specific controls (Global Illum, LOD) as they don't apply to raster mesh
-                ImGui::InputInt("Rays per pixel", &config.rayCount);
+            // ImGui::Checkbox("update Preview", &worldPreview);
-                ImGui::InputInt("Max reflections", &config.reflectCount);
+            // ImGui::Checkbox("Use Slower renderer", &config.slowRender);
-            }
+            // if (config.slowRender) {
-            ImGui::InputFloat("Lod dropoff", &config.lodDropoff);
+            //     ImGui::InputInt("Rays per pixel", &config.rayCount);
-            ImGui::InputInt("lod minimum Distance", &config.lodDist);
+            //     ImGui::InputInt("Max reflections", &config.reflectCount);
-            ImGui::Checkbox("use Global illumination", &config.globalIllumination);
+            // }
-            ImGui::Checkbox("use Lod", &config.useLod);
+            // ImGui::InputFloat("Lod dropoff", &config.lodDropoff);
+            // ImGui::InputInt("lod minimum Distance", &config.lodDist);
+            // ImGui::Checkbox("use Global illumination", &config.globalIllumination);
+            // ImGui::Checkbox("use Lod", &config.useLod);
 
             ImGui::End();
         }

util/grid/grid3eigen.hpp

@@ -211,7 +211,6 @@ private:
 
         auto accumulate = [&](const std::shared_ptr<NodeData>& item) {
             if (!item || !item->active || !item->visible) return;
-            avgPos += item->position;
             avgColor += item->color;
             avgEmittance += item->emittance;
             avgRefl += item->reflection;
@@ -236,7 +235,7 @@ private:
             
             auto lod = std::make_shared<NodeData>();
             
-            lod->position = avgPos * invCount;
+            lod->position = node->center;
             lod->color = avgColor * invCount;
             
             PointType nodeDims = node->bounds.second - node->bounds.first;
@@ -825,8 +824,6 @@ public:
     std::vector<std::shared_ptr<NodeData>> voxelTraverse(const PointType& origin, const PointType& direction,
                                          float maxDist, bool stopAtFirstHit, bool enableLOD = false) const {
         std::vector<std::shared_ptr<NodeData>> hits;
-        
-        if (empty()) return hits;
         float invLodf = 1.0f / lodFalloffRate_;
         uint8_t raySignMask = (direction.x() < 0 ? 1 : 0) | (direction.y() < 0 ? 2 : 0) | (direction.z() < 0 ? 4 : 0);
         Ray oray(origin, direction);
@@ -1014,30 +1011,27 @@ public:
         return outFrame;
     }
 
-    std::shared_ptr<NodeData> fastVoxelTraverse(const PointType& origin, const PointType& direction,
+    std::shared_ptr<NodeData> fastVoxelTraverse(const Ray& ray, float maxDist, bool enableLOD = false) const {
-                                         float maxDist, bool enableLOD = false) const {
         std::shared_ptr<NodeData> hit;
-        Ray oray(origin, direction);
         if (empty()) return hit;
+        float lodRatio = 1.0f / lodFalloffRate_;
 
         std::function<void(OctreeNode*, const PointType&, const PointType&, float, float)> traverseNode = 
                         [&](OctreeNode* node, const PointType& origin, const PointType& dir, float tMin, float tMax) {
             if (!node || tMin > tMax) return;
-            Ray ray(origin, dir);
 
             // LOD Check for fast traverse
             if (enableLOD && !node->isLeaf) {
                 float dist = (node->center - origin).norm();
                 float nodeSize = (node->bounds.second - node->bounds.first).norm();
-                if (dist > lodMinDistance_) {
+                if (dist > lodMinDistance_ && dist <= maxDist) {
                     float ratio = dist / (nodeSize + EPSILON);
-                    if (ratio > (1.0f / lodFalloffRate_)) {
+                    if (ratio > lodRatio) {
                         ensureLOD(node);
                         if (node->lodData) {
-                             // Project LOD
                             PointType toPoint = node->lodData->position - origin;
                             float projection = toPoint.dot(dir);
-                            if (projection >= 0 && projection <= maxDist) {
+                            if (projection >= 0) {
                                 PointType closestPoint = origin + dir * projection;
                                 float distSq = (node->lodData->position - closestPoint).squaredNorm();
                                 if (distSq < node->lodData->size * node->lodData->size) {
@@ -1096,9 +1090,9 @@ public:
             }
         };
         float tMin, tMax;
-        if (rayBoxIntersect(oray, root_->bounds, tMin, tMax)) {
+        if (rayBoxIntersect(ray, root_->bounds, tMin, tMax)) {
             tMax = std::min(tMax, maxDist);
-            traverseNode(root_.get(), origin, direction, tMin, tMax);
+            traverseNode(root_.get(), ray.origin, ray.dir, tMin, tMax);
         }
         return hit;
     }
@@ -1134,8 +1128,8 @@ public:
                 rayDir.normalize();
 
                 Eigen::Vector3f color = backgroundColor_;
-                
+                Ray ray(origin, rayDir);
-                auto hit = fastVoxelTraverse(origin, rayDir, std::numeric_limits<float>::max(), true);
+                auto hit = fastVoxelTraverse(ray, std::numeric_limits<float>::max(), true);
                 if (hit != nullptr) {
                     auto obj = hit;
                     

util/grid/gridmesh.hpp

@@ -0,0 +1,253 @@
+#ifndef GRIDMESH_HPP
+#define GRIDMESH_HPP
+
+#include "grid3eigen.hpp"
+#include "mesh.hpp"
+#include "../basicdefines.hpp"
+#include <vector>
+#include <limits>
+#include <memory>
+#include <unordered_map>
+#include <map>
+#include <cmath>
+
+struct GridPoint {
+    Eigen::Vector3f pos;
+    float val;
+    Eigen::Vector3f color; // Interpolated color field
+};
+
+struct GridCell {
+    GridPoint p[8];
+};
+
+template<typename T>
+class MetaballMesher {
+private:
+    float _resolution; 
+    float _isoLevel;
+    float _influenceRadius;
+
+    void vertexInterp(float isolevel, GridPoint p1, GridPoint p2, Eigen::Vector3f& outPos, Eigen::Vector3f& outColor) {
+        if (std::abs(p1.val - p2.val) < 0.00001f) {
+            outPos = p1.pos;
+            outColor = p1.color;
+            return;
+        }
+
+        float mu = (isolevel - p1.val) / (p2.val - p1.val);
+        outPos = p1.pos + mu * (p2.pos - p1.pos);
+        outColor = p1.color + mu * (p2.color - p1.color);
+    }
+
+    void sampleField(const Eigen::Vector3f& pos, const std::vector<std::shared_ptr<typename Octree<T>::NodeData>>& nodes, float& outVal, Eigen::Vector3f& outColor) {
+        float sumVal = 0.0f;
+        Eigen::Vector3f sumColor = Eigen::Vector3f::Zero();
+        float totalWeight = 0.0f;
+        
+        for (const auto& node : nodes) {
+            float r = node->size * _influenceRadius;
+            float distSq = (pos - node->position).squaredNorm();
+            
+            if (distSq < r * r) {
+                float dist = std::sqrt(distSq);
+                float x = dist / r;
+                float val = std::pow(1.0f - x*x, 3.0f);
+                
+                sumVal += val;
+                sumColor += node->color * val;
+                totalWeight += val;
+            }
+        }
+
+        outVal = sumVal;
+        if (totalWeight > 0.0001f) {
+            outColor = sumColor / totalWeight;
+        } else {
+            outColor = Eigen::Vector3f(1.0f, 0.0f, 1.0f);
+        }
+    }
+
+public:
+    MetaballMesher(float resolution = 0.5f, float isoLevel = 0.5f, float influenceMult = 2.5f) 
+        : _resolution(resolution), _isoLevel(isoLevel), _influenceRadius(influenceMult) {}
+
+    void polygonize(const std::vector<std::shared_ptr<typename Octree<T>::NodeData>>& nodes, 
+                   std::vector<Eigen::Vector3f>& outVerts, 
+                   std::vector<Eigen::Vector3f>& outColors,
+                   std::vector<int>& outIndices,
+                   int& startIndex) {
+        
+        if (nodes.empty()) return;
+
+        Eigen::Vector3f minBounds(1e9, 1e9, 1e9);
+        Eigen::Vector3f maxBounds(-1e9, -1e9, -1e9);
+        float maxNodeSize = 0;
+
+        for (const auto& node : nodes) {
+            minBounds = minBounds.cwiseMin(node->position);
+            maxBounds = maxBounds.cwiseMax(node->position);
+            if (node->size > maxNodeSize) maxNodeSize = node->size;
+        }
+
+        float padding = maxNodeSize * _influenceRadius;
+        minBounds -= Eigen::Vector3f(padding, padding, padding);
+        maxBounds += Eigen::Vector3f(padding, padding, padding);
+
+        float step = maxNodeSize * _resolution;
+        if (step <= 0.0001f) step = 0.1f;
+
+        int stepsX = static_cast<int>((maxBounds.x() - minBounds.x()) / step);
+        int stepsY = static_cast<int>((maxBounds.y() - minBounds.y()) / step);
+        int stepsZ = static_cast<int>((maxBounds.z() - minBounds.z()) / step);
+        
+        for (int x = 0; x < stepsX; ++x) {
+            for (int y = 0; y < stepsY; ++y) {
+                for (int z = 0; z < stepsZ; ++z) {
+                    
+                    Eigen::Vector3f basePos = minBounds + Eigen::Vector3f(x*step, y*step, z*step);
+                    GridCell cell;
+                    
+                    Eigen::Vector3f offsets[8] = {
+                        {0,0,0}, {step,0,0}, {step,0,step}, {0,0,step},
+                        {0,step,0}, {step,step,0}, {step,step,step}, {0,step,step}
+                    };
+
+                    int cubeIndex = 0;
+                    for (int i = 0; i < 8; ++i) {
+                        cell.p[i].pos = basePos + offsets[i];
+                        sampleField(cell.p[i].pos, nodes, cell.p[i].val, cell.p[i].color);
+                        if (cell.p[i].val < _isoLevel) cubeIndex |= (1 << i);
+                    }
+
+                    // Look up edges
+                    if (edgeTable[cubeIndex] == 0) continue;
+                    if (edgeTable[cubeIndex] == 255) continue; // Inside or Outside completely
+
+                    Eigen::Vector3f vertList[12];
+                    Eigen::Vector3f colorList[12];
+
+                    if (edgeTable[cubeIndex] & 1) vertexInterp(_isoLevel, cell.p[0], cell.p[1], vertList[0], colorList[0]);
+                    if (edgeTable[cubeIndex] & 2) vertexInterp(_isoLevel, cell.p[1], cell.p[2], vertList[1], colorList[1]);
+                    if (edgeTable[cubeIndex] & 4) vertexInterp(_isoLevel, cell.p[2], cell.p[3], vertList[2], colorList[2]);
+                    if (edgeTable[cubeIndex] & 8) vertexInterp(_isoLevel, cell.p[3], cell.p[0], vertList[3], colorList[3]);
+                    if (edgeTable[cubeIndex] & 16) vertexInterp(_isoLevel, cell.p[4], cell.p[5], vertList[4], colorList[4]);
+                    if (edgeTable[cubeIndex] & 32) vertexInterp(_isoLevel, cell.p[5], cell.p[6], vertList[5], colorList[5]);
+                    if (edgeTable[cubeIndex] & 64) vertexInterp(_isoLevel, cell.p[6], cell.p[7], vertList[6], colorList[6]);
+                    if (edgeTable[cubeIndex] & 128) vertexInterp(_isoLevel, cell.p[7], cell.p[4], vertList[7], colorList[7]);
+                    if (edgeTable[cubeIndex] & 256) vertexInterp(_isoLevel, cell.p[0], cell.p[4], vertList[8], colorList[8]);
+                    if (edgeTable[cubeIndex] & 512) vertexInterp(_isoLevel, cell.p[1], cell.p[5], vertList[9], colorList[9]);
+                    if (edgeTable[cubeIndex] & 1024) vertexInterp(_isoLevel, cell.p[2], cell.p[6], vertList[10], colorList[10]);
+                    if (edgeTable[cubeIndex] & 2048) vertexInterp(_isoLevel, cell.p[3], cell.p[7], vertList[11], colorList[11]);
+
+                    for (int i = 0; triTable[cubeIndex][i] != -1; i += 3) {
+                        outVerts.push_back(vertList[triTable[cubeIndex][i]]);
+                        outVerts.push_back(vertList[triTable[cubeIndex][i+1]]);
+                        outVerts.push_back(vertList[triTable[cubeIndex][i+2]]);
+                        
+                        outColors.push_back(colorList[triTable[cubeIndex][i]]);
+                        outColors.push_back(colorList[triTable[cubeIndex][i+1]]);
+                        outColors.push_back(colorList[triTable[cubeIndex][i+2]]);
+
+                        outIndices.push_back(startIndex++);
+                        outIndices.push_back(startIndex++);
+                        outIndices.push_back(startIndex++);
+                    }
+                }
+            }
+        }
+    }
+};
+
+template<typename T>
+class CachedGridMesh {
+private:
+    std::unique_ptr<Mesh> mesh_;
+    int id_;
+    
+    float _resolution = 0.5f;
+    float _isoLevel = 0.4f;
+
+public:
+    CachedGridMesh(int id = 0) : id_(id) {
+        mesh_ = std::make_unique<Mesh>(id_, std::vector<Vector3f>{}, std::vector<std::vector<int>>{}, std::vector<Color>{});
+    }
+
+    void setResolution(float res) { _resolution = res; }
+    void setIsoLevel(float iso) { _isoLevel = iso; }
+
+    void update(Octree<T>& grid, const Eigen::Vector3f& center = Eigen::Vector3f::Zero(), float radius = std::numeric_limits<float>::max()) {
+        std::vector<Vector3f> allVertices;
+        std::vector<std::vector<int>> allPolys;
+        std::vector<Color> allColors;
+        
+        auto nodes = grid.findInRadius(center, radius);
+        if (nodes.empty()) {
+            mesh_ = std::make_unique<Mesh>(id_, allVertices, allPolys, allColors);
+            return;
+        }
+        std::unordered_map<int, std::vector<std::shared_ptr<typename Octree<T>::NodeData>>> objectGroups;
+        
+        for (const auto& node : nodes) {
+            if (!node->active || !node->visible) continue;
+            objectGroups[node->objectId].push_back(node);
+        }
+        std::cout << "object map made" << std::endl;
+
+        MetaballMesher<T> mesher(_resolution, _isoLevel);
+        
+        int globalVertIndex = 0;
+        
+        for (auto& [objId, groupNodes] : objectGroups) {
+            std::vector<Eigen::Vector3f> objVerts;
+            std::vector<Eigen::Vector3f> objColors;
+            std::vector<int> objIndices;
+            int startIndex = globalVertIndex;
+            int localVertCounter = 0; 
+            
+            std::vector<Eigen::Vector3f> mVerts;
+            std::vector<Eigen::Vector3f> mColors;
+            std::vector<int> mIndices;
+            int mIdxStart = 0;
+
+            mesher.polygonize(groupNodes, mVerts, mColors, mIndices, mIdxStart);
+
+            for (size_t i = 0; i < mVerts.size(); ++i) {
+                allVertices.push_back(mVerts[i]);
+                allColors.push_back(mColors[i]);
+            }
+
+            for (size_t i = 0; i < mIndices.size(); i += 3) {
+                std::vector<int> tri = {
+                    mIndices[i] + globalVertIndex, 
+                    mIndices[i+1] + globalVertIndex, 
+                    mIndices[i+2] + globalVertIndex
+                };
+                allPolys.push_back(tri);
+            }
+
+            globalVertIndex += mVerts.size();
+            std::cout << "object done" << std::endl;
+        }
+
+        mesh_ = std::make_unique<Mesh>(id_, allVertices, allPolys, allColors);
+        mesh_->triangulate(); 
+    }
+
+    frame render(Camera cam, int height, int width, float near = 0.1f, float far = 1000.0f, frame::colormap format = frame::colormap::RGB) {
+        if (!mesh_) {
+            return frame(width, height, format);
+        }
+        return mesh_->renderFrame(cam, height, width, near, far, format);
+    }
+
+    Mesh* getMesh() {
+        return mesh_.get();
+    }
+    
+    size_t getVertexCount() const {
+        return mesh_ ? mesh_->vertices().size() : 0;
+    }
+};
+
+#endif

util/grid/mesh.hpp

@@ -12,6 +12,7 @@
 #include "../../eigen/Eigen/Dense"
 #include "../../eigen/Eigen/Geometry"
 #include "./camera.hpp"
+#include "../output/frame.hpp"
 
 using Vector2f = Eigen::Vector2f;
 using Vector3f = Eigen::Vector3f;
@@ -19,6 +20,12 @@ using Vector4f = Eigen::Vector4f;
 using Matrix4f = Eigen::Matrix4f;
 using Color = Eigen::Vector3f;
 
+struct Triangle2D {
+    Vector2f a, b, c;
+    Color color;
+    float depth;
+};
+
 class Mesh {
 private:
     int id;
@@ -27,6 +34,10 @@ private:
     std::vector<Color> _colors;
     mutable std::vector<Eigen::Vector3i> _tris;
     mutable bool _needs_triangulation = true;
+    
+    inline static float edgeFunction(const Vector2f& a, const Vector2f& b, const Vector2f& c) {
+        return (c.x() - a.x()) * (b.y() - a.y()) - (c.y() - a.y()) * (b.x() - a.x());
+    }
 public:
     Mesh(int id, const std::vector<Vector3f>& verts, const std::vector<std::vector<int>>& polys, const std::vector<Color>& colors)
         : id(id), _vertices(verts), _polys(polys), _colors(colors) {}
@@ -53,12 +64,12 @@ public:
     }
 
     void triangulate() {
-        std::vector<Eigen::Vector3i> newPols;
         if (!_needs_triangulation) return;
+        std::vector<Eigen::Vector3i> newPols;
         for (auto& pol : _polys) {
             if (pol.size() > 3) {
                 auto v0 = pol[0];
-                for (int i = 0; i < pol.size() - 1; i++) {
+                for (size_t i = 1; i < pol.size() - 1; i++) {
                     newPols.emplace_back(Eigen::Vector3i{v0, pol[i], pol[i+1]});
                 }
             } else if (pol.size() == 3) {
@@ -75,45 +86,188 @@ public:
     }
 
     bool colors(std::vector<Color> colorlist) {
-        if (colorlist.size() == 1) {
+        if (colorlist.size() == 1 || colorlist.size() == _vertices.size()) {
-            _colors = colorlist;
-        } else if (colorlist.size() == _vertices.size()) {
             _colors = colorlist;
+            return true;
         }
+        return false;
     }
 
-    void project_2d(Camera cam, int height, int width, float near, float far) {
+    std::vector<Triangle2D> project_2d(Camera cam, int height, int width, float near, float far) {
-        Vector3f zaxis = cam.direction - cam.origin;
+        triangulate();
-        zaxis = zaxis / zaxis.norm();
-        Vector3f xAxis = cam.up.cross(zaxis);
-        xAxis = xAxis / xAxis.norm();
-        Vector3f yAxis = zaxis.cross(xAxis);
 
-        Matrix4f viewMatrix = Matrix4f::Identity();
+        std::vector<Triangle2D> renderList;
-        viewMatrix.block<3,1>(0,0) = xAxis;
-        viewMatrix.block<3,1>(0,1) = yAxis;
-        viewMatrix.block<3,1>(0,2) = -zaxis;
-        viewMatrix.block<3,1>(0,3) = cam.origin;
 
-        viewMatrix = viewMatrix.inverse();
+        Vector3f forward = cam.forward();
+        Vector3f right = cam.right();
+        Vector3f up = cam.up;
 
-        float aspect = float(height) / float(width);
+        Matrix4f viewMatrixa = Matrix4f::Identity();
+        viewMatrixa.block<3,1>(0,0) = right;
+        viewMatrixa.block<3,1>(0,1) = up;
+        viewMatrixa.block<3,1>(0,2) = -forward; 
+        viewMatrixa.block<3,1>(0,3) = cam.origin;
+        
+        Matrix4f viewMatrix = viewMatrixa.inverse();
+
+        float aspect = float(width) / float(height);
         float fovrad = cam.fovRad();
-        float tanh = std::tan(fovrad / 2);
+        float tanHalfFov = std::tan(fovrad / 2.0f);
         
         Matrix4f projMatrix = Matrix4f::Zero();
-        projMatrix(0,0) = 1.0 / (aspect * tanh);
+        projMatrix(0,0) = 1.0f / (aspect * tanHalfFov);
-        projMatrix(1,1) = 1.0 / tanh;
+        projMatrix(1,1) = 1.0f / tanHalfFov;
-        projMatrix(2,2) = -(far + near) / (near - far);
+        projMatrix(2,2) = -(far + near) / (far - near);
-        projMatrix(2,3) = (-2.0f * far * near) / (near - far);
+        projMatrix(2,3) = -(2.0f * far * near) / (far - near);
         projMatrix(3,2) = -1.0f;
 
-        Vector3f viewDir = (cam.direction - cam.origin).normalized();
+        int n = _vertices.size();
+        
+        std::vector<Vector2f> screenCoords(n);
+        std::vector<float> linearDepths(n);
+        std::vector<bool> validVerts(n, false);
 
-        Vector3f mesh_center = _vertices.colwise().mean(); 
+        Eigen::MatrixXf homogeneousVerts(n, 4);
+        for (int i = 0; i < n; ++i) {
+            homogeneousVerts.row(i).head<3>() = _vertices[i];
+            homogeneousVerts(i, 3) = 1.0f;
+        }
 
+        Eigen::MatrixXf viewVerts = homogeneousVerts * viewMatrix.transpose();
+        Eigen::MatrixXf clipVerts = viewVerts * projMatrix.transpose();
+
+        for (int i = 0; i < n; ++i) {
+            float w = clipVerts(i, 3);
+            
+            if (w <= near) {
+                validVerts[i] = false;
+                continue; 
+            }
+
+            float x_ndc = clipVerts(i, 0) / w;
+            float y_ndc = clipVerts(i, 1) / w;
+
+            screenCoords[i].x() = (x_ndc + 1.0f) * 0.5f * width;
+            screenCoords[i].y() = (1.0f - (y_ndc + 1.0f) * 0.5f) * height;
+            
+            linearDepths[i] = w; 
+            validVerts[i] = true;
+        }
+
+        for (const auto& triIdx : _tris) {
+            int i0 = triIdx.x();
+            int i1 = triIdx.y();
+            int i2 = triIdx.z();
+
+            if (!validVerts[i0] || !validVerts[i1] || !validVerts[i2]) {
+                continue;
+            }
+
+            Triangle2D t2d;
+            t2d.a = screenCoords[i0];
+            t2d.b = screenCoords[i1];
+            t2d.c = screenCoords[i2];
+
+            t2d.depth = (linearDepths[i0] + linearDepths[i1] + linearDepths[i2]) / 3.0f;
+
+            // Handle Coloring
+            if (_colors.size() == n) {
+                t2d.color = (_colors[i0] + _colors[i1] + _colors[i2]) / 3.0f;
+            } else if (_colors.size() == 1) {
+                t2d.color = _colors[0];
+            } else {
+                t2d.color = {1.0f, 0.0f, 1.0f};
+            }
+
+            renderList.push_back(t2d);
+        }
+
+        std::sort(renderList.begin(), renderList.end(), [](const Triangle2D& a, const Triangle2D& b) {
+            return a.depth > b.depth;
+        });
+
+        return renderList;
+    }
+    
+    frame renderFrame(Camera cam, int height, int width, float near, float far, frame::colormap colorformat = frame::colormap::RGB) {
+        std::vector<Triangle2D> triangles = project_2d(cam, height, width, near, far);
+        return rasterizeTriangles(triangles, width, height, colorformat);
     }
 
+    static frame rasterizeTriangles(const std::vector<Triangle2D>& triangles, int width, int height, frame::colormap colorformat) {
+        frame outFrame(width, height, colorformat);
+        std::vector<float> buffer(width * height * 3);
+        
+        for (const auto& tri : triangles) {
+            int minX = std::max(0, (int)std::floor(std::min({tri.a.x(), tri.b.x(), tri.c.x()})));
+            int maxX = std::min(width - 1, (int)std::ceil(std::max({tri.a.x(), tri.b.x(), tri.c.x()})));
+            int minY = std::max(0, (int)std::floor(std::min({tri.a.y(), tri.b.y(), tri.c.y()})));
+            int maxY = std::min(height - 1, (int)std::ceil(std::max({tri.a.y(), tri.b.y(), tri.c.y()})));
+
+            float area = edgeFunction(tri.a, tri.b, tri.c);
+
+            for (int y = minY; y <= maxY; ++y) {
+                for (int x = minX; x <= maxX; ++x) {
+                    Vector2f p(x + 0.5f, y + 0.5f);
+
+                    float w0 = edgeFunction(tri.b, tri.c, p);
+                    float w1 = edgeFunction(tri.c, tri.a, p);
+                    float w2 = edgeFunction(tri.a, tri.b, p);
+                    bool inside = false;
+                    
+                    if (area >= 0) {
+                        if (w0 >= 0 && w1 >= 0 && w2 >= 0) inside = true;
+                    } else {
+                        if (w0 <= 0 && w1 <= 0 && w2 <= 0) inside = true;
+                    }
+
+                    if (inside) {
+                        int index = (y * width + x) * 3;
+                        if (index >= 0 && index < buffer.size() - 2) {
+                            buffer[index    ] = tri.color.x();
+                            buffer[index + 1] = tri.color.y();
+                            buffer[index + 2] = tri.color.z();
+                        }
+                    }
+                }
+            }
+        }
+
+        outFrame.setData(buffer, colorformat);
+        return outFrame;
+    }
 };
 
-#endif
+
+class Scene {
+private:
+    std::vector<std::shared_ptr<Mesh>> _meshes;
+
+public:
+    Scene() {}
+
+    void addMesh(std::shared_ptr<Mesh> mesh) {
+        _meshes.push_back(mesh);
+    }
+
+    void clear() {
+        _meshes.clear();
+    }
+
+    frame render(Camera cam, int height, int width, float near, float far, frame::colormap colorformat = frame::colormap::RGB) {
+        std::vector<Triangle2D> allTriangles;
+
+        for (auto& mesh : _meshes) {
+            std::vector<Triangle2D> meshTris = mesh->project_2d(cam, height, width, near, far);
+            allTriangles.insert(allTriangles.end(), meshTris.begin(), meshTris.end());
+        }
+
+        std::sort(allTriangles.begin(), allTriangles.end(), [](const Triangle2D& a, const Triangle2D& b) {
+            return a.depth > b.depth;
+        });
+
+        return Mesh::rasterizeTriangles(allTriangles, width, height, colorformat);
+    }
+};
+
+#endif