fix bugs

.vscode/tasks.json

@@ -30,6 +30,12 @@
             "type": "shell",
             "command": "xmake project -k compile_commands .vscode",
             "problemMatcher": []
+        },
+        {
+            "label": "debug frontend application",
+            "type": "shell",
+            "command": "xmake run -d implicit_surface_integrator",
+            "problemMatcher": []
         }
     ]
 }

3rdparty/xmake.lua

@@ -5,7 +5,6 @@ target("indirect_predicates")
     add_vectorexts("neon")
     add_vectorexts("avx")
     add_vectorexts("avx2")
-    add_vectorexts("avx512")
     add_vectorexts("sse")
     add_vectorexts("sse2")
     add_vectorexts("sse3")
@@ -27,7 +26,6 @@ rule("config.indirect_predicates.flags")
         target:add("vectorexts", "neon")
         target:add("vectorexts", "avx")
         target:add("vectorexts", "avx2")
-        target:add("vectorexts", "avx512")
         target:add("vectorexts", "sse")
         target:add("vectorexts", "sse2")
         target:add("vectorexts", "sse3")

application/implicit_surface_integrator.hpp

@@ -0,0 +1,41 @@
+#pragma once
+
+#include <implicit_surface_network_processor.hpp>
+
+class ImplicitSurfaceIntegrator : public ImplicitSurfaceNetworkProcessor
+{
+public:
+    using ImplicitSurfaceNetworkProcessor::ImplicitSurfaceNetworkProcessor;
+
+    void clear_results()
+    {
+        m_per_face_surf_int.clear();
+        m_per_face_vol_int.clear();
+    }
+
+    void compute()
+    {
+        m_per_face_surf_int.reserve(iso_faces.size());
+        m_per_face_vol_int.reserve(iso_faces.size());
+
+        Eigen::Vector3d temp_area_vector{};
+        for (const auto& face : iso_faces) {
+            const auto&     vertices = face.vertex_indices;
+            const auto&     v0       = iso_vertices[vertices[0]];
+            double          surf_int{};
+            Eigen::Vector3d area_vector_sum{};
+            for (auto iter = vertices.begin() + 2; iter != vertices.end(); ++iter) {
+                const auto &v1 = iso_vertices[*(iter - 1)], v2 = iso_vertices[*iter];
+                temp_area_vector  = (v1 - v0).cross(v2 - v0);
+                area_vector_sum  += temp_area_vector;
+                surf_int         += temp_area_vector.norm();
+            }
+            m_per_face_surf_int.emplace_back(surf_int / 2.0);
+            m_per_face_vol_int.emplace_back(area_vector_sum.dot(v0) / 6.0);
+        }
+    }
+
+private:
+    stl_vector_mp<double> m_per_face_surf_int{};
+    stl_vector_mp<double> m_per_face_vol_int{};
+};

application/main.cpp

