[Data/MeshDistanceField] Added an MDF structure

- This computes signed distances to the nearest mesh geometry within a 3D grid using a BVH

- This can be used for some algorithms, such as rendering ones that can use sphere marching through the 3D volume created from this field

include/RaZ/Data/MeshDistanceField.hpp

@@ -0,0 +1,56 @@
+#pragma once
+
+#ifndef RAZ_MESHDISTANCEFIELD_HPP
+#define RAZ_MESHDISTANCEFIELD_HPP
+
+#include "RaZ/Utils/Shape.hpp"
+
+#include <vector>
+
+namespace Raz {
+
+class BoundingVolumeHierarchy;
+class Image;
+
+/// 3-dimensional structure of signed distances to the closest mesh geometry in a specific area. Distances inside a mesh will be negative.
+class MeshDistanceField {
+public:
+  /// Creates a mesh distance field.
+  /// \param area Area inside which the distances will be computed.
+  /// \param width Number of divisions along the width; must be equal to or greater than 2.
+  /// \param height Number of divisions along the height; must be equal to or greater than 2.
+  /// \param depth Number of divisions along the depth; must be equal to or greater than 2.
+  MeshDistanceField(const AABB& area, unsigned int width, unsigned int height, unsigned int depth);
+
+  float getDistance(std::size_t widthIndex, std::size_t heightIndex, std::size_t depthIndex) const;
+
+  void setBvh(const BoundingVolumeHierarchy& bvh) { m_bvh = &bvh; }
+
+  /// Computes the distance field's values for each point within the grid.
+  /// \param sampleCount Number of directions to sample around each point; a higher count will result in a better definition.
+  /// \note This requires a BVH to have been set.
+  /// \see setBvh()
+  void compute(std::size_t sampleCount);
+  /// Recovers the distance field's values in a list of 2D floating-point images.
+  /// \return Images of each slice of the field along the depth.
+  std::vector<Image> recoverSlices() const;
+
+private:
+  constexpr std::size_t computeIndex(std::size_t widthIndex, std::size_t heightIndex, std::size_t depthIndex) const noexcept {
+    assert("Error: The given width index is invalid." && widthIndex < m_width);
+    assert("Error: The given height index is invalid." && heightIndex < m_height);
+    assert("Error: The given channel depth is invalid." && depthIndex < m_depth);
+    return depthIndex * m_height * m_width + heightIndex * m_width + widthIndex;
+  }
+
+  AABB m_area = AABB(Vec3f(0.f), Vec3f(0.f));
+  unsigned int m_width {};
+  unsigned int m_height {};
+  unsigned int m_depth {};
+  std::vector<float> m_distanceField {};
+  const BoundingVolumeHierarchy* m_bvh = nullptr;
+};
+
+} // namespace Raz
+
+#endif // RAZ_MESHDISTANCEFIELD_HPP

include/RaZ/RaZ.hpp

@@ -26,6 +26,7 @@
 #include "Data/Image.hpp"
 #include "Data/ImageFormat.hpp"
 #include "Data/Mesh.hpp"
+#include "Data/MeshDistanceField.hpp"
 #include "Data/MeshFormat.hpp"
 #include "Data/ObjFormat.hpp"
 #include "Data/OffFormat.hpp"

src/RaZ/Data/MeshDistanceField.cpp

