Added VoxelMeshSDF and VoxelToolLodTerrain.stamp_sdf()

config.py

@@ -69,6 +69,7 @@ def get_doc_classes():
         "VoxelBlockSerializer",
         "VoxelVoxLoader",
         "VoxelDataBlockEnterInfo",
+        "VoxelMeshSDF",
 
         "ZN_FastNoiseLite",
         "ZN_FastNoiseLiteGradient",

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.

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

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<const float> 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 {

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 <atomic>
+#include <memory>
+#include <vector>
+
+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<const Chunk *> near_chunks;
+	std::vector<const Triangle *> triangles;
+};
+
+struct ChunkGrid {
+	std::vector<Chunk> chunks;
+	Vector3i size;
+	Vector3f min_pos;
+	float chunk_size;
+};
+
+class GenMeshSDFSubBoxTask : public IThreadedTask {
+public:
+	struct SharedData {
+		std::vector<Triangle> 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<SharedData> 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<const Vector3> vertices, Span<const int> indices, std::vector<Triangle> &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<const Triangle> 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<float> sdf_grid, const Vector3i res, Span<const Triangle> 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<float> sdf_grid, const Vector3i res, Span<const Triangle> 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<float> sdf_grid, const Vector3i res, Span<const Triangle> 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<float> sdf_grid, const Vector3i res, Span<const Triangle> 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<const float> sdf_grid, Vector3i res, Span<const Triangle> 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<float> sdf_grid, Vector3i res, Vector3f min_pos, Vector3f max_pos);
+
+} // namespace zylann::voxel::mesh_sdf
+
+#endif // VOXEL_MESH_SDF_H

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 <scene/resources/mesh.h>
+
+namespace zylann::voxel {
+
+static bool prepare_triangles(
+		Mesh &mesh, std::vector<mesh_sdf::Triangle> &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> mesh) {
+	ZN_PROFILE_SCOPE();
+
+	ERR_FAIL_COND(mesh.is_null());
+
+	std::vector<mesh_sdf::Triangle> 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<gd::VoxelBuffer> 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<float> 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> 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<gd::VoxelBuffer> 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<VoxelMeshSDF> 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<VoxelMeshSDF> obj_to_notify;
+
+		void run(ThreadedTaskContext ctx) override {
+			ZN_PROFILE_SCOPE();
+			ZN_ASSERT(obj_to_notify.is_valid());
+
+			std::shared_ptr<mesh_sdf::GenMeshSDFSubBoxTask::SharedData> shared_data =
+					make_shared_instance<mesh_sdf::GenMeshSDFSubBoxTask::SharedData>();
+			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<float> 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<gd::VoxelBuffer>(), 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<gd::VoxelBuffer> 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<gd::VoxelBuffer> 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> mesh) {
+	Array result;
+
+	ZN_ASSERT_RETURN_V(is_baked(), result);
+	ZN_ASSERT(_voxel_buffer.is_valid());
+	const VoxelBufferInternal &buffer = _voxel_buffer->get_buffer();
+	Span<const float> 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<mesh_sdf::Triangle> 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

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 <core/object/ref_counted.h>
+
+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> 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> mesh, SceneTree *scene_tree);
+
+	Ref<gd::VoxelBuffer> get_voxel_buffer() const;
+	AABB get_aabb() const;
+
+	Array debug_check_sdf(Ref<Mesh> mesh);
+
+private:
+	void _on_bake_async_completed(Ref<gd::VoxelBuffer> buffer, Vector3 min_pos, Vector3 max_pos);
+
+	static void _bind_methods();
+
+	// Data
+	Ref<gd::VoxelBuffer> _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

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 <scene/3d/collision_shape_3d.h>
 #include <scene/3d/mesh_instance_3d.h>
@@ -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<VoxelMeshSDF> 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<gd::VoxelBuffer> 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<VoxelDataLodMap> 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<const float> 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<ops::SdfUnion, SdfBufferShape> 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

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<VoxelMeshSDF> mesh_sdf, Transform3D transform, float isolevel, float sdf_scale);
 
 protected:
 	uint64_t _get_voxel(Vector3i pos) const override;

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<FastNoise2>();
 #endif
+	ClassDB::register_class<VoxelMeshSDF>();
 
 	// Meshers
 	ClassDB::register_abstract_class<VoxelMesher>();

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;
 	}
 

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<uint8_t> &slice) const;
 
+	template <typename T>
+	bool get_channel_data(unsigned int channel_index, Span<T> &dst) const {
+		Span<uint8_t> dst8;
+		ZN_ASSERT_RETURN_V(get_channel_raw(channel_index, dst8), false);
+		dst = dst8.reinterpret_cast_to<T>();
+		return true;
+	}
+
 	void downscale_to(VoxelBufferInternal &dst, Vector3i src_min, Vector3i src_max, Vector3i dst_min) const;
 
 	bool equals(const VoxelBufferInternal &p_other) const;

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 {

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);
 }

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

util/math/vector3f.h

@@ -2,7 +2,7 @@
 #define ZYLANN_VECTOR3F_H
 
 #include "../errors.h"
-#include <core/math/math_funcs.h>
+#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