Merge pull request #1035 from MatPont/sdf

[SignedDistanceField] add module
topology-tool-kit · Jul 4, 2024 · cb7c18e · cb7c18e
2 parents 81a959b + 25c0929
commit cb7c18e
Show file tree

Hide file tree

Showing 14 changed files with 1,819 additions and 472 deletions.
diff --git a/core/base/boundingVolumeHierarchy/BoundingVolumeHierarchy.h b/core/base/boundingVolumeHierarchy/BoundingVolumeHierarchy.h
@@ -281,7 +281,10 @@ namespace ttk {
                    const IT *connectivityList,
                    const float *vertexCoords,
                    int *triangleIndex,
-                   float *distance) const {
+                   float *distance,
+                   std::vector<int> &triangles,
+                   std::vector<float> &distances,
+                   bool segmentIntersection = false) const {
       bool wasHit = false;
       float nearestTriangle = std::numeric_limits<float>::max();
       std::stack<Node *> stack;
@@ -308,11 +311,17 @@ namespace ttk {
               v1 *= 3;
               v2 *= 3;
               hasHit = MollerTrumbore(r, v0, v1, v2, vertexCoords);
-              if(hasHit && r.distance < nearestTriangle) {
-                *triangleIndex = triIdx;
-                nearestTriangle = r.distance;
+              if(hasHit
+                 and (not segmentIntersection
+                      or (r.distance >= 0 and r.distance <= 1))) {
                 wasHit = true;
-                *distance = r.distance;
+                triangles.emplace_back(triIdx);
+                distances.emplace_back(r.distance);
+                if(r.distance < nearestTriangle) {
+                  *triangleIndex = triIdx;
+                  nearestTriangle = r.distance;
+                  *distance = r.distance;
+                }
               }
             }
           } else {
@@ -329,6 +338,30 @@ namespace ttk {
       return wasHit;
     }
 
+    bool intersect(Ray &r,
+                   const IT *connectivityList,
+                   const float *vertexCoords,
+                   int *triangleIndex,
+                   float *distance,
+                   bool segmentIntersection = false) const {
+      std::vector<int> triangles;
+      std::vector<float> distances;
+      return intersect(r, connectivityList, vertexCoords, triangleIndex,
+                       distance, triangles, distances, segmentIntersection);
+    }
+
+    bool intersect(Ray &r,
+                   const IT *connectivityList,
+                   const float *vertexCoords,
+                   std::vector<int> &triangles,
+                   std::vector<float> &distances,
+                   bool segmentIntersection = false) const {
+      int triangleIndex;
+      float distance;
+      return intersect(r, connectivityList, vertexCoords, &triangleIndex,
+                       &distance, triangles, distances, segmentIntersection);
+    }
+
     bool wasNodeHit(const Ray &r, Node *n) const {
       float tmin = (n->m_minX - r.m_origin[0]) / r.m_direction[0];
       float tmax = (n->m_maxX - r.m_origin[0]) / r.m_direction[0];

diff --git a/core/base/geometry/Geometry.cpp b/core/base/geometry/Geometry.cpp
@@ -534,6 +534,60 @@ T Geometry::magnitude(const T *o, const T *d) {
   return sqrt(mag);
 }
 
+template <typename T>
+void Geometry::projectOnTrianglePlane(const T *p,
+                                      const T *a,
+                                      const T *normTri,
+                                      T *out) {
+  std::array<T, 3> ap{};
+  subtractVectors(a, p, ap.data());
+  std::array<T, 3> normTriScaled{};
+  scaleVector(normTri, dotProduct(normTri, ap.data()), normTriScaled.data());
+  subtractVectors(normTriScaled.data(), p, out);
+}
+
+template <typename T>
+void Geometry::projectOnEdge(const T *p, const T *a, const T *b, T *out) {
+  std::array<T, 3> ab{};
+  subtractVectors(a, b, ab.data());
+  std::array<T, 3> ap{};
+  subtractVectors(a, p, ap.data());
+  std::array<T, 3> abScaled{};
+  scaleVector(
+    ab.data(),
+    dotProduct(ap.data(), ab.data()) / dotProduct(ab.data(), ab.data()),
+    abScaled.data());
+  addVectors(a, abScaled.data(), out);
+}
+
+template <typename T>
+void Geometry::computeTriangleNormal(const T *a,
+                                     const T *b,
+                                     const T *c,
+                                     T *out) {
+
+  // triangle normal: cross product of two edges
+  // ab, ac vectors
+  std::array<T, 3> ab{};
+  subtractVectors(a, b, ab.data());
+  std::array<T, 3> ac{};
+  subtractVectors(a, c, ac.data());
+  // compute ab ^ ac
+  crossProduct(ab.data(), ac.data(), out);
+  // magnitude
+  const auto mag = magnitude(out);
+
+  if(mag > powf(10, -FLT_DIG)) {
+    // unitary normal vector
+    scaleVector(out, 1 / mag, out);
+    return;
+  }
+
+  out[0] = -1.0F;
+  out[1] = -1.0F;
+  out[2] = -1.0F;
+}
+
 template <typename T>
 int Geometry::subtractVectors(const T *a,
                               const T *b,
@@ -831,6 +885,12 @@ void Geometry::transposeMatrix(const std::vector<std::vector<T>> &a,
   template TYPE Geometry::magnitudeFlatten<TYPE>(                              \
     std::vector<std::vector<TYPE>> const &);                                   \
   template TYPE Geometry::magnitude<TYPE>(TYPE const *, TYPE const *);         \
+  template void Geometry::projectOnTrianglePlane<TYPE>(                        \
+    TYPE const *, TYPE const *, TYPE const *, TYPE *);                         \
+  template void Geometry::projectOnEdge<TYPE>(                                 \
+    TYPE const *, TYPE const *, TYPE const *, TYPE *);                         \
+  template void Geometry::computeTriangleNormal<TYPE>(                         \
+    TYPE const *, TYPE const *, TYPE const *, TYPE *);                         \
   template int Geometry::subtractVectors<TYPE>(                                \
     TYPE const *, TYPE const *, TYPE *, int const &);                          \
   template int Geometry::subtractVectors<TYPE>(std::vector<TYPE> const &,      \