From 868ae90fbeb781a125d4670cafe69e2d1927e061 Mon Sep 17 00:00:00 2001 From: Marc Gilleron Date: Mon, 2 May 2022 19:14:12 +0100 Subject: [PATCH] Added VoxelMeshSDF and VoxelToolLodTerrain.stamp_sdf() --- config.py | 1 + constants/voxel_constants.h | 1 + doc/source/changelog.md | 2 + edition/funcs.h | 28 + edition/mesh_sdf.cpp | 1001 ++++++++++++++++++++++++++++ edition/mesh_sdf.h | 142 ++++ edition/voxel_mesh_sdf_gd.cpp | 383 +++++++++++ edition/voxel_mesh_sdf_gd.h | 87 +++ edition/voxel_tool_lod_terrain.cpp | 97 +++ edition/voxel_tool_lod_terrain.h | 2 + register_types.cpp | 2 + storage/voxel_buffer_internal.cpp | 1 + storage/voxel_buffer_internal.h | 9 +- util/math/color8.h | 1 + util/math/conv.h | 2 + util/math/funcs.h | 4 + util/math/vector3f.h | 22 +- 17 files changed, 1783 insertions(+), 2 deletions(-) create mode 100644 edition/mesh_sdf.cpp create mode 100644 edition/mesh_sdf.h create mode 100644 edition/voxel_mesh_sdf_gd.cpp create mode 100644 edition/voxel_mesh_sdf_gd.h diff --git a/config.py b/config.py index 28c4a020..967b15ba 100644 --- a/config.py +++ b/config.py @@ -69,6 +69,7 @@ def get_doc_classes(): "VoxelBlockSerializer", "VoxelVoxLoader", "VoxelDataBlockEnterInfo", + "VoxelMeshSDF", "ZN_FastNoiseLite", "ZN_FastNoiseLiteGradient", diff --git a/constants/voxel_constants.h b/constants/voxel_constants.h index 6171d8d1..81a58c42 100644 --- a/constants/voxel_constants.h +++ b/constants/voxel_constants.h @@ -26,6 +26,7 @@ static const unsigned int MAX_VOLUME_SIZE = 2 * MAX_VOLUME_EXTENT; // 1,073,741, static const float INV_0x7f = 1.f / 0x7f; static const float INV_0x7fff = 1.f / 0x7fff; static const float INV_TAU = 1.f / Math_TAU; +static const float SQRT3 = 1.73205080757; // Below 32 bits, channels are normalized in -1..1, and can represent a limited number of values. // For storing SDF, we need a range of values that extends beyond that, in particular for better LOD. diff --git a/doc/source/changelog.md b/doc/source/changelog.md index 85191c84..2b97f270 100644 --- a/doc/source/changelog.md +++ b/doc/source/changelog.md @@ -19,6 +19,7 @@ Godot 4 is required from this version. - Added experimental support functions to help setting up basic multiplayer with `VoxelTerrain` (might change in the future) - Improved support for 64-bit floats - Added `ZN_ThreadedTask` to allow running custom tasks using the thread pool system + - Added `VoxelMeshSDF` to bake SDF from meshes, which can be used in voxel sculpting. - `VoxelGeneratorGraph`: added support for outputting to the TYPE channel, allowing use with `VoxelMesherBlocky` - `VoxelGeneratorGraph`: editor: unconnected inputs show their default value directly on the node - `VoxelGeneratorGraph`: editor: allow to change the axes on preview nodes 3D slices @@ -36,6 +37,7 @@ Godot 4 is required from this version. - `VoxelLodTerrain`: Editor: added option to show octree nodes in editor - `VoxelLodTerrain`: Added option to run a major part of the process logic into another thread - `VoxelToolLodTerrain`: added *experimental* `do_sphere_async`, an alternative version of `do_sphere` which defers the task on threads to reduce stutter if the affected area is big. + - `VoxelToolLodTerrain`: added `stamp_sdf` function to place a baked mesh SDF on the terrain - `VoxelInstancer`: Allow to dump VoxelInstancer as scene for debug inspection - `VoxelInstancer`: Editor: instance chunks are shown when the node is selected - `VoxelInstanceLibraryMultiMeshItem`: Support setting up mesh LODs from a scene with name `LODx` suffixes diff --git a/edition/funcs.h b/edition/funcs.h index 89f94960..35e80cb0 100644 --- a/edition/funcs.h +++ b/edition/funcs.h @@ -3,7 +3,9 @@ #include "../storage/funcs.h" #include "../util/fixed_array.h" +#include "../util/math/sdf.h" #include "../util/math/vector3.h" +#include "../util/math/vector3f.h" namespace zylann::voxel { @@ -109,6 +111,32 @@ inline void blend_texture_packed_u16( } } +// Interpolates values from a 3D grid at a given position, using trilinear interpolation. +// If the position is outside the grid, values are clamped. +inline float interpolate_trilinear(Span grid, const Vector3i res, const Vector3f pos) { + const Vector3f pfi = math::floor(pos - Vector3f(0.5f)); + // TODO Clamp pf too somehow? + const Vector3f pf = pos - pfi; + const Vector3i max_pos = math::max(res - Vector3i(2, 2, 2), Vector3i()); + const Vector3i pi = math::clamp(Vector3i(pfi.x, pfi.y, pfi.z), Vector3i(), max_pos); + + const unsigned int n010 = 1; + const unsigned int n100 = res.y; + const unsigned int n001 = res.y * res.x; + + const unsigned int i000 = pi.x * n100 + pi.y * n010 + pi.z * n001; + const unsigned int i010 = i000 + n010; + const unsigned int i100 = i000 + n100; + const unsigned int i001 = i000 + n001; + const unsigned int i110 = i010 + n100; + const unsigned int i011 = i010 + n001; + const unsigned int i101 = i100 + n001; + const unsigned int i111 = i110 + n001; + + return math::interpolate_trilinear( + grid[i000], grid[i100], grid[i101], grid[i001], grid[i010], grid[i110], grid[i111], grid[i011], pf); +} + } // namespace zylann::voxel namespace zylann::voxel::ops { diff --git a/edition/mesh_sdf.cpp b/edition/mesh_sdf.cpp new file mode 100644 index 00000000..b10887f8 --- /dev/null +++ b/edition/mesh_sdf.cpp @@ -0,0 +1,1001 @@ +#include "mesh_sdf.h" +#include "../util/godot/funcs.h" +#include "../util/math/box3i.h" +#include "../util/math/conv.h" +#include "../util/profiling.h" +#include "../util/string_funcs.h" // Debug + +namespace zylann::voxel::mesh_sdf { + +// Some papers for eventual improvements +// Jump flood +// https://www.comp.nus.edu.sg/%7Etants/jfa/i3d06.pdf +// GPU-based technique +// https://www.researchgate.net/profile/Bastian-Krayer/publication/332921884_Generating_signed_distance_fields_on_the_GPU_with_ray_maps/links/5d63d921299bf1f70b0de26b/Generating-signed-distance-fields-on-the-GPU-with-ray-maps.pdf +// In the future, maybe also gather gradients so it can be used efficiently with Dual Contouring? + +//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +float get_max_sdf_variation(Vector3f min_pos, Vector3f max_pos, Vector3i res) { + const Vector3f box_size = max_pos - min_pos; + const Vector3f max_variation_v = box_size / to_vec3f(res); + return math::max(max_variation_v.x, math::max(max_variation_v.y, max_variation_v.z)); +} + +void fix_sdf_sign_from_boundary(Span sdf_grid, Vector3i res, Vector3f min_pos, Vector3f max_pos) { + ZN_PROFILE_SCOPE(); + + if (res.x == 0 || res.y == 0 || res.z == 0) { + return; + } + + std::vector seeds; + + FixedArray dirs; + dirs[0] = Vector3i(-1, 0, 0); + dirs[1] = Vector3i(1, 0, 0); + dirs[2] = Vector3i(0, -1, 0); + dirs[3] = Vector3i(0, 1, 0); + dirs[4] = Vector3i(0, 0, -1); + dirs[5] = Vector3i(0, 0, 1); + + std::vector visited; + visited.resize(sdf_grid.size(), 0); + + const float tolerance = 1.1f; + const float max_variation = tolerance * get_max_sdf_variation(min_pos, max_pos, res); + + seeds.push_back(Vector3i()); + visited[0] = 1; + + // Spread positive sign from boundary. + // Has a bit of room for optimization, but profiling shows it accounts for a very small portion of time compared to + // calculating the distance field. + while (seeds.size() > 0) { + const Vector3i pos = seeds.back(); + seeds.pop_back(); + const unsigned int loc = Vector3iUtil::get_zxy_index(pos, res); + const float v = sdf_grid[loc]; + + for (unsigned int dir = 0; dir < dirs.size(); ++dir) { + const Vector3i npos = pos + dirs[dir]; + + if (npos.x < 0 || npos.y < 0 || npos.z < 0 || npos.x >= res.x || npos.y >= res.y || npos.z >= res.z) { + continue; + } + + const unsigned int nloc = Vector3iUtil::get_zxy_index(npos, res); + + ZN_ASSERT(nloc < visited.size()); + if (visited[nloc] == 1) { + continue; + } + + visited[nloc] = 1; + + ZN_ASSERT(nloc < sdf_grid.size()); + const float nv = sdf_grid[nloc]; + + if ((nv > 0.f) != (v > 0.f) && Math::abs(nv - v) < max_variation) { + continue; + } + + if (nv < 0.f) { + sdf_grid[nloc] = -nv; + } + seeds.push_back(npos); + } + } +} + +void partition_triangles( + int subdiv, Span triangles, Vector3f min_pos, Vector3f max_pos, ChunkGrid &chunk_grid) { + ZN_PROFILE_SCOPE(); + + const Vector3f mesh_size = max_pos - min_pos; + const float cs = math::max(mesh_size.x, math::max(mesh_size.y, mesh_size.z)) / subdiv; + const Vector3i grid_min = to_vec3i(math::floor(min_pos / cs)); + const Vector3i grid_max = to_vec3i(math::ceil(max_pos / cs)); + + chunk_grid.size = grid_max - grid_min; + chunk_grid.chunks.resize(Vector3iUtil::get_volume(chunk_grid.size)); + chunk_grid.chunk_size = cs; + chunk_grid.min_pos = min_pos; + + // Initialize chunk positions + { + Vector3i cpos; + for (cpos.z = 0; cpos.z < chunk_grid.size.z; ++cpos.z) { + for (cpos.x = 0; cpos.x < chunk_grid.size.x; ++cpos.x) { + cpos.y = 0; + unsigned int ci = Vector3iUtil::get_zxy_index(cpos, chunk_grid.size); + + for (; cpos.y < chunk_grid.size.y; ++cpos.y) { + ZN_ASSERT(ci < chunk_grid.chunks.size()); + Chunk &chunk = chunk_grid.chunks[ci]; + + chunk.pos = cpos; + ++ci; + } + } + } + } + + // Group triangles overlapping chunks + { + ZN_PROFILE_SCOPE_NAMED("Group triangles"); + + for (unsigned int triangle_index = 0; triangle_index < triangles.size(); ++triangle_index) { + const Triangle &t = triangles[triangle_index]; + + const Vector3f tri_min_pos = math::min(t.v1, math::min(t.v2, t.v3)); + const Vector3f tri_max_pos = math::max(t.v1, math::max(t.v2, t.v3)); + + const Vector3i tri_min_pos_i = to_vec3i(math::floor(tri_min_pos / cs)) - grid_min; + const Vector3i tri_max_pos_i = to_vec3i(math::floor(tri_max_pos / cs)) - grid_min; + + Vector3i cpos; + for (cpos.z = tri_min_pos_i.z; cpos.z <= tri_max_pos_i.z; ++cpos.z) { + for (cpos.x = tri_min_pos_i.x; cpos.x <= tri_max_pos_i.x; ++cpos.x) { + cpos.y = tri_min_pos_i.y; + unsigned int ci = Vector3iUtil::get_zxy_index(cpos, chunk_grid.size); + + for (; cpos.y <= tri_max_pos_i.y; ++cpos.y) { + ZN_ASSERT(ci < chunk_grid.chunks.size()); + Chunk &chunk = chunk_grid.chunks[ci]; + chunk.triangles.push_back(&t); + ++ci; + } + } + } + } + } + + // Gather chunks containing triangles + std::vector nonempty_chunks; + for (auto it = chunk_grid.chunks.begin(); it != chunk_grid.chunks.end(); ++it) { + const Chunk &chunk = *it; + if (chunk.triangles.size() > 0) { + nonempty_chunks.push_back(&chunk); + } + } + + // Find close chunks + { + ZN_PROFILE_SCOPE_NAMED("Group chunks"); + + Vector3i cpos; + for (cpos.z = 0; cpos.z < chunk_grid.size.z; ++cpos.z) { + for (cpos.x = 0; cpos.x < chunk_grid.size.x; ++cpos.x) { + cpos.y = 0; + unsigned int ci = Vector3iUtil::get_zxy_index(cpos, chunk_grid.size); + + for (; cpos.y < chunk_grid.size.y; ++cpos.y, ++ci) { + ZN_ASSERT(ci < chunk_grid.chunks.size()); + Chunk &chunk = chunk_grid.chunks[ci]; + + // Find closest chunk + + const Chunk *closest_chunk = nullptr; + // Distance is in chunks + int closest_chunk_distance_squared = 0x0fffffff; + + if (chunk.triangles.size() > 0) { + closest_chunk = &chunk; + closest_chunk_distance_squared = 0; + + } else { + for (auto it = nonempty_chunks.begin(); it != nonempty_chunks.end(); ++it) { + const Chunk &nchunk = **it; + const int distance_squared = (nchunk.pos - chunk.pos).length_squared(); + if (distance_squared < closest_chunk_distance_squared) { + closest_chunk = &nchunk; + closest_chunk_distance_squared = distance_squared; + } + } + } + + ZN_ASSERT(closest_chunk != nullptr); + + // Find other close chunks slightly beyond the closest chunk. + // This is to account for the fact the closest chunk might contain a triangle further away than + // a closer triangle found in a farther chunk. + + const int margin_distance_squared = + math::squared(sqrtf(closest_chunk_distance_squared) + constants::SQRT3); + + for (auto it = nonempty_chunks.begin(); it != nonempty_chunks.end(); ++it) { + const Chunk &nchunk = **it; + const int distance_squared = (nchunk.pos - chunk.pos).length_squared(); + // Note, this will include the closest chunk + if (distance_squared <= margin_distance_squared) { + chunk.near_chunks.push_back(&nchunk); + } + } + + // println(format("Chunk {} has {} / {} near non-empty chunks", chunk.pos, chunk.near_chunks.size(), + // nonempty_chunks.size())); + } + } + } + } +} + +/* +// Non-optimized version, suitable for single queries +float get_distance_to_triangle_squared(const Vector3f v1, const Vector3f v2, const Vector3f v3, const Vector3f p) { + // https://iquilezles.org/articles/triangledistance/ + + const Vector3f v21 = v2 - v1; + const Vector3f v32 = v3 - v2; + const Vector3f v13 = v1 - v3; + + const Vector3f p1 = p - v1; + const Vector3f p2 = p - v2; + const Vector3f p3 = p - v3; + + const Vector3f nor = v21.cross(v13); + + const float det = // + signf(v21.cross(nor).dot(p1)) + // + signf(v32.cross(nor).dot(p2)) + // + signf(v13.cross(nor).dot(p3)); + + if (det < 2.f) { + // Outside of the prism: get distance to closest edge + return math::min( + math::min( // + (v21 * math::clamp(v21.dot(p1) / v21.length_squared(), 0.f, 1.f) - p1).length_squared(), + (v32 * math::clamp(v32.dot(p2) / v32.length_squared(), 0.f, 1.f) - p2).length_squared()), + (v13 * math::clamp(v13.dot(p3) / v13.length_squared(), 0.f, 1.f) - p3).length_squared()); + } else { + // Inside the prism: get distance to plane + return math::squared(nor.dot(p1)) / nor.length_squared(); + } +} +*/ + +// Returns the distance from a point to a triangle, where some terms are precalculated. +// This may be preferred if the same triangles have to be queried many times. +float get_distance_to_triangle_squared_precalc(const Triangle &t, const Vector3f p) { + // https://iquilezles.org/articles/triangledistance/ + + const Vector3f p1 = p - t.v1; + const Vector3f p2 = p - t.v2; + const Vector3f p3 = p - t.v3; + + const float det = // + math::sign(t.v21_cross_nor.dot(p1)) + // + math::sign(t.v32_cross_nor.dot(p2)) + // + math::sign(t.v13_cross_nor.dot(p3)); + + if (det < 2.f) { + // Outside of the prism: get distance to closest edge + return math::min( + math::min( // + (t.v21 * math::clamp(t.v21.dot(p1) * t.inv_v21_length_squared, 0.f, 1.f) - p1).length_squared(), + (t.v32 * math::clamp(t.v32.dot(p2) * t.inv_v32_length_squared, 0.f, 1.f) - p2) + .length_squared()), + (t.v13 * math::clamp(t.v13.dot(p3) * t.inv_v13_length_squared, 0.f, 1.f) - p3).length_squared()); + } else { + // Inside the prism: get distance to plane + return math::squared(t.nor.dot(p1)) * t.inv_nor_length_squared; + } +} + +void precalc_triangles(Span triangles) { + ZN_PROFILE_SCOPE(); + for (size_t i = 0; i < triangles.size(); ++i) { + Triangle &t = triangles[i]; + t.v21 = t.v2 - t.v1; + t.v32 = t.v3 - t.v2; + t.v13 = t.v1 - t.v3; + t.nor = t.v21.cross(t.v13); + t.v21_cross_nor = t.v21.cross(t.nor); + t.v32_cross_nor = t.v32.cross(t.nor); + t.v13_cross_nor = t.v13.cross(t.nor); + t.inv_v21_length_squared = 1.f / t.v21.length_squared(); + t.inv_v32_length_squared = 1.f / t.v32.length_squared(); + t.inv_v13_length_squared = 1.f / t.v13.length_squared(); + t.inv_nor_length_squared = 1.f / t.nor.length_squared(); + } +} + +unsigned int get_closest_triangle_precalc(const Vector3f pos, Span triangles) { + float min_distance_squared = 9999999.f; + size_t closest_tri_index = triangles.size(); + + for (size_t i = 0; i < triangles.size(); ++i) { + const Triangle &t = triangles[i]; + //const float sqd = get_distance_to_triangle_squared(t.v1, t.v2, t.v3, pos); + const float sqd = get_distance_to_triangle_squared_precalc(t, pos); + + if (sqd < min_distance_squared) { + min_distance_squared = sqd; + closest_tri_index = i; + } + } + + return closest_tri_index; +} + +float get_mesh_signed_distance_at(const Vector3f pos, Span triangles /*, bool debug = false*/) { + float min_distance_squared = 9999999.f; + size_t closest_tri_index = 0; + + for (size_t i = 0; i < triangles.size(); ++i) { + const Triangle &t = triangles[i]; + //const float sqd = get_distance_to_triangle_squared(t.v1, t.v2, t.v3, pos); + const float sqd = get_distance_to_triangle_squared_precalc(t, pos); + + // TODO What if two triangles of opposite directions share the same point? + // If the distance comes from that point, it makes finding the sign ambiguous. + // sometimes two triangles of opposite directions share an edge or point, and there is no quick way + // to figure out which one must be taken. + // For now this is worked around in a later pass with a floodfill. + + /*if (debug) { + if (sqd < min_distance_squared + 0.01f) { + const Vector3f plane_normal = t.v13.cross(t.v1 - t.v2).normalized(); + const float plane_d = plane_normal.dot(t.v1); + float s; + if (plane_normal.dot(pos) > plane_d) { + s = 1.f; + } else { + s = -1.f; + } + println(format("Candidate tri {} with sqd {} sign {}", i, sqd, s)); + } + }*/ + + if (sqd < min_distance_squared) { + /*if (debug) { + println(format("Better tri {} with sqd {}", i, sqd)); + }*/ + min_distance_squared = sqd; + closest_tri_index = i; + } + } + + const Triangle &ct = triangles[closest_tri_index]; + const float d = Math::sqrt(min_distance_squared); + + //if (p_dir == CLOCKWISE) { + //const Vector3f plane_normal = (ct.v1 - ct.v3).cross(ct.v1 - ct.v2).normalized(); + const Vector3f plane_normal = ct.v13.cross(ct.v1 - ct.v2).normalized(); + //} else { + // normal = (p_point1 - p_point2).cross(p_point1 - p_point3); + //} + const float plane_d = plane_normal.dot(ct.v1); + + if (plane_normal.dot(pos) > plane_d) { + return d; + } + return -d; +} + +float get_mesh_signed_distance_at(const Vector3f pos, const ChunkGrid &chunk_grid) { + float min_distance_squared = 9999999.f; + const Triangle *closest_tri = nullptr; + + const Vector3i chunk_pos = to_vec3i(math::floor((pos - chunk_grid.min_pos) / chunk_grid.chunk_size)); + const unsigned chunk_index = Vector3iUtil::get_zxy_index(chunk_pos, chunk_grid.size); + ZN_ASSERT(chunk_index < chunk_grid.chunks.size()); + const Chunk &chunk = chunk_grid.chunks[chunk_index]; + + for (auto near_chunk_it = chunk.near_chunks.begin(); near_chunk_it != chunk.near_chunks.end(); ++near_chunk_it) { + const Chunk &near_chunk = **near_chunk_it; + + for (auto tri_it = near_chunk.triangles.begin(); tri_it != near_chunk.triangles.end(); ++tri_it) { + const Triangle &t = **tri_it; + + // TODO What if two triangles of opposite directions share the same point? + // If the distance comes from that point, it makes finding the sign ambiguous. + // sometimes two triangles of opposite directions share an edge or point, and there is no quick way + // to figure out which one must be taken. + // For now this is worked around in a later pass with a floodfill. + + //const float sqd = get_distance_to_triangle_squared(t.v1, t.v2, t.v3, pos); + const float sqd = get_distance_to_triangle_squared_precalc(t, pos); + + if (sqd < min_distance_squared) { + min_distance_squared = sqd; + closest_tri = &t; + } + } + } + + const float d = Math::sqrt(min_distance_squared); + ZN_ASSERT(closest_tri != nullptr); + + //if (p_dir == CLOCKWISE) { + //const Vector3f plane_normal = (ct.v1 - ct.v3).cross(ct.v1 - ct.v2).normalized(); + const Vector3f plane_normal = closest_tri->v13.cross(closest_tri->v1 - closest_tri->v2).normalized(); + //} else { + // normal = (p_point1 - p_point2).cross(p_point1 - p_point3); + //} + const float plane_d = plane_normal.dot(closest_tri->v1); + + if (plane_normal.dot(pos) > plane_d) { + return d; + } + return -d; +} + +struct Evaluator { + Span triangles; + const Vector3i res; + const Vector3f min_pos; + const Vector3f mesh_size; + const Vector3f hcs; + + inline Vector3f get_pos(const Vector3i grid_pos) const { + const Vector3f pf(float(grid_pos.x) / res.x, float(grid_pos.y) / res.y, float(grid_pos.z) / res.z); + const Vector3f pos = min_pos + mesh_size * pf + hcs; + return pos; + } + + inline float operator()(const Vector3i grid_pos) const { + return get_mesh_signed_distance_at(get_pos(grid_pos), triangles); + } +}; + +void generate_mesh_sdf_approx_interp(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos) { + ZN_PROFILE_SCOPE(); + + static const float FAR_SD = 9999999.f; + + const Vector3f mesh_size = max_pos - min_pos; + const Vector3f cell_size = mesh_size / Vector3f(res.x, res.y, res.z); + const Vector3f hcs(cell_size * 0.5f); + const Evaluator eval{ triangles, res, min_pos, mesh_size, hcs }; + + const unsigned int node_size_po2 = 2; + const unsigned int node_size_cells = 1 << node_size_po2; + + const Vector3i node_grid_size = math::ceildiv(res, node_size_cells) + Vector3i(1, 1, 1); + const float node_size = node_size_cells * cell_size.length(); + const float node_subdiv_threshold = 0.6f * node_size; + + std::vector node_grid; + node_grid.resize(Vector3iUtil::get_volume(node_grid_size)); + + // Fill SDF grid with far distances as "infinity", we'll use that to check if we computed it already + sdf_grid.fill(FAR_SD); + + // Evaluate SDF at the corners of nodes + Vector3i node_pos; + for (node_pos.z = 0; node_pos.z < node_grid_size.z; ++node_pos.z) { + for (node_pos.x = 0; node_pos.x < node_grid_size.x; ++node_pos.x) { + for (node_pos.y = 0; node_pos.y < node_grid_size.y; ++node_pos.y) { + const Vector3i gp000 = node_pos << node_size_po2; + + const size_t ni = Vector3iUtil::get_zxy_index(node_pos, node_grid_size); + ZN_ASSERT(ni < node_grid.size()); + const float sd = eval(gp000); + node_grid[ni] = sd; + + if (Box3i(Vector3i(), res).contains(gp000)) { + const size_t i = Vector3iUtil::get_zxy_index(gp000, res); + ZN_ASSERT(i < sdf_grid.size()); + sdf_grid[i] = sd; + } + } + } + } + + // Precompute flat-grid neighbor offsets + const unsigned int ni000 = 0; + const unsigned int ni100 = Vector3iUtil::get_zxy_index(Vector3i(1, 0, 0), node_grid_size); + const unsigned int ni010 = Vector3iUtil::get_zxy_index(Vector3i(0, 1, 0), node_grid_size); + const unsigned int ni110 = Vector3iUtil::get_zxy_index(Vector3i(1, 1, 0), node_grid_size); + const unsigned int ni001 = Vector3iUtil::get_zxy_index(Vector3i(0, 0, 1), node_grid_size); + const unsigned int ni101 = Vector3iUtil::get_zxy_index(Vector3i(1, 0, 1), node_grid_size); + const unsigned int ni011 = Vector3iUtil::get_zxy_index(Vector3i(0, 1, 1), node_grid_size); + const unsigned int ni111 = Vector3iUtil::get_zxy_index(Vector3i(1, 1, 1), node_grid_size); + + // Then for every node + for (node_pos.z = 0; node_pos.z < node_grid_size.z - 1; ++node_pos.z) { + for (node_pos.x = 0; node_pos.x < node_grid_size.x - 1; ++node_pos.x) { + for (node_pos.y = 0; node_pos.y < node_grid_size.y - 1; ++node_pos.y) { + float ud = FAR_SD; + + const unsigned int ni = Vector3iUtil::get_zxy_index(node_pos, node_grid_size); + + // Get signed distance at each corner we computed earlier + const float sd000 = node_grid[ni]; + const float sd100 = node_grid[ni + ni100]; + const float sd010 = node_grid[ni + ni010]; + const float sd110 = node_grid[ni + ni110]; + const float sd001 = node_grid[ni + ni001]; + const float sd101 = node_grid[ni + ni101]; + const float sd011 = node_grid[ni + ni011]; + const float sd111 = node_grid[ni + ni111]; + + // Get smallest one + ud = math::min(ud, Math::abs(sd000)); + ud = math::min(ud, Math::abs(sd100)); + ud = math::min(ud, Math::abs(sd010)); + ud = math::min(ud, Math::abs(sd110)); + ud = math::min(ud, Math::abs(sd001)); + ud = math::min(ud, Math::abs(sd101)); + ud = math::min(ud, Math::abs(sd011)); + ud = math::min(ud, Math::abs(sd111)); + + const Box3i cell_box = Box3i(node_pos * node_size_cells, Vector3iUtil::create(node_size_cells)) + .clipped(Box3i(Vector3i(), res)); + + // If the minimum distance at the corners of the node is lower than the threshold, + // subdivide the node. + if (ud < node_subdiv_threshold) { + // Full-res SDF + cell_box.for_each_cell_zxy([&sdf_grid, eval](const Vector3i grid_pos) { + const size_t i = Vector3iUtil::get_zxy_index(grid_pos, eval.res); + if (sdf_grid[i] != FAR_SD) { + // Already computed + return; + } + ZN_ASSERT(i < sdf_grid.size()); + sdf_grid[i] = eval(grid_pos); + }); + } else { + // We are far enough from the surface, approximate by interpolating corners + const Vector3i cell_box_end = cell_box.pos + cell_box.size; + Vector3i grid_pos; + for (grid_pos.z = cell_box.pos.z; grid_pos.z < cell_box_end.z; ++grid_pos.z) { + for (grid_pos.x = cell_box.pos.x; grid_pos.x < cell_box_end.x; ++grid_pos.x) { + for (grid_pos.y = cell_box.pos.y; grid_pos.y < cell_box_end.y; ++grid_pos.y) { + const size_t i = Vector3iUtil::get_zxy_index(grid_pos, eval.res); + if (sdf_grid[i] != FAR_SD) { + // Already computed + continue; + } + const Vector3f ipf = to_vec3f(grid_pos - cell_box.pos) / float(node_size_cells); + const float sd = math::interpolate_trilinear( + sd000, sd100, sd101, sd001, sd010, sd110, sd111, sd011, ipf); + ZN_ASSERT(i < sdf_grid.size()); + sdf_grid[i] = sd; + } + } + } + } + } + } + } +} + +void generate_mesh_sdf(Span sdf_grid, const Vector3i res, const Box3i sub_box, Span triangles, + const Vector3f min_pos, const Vector3f max_pos) { + ZN_PROFILE_SCOPE(); + ZN_ASSERT(Box3i(Vector3i(), res).contains(sub_box)); + ZN_ASSERT(sdf_grid.size() == Vector3iUtil::get_volume(res)); + + const Vector3f mesh_size = max_pos - min_pos; + const Vector3f cell_size = mesh_size / Vector3f(res.x, res.y, res.z); + + // ProfilingClock profiling_clock; + + Vector3i grid_pos; + const Vector3f hcs(cell_size * 0.5f); + + const Vector3i sub_box_end = sub_box.pos + sub_box.size; + + for (grid_pos.z = sub_box.pos.z; grid_pos.z < sub_box_end.z; ++grid_pos.z) { + for (grid_pos.x = sub_box.pos.x; grid_pos.x < sub_box_end.x; ++grid_pos.x) { + grid_pos.y = sub_box.pos.y; + size_t grid_index = Vector3iUtil::get_zxy_index(grid_pos, res); + + for (; grid_pos.y < sub_box_end.y; ++grid_pos.y) { + const Vector3f pf(float(grid_pos.x) / res.x, float(grid_pos.y) / res.y, float(grid_pos.z) / res.z); + const Vector3f pos = min_pos + mesh_size * pf + hcs; + + const float sd = get_mesh_signed_distance_at(pos, triangles); + + ZN_ASSERT(grid_index < sdf_grid.size()); + sdf_grid[grid_index] = sd; + + ++grid_index; + } + } + } + + // const uint64_t usec = profiling_clock.restart(); + // println(format("Spent {} usec", usec)); +} + +void generate_mesh_sdf_partitioned(Span sdf_grid, const Vector3i res, const Box3i sub_box, + const Vector3f min_pos, const Vector3f max_pos, const ChunkGrid &chunk_grid) { + ZN_PROFILE_SCOPE(); + ZN_ASSERT(Box3i(Vector3i(), res).contains(sub_box)); + ZN_ASSERT(sdf_grid.size() == Vector3iUtil::get_volume(res)); + + const Vector3f mesh_size = max_pos - min_pos; + const Vector3f cell_size = mesh_size / Vector3f(res.x, res.y, res.z); + + // ProfilingClock profiling_clock; + + Vector3i grid_pos; + const Vector3f hcs(cell_size * 0.5f); + + const Vector3i sub_box_end = sub_box.pos + sub_box.size; + + for (grid_pos.z = sub_box.pos.z; grid_pos.z < sub_box_end.z; ++grid_pos.z) { + for (grid_pos.x = sub_box.pos.x; grid_pos.x < sub_box_end.x; ++grid_pos.x) { + grid_pos.y = sub_box.pos.y; + size_t grid_index = Vector3iUtil::get_zxy_index(grid_pos, res); + + for (; grid_pos.y < sub_box_end.y; ++grid_pos.y) { + const Vector3f pf(float(grid_pos.x) / res.x, float(grid_pos.y) / res.y, float(grid_pos.z) / res.z); + const Vector3f pos = min_pos + mesh_size * pf + hcs; + + const float sd = get_mesh_signed_distance_at(pos, chunk_grid); + + ZN_ASSERT(grid_index < sdf_grid.size()); + sdf_grid[grid_index] = sd; + + ++grid_index; + } + } + } +} + +void generate_mesh_sdf_partitioned(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos, int subdiv) { + // TODO Make this thread-local? + ChunkGrid chunk_grid; + partition_triangles(subdiv, triangles, min_pos, max_pos, chunk_grid); + generate_mesh_sdf_partitioned(sdf_grid, res, Box3i(Vector3i(), res), min_pos, max_pos, chunk_grid); +} + +CheckResult check_sdf( + Span sdf_grid, Vector3i res, Span triangles, Vector3f min_pos, Vector3f max_pos) { + CheckResult result; + result.ok = false; + + ZN_ASSERT_RETURN_V(math::is_valid_size(res), result); + + if (res.x == 0 || res.y == 0 || res.z == 0) { + // Empty or incomparable, but ok + result.ok = true; + return result; + } + + const float tolerance = 1.1f; + const float max_variation = tolerance * get_max_sdf_variation(min_pos, max_pos, res); + + FixedArray dirs; + dirs[0] = Vector3i(-1, 0, 0); + dirs[1] = Vector3i(1, 0, 0); + dirs[2] = Vector3i(0, -1, 0); + dirs[3] = Vector3i(0, 1, 0); + dirs[4] = Vector3i(0, 0, -1); + dirs[5] = Vector3i(0, 0, 1); + + const Vector3f mesh_size = max_pos - min_pos; + const Vector3f cell_size = mesh_size / Vector3f(res.x, res.y, res.z); + const Vector3f hcs(cell_size * 0.5f); + const Evaluator eval{ triangles, res, min_pos, mesh_size, hcs }; + + Vector3i pos; + for (pos.z = 0; pos.z < res.z; ++pos.z) { + for (pos.x = 0; pos.x < res.x; ++pos.x) { + for (pos.y = 0; pos.y < res.y; ++pos.y) { + // + for (unsigned int dir = 0; dir < dirs.size(); ++dir) { + const Vector3i npos = pos + dirs[dir]; + + if (npos.x < 0 || npos.y < 0 || npos.z < 0 || npos.x >= res.x || npos.y >= res.y || + npos.z >= res.z) { + continue; + } + + const unsigned int loc = Vector3iUtil::get_zxy_index(pos, res); + const unsigned int nloc = Vector3iUtil::get_zxy_index(npos, res); + + const float v = sdf_grid[loc]; + const float nv = sdf_grid[nloc]; + + const float variation = Math::abs(v - nv); + + if (variation > max_variation) { + ZN_PRINT_VERBOSE(format("Found variation of {} > {}, {} and {}, at cell {} and {}", variation, + max_variation, v, nv, pos, npos)); + + const Vector3f posf0 = eval.get_pos(pos); + const Vector3f posf1 = eval.get_pos(npos); + + result.ok = false; + + result.cell0.grid_pos = pos; + result.cell0.mesh_pos = posf0; + result.cell0.closest_triangle_index = get_closest_triangle_precalc(posf0, triangles); + + result.cell1.grid_pos = npos; + result.cell1.mesh_pos = posf1; + result.cell1.closest_triangle_index = get_closest_triangle_precalc(posf1, triangles); + + return result; + } + } + } + } + } + + result.ok = true; + return result; +} + +static const float FAR_SD = 9999999.f; + +// Generates accurate SDF only very close to each triangle of the mesh, forming a "shell", and sets every other cells to +// `FAR_SD`. +void generate_mesh_sdf_seeds(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos) { + ZN_PROFILE_SCOPE(); + + // Fill SDF grid with far distances as "infinity", we'll use that to check if we computed it already + sdf_grid.fill(FAR_SD); + + const Vector3f mesh_size = max_pos - min_pos; + const Vector3f cell_size = mesh_size / Vector3f(res.x, res.y, res.z); + const Vector3f hcs(cell_size * 0.5f); + + const Evaluator eval{ triangles, res, min_pos, mesh_size, hcs }; + + for (unsigned int i = 0; i < triangles.size(); ++i) { + const Triangle &t = triangles[i]; + + const Vector3f aabb_min = math::min(t.v1, math::min(t.v2, t.v3)); + const Vector3f aabb_max = math::max(t.v1, math::max(t.v2, t.v3)); + + const Box3i tbox = + Box3i::from_min_max(to_vec3i(math::floor(aabb_min)), to_vec3i(math::ceil(aabb_max))).padded(1); + + // TODO: we could restrict the evaluation to use only the triangle and not all of them? + + tbox.for_each_cell_zxy([&sdf_grid, eval](const Vector3i &grid_pos) { + const size_t i = Vector3iUtil::get_zxy_index(grid_pos, eval.res); + if (sdf_grid[i] != FAR_SD) { + // Already computed + return; + } + ZN_ASSERT(i < sdf_grid.size()); + sdf_grid[i] = eval(grid_pos); + }); + } +} + +// Propagates SDF from a grid that contains a "shell" of accurate SDF values and `FAR_SD` everywhere else. +// The shell must be closed, and must contain all transitions from negative to positive. +void propagate_sdf_across_axis(Span sdf_grid, const Vector3i res, uint8_t axis, bool negative_axis) { + ZN_PROFILE_SCOPE(); + + struct L { + static inline void propagate(Span sdf_grid, unsigned int i, float &prev) { + const float v = sdf_grid[i]; + if (prev != FAR_SD && (v == FAR_SD || ((v < 0.f) == (prev < 0.f)))) { + if (prev < 0.f) { + const float v2 = prev - 1.f; + if (v2 > v) { + // TODO we could affine distances by propagating to neighbor cells as well? + sdf_grid[i] = v; + } + } else { + const float v2 = prev + 1.f; + if (v2 < v) { + sdf_grid[i] = v; + } + } + } + prev = v; + } + }; + + Vector3i pos; + // TODO Each sweep could start from the boundaries of the mesh instead of boundaries of the grid? + // const Vector3i start_min = Vector3i(); + // const Vector3i start_max = res; + + const int z_increment = res.x * res.y; + const int x_increment = res.x; + const int y_increment = 1; + + switch (axis) { + case Vector3i::AXIS_Y: + if (negative_axis) { + // -Y + for (pos.z = res.z - 1; pos.z >= 0; --pos.z) { + for (pos.x = res.x - 1; pos.x >= 0; --pos.x) { + unsigned int i = Vector3iUtil::get_zxy_index(pos, res); + float prev = sdf_grid[i]; + i -= y_increment; + + for (pos.y = res.y - 2; pos.y >= 0; --pos.y, i -= y_increment) { + L::propagate(sdf_grid, i, prev); + } + } + } + } else { + // +Y + for (pos.z = 0; pos.z < res.z; ++pos.z) { + for (pos.x = 0; pos.x < res.x; ++pos.x) { + unsigned int i = Vector3iUtil::get_zxy_index(pos, res); + float prev = sdf_grid[i]; + i += y_increment; + + for (pos.y = 1; pos.y < res.y; ++pos.y, i += y_increment) { + L::propagate(sdf_grid, i, prev); + } + } + } + } + break; + + case Vector3i::AXIS_X: + if (negative_axis) { + // -X + for (pos.z = res.z - 1; pos.z >= 0; --pos.z) { + for (pos.y = res.y - 1; pos.y >= 0; --pos.y) { + unsigned int i = Vector3iUtil::get_zxy_index(pos, res); + float prev = sdf_grid[i]; + i -= x_increment; + + for (pos.x = res.x - 2; pos.x >= 0; --pos.x, i -= x_increment) { + L::propagate(sdf_grid, i, prev); + } + } + } + } else { + // +X + for (pos.z = 0; pos.z < res.z; ++pos.z) { + for (pos.y = 0; pos.y < res.y; ++pos.y) { + unsigned int i = Vector3iUtil::get_zxy_index(pos, res); + float prev = sdf_grid[i]; + i += x_increment; + + for (pos.x = 1; pos.x < res.x; ++pos.x, i += x_increment) { + L::propagate(sdf_grid, i, prev); + } + } + } + } + break; + + case Vector3i::AXIS_Z: + if (negative_axis) { + // -Z + for (pos.x = res.x - 1; pos.x >= 0; --pos.x) { + for (pos.y = res.y - 1; pos.y >= 0; --pos.y) { + unsigned int i = Vector3iUtil::get_zxy_index(pos, res); + float prev = sdf_grid[i]; + i -= z_increment; + + for (pos.z = res.z - 2; pos.z >= 0; --pos.z, i -= z_increment) { + L::propagate(sdf_grid, i, prev); + } + } + } + } else { + // +Z + for (pos.x = 0; pos.x < res.x; ++pos.x) { + for (pos.y = 0; pos.y < res.y; ++pos.y) { + unsigned int i = Vector3iUtil::get_zxy_index(pos, res); + float prev = sdf_grid[i]; + i += z_increment; + + for (pos.z = 1; pos.z < res.z; ++pos.z, i += z_increment) { + L::propagate(sdf_grid, i, prev); + } + } + } + } + break; + default: + ZN_CRASH(); + } +} + +void generate_mesh_sdf_approx_sweep(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos) { + ZN_PROFILE_SCOPE(); + + generate_mesh_sdf_seeds(sdf_grid, res, triangles, min_pos, max_pos); + + propagate_sdf_across_axis(sdf_grid, res, Vector3i::AXIS_X, false); + propagate_sdf_across_axis(sdf_grid, res, Vector3i::AXIS_Y, false); + propagate_sdf_across_axis(sdf_grid, res, Vector3i::AXIS_Z, false); + + propagate_sdf_across_axis(sdf_grid, res, Vector3i::AXIS_X, true); + propagate_sdf_across_axis(sdf_grid, res, Vector3i::AXIS_Y, true); + propagate_sdf_across_axis(sdf_grid, res, Vector3i::AXIS_Z, true); +} + +void generate_mesh_sdf_full(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos) { + generate_mesh_sdf(sdf_grid, res, Box3i(Vector3i(), res), triangles, min_pos, max_pos); +} + +bool prepare_triangles(Span vertices, Span indices, std::vector &triangles, + Vector3f &out_min_pos, Vector3f &out_max_pos) { + ZN_PROFILE_SCOPE(); + + // The mesh can't be closed if it has less than 4 vertices + ZN_ASSERT_RETURN_V(vertices.size() >= 4, false); + + // The mesh can't be closed if it has less than 4 triangles + ZN_ASSERT_RETURN_V(indices.size() >= 12, false); + ZN_ASSERT_RETURN_V(indices.size() % 3 == 0, false); + + triangles.resize(indices.size() / 3); + + for (size_t ti = 0; ti < triangles.size(); ++ti) { + const int ii = ti * 3; + const int i0 = indices[ii]; + const int i1 = indices[ii + 1]; + const int i2 = indices[ii + 2]; + + Triangle &t = triangles[ti]; + t.v1 = to_vec3f(vertices[i0]); + t.v2 = to_vec3f(vertices[i1]); + t.v3 = to_vec3f(vertices[i2]); + + // Hack to make sure all points are distinct + // const Vector3f midp = (t.v1 + t.v2 + t.v3) / 3.f; + // const float shrink_amount = 0.0001f; + // t.v1 = math::lerp(t.v1, midp, shrink_amount); + // t.v2 = math::lerp(t.v2, midp, shrink_amount); + // t.v3 = math::lerp(t.v3, midp, shrink_amount); + } + + Vector3f min_pos = to_vec3f(vertices[0]); + Vector3f max_pos = min_pos; + + for (int i = 0; i < vertices.size(); ++i) { + const Vector3f p = to_vec3f(vertices[i]); + min_pos = math::min(p, min_pos); + max_pos = math::max(p, max_pos); + } + + out_min_pos = min_pos; + out_max_pos = max_pos; + + precalc_triangles(zylann::to_span(triangles)); + + return true; +} + +Vector3i auto_compute_grid_resolution(const Vector3f box_size, int cell_count) { + const float cs = math::max(box_size.x, math::max(box_size.y, box_size.z)) / float(cell_count); + return Vector3i(box_size.x / cs, box_size.y / cs, box_size.z / cs); +} + +// Called from within the thread pool +void GenMeshSDFSubBoxTask::run(ThreadedTaskContext ctx) { + ZN_PROFILE_SCOPE(); + ZN_ASSERT(shared_data != nullptr); + + VoxelBufferInternal &buffer = shared_data->buffer; + const VoxelBufferInternal::ChannelId channel = VoxelBufferInternal::CHANNEL_SDF; + //ZN_ASSERT(!buffer.get_channel_compression(channel) == VoxelBufferInternal::COMPRESSION_NONE); + Span sdf_grid; + ZN_ASSERT(buffer.get_channel_data(channel, sdf_grid)); + + if (shared_data->use_chunk_grid) { + generate_mesh_sdf_partitioned( + sdf_grid, buffer.get_size(), box, shared_data->min_pos, shared_data->max_pos, shared_data->chunk_grid); + } else { + generate_mesh_sdf(sdf_grid, buffer.get_size(), box, to_span(shared_data->triangles), shared_data->min_pos, + shared_data->max_pos); + } + + if (shared_data->pending_jobs.fetch_sub(1, std::memory_order_acq_rel) == 1) { + if (shared_data->boundary_sign_fix) { + fix_sdf_sign_from_boundary(sdf_grid, buffer.get_size(), shared_data->min_pos, shared_data->max_pos); + } + // That was the last job + on_complete(); + } +} + +} // namespace zylann::voxel::mesh_sdf diff --git a/edition/mesh_sdf.h b/edition/mesh_sdf.h new file mode 100644 index 00000000..7b890c3c --- /dev/null +++ b/edition/mesh_sdf.h @@ -0,0 +1,142 @@ +#ifndef VOXEL_MESH_SDF_H +#define VOXEL_MESH_SDF_H + +#include "../storage/voxel_buffer_internal.h" +#include "../util/math/vector3f.h" +#include "../util/math/vector3i.h" +#include "../util/span.h" +#include "../util/tasks/threaded_task.h" + +#include +#include +#include + +namespace zylann::voxel::mesh_sdf { + +// Utilities to generate a signed distance field from a 3D triangle mesh. + +struct Triangle { + // Vertices to provide from the mesh. + Vector3f v1; + Vector3f v2; + Vector3f v3; + + // Values precomputed with `prepare_triangles()`. + Vector3f v21; + Vector3f v32; + Vector3f v13; + Vector3f nor; + Vector3f v21_cross_nor; + Vector3f v32_cross_nor; + Vector3f v13_cross_nor; + float inv_v21_length_squared; + float inv_v32_length_squared; + float inv_v13_length_squared; + float inv_nor_length_squared; +}; + +struct Chunk { + Vector3i pos; + std::vector near_chunks; + std::vector triangles; +}; + +struct ChunkGrid { + std::vector chunks; + Vector3i size; + Vector3f min_pos; + float chunk_size; +}; + +class GenMeshSDFSubBoxTask : public IThreadedTask { +public: + struct SharedData { + std::vector triangles; + std::atomic_int pending_jobs; + VoxelBufferInternal buffer; + Vector3f min_pos; + Vector3f max_pos; + ChunkGrid chunk_grid; + bool use_chunk_grid = false; + bool boundary_sign_fix = false; + }; + + std::shared_ptr shared_data; + Box3i box; + + void run(ThreadedTaskContext ctx) override; + + // Called when `pending_jobs` reaches zero. + virtual void on_complete() {} + + void apply_result() override {} +}; + +// Computes a representation of the mesh that's more optimal to compute distance to triangles. +bool prepare_triangles(Span vertices, Span indices, std::vector &triangles, + Vector3f &out_min_pos, Vector3f &out_max_pos); + +// Partitions triangles of the mesh such that we can reduce the number of triangles to check when evaluating the SDF. +// Space is subdivided in a grid of chunks. Triangles overlapping chunks are listed, then each chunk finds which other +// non-empty chunks are close to it. The amount of subvidisions should be carefully chosen: too low will cause less +// triangles to be skipped, too high will make partitionning slower. +void partition_triangles( + int subdiv, Span triangles, Vector3f min_pos, Vector3f max_pos, ChunkGrid &chunk_grid); + +// A naive method to get a sampled SDF from a mesh, by checking every triangle at every cell. It's accurate, but much +// slower than other techniques, but could be used as a CPU-based alternative, for less +// realtime-intensive tasks. The mesh must be closed, otherwise the SDF will contain errors. +void generate_mesh_sdf_full(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos); + +// Compute the SDF faster by partitionning triangles, while retaining the same accuracy as if all triangles +// were checked. With Suzanne mesh subdivided once with 3900 triangles and `subdiv = 32`, it's about 8 times fasterthan +// checking every triangle on every cell. +void generate_mesh_sdf_partitioned(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos, int subdiv); + +// Generates an approximation. +// Subdivides the grid into nodes spanning 4*4*4 cells each. +// If a node's corner distances are close to the surface, the SDF is fully evaluated. Otherwise, it is interpolated. +// Tests with Suzanne show it is 2 to 3 times faster than the basic naive method, with only minor quality decrease. +// It's still quite slow though. +void generate_mesh_sdf_approx_interp(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos); + +// TODO Untested +// Generates an approximation. +// Calculates a thin hull of accurate SDF values, then propagates it across each axis using manhattan distance. +void generate_mesh_sdf_approx_sweep(Span sdf_grid, const Vector3i res, Span triangles, + const Vector3f min_pos, const Vector3f max_pos); + +Vector3i auto_compute_grid_resolution(const Vector3f box_size, int cell_count); + +struct CheckResult { + bool ok; + struct BadCell { + Vector3i grid_pos; + Vector3f mesh_pos; + unsigned int closest_triangle_index; + }; + BadCell cell0; + BadCell cell1; +}; + +// Checks if SDF variations are legit. The difference between two neighboring cells cannot be higher than the distance +// between those two cells. This is intented at proper SDF, not approximation or scaled ones. +CheckResult check_sdf( + Span sdf_grid, Vector3i res, Span triangles, Vector3f min_pos, Vector3f max_pos); + +// The current method provides imperfect signs. Due to ambiguities, sometimes patches of cells get the wrong sign. +// This function attemps to correct these. +// Assumes the sign on the edge of the box is positive and use a floodfill. +// If we start from the rough SDF we had, we could do a floodfill that considers unexpected sign change as fillable, +// while an expected sign change would properly stop the fill. +// However, this workaround won't fix signs inside the volume. +// I thought of using this with an completely unsigned distance field instead, however I'm not sure if it's possible to +// accurately tell when when the sign is supposed to flip (i.e when we cross the surface). +void fix_sdf_sign_from_boundary(Span sdf_grid, Vector3i res, Vector3f min_pos, Vector3f max_pos); + +} // namespace zylann::voxel::mesh_sdf + +#endif // VOXEL_MESH_SDF_H diff --git a/edition/voxel_mesh_sdf_gd.cpp b/edition/voxel_mesh_sdf_gd.cpp new file mode 100644 index 00000000..778a9aea --- /dev/null +++ b/edition/voxel_mesh_sdf_gd.cpp @@ -0,0 +1,383 @@ +#include "voxel_mesh_sdf_gd.h" +#include "../server/voxel_server.h" +#include "../server/voxel_server_updater.h" +#include "../storage/voxel_buffer_gd.h" +#include "../util/godot/funcs.h" +#include "../util/math/color.h" +#include "../util/math/conv.h" +#include "../util/profiling.h" +#include "../util/string_funcs.h" +#include "mesh_sdf.h" + +#include + +namespace zylann::voxel { + +static bool prepare_triangles( + Mesh &mesh, std::vector &triangles, Vector3f &out_min_pos, Vector3f &out_max_pos) { + ZN_PROFILE_SCOPE(); + ERR_FAIL_COND_V(mesh.get_surface_count() == 0, false); + Array surface; + { + ZN_PROFILE_SCOPE_NAMED("Get surface from Godot") + surface = mesh.surface_get_arrays(0); + } + PackedVector3Array positions = surface[Mesh::ARRAY_VERTEX]; + PackedInt32Array indices = surface[Mesh::ARRAY_INDEX]; + ERR_FAIL_COND_V( + !mesh_sdf::prepare_triangles(to_span(positions), to_span(indices), triangles, out_min_pos, out_max_pos), + false); + return true; +} + +bool VoxelMeshSDF::is_baked() const { + return _voxel_buffer.is_valid(); +} + +bool VoxelMeshSDF::is_baking() const { + return _is_baking; +} + +int VoxelMeshSDF::get_cell_count() const { + return _cell_count; +} + +void VoxelMeshSDF::set_cell_count(int cc) { + _cell_count = math::clamp(cc, 2, 256); +} + +float VoxelMeshSDF::get_margin_ratio() const { + return _margin_ratio; +} + +void VoxelMeshSDF::set_margin_ratio(float mr) { + _margin_ratio = math::clamp(mr, 0.f, 1.f); +} + +VoxelMeshSDF::BakeMode VoxelMeshSDF::get_bake_mode() const { + return _bake_mode; +} + +void VoxelMeshSDF::set_bake_mode(BakeMode mode) { + ERR_FAIL_INDEX(mode, BAKE_MODE_COUNT); + _bake_mode = mode; +} + +int VoxelMeshSDF::get_partition_subdiv() const { + return _partition_subdiv; +} + +void VoxelMeshSDF::set_partition_subdiv(int subdiv) { + ERR_FAIL_COND(subdiv < 2 || subdiv > 255); + _partition_subdiv = subdiv; +} + +void VoxelMeshSDF::set_boundary_sign_fix_enabled(bool enable) { + _boundary_sign_fix = enable; +} + +bool VoxelMeshSDF::is_boundary_sign_fix_enabled() const { + return _boundary_sign_fix; +} + +void VoxelMeshSDF::bake(Ref mesh) { + ZN_PROFILE_SCOPE(); + + ERR_FAIL_COND(mesh.is_null()); + + std::vector triangles; + Vector3f min_pos; + Vector3f max_pos; + ERR_FAIL_COND(!prepare_triangles(**mesh, triangles, min_pos, max_pos)); + + const Vector3f mesh_size = max_pos - min_pos; + + const Vector3f box_min_pos = min_pos - mesh_size * _margin_ratio; + const Vector3f box_max_pos = max_pos + mesh_size * _margin_ratio; + const Vector3f box_size = box_max_pos - box_min_pos; + + const Vector3i res = mesh_sdf::auto_compute_grid_resolution(box_size, _cell_count); + const uint64_t volume = Vector3iUtil::get_volume(res); + const VoxelBufferInternal::ChannelId channel = VoxelBufferInternal::CHANNEL_SDF; + Ref vbgd; + vbgd.instantiate(); + VoxelBufferInternal &buffer = vbgd->get_buffer(); + buffer.set_channel_depth(channel, VoxelBufferInternal::DEPTH_32_BIT); + buffer.create(res); + buffer.decompress_channel(channel); + Span sdf_grid; + ERR_FAIL_COND(!buffer.get_channel_data(channel, sdf_grid)); + + switch (_bake_mode) { + case BAKE_MODE_ACCURATE_NAIVE: + generate_mesh_sdf_full(sdf_grid, res, to_span(triangles), box_min_pos, box_max_pos); + break; + case BAKE_MODE_ACCURATE_PARTITIONED: + generate_mesh_sdf_partitioned( + sdf_grid, res, to_span(triangles), box_min_pos, box_max_pos, _partition_subdiv); + break; + case BAKE_MODE_APPROX_INTERP: + generate_mesh_sdf_approx_interp(sdf_grid, res, to_span(triangles), box_min_pos, box_max_pos); + break; + case BAKE_MODE_APPROX_SWEEP: + generate_mesh_sdf_approx_sweep(sdf_grid, res, to_span(triangles), box_min_pos, box_max_pos); + break; + default: + ZN_CRASH(); + } + + if (_boundary_sign_fix) { + mesh_sdf::fix_sdf_sign_from_boundary(sdf_grid, res, min_pos, max_pos); + } + + _voxel_buffer = vbgd; + _min_pos = box_min_pos; + _max_pos = box_max_pos; +} + +void VoxelMeshSDF::bake_async(Ref mesh, SceneTree *scene_tree) { + ZN_ASSERT_RETURN(scene_tree != nullptr); + VoxelServerUpdater::ensure_existence(scene_tree); + + //ZN_ASSERT_RETURN_MSG(!_is_baking, "Already baking"); + + struct L { + static void notify_on_complete(VoxelMeshSDF &obj, mesh_sdf::GenMeshSDFSubBoxTask::SharedData &shared_data) { + Ref vbgd; + vbgd.instantiate(); + shared_data.buffer.move_to(vbgd->get_buffer()); + obj.call_deferred( + "_on_generate_async_completed", vbgd, to_godot(shared_data.min_pos), to_godot(shared_data.max_pos)); + } + }; + + class GenMeshSDFSubBoxTaskGD : public mesh_sdf::GenMeshSDFSubBoxTask { + public: + Ref obj_to_notify; + + void on_complete() override { + ZN_ASSERT(obj_to_notify.is_valid()); + L::notify_on_complete(**obj_to_notify, *shared_data); + } + }; + + class GenMeshSDFFirstPassTask : public IThreadedTask { + public: + float margin_ratio; + int cell_count; + BakeMode bake_mode; + uint8_t partition_subdiv; + bool boundary_sign_fix; + Array surface; + Ref obj_to_notify; + + void run(ThreadedTaskContext ctx) override { + ZN_PROFILE_SCOPE(); + ZN_ASSERT(obj_to_notify.is_valid()); + + std::shared_ptr shared_data = + make_shared_instance(); + Vector3f min_pos; + Vector3f max_pos; + + PackedVector3Array positions = surface[Mesh::ARRAY_VERTEX]; + PackedInt32Array indices = surface[Mesh::ARRAY_INDEX]; + if (!mesh_sdf::prepare_triangles( + to_span(positions), to_span(indices), shared_data->triangles, min_pos, max_pos)) { + ZN_PRINT_ERROR("Failed preparing triangles in threaded task"); + report_error(); + return; + } + + const Vector3f mesh_size = max_pos - min_pos; + + const Vector3f box_min_pos = min_pos - mesh_size * margin_ratio; + const Vector3f box_max_pos = max_pos + mesh_size * margin_ratio; + const Vector3f box_size = box_max_pos - box_min_pos; + + const Vector3i res = mesh_sdf::auto_compute_grid_resolution(box_size, cell_count); + const VoxelBufferInternal::ChannelId channel = VoxelBufferInternal::CHANNEL_SDF; + shared_data->buffer.set_channel_depth(channel, VoxelBufferInternal::DEPTH_32_BIT); + shared_data->buffer.create(res); + shared_data->buffer.decompress_channel(channel); + + shared_data->min_pos = box_min_pos; + shared_data->max_pos = box_max_pos; + + switch (bake_mode) { + case BAKE_MODE_ACCURATE_NAIVE: + case BAKE_MODE_ACCURATE_PARTITIONED: { + // These two approaches are better parallelized + + const bool partitioned = bake_mode == BAKE_MODE_ACCURATE_PARTITIONED; + if (partitioned) { + mesh_sdf::partition_triangles(partition_subdiv, to_span(shared_data->triangles), + shared_data->min_pos, shared_data->max_pos, shared_data->chunk_grid); + } + shared_data->use_chunk_grid = partitioned; + + shared_data->boundary_sign_fix = boundary_sign_fix; + + // Spawn a parallel task for every Z slice of the grid. + // Indexing is ZXY so each thread accesses a contiguous part of memory. + shared_data->pending_jobs = res.z; + + for (int z = 0; z < res.z; ++z) { + GenMeshSDFSubBoxTaskGD *task = ZN_NEW(GenMeshSDFSubBoxTaskGD); + task->shared_data = shared_data; + task->box = Box3i(Vector3i(0, 0, z), Vector3i(res.x, res.y, 1)); + task->obj_to_notify = obj_to_notify; + + VoxelServer::get_singleton().push_async_task(task); + } + } break; + + case BAKE_MODE_APPROX_INTERP: + case BAKE_MODE_APPROX_SWEEP: { + VoxelBufferInternal &buffer = shared_data->buffer; + const VoxelBufferInternal::ChannelId channel = VoxelBufferInternal::CHANNEL_SDF; + Span sdf_grid; + ZN_ASSERT(buffer.get_channel_data(channel, sdf_grid)); + + if (bake_mode == BAKE_MODE_APPROX_INTERP) { + generate_mesh_sdf_approx_interp( + sdf_grid, res, to_span(shared_data->triangles), box_min_pos, box_max_pos); + } else { + generate_mesh_sdf_approx_sweep( + sdf_grid, res, to_span(shared_data->triangles), box_min_pos, box_max_pos); + } + + if (boundary_sign_fix) { + mesh_sdf::fix_sdf_sign_from_boundary(sdf_grid, res, box_min_pos, box_max_pos); + } + L::notify_on_complete(**obj_to_notify, *shared_data); + } break; + + default: + ZN_PRINT_ERROR(format("Invalid bake mode {}", bake_mode)); + report_error(); + break; + } + } + + void apply_result() override {} + + private: + void report_error() { + obj_to_notify->call_deferred("_on_generate_async_completed", Ref(), Vector3(), Vector3()); + } + }; + + ERR_FAIL_COND(mesh.is_null()); + + ERR_FAIL_COND(mesh->get_surface_count() == 0); + Array surface; + { + ZN_PROFILE_SCOPE_NAMED("Get surface from Godot") + surface = mesh->surface_get_arrays(0); + } + + _is_baking = true; + + GenMeshSDFFirstPassTask *task = ZN_NEW(GenMeshSDFFirstPassTask); + task->cell_count = _cell_count; + task->margin_ratio = _margin_ratio; + task->bake_mode = _bake_mode; + task->partition_subdiv = _partition_subdiv; + task->surface = surface; + task->obj_to_notify.reference_ptr(this); + task->boundary_sign_fix = _boundary_sign_fix; + VoxelServer::get_singleton().push_async_task(task); +} + +void VoxelMeshSDF::_on_bake_async_completed(Ref buffer, Vector3 min_pos, Vector3 max_pos) { + _is_baking = false; + + // This can mean an error occurred during one of the tasks + ZN_ASSERT_RETURN(buffer.is_valid()); + + _voxel_buffer = buffer; + _min_pos = to_vec3f(min_pos); + _max_pos = to_vec3f(max_pos); + emit_signal("baked"); +} + +Ref VoxelMeshSDF::get_voxel_buffer() const { + return _voxel_buffer; +} + +AABB VoxelMeshSDF::get_aabb() const { + return AABB(to_godot(_min_pos), to_godot(_max_pos - _min_pos)); +} + +Array VoxelMeshSDF::debug_check_sdf(Ref mesh) { + Array result; + + ZN_ASSERT_RETURN_V(is_baked(), result); + ZN_ASSERT(_voxel_buffer.is_valid()); + const VoxelBufferInternal &buffer = _voxel_buffer->get_buffer(); + Span sdf_grid; + ZN_ASSERT_RETURN_V(buffer.get_channel_data(VoxelBufferInternal::CHANNEL_SDF, sdf_grid), result); + + ZN_ASSERT_RETURN_V(mesh.is_valid(), result); + std::vector triangles; + Vector3f min_pos; + Vector3f max_pos; + ZN_ASSERT_RETURN_V(prepare_triangles(**mesh, triangles, min_pos, max_pos), result); + + const mesh_sdf::CheckResult cr = + mesh_sdf::check_sdf(sdf_grid, buffer.get_size(), to_span(triangles), _min_pos, _max_pos); + + if (cr.ok) { + return result; + } + + const mesh_sdf::Triangle &ct0 = triangles[cr.cell0.closest_triangle_index]; + const mesh_sdf::Triangle &ct1 = triangles[cr.cell1.closest_triangle_index]; + + result.resize(8); + result[0] = to_godot(cr.cell0.mesh_pos); + result[1] = to_godot(ct0.v1); + result[2] = to_godot(ct0.v2); + result[3] = to_godot(ct0.v3); + result[4] = to_godot(cr.cell1.mesh_pos); + result[5] = to_godot(ct1.v1); + result[6] = to_godot(ct1.v2); + result[7] = to_godot(ct1.v3); + return result; +} + +void VoxelMeshSDF::_bind_methods() { + ClassDB::bind_method(D_METHOD("bake", "mesh"), &VoxelMeshSDF::bake); + ClassDB::bind_method(D_METHOD("bake_async", "mesh", "scene_tree"), &VoxelMeshSDF::bake_async); + + ClassDB::bind_method(D_METHOD("get_cell_count"), &VoxelMeshSDF::get_cell_count); + ClassDB::bind_method(D_METHOD("set_cell_count", "cell_count"), &VoxelMeshSDF::set_cell_count); + + ClassDB::bind_method(D_METHOD("get_margin_ratio"), &VoxelMeshSDF::get_margin_ratio); + ClassDB::bind_method(D_METHOD("set_margin_ratio", "ratio"), &VoxelMeshSDF::set_margin_ratio); + + ClassDB::bind_method(D_METHOD("get_bake_mode"), &VoxelMeshSDF::get_bake_mode); + ClassDB::bind_method(D_METHOD("set_bake_mode", "mode"), &VoxelMeshSDF::set_bake_mode); + + ClassDB::bind_method(D_METHOD("get_partition_subdiv"), &VoxelMeshSDF::get_partition_subdiv); + ClassDB::bind_method(D_METHOD("set_partition_subdiv", "subdiv"), &VoxelMeshSDF::set_partition_subdiv); + + ClassDB::bind_method(D_METHOD("get_voxel_buffer"), &VoxelMeshSDF::get_voxel_buffer); + ClassDB::bind_method(D_METHOD("debug_check_sdf", "mesh"), &VoxelMeshSDF::debug_check_sdf); + + // Internal + ClassDB::bind_method(D_METHOD("_on_bake_async_completed", "buffer", "min_pos", "max_pos"), + &VoxelMeshSDF::_on_bake_async_completed); + + ADD_SIGNAL(MethodInfo("baked")); + + // These modes are mostly for experimentation, I'm not sure if they will remain + BIND_ENUM_CONSTANT(BAKE_MODE_ACCURATE_NAIVE); + BIND_ENUM_CONSTANT(BAKE_MODE_ACCURATE_PARTITIONED); + BIND_ENUM_CONSTANT(BAKE_MODE_APPROX_INTERP); + BIND_ENUM_CONSTANT(BAKE_MODE_APPROX_SWEEP); + BIND_ENUM_CONSTANT(BAKE_MODE_COUNT); +} + +} // namespace zylann::voxel diff --git a/edition/voxel_mesh_sdf_gd.h b/edition/voxel_mesh_sdf_gd.h new file mode 100644 index 00000000..ad66b33c --- /dev/null +++ b/edition/voxel_mesh_sdf_gd.h @@ -0,0 +1,87 @@ +#ifndef VOXEL_MESH_SDF_GD_H +#define VOXEL_MESH_SDF_GD_H + +#include "../storage/voxel_buffer_gd.h" +#include "../util/math/vector3f.h" + +#include + +class Mesh; +class SceneTree; + +namespace zylann::voxel { + +// Contains the baked signed distance field of a mesh, which can be used to sculpt terrain. +// TODO Make it a resource so we can pre-build, save and load the baked data more easily +class VoxelMeshSDF : public RefCounted { + GDCLASS(VoxelMeshSDF, RefCounted) +public: + enum BakeMode { // + BAKE_MODE_ACCURATE_NAIVE, + BAKE_MODE_ACCURATE_PARTITIONED, + BAKE_MODE_APPROX_INTERP, + BAKE_MODE_APPROX_SWEEP, + BAKE_MODE_COUNT + }; + + // The data cannot be used until baked + bool is_baked() const; + bool is_baking() const; + + int get_cell_count() const; + void set_cell_count(int cc); + + float get_margin_ratio() const; + void set_margin_ratio(float mr); + + BakeMode get_bake_mode() const; + void set_bake_mode(BakeMode mode); + + int get_partition_subdiv() const; + void set_partition_subdiv(int subdiv); + + void set_boundary_sign_fix_enabled(bool enable); + bool is_boundary_sign_fix_enabled() const; + + // Note, the mesh is not referenced because we don't want it to be a dependency. + // An SDF should be usable even without the original mesh. + // The mesh is only necessary when baking. + + void bake(Ref mesh); + + // Bakes the SDF asynchronously using threads of the job system. + // TODO A reference to the SceneTree should not be necessary! + // It is currently needed to ensure `VoxelServerUpdater` gets created so it can tick the task system... + void bake_async(Ref mesh, SceneTree *scene_tree); + + Ref get_voxel_buffer() const; + AABB get_aabb() const; + + Array debug_check_sdf(Ref mesh); + +private: + void _on_bake_async_completed(Ref buffer, Vector3 min_pos, Vector3 max_pos); + + static void _bind_methods(); + + // Data + Ref _voxel_buffer; + Vector3f _min_pos; + Vector3f _max_pos; + + // States + bool _is_baking = false; + + // Baking options + int _cell_count = 32; + float _margin_ratio = 0.25; + BakeMode _bake_mode = BAKE_MODE_ACCURATE_PARTITIONED; + uint8_t _partition_subdiv = 32; + bool _boundary_sign_fix = true; +}; + +} // namespace zylann::voxel + +VARIANT_ENUM_CAST(zylann::voxel::VoxelMeshSDF::BakeMode); + +#endif // VOXEL_MESH_SDF_GD_H diff --git a/edition/voxel_tool_lod_terrain.cpp b/edition/voxel_tool_lod_terrain.cpp index 6a66f926..568567f1 100644 --- a/edition/voxel_tool_lod_terrain.cpp +++ b/edition/voxel_tool_lod_terrain.cpp @@ -9,6 +9,7 @@ #include "../util/tasks/async_dependency_tracker.h" #include "../util/voxel_raycast.h" #include "funcs.h" +#include "voxel_mesh_sdf_gd.h" #include #include @@ -722,6 +723,100 @@ Array VoxelToolLodTerrain::separate_floating_chunks(AABB world_box, Node *parent *this, int_world_box, parent_node, _terrain->get_global_transform(), mesher, materials); } +// Combines a precalculated SDF with the terrain at a specific position, rotation and scale. +// +// `transform` is where the buffer should be applied on the terrain. +// +// `isolevel` alters the shape of the SDF: positive "puffs" it, negative "erodes" it. This is a applied after +// `sdf_scale`. +// +// `sdf_scale` scales SDF values (it doesnt make the shape bigger or smaller). Usually defaults to 1 but may be lower if +// artifacts show up due to scaling used in terrain SDF. +// +void VoxelToolLodTerrain::stamp_sdf( + Ref mesh_sdf, Transform3D transform, float isolevel, float sdf_scale) { + // TODO Asynchronous version + ZN_PROFILE_SCOPE(); + + ERR_FAIL_COND(_terrain == nullptr); + ERR_FAIL_COND(mesh_sdf.is_null()); + ERR_FAIL_COND(mesh_sdf->is_baked()); + Ref buffer_ref = mesh_sdf->get_voxel_buffer(); + ERR_FAIL_COND(buffer_ref.is_null()); + const VoxelBufferInternal &buffer = buffer_ref->get_buffer(); + const VoxelBufferInternal::ChannelId channel = VoxelBufferInternal::CHANNEL_SDF; + ERR_FAIL_COND(buffer.get_channel_compression(channel) == VoxelBufferInternal::COMPRESSION_UNIFORM); + ERR_FAIL_COND(buffer.get_channel_depth(channel) != VoxelBufferInternal::DEPTH_32_BIT); + + const Transform3D &box_to_world = transform; + const AABB local_aabb = mesh_sdf->get_aabb(); + + // Note, transform is local to the terrain + const AABB aabb = box_to_world.xform(aabb); + const Box3i voxel_box = Box3i::from_min_max(aabb.position.floor(), (aabb.position + aabb.size).ceil()); + + // TODO Sometimes it will fail near unloaded blocks, even though the transformed box does not intersect them. + // This could be avoided with a box/transformed-box intersection algorithm. Might investigate if the use case + // occurs. It won't happen with full load mode. This also affects other shapes. + if (!is_area_editable(voxel_box)) { + ZN_PRINT_VERBOSE("Area not editable"); + return; + } + + std::shared_ptr data = _terrain->get_storage(); + ERR_FAIL_COND(data == nullptr); + VoxelDataLodMap::Lod &data_lod = data->lods[0]; + + if (_terrain->is_full_load_mode_enabled()) { + preload_box(*data, voxel_box, _terrain->get_generator().ptr()); + } + + // TODO Maybe more efficient to "rasterize" the box? We're going to iterate voxels the box doesnt intersect + // TODO Maybe we should scale SDF values based on the scale of the transform too + + struct SdfBufferShape { + Span buffer; + Vector3i buffer_size; + Transform3D world_to_buffer; + float isolevel; + float sdf_scale; + + inline real_t operator()(const Vector3 &wpos) const { + // Transform terrain-space position to buffer-space + const Vector3f lpos = to_vec3f(world_to_buffer.xform(wpos)); + if (lpos.x < 0 || lpos.y < 0 || lpos.z < 0 || lpos.x >= buffer_size.x || lpos.y >= buffer_size.y || + lpos.z >= buffer_size.z) { + // Outside the buffer + return 100; + } + return interpolate_trilinear(buffer, buffer_size, lpos) * sdf_scale - isolevel; + } + }; + + const Transform3D buffer_to_box = + Transform3D(Basis().scaled(Vector3(local_aabb.size / buffer.get_size())), local_aabb.position); + const Transform3D buffer_to_world = box_to_world * buffer_to_box; + + // TODO Support other depths, format should be accessible from the volume + ops::SdfOperation16bit op; + op.shape.world_to_buffer = buffer_to_world.affine_inverse(); + op.shape.buffer_size = buffer.get_size(); + op.shape.isolevel = isolevel; + op.shape.sdf_scale = sdf_scale; + // Note, the passed buffer must not be shared with another thread. + //buffer.decompress_channel(channel); + ZN_ASSERT_RETURN(buffer.get_channel_data(channel, op.shape.buffer)); + + VoxelDataGrid grid; + { + RWLockRead rlock(data_lod.map_lock); + grid.reference_area(data_lod.map, voxel_box); + grid.write_box(voxel_box, VoxelBufferInternal::CHANNEL_SDF, op); + } + + _post_edit(voxel_box); +} + void VoxelToolLodTerrain::_bind_methods() { ClassDB::bind_method(D_METHOD("set_raycast_binary_search_iterations", "iterations"), &VoxelToolLodTerrain::set_raycast_binary_search_iterations); @@ -732,6 +827,8 @@ void VoxelToolLodTerrain::_bind_methods() { ClassDB::bind_method( D_METHOD("separate_floating_chunks", "box", "parent_node"), &VoxelToolLodTerrain::separate_floating_chunks); ClassDB::bind_method(D_METHOD("do_sphere_async", "center", "radius"), &VoxelToolLodTerrain::do_sphere_async); + ClassDB::bind_method( + D_METHOD("stamp_sdf", "mesh_sdf", "transform", "isolevel", "sdf_scale"), &VoxelToolLodTerrain::stamp_sdf); } } // namespace zylann::voxel diff --git a/edition/voxel_tool_lod_terrain.h b/edition/voxel_tool_lod_terrain.h index 762256e0..e93cb20c 100644 --- a/edition/voxel_tool_lod_terrain.h +++ b/edition/voxel_tool_lod_terrain.h @@ -9,6 +9,7 @@ namespace zylann::voxel { class VoxelLodTerrain; class VoxelDataMap; +class VoxelMeshSDF; class VoxelToolLodTerrain : public VoxelTool { GDCLASS(VoxelToolLodTerrain, VoxelTool) @@ -31,6 +32,7 @@ public: float get_voxel_f_interpolated(Vector3 position) const; Array separate_floating_chunks(AABB world_box, Node *parent_node); + void stamp_sdf(Ref mesh_sdf, Transform3D transform, float isolevel, float sdf_scale); protected: uint64_t _get_voxel(Vector3i pos) const override; diff --git a/register_types.cpp b/register_types.cpp index faa4c546..5334c93f 100644 --- a/register_types.cpp +++ b/register_types.cpp @@ -1,5 +1,6 @@ #include "register_types.h" #include "constants/voxel_string_names.h" +#include "edition/voxel_mesh_sdf_gd.h" #include "edition/voxel_tool.h" #include "edition/voxel_tool_buffer.h" #include "edition/voxel_tool_lod_terrain.h" @@ -141,6 +142,7 @@ void register_voxel_types() { #ifdef VOXEL_ENABLE_FAST_NOISE_2 ClassDB::register_class(); #endif + ClassDB::register_class(); // Meshers ClassDB::register_abstract_class(); diff --git a/storage/voxel_buffer_internal.cpp b/storage/voxel_buffer_internal.cpp index 7c7ba89d..b6510a80 100644 --- a/storage/voxel_buffer_internal.cpp +++ b/storage/voxel_buffer_internal.cpp @@ -589,6 +589,7 @@ void VoxelBufferInternal::duplicate_to(VoxelBufferInternal &dst, bool include_me void VoxelBufferInternal::move_to(VoxelBufferInternal &dst) { if (this == &dst) { + ZN_PRINT_VERBOSE("Moving VoxelBufferInternal to itself?"); return; } diff --git a/storage/voxel_buffer_internal.h b/storage/voxel_buffer_internal.h index dc09e944..a45d3a68 100644 --- a/storage/voxel_buffer_internal.h +++ b/storage/voxel_buffer_internal.h @@ -399,9 +399,16 @@ public: return Vector3iUtil::get_volume(_size); } - // TODO Have a template version based on channel depth bool get_channel_raw(unsigned int channel_index, Span &slice) const; + template + bool get_channel_data(unsigned int channel_index, Span &dst) const { + Span dst8; + ZN_ASSERT_RETURN_V(get_channel_raw(channel_index, dst8), false); + dst = dst8.reinterpret_cast_to(); + return true; + } + void downscale_to(VoxelBufferInternal &dst, Vector3i src_min, Vector3i src_max, Vector3i dst_min) const; bool equals(const VoxelBufferInternal &p_other) const; diff --git a/util/math/color8.h b/util/math/color8.h index a88f86ff..bfb373b3 100644 --- a/util/math/color8.h +++ b/util/math/color8.h @@ -5,6 +5,7 @@ namespace zylann { +// Color with 8-bit components. Lighter to store than its floating-point counterpart. struct Color8 { union { struct { diff --git a/util/math/conv.h b/util/math/conv.h index d5746b80..3ffcd92d 100644 --- a/util/math/conv.h +++ b/util/math/conv.h @@ -8,6 +8,8 @@ namespace zylann { +// Explicit conversion methods. Not in respective files because it would cause circular dependencies. + inline Vector3i to_vec3i(Vector3f v) { return Vector3i(v.x, v.y, v.z); } diff --git a/util/math/funcs.h b/util/math/funcs.h index b28dea9f..375efb4c 100644 --- a/util/math/funcs.h +++ b/util/math/funcs.h @@ -237,6 +237,10 @@ inline void sort(T &a, T &b, T &c, T &d) { sort(b, c); } +inline float sign(float x) { + return x < 0.f ? -1.f : 1.f; +} + } // namespace zylann::math #endif // VOXEL_MATH_FUNCS_H diff --git a/util/math/vector3f.h b/util/math/vector3f.h index 23c519ef..2097d446 100644 --- a/util/math/vector3f.h +++ b/util/math/vector3f.h @@ -2,7 +2,7 @@ #define ZYLANN_VECTOR3F_H #include "../errors.h" -#include +#include "funcs.h" namespace zylann { @@ -26,6 +26,7 @@ struct Vector3f { }; Vector3f() : x(0), y(0), z(0) {} + explicit Vector3f(float p_v) : x(p_v), y(p_v), z(p_v) {} Vector3f(float p_x, float p_y, float p_z) : x(p_x), y(p_y), z(p_z) {} inline float length_squared() const { @@ -212,6 +213,25 @@ inline T interpolate_trilinear(const T v000, const T v100, const T v101, const T return v; } + +inline Vector3f min(const Vector3f a, const Vector3f b) { + return Vector3f(min(a.x, b.x), min(a.y, b.y), min(a.z, b.z)); +} + +inline Vector3f max(const Vector3f a, const Vector3f b) { + return Vector3f(max(a.x, b.x), max(a.y, b.y), max(a.z, b.z)); +} + +inline Vector3f floor(const Vector3f a) { + return Vector3f(Math::floor(a.x), Math::floor(a.y), Math::floor(a.z)); +} + +inline Vector3f ceil(const Vector3f a) { + return Vector3f(Math::ceil(a.x), Math::ceil(a.y), Math::ceil(a.z)); +} + +inline Vector3f lerp(const Vector3f a, const Vector3f b, const float t) { + return Vector3f(Math::lerp(a.x, b.x, t), Math::lerp(a.y, b.y, t), Math::lerp(a.z, b.z, t)); } } // namespace math