Simplification of non-manifold input

cmake/wmtk_data.cmake

@@ -16,7 +16,7 @@ ExternalProject_Add(
     SOURCE_DIR ${WMTK_DATA_ROOT}
 
     GIT_REPOSITORY https://github.com/wildmeshing/data2.git
-    GIT_TAG 5cf9136e9004c15605227001a67df971af025fd7
+    GIT_TAG 36bb3968aa1e685d3da2e38b87bc6fa5ab328a13
 
     CONFIGURE_COMMAND ""
     BUILD_COMMAND ""

components/image_simulation/wmtk/components/image_simulation/EmbedSurface.cpp

@@ -19,6 +19,7 @@
 #include <bitset>
 #include <filesystem>
 #include <paraviewo/VTUWriter.hpp>
+#include <wmtk/io/read_triangle_mesh.hpp>
 #include <wmtk/utils/InsertTriangleUtils.hpp>
 #include <wmtk/utils/Logger.hpp>
 #include <wmtk/utils/ManifoldUtils.hpp>
@@ -856,7 +857,9 @@ EmbedSurface::EmbedSurface(const std::vector<std::string>& img_filenames, const
 
 EmbedSurface::EmbedSurface(
     const std::vector<std::string>& img_filenames,
-    const std::vector<Matrix4d>& img_transform)
+    const std::vector<Matrix4d>& img_transform,
+    const double tol_rel,
+    const double tol_abs)
     : m_img_filenames(img_filenames)
 {
     assert(img_filenames.size() == img_transform.size());
@@ -869,7 +872,7 @@ EmbedSurface::EmbedSurface(
             log_and_throw_error("Input file {} does not exist", m_img_filenames[i]);
         }
         MatrixXi F_single;
-        igl::read_triangle_mesh(m_img_filenames[i], Vs[i], F_single);
+        io::read_triangle_mesh(m_img_filenames[i], Vs[i], F_single, tol_rel, tol_abs);
 
         assert(Vs[i].cols() == 3);
         assert(F_single.cols() == 3);
@@ -920,13 +923,13 @@ EmbedSurface::EmbedSurface(
     F_surf_from_vector(tris);
 }
 
-void EmbedSurface::simplify_surface(const double eps)
+void EmbedSurface::simplify_surface(const double eps, const int num_threads)
 {
     // convert to STL vectors
     std::vector<Eigen::Vector3d> verts = V_surf_to_vector();
     std::vector<std::array<size_t, 3>> tris = F_surf_to_vector();
 
-    shortest_edge_collapse::ShortestEdgeCollapse surf_mesh(verts, 0, false);
+    shortest_edge_collapse::ShortestEdgeCollapse surf_mesh(verts, num_threads, false);
 
     // must be a small envelope to ensure correct tet tags later on
     surf_mesh.create_mesh(verts.size(), tris, modified_nonmanifold_v, eps);
@@ -1375,6 +1378,14 @@ std::pair<Vector3d, Vector3d> EmbedSurface::bbox_minmax() const
     return std::make_pair(bbox_min, bbox_max);
 }
 
+std::pair<Vector3d, Vector3d> EmbedSurface::bbox_surf_minmax() const
+{
+    const Vector3d bbox_max = m_V_surface.colwise().maxCoeff();
+    const Vector3d bbox_min = m_V_surface.colwise().minCoeff();
+
+    return std::make_pair(bbox_min, bbox_max);
+}
+
 std::vector<Eigen::Vector3d> EmbedSurface::V_surf_to_vector() const
 {
     std::vector<Eigen::Vector3d> verts;

components/image_simulation/wmtk/components/image_simulation/EmbedSurface.hpp

@@ -101,14 +101,16 @@ class EmbedSurface
      */
     EmbedSurface(
         const std::vector<std::string>& img_filenames,
-        const std::vector<Matrix4d>& img_transform = {});
+        const std::vector<Matrix4d>& img_transform = {},
+        const double tol_rel = -1,
+        const double tol_abs = -1);
 
     /**
      * @brief Simplify the input surface while staying within the eps envelope.
      *
      * @param eps The absolute envelope thickness.
      */
-    void simplify_surface(const double eps);
+    void simplify_surface(const double eps, const int num_threads = 0);
 
     /**
      * @brief Merge vertices that are closer than eps.
@@ -151,6 +153,7 @@ class EmbedSurface
     void write_emb_vtu(const std::string& filename) const;
 
     std::pair<Vector3d, Vector3d> bbox_minmax() const;
+    std::pair<Vector3d, Vector3d> bbox_surf_minmax() const;
 
     std::vector<Eigen::Vector3d> V_surf_to_vector() const;
     std::vector<std::array<size_t, 3>> F_surf_to_vector() const;

components/image_simulation/wmtk/components/image_simulation/image_simulation.cpp

@@ -99,8 +99,9 @@ void run_3D(const nlohmann::json& json_params, const InputData& input_data)
     write_unique_vtu();
 
 
-    if (!skip_simplify) {
+    if (params.preserve_topology && !skip_simplify) {
         // collapse for getting the right edge length
+        logger().info("Simplify after insertion");
         mesh.simplify();
     }
 

components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.hpp

@@ -147,7 +147,7 @@ nlohmann::json image_simulation_spec = R"(
   {
     "pointer": "/skip_simplify",
     "type": "bool",
-    "default": true,
+    "default": false,
     "doc": "If true, input simplification will be skipped."
   },
   {
@@ -189,7 +189,7 @@ nlohmann::json image_simulation_spec = R"(
   {
     "pointer": "/eps_simplify_rel",
     "type": "float",
-    "default": 2e-3,
+    "default": 2e-4,
     "doc": "Envelope thickness relative to the bounding box for the initial simplification."
   },
   {

components/image_simulation/wmtk/components/image_simulation/image_simulation_spec.json

@@ -142,7 +142,7 @@
   {
     "pointer": "/skip_simplify",
     "type": "bool",
-    "default": true,
+    "default": false,
     "doc": "If true, input simplification will be skipped."
   },
   {
@@ -184,7 +184,7 @@
   {
     "pointer": "/eps_simplify_rel",
     "type": "float",
-    "default": 2e-3,
+    "default": 2e-4,
     "doc": "Envelope thickness relative to the bounding box for the initial simplification."
   },
   {

components/image_simulation/wmtk/components/image_simulation/read_image_msh.cpp

@@ -475,17 +475,47 @@ InputData read_mesh(
     const bool use_sample_envelope = json_params["use_sample_envelope"];
     const int NUM_THREADS = json_params["num_threads"];
     const int max_its = json_params["max_iterations"];
-    const bool write_vtu = json_params["write_vtu"];
+    const bool debug_output = json_params["DEBUG_output"];
     const bool preserve_topology = json_params["preserve_topology"];
+    const bool skip_simplify = json_params["skip_simplify"];
+    const double epsr_simplify = json_params["eps_simplify_rel"];
+    double eps_simplify = json_params["eps_simplify"];
     const std::vector<std::string> input_names = json_params["input_names"];
 
+    double tol_rel = -1;
+    double tol_abs = -1;
+    if (!preserve_topology) {
+        tol_rel = 0.1 * epsr_simplify;
+        tol_abs = 0.1 * eps_simplify;
+    }
+
     // convert mesh into tet mesh
-    EmbedSurface image_mesh(input_paths, input_transforms);
+    EmbedSurface image_mesh(input_paths, input_transforms, tol_rel, tol_abs);
 
-    if (write_vtu) {
+    if (debug_output) {
         image_mesh.write_surf_off(output_filename + "_input.off");
     }
 
+    if (!preserve_topology && !skip_simplify) {
+        /**
+         * Simplify the input surface only if topology preservation is not required, since
+         * simplification can change the topology.
+         * If topology preservation is required, simplification is performed after insertion.
+         */
+        auto [bbox_min, bbox_max] = image_mesh.bbox_surf_minmax();
+        double diag = (bbox_max - bbox_min).norm();
+        if (epsr_simplify > 0) {
+            eps_simplify = diag * epsr_simplify;
+        }
+        logger().info("Simplify before insertion with eps = {:.4}", eps_simplify);
+
+        image_mesh.simplify_surface(eps_simplify, NUM_THREADS);
+
+        if (debug_output) {
+            image_mesh.write_surf_off(output_filename + "_input_simplified.off");
+        }
+    }
+
     input_data.V_envelope = image_mesh.V_surface();
     input_data.F_envelope = image_mesh.F_surface();
 
@@ -494,7 +524,7 @@ InputData read_mesh(
                                  : image_mesh.embed_surface(preserve_topology);
     image_mesh.consolidate();
 
-    if (write_vtu) {
+    if (debug_output) {
         image_mesh.write_emb_surf_off(output_filename + "_input_emb.off");
     }
 

components/shortest_edge_collapse/wmtk/components/shortest_edge_collapse/ShortestEdgeCollapse.cpp

@@ -6,6 +6,7 @@
 
 #include <Eigen/Core>
 #include <Eigen/Geometry>
+#include <paraviewo/VTUWriter.hpp>
 
 namespace wmtk::components::shortest_edge_collapse {
 
@@ -24,23 +25,16 @@ ShortestEdgeCollapse::ShortestEdgeCollapse(
         vertex_attrs[i] = {_m_vertex_positions[i], 0, false};
 }
 
-void ShortestEdgeCollapse::set_freeze(const TriMesh::Tuple& v)
+void ShortestEdgeCollapse::freeze_boundary()
 {
-    for (const Tuple& e : get_one_ring_edges_for_vertex(v)) {
+    for (const Tuple& e : get_edges()) {
         if (is_boundary_edge(e)) {
-            vertex_attrs[v.vid(*this)].freeze = true;
-            continue;
+            vertex_attrs[e.vid(*this)].freeze = true;
+            vertex_attrs[e.switch_vertex(*this).vid(*this)].freeze = true;
         }
     }
 }
 
-void ShortestEdgeCollapse::create_mesh_nofreeze(
-    size_t n_vertices,
-    const std::vector<std::array<size_t, 3>>& tris)
-{
-    wmtk::TriMesh::init(n_vertices, tris);
-}
-
 void ShortestEdgeCollapse::create_mesh(
     size_t n_vertices,
     const std::vector<std::array<size_t, 3>>& tris,
@@ -65,9 +59,7 @@ void ShortestEdgeCollapse::create_mesh(
     for (size_t v : frozen_verts) {
         vertex_attrs[v].freeze = true;
     }
-    for (const Tuple& v : get_vertices()) { // the better way is to iterate through edges.
-        set_freeze(v);
-    }
+    freeze_boundary();
 }
 
 void ShortestEdgeCollapse::partition_mesh()
@@ -113,6 +105,36 @@ bool ShortestEdgeCollapse::write_triangle_mesh(std::string path)
     return igl::write_triangle_mesh(path, V, F);
 }
 
+void ShortestEdgeCollapse::write_vtu(const std::string& path)
+{
+    const std::string out_path = path + ".vtu";
+    logger().info("Write {}", out_path);
+
+    MatrixXd V = MatrixXd::Zero(vert_capacity(), 3);
+    MatrixXi F = MatrixXi::Zero(tri_capacity(), 3);
+
+    VectorXd freeze(vertex_attrs.size());
+    freeze.setZero();
+
+    for (Tuple& t : get_vertices()) {
+        const size_t i = t.vid(*this);
+        V.row(i) = vertex_attrs[i].pos;
+        freeze(i) = vertex_attrs[i].freeze ? 1 : 0;
+    }
+
+    for (Tuple& t : get_faces()) {
+        const size_t i = t.fid(*this);
+        const auto vs = oriented_tri_vertices(t);
+        for (int j = 0; j < 3; j++) {
+            F(i, j) = (int)vs[j].vid(*this);
+        }
+    }
+
+    paraviewo::VTUWriter writer;
+    writer.add_field("freeze", freeze);
+    writer.write_mesh(out_path, V, F, paraviewo::CellType::Triangle);
+}
+
 bool ShortestEdgeCollapse::collapse_edge_before(const Tuple& t)
 {
     if (!TriMesh::collapse_edge_before(t)) return false;

components/shortest_edge_collapse/wmtk/components/shortest_edge_collapse/ShortestEdgeCollapse.h

@@ -45,14 +45,12 @@ class ShortestEdgeCollapse : public wmtk::TriMesh
         int num_threads = 1,
         bool use_exact_envelope = true);
 
-    void set_freeze(const TriMesh::Tuple& v);
-
-    void create_mesh_nofreeze(size_t n_vertices, const std::vector<std::array<size_t, 3>>& tris);
+    void freeze_boundary();
 
     void create_mesh(
         size_t n_vertices,
         const std::vector<std::array<size_t, 3>>& tris,
-        const std::vector<size_t>& frozen_verts = std::vector<size_t>(),
+        const std::vector<size_t>& frozen_verts = {},
         double eps = 0);
 
     ~ShortestEdgeCollapse() {}
@@ -64,6 +62,8 @@ class ShortestEdgeCollapse : public wmtk::TriMesh
         return vertex_attrs[loc.vid(*this)].partition_id;
     }
 
+    void write_vtu(const std::string& path);
+
 public:
     bool collapse_edge_before(const Tuple& t) override;
     bool collapse_edge_after(const Tuple& t) override;

components/shortest_edge_collapse/wmtk/components/shortest_edge_collapse/shortest_edge_collapse.cpp

@@ -4,12 +4,6 @@
 #include <CLI/CLI.hpp>
 
 #include <igl/Timer.h>
-#include <igl/is_edge_manifold.h>
-#include <igl/is_vertex_manifold.h>
-#include <igl/readOFF.h>
-#include <igl/read_triangle_mesh.h>
-#include <igl/remove_duplicate_vertices.h>
-#include <igl/writeDMAT.h>
 
 #include <stdlib.h>
 #include <chrono>
@@ -19,6 +13,7 @@
 #include <jse/jse.h>
 #include <nlohmann/json.hpp>
 #include <wmtk/Types.hpp>
+#include <wmtk/io/read_triangle_mesh.hpp>
 #include <wmtk/utils/Logger.hpp>
 #include <wmtk/utils/resolve_path.hpp>
 
@@ -45,6 +40,7 @@ void run_shortest_collapse(
     logger().info("runtime {:.4}s", timer.getElapsedTimeInSec());
     m.consolidate_mesh();
     m.write_triangle_mesh(output);
+    // m.write_vtu(std::filesystem::path(output).stem().string());
     logger().info("After collapse: #V = {}, #F = {}", m.vert_capacity(), m.tri_capacity());
 }
 
@@ -73,20 +69,7 @@ void shortest_edge_collapse(nlohmann::json json_params)
 
     MatrixXd V;
     MatrixXi F;
-    // read input and remove duplicates
-    {
-        if (!igl::read_triangle_mesh(path, V, F)) {
-            log_and_throw_error("Could not read mesh {}", path);
-        }
-        VectorXi SVI, SVJ;
-        MatrixXd temp_V = V;
-        igl::remove_duplicate_vertices(temp_V, 0, V, SVI, SVJ);
-        for (int i = 0; i < F.rows(); i++) {
-            for (int j = 0; j < 3; ++j) {
-                F(i, j) = SVJ[F(i, j)];
-            }
-        }
-    }
+    io::read_triangle_mesh(path, V, F, -1, -1);
     logger().info("Input: #V = {}, #F = {}", V.rows(), F.rows());
 
 
@@ -108,11 +91,11 @@ void shortest_edge_collapse(nlohmann::json json_params)
     const double envelope_size = env_rel * diag;
     VectorXi dummy;
     std::vector<size_t> modified_v;
-    if (!igl::is_edge_manifold(F) || !igl::is_vertex_manifold(F, dummy)) {
-        auto v1 = v;
-        auto tri1 = tri;
-        wmtk::separate_to_manifold(v1, tri1, v, tri, modified_v);
-    }
+    // if (!igl::is_edge_manifold(F) || !igl::is_vertex_manifold(F, dummy)) {
+    //     auto v1 = v;
+    //     auto tri1 = tri;
+    //     wmtk::separate_to_manifold(v1, tri1, v, tri, modified_v);
+    // }
 
     ShortestEdgeCollapse m(v, num_threads, !use_sample_envelope);
     m.create_mesh(v.size(), tri, modified_v, envelope_size);