Added voxel raycast to VoxelToolLodTerrain

- Made more precise with an option to use binary search (unoptimized)
- Distance along ray is now available in blocky terrains as well

doc/source/changelog.md

@@ -34,12 +34,14 @@ Ongoing development - `master`
     - Added `SdfSphereHeightmap` and `Normalize` nodes to voxel graph, which can help making planets
     - Added `SdfSmoothUnion` and `SdfSmoothSubtract` nodes to voxel graph
     - Added `VoxelInstancer` to instantiate items on top of `VoxelLodTerrain`, aimed at spawning natural elements such as rocks and foliage
+    - Implemented `VoxelToolLodterrain.raycast()`
 
 - Blocky voxels
     - Introduced a second blocky mesher dedicated to colored cubes, with greedy meshing and palette support
     - Replaced `transparent` property with `transparency_index` for more control on the culling of transparent faces
     - The TYPE channel is now 16-bit by default instead of 8-bit, allowing to store up to 65,536 types (part of this channel might actually be used to store rotation in the future)
     - Added normalmaps support
+    - `VoxelRaycastResult` now also contains hit distance
 
 - Breaking changes
     - `VoxelViewer` now replaces the `viewer_path` property on `VoxelTerrain`, and allows multiple loading points

edition/voxel_tool.cpp

@@ -12,15 +12,21 @@ Vector3 VoxelRaycastResult::_b_get_previous_position() const {
 	return previous_position.to_vec3();
 }
 
+float VoxelRaycastResult::_b_get_distance() const {
+	return distance_along_ray;
+}
+
 void VoxelRaycastResult::_bind_methods() {
 	ClassDB::bind_method(D_METHOD("get_position"), &VoxelRaycastResult::_b_get_position);
 	ClassDB::bind_method(D_METHOD("get_previous_position"), &VoxelRaycastResult::_b_get_previous_position);
+	ClassDB::bind_method(D_METHOD("get_distance"), &VoxelRaycastResult::_b_get_distance);
 
 	ADD_PROPERTY(PropertyInfo(Variant::VECTOR3, "position"), "", "get_position");
 	ADD_PROPERTY(PropertyInfo(Variant::VECTOR3, "previous_position"), "", "get_previous_position");
+	ADD_PROPERTY(PropertyInfo(Variant::REAL, "distance"), "", "get_distance");
 }
 
