From 4f0dba5eb41686a8b63ae858f7b47d7c03641dcd Mon Sep 17 00:00:00 2001 From: Yggdrasil75 Date: Mon, 2 Mar 2026 13:53:13 -0500 Subject: [PATCH] pushing this, lots of fun things to try. --- makefile | 2 +- tests/materialtest.cpp | 195 +++++++++++ tests/planet.cpp | 10 +- util/grid/grid3eigen.hpp | 679 ++++++++++++++++++++++++++++-------- util/sim/elementcontent.hpp | 383 ++++++++++++++++++++ util/sim/planet.hpp | 54 +-- 6 files changed, 1141 insertions(+), 182 deletions(-) create mode 100644 tests/materialtest.cpp create mode 100644 util/sim/elementcontent.hpp diff --git a/makefile b/makefile index 55510ef..d18a87e 100644 --- a/makefile +++ b/makefile @@ -7,7 +7,7 @@ STB_DIR := ./stb # Compiler and flags CXX := g++ -BASE_CXXFLAGS = -std=c++23 -O3 -fopenmp -march=native -I$(IMGUI_DIR) -I$(IMGUI_DIR)/backends -I$(STB_DIR) +BASE_CXXFLAGS = -std=c++23 -O3 -fopenmp -march=native -I$(IMGUI_DIR) -I$(IMGUI_DIR)/backends -I$(STB_DIR) -g BASE_CXXFLAGS += `pkg-config --cflags glfw3` CFLAGS = $(BASE_CXXFLAGS) LDFLAGS := -L./imgui -limgui -lGL diff --git a/tests/materialtest.cpp b/tests/materialtest.cpp new file mode 100644 index 0000000..1166517 --- /dev/null +++ b/tests/materialtest.cpp @@ -0,0 +1,195 @@ +#include +#include +#include +#include +#include + +// Include Eigen and project headers +#include "../eigen/Eigen/Dense" +#include "../util/grid/camera.hpp" +#include "../util/grid/grid3eigen.hpp" +#include "../util/output/frame.hpp" +#include "../util/output/bmpwriter.hpp" + +// Helper function to create a solid volume of voxels with material properties +void createBox(Octree& octree, const Eigen::Vector3f& center, const Eigen::Vector3f& size, + const Eigen::Vector3f& albedo, float emission = 0.0f, + float roughness = 0.8f, float metallic = 0.0f, float transmission = 0.0f, float ior = 1.45f) { + float step = 0.1f; // Voxel spacing + Eigen::Vector3f halfSize = size / 2.0f; + Eigen::Vector3f minB = center - halfSize; + Eigen::Vector3f maxB = center + halfSize; + + for (float x = minB.x(); x <= maxB.x(); x += step) { + for (float y = minB.y(); y <= maxB.y(); y += step) { + for (float z = minB.z(); z <= maxB.z(); z += step) { + Eigen::Vector3f pos(x, y, z); + + // .set(data, pos, visible, albedo, size, active, objectId, subId, emission, roughness, metallic, transmission, ior) + octree.set(1, pos, true, albedo, step, true, -1, 0, emission, roughness, metallic, transmission, ior); + } + } + } +} + +// Helper function to create a checkerboard pattern volume +void createCheckerBox(Octree& octree, const Eigen::Vector3f& center, const Eigen::Vector3f& size, + const Eigen::Vector3f& color1, const Eigen::Vector3f& color2, float checkerSize) { + float step = 0.1f; + Eigen::Vector3f halfSize = size / 2.0f; + Eigen::Vector3f minB = center - halfSize; + Eigen::Vector3f maxB = center + halfSize; + + for (float x = minB.x(); x <= maxB.x(); x += step) { + for (float y = minB.y(); y <= maxB.y(); y += step) { + for (float z = minB.z(); z <= maxB.z(); z += step) { + Eigen::Vector3f pos(x, y, z); + + // Use floor to correctly handle negative coordinates for the repeating pattern + int cx = static_cast(std::floor(x / checkerSize)); + int cy = static_cast(std::floor(y / checkerSize)); + int cz = static_cast(std::floor(z / checkerSize)); + + // 3D Checkerboard logic + bool isEven = ((cx + cy + cz) % 2 == 0); + Eigen::Vector3f albedo = isEven ? color1 : color2; + + octree.set(1, pos, true, albedo, step, true, -1, 0, 0.0f, 0.8f, 0.1f, 0.0f, 1.0f); + } + } + } +} + +int main() { + std::cout << "Initializing Octree..." << std::endl; + + // 1. Initialize Octree bounds + Eigen::Vector3f minBound(-10.0f, -10.0f, -10.0f); + Eigen::Vector3f maxBound(10.0f, 10.0f, 10.0f); + Octree octree(minBound, maxBound, 8, 16); + + // Set a dark background to emphasize the PBR light emission + octree.setBackgroundColor(Eigen::Vector3f(0.02f, 0.02f, 0.02f)); + octree.setSkylight(Eigen::Vector3f(0.01f, 0.01f, 0.01f)); + + std::cout << "Building scene..." << std::endl; + + // 2a. Build Room (Floor and 4 Walls) + Eigen::Vector3f cLightGray(0.8f, 0.8f, 0.8f); + Eigen::Vector3f cDarkGray(0.2f, 0.2f, 0.2f); + float chkSize = 1.0f; + + // Floor (Bounds: Z from -0.7 to -0.5) + // The boxes sit exactly on Z = -0.5 + createCheckerBox(octree, Eigen::Vector3f(0.0f, 0.0f, -0.6f), Eigen::Vector3f(14.4f, 14.4f, 0.2f), cLightGray, cDarkGray, chkSize); + + // Walls (Bounds: X/Y inner boundaries at +/- 7.0, rising from Z=-0.5 up to Z=7.5) + createCheckerBox(octree, Eigen::Vector3f( 7.1f, 0.0f, 3.5f), Eigen::Vector3f(0.2f, 14.4f, 8.0f), cLightGray, cDarkGray, chkSize); // +X + createCheckerBox(octree, Eigen::Vector3f(-7.1f, 0.0f, 3.5f), Eigen::Vector3f(0.2f, 14.4f, 8.0f), cLightGray, cDarkGray, chkSize); // -X + createCheckerBox(octree, Eigen::Vector3f( 0.0f, 7.1f, 3.5f), Eigen::Vector3f(14.0f, 0.2f, 8.0f), cLightGray, cDarkGray, chkSize); // +Y + createCheckerBox(octree, Eigen::Vector3f( 0.0f, -7.1f, 3.5f), Eigen::Vector3f(14.0f, 0.2f, 8.0f), cLightGray, cDarkGray, chkSize); // -Y + + // 2b. Create the 3x3 material sampler grid inside the room + Eigen::Vector3f cRed(1.0f, 0.1f, 0.1f); + Eigen::Vector3f cBlue(0.1f, 0.1f, 1.0f); + Eigen::Vector3f cPurple(0.6f, 0.1f, 0.8f); + Eigen::Vector3f size(1.0f, 1.0f, 1.0f); + + float sp = 2.0f; // spacing between cubes + + // --- LAYER 1: Metals --- + // (metallic = 1.0, slight roughness for blurry reflections, transmission = 0.0) + createBox(octree, Eigen::Vector3f(-sp, -sp, 0.0f), size, cRed, 0.0f, 0.15f, 1.0f, 0.0f, 1.45f); + createBox(octree, Eigen::Vector3f( 0, -sp, 0.0f), size, cBlue, 0.0f, 0.15f, 1.0f, 0.0f, 1.45f); + createBox(octree, Eigen::Vector3f( sp, -sp, 0.0f), size, cPurple, 0.0f, 0.15f, 1.0f, 0.0f, 1.45f); + + // --- LAYER 2: Opaque & Highly Refractive --- + // (metallic = 0.0, very low roughness. transmission = 0.0 for opacity, ior = 2.4 for extreme diamond-like reflection) + createBox(octree, Eigen::Vector3f(-sp, 0, 0.0f), size, cRed, 0.0f, 0.05f, 0.0f, 0.0f, 2.4f); + createBox(octree, Eigen::Vector3f( 0, 0, 0.0f), size, cBlue, 0.0f, 0.05f, 0.0f, 0.0f, 2.4f); + createBox(octree, Eigen::Vector3f( sp, 0, 0.0f), size, cPurple, 0.0f, 0.05f, 0.0f, 0.0f, 2.4f); + + // --- LAYER 3: Clear Glass --- + // (metallic = 0.0, near-zero roughness, transmission = 1.0 for full transparency, ior = 1.5 for glass) + createBox(octree, Eigen::Vector3f(-sp, sp, 0.0f), size, cRed, 0.0f, 0.01f, 0.0f, 1.0f, 1.5f); + createBox(octree, Eigen::Vector3f( 0, sp, 0.0f), size, cBlue, 0.0f, 0.01f, 0.0f, 1.0f, 1.5f); + createBox(octree, Eigen::Vector3f( sp, sp, 0.0f), size, cPurple, 0.0f, 0.01f, 0.0f, 1.0f, 1.5f); + + // White Light Box (Above) + // Placed near the ceiling (Z=7.4), made large (8x8) to cast soft shadows evenly over the whole 3x3 grid + createBox(octree, Eigen::Vector3f(0.0f, 0.0f, 7.4f), Eigen::Vector3f(8.0f, 8.0f, 0.2f), Eigen::Vector3f(1.0f, 1.0f, 1.0f), 15.0f); + + std::cout << "Optimizing and Generating LODs..." << std::endl; + octree.generateLODs(); + octree.printStats(); + + // 3. Setup rendering loop + int width = 512; + int height = 512; + int samples = 400; + int bounces = 5; + + struct View { + std::string name; + Eigen::Vector3f origin; + Eigen::Vector3f up; + }; + + // The walls are set perfectly at +/- 7.0 inner edges. + // Placing camera at +/- 6.8 will put it "just barely inside". + // Floor is at Z = -0.5, Wall top is at Z = 7.5 + std::vector views = { + {"+X", Eigen::Vector3f( 6.8f, 0.0f, 1.0f), Eigen::Vector3f(0.0f, 0.0f, 1.0f)}, + {"-X", Eigen::Vector3f(-6.8f, 0.0f, 1.0f), Eigen::Vector3f(0.0f, 0.0f, 1.0f)}, + {"+Y", Eigen::Vector3f( 0.0f, 6.8f, 1.0f), Eigen::Vector3f(0.0f, 0.0f, 1.0f)}, + {"-Y", Eigen::Vector3f( 0.0f, -6.8f, 1.0f), Eigen::Vector3f(0.0f, 0.0f, 1.0f)}, + {"+Z", Eigen::Vector3f( 0.0f, 0.0f, 7.3f), Eigen::Vector3f(0.0f, 1.0f, 0.0f)} // Looking down from just beneath wall top + }; + + Eigen::Vector3f target(0.0f, 0.0f, 0.5f); + + for (const auto& view : views) { + std::cout << "\nRendering view from " << view.name << " direction (Fast Pass)..." << std::endl; + + Camera cam; + cam.origin = view.origin; + cam.direction = (target - view.origin).normalized(); + cam.up = view.up; + + frame out = octree.fastRenderFrame(cam, height, width, frame::colormap::RGB); + + std::string filename = "output/fast/render_" + view.name + ".bmp"; + BMPWriter::saveBMP(filename, out); + } + + for (const auto& view : views) { + std::cout << "\nRendering view from " << view.name << " direction (Medium 60s Pass)..." << std::endl; + + Camera cam; + cam.origin = view.origin; + cam.direction = (target - view.origin).normalized(); + cam.up = view.up; + + frame out = octree.renderFrameTimed(cam, height, width, frame::colormap::RGB, 60, bounces, false, true); + + std::string filename = "output/medium/render_" + view.name + ".bmp"; + BMPWriter::saveBMP(filename, out); + } + + for (const auto& view : views) { + std::cout << "\nRendering view from " << view.name << " direction (Slow 400 Samples Pass)..." << std::endl; + + Camera cam; + cam.origin = view.origin; + cam.direction = (target - view.origin).normalized(); + cam.up = view.up; + + frame out = octree.renderFrame(cam, height, width, frame::colormap::RGB, samples, bounces, false, true); + + std::string filename = "output/slow/render_" + view.name + ".bmp"; + BMPWriter::saveBMP(filename, out); + } + + std::cout << "\nAll renders complete!" << std::endl; + return 0; +} \ No newline at end of file diff --git a/tests/planet.cpp b/tests/planet.cpp index c16f46d..ca33ac1 100644 --- a/tests/planet.cpp +++ b/tests/planet.cpp @@ -331,7 +331,10 @@ public: break; } - sim.grid.update(p.currentPos, p.currentPos, p, true, color, sim.config.voxelSize, true, -2, false, 0.0f, 0.0f, 0.0f); + if (!sim.grid.setColor(p.currentPos, color)) { + std::cout << "node failed to set" << std::endl; + } + // sim.grid.update(p.currentPos, p.currentPos, p, true, color, sim.config.voxelSize, true, -2, false, 0.0f, 0.0f, 0.0f); } for (auto& p : sim.config.interpolatedNodes) { v3 color = p.originColor; @@ -356,7 +359,10 @@ public: break; } - sim.grid.update(p.currentPos, p.currentPos, p, true, color, sim.config.voxelSize, true, -2, false, 0.0f, 0.0f, 0.0f); + if (!sim.grid.setColor(p.currentPos, color)) { + std::cout << "node failed to set" << std::endl; + } + // sim.grid.update(p.currentPos, p.currentPos, p, true, color, sim.config.voxelSize, true, -2, false, 0.0f, 0.0f, 0.0f); } } diff --git a/util/grid/grid3eigen.hpp b/util/grid/grid3eigen.hpp index 1c9373c..10cec9c 100644 --- a/util/grid/grid3eigen.hpp +++ b/util/grid/grid3eigen.hpp @@ -19,6 +19,8 @@ #include #include #include +#include +#include #ifdef SSE #include @@ -31,6 +33,12 @@ class Octree { public: using PointType = Eigen::Matrix; using BoundingBox = std::pair; + + struct Vector3fCompare { + bool operator()(const Eigen::Vector3f& a, const Eigen::Vector3f& b) const { + return std::tie(a.x(), a.y(), a.z()) < std::tie(b.x(), b.y(), b.z()); + } + }; struct NodeData { T data; @@ -41,20 +49,21 @@ public: bool visible; float size; Eigen::Vector3f color; - bool light; float emittance; - float refraction; - float reflection; + float roughness; + float metallic; + float transmission; + float ior; NodeData(const T& data, const PointType& pos, bool visible, Eigen::Vector3f color, float size = 0.01f, - bool active = true, int objectId = -1, int subId = 0, bool light = false, float emittance = 0.0f, - float refraction = 0.0f, float reflection = 0.0f) + bool active = true, int objectId = -1, int subId = 0, float emittance = 0.0f, + float roughness = 1.0f, float metallic = 0.0f, float transmission = 0.0f, float ior = 1.45f) : data(data), position(pos), objectId(objectId), subId(subId), active(active), visible(visible), - color(color), size(size), light(light), emittance(emittance), refraction(refraction), - reflection(reflection) {} + color(color), size(size), emittance(emittance), roughness(roughness), metallic(metallic), + transmission(transmission), ior(ior) {} - NodeData() : objectId(-1), subId(0), active(false), visible(false), size(0.0f), light(false), - emittance(0.0f), refraction(0.0f), reflection(0.0f) {} + NodeData() : objectId(-1), subId(0), active(false), visible(false), size(0.0f), emittance(0.0f), + roughness(1.0f), metallic(0.0f), transmission(0.0f), ior(1.45f) {} // Helper method to get half-size for cube PointType getHalfSize() const { @@ -110,8 +119,6 @@ private: std::map, std::shared_ptr> meshCache_; std::set> dirtyMeshes_; int nextSubIdGenerator_ = 1; - constexpr static float ior = 1.45f; - constexpr static float η = 1.0f / ior; void invalidateMesh(int objectId, int subId) { if (objectId < 0) return; @@ -232,21 +239,38 @@ private: } return true; } else { - // Keep in parent if size is larger than or equal to the node size - if (pointData->size >= node->nodeSize) { - node->points.emplace_back(pointData); - return true; - } - bool inserted = false; for (int i = 0; i < 8; ++i) { if (node->children[i] && boxIntersectsBox(node->children[i]->bounds, cubeBounds)) { + size++; inserted |= insertRecursive(node->children[i].get(), pointData, depth + 1); } } return inserted; } } + + bool invalidateNodeLODRecursive(OctreeNode* node, const BoundingBox& bounds) { + if (!boxIntersectsBox(node->bounds, bounds)) return false; + { + std::lock_guard lock(node->lodMutex); + node->lodData = nullptr; + } + if (!node->isLeaf) { + for (int i = 0; i < 8; ++i) { + if (node->children[i]) { + invalidateNodeLODRecursive(node->children[i].get(), bounds); + } + } + } + return true; + } + + void invalidateLODForPoint(const std::shared_ptr& pointData) { + if (root_ && pointData) { + invalidateNodeLODRecursive(root_.get(), pointData->getCubeBounds()); + } + } void ensureLOD(const OctreeNode* node) const { std::lock_guard lock(node->lodMutex); @@ -255,9 +279,10 @@ private: PointType avgPos = PointType::Zero(); Eigen::Vector3f avgColor = Eigen::Vector3f::Zero(); float avgEmittance = 0.0f; - float avgRefl = 0.0f; - float avgRefr = 0.0f; - bool anyLight = false; + float avgRoughness = 0.0f; + float avgMetallic = 0.0f; + float avgTransmission = 0.0f; + float avgIor = 0.0f; int count = 0; if (node->isLeaf && node->points.size() == 1) { @@ -271,9 +296,10 @@ private: if (!item || !item->active || !item->visible) return; avgColor += item->color; avgEmittance += item->emittance; - avgRefl += item->reflection; - avgRefr += item->refraction; - if (item->light) anyLight = true; + avgRoughness += item->roughness; + avgMetallic += item->metallic; + avgTransmission += item->transmission; + avgIor += item->ior; count++; }; @@ -300,9 +326,10 @@ private: lod->size = nodeDims.maxCoeff(); lod->emittance = avgEmittance * invCount; - lod->reflection = avgRefl * invCount; - lod->refraction = avgRefr * invCount; - lod->light = anyLight; + lod->roughness = avgRoughness * invCount; + lod->metallic = avgMetallic * invCount; + lod->transmission = avgTransmission * invCount; + lod->ior = avgIor * invCount; lod->active = true; lod->visible = true; lod->objectId = -1; @@ -333,7 +360,7 @@ private: return nullptr; } - bool removeRecursive(OctreeNode* node, const BoundingBox& bounds, const PointType& pos, float tolerance, bool& sizeDecremented) { + bool removeRecursive(OctreeNode* node, const BoundingBox& bounds, const PointType& pos, float tolerance) { { std::lock_guard lock(node->lodMutex); node->lodData = nullptr; @@ -352,17 +379,14 @@ private: if (it != node->points.end()) { node->points.erase(it, node->points.end()); - if (!sizeDecremented) { - size--; - sizeDecremented = true; - } + size--; foundAny = true; } if (!node->isLeaf) { for (int i = 0; i < 8; ++i) { if (node->children[i]) { - foundAny |= removeRecursive(node->children[i].get(), bounds, pos, tolerance, sizeDecremented); + foundAny |= removeRecursive(node->children[i].get(), bounds, pos, tolerance); } } } @@ -465,8 +489,47 @@ private: } } + PointType sampleGGX(const PointType& n, float roughness, uint32_t& state) const { + float alpha = std::max(EPSILON, roughness * roughness); + float r1 = float(rand_r(&state)) / float(RAND_MAX); + float r2 = float(rand_r(&state)) / float(RAND_MAX); + + float phi = 2.0f * M_PI * r1; + float denom = 1.0f + (alpha * alpha - 1.0f) * r2; + denom = std::max(denom, EPSILON); + float cosTheta = std::sqrt(std::max(0.0f, (1.0f - r2) / denom)); + float sinTheta = std::sqrt(std::max(0.0f, 1.0f - cosTheta * cosTheta)); + + PointType h; + h[0] = sinTheta * std::cos(phi); + h[1] = sinTheta * std::sin(phi); + h[2] = cosTheta; + + PointType up = std::abs(n.z()) < 0.999f ? PointType(0,0,1) : PointType(1,0,0); + PointType tangent = up.cross(n).normalized(); + PointType bitangent = n.cross(tangent); + + return (tangent * h[0] + bitangent * h[1] + n * h[2]).normalized(); + } + + PointType sampleCosineHemisphere(const PointType& n, uint32_t& state) const { + float r1 = float(rand_r(&state)) / float(RAND_MAX); + float r2 = float(rand_r(&state)) / float(RAND_MAX); + float phi = 2.0f * M_PI * r1; + float r = std::sqrt(r2); + float x = r * std::cos(phi); + float y = r * std::sin(phi); + float z = std::sqrt(std::max(0.0f, 1.0f - x*x - y*y)); + + PointType up = std::abs(n.z()) < 0.999f ? PointType(0,0,1) : PointType(1,0,0); + PointType tangent = up.cross(n).normalized(); + PointType bitangent = n.cross(tangent); + + return (tangent * x + bitangent * y + n * z).normalized(); + } + Eigen::Vector3f traceRay(const PointType& rayOrig, const PointType& rayDir, int bounces, uint32_t& rngState, - int maxBounces, bool globalIllumination, bool useLod) { + int maxBounces, bool globalIllumination, bool useLod) { if (bounces > maxBounces) return globalIllumination ? skylight_ : Eigen::Vector3f::Zero(); auto hit = voxelTraverse(rayOrig, rayDir, std::numeric_limits::max(), useLod); @@ -483,58 +546,110 @@ private: Ray ray(rayOrig, rayDir); rayCubeIntersect(ray, obj.get(), t, normal, hitPoint); - if (obj->light) { - return obj->color * obj->emittance; - } - Eigen::Vector3f finalColor = globalIllumination ? skylight_ : Eigen::Vector3f::Zero(); - if (obj->light) return finalColor; - - float refl = obj->reflection; - float refr = obj->refraction; - float diffuseProb = 1.0f - refl - refr; - float sum = refr + refl + diffuseProb; - if (sum <= 0.001f) return finalColor; - - float roll = float(rand_r(&rngState)) / float(RAND_MAX); - Eigen::Vector3f npc = Eigen::Vector3f::Zero(); - PointType nextDir; - PointType bias = normal * 0.002f; - - if (roll < diffuseProb) { - nextDir = randomInHemisphere(normal, rngState); - Eigen::Vector3f incoming = traceRay(hitPoint + bias, nextDir, bounces + 1, rngState, maxBounces, globalIllumination, useLod); - npc = obj->color.cwiseProduct(incoming); - } else if (roll < diffuseProb + refl) { - nextDir = (rayDir - 2.0f * normal.dot(rayDir) * normal).normalized(); - Eigen::Vector3f incoming = traceRay(hitPoint + bias, nextDir, bounces + 1, rngState, maxBounces, globalIllumination, useLod); - npc = obj->color.cwiseProduct(incoming); - npc /= refl; - } else { - float lη = η; - PointType n_eff = normal; - float cosI = -normal.dot(rayDir); - if (cosI < 0) { - cosI = -cosI; - n_eff = -normal; - lη = ior; - } - float k = 1.0f - lη * lη * (1.0f - cosI * cosI); - - if (k < 0.0f) { - nextDir = (rayDir - 2.0f * rayDir.dot(n_eff) * n_eff).normalized(); - } else { - nextDir = (lη * rayDir + (lη * cosI - std::sqrt(k)) * n_eff).normalized(); - } - Eigen::Vector3f incoming = traceRay(hitPoint - (n_eff * 0.002f), nextDir, bounces + 1, rngState, maxBounces, globalIllumination, useLod); - npc = obj->color.cwiseProduct(incoming); - - npc /= refr; + if (obj->emittance > 0.0f) { + finalColor += obj->color * obj->emittance; } - return finalColor + npc; - } + float roughness = std::clamp(obj->roughness, 0.01f, 1.0f); + float metallic = std::clamp(obj->metallic, 0.0f, 1.0f); + float transmission = std::clamp(obj->transmission, 0.0f, 1.0f); + + PointType V = -rayDir; + float cosThetaI = normal.dot(V); + bool isInside = cosThetaI < 0.0f; + PointType n_eff = isInside ? -normal : normal; + cosThetaI = std::max(0.001f, n_eff.dot(V)); + float coordMax = hitPoint.cwiseAbs().maxCoeff(); + float rayOffset = std::max(1e-4f, 1e-5f * coordMax); + + Eigen::Vector3f F0 = Eigen::Vector3f::Constant(0.04f); + F0 = F0 * (1.0f - metallic) + obj->color * metallic; + + PointType H = sampleGGX(n_eff, roughness, rngState); + float VdotH = std::max(0.001f, V.dot(H)); + + Eigen::Vector3f F_spec = F0 + (Eigen::Vector3f::Constant(1.0f) - F0) * std::pow(std::max(0.0f, 1.0f - VdotH), 5.0f); + + PointType specDir = (2.0f * VdotH * H - V).normalized(); + Eigen::Vector3f W_spec = Eigen::Vector3f::Zero(); + + if (specDir.dot(n_eff) > 0.0f) { + float NdotV = cosThetaI; + float NdotL = std::max(0.001f, n_eff.dot(specDir)); + float NdotH = std::max(0.001f, n_eff.dot(H)); + + float k_smith = (roughness * roughness) / 2.0f; + float G = (NdotV / (NdotV * (1.0f - k_smith) + k_smith)) * (NdotL / (NdotL * (1.0f - k_smith) + k_smith)); + + W_spec = F_spec * G * VdotH / (NdotV * NdotH); + } + + Eigen::Vector3f W_second = Eigen::Vector3f::Zero(); + PointType secondDir; + PointType secondOrigin; + + float transmissionWeight = transmission * (1.0f - metallic); + float diffuseWeight = (1.0f - transmission) * (1.0f - metallic); + + if (transmissionWeight > 0.0f) { + float eta = isInside ? obj->ior : (1.0f / obj->ior); + float k = 1.0f - eta * eta * (1.0f - VdotH * VdotH); + + if (k >= 0.0f) { + secondDir = ((eta * VdotH - std::sqrt(k)) * H - eta * V).normalized(); + secondOrigin = hitPoint - n_eff * rayOffset; + W_second = (Eigen::Vector3f::Constant(1.0f) - F_spec) * transmissionWeight; + W_second = W_second.cwiseProduct(obj->color); + } else { + Eigen::Vector3f tirWeight = (Eigen::Vector3f::Constant(1.0f) - F_spec) * transmissionWeight; + W_spec += tirWeight.cwiseProduct(obj->color); + } + } else if (diffuseWeight > 0.0f) { + secondDir = sampleCosineHemisphere(n_eff, rngState); + secondOrigin = hitPoint + n_eff * rayOffset; + W_second = (Eigen::Vector3f::Constant(1.0f) - F_spec) * diffuseWeight; + W_second = W_second.cwiseProduct(obj->color); + } + + W_spec = W_spec.cwiseMin(Eigen::Vector3f::Constant(4.0f)); + W_second = W_second.cwiseMin(Eigen::Vector3f::Constant(4.0f)); + + float lumSpec = W_spec.maxCoeff(); + float lumSecond = W_second.maxCoeff(); + + bool doSplit = (bounces <= 1); + + if (doSplit) { + Eigen::Vector3f specColor = Eigen::Vector3f::Zero(); + if (lumSpec > 0.001f) { + specColor = W_spec.cwiseProduct(traceRay(hitPoint + n_eff * rayOffset, specDir, bounces + 1, rngState, maxBounces, globalIllumination, useLod)); + } + + Eigen::Vector3f secondColor = Eigen::Vector3f::Zero(); + if (lumSecond > 0.001f) { + secondColor = W_second.cwiseProduct(traceRay(secondOrigin, secondDir, bounces + 1, rngState, maxBounces, globalIllumination, useLod)); + } + + return finalColor + specColor + secondColor; + } else { + float totalLum = lumSpec + lumSecond; + if (totalLum < 0.0001f) return finalColor; + + float pSpec = lumSpec / totalLum; + float roll = float(rand_r(&rngState)) / float(RAND_MAX); + + if (roll < pSpec) { + Eigen::Vector3f sample = traceRay(hitPoint + n_eff * rayOffset, specDir, bounces + 1, rngState, maxBounces, globalIllumination, useLod); + return finalColor + (W_spec / std::max(EPSILON, pSpec)).cwiseProduct(sample); + } else { + Eigen::Vector3f sample = traceRay(secondOrigin, secondDir, bounces + 1, rngState, maxBounces, globalIllumination, useLod); + return finalColor + (W_second / std::max(EPSILON, 1.0f - pSpec)).cwiseProduct(sample); + } + } + } + void clearNode(OctreeNode* node) { if (!node) return; @@ -575,6 +690,53 @@ private: } } + void optimizeRecursive(OctreeNode* node) { + if (!node) return; + + if (node->isLeaf) { + return; + } + + bool childrenAreLeaves = true; + for (int i = 0; i < 8; ++i) { + if (node->children[i]) { + optimizeRecursive(node->children[i].get()); + if (!node->children[i]->isLeaf) { + childrenAreLeaves = false; + } + } + } + + if (childrenAreLeaves) { + std::vector> allPoints = node->points; + for (int i = 0; i < 8; ++i) { + if (node->children[i]) { + allPoints.insert(allPoints.end(), node->children[i]->points.begin(), node->children[i]->points.end()); + } + } + + std::sort(allPoints.begin(), allPoints.end(), [](const std::shared_ptr& a, const std::shared_ptr& b) { + return a.get() < b.get(); + }); + allPoints.erase(std::unique(allPoints.begin(), allPoints.end(), [](const std::shared_ptr& a, const std::shared_ptr& b) { + return a.get() == b.get(); + }), allPoints.end()); + + if (allPoints.size() <= maxPointsPerNode) { + node->points = std::move(allPoints); + for (int i = 0; i < 8; ++i) { + node->children[i].reset(nullptr); + } + node->isLeaf = true; + + { + std::lock_guard lock(node->lodMutex); + node->lodData = nullptr; + } + } + } + } + template void writeVal(std::ofstream& out, const V& val) const { out.write(reinterpret_cast(&val), sizeof(V)); @@ -613,10 +775,11 @@ private: writeVal(out, pt->visible); writeVal(out, pt->size); writeVec3(out, pt->color); - writeVal(out, pt->light); writeVal(out, pt->emittance); - writeVal(out, pt->refraction); - writeVal(out, pt->reflection); + writeVal(out, pt->roughness); + writeVal(out, pt->metallic); + writeVal(out, pt->transmission); + writeVal(out, pt->ior); } } else { // Write bitmask of active children @@ -656,10 +819,11 @@ private: readVal(in, pt->visible); readVal(in, pt->size); readVec3(in, pt->color); - readVal(in, pt->light); readVal(in, pt->emittance); - readVal(in, pt->refraction); - readVal(in, pt->reflection); + readVal(in, pt->roughness); + readVal(in, pt->metallic); + readVal(in, pt->transmission); + readVal(in, pt->ior); node->points.push_back(pt); } @@ -751,13 +915,29 @@ private: bool rayCubeIntersect(const Ray& ray, const NodeData* cube, float& t, PointType& normal, PointType& hitPoint) const { BoundingBox bounds = cube->getCubeBounds(); - float tMin, tMax; - if (!rayBoxIntersect(ray, bounds, tMin, tMax)) { + float t0x = (bounds.first[0] - ray.origin[0]) * ray.invDir[0]; + float t1x = (bounds.second[0] - ray.origin[0]) * ray.invDir[0]; + if (ray.invDir[0] < 0.0f) std::swap(t0x, t1x); + + float t0y = (bounds.first[1] - ray.origin[1]) * ray.invDir[1]; + float t1y = (bounds.second[1] - ray.origin[1]) * ray.invDir[1]; + if (ray.invDir[1] < 0.0f) std::swap(t0y, t1y); + + float tMin = std::max(t0x, t0y); + float tMax = std::min(t1x, t1y); + + float t0z = (bounds.first[2] - ray.origin[2]) * ray.invDir[2]; + float t1z = (bounds.second[2] - ray.origin[2]) * ray.invDir[2]; + if (ray.invDir[2] < 0.0f) std::swap(t0z, t1z); + + tMin = std::max(tMin, t0z); + tMax = std::min(tMax, t1z); + + if (tMax < std::max(0.0f, tMin) || tMax < 0.0f) { return false; } - - if (tMin < EPSILON) { - if (tMax < EPSILON) return false; + + if (tMin < 0.0f) { t = tMax; } else { t = tMin; @@ -765,14 +945,28 @@ private: hitPoint = ray.origin + ray.dir * t; - normal = PointType::Zero(); PointType dMin = (hitPoint - bounds.first).cwiseAbs(); PointType dMax = (hitPoint - bounds.second).cwiseAbs(); - // float thresh = EPSILON * 1.00001f; - PointType nMin = (dMin.array() < EPSILON).cast(); - PointType nMax = (dMax.array() < EPSILON).cast(); - normal = -nMin + nMax; + float minDist = std::numeric_limits::max(); + int minAxis = 0; + float sign = 1.0f; + + for (int i = 0; i < Dim; ++i) { + if (dMin[i] < minDist) { + minDist = dMin[i]; + minAxis = i; + sign = -1.0f; + } + if (dMax[i] < minDist) { + minDist = dMax[i]; + minAxis = i; + sign = 1.0f; + } + } + + normal = PointType::Zero(); + normal[minAxis] = sign; return true; } @@ -822,12 +1016,6 @@ private: std::array color; }; - struct Vector3fCompare { - bool operator()(const Eigen::Vector3f& a, const Eigen::Vector3f& b) const { - return std::tie(a.x(), a.y(), a.z()) < std::tie(b.x(), b.y(), b.z()); - } - }; - std::pair vertexInterpolate(float isolevel, const PointType& p1, const PointType& p2, float val1, float val2, const PointType& c1, const PointType& c2) { if (std::abs(isolevel - val1) < EPSILON) return {p1, c1}; if (std::abs(isolevel - val2) < EPSILON) return {p2, c2}; @@ -986,13 +1174,12 @@ public: ensureLOD(root_.get()); } - bool set(const T& data, const PointType& pos, bool visible, Eigen::Vector3f color, float size, bool active, - int objectId = -1, int subId = 0, bool light = false, float emittance = 0.0f, float refraction = 0.0f, - float reflection = 0.0f) { + bool set(const T& data, const PointType& pos, bool visible, Eigen::Vector3f color, float size = 0.01f, bool active = true, + int objectId = -1, int subId = 0, float emittance = 0.0f, float roughness = 1.0f, + float metallic = 0.0f, float transmission = 0.0f, float ior = 1.45f) { auto pointData = std::make_shared(data, pos, visible, color, size, active, objectId, subId, - light, emittance, refraction, reflection); + emittance, roughness, metallic, transmission, ior); if (insertRecursive(root_.get(), pointData, 0)) { - this->size++; return true; } return false; @@ -1038,7 +1225,6 @@ public: readVal(in, maxPointsPerNode); readVal(in, size); - // Load global settings readVec3(in, skylight_); readVec3(in, backgroundColor_); @@ -1047,7 +1233,8 @@ public: readVec3(in, maxBound); root_ = std::make_unique(minBound, maxBound); - deserializeNode(in, root_.get()); + std::multimap, Vector3fCompare> pointMap; + deserializeNode(in, root_.get(), pointMap); in.close(); std::cout << "successfully loaded grid from " << filename << std::endl; @@ -1065,8 +1252,7 @@ public: bool remove(const PointType& pos, float tolerance = EPSILON) { auto pt = find(pos, tolerance); if (!pt) return false; - bool sizeDecremented = false; - return removeRecursive(root_.get(), pt->getCubeBounds(), pos, tolerance, sizeDecremented); + return removeRecursive(root_.get(), pt->getCubeBounds(), pos, tolerance); } std::vector> findInRadius(const PointType& center, float radius, int objectid = -1) const { @@ -1079,22 +1265,27 @@ public: return results; } + bool update(const PointType& pos, const T& newData) { + auto pointData = find(pos); + if (!pointData) return false; + else pointData->data = newData; + return true; + } + bool update(const PointType& oldPos, const PointType& newPos, const T& newData, bool newVisible = true, Eigen::Vector3f newColor = Eigen::Vector3f(1.0f, 1.0f, 1.0f), float newSize = 0.01f, bool newActive = true, - int newObjectId = -2, bool newLight = false, float newEmittance = 0.0f, float newRefraction = 0.0f, - float newReflection = 0.0f, float tolerance = EPSILON) { + int newObjectId = -2, float newEmittance = -1.0f, float newRoughness = -1.0f, + float newMetallic = -1.0f, float newTransmission = -1.0f, float newIor = -1.0f, float tolerance = EPSILON) { auto pointData = find(oldPos, tolerance); if (!pointData) return false; int oldSubId = pointData->subId; int targetObjId = (newObjectId != -2) ? newObjectId : pointData->objectId; - int finalSubId = oldSubId; if (oldSubId == 0 && targetObjId == pointData->objectId) { finalSubId = nextSubIdGenerator_++; - auto neighbors = findInRadius(oldPos, newSize * 3.0f, targetObjId); for(auto& n : neighbors) { if(n->subId == 0) { @@ -1104,18 +1295,26 @@ public: } } - if (!remove(oldPos, tolerance)) return false; + removeRecursive(root_.get(), pointData->getCubeBounds(), oldPos, tolerance); - bool res = set(newData, newPos, newVisible, - newColor != Eigen::Vector3f(1.0f, 1.0f, 1.0f) ? newColor : pointData->color, - newSize > 0 ? newSize : pointData->size, - newActive, targetObjId, finalSubId, newLight, - newEmittance > 0 ? newEmittance : pointData->emittance, - newRefraction >= 0 ? newRefraction : pointData->refraction, - newReflection >= 0 ? newReflection : pointData->reflection); + pointData->data = newData; + pointData->position = newPos; + pointData->visible = newVisible; + if (newColor != Eigen::Vector3f(1.0f, 1.0f, 1.0f)) pointData->color = newColor; + if (newSize > 0) pointData->size = newSize; + pointData->active = newActive; + pointData->objectId = targetObjId; + pointData->subId = finalSubId; + + if (newEmittance >= 0) pointData->emittance = newEmittance; + if (newRoughness >= 0) pointData->roughness = newRoughness; + if (newMetallic >= 0) pointData->metallic = newMetallic; + if (newTransmission >= 0) pointData->transmission = newTransmission; + if (newIor >= 0) pointData->ior = newIor; + + bool res = insertRecursive(root_.get(), pointData, 0); if(res) { - // Mark relevant meshes dirty invalidateMesh(targetObjId, finalSubId); if(finalSubId != oldSubId) invalidateMesh(targetObjId, oldSubId); } @@ -1126,11 +1325,13 @@ public: bool move(const PointType& pos, const PointType& newPos) { auto pointData = find(pos); if (!pointData) return false; - bool success = set(pointData->data, newPos, pointData->visible, pointData->color, pointData->size, - pointData->active, pointData->objectId, pointData->light, pointData->emittance, - pointData->refraction, pointData->reflection); - if (success) { - remove(pos); + + bool sizeDecremented = false; + removeRecursive(root_.get(), pointData->getCubeBounds(), pos, EPSILON); + + pointData->position = newPos; + + if (insertRecursive(root_.get(), pointData, 0)) { return true; } return false; @@ -1140,62 +1341,71 @@ public: auto pointData = find(pos, tolerance); if (!pointData) return false; pointData->objectId = objectId; + invalidateLODForPoint(pointData); return true; } bool updateData(const PointType& pos, const T& newData, float tolerance = EPSILON) { auto pointData = find(pos, tolerance); if (!pointData) return false; - pointData->data = newData; + invalidateLODForPoint(pointData); return true; } bool setActive(const PointType& pos, bool active, float tolerance = EPSILON) { auto pointData = find(pos, tolerance); if (!pointData) return false; - pointData->active = active; + invalidateLODForPoint(pointData); return true; } bool setVisible(const PointType& pos, bool visible, float tolerance = EPSILON) { auto pointData = find(pos, tolerance); if (!pointData) return false; - pointData->visible = visible; - return true; - } - - bool setLight(const PointType& pos, bool light, float tolerance = EPSILON) { - auto pointData = find(pos, tolerance); - if (!pointData) return false; - - pointData->light = light; + invalidateLODForPoint(pointData); return true; } bool setColor(const PointType& pos, Eigen::Vector3f color, float tolerance = EPSILON) { auto pointData = find(pos, tolerance); if (!pointData) return false; - pointData->color = color; - return true; - } - - bool setReflection(const PointType& pos, float reflection, float tolerance = EPSILON) { - auto pointData = find(pos, tolerance); - if (!pointData) return false; - - pointData->reflection = reflection; + invalidateLODForPoint(pointData); return true; } bool setEmittance(const PointType& pos, float emittance, float tolerance = EPSILON) { auto pointData = find(pos, tolerance); if (!pointData) return false; - pointData->emittance = emittance; + invalidateLODForPoint(pointData); + return true; + } + + bool setRoughness(const PointType& pos, float roughness, float tolerance = EPSILON) { + auto pointData = find(pos, tolerance); + if (!pointData) return false; + pointData->roughness = roughness; + invalidateLODForPoint(pointData); + return true; + } + + bool setMetallic(const PointType& pos, float metallic, float tolerance = EPSILON) { + auto pointData = find(pos, tolerance); + if (!pointData) return false; + pointData->metallic = metallic; + invalidateLODForPoint(pointData); + return true; + } + + bool setTransmission(const PointType& pos, float transmission, float tolerance = EPSILON) { + auto pointData = find(pos, tolerance); + if (!pointData) return false; + pointData->transmission = transmission; + invalidateLODForPoint(pointData); return true; } @@ -1325,7 +1535,7 @@ public: rayCubeIntersect(ray, obj.get(), t, normal, hitPoint); color = obj->color; - if (obj->light) { + if (obj->emittance > 0.0f) { color = color * obj->emittance; } else { float diffuse = std::max(0.0f, normal.dot(globalLightDir)); @@ -1351,6 +1561,163 @@ public: return outFrame; } + frame renderFrameTimed(const Camera& cam, int height, int width, frame::colormap colorformat = frame::colormap::RGB, + double maxTimeSeconds = 0.16, int maxBounces = 4, bool globalIllumination = false, bool useLod = true) { + auto startTime = std::chrono::high_resolution_clock::now(); + + generateLODs(); + PointType origin = cam.origin; + PointType dir = cam.direction.normalized(); + PointType up = cam.up.normalized(); + PointType right = cam.right(); + + frame outFrame(width, height, colorformat); + std::vector colorBuffer(width * height * 3, 0.0f); + std::vector accumColor(width * height, Eigen::Vector3f::Zero()); + std::vector sampleCount(width * height, 0); + + const float aspect = static_cast(width) / height; + const float fovRad = cam.fovRad(); + const float tanHalfFov = std::tan(fovRad * 0.5f); + const float tanfovy = tanHalfFov; + const float tanfovx = tanHalfFov * aspect; + + const PointType globalLightDir = (-cam.direction * 0.2f).normalized(); + const float fogStart = 1000.0f; + const float minVisibility = 0.2f; + + #pragma omp parallel for schedule(dynamic, 128) collapse(2) + for (int y = 0; y < height; ++y) { + for (int x = 0; x < width; ++x) { + int pidx = (y * width + x); + + float px = (2.0f * (x + 0.5f) / width - 1.0f) * tanfovx; + float py = (1.0f - 2.0f * (y + 0.5f) / height) * tanfovy; + + PointType rayDir = dir + (right * px) + (up * py); + rayDir.normalize(); + + Eigen::Vector3f color = backgroundColor_; + Ray ray(origin, rayDir); + auto hit = fastVoxelTraverse(ray, maxDistance_, true); + if (hit != nullptr) { + auto obj = hit; + + float t = 0.0f; + PointType normal, hitPoint; + + rayCubeIntersect(ray, obj.get(), t, normal, hitPoint); + color = obj->color; + + if (obj->emittance > 0.0f) { + color = color * obj->emittance; + } else { + float diffuse = std::max(0.0f, normal.dot(globalLightDir)); + float ambient = 0.35f; + float intensity = std::min(1.0f, ambient + diffuse * 0.65f); + color = color * intensity; + } + + float fogFactor = std::clamp((maxDistance_ - t) / (maxDistance_ - fogStart), minVisibility, 1.0f); + color = color * fogFactor + backgroundColor_ * (1.0f - fogFactor); + } + + accumColor[pidx] = color; + sampleCount[pidx] = 0; + } + } + + auto now = std::chrono::high_resolution_clock::now(); + std::chrono::duration elapsed = now - startTime; + + if (elapsed.count() < maxTimeSeconds) { + std::atomic timeUp(false); + std::atomic counter(0); + uint64_t totalPixels = static_cast(width) * height; + + uint64_t stride = 15485863; + auto gcd = [](uint64_t a, uint64_t b) { + while (b != 0) { + uint64_t temp = b; + b = a % b; + a = temp; + } + return a; + }; + while (gcd(stride, totalPixels) != 1) { + stride += 2; + } + + #pragma omp parallel + { + uint32_t localSeed = omp_get_thread_num() * 1973 + 9277; + int chunkSize = 64; + + while (!timeUp.load(std::memory_order_relaxed)) { + uint64_t startIdx = counter.fetch_add(chunkSize, std::memory_order_relaxed); + + if (omp_get_thread_num() == 0) { + auto checkTime = std::chrono::high_resolution_clock::now(); + std::chrono::duration checkElapsed = checkTime - startTime; + if (checkElapsed.count() >= maxTimeSeconds) { + timeUp.store(true, std::memory_order_relaxed); + break; + } + } + + if (timeUp.load(std::memory_order_relaxed)) break; + + for (int i = 0; i < chunkSize; ++i) { + uint64_t currentOffset = startIdx + i; + uint64_t pidx = (currentOffset * stride) % totalPixels; + + int y = pidx / width; + int x = pidx % width; + + float px = (2.0f * (x + 0.5f) / width - 1.0f) * tanfovx; + float py = (1.0f - 2.0f * (y + 0.5f) / height) * tanfovy; + + PointType rayDir = dir + (right * px) + (up * py); + rayDir.normalize(); + + uint32_t pass = currentOffset / totalPixels; + uint32_t seed = pidx * 1973 + pass * 12345 + localSeed; + + Eigen::Vector3f pbrColor = traceRay(origin, rayDir, 0, seed, maxBounces, globalIllumination, useLod); + + if (sampleCount[pidx] == 0) { + accumColor[pidx] = pbrColor; + sampleCount[pidx] = 1; + } else { + accumColor[pidx] += pbrColor; + sampleCount[pidx] += 1; + } + } + } + } + } + + #pragma omp parallel for schedule(static) + for (int pidx = 0; pidx < width * height; ++pidx) { + Eigen::Vector3f finalColor = accumColor[pidx]; + int count = sampleCount[pidx]; + + if (count > 0) { + finalColor /= static_cast(count); + } + + finalColor = finalColor.cwiseMax(0.0f).cwiseMin(1.0f); + + int idx = pidx * 3; + colorBuffer[idx] = finalColor[0]; + colorBuffer[idx + 1] = finalColor[1]; + colorBuffer[idx + 2] = finalColor[2]; + } + + outFrame.setData(colorBuffer, frame::colormap::RGB); + return outFrame; + } + std::vector> getExternalNodes(int targetObjectId) { std::vector> candidates; std::vector> surfaceNodes; @@ -1562,7 +1929,7 @@ public: PointType dirs[6] = {{1,0,0}, {-1,0,0}, {0,1,0}, {0,-1,0}, {0,0,1}, {0,0,-1}}; for(int i=0; i<6; ++i) { auto neighbor = find(node->position + dirs[i] * node->size, checkRad); - if(neighbor && neighbor->objectId == objectId && neighbor->active && neighbor->refraction < 0.01f) { + if(neighbor && neighbor->objectId == objectId && neighbor->active && neighbor->transmission < 0.01f) { hiddenSides++; } } @@ -1597,6 +1964,12 @@ public: return true; } + void optimize() { + if (root_) { + optimizeRecursive(root_.get()); + } + } + void printStats(std::ostream& os = std::cout) const { if (!root_) { os << "[Octree Stats] Tree is null/empty." << std::endl; diff --git a/util/sim/elementcontent.hpp b/util/sim/elementcontent.hpp new file mode 100644 index 0000000..afee099 --- /dev/null +++ b/util/sim/elementcontent.hpp @@ -0,0 +1,383 @@ +#ifndef ELEMENTS +#define ELEMENTS + +#include + +struct BaseElementProps { + float density; // kg/m^3 + float meltingPoint; // Kelvin + float boilingPoint; // Kelvin + float specificHeat; // J/(kg*K) + float electronegativity; // Pauling scale +}; + +static const std::array ELEMENT_DB = {{ + // 1: Hydrogen + {0.08988f, 14.01f, 20.28f, 14304.0f, 2.20f}, + // 2: Helium + {0.1785f, 0.95f, 4.22f, 5193.0f, 0.0f}, // No electronegativity, using 0 + // 3: Lithium + {534.0f, 453.69f, 1560.0f, 3582.0f, 0.98f}, + // 4: Beryllium + {1850.0f, 1560.0f, 2742.0f, 1825.0f, 1.57f}, + // 5: Boron + {2340.0f, 2349.0f, 4200.0f, 1026.0f, 2.04f}, + // 6: Carbon + {2267.0f, 4000.0f, 4300.0f, 709.0f, 2.55f}, + // 7: Nitrogen + {1.2506f, 63.15f, 77.36f, 1040.0f, 3.04f}, + // 8: Oxygen + {1.429f, 54.36f, 90.2f, 918.0f, 3.44f}, + // 9: Fluorine + {1.696f, 53.53f, 85.03f, 824.0f, 3.98f}, + // 10: Neon + {0.9002f, 24.56f, 27.07f, 1030.0f, 0.0f}, // No electronegativity + // 11: Sodium + {968.0f, 370.87f, 1156.0f, 1228.0f, 0.93f}, + // 12: Magnesium + {1738.0f, 923.0f, 1363.0f, 1023.0f, 1.31f}, + // 13: Aluminium + {2700.0f, 933.47f, 2792.0f, 897.0f, 1.61f}, + // 14: Silicon + {2329.0f, 1687.0f, 3538.0f, 705.0f, 1.9f}, + // 15: Phosphorus + {1823.0f, 317.3f, 550.0f, 769.0f, 2.19f}, + // 16: Sulfur + {2070.0f, 388.36f, 717.87f, 710.0f, 2.58f}, + // 17: Chlorine + {3.2f, 171.6f, 239.11f, 479.0f, 3.16f}, + // 18: Argon + {1.784f, 83.8f, 87.3f, 520.0f, 0.0f}, // No electronegativity + // 19: Potassium + {890.0f, 336.53f, 1032.0f, 757.0f, 0.82f}, + // 20: Calcium + {1550.0f, 1115.0f, 1757.0f, 647.0f, 1.0f}, + // 21: Scandium + {2985.0f, 1814.0f, 3109.0f, 568.0f, 1.36f}, + // 22: Titanium + {4506.0f, 1941.0f, 3560.0f, 523.0f, 1.54f}, + // 23: Vanadium + {6110.0f, 2183.0f, 3680.0f, 489.0f, 1.63f}, + // 24: Chromium + {7150.0f, 2180.0f, 2944.0f, 449.0f, 1.66f}, + // 25: Manganese + {7210.0f, 1519.0f, 2334.0f, 479.0f, 1.55f}, + // 26: Iron + {7874.0f, 1811.0f, 3134.0f, 449.0f, 1.83f}, + // 27: Cobalt + {8900.0f, 1768.0f, 3200.0f, 421.0f, 1.88f}, + // 28: Nickel + {8908.0f, 1728.0f, 3186.0f, 444.0f, 1.91f}, + // 29: Copper + {8960.0f, 1357.77f, 2835.0f, 385.0f, 1.9f}, + // 30: Zinc + {7140.0f, 692.88f, 1180.0f, 388.0f, 1.65f}, + // 31: Gallium + {5910.0f, 302.9146f, 2673.0f, 371.0f, 1.81f}, + // 32: Germanium + {5323.0f, 1211.4f, 3106.0f, 320.0f, 2.01f}, + // 33: Arsenic + {5727.0f, 1090.0f, 887.0f, 329.0f, 2.18f}, // Sublimes at 887K + // 34: Selenium + {4810.0f, 453.0f, 958.0f, 321.0f, 2.55f}, + // 35: Bromine + {3102.8f, 265.8f, 332.0f, 474.0f, 2.96f}, + // 36: Krypton + {3.749f, 115.79f, 119.93f, 248.0f, 3.0f}, + // 37: Rubidium + {1532.0f, 312.46f, 961.0f, 363.0f, 0.82f}, + // 38: Strontium + {2640.0f, 1050.0f, 1655.0f, 301.0f, 0.95f}, + // 39: Yttrium + {4472.0f, 1799.0f, 3609.0f, 298.0f, 1.22f}, + // 40: Zirconium + {6520.0f, 2128.0f, 4682.0f, 278.0f, 1.33f}, + // 41: Niobium + {8570.0f, 2750.0f, 5017.0f, 265.0f, 1.6f}, + // 42: Molybdenum + {10280.0f, 2896.0f, 4912.0f, 251.0f, 2.16f}, + // 43: Technetium + {11000.0f, 2430.0f, 4538.0f, 0.0f, 1.9f}, // No specific heat data + // 44: Ruthenium + {12450.0f, 2607.0f, 4423.0f, 238.0f, 2.2f}, + // 45: Rhodium + {12410.0f, 2237.0f, 3968.0f, 243.0f, 2.28f}, + // 46: Palladium + {12023.0f, 1828.05f, 3236.0f, 244.0f, 2.2f}, + // 47: Silver + {10490.0f, 1234.93f, 2435.0f, 235.0f, 1.93f}, + // 48: Cadmium + {8650.0f, 594.22f, 1040.0f, 232.0f, 1.69f}, + // 49: Indium + {7310.0f, 429.75f, 2345.0f, 233.0f, 1.78f}, + // 50: Tin + {7265.0f, 505.08f, 2875.0f, 228.0f, 1.96f}, + // 51: Antimony + {6697.0f, 903.78f, 1860.0f, 207.0f, 2.05f}, + // 52: Tellurium + {6240.0f, 722.66f, 1261.0f, 202.0f, 2.1f}, + // 53: Iodine + {4933.0f, 386.85f, 457.4f, 214.0f, 2.66f}, + // 54: Xenon + {5.894f, 161.4f, 165.03f, 158.0f, 2.6f}, + // 55: Caesium + {1930.0f, 301.59f, 944.0f, 242.0f, 0.79f}, + // 56: Barium + {3510.0f, 1000.0f, 2170.0f, 204.0f, 0.89f}, + // 57: Lanthanum + {6162.0f, 1193.0f, 3737.0f, 195.0f, 1.1f}, + // 58: Cerium + {6770.0f, 1068.0f, 3716.0f, 192.0f, 1.12f}, + // 59: Praseodymium + {6770.0f, 1208.0f, 3793.0f, 193.0f, 1.13f}, + // 60: Neodymium + {7010.0f, 1297.0f, 3347.0f, 190.0f, 1.14f}, + // 61: Promethium + {7260.0f, 1315.0f, 3273.0f, 0.0f, 1.13f}, // No specific heat data + // 62: Samarium + {7520.0f, 1345.0f, 2067.0f, 197.0f, 1.17f}, + // 63: Europium + {5244.0f, 1099.0f, 1802.0f, 182.0f, 1.2f}, + // 64: Gadolinium + {7900.0f, 1585.0f, 3546.0f, 236.0f, 1.2f}, + // 65: Terbium + {8230.0f, 1629.0f, 3503.0f, 182.0f, 1.2f}, + // 66: Dysprosium + {8540.0f, 1680.0f, 2840.0f, 170.0f, 1.22f}, + // 67: Holmium + {8790.0f, 1734.0f, 2993.0f, 165.0f, 1.23f}, + // 68: Erbium + {9066.0f, 1802.0f, 3141.0f, 168.0f, 1.24f}, + // 69: Thulium + {9320.0f, 1818.0f, 2223.0f, 160.0f, 1.25f}, + // 70: Ytterbium + {6900.0f, 1097.0f, 1469.0f, 155.0f, 1.1f}, + // 71: Lutetium + {9841.0f, 1925.0f, 3675.0f, 154.0f, 1.27f}, + // 72: Hafnium + {13310.0f, 2506.0f, 4876.0f, 144.0f, 1.3f}, + // 73: Tantalum + {16690.0f, 3290.0f, 5731.0f, 140.0f, 1.5f}, + // 74: Tungsten + {19250.0f, 3695.0f, 6203.0f, 132.0f, 2.36f}, + // 75: Rhenium + {21020.0f, 3459.0f, 5869.0f, 137.0f, 1.9f}, + // 76: Osmium + {22590.0f, 3306.0f, 5285.0f, 130.0f, 2.2f}, + // 77: Iridium + {22560.0f, 2719.0f, 4701.0f, 131.0f, 2.2f}, + // 78: Platinum + {21450.0f, 2041.4f, 4098.0f, 133.0f, 2.28f}, + // 79: Gold + {19300.0f, 1337.33f, 3129.0f, 129.0f, 2.54f}, + // 80: Mercury + {13534.0f, 234.43f, 629.88f, 140.0f, 2.0f}, + // 81: Thallium + {11850.0f, 577.0f, 1746.0f, 129.0f, 1.62f}, + // 82: Lead + {11340.0f, 600.61f, 2022.0f, 129.0f, 2.33f}, // Using 4+ value + // 83: Bismuth + {9780.0f, 544.7f, 1837.0f, 122.0f, 2.02f}, + // 84: Polonium + {9196.0f, 527.0f, 1235.0f, 0.0f, 2.0f}, // No specific heat data + // 85: Astatine + {8930.0f, 575.0f, 610.0f, 0.0f, 2.2f}, // Approx density, no specific heat + // 86: Radon + {9.73f, 202.0f, 211.3f, 94.0f, 2.2f}, + // 87: Francium + {2480.0f, 281.0f, 890.0f, 0.0f, 0.79f}, // Approx values + // 88: Radium + {5500.0f, 973.0f, 2010.0f, 94.0f, 0.9f}, + // 89: Actinium + {10000.0f, 1323.0f, 3471.0f, 120.0f, 1.1f}, + // 90: Thorium + {11700.0f, 2115.0f, 5061.0f, 113.0f, 1.3f}, + // 91: Protactinium + {15370.0f, 1841.0f, 4300.0f, 0.0f, 1.5f}, // No specific heat data + // 92: Uranium + {19100.0f, 1405.3f, 4404.0f, 116.0f, 1.38f}, + // 93: Neptunium + {20450.0f, 917.0f, 4273.0f, 0.0f, 1.36f}, // No specific heat data + // 94: Plutonium + {19850.0f, 912.5f, 3501.0f, 0.0f, 1.28f}, // No specific heat data + // 95: Americium + {12000.0f, 1449.0f, 2880.0f, 0.0f, 1.13f}, // No specific heat data + // 96: Curium + {13510.0f, 1613.0f, 3383.0f, 0.0f, 1.28f}, // No specific heat data + // 97: Berkelium + {14780.0f, 1259.0f, 2900.0f, 0.0f, 1.3f}, // No specific heat data + // 98: Californium + {15100.0f, 1173.0f, 1743.0f, 0.0f, 1.3f}, // No specific heat data + // 99: Einsteinium + {8840.0f, 1133.0f, 1269.0f, 0.0f, 1.3f}, // No specific heat data + // 100: Fermium + {9700.0f, 1125.0f, 1800.0f, 0.0f, 1.3f}, // Estimated values + // 101: Mendelevium + {10300.0f, 1100.0f, 0.0f, 0.0f, 1.3f}, // Estimated + // 102: Nobelium + {9900.0f, 1100.0f, 0.0f, 0.0f, 1.3f}, // Estimated + // 103: Lawrencium + {14400.0f, 1900.0f, 0.0f, 0.0f, 1.3f}, // Estimated + // 104: Rutherfordium + {17000.0f, 2400.0f, 5800.0f, 0.0f, 0.0f}, // Estimated + // 105: Dubnium + {21600.0f, 0.0f, 0.0f, 0.0f, 0.0f}, // Estimated + // 106: Seaborgium + {23500.0f, 0.0f, 0.0f, 0.0f, 0.0f}, // Estimated + // 107: Bohrium + {26500.0f, 0.0f, 0.0f, 0.0f, 0.0f}, // Estimated + // 108: Hassium + {28000.0f, 0.0f, 0.0f, 0.0f, 0.0f}, // Estimated + // 109: Meitnerium + {27500.0f, 0.0f, 0.0f, 0.0f, 0.0f}, // Estimated + // 110: Darmstadtium + {26500.0f, 0.0f, 0.0f, 0.0f, 0.0f}, // Estimated + // 111: Roentgenium + {23000.0f, 0.0f, 0.0f, 0.0f, 0.0f}, // Estimated + // 112: Copernicium + {14000.0f, 283.0f, 340.0f, 0.0f, 0.0f}, // Estimated + // 113: Nihonium + {16000.0f, 700.0f, 1400.0f, 0.0f, 0.0f}, // Estimated + // 114: Flerovium + {11400.0f, 284.0f, 0.0f, 0.0f, 0.0f}, // Estimated + // 115: Moscovium + {13500.0f, 700.0f, 1400.0f, 0.0f, 0.0f}, // Estimated + // 116: Livermorium + {12900.0f, 700.0f, 1100.0f, 0.0f, 0.0f}, // Estimated + // 117: Tennessine + {7200.0f, 700.0f, 883.0f, 0.0f, 0.0f}, // Estimated + // 118: Oganesson + {7000.0f, 325.0f, 450.0f, 0.0f, 0.0f} // Estimated +}}; + +struct PointProperties { + float weight = 0.0f; // Total mass + float density = 0.0f; // Mass / Volume + float meltingPoint = 0.0f; + float boilingPoint = 0.0f; + float specificHeat = 0.0f; + float electronegativity = 0.0f; +}; + +struct elementContent { + float hydrogen = 0.0f; + float helium = 0.0f; + float lithium = 0.0f; + float beryllium = 0.0f; + float boron = 0.0f; + float carbon = 0.0f; + float nitrogen = 0.0f; + float oxygen = 0.0f; + float fluorine = 0.0f; + float neon = 0.0f; + float sodium = 0.0f; + float magnesium = 0.0f; + float aluminum = 0.0f; + float silicon = 0.0f; + float phosporus = 0.0f; + float sulfur = 0.0f; + float chlorine = 0.0f; + float argon = 0.0f; + float potassium = 0.0f; + float calcium = 0.0f; + float scandium = 0.0f; + float titanium = 0.0f; + float vanadium = 0.0f; + float chromium = 0.0f; + float manganese = 0.0f; + float iron = 0.0f; + float cobalt = 0.0f; + float nickel = 0.0f; + float copper = 0.0f; + float zinc = 0.0f; + float gallium = 0.0f; + float germanium = 0.0f; + float arsenic = 0.0f; + float selenium = 0.0f; + float bromine = 0.0f; + float krypton = 0.0f; + float rubidium = 0.0f; + float strontium = 0.0f; + float yttrium = 0.0f; + float zirconium = 0.0f; + float niobium = 0.0f; + float molybdenum = 0.0f; + float technetium = 0.0f; + float ruthenium = 0.0f; + float rhodium = 0.0f; + float palladium = 0.0f; + float silver = 0.0f; + float cadmium = 0.0f; + float indium = 0.0f; + float tin = 0.0f; + float antimony = 0.0f; + float tellurium = 0.0f; + float iodine = 0.0f; + float xenon = 0.0f; + float caesium = 0.0f; + float barium = 0.0f; + float lanthanum = 0.0f; + float cerium = 0.0f; + float praseodymium = 0.0f; + float neodymium = 0.0f; + float promethium = 0.0f; + float samarium = 0.0f; + float europium = 0.0f; + float gadolinium = 0.0f; + float terbium = 0.0f; + float dysprosium = 0.0f; + float holmium = 0.0f; + float erbium = 0.0f; + float thulium = 0.0f; + float ytterbium = 0.0f; + float lutetium = 0.0f; + float hafnium = 0.0f; + float tantalum = 0.0f; + float tungsten = 0.0f; + float rhenium = 0.0f; + float osmium = 0.0f; + float iridium = 0.0f; + float platinum = 0.0f; + float gold = 0.0f; + float mercury = 0.0f; + float thallium = 0.0f; + float lead = 0.0f; + float bismuth = 0.0f; + float polonium = 0.0f; + float astatine = 0.0f; + float radon = 0.0f; + float francium = 0.0f; + float radium = 0.0f; + float actinium = 0.0f; + float thorium = 0.0f; + float protactinium = 0.0f; + float uranium = 0.0f; + float neptunium = 0.0f; + float plutonium = 0.0f; + float americium = 0.0f; + float curium = 0.0f; + float berkelium = 0.0f; + float californium = 0.0f; + float einsteinium = 0.0f; + float fermium = 0.0f; + float mendelevium = 0.0f; + float nobelium = 0.0f; + float lawrencium = 0.0f; + float rutherfordium = 0.0f; + float dubnium = 0.0f; + float seaborgium = 0.0f; + float bohrium = 0.0f; + float hassium = 0.0f; + float meitnerium = 0.0f; + float darmstadtium = 0.0f; + float roentgenium = 0.0f; + float cpernicium = 0.0f; + float nihnium = 0.0f; + float flerovium = 0.0f; + float moscovium = 0.0f; + float livermorium = 0.0f; + float tennessine = 0.0f; + float oganesson = 0.0f; +}; + +#endif \ No newline at end of file diff --git a/util/sim/planet.hpp b/util/sim/planet.hpp index 3ced496..b349249 100644 --- a/util/sim/planet.hpp +++ b/util/sim/planet.hpp @@ -41,29 +41,29 @@ struct Particle { Eigen::Vector3f currentPos; float plateDisplacement = 0.0f; - float temperature = -1; - float water = -1; + // float temperature = -1; + // float water = -1; Eigen::Vector3f originColor; bool surface = false; //gravity factors: - Eigen::Matrix velocity; - Eigen::Matrix acceleration; - Eigen::Matrix forceAccumulator; - float density = 0.0f; - float pressure = 0.0f; - Eigen::Matrix pressureForce; - float viscosity = 0.5f; - Eigen::Matrix viscosityForce; - float restitution = 5.0f; - float mass; - bool isStatic = false; - float soundSpeed = 100.0f; - float sandcontent = 0.0f; - float siltcontent = 0.0f; - float claycontent = 0.0f; - float rockcontent = 0.0f; - float metalcontent = 0.0f; + // Eigen::Matrix velocity; + // Eigen::Matrix acceleration; + // Eigen::Matrix forceAccumulator; + // float density = 0.0f; + // float pressure = 0.0f; + // Eigen::Matrix pressureForce; + // float viscosity = 0.5f; + // Eigen::Matrix viscosityForce; + // float restitution = 5.0f; + // float mass; + // bool isStatic = false; + // float soundSpeed = 100.0f; + // float sandcontent = 0.0f; + // float siltcontent = 0.0f; + // float claycontent = 0.0f; + // float rockcontent = 0.0f; + // float metalcontent = 0.0f; std::unordered_map neighbors; @@ -83,7 +83,7 @@ struct planetConfig { float displacementStrength = 200.0f; std::vector surfaceNodes; std::vector interpolatedNodes; - float noiseStrength = 1.0f; + float noiseStrength = 10.0f; int numPlates = 15; int smoothingPasses = 3; float mountHeight = 250.0f; @@ -207,9 +207,10 @@ public: Eigen::Vector3f normal = p.originalPos.normalized(); p.noisePos = p.originalPos + (normal * displacementValue * config.noiseStrength); p.currentPos = p.noisePos; - - grid.update(oldPos, p.currentPos, p, true, p.originColor, config.voxelSize, true, -2, false, 0.0f, 0.0f, 0.0f); + grid.move(oldPos, p.currentPos); + grid.update(p.currentPos, p); } + grid.optimize(); } void assignSeeds() { @@ -698,8 +699,10 @@ public: for (auto& p : config.surfaceNodes) { Eigen::Vector3f oldPos = p.currentPos; p.currentPos = p.tectonicPos; - grid.update(oldPos, p.currentPos, p, true, p.originColor, config.voxelSize, true, -2, false, 0.0f, 0.0f, 0.0f); + grid.move(oldPos, p.currentPos); + grid.update(p.currentPos, p); } + grid.optimize(); std::cout << "Finalize apply results completed." << std::endl; } @@ -787,9 +790,7 @@ public: newPt.originColor = p3.originColor; } - v3 interpNormal = (p1.originalPos.normalized() * w1 + - p2.originalPos.normalized() * w2 + - p3.originalPos.normalized() * w3); + v3 interpNormal = (p1.originalPos.normalized() * w1 + p2.originalPos.normalized() * w2 + p3.originalPos.normalized() * w3); interpNormal.normalize(); float r1 = p1.currentPos.norm(); @@ -813,6 +814,7 @@ public: } } } + grid.optimize(); std::cout << "Interpolated " << counter << " surface gaps." << std::endl; }