@@ -1,6 +1,37 @@
-#include <implicit_functions.hpp>
-#include <implicit_surface_network_processor.hpp>
+#include <load_functions.hpp>
+#include "implicit_surface_integrator.hpp"
+
+int main()
+{
+    std::cout << "generating tetrahedron background mesh..." << std::endl;
+    raw_point_t aabb_min{-1.0, -1.0, -1.0};
+    raw_point_t aabb_max{1.0, 1.0, 1.0};
+    auto        background_mesh = generate_tetrahedron_background_mesh(21, aabb_min, aabb_max);
+    std::cout << "loading lookup tables..." << std::endl;
+    load_lut();
+    std::cout << "reading implicit functions..." << std::endl;
+    auto implicit_functions = load_functions("functions.json");
+    if (!implicit_functions.success) {
+        std::cout << "Failed to load implicit functions from file." << std::endl;
+        return -1;
+    }
+
+    // compute SDF scalar field for the background mesh
+    std::cout << "computing SDF scalar field for the background mesh..." << std::endl;
+    Eigen::MatrixXd sdf_scalar_field(implicit_functions.functions.size(), background_mesh.vertices.size());
+    for (size_t i = 0; i < background_mesh.vertices.size(); ++i) {
+        const auto& vertex = background_mesh.vertices[i];
+        for (size_t j = 0; j < implicit_functions.functions.size(); ++j) {
+            std::visit([&](auto& function) { sdf_scalar_field(j, i) = function.evaluate_scalar(vertex); },
+                       implicit_functions.functions[j]);
+        }
+    }
+
+    ImplicitSurfaceIntegrator integrator(background_mesh, sdf_scalar_field);
+    labelled_timers_manager   timers_manager{};
+    std::cout << "running implicit surface network integrator..." << std::endl;
+    integrator.run(timers_manager);
+    integrator.compute();
 
-int main(){
     return 0;
 }

application/xmake.lua

@@ -1,6 +1,7 @@
 target("implicit_surface_integrator")
     set_kind("binary")
-    add_rules("config.indirect_predicates.flags")
+    -- add_rules("config.indirect_predicates.flags")
+    add_defines("SHARED_MODULE=0")
     add_deps("implicit_surface_network_process", "implicit_functions", "shared_module")
     add_files("main.cpp")
 target_end()

blobtree_structure/include/utils/eigen_alias.hpp

@@ -0,0 +1,27 @@
+#pragma once
+
+#include <tbb/tbbmalloc_proxy.h>
+
+#define EIGEN_NO_CUDA                     1
+#define EIGEN_USE_MKL                     1
+#define EIGEN_USE_MKL_ALL                 1
+#define EIGEN_USE_THREADS                 1
+// just an experiment, since it seems OKay to mix oneTBB and OpenMP, see
+// https://oneapi-src.github.io/oneTBB/main/tbb_userguide/appendix_B.html
+// #define EIGEN_DONT_PARALLELIZE \
+//     1 // use this macro to disable multi-threading in Eigen, since Eigen's multi-threading is
+//       // restricted to use OpenMP
+#define EIGEN_INITIALIZE_MATRICES_BY_ZERO 1
+#define EIGEN_NO_AUTOMATIC_RESIZING       1
+// #define EIGEN_MAX_CPP_VER                 17
+#include <Eigen/Dense>
+
+// namespace BlobtreeIntegral::IMath
+// {
+// using namespace Eigen;
+
+// static aligned_allocator<double> g_aligned_allocator{};
+
+// static ThreadPool       g_thread_pool(nbThreads());
+// static ThreadPoolDevice g_tensor_thread_pool_handle(&g_thread_pool, nbThreads());
+// } // namespace BlobtreeIntegral::IMath

blobtree_structure/interface/blobtree.h

@@ -0,0 +1,25 @@
+#pragma once
+
+#include "primitive_descriptor.h"
+
+typedef struct Blobtree    blobtree_t;     // fixed complete binary tree of 16 layers
+typedef struct Node        node_t;         // real node in the tree
+typedef struct VirtualNode virtual_node_t; // almost same as node_t, but has parent's and children's pointers to indicate the
+                                           // hierarchy, and it is outside of the tree
+
+EXTERN_C_BEGIN
+
+API blobtree_t* create_blobtree();
+API void        free_blobtree(blobtree_t* blobtree);
+
+API virtual_node_t* blobtree_new_virtual_node(const constant_descriptor_t* desc);
+API void            blobtree_free_virtual_node(virtual_node_t* node);
+
+API void virtual_node_set_parent(virtual_node_t* node, virtual_node_t* parent);
+API void virtual_node_set_left_child(virtual_node_t* node, virtual_node_t* child);
+API void virtual_node_set_right_child(virtual_node_t* node, virtual_node_t* child);
+API void virtual_node_add_child(virtual_node_t* node, virtual_node_t* child);
+API void virtual_node_remove_child(virtual_node_t* node, virtual_node_t* child);
+API void virtual_node_replace_primitive(virtual_node_t* node, const constant_descriptor_t* desc);
+
+EXTERN_C_END

blobtree_structure/interface/primitive_descriptor.h

@@ -0,0 +1,31 @@
+#pragma once
+
+#include <macros.h>
+
+typedef struct _raw_vector3d_t {
+    double x, y, z;
+} raw_vector3d_t;
+
+typedef struct _constant_descriptor_t {
+    double value;
+} constant_descriptor_t;
+
+typedef struct _plane_descriptor_t {
+    raw_vector3d_t point;
+    raw_vector3d_t normal;
+} plane_descriptor_t;
+
+typedef struct _sphere_descriptor_t {
+    raw_vector3d_t center;
+    double         radius;
+} sphere_descriptor_t;
+
+EXTERN_C API double evaluate_constant(constant_descriptor_t* desc, raw_vector3d_t point);
+EXTERN_C API double evaluate_plane(plane_descriptor_t* desc, raw_vector3d_t point);
+EXTERN_C API double evaluate_sphere(sphere_descriptor_t* desc, raw_vector3d_t point);
+
+#define evaluate(desc, point)                       \
+    _Generic((desc),                                \
+        constant_descriptor_t *: evaluate_constant, \
+        plane_descriptor_t *: evaluate_plane,       \
+        sphere_descriptor_t *: evaluate_sphere)(desc, point)

blobtree_structure/xmake.lua

@@ -0,0 +1,10 @@
+add_requires("eigen-latest")
+
+target("blobtree_structure")
+    add_rules("library.shared")
+    add_rules("config.indirect_predicates.flags")
+    add_deps("shared_module")
+    add_includedirs("include")
+    add_includedirs("interface", {public = true})
+    add_files("src/*.cpp")
+target_end()

implicit_arrangements/include/arrangement_builder.hpp

@@ -40,7 +40,7 @@ struct deduce_arrangement_result_type<3> {
 template <size_t N>
 IAComplex<N> init_ia_complex(uint32_t num_planes)
 {
-    IAComplex<N> ia_complex;
+    IAComplex<N> ia_complex{};
     ia_complex.vertices.resize(N + 1);
 
     if constexpr (N == 2) {
@@ -64,7 +64,7 @@ IAComplex<N> init_ia_complex(uint32_t num_planes)
 
         ia_complex.faces.resize(1);
         ia_complex.faces[0].edges = {0, 1, 2};
-        ia_complex.faces[0].signs = {true, true, true};
+        ia_complex.faces[0].signs.resize(4, true);
         ia_complex.faces[0].signs.resize(num_planes, false);
     } else {
         ia_complex.vertices[0] = {1, 2, 3};
@@ -87,10 +87,10 @@ IAComplex<N> init_ia_complex(uint32_t num_planes)
         ia_complex.edges[5].supporting_planes = {0, 1};
 
         ia_complex.faces.resize(4);
-        ia_complex.faces[0].edges            = {5, 3, 4};
-        ia_complex.faces[1].edges            = {2, 1, 5};
-        ia_complex.faces[2].edges            = {4, 0, 2};
-        ia_complex.faces[3].edges            = {1, 0, 3};
+        ia_complex.faces[0].edges            = small_vector_mp<uint32_t>{5, 3, 4};
+        ia_complex.faces[1].edges            = small_vector_mp<uint32_t>{2, 1, 5};
+        ia_complex.faces[2].edges            = small_vector_mp<uint32_t>{4, 0, 2};
+        ia_complex.faces[3].edges            = small_vector_mp<uint32_t>{1, 0, 3};
         ia_complex.faces[0].supporting_plane = 0;
         ia_complex.faces[1].supporting_plane = 1;
         ia_complex.faces[2].supporting_plane = 2;
@@ -105,8 +105,8 @@ IAComplex<N> init_ia_complex(uint32_t num_planes)
         ia_complex.faces[3].negative_cell    = INVALID_INDEX;
 
         ia_complex.cells.resize(1);
-        ia_complex.cells[0].faces = {0, 1, 2, 3};
-        ia_complex.cells[0].signs = {true, true, true, true};
+        ia_complex.cells[0].faces = small_vector_mp<uint32_t>{0, 1, 2, 3};
+        ia_complex.cells[0].signs.resize(4, true);
         ia_complex.cells[0].signs.resize(num_planes, false);
     }
     return ia_complex;
@@ -121,14 +121,15 @@ public:
     using arrangement_result_type = typename detail::deduce_arrangement_result_type<N>::type;
 
     ArrangementBuilder(const plane_type* planes, const uint32_t num_planes)
-        : m_planes(planes, planes + num_planes), m_coplanar_planes(num_planes + N + 1)
     {
         const auto* ia = lookup(planes, num_planes);
         if (ia != nullptr) {
             m_arrangement.arrangement         = *ia;
             m_arrangement.is_runtime_computed = false;
         } else {
-            auto     ia_complex         = init_ia_complex<N>(num_planes);
+            auto ia_complex = init_ia_complex<N>(num_planes + static_cast<uint32_t>(N) + 1);
+            m_planes        = PlaneGroup<N>(planes, planes + num_planes);
+            m_coplanar_planes.init(num_planes + static_cast<uint32_t>(N) + 1);
             uint32_t unique_plane_count = 0;
             for (size_t i = 0; i < num_planes; i++) {
                 auto plane_id          = static_cast<uint32_t>(i + N + 1);
@@ -139,6 +140,7 @@ public:
                     unique_plane_count++;
                 }
             }
+            std::cout << "test" << std::endl;
 
             m_arrangement.arrangement         = extract_arrangement(std::move(ia_complex));
             m_arrangement.is_runtime_computed = true;

implicit_arrangements/include/plane.hpp

@@ -27,17 +27,19 @@ template <size_t N>
 struct PlaneGroup {
     using PlaneType = typename detail::deduce_plane_type<N>::type;
 
+    PlaneGroup() = default;
+
     template <typename InputIt>
     PlaneGroup(InputIt first, InputIt last)
     {
         if constexpr (N == 2) {
-            planes = {
+            planes = small_vector_mp<PlaneType, N + 1>{
                 {1, 0, 0},
                 {0, 1, 0},
                 {0, 0, 1}
             };
         } else {
-            planes = {
+            planes = small_vector_mp<PlaneType, N + 1>{
                 {1, 0, 0, 0},
                 {0, 1, 0, 0},
                 {0, 0, 1, 0},

implicit_arrangements/src/add_plane.cpp

@@ -1,4 +1,5 @@
 #include <algorithm/glue_algorithm.hpp>
+#include <iostream>
 
 #include "ia_cut_face.hpp"
 #include "robust_assert.hpp"
@@ -51,19 +52,20 @@ inline void remove_unused_geometry(IAComplex<N>& data)
 
     // Shrink edges.
     {
+        active_geometries.clear();
         active_geometries.resize(data.edges.size(), false);
         for (auto& f : data.faces) {
             for (auto eid : f.edges) { active_geometries[eid] = true; }
         }
 
-        shrink(data.faces, index_map, active_geometries);
+        shrink(data.edges, index_map, active_geometries);
 
         for (auto& f : data.faces) {
             algorithm::transform<algorithm::ExecutionPolicySelector::simd_only>(
                 f.edges.begin(),
                 f.edges.end(),
                 f.edges.begin(),
-                [&](size_t i) {
+                [&](uint32_t i) {
                     ROBUST_ASSERT(index_map[i] != INVALID_INDEX);
                     return index_map[i];
                 });
@@ -72,6 +74,7 @@ inline void remove_unused_geometry(IAComplex<N>& data)
 
     // Shrink vertices.
     {
+        active_geometries.clear();
         active_geometries.resize(data.vertices.size(), false);
         for (auto& e : data.edges) {
             for (auto vid : e.vertices) {
@@ -80,7 +83,7 @@ inline void remove_unused_geometry(IAComplex<N>& data)
             }
         }
 
-        shrink(data.faces, index_map, active_geometries);
+        shrink(data.vertices, index_map, active_geometries);
 
         for (auto& e : data.edges) {
             e.vertices[0] = index_map[e.vertices[0]];
@@ -189,14 +192,14 @@ uint32_t add_plane(const PlaneGroup<3>& repo, IAComplex<3>& ia_complex, uint32_t
     }
 
     // Step 3: handle 2-faces.
-    small_vector_mp<std::array<uint32_t, 3>> subfaces;
+    small_vector_mp<std::array<uint32_t, 3>> subfaces{};
     subfaces.reserve(num_faces);
     for (uint32_t i = 0; i < num_faces; i++) {
         subfaces.emplace_back(ia_cut_2_face(ia_complex, i, plane_index, orientations.data(), subedges.data()));
     }
 
     // Step 4: handle 3-faces.
-    small_vector_mp<std::array<uint32_t, 3>> subcells;
+    small_vector_mp<std::array<uint32_t, 3>> subcells{};
     subcells.reserve(num_cells);
     for (uint32_t i = 0; i < num_cells; i++) {
         subcells.emplace_back(ia_cut_3_face(ia_complex, i, plane_index, subfaces.data()));
@@ -206,8 +209,8 @@ uint32_t add_plane(const PlaneGroup<3>& repo, IAComplex<3>& ia_complex, uint32_t
     {
         small_dynamic_bitset_mp<> to_keep(cells.size(), false);
         for (const auto& subcell : subcells) {
-            to_keep[subcell[0]] = subcell[0] != INVALID_INDEX;
-            to_keep[subcell[1]] = subcell[1] != INVALID_INDEX;
+            if (subcell[0] != INVALID_INDEX) to_keep[subcell[0]] = true;
+            if (subcell[1] != INVALID_INDEX) to_keep[subcell[1]] = true;
         }
 
         small_vector_mp<uint32_t> index_map{};
@@ -221,8 +224,8 @@ uint32_t add_plane(const PlaneGroup<3>& repo, IAComplex<3>& ia_complex, uint32_t
     }
 
     // Step 6: check for coplanar planes.
-    size_t coplanar_plane = INVALID_INDEX;
-    for (size_t i = 0; i < num_faces; i++) {
+    uint32_t coplanar_plane = INVALID_INDEX;
+    for (uint32_t i = 0; i < num_faces; i++) {
         const auto& subface = subfaces[i];
         if (subface[0] == INVALID_INDEX && subface[1] == INVALID_INDEX) {
             const auto& f  = faces[i];
@@ -232,6 +235,7 @@ uint32_t add_plane(const PlaneGroup<3>& repo, IAComplex<3>& ia_complex, uint32_t
 
     // Step 7: remove unused geometries.
     remove_unused_geometry(ia_complex);
+    std::cout << "test" << std::endl;
 
     return coplanar_plane;
 }

implicit_functions/include/implicit_function.hpp

@@ -2,15 +2,11 @@
 
 #include "utils/eigen_alias.hpp"
 
+template <typename Scalar, size_t Dim>
 class ImplicitFunction
 {
 public:
-    template <size_t Dim>
-    double evaluate_scalar(const Eigen::Ref<const Eigen::Vector<double, Dim>>& pos) const;
-
-    template <size_t Dim>
-    Eigen::Vector<double, Dim> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<double, Dim>>& pos) const;
-
-    template <size_t Dim>
-    Eigen::Vector<double, Dim + 1> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<double, Dim>>& pos) const;
+    Scalar                         evaluate_scalar(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>>& pos) const;
+    Eigen::Vector<Scalar, Dim>     evaluate_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>>& pos) const;
+    Eigen::Vector<Scalar, Dim + 1> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>>& pos) const;
 };

implicit_functions/include/implicit_functions.hpp

@@ -0,0 +1,12 @@
+#pragma once
+
+#include <variant>
+
+#include <primitive_functions.hpp>
+
+template <typename Scalar, size_t Dim>
+using implicit_function_t = std::variant<ConstantFunction<Scalar, Dim>,
+                                         PlaneDistanceFunction<Scalar, Dim>,
+                                         CylinderDistanceFunction<Scalar, Dim>,
+                                         SphereDistanceFunction<Scalar, Dim>,
+                                         ConeDistanceFunction<Scalar, Dim>>;

implicit_functions/include/load_functions.hpp

@@ -1,14 +1,14 @@
 #pragma once
 
-#include <string_view>
+#include <string>
 
 #include <container/small_vector.hpp>
 
-#include "implicit_function.hpp"
+#include "implicit_functions.hpp"
 
 struct load_functions_result_t {
-    small_vector_mp<std::unique_ptr<ImplicitFunction>> functions{};
-    bool                                               success{};
+    small_vector_mp<implicit_function_t<double, 3>> functions{};
+    bool                                            success{};
 };
 
-load_functions_result_t load_functions(const std::string_view& filename);
+load_functions_result_t load_functions(const std::string& filename);

implicit_functions/include/primitive_functions.hpp

@@ -4,94 +4,89 @@
 
 static constexpr double kEpsilon = 1e-6;
 
-class ConstantFunction : public ImplicitFunction
+template <typename Scalar, size_t Dim>
+class ConstantFunction : public ImplicitFunction<Scalar, Dim>
 {
 public:
-    explicit ConstantFunction(double value) : value_(value) {}
+    explicit ConstantFunction(Scalar value) : value_(value) {}
 
-    template <size_t Dim>
-    double evaluate_scalar(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const
-    {
-        return value_;
-    }
+    Scalar evaluate_scalar(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const { return value_; }
 
-    template <size_t Dim>
-    Eigen::Vector<double, Dim> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const
+    Eigen::Vector<Scalar, Dim> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const
     {
-        return Eigen::Vector<double, Dim>::Zero();
+        return Eigen::Vector<Scalar, Dim>::Zero();
     }
 
-    template <size_t Dim>
-    Eigen::Vector<double, Dim + 1> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const
+    Eigen::Vector<Scalar, Dim + 1> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const
     {
-        auto res = Eigen::Vector<double, Dim + 1>::Zero();
+        auto res = Eigen::Vector<Scalar, Dim + 1>::Zero();
         res[0]   = value_;
         return res;
     }
 
 private:
-    double value_{};
+    Scalar value_{};
 };
 
-template <size_t Dim>
-class PlaneDistanceFunction : public ImplicitFunction
+template <typename Scalar, size_t Dim>
+class PlaneDistanceFunction : public ImplicitFunction<Scalar, Dim>
 {
 public:
-    PlaneDistanceFunction(const Eigen::Ref<const Eigen::Vector<double, Dim>> &point,
-                          const Eigen::Ref<const Eigen::Vector<double, Dim>> &normal)
+    PlaneDistanceFunction(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &point,
+                          const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &normal)
         : point_(point), normal_(normal)
     {
         normal_.normalize();
     }
 
-    double evaluate_scalar(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const { return normal_.dot(pos - point_); }
+    Scalar evaluate_scalar(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const { return normal_.dot(pos - point_); }
 
-    Eigen::Vector<double, Dim> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const
+    Eigen::Vector<Scalar, Dim> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const
     {
         return normal_;
     }
 
-    Eigen::Vector<double, Dim + 1> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const
+    Eigen::Vector<Scalar, Dim + 1> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const
     {
-        Eigen::Vector<double, Dim + 1> res{};
+        Eigen::Vector<Scalar, Dim + 1> res{};
         res.template head<Dim>().array() = normal_.array();
         res[Dim]                         = evaluate_scalar(pos);
         return res;
     }
 
 private:
-    Eigen::Vector<double, Dim> point_{};
-    Eigen::Vector<double, Dim> normal_{};
+    Eigen::Vector<Scalar, Dim> point_{};
+    Eigen::Vector<Scalar, Dim> normal_{};
 };
 
-template <size_t Dim>
-class SphereDistanceFunction : public ImplicitFunction
+template <typename Scalar, size_t Dim>
+class SphereDistanceFunction : public ImplicitFunction<Scalar, Dim>
 {
 public:
-    SphereDistanceFunction(const Eigen::Ref<const Eigen::Vector<double, Dim>> &center, double radius)
+    SphereDistanceFunction(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &center, Scalar radius)
         : center_(center), radius_(radius)
     {
     }
 
-    double evaluate_scalar(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const
+    Scalar evaluate_scalar(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const
     {
         return (pos - center_).norm() - radius_;
     }
 
-    Eigen::Vector<double, Dim> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const
+    Eigen::Vector<Scalar, Dim> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const
     {
         auto   vec  = pos - center_;
-        double norm = vec.norm();
-        if (norm <= kEpsilon) return Eigen::Vector<double, Dim>::Zero();
+        Scalar norm = vec.norm();
+        if (norm <= kEpsilon) return Eigen::Vector<Scalar, Dim>::Zero();
         return vec / norm;
     }
 
-    Eigen::Vector<double, Dim + 1> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<double, Dim>> &pos) const
+    Eigen::Vector<Scalar, Dim + 1> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, Dim>> &pos) const
     {
-        Eigen::Vector<double, Dim + 1> res{};
+        Eigen::Vector<Scalar, Dim + 1> res{};
 
         auto   vec  = pos - center_;
-        double norm = vec.norm();
+        Scalar norm = vec.norm();
         if (norm <= kEpsilon) {
             res[Dim] = -radius_;
             return res;
@@ -103,52 +98,52 @@ public:
     }
 
 private:
-    Eigen::Vector<double, Dim> center_{};
-    double                     radius_{};
+    Eigen::Vector<Scalar, Dim> center_{};
+    Scalar                     radius_{};
 };
 
-template <size_t Dim>
+template <typename Scalar, size_t Dim>
 class CylinderDistanceFunction;
 
-template <>
-class CylinderDistanceFunction<3> : public ImplicitFunction
+template <typename Scalar>
+class CylinderDistanceFunction<Scalar, 3> : public ImplicitFunction<Scalar, 3>
 {
 public:
-    CylinderDistanceFunction(const Eigen::Ref<const Eigen::Vector3d> &axis_point,
-                             const Eigen::Ref<const Eigen::Vector3d> &axis_dir,
-                             double                                   radius)
+    CylinderDistanceFunction(const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &axis_point,
+                             const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &axis_dir,
+                             Scalar                                            radius)
         : axis_point_(axis_point), axis_dir_(axis_dir), radius_(radius)
     {
         axis_dir_.normalize();
     }
 
-    double evaluate_scalar(const Eigen::Ref<const Eigen::Vector3d> &pos) const
+    Scalar evaluate_scalar(const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &pos) const
     {
         auto vec  = pos - axis_point_;
         auto dist = vec.dot(axis_dir_);
         return (vec - dist * axis_dir_).norm() - radius_;
     }
 
-    Eigen::Vector3d evaluate_gradient(const Eigen::Ref<const Eigen::Vector3d> &pos) const
+    Eigen::Vector<Scalar, 3> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &pos) const
     {
-        auto            vec      = pos - axis_point_;
-        auto            dist     = vec.dot(axis_dir_);
-        auto            vec_para = dist * axis_dir_;
-        Eigen::Vector3d vec_perp = vec - vec_para;
-        double          norm     = vec_perp.norm();
-        if (norm <= kEpsilon) return Eigen::Vector3d::Zero();
+        auto                     vec      = pos - axis_point_;
+        auto                     dist     = vec.dot(axis_dir_);
+        auto                     vec_para = dist * axis_dir_;
+        Eigen::Vector<Scalar, 3> vec_perp = vec - vec_para;
+        Scalar                   norm     = vec_perp.norm();
+        if (norm <= kEpsilon) return Eigen::Vector<Scalar, 3>::Zero();
         return vec_perp / norm;
     }
 
-    Eigen::Vector4d evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector3d> &pos) const
+    Eigen::Vector<Scalar, 4> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &pos) const
     {
-        Eigen::Vector4d res{};
+        Eigen::Vector<Scalar, 4> res{};
 
-        auto            vec      = pos - axis_point_;
-        auto            dist     = vec.dot(axis_dir_);
-        auto            vec_para = dist * axis_dir_;
-        Eigen::Vector3d vec_perp = vec - vec_para;
-        double          norm     = vec_perp.norm();
+        auto                     vec      = pos - axis_point_;
+        auto                     dist     = vec.dot(axis_dir_);
+        auto                     vec_para = dist * axis_dir_;
+        Eigen::Vector<Scalar, 3> vec_perp = vec - vec_para;
+        Scalar                   norm     = vec_perp.norm();
         if (norm <= kEpsilon) {
             res[3] = -radius_;
             return res;
@@ -160,56 +155,56 @@ public:
     }
 
 private:
-    Eigen::Vector3d axis_point_{};
-    Eigen::Vector3d axis_dir_{};
-    double          radius_{};
+    Eigen::Vector<Scalar, 3> axis_point_{};
+    Eigen::Vector<Scalar, 3> axis_dir_{};
+    Scalar                   radius_{};
 };
 
-template <size_t Dim>
+template <typename Scalar, size_t Dim>
 class ConeDistanceFunction;
 
-template <>
-class ConeDistanceFunction<3> : public ImplicitFunction
+template <typename Scalar>
+class ConeDistanceFunction<Scalar, 3> : public ImplicitFunction<Scalar, 3>
 {
 public:
-    ConeDistanceFunction(const Eigen::Ref<const Eigen::Vector3d> &apex_point,
-                         const Eigen::Ref<const Eigen::Vector3d> &axis_dir,
-                         double                                   apex_angle)
+    ConeDistanceFunction(const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &apex_point,
+                         const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &axis_dir,
+                         Scalar                                            apex_angle)
         : apex_point_(apex_point), axis_dir_(axis_dir), apex_angle_cosine_(std::cos(apex_angle))
     {
         axis_dir_.normalize();
     }
 
-    double evaluate_scalar(const Eigen::Ref<const Eigen::Vector3d> &pos) const
+    Scalar evaluate_scalar(const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &pos) const
     {
         return apex_angle_cosine_ * (pos - apex_point_).norm() - (pos - apex_point_).dot(axis_dir_);
     }
 
-    Eigen::Vector3d evaluate_gradient(const Eigen::Ref<const Eigen::Vector3d> &pos) const
+    Eigen::Vector<Scalar, 3> evaluate_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &pos) const
     {
         auto vec_apex = pos - apex_point_;
         auto vec_norm = vec_apex.norm();
-        if (vec_norm <= kEpsilon) return Eigen::Vector3d::Zero();
+        if (vec_norm <= kEpsilon) return Eigen::Vector<Scalar, 3>::Zero();
         return vec_apex * apex_angle_cosine_ / vec_norm - axis_dir_;
     }
 
-    Eigen::Vector4d evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector3d> &pos) const
+    Eigen::Vector<Scalar, 4> evaluate_scalar_gradient(const Eigen::Ref<const Eigen::Vector<Scalar, 3>> &pos) const
     {
         auto vec_apex = pos - apex_point_;
         auto vec_norm = vec_apex.norm();
-        if (vec_norm <= kEpsilon) return Eigen::Vector4d::Zero();
+        if (vec_norm <= kEpsilon) return Eigen::Vector<Scalar, 4>::Zero();
 
-        Eigen::Vector4d res{};
+        Eigen::Vector<Scalar, 4> res{};
         res.template head<3>().array() = vec_apex.array() * apex_angle_cosine_ / vec_norm - axis_dir_.array();
         res[3]                         = apex_angle_cosine_ * vec_norm - vec_apex.dot(axis_dir_);
         return res;
     }
 
 private:
-    Eigen::Vector3d apex_point_{};
-    Eigen::Vector3d axis_dir_{};
-    double          apex_angle_cosine_{};
+    Eigen::Vector<Scalar, 3> apex_point_{};
+    Eigen::Vector<Scalar, 3> axis_dir_{};
+    Scalar                   apex_angle_cosine_{};
 };
 
-template <size_t Dim>
+template <typename Scalar, size_t Dim>
 class TorusDistanceFunction;

implicit_functions/interface/implicit_functions.hpp

@@ -1,4 +0,0 @@
-#pragma once
-
-#include <primitive_functions.hpp>
-#include <load_functions.hpp>

implicit_functions/src/main.cpp

@@ -3,22 +3,21 @@
 
 #include <nlohmann/json.hpp>
 
-#include <implicit_functions.hpp>
-#include "primitive_functions.hpp"
+#include <load_functions.hpp>
 
-load_functions_result_t load_functions(const std::string_view& filename)
+load_functions_result_t load_functions(const std::string& filename)
 {
     using json = nlohmann::json;
 
-    std::ifstream file(filename.data());
-    if (!file) {
+    std::ifstream fin(filename.c_str());
+    if (!fin) {
         std::cout << "Failed to open function file: " << filename << std::endl;
         return {};
     }
 
     json data;
-    file >> data;
-    file.close();
+    fin >> data;
+    fin.close();
 
     load_functions_result_t result{};
     const auto              num_functions = data.size();
@@ -26,23 +25,23 @@ load_functions_result_t load_functions(const std::string_view& filename)
     for (size_t j = 0; j < num_functions; ++j) {
         const auto type = data[j]["type"].get<std::string>();
         if (type == "plane") {
-            Eigen::Vector3d point, normal;
+            Eigen::Vector3d point{}, normal{};
             for (auto i = 0; i < 3; ++i) {
                 point[i]  = data[j]["point"][i].get<double>();
                 normal[i] = data[j]["normal"][i].get<double>();
             }
-            result.functions.emplace_back(std::make_unique<PlaneDistanceFunction<3>>(point, normal));
+            result.functions.emplace_back(PlaneDistanceFunction<double, 3>{point, normal});
         } else if (type == "line") {
-            Eigen::Vector3d point, direction;
+            Eigen::Vector3d point{}, direction{};
             for (auto i = 0; i < 3; ++i) {
                 point[i]     = data[j]["point"][i].get<double>();
                 direction[i] = data[j]["direction"][i].get<double>();
             }
-            result.functions.emplace_back(std::make_unique<CylinderDistanceFunction<3>>(point, direction, 0));
+            result.functions.emplace_back(CylinderDistanceFunction<double, 3>{point, direction, 0});
         } else if (type == "cylinder") {
-            Eigen::Vector3d axis_point, axis_direction;
+            Eigen::Vector3d axis_point{}, axis_direction{};
             if (data[j].contains("axis_point1") && data[j].contains("axis_point2")) {
-                Eigen::Vector3d axis_point_;
+                Eigen::Vector3d axis_point_{};
                 for (auto i = 0; i < 3; ++i) {
                     axis_point[i]  = data[j]["axis_point1"][i].get<double>();
                     axis_point_[i] = data[j]["axis_point2"][i].get<double>();
@@ -55,28 +54,29 @@ load_functions_result_t load_functions(const std::string_view& filename)
                 }
             }
             double radius = data[j]["radius"].get<double>();
-            result.functions.emplace_back(std::make_unique<CylinderDistanceFunction<3>>(axis_point, axis_direction, radius));
+            result.functions.emplace_back(CylinderDistanceFunction<double, 3>{axis_point, axis_direction, radius});
         } else if (type == "sphere") {
-            Eigen::Vector3d center;
+            Eigen::Vector3d center{};
             for (auto i = 0; i < 3; ++i) { center[i] = data[j]["center"][i].get<double>(); }
             double radius = data[j]["radius"].get<double>();
-            result.functions.emplace_back(std::make_unique<SphereDistanceFunction<3>>(center, radius));
+            result.functions.emplace_back(SphereDistanceFunction<double, 3>{center, radius});
         } else if (type == "cone") {
-            Eigen::Vector3d apex_point, axis_direction;
+            Eigen::Vector3d apex_point{}, axis_direction{};
             for (auto i = 0; i < 3; ++i) {
                 apex_point[i]     = data[j]["apex_point"][i].get<double>();
                 axis_direction[i] = data[j]["axis_direction"][i].get<double>();
             }
             double apex_angle = data[j]["apex_angle"].get<double>();
-            result.functions.emplace_back(std::make_unique<ConeDistanceFunction<3>>(apex_point, axis_direction, apex_angle));
+            result.functions.emplace_back(ConeDistanceFunction<double, 3>{apex_point, axis_direction, apex_angle});
         } else if (type == "constant") {
             double value = data[j]["value"].get<double>();
-            result.functions.emplace_back(std::make_unique<ConstantFunction>(value));
+            result.functions.emplace_back(ConstantFunction<double, 3>{value});
         } else {
             std::cout << "Unsupported function type: " << type << std::endl;
             return {};
         }
     }
 
+    result.success = true;
     return result;
 }

implicit_functions/xmake.lua

@@ -7,15 +7,13 @@ target("implicit_functions")
     add_vectorexts("neon")
     add_vectorexts("avx")
     add_vectorexts("avx2")
-    add_vectorexts("avx512")
     add_vectorexts("sse")
     add_vectorexts("sse2")
     add_vectorexts("sse3")
     add_vectorexts("ssse3")
     add_vectorexts("sse4.2")
     add_defines("SHARED_MODULE=0")
-    add_includedirs("include")
-    add_includedirs("interface", {public = true})
+    add_includedirs("include", {public = true})
     add_files("src/*.cpp")
     add_packages("eigen-latest", {public = true})
     add_packages("nlohmann_json")

network_process/include/background_mesh.hpp

@@ -2,6 +2,6 @@
 
 #include "utils/fwd_types.hpp"
 
-tetrahedron_mesh_t generate_tetrahedron_background_mesh(size_t             resolution,
-                                                        const raw_point_t& aabb_min,
-                                                        const raw_point_t& aabb_max);
+tetrahedron_mesh_t generate_tetrahedron_background_mesh(uint32_t                             resolution,
+                                                        const Eigen::Ref<const raw_point_t>& aabb_min,
+                                                        const Eigen::Ref<const raw_point_t>& aabb_max);

network_process/src/background_mesh.cpp

@@ -4,30 +4,30 @@
 #include <algorithm/glue_algorithm.hpp>
 
 struct TetrahedronVertexIndexGroup {
-    size_t v00, v01, v02, v03;
-    size_t v10, v11, v12, v13;
-    size_t v20, v21, v22, v23;
-    size_t v30, v31, v32, v33;
-    size_t v40, v41, v42, v43;
+    uint32_t v00, v01, v02, v03;
+    uint32_t v10, v11, v12, v13;
+    uint32_t v20, v21, v22, v23;
+    uint32_t v30, v31, v32, v33;
+    uint32_t v40, v41, v42, v43;
 };
 
-tetrahedron_mesh_t generate_tetrahedron_background_mesh(size_t             resolution,
-                                                        const raw_point_t& aabb_min,
-                                                        const raw_point_t& aabb_max)
+tetrahedron_mesh_t generate_tetrahedron_background_mesh(uint32_t                             resolution,
+                                                        const Eigen::Ref<const raw_point_t>& aabb_min,
+                                                        const Eigen::Ref<const raw_point_t>& aabb_max)
 {
     assert(resolution > 0);
 
     const auto         N = resolution + 1;
-    tetrahedron_mesh_t mesh;
+    tetrahedron_mesh_t mesh{};
     mesh.vertices.resize(N * N * N);
     mesh.indices.resize(resolution * resolution * resolution * 5);
 
-    for (auto i = size_t{0}; i < N; ++i) {
-        const auto x = (aabb_max.x() - aabb_min.x()) * i / resolution + aabb_min.x();
-        for (auto j = size_t{0}; j < N; ++j) {
-            const auto y = (aabb_max.y() - aabb_min.y()) * j / resolution + aabb_min.y();
-            algorithm::for_loop<algorithm::ExecutionPolicySelector::simd_only>(size_t{0}, N, [&](size_t k) {
-                const auto z      = (aabb_max.z() - aabb_min.z()) * k / resolution + aabb_min.z();
+    for (uint32_t i = 0; i < N; ++i) {
+        const auto x = (aabb_max[0] - aabb_min[0]) * i / resolution + aabb_min[0];
+        for (uint32_t j = 0; j < N; ++j) {
+            const auto y = (aabb_max[1] - aabb_min[1]) * j / resolution + aabb_min[1];
+            algorithm::for_loop<algorithm::ExecutionPolicySelector::simd_only>(0u, N, [&](uint32_t k) {
+                const auto z      = (aabb_max[2] - aabb_min[2]) * k / resolution + aabb_min[2];
                 const auto v0     = i * N * N + j * N + k;
                 mesh.vertices[v0] = {x, y, z};
 

network_process/src/implicit_surface_network_processor.cpp

@@ -12,11 +12,11 @@ bool ImplicitSurfaceNetworkProcessor::run(labelled_timers_manager& timers_manage
 
     const auto num_vert  = background_mesh.vertices.size();
     const auto num_tets  = background_mesh.indices.size();
-    const auto num_funcs = sdf_scalar_field.cols();
+    const auto num_funcs = sdf_scalar_field.rows();
 
     // compute function signs at vertices
     auto                scalar_field_sign = [](double x) -> int8_t { return (x > 0) ? 1 : ((x < 0) ? -1 : 0); };
-    Eigen::MatrixXd     scalar_field_signs(sdf_scalar_field.rows(), sdf_scalar_field.cols());
+    Eigen::MatrixXd     scalar_field_signs(num_funcs, num_vert);
     dynamic_bitset_mp<> is_degenerate_vertex(num_vert, false);
     {
         timers_manager.push_timer("identify sdf signs");
@@ -38,15 +38,14 @@ bool ImplicitSurfaceNetworkProcessor::run(labelled_timers_manager& timers_manage
         func_in_tet.reserve(num_tets);
         start_index_of_tet.reserve(num_tets + 1);
         start_index_of_tet.emplace_back(0);
-        for (Eigen::Index i = 0; i < num_funcs; ++i) {
-            for (const auto& func_signs : scalar_field_signs.colwise()) {
-                uint32_t   pos_count{};
-                const auto tet_ptr = &background_mesh.indices[i].v1;
-                for (uint32_t j = 0; j < 4; ++j) {
-                    if (func_signs(tet_ptr[j]) == 1) pos_count++;
-                }
+        for (Eigen::Index i = 0; i < num_tets; ++i) {
+            const auto tet_ptr = &background_mesh.indices[i].v1;
+            for (Eigen::Index j = 0; j < num_funcs; ++j) {
+                uint32_t pos_count{};
+                for (uint32_t k = 0; k < 4; ++k)
+                    if (scalar_field_signs(j, tet_ptr[k]) == 1) pos_count++;
                 // tets[i].size() == 4, this means that the function is active in this tet
-                if (0 < pos_count && pos_count < 4) func_in_tet.emplace_back(i);
+                if (0 < pos_count && pos_count < 4) func_in_tet.emplace_back(j);
             }
             if (func_in_tet.size() > start_index_of_tet.back()) { ++num_intersecting_tet; }
             start_index_of_tet.emplace_back(static_cast<uint32_t>(func_in_tet.size()));
@@ -72,7 +71,7 @@ bool ImplicitSurfaceNetworkProcessor::run(labelled_timers_manager& timers_manage
             for (uint32_t i = 0; i < num_tets; ++i) {
                 const auto start_index              = start_index_of_tet[i];
                 const auto active_funcs_in_curr_tet = start_index_of_tet[i + 1] - start_index;
-                if (num_funcs == 0) {
+                if (active_funcs_in_curr_tet == 0) {
                     cut_result_index.emplace_back(invalid_index);
                     continue;
                 }
@@ -115,6 +114,7 @@ bool ImplicitSurfaceNetworkProcessor::run(labelled_timers_manager& timers_manage
         }
         timers_manager.pop_timer("implicit arrangements calculation in total");
     }
+    std::cout << "test" << std::endl;
 
     // extract arrangement mesh: combining results from all tets to produce a mesh
     stl_vector_mp<IsoVertex> iso_verts{};

network_process/src/patch_connectivity.cpp

@@ -127,7 +127,7 @@ void compute_chains(const stl_vector_mp<IsoEdge>&                   patch_edges,
 
 void compute_shells_and_components(uint32_t                                                 num_patch,
                                    const stl_vector_mp<small_vector_mp<half_patch_pair_t>>& half_patch_pair_list,
-                                   stl_vector_mp<stl_vector_mp<uint32_t>>&                  shells,
+                                   stl_vector_mp<small_vector_mp<uint32_t>>&                shells,
                                    stl_vector_mp<uint32_t>&                                 shell_of_half_patch,
                                    stl_vector_mp<small_vector_mp<uint32_t>>&                components,
                                    stl_vector_mp<uint32_t>&                                 component_of_patch)

network_process/xmake.lua

@@ -3,7 +3,7 @@ add_requires("eigen-latest")
 target("implicit_surface_network_process")
     set_kind("static")
     add_defines("SHARED_MODULE=0")
-    add_rules("config.indirect_predicates.flags")
+    -- add_rules("config.indirect_predicates.flags")
     add_includedirs("include", {public = true})
     add_files("src/*.cpp")
     add_packages("eigen-latest", {public = true})

shared_module/container/detail/small_vector_base.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include <assert.h>
 #include <optional>
 #include <initializer_list>
 

shared_module/container/static_vector.hpp

@@ -406,7 +406,7 @@ public:
 
     void resize(size_type count)
     {
-        assert(count < max_size());
+        assert(count <= max_size());
 
         if (count > m_size)
             std::uninitialized_default_construct(end(), begin() + count);
@@ -418,7 +418,7 @@ public:
 
     void resize(size_type count, const T& value)
     {
-        assert(count < max_size());
+        assert(count <= max_size());
 
         if (count > m_size)
             std::uninitialized_fill(end(), begin() + count, value);