mesh_decimation

cmake/recipes/tests/wmtk_data.cmake

@@ -16,7 +16,7 @@ ExternalProject_Add(
     SOURCE_DIR ${WMTK_DATA_ROOT}
 
     GIT_REPOSITORY https://github.com/wildmeshing/data.git
-    GIT_TAG dfd3dc188aeb7e1c313e472db2ec56c5d14ba4b8
+    GIT_TAG 964b4918ceafd14dd2b0dfa18358b6b985fedfc5
 
     CONFIGURE_COMMAND ""
     BUILD_COMMAND ""

components/CMakeLists.txt

@@ -57,6 +57,7 @@ add_component(${WMTK_COMPONENT_PREFIX} "procedural")
 add_component(${WMTK_COMPONENT_PREFIX} "fusion")
 add_component(${WMTK_COMPONENT_PREFIX} "triangle_insertion")
 add_component(${WMTK_COMPONENT_PREFIX} "to_points")
+add_component(${WMTK_COMPONENT_PREFIX} "mesh_decimation")
 
 string(LENGTH ${json_components} json_components_length)
 math(EXPR json_components_length "${json_components_length}-2")

components/tests/CMakeLists.txt

@@ -17,7 +17,6 @@ wmtk::components)
 target_include_directories(${WMTK_COMPONENT_TEST_TARGET} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/..)
 target_compile_definitions(${WMTK_COMPONENT_TEST_TARGET} PRIVATE WMTK_TEST_DIR=\"${CMAKE_CURRENT_SOURCE_DIR}\")
 
-
 FetchContent_GetProperties(catch2)
 list(APPEND CMAKE_MODULE_PATH ${catch2_SOURCE_DIR}/extras)
 include(Catch)

components/tests/integration_test.cpp

@@ -58,16 +58,17 @@ int authenticate_json(const std::string& json_file, const bool compute_validatio
     in_args["root_path"] = json_file;
 
 
-    if (compute_validation && !contains_results(in_args)) {
-        spdlog::error("JSON file missing vertices edges faces or tetrahedra or meshes key. Add a * "
-                      "to the beginning of filename to allow appends.");
-        return 2;
-    }
-
-    //in_args["settings"] = R"({
-    //    "log_level": 5,
-    //    "opt_log_level": 5
-    //    })"_json;
+    // if (compute_validation && !contains_results(in_args)) {
+    //     spdlog::error("JSON file missing vertices edges faces or tetrahedra or meshes key. Add a
+    //     * "
+    //                   "to the beginning of filename to allow appends.");
+    //     return 2;
+    // }
+
+    // in_args["settings"] = R"({
+    //     "log_level": 5,
+    //     "opt_log_level": 5
+    //     })"_json;
 
     utils::set_random_seed(0);
     auto cache = wmtk::components::run_components(in_args, true);
@@ -88,9 +89,10 @@ int authenticate_json(const std::string& json_file, const bool compute_validatio
         auto tetrahedra = in_args["tests"]["tetrahedra"];
         if (meshes.size() != vertices.size() || meshes.size() != edges.size() ||
             meshes.size() != faces.size() || meshes.size() != tetrahedra.size()) {
-            spdlog::error("JSON size missmatch between meshes and vertices edges faces or "
-                          "tetrahedra or meshes key. Add a * "
-                          "to the beginning of filename to allow appends.");
+            spdlog::error(
+                "JSON size missmatch between meshes and vertices edges faces or "
+                "tetrahedra or meshes key. Set true for the parameter \"compute_validation\""
+                "to allow appends.");
             return 2;
         }
 
@@ -174,29 +176,29 @@ std::string tagsrun = "[.][integration]";
 #endif
 } // namespace
 
-#define WMTK_INTEGRATION(NAME, DO_VALIDATION)\
-TEST_CASE(std::string("integration_") + NAME, tagsrun) \
-{ \
-    std::string path = std::string("unit_test/") + NAME + ".json"; \
-    bool compute_validation = DO_VALIDATION; \
-    spdlog::info("Processing {}", NAME); \
-    auto flag = authenticate_json(WMTK_DATA_DIR "/" + path, compute_validation); \
-    REQUIRE(flag == 0); \
-}
-
+#define WMTK_INTEGRATION(NAME, DO_VALIDATION)                                        \
+    TEST_CASE(std::string("integration_") + NAME, tagsrun)                           \
+    {                                                                                \
+        std::string path = std::string("unit_test/") + NAME + ".json";               \
+        bool compute_validation = DO_VALIDATION;                                     \
+        spdlog::info("Processing {}", NAME);                                         \
+        auto flag = authenticate_json(WMTK_DATA_DIR "/" + path, compute_validation); \
+        REQUIRE(flag == 0);                                                          \
+    }
 
 
-WMTK_INTEGRATION("input",false);
-WMTK_INTEGRATION("to_points",false);
-WMTK_INTEGRATION("delaunay",false);
-WMTK_INTEGRATION("insertion",false);
-WMTK_INTEGRATION("insertion_open",false);
-WMTK_INTEGRATION("multimesh",false);
-WMTK_INTEGRATION("multimesh_boundary_2d",false);
-WMTK_INTEGRATION("multimesh_boundary_3d",false);
-WMTK_INTEGRATION("isotropic_remeshing",false);
-WMTK_INTEGRATION("isotropic_remeshing_mm",false);
-WMTK_INTEGRATION("disk_fan_mm",false);
-//WMTK_INTEGRATION("grid",false);
-WMTK_INTEGRATION("wildmeshing_2d",false);
-WMTK_INTEGRATION("wildmeshing_3d",false);
+WMTK_INTEGRATION("input", false);
+WMTK_INTEGRATION("to_points", false);
+WMTK_INTEGRATION("delaunay", false);
+WMTK_INTEGRATION("insertion", false);
+WMTK_INTEGRATION("insertion_open", false);
+WMTK_INTEGRATION("multimesh", false);
+WMTK_INTEGRATION("multimesh_boundary_2d", false);
+WMTK_INTEGRATION("multimesh_boundary_3d", false);
+WMTK_INTEGRATION("isotropic_remeshing", false);
+WMTK_INTEGRATION("isotropic_remeshing_mm", false);
+WMTK_INTEGRATION("disk_fan_mm", false);
+// WMTK_INTEGRATION("grid",false);
+WMTK_INTEGRATION("wildmeshing_2d", false);
+WMTK_INTEGRATION("wildmeshing_3d", false);
+WMTK_INTEGRATION("mesh_decimation", false);

components/tests/test_component_mesh_decimation.cpp

@@ -0,0 +1,64 @@
+#include <catch2/catch_test_macros.hpp>
+#include <nlohmann/json.hpp>
+#include <tools/DEBUG_TetMesh.hpp>
+#include <tools/DEBUG_TriMesh.hpp>
+#include <tools/TetMesh_examples.hpp>
+#include <tools/TriMesh_examples.hpp>
+#include <wmtk/Mesh.hpp>
+#include <wmtk/TriMesh.hpp>
+#include <wmtk/components/mesh_decimation/internal/MeshDecimation.hpp>
+#include <wmtk/components/mesh_decimation/internal/MeshDecimationOptions.hpp>
+#include <wmtk/components/mesh_decimation/mesh_decimation.hpp>
+#include <wmtk/function/simplex/AMIPS.hpp>
+#include <wmtk/invariants/TodoInvariant.hpp>
+#include <wmtk/io/MeshReader.hpp>
+#include <wmtk/io/ParaviewWriter.hpp>
+#include <wmtk/operations/EdgeCollapse.hpp>
+#include <wmtk/operations/attribute_new/CollapseNewAttributeStrategy.hpp>
+#include <wmtk/simplex/link.hpp>
+#include <wmtk/utils/mesh_utils.hpp>
+
+using json = nlohmann::json;
+using namespace wmtk;
+
+const std::filesystem::path data_dir = WMTK_DATA_DIR;
+
+TEST_CASE("component_mesh_decimation_options", "[components][mesh_decimation]")
+{
+    using namespace components::internal;
+
+    json o = {
+        {"input", "input_mesh"},
+        {"output", "output_mesh"},
+        {"constrait_value", 1},
+        {"target_len", 1.0},
+        {"cell_constrait_tag_name", "tag"},
+        {"pass_through", {"vertices"}}};
+
+    CHECK_NOTHROW(o.get<MeshDecimationOptions>());
+}
+
+TEST_CASE("preprocess", "[.]")
+{
+    using namespace wmtk::components;
+
+    auto mesh_in = wmtk::read_mesh(data_dir / "3d_images/sphere_regularized.hdf5");
+    Mesh& mesh = static_cast<Mesh&>(*mesh_in);
+
+    std::vector<wmtk::attribute::MeshAttributeHandle> pass_though;
+    wmtk::attribute::MeshAttributeHandle vertex =
+        mesh.get_attribute_handle<int64_t>("vertex_tag", PrimitiveType::Vertex);
+    wmtk::attribute::MeshAttributeHandle edge =
+        mesh.get_attribute_handle<int64_t>("edge_tag", PrimitiveType::Edge);
+    wmtk::attribute::MeshAttributeHandle face =
+        mesh.get_attribute_handle<int64_t>("face_tag", PrimitiveType::Triangle);
+    pass_though.push_back(vertex);
+    pass_though.push_back(edge);
+    pass_though.push_back(face);
+    internal::MeshDecimation MD(mesh, "tag", 1, 5, pass_though);
+    MD.process();
+
+    wmtk::io::ParaviewWriter
+        writer(data_dir / "3d_images/out.hdf", "vertices", mesh, false, true, true, false);
+    mesh.serialize(writer);
+}

components/wmtk_components/mesh_decimation/wmtk/components/mesh_decimation/CMakeLists.txt

@@ -0,0 +1,13 @@
+
+set(SRC_FILES
+    internal/MeshDecimationOptions.hpp
+    internal/MeshDecimationOptions.cpp
+    internal/MeshDecimation.hpp
+    internal/MeshDecimation.cpp
+	mesh_decimation.hpp
+    mesh_decimation.cpp
+    )
+
+
+#CURRENT_COMPONENT_LIB_NAME is set from the main cmake
+target_sources(${CURRENT_COMPONENT_LIB_NAME} PRIVATE ${SRC_FILES})

components/wmtk_components/mesh_decimation/wmtk/components/mesh_decimation/internal/MeshDecimation.cpp

@@ -0,0 +1,172 @@
+#include "MeshDecimation.hpp"
+#include "wmtk/Scheduler.hpp"
+#include "wmtk/function/simplex/AMIPS.hpp"
+#include "wmtk/invariants/MultiMeshLinkConditionInvariant.hpp"
+#include "wmtk/invariants/SimplexInversionInvariant.hpp"
+#include "wmtk/invariants/SmallerFunctionInvariant.hpp"
+#include "wmtk/invariants/TodoInvariant.hpp"
+#include "wmtk/operations/EdgeCollapse.hpp"
+#include "wmtk/operations/attribute_new/CollapseNewAttributeStrategy.hpp"
+#include "wmtk/operations/attribute_update/AttributeTransferStrategy.hpp"
+
+namespace wmtk::components::internal {
+
+MeshDecimation::MeshDecimation(
+    Mesh& mesh,
+    std::string constriant_name,
+    int64_t constrait_value,
+    double target_len,
+    const std::vector<attribute::MeshAttributeHandle>& pass_through_attributes)
+    : m_mesh(mesh)
+    , m_constrait_value(constrait_value)
+    , m_constriant_name(constriant_name)
+    , m_target_len(target_len)
+    , m_pass_through_attributes(pass_through_attributes)
+{}
+
+void MeshDecimation::process()
+{
+    using namespace wmtk::attribute;
+    using namespace wmtk::invariants;
+    PrimitiveType PV = PrimitiveType::Vertex;
+    PrimitiveType PE = PrimitiveType::Edge;
+    PrimitiveType PF = PrimitiveType::Triangle;
+    PrimitiveType PT = PrimitiveType::Tetrahedron;
+
+    volatile PrimitiveType x = m_mesh.top_simplex_type();
+
+    MeshAttributeHandle cell_tag_handle =
+        m_mesh.get_attribute_handle<int64_t>(m_constriant_name, m_mesh.top_simplex_type());
+    MeshAttributeHandle position = m_mesh.get_attribute_handle<double>("vertices", PV);
+    MeshAttributeHandle edge_handle = m_mesh.register_attribute<int64_t>("todo_edge_", PE, 1);
+    MeshAttributeHandle vertex_handle = m_mesh.register_attribute<int64_t>("todo_vertex_", PV, 1);
+    MeshAttributeHandle edge_len_handle = m_mesh.register_attribute<double>("len_edge_", PE, 1);
+
+    Accessor<int64_t> acc_cell = m_mesh.create_accessor<int64_t>(cell_tag_handle);
+    Accessor<double> acc_pos = m_mesh.create_accessor<double>(position);
+    Accessor<int64_t> acc_edge = m_mesh.create_accessor<int64_t>(edge_handle);
+    Accessor<int64_t> acc_vertex = m_mesh.create_accessor<int64_t>(vertex_handle);
+    Accessor<double> acc_len = m_mesh.create_accessor<double>(edge_len_handle);
+
+    switch (m_mesh.top_cell_dimension()) {
+    case 2:
+        for (const Tuple& edge : m_mesh.get_all(PE)) {
+            if (m_mesh.is_boundary(PE, edge)) {
+                acc_vertex.scalar_attribute(edge) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuple(edge, PV)) = 1;
+
+                acc_edge.scalar_attribute(edge) = 1;
+            } else if (
+                acc_cell.scalar_attribute(edge) !=
+                acc_cell.scalar_attribute(m_mesh.switch_tuple(edge, PF))) {
+                acc_vertex.scalar_attribute(edge) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuple(edge, PV)) = 1;
+
+                acc_edge.scalar_attribute(edge) = 1;
+            }
+        }
+        break;
+    case 3: {
+        for (const Tuple& face : m_mesh.get_all(PF)) {
+            if (m_mesh.is_boundary(PF, face)) {
+                acc_vertex.scalar_attribute(face) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuple(face, PV)) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuples(face, {PE, PV})) = 1;
+
+                acc_edge.scalar_attribute(face) = 1;
+                acc_edge.scalar_attribute(m_mesh.switch_tuple(face, PE)) = 1;
+                acc_edge.scalar_attribute(m_mesh.switch_tuples(face, {PV, PE})) = 1;
+            } else if (
+                acc_cell.scalar_attribute(face) !=
+                acc_cell.scalar_attribute(m_mesh.switch_tuple(face, PT))) {
+                acc_vertex.scalar_attribute(face) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuple(face, PV)) = 1;
+                acc_vertex.scalar_attribute(m_mesh.switch_tuples(face, {PE, PV})) = 1;
+
+                acc_edge.scalar_attribute(face) = 1;
+                acc_edge.scalar_attribute(m_mesh.switch_tuple(face, PE)) = 1;
+                acc_edge.scalar_attribute(m_mesh.switch_tuples(face, {PV, PE})) = 1;
+            }
+        }
+        break;
+    }
+    default:
+        std::runtime_error("MeshDecimation.cpp: mesh_decimation component only supports tetmesh "
+                           "and trimesh for now!");
+    }
+
+    for (const Tuple& edge : m_mesh.get_all(PE)) {
+        if (acc_vertex.scalar_attribute(edge) == 1 ||
+            acc_vertex.scalar_attribute(m_mesh.switch_tuple(edge, PV)) == 1) {
+            acc_edge.scalar_attribute(edge) = 1;
+        }
+        acc_len.scalar_attribute(edge) = (acc_pos.vector_attribute(edge) -
+                                          acc_pos.vector_attribute(m_mesh.switch_tuple(edge, PV)))
+                                             .norm();
+    }
+
+    auto op_scaffold = std::make_shared<operations::EdgeCollapse>(m_mesh);
+
+    op_scaffold->add_invariant(
+        std::make_shared<TodoInvariant>(m_mesh, edge_handle.as<int64_t>(), 0));
+    op_scaffold->add_invariant(std::make_shared<TodoSmallerInvariant>(
+        m_mesh,
+        edge_len_handle.as<double>(),
+        4.0 / 5.0 * m_target_len));
+
+    auto m_amips = std::make_shared<function::AMIPS>(m_mesh, position);
+    auto m_link_conditions = std::make_shared<InvariantCollection>(m_mesh);
+    m_link_conditions->add(std::make_shared<MultiMeshLinkConditionInvariant>(m_mesh));
+    auto m_function_invariant =
+        std::make_shared<SmallerFunctionInvariant>(m_mesh.top_simplex_type(), m_amips, 30);
+    auto m_inversion_invariant =
+        std::make_shared<SimplexInversionInvariant>(m_mesh, position.as<double>());
+
+    // Edge length update
+    auto compute_edge_length = [](const Eigen::MatrixXd& P) -> Eigen::VectorXd {
+        assert(P.cols() == 2);
+        assert(P.rows() == 2 || P.rows() == 3);
+        return Eigen::VectorXd::Constant(1, (P.col(0) - P.col(1)).norm());
+    };
+    auto m_edge_length_update =
+        std::make_shared<wmtk::operations::SingleAttributeTransferStrategy<double, double>>(
+            edge_len_handle,
+            position,
+            compute_edge_length);
+    auto m_prio_short_edges_first = [&](const simplex::Simplex& s) {
+        assert(s.primitive_type() == PrimitiveType::Edge);
+        auto acc = m_mesh.create_accessor<double>(edge_len_handle);
+        return std::vector<double>({acc.scalar_attribute(s.tuple())});
+    };
+
+
+    op_scaffold->add_invariant(m_link_conditions);
+    op_scaffold->add_invariant(m_function_invariant);
+    op_scaffold->add_invariant(m_inversion_invariant);
+
+    op_scaffold->add_transfer_strategy(m_edge_length_update);
+
+    op_scaffold->set_priority(m_prio_short_edges_first);
+
+    op_scaffold->set_new_attribute_strategy(
+        position,
+        wmtk::operations::CollapseBasicStrategy::Mean);
+    op_scaffold->set_new_attribute_strategy(vertex_handle);
+
+    op_scaffold->set_new_attribute_strategy(edge_handle);
+    op_scaffold->set_new_attribute_strategy(edge_len_handle);
+    op_scaffold->set_new_attribute_strategy(cell_tag_handle);
+
+    // pass_through
+    for (const auto& attr : m_pass_through_attributes) {
+        op_scaffold->set_new_attribute_strategy(attr);
+    }
+
+    while (true) {
+        Scheduler scheduler;
+        SchedulerStats pass_stats = scheduler.run_operation_on_all(*op_scaffold);
+        if (pass_stats.number_of_successful_operations() == 0) break;
+    }
+}
+
+} // namespace wmtk::components::internal