diff --git a/.vscode/tasks.json b/.vscode/tasks.json index 8de932d..97f30e2 100644 --- a/.vscode/tasks.json +++ b/.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": [] } ] } \ No newline at end of file diff --git a/3rdparty/xmake.lua b/3rdparty/xmake.lua index 5e52ac9..596ddbe 100644 --- a/3rdparty/xmake.lua +++ b/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") diff --git a/application/implicit_surface_integrator.hpp b/application/implicit_surface_integrator.hpp new file mode 100644 index 0000000..b42fec0 --- /dev/null +++ b/application/implicit_surface_integrator.hpp @@ -0,0 +1,41 @@ +#pragma once + +#include + +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 m_per_face_surf_int{}; + stl_vector_mp m_per_face_vol_int{}; +}; \ No newline at end of file diff --git a/application/main.cpp b/application/main.cpp index f60ad92..01cb6f4 100644 --- a/application/main.cpp +++ b/application/main.cpp @@ -1,6 +1,37 @@ -#include -#include +#include +#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; } \ No newline at end of file diff --git a/application/xmake.lua b/application/xmake.lua index 9bfeddd..8bda62e 100644 --- a/application/xmake.lua +++ b/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() \ No newline at end of file diff --git a/blobtree_structure/include/utils/eigen_alias.hpp b/blobtree_structure/include/utils/eigen_alias.hpp new file mode 100644 index 0000000..7e4bced --- /dev/null +++ b/blobtree_structure/include/utils/eigen_alias.hpp @@ -0,0 +1,27 @@ +#pragma once + +#include + +#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 + +// namespace BlobtreeIntegral::IMath +// { +// using namespace Eigen; + +// static aligned_allocator g_aligned_allocator{}; + +// static ThreadPool g_thread_pool(nbThreads()); +// static ThreadPoolDevice g_tensor_thread_pool_handle(&g_thread_pool, nbThreads()); +// } // namespace BlobtreeIntegral::IMath \ No newline at end of file diff --git a/blobtree_structure/interface/blobtree.h b/blobtree_structure/interface/blobtree.h new file mode 100644 index 0000000..644b533 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/blobtree_structure/interface/primitive_descriptor.h b/blobtree_structure/interface/primitive_descriptor.h new file mode 100644 index 0000000..ef432d7 --- /dev/null +++ b/blobtree_structure/interface/primitive_descriptor.h @@ -0,0 +1,31 @@ +#pragma once + +#include + +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) \ No newline at end of file diff --git a/blobtree_structure/xmake.lua b/blobtree_structure/xmake.lua new file mode 100644 index 0000000..2e5fd8d --- /dev/null +++ b/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() \ No newline at end of file diff --git a/implicit_arrangements/include/arrangement_builder.hpp b/implicit_arrangements/include/arrangement_builder.hpp index e410dee..d73e57e 100644 --- a/implicit_arrangements/include/arrangement_builder.hpp +++ b/implicit_arrangements/include/arrangement_builder.hpp @@ -40,7 +40,7 @@ struct deduce_arrangement_result_type<3> { template IAComplex init_ia_complex(uint32_t num_planes) { - IAComplex ia_complex; + IAComplex ia_complex{}; ia_complex.vertices.resize(N + 1); if constexpr (N == 2) { @@ -64,7 +64,7 @@ IAComplex 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 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{5, 3, 4}; + ia_complex.faces[1].edges = small_vector_mp{2, 1, 5}; + ia_complex.faces[2].edges = small_vector_mp{4, 0, 2}; + ia_complex.faces[3].edges = small_vector_mp{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 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{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::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(num_planes); + auto ia_complex = init_ia_complex(num_planes + static_cast(N) + 1); + m_planes = PlaneGroup(planes, planes + num_planes); + m_coplanar_planes.init(num_planes + static_cast(N) + 1); uint32_t unique_plane_count = 0; for (size_t i = 0; i < num_planes; i++) { auto plane_id = static_cast(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; diff --git a/implicit_arrangements/include/plane.hpp b/implicit_arrangements/include/plane.hpp index 09a611b..9933084 100644 --- a/implicit_arrangements/include/plane.hpp +++ b/implicit_arrangements/include/plane.hpp @@ -27,17 +27,19 @@ template struct PlaneGroup { using PlaneType = typename detail::deduce_plane_type::type; + PlaneGroup() = default; + template PlaneGroup(InputIt first, InputIt last) { if constexpr (N == 2) { - planes = { + planes = small_vector_mp{ {1, 0, 0}, {0, 1, 0}, {0, 0, 1} }; } else { - planes = { + planes = small_vector_mp{ {1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 0}, diff --git a/implicit_arrangements/src/add_plane.cpp b/implicit_arrangements/src/add_plane.cpp index a4e3135..389eb73 100644 --- a/implicit_arrangements/src/add_plane.cpp +++ b/implicit_arrangements/src/add_plane.cpp @@ -1,4 +1,5 @@ #include +#include #include "ia_cut_face.hpp" #include "robust_assert.hpp" @@ -51,19 +52,20 @@ inline void remove_unused_geometry(IAComplex& 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( 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& 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& 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> subfaces; + small_vector_mp> 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> subcells; + small_vector_mp> 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 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; } \ No newline at end of file diff --git a/implicit_functions/include/implicit_function.hpp b/implicit_functions/include/implicit_function.hpp index 9e8392a..76fa288 100644 --- a/implicit_functions/include/implicit_function.hpp +++ b/implicit_functions/include/implicit_function.hpp @@ -2,15 +2,11 @@ #include "utils/eigen_alias.hpp" +template class ImplicitFunction { public: - template - double evaluate_scalar(const Eigen::Ref>& pos) const; - - template - Eigen::Vector evaluate_gradient(const Eigen::Ref>& pos) const; - - template - Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref>& pos) const; + Scalar evaluate_scalar(const Eigen::Ref>& pos) const; + Eigen::Vector evaluate_gradient(const Eigen::Ref>& pos) const; + Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref>& pos) const; }; \ No newline at end of file diff --git a/implicit_functions/include/implicit_functions.hpp b/implicit_functions/include/implicit_functions.hpp new file mode 100644 index 0000000..3f7d6e9 --- /dev/null +++ b/implicit_functions/include/implicit_functions.hpp @@ -0,0 +1,12 @@ +#pragma once + +#include + +#include + +template +using implicit_function_t = std::variant, + PlaneDistanceFunction, + CylinderDistanceFunction, + SphereDistanceFunction, + ConeDistanceFunction>; \ No newline at end of file diff --git a/implicit_functions/include/load_functions.hpp b/implicit_functions/include/load_functions.hpp index ccc65f1..49ed8ce 100644 --- a/implicit_functions/include/load_functions.hpp +++ b/implicit_functions/include/load_functions.hpp @@ -1,14 +1,14 @@ #pragma once -#include +#include #include -#include "implicit_function.hpp" +#include "implicit_functions.hpp" struct load_functions_result_t { - small_vector_mp> functions{}; - bool success{}; + small_vector_mp> functions{}; + bool success{}; }; -load_functions_result_t load_functions(const std::string_view& filename); \ No newline at end of file +load_functions_result_t load_functions(const std::string& filename); \ No newline at end of file diff --git a/implicit_functions/include/primitive_functions.hpp b/implicit_functions/include/primitive_functions.hpp index 6c84854..1764120 100644 --- a/implicit_functions/include/primitive_functions.hpp +++ b/implicit_functions/include/primitive_functions.hpp @@ -4,94 +4,89 @@ static constexpr double kEpsilon = 1e-6; -class ConstantFunction : public ImplicitFunction +template +class ConstantFunction : public ImplicitFunction { public: - explicit ConstantFunction(double value) : value_(value) {} + explicit ConstantFunction(Scalar value) : value_(value) {} - template - double evaluate_scalar(const Eigen::Ref> &pos) const - { - return value_; - } + Scalar evaluate_scalar(const Eigen::Ref> &pos) const { return value_; } - template - Eigen::Vector evaluate_gradient(const Eigen::Ref> &pos) const + Eigen::Vector evaluate_gradient(const Eigen::Ref> &pos) const { - return Eigen::Vector::Zero(); + return Eigen::Vector::Zero(); } - template - Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref> &pos) const + Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref> &pos) const { - auto res = Eigen::Vector::Zero(); + auto res = Eigen::Vector::Zero(); res[0] = value_; return res; } private: - double value_{}; + Scalar value_{}; }; -template -class PlaneDistanceFunction : public ImplicitFunction +template +class PlaneDistanceFunction : public ImplicitFunction { public: - PlaneDistanceFunction(const Eigen::Ref> &point, - const Eigen::Ref> &normal) + PlaneDistanceFunction(const Eigen::Ref> &point, + const Eigen::Ref> &normal) : point_(point), normal_(normal) { normal_.normalize(); } - double evaluate_scalar(const Eigen::Ref> &pos) const { return normal_.dot(pos - point_); } + Scalar evaluate_scalar(const Eigen::Ref> &pos) const { return normal_.dot(pos - point_); } - Eigen::Vector evaluate_gradient(const Eigen::Ref> &pos) const + Eigen::Vector evaluate_gradient(const Eigen::Ref> &pos) const { return normal_; } - Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref> &pos) const + Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref> &pos) const { - Eigen::Vector res{}; + Eigen::Vector res{}; res.template head().array() = normal_.array(); res[Dim] = evaluate_scalar(pos); return res; } private: - Eigen::Vector point_{}; - Eigen::Vector normal_{}; + Eigen::Vector point_{}; + Eigen::Vector normal_{}; }; -template -class SphereDistanceFunction : public ImplicitFunction +template +class SphereDistanceFunction : public ImplicitFunction { public: - SphereDistanceFunction(const Eigen::Ref> ¢er, double radius) + SphereDistanceFunction(const Eigen::Ref> ¢er, Scalar radius) : center_(center), radius_(radius) { } - double evaluate_scalar(const Eigen::Ref> &pos) const + Scalar evaluate_scalar(const Eigen::Ref> &pos) const { return (pos - center_).norm() - radius_; } - Eigen::Vector evaluate_gradient(const Eigen::Ref> &pos) const + Eigen::Vector evaluate_gradient(const Eigen::Ref> &pos) const { auto vec = pos - center_; - double norm = vec.norm(); - if (norm <= kEpsilon) return Eigen::Vector::Zero(); + Scalar norm = vec.norm(); + if (norm <= kEpsilon) return Eigen::Vector::Zero(); return vec / norm; } - Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref> &pos) const + Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref> &pos) const { - Eigen::Vector res{}; + Eigen::Vector 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 center_{}; - double radius_{}; + Eigen::Vector center_{}; + Scalar radius_{}; }; -template +template class CylinderDistanceFunction; -template <> -class CylinderDistanceFunction<3> : public ImplicitFunction +template +class CylinderDistanceFunction : public ImplicitFunction { public: - CylinderDistanceFunction(const Eigen::Ref &axis_point, - const Eigen::Ref &axis_dir, - double radius) + CylinderDistanceFunction(const Eigen::Ref> &axis_point, + const Eigen::Ref> &axis_dir, + Scalar radius) : axis_point_(axis_point), axis_dir_(axis_dir), radius_(radius) { axis_dir_.normalize(); } - double evaluate_scalar(const Eigen::Ref &pos) const + Scalar evaluate_scalar(const Eigen::Ref> &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 &pos) const + Eigen::Vector evaluate_gradient(const Eigen::Ref> &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 vec_perp = vec - vec_para; + Scalar norm = vec_perp.norm(); + if (norm <= kEpsilon) return Eigen::Vector::Zero(); return vec_perp / norm; } - Eigen::Vector4d evaluate_scalar_gradient(const Eigen::Ref &pos) const + Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref> &pos) const { - Eigen::Vector4d res{}; + Eigen::Vector 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 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 axis_point_{}; + Eigen::Vector axis_dir_{}; + Scalar radius_{}; }; -template +template class ConeDistanceFunction; -template <> -class ConeDistanceFunction<3> : public ImplicitFunction +template +class ConeDistanceFunction : public ImplicitFunction { public: - ConeDistanceFunction(const Eigen::Ref &apex_point, - const Eigen::Ref &axis_dir, - double apex_angle) + ConeDistanceFunction(const Eigen::Ref> &apex_point, + const Eigen::Ref> &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 &pos) const + Scalar evaluate_scalar(const Eigen::Ref> &pos) const { return apex_angle_cosine_ * (pos - apex_point_).norm() - (pos - apex_point_).dot(axis_dir_); } - Eigen::Vector3d evaluate_gradient(const Eigen::Ref &pos) const + Eigen::Vector evaluate_gradient(const Eigen::Ref> &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::Zero(); return vec_apex * apex_angle_cosine_ / vec_norm - axis_dir_; } - Eigen::Vector4d evaluate_scalar_gradient(const Eigen::Ref &pos) const + Eigen::Vector evaluate_scalar_gradient(const Eigen::Ref> &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::Zero(); - Eigen::Vector4d res{}; + Eigen::Vector 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 apex_point_{}; + Eigen::Vector axis_dir_{}; + Scalar apex_angle_cosine_{}; }; -template +template class TorusDistanceFunction; \ No newline at end of file diff --git a/implicit_functions/interface/implicit_functions.hpp b/implicit_functions/interface/implicit_functions.hpp deleted file mode 100644 index 764bf6d..0000000 --- a/implicit_functions/interface/implicit_functions.hpp +++ /dev/null @@ -1,4 +0,0 @@ -#pragma once - -#include -#include \ No newline at end of file diff --git a/implicit_functions/src/main.cpp b/implicit_functions/src/main.cpp index 2419863..3dc8e7b 100644 --- a/implicit_functions/src/main.cpp +++ b/implicit_functions/src/main.cpp @@ -3,22 +3,21 @@ #include -#include -#include "primitive_functions.hpp" +#include -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(); 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(); normal[i] = data[j]["normal"][i].get(); } - result.functions.emplace_back(std::make_unique>(point, normal)); + result.functions.emplace_back(PlaneDistanceFunction{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(); direction[i] = data[j]["direction"][i].get(); } - result.functions.emplace_back(std::make_unique>(point, direction, 0)); + result.functions.emplace_back(CylinderDistanceFunction{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(); axis_point_[i] = data[j]["axis_point2"][i].get(); @@ -55,28 +54,29 @@ load_functions_result_t load_functions(const std::string_view& filename) } } double radius = data[j]["radius"].get(); - result.functions.emplace_back(std::make_unique>(axis_point, axis_direction, radius)); + result.functions.emplace_back(CylinderDistanceFunction{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 radius = data[j]["radius"].get(); - result.functions.emplace_back(std::make_unique>(center, radius)); + result.functions.emplace_back(SphereDistanceFunction{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(); axis_direction[i] = data[j]["axis_direction"][i].get(); } double apex_angle = data[j]["apex_angle"].get(); - result.functions.emplace_back(std::make_unique>(apex_point, axis_direction, apex_angle)); + result.functions.emplace_back(ConeDistanceFunction{apex_point, axis_direction, apex_angle}); } else if (type == "constant") { double value = data[j]["value"].get(); - result.functions.emplace_back(std::make_unique(value)); + result.functions.emplace_back(ConstantFunction{value}); } else { std::cout << "Unsupported function type: " << type << std::endl; return {}; } } + result.success = true; return result; } \ No newline at end of file diff --git a/implicit_functions/xmake.lua b/implicit_functions/xmake.lua index 529d4ea..abfeaae 100644 --- a/implicit_functions/xmake.lua +++ b/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") diff --git a/network_process/include/background_mesh.hpp b/network_process/include/background_mesh.hpp index b88728e..ea3326e 100644 --- a/network_process/include/background_mesh.hpp +++ b/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); \ No newline at end of file +tetrahedron_mesh_t generate_tetrahedron_background_mesh(uint32_t resolution, + const Eigen::Ref& aabb_min, + const Eigen::Ref& aabb_max); \ No newline at end of file diff --git a/network_process/src/background_mesh.cpp b/network_process/src/background_mesh.cpp index 21b26cf..4fedc6c 100644 --- a/network_process/src/background_mesh.cpp +++ b/network_process/src/background_mesh.cpp @@ -4,30 +4,30 @@ #include 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& aabb_min, + const Eigen::Ref& 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(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(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}; diff --git a/network_process/src/implicit_surface_network_processor.cpp b/network_process/src/implicit_surface_network_processor.cpp index a4baaef..af1878d 100644 --- a/network_process/src/implicit_surface_network_processor.cpp +++ b/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(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 iso_verts{}; diff --git a/network_process/src/patch_connectivity.cpp b/network_process/src/patch_connectivity.cpp index 4634ba8..9ab0d49 100644 --- a/network_process/src/patch_connectivity.cpp +++ b/network_process/src/patch_connectivity.cpp @@ -127,7 +127,7 @@ void compute_chains(const stl_vector_mp& patch_edges, void compute_shells_and_components(uint32_t num_patch, const stl_vector_mp>& half_patch_pair_list, - stl_vector_mp>& shells, + stl_vector_mp>& shells, stl_vector_mp& shell_of_half_patch, stl_vector_mp>& components, stl_vector_mp& component_of_patch) diff --git a/network_process/xmake.lua b/network_process/xmake.lua index 6c0e055..074c315 100644 --- a/network_process/xmake.lua +++ b/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}) diff --git a/shared_module/container/detail/small_vector_base.hpp b/shared_module/container/detail/small_vector_base.hpp index f74d743..c90a172 100644 --- a/shared_module/container/detail/small_vector_base.hpp +++ b/shared_module/container/detail/small_vector_base.hpp @@ -1,5 +1,6 @@ #pragma once +#include #include #include diff --git a/shared_module/container/static_vector.hpp b/shared_module/container/static_vector.hpp index 6748244..acd376d 100644 --- a/shared_module/container/static_vector.hpp +++ b/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);