Merge branch 'main' into spenbert/offsets

* main:
  Simplify before embedding if topology preservation is not required.
  Make shortest edge collapse work for non-manifold meshes.
  Fix the issue with the empty envelope after simplification.

Author: Spenbert02

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/EdgeCollapsing.cpp

@@ -389,7 +389,11 @@ void ImageSimulationMesh::simplify()
         m_envelope = nullptr;
         m_V_envelope.clear();
         m_F_envelope.clear();
-        init_envelope(V, F);
+        if (V.size() > 0 && F.size() > 0) {
+            init_envelope(V, F);
+        } else {
+            logger().warn("No surface faces left after simplification, skip re-building envelope");
+        }
     }
     logger().info("===== Simplification done =====");
 }

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/ImageSimulationMesh.cpp

@@ -586,8 +586,11 @@ void ImageSimulationMesh::init_envelope(const MatrixXd& V, const MatrixXi& F)
     if (m_envelope) {
         log_and_throw_error("Envelope was already initialized once.");
     }
+    if (V.size() == 0 || F.size() == 0) {
+        log_and_throw_error("Envelope vertices and faces cannot be empty.");
+    }
+
     assert(m_V_envelope.empty() && m_F_envelope.empty());
-    assert(V.size() != 0 && F.size() != 0);
     assert(V.cols() == 3); // vertices must be in 3D
     assert(F.cols() == 3); // envelope must be triangles
 

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;