Merge pull request #868 from wildmeshing/dzint/tet_remeshing

Dzint/tet remeshing

Author: daniel-zint

CMakeLists.txt

@@ -126,13 +126,13 @@ target_link_libraries(wildmeshing_toolkit PUBLIC
 
 if(MSVC)
     target_compile_options(wildmeshing_toolkit PUBLIC /MP)
-    target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELWITHDEBINFO>:/GL>)
-    target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELWITHDEBINFO>:/LTCG>)
-    target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELWITHDEBINFO>:/Ob2>) # inline whenever suitable
-    target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELWITHDEBINFO>:/Ot>) # favor faster code over binary size
-    target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELEASE>:/GL>)
-    target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELEASE>:/LTCG>)
-    target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELEASE>:/Ot>)
+    # target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELWITHDEBINFO>:/GL>)
+    # target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELWITHDEBINFO>:/LTCG>)
+    # target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELWITHDEBINFO>:/Ob2>) # inline whenever suitable
+    # target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELWITHDEBINFO>:/Ot>) # favor faster code over binary size
+    # target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELEASE>:/GL>)
+    # target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELEASE>:/LTCG>)
+    # target_compile_options(wildmeshing_toolkit PUBLIC $<$<CONFIG:RELEASE>:/Ot>)
     target_compile_definitions(wildmeshing_toolkit PUBLIC _SILENCE_ALL_MS_EXT_DEPRECATION_WARNINGS)
 endif()
 

cmake/recipes/libigl.cmake

@@ -12,6 +12,7 @@ CPMAddPackage(
     GIT_TAG 3ea7f9480967fcf6bf02ce9b993c0ea6d2fc45f6
     OPTIONS
         LIBIGL_PREDICATES ON
+        LIBIGL_COPYLEFT_TETGEN ON
 )
 # include(eigen)
 FetchContent_MakeAvailable(libigl)

components/image_simulation/wmtk/components/image_simulation/CMakeLists.txt

@@ -45,6 +45,7 @@ target_link_libraries(wmtk_${COMPONENT_NAME} PUBLIC
 	wmtk::sec_lib
 	VolumeRemesher::VolumeRemesher
 	jse::jse
+	igl_copyleft::tetgen
 )
 
 ######################

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

@@ -5,6 +5,7 @@
 #include <igl/is_vertex_manifold.h>
 #include <igl/remove_duplicate_vertices.h>
 #include <igl/remove_unreferenced.h>
+#include <igl/copyleft/tetgen/tetrahedralize.h>
 // igl must be included BEFORE VolumeRemesher
 #include <VolumeRemesher/embed.h>
 #include <VolumeRemesher/numerics.h>