-//----------------------------------------
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 VoxelTool::VoxelTool() {
 	_sdf_scale = VoxelBuffer::get_sdf_quantization_scale(VoxelBuffer::DEFAULT_SDF_CHANNEL_DEPTH);

edition/voxel_tool.h

@@ -12,10 +12,12 @@ class VoxelRaycastResult : public Reference {
 public:
 	Vector3i position;
 	Vector3i previous_position;
+	float distance_along_ray;
 
 private:
 	Vector3 _b_get_position() const;
 	Vector3 _b_get_previous_position() const;
+	float _b_get_distance() const;
 
 	static void _bind_methods();
 };

edition/voxel_tool_lod_terrain.cpp

@@ -1,9 +1,10 @@
 #include "voxel_tool_lod_terrain.h"
 #include "../terrain/voxel_lod_terrain.h"
 #include "../terrain/voxel_map.h"
+#include "../util/voxel_raycast.h"
 
 VoxelToolLodTerrain::VoxelToolLodTerrain(VoxelLodTerrain *terrain, VoxelMap &map) :
-		_terrain(terrain), _map(map) {
+		_terrain(terrain), _map(&map) {
 	ERR_FAIL_COND(terrain == nullptr);
 	// At the moment, only LOD0 is supported.
 	// Don't destroy the terrain while a voxel tool still references it
@@ -12,30 +13,178 @@ VoxelToolLodTerrain::VoxelToolLodTerrain(VoxelLodTerrain *terrain, VoxelMap &map
 bool VoxelToolLodTerrain::is_area_editable(const Rect3i &box) const {
 	ERR_FAIL_COND_V(_terrain == nullptr, false);
 	// TODO Take volume bounds into account
-	return _map.is_area_fully_loaded(box.padded(1));
+	return _map->is_area_fully_loaded(box.padded(1));
+}
+
+template <typename Volume_F>
+float get_sdf_interpolated(Volume_F &f, Vector3 pos) {
+	const Vector3i c = Vector3i::from_floored(pos);
+
+	const float s000 = f(Vector3i(c.x, c.y, c.z));
+	const float s100 = f(Vector3i(c.x + 1, c.y, c.z));
+	const float s010 = f(Vector3i(c.x, c.y + 1, c.z));
+	const float s110 = f(Vector3i(c.x + 1, c.y + 1, c.z));
+	const float s001 = f(Vector3i(c.x, c.y, c.z + 1));
+	const float s101 = f(Vector3i(c.x + 1, c.y, c.z + 1));
+	const float s011 = f(Vector3i(c.x, c.y + 1, c.z + 1));
+	const float s111 = f(Vector3i(c.x + 1, c.y + 1, c.z + 1));
+
+	return interpolate(s000, s100, s101, s001, s010, s110, s111, s011, fract(pos));
+}
+
+// Binary search can be more accurate than linear regression because the SDF can be inaccurate in the first place.
+// An alternative would be to polygonize a tiny area around the middle-phase hit position.
+// `d1` is how far from `pos0` along `dir` the binary search will take place.
+// The segment may be adjusted internally if it does not contain a zero-crossing of the
+template <typename Volume_F>
+float approximate_distance_to_isosurface_binary_search(
+		Volume_F &f, Vector3 pos0, Vector3 dir, float d1, int iterations) {
+
+	float d0 = 0.f;
+	float sdf0 = get_sdf_interpolated(f, pos0);
+	// The position given as argument may be a rough approximation coming from the middle-phase,
+	// so it can be slightly below the surface. We can adjust it a little so it is above.
+	for (int i = 0; i < 4 && sdf0 < 0.f; ++i) {
+		d0 -= 0.5f;
+		sdf0 = get_sdf_interpolated(f, pos0 + dir * d0);
+	}
+
+	float sdf1 = get_sdf_interpolated(f, pos0 + dir * d1);
+	for (int i = 0; i < 4 && sdf1 > 0.f; ++i) {
+		d1 += 0.5f;
+		sdf1 = get_sdf_interpolated(f, pos0 + dir * d1);
+	}
+
+	if ((sdf0 > 0) != (sdf1 > 0)) {
+		// Binary search
+		for (int i = 0; i < iterations; ++i) {
+			const float dm = 0.5f * (d0 + d1);
+			const float sdf_mid = get_sdf_interpolated(f, pos0 + dir * dm);
+
+			if ((sdf_mid > 0) != (sdf0 > 0)) {
+				sdf1 = sdf_mid;
+				d1 = dm;
+			} else {
+				sdf0 = sdf_mid;
+				d0 = dm;
+			}
+		}
+	}
+
+	// Pick distance closest to the surface
+	if (Math::abs(sdf0) < Math::abs(sdf1)) {
+		return d0;
+	} else {
+		return d1;
+	}
+}
+
+Ref<VoxelRaycastResult> VoxelToolLodTerrain::raycast(
+		Vector3 pos, Vector3 dir, float max_distance, uint32_t collision_mask) {
+	// TODO Transform input if the terrain is rotated
+	// TODO Implement broad-phase on blocks to minimize locking and increase performance
+
+	struct RaycastPredicate {
+		const VoxelMap &map;
+
+		bool operator()(Vector3i pos) {
+			// This is not particularly optimized, but runs fast enough for player raycasts
+			const float sdf = map.get_voxel_f(pos, VoxelBuffer::CHANNEL_SDF);
+			return sdf < 0;
+		}
+	};
+
+	Ref<VoxelRaycastResult> res;
+
+	// We use grid-raycast as a middle-phase to roughly detect where the hit will be
+	RaycastPredicate predicate = { *_map };
+	Vector3i hit_pos;
+	Vector3i prev_pos;
+	float hit_distance;
+	float hit_distance_prev;
+	// Voxels polygonized using marching cubes influence a region centered on their lower corner,
+	// and extend up to 0.5 units in all directions.
+	//
+	//   o--------o--------o
+	//   | A      |     B  |  Here voxel B is full, voxels A, C and D are empty.
+	//   |       xxx       |  Matter will show up at the lower corner of B due to interpolation.
+	//   |     xxxxxxx     |
+	//   o---xxxxxoxxxxx---o
+	//   |     xxxxxxx     |
+	//   |       xxx       |
+	//   | C      |     D  |
+	//   o--------o--------o
+	//
+	// `voxel_raycast` operates on a discrete grid of cubic voxels, so to account for the smooth interpolation,
+	// we may offset the ray so that cubes act as if they were centered on the filtered result.
+	const Vector3 offset(0.5, 0.5, 0.5);
+	if (voxel_raycast(pos + offset, dir, predicate, max_distance, hit_pos, prev_pos, hit_distance, hit_distance_prev)) {
+		// Approximate surface
+
+		float d = hit_distance;
+
+		if (_raycast_binary_search_iterations > 0) {
+			// This is not particularly optimized, but runs fast enough for player raycasts
+			struct VolumeSampler {
+				const VoxelMap &map;
+
+				inline float operator()(const Vector3i &pos) {
+					return map.get_voxel_f(pos, VoxelBuffer::CHANNEL_SDF);
+				}
+			};
+
+			VolumeSampler sampler{ *_map };
+			d = hit_distance_prev + approximate_distance_to_isosurface_binary_search(sampler,
+											pos + dir * hit_distance_prev,
+											dir, hit_distance - hit_distance_prev,
+											_raycast_binary_search_iterations);
+		}
+
+		res.instance();
+		res->position = hit_pos;
+		res->previous_position = prev_pos;
+		res->distance_along_ray = d;
+	}
+
+	return res;
 }
 
 uint64_t VoxelToolLodTerrain::_get_voxel(Vector3i pos) const {
 	ERR_FAIL_COND_V(_terrain == nullptr, 0);
-	return _map.get_voxel(pos, _channel);
+	return _map->get_voxel(pos, _channel);
 }
 
 float VoxelToolLodTerrain::_get_voxel_f(Vector3i pos) const {
 	ERR_FAIL_COND_V(_terrain == nullptr, 0);
-	return _map.get_voxel_f(pos, _channel);
+	return _map->get_voxel_f(pos, _channel);
 }
 
 void VoxelToolLodTerrain::_set_voxel(Vector3i pos, uint64_t v) {
 	ERR_FAIL_COND(_terrain == nullptr);
-	_map.set_voxel(v, pos, _channel);
+	_map->set_voxel(v, pos, _channel);
 }
 
 void VoxelToolLodTerrain::_set_voxel_f(Vector3i pos, float v) {
 	ERR_FAIL_COND(_terrain == nullptr);
-	_map.set_voxel_f(v, pos, _channel);
+	_map->set_voxel_f(v, pos, _channel);
 }
 
 void VoxelToolLodTerrain::_post_edit(const Rect3i &box) {
 	ERR_FAIL_COND(_terrain == nullptr);
 	_terrain->post_edit_area(box);
 }
+
+int VoxelToolLodTerrain::get_raycast_binary_search_iterations() const {
+	return _raycast_binary_search_iterations;
+}
+
+void VoxelToolLodTerrain::set_raycast_binary_search_iterations(int iterations) {
+	_raycast_binary_search_iterations = clamp(iterations, 0, 16);
+}
+
+void VoxelToolLodTerrain::_bind_methods() {
+	ClassDB::bind_method(D_METHOD("set_raycast_binary_search_iterations", "iterations"),
+			&VoxelToolLodTerrain::set_raycast_binary_search_iterations);
+	ClassDB::bind_method(D_METHOD("get_raycast_binary_search_iterations"),
+			&VoxelToolLodTerrain::get_raycast_binary_search_iterations);
+}

edition/voxel_tool_lod_terrain.h

@@ -9,9 +9,14 @@ class VoxelMap;
 class VoxelToolLodTerrain : public VoxelTool {
 	GDCLASS(VoxelToolLodTerrain, VoxelTool)
 public:
+	VoxelToolLodTerrain() {}
 	VoxelToolLodTerrain(VoxelLodTerrain *terrain, VoxelMap &map);
 
 	bool is_area_editable(const Rect3i &box) const override;
+	Ref<VoxelRaycastResult> raycast(Vector3 pos, Vector3 dir, float max_distance, uint32_t collision_mask) override;
+
+	int get_raycast_binary_search_iterations() const;
+	void set_raycast_binary_search_iterations(int iterations);
 
 protected:
 	uint64_t _get_voxel(Vector3i pos) const override;
@@ -21,8 +26,11 @@ protected:
 	void _post_edit(const Rect3i &box) override;
 
 private:
+	static void _bind_methods();
+
 	VoxelLodTerrain *_terrain = nullptr;
-	VoxelMap &_map;
+	VoxelMap *_map = nullptr;
+	int _raycast_binary_search_iterations = 0;
 };
 
 #endif // VOXEL_TOOL_LOD_TERRAIN_H

edition/voxel_tool_terrain.cpp

@@ -71,10 +71,13 @@ Ref<VoxelRaycastResult> VoxelToolTerrain::raycast(Vector3 pos, Vector3 dir, floa
 	Vector3i prev_pos;
 
 	RaycastPredicate predicate = { *_terrain, **library_ref, collision_mask };
-	if (voxel_raycast(pos, dir, predicate, max_distance, hit_pos, prev_pos)) {
+	float hit_distance;
+	float hit_distance_prev;
+	if (voxel_raycast(pos, dir, predicate, max_distance, hit_pos, prev_pos, hit_distance, hit_distance_prev)) {
 		res.instance();
 		res->position = hit_pos;
 		res->previous_position = prev_pos;
+		res->distance_along_ray = hit_distance;
 	}
 
 	return res;

register_types.cpp

@@ -1,6 +1,7 @@
 #include "register_types.h"
 #include "edition/voxel_tool.h"
 #include "edition/voxel_tool_buffer.h"
+#include "edition/voxel_tool_lod_terrain.h"
 #include "edition/voxel_tool_terrain.h"
 #include "editor/editor_plugin.h"
 #include "editor/fast_noise_lite/fast_noise_lite_editor_plugin.h"
@@ -100,6 +101,7 @@ void register_voxel_types() {
 	ClassDB::register_class<VoxelRaycastResult>();
 	ClassDB::register_class<VoxelTool>();
 	ClassDB::register_class<VoxelToolTerrain>();
+	ClassDB::register_class<VoxelToolLodTerrain>();
 	// I had to bind this one despite it being useless as-is because otherwise Godot lazily initializes its class.
 	// And this can happen in a thread, causing crashes due to the concurrent access
 	ClassDB::register_class<VoxelToolBuffer>();

terrain/voxel_map.cpp

@@ -125,7 +125,7 @@ VoxelBlock *VoxelMap::get_block(Vector3i bpos) {
 }
 
 const VoxelBlock *VoxelMap::get_block(Vector3i bpos) const {
-	if (_last_accessed_block && _last_accessed_block->position == bpos) {
+	if (_last_accessed_block != nullptr && _last_accessed_block->position == bpos) {
 		return _last_accessed_block;
 	}
 	const unsigned int *iptr = _blocks_map.getptr(bpos);

util/math/funcs.h

@@ -4,7 +4,19 @@
 #include <core/math/vector3.h>
 
 // Trilinear interpolation between corner values of a cube.
-// Cube points respect the same position as in octree_tables.h
+//
+//      6---------------7
+//     /|              /|
+//    / |             / |
+//   5---------------4  |
+//   |  |            |  |
+//   |  |            |  |
+//   |  |            |  |
+//   |  2------------|--3        Y
+//   | /             | /         | Z
+//   |/              |/          |/
+//   1---------------0      X----o
+//
 template <typename T>
 inline T interpolate(const T v0, const T v1, const T v2, const T v3, const T v4, const T v5, const T v6, const T v7,
 		Vector3 position) {
@@ -119,4 +131,12 @@ inline float smoothstep(float p_from, float p_to, float p_weight) {
 	return x * x * (3.0f - 2.0f * x);
 }
 
+inline float fract(float x) {
+	return x - Math::floor(x);
+}
+
+inline Vector3 fract(const Vector3 &p) {
+	return Vector3(fract(p.x), fract(p.y), fract(p.z));
+}
+
 #endif // VOXEL_MATH_FUNCS_H

util/math/vector3i.h

@@ -46,6 +46,13 @@ struct Vector3i {
 		z = Math::floor(f.z);
 	}
 
+	static inline Vector3i from_floored(const Vector3 &f) {
+		return Vector3i(
+				Math::floor(f.x),
+				Math::floor(f.y),
+				Math::floor(f.z));
+	}
+
 	_FORCE_INLINE_ Vector3 to_vec3() const {
 		return Vector3(x, y, z);
 	}

util/voxel_raycast.h

@@ -9,7 +9,9 @@ bool voxel_raycast(
 		Predicate_F predicate,
 		real_t max_distance,
 		Vector3i &out_hit_pos,
-		Vector3i &out_prev_pos) {
+		Vector3i &out_prev_pos,
+		float &out_distance_along_ray,
+		float &out_distance_along_ray_prev) {
 
 	VOXEL_PROFILE_SCOPE();
 
@@ -112,8 +114,12 @@ bool voxel_raycast(
 
 	/* Iteration */
 
+	float t = 0.f;
+	float t_prev = 0.f;
+
 	do {
 		hit_prev_pos = hit_pos;
+		t_prev = t;
 		if (tcross_x < tcross_y) {
 			if (tcross_x < tcross_z) {
 				// X collision
@@ -122,6 +128,7 @@ bool voxel_raycast(
 				if (tcross_x > max_distance) {
 					return false;
 				}
+				t = tcross_x;
 				tcross_x += tdelta_x;
 			} else {
 				// Z collision (duplicate code)
@@ -130,6 +137,7 @@ bool voxel_raycast(
 				if (tcross_z > max_distance) {
 					return false;
 				}
+				t = tcross_z;
 				tcross_z += tdelta_z;
 			}
 		} else {
@@ -140,6 +148,7 @@ bool voxel_raycast(
 				if (tcross_y > max_distance) {
 					return false;
 				}
+				t = tcross_y;
 				tcross_y += tdelta_y;
 			} else {
 				// Z collision (duplicate code)
@@ -148,6 +157,7 @@ bool voxel_raycast(
 				if (tcross_z > max_distance) {
 					return false;
 				}
+				t = tcross_z;
 				tcross_z += tdelta_z;
 			}
 		}
@@ -156,6 +166,8 @@ bool voxel_raycast(
 
 	out_hit_pos = hit_pos;
 	out_prev_pos = hit_prev_pos;
+	out_distance_along_ray = t;
+	out_distance_along_ray_prev = t_prev;
 
 	return true;
 }