@@ -0,0 +1,76 @@
+#include "RaZ/Data/BoundingVolumeHierarchy.hpp"
+#include "RaZ/Data/Image.hpp"
+#include "RaZ/Data/MeshDistanceField.hpp"
+#include "RaZ/Math/MathUtils.hpp"
+#include "RaZ/Utils/Ray.hpp"
+#include "RaZ/Utils/Threading.hpp"
+
+namespace Raz {
+
+MeshDistanceField::MeshDistanceField(const AABB& area, unsigned int width, unsigned int height, unsigned int depth)
+  : m_area{ area }, m_width{ width }, m_height{ height }, m_depth{ depth }, m_distanceField(width * height * depth, std::numeric_limits<float>::max()) {
+  if (m_width < 2 || m_height < 2 || m_depth < 2)
+    throw std::invalid_argument("[MeshDistanceField] The width, height & depth must all be equal to or greater than 2.");
+}
+
+float MeshDistanceField::getDistance(std::size_t widthIndex, std::size_t heightIndex, std::size_t depthIndex) const {
+  return m_distanceField[computeIndex(widthIndex, heightIndex, depthIndex)];
+}
+
+void MeshDistanceField::compute(std::size_t sampleCount) {
+  if (m_bvh == nullptr)
+    throw std::runtime_error("[MeshDistanceField] Computing a mesh distance field requires having given a BVH.");
+
+  std::fill(m_distanceField.begin(), m_distanceField.end(), std::numeric_limits<float>::max());
+
+  const Vec3f areaExtents = m_area.getMaxPosition() - m_area.getMinPosition();
+  const float widthStep   = areaExtents.x() / static_cast<float>(m_width - 1);
+  const float heightStep  = areaExtents.y() / static_cast<float>(m_height - 1);
+  const float depthStep   = areaExtents.z() / static_cast<float>(m_depth - 1);
+
+  Threading::parallelize(0, m_depth, [this, widthStep, heightStep, depthStep, sampleCount] (const Threading::IndexRange& range) {
+    for (std::size_t depthIndex = range.beginIndex; depthIndex < range.endIndex; ++depthIndex) {
+      for (std::size_t heightIndex = 0; heightIndex < m_height; ++heightIndex) {
+        for (std::size_t widthIndex = 0; widthIndex < m_width; ++widthIndex) {
+          const Vec3f rayPos = m_area.getMinPosition() + Vec3f(static_cast<float>(widthIndex) * widthStep,
+                                                               static_cast<float>(heightIndex) * heightStep,
+                                                               static_cast<float>(depthIndex) * depthStep);
+          float& distance = m_distanceField[computeIndex(widthIndex, heightIndex, depthIndex)];
+
+          for (const Vec3f& rayDir : MathUtils::computeFibonacciSpherePoints(sampleCount)) {
+            RayHit hit {};
+
+            if (!m_bvh->query(Ray(rayPos, rayDir), &hit))
+              continue;
+
+            if (rayDir.dot(hit.normal) > 0.f)
+              hit.distance = -hit.distance;
+
+            if (std::abs(hit.distance) < std::abs(distance))
+              distance = hit.distance;
+          }
+        }
+      }
+    }
+  });
+}
+
+std::vector<Image> MeshDistanceField::recoverSlices() const {
+  std::vector<Image> slices;
+  slices.reserve(m_depth);
+
+  for (std::size_t depthIndex = 0; depthIndex < m_depth; ++depthIndex) {
+    Image& slice = slices.emplace_back(m_width, m_height, ImageColorspace::GRAY, ImageDataType::FLOAT);
+
+    for (std::size_t heightIndex = 0; heightIndex < m_height; ++heightIndex) {
+      for (std::size_t widthIndex = 0; widthIndex < m_width; ++widthIndex) {
+        const float distance = getDistance(widthIndex, heightIndex, depthIndex);
+        slice.setPixel(widthIndex, heightIndex, distance);
+      }
+    }
+  }
+
+  return slices;
+}
+
+} // namespace Raz

src/RaZ/Script/LuaData.cpp

@@ -2,6 +2,8 @@
 #include "RaZ/Data/BoundingVolumeHierarchy.hpp"
 #include "RaZ/Data/BoundingVolumeHierarchySystem.hpp"
 #include "RaZ/Data/Color.hpp"
+#include "RaZ/Data/Image.hpp"
+#include "RaZ/Data/MeshDistanceField.hpp"
 #include "RaZ/Script/LuaWrapper.hpp"
 #include "RaZ/Utils/TypeUtils.hpp"
 
@@ -99,6 +101,17 @@ void LuaWrapper::registerDataTypes() {
     colorPreset["Yellow"]      = ColorPreset::Yellow;
     colorPreset["White"]       = ColorPreset::White;
   }
+
+  {
+    sol::usertype<MeshDistanceField> mdf = state.new_usertype<MeshDistanceField>("MeshDistanceField",
+                                                                                 sol::constructors<
+                                                                                   MeshDistanceField(const Raz::AABB&, unsigned int, unsigned int, unsigned int)
+                                                                                 >());
+    mdf["getDistance"]   = &MeshDistanceField::getDistance;
+    mdf["setBvh"]        = &MeshDistanceField::setBvh;
+    mdf["compute"]       = &MeshDistanceField::compute;
+    mdf["recoverSlices"] = &MeshDistanceField::recoverSlices;
+  }
 }
 
 } // namespace Raz

tests/src/RaZ/Data/MeshDistanceField.cpp