@@ -784,6 +785,8 @@ void EmbedSurface::remove_duplicates(const double eps)
 
 bool EmbedSurface::embed_surface()
 {
+    logger().info("Embed with VolumeInsertion");
+
     std::shared_ptr<SampleEnvelope> ptr_env;
     {
         const auto v_simplified = V_surf_to_vector();
@@ -838,6 +841,142 @@ bool EmbedSurface::embed_surface()
     return all_rounded;
 }
 
+bool EmbedSurface::embed_surface_tetgen()
+{
+    logger().info("Embed with TetGen");
+
+    std::shared_ptr<SampleEnvelope> ptr_env;
+    {
+        const auto v_simplified = V_surf_to_vector();
+
+        std::vector<Eigen::Vector3i> tempF(m_F_surface.rows());
+        for (size_t i = 0; i < tempF.size(); ++i) {
+            tempF[i] = m_F_surface.row(i);
+        }
+        ptr_env = std::make_shared<SampleEnvelope>();
+        ptr_env->init(v_simplified, tempF, 0.5);
+    }
+
+    Eigen::MatrixXd V_in;
+    Eigen::MatrixXi F_in;
+    {
+        const Vector3d vertices_max = m_V_surface.colwise().maxCoeff();
+        const Vector3d vertices_min = m_V_surface.colwise().minCoeff();
+        const double diag = (vertices_max - vertices_min).norm();
+
+        // bbox
+        double delta = diag / 15.0;
+        const Vector3d box_min(
+            vertices_min[0] - delta,
+            vertices_min[1] - delta,
+            vertices_min[2] - delta);
+        const Vector3d box_max(
+            vertices_max[0] + delta,
+            vertices_max[1] + delta,
+            vertices_max[2] + delta);
+
+
+        // add corners of domain
+        std::vector<Vector3d> points(8);
+        for (int i = 0; i < 8; i++) {
+            Vector3d& p = points[i];
+            std::bitset<sizeof(int) * 8> a(i);
+            for (int j = 0; j < 3; j++) {
+                if (a.test(j)) {
+                    p[j] = box_max[j];
+                } else {
+                    p[j] = box_min[j];
+                }
+            }
+        }
+
+        const double voxel_resolution = diag / 20.0;
+        std::array<int, 3> N; // number of grid points per dimension
+        std::array<double, 3> h; // distance between grid points per dimension
+        for (int i = 0; i < 3; i++) {
+            const double D = box_max[i] - box_min[i];
+            N[i] = (D / voxel_resolution) + 1;
+            h[i] = D / N[i];
+        }
+
+        std::array<std::vector<double>, 3> ds;
+        for (int i = 0; i < 3; i++) {
+            ds[i].push_back(box_min[i]);
+            for (int j = 0; j < N[i] - 1; j++) {
+                ds[i].push_back(box_min[i] + h[i] * (j + 1));
+            }
+            ds[i].push_back(box_max[i]);
+        }
+
+        const double min_dis = voxel_resolution * voxel_resolution / 4;
+        //    double min_dis = state.target_edge_len * state.target_edge_len;//epsilon*2
+        for (int i = 0; i < ds[0].size(); i++) {
+            for (int j = 0; j < ds[1].size(); j++) {
+                for (int k = 0; k < ds[2].size(); k++) {
+                    if ((i == 0 || i == ds[0].size() - 1) && (j == 0 || j == ds[1].size() - 1) &&
+                        (k == 0 || k == ds[2].size() - 1)) {
+                        continue;
+                    }
+                    const Vector3d p(ds[0][i], ds[1][j], ds[2][k]);
+
+                    Eigen::Vector3d n;
+                    const double sqd = ptr_env->nearest_point(p, n);
+
+                    if (sqd < min_dis) {
+                        continue;
+                    }
+                    points.emplace_back(ds[0][i], ds[1][j], ds[2][k]);
+                }
+            }
+        }
+
+
+        V_in.resize(m_V_surface.rows() + points.size(), 3);
+        V_in.block(0, 0, m_V_surface.rows(), 3) = m_V_surface;
+        for (size_t i = 0; i < points.size(); ++i) {
+            V_in.row(m_V_surface.rows() + i) = points[i];
+        }
+        std::array<size_t, 8> p;
+        for (size_t i = 0; i < p.size(); ++i) {
+            p[i] = m_V_surface.rows() + i;
+        }
+
+
+        F_in = m_F_surface;
+        // F_in.resize(m_F_surface.rows() + 12, 3);
+        // F_in.block(0, 0, m_F_surface.rows(), 3) = m_F_surface;
+        // F_in.row(m_F_surface.rows() + 0) = Vector3i(p[0], p[1], p[2]);
+        // F_in.row(m_F_surface.rows() + 1) = Vector3i(p[1], p[3], p[2]);
+        // F_in.row(m_F_surface.rows() + 2) = Vector3i(p[0], p[5], p[1]);
+        // F_in.row(m_F_surface.rows() + 3) = Vector3i(p[0], p[4], p[5]);
+        // F_in.row(m_F_surface.rows() + 4) = Vector3i(p[0], p[6], p[4]);
+        // F_in.row(m_F_surface.rows() + 5) = Vector3i(p[0], p[2], p[6]);
+        // F_in.row(m_F_surface.rows() + 6) = Vector3i(p[7], p[3], p[2]);
+        // F_in.row(m_F_surface.rows() + 7) = Vector3i(p[7], p[2], p[6]);
+        // F_in.row(m_F_surface.rows() + 8) = Vector3i(p[7], p[6], p[4]);
+        // F_in.row(m_F_surface.rows() + 9) = Vector3i(p[7], p[4], p[5]);
+        // F_in.row(m_F_surface.rows() + 10) = Vector3i(p[7], p[5], p[1]);
+        // F_in.row(m_F_surface.rows() + 11) = Vector3i(p[7], p[1], p[3]);
+    }
+
+    MatrixXi TF;
+
+    int ret = igl::copyleft::tetgen::tetrahedralize(V_in, F_in, "pc", m_V_emb, m_T_emb, TF);
+    if (ret != 0) {
+        log_and_throw_error("TetGen returned with {}", ret);
+    }
+
+    m_V_emb_r.resizeLike(m_V_emb);
+    for (int i = 0; i < m_V_emb.rows(); ++i) {
+        m_V_emb_r.row(i) = to_rational(m_V_emb.row(i));
+    }
+
+    // add tags
+    tag_tets_from_images(m_img_datas, m_xyz2ijk, m_V_emb, m_T_emb, m_T_tags);
+
+    return true;
+}
+
 void EmbedSurface::consolidate()
 {
     std::map<size_t, size_t> old2new;

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

@@ -110,6 +110,8 @@ class EmbedSurface
 
     bool embed_surface();
 
+    bool embed_surface_tetgen();
+
     /**
      * @brief Remove unreferenced vertices.
      */

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

@@ -121,6 +121,9 @@ std::tuple<double, double> ImageSimulationMesh::local_operations(
     std::tuple<double, double> energy;
 
     auto sanity_checks = [this]() {
+        if (!m_params.perform_sanity_checks) {
+            return;
+        }
         logger().info("Perform sanity checks...");
         const auto faces = get_faces_by_condition([](auto& f) { return f.m_is_surface_fs; });
         for (const auto& verts : faces) {
@@ -694,9 +697,8 @@ std::vector<std::array<size_t, 3>> ImageSimulationMesh::get_faces_by_condition(
         if (cond(m_face_attribute[fid])) {
             auto tid = fid / 4, lid = fid % 4;
             auto verts = get_face_vertices(f);
-            res.emplace_back(
-                std::array<size_t, 3>{
-                    {verts[0].vid(*this), verts[1].vid(*this), verts[2].vid(*this)}});
+            res.emplace_back(std::array<size_t, 3>{
+                {verts[0].vid(*this), verts[1].vid(*this), verts[2].vid(*this)}});
         }
     }
     return res;
@@ -1133,6 +1135,9 @@ void ImageSimulationMesh::write_vtu(const std::string& path)
     Eigen::MatrixXd V(vs.size(), 3);
     Eigen::MatrixXi T(tets.size(), 4);
 
+    Eigen::VectorXd v_sizing_field(vert_capacity());
+    v_sizing_field.setZero();
+
     Eigen::MatrixXd parts(tets.size(), 1);
     std::vector<MatrixXd> tags(m_tags_count, MatrixXd(tets.size(), 1));
     Eigen::MatrixXd amips(tets.size(), 1);
@@ -1156,6 +1161,7 @@ void ImageSimulationMesh::write_vtu(const std::string& path)
     for (const Tuple& v : vs) {
         const size_t vid = v.vid(*this);
         V.row(vid) = m_vertex_attribute[vid].m_posf;
+        v_sizing_field[vid] = m_vertex_attribute[vid].m_sizing_scalar;
     }
 
     std::shared_ptr<paraviewo::ParaviewWriter> writer;
@@ -1166,6 +1172,7 @@ void ImageSimulationMesh::write_vtu(const std::string& path)
         writer->add_cell_field(fmt::format("tag_{}", j), tags[j]);
     }
     writer->add_cell_field("quality", amips);
+    writer->add_field("sizing_field", v_sizing_field);
     writer->write_mesh(path, V, T);
 }
 

components/image_simulation/wmtk/components/image_simulation/Parameters.h

@@ -22,6 +22,8 @@ struct Parameters
 
     double stop_energy = 10;
 
+    bool perform_sanity_checks = false;
+
     void init(const Vector3d& min_, const Vector3d& max_)
     {
         box_min = min_;

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

@@ -6,20 +6,10 @@
 
 #include <jse/jse.h>
 #include <wmtk/TetMesh.h>
-#include <wmtk/utils/Partitioning.h>
-#include <paraviewo/VTUWriter.hpp>
 #include <wmtk/Types.hpp>
-#include <wmtk/envelope/Envelope.hpp>
-#include <wmtk/utils/InsertTriangleUtils.hpp>
 #include <wmtk/utils/Logger.hpp>
-#include <wmtk/utils/ManifoldUtils.hpp>
-#include <wmtk/utils/Reader.hpp>
-#include <wmtk/utils/io.hpp>
-#include <wmtk/utils/partition_utils.hpp>
 #include <wmtk/utils/resolve_path.hpp>
 
-#include <sec/ShortestEdgeCollapse.h>
-
 #include "EmbedSurface.hpp"
 #include "ImageSimulationMesh.h"
 #include "Parameters.h"
@@ -164,7 +154,9 @@ void image_simulation(nlohmann::json json_params)
         V_envelope = image_mesh.V_surface();
         F_envelope = image_mesh.F_surface();
 
-        const bool all_rounded = image_mesh.embed_surface();
+        // const bool all_rounded = image_mesh.embed_surface();
+        const bool all_rounded = json_params["use_tetgen"] ? image_mesh.embed_surface_tetgen()
+                                                           : image_mesh.embed_surface();
         image_mesh.consolidate();
 
         if (write_vtu) {