@@ -0,0 +1,175 @@
+#include "RaZ/Entity.hpp"
+#include "RaZ/Data/BoundingVolumeHierarchy.hpp"
+#include "RaZ/Data/Image.hpp"
+#include "RaZ/Data/Mesh.hpp"
+#include "RaZ/Data/MeshDistanceField.hpp"
+
+#include "CatchCustomMatchers.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+
+TEST_CASE("MeshDistanceField computation", "[data]") {
+  // See: https://www.geogebra.org/m/nn5jtkrt
+
+  // The MDF requires a definition of at least 2 on each axis
+  CHECK_THROWS(Raz::MeshDistanceField(Raz::AABB({}, {}), 0, 0, 0));
+  CHECK_THROWS(Raz::MeshDistanceField(Raz::AABB({}, {}), 1, 1, 1));
+  CHECK_NOTHROW(Raz::MeshDistanceField(Raz::AABB({}, {}), 2, 2, 2));
+
+  const Raz::AABB fieldBox(Raz::Vec3f(-1.f), Raz::Vec3f(1.f));
+
+  Raz::MeshDistanceField mdf(fieldBox, 9, 9, 9); // Choosing odd numbers to get the exact center position
+  CHECK(mdf.getDistance(0, 0, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(8, 8, 8) == std::numeric_limits<float>::max());
+
+  CHECK_THROWS(mdf.compute(1)); // No BVH set
+
+  Raz::BoundingVolumeHierarchy bvh;
+  mdf.setBvh(bvh);
+
+  mdf.compute(1); // The BVH is empty, nothing will be computed
+  CHECK(mdf.getDistance(0, 0, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(8, 8, 8) == std::numeric_limits<float>::max());
+
+  Raz::Entity mesh(0);
+  mesh.addComponent<Raz::Mesh>(Raz::AABB(fieldBox.getMinPosition() * 0.5f, fieldBox.getMaxPosition() * 0.5f));
+  bvh.build({ &mesh });
+
+  mdf.compute(1); // A sample count too low isn't enough to find much intersection
+  // Corners
+  CHECK(mdf.getDistance(0, 0, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(8, 8, 8) == std::numeric_limits<float>::max());
+  // Edges
+  CHECK(mdf.getDistance(4, 0, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(0, 4, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(0, 0, 4) == std::numeric_limits<float>::max());
+  // Faces
+  CHECK(mdf.getDistance(0, 4, 4) == 0.5f); // The first sample direction should be +X
+  CHECK(mdf.getDistance(4, 0, 4) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(4, 4, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(4, 4, 8) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(4, 8, 4) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(8, 4, 4) == std::numeric_limits<float>::max());
+  // Center
+  CHECK(mdf.getDistance(4, 4, 4) == -0.5f);
+
+  mdf.compute(2);
+  // Corners
+  CHECK(mdf.getDistance(0, 0, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(8, 8, 8) == 1.f);
+  // Edges
+  CHECK(mdf.getDistance(4, 0, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(0, 4, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(0, 0, 4) == 1.f);
+  // Faces
+  CHECK_THAT(mdf.getDistance(0, 4, 4), IsNearlyEqualTo(0.577350259f));
+  CHECK(mdf.getDistance(4, 0, 4) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(4, 4, 0) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(4, 4, 8) == std::numeric_limits<float>::max());
+  CHECK(mdf.getDistance(4, 8, 4) == std::numeric_limits<float>::max());
+  CHECK_THAT(mdf.getDistance(8, 4, 4), IsNearlyEqualTo(0.782987058f));
+  // Center
+  CHECK_THAT(mdf.getDistance(4, 4, 4), IsNearlyEqualTo(-0.577350259f));
+
+  mdf.compute(100);
+  // Corner to corner, should tend toward 0.87
+  CHECK_THAT(mdf.getDistance(0, 0, 0), IsNearlyEqualTo(0.912753105f));
+  CHECK_THAT(mdf.getDistance(8, 0, 0), IsNearlyEqualTo(0.966534495f));
+  CHECK_THAT(mdf.getDistance(8, 8, 0), IsNearlyEqualTo(0.980392158f));
+  CHECK_THAT(mdf.getDistance(0, 8, 0), IsNearlyEqualTo(1.06382978f));
+  CHECK_THAT(mdf.getDistance(0, 8, 8), IsNearlyEqualTo(1.162790895f));
+  CHECK_THAT(mdf.getDistance(0, 0, 8), IsNearlyEqualTo(1.045123935f));
+  CHECK_THAT(mdf.getDistance(8, 0, 8), IsNearlyEqualTo(1.048231006f));
+  CHECK_THAT(mdf.getDistance(8, 8, 8), IsNearlyEqualTo(1.132885695f));
+  // Edge to edge, should tend toward 0.71
+  CHECK_THAT(mdf.getDistance(4, 0, 0), IsNearlyEqualTo(0.724637687f));
+  CHECK_THAT(mdf.getDistance(0, 4, 0), IsNearlyEqualTo(0.772372246f));
+  CHECK_THAT(mdf.getDistance(0, 0, 4), IsNearlyEqualTo(0.749056816f));
+  CHECK_THAT(mdf.getDistance(4, 8, 8), IsNearlyEqualTo(0.761324167f));
+  CHECK_THAT(mdf.getDistance(8, 4, 8), IsNearlyEqualTo(0.76935935f));
+  CHECK_THAT(mdf.getDistance(8, 8, 4), IsNearlyEqualTo(0.725454986f));
+  // Side to side, should tend toward 0.5
+  CHECK_THAT(mdf.getDistance(4, 4, 0), IsNearlyEqualTo(0.510620058f));
+  CHECK_THAT(mdf.getDistance(0, 4, 4), IsNearlyEqualTo(0.503709972f));
+  CHECK_THAT(mdf.getDistance(4, 0, 4), IsNearlyEqualTo(0.50505048f));
+  CHECK_THAT(mdf.getDistance(4, 4, 8), IsNearlyEqualTo(0.501562536f));
+  CHECK_THAT(mdf.getDistance(8, 4, 4), IsNearlyEqualTo(0.504095972f));
+  CHECK_THAT(mdf.getDistance(4, 8, 4), IsNearlyEqualTo(0.50505048f));
+  // Inside at mid-distance between center & corners, should tend toward -0.25
+  CHECK_THAT(mdf.getDistance(3, 3, 3), IsNearlyEqualTo(-0.250781268f));
+  CHECK_THAT(mdf.getDistance(3, 3, 5), IsNearlyEqualTo(-0.252047986f));
+  CHECK_THAT(mdf.getDistance(3, 5, 3), IsNearlyEqualTo(-0.250781268f));
+  CHECK_THAT(mdf.getDistance(3, 5, 5), IsNearlyEqualTo(-0.252047986f));
+  CHECK_THAT(mdf.getDistance(5, 3, 3), IsNearlyEqualTo(-0.250781268f));
+  CHECK_THAT(mdf.getDistance(5, 3, 5), IsNearlyEqualTo(-0.251854986f));
+  CHECK_THAT(mdf.getDistance(5, 5, 3), IsNearlyEqualTo(-0.250781268f));
+  CHECK_THAT(mdf.getDistance(5, 5, 5), IsNearlyEqualTo(-0.251854986f));
+  // Center to side, should tend toward -0.5
+  CHECK_THAT(mdf.getDistance(4, 4, 4), IsNearlyEqualTo(-0.501562536f));
+
+  // Degenerate cases due to a grid definition too low, should tend toward 0 but fail to find an intersection other than the closest opposite inner face:
+
+  // Right on each corner
+  CHECK_THAT(mdf.getDistance(2, 2, 2), IsNearlyEqualTo(-1.010100961f));
+  CHECK_THAT(mdf.getDistance(2, 2, 6), IsNearlyEqualTo(-1.048736811f));
+  CHECK_THAT(mdf.getDistance(2, 6, 2), IsNearlyEqualTo(-1.010100961f));
+  CHECK_THAT(mdf.getDistance(2, 6, 6), IsNearlyEqualTo(-1.003125072f));
+  CHECK_THAT(mdf.getDistance(6, 2, 2), IsNearlyEqualTo(-1.022850156f));
+  CHECK_THAT(mdf.getDistance(6, 2, 6), IsNearlyEqualTo(-1.030927777f));
+  CHECK_THAT(mdf.getDistance(6, 6, 2), IsNearlyEqualTo(-1.045168638f));
+  CHECK_THAT(mdf.getDistance(6, 6, 6), IsNearlyEqualTo(-1.008191943f));
+  // Right in the middle of each edge
+  CHECK_THAT(mdf.getDistance(2, 4, 2), IsNearlyEqualTo(-0.50505048f));
+  CHECK_THAT(mdf.getDistance(2, 4, 6), IsNearlyEqualTo(-0.526315749f));
+  CHECK_THAT(mdf.getDistance(4, 2, 2), IsNearlyEqualTo(-0.515975773f));
+  CHECK_THAT(mdf.getDistance(4, 2, 6), IsNearlyEqualTo(-0.524368405f));
+  CHECK_THAT(mdf.getDistance(4, 6, 2), IsNearlyEqualTo(-0.522584319f));
+  CHECK_THAT(mdf.getDistance(4, 6, 6), IsNearlyEqualTo(-0.503709972f));
+  CHECK_THAT(mdf.getDistance(6, 4, 2), IsNearlyEqualTo(-0.537634432f));
+  CHECK_THAT(mdf.getDistance(6, 4, 6), IsNearlyEqualTo(-0.515463889f));
+  // Right in the middle of each face
+  CHECK_THAT(mdf.getDistance(2, 4, 4), IsNearlyEqualTo(-0.501562536f));
+  CHECK_THAT(mdf.getDistance(6, 4, 4), IsNearlyEqualTo(-0.511425078f));
+  CHECK_THAT(mdf.getDistance(4, 2, 4), IsNearlyEqualTo(-0.511425078f));
+  CHECK_THAT(mdf.getDistance(4, 6, 4), IsNearlyEqualTo(-0.501562536f));
+  CHECK_THAT(mdf.getDistance(4, 4, 2), IsNearlyEqualTo(-0.50505048f));
+  CHECK_THAT(mdf.getDistance(4, 4, 6), IsNearlyEqualTo(-0.503709972f));
+}
+
+TEST_CASE("MeshDistanceField slices", "[data]") {
+  // Creating a distance field with a single triangle inside
+  //
+  // -----^-----
+  // |  /   \  |
+  // |/_______\|
+
+  Raz::Entity mesh(0);
+  mesh.addComponent<Raz::Mesh>(Raz::Triangle(Raz::Vec3f(-1.f, -0.5f, 0.f), Raz::Vec3f(1.f, -0.5f, 0.f), Raz::Vec3f(0.f, 0.5f, 0.f)),
+                               Raz::Vec2f(), Raz::Vec2f(), Raz::Vec2f());
+
+  Raz::BoundingVolumeHierarchy bvh;
+  bvh.build({ &mesh });
+
+  Raz::MeshDistanceField mdf(Raz::AABB(Raz::Vec3f(-1.f, -0.5f, -0.25f), Raz::Vec3f(1.f, 0.5f, 0.25f)), 4, 3, 2);
+  mdf.setBvh(bvh);
+  mdf.compute(10);
+
+  constexpr float backDistance  = -0.289784729f;
+  constexpr float frontDistance = 0.271375328f;
+
+  CHECK_THAT(mdf.getDistance(1, 1, 0), IsNearlyEqualTo(backDistance));
+  CHECK_THAT(mdf.getDistance(1, 1, 1), IsNearlyEqualTo(frontDistance));
+
+  const std::vector<Raz::Image> imageSlices = mdf.recoverSlices();
+  REQUIRE(imageSlices.size() == 2); // The number of slices is equal to the MDF's depth definition
+
+  for (const Raz::Image& depthSlice : imageSlices) {
+    CHECK(depthSlice.getWidth() == 4);
+    CHECK(depthSlice.getHeight() == 3);
+    CHECK(depthSlice.getColorspace() == Raz::ImageColorspace::GRAY);
+    CHECK(depthSlice.getDataType() == Raz::ImageDataType::FLOAT);
+  }
+
+  CHECK_THAT(imageSlices[0].recoverPixel<float>(1, 1), IsNearlyEqualTo(backDistance));
+  CHECK_THAT(imageSlices[1].recoverPixel<float>(1, 1), IsNearlyEqualTo(frontDistance));
+}

tests/src/RaZ/Script/LuaData.cpp

@@ -156,6 +156,17 @@ TEST_CASE("LuaData Mesh", "[script][lua][data]") {
 #endif
 }
 
+TEST_CASE("LuaData MeshDistanceField", "[script][lua][data]") {
+  CHECK(Raz::LuaWrapper::execute(R"(
+    local mdf = MeshDistanceField.new(AABB.new(Vec3f.new(), Vec3f.new()), 2, 2, 2)
+
+    assert(mdf:getDistance(0, 0, 0) ~= 0)
+    mdf:setBvh(BoundingVolumeHierarchy.new())
+    mdf:compute(1)
+    assert(mdf:recoverSlices():size() == 2)
+  )"));
+}
+
 TEST_CASE("LuaData Submesh", "[script][lua][data]") {
   CHECK(Raz::LuaWrapper::execute(R"(
     local submesh = Submesh.new()